examples/seismic/tutorials/04_dask_pickling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 04 - Full waveform inversion with Dask and Devito pickling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "Here, we revisit [04_dask.ipynb: Full Waveform Inversion with Devito and Dask](https://github.com/devitocodes/devito/blob/master/examples/seismic/tutorials/04_dask.ipynb), but with a twist: we now want to show that it is possible to use pickle to serialize (deserialize) a Devito object structure into (from) a byte stream. This is specially useful in our example as the geometry of all source experiments remains essentially the same; only the source location changes. In other words, we can convert a `solver` object (built on top of generic Devito objects) into a byte stream to store it. Later on, this byte stream can then be retrieved and de-serialized back to an instance of the original `solver` object by the dask workers, and then be populated with the correct geometry for the i-th source location. We can still benefit from the simplicity of the example and create **only one `solver`** object which can be used to both generate the observed data set and to compute the predicted data and gradient in the FWI process. Further examples of pickling can be found [here](https://github.com/devitocodes/devito/blob/master/tests/test_pickle.py). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tutorial roughly follows the structure of [04_dask.ipynb](https://github.com/devitocodes/devito/blob/master/examples/seismic/tutorials/04_dask.ipynb). Technical details about [Dask](https://dask.pydata.org/en/latest/#dask) and [scipy.optimize.minimize](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html) will therefore treated only superficially."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What is different from 04_dask.ipynb\n",
    "\n",
    "* **The big difference between [04_dask.ipynb](https://github.com/devitocodes/devito/blob/master/examples/seismic/tutorials/04_dask.ipynb) and this tutorial is that in the former is created a `solver` object for each source in both forward modeling and FWI gradient kernels. While here only one `solver` object is created and reused along all the optimization process. This is done through pickling and unpickling respectively.**\n",
    "\n",
    "\n",
    "* Another difference between the tutorials is that the in [04_dask.ipynb](https://github.com/devitocodes/devito/blob/master/examples/seismic/tutorials/04_dask.ipynb) is created a list with the observed shots, and then each observed shot record of the list is passed as parameter to a single-shot FWI objective function executed in parallel using the `submit()` method. Here, a single observed shot record along information of its source location is stored in a dictionary, which is saved into a pickle file. Later, dask workers retrieve the corresponding pickled data when computing the gradient for a single shot. The same applies for the `model` object in the optimization process. It is serialized each time the model's velocity is updated. Then, dask workers unpickle data from file back to `model` object.\n",
    "\n",
    "\n",
    "* Moreover, there is a difference in the way that the global functional-gradient is obtained. In [04_dask.ipynb](https://github.com/devitocodes/devito/blob/master/examples/seismic/tutorials/04_dask.ipynb) we had to wait for all computations to finish via `wait(futures)` and then we sum the function values and gradients from all workers. Here, it is defined a type `fg_pair` so that a reduce function `sum` can be used, such function takes all the futures given to it and after they are completed, combine them to get the estimate of the global functional-gradient. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## scipy.optimize.minimize \n",
    "\n",
    "As in [04_dask.ipynb](https://github.com/devitocodes/devito/blob/master/examples/seismic/tutorials/04_dask.ipynb), 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": [
    "## Dask\n",
    "\n",
    "[Dask](https://dask.pydata.org/en/latest/#dask) is task-based parallelization framework for Python. It allows us to distribute our work among a collection of workers controlled by a central scheduler. Dask is [well-documented](https://docs.dask.org/en/latest/), flexible, an currently under active development.\n",
    "\n",
    "In the same way as in [04_dask.ipynb](https://github.com/devitocodes/devito/blob/master/examples/seismic/tutorials/04_dask.ipynb), 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": [
    "## Forward modeling\n",
    "\n",
    "We define the functions used for the forward modeling, as well as the other functions used in constructing and deconstructing Python/Devito objects to/from binary data as follows:"
   ]
  },
  {
   "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",
    "         '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",
    "         'nbl': 40}                # nbl 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, AcquisitionGeometry, Receiver\n",
    "from examples.seismic.acoustic import AcousticWaveSolver\n",
    "from examples.seismic import AcquisitionGeometry\n",
    "\n",
    "# Import convenience function for plotting results\n",
    "from examples.seismic import plot_image\n",
    "from examples.seismic import plot_shotrecord\n",
    "\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_circle=3.0, vp_background=2.5, \n",
    "                      origin=param['origin'], shape=param['shape'],\n",
    "                      spacing=param['spacing'], nbl=param['nbl'])\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_circle=2.5, vp_background=2.5, \n",
    "                      origin=param['origin'], shape=param['shape'],\n",
    "                      spacing=param['spacing'], nbl=param['nbl'],\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",
    "    v_curr = 1.0/np.sqrt(x.reshape(model.shape))\n",
    "    \n",
    "    if astype:\n",
    "        model.update('vp', v_curr.astype(astype).reshape(model.shape))\n",
    "    else:\n",
    "        model.update('vp', v_curr.reshape(model.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['geometry'], pkl['rec'].resample(dt)\n",
    "\n",
    "def dump_shot_data(shot_id, rec, geometry):\n",
    "    ''' Dump shot data to disk.\n",
    "    '''\n",
    "    pickle.dump({'rec':rec, 'geometry': geometry}, 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",
    "    # Reconstruct objects\n",
    "    with open(\"arguments.pkl\", \"rb\") as cp_file:\n",
    "        cp = pickle.load(cp_file)\n",
    "        \n",
    "    solver = cp['solver']\n",
    "\n",
    "    # source position changes according to the index\n",
    "    shot_id=param['shot_id']\n",
    "    \n",
    "    solver.geometry.src_positions[0,:]=[20, shot_id*1000./(param['nshots']-1)]\n",
    "    true_d = solver.forward()[0]\n",
    "    dump_shot_data(shot_id, true_d.resample(4.0), solver.geometry.src_positions)\n",
    "\n",
    "def generate_shotdata(solver):\n",
    "    # Pick devito objects (save on disk)\n",
    "    cp = {'solver': solver}\n",
    "    with open(\"arguments.pkl\", \"wb\") as cp_file:\n",
    "        pickle.dump(cp, cp_file) \n",
    "\n",
    "    work = [dict(param) for i in range(param['nshots'])]\n",
    "    # synthetic data is generated here twice: serial(loop below) and parallel (via dask map functionality) \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, We pass our function and the dictionary to the map() function of the client\n",
    "    # This returns a list of futures that represents each task\n",
    "    futures = c.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": [
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `Forward` run in 0.02 s\n",
      "Operator `Forward` run in 0.03 s\n",
      "Operator `Forward` run in 0.02 s\n",
      "Operator `Forward` run in 0.02 s\n",
      "Operator `Forward` run in 0.03 s\n"
     ]
    }
   ],
   "source": [
    "#NBVAL_IGNORE_OUTPUT\n",
    "from examples.seismic import plot_shotrecord\n",
    "\n",
    "# Client setup\n",
    "cluster = LocalCluster(n_workers=2, death_timeout=600)\n",
    "c = Client(cluster)\n",
    "\n",
    "# Generate shot data.\n",
    "true_model = get_true_model()\n",
    "# Source coords definition\n",
    "src_coordinates = np.empty((1, len(param['shape'])))\n",
    "# Number of receiver locations per shot.\n",
    "nreceivers = 101\n",
    "# Set up receiver data and geometry.\n",
    "rec_coordinates = np.empty((nreceivers, len(param['shape'])))\n",
    "rec_coordinates[:, 1] = np.linspace(param['spacing'][0], true_model.domain_size[0] - param['spacing'][0], num=nreceivers)\n",
    "rec_coordinates[:, 0] = 980. # 20m from the right end\n",
    "# Geometry \n",
    "geometry = AcquisitionGeometry(true_model, rec_coordinates, src_coordinates,\n",
    "                               param['t0'], param['tn'], src_type='Ricker',\n",
    "                               f0=param['f0'])\n",
    "# Set up solver\n",
    "solver = AcousticWaveSolver(true_model, geometry, space_order=4)\n",
    "generate_shotdata(solver)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dask specifics\n",
    "\n",
    "Previously in [03_fwi.ipynb](https://github.com/devitocodes/devito/blob/master/examples/seismic/tutorials/03_fwi.ipynb), 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 `__radd__` 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": [
    "#NBVAL_IGNORE_OUTPUT\n",
    "from devito import Function\n",
    "\n",
    "# Create FWI gradient kernel for a single shot\n",
    "def fwi_gradient_i(param):\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 reliably 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",
    "    nbl = model0.nbl\n",
    "\n",
    "    # Get src_position and data\n",
    "    src_positions, rec = load_shot_data(param['shot_id'], dt)\n",
    "\n",
    "    # Set up solver -- load the solver used above in the generation of the syntethic data.    \n",
    "    with open(\"arguments.pkl\", \"rb\") as cp_file:\n",
    "        cp = pickle.load(cp_file)\n",
    "    solver = cp['solver']\n",
    "    \n",
    "    # Set attributes to solver\n",
    "    solver.geometry.src_positions=src_positions\n",
    "    solver.geometry.resample(dt)\n",
    "\n",
    "    # Compute simulated data and full forward wavefield u0\n",
    "    d, u0 = solver.forward(vp=model0.vp, dt=dt, save=True)[0:2]\n",
    "        \n",
    "    # Compute the data misfit (residual) and objective function\n",
    "    residual = Receiver(name='rec', grid=model0.grid,\n",
    "                        time_range=solver.geometry.time_axis,\n",
    "                        coordinates=solver.geometry.rec_positions)\n",
    "\n",
    "    #residual.data[:] = d.data[:residual.shape[0], :] - rec.data[:residual.shape[0], :]\n",
    "    residual.data[:] = d.data[:] - rec.data[0:d.data.shape[0], :]\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, vp=model0.vp, dt=dt, 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[:])[nbl:-nbl, nbl:-nbl]    \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 = c.map(fwi_gradient_i, work, retries=1)\n",
    "    \n",
    "    # Perform reduction.\n",
    "    fg = c.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",
    "# 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 monitor the true relative solution error.\n",
    "relative_error = []\n",
    "def fwi_callbacks(x):    \n",
    "    # Calculate true relative error\n",
    "    true_vp = get_true_model().vp.data[param['nbl']:-param['nbl'], param['nbl']:-param['nbl']]\n",
    "    true_m = 1.0 / (true_vp.reshape(-1).astype(np.float64))**2\n",
    "    relative_error.append(np.linalg.norm((x-true_m)/true_m))\n",
    "\n",
    "# FWI with L-BFGS\n",
    "ftol = 0.1\n",
    "maxiter = 5\n",
    "\n",
    "def fwi(model, param, ftol=ftol, maxiter=maxiter):\n",
    "    # Initial guess\n",
    "    v0 = model.vp.data[param['nbl']:-param['nbl'], param['nbl']:-param['nbl']]\n",
    "    m0 = 1.0 / (v0.reshape(-1).astype(np.float64))**2\n",
    "    \n",
    "    # Define bounding box constraints on the solution.\n",
    "    vmin = 1.4    # do not allow velocities slower than water\n",
    "    vmax = 4.0\n",
    "    bounds = [(1.0/vmax**2, 1.0/vmin**2) for _ in range(np.prod(model.shape))]    # in [s^2/km^2]\n",
    "    \n",
    "    result = optimize.minimize(fwi_gradient,\n",
    "                               m0, args=(param, ), method='L-BFGS-B', jac=True,\n",
    "                               bounds=bounds, 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": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n",
      "Operator `initdamp` run in 0.01 s\n",
      "Operator `padfunc` run in 0.01 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      fun: 149.96025455476786\n",
      " hess_inv: <10201x10201 LbfgsInvHessProduct with dtype=float64>\n",
      "      jac: array([39.87949753, 37.44005203, 30.17471313, ...,  0.43074584,\n",
      "       -0.77040124, -2.02412486])\n",
      "  message: b'STOP: TOTAL NO. of ITERATIONS REACHED LIMIT'\n",
      "     nfev: 8\n",
      "      nit: 5\n",
      "     njev: 8\n",
      "   status: 1\n",
      "  success: False\n",
      "        x: array([0.15686476, 0.15623097, 0.15602988, ..., 0.15440663, 0.1551215 ,\n",
      "       0.17023674])\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": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAFwCAYAAADKcW5qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de5QdZZnv8e8vIQnBgCFivAQFIi4FDqNw4gwjOCHxYAIrCTjc5TYwEkVkGF3IMURd0MKEEQ86HhHJ6ByZk3g5wDAQhmtMwBkPGYkoQsR4YBJCqwyEhFsgSJLn/PFWS1HZ3bu69+7u9Nu/z1q1Kl31Vu23NnSePG+9F0UEZmZmuRkx2BUwMzPrDw5wZmaWJQc4MzPLkgOcmZllyQHOzMyy5ABnZmZZcoAzM7Nek3ScpBskPSbpJUmrJS2QtGuNa98u6VpJ64prfy3pUkmva1D2bEm/kvRy8Rkfr1vHnXr7UGZmZsAFwDrgIqATOAi4GJgm6f0Rsa3RRUUQWwqMAj5f3ON9wCXAO4ETS2XPBq4BFhTXfBD4hiRFxNXNKqjBGugt6W3AV4AjAJEq/9cRsW5QKmRmZrVJemNEPFU5djpwLfDBiFjWzXUfAu4AZkTEnaXjl5OC5m4R8aKknYDfArdFxBmlcv8AzAHeEhGv9FTHQWmilLQLsAx4N3AGcBopci9vlKKamdmOpRrcCvcV+0k9XDq62D9XOf4MKSap+PlPgTcCiyrl/jfwBuCwZnUcrHdwZwOTgWMi4p8j4iZSRN4L+Ngg1cnMzFoztdg/3EOZpcD/A/5W0v6SxkmaDpwPfDMiNhXlDij2D1WuX1Xs929WmcEKcHOAFRHxSNeBiFgD/Bg4epDqZGZmfSRpEtABLI2Ild2Vi4jNpOxrBClYPQ/8ELgF+GSp6IRiv7Fyiw2V890arE4mBwA3NTi+Cji+2cWjR4+NnXdp2lHHzGzY2Pzi8/z+9y8JYOb08bH+6S0t3e+nv9i0CthcOrQwIhY2KitpHOnv9C3AmT3dV9LOwA+AiaTXU+uAPwa+UFx/TksVLxmsADeB7aMypMi8e6MLJM0F5gKMGTuO9/3ZCf1XOzOzIea+H/2fP/x5/dOvcN/t/6Wl+4146082R8SUZuUkjQWWkF47TY2IziaX/CVwOLBvRDxaHPuRpGeBhZK+GREP8GqM2B34Xen6rsxtA00MmXFwEbEwIqZExJTRo8cOdnXMzIY9SaOA64EpwFER8WCNyw4ENpaCW5efFPv9in3Xu7YDKuW63r39stkHDVaA20jjTK27zM7MzHojorWtCUkjgMXAdFKHwRU1a/YEsLukfSvH/6TY/6bY3wusB06plDuVlL39uNkHDVYT5Sq2j8qQInPTqGxmZk30/xjnq0h9Ji4DNkk6pHSuMyI6Je0FPAp0RERHce47wKeBWyVdRnoHN4U06PunFIErIl6R9HnSwO7fkHpfTgfOAs6LiN83q+BgZXA3A4dImtx1QNLewKHFOTMz66sWk7easfHIYj+flG2Vt48W5wSMpBRrImItcAjwc+BS4FbS0LGFwBHlGVAi4pukTicnkAaHnwx8MiKuqlPBwcrg/p7UHfQmSZ8DAvgi8DhpWhYzM9uBRcTeNcqs5dWB2+XjvyQFrTqfcw19jAuDksEVA/mmA78mjUpfDKwBpkfEC4NRJzOznESMaGnLwaBNtlzMOXnsYH2+mVmuAmUTpFrh1QTMzDK0zQHOAc7MLD/O4GAIDfQ2MzPrDWdwZma5CZzB4QBnZpadACK2650/7DjAmZllx+/gwAHOzCxL25zBuZOJmZnlyRmcmVmG3ETpAGdmlp00YbKbKB3gzMyyIwc4/A7OzMwy5QzOzCxDfgfnAGdmliU3UTrAmZllxzOZJA5wZma5CXmgN+5kYmZmmXIGZ2aWITdROsCZmWUp3EDnAGdmlht3Mkkc4MzMsuOZTMCdTMzMLFPO4MzMMuQMzgHOzCw/Xk0AcIAzM8uSA5zfwZmZWaYc4MzMMhOkqbpa2ZqRdJykGyQ9JuklSaslLZC0a5PrLpYU3WybK2XXdlPumDrfg5sozcwyNABNlBcA64CLgE7gIOBiYJqk90fEtm6u+xZwe+XY64pjNzcof0dx37LVdSroAGdmlqEBCHCzI+Kp0s/3SNoAXAscDixrXK/oJAXEP5B0GikeXdvgkvURsaIvFXSAMzPLzQD0oqwEty73FftJvbzdGcB/krK1tvE7ODMza5epxf7huhdIehswDVgcEVsaFJkt6UVJL0taUff9GziDMzPLTpvmotxD0srSzwsjYmF3hSVNAjqApRGxsrtyDZxKSrYaNU8uIWWFa4A3AZ8EbpR0WkQsanZjBzgzswxFtHyL9RExpU5BSeOAm4AtwJm9/JzTgZ9FxC+qJyLivMrn3AisABYATQOcmyjNzLIjosWt9idJY0mZ1mRgRtGJpO61fwy8m8bZ23YiYitwHbCnpLc0K+8MzswsQwMxk4mkUcD1wBTgiIh4sJe3OAN4BfhuHz6+aY7qAGdmZr0maQSwGJgOzOptV35Jo4GTgNu66ZHZ6JqdgBOBdRHxRLPyDnBmZhkagAzuKuB44DJgk6RDSuc6I6JT0l7Ao0BHRHRUrp8FTKCb5klJJwNHA7cCj5M6mZwLHAycXKeCDnBmZpmJgVlN4MhiP7/Yyi4hzT4iYCSN+3ucAWwAbunm/muAicAVpEC4CVgJzIyIWuPlHODMzDK0rfVelD2KiL1rlFkLjXusRMTRTa5dQWr+7DP3ojQzsyw5gzMzy5DXg3OAMzPLkBzgcIAzM8tOm6bqGvIc4MzMchNtmapryHMnEzMzy5IzODOzDPVmPslcOcCZmWXI7+Ac4MzMsuQA5wBnZpad1ItysGsx+NzJxMzMsuQMzswsOx7oDW3O4CQdJ+kGSY9JeknSakkLJO1aKbe7pG9JWi9pk6Slkg5sZ13MzIat6FpRoO9bDtrdRHkBsBW4CJgJXA2cA9xVLI6HJJGWN58JnAccC4wClkvas831MTMbliLU0paDdjdRzq6szHqPpA2kBe0OB5YBc4BDgekRsRxA0r2ktX8uBP6qzXUyM7NhqK0ZXDfLjt9X7CcV+znAb7uCW3Hds6Ssrsf1gczMrB5ncAPTi3JqsX+42B8APNSg3Crg7ZLGDUCdzMyyFcC2Frcc9GsvSkmTgA5gaUSsLA5PANY2KL6h2O8OvNCf9TIzy10uWVgr+i3AFZnYTcAW4Mw23G8uMBdgzFgneWZmPXGA66cmSkljSe/UJgMzIqKzdHojKUurmlA6v52IWBgRUyJiyujRY9taXzMzy0/bMzhJo4DrgSnAERHxYKXIKuBDDS7dH1gXEW6eNDNrRUZj2VrR7oHeI4DFwHTgmIhY0aDYzcAkSVNL1+0GzC7OmZlZC4LWelDm0rzZ7gzuKuB44DJgk6RDSuc6i6bKm4F7gUWSPkNqkpwHCPhSm+tjZjYsOYNr/zu4I4v9fFIQK28fBYiIbcAs4C7gG8CNpNlPpkXE422uj5nZsOQMrs0ZXETsXbPcBuCsYjMzM2s7ryZgZpYht1A6wJmZZSmXZsZWOMCZmWUmpyVvWuEVvc3MrNfqrv/Z4LqLJUU32+ZK2RGS5klaK2mzpAckHVu3js7gzMwyNABNlBcA60jrf3YCBwEXA9Mkvb/oMd/It4DbK8deVxyrjoX+YvE584GfAicB10maFRG3NqugA5yZWYYGoImyzvqfDeoVnaSA+AeSTiPFo2tLxyaSgtvlEfHl4vBySfsClwNNA5ybKM3MMtTf4+Bqrv9Z1xnAfwJ3lI7NAEYDiyplFwEHStqn2U0d4MzMMhQtbn1UXf+zKUlvA6YBiyNiS+nUAcDLwCOVS1YV+/2b3dtNlGZm1sgeklaWfl4YEQu7K9zN+p91nEpKtq6tHJ8APBOxXWPrhtL5HjnAmZllpmuy5Ratj4gpdQq2uP7n6cDPIuIXvbyuKQc4M7PcDOA4uMr6n1Mr6382u/aPgXcDf93g9EZgvCRVsriuzG1Dg2tew+/gzMwy1DXYu69bHZX1P49qsP5nM2cArwDfbXBuFTAGeEfleNe7t182u7kDnJmZ9VrN9T97un40aVzbbd30yLydFPxOqRw/FXgoItY0+ww3UZqZZWgABno3Xf9T0l7Ao0BHRHRUrp9Fam6sdi4BICKelHQlME/S88D9wImkgDqnTgUd4MzMMjQA7+DK63/Or5y7hDSriYCRNG4tPIP0Hu2WHj5jPvACcD7wZmA1cEJE9HTNHzjAmZllJo1l698Mrs76nxGxFhpXJCKOrnH9VuDSYus1Bzgzswx5NQF3MjEzs0w5gzMzy43XgwMc4MzMsuQVvR3gzMyyE8A2Z3AOcGZmOervXpRDgTuZmJlZlpzBmZllyJ1MHODMzLLkAOcAZ2aWn2jLenBDnt/BmZlZlpzBmZllJnATJTjAmZllyQHOAc7MLEseB+cAZ2aWJWdw7mRiZmaZcgZnZpYbryYAOMCZmWXHvSgTBzgzswx5oLcDnJlZlpzAuZOJmZllyhmcmVmG/A7OAc7MLEsOcA5wZmbZiXAnE/A7ODMzy5QzODOzDLmF0gHOzCxLfgfnAGe2nbsWrhyQzzli7pQB+Rwbnhzg/A7OzCw7aaoutbQ1I+k4STdIekzSS5JWS1ogadc6dZS0n6TrJK0vXX9+pcxaSdFgO6bOZziDMzOzvrgAWAdcBHQCBwEXA9MkvT8itnV3oaQpwDLgbuCjwLPAO4FxDYrfUdy3bHWdCjrA2bAxUE2PddWtj5syrS8GoIlydkQ8Vfr5HkkbgGuBw0kBbDuSRgD/CPwwIj5cOrW8m89ZHxEr+lJBBzgzswz1d3yrBLcu9xX7ST1cejiwH/Cxdtepyu/gzMxyE12Dvfu+9dHUYv9wD2UOK/Y7S1oh6RVJT0r6mqSxDcrPlvSipJeL8rXev4EzOMvcjtYs2RflZ3BzpdXVhibKPSSVf4EWRsTC7gpLmgR0AEsjoqdfvLcW+x8AXwc+C0wprn0bUG62XELKCtcAbwI+Cdwo6bSIWNTsARzgzMyskfURUetfVJLGATcBW4AzmxTvajlcFBFfKP58t6SRwOWS9ouIhwEi4rzK59wIrAAWAE0DnJsozcwyE7Q2RKA381gWzYpLgMnAjIjobHLJ08X+rsrxO4v9Qd0+V8RW4DpgT0lvaVa3fs3gJN0OzAAui4jPlY7vDlwBHAOMBe4FPhURD/ZnfSxfOTRF1tHdc7rp0qoGYpy3pFHA9aQmxiNq/h2+qsn5bocXVDR9xH7L4CSdDLynwXGRov1M4DzgWGAUsFzSnv1VHzOz4aS/O5kU3f0XA9OBY3rRlf824GVS8lM2s9h3+69VSTsBJwLrIuKJZh/ULxlckaF9BfgU8N3K6TnAocD0iFhelL+X9BLxQuCv+qNOZmbWVlcBxwOXAZskHVI61xkRnZL2Ah4FOiKiAyAinpa0APi8pOdI4+WmAF8Aro2IR+APSdLRwK3A46ROJucCBwMn16lgfzVR/i3wUER8T1KjAPfbruAGEBHPSlpCehgHOGtquDRJ1lX9PtxkaQMw0PvIYj+/2MouIc0+ImAk27cWdgDPA58gzYjyO9Jrqy+WyqwBJhbHJwCbSNndzIi4o04F2x7gJB0GnE6D5snCAcBDDY6vAk6XNC4iXmh3vczMho3WxrLV+4iIvWuUWUsKctXjAVxZbN1du4LU/NlnbQ1wkkYD1wBfjoju5gqbAKxtcHxDsd8d2C7ASZoLzAUYM7bRdGVmZgZdky0Pdi0GX7szuAtJvSIva/N9KQYYLgTYbfxE/6czM+tBbJ84DTttC3CS3k5qh/0oMEbSmNLpMZLGk9pcN5KytKoJxX5ju+pkefF7t/o8+4lZe4cJTAZ2Jo0u31jaIL1E3AgcSHrXdkCD6/cndf30+zczsxYN0lyUO5R2NlH+HJjW4PhyUtD7NvAIcDNwpqSpEXEPgKTdgNlsP6TAzMz6IpMg1Yq2BbiIeIa0eN1rpHHdPBYRdxc/30yauWSRpM+QMrt5pJ42X2pXfczMhrNcsrBWDPhclMUqr7NI85B9A7gR2ApMi4jHB7o+ZmY5iha3HPT7agLRYNbOiNgAnFVsZmZmbeflcmyH5p6TrXOPyuEnp44irXCAMzPLkAOcA5yZWZYc4BzgbAfkZsn+4+ZKG04c4MzMMuQEzgHOzCxLbqJ0gDMzy45XE0gc4MzMcuNhAsAgzGRiZmY2EJzBmZllyAmcA5yZWZbcROkAZ2aWpQbTAA87DnBmZplxL8rEAc52CJ69ZOBVv3PPbGK5cYAzM8uQEzgHODOzLLmJ0gHOzCw/HugNeKC3mZllyhmcmVmGwimcMzgzs9xEG7ZmJB0n6QZJj0l6SdJqSQsk7VqnjpL2k3SdpPWl68+vlBkhaZ6ktZI2S3pA0rG1vgQc4MzMshTR2lbDBcBW4CJgJnA1cA5wl6QeY4ukKcC/A2OAjwJHAf8DGFkp+kXgYuDrwJHACuA6SUfVqaCbKM3MMjQALZSzI+Kp0s/3SNoAXAscDixrdFER/P4R+GFEfLh0anml3ERSEL08Ir7cVUbSvsDlwK3NKugMzszMeq0S3LrcV+wn9XDp4cB+wJVNPmIGMBpYVDm+CDhQ0j7N6ugAZ2aWo/5+CdfY1GL/cA9lDiv2O0taIekVSU9K+pqksaVyBwAvA49Url9V7PdvVhkHODOzDA10fJM0CegAlkZET3PvvbXY/wC4EzgC+BLpXdx3S+UmAM/E9t1BN5TO98jv4MzMctOegd57SCoHqoURsbBRQUnjgJuALcCZTe7blVgtiogvFH++W9JI4HJJ+0VETxlgbQ5wZmaZadNqAusjoukM3EWz4hJgMjA1IjqbXPJ0sb+rcvxOUueRg0hNnBuB8ZJUyeK6MrcNNOEmSjMz6xNJo4DrgSnAURHxYI3LVjU5v61Ubgzwjsr5rndvv2z2QQ5wZmYZGoCB3iOAxcB04JiIWFGzareROo/MqByfWey7mkVvB14BTqmUOxV4KCLWNPsgN1GamWVoAKbqugo4HrgM2CTpkNK5zojolLQX8CjQEREdRb2elrQA+Lyk50jj5aYAXwCujYhHinJPSroSmCfpeeB+4ERSQJ1Tp4IOcGZmGRqAgd5HFvv5xVZ2CWkGEpFmJ6m2FnYAzwOfIA3m/h1wBWnmkrL5wAvA+cCbgdXACRFxS50KOsCZmVmvRcTeNcqsJQW56vEgDfTucbB3RGwFLi22XnOAMzPLkBcTcIAzM8tQa9OR5MIBzswsM71YESBrDnBmZhnygqceB2dmZplyBmdmliEncA5wZmZZcnxzgLMdxBFzX53T9a6FPa20Ye1S/s4tP34H53dwZmaWKWdwZmY5cgLnAGdmlhsP804c4MzMchN+BwcOcGZmWXJ8cycTMzPLlDM4M7MMOYNzgDMzy1K4m4kDnJlZbgJncOAAZzsgz2rSfzx7yTDiAOdOJmZmlidncGZmGXIC108BTtJRwGeBg4FtwK+BCyNiWXF+d+AK4BhgLHAv8KmIeLA/6mNDl5srW+dmyWEowgO96YcmSkkfA24Cfgp8GDgeuA7YpTgvYAkwEzgPOBYYBSyXtGe762NmNhxFtLbloK0ZnKS9ga8Cn4mIr5ZO3VH68xzgUGB6RCwvrrsXWANcCPxVO+tkZmbDU7szuLNITZLf7KHMHOC3XcENICKeJWV1R7e5PmZmw5IzuPYHuMOAXwEnSXpU0hZJj0g6t1TmAOChBteuAt4uaVyb62RmNqxEG7YctLuTyVuL7QrgIuBR0ju4r0vaKSL+DpgArG1w7YZivzvwQvWkpLnAXIAxYx0Dzcx64k4m7Q9wI4Bdgb+IiH8qji0r3s3Nk/S1vt44IhYCCwF2Gz/R/+WGIfeorM89J83xrf1NlE8X+7sqx+8E3gS8BdhIytKqJhT7jW2uk5mZDUPtDnCrmpzfVpQ5oMG5/YF1EbFd86SZme1YJB0n6QZJj0l6SdJqSQsk7Vrj2uhme2+l3Npuyh1Tp47tDnA3FvsZleMzgc6IeAK4GZgkaWrXSUm7AbOLc2Zm1ooWe1DWbN68ANhK6m8xE7gaOAe4S1Kd2PId4E8r268blLujQbl76lSw3e/gbgWWA9dI2gP4D1Inkw8BZxZlbibNXLJI0mdITZLzAAFfanN9LFPVd0zD/Z2c37lZ1QC8g5sdEU+Vfr5H0gbgWuBwYFmT638TEStqfM76muW209YAFxFdqeMC4BLSu7ZfAadExHeLMtskzQK+DHwD2JkU8KZFxOPtrI+Z2XCUlsvp3whXCW5d7iv2k/r1w2tq+1RdEfFcRJwbEW+KiNER8Uddwa1UZkNEnBUREyJil4j4YEQ80O66mJnZgOp69fRwjbLnSHpZ0ouSlkn6QDflZhdlXpa0ou77N/BqApaJ7procmu6dFOk1dWG/G0PSeVfoIXFcK2GJE0COoClEdHsF28RcAvwW2Av4DOkIWVHRMTdpXJLSFnhGlJP/E8CN0o6LSIWNXsABzgzswy1oYVyfUTU+hdVMQPVTcAWXu1v0a2IOK30479Kuok0w9WlpBmxusqdV/mcG4EVpNdgTQOcFzw1M8vQQM1FKWksKdOaDMyIiM7e1zWeB/4FeF+TcltJq9PsKektze7rDM6ylsPsJ26WtB2VpFHA9cAU4Ig2rOnZm7yzaVkHODOzDPX3KIFirNtiYDowq69d+Yt77QbMAn7SpNxOwImkSUGeaHZfBzgzs9ykcQL9/SlXkcY5XwZsknRI6VxnRHRK2os06X5HRHQASLoAeBdpzHRXJ5MLgDcDp3TdQNLJpCXUbgUeJ3UyORc4GDi5TgUd4GzYqNvUN1BNmW56tP4yMPGNI4v9/GIruwS4mDSBx0he299jNfDhYns98BzwY+AvI6Kcwa0BJpJWp5kAbAJWAjMjoryIdrcc4MzMMtTf8S0i9q5RZi0pyJWPLSF1Sml27QpS82efuRelmZllyRmcmVmGvB6cA5zZdvxuzHLgAOcAZ2aWJQc4Bzgzs+wEEP3ezWTH504mZmaWJWdwZma56eV8krlygDMzy5ADnAOcmVmWHN/8Ds7MzDLlDM7MLEdO4RzgzMxy5HdwDnBmZtkJnMCBA5yZWZacwbmTiZmZZcoZnJlZbjzQG3CAMzPLUBCOcA5wZma5cSeTxAHOzCxDTuDcycTMzDLlDM7MLEPO4BzgzMyy5PjmAGdmliVncH4HZ2ZmmXIGZ2aWmfBAb8ABzswsS45vDnBmZllyBucAZ2aWJQc4dzIxM7M+kHScpBskPSbpJUmrJS2QtGuNa6Ob7b2VciMkzZO0VtJmSQ9IOrZuHZ3BmZllaAASuAuAdcBFQCdwEHAxME3S+yNiW5PrvwNcUzn268rPXyw+Zz7wU+Ak4DpJsyLi1mYVdIAzM8tMMCBNlLMj4qnSz/dI2gBcCxwOLGty/W8iYkV3JyVNJAW3yyPiy8Xh5ZL2BS4HmgY4N1GameUmXh0q0Net6Ue8Nrh1ua/YT2rDU8wARgOLKscXAQdK2qfZDRzgzMwy1N8BrhtTi/3DNcqeI+llSS9KWibpA5XzBwAvA49Ujq8q9vs3+wAHODMza2QPSStL29yeCkuaBHQASyNiZZN7LwI+Afw3YC7wBmCZpMNLZSYAz8T2K7duKJ3vkd/BmZllqA2v4NZHxJQ6BSWNA24CtgBnNisfEaeVfvxXSTcBDwGXAof1oa4NOYMzM8vQQDVRShoLLAEmAzMiorP3dY3ngX8B3lc6vBEYL0mV4l2Z2waacIAzM8tMVy/K/g5wkkYB1wNTgKMi4sE2VL3LKmAM8I5Kma53b79sdjMHODMz6zVJI4DFwHTgmJ66/Ne4127ALOAnpcO3A68Ap1SKnwo8FBFrmt3X7+DMzDI0AAO9rwKOBy4DNkk6pHSuMyI6Je0FPAp0REQHgKQLgHcBy4HfAnuRxru9mVIwi4gnJV0JzJP0PHA/cCIpoM6pU0EHODOzDA3AQO8ji/38Yiu7hDSriYCRvLa1cDXw4WJ7PfAc8GPgLyPiJ6+9DfOBF4DzSQFwNXBCRNxSp4IOcGZmuRmA9eAiYu8aZdaSglz52BJSp5Q6n7GV1LPy0t7X0AHOzCxLXkzAnUzMzCxTbQ1wkg6VdKekJyU9L+l+SWdVyuws6QpJvyuWWLhX0p+1sx5mZsPZQA0T2NG1LcBJ+iNgKTAKOBv4c9LEm9+WdE6p6LeL818gdQv9HXBHdR0gMzPrOwe49r6DO4nUW2Z2RLxQHLurCHynA1dLeg/wEeCsiPhfAJLuIQ3o66Bm108zM+tZJjGqJe1sohxNGpT3UuX4s6XPmVOU+UHXyYjYAnwfmCFpTBvrY2Y2bDmDa2+A+06x/5qkt0oaL+ls4IPAV4pzBwBrIuLFyrWrSAFy3zbWx8zMhrG2NVFGxEPFUgc3kpZBgJStfTwivl/8PIE0gWZV7eUPzMysZwFsyyQLa0XbApykdwI3kLKxj5OaKo8Gvilpc0QsbvH+c0nrBjFm7LgWa2tmlrGMmhlb0c5OJn9DythmRcQrxbEfSnoD8HeSvkfK3vZqcG3T5Q8iYiGwEGC38RP9n87MrAf+S7K97+AOBB4oBbcuPyGt1jqRlN3tI2mXSpn9gd+z/dLkZmZmfdLOAPcE8F5JoyvH/wTYTMrOlpDGyR3fdVLSTqQZou+MiJfbWB8zs2FKRLS25aCdTZRfB64Dlkj6Bukd3BzgZOArEfF74GeSfgB8tVgobw1wDrAP26/5Y2ZmfeR3cO3tRXm9pKOA/w58C9iZtA7QucA1paJnktYPuhQYDzwAzIyI+9tVFzOz4SzwOzho82oCEXEbcFuTMi8Bny42MzPrBx4m4NUEzMwsU14PzswsQ34H5wBnZpYlxzcHODOz7ET4HRz4HZyZmWXKGZyZWYb8Ds4BzswsS0Ees5G0wgHOzCwzXi4ncYAzM8uQmyjdycTMzDLlDM7MLENuonSAMzPLkuObA5yZWXY80DvxOzgzswz194Knko6TdIOkxyS9JGm1pAWSdu1NPSV9VlJI+rcG59YW56rbMXXu7QzOzMz64gJgHXAR0AkcBFwMTJP0/ojY1uwGkiYDn/vFSAEAAAvqSURBVAOe7KHYHcV9y1bXqaADnJlZhppGl9bNjoinSj/fI2kDcC1wOLCsxj2uBhYD76L7eLQ+Ilb0pYJuojQzy0zXQO9Wtqaf8drg1uW+Yj+p2fWSPgIcDMzrxaP1igOcmVmG+jvAdWNqsX+4p0KSdge+AlwYERua3HO2pBclvSxpRd33b+AAZ2Zmje0haWVpm9tTYUmTgA5gaUSsbHLvK4BfA99pUm4JcB4wAzgF2AzcKOnUOg/gd3BmZplp01yU6yNiSp2CksYBNwFbgDOblP0AcDpwcETPE4pFxHmVa28EVgALgEXN6uUMzswsQ9tQS1tdksaSMq3JwIyI6GxyyTXAt4FOSeMljSclWyOLn8d0d2FEbAWuA/aU9JZmdXMGZ2aWmwEa6C1pFHA9MAU4IiIerHHZfsX28QbnNgKfAr5a4z5Nn9ABzswsMwOxXI6kEaQu/tOBWb3oyj+twbGvAiNJ79se6eEzdwJOBNZFxBPNPsgBzszM+uIq4HjgMmCTpENK5zojolPSXsCjQEdEdABExN3VG0l6BtipfE7SycDRwK3A48CbgHNJQwtOrlNBBzgzswxt7f+POLLYzy+2sktIs4+IlJn1pb/HGmAiqcflBGATsBKYGRF31LmBA5yZWWYC2NrPTZQRsXeNMmuheY+ViDi8wbEVpObPPnOAMzPLTABbvJqAA5yZWY629qKrf648Ds7MzLLkDM7MLDNuokwc4MzMchNpVe/hzgHOzCxLjnB+B2dmZllyBmdmliMncA5wZmZ5coRzgDMzy1FsG+waDDoHODOz7LgbJbiTiZmZZcoZnJlZltxE6QBnZpYjv4NzgDMzy06EAxwOcGZmmXKAcycTMzPLkjM4M7PsuIkSHODMzPLkcXAOcGZmeXIG53dwZmaWJWdwZmbZ8Ts4cIAzM8uTA5wDnJlZnhzgHODMzHITXk0A3MnEzMwyVSvASdpT0v+UdK+kFyWFpL0blNtZ0hWSfifppaL8nzUoN0LSPElrJW2W9ICkY1t/HDMzA9I7uFa2DNTN4PYFTgA2Av/aQ7lvA2cDXwBmAb8D7pD03kq5LwIXA18HjgRWANdJOqp2zc3MrHv9HOAkHSfpBkmPFQnNakkLJO3am2pK+myRNP1bg3MtJUN1A9yPIuJNEXEUcF03lXwP8BHgUxHx9xHxQ1JQXAd0lMpNBC4ALo+IL0fE8oj4GLAcuLxuxc3MrDvRhq2pC4CtwEXATOBq4BzgLkl1WwcnA58DnuymSEvJUK1OJhG18tU5wCvAD0rXbZH0feCzksZExMvADGA0sKhy/SLgHyTtExFr6tTLzMy60f/NjLMj4qnSz/dI2gBcCxwOLKtxj6uBxcC7qMSjajJUHF4uaV9SMnRrs5u3s5PJAcCaiHixcnwVKaDtWyr3MvBIg3IA+7exTmZm1g8qwa3LfcV+UrPrJX0EOBiY102RnpKhAyXt0+wz2jlMYALpHV3VhtL5rv0zEdv1Ya2Wew1Jc4G5AGPGjmutpmZmWRu0mUymFvuHeyokaXfgK8CFEbFBUqNidZKhHlv7hsw4uIhYCCwE2G38RA/wMDPrTjDg4+AkTSL1t1gaESubFL8C+DXwnR7K9CkZKmtngNsI7NXgeFclNpTKjZekSsWr5br1/LNPrV+25KpNwPq+VjYTezC8v4Ph/vzg7wD8HXQ9/6t///7+2Tv4j5v3aPG+O0sqB6qFRaKxHUnjgJuALcCZPd1U0geA04GDGwSvtmpngFsFfFjSLpX3cPsDv+fVNHMVMAZ4B69NPbvevf2y2QdFxBslrYyIKa1Xe+ga7t/BcH9+8HcA/g4aPX9EzBzAzx8LLAEmA1MjorPJJdeQhpR1ShpfHNsJGFn8/FLRIbHlZKidnUyWAKOA47sOSNoJOBG4s6gwwO2k3panVK4/FXjIPSjNzIYGSaOA64EpwFER8WCNy/YDPk4KYF3bocAhxZ/PKcqVk6Gy2slQ7QxO0nHFH/9rsT9S0lPAUxFxT0T8TNIPgK8WD72mqOg+lIJZRDwp6UpgnqTngftJQXA6aaiBmZnt4IqxbotJf3fPiogVNS+d1uDYV4GRwHm82rJXToYuKZWtnQz1pomyOsD7G8X+HtKYB0htr5cBlwLjgQeAmRFxf+Xa+cALwPnAm4HVwAkRcUsv6tOwLXiYGe7fwXB/fvB3AP4OBuv5ryK12F0GbJJ0SOlcZ0R0StoLeBToiIgOgIi4u3ojSc8AO5XPtSMZUj+/4zMzswxJWkvjjoUAl0TExcWcxWu6fu7hXneTAtxhleMjSePkzubVZKgjIq6vVUcHODMzy9GQWi5H0tskXS/pWUnPSfonSW8f7Hr1h7oTmUraXdK3JK2XtEnSUkkHDla9+5Ok24tJWS+tHM/6O5B0lKQfSXqh+P9+paTppfO5P/+hku6U9KSk5yXdL+msSplaK5ns6OSVW9pqyAQ4SbuQ5jZ7N3AGcBrwTtLcZK8bzLr1k6YTmSoN/19SnD8POJbUk3W5pD0Ho9L9RdLJwHsaHM/6O5D0MdL4op8CHya987gO2KU4n/vz/xGwlPRMZwN/TpoO6tuSzikVrbuSyY7OK7e0U0QMiY3UIWUrsG/p2D6kgYWfHuz69cPzvrHBsdNJcxRML34+uvh5WqnM60njQ7422M/Qxu9id+AJ4OTieS8tncv2OwD2Bl4C/rqHMtk+f/Esf0MaRzuucvxe4N7iz+8pvoMzS+d3Ir2vuXmwn6GXzzui9OePFs+1d6VMrecFJpKmurqkcv0PgV8M9rMOxDZkMjhSr5kVEfGHweGRuon+mPRLnpWoN5HpHOC3EbG8dN2zpH/R5/Sd/C2pW/D3GpzL+Ts4C9gGfLOHMjk/P6TJdl8hBfqyZ3m1BarhSibA94EZksYMQD3bIlpYuYXtn7flyYqHuqEU4A4AHmpwfBXDZwWC6kSmPX0nby+mzxnSJB1GylzP7aZIzt/BYcCvgJMkPSppi6RHJJW/i5yfH16dq/Brkt4qabyks4EPkibrhformeTCK7fUNJQCXE+rFew+wHUZcGo8kWmzFRyG9PciaTRpWp8vR8Tqborl/B28lfSe+QrS+lcfAu4Cvi7p/KJMzs9PRDxEGmd7NPAb0rNeBXw8Ir5fFKu7kkku+nXllpwMmdUEhjP1YiLTzFwIjCUNJB2ORgC7An8REf9UHFtW9KqbJ+lrg1WxgSLpncANpKzj46SmyqOBb0raHBGLB7N+tmMbSgFuI43/Ndrdv2ayoJ4nMu3pO+k6PyQVwz/mk160j6m8RxmjNCnr82T8HQBPkzK4uyrH7yT1mnwLeT8/pE4mr5CmgnqlOPZDSW8A/k7S96i/kkkuBmzllqFuKDVRriK1KVftT41JN4ciNZ/ItKfvZF1EvNDPVexPk4GdSS/Ey5OyQhpCsRE4kLy/g1VNzm8j7+eH9N/4gVJw6/IT4A2knoKrgH2KoURl1ZVMclH3eVuerHioG0oB7mbgEEmTuw4UTTWHFueyotdOZHpMNJ7I9GZgkqSppet2A2Yz9L+Tn5MmZa1ukILeNNIvcs7fwY3Ffkbl+EzSXH9PkPfzQxoe8t7ifWzZnwCbSVlI3ZVMcuGVW+oa7HEKdTfgdaS/0B4ktcHPIU3m/B9UxsjksJEGdgdp4upDKtueRZkRwP8FHgdOIv1FeDfpl/5tg/0M/fS9VMfBZfsdACJNbvA06f3Th4C/L76Dv8j9+YvnO6543juK3/sPkQYtB3Blqdz3SVn9R0k9LK8nBcCDB/sZ+vjMx5X+Djin+Hlqb5+X1DlpM/BpUmedq0mZ/6zBfs4B+S4HuwK9/A//dtIL5+dI71/+mcogyFw2YG3xP3ej7eJSuQnAPxR/ob1IGsT5nsGufz9+L68JcLl/B8BupF6D/0lqfvoF8JHh8vzF8x1ZBO2nit/7nwOfAEaWyowFriRlfJuBfwcOH+y69/F5u/u9v7u3z0taguZzwGOkIQO/AI4b7GccqM2TLZuZWZaG0js4MzOz2hzgzMwsSw5wZmaWJQc4MzPLkgOcmZllyQHOzMyy5ABnZmZZcoAzM7MsOcCZmVmW/j/Gic/fGteXGQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAFwCAYAAADKcW5qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOy9f5QsV3Xf+93dPTP3SlcyEpgfFrYw4BUjgUAgHGI7Tz8ISNiSgAAGgoEHsXEcTJzkgWNQ7AcyBGJ7OU5WsI2ek2flCScswBiE+SlL8PK80AOBA5JMxEOWENeAQegK/bp3Zrp7vz+qZ+q7d9U+XdPTd6667/6sNWtOVZ1Tdbp+9On67h9HVBVJkiRJsmz0jnUHkiRJkuRokANckiRJspTkAJckSZIsJTnAJUmSJEtJDnBJkiTJUpIDXJIkSbKU5ACXJEmS7BgReaGIvF9EviYih0XkFhF5u4ic1KHtD4nIlSJyx6TtV0TkrSJyYkvdnxeR/yki65Nj/JOufRzs9EMlSZIkCYDXA7gDwJsAHARwNoA3AzhfRH5cVcdtjSaD2DUAVgD82mQfTwfwFgA/AuDFVPfnAbwLwNsnbZ4J4PdERFT196d1UI5VoLeI/CCAfwfgWQAEVef/uarecUw6lCRJknRGRL5fVb/j1r0CwJUAnqmq1wbtng3g4wAuVNVP0Pp3oBo0T1bVB0RkAOAbAD6qqq+kev8ZwKUAHqWqm6U+HhOJUkROAHAtgB8F8EoAL0c1cl/X9oqaJEmSPLjwg9uEz03+n1Zoujr5f49bfzeqMUkmy38PwPcDuMrV+78APBTAT07r47Gywf08gMcCeJ6q/qmqfhDViHw6gF84Rn1KkiRJdse5k/9fLtS5BsD/B+DfisgZInJARC4A8MsA/kBV75/UO3Py/ybX/ubJ/zOmdeZYDXCXArheVb+6tUJVbwPwFwCee4z6lCRJksyIiJwG4HIA16jqDVE9VT2C6u2rh2qwuhfAnwP4MIBfoqqnTv4fcru4y20POVZOJmcC+GDL+psBvGha47W1/XrCiZWjjohsr6ciAMCaF7Wl1MIcTZJ7Zd2UjhvufoB6ZE5WuAdTrddhvd8W7VnEnp24BzuthPDEl66HavvOfZvIZM2rV3q20gn7OnS8483S9Z7a2Kxr3r/Z9cTVdLs74kY7P2J8/EXl+06Y375K9+F375v4cmzeCx0eEQC46IKH6J3fHe7qmJ//0v03AzhCq65Q1Sva6orIAVTf6UMAryrtV0T2AXgPgIejMk/dAeDHAPz6pP0v7qrjxLEa4E5Fc1QGqpH5lLYGIvIaAK8BgP0nHMD5F1aONmsr9aMwGNjHYnNY3wWjUfsA57+wunyBdWWv/HeiAaXnRp73f2mDKvap4iqV7Tlco33spy9uXr/Wtx90hbYNqM0gWO/7Gg2k/gdMBJ93Lo/c9RjztnG98yH5fvk2XI+38b4ecaJ1Hnvaj9SPWXTveWevnd6Hvv7Xv11/uX3m4Eprm37hfPbpB0i/8GOG4VvHlAtPT3RNe4W+lbbtBeOOz/VPP7U/vRK6XWt1Bz2yUd9jV/6/kzHoq3+yve7O727icx97YqfjR/R+4LNHVPWcafVEZD+Aq1GZnc5V1YNTmvxjAOcBeLyq3jpZ93+LyPcAXCEif6CqX0Q9RpwC4JvUfuvN7S5MYWHCBCa/HK4AgFMf+ggdTJ66Pj19ff/00YA3oHr8ZdJ5DAq+NAvVTL2xu0HjL7r2OsV6QX9Gzkn3uWfWAxl/sfzpl75XL4i9JY7QQHikTwMhrR+4b5yTaMDbT+ed32xW3bMfDYq8ax4g/RejBF+uPOL7rxu+JErHjAYuwN470WB377rt3LU3jah9vY33/XcfZw+0j35B3PG39WD1le/Un8LcX413nvo6ft/K7n5plQYopusPkN3SdYDpwiyDJbf56afa56XLOWjUCRrxcfqNjtb3xEuftgYA+Pg3997iJCIrAN4H4BwAz1LVGzs0exKAQzS4bfHZyf8nAPgialvbmbAD3Jbt7a+mHehY2eAOof1NLXqzS5IkSXaC6u7+piAiPQDvBnABKofB6zv27FsAThGRx7v1f3fy/28m/z8D4E4AL3P1fhbV29tfTDvQsXqDuxm1hwxzBjqMykmSJMkUjr6N5J2ofCbeBuB+EXkGbTuoqgdF5HQAtwK4XFUvn2z7IwD/EsBHRORtqGxw56AK+v48JgOXqm6KyK+hCuz+G1TelxcAeDWA16kq2VzaOVYD3IcA/LaIPFZV/xoAROQxAH4CwK922UHbW71ftxIaGXauSyg7qRTsdiyfsCypTmKw9Xi9tq6v+hD0raPEGfHTZxzYLvsz8+EvUahLj2w5vfrWGbIND8Ch/lpd7tWS2j46Byc7ux3b8dboug2MrEl2Iac9RPJlyYYXyReDgvxrZE2+VlyncN1MG5Irv/Q132ZM9eqenjBoP2bpuyyS9I7m198sMuI4eI6OJnycrnY/liVLUjkvBKurfffa6/H6gbvhVwb11e/3WhKGdHsJ2y3Pmfy/bPLHvAVVVhNBZR3Y/gCqevtkMHwzgLcCeBiAr6MyQb2NM6Co6h9I5ZH2vwF4A6rB8JdU9fe6dPBYDXD/Byp30A+KyL9G9az9BqoP+a5j1KckSZKkI6r6mA51bkfLG4Wq/hWAn+l4nHdhxnHhmNjgJoF8FwD4Cqqo9HcDuA3ABap637HoU5IkyTKh2tvV3zJwzLwoJzknX3Csjp8kSbKsKGRpBqndsDBhArMgHFtlbDFxcDhjNWxy73YGAqPjs42Fdt6wKdAK0wVyPddCMHRkd+tqw+uqz1945sO2yx//ynq9YZ1CUHouzsrY6mr73BGyzXEZAPaR3e2ACTOo67Ar/qp7dpXasD2uXwgt6ElwEjreE0oVzXrft6h9YMNr1qvL1k7VHvrij8mMg/026gXro4B4v++SPauLra3ct93FI0RhDz7mUYzdjcIzTMyQ7Uv0fWLCWNzJ6QXfR2x367kYl7bwKH/scQ5wyz3AJUmSHJ/kGxyQE54mSZIkS8rCvsF1yhhAZZYFbKaLbnoUSybiwg+EdBbObtEjF3cdmSZGlhiTVFZMTbVLt99I6uoqVz7zcbWsKPKo7fI1N7op/DgNWJ/kyjFJnCMrUUby5Yl0rk+iTu8f2E6vsXzJciXVadzs5p5oz57SwGxjl/+4SSdZ04eeBO1Hpk0hdCVob2XNuI0EfRs7WbckWdbHjJfL2VjifUxb7+FrGh3HX/dLn1r//o+ekdLXR/Q102gSyZemGEuhrX1Q5BscFniAS5IkSdpRdPvxsezkAJckSbJ0pA0OWNABTqSWE6Tra/wcj991BoJZKPWZ1QsNJAo+B35KGtagpH11UbaKNlxw5g/ZftJzZeRLodnl+y7LzpiWx3W9+3V/XR7X0ucpzs3vRFJF91GS7TXSsMYue4r1tmRpuV7vPS0j+bLoNUhlCaZtUqczdclIw+vH7vhje9C6nyyr+jtsjvdxyVOSP9swSD4deXECTW/HTlCbKMHRP3yquwbj4IFhmbkk2c7xS0fhjzO9zTjf4NLJJEmSJFlOFvINLkmSJCmTEuUCD3Bb3kYmB6lXXAKZxHg3FQJkw6DcguRi6pHOUlQUWNXgYFmve3HyZuMi2t6+qZEEn2HcWvRNOkllgP3cz3xiLV/yx/nkTW5OxDG5mbLLqdKsxD1K4kzSJQAcJo3uZONtWdfZ5z7cSp/L5HlpApZdUC6f6z2bs52OyTL1HL1qgYLnJa0fOS2U2/BEsZEMCcRSJEuPI9cmkiVH7asbc/8xL35ae5IFn8CBE0KYoaJwDbS9mv2O8Amag3q8s0bShoIH7da6dDJZ4AEuSZIkiZAc4JA2uCRJkmRJyTe4JEmSJSRtcAs6wIkAK4Otcpw42ZiZgpQHfm2kgXed8LTLeo/pNruoe4OYscmRHaFD8tqqeeCSTTtoZJ0I7IilzxZt430/64mPtn2jz/3xL/1NvWBscFR2qWGOjCkTyrjOnnIy2YxOcGECnA1llXbHtjmecNX3sx/ce14Y6pJ1p3g+g3qlCU/DTCiczNsdZ0Q7HAXu+94Wxm0il/9N34bqbRhbXXt5HrzyqTSBLPXZhNW48YBDKjiMg+326q92YFBr1GOC75kwSgFukuRx43CT5ZQoF3KAS5IkSWIyk0lFDnBJkiTLhkoGemNBBziBzTwRESZ1LSWp7eAWPw9sstRoXjHntkxllhutFEtl195sM8eR1jp+H5Fc6dsYN2q0UwotePaTTtsu82X+2E1/SztwDuImE8q+7eI9JF3eM7LO4yeP6p0fILmSs5/4eeds9pO6HvezVzB9mKiWbnm+DZEsWQpdYVnSuuL7Nu31hoF06bexFLlBx1x3xzHyY9eHrIPe/7Kz4zospdq5IekauhCIcRB+U4rEsaaMdt3az93H17RXmHeSGXGGnr2PVlkYFnKAS5IkScqkRJkDXJIkyVKiGQW2oAOctGT5QMsrvZFp2l/pSzJAeHinLXVJ6tzFk665gzgB7ziQPMy8Vz6DhNnWLj36fvbYMy6QHptzhLXXK81/FrVn+eaiJz5iu/yxL3/PNhodph2QfMnSZc/OQXcPSZkP0Ic7iZTM/S4z7wrJkit0stnb0ifz5XPa77V/8K63R+hRWcoWEsiSPisJb9ukbZEMCcRSZFGGNBefde+g7HjxU+qvLf4EI7rs3iPS3KKcWJu71UhO3i4X8nyQjVxBUSqTqE61pqVkd9CQNc33WfOeSieTisUc4JIkSZICmckEyEwmSZIkyZKysG9wW7JPOUC23Vuya1CtScpcCCg33o3seVX4AeXnrmvvSyz9Ga+uIMDWf9BIvoykSyCWL01fvJRq2rfLL/7HZXR9Irny2T/6faa98ba88Zu0gSVKmo8OMPPODcnb8hCV73XOmixfrtG1XgnKgPO2JM/NvkncvHPMeS55UXbwjgScLBl4RD7g2oRSpJGJKUAfcPKjtq93EuULn3LCdnlEHyJ6LsX109z7fLLMBtsmSuSuhWcMgUd0KTjcPFdm19paBpwXdAZ6hyzsAJckSZIEaA5wQA5wSZIkS0kOcGmDS5IkSZaUhX2D29KbvTbNsD7fp0S7rM+XfuXYbCPtZcDZoILdzZK1QhozI7ZXNFo/dyX2enYTq1Il5zqOXrsRzGRfcD+TzOcxWVIKaWI6uFSXks+ybenZT3xU2+GtbQ6wyZtNJhSyzfVtaMEhSuTMWU3YNrfiOrdqsp/Q4ekkNkMLdha+4u/jLplIGjY4asMTyB4JTGbVcmBrM5PW+tk6h63bnv/kk8PDcF+tbZztcfV6P1GtSWJOq437vws9Mjb8wGXfn/cuYQKz4L8jxlMymSiOfqouEXkhgJcCOAfAwwHcAeBPAPwbVb230O7NAP73YPO6qu6jurcDOL2l3vNV9U+n9XFhB7gkSZIkZg8kytejGtTeBOAggLMBvBnA+SLy46phMOMfAviYW3fiZN2HWup/fLJf5pYuHcwBLkmSZAnZgwHuElX9Di1/WkTuAnAlgPMAXNveLz2IakDcRkRejmo8urKlyZ2qev0sHVzIAU61zkxi5Qpbz8o+nE7iaPUslhuLyRyMctce2uCXo0TSVj6J28cZRrpJY1HIQKNe0MbTde68LUoJaPgwvF+WLgHgEzfyvHOjoOxc3Hur20WWLw+N65tq4CTbE0j6WiXJd2AkNdvGS5bTaEhYnLQ3yDByeGQPEmYiMTebOx8k59pwgHGwHuC779KzTtkuj+LkJYYoc5DNVuJPCIUTRO7/vpdRBAO1b7Qx5aOkV/pjtj0Le+BF6Qa3LT43+X9ay7YSrwTwt6je1uZGOpkkSZIk8+Lcyf8vd20gIj8I4HwA71b1v6AAAJeIyAMisi4i14vI87rueyHf4JIkSZKYOeWifJiI3EDLV6jqFVFlETkNwOUArlHVG6J6LfwsqpetNnnyalRvhbcBeASAXwLwARF5uapeNW3HCznAffcB4L/cUF28V55Tv5/7OeL6pPP0CplIIqKExKWsApF3kze3Rvu20mOcyWSWhMZGvkQ7XWXNxhxw0f46HqdL+64JmkOvUtfmWU+sVZRP3ngH7aCQ9Nd4XlK5V3tXDvurYO4h+ZI7tI/0k9WCRLlTuRKwXpSciYQTsxzx17BLJpKxywajwfkIZEigW9YaT/TMchv25m3cx+DnkryrC8c3WXz6/OzsbuBoeGHzd1PHNl3Yqezfwp2qek6XiiJyAMAHAQwBvGqHx3kFgL9U1S/5Dar6OnecDwC4HsDbAUwd4FKiTJIkWToEusu/zkcS2Y/qTeuxAC6cOJF0bftjAH4U7W9vDVR1BOC9AB4tIo+aVn8h3+CSJEmSMnuRyUREVgC8D1Us3LNU9cYd7uKVADYB/PEMh5/6jrqYA5zWXl48L9yKc19boYjbAUXY+vncGJYOR+TWxQlehy4B73jE8ke9ntt4DzEzPx1LlEau3J1E6dmtZNHl+NW2dk/Ocpv2cpe+NLYF9XzsdHgbsAzXCOWJEgUHQc4AICRRknx5ROv1R3wfqHODGb6nhtH5Kd0skdxY8ioN6l1y1kPrKj7OO+hakeD+iGTJxv0VeP1az2TbSoLEyaaNc+ftmePs/MJF300PxqRbItID8G4AFwC4eKeu/CKyCuAlAD4aeGS2tRkAeDGAO1T1W9PqL+YAlyRJkhTZgze4dwJ4EYC3AbhfRJ5B2w6q6kEROR3ArQAuV9XLXfuLAZyKQJ4UkZcCeC6AjwD4Oionk9cCeCqqDCpTyQEuSZJkydA9iIMD8JzJ/8smf8xbUGUfEVSRx23+Hq8EcBeADwf7vw1VCrDfQjUQ3g/gBgAXqWqneLkc4JIkSZaQUjKEeaCqj+lQ53YECquqPndK2+tRyZ8zs/AD3Ef/Ry3wn7xmbR8nrFDWCDKDPPvHDmyXvX2A7WZmPVX0Wr2xz5HxY5NtcM4sMwzsc1bft22s+STOWLIXdLWBRW2abtzT912ysXSh0YZWPPOJP7Rd/vObOGTAXTheNAY+usEaoQX0mPH+OKWGuPQ6tG1oUm90dXwO4jjG7YmOG30zoRIF+yLVY7tbVxttZwJX+vBeKUyoG2Uv8d/C4f1aDMXp0H4OlJK/JzULP8AlSZIkTXI+uBzgkiRJlhDJAQ4LOsD1BNg/UWpWSdlZG1gdYIU+3YCSr173+XqqoqGb/2yT1Jgjw3rbeU+qM1V4GXOTZMkNmlSLJcpN513NsuQoCDPwGvpuZY5ZpIyoSdeudM5YEklNHY/TZV/+82u4UNo5Z/hgiZLlvoL0JxSCwLKkeOmx176tUa8D0awlHedps5qcbdNVlnwwoV30StjsJzxYRPPEAfGchUdTrmxDkW9wwIIOcEmSJEkBPTa2+QcbmaorSZIkWUoW8g2uJ8BJEzlylZKglpLSGk9FEg98VpINWmaJ8n2frTd8a8Me6CVPJIly2C5LbrrUEkaiNH2r6ZrQuCuRRLlXXliL8ovygjNrj8prb74jrhhJfz5tjdBN1WNvyYL0GEqR0UV09cPJlE2lwqb2hNOXPPn7bbWZHCJnyPBh2ncoN465u+NHXr6lhObciKXLxlR1fKmCy+4zpjBhIuoHZf6TvWUhB7gkSZKkTNrgcoBLkiRZSnKAW9ABrgdgbSJNsizZiNHm3LH0us6SwObI3gSHSZZ8gMp3U/mlZ1ldkz0nWfKMgrn98jiQPOau6HUIoJ7jbovM8ujN5AW6V894SRI0AdR87xQ6F0mU4fqj90EvPusR2+Wu98pue1P6OF2k9sacax33vWMKns7Go7KQwME4ctJ3Bvez5FH9nr/cqAoP1DtWLI5J4GiSTiZJkiTJUrKQb3BJkiRJiQz0Bub8BiciLxSR94vI10TksIjcIiJvF5GTXL1TROQPReROEblfRK4RkSfNsy9JkiTHLbo1o8Dsf8vAvN/gXg/gDgBvAnAQwNmopkw4X0R+XFXHUvnoXg3gMQBeB+AQgDcCuE5EntJ1uvOtiQWH7G6/aX+xGD2c1nP2kk0fJkDb7iH73AvPqvfgs5J0yZzR817cXG+XNrhZbsbOGUbmeExDycbS4ZjzNjnt2fPcyX0fzaTGu6Fz9pP2k8rJxUsu9tGWrtdqlnpRaEDDBsfJmvm5NI0KbWh9lAgFKIUTtGcrAmASeIcu/z60oMMNm29w8x/gLnEzs35aRO5CNaHdeQCuBXApgJ8AcIGqXgcAIvIZVHP//AqAfzbnPiVJkiTHIXOVKINpxz83+X/a5P+lAL6xNbhN2n0P1VtdcX6gJEmSpBuqsqu/ZWAvnEzOnfz/8uT/mQBuaql3M4BXiMgBVb1v2k63XtGPkIzoXfFZbuRNHE4wchfyBU9h1/4xles6/tKzzNE3PxkkWO9DAzipK1rLQCx/IFjvq4RzVZX2G0khhfv/waTfP5j6cszoKosSzz6rzubCbu0i9oRGGUJmma+sa2hAl6wmPbczNhGYLR2zn3TtJ2NlydjnXzucq2aYQPnGVoQzAh5XHNUBTkROA3A5gGtU9YbJ6lMB3N5S/a7J/1MATB3gkiRJkphleQvbDUdtgBORAwA+iCrc+lVz2N9rALwGAPqrJ02pnSRJcnyTA9xRGuBEZD8qm9pjAZzrPCMPoXpL85xK2xuo6hUArgCA1QOP0M2J/HiYPJDudVlJOL/xK86uX9itDODkAvM56jJLjN6TjOWPKJNBwwuKNnJS1Ki93xaVI+mzOub0vhWP2aEvgD1vXRXCLo9id2moY71dtl82zqck0yYBeMEbmBe1IPdFzJYsub0eP5e+n6ZeIFf2evaos2Q/6fK8NBKZGC/V6fsCgJFxsj1Ob9gOzH2AE5EVAO8DcA6AZ6nqja7KzQCe3dL0DAB3dLG/JUmSJAWWKJZtN8w70LsH4N0ALgDwPFW9vqXahwCcJiLnUruTAVwy2ZYkSZLsAsXuPCiXRd6c9xvcOwG8CMDbANwvIs+gbQcnUuWHAHwGwFUi8gbUgd4C4De7HGSstTR5iJIg/yOXBJl/wHASZJYivHejkDDBU3eVLnc071IU6FktT5cl/RxQUb0RuYWa9WPbLzXb2o/jJUqWNbtKlF28PWehtK8ucmNJfmUvt+v/ZmPnnVsQnv6407fL/t7nuQwHtM3XM7Dcx97A9Ew0g66peYe53YBYioy9KF0/A1mTvws6y6IzfPfb5989l8FCKTi8W6B35+4tLfNOtvycyf/LUA1i/PdzAKCqYwAXA/gkgN8D8AFUsfznq+rX59yfJEmS45J8g5vzG5yqPqZjvbsAvHrylyRJkiRzZyFnE1DUQdwve3KtoQ0L6ft6gXeShz2sBjTZHEsZXv6YRbIw0l8gRXiJMpIVWYocFeegq7f1zL4oIN63ocDerh6iXaTMknpSkj/bjl/qD0uPJfmVg/+f/PDV7fIXby90dEH40Uc/Zrt8eLNev9K3J2RAN/aY5PmVjl56pbyOTBfpzwdq7zT/ZNdclOX8lfG2LkTSY8Nzu8M93jRxTL8mqVAu6ACXJEmSlFkWmXE35ACXJEmyZCzTlDe7IWf0TpIkSXZM1/k/W9q9WUQ0+Dvi6vZE5I0icruIHBGRL4rIC7r2cSHf4E45Abj0qdXPk/WN+mfKcBRr230yEAyMXSnOXsB2N7ZP9Pu2DS9720GEtbu1283G7vMM2Z1/3P65R/STxdsk++P2bWxnc8lgzDmwmVC0dX21rS4be1zBDhq5PUf2vKZ9sd3WNgyytwDAkD4s19scL560c/KpjzPL+3r1h72PnpFVsq3t89egz88Fb2FXetuIb3c+v7P8ci7Zw6Jwgi7u/6V9d03wbPe185QtBRNctznk3P3e5e1sDyTKqfN/Bu3+EMDH3LoTJ+t8LPRvTI5zGYDPA3gJgPeKyMWq+pFpHVzIAS5JkiQpswcSZZf5P1v6pQdRDYjbiMjLUY1HV9K6h6Ma3N6hqr89WX2diDwewDsATB3gUqJMkiRZQo52HFzH+T+78koAfwvg47TuQgCrAK5yda8C8CQR+eFpO13INzhBnV2B3X/vP2zfiNnNfzDghKZxJpPoR08vCB8AgD7tpN/rKF8YV3aWKEmWcG7cPZYiSaPjQ1Jil0Y4g5EsSZbskSTXc7KolwLr5oVMFUFmlCgrimccSJEmHML1a2jq1eUNI0Pajq7T+WBZcp3aDx5ipT8+Du65tdn5o8G+h9ZloUdWar3xnk17Ru+hi38ilU8y2UYsRsId8P1Rrx66a83Phbm+Pd5Z/CCY0IBSKE6QINkkWC6GCcTb6v12lPR2qfz5895Flixl4el6nD3i3Mn/LxdrESLygwDOB/C7qjqkTWcCWAfwVdfk5sn/MwDcVtr3Qg5wSZIkyVHnYSJyAy1fMZnVpZVg/s8u/CwqNfFKt/5UAHdrM+jvLtpeJAe4JEmSJWMr2fIuuVNVz+lScZfzf74CwF+q6pd22G4qiznASS0F9kkK2b9mLyhLciZrhckCYndtZJZxN5mFZRqWUnpB9gTASgyskvZIUmPvxqoHNKedaRU5K/kfPu39GVI9L9NYaaf9vElDC6kb6ahdKGnOj1eXR8bbs12W3HQfeZ20M5YlN+geOOIkyvuo3hGWb8eb7WUAGJEX88qBuszKyti5ryoth45lDr6R1u+uy/2VutyrM65gvG7b99a2i/cP9m+XRx31NZb7eixnuxvZzGsYzItYzBzE5ZLEyH0LpMxy++mmg1mylXSlmCicZeMgW1Aj8frUA+5dHNyU+T+ntf0xAD8K4J+3bD4E4CEiIu4tbuvN7a6WNoZ0MkmSJFlCtoK9Z/3rgpv/86da5v+cxisBbAL445ZtNwNYA/A4t/6Myf+/mrbzHOCSJEmSHdNx/s9S+1VUcW0fDTwyP4Zq8HuZW/+zAG5S1aKDCbCgEqWgDq4eDGpdYW3Fawzt3nxGBvABwxwATV6MLGv23c8b3kcUhOox9dgjkSRX/+uDEwf3qA/9TnIlYIWN9s75QO+4TSzzlgJ22/YKuOBs6gRLkRskQ647FfAIbTtM7e+m9T4RAEaHgzLJfV6i9MtHi0jKZN1daN66/qqt129vfwS1XPmAuzicAIHng1sJZDO/3FWWjCjdN+ydbEwCRuukUuIAACAASURBVJYsePZGx+waqN3x83SqVvCINJ6TLP8Wd9z+IfYg0Hvq/J8icjqAWwFcrqqXu/YXo5IbvXMJAEBVvy0ivwPgjSJyL4AvAHgxqgH10i4dXMgBLkmSJCmzBzY4nv/zMrftLaiymgiAPtrVwleisqN9uHCMywDcB+CXATwSwC0AfkZVS222yQEuSZJkyVDYWdWPyjE6zP+pqrcjeMVU1ed2aD8C8NbJ347JAS5JkmQJydkEFnSAE5HtbCKrZIMbrdofCnx9G/aXrTreLMMhBJwFY8iu0j7hLNs72Fea6rgXdNb+Z3GpZndt5eTRvN4d0yZGmW5bK/XH2A0KWSOiI3rbJ2cZiexubGd7wKXU+N4mZbRh+9wwsLP55RHZs7q68h9ruJ/DI3E9vvnG9SP/wHjFVNtH12A1CNvw7uqliW/D7nA5sJv5ZMm9KPzGrG/fr9/3LMTPQZw5OXrcxgUbfnQ+o4xCAOqTcDTjHBaUhRzgkiRJkgJ7GAf3YCYHuCRJkiUkZ/Re0AFOBBgMmk45zYSk7QmJR5FbP2JZ0yR7lVgvYAlnQEmYG/NTBVkXZnG1Nvtid2p/PsKoxzh8wPSH6vU4tCGazM21Zylm6CXKUXs5kiXv2rT9PEISMob3B2UnUaqLNVgmOJvKeNhaHjqJkvM1j4ykRmV3GA0kudKta6XIYL1P6sxJlakNz8VYDE+R9npM8XGLso009jbdtb80F+I4yLLkpWHzfeLtH5NeFB7L44aFHOCSJEmSMkfbi3IRyEwmSZIkyVKysG9wW5kNeiv9sI5N0Fq/72+yYlPQAc28ZKTZbLh6LIX2yRONZQUvUXJmBuPxVfrR1UFysJ+51Fxa13uvMJY5WZYck+dmyXkrkq3Gzj7AkiUnS+aMJfexR6VXF0NZ8gHqwIJ4R84D81m7iYcjvr48L1lBNteCB24XwiTK7t6N5mOUDomXq+Ps/G2GJfnxuP0Z9fJ8ZGKIzmdpG39/+Dbm9GxLlLF54XhlYQe4JEmSJCYHuBzgkiRJlg+dy3xwC89CDnD3PDDGNV+oPOIufPoJ9QZxcqUEUgbJlcOhaRF6HtnA5FiWGAQSA3t7AcDIeIW1B7s2ZJYOgZzGI9NJFkYlJe3RyCo+wNZ4nJIsObUnzX3bRLK2nkm2bIK+yYuS2zSSIG+0bzueZEmmxbNuVkqevV28fpv3cXvZzEHnJf1++/MSJ2GO79DQa9nV41uHTycng5jtOYgDvdkUwl6U5VzLsZnmeGchB7gkSZIkRpESJZADXJIkyVKSA1wOcEmSJEtJxsEt8AC39euE9fhVbz9iO1eg9W+6dB+cVHkUTDjo7UcaCPkjzrjg2vR5YlPj6hy7Spfsc200qnB7zkRi9uuCCTiZbZABZpYcrw17R2T77PorVIKUGMcL3g7TW2kvcz13mvrSfrLnfTolWChmMgnsbtY2x+195iAKiwlutkYSZOoc28063/tBNphGJpMg2XLp3u8ysXK+wWWgd5IkSbKkLOwbXJIkSRKQswkAWNABTlWwOcl28YnP1lkrnvOMA6ZebyVI0GrkDqsdbtLyBoUQsPtuU6Ksy5ygmdVP760+Ngljqc/C7WM3fyNfzlFD6rynQvKVSMIpJ9OtyyyVrdLnPIEq3dOziYLR3x/0k27x0brd5kMNHqwYnTiQG/trtk1vtX0btT/R6TdrnC0kmGetq2xeuo/MHHDc5SB0ptEHkiWj8IHmPdmu6bFcKY24sfqhZXuWFEaOKHuJFqTHrlljGO7pvsk1WTchQjnAAQs6wCVJkiRlMtA7B7gkSZKlJF/gFnSAG2s9N9igX1/GT3z2PlOPJct+r31utsY8bSZjKyVoJrli5LzNRpT4l7eY+Zzc7dYLJAsj1XnHOJaNgjZd1UrtePsbYSeQoHzWCPbENJ/HyMT2+CskNUUJsHkevn3uut0zqmW4I2OS5KJ50QA7HxxvU1fPtOmQGaWURcRs6wXr3bLxfOT19Pi6a8AS4xpt2k/Pyz53DVb7XK63rVDZJeQxknp4rxTuySjRuG8Se0R39EA2ycVZOqzXj5x2aO/r9nuymXy6fZuRIb25opDhJ4Jvg/2T67gggvuespADXJIkSVImbXA5wCVJkiwlOcAt6AA3UuDujUo+4OBullUA4JOfu3e7/JxnnLRdNp5XLkA28vAybFodgT2xStPMmzZBECnLkr45y3rsCToOPN6K87QVksxGGC9I44lq98BS5MpgeoCt30c/uKb7VurycGQ/3CmcsJZlJyoPx9bzUrVeHh/FrA+94AyXro/3Mm1fz/KghvVYMmaJ0cuN9hrU63nKxUFDWub2dPwggXi1DeG2aesb20IvTreDIMlBV6meKSZojoK7OXGyew5mGYdYjt3yfvVJpNPJJAO9kyRJkiVlId/gkiRJkjKpUOYAlyRJspSkDW5BB7iRAoc2K335ANl49o1i+8A1ZI9jrXrk3HL/l6fUtjo1mQioPLbatkmQGtjd/Npo0kS2zfV88lttt2eZjClsa/AhEFF/Ck+Cdfdu75t3yWa7m5kkldoPRqYJVul8sK3NnlsqeztGYBfZ7UPeNeziaFo7unwef9t1+dwN21iHrDN9d615eRDY43ybKE9w5xAXc0La16tPGh7Y3aJ9AfZZ1nF0T/r7sH3fhcN0wj9jfE1WJl8A/vTlAJc2uCRJkqWjStUlu/qbhoi8UETeLyJfE5HDInKLiLxdRE6a2rhq/wQRea+I3Entf9nVuV1EtOXveV2OsZBvcEmSJMkx5/UA7gDwJgAHAZwN4M0AzheRH1eNsyKIyDkArgXwKQA/B+B7AH4EwIGW6h+f7Je5pUsHF3OA0zEwfAAAcN+wTrJ74tD+6thHn85IegXJ5bov1FLmBu3v75+1b7vsMx7YBKkkUZBo4NsgkC/MHHQNzUZbSq7M0pJrPYscZDO98JHirBGc2WUQZChpJJ828mt7X0xi3MIH6Op6Hp6PQkYNE0ayS11Sg/sGcNc0mi8skM38snYoty234a81E8mSfRePYDLycMhAoS9R9pFxpAOO42eHsYmO4+fa3ruF8x7I6KXzzoT3rlvum9CN9jp7IFFeoqrfoeVPi8hdAK4EcB6qAayBiPQA/BcAf66qz6dN1wXHuVNVr5+lg4s5wCVJkiRFjvb45ga3LT43+X9aoel5AJ4A4Bfm3SdP2uCSJEmWDd0K9p79b0bOnfz/cqHOT07+7xOR60VkU0S+LSL/QUTa5r26REQeEJH1Sf1O9jdgYd/gdDs57r2U0eIUl91ig7wqvRS5RVdPsmu+cGS7fN5ZNNcWbAaIESdyZimj+QnaMYqLl0yoQ7RD/gic1cTrJ0ZupE2lxLhGuTPZS3hXtpEG0nDJ1zCSHKPr4etH8/0Zz8+GJ1p7vdJxuiYRjoikqoZXqJHB6P5iCWzE6929Mm7fVsq0M46k0NK3XXB5+dQU53brKPmacxV5KrMkWNiZ9UptPzcAMKIVo+Bce6l9Fu9V6XQObRs+p4PJHJZ+v3OQKB8mIjfQ8hWqekVUWUROA3A5gGtU9YaoHoAfmPx/D4D/COBXAZwzafuDAFi2vBrVW+FtAB4B4JcAfEBEXq6qV037AAs6wCVJkiRHmTtV9ZwuFUXkAIAPAhgCeNWU6lvD9VWq+uuT8qdEpA/gHSLyBFX9MgCo6uvccT4A4HoAbwcwdYBLiTJJkmTJUOwuRGAneSwnsuLVAB4L4EJVPTilyXcn/z/p1n9i8v/s8HOpjgC8F8CjReRR0/p2VN/gRORjAC4E8DZV/de0/hQAvwXgeQD2A/gMgH+hqjfu9BhDeg13OZCNFyQ7cllFwHuvtXvzlWWr6QHQPXfDmCBSdCNS+4yUqYHGUaKQYDqS5DiZbUNyCdwTjTTVkK2kdRvLjT3jOWbbs1dZv2PAsU0O3O5l25hXbIbTGwVqc3nkkwfQMktnQ/JQZW/VoQuct/WU1o9by/44Rp4ryprTA5tL2FslPqPmsNEx+XyWLk4wN1vDizKQdkelNh0iuovPWLC+7+93Wl6JvCjbDz9XRGQFwPtQSYzP6vgdfvOU7R1nxZv+EY/aG5yIvBTAk1vWC6rR/iIArwPwAgArAK4TkUcfrf4kSZIcTxxtJ5OJu/+7AVwA4Hk7cOX/KIB1VC8/zEWT/6H9TkQGAF4M4A5V/da0Ax2VN7jJG9q/A/AvAPyx23wpgJ8AcIGqXjep/xlURsRfAfDPjkafkiRJkrnyTgAvAvA2APeLyDNo20FVPSgipwO4FcDlqno5AKjqd0Xk7QB+TUTuQRUvdw6AXwdwpap+Fdh+SXougI8A+DoqJ5PXAngqgJd26eDRkij/LYCbVPW/ikjbAPeNrcENAFT1eyJyNaoP02GAE6C30li76WSew8P21kPSEbxUxoyNLBHXs4GrJKnRr6C++0nEEl80bVzXwOSSXBi2534GXm1A4bMFZcDJioHcOLDT8Nk8hlReGXAdpfX2pPEy14vKVX/Iy5a2WYnSe6LWy7NJlHTdC/fXiDQ23sbz4G0O29f7bZv0HGwM65M7HNqLwPU2WdYkO8DI5XtlmTPy3Bw3PA3bJcYiJiC7Xs2nTcxx4h2b61EM2g6O2fGzMUUP4HC+wPZn1C9HX2F7EOj9nMn/yyZ/zFtQZR8RAH001cLLAdwL4J+iyojyTVRmq9+gOrcBePhk/akA7kf1dneRqn68SwfnPsCJyE8CeAVa5MkJZwK4qWX9zQBeISIHVPW+efcrSZLkuGF3sWzdDqH6mA51bkfLb0Gtfgn8zuQvans9KvlzZuY6wInIKoB3AfhtVY1yhZ0K4PaW9XdN/p8CoDHAichrALwGALDSlq4sSZIkAbaSLR/rXhx75v0G9yuovCLfNuf9YhJgeAUAyP7vz0uXJElSwCdgOB6Z2wAnIj+ESof9OQBrIrJGm9dE5CGoNNdDqN7SPKdO/h+afrAeMKgyuuwjZXfd6eE9tmPQkLhC6/tu3qgo4wm7BjcS49JipI2rn5uNwwkCV+eS+72xu7GdrOCKH9ndIpsZYBPlWhtar7VOta0usw2NbWur7s5jG9oqzQdn1ps69mKvDMatbXj9YODa9OvlPpd7ddnb4KwtpdvvLA0SSdtQAHvjjUb18iaVh0MqG3ucbb+x2Wsvk21uY1NcG5lab9PZtTfZpschCGS3a4QjBFlBTILoQkJzxtfbXj/LfIHeBhfYCqNE1i272KaUlSSKOSuFq9hMJtVRu96PxxPzDBN4LIB9qKLLD9EfUBkRDwF4Eipb25kt7c9A5fqZ9rckSZJdcoxyUT6omKdE+T8AnN+y/jpUg95/AvBVAB8C8CoROVdVPw0AInIygEvQDClIkiRJZmFJBqndMLcBTlXvRjV5nWHi6vo1Vf3UZPlDqDKXXCUib0D1ZvdGVG/xv9npYFLPgcQc9mECpIyskra0SrLTqpfkWDqkTRc/rZZifGLb6NdOPJdat/nPSi7/Nglxu/QoDfd9LrfLkgOn0RqJksorXPZyI0VwrBq5sF6/tmI/NEuOLFGumfVctqk71miZt61QedC3bVi+7PfqbZFcCQA9Wp5NomyXJVmSrLb1W7dFbv6bzuV/fbNetnJlvX59o+fa9Fq3rZNEub7hZE0OQTDHqddvunkah8P6HHI4gsme4qU/WmaJM8oc0kheHciKoVwJF86AjnR5ljXeFs1L2Mj8w/PoHbswgQc9e56LcjLL68Wo8pD9HoAPABgBOF9Vv77X/UmSJFlGdJd/y8BRn01AWyyoqnoXgFdP/pIkSZJk7izkdDk9AGuTYfN+VpD8zw56Rz/C6+md3kudr6DJIVjGY6+wrvNWMQ2Pxhl+IoVJWbtmJTEekaByvTAYxF6URpY0MqTtp5UY270j15wX5NrquHXb2uqI1lN51brzrRqJst62Mhi2rgeAQb9e7vfbpcx+37ZhWdJ4wopz4SXGWp9fm8mEZUgrMQ5H7duGo/pkb2zW5c2hvQgsZW5s1NuOUHl91cmaGyRfrvSoTV1ecQ8MS5b9Xvv9JWLb8G2pnFc3yPICxBlHoqwiDTOC2VewvuMzWUq4HX43RDKk20c0L2Ezkwl7gjf3s0yOIrthIQe4JEmSpEwOcDnAJUmSLCU5wC3oAKeovFIAWL1idNhWHNfy0gufcsJ2uRdIKYANauVkvL1ABgTiRKpdpYxY1nCyZrCpF8kaBc+rPn3wQeApCXSTJVedR+RaIFEaGXLVSZSBLLmPpEiWJb3cuMayJJVXB5vUZtO0GQxYoqzLg15d7vWs52Wftlm5stsUVixLsnTJ6wFgNK5PMMuSw2F9EVYGdXlzaJOPr5B8yQHtHOy+smEf/4HxHq23RXPlVcuRn1rsJazK3qPkRWnu48K3s5Eo6/KwkAQ5mvct8q6s+hCUqU7jGYvam/U+6UNUry77+eDYLFEymRzvLOQAlyRJkpTJF7gc4JIkSZaSlChzgEuSJFk6FDnAAQs6wKkCR7ZMI8P7t9c/70n7TL3xeHW7zBkX2KbgbXBkhrCu9P329dX+2D7XbrcrJU7mhVJS1mhiU6PHU31vT4smH7WJW+PEybZMGUrc5KU22XFUjpMlrw4oE8mA3ffHrfUBYGDakD1t0B4KUC2T3a6/QcepbXUlG5yJ9ShMnGvaj+t9s92N7VIAMBzVNjUOQQjDFDpmXNltgmgf1mpc9tmmyKZxl2GILhVG9DCxa//QnvbQzmQTNNN+XSqUkbHBUZsZQgPM8+/CEaQfPKOFyUut3ZzWIzDOIbDVuWwpOcAdg0wmSZIkSbIXLOQbXJIkSVImX+AWdoBTYLQOAHjuE2tZcjhs1GrFuAy7d9hRIF8OWIpxbVi6A0sU7I5ccvmncihRuG0IJA8rl7r20i5LmgwUTm40oQV8Pliy7dszbaVdPp/t5ape+xxstkx1nFu+keHoyvcK8lwp+0gnIlkydJ0HwJ+hDnZpuLVz8meWLwfG/b69XO2P5+urL9ZgzJKvS/BM9+4qZbRhiZHnoAOAlQFvq88vh54M3D3Fy5T7Gb1hfO/y06xcpmfMyqK2tZEog0wmnui55H2J+2xdEqf7j2YTJ+88zCgiJcqFHeCSJEmSEtFEqscTOcAlSZIsGelFWbGQA9xD9gMXPLHSBjZJlmzmWm7PWMD4Ge97JMdYKZIrFjIRGEmM5Up34A4/rhoyROBt2SUTAmClEFMuJHVlydJKoXEbm/miXW4seZV2kl/cCTRyHW3jOdd6TlvmOdf4urE3n5c12VOwZ7Qq9qgsSJ/GhW93v7BL3pG8zHPaDc16J9lG0jC3acjR0yVoLzeKuXcir0PTpMvjEs7z5rdxJhMvDZtj0kH5nmJlt7MXZrDfanm6LFl8XlrqJxULOcAlSZIkZfIFLge4JEmSpSQlygUd4FSbXlJAUyob00v7OJhPqrkP8gSj9hx4Ku63ER93HHg3qpOQWFKbJVlqnOC5vQy4QHEqGy+uRnA5qJ5OXd/sW0vnMWvAcew1aOdZo4Bj+kBjJ1Ga+dgCL0znmIs+94E8HTnZck9L3pnUTzNPnO1bFweB0vkIj16UNdvLpWTLYULioH21j+nB0F0fCjufW2ySsHPIta/3GInQ7Ls9CL4rXeTGxvqdfklkoDeADPROkiRJlpSFfINLkiRJymi+wuUbXJIkybKhc/ibhoi8UETeLyJfE5HDInKLiLxdRE7q0kcReYKIvFdE7qT2v+zq9ETkjSJyu4gcEZEvisgLOp0ELMEbXCkJMocQjMiGVtLdGXYX58wfvn3knmwzLuy9E28zzKB92ywTs5ZsaLN80tLEk/X62PYxHrfbo7g8cjY4DhsQabfboWOyE26vHTOkqLHBuUwkwWSo9rPFv0+j+20WW13Xel3scY1tHdY3t/HW9n42J1nlTChxvYjwnuxYL+7/7u1u6v5P68sceT2AOwC8CcBBAGcDeDOA80Xkx1VjY7SInAPgWgCfAvBzAL4H4EcAHHBVf2NynMsAfB7ASwC8V0QuVtWPTOvgwg9wSZIkSZM9GOAuUdXv0PKnReQuAFcCOA/VANZAql+D/wXAn6vq82nTda7ew1ENbu9Q1d/eqiMijwfwDgBTB7iUKJMkSZId4wa3LT43+X9aoel5AJ4A4HemHOJCAKsArnLrrwLwJBH54Wl9XMg3OEH9Kj+gIfrE/Xa8vvu++g05+jHTkDKCbaVfQzNJHkF51xRUjRm8sDvVmyXMoSuh/FuQ2sZBVpNSaIGaucxYuozDO0bj+vHheeJUbQbeSO4zoQ2NMIH2356+Xtu+drJtp3SVpqNyc3+77VFN1+eolPEk2mHX74LdsuNwgBLHxsfk3Mn/Lxfq/OTk/z4RuR7A0wAcAvDfAPwrVT082X4mgHUAX3Xtb578PwPAbaXO5BtckiTJEnK0nUw8InIagMsBXKOqNxSq/sDk/3sAfALAswD8Jipb3B9TvVMB3K1Nd9C7aHuRhXyDS5IkSQrMJ9D7YSLCA9UVqnpFW0UROQDgg6hyI7xqyn63XqyuUtVfn5Q/JSJ9AO8QkSeoaukNsDMLOcAp6os3JD+d+w9bpx2+wNGLf9ekv0dThmPMTblHxywldZ2lnmnT8Rx2kkK7HbJwjIwLWma63h/mnuwqV+4R/LKyG7mSvyN3wZ2qes60SiKyH8DVAB4L4FxVPTilyXcn/z/p1n8ClfPI2agkzkMAHiIi4t7itt7c7sIUUqJMkiRJZkJEVgC8D8A5AH5KVW/s0OzmKdu33lRuBrAG4HFu+xmT/3817UA5wCVJkiwhexDo3QPwbgAXAHieql7fsWsfReU8cqFbf9Hk/5Ys+jEAmwBe5ur9LICbVLXoYAIsqETJXpScOHXDhRXyNpP8tfDmb+c8a08K20g4y+WclGnXUuAsyZp5uRckTi61McmSqewTa9t51oat65t9jeJdaW45F1E+KgRx75R5SrOzeGQeTemvuyw5PTi8a5KDWejqhW16uUu5cg9Sdb0TwIsAvA3A/SLyDNp2UFUPisjpAG4FcLmqXj7p13dF5O0Afk1E7kEVL3cOgF8HcKWqfnVS79si8jsA3igi9wL4AoAXoxpQL+3SwYUc4JIkSZIyexDo/ZzJ/8smf8xbUGU1EQB9NNXCywHcC+Cfogrm/iaA30KVuYS5DMB9AH4ZwCMB3ALgZ1T1w106mANckiRJsmNU9TEd6tyOlhftidPI72BKsLeqjgC8dfK3Y3KAS5IkWUJyMoFFHeCkTqxsJi90pg6jbdNvCE7K7KVtnoyRJwLlNj1nhOPcvLzJTDA6g4bvb1Az6WJkfShp/VFmhsIxO+7aHaeDvaMRntFuKzO2tR7Zyfwkmr2oXl3u98auTXs9tsH1aFLTah/D1m2xna0rVsHh/bH9hfs2NnViW2HEPLKfmGwwvL7Ypr28W+xkv+58UNk8vnTaG8/bDMnJI+y58emT2s9hvIcu523WcO3lYjEHuCRJkiRE5xPovfDkAJckSbKE5ISnCzrA3X0YeP+NlTz0vDP7YT2bFaRdTGhIlIEs2Sddo+8Oydu4/bwzkURZTuK5rlwHOiSP9Y9EWG/Pks/WZSv/lsIEAunRSZT9fi03DgLpcdDfdMdhyZMmHOzqis8JlgvVWCpXiZM/h4cJ5sQrt2kv81x7pbkQI3NBw3TQQR4vfspgzjR7D7hjshTJO+APVMhq1DPHKXWOD0rF4DwBVl4fR2aNHKtmYiEHuCRJkqRMvsDlAJckSbKU5Pi2sAOcAP01AMCf3vTA9tp/+KRVUytSH0pZSdhDMvK27LtGLFlaj7eg+x5tlyj8Ldpj+YLbkIQk1Jex+wnH6pYEEpSfZH4cySwstTXaaOu2kqy50wwZXqrrBZ6X/V7sEWkylgTSo/eODGXJrhebTihnLyllLonkRp4zrjTX3WjcPtfdcGSPORq1bxvR/TEa2eOYbR3KQEHKnEH3jp7lnjud0d6kox2hF3wXlJK1m+OzGaFxrcwSFfkD+f3FfW075vFK5qJMkiRJlpIFfYNLkiRJiuQL3GIOcCLA2kSKO4ITttf/yU3rtqLWktSLz16r21MVL2WwRBlJHj7xaRT4WfK2MnIft2F5z3uCseck93kGb03u24hlVe/hRareqMdSldJ624ZlrFE/kK1co3G/XXIdB1LX2El641DGa5f0SvC+xck8vM0kSPbuhQzLuVpryNyf8di65g5HK9vl0XhA6+vyaFS32RzaR5nrDYd1vY3NPtXpuTb18uZwerlaltbykO6boVWGMRpzmaQ7Lhe8NcPnjb2eG5I3J+Bu31eJUArt+IxHjpuANQsYT2H2CnVfBq3ep84LNce3BR3gkiRJkgKaNjggB7gkSZKlJMe3dDJJkiRJlpSFfIPbmmAIAAY0RA+x1lYdAPAenkid5Ox9boh/2dPY3bxe7xMsR1hP5/gn1Diw2VgXZHtMdoU3xwkSxvoj8MVms4idFNQfsy5vBn1r2jGpPOR6ZL/ybv68rVfbiew1IPf/oQ8T4GTJ7SEDJcYUE9In2+147EILqG+9jgmWx8adv921n+1s1TLZyoa1PY5tbWxn29is6/h66xtcr309ABzZ6LeW17m8aS/2+ka9vEE3CJeHQ9MEm7TMttyxxjY4JrKh9wshNny/sjkrygJUwtrwYnu8wdjc4/AdjcKMXHJxdDhX+Qa3oANckiRJUqYxa8FxSA5wSZIkS4Yi3+CABR3gxgDWJxdvQK/xA58Emcqr9Iq/wm2cvPaJG+t6Z51WS1CnP7KWgJouzJyJBK31Sklq7cR1CDZ4F+K63Ofj0OfxF9cpRQT1WuwJ8RJM2/HLWaVZkiu5gbeHBoQhAy6jxtBk4aiv/ArdFKsr9gz0N2v5sd/nBMt12Sdo9tlQthAjH8fng8MBxnSxONsIYEMAWIq0EmW7+3+1HMiSm+3SY7Xca613ZJ3Xi2tTLx/ZYImS6th81RhSnMCIy1FWE0dsOuBQAJc9hWXJqeqxyQAAIABJREFUXvu1msdciGH76HkHbOJ0RM9IIRtMeNDu/VtW0skkSZIkWUoW8g0uSZIkKZMvcEdpgBORnwLwqwCeiuoN+isAfkVVr51sPwXAbwF4HoD9AD4D4F+o6o3te3QoMJy8l7Msud95Gq3R++ka6Xgr7Hnlp7YnJeDWb9cLX/3bWt569tnWY40zNUggwzWTGE/3gvLqoBqJkqUV8i7kfrn99W3+k/aDDl1HA68/3lMj2bLJPkK7Npku7IfjLBirQ/bMo/JqXV5xevTayoi21ddnZaVdeqzqkUQZzBvXc/cU14s8NEsSpUl8bM6T/TxR4uPNIXtXciYT59242S5f8vn0HpHs+cjbWIbcGHqJsi6zd+QGebkOncfrkO6xIXtRjuJnoossaebQc+YKfv6ie7chm5t6gbdlwfQQDTCN571Dlha/t7aE5rZfmoHeOAoSpYj8AoAPAvg8gOcDeBGA9wJVTi2pvpmvBnARgNcBeAGAFQDXicij592fJEmS4xHV3f0tA3N9gxORxwD4XQBvUNXfpU0fp/KlAH4CwAWqet2k3WcA3AbgVwD8s3n2KUmSJDk+mbdE+WpUkuQfFOpcCuAbW4MbAKjq90TkagDPRZcBTsfA8DAA4Ijs3159spMl9g/qnyEsUa5SuV94h2XPR5bUxAV9SyBfGInDJ+3tmKeXscmfA12jqLxzPdYV2z0d/d7GJjC5Lg/69nywZMuS1gpdnxXnWceS4wpdtxVyc101dawuatvU2wZc7tvPNuiPW8t8bn1AupggcnTCBv9z4uU4+fSQzu9wGMmVXPaSbxSATetdG/Z8tEmU2+sAwCbJiixFsncky5CAlfvCuQMRY56DfvtFaMiNfd7WLu+XkiCzk6vx9iw8vNH0dqU2DEuxXg2353D68Y9X5i1R/iSA/wngJSJyq4gMReSrIvJaqnMmgJta2t4M4IdE5MCc+5QkSXJcoXP4Wwbm/Qb3A5O/3wLwJgC3orLB/UcRGajqvwdwKoDbW9reNfl/CoD7/EYReQ2A1wAAVnIMTJIkKZFOJvMf4HoATgLwv6rqn0zWXTuxzb1RRP7DrDtW1SsAXAEAsv9hivHW3G+1ROmDtlmW3L9CEhbJFT7FJN8TI+NhxUGkpX5S2UgZtt4o8soi/HHGZlIpqtfnfXFgtJPXon3xzpwH4Jg6Pua53UiydbG/6FM9lln6QRmw3rDcntezxOi9KHlb1KbvJMo+57ZkWZJzafbicxjJlaUgduuZx9fKtmEJmOVxc96pPPQSJcnELDeyxOjnaWMpkuVGE5jt7qnIC3JckP46aWeNHI/Susl6HXbMF0tXga+BHxDsvHV8TJJi3TGH/GEDibLr6VBz4txzab6nZLIf53eZ49vcJcrvTv5/0q3/BIBHAHgUgEOo3tI8p07+H5pzn5IkSZLjkHkPcDdP2T6e1DmzZdsZAO5Q1YY8mSRJkjy4EJEXisj7ReRrInJYRG4RkbeLyEkd2mrw9xRX7/ag3vO69HHeA9wHJv8vdOsvAnBQVb8F4EMAThORc7c2isjJAC6ZbEuSJEl2wy5j4DrKm69HNfPWm1B9x/8+gF8E8EkR6TK2/BGAv+f+vtJS7+Mt9T7dpYPztsF9BMB1AN4lIg8D8NeonEyeDeBVkzofQpW55CoReQMqSfKNqETm3+x8JJ8aBM2sJKsmNIDWk0u5d/VmuwTbs6Tgwsx6eOgC7e0y4QKtdn2z9p/2Nl2xtqC4M8ZSN47Pm21E5619dWN+vV5gY+liz/PLHPrRY9uen7fOHIfXx58zzjTRjei8j93FHtF9OApsuXyvensa28ZGQXJj777Pc7OxrY2vu7fBlRKKR5j7ILgnev4+DK49r4+SkXtsqEa9fuzOR2Rj5QxB3m7XM/ZWba3nBw9jtqOyvdfia7WVGGYWU+cuuURVv0PLnxaRuwBcCeA8ANdOaf83qnp9h+Pc2bFeg7kOcKq69er4dgBvQWVr+58AXqaqfzypMxaRiwH8NoDfA7AP1YB3vqp+fZ79SZIkOR5RHH0vSje4bfG5yf/TjurBOzL3VF2qeo+qvlZVH6Gqq6p61tbgRnXuUtVXq+qpqnqCqj5TVb84774kSZIke8qW6enLHer+ooisi8gDInKtiPz9oN4lkzrrInJ9V/sbsMizCUjV9X0mobKtErmL+9CACJaD/sHZ+7bL6xtjV689G8M4CDMoYuQbL+NRmT43u1BHkk+JWK60MovdEEsu2l6t8zkIXcJNnbiNcfOXdjnL7y86n/4MStChYuhItL4QKhJlqhh3kBH9tijzxshda5NQOJDNSphrUriPWSrum0wk9Iy6DCXmtAdydN9cRNe31qPYazDybdqn/ismUY9kTZMcvSDzRmXfFX5OR0Fy7zm8vz1MRG6g5Ssm4VqtiMhpAC4HcI2q3hDVm3AVgA8D+AaA0wG8AVVI2bNU9VNU72pUb4W3ofLE/yUAHxCRl6vqVdM+wOIOcEmSJEnIHBTKO1X1nC4VJxmoPojKRPmqKdWhqi+nxf8uIh9EleHqragyYm3Ve507zgcAXI/KDDZ1gMsJT5MkSZaQvZpNQET2o3rTeiyAC1X14M77qvcC+DMAT59Sb4RqdppHi8ijpu13Md/gpAcMqgwmJ5P0uOYyVZh534Kh3Ht+XfC0OoRjfaMWBtY3yfuslM2Bpbug+4BTUCJPMqc82Awb7V5lRU+ygvS2RcMTK5LKAlkWKMhjfG4a8+PRMQNZc5YfpGVnz+nZMUr727UnK3/O0jbjqdh+3psS5/R9NY/ZfoZLkq29J+syP28u6YzZC0t6/cLN2wtkeJYlzbPT8af72BzH3pR8OvrBPTn2nrksHXa8P7oMJv57ip+xra+mY5G5RERWALwPwDkAntV5Ts+YnXyKqXUXc4BLkiRJihzt8W4S6/ZuABcAuHhWV/7Jvk4GcDGAz06pNwDwYlRJQb41bb85wCVJkiwbVZzA0T7KO1HFOb8NwP0i8gzadlBVD4rI6aiS7l+uqpcDgIi8HsDfQRUzveVk8noAjwTwsq0diMhLUU2h9hEAX0flZPJaAE8F8NIuHVzIAU5EsLayAgA4sFJrAhzADViZgq/1M59ey5CcYBYANkiK5PJmUAZiuc4GbRa8wmg9Bxz7NtZjbHfBroyRGAvzwUVee40gY5OotxmQCthAZl8vSnJrJVLb3mwLJTl/Qnb+BeADblvr7FK6BAoJeAPRdZbvsoYnKpVNsLyREW0b66nM5chv0d6jSgtdP4OZJ62DZ3F5X+w56p43Kptk6+NYCjVzCY7bpdhyEPq0Hlfwc7G5lWyZ97ODfe2C50z+Xzb5Y94C4M2obqs+rL/HLQCeP/n7PgD3APgLAP9YVfkN7jYAD0c1O82pAO4HcAOAi1SVJ9EOWcgBLkmSJClztMc3VX1Mhzq3w5luVfVqVE4p09pej0r+nJn0okySJEmWknyDS5IkWUJyPrgFHeAGAjxsYnvbT3Y3P+Hps57ePvP3+ia52Dp72noHG5xPODsyCZprSol5I5d/mwDYuUoH9XpBFo8S5ubvsf3LHlOD7BYmGbBLUrtJy2zj3CTb3KazwW0MuVz3ge17fKk2nQ82bxvyZKGFjA+R3a6Zxnt3RJeklGgm2uaTEG/h7TqDIKmHmeS10UaoXNfjCYJXXCjOamDvNAmzu4aEIGa3ds3ShLRhm2AHJXuaBPVKtvHos5X6yeEiG9sTnnZvf7ywkANckiRJUiYHuBzgkiRJlg5FIY/sccRCDnA9AU5cqS4ehwY8+8dONPWMWztLXUOWIa1WZsIESFJjGa6ZkLgmynTh5y8Tk+mBJcY4k0kUGtA1u4b9RddRGwqasMzkk/byMrv/86le37Rt1kmWXCe58shmvf4IyZKHh/aD3jPietxn1lLX7UHH1Amlg5o0K1732aWAaeaBlGC9W5ZBsJ4nsVt17VlurFefwGUnN3ImoH2067FGd7jPREJtTAYcJ3t3ufV8uArPs2ay27SHGTQCQqKwi8K93+W5LiQLsrIkmyScGYDDMEYd+8m34eHJvT9vaX0ZWMgBLkmSJCmgKVECOcAlSZIsJTnALegA1xPF2kSavOjptSzZyG5BMglLkcY7cmjvgs1IlizcLCw4xAmRbZtIloxkSGDnc72V5l/revNbOUjb17t9RXPijQJPR8DJlyQ/PkDl+6h8yEmUGJL8ODpM5fX2MmAlykXBaF0r7WUA6K9tF4e9unzPJEk5AGw46fBkuo7Gm5eu28BdbJYfx4V7IqLLXGh+3+MgQ4jZb+dEx3GfuyT69iaBaM7Gnsl+4jK7BJKt/z4zfaPy+mTBV8/xLQO9kyRJkiVlId/gkiRJkinkK9ziDnBb0gDLZiMnN0aJk7k8dG5LLKNFgasl78Ze4B3pg7Ztm/b2vUaC5h3Kkq5+JFkaachpsV3mEitKoR32BTivTFMmz0mWNccuUpzlxyFLlEfooEvgZ8afoav82ncZxSccwX6zvEr3y1owD+Asdp1SAHIkEXrPSyN7R7Ik76vw7ET369jLr/xdECQ8KCFkvGBPSZ+MoR95iBaC03l52yk8A70bLOwAlyRJkrSjyBc4IAe4JEmSpSTf4NLJJEmSJFlSFvIN7uQTenjW0yr7AU+uuTG0NpaNjdpOs77BGTXIblfw/7dJamN7GtvaonIxcXKYlSS2uYUuzIF9A4gnDI2SDpfazJIIpYQJtTBlypQR1iqwDHa3WWAbHIcQ6KhZd8LI3DtHoU/b+56eiWTsbkQ2v7JFsR/0edx4dKZ/IF9DjR2ymw3OJFvmyU85BMNndjHZYHhbHKZk+tDWoQz0BrCgA1ySJElSQovOX8cLOcAlSZIsGelkUrGwA9yWhMFu/uvrVn45shElTo73yxlHeoHEuDpwEuWAZcl6ByuU5dZnHumZ+brapUifDTx0s++YCYElF5ZmjRTjpaHAXbyr63iUpLbnXL17Zi6yej3PRbafDnq/16Aoc4dNnFxItpwYWEWP1PHGXGZ0x800ZVskm7udjaLQACqzDOgznETPmN3XLCEyU3db9Y3nbHR900Bt7/4C1i7D5wtcOpkkSZIkS8rCvsElSZIkMfkGt6ADnCownHhMsvS4vmlf1SPPSb7wfp42VjIGJI+tkAzJ0iMArKy0S5Rc7vVNk06y5LjhCjY9g8O4IDeOgm3DglfpyCRO5vZ8fN+f1m7aBL4++TQtj8g1bpWOs5/O4SnuGIewhql4fW14pL3eUhBMRsZzyLnTsUoaH18Pli695NNVymQiic9K4N6dtz3bh0lUXJJYZ9FPu2T+aWxrl2y5bz6TCdfUyGs4migOAMZbkryTPuMWxw0LOcAlSZIkZfINLm1wSZIkyZKykG9wqrotq22SLMmSJACsb7L0Vq9necxLGf1gbjYO2mZJEgBWSMtkj0qTULmjZFOa68rIgrw+kCgbciPJHJEs6T1MR6EXZccgVKKcsLour1BZB+07a7Sn5buF50KjW3y46nZCy6XExYuA+Ay+JNnSfHDo1595n5eJWUI2cjJJgq5NryQLBkSypCjfkwV5XtolPTZQNLw9Z5FSO7wBdY01s/e+n0QuSmbA+/aJ03lh3FipGegNYEEHuCRJkqRMjm85wCVJkiwl+QaXA1ySJMlSkgPcgg5witqexLakI84GNyR7kgS2goa7ep/LQeJkZ4iwtrZ2u5uX3aObj20KPk9wFAIQ2tZG/ny02934PHm7Hdv9bHvqs7fBoR2TKtmfd7Yf9NsNI2yTbNhOaXmV9n3fsL6g9/QO2EZjmvDT2OCOtK8HbBJjLh8LOIly34VJ9OmzDahM8SoH+vZKrXF4BpVNdh9vPgrsqrPY42xWkzgTCScn7gV1fFiAt2cH1fYGfw6Dz2ZOqE+cbpaOzUgmIi8E8FIA5wB4OIA7APwJgH+jqvdOaRt1+mxV/R9UrwfgXwH4BQCPBHALgMtV9f1d+riQA1ySJElSZg+GvdejGtTeBOAggLMBvBnA+SLy46pTp/L4IwDvcuu+4pZ/Y3KcywB8HsBLALxXRC5W1Y9M62AOcEmSJEuGYk8kyktU9Tu0/GkRuQvAlQDOA3DtlPZ/o6rXRxtF5OGoBrd3qOpvT1ZfJyKPB/AOAEs6wGktpbG8dnjdXlFOkGwyMwRzsfll4bLE8oeRFWn9yLj8uzaRy32QYQQARqQLdpElNxvt2yXGaH213L5tHJSBwlx1hJeGbJaTdvfoUrJmk6CZ2q/RPXBgZI96eFQ3Ojw+Ybt8ZFSXGzKkkSg36vKR72KusORoPjiFNrBE2fMS5T7aVrc/hcIu9juJch8tD6KsJj3bpheE3JTnMmxfz6qkD9Dl+5CfUTVmiFjCto8vhxlwpaDDjq7yqwTl5vdHjZHu2QzgvqfGZvK+lhelPQgTcIPbFp+b/D9tDoe4EMAqgKvc+qsA/GcR+WFVva20gwz0TpIkWUK2YuFm/ZuRcyf/v9yh7i+KyLqIPCAi14rI33fbzwSwDuCrbv3Nk/9nTDtADnBJkiRJGw8TkRvo7zWlyiJyGoDLAVyjqjdM2fdVAP4pgH8A4DUAHgrgWhE5j+qcCuBubUbT30XbiyykRKmo5TKW907YF0f7sycYe0p6L8pIfeBzPPI6Hv1OYLkudhSyUmQkV3qJkmWJ4S69KLt4R1b94T5H610/WT3hMmJYmuGbsmdkpzgzTJ/qGbmS+rzmJLkTSELepE5v0vrDI/uIHPruHa39nzthNhXy8Dz5cXXZnY+T6R4/gT43n4N93ouS5MtV401M5UImk16s4ocYOVsD6dBh5Eq6J4wk6HZg5dNuHpVxQvSYLvJlo0qgX/Z67d8RrftoYQ4K5Z2qek6XiiJyAMAHAQwBvGpafVV9OS3+dxH5IICbALwVwE/O0NdW8g0uSZJkCdkriVJE9gO4GsBjAVyoqgd33le9F8CfAXg6rT4E4CHSNOhuvbndhSnkAJckSbJkbHlRHu0BTkRWALwPVSzcT6nqjXPo+hY3A1gD8DhXZ8v29lfTdraQEuVdDwD/9QvVoP4Pn1ivHwQBwoCTUox3ZHwcZdlqWK/3XoMjn6G4befujom9KPk4zgsy9GhkibLUHlSPjhnIkI3+RP10D0MkS3Z/aOoy/wJbIUnNOUSaa9qnvq2wtOyOb89n3Z7l328cvL1bp48F99wab6Lyw3/gh7fLHMC94iTKlUDGXwmCvv1ylEyhBPfA3K9Og47kxzEvFBIaR+3L89a1e2WKmXiu1D44fiGPdBe5srGPbS/Ko+w26ZgEYb8bwAUALi65/HfY18kALgbwWVr9MQCbAF4G4C20/mcB3DTNgxJY0AEuSZIkKbMHw907AbwIwNsA3C8iz6BtB1X1oIicDuBWVNlHLgcAEXk9gL8D4DoA3wBwOqp4t0eiGsyq/qt+W0R+B8AbReReAF8A8GJUA+qlXTqYA1ySJMkSsgeB3s+Z/L9s8se8BVVWEwHQhxVjbgHw/Mnf96ESHP4CwD9W1c/a3eAyAPcB+GXUqbp+RlU/3KWDOcAlSZIsG7uLZet2CNXHdKhzO5zwqqpXo3JK6XKMESrPyrfuvIcLOsAJgNXJKfvjL9XGgpee5WxhgdZemnyUzQBspxK2czkRPZr0sOQiHyaZDexcvm+x3a4uNxInRy7/gZ2tuY263PGzdQ0TiLC2hrro3dXZXjIKTJ8DZzCJPvfnbv3aDD198PLVb9Smiqc+9vTtsj+HJntJEErjbXCRbbtEaJel5o1kyxwOEJjAjG3O26zMd8HOQwuinRl7HJw5raNN0mRmibIANVKzTH+a9tYi9+AkvSiTJEmSpWSuA5yI/ISIfEJEvi0i94rIF0Tk1a7OPhH5LRH5pogcFpHPiMj/Ms9+JEmSHM/sVZjAg525SZQichaAawBcD+DnATwA4IUA/pOIrKnq70+q/icAPw3gDQD+GsBrAXxcRP4ezwNUoid1dgaODHjvjX1Tb0gX6RVn13pUSTrQwK3cuPUXM3eQxKDtdfw2jdrbJuH+on01Jc7pffM3dpfjlJ6FeT4oRdfzQL7UqBIApXr++iwrX/hrkl9XTzbbznv8Q7bLNvl1Xfbz+HUNB2DMvcf7KrjvC7vsR3Ijy6U+2ZC1ArVvKLnvB5Q+f6hw+hCGIEm1fXZsZ0wYhGxdlDiT0/HKPG1wL0HlLXOJqt43WffJycD3CgC/LyJPBvCPALxaVf9PABCRT6MK6LscHV0/kyRJkjI5vs1XolxFFZR32K3/Hh3n0kmd92xtVNUhgP8G4EIRcXN+JEmSJLOQEuV83+D+CMAvAvgPIvI2VBLliwA8E8BWYs0zAdymqg+4tjejGiAfj3oqhBBBPecXZ1840WVmGJKH08duosS6VG3kvKBe8BTyTjSJj6lNIYnxqIPXIWDlQ+tVFt9ZkURo6gT1G+077re0v2nrZ2UW2Stq01HVLHrWLi0b95jFT/1Vvfzss35ou1zyQC5MkxgS3i9x4p9QvoykS9/R+Jd8IFe6xTAJu3tPihI0lxKFm3s3yH7kZU0/H2LSztwGOFW9aTLVwQdQTYMAVG9r/0RV/9tk+VRUCTQ9nac/SJIkScoojh+7col5Opn8CID3o3oD+yeopMrnAvgDETmiqu/e5f5fg2reIPRXT9plb5MkSZaYJZIZd8M8Jcp/g+qN7WJV3Zys+3MReSiAfy8i/xXV29vpLW2nTn+gqlcAuAIA1g48YnsG932FRLAsMbKqOArmAQOAP/tSve2nz2qPMva/jGwSYy63S5d+H12kx67M0rzrMbvIol0pSocFb7odH6fk5RZtO/HRdfn+Hc/8sbhI/QBFSZQbTTpcoIbcGEiRpfuwS3A3e3h6T2fWBc3XhHnEndwYeDR2vafCciGRdC/onLjvD+/N2kaOb/N1MnkSgC/S4LbFZ1HN1vpwVG93PywiJ7g6ZwDYQHNq8iRJkiSZiXkOcN8C8BQRWXXr/y6qaYjvQpV/bAWV8wkAQEQGqDJEf0JVo2mMkyRJks4IVHf3twzMU6L8jwDeC+BqEfk9VDa4SwG8FMC/U9UNAH8pIu8B8LuTifJuQ+V5+cOgaRKmoaglv/5KvX5tYF/KOQic3/BZvthwE4ut0LxvLFd+e7P+LfCSJwbzv6EQdF2QNUv1dsNe3aJH8zi7lStnaq/x9T1qiP+tGbmFBr9Jo/Uz8tEb/3a7/NNPfmRYL/QuLNzHXa5JMwFDe9nKku3rG+153/zFUGgTe+nGN1XkOem9KDloO5IeS20i0gY3Xy/K94nITwH4VwD+EMA+VPMAvRbAu6jqq1DNH/RWAA8B8EUAF6nqF+bVlyRJkuMZRdrggDnPJqCqHwXw0Sl1DgP4l5O/JEmS5CiQYQI5m0CSJEmypCzkfHCqwJGJq/+J9DPFy9I8pxVH/o9JhB+4eaOijBYs1a8M3LxixoZGCY1Jn/fzrGlgH2S353m735t6HW0fs/Bg1f5L7uqGB765uwOV7Gk9euRMBwq/NXl/87S1qc9I3L7tz770ne3yJU/+flutw7Xumoy74LEfutxHdrdmRh56Lunhm8f9zpgrGiRR9vazXr/dVhd9TsDdBsE98WB9DveShRzgkiRJkjI5vuUAlyRJsnSopg0OWNABbgjg7uFEoiS5cJ+7ohooO32TqNSFCdB7/T6qdtHZtd55ZMPqBdGNJDy1vTvOiBqZBM2B27OnQ77aIvOUKBtSSLCPouv4LMfZJbvOmGKkw377+tK2Yhte3qVEGUqRbj3X41AJWn/1F78DhiXLrnMERhTbhM8YVSkcX6N6MySLbp2KbUKv1y5/snzKkiTg5tsz89txOTal/MxTqklYPvmNdKnwLOQAlyRJkpRJG1wOcEmSJEuJ7lmqhwcviznAKTCcqCaHh/VFXN30SUzbPSef/WMn1rvyGUYoQfLmsJZmNjfr8oqTGJSyqURzZ206b01O+MxSpM1wYttEkuXR9LaM6OoZp4Ea3JCQuiTdLcwXNguzzGUWSow9LrtsdUaKZC/KYF+NbTNIlJEsOR5ypbiNUj1u445/9Ze+u12+5KyHtu56GV4k7L1C0qGr18VzspSVhBMxWx9b9902RVPP6XIqFnOAS5IkSYqkRJmB3kmSJMmSsvBvcBc/pS6vrtiPs0oeln2SFdfXY/dEM4ccyZXG467nvaBoYz+QGPxcdYEsyccfj50UGrTvyiy/6LrMAeel1Ehu1KBOY1skRcbOr50+W0nVufbmO4JG3k2O9eh++3ovUZrg7v709Y0+xJtiaH98brifPqk0S5EaeHGydOnoIld6jvVLhgRlIA4uN3H7fm63wHOSvzMabQJZ0sxHN8OZSolyCQa4JEmSpEmObznAJUmSLB0Z6F2RNrgkSZIl5GhPeCoiLxSR94vI10TksIjcIiJvF5GTdtJPEflVEVER+X9att0+2eb/ntdl34v7Bjc5/2x7YVd+wNqwTIYByuBQssuE9iP308i4A5M9bkDr+41wBLRuK7nLjzvUY3tYyRU/Wj+L+75/GOwEru12zKIJcQYX8y6ZSBou3dHPO7a7sc0KAPpr7dvY7ubbcGgA3Yecs3vNdY4jUfqy85/iI7omG9R8aMo+NIGWx5u8oS6qO2lBOMHVNx7aLl961im2SWRXLRDawzrWkQ71vOu9TZw8PSGy38f/3965B9tV1Xf887v3JvcmQEwQUQF5FWuFoUBgplSwPJzykogPVEBLxRGFOozVoVSE6QCC4OAgtSo+iqMzMqAgYKhWngE6StQQICGlOLzEWChgwisQSO799Y+1L+e31tlrn33PuY+cnd9n5szZZ++19l5r7X3OOuu7fr/fisz/M+dN6WT+n7vuDHE68DjwBWA1sA9wDnCIiLxDtSqad0BEdgXOBp6qSHZjcV7Lg3UK2L8dnOM4jpOlCzu0ibJIVW3stjtEZA3wA+Bg4LYa57gMuAJ4G/n+6BlVXdpNAV2idBzHaRjjjt69vDopEA9tAAAVmElEQVReI+7cxvlt8b59p/wicgKwEDhzAlWbEH05ghsUmDcU7kC0ZlIq6W3oRgsxmzVlgGgdK7tfcqnitepya1WlrgA5M/torTmTv9IUPyM9tq1bl01HLSTTOqnqlg1ekrsHFdfP5UkjSNy80rgGRHJjZhtgcKT8mJUuk+uMmL+Rc4yEPSsT9Sb9PNCFPZyVIm3QcPuVeDVxQ1lnn7fIHcL8TIy9Gl+ohjvB4hXPJqVrFSKVLydKLhpNlXSYM/lP80QSdiZdet/y8ifTzgwZmRxUvD9QlUhEFgBfBc5Q1TUdpNlFIvISwfflHuAiVb2+TmF8BOc4juOUsY2ILDOvT1YlFpHtgfOAW1R1WYdzXwz8Dvh+h3Q3AKcBhwMfAdYD14nIR+tUoC9HcI7jOE6eSYpF+Yyq7lcnoYhsCfyUsJrZSR3SvhM4EVioaZSIBFU9Lcl7HbAUuBD4Yady9WUHJwLDxdjTSgSjqURpo5JEQYztdt4iMid/pNZ3uUCs1eZS5Xki9WSwnlwoGVmyMsKI2W+VqTaDvUkMcNxr/rqST+5e3bjij3HCISs3WitIIzda6bEtXetKW1gZMrlvszNrEdpndzCVxyYoaaU/ZrPMHTYxwxkyz8qs5GYPmka1lpfrrZQpSXvUsbZsM6ZrfV688vnXtt+71zwmStY6Mk2XkRW7iWQiFQ/iVEmRbWUz27lnZWyaVhMQkTmEkdauwEGqurpDlm8DlwOrRWR+sW8IGCw+v6yqr5RlVNVREbka+LKIvFlVn6i6UF92cI7jOE4F0+ToLSKzgGuA/YC/VdWVNbK9vXidUnJsLfBZ4NIa5+lYQ+/gHMdxGsYkSZSViMgAwcT/UODoCZjyH1Ky71KCEclpwEMV1xwCPgw8rqpPdrpQX3ZwC+bC+/YLEsiYuYtt66dl8o9VeDZnJWErcSSOzbHjp7EUrHQ8nXr5oNJpu8Z+iC05c47mVdae8f0pP1d6zJJrpXYrudZ2JEXagMZDc+JMkdO2lSWtdBlfaJ6xfh0xFpGzrXVk4j+dc9qutMajM1E7p4EEzGcrPY6ahk4tzCLJ1NzTQVOademXKmdtmVtPLhw0m60TXr/iJZMmvtD799qydZmslWzGUpKkfTNyY9V3tMpCM5eHinP3ShSgWabmGjX4BvBB4AJgnYjsb46tVtXVIrIT8DBwnqqeB6Cqt6cnEpFngSF7TESOB44Bfg78AXgj8GmCa8HxdQrYlx2c4ziOU81o5yS9cmTxflbxspxLiD4ihJFZNxb7jwLbEiwutwbWAcuAI1T1xjon8A7OcRynYSjtRneTfg3VnWukeYwagoSqHlyybylB/uwa7+Acx3EahhI7+2+u9G0HN37v4gUGk7mx6AZHMUZe2xpLJi/qPBPtc30Zrd6a8rf9h5n405ebw4oCLFfMc41lF1bN58nNteXm2arLSU9EJv8rE+tgGyDZzrVF0UbM3Bokc3CtiTM7zzY3MfkfNp9nZcz/B9M1UqO5tswcLfXIR7OJz2CLbefdomcyuaqYm2Xns6w7wWDyIGfdCTQzNwfx/Jx1IaiIzXvtylYElQ/tPVyaptJNwMbPrhk4ORe9pMq1YDKpCgRdWe6C0WlyE9iU8UgmjuM4TiPp2xGc4ziOU45LlIG+7OCUlixWZb5rpSLJyBLpklhjuYgnXZaz7FyVeaoikdjtrFxotlMp1RzLRXZplzU7y591gzpXUSdqTCRLpkGQc1KklSEH4zw2+shWQ62CWhlydiJRWllyaMBKeuVlhrxMUtesO7eOnt0/kN7r6KEof/bbgzibdBl3gtSdYbY5xew6kVAgli+j/bn4OkQP/I9XlEuZx+/TKml7VJJMG2Sky5CnPH+v5vhV685lo6eQSpStz4NloUy09ymBJtCXHZzjOI7TCe/hfA7OcRzHaSR9OYJThQ2FwGxlyNSK0g7dIyUyWostzlNPkqsXBLlK4sxJkTkZMi2DLedYppxphJGerShr1i1HlZxsP958v4nVmltzLQ2CPDin9NiIudfzErlxjpElrRRpgyPPSvIMRJaTZn/NQNDdEEmRZn9039JrRmu72Wcl0ijTK5ksGUvD5AGJ0pn8ttkGk++YdUBeH0X6rpg7kORzCVfe19p+x5s2RMd2266VP44gkze9rBO9pCqiTjePQV1ZdMA8fOPPZFsZfQDXnx2c4ziO0wnv4byDcxzHaSIVfoWbC33ZwVmJ0upElcP4mhKSGsklkuTGMjIPMGplzRoyICRSpN2ukAvtdXKyYmwpmbeirCtr5iwiqyy08utotbZvvf/xJJPVmnOypJEh08DJxlF7gZEetzBa2UjytEdO22Z7yMqQqUWktcw1rdONLFnVhql0XpYnktpSaTlnlWrKrKlzeKR/2hPmZc2stWXFunMbMpaXNrRUGkcxZ/Ju9w+ZYv7mqdhSc/nTre3jF9ZztrdTHtGxLiwqU8vJOseyUip1fs/cjBLcyMRxHMdpKH05gnMcx3E64RKld3CO4zhNxOfg+rODU4WNhfgeLdLYJjnbObh2s9pO5MywkymFSAMftbNWY/lJkujZy5jc1w1inJt3S+fTRmvMwbVFJakh41eZSi9Zlcy1vZYoXRU04wKQmXcbSbLPN/Nuc2pGJRky58jNuw0kN7ssaMREqLvKcp35zlyEk6o8dYnnpoz5f5Xxe2YObyydPzJl22gKGs3BJfODs+13zKQbzhRnsKKY199bvv/YfasimZjt/Kl7j3JC+fxetTtCWSQT9Q6OPu3gHMdxnE54B+dGJo7jOE4j6dsR3LjcVne9Mc2sB9c29rdrZ+XWgEolSrudWyxqEix2s5FEMtttUUkoT1dFnQgOWRkSYvP/wdmt7YHZcbpMJBIbIDln/g+xC0DO/H9WImtG8vZA+X2vK0lWSY9x1Bcp319xryKpPHPe6gg0mWsm5exm7b7IBcF8iNX5vKZvf4Dsv+3UtWAyxyK5f/U33BNf8737lv9OdJQHS9Ccf0dyLP6O9hIWxSVK6OMOznEcx6nA/eC8g3Mcx2kmPoLryw5OURMlxMgvaTpryRVJQ1aOyv/LycuAiXVjJg9VElSNdHX/f+XStVlemYTWkjQqS83gsTevNLJkm0WkiSJhpciBqmDJrc9bmNNZKdJaR46k67QN2u3yqCSp3JizprXtMZret4zEmMsf8nQ+d9XzYQ+NjvUmN1b95KWRTUrPVZE/J9O2RRjKyHWRxJmcazKNBXKyc7o+3s+MZLlo3/Jgzd2Q/n7YE1o5OprhyM+kJNMvjqUvOzjHcRynCp+DA+/gHMdxmol3cP3bwdWZP42cnvNGTLWuEVk6JdfOrRUXO2AneXLruWXWoyv7PE7OijOVdQbsOng15aT/XPFE+QlmzzOZUolyxGwPl6bbIslipcjhaNtcMpIe4woMZtZpq7rXsSN8udRddQ9sdnuuNqk8Y9GYc7avuk42GDh5ctJhXafzunRj02BlwTHK5cpuzlWXqufDluFny1vhnxft2/rZrAzcbraj5yMNHhAVyKSLGjT+No8afXt0dDx9WgLv4Pq2g3Mcx3EyqLoVJe7o7TiO4zSUWh2ciOwgIv8mIneJyEsioiKyc0m6ERG5WESeEJGXi/R/U5JuQETOFJHHRGS9iNwnIh/ovTqO4zgOEObgenk1gLoS5W7Ah4C7gf8CDsukuxx4N/BPwCPAp4EbReSvVdWGOP0icDpwVnHO44CrReRoVf15x9JkRt+pHG7nvdTcsLoLd2bNwJPP0byZTVcxLxMvPlo+71Y1L2Ox8wCDlQsrlu9fvGKtSZRMjs2e39oeMI+LXZR0IF5c0prmb2Xn08y82XBq5m/y2Pm1OAhya386R5MLrlvVnnWigliz/DTdaGY7NbfPzZUduU+rrdvn+sqfiWdfbM0FLXkw//+015+nutJO1uR+CufQujHT76o8mTw/W77xte0D/zw+9rot7fxc+e9CWv7IPiCziG4aeH3DxtYd3pibTJ3iTkpEjgWOB/YDtgUeB64FvqSqL0zgPJ8HLgR+qaoHJscGgH8GPgW8CXgQOE9Vf1Ln3HWf4ztV9Y2qehRwdaaQewEnAJ9V1e+q6q2ETvFx4DyTbltC53aRqn5FVZeo6qeAJcBFNcvjOI7jZNFJeHXkdMLi618AjgAuA04Fbi46po6IyK7A2cBTmSRfBM4Bvg4cCSwlDIaOqnP+WiM41Vp/Bd4DbAB+ZPJtFJGrgM+LyLCqvgIcDswGfpjk/yHwPRHZRVUfrVMux3EcJ8PUy4yLVPVp8/kOEVkD/AA4GLitxjkuA64A3kbSH6WDoWL3EhHZjTAY6qj2TaYV5R7Ao6r6UrJ/FaFD263Y3gN4BXioJB3A7kBlB6e05EdrLjswHP9peOWV1g226XJrqbVdJ6MJVv23qS1/ZvLkzlWFjchRtW7VdSvMrbFy4yxj8p8GQTY6zYi5zhwjF85pkxvVbLf2W+kxlRRz8mPVmliWnNyokRQc57HyY05uTCXKDZHs3dreaGTJRftUBNM1eda/Wi6bhx12s/VhxCyMduSe9R6Q2x9obf/+5dYNGapoT3vMitZt903K76n14qiSB3PH6prv95KmW+xzdOeD6dGWfPnuhYPpwVLq/M5o8vDaZ+fKu18NG+umd94s6dzG+W3xvn2n/CJyArCQIHNeW5Kk58HQZHZwWwNrS/avMcfH35/V9ruaposQkU8CnwSYM3fL3krqOI7TaHSmDEUOKt4fqEokIguArwJnqOoaKf9H0/NgqG/84FT1O8B3AOZvva07eDiO4+RQpt0PTkS2J9hb3KKqyzokvxj4HfD9ijRdDYYsk9nBrQV2Ktk/Xog1Jt18EZGk4Gm6LM+tffqZxT/++jrgmW4L2xC2YRrbYH1mu2zYPk1Ma/3r8t1fT+vlum6DjV0e2wTZJJ+Db/5q2i41Xv/W7++rz93II4u36fG8IyJiO6rvFAONNkRkS+CnhEfnpKqTisg7gROBhSWd16QymR3cKuB9IjI3mYfbHXiV1jBzFTAM/Bnx0HP34v2/O11IVd8gIstUdb/ei92/bO5tsLnXH7wNwNugrP6qesQ0Xn8OcAOwK3CQqq7ukOXbBJey1SIy7oc0BAwWn18uDBJ7HgxNZiSTG4BZwAfHd4jIEPBh4KaiwAC/IFhbfiTJ/1HgfregdBzH6Q9EZBZwDcEX7ihVXVkj29uBUwgd2PjrAGD/YvvUIp0dDFlqD4Zqj+AKpz6AfYv3I0XkaeBpVb1DVe8RkR8BlxaVfrQo6C6YzkxVnxKRS4AzReQFYDmhEzyU4GrgOI7jbOIUvm5XEH67j1bVpTWzHlKy71KCwe5ptJQ9Oxg616StPRiaiESZOnh/s3i/g+DzAEF7vQA4H5gP3AccoarLk7xnAS8Cn6Hlnf4hVf2PCZSnVAvezNjc22Bzrz94G4C3wUzV/xsExe4CYJ2I7G+OrVbV1SKyE/AwIfrIeQCqent6IhF5FhiyxyZjMCRTPMfnOI7jNBAReYxyw0KAc1X1nCJm8aPjnyvOdTuhg0tDdQ0CZwInE4fquqZWGb2DcxzHcZpIXy2XIyJvEZFrROQ5EXleRK4VkR1nulxTgYgcKyI/EZHfFyszPCgiF4rIVkm6BSLy7yLyjIisE5FbRGTPmSr3VCIivyhWsjg/2d/oNhCRo0TkThF5sXjul4nIoeZ40+t/gIjcJCJPicgLIrJcRD6epKm1ksmmjvjKLZNK33RwIjKXENvsL4C/B/4OeCshNtkWM1m2KaJjIFMJ7v83FMdPAz5AsGRdIiI7zEShpwoROR7Yq2R/o9tARD5F8C+6G3gfYc7jamBucbzp9f9L4BZCnU4G3k8IB3W5iJxqkl5eHP8X4GjgCcJKJntPb4l7ZnzllrWElVty1K1vT8GK+x5V7YsXwSBlFNjN7NuF4Fj4uZku3xTU9w0l+04kxCg4tPh8TPH5EJPmdQT/kK/NdB0msS0WAE8SYtYpcL451tg2AHYGXgb+sSJNY+tf1OVLBD/aLZP9dwF3Fdt7FW1wkjk+RJivWTzTdZhgfQfM9ieKeu2cpKlVX8ISNq8Q5r9s/luBFTNd1+l49c0IjmA1s1RVX3MO12Am+kvCl7xRaL1Apu8B/ldVl5h8zxH+0TepTb5MMAu+suRYk9vg44Q40t+qSNPk+kMItruB0NFbnqOlQJWuZAJcBRwuIsPTUM5JQXtYuYX2+lYFK95TRHbpvcSbNv3Uwe0B3F+yfxUtx7+mkwYyrWqTHYvwOX2NiBxIGLl+OpOkyW1wIPA/wHEi8rCIbBSRh0TEtkWT6w+tWIVfE5HtRGS+iJwMvIsQrBfqrWTSJOrWt06w4kbTTx1c1WoFC6a5LNOOlAcy7bSCQ1+3i4jMJoT1+Yqqti1MUtDkNtiOMM98MWH9q8OAm4Gvi8hnijRNrj+qej/Bz/YY4I+Eun4DOEVVryqS1V3JpClM6cotTaJvVhPYnJEJBDJtGGcAcwiOpJsjA8BWwMdUdXy9rNsKq7ozReRrM1Ww6UJE3gr8hDDqOIUgVR4DfEtE1qvqFTNZPmfTpp86uLWU/xvN/ZtpBFIdyLSqTcaP9yWF+8dZhIn24WQeZVhCUNYXaHAbAH8ijOBuTvbfRLCafDPNrj8EI5MNhFBQG4p9t4rI64F/FZErqb+SSVOYtpVb+p1+kijHVwNP2Z0aQTf7EekcyLSqTR5X1RenuIhTya7ACGFC3AZlheBCsRbYk2a3waoOx8dodv0h3OP7TOc2zm+A1xMsBVcBuxSuRJZ0JZOmULe+PQcr7nf6qYNbDOwvIruO7yikmgOKY41C4kCm79XyQKaLge1F5CCTbx6wiP5vk3sJQVnTF4RO7xDCF7nJbXBd8X54sv8IQqy/J2l2/SG4h+xdzMda/oqwLOEa6q9k0hR85Za6zLSfQt0XsAXhB20lQYN/DyGY8yMkPjJNeBEcu5UQuHr/5LVDkWYA+BXwB+A4wg/h7YQv/Vtmug5T1C6pH1xj2wAQQnCDPxHmnw4Dvlu0wceaXv+ifscW9b2x+N4fRnBaVuASk+4qwqj+EwQLy2sIHeDCma5Dl3U+1vwGnFp8Pmii9SUYJ60HPkcw1rmMMPI/eqbrOS1tOdMFmOCN35Ew4fw8Yf7lehInyKa8gMeKh7vsdY5JtzXwveIH7SWCE+deM13+KWyXqINrehsA8whWg/9HkJ9WACdsLvUv6ndk0Wk/XXzv7wX+ARg0aeYAlxBGfOuBXwMHz3TZu6xv7nt/+0TrS1iC5mzg9wSXgRXAsTNdx+l6ebBlx3Ecp5H00xyc4ziO49TGOzjHcRynkXgH5ziO4zQS7+Acx3GcRuIdnOM4jtNIvINzHMdxGol3cI7jOE4j8Q7OcRzHaSTewTmO4ziN5P8BX/qlIkOFsR8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#NBVAL_SKIP\n",
    "\n",
    "# Plot FWI result\n",
    "from examples.seismic import plot_image\n",
    "\n",
    "slices = tuple(slice(param['nbl'],-param['nbl']) for _ in range(2))\n",
    "vp = 1.0/np.sqrt(result['x'].reshape(true_model.shape))\n",
    "plot_image(true_model.vp.data[slices], vmin=2.4, vmax=2.8, cmap=\"cividis\")\n",
    "plot_image(vp, vmin=2.4, vmax=2.8, cmap=\"cividis\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgoAAAF+CAYAAADjmjIfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd5hV1bnH8e87jaHDNNoMvUlHxgKoQEwsiVFRwG7UqLFETW40mhhjYotJbmyJ2KIxkUQEFY0aSzSCCFgGKQqCgLQBhKEjdcp7/zgH7ziZcmY4Z/aZM7/P85znzOy9197vfgjh59prrW3ujoiIiEhlkoIuQEREROKXgoKIiIhUSUFBREREqqSgICIiIlVSUBAREZEqKSiIiIhIleo9KJhZrpn90czmmNkeM3Mz61rhmHwze9TMloSPWWNmfzezbrW4zmXh9vvNbKmZXRHtexEREUl0KQFcsycwAZgLzAROqOSYs4H+wAPAIqATcAtQYGZD3H1tdRcws8uAR4DfAG8CxwMTzczc/aGaCszKyvKuXbtGfEMiIiIN2dy5cze7e3Zl+6y+F1wysyR3Lwv/fCnwGNDN3VeVOybb3YsqtOsCrATucPdfVnP+FGA98Kq7f6/c9ieAU4EO7l5cXY35+fleUFBQ63sTERFpiMxsrrvnV7av3h89HAwJNRxTVMm21UARod6F6gwHsoFJFbY/BWQCx0RWqYiIiDSYwYxmdhiQA3xaw6H9w9+fVNi+KPzdL5p1iYiIJLIGERTCjxMeJtSj8HgNh2eEv7dV2L61wn4RERGpQYMICsCfgBHA+e5eMQBEhZldbmYFZlZQVPRfTz5EREQapbgPCmZ2N3A5cIm7vxFBk4NBom2F7Qd7ErZSCXd/1N3z3T0/O7vSgZ8iIiKNTlwHBTO7GbgRuNbdn4qw2cGxCP0rbD84NmFxNGoTERFpDOI2KJjZtcAdwM3u/qdaNJ0DbAbOq7D9fEK9CbOiU6GIiEjiC2LBJcxsXPjHYeHvk82sCChy9xlmdjZwH/Aa8B8zO7pc853uvrjcuZYDq939eAB3LzazWwgtsLSO0IJL3wAuAa5x9wMxvTkREZEEEkhQAKZW+H1i+HsGMBo4CbDw90kVjj14zEEpQHL5A9z9YTNz4CfADcAa4IfuPhERERGJWL2vzNgQaGVGERFpTOJqZUYRERFpOBQUYmznvmL+vXhj0GWIiIjUiYJCjN3zxmdcOWkuC9ZuD7oUERGRWlNQiLEff7M3OS2bcN3keezeXxJ0OSIiIrWioBBjrZulcu9ZQ1izdQ+/+ueimhuIiIjEEQWFenBU90yuGt2TqXMLeXnh+qDLERERiZiCQj257pu9GJLXhp89/zHrtu8NuhwREZGIKCjUk9TkJO4/ewhlZc6PJ8+ntEzrV4iISPxTUKhHXTKbc9tpA/hg1VYemr486HJERERqpKBQz844vBPfHdyRe99cxrw122puICIiEiAFhXpmZtxx+gDat0rnusnz+VJTJkVEJI4pKASgddNU7jt7CIXb9nDri5oyKSIi8UtBISBHdM3gh2N68txHhfxzgaZMiohIfFJQCNC1x/diaOc23DztYwq37Qm6HBERkf+ioBCglOQk7j9rKO7w42fmU1JaFnRJIiIiX6OgELDOmc24/fT+fLhqGxOnrwi6HBERka9RUIgDY4fmctqQjtz/1jLmrtaUSRERiR8KCnHi9tMH0KF1Oj96Zh679hUHXY6IiAigoBA3WqWnct9ZQ1i3ba+mTIqISNxQUIgj+V0zuOYbvXh+3jpenL8u6HJEREQUFOLNNd/oybAubfnFtE9Yu1VTJkVEJFgKCnEmJTmJ+84aAsCPNGVSREQCpqAQh/IymnHH2AHMXb2NP72tt0yKiEhwFBTi1GlDOjF2aCceeGsZc1dvDbocERFppBQU4thtp/WnU9umXDd5Pjs1ZVJERAJQ70HBzHLN7I9mNsfM9piZm1nXSo67y8zeMLMt4WMuqsU1ngy3qfi5L4q3EnMt01O576yhbNixj1++8EnQ5YiISCMURI9CT2ACsA2YWc1x1wBNgZfreJ0iYHiFz711PFdghnVpy7Xf6MUL89fzwjxNmRQRkfqVEsA133H3dgBmdilwQhXHtXb3MjPrCVxYh+sccPf36lpkPLl6TA9mLiviFy98wrAubcnLaBZ0SSIi0kjUe4+Cu0c03y/S4xqDlOQk7j1rCAZcN3mepkyKiEi9SeTBjDlmttnMSszsMzO70cySgy6qrvIymnHnGQP5aM12HviPpkyKiEj9COLRQ32YD8wFFgHpwFjgN0Av4NLKGpjZ5cDlAJ07d66fKmvp1MEdmb50E3/6zzKO7ZXFEV0zgi5JREQSXEL2KLj7fe7+R3f/j7v/y90vA+4Hvm9mvapo86i757t7fnZ2dv0WXAu3nTaA3LbN+NHk+ezYqymTIiISWwkZFKrwdPg7P9AqDlGLJincf/YQvti5j1te+AR3D7okERFJYI0pKBzU4P9lHdq5LT86vhf/XLCeaZoyKSIiMdSYgsJ5hELCh0EXEg1XjenJkV0z+OWLi1izRW+ZFBGR2AgkKJjZODMbBwwLbzo5vG1UuWNGhY85Kbwpv1y78ud6y8yWl/u9i5m9Y2ZXmdkJZvZdM3uC0AJOj7j7itjeXf1ITjLuPXsIZnDt5HkUa8qkiIjEgAXxjNvMqrroDHcfHT5mOjCqsoPc3cqdazrQ1d27hn/PAJ4AhgLtgDJgSXjbxEjWZ8jPz/eCgoLIbiZgLy1YzzVPz+Oab/TkJyf0CbocERFpgMxsrrtXOoYvkOmR5f+hr+aY0RGea3SF37cCp9epsAbou4M7MuOzIh58eznH9srmyG6aMikiItHTmMYoJKxfndqfvIxm/PgZTZkUEZHoUlBIAKEpk0PZuHMfN0/7WFMmRUQkahQUEsSQvDb8+Fu9eXnhBp77SFMmRUQkOhQUEsgVo3pwVLcMbn3xE1Zt3h10OSIikgAUFBJIcpJx71lDSE4yrntmvqZMiojIIVNQSDAd2zTlN2cMYsHa7dz/5rKgyxERkQZOQSEBfWdQB8YPy+XB6ct57/MtQZcjIiINmIJCgvrVqf3pmtk8NGVyj6ZMiohI3SgoJKjmTVK476whFO3az881ZVJEROpIQSGBDc5rw/+c0JtXPt7A1LmFQZcjIiINkIJCgvvBcT04unsGv/rnIlZqyqSIiNSSgkKCOzhlMjU5iR/pLZMiIlJLCgqNQIfWTbn7jIEsKNzBvf/+LOhyRESkAVFQaCROHtiBs/LzeGjGCuas0JRJERGJjIJCI/LL7/ajW3jK5PY9B4IuR0REGgAFhUakefgtk1t27+dnz2vKpIiI1ExBoZEZmNuan5zQh1c/+YIpBWuDLkdEROKcgkIjdPmx3RnRI5Nf/XMxnxd9GXQ5IiISxxQUGqGkJOOeCUNokprEdZPnc6BEUyZFRKRyCgqNVPvW6dx9xiA+XreDezRlUkREqqCg0IidNKA95xyZxyPvrGD28s1BlyMiInFIQaGRu+WUfnTLas6Pp8xn225NmRQRka9TUGjkmqWl8MDZQ9m6+wA3Pb9QUyZFRORrFBSEAZ1ac8OJfXh90UYmf6gpkyIi8v8UFASAS4/pzjE9s7jtpcWs0JRJEREJq/egYGa5ZvZHM5tjZnvMzM2sayXH3WVmb5jZlvAxF9XyOqeb2Twz22dmq83sF2aWHKXbSDhJScYfJgwmPTWJ6ybP05RJEREBgulR6AlMALYBM6s57hqgKfBybS9gZicCzwEfAicD9wO/AO6q7bkak3at0vntmYP4ZN1O/vDG0qDLERGROJASwDXfcfd2AGZ2KXBCFce1dvcyM+sJXFjLa9wNvOvul4d/f9vMWgC/MLN73f2LOlXeCJzQvz3nHtWZR975nON6ZzOyZ1bQJYmISIDqvUfB3SPq0470uIrMLA8YAkyqsOspIJVQD4NU45bv9KNHdnP+Z8p8tmrKpIhIo5aIgxn7h78/Kb/R3VcCe4B+9V5RA9M0LZn7w1Mmb3xOUyZFRBqzRAwKGeHvbZXs21Zu/9eY2eVmVmBmBUVFRTErrqEY0Kk1N57Ul38v3sg/PlgTdDkiIhKQRAwKdeLuj7p7vrvnZ2dnB11OXLhkZDeO7ZXF7S8vZvmmXUGXIyIiAUjEoHCwJ6FtJfvaAlvrsZYGLSnJ+MP4wTRLS+Hap+ezv6Q06JJERKSeJWJQWBT+7l9+Y3ithmbA4nqup0HLCU+ZXLxhJ//7uqZMiog0NgkXFNx9DbAAOK/CrvOBYuDVei+qgftWv3acf3RnHpu5kpnLNH5DRKQxCWIdBcxsXPjHYeHvk82sCChy9xnhY0YB2UD78DH5ZvYlgLs/W+5cbwFd3L1nuUv8HHjZzB4BngaGElpw6X6toVA3N3+7H+99vpWfTFnAaz86jozmaUGXJCIi9cCCmPpmZlVddIa7jw4fMx0YVdlB7m7lzjUd6OruXStc4wzgVqAvsBH4M3Cnu9f4oD0/P98LCgpqOqzRWbx+J6c/OIvjemfz2IXDMLOaG4mISNwzs7nunl/ZvkB6FMr/Q1/NMaMjPFelx7n788DztSpMqtWvYytuPLkvt7+8mEnvr+GCo7sEXZKIiMRYwo1RkNi6eERXjuudzR0vL2bZRk2ZFBFJdAoKUitJScb/jh9EiyYpXDtZUyZFRBKdgoLUWk7LdH43bhCfbtjJ717TlEkRkUSmoCB1cvxh7bhweBcef3cl73ymKZMiIolKQUHq7OffPoze7Vrwk6kL2PLl/qDLERGRGFBQkDpLTw29ZXLH3mJ++qzeMikikogUFOSQHNahFTed1Je3lmxi0nurgy5HRESiTEFBDtnFI7syuk82d7zyKZ9pyqSISEJRUJBDZmb8ftxgWqancO3T89hXrCmTIiKJQkFBoiK7ZRN+P24wS77YxW9fWxJ0OSIiEiUKChI1Y/rmcNGIrvxl1iqmL90UdDkiIhIFCgoSVTed3Jc+7Vpy/dSFbNaUSRGRBk9BQaIqPTWZB84Zys59xdwwdYGmTIqINHAKChJ1fdq35Ocn9+XtpUX8bY6mTIqINGQKChIT3xvRlTF9srnzX5+y5IudQZcjIiJ1pKAgMWFm/H78YFqlp3Ld0/M1ZVJEpIGKKCiYWWszaxLrYiSxZLVowv+OH8TSjbu4+1VNmRQRaYhqDApmlgJsAU6IfTmSaEb3yeHikV15cvYq3l6iKZMiIg1NjUHB3UuAjYD6jqVObjypL33bt+SGZxdQtEtTJkVEGpJIxyhMAi6NZSGSuA5Omdy1r4QbntWUSRGRhiQlwuNWAeea2YfAi8AG4Gv/b+/uT0S3NEkkvdu15ObvHMYvX1zEk7NXcfHIbkGXJCIiEYg0KDwY/u4EDKtkvwMKClKtC47uwoylRfzm1SUc3T2Twzq0CrokERGpQaSPHrrV8Okek+okoZgZvxs3iNZNU/WWSRGRBiKioODuq2v6xLpQSQyZLZrwh/GDWbbpS+7616dBlyMiIjWI9NEDAGY2ABgFZABbgenuvigWhUniOq53Nt8/phuPv7uSUb2zOf6wdkGXJCIiVYh0waUUM5sELAD+CPw6/L3QzJ4ys+RIL2hmuWb2RzObY2Z7zMzNrGslx6Wb2e/NbIOZ7Q0ff1yE15gePm/Fz48irVNi66cn9eGwDq244dmFbNq1L+hyRESkCpGOUbgVmAD8ktCYhKbh718CZ4W/I9UzfK5twMxqjnscuCx87lMIzbR43cyGRHidhcDwCp/JtahTYqhJSjIPnD2E3ftLuH7qQsrKNGVSRCQeRRoUzgfucPc7w2MS9oe/7wTuAC6sxTXfcfd27v5tYGplB5jZYOBc4Mfu/pi7v0UoXKwBbovwOrvc/b0Kny9qUafEWK92LfnFKf1457Mi/jJ7VdDliIhIJSINCh2B2VXsmx3eHxF3L4vgsFOBYuCZcu1KCPUInKj3TiSO84/qzDcPa8dvX13C4vV6y6SISLyJNCisB0ZWsW9EeH809QdWuvueCtsXAWmEHl/UZKiZ7TCzYjNbaGbfj3KNEgVmxm/PHEjrZqlcN3keew9oyqSISDyJNCj8HbjZzG4xs+5m1tTMupnZz4CbgaeiXFcGoTEMFW0tt7867wA/ItQzMQ5YBvzZzH5RVQMzu9zMCsysoKioqA4lS11ltmjCPRNCUybv/NfioMsREZFyIp0e+StCiyr9OvzzQQY8TeTjBuqFu1ccXPmimU0jFHbuc/cvK2nzKPAoQH5+vkbW1bNje2Vz2bHdeGzmSkb1zuFb/TRlUkQkHkS64FKJu58LDAR+SGgmwg+Bge5+Xnj8QDRtA9pWsv1gT8LWSvbV5GkgndA9SBy6/sQ+9OvQihufW8imnZoyKSISD2oMCmaWZmZbzexUd1/k7g+FZz88FMPFlhYB3cysWYXt/YADwPJDOLd6C+JUk5TQWyb3HCjhJ1MXaMqkiEgcqDEouPsBoASoz//EewlIBcYf3GBmKYTWbHjD3ffX4ZznAXuBj6NSocREz5wW3HJKP2Yu28wTs1YGXY6ISKMX6RiFFwgNCnwjGhc1s3HhHw++ifJkMysCitx9hrvPM7NngPvMLBVYCVxJaJGn8yqcazmw2t2PD/9+LHAT8Dyh12O3Br5HaGDjTe6+Oxr3ILFz7pGdmb60iN+9tpThPTLp37F10CWJiDRa5l5z966ZjQUeAN4nFBo2UKEL393/E/FFzaq66Ax3Hx0+pilwJ6GFl9oQWj76RnefXuFcq4BV5dr1JLS89CAgi9B6DAuBP7r705HUl5+f7wUFBZHejsTA1t0HOOm+d2iZnsLL1xxL07SIVwkXEZFaMrO57p5f6b4Ig0JViyQ5oZkP7u4J8//kCgrxYdbyzZz/+Pucc2Rn7hqrMagiIrFSXVCI9NHDmCjWIxKRkT2zuPzY7jzyzueM6p3Nif3bB12SiEijU2NQMLM0YAjwlrt/EvuSRP7fT07ow6wVm7npuYUMyWtDu1bpQZckItKoRDrr4W5qXg1RJOrSUpK4/+yh7Csu43+mzNeUSRGRehbpEs6fElqZUaTe9chuwS+/249Zy7fw53c/D7ocEZFGJdKg8EvgFjPTiDIJxNlH5HFi/3b8/vWlfLJuR9DliIg0GpEGhRuBFsA8M1tuZjPN7J1ynxkxrFEEM+PuMwaR2bwJ106ex54D0V41XEREKhNpUCgFFgMzgbWEVmosLfepavqkSNS0bZ7GPRMGs3Lzbm5/+dOgyxERaRQimh55cDEjkaCN6JnFD47rwcMzVjCqdzYnDdCUSRGRWIq0R0EkbvzPt3ozsFNrbnp+IV/s0FsmRURiKeKgYGadzOweMysws5VmNiC8/UdmdlTsShT5utCUySHs15RJEZGYiygomFl/Qm9dvABYD3QG0sK7uwDXxaQ6kSp0z27Br07tx+wVW3h0pqZMiojESqQ9Cn8gtJZCN+AMQu93OGg2cHSU6xKp0YT8PE4e0J7/fX0pHxdqyqSISCxEGhSOAe529y+p8NZIYCOgEWVS78yM35wxkOyWTbhOUyZFRGIi0qBQ3fTHLGBvFGoRqbU2zdK4Z8IQVm7ZzW0vLQ66HBGRhBNpUPgAuLiKfROAWdEpR6T2hvfI5IpRPZj84Vpe/XhD0OWIiCSUSIPC7cB3zewNQgMaHfimmf0VGAvcGaP6RCLyP9/qzeDc1tz0/Mes364OLhGRaIkoKLj7DOB0QoMZnyA0mPFu4FjgdHd/P2YVikQgNTn0lsni0tCUyVJNmRQRiYqI11Fw91fcvRfQm9DgxsPcvbu7vxqz6kRqoWtWc351an/e+3wrj7yzIuhyREQSQq1XZnT35e4+292XxqIgkUMxflgu3xnYgXve+IwFa7cHXY6ISIOnJZwloZgZd40dSE7LJvzomfns3q8pkyIih0JBQRJO62ap3HPWEFZt2c2vX1oUdDkiIg2agoIkpKO7Z3LV6B5MKSjklYWaMikiUlcKCpKwfvTN3gzOa8PPnl/IOk2ZFBGpEwUFSVipyUk8cPYQSsucHz+jKZMiInWRUtUOMzuuNidy93cOvRyR6OqS2ZxfnzaA66cu4OEZK7h6TM+gSxIRaVCqDArAdP77BVCVsfBxyZFc0MxygRuBfGAw0BTo5u6rKhyXTmhFyPOBNsB84MZIA4mZXQb8hNAiUauAe9394UjaSmI58/BOTF+6iXv//Rkje2YxJK9N0CWJiDQY1QWFMTG6Zk9C74eYC8wETqjiuMeB7wA3AJ8DVwOvm9lwd59f3QXCIeER4DfAm8DxwEQzM3d/KCp3IQ2GmXHn2IHMW7Od6ybP45Vrj6VFk+r+py8iIgeZe/0+tzWzJHcvC/98KfAYFXoUzGwwoR6ES9z9L+FtKcAiYKm7n1rN+VOA9cCr7v69ctufAE4FOrh7cXU15ufne0FBQR3vUOLVByu3cvajczjj8Fz+d/zgoMsREYkbZjbX3fMr21erwYxmlmVmp5jZ98wsI7wt3cxqsxR0da+sPuhUoBh4ply7EmAycKKZNamm7XAgG5hUYftTQCah5aelETqyWwZXj+nJs3MLeXnh+qDLERFpECL6B95Cfg8UAv8k9GKoruHdLwI3R7mu/sBKd99TYfsiII3Q44vq2gJ8UklbgH6HXp40VNce34sheW342fMfs3Zrxf95iYhIRZH2BPwM+CFwG3AUoQGMB70EnBLlujKAbZVs31puf3VtqaR9tW3N7HIzKzCzgqKioogLlYYlNGVyKABnPjSbjwt3BFyRiEh8izQoXArc5u53AR9V2Lcc6BHVqgLg7o+6e76752dnZwddjsRQ58xmPHvFCFKTk5jwyBzeWPRF0CWJiMStSINCJ+C9KvYdAJpHp5yvbAPaVrL9YG/A1kr2lW9LJe0jaSuNRJ/2LZl29Qh6t2/JDybN5c8zP6e+B/aKiDQEkQaFdcCAKvYNBlZGp5yvLAK6mVmzCtv7EQomy2toC/8/VqF8W4DFh16eJIKclulMvuxoTurfnjte+ZRfvPAJJaWRjLUVEWk8Ig0KU4FfmtnIctvczHoTWtRocpTreglIBcYf3BCe9ngW8Ia776+m7RxgM3Behe3nE+pNmBXdUqUha5qWzIPnHs4Vo3rw9/fXcMlfC9i1r9rZsyIijUqkq878ChgBvAOsDm+bCuQBs4G7a3NRMxsX/nFY+PtkMysCitx9hrvPM7NngPvMLJVQj8WVhFZZPK/CuZYDq939eAB3LzazWwgtsLSO0IJL3wAuAa5x9wO1qVUSX1KScdPJfema2YxfvPAJ4x6aw+MX5ZPbtmKHlohI4xPxgktmlgycC5wI5ABbgNeAv4fXOIj8omZVXXSGu48OH9MUuDN8zTbAAkJLOE+vcK5VwKqD7cpt/wGh3o4uwBpCSzhPjKQ+LbjUeM1avpkrJs2lSUoyf/5evpZ7FpFGoboFl+p9ZcaGQEGhcVu+aRcXP/khRbv2c++EIZw8sEPQJYmIxFTUVmYUaQx65rRk2lUj6dehFVf+/SMenrFCMyJEpNGq7jXTK4ns7ZEAuHv3qFQkEgeyWjThH5cdzfVTF3D3q0tYtXk3t58+gNRkZWsRaVyqG8w4g68HheOBdoRmDWwM/zwS+AJ4K1YFigQlPTWZB84eStfM5vzp7eWs3baHiecNo3XT1KBLExGpN1UGBXe/6ODPZnY5oaWbR7h7YbnteYQGNM6JYY0igUlKMq4/sQ9ds5rzs+cXcuZDs/nLRUeQl6EZESLSOETaj3oDcGv5kADg7muBXwM3RrswkXgyblguf7vkKIp27ef0B2cxd3VlryIREUk8kQaFXGBfFfv2E1riWSShDe+RyfNXjaBFegrnPPYeLy3Qq6pFJPFFGhQWAzeYWXr5jeG1Dm5AyyJLI9EjuwXTrhrJ4NzWXPP0PB58e7lmRIhIQot0ZcafAq8Aa8zsX/z/YMZvA62Bk2NTnkj8yWiexqRLj+LGZxfy+9eXsnLzbu4aO5C0FM2IEJHEE1FQcPe3zGwo8AvgWKADsAF4A7jD3ZfErkSR+NMkJZl7zxpC16zm3PfmMgq37eHh84fRplla0KWJiESVVmashFZmlNp4Yd46fvrsQnLbNuWJi46ga1a037ouIhJbUVuZ0UL6m9mxZtbPzCw6JYo0XKcP7cSkS49i254DjJ04iw9XbQ26JBGRqIk4KJjZpYQeNywEpgMfA+vN7PuxKU2k4TiyWwbTrhpJ22ZpnPfY+7w4f13QJYmIREVEQcHMzgMeJRQOLiE0iPGS8O+Pmtk5MatQpIHomtWc568awdDObbhu8nzuf3OZZkSISIMX0RgFM1sALHT3CyrZ9xQw0N2HxKC+QGiMghyKAyVl3PT8Qp7/aB1jh3bi7jMH0iQlOeiyRESqFI0xCn2ASVXsmxTeLyJAWkoSfxg/mOtP6M20eeu44M8fsG33gaDLEhGpk0iDwi5CqzNWJje8X0TCzIwffqMXD5wzlPmF2xk7cRafF30ZdFkiIrUWaVB4FbjLzI4tv9HMhgN3hPeLSAWnDu7I05cdxc59JZzx0Gze+3xL0CWJiNRKpEHhp8AOYLqZrTGz981sNfAusDO8X0QqMaxLBi9cNZLM5mlc8Pj7PDe3sOZGIiJxIqKg4O5fAEOA6wi9Unon8B5wDTDU3TfGrEKRBNA5sxnPXzmSI7pm8JOpC7jnjaWaESEiDYJWZqyEZj1IrBwoKeMXL3zMlIJCvju4I78fN4j0VM2IEJFgVTfrIdKXQlU84X/1RLh7WV3OJdKYpKUk8dszB9EtqwW/fW0J67fv5dELhpHZoknQpYmIVCrSBZeamtndZrbCzPYDxRU+mvslEiEz48rRPZh43uF8sm4HYyfOZvkmzYgQkfgUaY/CROA84CVgMgoGIofs2wM70KF1Opf9rYAzJs7i4fOHMaJnVtBliYh8TaQrM24Bfu3uD8S+pOBpjILUp7Vb93DJkx+ycvNu7ho7kAlH5AVdkog0MtFYmXE/8Gn0ShKRg/IymvHcVSMY3iOTnz63kN++toSyMg0yFpH4EGlQeBI4O4Z1VMrMxpjZu2a218y2mtlTZtYuwrZexSdh3kkhiaNVeipPXHQE5x7VmcAC90QAACAASURBVIemr+Cap+exr7g06LJERCIeo3AL8JCZvQG8DmyreIC7PxHNwsKrQB683plAJqFVIN8ys2Huvj+C0zwJPFJh22fRrFMkWlKTk7jz9AF0y2zOXa9+yrrte3nswnyyW2pGhIgEJ9KgMAw4FcgBvlnJfgeiGhSAW4HVwOnuXgJgZp8CHwLfJzTAsibr3P29KNclEjNmxmXHdadzZjOumzyP0x+cxV8uPoLe7VoGXZqINFKRPnp4GNgCnEboTZHdKny6x6C2o4F/HwwJAO5eEK5jbAyuJxI3Tuzfnik/GM6B0jLOnDibmcuKgi5JRBqpSINCX+AGd3/J3Ze5++qKnxjUVkrl0zD3AwMiPMeVZrbfzPaY2X8qvtRKJJ4Nym3DC1ePpFPbplz0lw/5x/trgi5JRBqhSIPCUqB5LAup4ppHl99gZl2ADkBGBO0nAVcRelRyOaExDv8xs9GVHWxml5tZgZkVFBXpv94kPnRq05SpVwznmJ5Z/Hzax9z1r081I0JE6lWk6yicCPwOODVGvQeVXfM8Qv/Y3wk8QCgcPAqMAIrdvWktz9cS+ARY6+7HVHes1lGQeFNSWsavX1rMU++t5sT+7bjvrKE0TdM7IkQkOqpbRyHSoDAT6EnoH+vP+O9ZD+7uow610EqueztwPZBOaMDkM4R6Nga4e63HRZjZROD77l7tMHIFBYlH7s5fZq3i9lcWM7BTa/58YT45rdKDLktEEkA0FlwqBZYAs4HN4d/Lf2LyQih3vwXIAgYBHdz9HKAX8O6hnDYatYnUNzPjkmO68dgF+Szf9CWnPziLJV/sDLosEUlwDeo102Z2EvAqMNLdZ9eybStCjx5Wuftx1R2rHgWJd5+s28H3//ohu/eX8qdzhzK6T07QJYlIAxaNHoV6Z2ZDzeznZnZS+HMH8CLwu/Ihwcy6mFmJmf2y3LbrzewxMzvXzEab2feAWUB74OZ6vxmRKBvQqTUvXD2SzhnNuOTJD3lqzqqgSxKRBFXroGBmyWZWamaHx6Kgcg4A3yb0tspp4Z+vcPcbK5YEJPP1e1kK9CM0CPLfwD3ASuAYd58Z47pF6kWH1qEZEWP65HDLi4u47aXFlGpGhIhEWaQrM1ZkUa2iEu6+CKh2dkL4uFUV63H3lwi9ElskoTVvksKjF+ZzxyuLeWLWStZs3cP9Zw+heZO6/tUWEfm6uj560H+2iMSJ5CTj1u/257bT+vOfJRuZ8MgcvtixL+iyRCRB1DUoxLxHQURq58LhXXn8e0ewavNuTn9wFovW7wi6JBFJALUOCu5eCowhNA5AROLImL45TL1iBGYw/uE5vPXpxqBLEpEGrk49Cu4+w913R7sYETl0/Tq24sWrR9IjuwWX/a2Av8xaGXRJItKA1RgUzKyrmY00s+wq9meZ2YXRL01E6iqnVTrP/OBovnlYO3790mJuffETSkpjsi6aiCS4KoOCmaWa2TPACuAdYL2ZPWxmFV8O1QP4SwxrFJE6aJaWwkPnD+OyY7vx1zmruexvBXy5v6TmhiIi5VTXo3AdcCrwK+A7wB+AC4HZZtY+9qWJyKFKTjJu/k4/7jh9AO8s28y4h2azfvveoMsSkQakuqDwPeAOd7/d3V9z95sIvfa5FTDLzLrVS4UicsjOP7oLT1x0BIXb9nL6g7P4uFAzIkQkMtUFhW5UePmSuy8k9JrnvcC7ZtYvhrWJSBSN6p3Nc1eOIDU5iQmPzOGNRV8EXZKINADVBYXtQGbFje6+ATgO2ADMACp9iYSIxJ8+7Vsy7eoR9G7Xgh9MmsufZ35OQ3oxnIjUv+qCwgJC71f4L+6+ldBaCkuA+2NQl4jESE7LdCZfPpyT+rfnjlc+5eYXNCNCRKpWXVB4ATjBzP6rVwHA3XcBJxB66ZJWahRpQJqmJfPguYdzxage/OP9NVzy1wJ27SsOuiwRiUNVBgV3f8zdO7v7lmqO2evuJ7t73L6uWkQql5Rk3HRyX+4+YyCzl29m3ENzKNy2J+iyRCTOHPI/8Gb2TTNbEI1iRKT+nX1kZ/56yZGs37GX0x+czfy124MuSUTiSDR6AloDA6JwHhEJyMieWUy7agRN05I4+9E5vPrxhqBLEpE4oUcGIgJAz5yWTLtqJP06tOLKv3/EwzNWaEaEiCgoiMj/y2rRhH9cdjSnDOrA3a8u4WfPf0yxZkSINGopQRcgIvElPTWZB84eStfM5vzp7eWs3baHiecNo3XT1KBLE5EAVBkUzKx7hOfQex9EEkxSknH9iX3omtWcnz2/kDMfms1fLjqCvIxmQZcmIvWsuh6F5UAkDygtwuNEpIEZNyyXTm2acsWkuZz+4CwevTCfYV3aBl2WiNSj6oLCxfVWhYjEreE9Mnn+qhFc8uSHnPPYe/xh/GC+O7hj0GWJSD2pMii4+1/rsxARiV89slsw7aqR/OCpAq55eh5rtu7hqtE9MNOirCKJTrMeRCQiGc3TmHTpUZw2pCO/f30pNzy7kAMlmhEhkug060FEItYkJZn7zhpC18zm3P/WMgq37eHh84fRplla0KWJSIyoR0FEasXM+PG3enPvWYP5aPV2zpg4m1WbdwddlojESFwHBTMbY2bvmtleM9tqZk+ZWbsI26ab2e/NbEO4/RwzOy7WNYs0FmOH5jLp0qPYtucAYyfO4sNVW4MuSURiIG6DgpkdC7wBbAfOBK4DjgPeMrMmEZziceAy4JfAKcAG4HUzGxKbikUanyO7ZTDtqpG0aZbGeY+9z4vz1wVdkohEWdwGBeBWYDVwurv/y92fIhQY+gPfr66hmQ0GzgV+HH5d9lvABGANcFtsyxZpXLpmNWfaVSMY2rkN102ez31vfqZ3RIgkkHgOCkcD/3b3koMb3L0A2AKMraHtqUAx8Ey5tiXAZODECHskRCRCbZql8dT3j+KMwztx35vL+J8pC9hfUhp0WSISBfEcFEqBA5Vs30/Nr7XuD6x09z0Vti8C0oCeh16eiJSXlpLEH8YP5voTejNt3jou+PMHbNtd2V9hEWlI4jkoLCXUq/AVM+sCdAAyamibAWyrZPvWcvu/xswuN7MCMysoKiqqQ7kiYmb88Bu9eOCcocwv3M7YibP4vOjLoMsSkUMQz0HhfuBIM7vDzHLMrC/wFFAW/kSVuz/q7vnunp+dnR3t04s0KqcO7sjTlx3Fzn0lnPHQbN77fEvQJYlIHcVtUHD3vwN3AD8BNgKLgXXAvwjNYKjONqCyN9cc7EnQPC6RGBvWJYMXrhpJZvM0Lnj8fZ6bWxh0SSJSB3EbFADc/RYgCxgEdHD3c4BewLs1NF0EdDOziu/E7Udo3MPyaNcqIv+tc2Yznr9yJEd0zeAnUxdwzxtLNSNCpIGJ66AA4O673f1jd99oZicBfYGHa2j2EpAKjD+4wcxSgLOAN9x9f8wKFpGvad0slScvPpIJ+bk88J/lXDt5PvuKNSNCpKGI23c9mNlQ4GTgo/CmY4AbgN+5++xyx3UBVgC3ufttAO4+z8yeAe4zs1RgJXAl0A04r/7uQkQgNCPit2cOoltWC3772hLWb9/LoxcMI7OFZiqLxLt47lE4AHyb0NoH08I/X+HuN1Y4zoBk/vteLgb+QmicwytAHnCSu3+EiNQ7M+PK0T2YeN7hfLJuB2Mnzmb5Js2IEIl3pueF/y0/P98LCgqCLkMkYc1bs43L/lbAgZIyHj5/GCN6ZgVdkkijZmZz3T2/sn3x3KMgIglqaOe2TLtqJO1apXPhEx8w5cO1QZckIlVQUBCRQORlNOO5q0YwvEcmP31uIb99bQllZerhFIk3CgoiEphW6ak8cdERnHNkZx6avoKLn/yQmcuKKFVgEIkbcTvrQUQah9TkJO4aO4BeOS24/61lXPD4B3Rsnc64YbmMG5ZH58yKy6GISH3SYMZKaDCjSDD2FZfy5qcbmVJQyMxlRbjD8O6ZjM/P5eQBHWialhx0iSIJqbrBjAoKlVBQEAne+u17ef6jQqYUFLJm6x5aNknhlMEdmZCfy5C8NphZ0CWKJAwFhVpSUBCJH2VlzgertjKlYC2vfvwFe4tL6ZXTgvH5uYwdmkt2Sy3aJHKoFBRqSUFBJD7t2lfMKws3MKVgLR+t2U5KkjGmbw4T8vMY3Seb1GSNzxapCwWFWlJQEIl/yzftYmpBIc99tI7NX+4nq0UTzji8ExPyc+mZ0zLo8kQaFAWFWlJQEGk4ikvLmLG0iCkFa/nPkk2UlDlDO7dhQn4epwzqQMv01KBLFIl7Cgq1pKAg0jAV7drPC/PWMaVgLcs2fUl6ahLfHtCB8fl5HNUtg6QkDYAUqYyCQi0pKIg0bO7OgsIdTClYy0vz17NrfwmdM5oxflguZw7LpWObpkGXKBJXFBRqSUFBJHHsPVDKa4s2MLWgkNkrtmAGx/TMYkJ+Ht/q1470VK3NIKKgUEsKCiKJae3WPUydW8hzcwtZt30vrZumctqQjkzIz6N/x1Zam0EaLQWFWlJQEElsZWXO7BVbmFKwltcWfcGBkjIO69CKCfm5nD6kE22bpwVdoki9UlCoJQUFkcZjx55i/rlgHVMKCvl43Q7SkpP4Zr8cxufncVyvbJI1AFIaAQWFWlJQEGmcPt2wk6kFhUybV8i2PcW0b5XOmcM6MX5YHl2zmgddnkjMKCjUkoKCSON2oKSMtz7dyJSCtcz4rIgyhyO7ZjA+P5dvD+xA8yZ68a4kFgWFWlJQEJGDNu7cx3MfFTK1oJCVm3fTPC2ZUwZ1ZHx+LsO6tNUASEkICgq1pKAgIhW5OwWrtzG1YC0vL9zAngOldM9qzvj8PM48vBM5rdKDLlGkzhQUaklBQUSqs3t/Ca98vIGpBWv5cNU2kpOM0b2zGZ+fyzf6tiMtRS+nkoZFQaGWFBREJFKfF33Js3MLeXZuIZt27SejeRpjh3ZiQn4efdrr5VTSMCgo1JKCgojUVklpGTOXbWZKwVre/HQjxaXO4NzWjMvP49TBHWndVC+nkviloFBLCgoicii27j7w1cuplnyxiyYpSZw0oD0T8vMY3j1TL6eSuNNgg4KZjQRuBYYATYFlwJ/c/Yka2q0CulSya6y7v1DTdRUURCQa3J1P1u1kSsFaXpy/jp37SujUpinjhuUyblgueRnNgi5RBGigQcHMBgHvA+8B9wF7gHHA5cBV7v5QNW1XAUuAX1XYtdTdt9V0bQUFEYm2fcWlvLF4I1ML1vLu8s24w8iemUzIz+PE/u31cioJVEMNCncB1wMZ7v5lue1zANx9eDVtVwHvuvv5dbm2goKIxFLhtj08N3cdU+eupXDbXlqmp3Dq4NDLqQblttbaDFLvqgsK8by8WBpQDOytsH0H0Lb+yxERiY7cts247pu9uOYbPXlv5RamFoRmTfz9/TX0adeS8fm5jB3aicwWTYIuVSSuexQGEHr08CRwJ6FHD+OBicAF7j65mrarCIWJVCAZmAfcHcn4BFCPgojUv537inlpwXqmFhQyf+12UpKM4w/LYUJ+HqN6Z5OSrLUZJHYa5KMHADM7ApgGdApvKgaudPfHa2j3R+BDYCXQDvghMIpQwJhU03UVFEQkSJ9t3MXUgrU8/9E6tuw+QHbLJpxxeOjlVD1zWgRdniSgBhkUzKwX8BbwKfBHQo8gTgOuBC5y97/X4lzJhAZFtnf3vCqOuZzQQEk6d+48bPXq1Yd2AyIih6i4tIy3l2xiSkEhby/dRGmZM6xLWybk5/KdQR1poZdTSZQ01KAwFTgc6OvuxeW2/x04Echx97JanO+nwG+Bju6+obpj1aMgIvFm0659TPsotDbDiqLdNE1N5tsDOzAhP5cju2VoAKQckoY6mHEgsKB8SAj7ADgXyAG+qMN54zMZiYhUI6dlOj8Y1YPLj+vOR2u28+zctby0YAPPfVRI18xmjM/P44zDO9GhddOgS5UEE889CtOBzoR6FA6U2/4PYCzQuvz2Gs6VQmhgZJa7V7YQ09eoR0FEGoI9B0p49eMvmFKwlvdXbiXJ4Nhe2UzIz+Ob/XJokqK1GSQyDbVH4U/AVOAlM5tIaIzCqcA5wL0HQ4KZlQB/dffvh38/h9BYhn8BawkNZrya0GOMc+r7JkREYqVZWgpnDsvlzGG5rN6y+6uXU139j49o0yyV04d0Ynx+Lv07tg66VGnA4rZHAcDMTgZuBPoD6cAK4FHgEXcvDR/jhILCReHfjwbuCrfJAHYDBcDv3f31SK6rHgURaahKy5x3l29masFa3li0kQOlZfTv2IoJ+XmcNqQjbZqlBV2ixKEGOZgxSAoKIpIItu85wIvz1zOlYC2L1u8kLTmJE/q3Y3x+Hsf0zCJZL6eSMAWFWlJQEJFE88m6HTw7t5AX5q9j+55iOrRO/+rlVF0ymwddngRMQaGWFBREJFHtLynlzcWbmFKwlneWFeEOR3XLYEJ+HicPbE+ztHgeuiaxoqBQSwoKItIYrN++l+c/KmTq3EJWb9lDiyYpfHdwB8bn5zE0r43WZmhEFBRqSUFBRBoTd+eDlVuZUlDIvz7ewN7iUnpkN2dCfh5jD+9ETsv0oEuUGFNQqCUFBRFprL7cX8IrC9czpaCQuau3kZxkDOvSltF9shnTJ4e+7VuqpyEBKSjUkoKCiAgs3/QlL8xbx3+WbGLxhp0AtG+Vzpi+2YzqncMxvbL0vokEoaBQSwoKIiJf98WOfcz4bBPTlxYxc9lmvtxfQmqycUTXjK96G3rmtFBvQwOloFBLCgoiIlUrLi2jYNU2pn+2ielLili6cRcAndo0/So0jOiZqRkUDYiCQi0pKIiIRG7d9r3MWFrE20s3MWv5ZvYcKCUtOYmjumcwuk8OY/pk0y2ruXob4piCQi0pKIiI1M3+klI+XLmN6Us38fbSTawo2g1Al8xmjO6dzei+OQzvnkl6ql5YFU8UFGpJQUFEJDrWbt0TDg1FzF6xmX3FZTRJSWJ4j0zG9MlhTJ8cOmc2C7rMRk9BoZYUFEREom9fcSnvr9zK20s2MX3pJlZt2QNA96zmoUcUfbM5sluGXo8dAAWFWlJQEBGJvZWbd3/V2/De51s4UFJG09RkRvbMZHSfHEb3ySa3rXob6kN1QUFDUkVEJBDdsprTLasbF4/sxt4Dpcz5fDNvLwkNinzz000A9MppwZi+odCQ3yWDtJSkgKtufNSjUAn1KIiIBMfdWVG0+6sBkR+s3EpxqdM8LZljemUxpk8Oo/vk0L61lpaOFj16qCUFBRGR+PHl/hJmL9/M20uLmLF0E+t37AOgb/uWod6G3tkc3qUtqcnqbagrBYVaUlAQEYlP7s5nG7/k7aWhAZEFq7ZRUua0TE/huF7ZjOqTzeje2eS0Um9DbSgo1JKCgohIw7BzXzGzlm1menjBp0279gMwoFOr8COKbIbktSU5SYs9VUdBoZYUFEREGh53Z/GGnUxfWsT0pZuYu3obZQ5tmqVyXK9sRvfJZlTvbDJbNAm61LijoFBLCgoiIg3fjj3FzFxexNtLipjx2SY2f3kAMxiU24bRvbMZ0zeHQZ1ak6TeBgWF2lJQEBFJLGVlzqL1O3k7PJNi/trtuENm8zRG9Q6NbTiuVzZtm6cFXWogFBRqSUFBRCSxbd19gJnLinh7ySZmfFbEtj3FJBkM7dz2q96Gfh1aNZreBgWFWlJQEBFpPErLnIWF23k7PLZhYeEOALJbNmFU79Brs4/plUXrpqkBVxo7Cgq1pKAgItJ4Fe3azzufhWZRzFy2mR17i0lOMoZ1bsvovqHg0Ld9y4R6bbaCQi0pKIiICEBJaRnz124Pr9tQxKL1OwFo3yqd0X2yGd0nh5E9M2mZ3rB7GxpsUDCzkcCtwBCgKbAM+JO7P1FDuyTgRuAHQHtgKXCbuz8XyXUVFEREpDIbd+5jxtIipn+2iZmfbWbX/hJSkowjumYwJtzb0DOnRYPrbWiQQcHMBgHvA+8B9wF7gHHA5cBV7v5QNW3vBK4HbgbmAmcDlwGnuPu/arq2goKIiNSkuLSMuau3fbVuw5IvdgHQqU1TRvcJhYYRPTNplhb/719sqEHhLkL/2Ge4+5flts8BcPfhVbTLAdYCd7v7reW2vwVku/ugmq6toCAiIrW1fvver0LDu8s3s+dAKWnJSRzVPeOr12Z3z2oel70NDfU102lAMbC3wvYdQNtq2p0YbjupwvZJwBNm1s3dV0atShEREaBjm6ace1Rnzj2qM/tLSilYtS38Bswibn95Mbe/DJ0zmjGmTzaj++YwvHsm6anJQZddo3gOCk8CVwIPhB8l7AHGA8cDF1TTrj+wH1heYfui8Hc/QEFBRERipklKMiN7ZjGyZxY3fwfWbt3D9M+KmL5kE88UrOWvc1bTJCWJ4T0yv1q3oUtm86DLrlTcPnoAMLMjgGlAp/CmYuBKd3+8mjaPAqe6e/sK23sSGgx5obs/Vd119ehBRERiZV9xKR+s3PrVTIqVm3cD0D2r+VePKI7sllGvvQ0N8tGDmfUCniPUE3AFoUcQpwEPm9k+d/97lK93OaGBknTu3DmapxYREflKemoyx/XO5rje2dz6XVi1efdXjygmvb+aJ2atpGlqMiN7ZjKqTw5j+mST27ZZYPXGbY+CmU0FDgf6untxue1/JzQOIcfdyypp91vgOqCpl7s5MzuS0CyKU9z9lequrR4FEREJwt4Dpbz3+Zav3kmxdmtomF6vnBZfzaTI75pBWkpSVK/bIHsUgIHAgvIhIewD4FwgB/iiknaLgCZAD74+TqFf+HtxlOsUERGJiqZpyYzpm8OYvjm4O59v3s3bS0KPKP46ezWPzVxJ87RkTuzfnnvOGlIvNcVzUPgCGGJmae5+oNz2o4B9wNYq2r1GaCzDecCvy20/H/hEMx5ERKQhMDN6ZLegR3YLLj22O7v3lzB7Rai3oT4nWMZzUPgTMBV4ycwmEhqjcCpwDnDvwfBgZiXAX939+wDuvsnM7gF+Zma7gI+As4BvhNuLiIg0OM2bpPCtfu34Vr929XrduA0K7v6smX2b0FLMfwbSgRXA1cAj5Q5NDn/Kuxn4ktBYhYNLOE9w95djXbeIiEgiidugAODurwKv1nDMf/XAuHspcEf4IyIiInUU3WGTIiIiklAUFERERKRKCgoiIiJSJQUFERERqZKCgoiIiFRJQUFERESqpKAgIiIiVVJQEBERkSopKIiIiEiVFBRERESkSgoKIiIiUiVz96BriDtmVgSsjuIps4DNUTxfkHQv8SdR7gN0L/EqUe4lUe4Don8vXdw9u7IdCgr1wMwK3D0/6DqiQfcSfxLlPkD3Eq8S5V4S5T6gfu9Fjx5ERESkSgoKIiIiUiUFhfrxaNAFRJHuJf4kyn2A7iVeJcq9JMp9QD3ei8YoiIiISJXUoyAiIiJVUlCoIzPLNbM/mtkcM9tjZm5mXSNsm2RmPzOzVWa2z8wWmNmZsa242noO5V5WhY+v+Dk9tlVXWss4M3vOzFab2V4zW2pmvzGzlhG0TTez35vZhnDbOWZ2XH3UXUU9h3Ivlf15uJkNqY/aK6nnRDP7j5l9YWb7zazQzKaYWb8I2rY1sz+b2WYz221mb5rZwPqou5Ja6nQfZta1mj+TNvVVf03M7LVwTXdEcGxc/X0pr5b3EW9/V0ZXUc/2CNrG7M8kJRonaaR6AhOAucBM4IRatL0duB64Odz+bGCqmZ3i7v+KdqEROJR7AXgd+FWFbUsPvaxaux5YA/wcKASGEqprjJmNcPeyato+DnwHuAH4HLgaeN3Mhrv7/JhWXblDuReAJ4FHKmz7LMo1RiqD0P+2JgJFQGfgJuA9Mxvo7pWuWWJmBrwEdAWuAbYBPwPeNrMh7l5YD7WXV6f7KOc3wD8rbNsV9SrrwMzOAQbXokm8/X0B6nQfEF9/Vw66Fviw3O8lEbSJ3Z+Ju+tThw+QVO7nSwEHukbQLgfYD/y6wva3gIUN6V7Cx68CJgX95xGuJbuSbReG7+cb1bQbHD7m4nLbUgiFnX82pHsJH+fAHUH/edRQY59wnT+p5pjTwseMKbetNbAVeCDoe6jFfXQNH3Np0PVWUV9b4AvgnEj+txOPf1/qch/hNnH1dwUYHa7pm7VsF9M/Ez16qCOv+b/oqnIikAZMqrB9EjDQzLodUmF1cAj3ElfcvaiSzQdTeadqmp4KFAPPlDtXCTD5/9o781irqisOfz+egihFwHkoRYtapQ6x2mgCKtUEZ6y10dYJja1TmxI7aXDAocG2ommjMcYhaEFN1aJYWge0CEqxNmBFHBqNUFRABB4oIKKs/rH2bY/n3XvffZd337lP15ecnHf22Xuftc6+5521915nL2CEpF6dJmSNbIIu3YXlaV+tp3Qi8K6Z/a2UYGar8FGGkQ2UrSPUokez82vgZTO7r8b8Tfe8JDqqx+eJhrZJGApdzxB8ROGNXPr8tG933rYJOSH5NqyXNLsI/4QqHJ72r1bJMwR4y8zW5tLn40bd4EYIVge16FLiwtQea9O8+rBGClYLklok9ZS0Bz7UuwSo9k99CPBymfT5wEBJfRogZrvUoUeJcZI+kbRK0pSifC2ySBqKj1Rd3IFiTfe81KlHiaZ7VoBJkj6VtFzSvZIGtpO/oW0ShkLXMwBotTQ2lGFF5nx34lF8/ngEcDrwETBZ0hmFSgVI2gW4BphmZv+sknUAPv+dp2napAO6gI9OXQQcBfwQ2AZ4WtIRDRWyfZ7HjeR/A/vhUyjvVcnfXrv071zxaqajeqzHDYrzgeG4/8m+wCxJezdY1opI6pnkusHMOuJT1FTPyyboAc33rKwCxuNTwN/C/dmOAv4uafsq5RraJuHMGGwSZvbj7LGkycBs3HErP73SZaTe5iP4kPA5RcnRGXRUFzM7M3M4U9IjeM/8OmBoQ4SsjTOBvsDu+MvySUlDzWxBgTLVQ4f0MLPFwAWZpJmSHsN7e2OAoozqXwC9gV8VdP3Oom49sOmwDwAACRBJREFUmu1ZMbO5wNxM0jOSZgD/wB0cL+9qmSBGFIpgJdAveXRnKVl8K+jGmNmnwAPArpJ2KkIGSb3xkY7dgRHWvnf8Ssr3Tgtvkzp0aYOZfQBMBQ7uZPE6KserZvZ8mkM+EuiDfzVQifbapVwPquHUoUe5OhYBz1JQm6Sh7DHAFUAvSf0yn2qWjlsqFG+a52UT9WhDszwrWcxsDj56VU2mhrZJGApdz3ygF/DVXHrJN+GVrhWnoXT5sp+SNgceBA4CjjWzeTUUmw/sJmnLXPo+wMe09SfpEurUpRpNswyrmbXi97Xa3Ol8fO41zz7Af8zsw0bI1hFq1KNqFZ0oTkfYHdgCH/VbmdnAR0lW4tMj5Wim52VT9KhG0zwrGarJ1NA2CUOh63kM9049PZd+Bu6x+1bXi9R5SNoMOBX/R76ki6/dA5iEz+2dZGazayz6KLA58N1MXSU9njCz9Z0ta3tsgi7l6uoLHI8PXzYFknYAvga8WSXbFGAXSSUnzpIuJ9B2PYJCqFGPcuUG4kPbRbXJi7i/RH4Df+kOp/LLpZmel03Row1N+qwchH+GW02mxrZJ0d+NducNOCVtt+LW3oXp+PBMnk+AO3Plrsed/i7Bv5u9FdgIHN+ddMG/V74f9zYeji8cNTOVP60AHUqyXwccktt2TXm+kvS4Mlf2frz3cR4+nPxgaqMDC2qPunTBe1G3A99Pv62zgXl4r2JYQbpMxoeGR6bfyfnAa0ArsGfKc3jS5axMuR7ALGBR+m2NAKbjw6hf7kZ6jAduwhc1G477KyxM5fYqok2q6PiZdQW6y/NSjx5N+qxMSs/8yXgn4afA+/jia9sW1SaFN2h33tKPsdw2PZdnQq5cC+6UshD3iH4JOKW76YK/tJ4GluKjJK3ANHwuvQgdFlTRY2zKMyh7nCnbG7gR/9TtI9yz/YgC26MuXfDe9nPpn8sG/Dv/KcA3C9Tll/iKhq3AWnwRmNvILOrF/xeaGZUrOwC4CzcO1uILk+3fnfQAzsXXwFiZ2mQJcC9NZiQkWfMv2G7xvNSjR5M+K5fh74NVSaZFeJTInYpsk4geGQRBEARBRcJHIQiCIAiCioShEARBEARBRcJQCIIgCIKgImEoBEEQBEFQkTAUgiAIgiCoSBgKQRAEQRBUJAyFICgYSaMkmaTBmbTRkk4uUKZ+ksZKOrDMuemSphcgVlOR2uy6ouUIgkYT0SODoDkZjQcN+lNB1+8HXAW8DczJnbuo68UJgqAowlAIgi8IknpZJ6zDb2afp8BlTUuKMLu5mX1ctCzBF5uYegiCJkPSAnw999PT8LZJmpA5v7+kKZJWSlon6TlJw3J1TJD0tqRDJc2StA74TTp3mqSnJS2T9KGkuZLOzpQdBJSCk92ekWFUOt9m6kHSXpImS2pNMs2WdHQuz9hUzx6SpqZrL5R0ZQqCVe2eDEplz5d0jaTF6VqPSto1l9ckja1QflSZe3RQ6R5Jel3Scen8JZIWSFot6RFJ25UXTWNSPeskzZB0QJlMJ6d7sjbJ/UAKDJXNs0DSREnnSnoNjzlwXLX7EgRdQRgKQdB8fBtfr/1x4NC0XQuQfAZm4TEQfgB8B1+jfpqkb+Tq2RoPFHMfcAweXwA8NO+DeATTk/DIc3dIuiCdX4wHpQEYl5FhajlhJe2MT5PsD/wID37UCkyVdEyZIpPxGCEnAQ8DV+MBeWrhMjyk87nAT5JcE2ssW46+wD3AHfh9fw94SNJ4PIDTxfg00HDgljLlzwKOxfUeBewAPCVpQClDuq8P4SHkT8GDSX0deEbSl3L1DceDxV0NHI2v+x8ExVJ08I7YYvuib/gLxoDBmbQFwMQyeZ8CXgV6ZtJaUtrDmbQJqc6R7Vy7Bz4FeTvwr0z6oFT+vDJlpvPZYGE34NHsBudkeh2Yk0kbm+o8J1ffPDwUbjU5S/JMz6X/LKXvnEkrFzCnVH5UmXt0WCZtv5T2OtCSSb8RD9LTkrvO+8BWuetsAK5Nx33wAD935eTZDR8xGJ1r87XAjkX/JmOLLbvFiEIQdBMk9cbDGD8AbJS0WYo5Lzxq52G5IhuAP5epZw9J90l6J+XZgIem3atO0Q4DZpvZG6UEM/sUH8k4QFLfXP78yMTLwEBq4y+543lpX2v5PGvMbEbm+LW0n5Z0yKZvBuyUl8fM1pQOzGwBMBsf6SDt+wKTSu2V2mxRqjPfZrPNbEmdugRBQwhnxiDoPgzAe+pXpK0NknqY2cZ0uCz3skNSH+BJvOd6KfAm3rO9EB/Or1euuWXSl+BGTH9gdSZ9RS7femCLGq9VriwdKJ+nNXtgZh+7DyErc/lKDoX56ywtU+dSYEj6e/u0n1bh+vnrLK4oaRAURBgKQdB9aAU24nPl95TLkDESwIfG8xyKO0oOM7NnS4mpl1svK4Ady6TvmGTIvwwbzXqgZy5tmwZda4cKae+kv5en/Shgfpm8H+SOy7VZEBRKGApB0JysB3pnE8xsjaSZuNPgnJxRUCtbpv2GUoKk/sDIMtcnL0MFngFGSxqUht6R1AKcCsw1s9XVCjeAhbizYJZGfT1wrKStStMP6YuRQ4Dr0/lZuDEw2MzubpAMQdBQwlAIgubkFWCYpOPxIfz300v4EmAG8LikO/Gh6m2BA3FHu0vbqXcWPg1wi6SrgK2Ay3GnvK0z+ZbiveHTJL0ErAHeMrPltOUmvMf8ZKpzNb4o054U83nf/cDlksbg/gLDgO816FrrgCck/RbohX+tsBq/J5jZakk/x+/3dsBfcefGXXB/k+lmdm/ZmoOgSQhnxiBoTi7DPe//CLyAfzGAmc0BDsZf4r8HngB+B+yLGxBVMbNl+GeALfgnkuPwTwMn5vJtxB0c++Pz6y8AJ1So811gKD60fmuqdwBwnJk9VrPGncc44Gb8k8WHgb2BMxt0rXtw58ybgbuBZcCRZvY/Xwozuw04EXcW/QPukDkW76i92CC5gqDTkFlMiQVBEARBUJ4YUQiCIAiCoCJhKARBEARBUJEwFIIgCIIgqEgYCkEQBEEQVCQMhSAIgiAIKhKGQhAEQRAEFQlDIQiCIAiCioShEARBEARBRcJQCIIgCIKgIv8Fyn1KjsBCkCoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#NBVAL_SKIP\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Plot model error\n",
    "plt.plot(range(1, maxiter+1), relative_error); plt.xlabel('Iteration number'); plt.ylabel('L2-model error')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be observed in last figures, the results we obtain are exactly the same to the ones obtained in [04_dask.ipynb](https://github.com/devitocodes/devito/blob/master/examples/seismic/tutorials/04_dask.ipynb)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}