{ "cells": [ { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# 05 - FWI with total variation (TV) minimization as constraints" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "There is a lot of research on regularization to improve the quality of the final result beyond the simple box constraints implemented in the previous tutorials. In this tutorial we look at how more advanced FWI techniques such as [total variation denoising]( https://en.wikipedia.org/wiki/Total_variation_denoising) applied as a constraint can be implemented using [Devito](http://www.opesci.org/devito-public) and [scipy.optimize.minimize](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html) using [scikit-image](http://scikit-image.org/). This is a variant of the FWI with TV constrains algorithm described in [Peters and Herrmann 2017](https://doi.org/10.1190/tle36010094.1).\n", "\n", "[Dask](https://dask.pydata.org/en/latest/#dask) is also used here to speed up the examples.\n", "\n", "This tutorial uses the same synthetic datasets and model setup as the previous two tutorials, so check back if you get lost on parts of the code specific to Devito, SciPy.optimize or Dask. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Setting up (synthetic) data\n", "We are going to set up the same synthetic test case as for the previous tutorial (refer back for details)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "#NBVAL_IGNORE_OUTPUT\n", "from examples.seismic import Model, demo_model\n", "\n", "# Define the grid parameters\n", "def get_grid():\n", " shape = (101, 101) # Number of grid point (nx, nz)\n", " spacing = (10., 10.) # Grid spacing in m. The domain size is now 1km by 1km\n", " origin = (0., 0.) # Need origin to define relative source and receiver locations\n", "\n", " return shape, spacing, origin\n", "\n", "# Define the test phantom; in this case we are using a simple circle\n", "# so we can easily see what is going on.\n", "def get_true_model():\n", " shape, spacing, origin = get_grid()\n", " return demo_model('circle-isotropic', vp=3.0, vp_background=2.5, \n", " origin=origin, shape=shape, spacing=spacing, nbpml=40)\n", "\n", "# The initial guess for the subsurface model.\n", "def get_initial_model():\n", " shape, spacing, origin = get_grid()\n", " grid = get_true_model().grid\n", " return demo_model('circle-isotropic', vp=2.5, vp_background=2.5, \n", " origin=origin, shape=shape, spacing=spacing, nbpml=40,\n", " grid=grid)\n", "\n", "\n", "from examples.seismic.acoustic import AcousticWaveSolver\n", "from examples.seismic import TimeAxis, RickerSource, Receiver\n", "\n", "# This is used by the worker to get the current model. \n", "def get_current_model(param):\n", " \"\"\" Returns the current model.\n", " \"\"\"\n", " model = get_initial_model()\n", " model.m.data[:] = np.reshape(np.load(param['model']), model.m.data.shape)\n", " return model\n", "\n", "# Inversion crime alert! Here the worker is creating the 'observed' data\n", "# using the real model. For a real case the worker would be reading\n", "# seismic data from disk.\n", "def get_data(param):\n", " \"\"\" Returns source and receiver data for a single shot labeled 'shot_id'.\n", " \"\"\"\n", " true_model = get_true_model()\n", " dt = true_model.critical_dt # Time step from model grid spacing\n", "\n", " # Set up source data and geometry.\n", " src = RickerSource(name='src', grid=true_model.grid, f0=param['f0'],\n", " time_range=TimeAxis(start=param['t0'], stop=param['tn'], step=dt))\n", " src.coordinates.data[0, :] = [30, param['shot_id']*1000./(param['nshots']-1)]\n", "\n", " # Set up receiver data and geometry.\n", " nreceivers = 101 # Number of receiver locations per shot\n", " rec = Receiver(name='rec', grid=true_model.grid, npoint=nreceivers,\n", " time_range=src.time_range)\n", " rec.coordinates.data[:, 1] = np.linspace(0, true_model.domain_size[0], num=nreceivers)\n", " rec.coordinates.data[:, 0] = 980. # 20m from the right end\n", "\n", " # Set up solver - using model_in so that we have the same dt,\n", " # otherwise we should use pandas to resample the time series data. \n", " solver = AcousticWaveSolver(true_model, src, rec, space_order=4)\n", "\n", " # Generate synthetic receiver data from true model\n", " true_d, _, _ = solver.forward(src=src, m=true_model.m)\n", "\n", " return src, true_d, solver\n", "\n", "# Define a type to store the functional and gradient.\n", "class fg_pair:\n", " def __init__(self, f, g):\n", " self.f = f\n", " self.g = g\n", " \n", " def __add__(self, other):\n", " f = self.f + other.f\n", " g = self.g + other.g\n", " \n", " return fg_pair(f, g)\n", " \n", " def __radd__(self, other):\n", " if other == 0:\n", " return self\n", " else:\n", " return self.__add__(other)\n", "\n", "from devito import Function\n", "\n", "# Create FWI gradient kernel for a single shot\n", "def fwi_gradient_i(param):\n", " from devito import clear_cache\n", "\n", " # Need to clear the workers cache.\n", " clear_cache()\n", "\n", " # Communicating the model via a file.\n", " model0 = get_current_model(param)\n", " src, rec, solver = get_data(param)\n", " \n", " # Create symbols to hold the gradient and the misfit between\n", " # the 'measured' and simulated data.\n", " grad = Function(name=\"grad\", grid=model0.grid)\n", " residual = Receiver(name='rec', grid=model0.grid, time_range=rec.time_range,\n", " coordinates=rec.coordinates.data)\n", " \n", " # Compute simulated data and full forward wavefield u0\n", " d, u0, _ = solver.forward(src=src, m=model0.m, save=True)\n", " \n", " # Compute the data misfit (residual) and objective function \n", " residual.data[:] = d.data[:] - rec.data[:]\n", " f = .5*np.linalg.norm(residual.data.flatten())**2\n", " \n", " # Compute gradient using the adjoint-state method. Note, this\n", " # backpropagates the data misfit through the model.\n", " solver.gradient(rec=residual, u=u0, m=model0.m, grad=grad)\n", " \n", " # Copying here to avoid a (probably overzealous) destructor deleting\n", " # the gradient before Dask has had a chance to communicate it.\n", " g = np.array(grad.data[:])\n", " \n", " # return the objective functional and gradient.\n", " return fg_pair(f, g)\n", "\n", "import numpy as np\n", "from distributed import LocalCluster, Client\n", "\n", "# Dumps the model to disk; workers will pick this up when they need it.\n", "def dump_model(param, model):\n", " np.save(param['model'], model.astype(np.float32))\n", "\n", "def fwi_gradient(model, param):\n", " # Dump a copy of the current model for the workers\n", " # to pick up when they are ready.\n", " param['model'] = \"model_0.npy\"\n", " dump_model(param, model)\n", "\n", " # Define work list\n", " work = [dict(param) for i in range(param['nshots'])]\n", " for i in range(param['nshots']):\n", " work[i]['shot_id'] = i\n", " \n", " # Distribute worklist to workers.\n", " fgi = client.map(fwi_gradient_i, work)\n", " \n", " # Perform data reduction.\n", " fg = client.submit(sum, fgi).result()\n", " \n", " # L-BFGS in scipy expects a flat array in 64-bit floats.\n", " return fg.f, fg.g.flatten().astype(np.float64)\n", "\n", "# Start Dask cluster\n", "cluster = LocalCluster(n_workers=5, death_timeout=600)\n", "client = Client(cluster)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## FWI with L-BFGS-B\n", "Equipped with a function to calculate the functional and gradient, we are finally ready to define the optimization function." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from scipy import optimize\n", "from skimage.restoration import denoise_tv_chambolle as denoise\n", "\n", "# Define bounding box constraints on the solution.\n", "def apply_box_constraint(m):\n", " # Maximum possible 'realistic' velocity is 3.5 km/sec\n", " # Minimum possible 'realistic' velocity is 2 km/sec\n", " return np.clip(m, 1/3.5**2, 1/2**2)\n", "\n", "# Many optimization methods in scipy.optimize.minimize accept a callback\n", "# function that can operate on the solution after every iteration. Here\n", "# we use this to apply constraints and to monitor the true relative\n", "# solution error.\n", "relative_error = []\n", "\n", "def fwi_tv_callbacks(x):\n", " # Apply boundary constraint\n", " x.data[:] = denoise(x.reshape(181, 181), weight=5.0e-3).flatten()\n", " x.data[:] = apply_box_constraint(x)\n", " \n", " # Calculate true relative error\n", " true_x = get_true_model().m.data.flatten()\n", " relative_error.append(np.linalg.norm((x-true_x)/true_x))\n", "\n", "def fwi(model, param, ftol=1e-6, maxiter=20):\n", " result = optimize.minimize(fwi_gradient,\n", " model.m.data.flatten().astype(np.float64),\n", " args=(param, ), method='L-BFGS-B', jac=True,\n", " callback=fwi_tv_callbacks,\n", " options={'ftol':ftol,\n", " 'maxiter':maxiter,\n", " 'disp':True})\n", "\n", " return result" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now apply our FWI function and have a look at the result." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " fun: 368.47651052751235\n", " hess_inv: <32761x32761 LbfgsInvHessProduct with dtype=float64>\n", " jac: array([8.51090552e-12, 3.93718148e-11, 1.02249945e-10, ...,\n", " 1.89725249e-10, 7.54411256e-11, 1.66406802e-11])\n", " message: b'STOP: TOTAL NO. of ITERATIONS REACHED LIMIT'\n", " nfev: 6\n", " nit: 5\n", " status: 1\n", " success: False\n", " x: array([0.16, 0.16, 0.16, ..., 0.16, 0.16, 0.16])\n" ] } ], "source": [ "#NBVAL_SKIP\n", "\n", "# Change to the WARNING log level to reduce log output\n", "# as compared to the default DEBUG\n", "from devito import configuration\n", "configuration['log_level'] = 'WARNING'\n", "\n", "# Set up inversion parameters.\n", "param = {'t0': 0.,\n", " 'tn': 1000., # Simulation lasts 1 second (1000 ms)\n", " 'f0': 0.010, # Source peak frequency is 10Hz (0.010 kHz)\n", " 'nshots': 8} # Number of shots to create gradient from\n", "\n", "model0 = get_initial_model()\n", "\n", "# Apply FWI with TV.\n", "result = fwi(model0, param, ftol=1e-6, maxiter=5)\n", "\n", "# Print out results of optimizer.\n", "print(result)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUYAAAEKCAYAAABuTfznAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztvXuYZUV57/95p3ume2Z6YGBGQGBwvODteEPxdvAc7wY9McbEnKDRqPFINDFHFE285KhHT6LEROXnDSdqUIMaL6jEEBUTDNEEFBAdZDQhSJSL4IADM8x0z3T3+/ujqtau/e5aa9XevfuyZ+r7PP3UWrVq1aq19upa33qvoqoUFBQUFHSwarkHUFBQULDSUCbGgoKCAoMyMRYUFBQYlImxoKCgwKBMjAUFBQUGZWIsKCgoMFjUiVFEThWRH4nItSLyusTxCRH5G3/8MhHZupjjKSgoKMjBok2MIjIGvB94OvBA4Lki8kDT7CXAL1T1PsC7gbMWazwFBQUFuVhMxvgo4FpVvU5V9wOfBp5l2jwL+Jjf/hzwZBGRRRxTQUFBQSvGF7Hv44CfRvs3AI+ua6OqsyJyB7AJ2Bk3EpHTgdMBWLP+ERx9f1jnjq2f2u07uqFqfytHA3AX6wE4cGCNO3DAz7mzvuFcdJHgADRvSq053nYsVd/kZDSIA1LdJyRVLzVlP+fm1i8WBrleznOta5Oqb+tvkHNSsPdat78qUWePhXIscU6oO2DKG67Yqap3A7iPiO7NHPbN8FVVPTWz+YrFYk6MQ4OqbgO2AcgJJytnXg4nu2MPOeUfAXg7r6/av5c/AOAyPw/fcOMWd+Bnk64M0+6e6CLTLeWs2W86NttSb7dT+zmwv954TX3qWF3bpnNz6xcLg1wv57nWtUnVt/1uOefkoO23DeVkos2kKadqSoCNvgy84me+fI38Z2iyF/jdrEHDW2BzZtMVjcV8tW8EtkT7x/u6VJsbRGQcOBy4rbFXwf3g/sfdySagMxkC3OIZ45inhFMbHavcM+s/j5OrXZma5NomxPicumNtfaXq6soYdf9gw5gw+umjqe1IfGo9hjHJ9fOB6+f3m204FtfHx+2EWDdBbqSDMI0FkrCr91LSMIyDFYspY/wOcKKI3FNE1gCnAReYNhcAL/TbzwH+UUtUi4KCFYVVwNrMv4MFi/Yh8DLDVwBfxUkyPqqqPxCRtwKXq+oFwEeAT4jItcDtuMkzD55x3Trj2OH3JzqH5vxtzXkByprJ/QBMeua4f9I1ng8MEmDadzArXf03MsY2hrinpr6p34UwxxzUMZFUn+Pm2CBL6sWmGrnL4MVaFg9DJJJC29I5xRhtm7oldLyUnvJCxbCKSvxeAqzOGvTBg0V9bVX1QuBCU/emaHsa+I3FHENBQcHCcCgupUfvfhX3VfZs7I6dTmAydkxHxbxmbKbrlLExd2zd1D4AJjyDnJleU7WZ83VznkXOh2Oz/hFNGyYZbwdmWPeVb2Kbtm2d4iZGE5vsF01ssI1dDsIUF0uBUte2ifUtBmPsB02KrzamaOWHqbo2WSOw2v9PHCiMsQujNzEWFBQsKQ5Fxlh8pQsKChoRGGPOX2tfIltE5GIRuUZEfiAir0y0OVxE/lZEvufbvDg69kIR+Xf/90J77rAwmh+CWaKlplOa7I+WxWPr3dJ5bIB1TpdCBjrKmNTyq9/uc5ZKdukcL5XarjfIktMu4WPUHWtS3OQqbBaK3GVwjhJrGCY3OWizGU2Z3vSzlJ4ydbXKl85LNjHpxE4HbP8RglZ6SJgFzlTVK0VkA3CFiFykqtdEbX4fuEZVnykidwN+JCLn+Tt4M86KWf25F6jqL4Y3PIfCGAsKChoxTMaoqjer6pV+ezewA+cB19UM2ODdg6dwFiuzwC8BF6nq7X4yvAhYFC+b0WSMEH3tHaObme7Y64yNB0WMY5FBoTLry/2+7XzEMmuVLIMYeNd5vKTOrdtfKiyEOdrjqf7arpcztrr9uK4fxrgYRtoBOZ5DOV5Idcww1Ft2GG/XMUe/v9orGqFjytb2e/QxUWwWkcuj/W3ec60HPprWScBl5tD7cDbONwEbgN9U1XkRSbkZ20l1KBjdibGgoGBJ0KdWeqeqntzap8gU8HngDFW90xz+JeAq4EnAvYGLROSf84ewcIzuxBi+3N5U5sD4hurQHXuMRKSODS6E/aXqcsu6utxz+kGuG1qKyfUjW7RtBjm3DjlMsU2m2M+zHxb6ZYo5fs91JdTLFCvm6F7YIFeEjilbE4atlRaR1bhJ8TxVPT/R5MXAO7wX3LUi8mPg/jgX4idE7Y4HvjHEoVUoMsaCgoJGDNMl0MsNPwLsUNV31TT7CfBk3/5o4H7AdTgvuqeJyBEicgTwNF83dIweY5zHMbVgVB2igeyJ4zJ54l/HInIYY91+E/NYiKxqmMbFOcynjtE1aZgtFiKX7Af9aJZzGWTbsVzk/AfluvcNwhhj977ati78QJAtrolkjBPMdF8v8SyGbOB9CvACYLuIXOXr3gCcAKCq5wBvA84Vke3+8n+kqjsBRORtuDgMAG9V1duHN7QORm9iLCgoWFIMcymtqt+kJbKmqt6EY4OpYx8FPjqk4dRi9CbGeRxbDPHjQpikfrScTeyvn3BgC5FR5bYdlk3dINriHI113TnDQD9ywlyZbVu/begnDFi/THGQEGJNdoy+XDXlwswG2WKOXDFGcQksKCgoMDgUXQJH734P4OSK/Wj02gI1pOoGYYxt9YNgEHlhP20X4sWyWLaXOayv7lg/ssW6/YCF/nfkaqFT73CuNrrxHBdSLNj1VmV0wyGYc9PzKoyxoKCgwEA4uILQ5qBMjAUFBY0QYHXuTLFYq4glxuhNjGEpHX6AlNnCZnOsyfUqoF9jbbsd7w+y5FyuF6pJKZOrfOnnOoOc088Sup9zc68Pg7n12bZ17+NCYismz3HmOau8WU5QuqyZ2O8v26B8ScT/FIHxMjEWFBQUdCACq8fa2x1MGL2JMTDGgMAOj4nqjjHHApvMUb4MIwVqjsH3IAbedceGyeiGxRz7UfLk9tHPs2+6zmIoW1KMMZchDosxVv17hhjMc8brGWLIi9SofOmHMR4kOMRut6CgoF+IwOqJ9nYHE0ZvYgwG3iE37vG+fEzU5v6+jPPnxmgKCDFI6LBBjMMXItOsY75NbC8XqfZtTLFJFlfXbz+MeBCTm4XINPsx1q6TPeacM0jYsUYTH2ees8ozxPEGM50epP4nAg5BQ8ZD7HYLCgr6RpkYRwhBjvgUV5z4m9+rDm3A5Y/exzoA9vpyxgeu3Tfj9+MsgT402YFQFwLfTpvUBqlMf/ZYXTlshlqHJiPtQbAQA/KFnNtU388zqEMT28s9t8mtr9/QYSkZY5t8clyj6zlmaGWLTS6As1bGmMpKGV//EMEhdrsFBQV9Q4CilV7hGMPJDr0c8cjn3AjAR/mdqsmF/A8AruEBANzm1dN7vf3+vgnPGCc6EuWgnQuscv+cOxbSIgR2uT9KoTC/x/VTGwB3D+n6VJ3NM23347o62WYKw2BWw2B9w75e7pgGecMHkTE2uaYOcm6rbNEzxSjobAgrZmWLFrPRLDcXOky9owFlKV1QUFBgIEDRSq9wrMLZJXobxTVj7ot5BZ00E/s8M5zAfUHX4cIuBYf5UB/YIXS+nPt93cyYZ4xjnjFOePnk5LrqnBn/tZ6umKN/e0LQ3PB0A3Mc5GnHX/C6FKs5csQmG866c3Nlfv1omHPO6afNQhhh2/4gIcT6CR1WxwZT/fe0dQ9nVcQKx3xdk92ixdxchoyxMMaCgoICgyFOjCKyBfg4cDQuTeo2VT3btHkt8Ft+dxx4AHA3Vb1dRK4HdgNzwGxO4q1BUCbGgoKCdgxP+TILnKmqV4rIBuAKEblIVa8JDVT1ncA7AUTkmcCrTAqDJ4ZUB4uF0ZsYw9fLU/+f3XgsAJ877jlVk7B0DpjxApI58+umnOnHfJuQCyOcG9rGy5Qg5A7LmXm/lGF8dbhAumw61nROnYH3IO52tm1OoIZ+zmkby7CCDeSa3OS86TnG2oMEgujXbKfxOt6Iu8rf0lG+9Bh0j6UNu+eiG9tfmafRXcYYImNU1ZuBm/32bhHZgcsNfU3NKc8FPjWcq+ejZAksKChoRpgYc/766VZkK3AScFnN8XXAqbhUqwEKfE1ErhCR0/u7Yj5GnjFyvfu0bt/84KrJURO3AB0ly5hhhoE5xgxytqob7zpmzXf2R0bhM17ZMl99dU12wqYMg4uJFMvst4wxjHNs/SBIva05rnn99p/D8HNC2uW2zQoi4c1zDCsc71K+NBt0V+/2XPTee3O0HhOzGP1ppTeLyOXR/jZV3dbTpcgUbsI7Q1XvrOnrmcC3zDL6cap6o4gcBVwkIj9U1UuyR5eJ0ZsYCwoKlhb9LaV3tilERGQ1blI8T1XPb2h6GmYZrao3+vJWEfkC8CigTIwVwhfOi2DvuL4Td2zmGMfg1k3tA3plLXOJ2w5f0WDQ3THs7jb0no8YY8UQ+3UNTNXlsLDFRFPQikECXQyDMebI/EJ/kzQjh202sb9+Q4g1Hasz00ma+nQbcgcj7pRpjjXTse+7fcchclhYIhmjiAjwEWCHqr6rod3hwOOB50d164FVXja5Hpdi9a3DGVk3RndiLCgoWBoM1yXwFOAFwHYRucrXvQE4AUBVz/F1zwa+pqp3ReceDXzBza2MA59U1a8MbWQRFnViFJFTgbNxj/XDqvoOc/zVwP/Cffd/DvyOqv5nc6e4URvGWOWZBqanj3RlytG+DbMmF3idvDDebgsAkRN2rC7s00LCZ+W0yWF//dzXYjDGHJlfGJNljv1oo/uRF+awvlx22RR01hhyB6aYcvuzssXUygi6g6fML71W+pu+x7Z25wLnmrrrgIcOZyTNWDSttIiMAe8Hng48EHiuiDzQNPsucLKqPgT4HPBnizWegoKCARGULzl/BwkWkzE+CrjWz/KIyKeBZxHZK6nqxVH7S4nkCa2wjDFGqKu+5sZFL6Cfux+G1nYQ+V3TWAZBm01iitUuJDTaIDLTQRicvX7Tb9vUb9v1FqJZzu0DeuwVrRa6SQNdJ1sMFhUhxB4AYbsEkejCYt7uccBPo/0bgEc3tH8J8PepA95eydksrT1hSMMrKCjIQpkYlwci8nzgZJwWqgfeDmobgGw8WZkFdplGMXMMya/6CQtfxx4scphcjqfIQrTPViObg1ymmBMarU773u+Y6lD3bFJszMoWF/Jch52cqu8gEr1BZ8eMbHGsJaQY9MoWgxY62Cwe6LKsqAmZZ7EiZoqlw2Le7o3Almj/eF/XBRF5CvBG4PGqOmOPFxQULDNKoNqh4jvAiSJyT9yEeBrwvLiBiJwEfAg4VVVvzep1DscWrYwx/lJPmXJjTX2MOplVk1zSspKl+qouxEe6TQaYwxjb2AXgPLfijg9kDM7bhVY3ZpSXqetZpmjR9Ez6sWPMlS0uJFDteGewbVropnQFAb22uf5CkR1jY/qNgLKUHh5UdVZEXgF8Ffe9+aiq/kBE3gpcrqoX4CJoTAGf9bZJP1HVX1msMRUUFAyAEqh2uFDVC4ELTd2bou2nLOb1CwoKhoDCGEcAQfFic6LES5jNpgzHwhI6aR5hjuUqY8KYUvtNxs1txuGppW1AP0vo3KVzk9uizV3Tc714mbzP1NnGqSW1XULb/bW+jJbWbVHMc5bQdr/RfKamzHEjbC39M0m497Xlb0nBLqF7lC7BlRXyxCNlYiwoKCgwKBPjCGAW+BkdFhO+upujNmHbKl8sk4yVMJZNDmL8bfdzjKetIqWJKdrrWNbcNOY2pmhZYVzXE4oqsL59pgSqIMF1SpcmC+zAZNaa/XDOhqityfe9GIbk/bgE9mOuU/XfHUosBIiAjtKls9+/sqUns6UNGBFvNylfoGilCwoKCrpQGOMIYJZuY+5wBxujuuN9eR9fbvXlZtM2ZoxGzrMq4ws9P2s+o2G/crMyrCZlCN2WXzr+hepYZhNyZYqhjA3nq/EG05vADHf7cq+pj7ctY8yxbQpMcZ/ZT8klPXucNQbKltkNg0nG/eUyx1RddW53KDFrmgP1ssUsMx1jntMJFGGeVbzd9JyKVrqgoKDAoDDGEcRWXz4mqgvxgwNjrGSK7vM4OeWYzkSXTKc5gVAKdQFAZ4yzfkcbGH127dfb5p7uR144DNniLlMPdJhiiCwfItBb5hgzRit/7Mew28oWg0yxqQ/fZtoECqkLQwb15LUfGWOOgXedy59hiiGhVTJNQVu6gihNgdVCz9kVTM7KpRh4A4fc7RYUFPSN4hI4IhinI0cMJuKv6DjgH7n1JqCe/VWpUKMkWWN0t02lVrWoEmiNeeY45hnjRDqB1sx0h6FWLHLPZOeeoFcDnLJ9bGOVqXPqmGKPNjpmZ4EhWsZomWKcrrZOxtiENtniYQ3nhoewzhULSRUxFJvE1DnpEGKNQWczXf/iNAXZssUmGWMKhyBjLOlTCwoKmiG4ST7nr60rkS0icrGIXCMiPxCRVybavFZErvJ/V4vInIgc6Y+dKiI/EpFrReR1w7pFi9H7DozjZIYP8/vPd8VZ9/7fVZMrvJDxei+A3O3lUDNGtRYzScsiLYNMwaZaDWU4dzww1vW9jGDGs4XOx9u/VU1Mp8nGMfecOuZYyRPjTJaWKYb9faa+SSsdENjfanpRxxSb2GboZ3f3/qxJUJaTksIiJyydZYhVqLsodFhNAqvxSrbYWUFANzuse+9SKVADaj1cmp5FYyAQj+EupWeBM1X1ShHZAFwhIhepahzA+p24OAqIyDOBV6nq7VFWgKfi4rt+R0QuiM8dFgpjLCgoaEZYSuf8tUBVb1bVK/32bmAHLqh1HZ5LJ4VqlRVAVfcDISvA0FEmxoKCgnbkT4ybReTy6O/0ui5FZCtwEnBZzfF1wKm4HNSQzgrQNKkOjNFbSq8GjqEy05m6/88B2BI9r6s4CegoR8ISej9R5GJgrGt9MGOOtStf6pY7E9TE200YyVYmPZVpRVgS0l1CvgtgTrCKHje/sCSNl9J1yha7hM4x8E4NLsAqXwJWm/q10bG65XdLru/UEPpRzNQqXdwSetVURxFlzXDC0tkqBftR9Fl0K19MvMWcQCE5yqn+ltI7VfXktkYiMoWb8M5Q1Ttrmj0T+Jaq3l5zfNEwehNjQUHB0mLIWmkRWY2bFM9T1fMbmp5GZxkNmVkBhoHRmxjHcIJu/6Xec+3dAPjTh76xajLjmeEujgBg9x2pkN3dBt5rJtx2YHtr8F93/ymdYD/9Ipw7VillOp/lsfExPwZ3vfDln6+++gn3rTqm2I+Bdw9DCMqClOueZX22TQ7dSClboDmIRB1S16lzOWyubjzWj7nOVDdTXNfFGN07s2bMM8caBV8T5gxVq5QvxpgbMpQug+boGaJLoLiI1B8BdqjquxraHY7LAfX8qLo1K8CwMHoTY0FBwdJiuIzxFOAFwHYRucrXvQE4AUBVz/F1zwa+pqp3hRPrsgIMbWQRRm9iDDZV4av3Q1dczSM7bdoCCIQYDwnTCusuaJlkk1G4hc3WFsuS5sa6TXiCPGo6MIBpY/gdbw8SGq0WJq9K8oKBya2raRuf0xaY1gahTdUFueE6sx/LGI1h9yAPY5CwY6GszHPc7xaYYmyCs27M15nVR8rBALrZoZUpWuZY1XfliA7vDt1ljgnYEhl4q+o36X3pUu3OBc5N1PdkBVgMjN7EWFBQsLQ4BD1fRvd2w9fwBl/GXzwrUqxzpZuKPlyTjqFNT4Wy240rMLomuaRlAE37lfzRyxrHTbiz+So0VTTGXIaTwzIDq660001ud1Zr3KR5zo3v1cQYbVCJFGMM22Hcq7u7GEb2vtQ5JiDEpGGKgSW6EfrVR8UY81cYQS49Y/artjZQBOTLFvtJ+9AZ0CGF0Z0YCwoKlgaFMY4AFPfFC0wnBK2Nv4J1uaGbGIHNOT3lvr7zk66cnnKyrP0JO7W5Kfc5tdrHJlRt/DlznjlWssagYRyPNLWDBEsN91MX5LZijoFtbIpODts2UG2bHLEJKc2zfQ1b8kvHqHXN82UqKLFtM1Biq+4VRPjt10YBNQJTXON5n2WOAR2X0shl1KiBqzZz3droKlAE1NtuNgWjzQlMUgLVFhQUFBgUxjgCUNwX0DLG2DarjQGkmFXdVzawC29XOB9pZoO0sQowGhTJY3EYrjx7tSCvDDKjINsMjLVrvJYFVmNMdFynTbUeME0eIsGmskcrnYG6NyznzWtq0/bb9miPo3NtkjRb2tVDfL6X/a42IcQCG4xXC5Yprmmxha3TPHe1MaHF6NJK+3IhoddSGOGJUUTWA9Oqmp9/lpG93YKCgiXDCE2MIrIKZ/j9W8AjcfqrCRHZCfwd8CFVvbatnxJEoqCgoBU6lve3AnAxcG/g9cAxqrpFVY8CHgdcCpwlIs9v6gBG5jsQYR63VLBKlzhLYJ0gvsnw2wZoqMs6N91RBMxPdgvCu0NUdJCTPya0CWYfoc/pyajXqZoAEzlKpZD3ZpgRrm1927F+MUiwB3vdJkVKm6ImqbAJShfvOjrWbaoVL5d7XQDrzHS643imUEWL74m5GCmm2sxzcn7rxO+lq2B/RhDaFYKnqGqPNtAHovg88Hnvq92I0ZsYCwoKlhQqMDuWu7icX9SxtCFMiiJyb+AGVZ0RkScADwE+rqq7UhOnxehNjHO4jHY2V0msNAjsqE4g32QIbff7YEA2P4dlAE1mPEFAH/LHBDOQualOKK7Or+k/eNZIOzCeprBSdRnyUvt1z6/ueaaOVX0ptQjKHTu2JjOTNjaUw4rqmHZSYeOevjXoHm9gg71mOeNdba3LaKx86ZjwjHeVVdAIG1os3h7GM4mgIsyN504V/QdbWSR8HjhZRO4DbAO+BHwSeEbOyaM3MRYUFCw5QqK3EcK8DzrxbOC9qvpeEflu7smjNzHO0pEvQq+5CXS++HVMIMUIcrO/RcwnuO9N2NwdGUyxjj1WQW4TBrUhT8zclLuxSs5k8wen4M/FuB5aV8dgNgSwjhAEodvsJNSn5GrWPbKf0FqhDPK0/SbAcGz0vNe7BIY2IdRcb2ZGtx+7zs1b1hVgnk2co6eSKRqD7jUJMx17X233G8r4/qr7CeWML0Nucuv+F28vUKZooUiWKdEKwwEReS7wQlzAW2iPa1dh9CbGgoKCJYUitRHEVzBeDLwM+BNV/bGP4fiJ3JNHb2Kcpzssf/gabo7qwl1tNMdC2RRYoNZAuD50fZUDuObTHF6qHFfBNcZtbGKikyZh7UTzyznewOBCXdhf6937Avvb4NMXxC5t60ybcMzWr4lSOfQGTMhnjJYZ2v29URCJfd7YPGSA3Ov39425NnvXr+sq90Xn1vVvGVwKdVkkU/fZFjrMMuOZyK4h1O2dc+Of8Qy4kyu66iTusLcuF00yRqQaz0qHiGwD/h74uqpWqUNV9cfAWbn9NE6MInI8zljyvwHH4pxlr8YZSv69qi6vCqqgoGDRMWJL6Y8ATwdeLSL7ga8BX1HV7/XTSe3EKCJ/hcvA9WXcTHsrjk/dF5e5640i8jpVvaShj1OBs3FBiz6squ+oaffrwOeAR6rq5c1DVrqDFnixQWzHeLwpj6G7TUpzWqupTIcfg45czmqj69C0HLFsM5UoybIVywYts3N1jt0FRlhXbmRX1368PeXLNgYZj3tixj+b2fS3c258VbTtZYtewB/YYGBQYX9v5JIYmGJvOZWs3+PL1LG9Vf/dcss4eVplR1hpiZvlh/F2nZ423F/oM75eYIpBVjq9x9/7tLFl7ccCIWCAdeKwJkYR2QJ8HDga98+8TVXPTrR7AvAe3D/4TlV9vK+/HpedbQ6YtYm3VPUyXNbBt4jIJuBpwJki8hDgStwk+Zm2cTY9or9Q1asT9VcD54vIGnw48hRyk2P7pNuvpCaFYkFBwfJiyDLGWeBMVb3S/+9fISIXxfOCiGwEPgCcqqo/EZGjTB9PVNWdtEBVb8Ml0/qU7/cROFLXitqJsWZSjI/vB5p8Dqvk2H5QITn2Nabd23CM9LU5A+7Afznv43fj78ZjfBmYYkoLDYYxem1zi/Z2fDxfZmb3U3KounD3VhPs6rqZYSinatgfwBF+e2NNucmr+MP+EXOdczfc4diy3OErbBmyccTZYgODCbda56UTGwJP+O1xvxJY7zsJv1eIRRt7ovi66cNduWu929jllwW2/EW0pOgc88nSPHMM9YE5xixzr2Gt+yu2161RTjHGOrZpzw2aZ+jIFDtMMWij6S5TMsY65ATlSMAtpYejjlDVm4Gb/fZuEdmBW5nG88LzgPNV9Se+3a39XsdPrr+NS7RcDT6WOzah1ZxdRH5ZRL4rIreLyJ0isltE6vLAxmhNji0iDwe2qOrftYzh9JDAu9tWp6CgYLHhlC9rsv6AzeF/1f+dXteviGwFTqJ3tXhf4AgR+YaIXCEiv901HPiar6/tG5cXZiuwHbgi+stCzmfgPcCvAdtVtcF9oT/4KBjvAl7U1lZVt+Gs1xE5WWF1hw0+xZcv67SfOv7nQORvbL+6KdQwxbEchjhn7NIyEqrXpT0ITNHKBuNtywwt+9vMbdU5m/y2PbbpLpfDfDJ8i0PZORVCmvNBGGNcBwnGGGHCHFvvS8sY19OBj6E76RnjMUe6QR2zyQ/uqP90t3CU6+S2yGwhMMPbfCc7/bGwb9kmpFhlt1Y8IJaD1mmdK5niTLC1NDaKxNpnY69oPb768XEfEC42dPZSeqeV+6UgIlM475QzVNUSrXHgEcCTcTks/lVELlXVfwMep6o3+uX1RSLywxo9x6Sqvjp30BY5DpA/Ba4eYFJsS469AXgQ8A0vUH0McIGItD7UgoKCpYRbSuf8ZfXmgjh8HjhPVc9PNLkB+Kqq3uVliZcADwVQ1Rt9eSvwBZzILoVPiMhLReTuInJk+Mu945w7+UPgQhH5JyIe0JQs26MxObaq3kFkfSgi3wBe066VLigoWEoM01xHRARnUrOjYQ75EvA+ERnHBa16NPBuH3R2lZdNrsdpnN9a08d+4J3AG+nk5lDgXjnjzJkY/wRH4Cepj6zVg7rk2CKodrDRAAAgAElEQVTyVuByVb0gt68urMIFkn6Y33+Fu+ez7t2RqV7hNTHf9Y1+6onrdCXATuQQGe8OzFBnoBnn8R3zy2+rkKkMvsfTS2voKFU6JjfdS+iU+YxVlIRl8VHcAsDRfj28KZLDhrqjb3dLTLnJHwhlWELf4suwfIbOstouqe+qKaGztKtbSgfEb174nwtL57CiXG/Kw6NzwnZIS3Ok2fd6zCOPcmvOI4+9oTp1+ii3fcv6uwFwK0cDnaX0LWYfOsvtoACLj0FnSW1ztUBktmPcFfd6EU/PshnqlSx1ocXiNgHjLfWpupoZYYh2jKcALwC2i8hVvu4NeAsXVT1HVXeIyFeA7+P+FT+sqleLyL2AL7i5lXHgk6r6lZrrnAncJ0d7nULOxHisqj5okM5TybFV9U01bZ8wyDUKCgoWF8NkjKr6TRqzm1Xt3oljfHHddfgldQauBfa2tqpBzsR4oYg8TVW/NuhFhooxnKG2N9N57L0vBuAPP/u+qsn//Y0/BGA7DwYiofYe/3vEX9mA6usa8hP7LIHejGfeKGcADngCbRU1lkmuqQzBO5/uXmNt18YqWGL2V8cQj/X0L9Qfe9fPqnMmLUMMdgKBId5s9jMY417PEHf7Mn77mjJOx4i9+cOjD057QX2x1v9shwV2GCtfQp1limH/aF8e68vI4GPS193jWKek27Cl+5kHd8nY1dG6ANqJIihWxrpyXztUShjLFCulYOK9bGOMTRn/+gktl6G8USTJhFc47gKuEpGL6RYBZpnr5EyMLwdeIyIzuPdcXP/alKG9oKDgIMGIuQQGfNH/DYTWiVFVN9g6L0BdHqzCyaL8l+6yWx4NwO/+xnuqJj/ifgDcMudow/xOTzWC7XJTNj2b0mDc3+p4N5OM21hWecAzRZtJLn63qlQGVVCH7mAOVo4IESP0NO9ov393TwePvcMxoNWx9ehPfGmZo5Ux2hIqhnibH8Kd/jYCgdxnyng7MMUcr7TVpgyca63/zm/wY4qlekd6xrjastrAGIMBSJB/Bvko9LCxI2fdxsSx7sGNTbSbaFkTnHUhiEWX26KDDTI7b9MSWBOcxBiz3P0GMfA+eMOOXa2qXXaLIvLLuSfnGHi/1eyvAv46e3gFBQUjj1nGsv5WEP5SRCrdiI/N+H9yT85ZSm8Rkder6ttFZAL4DJAdCXfR4JWM8192bHDbU36vcywwtF1ezRnEdDmMMbfsqvNf/pADetK7fhltdWzoHbY7gRm6Dbqt8TZExtn+hipZ412eKQZ54U10UMcUaxjjgUjGeHtgjGHfl4EJBVIWM0bLFFuTa9BhilbWGMogs+no5+FOP7YjPdvaFKRIbS6JMQwxXO9tEY4+wT3X/WNxGLAQODaEQOsOPGHDoEEcgs1cPFg2WFYYh9MbZm7onPQcS+QSuIR4DvA5EXkeLjrYb+PMe7KQc7e/A5wnIq8HnghcqKrvaTmnoKDgIMEoLqVV9ToROQ0nZ/wJ8DRV3ddyWoWmsGMPj3bPBj4EfAu4REQerqpXDjjmhSGkTw2K10t9GduAHRPkgX4/56vb9lVtYow9wSnSItjYDTDYwwUZY6+73y982WGMVu64ec6Vk9YWMZYT3l5TBrrnmZf68s6IteQyxUFkjDHCYwy/YF+5nMIj9YPdZH+n8P8cK1UnTGnsJQ9b70a/6aiOfNe6AHZ+r+4QZhMJTXaPS6hN/pWySUzVxecE9MP+Uu9wtlY624R5WSEi2+kYdIOTOo8Bl4kIqvqQnH4aw46Z/V8AD/T1Cjwpf7gFBQWjihFbSmcrWJrQFHbsicO4wNDhPNo78sLrfRnfiU0nmpMsPFejl0oVatsE7xmTbClmjBPG0yUnhJj1hglhwWqDPKTqarxW9nmisy8iN5apHTClZYfxdj8yxjaExxrbSwb5Y2Cva/24D/j7WR3YYEorbT1q7jD7XuO94fAogMdE9+9iQ79NVCuA3rQStWiixE32irkY4lw2Qkvp21R1T1MDEZlqa1OrlRaR5zeZ5YjIvUXkce3jLCgoGGUEGWPO3wrAl0TkL0Tkv3t/agBE5F4i8hIR+SoZwWqbvimbcJbjIY7Zz3Fc6D7A43G63tct5A4KCgpWPkZJ+aKqTxaRZwC/C5wiIkfg+PaPcLmqXqiqP2vqA5qX0meLyPtwssRTgIfgVk87gBeE6LrLglk6S41gihNHd7ZL3Jwlde6SpUlwPZmuT2UR7M3aZzPwdS+x4+3QVuqCOcQC+xlT2iWaX4LO+v3UstiWFoOs8mLYfm3y39RPY5fslSLI399qe9/xqrbtufn9dXd1Qomsneh+9mtNhsROxPWUG6FXwtTF9sxx71usuIuZBt6j5BKYitHQLxofi6rOARf5v4KCgkMQo8QYh4WRUTX1oMk4dpAoxxa1uUoa2pjrVMElxnrNNjpMIzDHdBiyWJhfscq7/I3ZCNopdmTHZg2fM2Bd9gI7S+mjBlG2WIZo++oNy1A//AM5tj5zprRRx30pEfOe6Pmd0mWsYKvNrR3yC9kgM01ufjlmOouIMjEWFBQURBhylsCRwOhNjMFcJ6Cf8EsBdQE8Y9hjCxWkUZclcC5ZBnaZlFn1I38awrj7gTXSDvs5TNK+jHVMsi9YphzX1bHoxDl1v1Nd2V03wI/QdkqOk8KQMGJ2jACIyF/gg2MPcn7r3Xr/6F+nNw1hXUjxgoKCgwxDTG2wBfg4LmKmAttU9exEuyfgEvGtxiXYeryvPxXniTeGi+z9jppL7QC2+fQIfwV8yqdTyULOZ+BLOBPYK+gNWL/8aPqyLuQjZ1llU+j3pkATpBlDKnNgwcGNnIyTFXJWNRZN8vAFIKRPHRJmgTNV9UoR2QBcISIXqWqVV9rnhP4AcKqq/sRnBERExoD3A0/FhZH5johcEJ9bjVn1w8CHReR+wIuB74vIt4C/VNWL2waZ8wiPV9VWg8iCgoKDE8OUMarqzfi48T6p1Q5cvvl4cnsecH4wCfQZAcFlBLzWpzhARD4NPMucW8FPpPf3fzuB7wGvFpHfVdXTmsaZMzH+i4g8WFW3Z7RdfAQZow3gENsoTpm6Nqf6Jtg2k4ntOqbYwBDsi1blpPblrCm72vjrrG5JYNR6bBFgZYnDcAlcEGwwibjOHqurp/f3qStj2Lq52ZrJJed3zDlnZcgYN4tInOlzm88L3wMR2QqcBFxmDt0XWO0zh24AzlbVj+Mm0DgM8w24DIKpvt+N85v+R+BPVfXb/tBZIvKjtptoiq4TolSMAy8WketwS+mQ2iArSkVBQcHoow8Z405Vbc0NLyJTuNzSZ6jqnebwOPAI4Mk4a61/FZFL6Q/fB/5YVe9KHKvLRd01gDoMJUrFoiGwteDVckx0bLM51sYcU8dywo5Ztlrz5U99bUNdkN2EFy94GOw3ZVedD4g7OeH5mA2tFTsptLGjPrDc7C+HHPWw6NRvMmZK+/v556fR6iCE3Zoxv1cec6wZbdM7NwiLHASZYceGaccoIqtxk+J5qnp+oskNuGAQdwF3icgluOyAN4DPhexwPC5nfQrPV9W/Mtf9B1V9co4Spskl8D99Z59Q1ReYC3wClxu2oKDgIMcwZYw+MM1HgB2q+q6aZl8C3uc1ymtwy+V3Az8EThSRe+ImxNNw8si4/0lcosnN3k86WNEfhluKZyHne/RfzIXHcDS3oKDgEIDTSg/NV/oUHKnaLiJX+bo3ACcAqOo5qrpDRL6CWw7P48xyrgYQkVcAX8Xx/ZSd4u8CZ+AS58bBtO8E3kcmmmSMr/cDXisid9KZefcDSWHqkkBwow7L5K2+fFDU5nhf2qAO/TjiD7KUrpZkcQDhDuYSipSwNAsRovdX+84RLs4hErb3jrm2h633KwIbXzBWENnltV02+iGNB4VOZJBlDaytS2AO+jHwDqhbBcfj6ckoGMpwf/a+4/9r+7xqInnvXd+Jyrev+n1szhdXBhFIHOk6e/nZpEhpoy5LoFwb5lJaVb9Jjy9kst07gXcm6hsDRHibyLNF5A9U9b2DjrNpKf124O0i8nZVff2gFygoKBh9jIqvtIg8SVX/EbhRRH7NHq+RafYg53vzBn+Bx+G01P+sqgMnsl4wBMd6grIlhMp9StRmynv/2yyBIbhEikHmCO3tvmVf1b7reGy8m6LGcpq6rHM2p0iKMe7ydHndJhdU4rC7PB8L+rdYx1cXYsuwo3UhgncUOCFYGwU2FphiXW6WGG2sMn6MNktgbZ7p6JyeTIKe9a22bNBE5e7aDuWmdP3uic6zt79PYPQz1e/Y/XvGx3omlfBejJvcRE3R4XNcVOv+mxfoFjpivtKPx5noPDNxTIGhTYzvxwWn/ZTff5mIPFVVfz/nAgUFBaONUfKVVtU3+/LFC+kn526fBDxAVRVARD4GDOSYPRQEGeNWv/+rjuKcftyHqibXcm8ArjrmJABu/6FRRtnMa9DOGHPkQJOOua2adCGobGDS+OWaMUwxBKENjOQXnhWujQLVVkwxBEsd8+WRPq+0DUMGveHZ6ko/1CMjdjEb5ZiOEe7iTrMPvawyoCkYrWWKdXmlO/zNpX4D2OQJ2mGb7AFfHm324+2jzDl+/85NbiS7KkF2Z3u3yQq4p2KQQdYYyxjHfWkZo3/YVR5yX58TDMS6/eWc0/RfniF/H7JL4JJARP4U+DNV3eX3j8C5Iv5xzvm1OV8iXIvXGHls8XUFBQWHAMJSOudvBeHpYVIEUNVfAM/IPTmHMW4AdojIt3Fr9EcBl4vIBf6Cv9LfeBeI4BLov7KPPe5fAPjQuWdUTf7fi84E4NYxRxd2bXZf+/lZL3BK5QfLZoyRxrmSFTkGEALTrpl0lM26BKa00vuNVnqt51yBidxWWat3Ak/0hCg73JVH+QzQXWytz+yHEr3bR/u6tV75fZtnopbJxSLNuoyCTajVMJvrHBmN7TDv9rnaMkTD/irGePfogsea0n/27zrW8YRb/HtzS3Uy3OYvEMpdHAH0yhxjGaNlimPmPZmf9O/SVEJJa987m4N6iUPmjcpSOsKYiEyo6gyAiKyFfJujnLt906AjKygoGH2MaGqD84B/EJHg/fJi4GO5J7dOjKr6TyJyD+BEVf26n3nHVXV327mLgjlcTunr3e6//odLf332i06vmlzDA4GOvGes+kI39FvHDGtYYYwq2ZVvO17t1wcqrbNjDCkNwjmpwKe1fR3u+jp2vJMEbbLGRq9HWxvKwLgAT0A57DZfemoYcjff6Zl3nIva5ppuixUMCcbo/weDTeK6MLb10Ul2vHWM8VizH9f58vZj3YVupZsp3hqdFOp2VozRrULC7xbsHFOrgorhmxXFdBVUIvxIEXO07+N0zX4qF3UdBmSOozgxqupZIvI9OvYqb1PVr+ae3zoxishLgdNxr9y9cebT5+AcvAsKCg4BrDD5YS6+i/vuqt/ORs5S+vdxcsXLAFT130PgyGWB4r6MwZnoHPeVPeNl51RNDj/+FgBmph2DOrDHpFNK3XUmU4zlhrUMsSYJVorx2WASu7t0r7Ztt1yyo9le60tva7e+08fRJ7pQdkdvcoJCCWKzUN7iyxDxLtZE32bqvKwxaL83NaVrtewl3Hrq/2vSHAv7ObaIljFaTXO4z8ASgbuO8rLECdc4yHFvoXv/tkiVHawEwjGrpQ6/RYpZhd99zZh/OFPdx/eHFc14pPmd9jw6yMMtUySxb9/rfjy+GtrMs2qYLoFLAhH5nzjPmW/gqPh7ReS1qvq5nPNzJsYZVd3vfL/BO3anfd56B9cahtzfwFt8n99T1efZNgUFBcuLUVtKA28EHhmC3IrI3YCvA0ObGP9JRILP9FOB3wP+tu2knDDkInIi8HrgFFX9xbIy0YKCgiRGUcYIrIoif4Nb/+SYJwJ5E+PrgJcA23GRKy4EPpxxXk4Y8pcC7/c2RpgbqYECB+BnfqnxlVDfEVzfcbL3Fwxug2bpkjTWrpYS0n2wIQp33RI6LJmaltABVgnTe7wzyLoABtZVcHO0BAwKhZuOdOviTUc6/8jN93T762+ddw3Dkjq1lA72OHeYsmkpncrOFyMVUdsGuhhkKe331e/vPNL9+PGyODynnZUJTvfy2JbQa9jd+Q2CmY77/WJZnM0SGHJPj3sRy5g3s9rvHQKC6Adgeo/rn3H/UPb497JuuRwf6yfDYMYM4CzkRm5i/IqIfJWOx95v0hB8wiJHKz0vIl8EvqiqP+9jYDlhyO8L4JPUjAFvUdWvmDaIyOk4BRDdtuYFBQWLj9FxCQxQ1deKyK/jwpyBS7Hwhdzzm8KOCfBm4BV4Cioic8B7h5g6dRw4EXgCTtt9ic8vsytu5HNGbHNjeISXb3qDkOs9c4yzTIQv4v19aZljKqJ3j7lOYI7uiz0flDH+6w4w680t2pQtOfmFw4s3U+0HRUtH6G1Dke0yboMh8MRGOo8v1B3h68KxjRO+3GLK6NyNdzjJ/2rLGG1AithgfokYo4ZAD4e73/8XYxv9ULpdKjvs74jqXMsILRsM5b4obIV1+dtvlC1NSpeKKfr9wL7W+Pr9E67PNROd3zqsRvbu8YoZf/1qRdOUT93u2/rUSqlBUTOiS2lU9fO4SOF9o+kz8CrcbPtIVf0xgIjcC/igiLxKVd/d0veNtIchvwG4TFUPAD8WkX/DTZTf6eMeCgoKFhGK1Ip6VhpEZDdp5XDIVXVY4lgPmibGFwBPVdUQtAtVvU5Eng98DRdqvAnfoSUMOfBF4LnAX4nIZtzS+rrmboUup7fwhdsZNbnBl9Y5P4iMLHOERECIcK7/QnuH//lI5jg3Gxp3WGQMyxRTuaQt4xirmKLD3oopdPrbZfoLzCMEl1gTjScwxnWGVU75ckNP2aF/Gw7f3V1u6e4rMNWJ6Hrh2tU9z/lnMOvKufFe5jE71s2Owz/hPiNL3RsxuD01LK+tjPvZY0xtrFtfHDhhEBlbx0g/vVII15nwNDv+3cbWh9B1c378DhVzDO9lylynzp2Qmv0WDDO6johsAT6OM6RS3BL3bNPmCbj0Bj/2VeeHVaqIXI97HHPArE28par19m59oOluV8eTYnThn/tkNo1Q1dlUGHIReStwuape4I89TUSuwd3oa1X1tvpeCwoKlgNDXErP4qLcXCkiG4ArROSi2FrF459VtS4h3xNTc5OFiDwO57EXiNeGsPptQ9PEmKZB7ccqpMKQq+qbom0FXu3/8iB0u4IHFpiSmwTys9PUh69sU47owCp75DNRkAD/Na/yBQ/BBrYT3LQ3ZFXby2llmm67WyNqWaVlfzFjDMc6bfaZPma6+o7rKjmrl7s2Dd3m0rbadxtow207dmdZZV1Qh5htWmZqn7V15YtRp2nuyBNnEm3TgtZwbuW6mmrn36m5KTemPeFdm/QvbPwOW7fBQbDIMkZVvRm42W/vFpEdOEWtnRgXBBF5M3AycD/gr3BJtf6ajjKmEU2P8KE+10vPNemNL1xQUHCQQhHm5rMnxs0iEqtCt3nlaQ9EZCtwEt6rzuCx3tf5JuA1UdIrBb4mIgp8qK5v4Nm+7ysBVPUmz1Cz0JTzZWWqocZxeaPD1ByY3UbTBnoZok1tENs3hrq6Kb+SPXbsJedn8x5RSrZoYdlSpf2ci/JKezu3oA3fPz3RPY6qbPjeee16CKYbZFgTPrDBmkjrHuwx13mmaGVhSdlYhlzVYrZiat3Mzbo+xjK/Hrmgf07hGc34ZxPYfNdvVfe71YSPA5jwz2XNRDfTDliXSObQZpUQ7rcKH9cggw62jqt9eSBEB4nex+rdtczR/j+kZIxNgWrnpXqeGdhp5X4piMgUTmN8hqpaAnYlcA9V3SMiz8DpIk70xx6nqjd6Z5CLROSHqnpJ4hL7VVX9BIqIrE+0qUW2JXhBQcGhCVVhbnYs6y8HXkfxeeC8VHIqVb1TVff47QuB1V5GiKre6MtbgS/gHElS+IyIfAjY6APhfB34y9x7Hi2rTXCyqo30MsWYMaa0ztD7xWxKcZCBVAgyN8R2u8XOkLrZkmWKe6MAGJU3xB5/Y4EBW0Yc31ePnZrXrhNKhyqgbMq208pf6+pTbVL9WrTZ36V+r7o2ObZ7AT1yZf9MfP30VIdkTPsEa1MbvX7YHwq/dZMMro49W7vGuI81FRv3No5jrpyY9IFRfMl49JK3BVkeNJiEkj3ptcHbR38E2KGq76ppcwxwi2d8j8IRuNs861vlZZPrgacBSZtqVf1z78J8J07O+CZVvSh3nKM3MRYUFCwpVIXZA0OTrJ2CMwXcLiIhRtYb8C5tqnoO8Bzg5SIyiwvxeZqfJI8GvuAD2owDn7SeciLyfl//LT8RZk+GMUZvYlyF+9IHhhgi/6cYYz+JrSZzy06w/p4wYwNEAu3ViPowZF5Wtj+W7QSmGJxTQrmnZj/etmyyjjU33UIOC6xjjE3nWAwjNH8du423p2rK8C51XdedtNf85sFrJchbZ7tYXxp1LDPWSo8buWPQYPfKGqMbs+9qnawxRlZoMmF+bjhThap+k66IvMk27wPel6i/DnhoyyX+DfhzEbk78BngU6raVyxGKDLGgoKCNihOYZXzt9xDVT1bVR+Lyy99G/BREfmhiLxZRO6b20+ZGAsKCpoxLzA9nve3QqCq/6mqZ6nqSTjvul8FduSev3LuJBdhKR2WPZtNCWmlQLxvlxxQv6yqSreEXj3VMctYN+UNnye6XeOskXMK1sWsWkrP+TBkwRQnCkXVsyy2Buy7zH5cZ5fbdQqbphwiC8k2l/Om9bP8rhOL1IlAYtMs6xqaUuDZ6/rt+clgDuTKtRPdxumx65wVjwyCHpOfELLMm10diEQ7PXmqrUtsjhKm7tgQMg0uJXxA7afjXJGfjIvk/Zbc80dvYiwoKFhahJTFIwCviX4uLof0t4FPA6er6l2NJxqM5sQ4Tudrf4wpoZcx2i9mFmPsZojBAHrt+g5jtAbO1izDunilTDo6gWq7DZQr84jpyC192pSWOVqlDPSyyTrFTcrUp4LNEm2NmReiHYmx2hyz+w35l62JlmWDMWMM9xivMuK+Uu9HpdBwYwoBRCrD8rHu0GLQq4wbqxhk+r2I34+6oBWVK+JktzIGYL6NMdr/g9SxGnOdUZkYcRkBPonzx/7FoJ2M5sRYUFCwdPBB80cBqvqkYfQzuhNjYADH+/L+0TH7pbQsKMkYXaNJLzdc55licIuz8kOoD0hrkWaKxkzHyxaDu18lW0wZNVvmWFfG23WssmKXIYRdnNsgeGoFhhiCX1nmGNOJ3P+gOEBTHUNcW1MCVfgt7/66K6QC8IebciyH332P2W9iT+bZH6jMqXyAi/XrqEPH9S+9omh6PyyqcxNZK+eDvNEyxxxW2ASF6LU/JDC6E2NBQcHSYLSW0kPB6E2MIa90+Bp6pvigR3SCfocQVLfNuGRHe70rXZDb2cAJ0AmesG6sOyxXnRwR6rXOKQ1lXO+OGW3mbCj9OVXZ1UF/ZbzdxCqBDhuM/flvrzm2z5QxS2yTPzbJFusYYijj4MshUIphqLOeueWw6LbnlzqnKoPLpmFuiVgL4TfOXWFAc8oE6Gin43c4sNgexthk6G1/jtGXMQ4FozcxFhQULC3KxDgCmMPJhvwPNXV/l7jwTZEv+Rd4NgDXTDwQgNsmHHMMmt+UU38dQwz7TV95+1WvtVFM2LgFrWawi6u++tMm6RHUM50ctLnqVcwrRTfrtNK2jLfbZI0pGaO9ju0jdU5gkzUPo59nlMO8qyDI7veZH/cyRrMqAZib9EmvQgAIE9Q2x4U0FbA4RixjrNwEq2C2LbJGaJbFBpSJsaCgoCCBMjGucMwCPwOudbt7fng3AH76iE5CwvB1DWwvhO23bC+lYbayxH5sEa1ssSnQ6t45H7bfBouoSrrLcO9xaeub0BRAA0hrfsN2YIFWa9x0oTbGGB/PiaRqsdqUZkxN99v21jcxKyuv2+OuH+4mZozBwiAwx/0+LUawdKizd81B9X5GycUqb5jAIusYY1MCrRTmyWOWBxFGb2IsKChYWpSldEFBQYFBmRhHALM4N7cQ4vKvXXHmxvdXTY7cehMQZajzsEqX/V22FWFZHTK2dS+Lw7kpw1trWhGWzp19nyd5pmMEPGOW0PNVdG6TLzjHzMQiFWeyzt2t6jOMLVak1JjENF7QLrtzDL7Dkj3XbCe+pl1Sm8NNS8U6U5Ucs6ee/nz07yjfzrQxwQpL3bC0roy0x3qVMTm5cizGfX+1SphUXqO6gCsxDsGJsYQdKygoaMds5l8LRGSLiFwsIteIyA9E5JWJNk8QkTtE5Cr/96bo2Kki8iMRuVZEXjeUe0tg9BhjwA2+/LovZzsBBm5/zHFuY6uv8MECVm10ATZC9rd1UQixsTHHmKx5Tg568iLPdWeos+wQYpc//zWvy98SR+O2LDKHMdYxRZs3u3INPLKm0xQCS4sz5qVMeOILpQZZ5wp4mNmPs18eZupWd3dh73dYypc2ZUX0HjLrLn6gUsI4JhcYZDDO7rj3RaY+CRaZQrwq6vQTmKm73nwTY8yZAYbLGGdxAR6u9OlMrxCRi1TV5pX+Z1X95bhCRMaA9wNPxc0A3xGRCxLnLhijOzEWFBQsDebp/c4NCFW9GbjZb+8WkR3AcUDO5PYo4Fqf4gAR+TTwrMxz+8LoT4zX+/LSqC583UNgia2umN/sUrtNb/TllHbO8V/v1VW+5e5PZJDfzCbCt8+ZPM6VkbbN8zwdsYncQBD9uLKl0BbI1TLHPXFor8AerRwvMLjdZh96DbxThuPxgJr6t0wxvo5lkx51GQ1zGGOTS2COzLIW3fLHef+uhRVFWMGMR8baVX7vsVC0U7bAGEM/c76sgkvMGvOdeLsunzo4xpi/gNosIpdH+9tUdVuqoYhsBU4CLkscfqyIfA+4CXiNqv4AN4H+NGpzA/Do7JH1gdGfGAsKChYf+Uvpnap6clsjEZnC5Zv/J/IAABc+SURBVJY+Q1XvNIevBO6hqntE5BnAF4ET+xjtgjF6E+MYjkBYg9M4OOv1vrQhqEKIsiqUfcSOfLa1kHXtQF1ahBhthtaWeTRpOXOCPeTKGGNYphjOsUwxDuQaULHHwM6sJtgagMfbqZBkbYO0/Yb6FDsM236MdTLFnDd8kGAcg/znjIfn6S0RPFOsZM9RQIgxuzKpYY5d+77Nmih4LXSYaaUxn4zf+zA2U8YYslZaRFbjJsXzVPX8nstFE6WqXigiHxCRzcCNwJao6fG+bugYvYmxoKBgaTHEiVFcUuiPADtU9V01bY4BbvG5pB+Fs565DUd/ThSRe+ImxNOA5w1nZN0YvYlxHKdlno32oYbx+DKwySbXqLZc1DmoY4w5QQnq5IeDyBibbPZsmXpuFpVmPNg6NjHGoKG2AShyBKFWOx2ul7Jj7NNucSHMEXpXKIP851S2o2nmGMuqZ/z1Jgz7a5I5VnWhjZcxBhnmdGXX2JCLOnVfw3UJPAV4AbBdRII18huAEwBU9RzgOcDLRWQW93KdpqoKzIrIK4Cv4u7yo172OHSM3sRYUFCw9BgSY1TVb5JM3tPV5n3A+2qOXQhcOJzR1GP0JsbAGO0XLmY+NrFVXYyC1FdwEMaYK2tsYox1ZVNqgxztdB1zsuzZ7qdQXbdOixxvW6aY4wFTZ89Y490C7bLFJtlwP88+wMqtB0H1rENoue5gxRDZIprgyjmwASZ6PGImI4+vIG9cQhnjKGD0JsaCgoKlxQglwxoWysRYUFDQjP7sGA8KLOrEKCKnAmfjBKUfVtV3mOMnAB/DGdCMAa/zMoR6rMItkUNO4FTeYJtj2BozpwTN1sC1zQi46dgwl9Kpc9rG0rSkrjP4buszPqfHlChe4tYFj2hSwtiHnZFPepAltEXbEno80XYx4hKOe8Pv2c7sY014rPnO2FjvcwyBJ4JranAXDOY7wTnhQGQWVClimpQvZSk9PGT6Nf4x8BlV/aCIPBAnVN26WGMqKCgYAMrQXAJHBYvJGHP8GpWO9e7hOPefZghu1Dav9OaoTWCRbWyiHzex1H6bS94gX9kcN79c5phC27Nocg2r6yOlIOpRlCQUJ/0iFfwglyn2w/T7YYd1ir2mfnuu022+A70mPHOz3dkIbTi87iGZda8nmxMmiAV0HBoaXQPLUnqoyPFrfAvwNRH5A2A98JRURyJyOnA6AJMnDHucBQUFTShL6SXHc4FzVfUvROSxwCdE5EGqOh838k7o2wDk8JNd5IfAGANT3BqdEOoqcx0fLGK8j1+3yutssvWljH7bjLVTLnyWjViDdVvfD1KyMctsbNscxmgZYoqltZkQ9eMhmMPwB5ExLkS+OwirrEMP24xDlnWb8AT5YI6ssdO9CdTszXe6clGHABM2P0yMMjEOFTl+jS8BTgVQ1X8VkUnctHbrIo6roKCgHxRznaHiO7T7Nf4EeDJwrog8APe9+nlW72HkCcY4ecztAGw43Pmy9ZW/d677S90JNuvkP1VIMai+6lWw2boAEE1GwaGuJ9VA65AXhhyXuTqNdR1zjNvkGr33M7YUY1yIK2A/jHEY2mjLEJv6NJpqK2sMmBvrvUHLFKt9n50w/l1XhWC2Uw2MEYqMcVhQ1aRfo4i8FbhcVS8AzgT+UkRehfsuvcj7RBYUFKwUlPSpw0XKr1FV3xRtX4NzKu+jU9wXN3zZPGM8fOvPqibHTjjlts0RbREntqoSWo15hjjmc0FP+NzQ6/3+TIcx1ia0GjduVjmok9ul2NhCUCfDbBqrdRdsklvm2HD2i36CYvRjxxiQIz+sewaNrC+zbeo5Gk21XckGOeFclA4haJ/rEmmF+rmxOCCuDzAx6YI318oYy1K6oKCgIEIx1xkhhC/bfdxn+GETV1WHNviQ+zat6Wy1P95VD7DfpzitvAboTkZU7U/UU58qrFOwR5vNcNBfDAxi15jaz5G92f2FyBb7QZscMke2WMeerdw3dQzTZiFyycb3IrxD7Vrq4OkS3uHx6h0Oto/h3e4MMtg2TlvvMIuilS4oKCiIUMx1RgCK+xL7kR9/nLMh/y3Oq5pcy30AuIljAdjtP4P7fODTFGMMmKt5JHOGbXYdq2weQ/KrGtvHYflK5zK5JvTzy9fZQKaul8sQh/2PVnc/TQyybgxJmZ8vrfXAQtBkrdDTxifU8rv7k40d1kyEU40dY2I9vGbC99Tk+XIIKl9WLfcACgoKVjgCY8z5a4GIbBGRi0XkGhH5gYi8sqHtI0VkVkSeE9XNichV/u+ChdxWE0aPMRYUFCw9hsfwZ4EzVfVKEdkAXCEiF5ngMiEIzVnA18z5+1T1YUMbTQ1Gb2Kcw+Ug8dY5M17R8V/5l6rJLRztm7qlbVhC7/b5iWe8omV/5LQf6sJSOShjQv/BTGcmMvCeDuY5IURUMLEIOVLC8sPux9v95JVeyNI5wC4Tm96A3LY5ARqa0NZm2KY+bWi6nl1SD+IamLOEtm0r+CX1eLfht9sO+aS9gmbMuA8aZYzb9jfQFkRiSOY6qnozcLPf3i0iO3BxFa4xTf8Al0nwkcO5cn8oS+mCgoJmBHOdnL8+ICJbgZOAy0z9ccCzgQ8mTpsUkctF5FIR+dX+rpiP0WSMu4Br3e7P/8ZF2/mN3/xs1WSjTwsYWGBgjHt9WTHGuQ5j3D8dGKGrq3LxBoZo3f7ibcsacthfXY7ouv24bhCFTRtymGMOBmF3gyhoFqLUGUagjjrTnmEFnqhjlVV9r+F3yBMTTHiCu2CTEmYCo3xJmev0p5XeLCKXR/vbfBCYLojIFI4RnhHnkfZ4D/BHqjrvsq124R6qeqOI3Av4RxHZrqr/kT26TIzexFhQULC0mKefQLU7VfXkpgYisho3KZ6nqucnmpwMfNpPipuBZ4jIrKp+UVVvBFDV60TkGzjGWSbG6ut1vd//iit2TD+80+b+vtxcE24syGWmo69RLnNLmXDUtWk6t+2cfsx16o4PgoW+EW2BbgeROfZjFlRnWtQPFjLGYZm1tJkfGcNv6LimhqyAVtaYMtep6ppkjDA0zxdxs91HgB2q+q5UG1W9Z9T+XODLqvpFETkC2KuqMyKyGedO/GfDGVk3Rm9iLCgoWHoML7TLKcALgO0iEtzV3gCcAKCq5zSc+wDgQyIyj9OPvMNqs4eF0ZsYx+iWg1zvy29GdSGexDH+qxpCKjXJltpYXo7Mr66+icktl7vdIAbe/fQxDLlkzrOoY4pNzHEYpie2j37cCOsQP7O2gBPVfidlRDD+nhk3aRCq0GQNZuE5gYqHAFX9JsnsZrXtXxRt/wvw4EUYVg+KVrqgoKDAYPQY4zhwTLQfvnB7orobfBm+ujZtakA/GsQc1pe7P6y2Oef0i35YYFM4sJz+LPp5jnW/1zBkjCnk9tfkRti2khjErrEreK/bCcGUZ4KWumKOvbLGyo5xqQOdrHCUx1BQUNCC/tTSBwNGb2Jcg8seY2Ui8Z1YW6zwRd5DPXK1wzk2gv2wlcX0+ugHOUwxJ+VAW3854cBy9vuR47ZhmM+4aUUxzMATqf3x7uC2+z1TDFrqNetdUNooOUdvP8nf59CLVDt6E2NBQcES49CLO1YmxoKCghYUxrjysRqnfNno96vc0VGbtrtKLbusC1adYLwpwrXtaxhGxjlxBBdyvWEsoYelhKkzp0o9i1yRxzCM3lNYiLhkIUvrrOfqltTz4z4qvc/rElxgJ8Y6Zjt5ypcyMRYUFBQYKEX5stIxhmOLx/v9kFc6VrhYBlVnepMTECLHwNvuL8Slra2vFJrYZb9oYt51TLGJMfajhGljwsNifcNg9AsxQ1qMaOAxqmfvjL9DeLyghJlZ31G/9LgJ1jLGImMsKCgoiFCW0isfq3Bf2iBjDMbemzs/3Cqf+Wy+ysHiSxs6rJ8QYjkugYsdEGIQI+M29GO0PQhTHMRsJwd1bHIQM6thGoX308diM8eAcSdbnJl2/xdrJjsyxlje6NqmOiiMsaCgoMCgMMaVD8GNOnxlN7qwH1Obd1VN1q7vFhTPzTnGGILRhpy8IUwTNASkrWOQqbphss0c1LVNaW9Tx9rqc5liSi45DK30MJGj0V4ud8KUxYM9dxA5a+hvj9NSHxhfC8B+r6WGSN7Y2G9hjAUFBQUGxSVwdFC5BLqvX8wS17G3q2lwng8uUVXCq0jWElIZhNQGByoGaRJd5SS0sqwpJTuyDHGx2ErbL9wP+8txCcxliv24Bg7yluZotOva5PwWg1gLtLmQxkFOmlhk7jiCC2z1HrrV0Ey0Uqrkja2MsSylCwoKCgzKUnrlI5IxTk45dhizxDV0M8OAWbrTSaYQEgqFsPDzIS3CeCLYbb8eNv0EQejHbrLuuk1j7If99RM8YjGRkhMOg03W1TfJauvQD5O0iH/rYT5j09eBybXVoUre2CjjPvQYYwlUW1BQ0IIwMeb8NUNEtojIxSJyjYj8QERe2dD2kSIyKyLPiepeKCL/7v9euKDbasBoMsaCgoIlxFC10rPAmap6pYhsAK4QkYts7hYRGQPOAr4W1R0JvBmXRVD9uReo6i+GNbiA0ZsYg4G3dwHceLgz09nEzp6mu9kAwH6TR9rmkIZI2VJnDJ5a2u4xddMt9YPklc5xJ2yCXRYOEhDCts2tzxnXsM/JGWPuEnohz7mfc3JgFXj9mO/Y33aykydm7+S67v6T/Q1PK62qNwM3++3dIrIDOA6wSa3+AJdi9ZFR3S8BF6nq7QAichFwKvCpoQwuwuhNjAUFBUuMvmSMm0Xk8mh/m6puSzUUka24vNCXmfrjgGcDT6R7YjwO+Gm0f4OvGzpGd2L0ypcT/HN6GFdVh37KFgD2+0jGuzgCgN27HIM84MtkXum2Mv6i5p6T44LYjyFvGzvpx5g6hyG2meA0KXuWSpmZq1CJ64bRfxP6vfec9nV5ZOJjAZWBt9mP7n8+MEa7yulCX0vpnap6clsjEZnCMcIzVPVOc/g9wB+p6rxLQ730WLSJUUQ+CvwycKuqPihxXICzgWcAe4EXqeqVizWegoKCQTFcrbSIrMZNiuep6vmJJicDn/aT4mbgGSIyC9wIPCFqdzzwjaENLMJiMsZzgfcBH685/nTgRP/3aOCDvmzGPO7r5r9sm7gNgNfw51WTs/gjAK5nKwC773ACyQM/O8w1COLIOAdMm1wwhzHmBpNoOzYoBgn/1W87aGasK9XcbSWsjYYpV03JvAeCZ2TBo3bhjLH5am62+wiwQ1Xflbya6j2j9ucCX1bVL3rly5+KyBH+8NOA1w9lYAaL9rqo6iVehlCHZwEfV1UFLhWRjSJydy+cLSgoWDEYqkvgKcALgO0iEuRfbwBOAFDVc+pOVNXbReRtwHd81VuDImbYEDcvLQ78xPjlmqX0l4F3qOo3/f4/4OQKlyfang6c7ncfBFy9WGNeBGyGhMp8ZWKUxgqjNd5RGivA/VR1A4CIfIVOSOg27FTVUxdvWEuDlbDAaIXXam0DEJHLc4S7KwWjNN5RGiuM1nhHaazgxhu2D4aJrl8sp+fLjeDVxw7H+7qCgoKCZcVyTowXAL8tDo8B7ijyxYKCgpWAxTTX+RROtb5ZRG7AufKshkrAeiHOVOdanLnOizO7ThqLrmCM0nhHaawwWuMdpbHC6I13qFhU5UtBQUHBKKJE1ykoKCgwKBNjQUFBgcGKnRhF5FQR+ZGIXCsir0scnxCRv/HHL2sxJl9UZIz11T7+3PdF5B9E5B7LMc5oPI3jjdr9uoioiCybmUnOWEXkf0bx/T651GM0Y2l7F07w8Qi/69+HZyzHOP1YPioit4pI0i7YK0b/P38v3xeRhy/1GJcNqrri/oAx4D+AewFrgO8BDzRtfg84x2+fBvzNCh7rE4F1fvvlyzXW3PH6dhuAS4BLgZNX6lhxLqXfBY7w+0et5GeLU2q83G8/ELh+Gcf734GHA1fXHH8G8Pc4v8HHAJct11iX+m+lMsZHAdeq6nWquh/4NM6FMMazgI/57c8BT5blCcXROlZVvVhVQ+6FS3E2m8uFnGcL8DZcoNDFSgOfg5yxvhR4v/pgpap66xKPMUbOeBXwTvscDty0hOPrHojqJUCTS13ltquqlwIbReTuSzO65cVKnRhz4q5VbVR1FrgD2LQko6sZh0dbjLiX4L7Cy4XW8fol0xZV/bulHFgCOc/2vsB9ReRbInKpiCynl0bOeN8CPN+bsF2IC8i6UrFk8Q9XGkbCJfBggYg8HxdS6fHLPZY6iMgq4F3Ai5Z5KLkYxy2nn4Bj4peIyINVdVfjWcuH5wLnqupfiMhjgU+IyINUdX65B1bQwUpljDnuglUbERnHLUtuW5LR1YzDI+naKCJPAd4I/IqqzizR2FJoG+8GXKCOb4jI9TjZ0gXLpIDJebY3ABeo6gFV/THwb7iJcjmQM96XAJ8BUNV/xYWczQ3QsNQ4dN12l1vIWSP0HQeuA+5JR4j9X0yb36db+fKZFTzWk3BC+RNH4dma9t9g+ZQvOc/2VOBjfnszbum3aQWP9+9xQZkBHoCTMcoyvg9bqVe+/A+6lS/fXq5xLvlzWe4BNPxgz8B9/f8DeKOveyuOcYH70n4W51L4beBeK3isXwduAa7yfxes5Gdr2i7bxJj5bAW39L8G2A6ctpKfLU4T/S0/aV4FPG0Zx/opXGKqAzjm/RLgZcDLomf7fn8v25fzPVjqv+ISWFBQUGCwUmWMBQUFBcuGMjEWFBQUGJSJsaCgoMCgTIwFBQUFBmViLCgoKDAoE+NBBBHZIiI/9vl3EZEj/P7WRbrey0Tkt/32i0Tk2OjYh0XkgUO6zq+KyJv89rki8pwB+7mbz3hXUNCIMjEeRFDVnwIfBN7hq94BbFPV6xfpeueo6sf97ouAY6Nj/0tVrxnSpf4Q+MBCO1HVnwM3i8gpCx9SwcGMMjEefHg38BgROQN4HPDntoGIbBWRH4rIeSKyQ0Q+JyLr/LEn+1iB2328vglf/44opuSf+7q3iMhrPIM7GThPRK4SkbUi8o3gRigiz/X9XS0iZ0Xj2CMifyIi3/MBII5OjPW+wIyq9uRkFpG3eQY5JiLXi8jb/fUvF5GHi8hXReQ/RORl0WlfBH5r8MdbcCigTIwHGVT1APBa3AR5ht9P4X7AB1T1AcCdwO+JyCRwLvCbqvpgnIvby0VkE/BsnHvbQ4D/Z675OeBy4LdU9WGqui8c88vrs4AnAQ8DHikiv+oPrwcuVdWH4mI/vjQxzlOAK22liLwTuBvwYlWd89U/UdWHAf/s7+M5OFe2/xudejnw32qeSUEBUCbGgxVPx7l6PaihzU9V9Vt++69x7PJ+wI9V9d98/cdwwUzvwMVl/IiI/Bouq2MuHgl8Q1V/ri483Hm+T4D9wJf99hU4v12LuwM/N3X/BzhcVV+m3a5bF/hyOy6o6m6/fJ4RkY3+2K1ES/6CghTKxHiQQUQeBjwVx5Re1RBY1PqC1vqG+gntUbiAwL8MDEuBcSCa2OZIh8Hbh/OLj/Ed4BFByRQhRC2aj7bDfuh70vdZUFCLMjEeRPARzD+IW0L/BHgnCRmjxwk+HiDA84BvAj8CtorIfXz9C4B/EpEpHEO7EHgV8NBEf7txIcssvg08XkQ2i8gYLh7hP/VxWzuA+5i6r+AUS38nIqlrNuG+QDLHSUFBQJkYDy68FCdnu8jvfwB4gIikAuP+CPh9EdkBHAF8UFWngRcDnxWR7TimdQ5uwvuyiHwfN4G+OtHfucA5QfkSKlX1ZuB1wMW4iDJXqOqX+rinS4CTbNoKVf0s8Je4WJFrk2em8URguSOTF6xwlOg6hyC8XeOXVbVJBrliICJnA3+rql8fQl+XAM9SnyOmoCCFwhgLRgF/CqxbaCcicjfgXWVSLGhDYYwFBQUFBoUxFhQUFBiUibGgoKDAoEyMBQUFBQZlYiwoKCgwKBNjQUFBgcH/Dy/GyjkgfNiXAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#NBVAL_SKIP\n", "\n", "# Show what the update does to the model\n", "from examples.seismic import plot_image, plot_velocity\n", "\n", "model0.m.data[:] = result.x.astype(np.float32).reshape(model0.m.data.shape)\n", "model0.vp = np.sqrt(1. / model0.m.data[40:-40, 40:-40])\n", "plot_velocity(model0)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAASMAAAD8CAYAAAA8P8JjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztfXvQZVV152/1i6abR9PQaREYGyeUKaLlaKGSmJmxxEwYouJMWShJHCZhqitVZkTjjIKZFGNNMqMzlI9UUs50BSNGCmSUDMhYUYMQkzESG58IokYQmneDPAR5dPeaP85e3P2tu9be+5x7v++e7l6/qq/2Pefss8+++55v799eT2JmBAKBwKKxatEdCAQCASAmo0AgMBLEZBQIBEaBmIwCgcAoEJNRIBAYBWIyCgQCo0BMRoFAYBRYlsmIiE4noluJ6AdEdP5yPCMQCBxYoHkbPRLRagDfA/DLAHYB+CqAs5n55rk+KBAIHFBYswxtvhzAD5j5hwBARJcDOBOAOxmtIuLVANan4yNT+Rya1HkwzZlPpOOnU7knlXtTuS9rlxvLHOSU3nXrmnfcAm9psM7rc33uHVJnf0TLb+DV6fP7eXVL4+r9flyoo9+/VaoEgNWpfCaVeyflbmbeAgCnn3467969u9C7CW688cbPMfPpTZVnwHJMRscBuDM73gXgFboSEW0HsB3oBvIoAC9I116bynevm9T/86e6cmc6lgc8lMpHU/nE5JZnf4wn1fHT6jjH2lSuU8er1Hk5lh8+/2xdy8/n2KeO9zrX9xp19LWWe/XzvOf2wSz3zgt6rL3z1m/g1bHa9J7jyTus8da/l/6d8vdSTyiyYMt7uTGVG7J75NwDqXwwlY8AP5I6u3fvxs6dO9ECIjqmqeKMWI7JqAnMvAPADgBYR8RrMRlgmWBkAgKA+1IpP/rhqZQf6ZBUZrdMTTq6tCYlXedxdawnq/wl9K6tVqUFqetNFqV7NEqTg36OV7dlgunT1yGoCTRbJgvdhnWPN/mU6rYir19bCEoLz1pVrlflRkwgOwt5dx8zn8qY7CvGgeWYjO4CcEJ2fHw6FwgERgPGZM8wDizHZPRVACcR0YnoJqE3A/i1lhuFfQi9vKlQR1YIYUjr1XVgmvb2YUZ6a+cd56tYbYtnMaTSlmAoSm1pxrNanbcY0XJs7fTzLchzW9hN61bLkq3U2rDg9dvbKufXNBOS0nof5TnyTmlGtEGVwIQl6S3dUhwEzIiZ9xDR7wD4HLqx/Cgzf2fezwkEArPgIJiMAICZPwvgs8310a0WwjYeTmU+o+vZfZ1z/HR2TlYcLQTUzMiSM61SdWvH+WfZq3sypPy7yMqn2ccQAzDNclra0s/12sgxT4F13larvKeP3KePMLqFgXkYonjwGHr+XGE3IhP1ZEbrMYFmRPbvf5BMRoFAYOyIySgQCIwGMRmZ2IfpbZQl659F0OsJEHU/LGhDMmt75bUvNFy2a/nWTtP5vnYsFnQb1raqVQCbP3u5Vfn6eQK9PSttvWpq+iFCaQslQXV+Pa/3jKrrbc/yLdch6lytBKbFA/b32oelAorFYzSTUSAQWEnENs2ECLBbVJyGifuS85ZAudX4saWOJ3zMr6009Gov/WhRm7ewgZo5wCxoMXPwGJGlpu8jwJ4HPCt4S4DtKVI0kzkEE3gq/RZmVGfT85mMiOij6Bwn7mfmF6pr7wRwEYAtzFz0P4kQIoHAQQlhRi1/VXwMwJTvGhGdAOBfALijpZFRMCOBrBiyk83N2LWRoceihjCjltWrtOIJZPVdq45Lso7lWLm1jKeFIel7S9Dt9mFILUaINSZUMh71njNkfK3vVfML9MxJrHMeI8rNVrSa3lP15+YiZVmRYH7bNGb+EhFtMy59EMC7AFzV0s6oJqNAILBS2IfldAchojMB3MXM3yRqi38wismI0a0a4nEvM3ppqDy2Y61EXh3rHkFtNV6PadRW51kM8Ep1NeT7WFqw1nv7PK8PWsaglRH1kRn1QSmSQl9GZLkn1RjREG1azox0WBGfuTYzo2OIKHfx35Ec3U0Q0QYA70G3RWvGKCajQCCw0ui1TdvNzKf0aPwfAzgRgLCi4wF8jYhezsz3ejeNYjIidLO4xF0RWZHt4LcUngwpP6dXJM1urBV2SPiKIXW958zCRrQsxwpjsRyOuSXZUYvdTysTGiIzKmGIVrCVEVmO1JrNaEaUa9M8BqSdsXM5k/7utpZ3+VT7zPxtAD8jx0R0O4BTQpsWCAQMzE+bRkSXAfg7AC8gol1EdO6QHo2CGa0GsAn+agP40RlLq7JmQHqFtZhXSyC0VsxiNd33OtBmD+RZU/f5vkPa8JhLH21aibUOkRW12oVZVvae7MjTmAE+m9GMKGc5+h7NkCyborYAenPVpp1dub6tpZ1RTEaBQGClcXAEVwsEAqNHuIOYWANgS3as4wIBE6G2ptaarrYYmWknV0vA6wk1S1uRmuCwZXvR97pV19uu5dfmISwcsrXT97YIoVvcQWp96xPTW2+9rJhEnhuSp74HfBW+t23L7/G2ZyWxQskhfOKENR6MYjIKBAIrjWBGJtYA2IzJDC4s6MGsjhhESuYDLZQWtMTAlvNWHGtBLVxFSQALdU0zlj7squUer43S9/IMI62VdkgfvOcKSqmevDqzqO/7oMXZ1WNRnurdOlcr83v071JiRF5s7WnEZBQIBBaO5XUHGYJRTEaEjukIIxITzXzWf3EqD4cNy1G25jJi3VNzkNVtWjGc9Yqu9/fWPd4qX7qnhpIcwWN6peBqfVD7PlbbNfX8EEap0aKm947zz15gNP2bW3Hc+zAjbYbSHh7EdkmZILZpgUBgFIjJyMU+TDLJyuz/r7PrwoiEWIqmTcKNPK2uA+350vKV70mnzpBstLrUDCk/pxlLKQTGLDKceWrvWlbnWRyAl0smpOFpz0raNC1ja2E53jUtS8q/t3etpFHUjMgfx5iMAoHAwhHMyISI0u5Mx29M5fsPm9T5fz/pyh+kY9G0iZZNZ3wF6uymdI+XQdY7X2pXn7fcUIa4jngMaQhjKWGIM2+t/Xk7RQ7JY1+TFZUC9mntWY31WOc89xCLTdW0u5bGT97lgzajbCAQ2B8Q2jQT+9DJgPRq872fTOqIbEhWBh3cTO4prWJ6xVivjoHJz7NOHeu9eksaJY8dWNoNK41RDUO0T0PCs7Yyo3nJm2osxzqupVOyNGOtIYxbAqTVWI91ztOQWdq0mtbV+l5aVjmNYEaBQGDhiG1aIBAYBWIyMiExsAV3p/Lq7NyGVJYM0gA7rotAC3xLLgk6jrBHl3MK7GXN0MclIWqfzB4easLOljp9oigOidddctH0YlCXtmLeFq4Uz9ozbPVcPYD69qyPal9vy6ysHjWn4JLLipRWvPaYjAKBwIgQXvsuZCYXd5DvZteOTqWOW11a+TzXjpJqv2a46GWF0M+2+taCGruyMCS/WJ/wHJ7wdN5xp73wLfKcFuF+jSGV3g/PoNViKjWVviWM1vd6wu6Sar+FabYJsEObFggERoHYppmQME8ygz+SyruzOjKHb0xlnxWiZsyWq/Zb3T9KrErqlljUPFBzRrVYzjp1znPm7RObeghKzNJjvSUWrOGF0bCet1fV0YxoiLOrZj3WNc2iSjLM2vuefy/9jm427xzfZBTZQQKBgxJzzQ7yUSK6n4huys79DyL6LhF9i4j+gog21doZzIyI6AQAHwewFd0328HMHyaizQA+CWAbgNsBnMXMP25pU1aMR1P5QHZNZnkv04egpGlpCRGhV1B9rNmPdW8tAFcO3f8hcpia3MdanfU177jUfk3Tk8NjNZYWSNrVQcFK41gLGWsFGtPstxbS2DqnNWAt8p8am+oTHkSzOsB30J7G3JjRxwD8Mbr5QPAFABcw8x4iej+ACwC8u9TILMxoD4B3MvPJAE4F8FYiOhnA+QCuZeaTAFybjgOBwKggAuyWvzKY+UuYBN2Qc59nZpntvoIuq2wRg5kRM98D4J70+TEiugXAcQDOBPCqVO0SANejMiMC3dBomVG+UomdkbffFrQwiZLWq5YXy5M/leqstKyoxHL0NW+VtsKcDGFEAs+eqhQ2VSCrvtRt0arVfr+8nVbnV6CuASsFQfPGscXuzYOWjQHT7D1nxhOsqMzot9DtloqYiwCbiLYBeAmAGwBsTRMV0Gnptzr3bAewHWhLYx0IBOaJXpPRMUS0Mzvewcw7Wm4kot9LD7q0VnfmyYiIDgPwaQBvZ+ZHiejZa8zMRMTWfenL7ACAjUS8FtMhZXOC+Ji61hKyYYjWx2NNJSfN1rZa0MfyumYdbtm6eGly+oxjq/Mm0I9ZCmRFrzno9mGlJafXVufX/JwXDtazJcqv1X63PuNpBfnT52wLbKDHZLSbmU9prSwgon8L4LUATmNmcx7IMdNkRERr0U1ElzLzlen0fUR0LDPfQ0THArh/lmcEAoHlwPJu04jodADvAvDPmfmJWn1gNm0aAbgYwC3M/IHs0tUAzgHwvlReVWtrFbrZW1YTYUj5bK/Dzcq3E8Yks3++d/a2fy1BwbRsYwhkVfS0QoCfr75FtqLh2Q7lq7NmRDWGlN/fhxEJWrRaHjRL1Cyn5XktWidPE1YKdtYaOrZ0bx/Zm/4eLbZy5eD985uMiOgydHLiY4hoF4AL0WnPDgHwhbRb+goz/3apnVmY0SsBvAXAt4noG+nce9BNQlcQ0bkAfgTgrBmeEQgElgXzcwdh5rON0xf3bWcWbdrfossyZOG0oe0GAoGVQjjKTmE1gE2YppMlAbYXVzoX1pXU1cieV1K91wTa1pZLG0SWHD09wXifONba8E5vHazc7fqaPj8kDrMFT5BcasPbovZR6evnWSrwVqdXayy8MS+9c56gukVgLagJ6IHpaKURAzsQCIwYMRmZWA3gSEwLoR/O6ggzelIdi8OLCL03YALtOtIiMKy5LbQYPa516lgroVbDtrAOT0WsV2XNevJzXlkSYPcRvHrjpPucwwtyphln6V7vvHZeBurC+5ax0ELi0rtWCtNi9TmH57Bdyt/n7Qg6xGQUCARGg5iMpiAyI5nJH1MlANyXSnGAeW4qNTOyZEZ9HBBrbgSlbLRedgmdYcTqh977a3V2Dr3C6tVaM6J8TDaqc8Ik9b2WoV8f2ZunipY2WvQ4mtW0sEV9r2ZXLYHSWowea/dYrKoU2teD9x62hLIps+wIrhYIBEaB2KaZIHQriMzowoJuzurICvPCVAojkpVea4Pye1qDmufwQomUcmpJnafU8Sp13np2KTxG3gbgr+ieHEjGCJgwIc2QSmzKM3oshW+pGeeVNEl69dc55eTYGk8P0p8Sc5aypFnUcqYh2rQaQy8lbKhlR87PVXcEHKr9QCAwBsw7nMSMGMVkxOhWOQmm9tVUviyrI67/OhxCabX2VvQSakH8S3v1FtmQ9zzvuQIrqL4nK9LsJ3dArjEjSwM3ZBz199D2YCX7Io8NeGFarfY8pmmFRvHkP5oFlerWgtZZffT6bJ3TzNLTquWfi7+XxHoeEUYxGQUCgRWGTlY4AoxiMtqLzqboH9Lxyal8X1bnm6m8PZWiadPjWQoK1rKyy2KhnVzlWLOffK/eyhwsZ17vvbDa9BiRLjeoY2DCiGoMqcUCW6P0vazxAuwEClqLJn3SspC8j6Wg9vlxH9lbi1W/p7EtycRqlv+lYH+1VFr5Zz90CIIZBQKBESFkRoFAYOEIZmRjD7ptmuRJe1EqT8jqSA6UUuQ+oG18+xjP6e1aS4hc6aPQ5JJTo+faoWG5FbSq9i0XmVa3kPxzzdhxr3FNQ2/FcoWE3oZplX5LfrEaSgJsb3tWcjT2BNZ9oouWzDs8t4+SO4jcX+1DTEaBQGDhYMQ2zcI+dAJpmd2FIf33rI5cE+dZifiow3RYRo+eQ+SQEA6lSJDaQVav/lbu9rIzo31e97uWdcIyC2hlZCXM410u5WfzMr+UfjdPXV4yB6iZSJQMJT1H2RJK2Wn0dW8noE1LStFDTTDawm2uIEYxGQUCgQUgmNE0RJYmq8ztRh0vLIe+bsV79gKKWTGC+8R31vW0M2bNVcCqq5nXPLb1pZVWu1RY8h8ts9EY4k7zpCrza16G1xZ442X9rt7v4snR8s815+GWvunvZwWAq4UMKWU9WaWOlyAE2IFAYDQIZuRDVh1xC9lnXPNgGah5zp8lozatFdHaM0+jZN3jyY5yNuJpiIQxWKullzNLjuXellAl2k2j5OBZg2XA6IXiLbnT6PIJdWy5QNTyppUMJfto0/oyI4t86N/AcwnKP9dKi/0W55o5MiMi+ii6/Gj3M/ML07nN6LLIbkO32TmLmX9cameWTDyBQGB/hUxGLX91fAzA6erc+QCuZeaTAFybjosYBTMidCuLDhP7QFanZutSkhlprYgu85VIr4r71HmBlu3k5/Q9WuZihTH1QrtaISI87ZnHevIMeq1uDENQ0ui0OATXNEaPq/P599LsSdinDgWcvx9eELUSc25lRpa21ctOrO+xcqB5LFgzQAvmtTn6pjHzl1KK+xxnosulBgCXALgewLtL7YxiMgoEAgvA8gqwtzLzPenzvZgE3nAxislIgqvJqiOB0x7P6pTsKgB7FdOrR2mfLagF8fdsX6y+6FXZYmJyTvq6QdXJx0BDh+fQ2XWtVbsmAyuhtqJb8FbsUriMWkAxSxOn2ZKMm3YEtgKl1ULGWtbocq/HKC35R80q3dIiejKiktW2HltXm9YuwD6GiHZmxzuYeUfrzczMRMS1eqOYjAKBwALQzox2M/MpPVu/j4iOZeZ7iOhYAPfXbggBdiBwMEKYUcvfMFwN4Jz0+RwAV9VuGAUzInS0VyISWoJeTc092mo5oWo3DN1+LjD01PIt2xs5p51qdVv5VlKf07+9tJFv1+5V1/TWch6uCaU4zN69JbRkxqitjC056/RWX8dtyp2GxeSjFqPIcqfx4ha1jIn3PbQZRP7Zi/TYkg3Z3abNyR2EiC5DJ6w+hoh2AbgQXTiyK4joXAA/AnBWrZ1RTEaBQGABmJPRIzOf7Vw6rU87o5iMVqNbuUrGZo+k0lMNW8LGIQ6X3gpeyhSqodvwMqRa0OplYYuWIZyn1i1FG9SsULusyPOv56q8sQm/QrSkj54xJOB/rxb1tcAz9bBie2u25EW0LMWxFgxhmPr7WuFANOPX75A1Nk1KiXAHCQQCo0FMRtNYDeDI7FjUs/mKpGMny3FL3rTWUCL5M7W8oI8cRtfR8qaSC4nFGPLzVjueQ65mBcD0OL2fb0uftFTuK9lde1S5Vx0Lpl+nz/HfqGtdb/4rvQTAUtcY7c6i88+1OId60GOUf/beDws6cJkOYeM5BgO+vLMlF98Qx+kio4x4RoFAYDQIZjSNVehkI1rG8VhW5xl17nBVWquaJw8pyaZKubn6oua6op9daqNkeFcLN2tllD2dL0qfLkmldrawTAo1QyphjSqX9u49fAEA4Cv036aeog0XtauH7ings6gWeZPHMDX7yT97iigt48nr6ZH1GF9JzuTNH73f00hVFAgERoEDUYBNRKsB7ARwFzO/lohOBHA5gKMB3AjgLcxctGhYhaWrt4SWzZmROM1Kvc3quCT/qWmOrOypfUJC1GCxG++aJ+/J7WO0FkizxMNVveP4ddndopf8eCp/4pQZ79iT1ueap2cO/cVIvtFhqey+2al8cnZTd+1p+nsAk99fysfVcf5+aBYlx8KYNBsBhmnrPDLhOQKXQqTovg2RhZVQ/blGJjOahwX2eQBuyY7fD+CDzPyzAH4M4Nw5PCMQCMwT8w0hMhfMxIyI6HgAvwrgDwH8LhERgFcD+LVU5RIA/xnAR1rakxXivlTenF0TbduWVLZYHXvXSgH59ezs2QhZ6WA8GZEO02GFpPDkPJrtABPrYillbA49RF2494iu3POZyc01gYwlMip5cuYopmtNnGT9Q6lM59ffPbknffl1qU9HvzyVQoESZeZHlhwCmLAk4X2PqvNS5mFHPKt+LcuxAtt5gdG8MC55u5oRadbWJ1Da4JAvB+A27UMA3oXJ/8rRAB5mZpFw7gJw3IzPCAQC88aBJMAmIgkzeSMRvWrA/dsBbAe6mexJTHyuvprK3E1YGJFM5lp70mKH02JF3bqNLuVL97R3Je2WzOaa9Ui5GRMcsVadlPLoVF6RylsSP8iFK0IN5Jynoupj8i3ok8jeUmVqOnhpKv8wlYkKUSJXRz04ufWoRIk4XXto6S3Plo9MbpliT5o0CnJbqNbA+C12RjW7I6CdvFjvclUWNjKZ0SzM6JUAXk9EZ6B7pY4A8GEAm4hoTWJHxwO4y7o5xUPZAQDPaYh1EggE5ogRbtMGC7CZ+QJmPp6ZtwF4M4AvMvOvA7gOwBtTtabQAYFAYAE4kATYDt4N4HIi+gMAXwdwce2Gveho9K3p+GWpfF9W55uqFLplqWwFNeGetdXSgmkv/EgpBIbemeit2BHZPSJ8lh2W3nkdJY1swQRbnFJukoG0bCS0RFdv0yzJa81HxVrSanYU1jZN2yTI8YtTKd9H9mB5kPR0jtK5o1O5+aGlVfMtsnaN8RycLdGKt9XyjvNzta2e9T5qtEQcLeJAdQdh5uvRBdwGM/8QwMvn0W4gEFhGjGybNgoL7L3oFmpR8soe7+ez5CZPJJok6n5vJcqhM2p6AuwW9bxnuGi5n2jBtazGwohyYbR81iTnULENfI4q80pb1bEwiB+mUqnEAUxbDmofDMuPwQu8rGHZSOiBKwmwPbsGGSQZgwfU+fycsoWgVP5M0o6sz6TRtcweFrsReCp8z7AR8F1VSsHVSqYkOfLzTU61B5I2LRAI7McYoQB7FJPRPnQriaxUX0vlZzKhkSz2ov7X4hBZTUpGj16OMEsj7d2r61nnNCOSBV/kQ5swwRQjkkrCAraqY+va81Mp1qJan22p9j2jR4sOeNHOSvBU+yUvZT1gj6tSBu55ql7ejpf8LP2QR9w3uUW+u8eEdLbdrJln4d2riWZ+rZQ7DmiTGQlKZirVn+lAlBkFAoH9DMGMbEhGWREBfDmV+eIsC55e7GXRbDF6fFIdl5iRFcZW19XQi77nyJqLOkSL9qyMSMuDLJmRlhXVLPxKHqWesMOyvBsS/1XTUTmW5+TMqFXNJLKxkrBPUPBcPSL14cnUrnbTkKGy5IICLTvycrzln4cQTY2S9qzJXmeEk1GkKgoEDlbMMVUREb2DiL5DRDcR0WVEVAvVNYXRMKO1mLAEWehvyOqIuECnhCkt1rVA/JY2TTuz7lN1tOyoFPRMa9G0iwcAbFyrTm52KudxefU1rT3zWI91ruYtmp/ro03T9+q4uqUfTrevfyj53lbqVU/2ZY1FGqfD0wsnw6eVejmx1MH3dJdL7iD6Wh+CKahpe/O+FNufozaNiI4D8DYAJzPzT4noCnSG0B/r084oJqNAILDCmP82bQ2AQ4noGXTz+d2V+mYDC4fIjESm8txU3pnVkW/mpZWxoMdaLwSWzEgzH8Eh6thakdaqulKWEglOCZSOcM7nMUSkAVEteoKKEstpjYWRf64ZsOTn9Q+jB3+VcV4L9eRYy5k+mcozs3u1gO4JdazjFWefD03XNqa+aEWfxX5rKataEjJ6KaYs9JGnNFtjz2kyYua7iOgiAHcA+CmAzzPz5/u2EzKjQOBgRL/01scQ0c7sb3veFBEdhW5pOBEdl9hIRL/Rt0ujYEaBQGABaGdGu5n5lML11wC4jZkfAAAiuhLALwL4RJ/ujGYyWoXpLKpHZ9eF2gr7XrLVgU1NaypUi2pr38916rxGTuG1vd1GVUqfD80fqPdwniNprpuQB3nS0hbpfstWaznRknJV73X097TSrOjUsdrOwgqzmepuUNs07eML1B2mW3a7NT3AkC2ZJcAuYr7uIHcAOJWINqDbpp2GLi5+L4xmMgoEAiuIOQqwmfkGIvoUOueJPeiidezo284oJiMRYAusjBhyTuSQfSzZPYJgRcDQubJqWmxLgK0dZbXK30zU5jEhS1Iv52o63Fms6fI2PY/j5WJR3uBramFRltqPYEmj15WrWkaPmhG1ZJPRrKnP8C2LcHeO7iDMfCGAC2dpYxSTUSAQWGGM0AJ7VJORF6sa8FX5um7L+Ho2daU6uq7VR0+eMCVfsGwJ+gTsLgXgbkVNN126Z6VfYv09refr7+GNp8WMVtlVS9ljPIa0XNDGt3NBTEaBQGDhOFAjPc4DNZ/HeZAAb3UphQOZx0o01cZKLacrjQP1ex2IYNixmheI0UxGgUBghRHMyIdsYWXCtvw7BVbeeqC8OGsllCU20Y6QumyB5z3x7BbdCs9Ru8kyWJlFxbIc8oKSO8gsaIk+76lKBQ2GP3pIhmi9SqzbUlQuDCHADgQCo0DIjMoQwiC2RA9l1yRemFhntzjKeqtUydlRm6eUNCtA2Y5Evs9T6th0Dq2lKLUeVNOItWjMWl7IWaKA1WCpqjzzZn09l3l4MTxK4XTV2Hux3EpG4jU7NEsJ6oUfGUImZ/pJghkFAoGFI7ZpgUBgFIhURT72YhJzWGIXPZhdPyGVXsDA0iTvsX7LDk5HeqxtB/Pn1oINyvc7PHsJ1ul0Ejp7h+xLre2F5/rQYpfQx+jR21d4bZau6ef2Sc2iv6e13dUJy0qRHlVdXcWK1thqDFvKydcnjtEQ29RmwhPMKBAILBwhwLYhjFGSgorg+p9mdYSpPJpKSXzRklHWW10sZqTPOem3TE2xXkmF3IjDryY9ALBOpPU6zsjhqnIeiFnqSsYQaVh3WjqSO5RqKa1mI0Nyo1nwMsrqYOKlIOJeOJA3pTLXcMgY6LQxjzllVufpZ5Y2UUqU4uU6E5QyEHtuSEOcvltQZVPBjAKBwMIRzMiGCPaFGf1SKv9dVkcSYHw7lberNqxVwPORbMkO4uVPK6n2vYQUetFeEt8rySsOfVhdLEX40teERXkGk5ZZgIdSxDntPlAyvvQGv5RR1ks4p+OCy1jlQsUH1TUvl9zDmCDRa02m5PeaMsmAH7XF+7rWcEsdj2XNSlia5phwBwkEAqNBMKNpSHA1YRKvSOUZ/2VS5xu/35WSEEMWPh2B1PKDK77VAAAcEklEQVRIqDGiIczIiqYhv61mRF6cr/zzlrRKr2vppMbzK/dYKh2dLlVYiVYlAW25zjSkL57Gz2JGOgSvlyHlvlQKlc4/S3mvcz6TMz2aXh5Nmh5XpTUUAk/GWGJE3jtU0rTPmz2FnVEgEBgHYjKysQrdonhcOv5yKr/3+5M6t6ZSMyFZiSymMg9m5KWMt8QksoLK87yc7SWSsyWt3OsKOeJdoZTkb9CsQ5gGMK3a022UUqHWmJFlM6QH0NOQ5f3UMqM8BS8wERjm2jTNgByG9OhT07dIM1rhZmnTPI3YVGjhhPy31mnhrNRx+XnrmsZM80ls0wKBwMJxoDEjItoE4E8BvBDd1/stdCTmkwC2oVvDzmLmH5faWY1u8RPRx1+l8p1ZnS2ptAxpgX4GvR5jsup4dUuiHPmN9UpYYkSaCG1OMqQjtGU2MFm6tbBDNEpvUvVy2xox1NLtloxrdHJ4raWzvpgeME9WlGdd8JjRLamU7/eQOs7PKbqz9yH7VmB6+MR2TZhsyabIY0IC633UjFkfr1bHwDR56ePBUXS8nbM7iDUXMPPf9Wlj1qQDHwbwl8z8cwBejO61OR/Atcx8EoBr03EgEBgb9jb+tcGaC3ph8GREREcC+GcALgYAZn6amR9Gl+b2klTtEgBvGPqMQCCwTOiX3rqIwlzQC7Ns005ER4r/jIheDOBGAOcB2MrM96Q69wLYWmuI0NHd56ZjkcPekNU5MpWSZVbosecEa0HPvFZM7Nr2TLuFWJB2rRA6wFJVsXYdkV2T/JKbU4VNd0/uOUK2IptTuUkdvzeVLds0/WCr016mWg1rIGvGji3btEtTqfdV+Z4r7cP2/mTJ4dQWrM9QyBDk/4+lZCP5eSvp7VOqjn4vBKUYWVYdD9X/ifnJjMy5gJkfL9+2FLNs09YAeCmAjzDzS9D9lku2ZMzM6ObgKRDRdiLaSUQ7n7AqBAKB5YMIsNu2acfI/2r6265aq84FLZiFGe0CsIuZhcB8KnXgPiI6lpnvIaJjAdxv3czMO5BS4B5LxMBksRSG9IKs/vdTKSubMCRZRC1htGean3tWAHY0Ps8coMUW0VtwhBHlz5fVUjTdOqLlA+oYmIQg2XTf0nKjdPrn5eZElTijEHr5l2MdesOib1687uUSYL8wlUpgv/eJJYdLPnt+sfpr5p/1Vy4Jrr0QM/oei2CWoqfk91ouRkNU/Hr3MIV21f5uZj6lcN2bC3phMDNi5nsB3ElEMmecBuBmAFcDOCedOwfAVUOfEQgElgmiTWv5qzXlzwW9MKud0b8HcCkRrUPny/qb6Ca4K4joXAA/AnBWrZF96FYuzTZyWzdR7d+ZSpEJyGqmjR8Bf3EuRcnwAmB5Kn3LlEDg7fczu7tnV2xPxCIsKicQmlSIH+mG9OIcvjuV1DGif8T/anLzYYlDHCZPFk6RhC2lKGR71KiUooTJFyJt7SjlYaoE5Bd/nL4EoB4NJBdI1IieHFusA8Y1wCZ8Nf9f7VNshafR7WoZY94vfY/u8+BEMfO3M7Lmgl6YaTJi5m9gIm/Ocdos7QYCgRXAHCejwlzQjFFYYO/DUi2H9kwAJiuMiBVkPdfOqbk8Rq9aWnYkbVlsqvY7tWjxainR8mv6nlKYE82atGeF9gJZT3/x7L3Cpl7HH0qfhBE9pY5zZrSnK9Zo/c8edbzG+Cyl7l1X/g1NvKG1+4znsVKKIOsp/izxSIsRrD5uDUOsDRpLdaXP8jta772c08H9BC3fbwkinlEgEBgNDiR3kHlhLzpmpG1ucrYk17Q2Q+paYRhqzouWLYjHVFqcXPU9Wn5Q8nn1tDF9VnZPjpF/PxmLv6a3AwAuYpHCab6Rsx7rnHWcw2NG3fF76UVLWgbqCQ28McrPaXhjAkz77nrO16WwxJ6c0JIp6jqeL3TO7j2XopZwIyvpDjIPjGIyCgQCK4+REaNxTEb70MkDNCPKjSG186KGdV5rxoZskUsrXA2ezUnOBkRSo7VAfYLC176XlTFXzn2RuiRQesX/Mpu2qr3xaiIAvtzMYjeljK4e9O/kmTVZEUu0DFFbVeeoJXkovSea+WhmXmJvXh4FQekdsPo0Qqf9cUxGgUBg5TEy+fU4JqN96BiAjoyRMyNtJ+IFpLJCuraW1rkhyfO8KK3SZ8vOSL5rzV8KmA5B4SWctGymav5QUp6QGE1+f4sGx4OXJqoPSvHm5Jpn4C0W7HmcOd03/dtbtkLed9bB/axgf/N4H7W2t8RuSl4CwYwCgcBoEMwoEAgsHPswukxF45iMRMuo84vlqn1vC+K5TwDTVF0fb1Tn8886mYUW/PaJ2qhV+/lL4Bn66THIqbu4xmi/Ui/lWskloWX75NH5lpW1Fvpin3HN297KuGkjSOucbPVLWxEvuklJgC6fLadn7x4N/RuUssfoPtUE2dazvXc1mFEgEFg4QmZUwF5Mr25bsuterGGBZeDoMSPNiKyVqJZJ1oInuC45Qnpqf51RJA8hskmd08JZjynln73vV0oO68F6qb2xKKn29fh4Zg6WUay8M1pwLOet98N7nueCkbfvQT8/Zx+e87VOoJKzrZr63xp7bSDsISajQCCwcIzQNW1ck5GsTMKInpdd00Gs9ECW9ts1mYrlXFszdrTSz3ss4Gl1bLEBjz21pKT3ErDq6/k57WSr5RZWxhQY1wCbBXjyHy0/KwU78+LACSzWIW3oBLklVuqVuQmGhpbleLLElhCynsFm/tmTHVlt6j54Ro8j8wYZ12QUCARWDrFNM0DoVgCZqSVc3C9mdWSF0waR2sjMYkY1LVOJBQhatCSeUZ63muWfNXsStGS9nSXmfQtbbE16YIVG8ZyELdcYjwHp9rVsJ++jNlbVhp95H2t9K31fac+TKa5S9fJn1xhSyR1Efiedyq7EcsLoMRAIjBohMzKwCt2KLavJtlRekNW5JpXfTaWEmNcrQmm/7TGiPgHS9Pn8+ZoNePZFJZmRlT5pKGZhVaWwGR4sFqBlfdodpeQu4bHUebiS5J+1/ZIXKN9qx3uXrL7X/vmt38vTosnvpBma/uwhmFEgEBgFYjJysBpdkkZJUXR7Ku/M6mjZkJfEsbSq9GFCgloKmlKgtKfUccnBs8XxUaDlH7qPpWBuNVsoK9yKN16WM+g8oPvgjVHJvkmgw7SW0gB5zCiH3O9pHzUrGRJ6xkqdpWVH3jFgJ5/UGKM2bZYkjoFAYD/GnLJbAwCIaDURfZ2IrqnXtjEKZhQIBFYWy7BNOw/ALZhkzuqNUUxGq9G5OJyUjv86lf8xq7M1lV422BbVe5/IgVoYrY+txKteri5dlpK1en0rCVG1AHadOrZU4H1cBXTfvDjMlmpf91H3tSRQrrnVlFwtNKztpzYVkO+ljR3ze7TgWDvM6uMW95qSY6u0owXVLdEha+/5vCYjIjoewK8C+EMAvzu0nVFMRoFAYGUxZ3eQDwF4F5a6UPbGKCajVehWnOek45el8mtZHVmRhCGJs6gY8encYcD06tSi8hToVbgWFgTws1roMl959T1ehtAWwatmRC0syGMdh2R1Sq4iXh89A1BvTPLPXl3NmKznaWgWZzExqeNFwSyZLHjMxfoddV0vfvZe41yrIDvvwxyZ0TFEtDM73sHMOwCAiF4L4H5mvpGIXtXe5DRGMRkFAoGVRU9t2m5m9rLFvhLA64noDHRc4Agi+gQz/0bfPo1iMpKBkRld2M/JWZ0fpvLmVAqL2pxK7eYATK/ofQKJ1UJgWIHSPBbgnbeuSfu6r9bqXFNNt7hxeLIdi2H2iQvutav7agVI88bNMxAtoSSf0czIM1hsYRCarVluNTIGNUdW6/fSrEp+H4slennmcsxLgM3MFyDZJydm9B+GTETASCajQCCw8gh3EAN70QXK0jKdXBq2NasLTAwin1R183AZ2rVhSIjVVoZknfNYQSmMRalvXp/0c0qMqBbSVYeqAPq5PNT6WHKU9eRJ3tiXoL+7pQHU32NIJhhtBFlCrU5J86YZl36u9XuVYlwvhwU2M18P4Pqh949iMgoEAiuPcAcxsBddaBC9ClihJHToWB143Voh+riDWH2zjq1V2tO8eezAqtMnzK121mzVDuXwwsDOmxnVGFL+2ZMVlbRpGto9Q9sU5e21hkjpAytnXcmeKEc+rtpxumZvBEyzJWs7NkZ3kFFMRoFAYGURjrIOhBl5KyEwHXjLCy6Vr3yaxQxxxPNCh/QJteoxpPycxzqsPnvhOXQbFhOUcfScUUsyCI+1lWRTNRlVSX7mMaCWvPKeRjHvqx63WpslaFbfkpPMa9dylNXQ70Cu/dQO2l5fQoAdCAQWjjEyo5m89onoHUT0HSK6iYguI6L1RHQiEd1ARD8gok8SUYuiIRAIrCDEHWReXvvzwGBmRETHAXgbgJOZ+adEdAWANwM4A8AHmflyIvqfAM4F8JFSW/vQba8eScc6KwQwLYzVavtShMJZtmd9BNi17VlJ8Dqkj9pAUm8Dh6iqtVAc8LPp9olqWBPqA74Ae4iRo5d5uOSyotvwjnPId/faKLl29FGo6DHW32efca0kwNZ9GwNmjWe0BsChRLQGnfHzPQBeDeBT6folAN4w4zMCgcCcIdq0lr+VwmBmxMx3EdFFAO4A8FMAnwdwI4CHmXlPqrYLwHHVttB9aWFED6RyU1bnyFR6rgmleM8a3qqdf9ZC73msLlYb3io8RM1cMsisoSSY14LePiu7167Vx5oKv0WV7wnztcDX6kON5ZSeK5A2Stlo9RiUTAtqUTa1mwgwzZqs9+CAkhkR0VEAzgRwIrqIsRsBnN7j/u1EtJOIdu6pVw8EAnPG3sa/lcIs2rTXALiNmR8AACK6Ep0H7yYiWpPY0fEA7rJuTiEIdgDARiIGpl07tmb1hSVZsiFgmGuHtXLUgqq1qoOBaVmOFSJiFnjqeS9nGDCt5vXYR4vRo6CP42yJgekxr7l/WGp6LZcpyVa0ZsVjkqXvV5MvWe+JJztqYX4t4Vt0LjxLtT/G9NazyIzuAHAqEW0gIgJwGjqn+usAvDHVOQfAVbN1MRAILAcOGGbEzDcQ0afQxUDbA+Dr6JjO/wVwORH9QTp3ca0tQjfjy8q0RZX5Z3EDaQmUVsuiYbk+eKEutIbCkgl4zy+t8JrNtEAbSnrnrRXRc6IsadM065hFrtViNOqx0xZtpGaDVl4xr09aC6rbts558qaSY66u0yI70vDcQ/LPUlpuQWNkRjMZPTLzhQAuVKd/CODls7QbCASWF4w2K/GVxKgssGXFkIBpJ2TXhBnJPrgW1Dxvz2NClrm8Dg2rVxcdrD2Hp6UryXDmQYNrMinrpdNaQt3HZ4y6fRxla7+PtSrXNG99VnJPLlMKIatRkg/W5GYlNuW1YUHLtbzn5vW07Z0nCzugmFEgENg/MUbV/mgmo1WYDI6ElM3DzoqGzVsttXwhr+MxoZYMogJP/mNZ2Aq8lD5WnaHXc7SsdDr0qWZGFhvwZBlDtGqCUl9LciWrXqlPfay3PQwx/Ctp0WrH+b01mWIpRG1NZhSTUSAQGAVim2aA0M3mMsuLrOhNWZ3bU3l3Kh9LpU7AVwrpoVdaa9WsMaCSVbB3b4lN9Qkdou/RaFnpPM1RSZ7RyoDmGZwsR58QIhqlNFXSbs2TuxSkro8f2xBmWWOl1j1aZvSwUSeCqwUCgVFgjNu0WR1lA4HAfop5GT0S0QlEdB0R3ZxCCp03pD+jYEar0Ln8a1Xkcc+d1Lk37c+Eouu8W1budi8rh66bq76fUnVqWU6tEBhe+5ZwuBanu0+8aU+1n29rvDp6G2BR+JaYza33tF63UPoHqa2uVoZXT6hvGRT2DTdS2u62OEnr562Hjfx76XfKMkeZs9HjHgDvZOavEdHhAG4koi8w8821G3OMYjIKBAIrj3lt05j5HnThg8DMjxHRLeiidex/k9FqdG4ez0vHX0jlf7p7UkdU+7JSSLgRj7nk54bkMdPq/xaDSa+ONq7LBaZawOoxonkxCO8FbAmbURN4DnEotdpcTtlBnzAdJXjxpfsYNrYwI89os8SypI5mfjmWS2ZERNsAvATADX3vHcVkFAgEVhY9tWnHENHO7HhHirqxBER0GIBPA3g7Mz/at0+jmYxWYeLyIarIT2TXt6VSGJIXGqIPUyk5ynrGlHpFtEwJNKtZr45LK18LI/JW23nIAIYEZhO0hGct1W01rvTa7NsnD33GseZU2+IWos9b2UG083PJBETLjObgDrKbmU8pVSCitegmokuZ+cr2picYzWQUCARWDvPcpqUQQhcDuIWZPzC0nVFMRvvQacVk9hd3kDuzOl9LpTjRSigRvTKUgnVpWVFLVlhBS8hVjwG1hN7oywqGwnsBZzEsFJTYlMeMLJboORzPYyxKbbS4nXgynCFB9/Sx9Q542lWt8bMMN61ssznmKDN6JYC3APg2EX0jnXsPM3+2TyOjmIwCgcDKYp6qfWb+W3SOFDNhFJPRHgAPYlrDsjGrI+ckWL/IlUqm/F5wrpJcRrMbGHWs61adWe6dJyNq0aq1BIDTaHmZvXG0nme59ORttGQGnmXcdLut7iI5hjhFe+zHOlcLq1trL8fYLLBHMRkFAoGVxT6Eb5qJPegYTymNjXZW9BLUlcJK6MBsJe2WvtfDEBlB3/uHosRyvHGz7vFY0xBmVGJgMhba5sqDZVk+T5Q0mZ5TdB9m1Oe314xI2JoVJljfswIyo7lgFJNRIBBYWRxwMbADgcD+i2BGBvYAeAiTrZiVG03nghJ4MYPyz5o6a0G2Rfe9SIgaluBwFkgbNUfWWl9azufXSuOIhmveeS/yojWeeqX2FA8tbQzZtrUIelt/41mMRi0Bto7bJcqZoQLsMYYQGcVkFAgEVhYRXM0Bo5v5xR1EZ48FptW7tfjW+Tmvrpy3Ivl5oTtaDBd13SEsp9ZGC1rYTUsmDt1On2wdNQPG0lh4jMgygpyFEennabT0sZZhZF7wWM9ao46ghdWPAaOYjAKBwMoiBNgOCN3MLk6w4vJxeFZHG555YUBagqt5oUSAaRbluZK0sCnPYK3kXNuCViO2WXKTWbI33e4sK2uL0aPHhi0jxD5ypVlQU8+XMtd6LKnPOOp3y3KG1X0Io8dAIDBaBDMqYDUmM7q4gWzKrmuDxeViRq2B2CyHWn1tlh+7RU7hoWQI2sqIrKwnrc8vYR6rsQ4Tm38uyZWs4xbk9+j+awNNreXKmZl2dyplHfHgMaOcDVkaNgvBjAKBwMIR2jQHhG6m124buTZtg7pHz/4WM6pp3Kyg+nrV0lk5tfwnD3buuVSU5EKzhJXVaHHx6MOINGor6RA21/K8mizJqqu/n5WIoNWxuVRH2tNuGvp8CfJutTCVFmakE1tYCDujQCAwCsRkVMAqTGb3Q1KZsyFPm6YHNF+JWlfyfIXVQbI8621LG+SFvlinji2bkFbbEAs1RmQxiJpGrGRn5GFWa/S+/xyW5bzHhKxxnofGbRYZlCf/HNJG/k4J0+rzG48Bo5mMAoHAyiGYUSAQGA2CGRkQAbZW6R+d1ZFZXDLI6syyLTnQWnKftWaQtYTfWgBfyxYL9I+fBLTHJGrZStbaHIp5qP+9yI8lhYDeeuk6Q7aSsxikWqp9rf7XaNmu6XfpkOyafmet35Rh51NbJKrbZSL6KBHdT0Q3Zec2E9EXiOj7qTwqnSci+iMi+gERfYuIXrqcnQ8EAsMgRo8tfyuFFmb0MQB/DODj2bnzAVzLzO8jovPT8bsB/EsAJ6W/VwD4SCqryAXYz03li7Lrklz2zlRKDGxhSjrDLDBRu8/CjDQT0gLD3PxAMyAvYuG8hbslgXXrPSXM4qw7BN74lNT0NZTY1DxdcvT1vI9aoVFjSIDPkrSJSd6GtG+x9xxjkxlVmREzfwlduKEcZwK4JH2+BMAbsvMf5w5fAbCJiI6dV2cDgcB8IALslr+VwlCZ0VZmvid9vhfA1vT5OCxNd7YrnbsHBTA6FiMzuciKXnfipM5nbutKaVyY0IOpfCyVT0xuebaOJ+exXDo8Y7INznXLJcGTV7QEFGtxdm3FEDlQqY/LGa+7DyyVfM2osg8jGsKUBH1YpPd+WI7A2uSk5Eoiz5b/Ca/uASfAZmYmIu57HxFtB7AdWJ5g6oFAwMeB5A5yHxEdy8z3pG3Y/en8XQBOyOodn85NgZl3ANgBAET0wAPA49cBuwHgulTnt28b2LvlxzFIfd0PEH1dHiykr0/Uq1iQvj5PTuwDPvd4d74FK/I9h05GVwM4B8D7UnlVdv53iOhydILrR7LtnAtm3kJEO5n5lIH9WVFEX5cH0dflgdVXZj59Uf3xUJ2MiOgyAK8CcAwR7QJwIbpJ6AoiOhfAjwCclap/FsAZAH6AbhL/zWXocyAQOABRnYyY+Wzn0mlGXQbw1lk7FQgEDj6MSXa8Y9Ed6IHo6/Ig+ro82C/6Sh2ZCQQCgcViTMwoEAgcxBjFZEREpxPRrcmn7fxF9ycHEZ1ARNcR0c1E9B0iOi+dN/3zFg0iWk1EXyeia9LxiUR0QxrbTxKRZVO3EBDRJiL6FBF9l4huIaJfGPG4viP9/jcR0WVEtH4sY3ug+I8ufDIiotUA/gSdX9vJAM4mopMX26sl2APgncx8MoBTAbw19U/8804CcG06HgPOA3BLdvx+AB9k5p8F8GMA5y6kVzY+DOAvmfnnALwYXb9HN65EdByAtwE4hZlfiM4Q+s0Yz9h+DIBW1XvjmPuPbkfnPzoOMPNC/wD8AoDPZccXALhg0f0q9PcqAL8M4FYAx6ZzxwK4dQR9Ox7di/dqANegi86yG8Aaa6wX3NcjAdyGJLfMzo9xXMXNaTM6DfQ1AH5lTGMLYBuAm2rjCOB/ATjbqrfov4UzI/j+bKMDEW0D8BIAN8D3z1skPgTgXZi4HR0N4GFm3pOOxzS2JwJ4AMCfpW3lnxLRRoxwXJn5LgAXAbgDnZ/lIwBuxHjHFujvP7pwjGEy2i9ARIcB+DSAtzPzo/k17paYhaoliei1AO5n5hsX2Y8eWAPgpQA+wswvQefXvGRLNoZxBYAkbzkT3QT6XHRxAEdnwexhLONYwxgmo2Z/tkWBiNaim4guZeYr0+n7JDyK8s9bFF4J4PVEdDuAy9Ft1T6MLoyLGLeOaWx3AdjFzDek40+hm5zGNq4A8BoAtzHzA8z8DIAr0Y33WMcW8MdxtP9vY5iMvgrgpKSZWIdOMHj1gvv0LIiIAFwM4BZm/kB2SfzzgKX+eQsBM1/AzMcz8zZ0Y/hFZv51dH7Hb0zVFt5PATPfC+BOInpBOnUagJsxsnFNuAPAqUS0Ib0P0tdRjm2CN45XA/g3Sat2Khr9R1cEixZaJSHaGQC+B+AfAPzeovuj+vZL6CjutwB8I/2dgU4ecy2A7wP4KwCbF93XrM+vAnBN+vx8AH+Pzl/wfwM4ZNH9y/r5TwDsTGP7fwAcNdZxBfBeAN8FcBOAP0cXdnoUYwvgMnSyrGfQMc5zvXFEp9T4k/S/9m10GsKFjy8zhwV2IBAYB8awTQsEAoGYjAKBwDgQk1EgEBgFYjIKBAKjQExGgUBgFIjJKBAIjAIxGQUCgVEgJqNAIDAK/H9F+zKnJnzTtQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#NBVAL_SKIP\n", "\n", "# Plot percentage error\n", "plot_image(100*np.abs(model0.vp-get_true_model().vp.data)/get_true_model().vp.data, vmax=15, cmap=\"hot\")" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEaCAYAAADDgSq4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xd8VHW+//HXJwkkhBJ6CyWIgCDSawJXXXsDdW0IKFlX1CuL26+76+967929d7vrqlhXARsKCHbFsrpKQq+CgFIChKr0DiGf3x9z2I0xwEBmMpPM+/l4zIPMOWfO+cwk5JPv5/s936+5OyIiItGQFOsARESk6lKSERGRqFGSERGRqFGSERGRqFGSERGRqFGSERGRqFGSERGRqFGSEQHM7GYzm2tme81sk5m9Y2YDYh2XSGWnJCMJz8x+DDwI/B/QBGgFPAoMjmVcx5hZSqxjEDldSjKS0MwsA/gf4G53n+Lu+9z9iLu/4e4/M7NUM3vQzDYGjwfNLDV47XlmVmhmPzGzrUELKDfY19fMNptZcolrXWNmi4Ovk8zsXjNbZWbbzGyimdUP9mWZmZvZbWa2Dvh7sP0WM1sbHP//zKzAzC48hfPdambrzOxrM/tVibiSzeyXwWv3mNk8M2sZ7DvLzN43s+1mtsLMbqiAb4tUIUoykuj6A2nA1OPs/xXQD+gGdAX6APeV2N8UyAAygduAMWZWz91nAfuA75Q49mbgxeDrHwBXA+cCzYEdwJhS1z4X6AhcYmadCLWuhgLNSlzzmHDONwDoAFwA/KeZdQy2/xgYAlwO1AG+B+w3s5rA+0HMjYGbgEeDWETC4+566JGwD0K/tDefYP8q4PISzy8BCoKvzwMOACkl9m8F+gVf/wZ4Jvi6NqGk0zp4vgy4oMTrmgFHgBQgC3DgjBL7/xOYUOJ5OnAYuPAUzteixP7ZwE3B1yuAwWW89xuBT0ttewK4P9bfNz0qz0O1Xkl024CGZpbi7kVl7G8OrC3xfG2w7Z+vL/W6/UCt4OsXgXwzuwu4Fpjv7sfO1RqYambFJV57lFCf0DHrS8Xxz+fuvt/MtpXYH875Nh8nzpaEkmlprYG+ZrazxLYU4LkyjhUpk8plkuhmAIcIlZrKspHQL9tjWgXbTsrdPyeUlC7jm6UyCCWMy9y9bolHmrtvKHmKEl9vAloce2JmNYAGp3i+41kPtD3O9n+UOmctd78rjHOKAEoykuDcfRehUtQYM7vazNLNrJqZXWZmfwAmAPeZWSMzaxgc+/wpXOJF4B7g34BJJbY/DvyvmbUGCM5/otFsk4GrzCzbzKoD/wVYOc5X0t+AX5tZOwvpYmYNgDeB9mY2PPhMqplZ7xJ9OSInpSQjCc/d/0yo8/s+4CtCf8GPAl4l1K8yF1gMfAbMD7aFawKhzvi/u/vXJbb/FXgdeM/M9gAzgb4niHEpoc79lwi1avYS6v85dDrnK+UBYCLwHrAbeBqo4e57gIsJdfhvJFRu+z2QGuZ5RTB3LVomUtmYWS1gJ9DO3dfEOh6R41FLRqSSMLOrgnJeTeBPhFpWBbGNSuTElGREKo/BhMpWG4F2hIYgqxQhcU3lMhERiRq1ZEREJGqUZEREJGqq5B3/Qcfoo4Sm3fjY3V843rENGzb0rKysigpNRKRKmDdv3tfu3uhkx0U1yZjZj4DvE7pz+TMg190PnsZ5ngGuBLa6e+dS+y4ldI9AMvA3d/8doSk8Jrv7G2b2MnDcJJOVlcXcuXNPNSQRkYRmZmtPflQUy2VmlgmMBnoFiSGZ0E1dJY9pbGa1S207s4zTjQMuLeMayYRmmr0M6AQMCWaIbcG/5nk6Wr53IiIipyvafTIpQI1g0aV0vj3n07nAqyXW57gdeLj0Sdz9E2B7GefvA6x099XufpjQ3dCDgUL+Nc+T+p1ERGIkar+Ag4n5/gSsIzQNxi53f6/UMZOAacDLZjaU0DoW15/CZTL55ky1hcG2KcB3zewx4I2yXhjc2Pbkrl27TuFyIiJyKqJZLqtHqFXRhtA05TXNbFjp49z9D8BB4DFgkLvvLe+1PbS6Ya6733W8Tn8PrXw4MiMjo7yXExGR44hmKelCYI27f+XuRwi1LrJLH2RmA4HOhFYmvP8Ur7GB0FoYx7QItomISByIZpJZB/QL5loyQku+Lit5gJl1B54k1OLJBRqY2anMcDsHaGdmbYLpz28iNBOtiIjEgWj2ycwitAbGfELDl5MIJZSS0oEb3H2VuxcDt/DNVQgBMLMJhBaX6mBmhWZ2W3CNIkJTsk8jlMAmBlOii4hIHEj4uct69erluk9GROTUmNk8d+91suM0vFdERKJGSUZERKJGSUZERKJGSUZERKJGSUZERKJGSUZERKJGSUZERKJGSUZERKJGSUZERKJGSUZERKJGSUZERKImJdYBxNqKzXs4948fndZr69aoxg29W3JN90zSqyf8Ryki8i0J/5sxvXoy3VvWPa3Xfrl1L7+auoTfv7Ocm/q0Yni/1rSsnx7hCEVEKi/NwlyOWZjdnXlrdzA2v4B3l2zG3bmoUxNyc9rQt019QsvoiIhUPeHOwpzwLZnyMDN6ZdWnV1Z9Nu48wPMz1zJh9jqmLd1Cx2Z1yM3OYlC35qRVS451qCIiMaGWTITXkzl45CivLdzA2LwClm/eQ730agzp04rh/VvTLKNGxK4jIhJL4bZklGSitGiZuzNz9XbG5a/h/c+3YGZc2rkpudlZ9GxdT6U0EanUVC6LMTOjf9sG9G/bgPXb9/PczLW8NHsdby3exDmZGYzIzuLKrs1ITVEpTUSqLrVkKnD55f2Hi5gyfwPj8gtYuXUvDWtV5+Y+rRjWrzWN66RVSAwiIpGgclmYKjLJHOPu5K3cxti8Nfx9xVaSzbiiSzNyc9rQ7TSHU4uIVCSVy+KYmTGgXUMGtGtIwdf7GD+jgElzC3lt4Ua6taxLbk4Wl3VuRvUUTcggIpWbWjIxaMmUZe+hIl6ZV8i4/ALWfL2PxrVTGdavNTf3bUXDWqmxDk9E5BtULgtTvCSZY4qLnX98+RXj8gr4xxdfUT05iau6Nic3J4vOmRmxDk9EBFC5rNJKSjLO79CY8zs0ZuXWvTw7o4DJ8wp5ZX4hvbPqMSK7DZec3YSUZJXSRCT+qSUTZy2Zsuw+eISJc9bz7Iy1rNu+n2YZaQzv35qbereifs3qsQ5PRBKQymVhqgxJ5pijxc5Hy7cyNn8NeSu3kZqSxNXdMhmRk0XHZnViHZ6IJBAlmTBVpiRT0orNexiXX8DUBYUcPFJMvzPqk5vThgs7NiE5SbMJiEh0KcmEqbImmWN27j/My0EpbcPOA7SoV4Nb+rfmxl6tyEivFuvwRKSKUpIJU2VPMscUHS3mg2VbGJtXwKw126lRLZlre2QyIjuLdk1qxzo8EalilGTCVFWSTElLN+5ifH4Bry7cyOGiYga2a8iI7CzO79CYJJXSRCQClGTCVBWTzDHb9x1mwux1PDdjLZt3H6R1g3Ru7Z/Fdb1aUCdNpTQROX1KMmGqyknmmCNHi3l3yWbG5Rcwb+0OalZP5rqeLbg1O4szGtWKdXgiUgkpyYQpEZJMSYsLdzIur4A3Fm/kyFHnvA6NyM1pw8AzG6qUJiJhU5IJU6IlmWO27jnIhFnreX7WWr7ac4gzGtVkRHYW3+3RgpqpmghCRE5MSSZMiZpkjjlcVMzbn21ibN4aFhXuonZqCjf0bsmt/bNo1SA91uGJSJxSkglToieZkuav28G4vALe/mwTR9254Kwm5OZkkd22gZaLFpFvUJIJk5LMt23edZAXZq3lxVnr2LbvMO2b1GJEdhuu6Z5JjepaLlpElGTCpiRzfAePHOWNRRsZm1fA55t2k1GjGjf1bsnw/q1pUU+lNJFEpiQTJiWZk3N35q7dwdi8NUxbugV35+JOTcnNyaJPm/oqpYkkIK0nIxFjZvTOqk/vrPps2HmA52euZcLsdby7dDMdm9UhNzuLQd2ak1ZNpTQR+Sa1ZNSSOS0HDh/ltYUbGJtXwIote6hfszpD+rRkWL/WNMuoEevwRCTKVC4Lk5JM+bg7M1ZvY2xeAR8s20KSGZd1DpXSerSqp1KaSBWlcplUCDMju21Dsts2ZP32/Tw7o4CX5qznzcWbOCczg9ycLK7o0ozUFJXSRBKRWjJqyUTcvkNFTFmwgXF5a1j11T4a1qrOzX1bM6xvKxrXSYt1eCISASqXhUlJJnrcnekrv2ZsXgF/X76VasnGFec0IzenDV1b1o11eCJSDiqXScyZGQPbNWJgu0as+Xof4/MLmDyvkFcXbqR7q7qMyM7i8nOaUS05KdahikiUqCWjlkyF2nPwCJPnFTI+v4CCbftpUieVYX1bM6RvKxrWSo11eCISJpXLwqQkExvFxc4/vviKZ/LW8OmXX1M9JYlBXZszIjuLzpkZsQ5PRE5C5TKJa0lJxvlnNeb8sxqzcusexuev5ZX5hUyeV0jvrHrk5rTh4k5NSFEpTaRSU0tGLZm4sevAESbNXc/4GQWs336A5hlpDOvfmiG9W1GvZvVYhyciJahcFiYlmfhztNj5cNkWxuUXkL9qG6kpSVzTPZMROVmc1bROrMMTEZRkwqYkE9+Wb97N+PwCpszfwKGiYvqf0YAROVlc2LEJyVouWiRmlGTCpCRTOezYd5iX5qznuRkFbNx1kBb1anBr/yxu6NWSjPRqsQ5PJOEoyYRJSaZyKTpazPufb2FsXgGzC7ZTo1oy3+2ZyYjsLM5sXDvW4YkkjIgkGTNLBq5291ciGVw8UZKpvJZs2MX4/AJeW7SRw0XFDGzXkNycLM5r35gkldJEoipiLZngRD0jFlmcUZKp/LbtPcSE2et4buZatuw+RFaDdG7pn8X1vVpQO02lNJFoiGSS+S2wBXgZ2Hdsu7vvLm+Q8UBJpuo4crSYd5ZsZlzeGuav20nN6slc36slt2Zn0aZhzViHJ1KlRDLJrC9js7t7q9MNLp4oyVRNi9bvZFx+AW8u3siRo875HRqRm9OGge0aao0bkQhQx3+YlGSqtq17DvLCzHW8MGsdX+89RNtGNRmRncW1PVpQM1UTXoicrki2ZFKAkcC/BZs+Bv7m7kXlDTIeKMkkhkNFR3n7s02MzStgceEuaqelcGOvltzSP4tWDdJjHZ5IpRPJJPMEUBN4Ntg0DDjo7iPLHWUcUJJJLO7O/HU7GZu3hneWbKbYnQvOasL3crLo37aBSmkiYYrkBJn93L1riefvmdmi0w9NJHbMjJ6t69GzdT027TrACzPX8eLsdXywbAvtm9RiRHYbrumeSY3qWi5aJBLCmeK22Myyjj0Jvi6OTjgiFadZRg1+ekkH8u/9Dn+4rgvJSUn8cupn9Pvth/z2nWVs2Hkg1iGKVHrhlMsuBp4GVgAGnAnc5u4fRD+86FO5TI5xd2av2c64/AKmLd0MwCVnNyU3pw29s+qplCZSQkTKZWaWBOwC2gMdg83L3F1/4kmVY2b0PaMBfc9oQOGO/Tw3cy0vzV7PO0s206lZHUbkZDGoa3PSqqmUJhKucFoyC929WwXFU+HUkpETOXD4KFMXbGBc/hq+2LKX+jWrc3OfVgzr15qmGWmxDk8kZiI5uuwvwMfu/lqkgosnSjISDndnxqptPJNXwIfLt5BsxmXnNGNEdhY9WtVVKU0STiSTzA4gAzgEHCDUL+PuXj8SgcaakoycqnXb9vPsjAJenruePQeL6NIig9ycLC4/pxmpKSqlSWKI1CzMBrQENpTe5+5HyxVhnFCSkdO171ARU+YXMja/gNVf7aNhrVSG9m3F0H6taFxbpTSp2iLZklni7p0jFlmcUZKR8ioudj5d+TXj8tbw0YqvqJZsXNmlObk5WXRpUTfW4YlERSRvxlxoZt3dfUEE4hKpcpKSjHPbN+Lc9o1Y/dVenp2xlklz1zN1wQZ6tKrLiJw2XNa5KdWSw7ktTaRqCaclsxToAKwiNNX/sT6ZHtEPL/rUkpFo2H3wCJPnFjJ+RgFrt+2nSZ1UhvdrzZA+rWhQKzXW4YmUWyTLZW3L2u7uq04ztriiJCPRVFzsfLRiK+PyC/j0y6+pnpLE4K7NGZGTxdnNM2Idnshpi1i5zN1XmVk/oL27P2tmDQhNmCkiJ5GUZFzQsQkXdGzCl1v2MC6/gCnzNzBpXiF9suqTm5PFRZ2akKJSmlRR4bRk7gNygLbu3t7MMoGX3X1ARQQYbWrJSEXbtf8IE+euZ/yMAgp3HKB5RhrD+2cxpE9L6qZXj3V4ImGJZLlsIdAdmO/u3YNti929S0QijTElGYmVo8XOB8u2MDZvDTNXbyetWhLXdM9kRHYbOjStHevwRE4okqPLDrm7m5kHJ9YKTyIRkJxkXHJ2Uy45uynLNu1mfFBKmzB7PdltG5Cb04YLOzbWbAJSqYVTCJ5iZmOADDPLBd4DnoluWCKJpWOzOvzuu12Y8YsL+PmlHVjz9T5uf3Yuo15cwO6DR2IdnshpO2m5DMDMLgMuJjR8eZq7vxPtwCqKymUSj4qOFvPkp6v583tf0KJeDcbc3IPOmRqNJvEjYn0yVZ2SjMSzOQXb+cGLC9i+/zD/eWUnhvZtpfKZxIVwk4zGTYrEsd5Z9Xlr9AD6ndGA+15dwuiXFrL3UFGswxIJm5KMSJxrUCuVcSN687NLOvDW4o1c9fB0Pt+4O9ZhiYQlrCRjZtXN7MxoByMiZUtKMu4+/0xevL0f+w4Vcc2jeUyYvY5EL3dL/DtpkjGzK4DPgPeD593MbGq0AxORb+t3RgPeGj2Q3ln1+cWUz/jRywvZp/KZxLFwWjL/A/QFdgK4+0JArRqRGGlUO5Xx3+vDjy5sz2uLNjLokems2Lwn1mGJlCmcJHPE3XeW2qY2ukgMJScZ91zYjhdu68uuA0UMHjOdiXPXxzoskW8JJ8ksM7MbgCQza2NmfwFmRjkuEQlD9pkNefueAXRvWY+fT17MTyYuYv9hlc8kfoSTZEYBPYFiYApwCPhhNIMSkfA1rp3G89/vy+gL2jFlQSGDH8njyy0qn0l8CCfJtHX3/3D37sHjXnffH/XIRCRsyUnGjy9qz3Pf68uO/YcZ9Eger8wrjHVYImElmTFmtsTM7jezs6IekYictgHtGvLW6IF0aZHBTyYt4ueTF3Hg8NFYhyUJ7KRJxt0HApcAe4DxZrbAzO6NemQiclqa1Enjhe/3ZdT5ZzJpXiFXj8lj5da9sQ5LElRYN2O6+wZ3fwAYQeiemV9HMygRKZ+U5CR+ekkHxuX24au9hxj0yHReXbAh1mFJAgrnZsx2ZnZfsHjZU8AcoFXUIxORcju3fSPeHj2Qzs0z+OHLC/nFlMUcPKLymVSccFoyLwIHgEHuPsDdH3b3TVGOS0QipGlGGi/e3pe7zmvLhNnruXpMHqu/UvlMKkY4fTK93f3P7r6uIgISkchLSU7iPy49i7G5vdmy+yBXPTyd1xdtjHVYkgCOm2TMbELw7wIzm1/iscDM5ldciCISKed3aMxbowdyVrM6jJ6wgF9N/UzlM4mqlBPs+1nw73UVEYiIVIzmdWvw0sh+/GnaCp74ZDUL1+9kzM09yGpYM9ahSRV03JaMux+7k+s2d19V8gHcVjHhiUg0VEtO4heXd+TpW3tRuOMAVz48nbcWq6tVIi+cjv9Ly9h2RaQDEZGKd0HHJrx9z0DaNanF3S/O5z9fW8KhIpXPJHJO1Cdzh5ktADqU6pP5ElhecSGKSDRl1q3ByyP78/0BbXh2xlque2wG67Zp5iiJDDveynpmVg9oAPwWKHmH/x5331oBsVWIXr16+dy5c2MdhkhceG/pZn46aREO/PG6LlzauVmsQ5I4ZWbz3L3XyY47UZ/MDndf6e7XB/0wOwjdL5NiZs0jGKuIxImLz27KW6MHckbDmtz5/Hz+6/WlHC4qjnVYUomFc8f/5Wb2BVAIzALWA3+PdmAiEhst66cz6c5scnOyGJdfwPWP57N+u8pncnrC6fj/PyAHWOHuLQkNBPg0qlGJSExVT0ni/qvO5vFhPVj99T6ueOhT3lu6OdZhSSUUTpIpcvevCK2Mae7+PtAnynGJSBy4tHMz3vrBQFo3qMnI5+bx6zc/V/lMTkk4SWaXmdUCpgPPmtmfCfXNiEgCaNUgncl39efW/q15evoabnhiBoU7VD6T8ISTZK4GDhJacvljYANwVRRjEpE4k5qSzH8P7syYm3uwcuternhoOh8u2xLrsKQSCGeCzD3uXuTuR9z9aXd/ICifiUiCuaJLM978wQAy69bgtvFz+e3byzhyVOUzOb4T3Yy5w8y2l3jsKPlvRQYpIvEjq2FNpvx7NsP6teKJT1Zz05Mz2bhTFXQp24laMg2BRiUeDUv9KyIJKq1aMr+5+hweGtKd5Zt2c8VDn/LRiipzj7ZE0Iluxjx67AH0BoYGX2cAmRUVoIjEr0Fdm/PGDwbQpE4auWPn8Pt3l1Ok8pmUEM7NmPcB9wP3BZtqEFotU0SEMxrV4tW7cxjSpxWPfbyKm5+axeZdB2MdlsSJcEaXXQdcDuwDcPcNQJ1oBiUilUtatWR+e+05PHhjN5Zs3MXlD33KP77Q+CAJL8kc8tAsmg5gZunRDUlEKquru2fy+qgBNKqVyoixs/nTtBUqnyW4cJLMFDMbA2SYWS7wHvBMdMMSkcrqzMah8tn1PVvwyEcrGfq3WWzdrfJZojruVP/fOMjsMuBiwIBp7v5OtAOrKJrqXyR6XplXyH2vLqFmajIP3tidAe0axjokiZByT/UfnCTZzN5393fc/Ufu/sOqlGBEJLq+27MFr4/KoV56dYY/M4sH3v+Co8Un/8NWqo4TJplgyHKymamjX0ROS7smtXltVA7XdM/koQ+/ZPjTs9i6R+WzRBHWBJnAIjN7wsweOPaIdmAiUnWkV0/hgRu68YfrujB/3Q4u/+t08ld+HeuwpAKEk2TeBH4DzAaWlniIiJySG3q15LW7B5BRI4VhT8/ioQ+/VPmsigur478qU8e/SMXbd6iIX039jFcXbmRgu4b85cZuNKyVGuuw5BREpONfRCQaaqam8Jcbu/G7a89h9prtXP7XT5m5elusw5IoUJIRkZgwM27q04pX786hVmoKNz81kzEfraRY5bMqJewkY2Zqy4pIxHVsVofXfzCAK7o054/TVjBi3By27T0U67AkQsKZILOPmX0GfBk872pmD0c9MhFJGLVSU3jopm787zWdmbl6G1c8NJ05BVq2qioIpyXzEHAlsA3A3RcB50czKBFJPGbG0L6tmXJXNmnVkrjpyZk89vEqlc8quXCSTJK7ry217Wg0ghER6ZyZwRs/GMClZzfl9+8u57bxc9ix73Csw5LTFE6SWW9mfQAPppn5IfBFlOMSkQRWO60aj9zcnV8PPpu8ldu4/KFPmbdW5bPKKJwkcxfwY6AVsAXoF2wTEYkaM2N4/yxeuSubaslJ3PjETJ78ZBWJfm9fZaObMXUzpkjc233wCP8xeTHvLNnMhR0b86fru1I3vXqsw0po4d6MedIkY2ZPESxYVpK7jzz98OKHkoxI5eDujM8v4H/fXkbj2mk8fHN3erSqF+uwElYk7/j/APgweOQBjQENYheRCmVmjMhpw+Q7szGDGx6fwd8+Xa3yWZw75XKZmSUB0909OzohVSy1ZEQqn137j/DTyYt4//MtXNypCX+8risZ6dViHVZCiebcZW2AJqfxOhGRiMhIr8aTw3ty3xUd+fvyrVzx8KcsWr8z1mFJGcK543+HmW0PHjuB94FfRD80EZHjMzO+P/AMJt7ZH3e47vF8xuWtUfkszpxs+WUDugKNgkc9dz/D3SdWRHAiIifTo1U93ho9gHPbN+K/3vicf39hPrsPHol1WBI42fLLDrzt7keDh/5EEJG4Uze9Ok/d0otfXn4W732+hSsfms6SDbtiHZYQXp/MQjPrHvVIRETKwcwY+W9tmXhHP44cLebaR/N5bkaBymcxdtwkY2YpwZfdgTlmtsLM5pvZAjObXzHhiYicmp6t6/P26IHknNmA//faUkZNWMAelc9iJuUE+2YDPYBBFRSLiEhE1KtZnadv7c0Tn6zmT++tYOmGXYwZ2oOzm2fEOrSEc6JymQG4+6qyHhUUn4jIaUlKMu46ry0vjezHwSPFXPNoPi/MWqvyWQU7UUumkZn9+Hg73f2BKMQjIhJRvbPq89boAfxo4iJ+NXUJs1Zv5/+uPYdaqSf69SeRcqKWTDJQC6h9nIeISKXQoFYq40b05meXdODNxRsZ9PB0lm3aHeuwEsJxp5Uxs/nu3qOC46lwmlZGJLHMXL2N0RMWsOvAEf570Nnc2LsloVsC5VREYloZfeoiUuX0O6MBb98zkD5t6nPvlM/41atLOHK0ONZhVVknSjIXVFgUIiIVqGGtVMbl9uGu89ry4qx13PL0bC3xHCXHTTLurrVORaTKSk4y/uPSs3jghq7MW7uDqx/NY+XWPbEOq8o5nVmYRUSqjGt7tGDCyH7sO1TENWPy+ccXX8U6pCpFSUZEEl7P1vV49e4cWtRPJ3fsbJ6ZrtmcI0VJRkQEaFEvncl39ufCjk34nzc/55dTNSAgEpRkREQCNVNTeHxYT+4+vy0TZq9j+NOzNCCgnJRkRERKSEoyfnbJWTx4Yzfmr9upAQHlpCQjIlKGq7tn8tLIfuw7dJRrxuTz8YqtsQ6pUlKSERE5jh6t6vHaqBxa1k/ne+Pm8LQGBJwyJRkRkRPIrFuDyXf156JOTfj1m5/ziymfcbhIAwLCVSWTjJnVNLPxZvaUmQ2NdTwiUrmlV0/hsaE9GXX+mbw0Zz3Dn57Fdg0ICEvUkoyZdTCzhSUeu83sh6d5rmfMbKuZLSlj36XBqp0rzezeYPO1wGR3vx0tuiYiEZCUZPz0kg789aZuLFi/k6vH5PHlFg0IOJmoJRl3X+Hu3dy9G9AT2A9MLXmMmTU2s9qltp1ZxunGAZeW3mhmycAY4DKgEzDEzDoBLYD1wWFHy/lWRET+aXC3TF4e2Y8DR45yzaP5fLRcAwJOpKLKZRcAq9x9bant5wKvmlkqgJndDjxc+sXu/glQ1lxqfYCV7r7a3Q8DLwGDgUJCiQaqaElQRGJU5oF6AAAMcElEQVSne6t6vHZ3Dq3qp3Pb+Dn87dPVGhBwHBX1C/gmYELpje4+CZgGvBz0nXwPuP4UzpvJv1osEEoumcAU4Ltm9hjwRlkvNLOrzOzJXbt2ncLlRERCmgcDAi7u1JTfvLWMe1/RgICyRD3JmFl1Qv0ik8ra7+5/AA4CjwGD3H1vea/p7vvcPdfd73L3F45zzBvuPjIjI6O8lxORBJVePYVHh/Zg9HfO5OW56xmmAQHfUhEtmcuA+e6+paydZjYQ6Eyov+b+Uzz3BqBliectgm0iIhUiKcn48cWhAQEL1+9k8JjpfKEBAf9UEUlmCGWUygDMrDvwJKF+lFyggZn95hTOPQdoZ2ZtghbTTcDr5YxXROSUDe6WycQ7+nPwSDHXPprP35eX+Xd1wolqkjGzmsBFhPpIypIO3ODuq9y9GLgFKD04ADObAMwAOphZoZndBuDuRcAoQv06y4CJ7r408u9EROTkurWsy+ujcshqmM5t4+fy1CcaEGCJ/gH06tXL586dG+swRKQK2X+4iJ9OWsTbn23m+p4t+M01nUlNSY51WBFlZvPcvdfJjtPwXhGRCEuvnsIjQ3ow+oJ2TJpXyLC/zWLb3kOxDismlGRERKIgKcn48UXteXhIdxYX7mLwmDxWbE68AQFKMiIiUXRV1+ZMvKM/h4uKufbRPD5cllgDApRkRESirGvLurw+agBnNKrF95+dy5OfrEqYAQFKMiIiFaBpRhoT7+jP5Z2b8X9vL+dnkxdzqKjqT62YEusAREQSRY3qyTw8pDtnNq7FXz/8krXb9vHYsJ40rJUa69CiRi0ZEZEKlJRk/Oii9jxyczAg4JE8lm/eHeuwokZJRkQkBq7s0pxJd/anqLiY7z6azwefV80BAUoyIiIx0qVFXV67ewBtG9fi9ufm8vg/qt6AACUZEZEYapqRxssj+3P5Oc343TvL+cmkRVVqQIA6/kVEYqxG9WQeGdKd9o1r85cPvmDttv08MbxqDAhQS0ZEJA6YGfdc2I4xN/dg6cbQgIBlmyr/gAAlGRGROHJFl2ZMuiM7NCDgsXzeW7o51iGVi5KMiEicOadFBq+PGkC7xrW44/l5PPZx5R0QoCQjIhKHmtRJ4+U7+nNll+b8/t3QgICDRyrfgAB1/IuIxKm0ask8dFM32jWuxQPvf0HB1/t4YngvGtWuPAMC1JIREYljZsboC9rx6NAefL5pN4Mfmc7nGyvPgAAlGRGRSuDyc5ox+c5sih2uezyfaZVkQICSjIhIJdE5M4PXR+WEBgQ8N48xH62M+wEBSjIiIpVI42BAwKCuzfnjtBX8eGJ8DwhQx7+ISCWTVi2Zv97UjfZNavGn976gYNs+nhjek8a102Id2reoJSMiUgmZGaO+047HhvZg+aY9XP1IHks37op1WN+iJCMiUolddk4zJt3ZHweue2wG7y6JrwEBSjIiIpVc58wMXrs7hw5Na3Pn8/E1IEBJRkSkCmhcJ42XRvZjcLfQgIAfvrwwLgYEqONfRKSKSKuWzIM3dqN9k9r8cdoK1m7bz5O3xHZAgFoyIiJViJlx9/ln8viwnqzYvIfBj+SxZEPsBgQoyYiIVEGXdm7K5Lv6Y8D1j8/g3SWbYhKHkoyISBV1dvMMXh11bEDAfB7+8MsKHxCgJCMiUoU1rh0aEHB1t+b8+f0vuOelih0QoI5/EZEqLq1aMn+5sRvtmwYDArbv56nhPWlcJ/oDAtSSERFJAGbGv58XGhDw5ZY9DKqgAQFKMiIiCeSSs5sy+c5skpOMDTsPRP16KpeJiCSYTs3r8OFPziWtWnLUr6WWjIhIAqqIBANKMiIiEkVKMiIiEjVKMiIiEjVKMiIiEjVKMiIiEjVKMiIiEjVKMiIiEjUWL0t0xoqZfQWsPc2XNwS+jmA4IiVlALFbCCQxJPpnXJ7339rdG53soIRPMuVhZnPdvVes45CqycyedPeRsY6jKkv0z7gi3r/KZSLx641YB5AAEv0zjvr7V0umHNSSERE5MbVkyufJWAcgIhLP1JIREZGoUUtGRESiRklGRESiRouWiVQxZlYTeBQ4DHzs7i/EOKQqJ5E/41N972rJRJCZ1TSz8Wb2lJkNjXU8Ejtm1tLMPjKzz81sqZndU45zPWNmW81sSRn7LjWzFWa20szuDTZfC0x299uBQad73XhnZmlmNtvMFgWf8X+X41yV8jM2s2QzW2Bmb5bjHFF970oyJ3G8b0A8/+BJXCgCfuLunYB+wN1m1qnkAWbW2Mxql9p2ZhnnGgdcWnqjmSUDY4DLgE7AkOAaLYD1wWFHy/k+4tkh4Dvu3hXoBlxqZv1KHpAAn/E9wLKydsTLe1eSOblxlPoGVIIfPIkxd9/k7vODr/cQ+kWQWeqwc4FXzSwVwMxuBx4u41yfANvLuEwfYKW7r3b3w8BLwGCgkNDPIlTh/+Mesjd4Wi14lB4uW2U/YzNrAVwB/O04h8TFe6+yP4CRcpxvQNz+4En8MbMsoDswq+R2d58ETANeDsqr3wOuP4VTZ/KvP2og9POXCUwBvmtmj1HF72gPykULga3A++6eSJ/xg8DPgeKydsbLe1fH/+kp68PvCzwEPGJmV1DF/3NLeMysFvAK8EN33116v7v/wcxeAh4D2pb4y/y0ufs+ILe856kM3P0o0M3M6gJTzayzuy8pdUyV+4zN7Epgq7vPM7PzjndcPLx3/bUdQe6+z91z3f2uRBptImUzs2qEEswL7j7lOMcMBDoDU4H7T/ESG4CWJZ63CLYlHHffCXxE2X0LVfEzzgEGmVkBoUrKd8zs+dIHxcN7V5I5PfH6gydxwswMeBpY5u4PHOeY7oSmJhpM6C/DBmb2m1O4zBygnZm1MbPqwE3A6+WLvPIws0ZBCwYzqwFcBCwvdUyV/Izd/Rfu3sLdswjF9Hd3H1bymHh570oypycuf/AkruQAwwn9hbkweFxe6ph04AZ3X+XuxcAtlLG2kZlNAGYAHcys0MxuA3D3ImAUobr7MmCiuy+N3luKO82Aj8xsMaH/k++7e+mhvIn8GcfFe9fcZScRfAPOI7RA2Rbgfnd/OviF8SCQDDzj7v8buyhFROKTkoyIiESNymUiIhI1SjIiIhI1SjIiIhI1SjIiIhI1SjIiIhI1SjIiIhI1SjKSkMxsb/BvlpndHOFz/7LU8/xInj/SzGyEmT0S6zikalKSkUSXBZxSkjGzk00s+40k4+7ZpxhTpRIsfSFSJiUZSXS/AwYG0778KJg6/o9mNsfMFpvZHQBmdp6ZfWpmrwOfB9teNbN5FlqVcWSw7XdAjeB8LwTbjrWaLDj3EjP7zMxuLHHuj81sspktN7MXgrnPviE45vcWWg3yi2Dyw2+1RMzszWMz85rZ3uCaS83sAzPrE5xntZmVXFivZbD9SzO7v8S5hgXXW2hmTxxLKMF5/2xmi4D+kfpmSNWjqf4l0d0L/NTdrwQIksUud+9tocWe8szsveDYHkBnd18TPP+eu28PJmecY2avuPu9ZjbK3buVca1rCa3g2JXQNEVzzOyTYF934GxgI5BHaO6z6WWcI8Xd+wTTGt0PXHiS91eT0OSJPzOzqcBvCE0k2QkYz7/m3OtDaLbe/UFcbwH7gBuBHHc/YmaPAkOBZ4PzznL3n5zk+pLglGREvulioIuZXRc8zwDaAYeB2SUSDMBoM7sm+LplcNy2E5x7ADAhWANli5n9A+gN7A7OXQhgoUW4sig7yRxbMmBecMzJHAbeDb7+DDgUJIzPSr3+fXffFlx/ShBrEdCTUNIBqEFocTAIrfz6ShjXlwSnJCPyTQb8wN2nfWNjqPy0r9TzC4H+7r7fzD4G0spx3UMlvj7K8f9vHirjmCK+WfouGccR/9cEhcXHXu/uxaX6lkpPYuiEPovx7v6LMuI4GCRLkRNSn4wkuj1A7RLPpwF3WWjBMcysvZnVLON1GcCOIMGcBfQrse/IsdeX8ilwY9Dv0wj4N2B2BN5DAaHVIZPMrCWh0tepusjM6gelv6sJlew+BK4zs8YAwf7WEYhXEohaMpLoFgNHgw7sccBfCZWR5ged718R+qVb2rvAnWa2DFgBzCyx70lgsZnNd/ehJbZPJdRJvohQS+Hn7r45SFLlkQesITQgYRkw/zTOMZtQ+asF8Ly7zwUws/uA98wsCTgC3E0Za5KIHI+m+hcRkahRuUxERKJGSUZERKJGSUZERKJGSUZERKJGSUZERKJGSUZERKJGSUZERKJGSUZERKLm/wPa53G7OcehxQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#NBVAL_SKIP\n", "import matplotlib.pyplot as plt\n", "\n", "# Plot objective function decrease\n", "plt.figure()\n", "plt.loglog(relative_error)\n", "plt.xlabel('Iteration number')\n", "plt.ylabel('True relative error')\n", "plt.title('Convergence')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook is part of the tutorial \"Optimised Symbolic Finite Difference Computation with Devito\" presented at the IntelĀ® HPC Developer Conference 2017." ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.5" }, "widgets": { "state": {}, "version": "1.1.2" } }, "nbformat": 4, "nbformat_minor": 1 }