examples/seismic/tutorials/12_time_blocking.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementation of time blocking for compression/serialization and de-serialization/de-compression of wavefields with Devito operators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction \n",
    "\n",
    "The goal of this tutorial is to prototype the compression/serialization and de-serialization/de-compression for wavefields in Devito. The motivation is using seismic modeling operators for full waveform inversion (FWI). Some of the steps in FWI require the use of previously computed wavefields, and of particular interest the adjoint of the Jacobian linearized operator -- an operator that maps data perturbation into velocity perturbation and is used to build FWI gradients -- requires a zero-lag temporal correlation with the wavefield that is computed with the nonlinear source. \n",
    "\n",
    "There are implemented alternatives to serialization/de-serialization like checkpointing, but we investigate the serialization option here. For more information on checkpointing, please see the details for ```pyrevolve```, a python implementation of `optimal checkpointing` for Devito (https://github.com/devitocodes/pyrevolve).\n",
    "\n",
    "We aim to provide a *proof of concept* for compression/serialization and de-serialization/de-compression of the nonlinear wavefield. We will achieve this via _time blocking_: we will run a number of time steps in the generated c kernel, and then return control to Python for compression/serialization (for the nonlinear forward operator), and de-serialization/de-compression (for the linearized Jacobian forward and adjoint operators).\n",
    "\n",
    "In order to illustrate the use case for serialization, we outline the workflow for computing the gradient of the FWI objective function, ignoring a lot of details, as follows:\n",
    "\n",
    "1. Generate the nonlinear forward modeled data at the receivers $d_{mod}$\n",
    "$$ \n",
    "d_{mod} = F m \n",
    "$$\n",
    "\n",
    "1. Compress and serialize the nonlinear source wavefield to disk in *time blocks* during computation of 1. The entire nonlinear wavefield is of size $[nt,nx,nz]$ in 2D, but we deal with a block of $M$ time steps at a time, so the size of the chunk to be compressed and serialized is $[M,nx,nz]$. \n",
    "\n",
    "1. Compute the data residual $\\delta r$ by differencing observed and modeled data at the receivers \n",
    "$$\n",
    "\\delta r = d_{obs} - d_{mod}\n",
    "$$\n",
    "\n",
    "1. Backproject the data residual $\\delta r$ via time reversal with the adjoint linearized Jacobian operator.\n",
    "\n",
    "1. De-serialize and de-compress the nonlinear source wavefield from disk during computation in step 4, synchronizing time step between the nonlinear wavefield computed forward in time, and time reversed adjoint wavefield. We will deal with de-serialization and de-compression of chunks of $M$ time steps of size $[M,nx,nz]$.\n",
    "\n",
    "1. Increment the model perturbation via zero lag correlation of the de-serialized nonlinear source wavefield and the backprojected receiver adjoint wavefield. Note that this computed model perturbation is the  _gradient_ of the FWI objective function.\n",
    "$$\n",
    "\\delta m = \\bigl( \\nabla F\\bigr)^\\top\\ \\delta r\n",
    "$$\n",
    "\n",
    "Please see other notebooks in the `seismic/self=adjoint` directory for more details, in particular the notebooks describing the self-adjoint modeling operators. \n",
    "\n",
    "| Self-adjoint notebook description | name | \n",
    "|:---|:---|\n",
    "| Nonlinear operator | `sa_01_iso_implementation1.ipynb` |\n",
    "| Jacobian linearized operators | `sa_02_iso_implementation2.ipynb` |\n",
    "| Correctness tests | `sa_03_iso_correctness.ipynb` |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outline \n",
    "1. Definition of symbols\n",
    "1. Description of tests to verify correctness\n",
    "1. Description of time blocking implementation\n",
    "1. Compression note -- use of ```blosc```\n",
    "1. Create small 2D test model\n",
    "1. Implement and test the Nonlinear forward operation\n",
    "  - Save all time steps \n",
    "  - Time blocking plus compression/serialization\n",
    "  - Ensure differences are at machine epsilon \n",
    "1. Implement and test the Jacobian linearized forward operation\n",
    "  - Save all time steps \n",
    "  - Time blocking plus compression/serialization\n",
    "  - Ensure differences are at machine epsilon \n",
    "1. Implement and test the Jacobian linearized adjoint operation\n",
    "  - Save all time steps \n",
    "  - Time blocking plus compression/serialization\n",
    "  - Ensure differences are at machine epsilon \n",
    "1. Discussion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of symbols\n",
    "\n",
    "We show the symbols here relevant to the implementation of the linearized operators.\n",
    "\n",
    "| Symbol &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; | Description  | Dimensionality | \n",
    "|:---|:---|:---|\n",
    "| $m_0(x,y,z)$ | Reference P wave velocity    | function of space |\n",
    "| $\\delta m(x,y,z)$ | Perturbation to P wave velocity    | function of space |\n",
    "| $u_0(t,x,y,z)$ | Reference pressure wavefield | function of time and space |\n",
    "| $\\delta u(t,x,y,z)$ | Perturbation to pressure wavefield | function of time and space |\n",
    "| $q(t,x,y,z)$ | Source wavefield | function of time, localized in space to source location |\n",
    "| $r(t,x,y,z)$ | Receiver wavefield | function of time, localized in space to receiver locations |\n",
    "| $\\delta r(t,x,y,z)$ | Receiver wavefield perturbation | function of time, localized in space to receiver locations |\n",
    "| $F[m; q]$ | Forward nonlinear modeling operator | Nonlinear in $m$, linear in $q$: $\\quad$ maps $m \\rightarrow r$ |\n",
    "| $\\nabla F[m; q]\\ \\delta m$ | Forward Jacobian modeling operator | Linearized at $[m; q]$: $\\quad$ maps $\\delta m \\rightarrow \\delta r$ |\n",
    "| $\\bigl( \\nabla F[m; q] \\bigr)^\\top\\ \\delta r$ | Adjoint Jacobian modeling operator | Linearized at $[m; q]$: $\\quad$ maps $\\delta r \\rightarrow \\delta m$ |"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Description of tests to verify correctness\n",
    "\n",
    "In order to make sure we have implemented the time blocking correctly, we numerically compare the output from two runs:\n",
    "1. all time steps saved implementation -- requires a lot of memory to hold wavefields at each time step\n",
    "1. time blocking plus compression/serialization implementation -- requires only enough memory to hold the wavefields in a time block  \n",
    "\n",
    "We perform these tests for three phases of FWI modeling:\n",
    "1. nonlinear forward: maps model to data, *forward in time*\n",
    "1. Jacobian linearized forward: maps model perturbation to data perturbation, *forward in time*\n",
    "1. Jacobian linearized adjoint: maps data perturbation to model perturbation, *backward in time*\n",
    "\n",
    "We will design a small 2D test experiment with a source in the middle of the model and short enough elapsed modeling time that we do not need to worry about boundary reflections for these tests, or runnin out of memory saving all time steps."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Description of time blocking implementation\n",
    "\n",
    "We gratefully acknowledge Downunder Geosolutions (DUG) for in depth discussions about their production time blocking implementation. Those discussions shaped the implementation shown here. The most important idea is to separate the programmatic containers used for propagation and serialization. To do this we utilize two distinct ```TimeFunction```'s.\n",
    "\n",
    "### Propagation uses ```TimeFunction(..., save=None)```\n",
    "\n",
    "We use a default constructed ```TimeFunction``` for propagation. This can be specified in the constructor via either ```save=None``` or no ```save``` argument at all. Devito backs such a default ```TimeFunction``` by a ```Buffer``` of size ```time_order+1```, or 3 for second order in time. We show below the mapping from the monotonic *ordinary time indices* to the buffered *modulo time indices* as used by a ```Buffer``` in a ```TimeFuntion``` with ```time_order=2```.\n",
    "\n",
    "### Modulo indexing for ```Buffer``` of size 3\n",
    "```\n",
    "Ordinary time indices:  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15\n",
    "Modulo   time indices:  0  1  2  0  1  2  0  1  2  0  1  2  0  1  2  0\n",
    "```\n",
    "\n",
    "*Important note:* the modulo indexing of ```Buffer``` is the reason we will separate propagation from serialization. If we use a larger ```Bufffer``` as the ```TimeFunction``` for propagation, we would have to deal with the modulo indexing not just for the current time index, but also previous and next time indices (assuming second order in time). This means that the previous and next time steps can overwrite the locations of the ordinary time indices when you propagate for a block of time steps. **This is the reason we do not use the same ```TimeFunction``` for both propagation and serialization.**\n",
    "\n",
    "### Generated code for a second order in time PDE\n",
    "We now show an excerpt from Devito generated code for a second order in time operator. A second order in time PDE requires two wavefields in order to advance in time: the wavefield at the next time step $u(t+\\Delta t)$ is a function of the wavefield at previous time step $u(t-\\Delta t)$ and the wavefield at the current time step $u(t)$. Remember that Devito uses a ```Buffer``` of size 3 to handle this. \n",
    "\n",
    "In the generated code there are three modulo time indices that are a function of ordinary time and cycle through the values $[0,1,2]$:\n",
    "* ```t0``` -- current time step\n",
    "* ```t1``` -- next time step\n",
    "* ```t2``` -- previous time step\n",
    "\n",
    "We show an excerpt at the beginning of the time loop from the generated code below, with ordinary time loop index ```time```. Note that we have simplified the generated code by breaking a single for loop specification line into multiple lines for clarity. We have also added comments to help understand the mapping from ordinary to modulo time indices.\n",
    "\n",
    "```\n",
    "for (int time = time_m; time <= time_M; time += 1) {\n",
    "  t0 = (time + 0)%(3);    // time index for the current time step \n",
    "  t1 = (time + 1)%(3);    // time index for the next time step\n",
    "  t2 = (time + 2)%(3);    // time index for the  previous time step\n",
    "  \n",
    "  // ... PIE: propagation, source injection, receiver extraction ...\n",
    "}\n",
    "\n",
    "```\n",
    "\n",
    "It should be obvious that using a single container for both propagation and serialization is problematic because the loop runs over ordinary time indices ```time_m``` through ```time_M```, but will access stored previous time step wavefields at indices ```time_m-1``` through ```time_M-1``` and store computed next time step wavefields in indices ```time_m+1``` through ```time_M+1```. \n",
    "\n",
    "### Serialization uses ```TimeFunction(..., save=Buffer(M))``` \n",
    "\n",
    "We will use an independent second ```Buffer``` of size $M$ for serialization in our time blocking implementation. This second ```TimeFunction``` will also use modulo indexing, but by design we only access indices ```time_m``` through ```time_M``` in each time block. This means we do not need to worry about wrapping indices from previous time step or next time step wavefields.\n",
    "\n",
    "### Minimum and maximum allowed time index for second order PDE\n",
    "\n",
    "It is important to note that for a second order in time system the minimum allowed time index ```time_m``` will be $1$, because time index $0$ would imply that the previous time step wavefield $u(t-\\Delta t)$ exists at time index $-1$, and $0$ is the minimum array location.\n",
    "\n",
    "Similarly, the maximum allowed time index ```time_M``` will be $nt-2$, because time index $nt-1$ would imply that the next time step wavefield $u(t+\\Delta t)$ exists at time index $nt$, and $nt-1$ is the maximum array location.\n",
    "\n",
    "### Flow charts for time blocking\n",
    "\n",
    "Putting this all together, here are flow charts outlining control flow the first two time blocks with $M=5$.\n",
    "\n",
    "#### Time blocking for the nonlinear forward\n",
    "\n",
    "```\n",
    "Time block 1\n",
    "    Call generated code Operator(time_m=1, time_M=5)\n",
    "\tReturn control to Python\n",
    "    Compress indices 1,2,3,4,5\n",
    "    Serialize indices 1,2,3,4,5 \n",
    "        (access modulo indices 1%5,2%5,3%5,4%5,5%5)\n",
    "\n",
    "Time block 2\n",
    "    Call generated code Operator(time_m=6, time_M=10)\n",
    "\tReturn control to Python\n",
    "    Compress indices 6,7,8,9,10\n",
    "    Serialize indices 6,7,8,9,10\n",
    "        (access modulo indices 6%5,7%5,8%5,9%5,10%5)\n",
    "```\n",
    "\n",
    "#### Time blocking for the linearized Jacobian adjoint (time reversed)\n",
    "\n",
    "```\n",
    "Time block 2\n",
    "    De-serialize indices 6,7,8,9,10\n",
    "    De-compress indices 6,7,8,9,10\n",
    "        (access modulo indices 6%5,7%5,8%5,9%5,10%5)\n",
    "    Call generated code Operator(time_m=6, time_M=10)\n",
    "\tReturn control to Python\n",
    "\n",
    "Time block 1\n",
    "    De-serialize indices 1,2,3,4,5 \n",
    "    De-compress indices 1,2,3,4,5\n",
    "        (access modulo indices 1%5,2%5,3%5,4%5,5%5)\n",
    "    Call generated code Operator(time_m=1, time_M=5)\n",
    "\tReturn control to Python\n",
    "```\n",
    "\n",
    "## Arrays used to save file offsets and compressed sizes\n",
    "\n",
    "We use two arrays the length of the total number of time steps to save bookeeping information used for the serialization and compression. During de-serialization these offsets and lengths will be used to seek the correct location and read the correct length from the binary file saving the compressed data.\n",
    "\n",
    "| Array | Description |\n",
    "|:---|:---|\n",
    "| file_offset | stores the location of the start of the compressed block for each time step | \n",
    "| file_length | stores the length of the compressed block for each time step | \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compression note -- use of ```blosc```\n",
    "\n",
    "In this notebook we use ```blosc``` compression that is not loaded by default in devito. The first operational cell immediately below ensures that  ```blosc``` compression and Python wrapper are installed in this jupyter kernel. \n",
    "\n",
    "Note that ```blosc``` provides *lossless compression*, and in practice one would use *lossy compression* to achieve significantly better compression ratios. Consider the use of ```blosc``` here as a placeholder for your compression method of choice, providing all the essential characteristics of what might be used at scale. \n",
    "\n",
    "We will use the low level interface to blosc compression because it allows the easiest use of numpy arrays. A synopsis of the compression and decompression calls is shown below for devito ```TimeFunction``` $u$, employing compression level $9$ of the ```zstd``` type with the optional shuffle, at modulo time index ```kt%M```.\n",
    "\n",
    "```\n",
    "    c = blosc.compress_ptr(u._data[kt%M,:,:].__array_interface__['data'][0], \n",
    "                           np.prod(u2._data[kt%M,:,:].shape), \n",
    "                           u._data[kt%M,:,:].dtype.itemsize, 9, True, 'zstd')\n",
    "    \n",
    "    blosc.decompress_ptr(u._data[kt,:,:], d.__array_interface__['data'][0])\n",
    "```\n",
    "\n",
    "```blosc``` Reference:\n",
    "* c project and library: https://blosc.org\n",
    "* Python wrapper: https://github.com/Blosc/python-blosc\n",
    "* Python wrapper documentation: http://python-blosc.blosc.org\n",
    "* low level interface to compression call: http://python-blosc.blosc.org/tutorial.html#compressing-from-a-data-pointer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: blosc in /cpfs/lfs03/ESDRD/wask/local/lib/python3.8/site-packages (1.9.1)\n",
      "\u001b[33mWARNING: You are using pip version 20.2.1; however, version 20.2.2 is available.\n",
      "You should consider upgrading via the '/cpfs/lfs03/ESDRD/wask/local/bin/python3.8 -m pip install --upgrade pip' command.\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "# Install pyzfp package in the current Jupyter kernel\n",
    "import sys\n",
    "!{sys.executable} -m pip install blosc\n",
    "import blosc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports \n",
    "\n",
    "We have grouped all imports used in this notebook here for consistency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from examples.seismic import RickerSource, Receiver, TimeAxis\n",
    "from devito import (Grid, Function, TimeFunction, SpaceDimension, Constant, \n",
    "                    Eq, Operator, configuration, norm, Buffer)\n",
    "from examples.seismic.self_adjoint import setup_w_over_q\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm\n",
    "import copy\n",
    "import os\n",
    "\n",
    "# These lines force images to be displayed in the notebook, and scale up fonts \n",
    "%matplotlib inline\n",
    "mpl.rc('font', size=14)\n",
    "\n",
    "# Make white background for plots, not transparent\n",
    "plt.rcParams['figure.facecolor'] = 'white'\n",
    "\n",
    "# Set logging to debug, captures statistics on the performance of operators\n",
    "# configuration['log-level'] = 'DEBUG'\n",
    "configuration['log-level'] = 'INFO'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Instantiate the model for a two dimensional problem\n",
    "We are aiming at a small model as this is a POC.\n",
    " - 101 x 101 cell model \n",
    " - 20x20 m discretization \n",
    " - Modeling sample rate explicitly chosen: 2.5 msec\n",
    " - Time range 250 milliseconds (101 time steps)\n",
    " - Wholespace model\n",
    "   - velocity: 1500 m/s\n",
    "   - density: 1 g/cm^3\n",
    " - Source to left of center \n",
    " - Vertical line of receivers to right of center \n",
    " - Velocity perturbation box for linearized ops in center of model\n",
    " - Visco-acoustic absorbing boundary from the self-adjoint operators linked above, 10 points on exterior boundaries.\n",
    " - We generate a velocity perturbation for the linearized forward Jacobian operator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shape;            (101, 101)\n",
      "origin;           (0.0, 0.0)\n",
      "spacing;          (20.0, 20.0)\n",
      "extent;           (2000.0, 2000.0)\n",
      "\n",
      "grid.shape;       (101, 101)\n",
      "grid.extent;      (2000.0, 2000.0)\n",
      "grid.spacing_map; {h_x: 20.0, h_z: 20.0}\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `WOverQ_Operator` run in 0.01 s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "time_range;  TimeAxis: start=0, stop=250, step=2.5, num=101\n",
      "\n",
      "src_coordinate  X;            +800.0000\n",
      "src_coordinate  Z;           +1000.0000\n",
      "rec_coordinates X min/max;   +1200.0000   +1200.0000\n",
      "rec_coordinates Z min/max;    +200.0000   +1800.0000\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "# Define dimensions for the interior of the model\n",
    "nx,nz = 101,101\n",
    "npad = 10\n",
    "dx,dz = 20.0,20.0  # Grid spacing in m\n",
    "shape = (nx, nz)   # Number of grid points\n",
    "spacing = (dx, dz) # Domain size is now 5 km by 5 km\n",
    "origin = (0., 0.)  # Origin of coordinate system, specified in m.\n",
    "extent = tuple([s*(n-1) for s, n in zip(spacing, shape)])\n",
    "\n",
    "# Define the dimensions \n",
    "x = SpaceDimension(name='x', spacing=Constant(name='h_x', value=extent[0]/(shape[0]-1)))\n",
    "z = SpaceDimension(name='z', spacing=Constant(name='h_z', value=extent[1]/(shape[1]-1)))\n",
    "\n",
    "# Initialize the Devito grid \n",
    "dtype = np.float32\n",
    "grid = Grid(extent=extent, shape=shape, origin=origin, dimensions=(x, z), dtype=dtype)\n",
    "\n",
    "print(\"shape;           \", shape)\n",
    "print(\"origin;          \", origin)\n",
    "print(\"spacing;         \", spacing)\n",
    "print(\"extent;          \", extent)\n",
    "\n",
    "print(\"\")\n",
    "print(\"grid.shape;      \", grid.shape)\n",
    "print(\"grid.extent;     \", grid.extent)\n",
    "print(\"grid.spacing_map;\", grid.spacing_map)\n",
    "\n",
    "# Create velocity and buoyancy fields. \n",
    "space_order = 8\n",
    "m0 = Function(name='m0', grid=grid, space_order=space_order)\n",
    "b = Function(name='b', grid=grid, space_order=space_order)\n",
    "m0.data[:] = 1.5\n",
    "b.data[:,:] = 1.0 / 1.0\n",
    "\n",
    "# Perturbation to velocity: a square offset from the center of the model\n",
    "dm = Function(name='dm', grid=grid, space_order=space_order)\n",
    "size = 4\n",
    "x0 = (nx-1)//2 + 1\n",
    "z0 = (nz-1)//2 + 1\n",
    "dm.data[:] = 0.0\n",
    "dm.data[(x0-size):(x0+size+1), (z0-size):(z0+size+1)] = 1.0\n",
    "\n",
    "# Initialize the attenuation profile for the absorbing boundary\n",
    "fpeak = 0.001\n",
    "w = 2.0 * np.pi * fpeak\n",
    "qmin = 0.1\n",
    "wOverQ = Function(name='wOverQ_025', grid=grid, space_order=space_order)\n",
    "setup_w_over_q(wOverQ, w, qmin, 100.0, npad)\n",
    "\n",
    "# Time sampling\n",
    "t0 = 0    # Simulation time start\n",
    "tn = 250  # Simulation time end\n",
    "dt = 2.5  # Simulation time step interval\n",
    "time_range = TimeAxis(start=t0, stop=tn, step=dt)\n",
    "nt = time_range.num\n",
    "print(\"\")\n",
    "print(\"time_range; \", time_range)\n",
    "\n",
    "# Source 10 Hz center frequency\n",
    "src = RickerSource(name='src', grid=grid, f0=fpeak, npoint=1, time_range=time_range)\n",
    "src.coordinates.data[0,:] = [dx * ((nx-1) / 2 - 10), dz * (nz-1) / 2]\n",
    "\n",
    "# Receivers: for nonlinear forward and linearized forward\n",
    "# one copy each for save all and time blocking implementations\n",
    "nr = 51\n",
    "z1 = dz * ((nz - 1) / 2 - 40)\n",
    "z2 = dz * ((nz - 1) / 2 + 40)\n",
    "\n",
    "nl_rec1 = Receiver(name='nl_rec1', grid=grid, npoint=nr, time_range=time_range)\n",
    "nl_rec2 = Receiver(name='nl_rec2', grid=grid, npoint=nr, time_range=time_range)\n",
    "ln_rec1 = Receiver(name='ln_rec1', grid=grid, npoint=nr, time_range=time_range)\n",
    "ln_rec2 = Receiver(name='ln_rec2', grid=grid, npoint=nr, time_range=time_range)\n",
    "nl_rec1.coordinates.data[:,0] = nl_rec2.coordinates.data[:,0] = \\\n",
    "    ln_rec1.coordinates.data[:,0] = ln_rec2.coordinates.data[:,0] = dx * ((nx-1) / 2 + 10)\n",
    "nl_rec1.coordinates.data[:,1] = nl_rec2.coordinates.data[:,1] = \\\n",
    "    ln_rec1.coordinates.data[:,1] = ln_rec2.coordinates.data[:,1] = np.linspace(z1, z2, nr)\n",
    "\n",
    "print(\"\")\n",
    "print(\"src_coordinate  X;         %+12.4f\" % (src.coordinates.data[0,0]))\n",
    "print(\"src_coordinate  Z;         %+12.4f\" % (src.coordinates.data[0,1]))\n",
    "print(\"rec_coordinates X min/max; %+12.4f %+12.4f\" % \\\n",
    "      (np.min(nl_rec1.coordinates.data[:,0]), np.max(nl_rec1.coordinates.data[:,0])))\n",
    "print(\"rec_coordinates Z min/max; %+12.4f %+12.4f\" % \\\n",
    "      (np.min(nl_rec1.coordinates.data[:,1]), np.max(nl_rec1.coordinates.data[:,1])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot velocity and density models\n",
    "\n",
    "Next we plot the velocity and density models for illustration, with source location shown as a large red asterisk and receiver line shown as a black line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1kAAAPSCAYAAACJdhAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzde3zP9f//8fvbMLMxh50dZmLsO3xoyzFROTWHPihGzeQwlZQQqT4ZH4fqk09EqZVMOYQOJCUTMo1iy4dEKafEDnzmsDlvr98fft4fbzu98dreb3O7Xi7vy/Z+Pp/v1+vxeueyR4/X6/l6viyGYRgCAAAAAJiijKMDAAAAAIDShCILAAAAAExEkQUAAAAAJqLIAgAAAAATUWQBAAAAgIkosgAAAADARBRZwG2uffv2at++/S27fQDArcNisSg2NtbRYahOnToaOHCgo8NAKUaRBYfIyMjQ888/r9DQULm7u6tixYpq0qSJnn/+eR05csShsf3yyy+KjY3VgQMHSmR/U6dO1fLly0tkX87gzJkzio2N1YYNGxwdCgAUilz1P6U9Vx05ckSxsbHavn27o0NBKVHW0QHg9pOSkqIHHnhAJ0+eVL9+/fTUU0+pTJky2rFjh95//3199tln+u233xwW3y+//KKJEyeqffv2qlOnTrHvb+rUqXrooYf097//vdj35Qhr1qyxeX/mzBlNnDhRkrjCBcBpkatslbZc9euvv6pMmf9dazhy5IgmTpyoOnXqqGnTpg6MDKUFRRZK1MmTJ61/oJOTkxUaGmrTP3XqVL366quOCA3FpHz58o4OAQCuC7mq9HN1dXV0CCjtDKAEvfrqq4Yk46OPPrL7M5988okRFhZmVKhQwahWrZoRGRlpHDx40GZMdHS04erqahw+fNh48MEHDXd3d8PLy8sYPXq0cenSJZuxS5YsMcLDw41KlSoZHh4eRsOGDY1JkyYZhmEY8+bNMyTlec2bN88wDMPYuHGj0adPH6N27dpG+fLlDT8/P2PIkCHG8ePHbfYxYcIEQ5KxZ88eIzo62vD09DQqV65sDBw40MjOzraOy29f7dq1K/C7CA0NNe6+++58++rVq2fz2dzcXOPNN980GjVqZLi6uhre3t7GoEGDjIyMDJvPtWvXLs8+s7OzjTFjxhi1atUyypcvb9SrV8+YNm2akZOTk2e/ixcvNlq0aGFUrFjR8PT0NNq0aWMsX7483+3v378/32OOjo421qxZY0gyPv300zz7WLFihSHJ+Oqrrwr8bgDALOSqm8tVhmEY586dM0aOHGl4eXkZHh4eRvfu3Y0///zTkGRMmDDBZuyRI0eMwYMHG76+vkb58uWNhg0bGm+//bbNmPXr1xuSjEWLFhlTpkwxatSoYbi6uhr33XefsXfvXpuxe/fuNR5++GHDz8/PKF++vOHv72/06tXLOHLkiHVMYGCgER0dbbPta18TJkww4uLiDElGcnJynmOcOXOmIcn45ZdfCv0ucHviShZK1BdffKEKFSqoT58+do1fsGCBoqKiFBYWpmnTpikjI0NvvvmmNm3apJ9++kleXl7Wsbm5uerSpYuaN2+u119/XWvXrtX06dN1xx136IknnpAkrV27VpGRkbrvvvs0bdo0ubi46Ndff9WmTZskSffcc4+efvppvfnmm3rhhRcUEhIiSWrdurUkadmyZTp58qRiYmLk4+NjnTby888/KykpSRaLxSb+yMhI1a1bV9OmTVNKSoref/99+fj4WM+AfvTRRxoyZIiaN2+umJgYSZKvr2+B30dkZKRefvllHT58WDVr1rS2Jycn6/fff9fo0aOtbU888YTmzp2r6OhoPfXUU/rzzz81a9Ys/fjjj9q6dasqVKiQ7z4Mw9Df//53JSQkaNCgQQoLC9O3336r8ePH68CBA3rnnXesYydPnqx//OMfatmypSZMmCA3NzclJyfrm2++0YMPPphn297e3pozZ46eeOIJ9ezZU7169ZIk3XHHHWrRooVq1KihBQsWWNuv/nfg5+enTp06FfjdAIBZyFU3l6skaciQIVqwYIH69++v1q1ba8OGDeratWuecenp6WrZsqVycnL05JNPysfHR99++62efPJJHT9+XC+99JLN+Ndee00uLi4aM2aMTp48qddee02PPPKIfvjhB0nSxYsX1blzZ509e1bDhw+Xv7+/jh49qtWrV+vIkSPy9/fPE0NISIgmTZqkl19+WTExMWrbtq0kqUmTJgoMDNTTTz+tBQsW6M4778zz3z08PNz6/QM2HF3l4fZStWpV429/+5tdYy9cuGD4+voaISEhxpkzZ6ztV844jR492toWHR1tSDImTpxos41mzZoZYWFh1vcjR440KleunOeM4dWWLVtmSDLWr1+fp+/qM3tXLFy40JBkJCYmWtuunB187LHHbMb27NnTqF69uk2bu7u79WxaUfbu3WtIMv7973/btD/33HNG2bJlrVepvv/+e0OSMX/+fJtxiYmJhiTj3XfftbZdeyXrylWj2NhYm88OHDjQkGTs3LnTMAzD+P33340yZcoYPXr0yPN95ubmFrj9jIyMfM9kGoZhjBs3znB1dTX++9//WttOnjxpVKhQwXj22WcL+WYAwDzkqpvLVdu3bzckGU888YRNe1RUVJ6//0OHDjV8fX2N9PR0m7FDhgwx3NzcjMzMTMMw/vd9NmzY0Dh//rx13JWrSVdy05V9L1u2rNAYr76SZRiGsXXrVpurgVfr27ev4efnZ/Pf49dffzUkGTNnzix0P7h9sbogStSpU6dUqVIlu8Zu27ZNaWlpeuKJJ+Tm5mZtb9++vcLCwrRq1ao8nxk6dKjN+7Zt22rfvn3W956ensrOzs6zGIO9KlasKOny1Z5Tp07p2LFj1jOHycnJdsVz/PhxnTp16ob2X69ePYWFhWnJkiU27UuXLtX9999vPVu6dOlSeXh4qEuXLjp27Jj11bBhQ/n6+mr9+vUF7mPVqlUqU6aMnnnmGZv2K1fJrnzvn3/+uXJzc/WPf/xDLi4uNmOvPUtqr+joaJ0/f15Lly61tn3yySc6d+6coqKibmibAHC9yFU3l6u++uorSdJTTz1l0/7000/bvDcMQ5988om6du0qi8Vik686deqks2fPWq9QXTFgwACbe32vXHW68v1VrlxZkvTNN98oOzv7huK/VnR0tFJTU7V27Vpr24IFC1S2bFlFRkaasg+UPhRZKFGVK1fW6dOn7Rp78OBBSVKDBg3y9IWEhORZtrZcuXJ5pgFUrVpVmZmZ1vdPPvmkGjRooIiICNWoUUPR0dH64osvZBiGXTH9+eefioyMlKenpzw9PeXt7a2goCBJl2+Uvlbt2rXzxCPJJqbrFRkZqR9++MF6/Js3b9bBgwdt/tD/9ttvysrKkq+vr7y9vW1eaWlpSk9PL3D7Bw8elK+vr6pUqWLT3qBBA5UpU8a63z/++EOS8twQfjNCQkJ01113acGCBda2BQsWKDQ0VM2aNTNtPwBQGHLVzeWqgwcPymKxqF69ejbtwcHBNu8zMjKUmZmpDz74IE+uujJV89p8VVSsQUFBGjVqlN5//315eXmpQ4cOmjlzpo4fP35DxyJJnTp1kp+fn01uWrhwoTp37iwfH58b3i5KN+7JQokKCQlRSkqKLly4YPqqc1cvxVoQHx8f/fTTT1q7dq2+/vprrV69Wh9++KG6deumL774otArMDk5OerUqZMyMjI0fvx4hYSEyN3d3Tq/Pjc3N89nrr3Cc4W9iTI/ffv21dixY7VkyRKNGzdOS5Yskaurq3r27Gkdk5ubq+rVq+vjjz/OdxtXkpIzGjBggJ5++mkdOHBAZcuW1XfffaepU6c6OiwAtxFy1WU3k6vscSWWfv36adCgQfmOufZEnj2xTp8+XYMGDdIXX3yhNWvWaPTo0Zo8ebK+++47/d///d91x+ni4qJHHnlE77zzjrKzs/Wf//xH+/btIzehUBRZKFE9evRQUlKSli1bpkceeaTQsYGBgZIuP8vi2gUP9uzZc8PPBSlfvrwiIiIUEREhwzA0fvx4vfrqq0pKSlKbNm0KTF47d+7Unj17FB8fr+joaGv73r17byiOK653al2tWrXUunVrLVmyRM8995yWLVumLl26yNPT0zrmjjvuUEJCglq2bCkPD4/r2n5gYKASEhJ08uRJm23+9ttvys3NtX7vd9xxhyRp165dCg8Pt3v7RR1vv379NHr0aC1YsMD6PzdF/VsBADORq/K6nlwVGBgowzD0+++/2xQ11z5XzNvbW5UqVdKlS5fUoUOHm4rvWqGhoQoNDdX48eO1Y8cOhYWF6Y033tB7772X7/iiji86OlrTp0/X559/rqSkJFWuXDnfBZ6AK5guiBI1bNgw1ahRQ6NHj9bu3bvz9J8+fVovvviiJCk8PFy+vr569913de7cOeuYxMREbdu2Td26dbvu/V87XcBisVinoZ04cUKS5O7uLinvNIkrZ8+uPbP3+uuvX3ccV3N3d7/uKRmRkZH66aef9MEHH+jIkSN55oT37dtXubm5mjRpUp7P5uTkFLq/bt26KTc3V2+++aZN+7///W9Jsq4O1bNnT5UpU0aTJk1STk6OzdjCzn5euVegoBiqV6+url27asGCBVq4cKHuvfdem5UUAaC4kavyup5c9cADD0iSZs+ebdM+a9Ysm/cuLi566KGHtHz5cv3nP//Js52MjIzrjvPUqVO6dOmSTVtISIjc3Nys311+Cvo+r2jcuLGaNWumefPmaenSpXr44YcLXKUXkLiShRJWpUoVLV++XBEREbrzzjvVv39/3XXXXSpTpox+/vlnLV68WNWqVdOUKVNUrlw5/etf/9KAAQPUtm1bPfroo9ZlcWvUqKFx48Zd9/6HDBmiY8eO6f7771etWrX0119/afbs2fL399c999wjSWrWrJlcXFw0bdo0nThxQm5ubmrRooUaNmyo+vXra/To0Tp8+LCqVaumr7/+WocPH76p7yQ8PFxr167V66+/rpo1a8rHx0f33XdfoZ95+OGHNXLkSI0aNUoVK1ZU9+7dbfrvueceDR8+XP/617+0Y8cOde7cWa6urvr999/1ySefaNKkSRo4cGC+2+7WrZs6duyoCRMm6ODBg7rzzju1bt06ffrppxo2bJgaNWok6fKVrJdfflmxsbG6++671atXL1WsWFEpKSmqUKGC3nrrrXy37+bmptDQUH388ccKDg5W9erVFRQUpBYtWljHDBgwwDr9MT4+3s5vEgDMQa7K63pyVdOmTdWvXz/NmTNHJ0+eVJs2bbR+/fo8V7Ik6ZVXXtGGDRvUqlUrDR06VKGhocrMzNT27dv1+eef2xSu9li3bp2GDx+uhx56SA0aNJBhGFqyZIlOnz6tvn37Fvi5O+64Q1WrVtWcOXPk4eGhSpUqqVGjRtacJ13OTc8++6wksRgTiuaAFQ0BIy0tzRg7dqwREhJiuLm5GRUqVDAaN25svPDCC0ZqaqrN2CsPeHR1dTWqVq1q9O3b1zhw4IDNmCsPeLzWleVpr95Wly5drA88rFWrlhEdHW3s27fP5nMffPCBUa9ePcPFxcVmSdc9e/YYnTt3NipXrmxUrVrV6Nevn5GamppnSdor+z169KjNdq88QHL//v3Wtt9++8249957DXd3d7se8HjF/fffb0gy+vbtW+CYuXPnGnfddZfh5uZmVKpUyQgNDTVGjx5t84DM/B5GnJWVZYwePdqoUaOGUa5cOeOOO+4wpk6dmu9ywvPnz7c+gLNKlSpGmzZtjBUrVhS6/S1bthjNmzc3XF1drQ8jvtqFCxcMLy8vo2LFisbp06ft+j4AwGzkqv3WtuvNVWfPnjWefvppo3r16oa7u3uhDyNOT083RowYYdSuXdsoV66c4evra7Rv396YPXu2dcyVJdwXL15s89krD7m/cuz79u0zBg8ebNSrV89wc3Mzqlatatx9993G8uXLbT537RLuhmEYX375pdG4cWOjXLly+caZlpZmlC1b1ggMDLR5VAmQH4thFPNdjQBwnXJyclSzZk3dd999WrhwoaPDAQBAmZmZ8vPz05gxYzRlyhRHhwMnxz1ZAJzOypUrlZqaanPTNgAAjjR//nxduHCB3AS7lOoi6+2331ZQUJAqVKigsLAwJSYmOjokAIX44Ycf9P7772vkyJFq3LixOnbs6OiQgGJBfgJuHevWrdNbb72lSZMmqVu3bnme9wXkp9QWWUuWLNEzzzyjF154QT/99JNat26tBx54QIcOHXJ0aAAKMGfOHD3++OOqVq2aFixYcN3L2wO3AvITcGuZNGmSnn32WTVu3Fhvv/22o8PBLaLU3pPVokULNWnSxOZ5CPXr19dDDz2kadOmOTAyAMDtjPwEAKVfqbySdeHCBSUnJ+d5KGCnTp2UlJTkoKgAALc78hMA3B5K5XOyjh07ppycHPn6+tq0+/r6au3atXnGx8XFKS4uTpKU/MseKahhicQJAE7lyAEZmcccHUWpdjP5aceeHfJqWL1E4gQAZ3PhwEUdO3br5KhSWWRdr5iYGMXExEiSLKHh0tJtDo4IABygT7ijI8A1rs5PAeH+GrztMQdHBACO8XX4GkeHcF1K5XRBLy8vubi4KC0tzaY9LS1Nfn5+DooKAHC7Iz8BwO2hVBZZ5cuXV1hYmBISEmzaExIS1Lp1awdFBQC43ZGfAOD2UGqnC44aNUpRUVFq3ry52rRpo3feeUdHjhzR448/7ujQAAC3MfITAJR+pbbI6tu3r44fP67Jkyfr6NGjatSokb766isFBgY6OjQAwG2M/AQApV+pfU7WjWLhCwC3rT7hMnbx989ZsfAFgNvZ1+FrtG3brZOjSuU9WQAAAADgKBRZAAAAAGAiiiwAAAAAMBFFFgAAAACYiCILAAAAAExEkQUAAAAAJqLIAgAAAAATUWQBAAAAgIkosgAAAADARBRZAAAAAGAiiiwAAAAAMBFFFgAAAACYiCILAAAAAExEkQUAAAAAJqLIAgAAAAATUWQBAAAAgIkosgAAAADARBRZAAAAAGAiiiwAAAAAMBFFFgAAAACYiCILAAAAAExEkQUAAAAAJqLIAgAAAAATUWQBAAAAgIkosgBc9srIyy8AAJzMmpEJWjMywdFhAHYr6+gAADiJPdsdHQEAAPlK3Z7u6BCA68KVLAAAAAAwEUUWAAAAAJiIIgsAAAAATORURda0adN01113qXLlyvL29lb37t31888/24wZOHCgLBaLzatly5Y2Y86fP68RI0bIy8tL7u7u6tGjhw4fPlyShwIAKGXIUQAAezlVkbVhwwY9+eSTSkpK0rp161S2bFl16NBB//3vf23GdejQQUePHrW+vvrqK5v+kSNH6tNPP9XixYuVmJioU6dOqVu3bsrJySnJwwEAlCLkKACAvZxqdcFvvvnG5v1HH30kT09Pff/99+revbu13dXVVX5+fvlu4+TJk5o7d67mzZunjh07WrcTGBiotWvXqnPnzsV3AACAUoscBQCwl1NdybrW6dOnlZubq6pVq9q0b9q0ST4+PgoODtbQoUOVnv6/ZT2Tk5N18eJFderUydpWq1YthYSEKCkpqcRiBwCUbuQoAEBBnOpK1rWeeeYZNW3aVK1atbK2denSRb169VJQUJAOHDigl156Sffdd5+Sk5Pl6uqq1NRUubi4yMvLy2Zbvr6+Sk1NzXc/cXFxiouLu/wmM6PYjgcAUHqURI66Oj9lZ5wp3gMCAJjGaYusUaNGadOmTdq0aZNcXFys7ZGRkdbfGzdurLCwMAUGBmrVqlXq1avXDe0rJiZGMTExkiRLaPjNBQ4AKPVKKkddnZ8Cwv1vPnAAQIlwyumCzz77rBYvXqx169apbt26hY4NCAhQzZo1tXfvXkmSn5+fcnJydOzYMZtxaWlpBc6RBwDAXuQoAEBRnK7IeuaZZ6zJq2HDhkWOP3bsmP766y/5+18+wxcWFqZy5copISHBOubw4cPavXu3WrduXWxxAwBKP3IUAMAeTjVdcPjw4froo4+0fPlyVa1a1To/3cPDQx4eHsrKylJsbKx69+4tf39/HThwQOPHj5ePj4969uwpSfL09NTgwYM1duxY+fj4qHr16ho1apSaNGmiDh06OPLwAAC3MHIUAMBeTlVkvf3225Kk+++/36Z9woQJio2NlYuLi3bu3KkPP/xQJ06ckL+/v+69914tXbpUlSpVso6fMWOGypYtq759++rs2bO6//779eGHH9rMmwcA4HqQowAA9nKqIsswjEL73dzc8jynJD+urq6aNWuWZs2aZVZoAIDbHDkKAGAvp7snCwAAAABuZRRZAAAAAGAiiiwAAAAAMBFFFgAAAACYiCILAAAAAExEkQUAAAAAJqLIAgAAAAATUWQBAAAAgIkosgAAAADARBRZAAAAAGAiiiwAAAAAMBFFFgAAAACYiCILAAAAAExEkQUAAAAAJqLIAgAAAAATUWQBAAAAgIkosgAAAADARBRZAAAAAGAiiiwAAAAAMBFFFgAAAACYiCILAAAAAExEkQUAAAAAJqLIAgAAAAATUWQBAAAAgIkosgAAAADARBRZAAAAAGAiiiwAAAAAMBFFFgAAAACYiCILAAAAAEzkVEVWbGysLBaLzcvPz8/abxiGYmNjFRAQIDc3N7Vv3167du2y2UZmZqaioqLk6ekpT09PRUVF6cSJEyV9KACAUoYcBQCwl1MVWZLUoEEDHT161PrauXOnte+1117T9OnTNWvWLG3dulU+Pj7q2LGjTp8+bR3Tv39/paSkaPXq1Vq9erVSUlIUFRXliEMBAJQy5CgAgD3KOjqAa5UtW9bmzOAVhmFoxowZev7559W7d29J0vz58+Xj46NFixZp2LBh2r17t1avXq1NmzapVatWkqR3331Xbdu21a+//qoGDRqU6LEAAEoXchQAwB5OdyVr3759CggIUFBQkCIjI7Vv3z5J0v79+5WamqpOnTpZx7q5uemee+5RUlKSJGnz5s3y8PBQ69atrWPatGkjd3d36xgAAG4UOQoAYA+nKrJatGih+Ph4rV69Wu+9955SU1PVunVrHT9+XKmpqZIkX19fm8/4+vpa+1JTU+Xt7S2LxWLtt1gs8vHxsY7JT1xcnMLDwxUeHi5lZhTDkQEAbnWOyFFX56fsjDPFdGQAALM51XTBBx54wOZ9y5YtVbduXc2fP18tW7Ystv3GxMQoJiZGkmQJDS+2/QAAbl2OyFFX56eAcP9i2QcAwHxOdSXrWh4eHgoNDdXevXutc+DT0tJsxqSlpVn7/Pz8lJGRIcMwrP2GYSg9PT3fOfQAANwochQAoCBOXWSdO3dOe/bskb+/v4KCguTn56eEhASb/sTEROv89latWikrK0ubN2+2jtm8ebOys7Nt5sADAHCzyFEAgII41XTBMWPGqHv37qpdu7bS09P1z3/+U9nZ2YqOjpbFYtHIkSM1depUNWzYUMHBwZo8ebI8PDzUv39/SVJISIi6dOmiYcOGKS4uTpI0bNgwdevWjVWbAAA3hRwFALCXUxVZhw8fVr9+/XTs2DF5e3urZcuW2rJliwIDAyVJY8eO1dmzZzV8+HBlZmaqRYsWWrNmjSpVqmTdxqJFizRixAh17txZktSjRw/Nnj3bIccDACg9yFEAAHtZjKsnh+PywhdLtzk6DKDkDWx/+Wf8BkdGAUfqEy5jF3//nFVAuL8Gb3vM0WEADvFh+4WSpAEbHnFwJHCUr8PXaNu2WydHOfU9WQAAAABwq6HIAgAAAAATUWQBAAAAgIkosgAAAADARBRZAAAAAGAiiiwAAAAAMBFFFgAAAACYiCILAAAAAExEkQUAAAAAJqLIAgAAAAATUWQBAAAAgIkosgAAAADARBRZAAAAAGAiiiwAAAAAMBFFFgAAAACYqKw9g3Jzc7VhwwZ99913OnDggM6ePStvb2/deeed6tSpk2rVqlXccQIAkAf5CQDgjAq9knX27FlNmTJFtWrVUteuXbVmzRplZWWpfPny2r9/vyZOnKigoCBFRERoy5YtJRUzAOA2R34CADizQq9k1a9fX61atVJcXJw6deqkcuXK5Rlz8OBBLVq0SH379tVLL72koUOHFluwAABI5CcAgHMrtMhavXq1GjVqVOgGAgMDNX78eI0ePVoHDx40NTgAAPJDfgIAOLNCpwsWlcCuVr58edWvX/+mAwIAoCjkJwCAM7Nr4YsrLly4oJ9//lnp6enKzc216YuIiDA1MAAA7EV+AgA4E7uLrISEBEVFRSk9PT1Pn8ViUU5OjqmBAQBgD/ITAMDZ2P2crOHDh6tbt27av3+/zpw5o7Nnz1pfZ86cKc4YAQAoEPkJAOBs7L6SdfToUb3wwgsKDAwszngAALgu5CcAgLOx+0pWt27dlJSUVJyxAABw3chPAABnY/eVrHfeeUePPPKIkpOT1ahRozzPJBkwYIDpwQEAUBTyEwDA2dhdZH3zzTf69ttv9dVXX6lixYqyWCzWPovFQhIDADgE+QkA4Gzsni44ZswYPfXUUzp9+rSysrJ0+vRp6+vUqVPFGSMAAAUiPwEAnI3dRdaJEyf0+OOPy93dvTjjAQDgupCfAADOxu4iq3fv3lq7dm1xxqI6derIYrHkeXXt2lWSFBsbm6fPz8/PZhuGYSg2NlYBAQFyc3NT+/bttWvXrmKNGwDgOCWRnyRyFADAfnbfk1W3bl29+OKL2rhxo5o0aZLnxuJRo0bddDBbt261eWjk0aNHFRYWpj59+ljbGjRooA0bNljfu7i42Gzjtdde0/Tp0xUfH68GDRpo0qRJ6tixo3799VdVqlTppmMEADiXkshPEjkKAGA/u4usDz74QJUqVVJSUlKepXItFospSczb29vm/dy5c1W5cmWbBFa2bNk8ZwavMAxDM2bM0PPPP6/evXtLkubPny8fHx8tWrRIw4YNu+kYAQDOpSTyk0SOAgDYz+4ia//+/cUZRx6GYWju3Ll69NFH5ebmZm3ft2+fAgIC5OrqqhYtWmjq1KmqW7euNcbU1FR16tTJOt7NzU333HOPkpKSSGAAUAqVdH6SyFEAgMLZfU9WSUtISND+/fs1dJe99A4AACAASURBVOhQa1uLFi0UHx+v1atX67333lNqaqpat26t48ePS5JSU1MlSb6+vjbb8vX1tfblJy4uTuHh4QoPD5cyM4rhaAAApUlJ5air81N2xpliOhoAgNkKLbImT56s7Oxsuzb0/fffa+XKlaYEJUnvvfee7rrrLv3tb3+ztj3wwAPq06ePmjRpog4dOujLL79Ubm6u5s+ff1P7iomJ0bZt27Rt2zapqnfRHwAAOJQj85NUcjnq6vzk7l3RjNABACWg0CLrjz/+UO3atRUTE6OVK1fq6NGj1r5z584pJSVFb775ppo3b66oqChVrVrVlKDS09O1YsUKmzOE+fHw8FBoaKj27t0rSdZ58GlpaTbj0tLSCpwjDwC49TgqP0nkKABA0QotsubNm6cNGzbIYrFowIABqlmzpsqWLSs3Nze5u7srPDxcH374oYYMGaLdu3fr7rvvNiWo+Ph4ubq6ql+/foWOO3funPbs2SN/f39JUlBQkPz8/JSQkGAzJjExUa1btzYlNgCA4zkqP0nkKABA0Ypc+KJx48Z69913NWfOHO3YsUMHDx7U2bNn5eXlpaZNm8rLy8vUgAzD0Pvvv6/IyEh5eHjY9I0ZM0bdu3dX7dq1lZ6ern/+85/Kzs5WdHS0pMurSI0cOVJTp05Vw4YNFRwcrMmTJ8vDw0P9+/c3NU4AgGOVdH6SyFEAAPvYvbpgmTJl1LRpUzVt2rQ449GGDRu0d+9eLViwIE/f4cOH1a9fPx07dkze3t5q2bKltmzZosDAQOuYsWPH6uzZsxo+fLgyMzPVokULrVmzhuePAEApVVL5SSJHAQDsYzEMw3B0EM7EEhouLd3m6DCAkjew/eWf8RscGQUcqU+4jF38/XNWAeH+GrztMUeHATjEh+0XSpIGbHjEwZHAUb4OX3N5kbpbhNMu4Q4AAAAAtyKKLAAAAAAwEUUWAAAAAJjohoqstLQ05ebmmh0LAAA3hfwEAHAGdhdZFy9e1NixY1WpUiXVqFFDBw4ckCSNGzdOb7/9dnHFBwBAochPAABnY3eRNXHiRK1cuVILFiyQq6urtb158+aKj48vjtgAACgS+QkA4Gzsfk7W4sWL9cEHH6hdu3YqU+Z/tVmjRo3022+/FUtwAAAUhfwEAHA2dl/JOnLkiM0DFa+4dOmSLl26ZGpQAADYi/wEAHA2dhdZoaGh2rhxY572pUuXKiwszNSgAACwF/kJAOBs7J4uOGHCBD366KP6888/lZOTo2XLlmnPnj1atGiRVq1aVZwxAgBQIPITAMDZ2H0lq3v37lq6dKnWrFmjMmXKaOLEidq7d69WrlypDh06FGeMAAAUiPwEAHA2dl/JkqTOnTurc+fOxRULAAA3hPwEAHAmdl/Jqlu3ro4fP56n/cSJE6pbt66pQQEAYC/yEwDA2dhdZB04cEA5OTl52s+fP6+//vrL1KAAALAX+QkA4GyKnC742WefWX9ftWqVPD09re9zcnL07bffqk6dOsUSHAAABSE/AQCcVZFF1kMPPSRJslgsGjx4sE1fuXLlVKdOHU2fPr14ogMAoADkJwCAsyqyyMrNzZUkBQUFaevWrfLy8ir2oAAAKAr5CQDgrOxeXXD//v3FGQcAADeE/AQAcDbXtYR7Zmamvv76ax06dEgXLlyw6Xv55ZdNDQwAAHuRnwAAzsTuImvLli3q2rWrXF1dlZGRoRo1aujo0aNydXVVnTp1SGIAAIcgPwEAnI3dS7g/99xzeuSRR/TXX3+pQoUKWrdunQ4dOqTw8HCNGzeuOGMEAKBA5CcAgLOxu8jasWOHnnrqKVksFrm4uOj8+fPy9fXVq6++qtjY2GIMEQCAgpGfAADOxu4iq3z58tbffX19dfDgQUmSh4eHjhw5Yn5kAADYgfwEAHA2dt+Tdeedd2rr1q0KDg5W+/bt9dJLLyktLU0LFixQkyZNijNGAAAKRH4CADgbu69kTZkyRQEBAZKkyZMny9vbWyNGjFBmZqbi4uKKLUAAAApDfgIAOBu7r2SFh4dbf/f29tbXX39dLAEBAHA9yE8AAGdj95UsoLSrduK4vh7WRdVOHHd0KAAAWLkdP6OoLh/L7fgZR4cCwE52F1n//e9/9cQTTyg4OFhVqlRR5cqVbV7ArS56ebw6f/+NBqyY7+hQAFwH8hNKu2bxO1X/m/1qNn+no0MBYCe7pwsOHjxYP/30k2JiYhQQECCLxVKccQElyzA06qM3ZJE06sM3NGPAsxL/xoFbAvkJpZphqPUbW2WR1PqNrUp6tjn5CbgF2H0l69tvv9WSJUv0wgsvaODAgYqOjrZ52WPjxo3q0aOHatSoIYvFovj4eJt+wzAUGxurgIAAubm5qX379tq1a5fNmMzMTEVFRcnT01Oenp6KiorSiRMnbMbs3LlT7dq1k5ubm2rUqKFJkybJMAx7DxW3obbJifI8fVKSVOX0Cd2dssnBEQGwlxn5SSJHwTkFJv6pCifPSZIqnDinwE2HHRwRAHvYXWT5+PjIw8PjpnaWlZWlRo0aaebMmXJzc8vT/9prr2n69OmaNWuWtm7dKh8fH3Xs2FGnT5+2junfv79SUlK0evVqrV69WikpKYqKirL2nzp1Sh07dpSvr6+2bt2qmTNn6l//+pf+/e9/31TsKN1GfjRDFc9mS5Iqns3Wsx++4eCIANjLjPwkkaPgnFrP2Kpy2RclSeWyL6r1Gz86OCIA9rAYdp4+W7JkiZYuXar58+ebksw8PDw0e/ZsDRw4UNLlM4QBAQF66qmn9OKLL0qSzp49Kx8fH73++usaNmyYdu/erf/7v//Tpk2b1KZNG0nSpk2b1LZtW+3Zs0cNGjTQnDlzNG7cOKWlpVmT5OTJkzVnzhwdPny4yGkkltBwaem2mz4+OK/lIx7Ug+u/sGk7X668XC9eKPC9JK24t4f+PmtFicToEAPbX/4Zv8GRUcCR+oTL2HXr/f0zOz9JzpmjAsL9NXjbY6YcH5xT/weXKeSL323aLpUvo7IXcgt8L0m7e9TXohUPlUiMjvJh+4WSpAEbHnFwJHCUr8PXaNu2WydH2X0la/LkyVqzZo18fHwUEhKiJk2a2Lxu1v79+5WamqpOnTpZ29zc3HTPPfcoKSlJkrR582Z5eHiodevW1jFt2rSRu7u7zZi2bdvanIXs3Lmzjhw5ogMHDtx0nLj1vfDMVB3wr62zrhWsbdcWVFe/P+NaQQcCAvXCM1NLLEYA9ivu/CSRo1Ay1k5tr8zalXWxwv9umb+2oLr6/cUKZZUZWFlrp7YrsRgB2MfuhS8eeqh4z5CkpqZKknx9fW3afX199ddff1nHeHt725zps1gs8vHxsX4+NTVVNWvWzLONK31BQUF59h0XF/e/B1ZmZphzQHBav9QLVeiKX/TBPwap23dfyv1cwUviZrlV1Mp23TVk0lydqeheglECsFdx5yfJcTnq6vyUncHy3aVdeqi3Zv0yVD0HrVKDL/9Q+TMXCxx7oWI57eleT8vnRuiie/kSjBKAPewusiZMmFCccThUTEyMYmJiJP3/6YIo9c5UdFfk9CWKWfqOZr4yUhUunM8z5lx5V41+brri+jzugAgB2Ot2yU8B4f4OjgYl4aJ7eS1d0lPh7/ykiJEJKnc+J8+YS64u+nr6/dr2eDMHRAjAHk7zMGI/Pz9JUlpamk17Wlqatc/Pz08ZGRk2qzAZhqH09HSbMflt4+p9AFf8FHKnzpd3zbfvfDlXpYSElXBEAJwROQol7eidvspxdcm375Kri46G+ebbB8A5FFpkVa5cWceOHZMkVapUKc8DHs182GNQUJD8/PyUkJBgbTt37pwSExOt89tbtWqlrKwsbd682Tpm8+bNys7OthmTmJioc+fOWcckJCQoICBAderUuek4UbqE79qmspcuT8fItViUXaGicv//VJ+yORcVfgsuAgDcDkoyP0nkKJS8GtuOqszFy/df5VqkCxXLKvf/z0QtczFXAdtSHRgdgKIUOl1w1qxZqlSpkiRp9uzZN72zrKws/f775VVzcnNzdejQIW3fvl3VqlVT7dq1NXLkSE2dOlUNGzZUcHCwJk+eLA8PD/Xv31+SFBISoi5dumjYsGHWOerDhg1Tt27d1KBBA0mXl8+dOHGiBg4cqJdeekm//fabXnnlFU2YMIEHVCKPtsmJcj93VmdcKyi9uq9GjpuhGa88I9//psv93FndnZKodyKfcHSYAK5hdn6SyFFwLoGJf6r82Uu6WKGssnwr6usZHfXAMwnySD+j8mcvKTDxT2194k5HhwmgAHYv4W6GDRs26N57783THh0drfj4eBmGoYkTJ+rdd99VZmamWrRoobfeekuNGjWyjs3MzNSIESP0xReXl+Du0aOHZs+erSpVqljH7Ny5U8OHD9ePP/6oqlWr6vHHH9fLL79sVwJjCffbyx+d66p26iEt6/iQdXGLimey9cE/Bql3wqc6GBCoeqv/cHSYJYMl3HGLLuFuFmfPUSzhfnt5tu7bqnLolH5+qKF1cYty2RfUc9AqhX76q04EeuqNP26fk4As4Y5bbQn3Ei2ybgUUWbeXL5/oqk869VZ8z0F5+gZ+/oEeWvOpus1Z5YDIHIAiC7d5keXsKLJuL492XapdvRvop0F/y9PX7IP/KPTTX7VgVR8HROYYFFm41YqsQqcLlilTxu7pCzk5eVe/AZxdYQVUfM9B+RZfAByP/ITSrrAC6qdBf8u3+ALgPAotspYuXWpNYmlpaXr55ZfVs2dPtWrVStLlG3qXL1+uiRMnFn+kAAD8f+QnAIAzK7TIuvoBjz169NC0adM0dOhQa9ugQYPUvHlzLV++XE8++WTxRQkAwFXITwAAZ2b3c7LWrVuX7w3B9957rzZs2GBmTAAA2I38BABwNnYXWV5eXvrkk0/ytH/yySfy9vY2NSgAAOxFfgIAOJtCpwtebdKkSXrssce0fv1665z3LVu2aO3atZo7d26xBQgAQGHITwAAZ2N3kTVgwAA1bNhQM2fOtD7/IyQkRN9//71atGhRbAECAFAY8hMAwNnYVWRdvHhRjz76qKZOnaqFCxcWd0wAANiF/AQAcEZ23ZNVrlw5rVmzxu5nkgAAUBLITwAAZ2T3whe9evXSZ599VpyxAABw3chPAABnY/c9WbVr19bkyZOVmJio8PBwubu72/SPGjXK9OAAACgK+QkA4GzsLrLi4+NVtWpV7dixQzt27LDps1gsJDEAgEOQnwAAzsbuImv//v3FGQcAADeE/AQAcDZ235N1taysLGVnZ5sdCwAAN4X8BABwBtdVZL311luqXbu2PD09VblyZQUGBurtt98urtgAALAL+QkA4Ezsni44depUTZs2TWPGjNHdd98tSUpMTNTzzz+vU6dO6fnnny+2IAEAKAj5CQDgbOwust555x3FxcWpX79+1rb7779f9evX1wsvvEASAwA4BPkJAOBs7J4umJ6errvuuitPe/PmzZWWlmZqUAAA2Iv8BABwNnYXWcHBwVq0aFGe9kWLFqlBgwamBgUAgL3ITwAAZ2P3dMHY2Fj16dNHGzduVJs2bSRJ33//vb777jstW7as2AIEAKAw5CcAgLOxu8jq1auXfvjhB73xxhv68ssvJUkhISH68ccf1axZs2ILEACAwpCfAADOxu4iS5LCwsK0YMGC4ooFAIAbQn4CADiT6yqyzp8/r4ULF+qXX36RxWJRaGio+vXrJ1dX1+KKDwCAIpGfAADOxO6FL3755RfVr19fo0aN0g8//KAtW7Zo5MiRCg4O1u7du4szRgAACkR+AgA4G7uLrGeeeUbNmjXToUOHlJiYqMTERB06dEh/+9vfNHLkyOKMEQCAApGfAADOxu7pgt9//722bt2qypUrW9sqV66sKVOmqGXLlsUSHAAARSE/AQCcjd1XsipUqKATJ07kaT958qQqVKhgalAAANiL/AQAcDZ2F1ndu3fX0KFD9f333ysnJ0c5OTnatGmThg0bph49ehRnjAAAFIj8BABwNnYXWTNnzlT9+vXVtm1bVahQQRUqVFC7du0UHBysGTNmFGeMAAAUiPwEAHA2dhdZVapU0YoVK/Tbb7/ps88+02effaZff/1Vn3/+uTw9Pe3axsaNG9WjRw/VqFFDFotF8fHx1r6LFy9q3LhxatKkidzd3eXv76/+/fvr0KFDNtto3769LBaLzSsyMtJmTGZmpqKiouTp6SlPT09FRUXlO5UEAHDrIz8BAJyNXUXWoUOHdOHCBUlSvXr11L17d3Xv3l21a9fOk2QKk5WVpUaNGmnmzJlyc3Oz6Ttz5oxSUlL04osvKiUlRStWrNCff/6pLl266NKlSzZjH3vsMR09etT6evfdd236+/fvr5SUFK1evVqrV69WSkqKoqKi7I4TAHBrID8BAJxRkasLfvzxx5o8ebJSUlLy9BmGoYiICE2YMEEPP/xwkTuLiIhQRESEJGngwIE2fZ6enkpISLBpe/fddxUaGqrdu3ercePG1vaKFSvKz88v333s3r1bq1ev1qZNm9SqVSvrdtq2batff/1VDRo0KDJOAIDzIz8BAJxVkVey4uLi9Nxzz6l8+fJ5+lxdXTVu3DjNmTOnWII7deqUJKlq1ao27R9//LG8vLwUGhqqMWPG6PTp09a+zZs3y8PDQ61bt7a2tWnTRu7u7kpKSiqWOAEAJY/8BABwVkVeydq9e7fuvvvuAvtbt26t5557ztSgJOnChQsaPXq0unfvrpo1a1rb+/fvr8DAQAUEBGjXrl0aP368duzYoTVr1kiSUlNT5e3tLYvFYv2MxWKRj4+PUlNT891XXFyc4uLiLr/JzDD9WAAA5rvd8lN2xhnTjwUAUDyKLLJOnjxpne+en/Pnz1vP6Jnl0qVLevTRR3XixAl98cUXNn0xMTHW3xs3bqy6deuqRYsWSklJ0Z133nlD+4uJibFu1xIafuOBAwBKzO2WnwLC/W88cABAiSpyumBQUJC2bt1aYP+PP/6oOnXqmBbQpUuX1K9fP+3YsUPffvutqlevXuj48PBwubi4aO/evZIkPz8/ZWRkyDAM6xjDMJSenl7gPHkAwK2H/AQAcFZFFlm9evXSiy++qKNHj+bpO3LkiP7xj3+od+/epgRz8eJF9e3bVzt27ND69evtSjo7d+5UTk6O/P0vn+Fr1aqVsrKytHnzZuuYzZs3Kzs722YePADg1kZ+AgA4qyKnC44bN06ff/65goOD9eijj6phw4aSLs+FX7hwoWrXrq2xY8fatbOsrCz9/vvvkqTc3FwdOnRI27dvV7Vq1RQQEKCHH35YW7du1cqVK2WxWKxz1D09PeXm5qY//vhDCxcuVEREhLy8vPTLL79o9OjRatasmdq0aSNJCgkJUZcuXTRs2DDrPPZhw4apW7durNwEAKUI+QkA4KwsxtXzFgpw8uRJjR8/XkuWLFFmZqakyw9/jIyM1JQpU/KsrlSQDRs26N57783THh0drdjYWAUFBeX7uXnz5mngwIH6888/9eijj+rnn39WVlaWatWqpa5du2rChAmqVq2adXxmZqZGjBhhnS/fo0cPzZ49W1WqVCkyRktouLR0m13HA5QqA9tf/hm/wZFRwJH6hMvYdWv9/bud8lNAuL8Gb3vMruMBSpsP2y+UJA3Y8IiDI4GjfB2+Rtu23To5yq4i6wrDMHTs2DEZhpFnhaTSgiILty2KLNyCRdYVt0N+osjC7YwiC7dakVXkdMGrWSwWeXt7F1csAADcEPITAMCZFLnwBQAAAADAfhRZAAAAAGAiiiwAAAAAMBFFFgAAAACYqMgi6/7779dnn31WYH9aWppcXFxMDQoAgKKQnwAAzqrIImv9+vXq06ePJkyYUOCY61gFHgAAU5CfAADOyq7pgnPmzNGMGTPUs2dPZWdn5+kvjc8jAQA4P/ITAMAZ2VVkPfjgg9qyZYt27dqlli1bat++fcUdFwAARSI/AQCckd0LX4SEhGjr1q2qVauW7rrrLq1du7Y44wIAwC7kJwCAs7mu1QU9PT21atUqDR06VBEREXrjjTeKKy4AAOxGfgIAOJOyRQ24dj67xWLRK6+8oqZNm2rIkCFat25dsQUHAEBByE8AAGdV5JWsglZmioyM1KZNm7Rz507TgwIAoCjkJwCAsyryStb69etVrVq1fPuaNm2q5ORkrVq1yvTAAAAoDPkJAOCsiiyy2rVrV2h/9erVNWDAANMCAgDAHuQnAICzuq6FLwAAAAAAhaPIAgAAAAATUWQBAAAAgIkosgAAAADARBRZAAAAAGAiiiwAAAAAMBFFFgAAAACYiCILAAAAAExEkQUAAAAAJqLIAgAAAAATUWQBAAAAgIkosgAAAADARBRZAAAAAGAiiiwAAAAAMFGJFlkbN25Ujx49VKNGDVksFsXHx9v0Dxw4UBaLxebVsmVLmzHnz5/XiBEj5OXlJXd3d/Xo0UOHDx+2GXPo0CF1795d7u7u8vLy0tNPP60LFy4U9+EBAG5R5CcAgJlKtMjKyspSo0aNNHPmTLm5ueU7pkOHDjp69Kj19dVXX9n0jxw5Up9++qkWL16sxMREnTp1St26dVNOTo4kKScnR127dtXp06eVmJioxYsX65NPPtHo0aOL/fgAALcm8hMAwExlS3JnERERioiIkHT5rGB+XF1d5efnl2/fyZMnNXfuXM2bN08dO3aUJH300UcKDAzU2rVr1blzZ61Zs0a7du3SwYMHVatWLUnSa6+9piFDhmjKlCmqXLmy+QcGALilkZ8AAGZyunuyNm3aJB8fHwUHB2vo0KFKT0+39iUnJ+vixYvq1KmTta1WrVoKCQlRUlKSJGnz5s0KCQmxJjBJ6ty5s86fP6/k5OSSOxAAQKlCfgIA2KtEr2QVpUuXLurVq5eCgoJ04MABvfTSS7rvvvuUnJwsV1dXpaamysXFRV5eXjaf8/X1VWpqqiQpNTVVvr6+Nv1eXl5ycXGxjrlWXFyc4uLiLr/JzDD/wAAAtzRnyE/ZGWeK4cgAAMXBqYqsyMhI6++NGzdWWFiYAgMDtWrVKvXq1avY9hsTE6OYmBhJkiU0vNj2AwC4NTlDfgoI9y+2/QAAzOV00wWvFhAQoJo1a2rv3r2SJD8/P+Xk5OjYsWM249LS0qzz5P38/JSWlmbTf+zYMeXk5BQ4lx4AgOtBfgIAFMapi6xjx47pr7/+kr//5bN3YWFhKleunBISEqxjDh8+rN27d6t169aSpFatWmn37t02y+YmJCTI1dVVYWFhJXsAAIBSifwEAChMiU4XzMrK0u+//y5Jys3N1aFDh7R9+3ZVq1ZN1apVU2xsrHr37i1/f38dOHBA48ePl4+Pj3r27ClJ8vT01ODBgzV27Fj5+PioevXqGjVqlJo0aaIOHTpIkjp16qTQ0FANGDBA06dP1/Hjx/Xcc89p6NChrNwEAMgX+QkAYKYSvZK1bds2NWvWTM2aNdPZs2c1YcIENWvWTC+//LJcXFy0c+dOPfjggwoODlZ0dLQaNGigzZs3q1KlStZtzJgxQz179lTfvn3Vpk0beXh4aOXKlXJxcZEkubi4aNWqVapYsaLatGmjvn37qnfv3nr99ddL8lABALcQ8hMAwEwWwzAMRwfhTCyh4dLSbY4OAyh5A9tf/hm/wZFRwJH6hMvYxd8/ZxUQ7q/B2x5zdBiAQ3zYfqEkacCGRxwcCRzl6/A12rbt1slRTn1PFgAAAADcaiiyAAAAAMBEFFkAAAAAYCKKLAAAAAAwEUUWAAAAAJiIIgsAAAAATESRBQAAAAAmosgCAAAAABNRZAEAAACAiSiyAAAAAMBEFFkAAAAAYCKKLAAAAAAwEUUWAAAAAJiIIgsAAAAATESRBQAAAAAmosgCAAAAABNRZAEAAACAiSiyAAAAAMBEFFkAAAAAYCKKLAAAAAAwEUUWAAAAAJiIIgsAAAAATESRBQAAAAAmosgCAAAAABNRZAEAAACAiSiyAAAAAMBEFFkAAAAAYCKKLAAAAAAwEUUWAAAAAJioRIusjRs3qkePHqpRo4YsFovi4+Nt+i0WS76v4cOHW8cMHDgwT3/Lli1ttnP+/HmNGDFCXl5ecnd3V48ePXT48OGSOEQAwC2I/AQAMFOJFllZWVlq1KiRZs6cKTc3tzz9R48etXmtXLlSktSnTx+bcR06dLAZ99VXX9n0jxw5Up9++qkWL16sxMREnTp1St26dVNOTk7xHRwA4JZFfgIAmKlsSe4sIiJCERERki6f8buWn5+fzfsVK1YoODhY7dq1s2l3dXXNM/aKkydPau7cuZo3b546duwoSfroo48UGBiotWvXqnPnziYcCQCgNCE/AQDM5LT3ZGVlZenjjz/W0KFD8/Rt2rRJPj4+Cg4O1tChQ5Wenm7tS05O1sWLF9WpUydrW61atRQSEqKkpKQSiR0AUHqRnwAARSnRK1nXY9GiRbpw4YKio6Nt2rt06aJevXopKChIBw4c0EsvvaT77rtPycnJcnV1VWpqqlxcXOTl5WXzOV9fX6Wmpua7r7i4OMXFxV1+k5lRLMcDACgdHJWfsjPOFM8BAQBM57RF1nvvvacHH3xQ3t7eNu2RkZHW3xs3bqywsDAFBgZq1apV6tWr1w3tKyYmRjExMZIkS2j4jQcNACj1HJWfAsL9bzxoAECJcsrpgtu3b9e2bdvynYpxrYCAANWsWVN79+6VdHnefE5Ojo4dO2YzLi0trcB58gAA2IP8BAD/j707j6/pWv8H/jkZnByZ0MwJCUKMNcRMRBtCG6HaahQRqqiaWmrqbSW0t7SG29YPvVREUXpLa/hWaWrOjVkQSdQYSkTEGCJBzvP7g+zryElk2Jk/79crr5y99tprr7X2PufJc84+O5QfZTLJWrx4MWrXro2uXbs+t25qaiouX74MZ+fH7/B5e3vD3NwckZGRSp1Lly4hISEBHTp0KLY+svOVDgAAIABJREFUExFRxcf4RERE+VGilwvevXsXZ86cAQDo9XpcvHgRR48eRY0aNVCrVi0AQHp6OlatWoVJkyZBo9Hk2D4sLAxvvPEGnJ2dkZiYiKlTp8LBwQF9+vQBANja2mLo0KGYNGkSHBwc8MILL2D8+PF48cUX8xUUiYio8mF8IiIiNZXoJ1mHDh1CixYt0KJFC9y/fx+hoaFo0aIFpk2bptT56aefcO/ePQwZMiTH9qampoiNjUXv3r1Rv359hISEwMvLC3v37oW1tbVS7+uvv0afPn0QFBSEjh07wsrKCps2bYKpqWmJjJOIiMoXxiciIlKTRkSktDtRlmgatwL+c6i0u0FU8gZ3efw7Ymdp9oJK01utIHF8/SurXFo5Y+ihnAkeUWXwQ5dVAIBBOweUck+otPze6g8cOlR+YlSZ/E4WERERERFReVVmb+FORCWsQfPS7gEREZFRTs0dSrsLRAXCJIuIHpvydWn3gIiIyCj/r7uVdheICoSXCxIREREREamISRYREREREZGKmGQRERERERGpiEkWERERERGRiphkERERERERqYhJFhERERERkYqYZBEREREREamISRYREREREZGKmGQRERERERGpiEkWERERERGRiphkERERERERqYhJFhERERERkYqYZBEREREREamISRYREREREZGKmGQRERERERGpiEkWERERERGRiphkERERERERqYhJFhERERERkYqYZBEREREREamISRYREREREZGKmGQRERERERGpiEkWERERERGRiphkERERERERqYhJFhERERERkYqYZBEREREREamISRYREREREZGKSjTJmjlzJlq3bg0bGxvY29sjMDAQJ06cMKgjIggLC4OLiwt0Oh26dOmCuLg4gzo3b95EcHAwbG1tYWtri+DgYNy6dcugTmxsLHx9faHT6eDq6ooZM2ZARIp9jEREVP4wPhERkZpKNMnauXMn3n//fURHR2P79u0wMzND165dcePGDaXOV199hblz52L+/Pk4ePAgHBwc0K1bN6SlpSl1+vfvjyNHjmDLli3YsmULjhw5guDgYGX9nTt30K1bNzg6OuLgwYP45ptvMHv2bMybN68kh0tEROUE4xMREalJI6X49tndu3dha2uL9evXIzAwECICFxcXjB49Gv/4xz8AAPfv34eDgwPmzJmDESNGICEhAY0aNUJUVBQ6duwIAIiKioKPjw9OnjwJLy8vLFq0CJMnT8bVq1eh0+kAAJ9//jkWLVqES5cuQaPR5NonTeNWwH8OFf/giYjKmrdaQeL4+geUzfjk0soZQw8NKf7BExGVQb+3+gOHDpWfGFWq38lKS0uDXq9H9erVAQDnz59HcnIy/P39lTo6nQ6dO3dGdHQ0AGDv3r2wsrJChw4dlDodO3aEpaWlQR0fHx8lgAFA9+7dkZSUhMTExBIYGRERlWeMT0REVBRmpbnzcePGoXnz5mjfvj0AIDk5GQDg6OhoUM/R0RGXL19W6tjb2xu826fRaODg4KBsn5ycDDc3txxtZK+rXbu2wbrFixdj8eLFAADt2RNoEtJKrSFWOteuXYO9vX1pd6Pc4vwVDeevaE5eOFnaXSgzymJ8unHiJn5v9YdaQ6x0+PpQNJy/ouMcFs3Jk+UrRpVakjV+/HhERUUhKioKpqampdUNAMDw4cMxfPhwAECrVq3K1UeRZQ3nr2g4f0XD+SuaVq34BhPA+FRRcf6KhvNXdJzDoilvMapULhf88MMPsXr1amzfvh116tRRyp2cnAAAV69eNah/9epVZZ2TkxOuXbtmcCcmEUFKSopBHWNtPL0PIiKiZzE+ERGRGko8yRo3bpwSwBo0aGCwrnbt2nByckJkZKRSlpGRgT179ijXuLdv3x53797F3r17lTp79+7FvXv3DOrs2bMHGRkZSp3IyEi4uLjAw8OjGEdHRETlFeMTERGpxTQsLCyspHY2atQoLF++HD///DNq1aqFu3fv4u7duwCAKlWqQKPRICsrC7NmzYKXlxeysrIwfvx4XLlyBYsXL4ZWq4W9vT3279+PH3/8ES1atMClS5cwYsQItGnTBmPGjAEA1K9fH//+979x9OhRNGjQAFFRUfjoo48wdepUgy8k58bb27tY56Gi4/wVDeevaDh/RVNZ54/xqXLg/BUN56/oOIdFU67mT0oQAKM/oaGhSh29Xi+hoaHi5OQkWq1WOnfuLLGxsQbt3LhxQwYMGCDW1tZibW0tAwYMkJs3bxrUOX78uPj4+IhWqxUnJycJCwsTvV5fEsMkIqJyhvGJiIjUVKr/J4uIiIiIiKiiKdX/k0VERERERFTRMMkiIiIiIiJSEZOsJxYuXIjatWvDwsIC3t7e2LNnT2l3qUwICwuDRqMx+Hn6NsMigrCwMLi4uECn06FLly6Ii4szaOPmzZsIDg6Gra0tbG1tERwcjFu3bpX0UErE7t270atXL7i6ukKj0SAiIsJgvVrzFRsbC19fX+h0Ori6umLGjBmoCFf+Pm/+Bg8enON8bNeunUGdzMxMjBkzBnZ2drC0tESvXr1w6dIlgzoXL15EYGAgLC0tYWdnh7Fjx+LBgwfFPbxiN3PmTLRu3Ro2Njawt7dHYGAgTpw4YVCH52D5xBiVE+NTwTA+FR1jVOFVyvhUOl8FK1vWrFkjZmZmsnjxYomPj5fRo0eLpaWlXLhwobS7VupCQ0PFy8tLrly5ovykpKQo62fNmiVWVlaydu1aiY2Nlb59+4qzs7PcuXNHqdOjRw9p1KiRREdHS3R0tDRq1Eh69uxZGsMpdr/99ptMnTpVfv75Z9HpdLJs2TKD9WrM1+3bt8XR0VH69u0rsbGx8vPPP4uVlZXMmTOnpIZZbJ43fyEhIdK1a1eD8/H69esGdd577z1xdnaWP/74Qw4fPiy+vr7SrFkzefTokYiIPHr0SJo0aSK+vr5y+PBh+eOPP8TZ2VlGjx5dUsMsNv7+/hIeHi6xsbFy/Phxee2118TR0dFgjngOlj+MUcYxPhUM41PRMUYVXmWMT0yyRKRNmzby7rvvGpR5enrKlClTSqlHZUdoaKg0btzY6Dq9Xi9OTk7y+eefK2Xp6eliZWUl3333nYiIxMfHCwCJiopS6uzZs0cAyMmTJ4u386XM0tLS4AVYrflauHChWFtbS3p6ulLns88+ExcXlwp1h7Jn50/kcQALCAjIdZtbt26Jubm5rFy5Uim7ePGiaDQa2bJli4iIbN68WTQajVy8eFGps2LFCtFqtXL79m11B1HK0tLSxMTERDZu3CgiPAfLK8Yo4xifCo/xqegYo4qmMsSnSn+54IMHD3D48GH4+/sblPv7+yM6OrqUelW2nDt3Di4uLqhduzb69euHc+fOAQDOnz+P5ORkg7nT6XTo3LmzMnd79+6FlZWVwf9/6dixIywtLSvd/Ko1X3v37oWPjw90Op1Sp3v37khKSkJiYmLJDKYURUVFwcHBAfXr18ewYcOQkpKirDt8+DAePnxoMMc1a9ZEw4YNDeavYcOGqFmzplKne/fuyMzMxOHDh0tuICUgLS0Ner0e1atXB8BzsDxijMob45M6+NqgHsao/KkM8anSJ1mpqanIysqCo6OjQbmjoyOSk5NLqVdlR9u2bREREYEtW7ZgyZIlSE5ORocOHXD9+nVlfvKau+TkZNjb20Oj0SjrNRoNHBwcKt38qjVfycnJRtt4eh8VVY8ePfDDDz9g27ZtmDt3Lg4cOICXX34ZmZmZAB6P39TUFHZ2dgbbPTvHz86fnZ0dTE1NK9z8jRs3Ds2bN0f79u0B8Bwsjxijcsf4pB6+NqiDMSr/KkN8MivRvVG588orrxgst2vXDnXq1MHy5ctzfJmTqLj169dPedy0aVN4e3vD3d0dv/32G15//fVS7FnZM378eERFRSEqKgqmpqal3R0i1TE+UVnDGJU/lSU+VfpPsrLfHbh69apB+dWrVw3uUkSPWVlZoXHjxjh9+rQyP3nNnZOTE65du2ZwVxcRQUpKSqWbX7Xmy8nJyWgbT++jsnBxcYGbmxtOnz4N4PH4s7KykJqaalDv2Tl+dv6yPy2oKPP34YcfYvXq1di+fTvq1KmjlPMcLH8Yo/KP8anw+NpQPBijcqpM8anSJ1lVqlSBt7c3IiMjDcojIyMNrvmkxzIyMnDy5Ek4Ozujdu3acHJyMpi7jIwM7NmzR5m79u3b4+7du9i7d69SZ+/evbh3716lm1+15qt9+/bYs2cPMjIylDqRkZFwcXGBh4dHyQymjEhNTcXly5fh7OwMAPD29oa5ubnBHF+6dAkJCQkG85eQkGBwy9zIyEhotVp4e3uX7ACKwbhx45QA1qBBA4N1PAfLH8ao/GN8Kjy+NhQPxihDlS4+lehtNsqoNWvWiLm5uSxZskTi4+Nl7NixYmlpKYmJiaXdtVI3YcIE2blzp5w7d0727dsnAQEBYm1trczNrFmzxMbGRtatWyexsbESFBRk9HabTZo0UW632aRJkwp7i9y0tDSJiYmRmJgY0el0Mn36dImJiVFutazGfN26dUscHR0lKChIYmNjZd26dWJtbV0hbpGb1/ylpaXJhAkTJDo6Ws6fPy87duyQdu3aiaurq8H8vffee+Lq6iqRkZFy5MgR6dKli9Hb47700kty5MgRiYyMFBcXl3J/e1wRkffff1+sra1l27ZtBrcQTktLU+rwHCx/GKOMY3wqGManomOMKrzKGJ+YZD2xYMECcXd3lypVqkjLli1l165dpd2lMiH7BDc3NxcXFxd5/fXXJS4uTlmv1+slNDRUnJycRKvVSufOnSU2NtagjRs3bsiAAQPE2tparK2tZcCAAXLz5s2SHkqJ2LFjhwDI8RMSEiIi6s3X8ePHxcfHR7RarTg5OUlYWFiFuD1uXvOXnp4u/v7+Ym9vL+bm5lKrVi0JCQkxuM2tiEhGRoaMHj1aatSoITqdTnr27JmjzoULFyQgIEB0Op3UqFFDxowZIxkZGSU51GJhbO4ASGhoqFKH52D5xBiVE+NTwTA+FR1jVOFVxvikEakg/4abiIiIiIioDKj038kiIiIiIiJSE5MsIiIiIiIiFTHJIiIiIiIiUhGTLCIiIiIiIhUxySIiIiIiIlIRkywiIiIiIiIVMckiKkFr166FRqNRliMiImBlZVWKPcqfIUOGYMaMGUVqIzY2Fq6urrh3755KvSIiIjUxRjFGkXqYZFGZptfr0blzZwQGBhqUp6enw8vLC++9916e26elpeHTTz9Fo0aNoNPp4OjoiC5dumD16tXQ6/XF2fV8CQoKwrlz51RvV6PRYO3ataq0FRsbi/Xr1+ODDz4oUjtNmzZFu3btMG/ePFX6RURU2hijCocxiioDJllUppmYmCAiIgI7duxAeHi4Uj558mRkZWVh7ty5uW5769YttG/fHuHh4Zg4cSIOHTqEqKgohISE4LPPPsPFixeLrd8PHjzIVz2dTgcHB4di64ca5s+fjzfeeAM2NjZFbmvIkCFYtGgRHj16pELPiIhKF2NU6WOMojJLiMqBRYsWiY2NjSQmJsqff/4ppqamsmfPnjy3GTlypFStWlX+/vvvHOvu378v9+/fFxGRGzduyKBBg6RatWpiYWEhfn5+cuLECYP669atkyZNmkiVKlXEzc1NPv/8c9Hr9cp6d3d3CQ0NlSFDhoitra28+eabIiKyfPlyqVWrluh0OgkICJD/9//+nzz9tFu2bJlYWloqy6GhodK4cWNZvXq11KlTR6ysrKR3795y7do1pc6BAwekW7du8sILL4i1tbV07NhRoqOjDfoCQPlxd3dX1m3cuFFatmwpWq1WPDw85OOPP5bMzMxc5/DRo0dia2sr69evNyh3d3eX6dOnS0hIiFhZWYmbm5usWbNGbt68KUFBQWJpaSmenp6ydetWg+0yMzNFq9VKZGRkrvskIipvGKMYo4iexSSLyg1/f3/x8fGRmjVrysSJE/Osm5WVJdWrV5dhw4Y9t91evXqJl5eX7Nq1S44fPy6BgYHi5uYm6enpIiJy6NAhMTExkWnTpslff/0lK1euFEtLS/n222+VNtzd3cXa2lq+/PJLOX36tJw6dUr27dsnGo1GPv/8c/nrr7/ku+++kxo1ajw3gFlaWsprr70mx44dk+joaKlVq5YMHz5cqbNt2zb54YcfJD4+XhISEmTUqFFSrVo1SU1NFRGRlJQUASBLliyRK1euSEpKioiIbNmyRaytrSU8PFzOnDkj27dvl/r168uECRNynZsjR44IALl06ZJBubu7u1SvXl0WLFggp06dkvHjx4tWq5VXXnlFli9fLqdPn5Z33nlH7O3tlT8UsrVt21Y++eST5x4XIqLyhDHqMcYooseYZFG5ce7cOdFoNOLp6SkZGRl51r169aoAkHnz5uVZ79SpUwJAdu3apZTdunVLbGxsZMmSJSIi0r9/f3nppZcMtgsNDRVXV1dl2d3dXXr27GlQ5+2335auXbsalA0dOvS5AUyr1cqtW7eUss8//1zq1q2b6xj0er04OTnJihUrlDIA8vPPPxvU8/HxkRkzZhiU/frrr2JpaWnwjuez6zUajWRlZRmUu7u7S79+/ZTltLQ0ASBjxoxRys6fPy8A5ODBgwbb9unTRwYOHJjreIiIyiPGKOMYo6iy4neyqNwIDw+HTqfDpUuXcP78+Tzriki+2kxISICJiQnat2+vlNna2qJp06aIj49X6nTs2NFgu06dOuHy5cu4c+eOUtaqVascbT/dLoAcy8a4u7vD1tZWWXZxcUFKSoqynJKSghEjRqB+/fqwtbWFtbU1UlJSnnv9/uHDh/HPf/4TVlZWyk///v1x7949JCcnG93m/v37MDc3h4lJzpeKF198UXlsZWWFqlWromnTpkqZo6Oj0t+n6XQ63L9/P8++EhGVN4xRjzFGET1mVtodIMqPgwcPYtasWdi4cSMWLVqEkJAQREdHw9TU1Gh9e3t7VKtWDQkJCYXe59O3sc1PHUtLy0Lv62nm5uY59vH0XaZCQkJw9epV/Otf/4KHhwe0Wi38/Pye+0VmvV6P0NBQ9O3bN8c6e3t7o9vY2dnhwYMHSE9PR9WqVZ/bz6fLsufm2Ttk3bhxAx4eHnn2lYioPGGMYowiehY/yaIyLyMjA4MGDcLgwYPxyiuvYPHixThz5gy++uqrXLcxMTFBv379sGrVKly6dMlomxkZGWjYsCH0ej327t2rrLtz5w5iY2PRqFEjAEDDhg3x3//+12D7qKgouLm5wdraOtc+NGzYEPv27TMoe3a5MKKiojBmzBgEBASgcePGsLa2xpUrVwzqmJubIysry6CsZcuWOHnyJDw9PXP8mJkZf7+lefPmAKC8Y6qGEydOoGXLlqq1R0RUmhijDDFGET3GJIvKvKlTpyIjI0P53xVOTk5YsGABwsLCEBcXl+t2//znP1GrVi20bdsWy5YtQ1xcHM6cOYMVK1bA29sbycnJqFevHnr37o0RI0Zgz549iI2NxcCBA2FjY4P+/fsDACZMmIBdu3YhLCwMp06dwqpVqzB37lxMmjQpz36PHTsWf/75J2bOnInTp09jyZIl+PXXX4s8H/Xr18fKlSsRHx+PgwcPol+/fqhSpYpBHQ8PD2zbtg3Jycm4efMmAGDatGn48ccfMW3aNJw4cQInT57E2rVr8xyHvb09WrZsiaioqCL3GwASExNx+fJl+Pv7q9IeEVFpY4wyxBhF9ERpfymMKC+7du0SU1NT2bFjR451b775pnh7e8vDhw9z3f7WrVvy8ccfi5eXl2i1WrG3txdfX19ZvXq18kXZgtwe19zcPNfb486ePTvH/sPDw6VmzZpiYWEhPXr0kPnz5+fr9rhPe7bO0aNHpU2bNmJhYSF16tSRH374QRo3biyhoaFKnY0bN4qnp6eYmZkZ3B5369at0qlTJ9HpdGJtbS3e3t4yf/78XOdPROS7776TVq1aGZQZG6+lpaUsW7ZMWb5//74AkE2bNillX3zxhXTv3j3P/RERlReMUYxRRLnRiOTz25dEVCllZGSgQYMGWLFiBXx8fArdTmZmJurVq4fVq1fn+JI2ERFRYTBGUVnFywWJKE8WFhb44YcfcOPGjSK1c+HCBfzjH/9g8CIiItUwRlFZxU+yiIiIiIiIVMRPsoiIiIiIiFTEJIuIiIiIiEhFTLKIiIiIiIhUxCSLiIiIiIhIRUyyiIiIiIiIVMQki4iIiIiISEVMsoiIiIiIiFTEJIuIiIiIiEhFTLKIiIiIiIhUxCSLiIiIiIhIRUyyiIiIiIiIVMQki4iIiIiISEVMsoiIiIiIiFTEJIuIiIiIiEhFTLKIiIiIiIhUxCSLiIiIiIhIRUyyiIiIiIiIVMQki4iIiIiISEVMsoiIiIiIiFTEJIuIiIiIiEhFTLKIiIiIiIhUxCSLiIiIiIhIRUyyiIiIiIiIVMQki4iIiIiISEVMsoiIiIiIiFTEJIuIiIiIiEhFTLKIiIiIiIhUxCSLiIiIiIhIRUyyiIiIiIiIVGRW2h0oLI3GE8B9AJrskqd+kM/Hua0r7HbG6uU2gDx2CyPrni0r6Lr8dN3Y9kXtU0H2+9w+yZPH8nhRWX78WKMRpb4GAs2Teo+LxOhvQJ7arWH95y3n3W7h21C7vcK0YWxeCtWePNlGAMiT39nkqd/Plhe0nrHy57WR2+/81CtM/wrS9/zWk6cWsx/LU4/z6OKz01SQ7ha0vRe7d8eWLVtQWfwvPgEFj035edHNqw0UsJ6xAeTRxLPNP/u4sOvyE9+Ko09FiasGbctTZc/EI41hzMrva6kS64zUzet3bu2p0Yba7RWmjbzmpcDtiSjx6fEy/qcsxouCtPFs+fPaKEjfC9AP5eHTsUny7np+41RuQyxoG8UVo8ptkgWkAxgFwPzJstmTx9lDyuvx8+o922Zuj59X7xkaAKZGdmv2pPzpJrLXGSsv7jaeLs9r3dPtPbv89H6NrTPWtrG+mwEwe/I0MMsCzLJgYvro8aJ5FkzNsh7vyiwLpmZPys2yYGqSBVM8WYcsmOLJOmQ9Wc5e9whmBvUMtzG+Lmf5s+3ntk3++pRXe3n3t2BjNNanwo3R6LqsLJg+elKepYfpI+DJoYMmC3iyyePf2Y+zclkuaL3nbZdbvSwj26hZL6/tnlfP2Jie1Hv4CHgy1XiY9dTjR4+rPXyy2cOnmsh+nL3ueY/zW+/Z9gFgU2oqKpfs+ASoH1fUinXPyE4Q1IgXz4tFareR1zZ5xRxj+zX2OK/YBzyOUU/iEwCYmD6CmfnTsSn78SMlPj3ubknGFbViXdHjSsnHOiPtZT15/OiREp8eLz+JT0DxxQEgZwxTIQ6oUk+t2PnUdg+fPH706HF8Uh4/MowXBY0rzz5+Xr3c2geKL0bxckEiIiIiIiIVMckiIiIiIiJSEZMsIiIiIiIiFTHJIiIiIiIiUhGTLCIiIiIiIhUxySIiIiIiIlIRkywiIiIiIiIVMckiIiIiIiJSEZMsIiIiIiIiFTHJIiIiIiIiUhGTLCIiIiIiIhUxySIiIiIiIlIRkywiIiIiIiIVMckiIiIiIiJSEZMsIiIiIiIiFTHJIiIiIiIiUpFZaXegsBo1qg6d7tfS7kaJunbtGuzt7QF5UvDgyU8FpYy3BGU987uklcaYDZk+87v4lf6YS9Zzx2uK/01/lZLoUfG7f/9+aXehRFXG+GRA8L/YVMwxqrK9fmShMscnwPAFsmIrG/NtRPYhqCDxCSi+GFVukyydTodDhw6VdjdKVKtWrSrVmCvbeAGOuTKobOMFHo+5MqmM8am0VMbnU2nhXJcsznfJKa4YxcsFiYiIiIiIVMQki4iIiIiISEWmYWFhYaXdicLy9vYu7S6UuMo25so2XoBjrgwq23iByjfmyjbe0sS5Ljmc65LF+S45xTHXGhGR51cjIiIiIiKi/ODlgkRERERERCpikkVERERERKQiJllEREREREQqKvUka/fu3ejVqxdcXV2h0WgQERGR721Pnz4Na2trWFlZ5Vi3a9cueHt7w8LCAnXq1MF3332nYq+LpjjGvHPnTmg0mhw/J0+eVLn3BVfQ8SYmJhody5YtWwzqVaRjnJ8xl+VjDBTuvBYRfP3112jQoAG0Wi2cnZ0xZcoUgzpl9TgXx3gr2jEOCwszOh6NRoOUlBSlXmxsLHx9faHT6eDq6ooZM2agLH1deOHChahduzYsLCzg7e2NPXv25Fl/wYIFaNiwIXQ6Hby8vPDDDz/kqLNu3To0atQIWq0WjRo1wq+/VuJ/XvwUtec6IiLC6PmXkZFRnMMo8wrz+pWf5ynP65yKY655XhtX0LnOyMjA4MGD8eKLL8Lc3BxdunQxWq+wf4eUepJ19+5dNGnSBN988w10Ol2+t3vw4AH69euHzp0751h3/vx5vPrqq+jQoQNiYmIwdepUjBkzBuvWrVOz64VWHGPOFhcXhytXrig/9erVU6PLRVLY8W7ZssVgLC+//LKyrqIe47zGnK0sHmOgcGOeMGECFi5ciC+//BIJCQnYvHmzwfldlo9zcYw3W0U5xh999JHBOK5cuQJfX1906dIFDg4OAIA7d+6gW7ducHR0xMGDB/HNN99g9uzZmDdvXnEPJ19++uknjBs3Dh9//DFiYmLQoUMHvPLKK7h48aLR+osWLcLkyZMxbdo0xMXFYfr06Rg1ahQ2bdqk1Nm7dy+CgoIwYMAAHD16FAMGDEDfvn2xf//+khpWmVQccw0AVatWzXEeWlhYlMSQyqyCPpfz8zzleW1cccw1wPPamILOdVZWFiwsLDB69GgEBAQYrVOkv0OkDLG0tJRly5blq+4HH3wggwcPlmXLlomlpaXBukmTJomnp6dB2dChQ6Vdu3ZqdVU1ao15x44dAkCuXbtWDL1UT37Ge/78eQEgBw8ezLVORTvG+RlzeTnGIvleJJmhAAAgAElEQVQb88mTJ8XMzEzi4+NzrVNejrNa461ox/hZFy9eFBMTE1m1apVStnDhQrG2tpb09HSl7LPPPhMXFxfR6/VqdbfQ2rRpI++++65Bmaenp0yZMsVo/fbt28sHH3xgUDZ+/Hjp2LGjsvzWW29J165dDer4+flJv379VOp1+VQcc20sXpKh/DyX8/M85Xn9fGrNNc/r5ytojBo1apT4+vrmKC/K3yGl/klWYfz222/4v//7P8yfP9/o+r1798Lf39+grHv37jh06BAePnxYEl1U3fPGnK1Vq1ZwdnaGn58fduzYUUK9Kx6vv/46HBwc0LFjR6xdu9ZgXUU8xkDeY85WUY7xhg0bUKdOHWzZsgV16tSBh4cHQkJCDC4jq0jHOT/jzVZRjvGzli5diurVq+ONN95Qyvbu3QsfHx+Ddx27d++OpKQkJCYmlkIv/+fBgwc4fPhwjnPQ398f0dHRRrfJzMzM8W6yTqfDgQMHlHM2t/M6tzYrg+KaawC4f/8+3N3d4ebmhp49eyImJkb9AVRw+Xme8rxWR35fE3lel4yi/B1S7pKspKQkDBs2DCtXrjT6XSwASE5OhqOjo0GZo6MjHj16hNTU1JLopqryM2ZnZ2csWrQI69atwy+//AIvLy/4+fk993r2ssjKygpz5szBf/7zH2zevBl+fn4ICgrCypUrlToV7RjnZ8wV6RgDwLlz53DhwgWsWbMGERERWLFiBU6ePInAwEDo9XoAFes452e8Fe0YPy0rKwvh4eEIDg6GVqtVynM7xtnrSlNqaiqysrKM9i+3vnXv3h3h4eE4ePAgRASHDh3C999/j4cPHyrnbG5jLu3xlqbimmsvLy+Eh4djw4YNWL16NSwsLNCxY0ecPn262MdUkeTnecrzWh35mWue1yWnKH+HmBVnx4pDcHAwRo4cibZt25Z2V0pMfsbs5eUFLy8vZbl9+/ZITEzE7Nmz4ePjUxLdVI2dnR0mTJigLLdq1Qqpqan46quvMHDgwFLsWfHJz5gr0jEGAL1ej8zMTKxYsQL169cHAKxYsQJeXl44ePBghXuO52e8Fe0YP23Lli34+++/MWzYsNLuSrH69NNPkZycjA4dOkBE4OjoiJCQEHz11VcwMSl372uWafmZ6/bt26N9+/bKNh06dEDz5s0xf/58fPvtt6XVdaIi4XldPpS7V/zt27dj+vTpMDMzg5mZGYYOHYp79+7BzMwMixcvBgA4OTnh6tWrBttdvXoVZmZmsLOzK41uF0l+xmxM27ZtK8y7Gs+OpaIdY2Pyc/zK8zF2dnaGmZmZknAAQL169WBqaqp80b0iHef8jNeY8nyMn7Z48WJ06NABjRo1MijP7RhnrytNdnZ2MDU1Ndq/3Pqm0+kQHh6O9PR0JCYm4uLFi/Dw8IC1tTXs7e0B5D7m0h5vaSquuX6WqakpWrVqVSGeUyUpP89TntfqKMxrIs/r4lOUv0PKXZIVGxuLo0ePKj8zZsyATqfD0aNH0bdvXwCPM/zIyEiD7SIjI9GqVSuYm5uXRreLJD9jNubo0aNwdnYuwZ4Wn2fHUtGOsTH5OX7l+Rh37NgRjx49wtmzZ5Wyc+fOISsrC+7u7gAq1nHOz3iNKc/HOFtSUhJ+++03o59itW/fHnv27DG49XBkZCRcXFzg4eFRgr3MqUqVKvD29jZ6Dnbo0CHPbc3NzeHm5gZTU1OsWbMGPXv2NPh0pTBtVmTFNdfPEhEcP3683D+nSlp+nqc8r9VRmNdEntfFp0h/h+T7thvFJC0tTWJiYiQmJkZ0Op1Mnz5dYmJi5MKFCyIiMmXKFHn55Zdz3d7YHVbOnTsnVatWlXHjxkl8fLwsWbJEzM3NZe3atcU6lvwqjjH/61//kl9//VVOnTolJ06ckClTpggAWbduXbGOJT8KOt6IiAhZtWqVxMfHy8mTJ2X27Nlibm4u8+bNU+pUtGOcnzGX5WMsUvAxZ2VlScuWLaVz585y5MgROXLkiHTu3Fnatm0rWVlZIlK2j3NxjLeiHeNsn332mdjY2Mi9e/dyrLt165Y4OjpKUFCQxMbGyrp168Ta2lrmzJlT7OPJjzVr1oi5ubksWbJE4uPjZezYsWJpaSmJiYkiIhIcHCzBwcFK/b/++kt++OEHOXXqlOzfv1+CgoKkRo0acv78eaXOf//7XzE1NZWZM2dKQkKCfPHFF2JmZib79u0r6eGVKcUx12FhYbJlyxY5e/asxMTEyJAhQ8TMzEz2799f0sMrUwr6XM7P85TntXHFMdc8r40rTIyKi4uTmJgYCQoKEm9vb2X7bEX5O6TUk6zsWxY/+xMSEiIiIiEhIeLu7p7r9rndxnLnzp3SokULqVKlinh4eMiiRYuKaQQFVxxj/vLLL8XT01MsLCykevXq0qlTJ/ntt9+KcRT5V9DxRkRESMOGDaVq1apibW0t3t7esmLFihztVqRjnJ8xl+VjLFK48zopKUnefPNNsbKyEnt7e+nfv78kJycb1Cmrx7k4xlsRj7FerxcPDw8ZOXJkru0eP35cfHx8RKvVipOTk4SFhZWJ27dnW7Bggbi7u0uVKlWkZcuWsmvXLmWdr6+vwW1/4+PjpXnz5qLT6cTGxkZ69+4tJ0+ezNHmzz//LF5eXmJubi4NGjQoM4l0aVN7rj/44AOpVauWVKlSRezt7cXf31+io6NLajhlVmGey/l5nvK8zqk45prntXGFmWt3d3ej2zytsH+HaESe+XfdREREREREVGjl7jtZREREREREZRmTLCIiIiIiIhUxySIiIiIiIlIRkywiIiIiIiIVMckiIiIiIiJSEZMsIiIiIiIiFTHJogorLCwMTZo0Ua29iIgIWFlZqdbes6ZPn4533nmn2NovS1JSUmBvb49Lly6VdleIiMo9Dw8PzJkzp1javnnzJhwdHXH27FnV2oyJiYGtrS2qV6+OY8eOqdbuzp07odFokJqaWug2YmNj4erqinv37qnWL6qcmGRRmdOrVy/4+fkZXZeQkACNRoM//vijhHsFBAUF4dy5c8qymklcSkoK5s6di08++USV9h48eABbW1scPXpUlfbU5uDggEGDBiE0NLS0u0JEVGCDBw+GRqOBRqOBubk5HBwc8NJLL2HBggV4+PBhiffn4MGDeP/995VljUaDtWvXqtL2F198gVdffRV169ZVpb1Tp06he/fu6Nu3LwIDA9G9e3ecOXMmR71PP/0UDRo0gKWlJapXrw4/Pz9ER0cbbbNBgwZYv369Kv1r2rQp2rVrh3nz5qnSHlVeTLKozBk6dCh27NiBxMTEHOuWLl0Kd3d3dO3atcT7pdPp4ODgUCxtf//992jTpg3q1KmjSns7duxA9erV0bx5c1XaKw5DhgzBqlWrcOPGjdLuChFRgXXt2hVXrlxBYmIi/vjjDwQGBiI0NBQ+Pj4l/imIvb09qlatqnq76enp+P777zF06FBV2vv777/RrVs39O/fH0uWLMHy5cvx+uuvo1u3brh8+bJBXS8vLyxYsACxsbGIiopC7dq10aNHD1y9etWg3l9//YWLFy/C399flT4Cj+PTokWL8OjRI9XapMqHSRaVOQEBAXB0dMSyZcsMyh8+fIgVK1bgnXfegYmJCeLj4xEQEABra2s4ODjg7bffRnJycq7t6vV6fPbZZ6hZsya0Wi2aNm2KDRs2GNRJSkrCgAED8MILL6Bq1apo3rw5duzYAcDwcsGIiAhMnz4dcXFxyruZEREReOedd9CzZ88c+61Vq1ae74r9+OOPCAwMNCjr0qULRo4ciQkTJqBGjRqwt7fHN998g8zMTIwaNQrVqlVDrVq1sGLFihztbdiwAb179zbo9++//44GDRqgatWq6NWrF27fvo21a9eiXr16sLW1RXBwMO7fv6+0sXv3brRr1w5WVlawtbVFmzZtcOLECWV9dHQ0fH19UbVqVbi6umLkyJG4c+eOsl5EMHfuXNSrVw9arRZubm6YOnWqsr5JkyZwcXHBL7/8kuu8EBGVVVqtFk5OTnB1dUXz5s0xfvx47Ny5E0eOHMFXX32l1Hvw4AEmT54MNzc3VK1aFa1bt8bWrVuV9dmXuG3btg1t27ZF1apV0apVKxw5ckSpc/v2bQQHB8PBwQEWFhaoU6cOvv76a2X905cLenh4AAD69u0LjUYDDw8PJCYmwsTEBIcOHTIYw5IlS2BnZ4cHDx4YHePmzZuh0WjQsWPHHP39/fff4e3tDZ1OBx8fH1y6dAm7du1Cs2bNYGVlhZ49e+L69evKdqmpqfD398egQYPw9ddfK7Fz4cKFCAoKgr+/v0H9gQMHws/PD3Xq1EHjxo0xb948pKWl5bhCY8OGDejWrZvRJDMzMxN9+vRBy5YtkZKSAgDYv38/WrZsCQsLC7Ro0UIZ486dO5Xt/P39cePGDYMyogITojJo8uTJUqtWLcnKylLK1q1bJyYmJnLx4kVJSkqSF154QSZNmiTx8fFy7Ngx6dmzp7Rp00bZJjQ0VBo3bqxsP2/ePLG2tpZVq1bJX3/9JZ9++qmYmJhITEyMiIjcvXtXPD09pUOHDrJ79245c+aMrFu3TrZv3y4iIsuWLRNLS0sREUlPT5cJEyaIl5eXXLlyRa5cuSLp6ekSHR0tpqamkpSUpOx3y5YtYm5uLikpKUbHev36ddFoNBIVFWVQ7uvrK9bW1hIaGiqnTp2SOXPmCADp0aOHfP3113L69Gn55JNPpEqVKgb70+v14urqKtu2bVP6bWZmJn5+fnLo0CGJjo4WZ2dn8fPzk549e8qxY8dk+/btUq1aNZkzZ46IiDx8+FCqVasmEyZMkDNnzkhCQoKsWrVK4uPjRUTk+PHjYmlpKXPmzJFTp07Jvn37pF27dvLGG28o/ZgyZYrY2trK0qVL5fTp0xIdHS0LFiwwGGNQUJAMHDjwuecDEVFZEhISIgEBAUbXBQYGGsSe/v37S9u2bWXXrl1y9uxZmT9/vpibm8vRo0dFRGTHjh0CQFq3bi3bt2+XhIQE8ff3lwYNGoherxcRkdGjR0uzZs1k//79kpiYKDt27JD//Oc/yj7c3d1l9uzZIiKSkpIiAGTJkiVy5coVJfb4+/vLyJEjDfrarl07+eCDD3Id59ixY6Vr164GZU/3d/fu3XLs2DFp3LixdOjQQV5++WXZt2+fHDx4UDw8PGT06NH5ndI8ZWZmyuzZs8XGxkaSk5MN1rVv316WLl1q0Ldr167J7du3pUuXLtK5c2e5ffu2iIikpaWJnZ2dvP3223LixAn5448/pFGjRgJAduzYYdBu27Zt5ZNPPlGl/1Q5McmiMunUqVMCQLZu3aqUvfrqq9KjRw8REfn000/l5ZdfNtjmxo0bAkD2798vIjmTLBcXF5k+fbrBNr6+vjJgwAAREVm8eLFYWVnJtWvXjPbp6STLWPvZGjduLDNnzlSW33rrLYPk41kxMTECQM6dO5ejb+3atVOW9Xq92NnZSWBgoFL24MEDMTc3l59//lkpO3DggFSvXl0ePnyo9BuAnDx5UqkzYcIEMTExMRjr0380XL9+XQDIzp07jfY5ODhY3nnnHaPjuHr1qqSlpYlWq5VFixblOm4RkQ8//FA6deqUZx0iorImryRr8uTJotPpRETkzJkzotFo5MKFCwZ1evfurSQ82YnBli1blPVRUVECQP7++28ReZy4DRkyJNf+PJ1kiYgAMIgLIiI///yzVKtWTe7fvy8iIvHx8QJAYmNjc223d+/eMmjQIIMyY/2dP3++AJDDhw8rZbnFyILYtGmTWFpaikajERcXFyW+Z0tOThYzMzMlkczuW1xcnLRs2VICAwOV8YqIfPfdd1K9enVJT09XylatWmU0yerTpw/fBKQi4eWCVCbVq1cPvr6+CA8PB/D4Mr6tW7cq14UfPnwYu3fvhpWVlfJTs2ZNADB6B6Q7d+4gKSnJ4JIHAOjUqRPi4+MBPL7b0Ysvvgg7O7si9X3YsGHKpY43btzAhg0b8ryePfsSPQsLixzrXnzxReWxRqOBg4MDmjZtqpSZm5ujevXqymUQwONLJwICAmBmZqaUabVaeHl5KcuOjo5wcnIyGKujo6PSTo0aNTB48GB0794dAQEBmDdvHi5evKjUPXz4MFauXGkw/9lze/bsWcTHxyMzMzPXG5hk0+l0BpcoEhGVdyICjUYDADhy5AhEBI0aNTJ4vfztt99yxKqnX+9dXFwAQHlNHjlyJH766Sc0a9YMH330EXbt2lXgfvXu3RtVqlRRLtEODw9HmzZt8ryB0/37943Gpmf76+joCAAG8enpmFJYL730Eo4ePYro6Gj06NEDb731Fq5cuaKs37RpE9q1awd7e3uD7bp37w43Nzf88ssvBv0/efIkmjRpAp1Op5S1bdvW6L4Zn6iomGRRmTV06FCsX78eN27cQEREBGrUqKF8z0iv1yMgIABHjx41+Dl9+nSO70Q9T3YwVEtwcDAuXLiAqKgorFq1Cvb29ujevXuu9bMTnZs3b+ZYZ25unqOvxsr0er2yvH79emWesj2dcOW3nWXLlmH//v3o3LkzNm7cCC8vL+V7BHq9Hu+++67B3B87dgynT58u0M02bty4kSM4EhGVZ/Hx8cpNjPR6PTQaDQ4ePGjwepmQkKC8iZjt6dfk7LiU/Zr8yiuv4MKFC/joo4+QmpqKgIAADBkypED9Mjc3x6BBgxAeHo5Hjx5hxYoVz72hhZ2dndHYlFt/ny17OqYUhqWlJTw9PdGuXTssXboU5ubm+P7775X1xuIdAPTs2RNRUVEG3yMuKMYnKiomWVRmvfnmm7CwsMDKlSsRHh6OQYMGKS/gLVu2RFxcHNzd3eHp6WnwY21tnaMtGxsbuLi44L///a9BeVRUFBo1agQAaNGiBY4fP57v/69RpUoVZGVl5SivUaMGXn/9dYSHhyM8PBwhISEwMcn9qVa3bl3Y2Ngon6gVxdmzZ3HmzBn06NGjyG0BQLNmzTB58mTs3LkTXbp0wfLlywH8b/6fnXtPT0/odDo0bNgQWq0W27Zty7P9EydOoGXLlqr0lYiotJ04cQJbtmzBm2++CeBxXBERJCcn53itdHV1LVDbdnZ2CA4ORkREBJYuXYrly5cjMzPTaF1zc3Oj8endd9/Fjh07sHDhQqSlpaFfv3557rNFixaqxCa16PV6Zcz37t3Dtm3b8Nprr+Wo99lnn+G9996Dn5+fwY0yGjRogBMnThh8QnXgwAGj+2J8oqJikkVllk6nQ//+/REWFoazZ88avOM2atQo3L59G0FBQdi/fz/OnTuHP//8E8OHD0daWprR9iZOnIg5c+Zg9erVOHXqFKZNm4Y9e/bgo48+AgD0798fDg4O6N27N/bs2YNz585h48aNyt0Fn+Xh4YELFy7gyJEjSE1NNQh2w4YNw6pVq3Ds2LHn/oNhExMTdO3aFVFRUQWdohw2bNgAPz+/Iv/T5PPnz2PKlCmIjo7GhQsXsGPHDhw/flxJSCdPnowDBw7gvffeQ0xMDM6cOYP/+7//w4gRIwAA1tbWGDduHKZOnYply5bh7NmzOHDgABYtWqTsIz09HYcPH1YtISQiKkmZmZlITk5GUlISjh07hnnz5qFLly7w9vZW4kr9+vUxYMAADB48GGvXrsW5c+dw6NAhzJkzp0B3Vp02bRrWr1+P06dPIyEhAb/88gvq1KkDrVZrtL6Hhwe2bduG5ORkg0+ivLy80KlTJ0ycOBFvvvkmbGxs8txv9+7dkZCQYHDXv5Jw584dfPLJJ9i/fz8uXryIw4cP45133sGlS5fw1ltvAQC2bt2KOnXqwNPT02gb//znPzFixAh07dpV+YfH/fv3h6mpKYYNG4b4+Hj8+eef+OKLLwAYXtWSmJiIy5cvq3pbeKp8mGRRmfbuu+/i5s2b6NChAxo2bKiUZ38qZWJigh49eqBx48YYNWoUtFptrkFn7NixmDhxIiZNmoQmTZrg119/xbp169CsWTMAjy9L2LVrF9zc3BAYGIgmTZogNDQ018sJ33jjDbz66qvw8/ODvb09Vq9erazr0qUL3Nzc0KVLl3z976vhw4fjp59+MvrOY0HkdulEQVWtWhWnTp1C3759Ub9+fYSEhGDAgAGYPHkygMfX4u/evRuJiYnw9fVFs2bNMHXqVOW6fACYOXMmJk+ejM8++wwNGzbEG2+8gUuXLinrN2zYgFq1asHHx6fI/SUiKml//vknnJ2dUatWLfj5+WHjxo0ICwvD7t27YWlpqdRbtmwZhgwZgkmTJqFBgwbo2bMndu/eDXd393zvS6vV4h//+AeaNWuGjh07Ii0tDZs2bcq1/ty5c7Fjxw7UrFkTLVq0MFg3dOhQPHjwIF//+6pp06Zo06YN1qxZk+++qsHMzAxxcXHo06cP6tWrh8DAQFy/fh27d+9WvguWn3j3xRdfYNiwYfDz88OxY8dgbW2NTZs2IS4uDi1atMDEiRMRFhYGwPB70atXr4a/v3+BjhHRszQiIqXdCaKK5v79+3B1dcX8+fMxYMCAfG3Tvn17vP/++wgODi7UPlNTU+Hs7Iy///4bTk5OhWqjJLVp0wYffPAB+vfvX9pdISKqNL788kssXboUp06dylf9LVu2YNy4cYiPj4epqWkx9y5/srKy4ODggN9//x1t2rQpUlsbNmxAnz59kJKSAjs7O2RmZqJevXpYvXp1jptlERWEaVh2Ck9ERabX63Ht2jXMnj0bcXFxWLJkSb6DUps2bXDlypUC3TjiaZcuXYKnp+dz7+hXFqSkpODBgwcYMWKE6jceISKinO7evYvTp09j1KhRmDhxYr6TE09PT4gIXFxcUK1atWLuZf6kpqaiSpUq6NevX4FjyPLly5GWlgYRQVRUFMaNG4fOnTsrn+ydO3cONWvWRJ8+fYqj61SJ8JMsIhUlJiaidu3acHNzw9KlS3k9NxERlQmDBw/G6tWr0atXL6xevTrHXWcri6+++goLFy7ElStX4OTkhICAAHz55ZdGb5pFVBRMsoiIiIiIiFTEG18QERERERGpiEkWERERERGRiphkERERERERqYhJFhERERERkYqYZBEREREREamISRYREREREZGKmGQRERERERGpiEkWERERERGRiphkERERERERqYhJFhERERERkYqYZBEREREREamISRYREREREZGKmGQRERERERGpiEkWERERERGRiphkERERERERqYhJFhERERERkYqYZBEREREREamISRYREREREZGKmGQRERERERGpiEkWERERERGRiphkERERERERqYhJFhERERERkYqYZBEREREREamISRYREREREZGKmGQRERERERGpiEkWERERERGRiphkERERERERqYhJFhERERERkYqYZBEREREREamISRYREREREZGKmGRRmeHh4YHBgweX2/YrorCwMGg0GiQnJ5fI/niMiKi0RUREQKPRIDExsbS7UuJiYmLg4+MDKysraDQaHD16tLS7VOyy41xhVOZzhZ6PSRYVWK9evWBhYYFbt27lWmfs2LHQaDQ4depUCfasYJKSkhAWFqZ6ENm5cyc0Go3yY2pqCkdHR/Tt2xcJCQmq7mvz5s0ICwtTtc3iVh77TESkprt37yI0NBSvvvoq7O3todFoMGvWrFzrZ2ZmYsqUKXB1dYVOp0ObNm2wdetWVfuk1+sRFBSEK1euYO7cuVixYgXc3d2fu92ZM2cwZswYeHp6wsLCAi+88AICAwOxb98+VftHVN4wyaICGzhwIDIzM7Fu3Tqj67OysvDTTz+hdevWqF+/fgn3Lnd//fUXlixZoiwnJSVh+vTpxfZO3ahRo7BixQp8//33CAoKwqZNm9CpUyckJSWpto/Nmzdj+vTpqrVXEvLq87PHiIioIkpNTcWMGTMQGxuLFi1aPLf+4MGDMXfuXLz99tv45ptvYG5ujoCAAOzatUu1PiUlJeH06dMYM2YMRowYgYEDB6J69ep5bhMREYEXX3wR69atw5tvvolvv/0WI0eOxIEDB+Dj44Nff/1Vtf4RlTdmpd0BKn969eoFGxsb/Pjjjxg6dGiO9ZGRkUhJScE//vGPUuhd7rRabYnur1OnTujXrx8AYMiQIahXrx7Gjh2L5cuXY+rUqUVq+969e7C0tFSjm6XSfm5K+hgREZUGZ2dnXL58GS4uLkhMTETt2rVzrXvgwAGsWbMGs2bNwuTJkwEAgwYNQpMmTTBx4kQcOHBAlT6lpKQAAGxtbfNV/6effsI777yD3r1748cff4ROp1PWjR49Gk2bNsWwYcPg7+9fKvGEqLTxkywqMAsLC7zxxhvYuXOn0U9lVq1aBVNTUyXByMzMxPTp01GvXj1otVq4urriww8/RHp6+nP3lZqaiuHDh8PJyQkWFhZo0qSJ0U86RAQLFixA8+bNodPpYGdnh27dumHPnj1Knae/77Nz5060bt0awOMEKPvSvrCwMCxZsgQajQZHjhzJsZ9vv/0WGo2mUJf9vfzyywCA8+fPK2V//PEHfH19YWVlBSsrK/To0SPHJ2uDBw+GhYUFEhMTlQQ3ICAAgwcPxoIFCwDA4PLExMREJCYmQqPRICIiIkc/sseZLft69BMnTiA4OBg1atRAkyZNDLa5ceMG+vfvD1tbW1SvXh0jRozA3bt3Deps3LgRgYGBcHNzg1arhbu7OyZOnIiMjAyDseTWZ8D4d7Lycw5kj3fWrFlYsmQJ6tatC61Wi9atW+PgwYO5HBEiovw5duwYXn31VdjY2MDS0hJdunQxiC/Zdu7ciVatWsHCwgJ169bFv//9b6Pf+dFqtXBxccnXvteuXQsTExMMHz5cKbOwsMDQoUNx8ODBfH0f6Hn9Hzx4MLy9vQH8LyZ26dIl1/auXLmCESNGoHHjxlizZo1BggUATk5OGDNmDK5fv/7cyxqz5ychIQEDBw6Era0t7Ozs8PHHH0NEcD4aBjgAACAASURBVPnyZbz22muwsbGBo6MjZs+eXeDxZYuKikLr1q0Njk9eY3z33Xfh5OQErVaLhg0bYtGiRXmOhehp/CSLCmXgwIFYtmwZ1qxZg/Hjxyvl6enpWL9+Pbp16wYHBweICPr06YNdu3Zh2LBhaNSoERISErBw4ULExcVh69atuX7hNCMjAy+99BJOnjyJUaNGoW7duli/fj2GDx+O69evY8qUKUrd4cOH4/vvv0f37t0xePBgiAiio6Oxe/du+Pj45Gi7YcOGmDFjBqZNm4bhw4crdV588UW4u7tj7NixWLlyJVq2bGmw3cqVK9GqVSs0bNiwwHN29uxZAMALL7wAAPjxxx8xcOBAdOvWDTNnzkRmZiYWL14MHx8fHDx4EA0aNFC21ev18Pf3R5s2bTB79myYmZmhUaNGSEpKQmRkJFasWKHUtbe3x7Vr1wrcv6CgINSuXRuff/45Hjx4YLCuX79+cHV1xRdffIGjR49i8eLF+Pvvv7F582alzrJly6DVajF27FjY2tpi3759+Ne//oW///4ba9asAQCMGDEi1z4bU5BzAHj8zurdu3cxYsQIaDQafPXVV3j99ddx7tw5mJubF3hOiIgSEhLg4+MDS0tLTJw4ERYWFliyZAm6du2KyMhIdO7cGcDjm0b06NEDTk5OCAsLg16vx4wZM2BnZ1ek/cfExKBu3bo5Lt1r06aNst7Dw6NI/R8xYgTq1q1rEBMdHR1zbXPOnDm4ffs21qxZk+sVCE2bNgWAfH83++2330aDBg0wa9YsbN68GTNnzkSNGjUQHh6Ozp0748svv8SqVaswadIkeHt7K29c5vf4xMbGwt/fH/b29ggLC0NWVhamT59uNP6kpKSgXbt2yMrKwvvvvw8HBwds27YN77//Pq5fv45PPvkkX2OiSk6ICiErK0vc3NzE29vboPzHH38UALJy5UoREVm1apVoNBrZuXOnQb2VK1cKANm6datS5u7uLiEhIcryN998IwAkIiJCKXv06JH4+f1/9u48rqo6/+P4+4qICEomuwtqCZrLqFAu5NKkYiiWWrkkoqPilOMvU7OamqRGrWmm0tGxojG1XBrHbDFNxS0XsHEZc0nNSmVc2AoXlFDh/P4gbl7Zrnou94qv5+NxH957vt97zucA8uFzv9/zPfcbHh4eRlZWlmEYhrFhwwZDkvHEE08Ui7OgoKDU/W/fvt2QZMydO7fY+wYMGGAEBgYaly9ftm47dOiQIcmYMWNGGV+ZX+NJTEw0MjMzjZMnTxorVqwwGjVqZFgsFmP79u1GTk6OUbt2bWP48OE27/3pp58MPz8/Y9CgQdZtcXFxhiTjqaeeKnasMWPGGCX9Nz5y5Eip5ybJmDx5svX15MmTDUlGv379ivUtauvZs6fN1/JPf/qTIclISkqybjt//nyx90+dOtWwWCxGampquTEbxvX/DBSdb506dYyffvrJ2vfTTz81JBnLly8v8XgAcLW5c+cakowjR44YhmEYffv2Ndzd3Y1vv/3W2iczM9OoU6eOTQ6MiYkxqlevbvzvf/+zbjt8+LBRtWrVUn/nGcavv79eeeWVEtubN29udO7cudj2/fv3G5KMWbNmlXk+9sZfVk68UkFBgeHr62s0bNjQJi9c7YsvvjAkGVOmTClzf0V55ne/+5112+XLl4169eoZFovF5v3Z2dmGp6en8dhjj13z+T300EOGh4eHcezYMeu2Q4cOlfj9GTVqlBEQEGBkZGTYbB85cqTh6elpZGdnG4ZR/GcFuBLTBXFdqlSpokGDBmnnzp02n1ItXLhQXl5eeuihhyRJS5YsUWhoqJo3b66srCzro0uXLrJYLNqwYUOpx1ixYoX8/Pw0ZMgQ6zY3NzeNGzdOeXl5Wrt2raTCqRSSSlxM4XqXZY2Li1NaWpr1GFLhKFbVqlWt0yDLEx8fLz8/PwUHB6tXr17Kz8/X4sWLFRERoaSkJGVnZ2vw4ME2X5f8/Hx16tSpxK/LE088cV3nYq/HH3+81LY//OEPNl/L//u//5Mkff7559ZtNWrUkFQ46nbmzBllZWXp3nvvlWEYJU69tIe9PwNF+vfvb/Npb9EI5Q8//HBdxwdwa8vPz9fq1asVExOjJk2aWLf7+vpq2LBh2rlzp9LT05Wfn6+1a9eqT58+qlevnrXfnXfeqQceeOCGYsjNzS1xtKh69erW9huN/1ocOHBAWVlZeuCBB8rMsUW/d+2dFjly5Ejrczc3N0VERMgwDJtrv2+77TaFhYVZ930t35/Vq1erT58+atCggbVfaGiooqKibOIwDENLly5Vr169ZLFYbHJ0jx49lJubq6+++squc8KtjemCuG5DhgzRX//6Vy1cuFAvvfSSsrKytHr1ag0YMMB6keu3336rQ4cOlTodrOhC25IcO3ZMd955p9zc3Gy2F03VK5qH/v333ysgIOCGp2RcqUePHgoMDNSCBQusv4AXLlyoqKgo+fv727WP559/Xl27dlW1atVUr149hYSEWM+lqDDt3r17ie+tUqVKsddlTQcxwx133FFq25XJSypMYLVr17a5FmDfvn2aNGmSNm7cWCzpnzlz5rpisvdnoMiVyVOSteDKzs6+ruMDuLVlZmbqwoULCgsLK9Z25e+hgoIC5ebm6s477yzWr6Rt18LT01N5eXnFthdd73r19VBXsjf+sqYGXu348eOSVO7y7uvXr5ckRUZG2rXfq39/+/j4yN3dXYGBgcW2FxWG9p6fYRjKzc0tlsukwkJrxYoV1teZmZnKzs7We++9p/fee6/EWMv62wUoQpGF69aqVSu1bNlSixcv1ksvvaQlS5bo8uXLNqMOBQUFuuuuuzRjxowS92HvJ1wVzc3NTY899pjefvttnT9/Xl9//bV++OEHTZs2ze59tGjRQt26dSuxraCgQFLh8rd169Ytd1/u7u6qWtX+/66lfbqYn59f6nvKStTlOXPmjO677z55eXlp6tSpuvPOO+Xp6akTJ05o2LBh1vN1tKuLsSKGYVTI8QHAbEFBQTp27Fix7adOnZLkvDxa1uJVJ0+e1Oeff6527drZfSuXkn5/X/2BYxFH/k4vyleDBg3S7373uxL7NG/e3GHHR+VBkYUbMmTIED3zzDPavn27Fi5cqICAAJvRmTvuuEM7d+7U/ffff81T90JCQvTf//5X+fn5Nr98Dx48KEnWkZ077rhDq1atUmZmZqkjZiUpL564uDi9/vrr+vjjj5WcnKxatWrpwQcfvKZzKE3RqJGfn1+phZg9SjuHohGcq28YXVKitsfhw4dtEmVWVpays7Ot34MNGzYoKytLS5cuVZcuXaz9kpKS7I65JPb+DACAI/j5+alGjRo6dOhQsbYrfw/5+vqqevXq+u6774r1K2nbtWjdurXWr1+v7Oxsm+nQRVPWWrdufcPxX4uiXLB3795S+4wdO1Z5eXmaOnXqNe37Wl3L98fT01OHDx8u1u/qhTn8/PxUs2ZNXb58+YbyM8A1WbghgwcPVpUqVTRlyhQlJydrwIABNn8MDxgwQOnp6SUue5qXl6dz586Vuu/evXsrMzNTixYtsm4rKCjQjBkz5OHhYf3l9/DDD0uSzbLkRcr6tKtoSmNpU8latmypNm3aaO7cuVqyZIkeeeQR6xz4GxUVFaXbbrtN06ZNK7aSnyS7Vwcs7Rxq1aolX19fbdq0yWb77NmzryveWbNm2Xwt//73v0uSevXqJenXTyCv7FNQUKA33njD7phLYu/PAAA4gpubm3r27Knly5dbV4iVCm9rMX/+fEVERCggIEBubm7q1q2bPvvsM+t0OqmwwPriiy9uKIaHH35YBQUFSkxMtG7Ly8vT3LlzFR4eXuY9tuyN/1o0bNhQ99xzjz777DPt2LHDps0wDD399NNatmyZJk6cqPvvv/+a9n2truX7ExUVpeXLlys1NdXa79tvvy22xLybm5sefvhhffLJJ/r666+LHfN6Vu/FrYmRLNyQevXqqUuXLvrss88kyWaqYNHrpUuXasyYMfryyy+tCyEcOnRIS5Ys0b///e9S78UxatQoJSYmasSIEfrvf/+rxo0b65NPPtG6dev0yiuvWJdC79q1q4YNG6bZs2fr+++/t15knJKSolatWumPf/xjifsvWhL3rbfekre3t2rWrKkWLVrY3CNq6NCheuqppyRJsbGxN/S1ulKtWrX09ttv67HHHlObNm00aNAgBQQEKDU1VatWrVLz5s1LvMfV1SIiIiQVLkzxwAMPqGrVqoqJiZGXl5dGjhypV199VSNHjlRERIQ2bdpk91K6Vztx4oSio6PVu3dvff3113r33XfVo0cP66hlZGSk6tSpo7i4OI0dO1bu7u5aunRpsXtplRfz1ez9GQAAR5kyZYrWrFmje++9V2PGjLEuEX769GnrwktS4Qd9Rf0ef/xxFRQUaNasWWrevHmJf6zPmjVLp0+fts442LBhgy5fviypcCSo6KbA7dq10yOPPKIXXnhBWVlZatKkid5//30dOXKkxNkC1xv/tUhMTFSXLl2sy783a9ZMaWlp+vDDD3Xw4EE9/fTT+stf/nJd+75W9p7fSy+9pFWrVqlTp04235+77rpLe/bssdnnq6++qo0bN6pDhw4aNWqUmjdvruzsbO3evVsff/yxzf0fgVI5aVVDVCJz5swxJBmhoaEltl+6dMn461//arRo0cLw8PAwbrvtNqNt27bGiy++aPz444/Wflcv320Yhcuwjhw50vD39zeqVatm3HXXXUZiYmKxY+Tn5xtvvPGGcddddxnVqlUz6tSpY3Tv3t3YsmVLmfv//PPPjZYtWxru7u7FljY3DMNIT083qlataoSEhJS5VO2VipZwX7x4cbl9N23aZPTs2dO47bbbjOrVqxuNGzc2YmNjjeTkZGufuLg4w8PDo8T35+fnG+PGjTMCAgIMi8Vis5TshQsXjBEjRhg+Pj5GzZo1jUcffdTIyMgodQn3U6dOFdt/Udu+ffuMQYMGGbVq1TJ8fHyMkSNHGmfOnLHpu23bNiMyMtKoUaOG4e/vbzz++OPGnj17ii0JXFbM1/szUNYSyCV9XwGgNCUty717927jgQceMLy9vY0aNWoYnTt3Nr788sti7123bp3Rpk0bo1q1asYdd9xh/POf/zQmTJhgVK9evVjfkJAQQ1KJj6uXBM/NzTWefvppIygoyPDw8DDCw8ONlStX2n1O9sRv7xLuRQ4fPmzExcUZQUFBRpUqVQxJRr169YytW7faHZdhlJ6DSst9Xbp0McLCwq75/AzDML788ksjPDzcqFatmtG4cWPj7bffth7/ahkZGcbYsWONBg0aGO7u7kZAQIDRtWtXmyXzWcIdZbEYBleEA6XJzs5WYGCgJk6c6PC55QCAyuehhx7S/v37S7weqDKJj4/XnDlztGHDBusNgIFbWaW+Jmv27Nlq1KiRqlevrvDwcG3evNnZIeEmM3/+fF28eFFxcXHODgVAJUJ+qpyuvn3F4cOHtXLlylKnxVcmb775pu644w7FxsZe9207gMqk0o5k/etf/9KQIUM0e/Zs3XvvvZo9e7bmzp2rb775pti9GICrrV+/XgcOHNCf/vQnRUZGavny5c4OCUAlQX6qvIKCgjRs2DA1btxYx44d01tvvaW8vDz997//LfEeTQAqr0pbZLVr106tWrXSu+++a93WpEkTPfzww3rllVecGBluBl27dlVycrI6dOigBQsWqH79+s4OCUAlQX6qvIYPH64NGzYoLS1NHh4e6tChg6ZNm6a2bds6OzQAFaxSri548eJF7dy5UxMnTrTZ3qNHDyUnJzspKtxMNm7c6OwQAFRC5KfKbe7cuc4OAYCLqJRFVlZWlvLz84vd+yEgIEBr164t1j8xMdF6/4mDBw+qadOmFRInALiSo0ePKisry9lhVGo3kp927twjybciwgQAl1OnzsWbKkdVyiLrWsXHxys+Pl5S4T18dvznP06OCAAqXsQ99zg7BFzlyvxksQRLinduQADgJA0bfu7sEK5JpVxd0NfXV25ubkpPT7fZnp6ersDAQCdFBQC41ZGfAODWUClHsqpVq6bw8HAlJSXpkUcesW5PSkpS//79nRgZAOBWZl5+Gm9+cADgct5wdgDXrVIWWZI0fvx4xcbG6p577lFkZKTefvttnTx5Ur///e+dHRoA4BZGfgKAyq/SFlkDBgzQjz/+qClTpujUqVNq0aKFVq5cqZCQEGeHBgC4hZGfAKDyq7RFliQ98cQTeuKJJ5wdBgAANkzNT71rmbMfAHAFn591dgSmqJQLXwAAAACAs1BkAQAAAICJKLIAAAAAwESV+posAABuOSOdHQAAXIN/OjsAx2AkCwAAAABMRJEFAAAAACZiuiAAAJVUswd3OTsEACjmwKdtnR2CwzGSBQAAAAAmosgCAAAAABMxXRAAgFtAfy11dggAbmEf6WFnh1ChGMkCAAAAABNRZAEAAACAiZguCADALaavPnF2CABuAR/rIWeH4DSMZAEAAACAiSiyAAAAAMBETBcEAOAW1nbrAWeHAKAS2RXZzNkhuARGsgAAAADARBRZAAAAAGAipgsCAIBCHzg7AAA3pVhnB+B6GMkCAAAAABNRZAEAAACAiSiyAAAAAMBEXJMFAABKdHaesyMA4IpqDXN2BK6PkSwAAAAAMBFFFgBJ0rinntK4p55ydhgAABTz7OXCB3CzYLogAEnS7q+/dnYIAFzYG3nOjgC3spW//Ouf79QwbmnjPZwdwc2FkSwAAAAAMBFFFgAAAACYiCILAAAAAEzkUkXWK6+8orvvvlu1atWSn5+fYmJitG/fPps+w4YNk8VisXm0b9/epk9eXp7Gjh0rX19feXl5qU+fPjp+/HhFngoAoJIhRwEA7OVSRdbGjRv1xBNPKDk5WevXr1fVqlXVrVs3/fTTTzb9unXrplOnTlkfK1eutGkfN26cPvroIy1evFibN2/W2bNn1bt3b+Xnc7UkAOD6kKMAAPZyqdUFV69ebfP6gw8+kI+Pj7Zu3aqYmBjrdg8PDwUGBpa4jzNnzmjOnDmaO3euunfvbt1PSEiI1q5dq6ioKMedAACg0iJHAQDs5VIjWVc7d+6cCgoKVLt2bZvtW7Zskb+/v0JDQzVq1ChlZGRY23bu3KlLly6pR48e1m3169dXs2bNlJycXGGxAwAqN3IUAKA0LjWSdbUnn3xSrVu3VocOHazbevbsqX79+qlRo0Y6evSoXnjhBf32t7/Vzp075eHhobS0NLm5ucnX19dmXwEBAUpLSyvxOImJiUpMTJQkZWZmOu6EAACVRkXkqCvzk3TBkacDADCRyxZZ48eP15YtW7Rlyxa5ublZtw8cOND6vGXLlgoPD1dISIhWrFihfv36Xdex4uPjFR8fL0mKiIi4scABAJVeReWoK/OTxRJ844EDACqES04XfOqpp7R48WKtX79ejRs3LrNvcHCw6tWrp8OHD0uSAgMDlZ+fr6ysLJt+6enppc6RBwDAXuQoAEB5XK7IevLJJ63Jq2nTpuX2z8rK0okTJxQUFCRJCg8Pl7u7u5KSkqx9jh8/rgMHDqhjx44OixsAUPmRowAA9nCp6YJjxozRBx98oE8++US1a9e2zk/39vaWt7e3cnJylJCQoP79+ysoKEhHjx7Vc889J39/f/Xt21eS5OPjoxEjRmjSpEny9/dXnTp1NH78eLVq1UrdunVz5ukBAG5i5CgAgL1cqsiaPXu2JOn++++32T558mQlJCTIzc1Ne/fu1fvvv6/Tp08rKChI9913n5YsWaKaNWta+0+fPl1Vq1bVgAEDlJubq/vvv1/vv/++zbx5AACuBTkKAGAvlyqyDMMos93T07PYfUpK4uHhoZkzZ2rmzJlmhQYAuMWRowAA9nK5a7IAAAAA4GZGkQUAAAAAJqLIAgAAAAATUWQBAAAAgIkosgAAAADARBRZAAAAAGAiiiwAAAAAMBFFFgAAAACYiCILAAAAAExEkQUAAAAAJqLIAgAAAAATUWQBAAAAgIkosgAAAADARBRZAAAAAGAiiiwAAAAAMBFFFgAAAACYiCILAAAAAExEkQUAAAAAJqLIAgAAAAATUWQBAAAAgIkosgAAAADARBRZAAAAAGAiiiwAAAAAMBFFFgAAAACYiCILAAAAAExEkQUAAAAAJqLIAgAAAAATUWQBAAAAgIkosgAAAADARC5VZCUkJMhisdg8AgMDre2GYSghIUHBwcHy9PRU165dtX//fpt9ZGdnKzY2Vj4+PvLx8VFsbKxOnz5d0acCAKhkyFEAAHu5VJElSWFhYTp16pT1sXfvXmvba6+9ptdff10zZ87U9u3b5e/vr+7du+vcuXPWPoMHD9auXbu0atUqrVq1Srt27VJsbKwzTgUAUMmQowAA9qjq7ACuVrVqVZtPBosYhqHp06fr2WefVf/+/SVJ8+fPl7+/vxYtWqTRo0frwIEDWrVqlbZs2aIOHTpIkt555x116tRJhw4dUlhYWIWeCwCgciFHAQDs4XIjWT/88IOCg4PVqFEjDRw4UD/88IMk6ciRI0pLS1OPHj2sfT09PdW5c2clJydLklJSUuTt7a2OHTta+0RGRsrLy8vaBwCA60WOAgDYw6WKrHbt2mnevHlatWqV3n33XaWlpaljx4768ccflZaWJkkKCAiweU9AQIC1LS0tTX5+frJYLNZ2i8Uif39/a5+SJCYmKiIiQhEREcrMzHTAmQEAbnbOyFFX5ifpgmNODABgOpeaLvjAAw/YvG7fvr0aN26s+fPnq3379g47bnx8vOLj4yXpl0QGAIAtZ+SoK/OTxRLskGMAAMznUiNZV/P29lbz5s11+PBh6xz49PR0mz7p6enWtsDAQGVmZsowDGu7YRjKyMgocQ49AADXixwFACiNSxdZP//8sw4ePKigoCA1atRIgYGBSkpKsmnfvHmzdX57hw4dlJOTo5SUFGuflJQUnT9/3mYOPAAAN4ocBQAojUtNF5w4caJiYmLUoEEDZWRk6M9//rPOnz+vuLg4WSwWjRs3TtOmTVPTpk0VGhqqKVOmyNvbW4MHD5YkNWvWTD179tTo0aOVmJgoSRo9erR69+7Nqk0AgBtCjgIA2Muliqzjx49r0KBBysrKkp+fn9q3b69t27YpJCREkjRp0iTl5uZqzJgxys7OVrt27bRmzRrVrFnTuo9FixZp7NixioqKkiT16dNHs2bNcsr5AAAqD3IUAMBeLlVkffjhh2W2WywWJSQkKCEhodQ+tWvX1oIFC0yODABwqyNHAQDs5dLXZAEAAADAzYYiCwAAAABMRJEFAAAAACaiyAIAAAAAE1FkAQAAAICJKLIAAAAAwEQUWQAAAABgIoosAAAAADARRRYAAAAAmIgiCwAAAABMRJEFAAAAACaiyAIAAAAAE1FkAQAAAICJKLIAAAAAwEQUWQAAAABgIoosAAAAADBRVXs6FRQUaOPGjfryyy919OhR5ebmys/PT23btlWPHj1Uv359R8cJAEAx5CcAgCsqcyQrNzdXU6dOVf369dWrVy+tWbNGOTk5qlatmo4cOaKXXnpJjRo1UnR0tLZt21ZRMQMAbnHkJwCAKytzJKtJkybq0KGDEhMT1aNHD7m7uxfrc+zYMS1atEgDBgzQCy+8oFGjRjksWAAAJPITAMC1lVlkrVq1Si1atChzByEhIXruuec0YcIEHTt2zNTgAAAoCfkJAODKypwuWF4Cu1K1atXUpEmTGw4IAIDykJ8AAK7MroUvily8eFH79u1TRkaGCgoKbNqio6NNDQwAAHuRnwAArsTuIispKUmxsbHKyMgo1maxWJSfn29qYAAA2IP8BABwNXbfJ2vMmDHq3bu3jhw5ogsXLig3N9f6uHDhgiNjBACgVOQnAICrsXsk69SpU/rjH/+okJAQR8YDAMA1IT8BAFyN3SNZvXv3VnJysiNjAQDgmpGfAACuxu6RrLfffluPPfaYdu7cqRYtWhS7J8nQoUNNDw4AgPKQnwAArsbuImv16tVat26dVq5cqRo1ashisVjbLBYLSQwA4BTkJwCAq7F7uuDEiRP1hz/8QefOnVNOTo7OnTtnfZw9e9aRMQIAUCryEwDA1dhdZJ0+fVq///3v5eXl5bBgGjZsKIvFUuzRq1cvSVJCQkKxtsDAQJt9GIahhIQEBQcHy9PTU127dtX+/fsdFjMAwLkqIj9J5CgAgP3sLrL69++vtWvXOjIWbd++XadOnbI+du3aJYvFokcffdTaJywszKbP3r17bfbx2muv6fXXX9fMmTO1fft2+fv7q3v37jp37pxDYwcAOEdF5CeJHAUAsJ/d12Q1btxYzz//vDZt2qRWrVoVu7B4/PjxNxyMn5+fzes5c+aoVq1aNgmsatWqxT4ZLGIYhqZPn65nn31W/fv3lyTNnz9f/v7+WrRokUaPHn3DMQIAXEtF5CeJHAUAsJ/dRdZ7772nmjVrKjk5udhSuRaLxbQkVsQwDM2ZM0dDhgyRp6endfsPP/yg4OBgeXh4qF27dpo2bZoaN24sSTpy5IjS0tLUo0cPa39PT0917txZycnJJDAAqIQqOj9J5CgAQNnsLrKOHDniyDiKSUpK0pEjRzRq1Cjrtnbt2mnevHlq2rSpMjIyNGXKFHXs2FH79+9XnTp1lJaWJkkKCAiw2VdAQIBOnDhRofEDACpGRecniRwFACib3UVWRXv33Xd199136ze/+Y112wMPPGDTp3379mrcuLHmz59/Q59UJiYmKjExUZKUmZl53fsBANwaKipHXZmfpAvXGy4AoIKVufDFlClTdP78ebt2tHXrVi1fvtyUoDIyMvTpp5/afEJYEm9vbzVv3lyHDx+WJOs8+PT0dJt+6enppc6Rl6T4+Hjt2LFDO3bsKDbnHgDgepyVn6SKzVFX5iepxo0HDwCoEGUWWd9//70a+rnm/QAAIABJREFUNGig+Ph4LV++XKdOnbK2/fzzz9q1a5f+/ve/65577lFsbKxq165tSlDz5s2Th4eHBg0aVGa/n3/+WQcPHlRQUJAkqVGjRgoMDFRSUpJNn82bN6tjx46mxAYAcD5n5SeJHAUAKF+Z0wXnzp2rvXv3atasWRo6dKjOnj0ri8Uid3d3Xbx4UYZhqG3btoqPj1dcXJw8PDxuOCDDMPTPf/5TAwcOlLe3t03bxIkTFRMTowYNGigjI0N//vOfdf78ecXFxUkqvMB53LhxmjZtmpo2barQ0FBNmTJF3t7eGjx48A3HBgBwDc7ITxI5CgBgn3KvyWrZsqXeeecdvfXWW9qzZ4+OHTum3Nxc+fr6qnXr1vL19TU1oI0bN+rw4cNasGBBsbbjx49r0KBBysrKkp+fn9q3b69t27YpJCTE2mfSpEnKzc3VmDFjlJ2drXbt2mnNmjWqWbOmqXECAJyrovOTRI4CANjH7oUvqlSpotatW6t169aOjEf33XefDMMose3DDz8s9/0Wi0UJCQlKSEgwOTIAgCuqqPwkkaMAAPYp85osAAAAAMC1ocgCAAAAABNRZAEAAACAiSiyAAAAAMBE11Vkpaenq6CgwOxYAAC4IeQnAIArsLvIunTpkiZNmqSaNWuqbt26Onr0qCTpmWee0ezZsx0VHwAAZSI/AQBcjd1F1ksvvaTly5drwYIFNjd1vOeeezRv3jxHxAYAQLnITwAAV2P3fbIWL16s9957T126dFGVKr/WZi1atNC3337rkOAAACgP+QkA4GrsHsk6efKkzV3ri1y+fFmXL182NSgAAOxFfgIAuBq7i6zmzZtr06ZNxbYvWbJE4eHhpgYFAIC9yE8AAFdj93TByZMna8iQIfrf//6n/Px8/fvf/9bBgwe1aNEirVixwpExAgBQKvITAMDV2D2SFRMToyVLlmjNmjWqUqWKXnrpJR0+fFjLly9Xt27dHBkjAAClIj8BAFyN3SNZkhQVFaWoqChHxQIAwHUhPwEAXIndI1mNGzfWjz/+WGz76dOn1bhxY1ODAgDAXuQnAICrsbvIOnr0qPLz84ttz8vL04kTJ0wNCgAAe5GfAACuptzpgsuWLbM+X7FihXx8fKyv8/PztW7dOjVs2NAhwQEAUBryEwDAVZVbZD388MOSJIvFohEjRti0ubu7q2HDhnr99dcdEx0AAKUgPwEAXFW5RVZBQYEkqVGjRtq+fbt8fX0dHhQAAOUhPwEAXJXdqwseOXLEkXEAAHBdyE8AAFdzTUu4Z2dn64svvlBqaqouXrxo0/biiy+aGhgAAPYiPwEAXIndRda2bdvUq1cveXh4KDMzU3Xr1tWpU6fk4eGhhg0bksQAAE5BfgIAuBq7l3B/+umn9dhjj+nEiROqXr261q9fr9TUVEVEROiZZ55xZIwAAJSK/AQAcDV2j2Tt2bNHc+bMkcVikZubm/Ly8tS4cWP95S9/0eDBg/XYY485Mk7gppXg5lbxxyzhnkFAZUV+AgC4GrtHsqpVq2Z9HhAQoGPHjkmSvL29dfLkSfMjAwDADuQnAICrsXskq23bttq+fbtCQ0PVtWtXvfDCC0pPT9eCBQvUqlUrR8YIAECpyE8AAFdj90jW1KlTFRwcLEmaMmWK/Pz8NHbsWGVnZysxMdFhAQIAUBbyEwDA1dg9khUREWF97ufnpy+++MIhAQEAcC3ITwAAV2P3SBZQ6f34oxQdXfgvAAAuwlPSY7/8C+DmYHeR9dNPP+nxxx9XaGiobrvtNtWqVcvmAdz05s+XVq+W3n/f2ZEAuAbkJ1R2rSXdKek3zg4EgN3sLrJGjBihL774QsOGDdP06dM1c+ZMm4c9Nm3apD59+qhu3bqyWCyaN2+eTbthGEpISFBwcLA8PT3VtWtX7d+/36ZPdna2YmNj5ePjIx8fH8XGxur06dM2ffbu3asuXbrI09NTdevW1csvvyzDMOw9VdyKDEOaPr3w+fTpha8B3BTMyE8SOQquq70ki6QOzg4EgN3sviZr3bp1SkpKUrt27a77YDk5OWrRooWGDh2qoUOHFmt/7bXX9Prrr2vevHkKCwvTyy+/rO7du+vQoUOqWbOmJGnw4MFKTU3VqlWrJEkjR45UbGysli9fLkk6e/asunfvrs6dO2v79u06ePCghg8fLi8vL02YMOG6Y0clt3mzdOZM4fPsbGnLFqlTJ+fGBMAuZuQniRwF1xQiqfovz6tLaiAp1XnhALCT3UWWv7+/vL29b+hg0dHRio6OliQNGzbMps0wDE2fPl3PPvus+vfvL0maP3++/P39tWjRIo0ePVoHDhzQqlWrtGXLFnXoUPh5zjvvvKNOnTrp0KFDCgsL08KFC3XhwgXNnz9fnp6eatGihQ4ePKg33nhD48ePl8ViuaFzQCU1Y4Z0/nzh8wsXCkezKLKAm4IZ+UkiR8E1tZfk/stzdxWOZlFkAa7vmpZwf/HFF5WTk+OQQI4cOaK0tDT16NHDus3T01OdO3dWcnKyJCklJUXe3t7q2LGjtU9kZKS8vLxs+nTq1Emenr9eHhoVFaWTJ0/q6NGjDokdN5mHHpLc3GwfK1f+OkXQMApfX93noYecGzeAEjk6P0nkKFSMgZISrno00a9/rFX55fXVfQZWWIQA7GX3SNaUKVN09OhR+fv7KyQkRO7u7jbte/bsuaFA0tLSJEkBAQE22wMCAnTixAlrHz8/P5tP+iwWi/z9/a3vT0tLU7169Yrto6itUaNGNxQnKoGpU6Wvv5YyMqSffy7cdvGibZ8rX1evLgUEFL4PgMtxdH6SyFGoGOskBUry0q+jV1f/oXbl60uSzv/yPgCuxe4i6+GHH3ZkHE6VmJhovWFlZmamk6OBwzVvLu3bJ40YIa1YUTg9sDQ1aki9e0v//Kfk5VVxMQKw262Sn6QyflehUsiU9A9JD0oKlVStjL4XJR2S9JkKiy0ArsXuImvy5MmOjEOBgYGSpPT0dDVo0MC6PT093doWGBiozMxMGYZh/aTQMAxlZGTY9ElPT7fZd9Hroj5Xi4+PV3x8vCTbm1qiEvPykj78UHrnHempp6S8vOJ9PDykv/1NGj264uMDYDdH5yfJeTnqyvxksQSbfFZwRZckLZUULukBlfyH2mVJqyXtrMC4AFwbl7kZcaNGjRQYGKikpCTrtp9//lmbN2+2zm/v0KGDcnJylJKSYu2TkpKi8+fP2/TZvHmzfi6aBiYpKSlJwcHBatiwYcWcDG4ebdoUFlMl8fCQ2rat2HgAuCRyFCpamgqLqZJclnSqAmMBcO3KLLJq1aqlrKwsSVLNmjWL3eDxWm/2mJOTo927d2v37t0qKChQamqqdu/erdTUVFksFo0bN05/+ctftGzZMu3bt0/Dhg2Tt7e3Bg8eLElq1qyZevbsqdGjRyslJUUpKSkaPXq0evfurbCwMEmFy+fWqFFDw4YN0759+7Rs2TK9+uqrrNqEku3cKV36ZaKFxVI4PbDo5+TSJWnHDufFBqBUZucniRwF1xKsX/9IK1Dh9MCCX15X+aUdgOsqc7rgzJkzrff+mDVr1g0fbMeOHbrvvvusrydPnqzJkycrLi5O8+bN06RJk5Sbm6sxY8YoOztb7dq105o1a6wxSNKiRYs0duxYRUVFSZL69OljE5uPj4+SkpI0ZswYRUREqHbt2powYYLGjx9/w/GjEtq8WcrN/XVxizfeKJw+mJFRuH3LFunxx50dJYCrmJ2fJHIUXEsDFV6TVbS4xSpJPVW4KEa1X9r5GBBwXWUWWXFxcSU+v15du3Yt8672FotFCQkJSkhIKLVP7dq1tWDBgjKP07JlS23atOl6w8St5KuvCpdn79Pn18UtuncvXBRj2bLCdgAux+z8JJGj4FrqqXDk6qB+XdziexUuinHXL+0AXJfLXJMFOEWzZoWLXyxe/OvqgVcuitG0qXPjAwDckjJVWFx9pF9XDyxaFOMzSVlOiguAfcocyapSpYrdc8Tz8/NNCQioUJ9/Xnrb8OGFDwAuh/yEym5RGW27f3kAcF1lFllLliyxJrH09HS9+OKL6tu3rzp06CCpcNWkTz75RC+99JLjIwUA4BfkJwCAKyuzyLryBo99+vTRK6+8olGjRlm3/e53v9M999yjTz75RE888YTjogQA4ArkJwCAK7P7mqz169fbrLpU5L777tPGjRvNjAkAALuRnwAArsbuIsvX11dLly4ttn3p0qXy8/MzNSgAAOxFfgIAuJoypwte6eWXX9bw4cO1YcMG65z3bdu2ae3atZozZ47DAgQAoCzkJwCAq7G7yBo6dKiaNm2qGTNm6LPPPpNUeHf7rVu3ql27dg4LEACAspCfAACuxq4i69KlSxoyZIimTZumhQsXOjomAADsQn4CALgiu67Jcnd315o1a+y+JwkAABWB/AQAcEV2L3zRr18/LVu2zJGxAABwzchPAABXY/c1WQ0aNNCUKVO0efNmRUREyMvLy6Z9/PjxpgcHVAYJ+fnODgGo1MhPAABXY3eRNW/ePNWuXVt79uzRnj17bNosFgtJDADgFOQnAICrsbvIOnLkiCPjAADgupCfAACuxu5rsq6Uk5Oj8+fPmx0LAAA3hPwEAHAF11Rk/eMf/1CDBg3k4+OjWrVqKSQkRLNnz3ZUbAAA2IX8BABwJXZPF5w2bZpeeeUVTZw4Uffee68kafPmzXr22Wd19uxZPfvssw4LEgCA0pCfAACuxu4i6+2331ZiYqIGDRpk3Xb//ferSZMm+uMf/0gSAwA4BfkJAOBq7J4umJGRobvvvrvY9nvuuUfp6emmBgUAgL3ITwAAV2N3kRUaGqpFixYV275o0SKFhYWZGhQAAPYiPwEAXI3d0wUTEhL06KOPatOmTYqMjJQkbd26VV9++aX+/e9/OyxAAADKQn4CALgau4usfv366auvvtKbb76pzz//XJLUrFkz/ec//1GbNm0cFiAAAGUhPwEAXI3dRZYkhYeHa8GCBY6KBQCA60J+AgC4kmsqsvLy8rRw4UJ98803slgsat68uQYNGiQPDw9HxQcAQLnITwAAV2L3whfffPONmjRpovHjx+urr77Stm3bNG7cOIWGhurAgQOOjBEAgFKRnwAArsbuIuvJJ59UmzZtlJqaqs2bN2vz5s1KTU3Vb37zG40bN86RMQIAUCryEwDA1dg9XXDr1q3avn27atWqZd1Wq1YtTZ06Ve3bt3dIcAAAlIf8BABwNXaPZFWvXl2nT58utv3MmTOqXr26qUEBAGAv8hMAwNXYXWTFxMRo1KhR2rp1q/Lz85Wfn68tW7Zo9OjR6tOnjyNjBACgVOQnAICrsbvImjFjhpo0aaJOnTqpevXqql69urp06aLQ0FBNnz7dkTECAFAq8hMAwNXYXWTddttt+vTTT/Xtt99q2bJlWrZsmQ4dOqSPP/5YPj4+du1j06ZN6tOnj+rWrSuLxaJ58+ZZ2y5duqRnnnlGrVq1kpeXl4KCgjR48GClpqba7KNr166yWCw2j4EDB9r0yc7OVmxsrHx8fOTj46PY2NgSp5IAAG5+5CcAgKuxq8hKTU3VxYsXJUl33nmnYmJiFBMTowYNGhRLMmXJyclRixYtNGPGDHl6etq0XbhwQbt27dLzzz+vXbt26dNPP9X//vc/9ezZU5cvX7bpO3z4cJ06dcr6eOedd2zaBw8erF27dmnVqlVatWqVdu3apdjYWLvjBADcHMhPAABXVO7qgh9++KGmTJmiXbt2FWszDEPR0dGaPHmyHnnkkXIPFh0drejoaEnSsGHDbNp8fHyUlJRks+2dd95R8+bNdeDAAbVs2dK6vUaNGgoMDCzxGAcOHNCqVau0ZcsWdejQwbqfTp066dChQwoLCys3TgCA6yM/AQBcVbkjWYmJiXr66adVrVq1Ym0eHh565pln9NZbbzkkuLNnz0qSateubbP9ww8/lK+vr5o3b66JEyfq3Llz1raUlBR5e3urY8eO1m2RkZHy8vJScnKyQ+IEAFQ88hMAwFWVO5J14MAB3XvvvaW2d+zYUU8//bSpQUnSxYsXNWHCBMXExKhevXrW7YMHD1ZISIiCg4O1f/9+Pffcc9qzZ4/WrFkjSUpLS5Ofn58sFov1PRaLRf7+/kpLSyvxWImJiUpMTJQkZWZmmn4uAADz3Wr5Sbpg+rkAAByj3CLrzJkz1vnuJcnLy7N+omeWy5cva8iQITp9+rQ+++wzm7b4+Hjr85YtW6px48Zq166ddu3apbZt217X8eLj4637jYiIuP7AAQAV5lbLTxZL8PUHDgCoUOVOF2zUqJG2b99eavt//vMfNWzY0LSALl++rEGDBmnPnj1at26d6tSpU2b/iIgIubm56fDhw5KkwMBAZWZmyjAMax/DMJSRkVHqPHkAwM2H/AQAcFXlFln9+vXT888/r1OnThVrO3nypP70pz+pf//+pgRz6dIlDRgwQHv27NGGDRvsSjp79+5Vfn6+goKCJEkdOnRQTk6OUlJSrH1SUlJ0/vx5m3nwAICbG/kJAOCqyp0u+Mwzz+jjjz9WaGiohgwZoqZNm0oqnAu/cOFCNWjQQJMmTbLrYDk5Ofruu+8kSQUFBUpNTdXu3bt1++23Kzg4WI888oi2b9+u5cuXy2KxWOeo+/j4yNPTU99//70WLlyo6Oho+fr66ptvvtGECRPUpk0bRUZGSpKaNWumnj17avTo0dZ57KNHj1bv3r1ZuQkAKhHyEwDAVZVbZHl7e2vr1q167rnn9K9//UvZ2dmSCm/+GBsbq6lTp6pmzZp2HWzHjh267777rK8nT56syZMnKy4uTgkJCfr0008lSeHh4Tbvmzt3roYNG6Zq1app3bp1mjFjhnJyclS/fn316tVLkydPlpubm7X/okWLNHbsWEVFRUmS+vTpo1mzZtkVIwDg5kB+AgC4qnKLLKnwk7rZs2frH//4h7KysmQYRrEVkuzRtWtXm7noVyurTZLq16+vL7/8stzj1K5dWwsWLLim2AAANx/yEwDAFdlVZBWxWCzy8/NzVCwAAFwX8hMAwJWUu/AFAAAAAMB+FFkAAAAAYCKKLAAAAAAwEUUWAAAAAJio3CLr/vvv17Jly0ptT09Pt1meFgCAikB+AgC4qnKLrA0bNujRRx/V5MmTS+1T3tK2AACYjfwEAHBVdk0XfOuttzR9+nT17dtX58+fL9Z+rfcjAQDADOQnAIArsqvIevDBB7Vt2zbt379f7du31w8//ODouAAAKBf5CQDgiuxe+KJZs2bavn276tevr7vvvltr1651ZFwAANiF/AQAcDXXtLqgj4+PVqxYoVGjRik6Olpvvvmmo+ICAMBu5CcAgCupWl6Hq+ezWywWvfrqq2rdurVGjhyp9evXOyw4AABKQ34CALiqckeySluZaeDAgdqyZYv27t1relAAAJSH/AQAcFXljmRt2LBBt99+e4ltrVu31s6dO7VixQrTAwMAoCzkJwCAqyq3yOrSpUuZ7XXq1NHQoUNNCwgAAHuQnwAAruqaFr4AAAAAAJSNIgsAAAAATESRBQAAAAAmosgCAAAAABNRZAEAAACAiSiyAAAAAMBEFFkAAAAAYCKKLAAAAAAwEUUWAAAAAJiIIgsAAAAATESRBQAAAAAmosgCAAAAABNRZAEAAACAiSiyAAAAAMBEFVpkbdq0SX369FHdunVlsVg0b948m/Zhw4bJYrHYPNq3b2/TJy8vT2PHjpWvr6+8vLzUp08fHT9+3KZPamqqYmJi5OXlJV9fX/3f//2fLl686OjTAwDcpMhPAAAzVWiRlZOToxYtWmjGjBny9PQssU+3bt106tQp62PlypU27ePGjdNHH32kxYsXa/PmzTp79qx69+6t/Px8SVJ+fr569eqlc+fOafPmzVq8eLGWLl2qCRMmOPz8AAA3J/ITAMBMVSvyYNHR0YqOjpZU+KlgSTw8PBQYGFhi25kzZzRnzhzNnTtX3bt3lyR98MEHCgkJ0dq1axUVFaU1a9Zo//79OnbsmOrXry9Jeu211zRy5EhNnTpVtWrVMv/EAAA3NfITAMBMLndN1pYtW+Tv76/Q0FCNGjVKGRkZ1radO3fq0qVL6tGjh3Vb/fr11axZMyUnJ0uSUlJS1KxZM2sCk6SoqCjl5eVp586dFXciAIBKhfwEALBXhY5kladnz57q16+fGjVqpKNHj+qFF17Qb3/7W+3cuVMeHh5KS0uTm5ubfH19bd4XEBCgtLQ0SVJaWpoCAgJs2n19feXm5mbtc7XExEQlJiZKkjIzMx1wZgCAm5kr5CfpgunnBQBwDJcqsgYOHGh93rJlS4WHhyskJEQrVqxQv379HHbc+Ph4xcfHS5IiIiIcdhwAwM3JFfKTxRLssOMAAMzlctMFrxQcHKx69erp8OHDkqTAwEDl5+crKyvLpl96erp1nnxgYKDS09Nt2rOyspSfn1/qXHoAAK4F+QkAUBaXLrKysrJ04sQJBQUFSZLCw8Pl7u6upKQka5/jx4/rwIED6tixoySpQ4cOOnDggM2yuUlJSfLw8FB4eHjFngAAoFIiPwEAylKh0wVzcnL03XffSZIKCgqUmpqq3bt36/bbb9ftt9+uhIQE9e/fX0FBQTp69Kiee+45+fv7q2/fvpIkHx8fjRgxQpMmTZK/v7/q1Kmj8ePHq1WrVurWrZskqUePHmrevLmGDh2q119/XT/++KOefvppjRo1ipWbAAAlIj8BAMxUoSNZO3bsUJs2bdSmTRvl5uZq8uTJatOmjV588UW5ublp7969evDBBxUaGqq4uDiFhYUpJSVFNWvWtO5j+vTp6tu3rwYMGKDIyEh5e3tr+fLlcnNzkyS5ublpxYoVqlGjhiIjIzVgwAD1799ff/vb3yryVAEANxHyEwDATBU6ktW1a1cZhlFq++rVq8vdh4eHh2bOnKmZM2eW2qdBgwb6/PPPrytGAMCth/wEADCTS1+TBQAAAAA3G4osAAAAADARRRYAAAAAmIgiCwAAAABMRJEFAAAAACaiyAIAAAAAE1FkAQAAAICJKLIAAAAAwEQUWQAAAABgIoosAAAAADARRRYAAAAAmIgiCwAAAABMRJEFAAAAACaiyAIAAAAAE1FkAQAAAICJKLIAAAAAwEQUWQAAAABgIoosAAAAADARRRYAAAAAmIgiCwAAAABMRJEFAAAAACaiyAIAAAAAE1FkAQAAAICJKLIAAAAAwEQUWQAAAABgIoosAAAAADARRRYAAAAAmIgiCwAAAABMRJEFAAAAACaq0CJr06ZN6tOnj+rWrSuLxaJ58+bZtFsslhIfY8aMsfYZNmxYsfb27dvb7CcvL09jx46Vr6+vvLy81KdPHx0/frwiThEAcBMiPwEAzFShRVZOTo5atGihGTNmyNPTs1j7qVOnbB7Lly+XJD366KM2/bp162bTb+XKlTbt48aN00cffaTFixdr8+bNOnv2rHr37q38/HzHnRwA4KZFfgIAmKlqRR4sOjpa0dHRkgo/8btaYGCgzetPP/1UoaGh6tKli812Dw+PYn2LnDlzRnPmzNHcuXPVvXt3SdIHH3ygkJAQrV27VlFRUSacCQCgMiE/AQDM5LLXZOXk5OjDDz/UqFGjirVt2bJF/v7+Cg0N1ahRo5SRkWFt27lzpy5duqQePXpYt9WvX1/NmjVTcnJyhcQOAKi8yE8AgPJU6EjWtVi0aJEuXryouLg4m+09e/ZUv3791KhRIx09elQvvPCCfvvb32rnzp3y8PBQWlqa3Nzc5Ovra/O+gIAApaWllXisxMREJSYmSpIyMzMdc0IAgErBWflJuuCI0wEAOIDLFlnvvvuuHnzwQfn5+dlsHzhwoPV5y5YtFR4erpCQEK1YsUL9+vW7rmPFx8crPj5ekhQREXH9QQMAKj1n5SeLJfj6gwYAVCiXnC64e/du7dixo8SpGFcLDg5WvXr1dPjwYUmF8+bz8/OVlZVl0y89Pb3UefIAANiD/AQAsIdLFlmJiYlq1KiRunXrVm7frKwsnThxQkFBQZKk8PBwubu7Kykpydrn+PHjOnDggDp27OiwmAEAlR/5CQBgjwqdLpiTk6PvvvtOklRQUKDU1FTt3r1bt99+uxo0aCBJunDhghYuXKhJkybJYrEUe39CQoL69++voKAgHT16VM8995z8/f3Vt29fSZKPj49GjBihSZMmyd/fX3Xq1NH48ePVqlUru5IiAODWQ34CAJipQkeyduzYoTZt2qhNmzbKzc3V5MmT1aZNG7344ovWPv/61790/vx5DR8+vNj73dzctHfvXj344IMKDQ1VXFycwsLClJKSopo1a1r7TZ8+XX379tWAAQMUGRkpb29vLV++XG5ubhVyngCAmwv5CQBgpgodyeratasMwyizz/Dhw0tMYJLk6emp1atXl3scDw8PzZw5UzNnzryuOAEAtxbyEwDATC55TRYAAAAA3Kxcdgl3ABWr9W9+4+wQALiw8R7OjgC3sozLhf+O5y9X3CT4UQUgSZr+5pvODgEAgBK9yl+suMkwXRAAAAAATMTnAgAAoES1hjk7AgC4OTGSBQAAAAAmosgCAAAAABNRZAEAAACAibgmCwAAFIp1dgAAUDkwkgUAAAAAJqLIAgAAAAATMV0QAIBb2K7IZs4OAQAqHUayAAAAAMBEFFkAAAAAYCKmCwIAcIv5WA85OwQAqNQYyQIAAAAAE1FkAQAAAICJmC4IAMAt4CM97OwQAOCWwUgWAAAAAJiIIgsAAAAATMR0QQAAKqkDn7Z1dggAcEtiJAsAAAAATESRBQAAAAAmYrogAACVyT+dHQAAgJEsAAAAADARRRYAAAAAmIgiCwAAAABMxDVZAADczD4/6+wIAABXYSQLAAAAAExUoUXWK6+8orvvvlu1atWSn5+fYmJitG/fPps+hmEoISFBwcHB8vT0VNeuXbV//36bPtkaT+2fAAAgAElEQVTZ2YqNjZWPj498fHwUGxur06dP2/TZu3evunTpIk9PT9WtW1cvv/yyDMNw+DkCAG4+5CcAgJkqdLrgxo0b9cQTT+juu++WYRh68cUX1a1bN33zzTe6/fbbJUmvvfaaXn/9dc2bN09hYWF6+eWX1b17dx06dEg1a9aUJA0ePFipqalatWqVJGnkyJGKjY3V8uXLJUlnz55V9+7d1blzZ23fvl0HDx7U8OHD5eXlpQkTJlTkKQMAbgI3X356w9TzBwCYy2I48eOznJwc+fj46JNPPlFMTIwMw1BwcLD+8Ic/6Pnnn5ck5ebmyt/fX3/72980evRoHThwQHfddZe2bNmiyMhISdKWLVvUqVMnHTx4UGFhYXrrrbf0zDPPKD09XZ6enpKkKVOm6K233tLx48dlsVhKjSkiIkI7/vMfx588ALiYiHvu0Y4dO5wdhktwxfxksQRLinf4uQOAKwoP//ymylFOvSbr3LlzKigoUO3atSVJR44cUVpamnr06GHt4+npqc6dOys5OVmSlJKSIm9vb3Xs2NHaJzIyUl5eXjZ9OnXqZE1gkhQVFaWTJ0/q6NGjFXBmAICbGfkJAHAjnLq64JNPPqnWrVurQ4cOkqS0tDRJUkBAgE2/gIAAnThxwtrHz8/P5tM+i8Wi/2/v3uOiqtY+gP8GGGAcxtvhJqAooXjPhBRKxbIXLbyUXbRM0cyy1DSt1E5HqCyzst6Ox0tagKnpSa20Yy9FhgZhmoSKomniJUtARBNFUJjn/UPZh839socZ5Pf9fObj3muvvfazFuN+WLBn4e7urpyfmZkJHx+fcm2UHOvQoYPq2IoVK7BixQoAwIEDBxDUp49WXWxyzp49Czc3N2uH0Whx/OqH41c/hw8ftnYINsMW85OTUy66d/+PVl1scnh/qB+OX/1xDOunseUoq02yZs6ciaSkJCQlJcHe3t5aYQAAnnrqKTz11PVHMIKCghrVryJtDcevfjh+9cPxq5+goCBrh2ATmJ9uThy/+uH41R/HsH4aW46yyuOCzz//PNatW4fvv/8efn5+SrmnpycAICsrS1U/KytLOebp6YmzZ8+qVmISEWRnZ6vqVNRG6WsQERGVxfxERERaaPBJ1vTp05UE1rlzZ9WxDh06wNPTE/Hx8UpZQUEBEhMTlWfcQ0JCcOnSJezcuVOps3PnTly+fFlVJzExEQUFBUqd+Ph4eHl5oX379hbsHRERNVbMT0REpBX7qKioqIa62JQpU7Bq1Sps2LAB7dq1w6VLl3Dp0iUAgKOjI3Q6HYqLi/HWW28hICAAxcXFmDlzJs6cOYMVK1bAyckJbm5u2LVrFz799FPcdtttOH36NJ5++mn06dMH06ZNAwB06tQJH374Ifbu3YvOnTsjKSkJL7zwAubOnav6QHJlAgMDLToONzuOX/1w/OqH41c/TXX8mJ+aBo5f/XD86o9jWD+NavykAQGo8BUZGanUMZvNEhkZKZ6enuLk5CQDBgyQtLQ0VTu5ubkyZswYMZlMYjKZZMyYMXL+/HlVnf3790v//v3FyclJPD09JSoqSsxmc0N0k4iIGhnmJyIi0pJV/04WERERERHRzcaqfyeLiIiIiIjoZsNJ1g1Lly5Fhw4d4OzsjMDAQCQmJlo7JJsQFRUFnU6nepVeAUtEEBUVBS8vLxgMBgwcOBAHDx5UtXH+/HmMHTsWLVq0QIsWLTB27FhcuHChobvSIH744QcMHz4c3t7e0Ol0iI2NVR3XarzS0tIQGhoKg8EAb29vvPbaa7gZfild3fiNHz++3PsxODhYVaewsBDTpk2Dq6srjEYjhg8fjtOnT6vqnDp1CsOGDYPRaISrqyuee+45XL161dLds7gFCxbg9ttvR/PmzeHm5oZhw4bhwIEDqjp8DzZOzFHlMT/VDvNT/TFH1V2TzE/WeUrRtqxfv14cHBxkxYoVkp6eLlOnThWj0SgnT560dmhWFxkZKQEBAXLmzBnllZ2drRx/6623xMXFRTZu3ChpaWny8MMPS5s2beTixYtKnSFDhkjXrl0lOTlZkpOTpWvXrjJ06FBrdMfitm7dKnPnzpUNGzaIwWCQmJgY1XEtxuuvv/4SDw8PefjhhyUtLU02bNggLi4u8u677zZUNy2muvGLiIiQe+65R/V+PHfunKrO5MmTpU2bNvLtt99KSkqKhIaGyq233ipFRUUiIlJUVCTdu3eX0NBQSUlJkW+//VbatGkjU6dObahuWkxYWJhER0dLWlqa7N+/X+6//37x8PBQjRHfg40Pc1TFmJ9qh/mp/pij6q4p5idOskSkT58+8uSTT6rK/P39Zc6cOVaKyHZERkZKt27dKjxmNpvF09NT5s+fr5Tl5+eLi4uLLF++XERE0tPTBYAkJSUpdRITEwWAHD582LLBW5nRaFTdgLUar6VLl4rJZJL8/Hylzuuvvy5eXl431Yfny46fyPUEFh4eXuk5Fy5cEL1eL2vWrFHKTp06JTqdTuLi4kRE5OuvvxadTienTp1S6qxevVqcnJzkr7/+0rYTVpaXlyd2dnayZcsWEeF7sLFijqoY81PdMT/VH3NU/TSF/NTkHxe8evUqUlJSEBYWpioPCwtDcnKylaKyLRkZGfDy8kKHDh0wevRoZGRkAACOHz+OzMxM1dgZDAYMGDBAGbudO3fCxcVFtTTxnXfeCaPR2OTGV6vx2rlzJ/r37w+DwaDUGTx4MP7880+cOHGiYTpjRUlJSXB3d0enTp0wadIkZGdnK8dSUlJw7do11Ri3bdsWXbp0UY1fly5d0LZtW6XO4MGDUVhYiJSUlIbrSAPIy8uD2WxGq1atAPA92BgxR1WN+UkbvDdohzmqZppCfmryk6ycnBwUFxfDw8NDVe7h4YHMzEwrRWU7+vbti9jYWMTFxWHlypXIzMzEHXfcgXPnzinjU9XYZWZmws3NDTqdTjmu0+ng7u7e5MZXq/HKzMyssI3S17hZDRkyBJ988gm2bduGRYsWYffu3bj77rtRWFgI4Hr/7e3t4erqqjqv7BiXHT9XV1fY29vfdOM3ffp09OrVCyEhIQD4HmyMmKMqx/ykHd4btMEcVXNNIT85NOjVqNG59957VfvBwcHw8/PDqlWryn2Yk8jSRo8erWz36NEDgYGB8PX1xdatWzFy5EgrRmZ7Zs6ciaSkJCQlJcHe3t7a4RBpjvmJbA1zVM00lfzU5H+TVfLTgaysLFV5VlaWapUius7FxQXdunXD0aNHlfGpauw8PT1x9uxZ1aouIoLs7OwmN75ajZenp2eFbZS+RlPh5eUFHx8fHD16FMD1/hcXFyMnJ0dVr+wYlx2/kt8W3Czj9/zzz2PdunX4/vvv4efnp5TzPdj4MEfVHPNT3fHeYBnMUeU1pfzU5CdZjo6OCAwMRHx8vKo8Pj5e9cwnXVdQUIDDhw+jTZs26NChAzw9PVVjV1BQgMTERGXsQkJCcOnSJezcuVOps3PnTly+fLnJja9W4xUSEoLExEQUFBQodeLj4+Hl5YX27ds3TGdsRE5ODv744w+0adMGABAYGAi9Xq8a49OnT+PQoUOq8Tt06JBqydz4+Hg4OTkhMDCwYTtgAdOnT1cSWOfOnVXH+B5sfJijao75qe54b7AM5ii1JpefGnSZDRu1fv160ev1snLlSklPT5fnnntOjEajnDhxwtqhWd2sWbNk+/btkpGRIT/99JOEh4eLyWRSxuatt96S5s2by6ZNmyQtLU1GjRpV4XKb3bt3V5bb7N69+027RG5eXp6kpqZKamqqGAwGefXVVyU1NVVZalmL8bpw4YJ4eHjIqFGjJC0tTTZt2iQmk+mmWCK3qvHLy8uTWbNmSXJyshw/flwSEhIkODhYvL29VeM3efJk8fb2lvj4ePnll19k4MCBFS6Pe9ddd8kvv/wi8fHx4uXl1eiXxxURefbZZ8VkMsm2bdtUSwjn5eUpdfgebHyYoyrG/FQ7zE/1xxxVd00xP3GSdcOSJUvE19dXHB0dpXfv3rJjxw5rh2QTSt7ger1evLy8ZOTIkXLw4EHluNlslsjISPH09BQnJycZMGCApKWlqdrIzc2VMWPGiMlkEpPJJGPGjJHz5883dFcaREJCggAo94qIiBAR7cZr//790r9/f3FychJPT0+Jioq6KZbHrWr88vPzJSwsTNzc3ESv10u7du0kIiJCtcytiEhBQYFMnTpVWrduLQaDQYYOHVquzsmTJyU8PFwMBoO0bt1apk2bJgUFBQ3ZVYuoaOwASGRkpFKH78HGiTmqPOan2mF+qj/mqLprivlJJ3KT/BluIiIiIiIiG9DkP5NFRERERESkJU6yiIiIiIiINMRJFhERERERkYY4ySIiIiIiItIQJ1lEREREREQa4iSLiIiIiIhIQ5xkERERERERaYiTLKIGtHHjRuh0OmU/NjYWLi4uVoyoZiZMmIDXXnutXm2kpaXB29sbly9f1igqIiLSEnMUcxRph5MssmlmsxkDBgzAsGHDVOX5+fkICAjA5MmTqzw/Ly8P//jHP9C1a1cYDAZ4eHhg4MCBWLduHcxmsyVDr5FRo0YhIyND83Z1Oh02btyoSVtpaWn48ssvMWPGjHq106NHDwQHB+O9997TJC4iImtjjqob5ihqCjjJIptmZ2eH2NhYJCQkIDo6WimfPXs2iouLsWjRokrPvXDhAkJCQhAdHY0XX3wRe/bsQVJSEiIiIvD666/j1KlTFov76tWrNapnMBjg7u5usTi0sHjxYjz44INo3rx5vduaMGECli1bhqKiIg0iIyKyLuYo62OOIpslRI3AsmXLpHnz5nLixAn57rvvxN7eXhITE6s855lnnpFmzZrJ77//Xu7YlStX5MqVKyIikpubK+PGjZOWLVuKs7OzDBo0SA4cOKCqv2nTJunevbs4OjqKj4+PzJ8/X8xms3Lc19dXIiMjZcKECdKiRQt56KGHRERk1apV0q5dOzEYDBIeHi7/+te/pPR/u5iYGDEajcp+ZGSkdOvWTdatWyd+fn7i4uIiI0aMkLNnzyp1du/eLf/zP/8jf/vb38RkMsmdd94pycnJqlgAKC9fX1/l2JYtW6R3797i5OQk7du3l5dfflkKCwsrHcOioiJp0aKFfPnll6pyX19fefXVVyUiIkJcXFzEx8dH1q9fL+fPn5dRo0aJ0WgUf39/+eabb1TnFRYWipOTk8THx1d6TSKixoY5ijmKqCxOsqjRCAsLk/79+0vbtm3lxRdfrLJucXGxtGrVSiZNmlRtu8OHD5eAgADZsWOH7N+/X4YNGyY+Pj6Sn58vIiJ79uwROzs7mTdvnvz666+yZs0aMRqN8s9//lNpw9fXV0wmkyxcuFCOHj0qR44ckZ9++kl0Op3Mnz9ffv31V1m+fLm0bt262gRmNBrl/vvvl3379klycrK0a9dOnnrqKaXOtm3b5JNPPpH09HQ5dOiQTJkyRVq2bCk5OTkiIpKdnS0AZOXKlXLmzBnJzs4WEZG4uDgxmUwSHR0tv/32m3z//ffSqVMnmTVrVqVj88svvwgAOX36tKrc19dXWrVqJUuWLJEjR47IzJkzxcnJSe69915ZtWqVHD16VJ544glxc3NTvlEo0bdvX3nllVeq/boQETUmzFHXMUcRXcdJFjUaGRkZotPpxN/fXwoKCqqsm5WVJQDkvffeq7LekSNHBIDs2LFDKbtw4YI0b95cVq5cKSIijz32mNx1112q8yIjI8Xb21vZ9/X1laFDh6rqPProo3LPPfeoyiZOnFhtAnNycpILFy4oZfPnz5dbbrml0j6YzWbx9PSU1atXK2UAZMOGDap6/fv3l9dee01V9sUXX4jRaFT9xLPscZ1OJ8XFxapyX19fGT16tLKfl5cnAGTatGlK2fHjxwWA/Pzzz6pzH3jgAXn88ccr7Q8RUWPEHFUx5ihqqviZLGo0oqOjYTAYcPr0aRw/frzKuiJSozYPHToEOzs7hISEKGUtWrRAjx49kJ6ertS58847Vef169cPf/zxBy5evKiUBQUFlWu7dLsAyu1XxNfXFy1atFD2vby8kJ2drexnZ2fj6aefRqdOndCiRQuYTCZkZ2dX+/x+SkoK3njjDbi4uCivxx57DJcvX0ZmZmaF51y5cgV6vR52duVvFT179lS2XVxc0KxZM/To0UMp8/DwUOItzWAw4MqVK1XGSkTU2DBHXcccRXSdg7UDIKqJn3/+GW+99Ra2bNmCZcuWISIiAsnJybC3t6+wvpubG1q2bIlDhw7V+Zqll7GtSR2j0Vjna5Wm1+vLXaP0KlMRERHIysrC+++/j/bt28PJyQmDBg2q9oPMZrMZkZGRePjhh8sdc3Nzq/AcV1dXXL16Ffn5+WjWrFm1cZYuKxmbsitk5ebmon379lXGSkTUmDBHMUcRlcXfZJHNKygowLhx4zB+/Hjce++9WLFiBX777Te8/fbblZ5jZ2eH0aNHY+3atTh9+nSFbRYUFKBLly4wm83YuXOncuzixYtIS0tD165dAQBdunTBjz/+qDo/KSkJPj4+MJlMlcbQpUsX/PTTT6qysvt1kZSUhGnTpiE8PBzdunWDyWTCmTNnVHX0ej2Ki4tVZb1798bhw4fh7+9f7uXgUPHPW3r16gUAyk9MtXDgwAH07t1bs/aIiKyJOUqNOYroOk6yyObNnTsXBQUFyt+u8PT0xJIlSxAVFYWDBw9Wet4bb7yBdu3aoW/fvoiJicHBgwfx22+/YfXq1QgMDERmZiY6duyIESNG4Omnn0ZiYiLS0tLw+OOPo3nz5njssccAALNmzcKOHTsQFRWFI0eOYO3atVi0aBFeeumlKuN+7rnn8N1332HBggU4evQoVq5ciS+++KLe49GpUyesWbMG6enp+PnnnzF69Gg4Ojqq6rRv3x7btm1DZmYmzp8/DwCYN28ePv30U8ybNw8HDhzA4cOHsXHjxir74ebmht69eyMpKanecQPAiRMn8McffyAsLEyT9oiIrI05So05iugGa38ojKgqO3bsEHt7e0lISCh37KGHHpLAwEC5du1apedfuHBBXn75ZQkICBAnJydxc3OT0NBQWbdunfJB2dosj6vX6ytdHvedd94pd/3o6Ghp27atODs7y5AhQ2Tx4sU1Wh63tLJ19u7dK3369BFnZ2fx8/OTTz75RLp16yaRkZFKnS1btoi/v784ODiolsf95ptvpF+/fmIwGMRkMklgYKAsXry40vETEVm+fLkEBQWpyirqr9FolJiYGGX/ypUrAkC++uorpezNN9+UwYMHV3k9IqLGgjmKOYqoMjqRGn76koiapIKCAnTu3BmrV69G//7969xOYWEhOnbsiHXr1pX7kDYREVFdMEeRreLjgkRUJWdnZ3zyySfIzc2tVzsnT57E3//+dyYvIiLSDHMU2Sr+JouIiIiIiEhD/E0WERERERGRhjjJIiIiIiIi0hAnWURERERERBriJIuIiIiIiEhDnGQRERERERFpiJMsIiIiIiIiDXGSRUREREREpCFOsoiIiIiIiDTESRYREREREZGGOMkiIiIiIiLSECdZREREREREGuIki4iIiIiISEOcZBEREREREWmIkywiIiIiIiINcZJFRERERESkIU6yiIiIiIiINMRJFhERERERkYY4ySIiIiIiItIQJ1lEREREREQa4iSLiIiIiIhIQ5xkERERERERaYiTLCIiIiIiIg1xkkVERERERKQhTrKIiIiIiIg0xEkWERERERGRhjjJIiIiIiIi0hAnWURERERERBriJIuIiIiIiEhDnGQRERERERFpyMHaAdSVTucO4CoAXUlJmX/rU1bRsYq2rVVP67arb0NXqlin++9+6W1LHKuoTlXHtL6uNWKCyPVX2e2aHCtdXtWxsmX1va6WMWl13epiKtmtZLtsWWX1yp5T0zZq25alj1my/Z6DByMuLg5NhU7nD+BKyV6pF2qwXVluqqxe2W3Usl5FHaiiibLNl92u67GK/q1ptyx93RrFJKXKBDqdKMdKtnU36uhu/M/QlfrfooOU24dSr3zdqv6trD0t2tC6vbq0UdW41Lo9Eej+2yR0/z2EUpeq/N+a1KvoplhReU2O1aaNsuXVtVGb2GsRh7JZNi1XEXrZ7arCrUm96tqwVI5qtJOs6xOsgQDsb+zb33jZldmvaNuuDscqKi97XmXXrS4muwqO1bW9ulyr8j7qdNfbc3C4/qpsu67HrFXPFmOyg/l6QVHR9Vdl21Uds1a9RhiTiMAMlIx6tduVHSu+sS01aKNs3equVYz/JoHqYqqq7YquVVxJ3cr6WFUc1V3rq5wcNC35AKbc2NbjeqrV39ivbrum9UrSd1XbVR0ro2SCYF/mUg74b1oo3URF5Q5VtFG6XOs2qjqnonpVXbei7cr6r2wL4FB8/QXAzr4IDvrr2/YOxbB3KNkugoNDMeztbuyjGPYo2S6CA9Tl9ii6canS9cqfoz6mPqeq9kufU/GximKqrr2KYqp5vNXHpFEfi29sFxXBodgM+xspwr4I0JXcTItuvEq2iyvYr2u94krOK12vqvMsVa+i82pSr2zbpc67Vir9XisutV3032rXbrxKmii7XbpeZdvV1ausfcByOYqPCxIREREREWmIkywiIiIiIiINcZJFRERERESkIU6yiIiIiIiINMRJFhERERERkYY4ySIiIiIiItIQJ1lEREREREQa4iSLiIiIiIhIQ5xkERERERERaYiTLCIiIiIiIg1xkkVERERERKQhTrKIiIiIiIg0xEkWERERERGRhjjJIiIiIiIi0hAnWURERERERBriJIuIiIiIiEhDDtYOoK66dvWGwXDa2mHUytmzZ+Hm5mbtMGpNq7jNZuDq1evbJf9aSlMf64bWYHHb3fi5kKPj9Vc9cKwbzpUrV6wdQoPq2rUVDIYvrB1Gg1K9LwVAyT3ewvd6a7LG/8XiGy9rsI17j/2NV8OwjT43rGr7XPIlqF8KtimWylGNdpJlMBiwZ88ea4dRK0FBQY0uZqBxxt0YYwYYd0NqjDEDjTPuoKAga4fQoBpjfqqvxvi+rK+m1uem1l+AfW4qLJWj+LggERERERGRhjjJIiIiIiIi0pB9VFRUlLWDqKvAwEBrh1BrjTFmoHHG3RhjBhh3Q2qMMQONM+7GGHN9NLX+AuxzU9DU+guwz02FJfqsExHRvFUiIiIiIqImio8LEhERERERaYiTLCIiIiIiIg3Z5CRrxYoVuOuuu9CyZUvodDqcOHGiRudt2rQJXbt2hZOTE7p27YovvlD/nRIRQVRUFLy8vGAwGDBw4EAcPHhQk5gLCwsxbdo0uLq6wmg0Yvjw4Th9uuq/49W+fXvodLpyr/DwcKVOVFRUueOenp6axFzXuGsSkyXHuq5xL1iwALfffjuaN28ONzc3DBs2DAcOHFDVGT9+fLm+BQcH1ynGpUuXokOHDnB2dkZgYCASExOrrL9jxw4EBgbC2dkZfn5+WL58eb3btHTcn3/+OcLCwuDm5gaTyYS+fftiy5YtqjqxsbEVvs8LCgqsFvf27dsrjOnw4cOqetXdUxoy5oremzqdDkajsdb9qqsffvgBw4cPh7e3N3Q6HWJjY6s9Jy0tDaGhoTAYDPD29sZrr72Gsk+pW3qcLc1SOcuWWSp32BJL3MNtmSXuobbMUvczW1Xb/p44caLCr3FcXFzDBKyBmnzfVxHNvs5ig95//31588035f333xcAcvz48WrPSU5OFnt7e5k/f76kp6fL/Pnzxd7eXn766SelzltvvSUuLi6yceNGSUtLk4cffljatGkjFy9erHfMkydPljZt2si3334rKSkpEhoaKrfeeqsUFRVVek52dracOXNGef3yyy+i0+kkNjZWqRMZGSkBAQGqetnZ2fWOtz5x1yQmS451XeMOCwuT6OhoSUtLk/3798v9998vHh4ecu7cOaVORESE3HPPPaq+lT5eU+vXrxcHBwdZsWKFpKeny9SpU8VoNMrJkycrrJ+RkSHNmjWTqVOnSnp6uqxYsUIcHBxk48aNdW6zLmp7jeeee04WLFggu3btkqNHj0pUVJTY2dnJDz/8oNSJiYmRZs2aqcb0zJkzmsVcl7gTEhIEgBw8eFAVU+n3T03uKQ0Z84ULF8qNoZ+fn4wfP75W/aqPrVu3yty5c2XDhg1iMBgkJiamyvp//fWXeHh4yMMPPyxpaWmyYcMGcXFxkXfffVepY+lxbgiWylm2zFK5w1ZY4h5uyyxxD7V1lrif2bLa9vf48eMCQOLi4lRf48LCwoYJWAM1+b6vLC2/zjY5ySrx888/1zhhPfLII3LPPfeoygYNGiSjR48WERGz2Syenp4yf/585Xh+fr64uLjI8uXL6xXnhQsXRK/Xy5o1a5SyU6dOiU6nk7i4uBq3M3/+fGnRooXk5+crZZGRkdKtW7d6xVeZusZdXUyWHOv6xF1WXl6e2NnZyZYtW5SyiIgICQ8Pr3eMffr0kSeffFJV5u/vL3PmzKmw/ksvvST+/v6qsokTJ0pwcHCd26wLLa5x++23y8yZM5X9mJgYMRqNmsVYkdrGXfINwtmzZytts7p7Sn3Vd6yTkpIEgPz4449KWU36pRWj0Vhtkl66dKmYTCbVPe31118XLy8vMZvNImL5cW5IWuYsW2ap3GFLLHEPt2WWuIc2JlrdzxqLmvS3ZJL1888/N0xQDaCi7/vK0vLrbJOPC9bFzp07ERYWpiobPHgwkpOTAQDHjx9HZmamqo7BYMCAAQOUOnWVkpKCa9euqdpu27YtunTpUuO2RQQff/wxHn/8cRgMBtWxjIwMeHl5oUOHDhg9ejQyMjLqFa8WcVcVkyXHur5xl5aXlwez2YxWrVqpypOSkuDu7o5OnTph0qRJyM7OrlV8V69eRUpKSrn3Y1hYWKXxVfb+3bNnD65du1anNmtLq2vk5eWVG9MrV6XQjrAAABdsSURBVK7A19cXPj4+GDp0KFJTUzWJGahf3EFBQWjTpg0GDRqEhIQE1bHq7inWirnEypUr0a1bN9xxxx3ljlXVr4a0c+dO9O/fX3VPGzx4MP7880/lkTpLjrMta8z9tlTusBWWuIfbMkvdQ282Nbmf3YxGjhwJd3d33Hnnndi4caO1w6mXyr7vK03Lr/NNM8nKzMyEh4eHqszDwwOZmZnK8ZKyyurU59r29vZwdXWtc9vx8fE4fvw4Jk2apCrv27cvYmNjERcXh5UrVyIzMxN33HEHzp07V6+Y6xN3dTFZcqzrE3dZ06dPR69evRASEqKUDRkyBJ988gm2bduGRYsWYffu3bj77rtRWFhY43ZzcnJQXFxcq/5X9v4tKipCTk5OndqsLS2usWTJEpw+fRpjx45VygICAhAdHY3Nmzdj3bp1cHZ2xp133omjR49aLe42bdpg2bJl2LRpEz7//HMEBARg0KBBqs8gVHdPaeiYS/vrr7/w2Weflbtf1KRfDamyMSw5VlUdrd7Xtqox99tSucNWWOIebsssdQ+92dTkfnYzcXFxwbvvvovPPvsMX3/9NQYNGoRRo0ZhzZo11g6tzir6vq8sLb/ODrUPsW5eeeUVvPHGG1XWSUhIwMCBAxsmoBqoacxaWLlyJW6//XbceuutqvJ7771XtR8cHAw/Pz+sWrUKM2fOrLAtS8ddl5hqoiHHe+bMmUhKSkJSUhLs7e2V8tGjRyvbPXr0QGBgIHx9fbF161aMHDlSk2vfrDZt2oQXX3wR//73v+Hr66uUh4SEqG5od9xxB3r16oXFixfjn//8pzVCRUBAAAICApT9kJAQnDhxAu+88w769+9vlZhqY82aNTCbzarJLND4+2VLGmPOqq/Gmjuo4fFec/NzdXXFrFmzlP2goCDk5OTg7bffxuOPP27FyOqmsu/7LKnBJlkzZsyo9ovSrl27Orfv6emJrKwsVVlWVpayclHJv1lZWarrlK5T15h/+uknFBcXIycnB25ubqq2a3Kzyc7OxubNm7FkyZJq67q4uKBbt25V/hagoeKuLKa6jHVDxv38889j/fr1SEhIgJ+fX5V1vby84OPjU6vfuri6usLe3r7K92NZlb1/HRwc4OrqChGpdZu1VZe4S2zcuBHjxo3DJ598gmHDhlVZ197eHkFBQZr9Jqs+cZfWt29frF+/Xtmv7p5SH/WNeeXKlXjwwQfRunXrauuW7VdDqmwMS45VVcfaq85ZO2dZg7Vzh62wxD3cllnqHnqzqcn97GbXt29fxMTEWDuMWqvN931afp0b7HFBV1dXdO7cucpXs2bN6tx+SEgI4uPjVWXx8fHK5xU6dOgAT09PVZ2CggIkJiZW+JmG2sQcGBgIvV6vavv06dM4dOhQpW2XFhsbCycnJzz66KPV1i0oKMDhw4fRpk2bSus0VNyVxVSXsW6ouKdPn45169bh+++/R+fOnavtW05ODv74448qx7ssR0dHBAYGVvl+LKuy929QUBD0en2d2qytul7js88+w9ixYxEbG4uHHnqo2uuICPbv31+rMa2KVmOzd+9eVUzV3VPqoz4x7969G/v27Sv3qGBlyvarIYWEhCAxMVG1XH98fDy8vLzQvn17pY4l39d1Ze2cZQ3Wzh22whL3cFtmqXvozaYm97ObXWP8Gtf2+z5Nv851W5/Dss6cOSOpqamydu1aASBbt26V1NRU1ZKLd999t2rVmx9//FHs7e1lwYIFcujQIXnzzTfFwcGh3BLuzZs3l02bNklaWpqMGjVK0yXcvb29JT4+Xn755RcZOHBgueVsAwICZPHixarzzGazdOzYsdyqPiVmzZol27dvl4yMDPnpp58kPDxcTCaTnDhxot4x1zXumsRkybGua9zPPvusmEwm2bZtm2o50ry8PBG5vurMrFmzJDk5WY4fPy4JCQkSHBws3t7etY57/fr1otfrZeXKlZKeni7PPfecGI1GZYzGjh0rY8eOVeqXLP87ffp0SU9Pl5UrV4pery+3hHtVbWqhtnGvW7dOHBwc5H//938rXfY+KipK4uLi5NixY5KamioTJkwQBwcH2bVrl9Xifv/99+WLL76QI0eOyIEDB2TOnDkCQDZt2qTUqck9pSFjLjFx4kTp2LFjhW3WpF/1kZeXJ6mpqZKamioGg0FeffVVSU1NVZZ5njNnjtx9991K/QsXLoiHh4eMGjVK0tLSZNOmTWIymVRL4Vp6nBuCpXKWLbNU7rAVlriH2zJL3ENtnSXuZ7astv2NjY2VtWvXSnp6uhw+fFjeeecd0ev18t5771mrC7VW3fd9Ipb9OtvkJCsyMlIAlHuVXm7S19dXIiIiVOdt2LBBAgICRK/XS+fOncv9ZzebzRIZGSmenp7i5OQkAwYMkLS0NE1iLigokKlTp0rr1q3FYDDI0KFD5dSpU6o6ACQyMlJV9v333wuASr/hLJmc6PV68fLykpEjR8rBgwc1ibmucdckJkuOdV3jrug9VbpOfn6+hIWFiZubm+j1emnXrp1ERESUa7emlixZIr6+vuLo6Ci9e/eWHTt2KMdCQ0MlNDRUVX/79u1y2223iaOjo7Rv316WLVtWqza1Upu4Q0NDKxzT0nVmzJgh7dq1E0dHR3Fzc5OwsDBJTk62atwLFy4Uf39/cXZ2llatWkm/fv1k69at5dqs7p7SkDGLiFy8eFGMRqMsXLiwwvZq2q+6Klm2ueyr5F4cEREhvr6+qnP2798v/fv3FycnJ/H09JSoqKhyy+BaepwtzVI5y5ZZKnfYEkvcw22ZJe6htsxS9zNbVdv+xsbGSpcuXaRZs2ZiMpkkMDBQVq9ebZ3g66i67/tELPt11t0IgoiIiIiIiDRw0yzhTkREREREZAs4ySIiIiIiItIQJ1lEREREREQa4iSLiIiIiIhIQ5xkERERERERaYiTLCIiIiIiIg1xkkU3raioKHTv3l2z9mJjY+Hi4qJZe2W9+uqreOKJJyzWvi3Jzs6Gm5sbTp8+be1QiIganfHjx2Po0KHWDqNC58+fh4eHB44dO1brc9PS0uDt7Y3Lly9bIDKihsVJFtmc4cOHY9CgQRUeO3ToEHQ6Hb799tsGjgoYNWoUMjIylH0tJ3HZ2dlYtGgRXnnlFU3au3r1Klq0aIG9e/dq0p7W3N3dMW7cOERGRlo7FCKiJq+goADjx49Hz549odfrMXDgwArr7dixA4GBgXB2doafnx+WL19ers6bb76J++67D7fccouq/NNPP0VISAhcXFxgNBrRt29frFmzRlWnR48eCA4OxnvvvadZ34ishZMssjkTJ05EQkICTpw4Ue7Yxx9/DF9fX9xzzz0NHpfBYIC7u7tF2v7oo4/Qp08f+Pn5adJeQkICWrVqhV69emnSniVMmDABa9euRW5urrVDISJq0oqLi+Hs7IypU6ciPDy8wjrHjx/HfffdhzvuuAOpqamYO3cupk2bhk2bNil18vPz8dFHH2HixImqc2fPno0JEyZgxIgRSElJQWpqKkaOHImJEydizpw5qroTJkzAsmXLUFRUpH1HiRoQJ1lkc8LDw+Hh4YGYmBhV+bVr17B69Wo88cQTsLOzQ3p6OsLDw2EymeDu7o5HH30UmZmZlbZrNpvx+uuvo23btnByckKPHj2wefNmVZ0///wTY8aMwd/+9jc0a9YMvXr1QkJCAgD144KxsbF49dVXcfDgQeh0Ouh0OsTGxuKJJ54o9wiH2WxGu3btqvzJ3Keffophw4apygYOHIhnnnkGs2bNQuvWreHm5oYPPvgAhYWFmDJlClq2bIl27dph9erV5drbvHkzRowYoYr7//7v/9C5c2c0a9YMw4cPx19//YWNGzeiY8eOaNGiBcaOHYsrV64obfzwww8IDg6Gi4sLWrRogT59+uDAgQPK8eTkZISGhqJZs2bw9vbGM888g4sXLyrHRQSLFi1Cx44d4eTkBB8fH8ydO1c53r17d3h5eeHzzz+vdFyIiKhqhYWFmDFjBjw8PODs7Izg4GAkJSWp6mzduhUBAQFwdnbGgAEDsH79euh0OuWHmUajEcuXL8dTTz0FHx+fCq+zfPlyeHl5YfHixejSpQsmTZqEiIgIvPvuu0qdr7/+GjqdDnfeeadStnv3brz99ttYuHAh5syZg4CAAHTq1AmzZ8/GwoULsXDhQuzevVupHxYWhtzcXGzfvl27QSKyAk6yyOY4ODggIiICsbGxMJvNSvlXX32FnJwcTJgwAWfOnMGAAQPQvXt37N69G9999x0uXbqEESNGqM4p7YMPPsA777yDhQsXIi0tDQ888ABGjhypPFJ3+fJlhIaG4sSJE/jyyy+RlpaGefPmVdjWqFGjMGvWLAQEBODMmTM4c+YMRo0ahUmTJiEuLg5nzpxR6sbHxyMzMxNjx46tsK3c3Fykp6cjKCio3LG1a9fCZDJh165dmDNnDmbMmIH7778fnTp1wp49exAREYEnn3xSdT0RwZYtW5RJFnA9CS9atAhr167Ftm3bsGfPHjz44INYtWoVNm3ahC+//BL/+c9/sHTpUgBAUVERRowYgX79+mHfvn3YtWsXZsyYAXt7ewDXn5sPCwvD8OHDsW/fPnz++efYu3ev6jNlL7/8Ml5//XXMnTsXBw8exIYNG9C2bVtV//r06YMdO3ZUOC5ERFS9l156Cf/+978RHR2N1NRU9OjRA0OGDFHywqlTpzBy5EiEh4dj3759eO655/DSSy/V+jo7d+5EWFiYqmzw4MHYs2cPrl27BgBITExEYGAgdDqdUmft2rVwcXHBs88+W67NZ555BkajEevWrVPKHB0d0atXL+YGavyEyAYdOXJEAMg333yjlN13330yZMgQERH5xz/+IXfffbfqnNzcXAEgu3btEhGRyMhI6datm3Lcy8tLXn31VdU5oaGhMmbMGBERWbFihbi4uMjZs2crjCkmJkaMRqOyX7b9Et26dZMFCxYo+4888og8+OCDlfY1NTVVAEhGRka52IKDg5V9s9ksrq6uMmzYMKXs6tWrotfrZcOGDUrZ7t27pVWrVnLt2jUlbgBy+PBhpc6sWbPEzs5O1deIiAgJDw8XEZFz584JANm+fXuFMY8dO1aeeOKJCvuRlZUleXl54uTkJMuWLau03yIizz//vPTr16/KOkREpFZyv7506ZLo9XpZtWqVcqyoqEj8/Pzk73//u4iIzJkzRzp37qw6/4033hAAcvz48XJtT5kyRUJDQ8uVd+zYsVwO3bFjhwCQP//8U0RERowYIePGjVPVGTJkiPTs2bPSvvTs2VPuvfdeVdkDDzwgjz/+eKXnEDUG/E0W2aSOHTsiNDQU0dHRAK4/xvfNN98oz3mnpKTghx9+gIuLi/Iq+S1JRSsaXbx4EX/++afqEQYA6NevH9LT0wEAqamp6NmzJ1xdXesV+6RJk5RHHXNzc7F58+Zyz6eXVvKInrOzc7ljPXv2VLZ1Oh3c3d3Ro0cPpUyv16NVq1bIzs5WyjZv3ozw8HA4ODgoZU5OTggICFD2PTw84Onpqeqrh4eH0k7r1q0xfvx4DB48GOHh4Xjvvfdw6tQppW5KSgrWrFmjGv+SsT127BjS09NRWFhY6QImJQwGg+oRRSIiqrljx47h2rVrqtxmb2+PkJAQJbcdPnwYt99+u+q8vn37WiSeK1euVJjLquPo6KjaZ26gmwEnWWSzJk6ciC+//BK5ubmIjY1F69atlUfgzGYzwsPDsXfvXtXr6NGjtV7WtvRjDVoYO3YsTp48iaSkJKxduxZubm4YPHhwpfVLJjrnz58vd0yv15eLtaKy0o9Ifvnll6pHBQGoJlw1bScmJga7du3CgAEDsGXLFgQEBOCbb74BcH38n3zySdXY79u3D0ePHq3VYhu5ublwc3OrcX0iIqoZrXObp6cnsrKyVGVZWVlwcHBQ8pirq2u5XNapUyccO3YMhYWF5dosLCzEsWPH0KlTJ1U5cwPdDDjJIpv10EMPwdnZGWvWrEF0dDTGjRunTAx69+6NgwcPwtfXF/7+/qqXyWQq11bz5s3h5eWFH3/8UVWelJSErl27AgBuu+027N+/Hzk5OTWKz9HREcXFxeXKW7dujZEjRyI6OhrR0dGIiIiAnV3l/9VuueUWNG/eXPmpY30cO3YMv/32G4YMGVLvtgDg1ltvxezZs7F9+3YMHDgQq1atAvDf8S879v7+/jAYDOjSpQucnJywbdu2Kts/cOAAevfurUmsRERNzS233AJHR0dVbisuLsbOnTuV3Na5c2fs2bNHdV7phSZqKiQkBPHx8aqy+Ph4BAUFKbn5tttuK5fLxowZg8uXL2PZsmXl2ly6dCny8/Mxbtw4VTlzA90MOMkim2UwGPDYY48hKioKx44dUz1yN2XKFPz1118YNWoUdu3ahYyMDHz33Xd46qmnkJeXV2F7L774It59912sW7cOR44cwbx585CYmIgXXngBAPDYY4/B3d0dI0aMQGJiIjIyMrBlyxZldcGy2rdvj5MnT+KXX35BTk6O6qd0kyZNwtq1a7Fv375q/8CwnZ0d7rnnnnKrQdXF5s2bMWjQoHr/0eTjx49jzpw5SE5OxsmTJ5GQkID9+/crSXv27NnYvXs3Jk+ejNTUVPz222/4z3/+g6effhoAYDKZMH36dMydOxcxMTE4duwYdu/erUqy+fn5SElJ0WxCSETU1BiNRjzzzDOYPXs2vv76axw6dAjPPPMMsrKylIUmJk+ejGPHjuGFF17Ar7/+is8//xwffvghAPVvu9LT07F3717k5OTg0qVLylMKJSZPnow//vgDM2bMwKFDh/DRRx8hNjZWyaHA9YUwDh06hHPnzillffr0wUsvvaSsJnjkyBEcPXoUb7/9NubMmYP58+er/ubkiRMn8Mcff5RbZIOo0bH2h8KIqpKSkiIA5I477ih37MiRI/Lggw9Ky5YtxdnZWTp16iRTp06VwsJCESm/MEVxcbG89tpr4uPjI3q9Xrp37y5ffPGFqs3ff/9dHnnkEWnRooUYDAbp1auXJCQkiEj5hS8KCgqU6wOQmJgY5ZjZbBY/Pz+56667atTPuLg48fT0lKKiIqUsNDRUpkyZoqrXrVs3iYyMVJV5eHjI4sWLRUSkf//+8uGHH6qOl41bROSdd94RX19fVdns2bMlMDBQREQyMzPlgQceEC8vL3F0dJS2bdvKiy++KFevXlXq//zzzzJ48GAxmUzSrFkz6d69u/zjH/9QjhcXF8uCBQukQ4cOotfrxcfHR15++WXl+KeffioBAQE1Gh8iIvqv0gsVFRQUyPTp08Xd3V0cHR2lb9++kpiYqKr/1VdfSceOHcXJyUn69esn0dHRAkAyMzOVOr6+vgKg3Ku07du3y2233SaOjo7Svn37Chc3Cg4Oln/961/lylevXi19+/YVg8GgtP3ZZ5+Vq/fmm2/K4MGD6zQuRLZEJyJirQke0c3qypUr8Pb2xuLFizFmzJganRMSEoJnn3220qXeq5OTk4M2bdrg999/h6enZ53aaEh9+vTBjBkz8Nhjj1k7FCKiJuWDDz7AvHnzcOHCBc0/uxUXF4fp06cjPT1d+bMfZWVnZ+Ouu+6Cu7s7vv76axgMBgDXP6PVsWNHrFu3rtxCVUSNjX1UVFSUtYMgulmYzWacPXsW77zzDg4ePIiVK1dWmmTK6tOnD86cOVOrhSNKO336NPz9/atd0c8WZGdn4+rVq3j66ac1T/BERKS2ZMkSAMDVq1cRFxeHOXPmYNy4cbj33ns1v5a/vz9EBF5eXmjZsmWFdYxGIx5++GHk5OTAzs4OHTp0AABkZGSgbdu2eOCBBzSPi6ih8TdZRBo6ceIEOnToAB8fH3z88cd8ppyIiKzu+eefx2effYZz587Bx8cHo0ePxrx588otnU5E2uEki4iIiIiISENcXZCIiIiIiEhDnGQRERERERFpiJMsIiIiIiIiDXGSRUREREREpCFOsoiIiIiIiDTESRYREREREZGG/h8vk/MuJmU66QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x1008 with 8 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# note: flip sense of second dimension to make the plot positive downwards\n",
    "plt_extent = [origin[0], origin[0] + extent[0], origin[1] + extent[1], origin[1]]\n",
    "\n",
    "vmin, vmax = 1.4, 1.7\n",
    "pmin, pmax = -1, +1\n",
    "dmin, dmax = 0.9, 1.1\n",
    "\n",
    "plt.figure(figsize=(12,14))\n",
    "\n",
    "# plot velocity \n",
    "plt.subplot(2, 2, 1)\n",
    "plt.imshow(np.transpose(m0.data), cmap=cm.jet, \n",
    "           vmin=vmin, vmax=vmax, extent=plt_extent)\n",
    "plt.colorbar(orientation='horizontal', label='Velocity (m/msec)')\n",
    "plt.plot(nl_rec1.coordinates.data[:, 0], nl_rec1.coordinates.data[:, 1], \\\n",
    "         'black', linestyle='-', label=\"Receiver\")\n",
    "plt.plot(src.coordinates.data[:, 0], src.coordinates.data[:, 1], \\\n",
    "         'red', linestyle='None', marker='*', markersize=15, label=\"Source\")\n",
    "plt.xlabel(\"X Coordinate (m)\")\n",
    "plt.ylabel(\"Z Coordinate (m)\")\n",
    "plt.title(\"Constant velocity\")\n",
    "\n",
    "# plot density\n",
    "plt.subplot(2, 2, 2)\n",
    "plt.imshow(np.transpose(1 / b.data), cmap=cm.jet,\n",
    "           vmin=dmin, vmax=dmax, extent=plt_extent)\n",
    "plt.colorbar(orientation='horizontal', label='Density (m^3/kg)')\n",
    "plt.plot(nl_rec1.coordinates.data[:, 0], nl_rec1.coordinates.data[:, 1], \\\n",
    "         'black', linestyle='-', label=\"Receiver\")\n",
    "plt.plot(src.coordinates.data[:, 0], src.coordinates.data[:, 1], \\\n",
    "         'red', linestyle='None', marker='*', markersize=15, label=\"Source\")\n",
    "plt.xlabel(\"X Coordinate (m)\")\n",
    "plt.ylabel(\"Z Coordinate (m)\")\n",
    "plt.title(\"Constant density\")\n",
    "\n",
    "# plot velocity perturbation\n",
    "plt.subplot(2, 2, 3)\n",
    "plt.imshow(np.transpose(dm.data), cmap=\"seismic\", \n",
    "           vmin=pmin, vmax=pmax, extent=plt_extent)\n",
    "plt.colorbar(orientation='horizontal', label='Velocity (m/msec)')\n",
    "plt.plot(nl_rec1.coordinates.data[:, 0], nl_rec1.coordinates.data[:, 1], \\\n",
    "         'black', linestyle='-', label=\"Receiver\")\n",
    "plt.plot(src.coordinates.data[:, 0], src.coordinates.data[:, 1], \\\n",
    "         'red', linestyle='None', marker='*', markersize=15, label=\"Source\")\n",
    "plt.xlabel(\"X Coordinate (m)\")\n",
    "plt.ylabel(\"Z Coordinate (m)\")\n",
    "plt.title(\"Velocity Perturbation\")\n",
    "\n",
    "# Plot the log of the generated Q profile\n",
    "q = np.log10(w / wOverQ.data)\n",
    "lmin, lmax = np.log10(qmin), np.log10(100)\n",
    "\n",
    "plt.subplot(2, 2, 4)\n",
    "plt.imshow(np.transpose(q.data), cmap=cm.jet, vmin=lmin, vmax=lmax, extent=plt_extent)\n",
    "plt.colorbar(orientation='horizontal', label='log10(Q)')\n",
    "plt.plot(nl_rec1.coordinates.data[:, 0], nl_rec1.coordinates.data[:, 1], \\\n",
    "         'black', linestyle='-', label=\"Receiver\")\n",
    "plt.plot(src.coordinates.data[:, 0], src.coordinates.data[:, 1], \\\n",
    "         'red', linestyle='None', marker='*', markersize=15, label=\"Source\")\n",
    "plt.xlabel(\"X Coordinate (m)\")\n",
    "plt.ylabel(\"Z Coordinate (m)\")\n",
    "plt.title(\"log10 of $Q$ model\")\n",
    "\n",
    "plt.tight_layout()\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation of the nonlinear forward\n",
    "\n",
    "We copy the nonlinear forward PDE described in the 1st self-adjoint notebook linked above:\n",
    "\n",
    "$$\n",
    "L_t[\\cdot] \\equiv \\frac{\\omega_c}{Q} \\overleftarrow{\\partial_t}[\\cdot] + \\partial_{tt}[\\cdot]\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\frac{b}{m^2} L_t[u] =\n",
    "    \\overleftarrow{\\partial_x}\\left(b\\ \\overrightarrow{\\partial_x}\\ u \\right) +\n",
    "    \\overleftarrow{\\partial_y}\\left(b\\ \\overrightarrow{\\partial_y}\\ u \\right) +\n",
    "    \\overleftarrow{\\partial_z}\\left(b\\ \\overrightarrow{\\partial_z}\\ u \\right) + q\n",
    "$$\n",
    "\n",
    "## Quantity that is serialized \n",
    "\n",
    "Recall that for Jacobian operators we need the scaled time derivatives of the reference wavefield for both the Born source in the Jacobian forward, and the imaging condition in the Jacobian adjoint. The quantity shown below is used in those expressions, and this is what we will serialize and compress during the nonlinear forward.\n",
    "\n",
    "$$\n",
    "\\frac{2\\ b}{m_0}\\ L_t[u_0] \n",
    "$$\n",
    "\n",
    "## The two implementations\n",
    "\n",
    "We borrow the stencil from the self-adjoint operators shown in the jupyter notebooks linked above, and make two operators. Here are details about the configuration.\n",
    "\n",
    "### Save all time steps implementation\n",
    "  * wavefield for propagation: ```u1 = TimeFunctions(..., save=None)```\n",
    "  * wavefield for serialization: ```v1 = TimeFunctions(..., save=nt)```\n",
    "  * Run the operator in a single execution from ```time_m=1``` to ```time_M=nt-1```\n",
    "\n",
    "### Time blocking implementation\n",
    "  * wavefield for propagation: ```u2 = TimeFunctions(..., save=None)```\n",
    "  * wavefield for serialization: ```v2 = TimeFunctions(..., save=Buffer(M))```\n",
    "  * Run the operator in a sequence of time blocks, each with $M$ time steps, from ```time_m=1``` to ```time_M=nt-1```\n",
    "  \n",
    "## Note on code duplication\n",
    "\n",
    "The stencils for the two operators you see below are exactly the same, the only significant difference is that we use two different ```TimeFunction```s. We could therefore reduce code duplication in two ways:\n",
    "\n",
    "1. Use the placeholder design pattern and ```stencil.subs``` to substitude the appropriate ```TimeFunction```. \n",
    "Please see the FAQ for more information [https://github.com/devitocodes/devito/wiki/FAQ#how-are-abstractions-used-in-the-seismic-examples](https://github.com/devitocodes/devito/wiki/FAQ#how-are-abstractions-used-in-the-seismic-examples)\n",
    "\n",
    "2. Write a function and use it to build the stencils.\n",
    "\n",
    "To increase the clarity of the exposition below, We do neither of these and duplicate the stencil code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Kernel` run in 0.01 s\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "# Define M: number of time steps in each time block\n",
    "M = 5\n",
    "\n",
    "# Create TimeFunctions\n",
    "u1 = TimeFunction(name=\"u1\", grid=grid, time_order=2, space_order=space_order, save=None)\n",
    "u2 = TimeFunction(name=\"u2\", grid=grid, time_order=2, space_order=space_order, save=None)\n",
    "v1 = TimeFunction(name=\"v1\", grid=grid, time_order=2, space_order=space_order, save=nt)\n",
    "v2 = TimeFunction(name=\"v2\", grid=grid, time_order=2, space_order=space_order, save=Buffer(M))\n",
    "\n",
    "# get time and space dimensions \n",
    "t,x,z = u1.dimensions\n",
    "\n",
    "# Source terms (see notebooks linked above for more detail)\n",
    "src1_term = src.inject(field=u1.forward, expr=src * t.spacing**2 * m0**2 / b)\n",
    "src2_term = src.inject(field=u2.forward, expr=src * t.spacing**2 * m0**2 / b)\n",
    "nl_rec1_term = nl_rec1.interpolate(expr=u1.forward)\n",
    "nl_rec2_term = nl_rec2.interpolate(expr=u2.forward)\n",
    "\n",
    "# The nonlinear forward time update equation\n",
    "update1 = (t.spacing**2 * m0**2 / b) * \\\n",
    "    ((b * u1.dx(x0=x+x.spacing/2)).dx(x0=x-x.spacing/2) + \\\n",
    "     (b * u1.dz(x0=z+z.spacing/2)).dz(x0=z-z.spacing/2)) + \\\n",
    "    (2 - t.spacing * wOverQ) * u1 + \\\n",
    "    (t.spacing * wOverQ - 1) * u1.backward\n",
    "\n",
    "update2 = (t.spacing**2 * m0**2 / b) * \\\n",
    "    ((b * u2.dx(x0=x+x.spacing/2)).dx(x0=x-x.spacing/2) + \\\n",
    "     (b * u2.dz(x0=z+z.spacing/2)).dz(x0=z-z.spacing/2)) + \\\n",
    "    (2 - t.spacing * wOverQ) * u2 + \\\n",
    "    (t.spacing * wOverQ - 1) * u2.backward\n",
    "\n",
    "stencil1 = Eq(u1.forward, update1)\n",
    "stencil2 = Eq(u2.forward, update2)\n",
    "\n",
    "# Equations for the Born term\n",
    "v1_term = Eq(v1, (2 * b * m0**-3) * (wOverQ * u1.dt(x0=t-t.spacing/2) + u1.dt2))\n",
    "v2_term = Eq(v2, (2 * b * m0**-3) * (wOverQ * u2.dt(x0=t-t.spacing/2) + u2.dt2))\n",
    "\n",
    "# Update spacing_map (see notebooks linked above for more detail)\n",
    "spacing_map = grid.spacing_map\n",
    "spacing_map.update({t.spacing : dt})\n",
    "\n",
    "# Build the Operators\n",
    "nl_op1 = Operator([stencil1, src1_term, nl_rec1_term, v1_term], subs=spacing_map)\n",
    "nl_op2 = Operator([stencil2, src2_term, nl_rec2_term, v2_term], subs=spacing_map)\n",
    "\n",
    "# Run operator 1 for all time samples\n",
    "u1.data[:] = 0\n",
    "v1.data[:] = 0\n",
    "nl_rec1.data[:] = 0\n",
    "nl_op1(time_m=1, time_M=nt-2)\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.145e+01\n",
      "2.669e-03\n",
      "1.381e-02\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "# Continuous integration hooks for the save all timesteps implementation\n",
    "# We ensure the norm of these computed wavefields is repeatable\n",
    "print(\"%.3e\" % norm(u1))\n",
    "print(\"%.3e\" % norm(nl_rec1))\n",
    "print(\"%.3e\" % norm(v1))\n",
    "assert np.isclose(norm(u1), 4.145e+01, atol=0, rtol=1e-3)\n",
    "assert np.isclose(norm(nl_rec1), 2.669e-03, atol=0, rtol=1e-3)\n",
    "assert np.isclose(norm(v1), 1.381e-02, atol=0, rtol=1e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the time blocking implementation over blocks of M time steps\n",
    "\n",
    "After each block of $M$ time steps, we return control to Python to extract the Born term and serialize/compress. \n",
    "\n",
    "The next cell exercises the time blocking operator over $N$ time blocks. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "# We make an array the full size for correctness testing\n",
    "v2_all = np.zeros(v1.data.shape, dtype=dtype)\n",
    "\n",
    "# Number of time blocks\n",
    "N = int((nt-1) / M) + 1\n",
    "\n",
    "# Open a binary file in append mode to save the wavefield chunks  \n",
    "filename = \"timeblocking.nonlinear.bin\"\n",
    "\n",
    "if os.path.exists(filename):\n",
    "    os.remove(filename)\n",
    "f = open(filename, \"ab\")\n",
    "\n",
    "# Arrays to save offset and length of compressed data\n",
    "file_offset = np.zeros(nt, dtype=np.int64)\n",
    "file_length = np.zeros(nt, dtype=np.int64)\n",
    "\n",
    "# The length of the data type, 4 bytes for float32\n",
    "itemsize = v2.data[0,:,:].dtype.itemsize \n",
    "\n",
    "# The length of a an uncompressed wavefield, used to compute compression ratio below\n",
    "len0 = 4.0 * np.prod(v2._data[0,:,:].shape)\n",
    "\n",
    "# Loop over time blocks\n",
    "v2_all[:] = 0\n",
    "u2.data[:] = 0\n",
    "v2.data[:] = 0\n",
    "nl_rec2.data[:] = 0\n",
    "for kN in range(0,N,1):\n",
    "    kt1 = max((kN + 0) * M, 1)\n",
    "    kt2 = min((kN + 1) * M - 1, nt-2)\n",
    "    nl_op2(time_m=kt1, time_M=kt2)\n",
    "\n",
    "    # Copy computed Born term for correctness testing\n",
    "    for kt in range(kt1,kt2+1):\n",
    "\n",
    "        # assign\n",
    "        v2_all[kt,:,:] = v2.data[(kt%M),:,:]\n",
    "    \n",
    "        # compression\n",
    "        c = blosc.compress_ptr(v2._data[(kt%M),:,:].__array_interface__['data'][0], \n",
    "                               np.prod(v2._data[(kt%M),:,:].shape), \n",
    "                               v2._data[(kt%M),:,:].dtype.itemsize, 9, True, 'zstd')\n",
    "\n",
    "        # compression ratio\n",
    "        cratio = len0 / (1.0 * len(c))\n",
    "        \n",
    "        # serialization\n",
    "        file_offset[kt] = f.tell()\n",
    "        f.write(c)\n",
    "        file_length[kt] = len(c)\n",
    "\n",
    "        # Uncomment these lines to see per time step output\n",
    "#         rms_v1 = np.linalg.norm(v1.data[kt,:,:].reshape(-1))\n",
    "#         rms_v2 = np.linalg.norm(v2_all[kt,:,:].reshape(-1))\n",
    "#         rms_12 = np.linalg.norm(v1.data[kt,:,:].reshape(-1) - v2_all[kt,:,:].reshape(-1))\n",
    "#         print(\"kt1,kt2,len,cratio,|u1|,|u2|,|v1-v2|; %3d %3d %3d %10.4f %12.6e %12.6e %12.6e\" % \n",
    "#               (kt1, kt2, kt2 - kt1 + 1, cratio, rms_v1, rms_v2, rms_12), flush=True)\n",
    "        \n",
    "# Close the binary file\n",
    "f.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.669e-03\n",
      "1.381e-02\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "# Continuous integration hooks for the time blocking implementation\n",
    "# We ensure the norm of these computed wavefields is repeatable\n",
    "# Note these are exactly the same norm values as the save all timesteps check above\n",
    "print(\"%.3e\" % norm(nl_rec1))\n",
    "print(\"%.3e\" % np.linalg.norm(v2_all))\n",
    "assert np.isclose(norm(nl_rec1), 2.669e-03, atol=0, rtol=1e-3)\n",
    "assert np.isclose(np.linalg.norm(v2_all), 1.381e-02, atol=0, rtol=1e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correctness test for nonlinear forward wavefield\n",
    "\n",
    "We now test correctness by measuring the maximum absolute difference between the wavefields computed with the```Buffer``` and save all time steps implementations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative norm of difference wavefield; +0.0000e+00\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "norm_v1 = np.linalg.norm(v1.data.reshape(-1))\n",
    "norm_v12 = np.linalg.norm(v1.data.reshape(-1) - v2_all.reshape(-1))\n",
    "\n",
    "print(\"Relative norm of difference wavefield; %+.4e\" % (norm_v12 / norm_v1))\n",
    "\n",
    "assert norm_v12 / norm_v1 < 1e-7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation of Jacobian linearized forward\n",
    "\n",
    "As before we have two implementations:\n",
    "1. operates on all time steps in a single implementation and consumes the all time steps saved version of the nonlinear forward wavefield\n",
    "1. operates in time blocks, de-serializes and de-compresses $M$ time steps at a time, and consumes the compressed and serialized time blocking version of the nonlinear forward wavefield \n",
    "\n",
    "One difference in the correctness testing for this case is that we will assign the propagated perturbed wavefields to two ```TimeFunction(..., save=nt)``` for comparison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Kernel` run in 0.01 s\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "# Create TimeFunctions\n",
    "duFwd1 = TimeFunction(name=\"duFwd1\", grid=grid, time_order=2, space_order=space_order, save=nt)\n",
    "duFwd2 = TimeFunction(name=\"duFwd2\", grid=grid, time_order=2, space_order=space_order, save=nt)\n",
    "\n",
    "ln_rec1_term = ln_rec1.interpolate(expr=duFwd1.forward)\n",
    "ln_rec2_term = ln_rec2.interpolate(expr=duFwd2.forward)\n",
    "\n",
    "# The Jacobian linearized forward time update equation\n",
    "update1 = (t.spacing**2 * m0**2 / b) * \\\n",
    "    ((b * duFwd1.dx(x0=x+x.spacing/2)).dx(x0=x-x.spacing/2) + \\\n",
    "     (b * duFwd1.dz(x0=z+z.spacing/2)).dz(x0=z-z.spacing/2) + \\\n",
    "     (dm * v1)) + (2 - t.spacing * wOverQ) * duFwd1 + \\\n",
    "    (t.spacing * wOverQ - 1) * duFwd1.backward\n",
    "\n",
    "update2 = (t.spacing**2 * m0**2 / b) * \\\n",
    "    ((b * duFwd2.dx(x0=x+x.spacing/2)).dx(x0=x-x.spacing/2) + \\\n",
    "     (b * duFwd2.dz(x0=z+z.spacing/2)).dz(x0=z-z.spacing/2) + \\\n",
    "     (dm * v2)) + (2 - t.spacing * wOverQ) * duFwd2 + \\\n",
    "    (t.spacing * wOverQ - 1) * duFwd2.backward\n",
    "\n",
    "stencil1 = Eq(duFwd1.forward, update1)\n",
    "stencil2 = Eq(duFwd2.forward, update2)\n",
    "\n",
    "# Build the Operators\n",
    "lf_op1 = Operator([stencil1, ln_rec1_term], subs=spacing_map)\n",
    "lf_op2 = Operator([stencil2, ln_rec2_term], subs=spacing_map)\n",
    "\n",
    "# Run operator 1 for all time samples\n",
    "duFwd1.data[:] = 0\n",
    "ln_rec1.data[:] = 0\n",
    "lf_op1(time_m=1, time_M=nt-2)\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.438e+00\n",
      "2.681e-02\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "# Continuous integration hooks for the save all timesteps implementation\n",
    "# We ensure the norm of these computed wavefields is repeatable\n",
    "print(\"%.3e\" % norm(duFwd1))\n",
    "print(\"%.3e\" % norm(ln_rec1))\n",
    "assert np.isclose(norm(duFwd1), 6.438e+00, atol=0, rtol=1e-3)\n",
    "assert np.isclose(norm(ln_rec1), 2.681e-02, atol=0, rtol=1e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the time blocking implementation over blocks of M time steps\n",
    "\n",
    "Before each block of $M$ time steps, we de-serialize and de-compress. \n",
    "\n",
    "### TODO\n",
    "* In the linearized op below figure out how to use a pre-allocated array to hold the compressed bytes, instead of returning a new bytearray each time we read. Note this will not be a problem when not using Python, and there surely must be a way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "# Open the binary file in read only mode\n",
    "f = open(filename, \"rb\")\n",
    "\n",
    "# Temporay nd array for decompression\n",
    "d = copy.copy(v2._data[0,:,:])\n",
    "\n",
    "# Array to hold compression ratio\n",
    "cratio = np.zeros(nt, dtype=dtype)\n",
    "\n",
    "# Loop over time blocks\n",
    "duFwd2.data[:] = 0\n",
    "ln_rec2.data[:] = 0\n",
    "for kN in range(0,N,1):\n",
    "    kt1 = max((kN + 0) * M, 1)\n",
    "    kt2 = min((kN + 1) * M - 1, nt-2)\n",
    "    \n",
    "    # 1. Seek to file_offset[kt] \n",
    "    # 2. Read file_length[kt1] bytes from file\n",
    "    # 3. Decompress wavefield and assign to v2 Buffer\n",
    "    for kt in range(kt1,kt2+1):\n",
    "        f.seek(file_offset[kt], 0)\n",
    "        c = f.read(file_length[kt])\n",
    "        blosc.decompress_ptr(c, v2._data[(kt%M),:,:].__array_interface__['data'][0])\n",
    "        cratio[kt] = len0 / (1.0 * len(c))\n",
    "\n",
    "    # Run the operator for this time block\n",
    "    lf_op2(time_m=kt1, time_M=kt2)\n",
    "    \n",
    "    # Uncomment these lines to see per time step outputs\n",
    "#     for kt in range(kt1,kt2+1):\n",
    "#         rms_du1 = np.linalg.norm(duFwd1.data[kt,:,:].reshape(-1))\n",
    "#         rms_du2 = np.linalg.norm(duFwd2.data[kt,:,:].reshape(-1))\n",
    "#         rms_d12 = np.linalg.norm(duFwd1.data[kt,:,:].reshape(-1) - duFwd2.data[kt,:,:].reshape(-1))\n",
    "#         print(\"kt1,kt2,len,cratio,|du1|,|du2|,|du1-du2|; %3d %3d %3d %10.4f %12.6e %12.6e %12.6e\" % \n",
    "#               (kt1, kt2, kt2 - kt1 + 1, cratio[kt], rms_du1, rms_du2, rms_d12), flush=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.438e+00\n",
      "2.681e-02\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "# Continuous integration hooks for the save all timesteps implementation\n",
    "# We ensure the norm of these computed wavefields is repeatable\n",
    "# Note these are exactly the same norm values as the save all timesteps check above\n",
    "print(\"%.3e\" % norm(duFwd2))\n",
    "print(\"%.3e\" % norm(ln_rec2))\n",
    "assert np.isclose(norm(duFwd2), 6.438e+00, atol=0, rtol=1e-3)\n",
    "assert np.isclose(norm(ln_rec2), 2.681e-02, atol=0, rtol=1e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correctness test for Jacobian forward wavefield\n",
    "\n",
    "We now test correctness by measuring the maximum absolute difference between the wavefields computed with the```Buffer``` and save all time steps implementations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative norm of difference wavefield; +0.0000e+00\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "norm_du1 = np.linalg.norm(duFwd1.data.reshape(-1))\n",
    "norm_du12 = np.linalg.norm(duFwd1.data.reshape(-1) - duFwd2.data.reshape(-1))\n",
    "\n",
    "print(\"Relative norm of difference wavefield; %+.4e\" % (norm_du12 / norm_du1))\n",
    "\n",
    "assert norm_du12 / norm_du1 < 1e-7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation of Jacobian linearized adjoint\n",
    "\n",
    "Again we have two implementations:\n",
    "1. operates on all time steps in a single implementation and consumes the all time steps saved version of the nonlinear forward wavefield\n",
    "1. operates in time blocks, de-serializes and de-compresses $M$ time steps at a time, and consumes the compressed and serialized time blocking version of the nonlinear forward wavefield \n",
    "\n",
    "For correctness testing here we will compare the final gradients computed via these two implementations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Kernel` run in 0.01 s\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "# Create TimeFunctions for adjoint wavefields\n",
    "duAdj1 = TimeFunction(name=\"duAdj1\", grid=grid, time_order=2, space_order=space_order, save=nt)\n",
    "duAdj2 = TimeFunction(name=\"duAdj2\", grid=grid, time_order=2, space_order=space_order, save=nt)\n",
    "\n",
    "# Create Functions to hold the computed gradients\n",
    "dm1 = Function(name='dm1', grid=grid, space_order=space_order)\n",
    "dm2 = Function(name='dm2', grid=grid, space_order=space_order)\n",
    "\n",
    "# The Jacobian linearized adjoint time update equation\n",
    "update1 = (t.spacing**2 * m0**2 / b) * \\\n",
    "    ((b * duAdj1.dx(x0=x+x.spacing/2)).dx(x0=x-x.spacing/2) +\n",
    "     (b * duAdj1.dz(x0=z+z.spacing/2)).dz(x0=z-z.spacing/2)) +\\\n",
    "    (2 - t.spacing * wOverQ) * duAdj1 + \\\n",
    "    (t.spacing * wOverQ - 1) * duAdj1.forward\n",
    "\n",
    "update2 = (t.spacing**2 * m0**2 / b) * \\\n",
    "    ((b * duAdj2.dx(x0=x+x.spacing/2)).dx(x0=x-x.spacing/2) +\n",
    "     (b * duAdj2.dz(x0=z+z.spacing/2)).dz(x0=z-z.spacing/2)) +\\\n",
    "    (2 - t.spacing * wOverQ) * duAdj2 + \\\n",
    "    (t.spacing * wOverQ - 1) * duAdj2.forward\n",
    "\n",
    "stencil1 = Eq(duAdj1.backward, update1)\n",
    "stencil2 = Eq(duAdj2.backward, update2)\n",
    "\n",
    "# Equations to sum the zero lag correlations\n",
    "dm1_update = Eq(dm1, dm1 + duAdj1 * v1)\n",
    "dm2_update = Eq(dm2, dm2 + duAdj2 * v2)\n",
    "\n",
    "# We will inject the Jacobian linearized forward receiver data, time reversed\n",
    "la_rec1_term = ln_rec1.inject(field=duAdj1.backward, expr=ln_rec1 * t.spacing**2 * m0**2 / b)\n",
    "la_rec2_term = ln_rec2.inject(field=duAdj2.backward, expr=ln_rec2 * t.spacing**2 * m0**2 / b)\n",
    "\n",
    "# Build the Operators\n",
    "la_op1 = Operator([dm1_update, stencil1, la_rec1_term], subs=spacing_map)\n",
    "la_op2 = Operator([dm2_update, stencil2, la_rec2_term], subs=spacing_map)\n",
    "\n",
    "# Run operator 1 for all time samples\n",
    "duAdj1.data[:] = 0\n",
    "dm1.data[:] = 0\n",
    "la_op1(time_m=1, time_M=nt-2)\n",
    "\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.626e+01\n",
      "1.426e-04\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "# Continuous integration hooks for the save all timesteps implementation\n",
    "# We ensure the norm of these computed wavefields is repeatable\n",
    "print(\"%.3e\" % norm(duAdj1))\n",
    "print(\"%.3e\" % norm(dm1))\n",
    "assert np.isclose(norm(duAdj1), 4.626e+01, atol=0, rtol=1e-3)\n",
    "assert np.isclose(norm(dm1), 1.426e-04, atol=0, rtol=1e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the time blocking implementation over blocks of M time steps\n",
    "\n",
    "Before each block of $M$ time steps, we de-serialize and de-compress. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n",
      "Operator `Kernel` run in 0.01 s\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "# Open the binary file in read only mode\n",
    "f = open(filename, \"rb\")\n",
    "\n",
    "# Temporay nd array for decompression\n",
    "d = copy.copy(v2._data[0,:,:])\n",
    "\n",
    "# Array to hold compression ratio\n",
    "cratio = np.zeros(nt, dtype=dtype)\n",
    "\n",
    "# Loop over time blocks\n",
    "duAdj2.data[:] = 0\n",
    "dm2.data[:] = 0\n",
    "for kN in range(N-1,-1,-1):\n",
    "    kt1 = max((kN + 0) * M, 1)\n",
    "    kt2 = min((kN + 1) * M - 1, nt-2)\n",
    "        \n",
    "    # 1. Seek to file_offset[kt] \n",
    "    # 2. Read file_length[kt1] bytes from file\n",
    "    # 3. Decompress wavefield and assign to v2 Buffer\n",
    "    for kt in range(kt1,kt2+1,+1):\n",
    "        f.seek(file_offset[kt], 0)\n",
    "        c = f.read(file_length[kt])\n",
    "        blosc.decompress_ptr(c, v2._data[(kt%M),:,:].__array_interface__['data'][0])\n",
    "        cratio[kt] = len0 / (1.0 * len(c))\n",
    "\n",
    "    # Run the operator for this time block\n",
    "    la_op2(time_m=kt1, time_M=kt2)\n",
    "    \n",
    "    # Uncomment these lines to see per time step outputs\n",
    "#     for kt in range(kt2,kt1-1,-1):\n",
    "#         rms_du1 = np.linalg.norm(duAdj1.data[kt,:,:].reshape(-1))\n",
    "#         rms_du2 = np.linalg.norm(duAdj2.data[kt,:,:].reshape(-1))\n",
    "#         rms_d12 = np.linalg.norm(duAdj1.data[kt,:,:].reshape(-1) - duAdj2.data[kt,:,:].reshape(-1))\n",
    "#         print(\"kt2,kt1,kt,cratio,|du1|,|du2|,|du1-du2|; %3d %3d %3d %10.4f %12.6e %12.6e %12.6e\" % \n",
    "#               (kt2, kt1, kt, cratio[kt], rms_du1, rms_du2, rms_d12), flush=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.626e+01\n",
      "1.426e-04\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "# Continuous integration hooks for the save all timesteps implementation\n",
    "# We ensure the norm of these computed wavefields is repeatable\n",
    "# Note these are exactly the same norm values as the save all timesteps check above\n",
    "print(\"%.3e\" % norm(duAdj2))\n",
    "print(\"%.3e\" % norm(dm2))\n",
    "assert np.isclose(norm(duAdj2), 4.626e+01, atol=0, rtol=1e-3)\n",
    "assert np.isclose(norm(dm2), 1.426e-04, atol=0, rtol=1e-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Correctness test for Jacobian adjoint wavefield and computed gradient\n",
    "\n",
    "We now test correctness by measuring the maximum absolute difference between the gradients computed with the```Buffer``` and save all time steps implementations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative norm of difference wavefield,gradient; +0.0000e+00 +0.0000e+00\n"
     ]
    }
   ],
   "source": [
    "# NBVAL_IGNORE_OUTPUT\n",
    "\n",
    "norm_du1 = np.linalg.norm(duAdj1.data.reshape(-1))\n",
    "norm_du12 = np.linalg.norm(duAdj1.data.reshape(-1) - duAdj2.data.reshape(-1))\n",
    "\n",
    "norm_dm1 = np.linalg.norm(dm1.data.reshape(-1))\n",
    "norm_dm12 = np.linalg.norm(dm1.data.reshape(-1) - dm2.data.reshape(-1))\n",
    "\n",
    "print(\"Relative norm of difference wavefield,gradient; %+.4e %+.4e\" % \n",
    "      (norm_du12 / norm_du1, norm_dm12 /norm_dm1))\n",
    "\n",
    "assert norm_du12 / norm_du1 < 1e-7\n",
    "assert norm_dm12 / norm_dm1 < 1e-7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Delete the file used for serialization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.remove(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}