tutorials/2D LS-RTM with devito and Automatic Differentiation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "19da947c",
   "metadata": {},
   "source": [
    "# Least-Squares Reverse Time Migration through Devito and Automatic Differentiation\n",
    "\n",
    "@Author: Francesco Picetti - picettifrancesco@gmail.com\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "329c7dfc",
   "metadata": {},
   "source": [
    "In this tutorial, we put together torch's Automatic Differentiation abilities and the OccamyPy physical operators (well, actually Devito's) for RTM. Basically, we can write a standard pytorch optimization loop with the modeling provided by OccamyPy; the tricky part is to insert a linear operator (i.e., Born) in the AD graph."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03115cc3",
   "metadata": {},
   "source": [
    "To use this notebook, we need devito installed. On your occamypy-ready env, run\n",
    "\n",
    "```pip install --user git+https://github.com/devitocodes/devito.git```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "2c4d250f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING! DATAPATH not found. The folder /tmp will be used to write binary files\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/nas/home/fpicetti/miniconda3/envs/occd/lib/python3.10/site-packages/dask_jobqueue/core.py:20: FutureWarning: tmpfile is deprecated and will be removed in a future release. Please use dask.utils.tmpfile instead.\n",
      "  from distributed.utils import tmpfile\n"
     ]
    }
   ],
   "source": [
    "import torch\n",
    "import occamypy as o\n",
    "import born_devito as b\n",
    "\n",
    "# Plotting\n",
    "from matplotlib import rcParams\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "import matplotlib.pyplot as plt\n",
    "rcParams.update({\n",
    "    'image.cmap'     : 'gray',\n",
    "    'image.aspect'   : 'auto',\n",
    "    'image.interpolation': None,\n",
    "    'axes.grid'      : False,\n",
    "    'figure.figsize' : (10, 6),\n",
    "    'savefig.dpi'    : 300,\n",
    "    'axes.labelsize' : 14,\n",
    "    'axes.titlesize' : 16,\n",
    "    'font.size'      : 14,\n",
    "    'legend.fontsize': 14,\n",
    "    'xtick.labelsize': 14,\n",
    "    'ytick.labelsize': 14,\n",
    "    'text.usetex'    : True,\n",
    "    'font.family'    : 'serif',\n",
    "    'font.serif'     : 'Latin Modern Roman',\n",
    "})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "206f7020",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5b44c66e",
   "metadata": {},
   "outputs": [],
   "source": [
    "args = dict(\n",
    "    filter_sigma=(1, 1),\n",
    "    spacing=(10., 10.),  # meters\n",
    "    shape=(101, 101),    # samples\n",
    "    nbl=20,              # samples\n",
    "    nreceivers=51,\n",
    "    src_x=[500],         # meters\n",
    "    src_depth=20.,       # meters\n",
    "    rec_depth=30.,       # meters\n",
    "    t0=0.,\n",
    "    tn=2000.,            # Simulation lasts 2 second (in ms)\n",
    "    f0=0.010,            # Source peak frequency is 10Hz (in kHz)\n",
    "    space_order=5,\n",
    "    kernel=\"OT2\",\n",
    "    src_type='Ricker',\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c26c4895",
   "metadata": {},
   "outputs": [],
   "source": [
    "def unpad(model, nbl: int = args[\"nbl\"]):\n",
    "    return model[nbl:-nbl, nbl:-nbl]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaed305a",
   "metadata": {},
   "source": [
    "create hard model and migration model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "308e7fdb",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "model_true, model_smooth, water = b.create_models(args)\n",
    "model = o.VectorNumpy(model_smooth.vp.data.__array__())\n",
    "model.ax_info = [\n",
    "    o.AxInfo(model_smooth.vp.shape[0], model_smooth.origin[0] - model_smooth.nbl * model_smooth.spacing[0],\n",
    "             model_smooth.spacing[0], \"x [m]\"),\n",
    "    o.AxInfo(model_smooth.vp.shape[1], model_smooth.origin[1] - model_smooth.nbl * model_smooth.spacing[1],\n",
    "             model_smooth.spacing[1], \"z [m]\")]\n",
    "model_extent = [0., 1000., 1000., 0.]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c38f42a",
   "metadata": {},
   "source": [
    "Define acquisition geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "401e2d94",
   "metadata": {},
   "outputs": [],
   "source": [
    "rec = b.build_rec_coordinates(model_true, args)\n",
    "src = [b.build_src_coordinates(x, args[\"src_depth\"]) for x in args[\"src_x\"]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eab9377d",
   "metadata": {},
   "source": [
    "### create the common shot gather (CSG) data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8538b9e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "d = o.VectorTorch(b.propagate_shots(model=model_true, src_pos=src, rec_pos=rec, param=args))\n",
    "assert d[:].requires_grad is False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7803a2d7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsAAAAGmCAYAAAByajI1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABvHUlEQVR4nO29a7BlV1X3/Z9JCLnnpMMlFwLJaYjKVU8a8YNg+dj9vPLJ8rE7vFZZWlRJB9FS8dItliJiaTyNypenhG60EMtLJWkp7+jbjZZoFUKSlkBShEufEEggIcnJyY3cSOb74ayT3r179zn7P84eY4891/9X1ZWcvedaa+41b/855phjllorhBBCCCGE6AunTDsDQgghhBBCRCIBLIQQQggheoUEsBBCCCGE6BUSwEIIIYQQoldIAAshhBBCiF4hASyEEEIIIXqFBLAQYiqUUuZLKYullAdLKUdHfL/23aFSysKEnrnQ3e+E580apZT93fvZPu28jENX3jeXUmop5YZp50cI0W8kgIUQU6HWulRr3QvgWgBbSimLQ9/vBXCg1rqj1npkQs88AmBx+PNObM+UKKu1XgNgib2ulLLbITsbUmtdArBrGs8WQohhJICFENNmBavCaE8pZX7oOw9L7fKIzw4BuM7hWRnZMcVnj3r3QggRzmnTzoAQQtRaD5dSDgO4AcBV03h+9DOnQWf9HZ5kCCFE75AAFkJkYReAB0spu2utB06WqJSyB8eW/udrrfu6z7fjmHvD27Aq9OYBXNi5U5zsfgvddfO11q3svQbyMw9gqdZ6sPt8DsDu7rsdAPavuXIM3H8JwP7ue6yXz4FnreAkltQxnrkDwHx3Hwy8u5Ne132/G8DetfczLt11u7r7rjDXCiGEK7VW/dM//dO/qf0DsHvg//cAeBDA3PB33d83ANg+8Pc8gEMDf2/HqtvEYJqjABYG/l4AcHTovqM+G+deNwDYOfD3obXvcUxUD147N/D3TgA3d89ZALC4wXtaHHpXcwDqUP42euZ2ADef5N7rXTcPYA9Zrotd/gbvu7v77IZp1zv90z/96/c/+QALIdJQVy2SywA+NPxdZ6ndXgfcFerqxqotA5EQlrEquAZdGtassyzr3qvzV95ZO4tvxw0Arun+fx6rgnPw2sG/V7Aqlg/XWo/U9a3Uc1gVoM9ZxmutKwCGNwdu9MyTse51dXXD4r4x7nNcfgEc6cpojV64mggh8iMXCCFENnYBuHlE6LNtGB31YG3Z/vDA34OsANhizMt699oOYGUoDNlWdAK51roLeE4MznfXDedj3CgO2zGGC8GYz5zYdeuwrfsvHaVCCCEikAVYCJGKuup7egCr1tRB5uJzsy5zWPX5PTzwb2+tdQfwXMzhGwBcjVUhOEoMrhDP2zCCwpjPHEw/b7lOCCFmHQlgIURG9mLVAjnoFnAYo10Z5gHcGJGpIY5gRH5KKXOdJfXjAK6ttR7o3BXmuu8t7hgjnzX8XMMzFxzyCgA3df+dM14vhBCuSAALIabNCZEFOhH2Ngwsw3eW4cODLgdrbhJDfrjDzE0qo4P36nyDbyql7BxKczVWxepcPf4Aj7XfQp9q1/nRHhg8xKITrgsDeRrnmcM+zEfGua47xW0Pkd8VrFrxtw+J6OF3JYQQU0ECWAgxFTpRdQNWD8DY3wm65+hE7eGhz3YB2FFK2d2JwbfUWq/q7rcA4F0YCPPV/XcbgGtKKTuH0iyuc92G9+ryswPA67v87Cyl7OysqEcA7OtOmNveifZd6MKddX/vXbv/OJbWunry21z3nO1dXo4AeFcpZftGz+zusSakF7G6oXBpnOuwKoSvAUGX370A9nflu3/g6+1DfwshRCil1jrtPAghhBBCCBGGLMBCCCGEEKJXSAALIYQQQoheIQEshBBCCCF6hQSwEEIIIYToFRLAQgghhBCiV0gACyGEEEKIXiEBLIQQQggheoUEsBBCCCGE6BUSwEIIIYQQoldIAAshhBBCiF4hASyEEEIIIXqFBLAQQgghhOgVEsBCCCGEEKJXSAALIYQQQoheIQEshBBCCCF6hQSwEEIIIYToFRLAQgghhBCiV0gACyGEEEKIXiEBLIQQQggheoUEsBBCCCGE6BUSwEIIIYQQoldIAAshhBBCiF4hASyEEEIIIXrFaVEPKqXMA9gJ4AiABQAHaq0rbFrmPkIIIYQQQgxTaq0xDyrlUK11R/f/8wD21lqvYdMy9xFCCCGEEGKYEBeITqg+R611CcDVbFrmPkIIIYQQQowiygd4AcDK8IfDgnaMtMx9hBBCCCGEOIEoAbxlxGfLAObItMx9hBBCCCGEOIGwTXCRlFJ2A9gNAM973vOuesELXjDlHAkhhBBCCAvf+MY37q+1vnCS94wSwKOstFswwp1hg7Rj3afWegDAAQC45JJL6jXXaI+cECIX7Abkb3/72/QzTjuN6+JLKfQzhBDCm/e85z13TvqeUQL4CEa4L3Sb2MZO23XO495HCCFMWKLjsOKRfca3vvUtKj0AnHvuuVR6799geYYQQngQIoAHxCuA5zatXT/093KtdWW9tBvdRwghRsEKNYtIY5/x9NNPU+lXVlao9ABwxhlnUOlPP/10Kn3Ee5JgFkJ4EOkDvKuUsgfdARZDsXsXARxC57awQdr1vhNCiKngbdG97bbbqPQAcN5551Hpn/e851HpJU6FELNKmADu3BT2dX8eHvpuF5H2pN8JIfqBtxXRsrT/zDPPUOkffPBBKv2nP/1pKj0AfOd3fieV/pxzzqHSWwSwd1lIlAshxqHJKBBCCBHNE088QaW/5ZZbnHJif8YLX8htsmYFsxBCZEECWAgxdbJZdJ999lkqPQA8+uijVPqMAvj7vu/7qPRnnXUWlR4ATjmFCz8vi7EQwgMJYCGEGMIScuzee+91yEks7G/YsmXU2UTrw260E0IIDySAhRATJWNoLNaiy1pzAeCTn/wkfU022N9w2WWX0c+Ym5uj0p966qlUeoVyE0KMgwSwEGLm8D5EYnl5mUoP9NMCbHlPrN+wt8uEEKKfSAALIdYlo08lmyd2g5ol5FgfsbynF73oRVR679BsimUsRD+RABZCTB3vTWqPPPIIlf4zn/kMlb6vWN7Ttm3bqPRnn302lT5CAAshZh8JYCF6RkbrlbdF9+jRo1R64QdbFqzPMBuZQrGMhegnEsBCiIli2VTEXvPYY49R6W+99VYqfUZa2azFlgV7mMeZZ55JpW/lvQohOCSAhZhhWhm8n3rqKSr9XXfdRaW/7777qPQRvPjFL073jIiNfGxZsGXNHv98xhlnUOktKDKFEPmQABZCrEvEIRKPP/44lf5LX/oS/Yxs3HPPPemekVFEsWV9xRVXUOktcYkVmUKI2UcCWIhEtOAraDlEghVqt99+O/0MMZuwZf3a176WSs9usgP8D/NQZAoh/JEAFkKsi/cGNQD46le/Sl+TDcuydTZaEFFsXbr00kvpZ3iHZhNC+CMBLIQjGQUFm6enn36aSv/AAw9Q6QHgv//7v+lrhBgFW5euvPJK+hms3zBrMVZkCiH8kQAWomd4C+C7776bSp+RFqy5EbSwWctSX7Md5iGE4JEAFoIgm9XEIkCeeeYZKv2DDz5Ipf/0pz9NpRdimljq6+WXX06lf/7zn0+lj7AAy2Is+o4EsBA9g7XofvOb36TSs6euRfCJT3xi2lnYNO94xzvoa/74j//YISebgy2LN73pTU45WcVSX9k2ceGFF1LpTztNQ7MQ3qiVid6ScZk7Y8ixpaUl+hnZeOMb3zjtLGyaD3zgA/Q1GQVwC2XBtglvizGQLzRbC+4xom0kgIUgyNZBs9ZcgA85dtttt9HP8Ob973//tLOwaV72spele8add97plBM7bFm/853vdMrJMdg28cpXvpJKbwnNZhHNDBGHeQgRiQSwEI54W3QtIce+9rWv0ddk4xd/8RennYVNExH6rYXwcmxZRwhgFrbNRYRmy2YxFiIaCWDRDC1s0mAPkbj//vvpZ2QMOfb2t7992lnYNL/0S7807SxsGvY3/NEf/ZFTTuywdemDH/ygU06Owba5l7/85fQzzjzzTCq9DvMQfae0vkxxySWX1GuuuWba2RABZOw8WYvuY489RqU/cuQIlR4A/v3f/52+xpsW+qE+DsYqNx9+8Ad/kL5mYWGBSs+6WbAWYwsZ+3CRg/e85z0311q3TfKesgCLtGTsDL1j6C4vL1Ppb775Zip9BG9+85unnYVN85GPfGTaWZgJLO/pp37qpxxyYsdSXz/2sY855OQYlnbNbrRjLcA6zEO0hgSwEARsh8u6NLDhlR5++GEqfQT//M//PO0sbJr3vOc9087CTGB5T9kEsKW+egspS7tm+44XvOAFVHod5iFaQwJYhJAxJE7EIRIPPfQQlf6zn/0slT6C173uddPOwqZh3T6+8pWv+GSkMSzviS0LizuAN2ybuOWWW5xycgy272A32mW0AEdEppCQbxcJYCEIWIvufffdR6XPGKHh3/7t36adhU3zu7/7u9POguhgyyKjAGbbBHsQhgW272D7pgsuuIBKr8M8RHZUQ4WJjL5XbJ5Yay7AL01+/vOfp5/hDWv52bJli1NO7Nx6661U+o9//ONOObHDCoqzzjqLfsa3vvUtKj177LUFtizYsn71q19NpbfAtgm2zd19991Uegts33TRRRdR6VmXCQA49dRTqfTeFmMg51gnJoMEsOgtrDUXAFZWVqj0t99+O/0Mbz760Y9OOwub5vd+7/emnYUTOOOMM6j0559/PpXectABK0LYUwItcahZ2LL+q7/6K6ec2GHb3Bve8AannByD7Zu+53u+h0rP1m+AF8BCbAYJYAEg5yzX+xAJNuQYAHzpS1+ir/HEYhX83u/9XoecbI677rqLSv/Xf/3XTjmxw24qOuecc6j0liVlVgCzv4EtNwtsWe/bt49K/5KXvIRKb4FtcxHWfha272PrEpDzMA9FpmgXCWDRDKxFl92gBgBf+MIX6Gs8+fCHPzztLEyEbD66rDgFeNHCCtoIn0r2N1je06OPPkpfw8DWpQ984ANOObFjaddvectbHHJyDLbvY49/Bvj65H2Yh2gbHYTRIBFlGrH7lrXoPvLII1R6yyESn/jEJ+hrPMnYfi0hnFifSov/NoPFKsj+BnbwtghgdlL41FNPUenZuNWAv9WYXUa3/IbzzjuPvsabbJbEN73pTfQ17GEe5557LpXecphHxFjHkq2sI5jpgzBKKfMAdgI4AmABwIFa68pJ0i4A2N79+XoAb1tLW0pZBHAUwPUArgawXGs96Jr5HpCxQXkL4C9/+ctU+ggyHi3L8ju/8zv0Nd6ClrVssv68gL9F1zJ4s89g25zlPbFlwS7ts3XJUl/f97730dd4w/Yd3kd3W/rXV7ziFVT6jKfZZRTMYjSRLhD7a607AKCUsgRgEcAJptlSyhyAbbXWfd3fOwF8HMBVg/fqrr+21nrAOd9Tp4UGYvkNrDXq/vvvp9J//etfp9JH8M53vnPaWdg011133bSzcAKspciy4YwVm9473i3XsL/B8p7YsvD2bbXU14wCmO07vAWwpX9l+3B21cWyyS6jcUh+xpMhRAB31t/nqLUulVKuxggBDGAbgL0A1oTtYQA3lFLmOivwjbVWleY6ZNyg1teQY3v27Jl2FjbNu9/9bip9RCxjVnixvoUWyya7gYcVmxEhnFgs74ktCzbyypNPPkmlt9RXtk28973vpZ/hDds3sZsLLXiHZrP4DHuHWsvYrvtClAV4AcDK8IellPla69LgZ7XWw6WUXQMfzXefrwxdu1Br5Z04E9DHymsJOcb65n3xi1+kn+HN4uLitLOwaf7yL/9y2lk4AdaKyAo1SwxTbxeIjAOl5T2xZcGWNSuALbBtIqMAZvumCAHM9uHsiXxzc3NUeqCNAz1kMR5NVMmOWqdYBjA3KvGQsH0LgMGWN19K2Q7gps4f+LphIVxK2Q1gN2CLRZiNFiy6lp3f2UKOAcBb3/rWaWdh07BxVZeWljZOtEnYQYZt12eeeSaV3rK0z1qXIlwg2HYd4SPJTobZsmYtxpbJOdsm2Db367/+61T6CNi+LyJCDTtGvOhFL6KfwU7yMro29dHoNg6ppzadP/DCmu8wAKz5Bnff7wdwCMDWwes6v+ADwGoUiJDM9gx20GAHJSDnJrWMIZNY/vzP/3zaWTgBb0HLWh0tS6UtWIAjJtvs5Jkta7YuPfDAA1R6C2ybyyiA2b4vQgCzY8SrXvUq+hnsCoQO85gdogTwKGvvFoxwixhiEcCgOwQGfIHXfInnR10YSbZBJiLkGLsx5Y477qDSA/7xQn/kR36EvsZiGfTkD//wD+lrvGMZWwYA9lhg1o80wgKc0VLE9gWsKH/66aep9BbYyTZblyyTc1bEs23O0q5/+Zd/mb6GgW0Tlv717/7u76j07BhhGYcuvPBCKr33YR6Av59xNk3jRZQAPoIRbhDD/r+DlFL2ANhba11ZE72d68Mijo8IISYAO8iwG9QiltFZLINMNv7sz/5s2lk4AYufnbegZS26Ft9W9hkR7gYsEdYrdnBly5qtS5b66m01trRrbwHMYulfWQHMYhmHrrzySio9G5pNh3lMjxAB3Flqn/u7s9peP/T38kCs350ADg5sfNsO4CCAmwBcO3Ddzu7zSeaVviabRZe15gLAE088QaW/8847qfT33HMPld7CD/zAD1Dpt27dunGiYP7v//2/VPpbb73VKSfHYIWa5QhU9mABNpYsK6IsAjjC8uMN23dYfoO3yGats5b6+uCDD1Lp2fdqadds3/FzP/dz9DMYLP0r24f/x3/8B5XeMg6xYx27AmHZZOd9BHSEj3EGq3GkD/Cuzqp7BKt+vYMh0Bax6st7oBPDNwDHvaAldIK4lLLU3WcFwNZa63EuEsImgNlDJL761a/Sz/DGEtA+G3/6p3867SycABtrk/WZA/x9ellxGuEDnDFcUsalTzZP7GqWpb6ybYKNb2uB7Tu8BbAFtg+3nDbHwo51rPBnJ/NAzsnzLBImgDt3h7UNbIeHvts1lO6kvWoX8YEKf5ZtB6T3oMRacwHg7rvvptIfPXqUfgbLd3/3d1Pp3/jGN/pkZBP8yZ/8CZX+M5/5jE9GBmAtcOzOaUvklWwWXYsFuI+bXyy/mR282WewBgBLFAi2TbAWY0vcdLbvYPumn/7pn6bSW2D7cHaMsPSv7FjHukxY+kvWzYIl48TWg9RRIKZBBrP8MGxn+Nhjj9HPiDi8gCXjTmiWD33oQ9POwgmwS77e7glAPotuxNJ+RiuOZfWIxbuPZeuSRQCzbYJtc/feey+V3gLbN0UIYBZ2jLj66qudcnIMdix92cteRj+DrePesYwzxigfBwngKcAWPHsksMXPyduXlJ0VA8CuXfm8W/7iL/6CSv/pT3/aKSerWDo29vQk7wgNAC+aWUHLvieLOM0oaFlYEW+xAHv7GWeMZfz4449T6S2b7Ng8sX0T2/f9xE/8BJXeAjtGWMYh9nAOdix9+ctfTqUH/EOzZTQEetALATzrm9TYkGOWM9i9+YVf+IVpZ2Ei7N+/f9pZOA5LYHd2yY0VpxEhxFhBGzEAeLtAzKqVZZhs/bGlvrJtgm1zlnbt3e+zfV+EAGaxjEM/+7M/65CTY1jK7dJLL6XSs3U8Yn9CBpHdCwHsjbdF97777qPS33bbbVR6Cy95yUuo9O94xzuccmLnb/7mb+hr/uu//sshJ8dgLZsXX3wx/Qxviy67PAf4n6IW4Z6QoUOPxvKeWANAxGDM4r0xz9Ku2Y127DjE9n2W/vXHfuzH6GsYLOPQtddeu3GiAe666y4qvWW8vvzyy6n0rM8w24fPat8nATyExWLivUmNdWmw+ACz7N692/0Z3nzwgx+cdhZO4MUvfjGV3hLD1NuiG7GBzNulYVYtGtFEhD/yfq+W+sqKeLbNWdo123d47/uw9K/eAtgCO9a9+93vptJbxmtWE7Bub2yfrzBoPYad3bcQcuw3fuM3pp2FE2CDnB8+fHjjRMGwpw5ZQjh5bziz+CV7C9oW/HMzknGi4L3hB+AFMNvmLO2a7Tu8BbClf2X78Pl5/4Ng2bGOFcAWvEOz9eUwj14IYGZ2EnGIBBty7Ctf+QqV3sI73/lOKn2G2dswGf2M2aUqdqbO7kYH+MHYe8MZ4O+jGyGAM7YJb/roY2yBHVcs7ZrtO9i+KWIcYvvwf/iHf3DKyTHY+sSOpe9///up9ABfFqzmYFcgIg7z8KAXApjBEhKHPY+c9RGK4Fd/9VennYUT+MY3vkGl/9jHPuaUk2N4n4rGWn4i/G29N5wB/oI2oyjKmKeMB2GwZc2KzYj4zWybs8QBZvsOtm9irY4WYxLbh7NjhMW3moUdSy0CmIXVHGxotohxyINeCGCmU2c3BgB8I/z85z9PP4PhrW99K31NRMfAwloDLIMGyxVXXEGlZ60yLWw4s4iibEd79pUW3pO3YAb4NsG2OYtlne072L6J7fsshyWxfTg7Rlx//fVUegvsWGoZrz/84Q9T6VnNwbqWWFYsLL73k6Z5AVxrpRoVG3IM4AWw97Lhz//8z7ve38JDDz1EX/NP//RPDjk5hmXZ5oUvfCGV3tuia+lEvMVmRAQFb6HWghDMSjYrc0QkC/YZlnbt7WfM9n133nknlR7gV2DZMcIyDllOamOwjNesAGbbHKtpLrvsMio9YDPeTJpeCGDGqms5ROILX/gCfQ3Dzp07qfTs8ZARWBq5ZTLCYNlAwVpNvAVwhO9VxiDqEqizS7aya2Ujn7cAZvs+S//KHjrBjhGWcegjH/kIfQ2DZbxmNcHBgwep9KymeelLX0qlB2wHJk2a5gXws88+S53C881vfpN+BusDzJIx5BjrKvKP//iPTjk5BrvMaAk2zw4abDiZjCeWSdCKWYKtSxEhnCI2/LB9B9s3sX2fpX9lN3dFjEPsMyJ8W1lNwApgVtNYdBN7mIcHzQvgb3/721SAcHYGauGHf/iHqfQ7duxwyokd1vdqeXnZKSfHYC0ObNxMgJ+1SgDnoIXf0Fe8XcYybuSz4C2A2b7P0r+yffjtt99OpbeMQ+xY94EPfIB+BgurCVjN8S//8i9UeotusrhNTJpeCOB777137PTsqWsWfvInf9L9Gd78/d//vfsz2CU9tsO1xNpkBw3vE8gyhJKZBBKo4mRkOzo5goh2zfZN3hZjgO/DWYsxG7IU4Me6CAHMwmoOVgBbdBOjy7xoXgA/+eSTdHBtlu///u+n0v/4j/+4U07s/Mqv/AqV3vvceYAP3s12nhYfJHZ5q4XoBhKnomVaqd/ebhls32fpX9k+nB0jLMcOs2MdO5b+wR/8AZXeAqs5/viP/5hKzx6TDfCHnnjQvAB++umn3cVaRkHL8rd/+7eu92ePAgV4HzL29JqII3u1gUwIkdEv2Tv+saV/Zftwdoy44447qPQAv9GOHUsjBDALq2ksAjjCiLYRJePyzyQppVA/cNu2bfQzbrzxRvoaT37rt36Lvua9732vQ06O8epXv9r9GrYzPPPMM6n0gL+PbkYLsPAhY9m1Ph60DFt2bHo29JvlUClmwzrAb7669dZbqfTWaxgsRyf/9m//tkNO7Lz+9a+nr7npppvYS26utfICbR2atwCz/OiP/ui0s7BpPvrRj7o/g52pWzZEsD5k3u4JlmskaHOg9zoeGd+TRPl4ZIuXbOlf2T6cHSMs4xBrNX7ssceo9JbxOpsAtugmgwCeOM1bgM8666x65ZVXjp3+M5/5jF9mjOzbt49Kv3fvXqecHOO1r30tlf5Vr3oV/Qy2s2ItupYlwD5GUIhA70lMi9bHwEnhbTEG+JPgWIuxZeMV6zf82c9+ln4Gy+LiIpV+z549Tjmxw8Y/vuWWW2QBZjn//PPx5je/edrZ2BQRFl3vmbRlRzDrQyZx6oPek2iZiPrdR5EdcSQ6O0ZERKZgn/HII49Q6QFeE2QUwKwuu+WWWyaeh+YF8IUXXpgu7NjHPvYxKn2Ej/Ell1xCpd+yZQuVvgV/W+s12WjhNwgxS0RsgvMmY7xkdoywjEPsWMeOpZaTZFlNwGqOCKMhq8t+//d/f+J5aF4An3HGGfiu7/quaWfjON73vvdR6dmlpPPOO49KDwAve9nLqPQXXHABlZ6NIQn4R1BoQQi28BuEEMdjadfZRHNE3+Qdyxjgxzp2LP3GN75BpQeAhx9+mErPao4IAZxBlzUvgCP4z//8Tyr9pz71KaecrMLOQAF+lsseUmHxt+2juOvjbxZCbJ4WrMwsEaHf2LHO22IM8AKY1RyspnnjG99Ipc+CBPAEYOP4sXEFzz//fCo9OwMF/C26EQI4o3jMmCchhMgomL37S8s4xI513hZjgLcaP/TQQ1R6VtNIADfC5z73OfqaT37ykw45OcbFF19MpWdnoAA/y40IiZNRPGbMkxBCeNOCYI4IfeltMQZ4TcAKYFbTWHTTa17zGvqaSSMBPITlVBb2HGx2hnj55ZdT6efm5qj0AD/LZTcfZBSOGfMkhBAtkNGP2ZIndqxjx1LLeM1qAjb8G6tpLLrpIx/5CH3NpGleAD/11FO48847x05vOdKP5dJLL6XSX3jhhVR6dgYK5DyylyVjnoQQQoxHC1Zmdiy1jNesJmA1x4MPPkilt+gmRpd50bwAvueee6ig0UtLS/Qz2CN42dkbazG2NKiMMXQlaIUQQpyMCMHs7WZhGa+9V5HZI6Ytuok9zMOD5gXwI488gk984hNjp7eILnYXJzt7izjhTGJTCCGE2BwRkSlYTcBqDlbTsC4TAChd5kXzAviJJ57A5z//+bHTX3TRRfQz2NkV6/QeEXGBRYJZCCHELJFx3IqITMFqDm8fYwCULvMiTACXUuYB7ARwBMACgAO11pWTpF0EcBTA9QCuBrBcaz3I3mcN5iAJ9thDwN9FIWLDWQsxIYUQQohpEhG+k72G1RysprHoJssBIJMm0gK8v9a6AwBKKUsAFgFcs176Ls21tdYD1vuceuqp1Mlol1122dhp1/A+RIL1KYoQsxLMQgghxPFECOBsodksuondBMduzBuHEAHcWW2fo9a6VEq5GicXrjfWWk+oFYb74LTTTqM2qbEzHwA4/fTTqfTeMXTZo5MtSAALEUvEBh4hxObIKIBZWE1j0U1s8ICZFcBYdVVYGf6wlDJfaz3p9sFSykKt9chm7nP66adTDt2MtXjwGQzeERcsjYMVzRnD1QgxKVqprxl/h0S5mCW8BW3EIVHsM1hNY9FN7Ea7L3zhC/QzNiJKAI+ypy8DmDtJ+vlSynYAN3X+wNd1Qnis+5RSdgPYDazulmQCTT/vec8bO+3A81Klj0ACWMwSqn95YMsiY/8n+oN3/YvwAfa+v0U3WQ4AmTQpo0DUWvet/X8pZT+AQwC2EtcfAHAAAC688MJ6zjnnjP1sy45MtkPv4wCgjXlikqhu9Ic+9pfCh77WDW+NYtFNjC7zIkoAj7L2bsEIdwYAKKXMrUV26Px813x/qfsAq6b/s846a+yMWhoI6z7g7W4Q4VPE/gYJFiGEELOGt7tBhHEom0YBQOkyL6IE8Ej3hVF+u53rwyKAqzZzn4H7UWHELELNu3JFNKhsSyoWJLJnF5WdmBSyGPeHjO6EGcc6b41i6b/Z8K4ehOSgs+I+93dn0b1+6O/lzup7E4BrB77bCeDgOPcZRSmFMs9nFMARszHvWa4lMoVE0WyichOzhKJr9IeI1dGI/THZNIflN0Qc2LURkRJ8VyllD7oDLGqtg6HLFrHq53ug1rpSSlnq0q4A2Fpr3TXmfUbCVGCLUGOveeaZZ6j0GSrKMC3MpCXUhBAiD3216LKwYxerObwFM+Afym0cwgRw56awtrnt8NB3u4b+PoJVgUvdZxSsBdhSkGzlYguevb+lwXqL7IjQbBK0Pui9CnE8crPIQUTIMRZvcWq5xju9hQyGvek7YQTgbQFuwR+nBRSaTQghZpc+ThQyukBEjI29sgBPi1IK9aIjwqCxlfHb3/42ld7yG7zdMjJ2bC2EZsuWHyHEiUS002x9bLb8ADF58nY3AHhN4H06rEVzSAAHkc0C7J3esnzBVsaIJUCFZtuYFkS8EOJ4MorHCFpwaYgYr73zxJLxOOdxaF4AR/gAe1fGCId070gTlhmiQrNN//5CiHgyRqbo6wa1bBEXLNd4G6wsYlY+wEEwhZNxs5a3YLZcw1beiA5dodmEECIHGa2CEf13hADOZtG1lJsswEF4+wBnE8ARLhARlTfDDHGYbBvt5AIhRHtkdIHImCeWiChPEREXsll05QOclFIKVfgRM5lsghnwn7VGbC7MSDbBLIQQQBuCtgWLbsYNkhEGsQz1r3kBDHCFE3GiiffBGRkFcERcwYhG672bVgJYCBGBt09vhIUvYhxqQQB7W3TlAtEILRwjbGlQ3m4WGd0ZMhJhMZZVWohYWtkQNutY+rIIAxSLt+bIuOHRAwngISIc973TZwyrErExjyViBiqxKYTISEbRkjHiQraQY0A+zZGxLo2DBPAQLQjgjKfZsYG7gZwb87zLooXQbEKIzRGxEsmSYcl6mIhxKFvIMSCf5pAA7jHZHMxbCOUGxMzuWbwHgYwW44x5EmKaZBzws+UpYtUvo3tgxhBifRG0LBLAEyDbZoKMm7sydoZ9tRhLoAqRmwgB0oJFN6ORhCVCAPfFp5dFAngKZKyMEW4WLBmtARlFNou3RVeiXMwS2qA2Hhn711YELUtGDTGLSABPAW8L8Gmn8cXK+ka14AJh6TyzRbOQz7AQ+emjYMm46pfxRFLLeC2XhskgATwFvDvDjBaNCN+rCAHsbSnPaImXD7AQx5NRUGQLrxkhgDPugYi4f0ZNMIv0QgAzjSRjRYlwqmctmxmXxSOCokccEemNBK0PfXxPGfvLjLTwnjIeOhExDrHXRBwi0UJ9ytBf9kIAM0QskXjfPyKsSkSeWohNnNGi2wIZOs9hMuapj7QgDixks+hGhONkyXjQVcaQdyyz2vf1QgDPugWYJUIAe1uMAf+yyOiXnNFizBLxXjMGp+/j5KWFDT99tcBl9LdliTDEeFt0W6hLFjKI5l4IYG8ilmG8aWGWG9FBs8t6GS3rLBkt6xkFc4YOPRpLWWcUtCwZ3YiyuYxFtKGIcaiPgrYvfZkE8BARlkpvwZxxE5zFsun9niIsdk8//TSVPuNkKuK0JXZw9Y5aknH5NiMZBQW7qz6ib4oQ5Wye2DYU4QLhPdZZyrqFTevefdOs9n0SwMJERCPP6ALhbQXJKIDZ38yK/ohnZLRit0DGIP7Pe97zqPQWUZTRtSmjTy9LRrHZgkVXjEYCeAJ4i5aMO10zznJZInYpZ4StTxECmL0mm2C2XJPRapJRUHgLWvb+lmvYPGXsL1kiJkcsGcdGC7LoTgYJYAGgjQofEaaHFVIZrSze78niAsFe450+owtExo2kES4QbJtjXSAsEzZvN4uMx7Sz6S0TC/Y9Zdw0nHEVT4ymeQFca00XBSKbxdhCxk1FEZbKJ554gkrvLdQyHjBiEcDeojzjBp4WiLCosfUpQmyyAjjjpmH2N1hOLGNhy47NU0S7bmG8jiBDnpoXwC2QsYFEWAWffPJJKj0rTr/1rW9R6S3XeAvmjKGGWiBjvNAIMrr4eE92IsjYTlnxeMYZZ1DpM9aljH1Zxjz1BQngITLMSjZLxOYuNv1TTz1FpQd48fj4449T6S0C+NFHH6XSs78hYtBghdfpp59OpbdYirwtahFWwRbCJWXcLMj2NRErFuw1bP+XsR+wvCcW7zaRMRpHRGQKMRoJ4CmQcQmavcbbOgvwAvWxxx5zvT/Ai2yL8GdgxanlmrPOOotKf/bZZ1PpAeDMM8+k0j//+c+n0rP+iBYR38JGJ++VIIDva1hXJbZvYts0wPc1LJY8efc13hMRwN9lzFJfWcs32zdlNBhk7Js8kACeAN6DRkbfVja9ZcBgBwE2PTtQAnxZeO9IZ4UjAJx33nlU+vPPP9/1/gBwzjnnUOm9BbPFKpNxUMo42WbbkLegZVd1AODhhx+m0j/00ENU+ogVCO/IK5b+lSViwua9kZkV2Ba8NxfOqmDuhQBmGkmE+4C3tdViPWCtod6uAABv0WDLwdJoWesm27mx95+bm6PSW6658MILqfQXXHABlR7gRTNrlY5w48h46p/3ZtUIKx/bD7B9GStmAeDBBx+k0j/wwANU+pWVFSq95RrWKMH24RZDDDvWeY+lAP+72cm85T15GwD64sYRJoBLKfMAdgI4AmABwIFa68pJ0t4A4G2jvi+lLAI4CuB6AFcDWK61Hlzv2UxhRizpefu2PvLII1R6gB80IlwB2EYYYT1lOzdWbLLpLWJzy5Ytrs9gLcaAvwWYrRsWC1yGDj0ai8GAFc2sQIiwAHu36+XlZSo9wItyVjCz6S3vlS07byu25ZqIVRTviW1Gq7QHkRbg/bXWHQBQSlkCsAjgmpOk3Qlg59CAsrfWum/tXt3119ZaD6z30For1UgslZFttOzMmxW0ER2PtysAwFv5WGH3ghe8gEpvueaFL3whlT5CbJ577rlUem9xCvAdrre7gQV24IuwtnpbpTO+V2+DBOAvsi1GDNbNghXM9913H5X+/vvvp9JbrmEnChGRf1grs8UqzQp/ts9nx3dLnx8RVm/DPEQ8pLP+PketdamUcjVGCOAu7Y5a6+GBz3YPCN0ba61j99K1VqpDjNgQwS65sZ2hpZGzgwa7pMw2QIAXmxdffDGV/pJLLqHSA8CLX/xiKj0ryllXAEvHk3EDmbeQitjA4/2MiIMwvKNrWJ7B/gY2T5Y25L3aZPGjZ/sa1rWJvb/FPYs1erDlYBHl7HjNCtoIv+SM0a0s7W7SREnwBQArwx+WUuZrrUtDHy8PflZK2YlVd4fhaxdqrUc2evCzzz5Lzb4tM2+2gbDpLaKche14WHFqEZsve9nLqPSXXXYZlf6iiy6i0gP872Z9eiMCu3vvwre4u3if1JbxNLsWBHDEDnbvEHYRET/YSaclkor3ASOsoYc1SAB8n/y1r32NSn/nnXdS6QHg61//OpWeFdkWgxVr7Wf7ZLbPt/gxZ4jXHSWAR00dlwHMDX846PdbSpkDsGXIF3i+lLIdwE2dP/B1w0K4lLIbwG5gtSNhGq5lQwRbGVlByw6U7JI1wC/Vv/SlL6XSX3HFFVR6ALj88sup9KzItrhAsBMF71ibERse2YHPkqdsRyFnPM0uQgBnPEUtW3qAt16xgjYilKD3b7CMQ+w1bH/MTkQA3srMtrlvfvObVHqAd6lh+/wIi3GGg4Cm74SxPu8CcO3gBwN+wCil7AdwCMDWoTQHABwAgPPPP78yvkusmAV4HzK24Fm/0Je85CVUegCYn5/fONEAV155JZWeFbMAb0Fg/WEtZ9WzosV7R3rEbnE2TxaLhveyoXcoI8s1GZclI3yA2Wu83TIsoogVXqxLQ0R0F+/IKxEhGlk3joh9H+w4tLQ0vAi+MXfddReVnvX3ZlfCLRZgyzWTJkoAj7L2bsEIt4ghttda9w5+UEqZW7MId77E6yq3Z599lipMy+DNdrhsA2EFLStOAeA7vuM7qPRbt27dONEAFncDtjNkB0rL5gN25s2KTTZcEtuxAXyeInzQI3Zzi37gHXsb4MVgxCqKd5gyts1ZLMDsZlh2L4plssPmKWLPBAs7OWfHCEuf3ycL8BGMcIMY4f/7HJ2bw/KIzxYBXDXug5999lmqM7GcrMXOQl/+8pdT6V/5yle6pgd4FwXWZcLS8bB+S6yl0hJqyHuXsncoI8B/+cwysfA+ZrWFzV0RRESmiNiQyGARIGx/xq4sRqzssOkjouywz2A3V7PjFsAbYtgVW1Y/ALy1n52MfPnLX6bSs4YbIGZv00aECODOUvvc353V9vqhv5eHfH1HbZy7CQMuEd0GuXVjAJdSqBmZJa4qK2hf/epXu6ZnrbMA31mxAsGyMYrtoNkwPRbfK1YAe1tbI07YYwWtZWmLFV7eS/UW0eUdazjjSXARm1+8XUssv8Hbp9ySJ7aPZfsBdiKc8fAjy/jOWoDZvSgWo1vGMJAsltXLSRPpA7yrlLIH3UEYtdbBEGiLWPXlHYzpuwLgxsEb1FpXSilL3X1WAGytte5a76GnnXYaXvSiF42dSctmLVagvva1r6XSswKbnUVbYDtoS/iZe+65xzW9JU/eweAjDhhhy44dvDO6J2QcACL8bb2tpxFkrE/e1ntLubHtmp3YsoLWIoDZa7wFM8BbjdmJMBta0/IMb8Fs6ZvuuOMO+ppJUzJuxJgkl156af2Zn/mZsdO/7nWvo5/xmte8hkpv2RDmDTu7Z0PDsE77lmewgtYS8YN9TxHRB1i8rYKWPsXbAhzhF+odfSBCAEfUV2/fU++6ZLkmo3uMd31l/aQB3t0gIhwnuweHfYblPXnzla98hUr/uc99jn7GLbfcQqX/zd/8zZtrrdvoB61D9igQm+bss8/GG97whrHTWwQwY2GOwGIVvPvuu6n0bDxF9v6Av6CNCNeVEVYMRgzG3hPxiPBe7DMiyoG9JsIHOMMRqINE1O8IFxyWjKEEvePPRri7sG2IjXcP2NwmGFgjnUXERxy3vBHNC+DzzjsPO3bsmHY2NgW7VPWlL32JfsbRo0ep9GwAcsuGM9Z9gO3cLINMtpO1LGLCe8C3iFnvAb+F8F4RPsARS58RkxHv+2dsQ6yl3NsXO+KIX9aIYdkz4b0vw/KeXvGKV1DpLZvQGSxGwAy6rHkB3AKs2GSXLyzPYK2zrJgF+M7Qe7C3XON96pVlIPa2wEUIYG+LcUaf4YxY3pN32UVMdrwFsMVXlc2T98lxGSOERMT3ZrH0x6z1lN1H1BckgKcAK1Bvv/12Kr3FuZyNiMDOciNcByKsp95+mxHWsYzWq4xL9d6weYqwAEeQzX82YxuK8Pf2fk+WPt/bim1xD2Rd6yJOeGRhx62M+5Q8kACeAGz4LVYA33vvvVR6y2l2bMeQ0doaIYCzCdoIv9AIvH0kIyy63lbsCAGc8QjUbILZeo032VYtLPGVva3Ylokw2ybYsdQyXrOagI0DzB57bYmvnIHmBXCtlVpKt1hP2WvY6AZsvDzL8lnE8aEs3p1hhP9stvQRRAi1PlqMWyFbRIRW2hCLdzlYJlPeLmYZJ3iW8ZrVBKxGYf2SLeFjLddMmuYF8GOPPYZPfepTY6dnxSnAb/DyDqVlEXbsue0ZQ2NFkG2wzOhu0Fe820REWUe4D2Qjon5nbHPedSNCALO0EvKO1QSsYGYFsMWKbTmIatI0L4Affvhh/Ou//uvY6S27RL3jI7Li1BKSxHspyTLL9T6AIWI5LKN4bCFP3hbdiA08Eb6CEasi3mSbdFrI2OZYMkZSYeurxS3D243NUjfYDeKsHzN7qNQXv/hFKj3Au1l40LwAfvzxx3HbbbeNnd4yAJx//vlUeraBsALYEiOQ/d3egz2Qc+DzFi0tDJQReL+njALYkqeMFrUIgZDp/lnJuDrgfWBNxNjoPZYCvABmT9hjV7UtFuAMJzw2L4CffPJJyv+FdRYH+A6UtdBa8sTCdiRs+hZOOLPA5ilDpzANvDfBZbQAe4f5A9qIVsCS8dCJjHhb+y11j7XQsoLWsnfFYjVmiBDArNvlysoKlZ7dlAfYQqNOmuYFcK2Vmv1YKjtbudjg3aw/juW0G7Yj8Y5va7mmhRBiLBmPHba8I28REjE58vYBjrAAR0wUMlrvWbxXdjJGpshoAc4Y+Yc19ljGa1YTsJqD1TSshdl6zaRpXgAD/jM4NuwJ64/DBr22/F62kbN5snSG7O9gZ9IWq7S3JbsFK3YLR/ZGTCwi8pRxY553vNqIzVoZ3TKyubtExFn3HtsBf3cDVg8AvIWWfYYlXjJLRNltmIdpZ8CbU045hXI5sHQirPBiZ2Ns5WXdEwD/5TDL0pN3rOEIa0DGsGbZQo5Z8BZqGY9zjiDjRCGj2PR+RsQkMtux7pZnRLiYsdZW9uhki/8sqwlYzcG+J4uYtWzWnzTNC+BTTz0V55577tjpI3ww2Rki26Ai3A1YLIOM9wwx4xJdxBJ0RrJFgcgYhiqjD3DERKGFUGsRE+FsBwdFlFvEqh/rq8oKWtaaC/CawNvdgF0RBnJEnDmpwiilnAdgi+Gey7VW3qbvxGmnnYYLLrhg7PTsTAngK5d3CBOLiPfuSCxCjQ2T4n1MseUZGc+2z2aBszwj4859b+tVxICRUWxm9FlniciTt6CN6F/Z/owdhyxhTtkYumwEBYsAZn8HOw6xdcMigNnoVh6sV6PfBeBGAGxL3NZdmwLWBcKCt3hkl2AirIIRy2HsIMCWc4SlnHVHiRB23qI8o09lBBmFF0vGPHmT8XCECP9Ztm+KmICxYyO7WSvC35ZNz1pzAd6I5j1uWcRsdheIm2qtH2VvWEqxWI3dKKVQhWMZiL13fbLpLVbsjKLI+xmWBujdkWRcKo1YZswmvCLea0ZrK0u2cgNi2pB3FByL2PSO5MO+J0s/4B2ui7XOAsADDzxApY/YcMb2HeweHNaiaxHAqS3Atda/Gf6slPIrAJZqrR8tpVwPoALYX2v9t4HrPuSSUyOlFKqhW/xOvWPoshZgS8fDPoP1i4p4r94DAOAf/SJiZ2y2SBaWa7KFfLJe431/WdY3JuNKkKUf8HZRYA0SFmHHCuCIDWfsWMeOpRH7Y7w1SkR99YDNwR2d+P1VAEdrre8qpfyQR8YmRSmF6qwsERTY2ZW3BTgivBfbUWVcZow4qtM7hF3Eb8gYmYIdjFnBknETXIQAbuF3RxyPy062IwRFNkFrOeiA9bdlrbMWAcyOdex7tZS1t0WXvX/GyFPjwL75tdq5HcDe7v+5c4CTE7Exiu082fQWAey9+crilsH+btYaYHHc9+4YIgZv7+gDlk2Y2aJfRERcaEEAR7ynjIcjZPS39W6nrGXTsuGM7cO9D6EC/KOWWEKEemuICENMBtjRc2sp5UEAW2utn3HIjwtMxxDhD8ZW+IyhrtjIF5bfwHag3hMRIJ+gtcy8WSI28mW0VHoTYYnPeLAKi/fhCBn7AYug8N4rkjG+LTtGWOq3d0QEy14UNkISqzm8/cOBHH04K4CvB7AbwI5SyvlYjfZw/8RzNWN4ixy28kZY4NjKa/EHY/PEinKLhYJ13GfTs51nhP9iRLgu73jJGcloAc5oycloAfZ+hqUPZwWwt0XX0r96G1Ys9ZudILF9eMRKpPemzVllvTjA3z1s5a21PgTgfQMf/dqI6/7X4Ka4WSOi4NnKG2E98H6GJU/e0TUsHbT3xryMoYkiNhV5W7Jb8AGOIOPvZvuOiDbE/gZW0LJ9GeAfEoy16Fr6V/Z3s3XDYu1njRisddYSDSGbRXdWBfN6b+X/LaUskfcrAP43gFQCmCnMiM1a2Xa8W4gYKFlrQMQgw+aJHZTYztAiHL0nLxETsBZOgusr2WLoWuqrd/2znNzF9jVseu8Dnyyw/Z/F2sq6KLB9uMUH2HuS10JM83FYTwBfB+At5P1qd10qZnV2soa3D5zlmohBhr0mIvyM9zHWbMcWsWTtvWJhucZ70pnx9LsIWjgVzftEPsD/MCPL4QisRdf7OF3Le2UttGzfFOFvy4rsiAgKsypQvVkvDvD/APifwLy4wIZBsz7DkwiLcbbfYL2GwWKhYC053n52ls7TeznMYtGIOCmLQQJ4dp8RsfGP3dPAbjizuA+w13ifMBrhCsVaWy0CmBW0rIjPaDBgydiXjcP0IxE7U0oJ2SmfCYsF2Ltzy3iyliUkjrflhx3ELBbgjAc2eFuZI1aBZnUQ2AwRgtn7JMKIAxtY66zlyF5vARxxqI+3v63FBcLbopuxD+8LvRDATENsoWJFHDuc8T1F+Hmyy4DegzE7EAP5ThGyXsMQEXLM2y0jgojf4B2vNuLEMm+LrqVdW34HA9tvWMSmd9QcS54yxtDNOP6yZOj/mhfAAGf9yeg+wJJRAGc8TSziPbHP8LYwA/zgzYpTiwsEew07KGVzsegzlhBfDBFtiBW0bJuz5In93WwdZ9ucJbqBtwC29E0Z98e0oDkkgANgXSD6KoDZQSljrM2MR/xmjK/MWpdYAWyxXrEDGTtQRoT1yRhD1xuL/2y2NmGpr6ygZZ9hadfe4eJYdwM2PcD76HofNQ/4u/tFhOHzRgJ4A0opCwA+VGu9aoN08wB2AjgCYAHAgVrrykbfrXM/qmOIiFbAEmHZ9D5HPsIHmMUyeHuHWotYvvU+AtWyuZC1XnlbESNCuWUbxIAcg9Iw3j7AlvrKtgnv+gr4uyhEuBt4bzjLGCM/wujm3ddETIQ9oAVwKeVtABZqrT9TSvkhADfWWtf12C+lbAewjFXRuhH7a607uuuWACwCuGaM704KU4EjZmMtDHwRO629T/mJsNh5v9eMx+NaRDl7DStyIiZsGQdK77KOaNfefvQRp1RGRFJhxeM555xDpWctuuz9AV7QRghg73bagh6wvNeISeFGUAK4lPL7AI5i1QKLWuvHSyn/B8BH17uu1nq4u36j+88PXbdUSrkawDXrfbfBPdNZgDM2EG/LT4S/bYQA9hYI7GBvgX0GW/8sosjbAhxhbfD2M44QwBGwZeFtAbbUDbYsWHeDiAgKrLuBt3sC0M8NZxnH9wgDV4a+iW1lN9Za/6aU8j0Dn61MMD8Lo+7Xid+TfldrXffEumwWYO84rBlnlBa8xaNlkMl2ipolegJr8YqwBrBlzS5bZxRFEZtlvCdsEZMdbzciy3tlLbTeS/uAv4U2QgB7H/rUip++t6D1NpJkgVUAF3T/HXybC5jc0cdbRny2DGBug+9OSinF/SjkbLtEM/rbWmAbLSsGLb5/2Rq6pdy8j5jOuLkwwmLs7b6SMTRbhB+9t/XeYvRgra3sMyKO7G3hiN+MG59ZIvoa7zZkWbnMUBasAL6jlPL/AXiw88FdwKofbipKKbsB7AaACy64gGokEdEKsgnmqGeweDfyjEdGe/vAAfwGHjZ9RGfo7RdqCUPlvXxrwXvTpuU9eft7s3UpwlLJikGLeGQFLVtfI474bUEAe0/Oo57B0IsoEJ3P7xJWIzFcCODXuiOTJ8Uoi+4WrLo+rPfdcD4PADgAAC996UsrU+kjxGM2wWy5JuMGHraRZ/T3jpiIeC8bRmwq8o71GiGAI0IJekf8sLwn7wlVxIll3jFxLaLcO5Y2m94igL3bRMaNpBk3kEVYmGfRAoxa6x0A3rf2dynlV2qtfzCh/BzBCFeHbsMbTvbdejfMuAnOO31ErMOMHU/EJjjvZWvvyZHlGREdFeuWwVoqLbFeWbz97CIEMPueHnnkESq95RlsWbPt2rKKknEDmfdGu4iNfC2MQxEuY9nGxozhY8eBjQLxNqy6PKy9/QLgfABmAdxtcFuuta4MCN3B764HjhPBJ3w3xjOs2RMzREQ5ewvaiJl3C4NMxk1w3uHiIgSw9wEPAG8BZusGK2gjQt5lPDgo40RYY3V/yFDW7BRtrtZ6nBW2E8Xr0sUBXuj+fxHAobXQaFgV1IfQuSwA2FVK2YPusIta62CYs/W+G0mtlRoELLOSiNBV3kTs/mbJGMPUewNPxjisEQexsG2ItRhHiM2HH143HPoJRMQw9Rbllo2kbBvyfk8R/osZ2zVbDhmPms/oA8y+V4s7A9tfeqePGEs9YAXw4RGfXbfRRZ3YPQxg34jvdg39vTSQ7vC4360HUziWDt27c/OOgwnktAawZDycw3v5zFLW3p2hpQ2x17BCjRXMERv5MvrRe9dvwH/zX0Q0mIxtiC0L7z4/4x6LCLz7fMs12dIDs3kS3IPdwReDfrfXAPiZyWVpsjz77LOmjRoM3n6hGf1xMnY8LBG7UL0FiKXj8RaPESfBsXl69NFHqfQRS/sZV468oxsAvD8sCyuALfWVrX8slnadcSMzSx/HlQhDTMZVvwywAvjtANaONV7jKiQWwLVWqrOKWA7zbuQZnMuHyZinjB16hPXAO9RVxLI4m54VOY899hiVHuBFNiuiLAMl2+68j9MFeJHtXTcs9dXbqBIhgL37JgsZhVQGS+Uw3nmK0DQZJjusAD5aa/21wQ9KKT82wfxMHNYCHGHKz9jIWTJU3mEyLtF5C/+Ik7hY8RjRhrwt5RarICuAWZFtsRizYpN9hiWCQjY/+ohVFBZLnrL1NRldzCLImCeWiLHUss9i0rACeFTIsaOTyIgXzz77LLWcmdEfJ+NyRMQzvBthRJiebIIZyLkxzxu2vkYsi7MhxCyWStYdgK3flveUTSBkjPgR4R6TcWLhvfkq47gV8Qzv6EURUUs8YBXAFaWUDwC4ufu7YPVQjP9normaMEyjsnQ83gI4opG30JGwDcrSQXsfSpJxM6J3zEnLNS28pwhB4R2Gz0K2w2Es9dXbKGEpa+9nZNyg621MshBhPfV+RoTRI8NE2OIDfBjABQOfXTi57EyeWiv1ojMe6ZfRqT5DCJNhMjSozZLRjSNipu59Oh17/4jVAfYZlnbtPWGLOE7cu6wzWq8y9uHekS+inuFNRH/pvVIYseqcYbxme6+3DR99XEoZOxzZNGBPgrNUXla0eMc6jAjhlHEDRYT7gLclMaMAjrCeZrMARxyPy0ZQsJQ16wLB/gbLe8pmAY4QwC24pbXQD1jIGLXJexKZMfKUB9RbHBa/HQ9OKC8ulFKoTj2io8roe5XNj9mC98BqucZ78LYIkGzWA8BflLPlEHE8rvfRyQBfP9jfYHlP2eLPZtzAE+FHz75Xts+POLmwhXElYhzKuBclw0b6dXvHUsr/qrX+28DfPz2cBMl9gEsplBUk4yYklojQWBFuGSzeVmzLNdkEs+Uatuws79W7fkRYgNkQYt4CBODrE2uVjrAAs0TUV29BERHyjq1PEZPzjBu+vSfn7CoN4D8Ba0EHjcNGNXRfKeWna62f6f5+O048+S21D/App5yCM888c+z0EZ0hSwtH/GY8FtOCt8iOEPEsrOUxwh/M2/fPMiix1lNWMEeIIja95T2xeNe/iBUL79PvAP/fnVGcRowR3isQFjGbzaLb5Ca4Wuu2oY9m0gfY2wLsffZ8xk1wGWZvwoZ32UVMwLwDtVuEXQaftmgi9kx4RzeI6MuyCRbhR0b/WW/DSoTm8IA+CGPtf0opVwC4AqNjA6eC6aQjOnSWiFluhtnYMBmtod5kPAkuY/xP79BsliN+s23kA3KKR28B7B1KC/DvkyOWxfvYX7Ywzlmu8Z5QWe6fIZIUK4CvBvAnAFBrvQPAHaWU/wPgKxPO18QopaRr6C00wghxmjFagffgnTEec8QKRDaLrqUusc+IyBP7XtnJkeVwDu+yzlhfWSI25nn3r33tmzKO7xFHG7Nk0GUbCuBSyvkAfqj786pSyvJQkh0APjrpjE0KNg5whFk+Y4PKKIC9Z7kZGuBmsdQN7/ifETN71qeSjVZgiW7AWo3ZZ1hEkfeJZcwx89ZnZNw0nPHAhowhxFi8BW3EJLKP47uFDBOFDUeRWutDpZT/AbAXwDxWIz+s8SCAfU55mxjZDsLI4PsyjPcSScaOx4K3b5R3ODrAP9h8xLK4d3xbZuOs9ZoWBPDjjz9OpbdcYzlumcFSX73bUEQoS28DQ8RYGjFG9PHo+L4Yk8Yyo3TuDm8vpfxQrfXjznmaOEzHEGG98m60EdZW1gJn6Xi8g81HNFq2PrFLyk888QSV3nKNt88w4C9ozz77bCr9OeecQ6UH/GPoWiIJsMKLFZvf+ta3qPQA8Oijj1LpH3vsMSp9hIXZu51GRCtg+8sIi3HG48TZZ2SMmpPRApxhBYI9CGPmxG+tlar0ET5C2WKeAnxlZBt5RBigjMth7Htil5QtAph9hsXPk4UVwKzYZAXt+eefT6W3PIN1mYg4WICtGxFh0LwNBhYLM9smIiIDeAtgtg/PuOoXYe2P2MiczQe9aQvwrMNU+ohwHhmPrPReDrPEbc3YCLNZdC1L0KzIYX9zxDHCrAA+77zzqPQRApiNAxwhgNn6Z8mTt0Bg25ylb/Ke2EbAlgM72WnBYhzxjIybML2j7GSheQEcYQHOEM5jkIxHUGY8Hcxyf9Za5L2kbLFesYN3RNB/b0F7wQUXuKYHgHPPPZdK34IAtmwW9D4BjxXAljbEXhPhA+z9jIiNpBmPac8YyacFzTETm+BawHsTXDYynsATYZ2N8AdjBz5WUHgLbMB/6dMy8LFikPXp9U5vuYb9zRGb4NiBLGISya5ysD7Dlvqa0Y+eJaPVLuIoeJaMG/OyETGx8KB5AVxrdXeByObSEHFkakaR7b3j3XINa9FlB0qLiM+24QzgLbqstZV1T7D8BlbQZhTAEfdn67i3X7zFPcFbAEfEV85YN7wjo2QU/RlDr7JkjGQxDs0LYGB2CyeSjO/I21cw4hQ17zxZOnRWALNCjXVnsFzD5ondcGZx42CvifBx93Zfsbwntiy865+lvnpHpogIg8bWDbbvi3DZiTAOCR8ylEUvBDBDxpAkGcN7sWQ8sMFiZWGv8Q41ZNmFz1o3va2zAJ8nNuZuxDGx3uGSIizAEa5NbFmwZc3WpYhIKhFHTHv3TREC2HtSaMlTK+NppvsDOd5rLwRwNh/gbII5gr4GIPeOr8xaxwDe4sUKCotFjf0d3tZWS331PjAkowCOOLHMu01Y6ivbJjKGHvTeSBURUSmiDbVAC4LZAwngTaS1IgvweGQLLwf4+5yxPnAWaysb4ou1AFv8Z71j4kb4h7MCuIV2GiFy2LJm65Klvnq7QkVMzr1P2LOQMeJCC+002/2z0LwArrX2pjA3QwtWZpaI5Vs2fUZ/W+8jfgH/mKHsQBlhgesr3n0HW5cs9ZVtE6wF2LIxz7vOevvntkIfx1ILGX538wIYyBelwfv+EZtlMoZBY7HkyTuCAmvRjfC3ZUW5xS/ZezNLxuD0GTaBDJNxYx4LmydLffUO22c5nMO7PkUcMe1NxkOiWogCkUHMWuiFAM6Gt6CNOIEnwqfSezC27GD3tuh6pwf4JWJWIFgEjnfED+8NakAby7cRAjib/3ZEJBXvyBeAv59xxBG/Getrxo153v7bsypoWSSAh4hoIN7BviMc/SMOnWCJ6AxZiy67VOrtngDwvyFiQuUtaL0jilie0YIAjhi8Mx6O4O1mYWnX3kdAR/jRZxReGcffbOHfMp6fMA7NC+BSClUhM/qFRoRwamH5NmLm7e0CwaaPsGKzRHSG3oLWIoC9BUVGAWypf94nZUX0lyxsnizv1buvyRjLOCMZrdIZ32uvBHApZQHAh2qtV42Rbnv35+sBvK3WutJ9twjgKIDrAVwNYLnWenCjZzOVJUJQZLRosA0k425xlohBhl36jBDA3ps0IgY+Vmyyg3eEBTjjoMT2TZb35H20Nnv/iBWOCPcs774pYxvK6NvfgstihNXbUj8mTYgALqVsB7AMYGGDdHMAttVa93V/7wTwcQCDonk/gEUA19ZaD4zxbMpqZ6mM2TrciEae8bzzCCsLawH2Plkro7uLpf55H2IScUhKxsHbe0NOxpWgiL7Ju91Z7u/d13j3fQDfTr0NN5ZrIsKWZnO7jNhc6EGIAK61HgbG+sHbAOwFsCZsDwO4oZQy11mBb6y1Um+tlEI1xIgjUL0tcBmXoCMGygj/RXYQ8N5AFnFiGUvEEdPegjajALbgfTxuRgEcgbcRI0IUefdlESfBRbS5CJHN4j1OSABPgVrr4VLKroGP5rvPVwbTlVIWaq1HxrlnKYVaGopwgfAWmxYB4j3zjrBeRcT/ZK9hB42IGJIZN5B5C1RvgQ34D5QRPsARy7fZQthZ3qv3ql9EWWfbzwD4+xlHCGC274gwxHj7xVtWLDKcypdKAAPAkLB9C4B9A3/Pd+4UN3X+wNdtJIQjXCBaiGEaIXJYslk0AH8XhYjwct4hwSzW02wCOMKPuQUBHLGk3EciLMDeGx4t/St7jffGU8sz2PcUsT/Ge9ySBXjCdP7AC7XWHWufrfkGd9/vB3AIwNYR1+4GsBsALrjgAkrURsQwbUGARPg5sY2W3dTBpgdilgEZIjp07w1nlmu8LboWAZzx6G4Wtp1a3pP3MdYZy4F9hmUlkoUtB+/Yx0DOlR227Ng8Rew7yma4yUJaAYzVjW6D7hAY8AVGrXWplDI/6sJuc9wBAHjpS19avWca2QRtxPJtxl3N3u4JgL+/d8TgzZZ1hP+styhn86ST4MYjY3SXiJMzvZ+RMU8RFmDvPj8ivGFGlzFvARyxcu5BSgFcStkDYG+tdWVN9HauD4s4PiLEWDADU4T/rPe58BYB4h1r0+IPxp6G5B1yDPC36LLlYOk8WWurd3ogn0tDhLBrQQBH+P61EKEmAu/JeUSf7z0Rjjhi2tvAAOTb1JYxkso4TF0Ad1bc5YFYvzsBHBzY+LYdwEEANwG4duC6nd3n61JrpSqkRQCzAtV7uTcivmjGDWcZY+iyRJyw572ikNHKErGRtAXhFWHFzhZBoZXjcb2tdhlX/di+zDIOeY+/EX14C4fDeBAZB3ih+/9FAIfWQqNh1ap7CMCBTgzf0KVbu3wJnSAupSx11uEVAFtrrce5SIyi1kpVFsts7IknnqDSewvmiF3NbEdiOdve2wKcMUxPhPUgYwgx72e0cHBLRiL2TGQMQ+Utsi1tyNt6HyFY2D48o7tBhMuid18T0YYsLjKTJjIO8GEcH9Fh7btdA/+/BOCkPUUX8WGs8GcD11AVMmIDj3eol4gg6hERF7w3nFk6dG8rX4T1wHvQyBhBIUJEZTzZLSPe7lYR/aV3qKuINuS90S6if804DrHvNWLPDqs5IsKgZWDqLhDe1FqpwmetuZZrvEOIRRzm0ceIC4D/hseM7gYRefIW/hmPTI24v/dSpkWoZRssIzacRUzOvZ8RYTHOGJmCFY+sHoiYgGWMjJKBXghgbxcIbwuZd0cF+EdQiDhghCXjIRIZLcARnae3ZT3jb4ggY57YsvDekNNCPxDxDPY3R0x02DHCMg55HxhiMbqxeK9OWXRTb3yAp8mzzz5L+dyy/rkAL4C9NxNYZrneERQiNpxFWPmybSCLsLZmtJ5mFHbeZLQAR5BtcmS5JsJXNZsF2IL35sKIjXkRq6MZQ7OxZOibmhfAAFeYEYO396w1o59TxGYZloiQd97pMx7Y0Ir1lMXbUhkxYGQYlDZLxvoaEd3F2y85IkSed/1rZX8Mmydvw03EgUweNC+Aa63uAtj7yF52hmiJuMDOcjP620ZEUPDe8JhRAEds7vLeVR9xwEMLRAhgbytfxG/w3mgXIYAzTna8jUkWvEN+WsbrbIJWAjgptVb3wczbousdDxfwt+hG+NllPDK6BQtwRjJairx9WzMSYeVr4T1FTDqznU4Xsbkw4v5sHfce3y3XsK6dGQ8Y8aB5AQxw4itiiYS16GZcUvH2zwX8Q71kjFebMeRYBNlETsbBOyMZ3TIylkPGNudNxnKIOLI3o8siqzkiXCAyuL31QgAzFdgigL1dFM4880wqvWVG6W3RtQwALRzZ6y1o+7rhrAVR1MJ7bYEIH2D2vWY8iTCjbzWLZWz09pW25InVBN6b4GQBTkophep8LJXRO0xKH+PhWq7JeGKZ4jWOR8alT28y5ikjGeMrZ3yGt99mxhP2Mm7Mizj2OpvmsOgmS8StSdO8AAa4JQlLQXof2ctWRsvAynZuEf62bANh00cIYO8TyCI2bWa0tvZRMLdCNkHbQtxqwL9dtxBLO8K1ydudEPB3gYg4zjmDW1DzAriUQlWWCIf0iBBiLN5LJBmtrRGHnmS0srBlHXEspvfgnZFsJ6IBOZYlN0vGOMAttOsWLMAWa6t33xSRp4wb7xUFIoBTTjmFmv1YglJHHArBEBFCx/t4SMszInyAsy0zRlg0+ig2W/jNWfEWqBEhnFqIpa1QguORsb/0tkp7H/4ByAIcQoQFOJtFNyImX8QSibcjfkRkCu/BPsKnLWPIu4zH43ovKUeEZmNpwX0g40EYEW0ool2zsH14hA8wO9ZFnLDn/YyIUG6WFdhJ0wsBzPjoWmYy2Sy6GQ94aCHigvUaTyJ8BSMsOREHBXjjbSmK+M0REwXvmLgtuEBEtLkWLLQZD27JuNrkLYBlAU7KKaecQs1OLGI2m0XXYm3NeGKZd0iwjIKCxdJ5su+JtTZkXL6NECAZl0pZMi7Ve98/Y4QGC9n6mojDPLxP5AP83QMtPsDZjrG26CaL1XjSNC+AWQuwJSi1d6ON2NzlbWXOKIAjyChyMm6+YsloUWvhvbK0IOwiiFh1ydjXsHgbMTK6QFjGa/YZ3uktuikDvRDAzOzEMhvzbrQRFmDvPGUUpy2QMa5lxECcURT1sY63EEM3o3DM6O/dAhHuLt4GLsA/T96CGcjRh/dCADOzk4glZe+YuxGb4DKGAYoIP+NNxh3EGQUCS0ahlmEAGKaFsmaJiKQSUTcybmrLRoTPOjuWWtwH2Gd4W2gtE7YMVuPmBTDAFU7GHcER1taMVpaI+LPZ6KM4jSCjAM5IRBg+lhbqeMZVkRbqa8b+MuMmzIwGqwzjdfMCuJRCvWiLr2o2F4hWlsK8LRotDKwZOpFhWqh/LYiDrGQULSxqdz70sQ/P6JbRlz0TzQtggGtUEbFhWxC0EdZZ72gFFvp4BC87KEX4g3nHDRbjEbGyw5IxDFXGUFfZonFYyBg3PQJvDREhmCWAA6i1ugvgbMsRGf1tLZ0IK6RaGMgyLtFlPDLVOzRbX3fhR5BNoGYU8Rn7gWz3B9ro8+V2OT2aF8AA5zCeMYJCxo5EG87yPIMl48CXLU8W60SGDj2aFtpQXwUwi/zoxyPj5CVjJIsMdbx5AVxrpQqzrwKYJdtyG5CjQc0CGd9TRiuLyEEf62tGMk52+jqutCCA5QIRQK2VOpnFEpQ6IkwZg6XBegdqz+gDJ8YjY1mzZMyT6A/qm3zIuD8mm4sP4B/KbVb7114IYG8XiGw+vRECeFYrvIhxB/CmhfqX8Te0YClv4TdkrBsZacG1zlLW3hoi4vCtDHW8FwKYsQBnFMAsGZeqLLQwkGUkQ8cTTSttIhsZ/UL7WA5iPDLWjYwW4IhNwxnG914IYEbUWuIAZyjIQTKebGQhY56ykbEcIlCe+kML/UBf26k3GX9zRuNQtvRZ6IUAZvx6M3ZUGd0TMjbyPpJxAIgg4+/OmCeRA9UNH1o4VCWCbKvUQA7/7V4IYItVl6GFRpgxTxmRKN+YFupSC78hggj/RRZNznOQ0RCTsV1nrEuyADdE3wRwxriWEbSwVCrGo4X62ldaCHmn+rcxGcehVsrNW6CymmlWx9IwAVxKWQDwoVrrVRukWwRwFMD1AK4GsFxrPdh9Nw9gJ4AjABYAHKi1rmz07GwWh77MrjaL93ttpTMUQhxD7ToHGS3ArZBNQ8yqRgkRwKWU7QCWsSpax2E/gEUA19ZaDwx+Xmvd0d1zqUtzzXo3qrVSviYZQ5KwWO7v3ZFYQmm1EH+2BWa1cxMiC33smzL+5ohxiCVjJJUITZNhXAkRwLXWw8DYDeLGWusJCTvr7+A9l0opV2MDAdylHTOnMWR0SM+I9+EcFlpYvmWRT6VomYxCLYIWfnfGGOURRIQpY5AAnjCllIVa65GBjxYArIxIN19rXZrUcy2FwlYuNnh3hCU0mx+z5RkRnWHGU/xYMnQ8w/RxYiHGowWhxtLCb47ojzPu+4gQp97jUF+MdBkF8HznMnFT5w98XSeEt4xIuwxgbvjDUspuALsB4Nxzz6UKJ6MLRMbliIwddERnmNEqzdJHsZlxoOwrGfsOlhZ+A0tGQZuRjON1NjfNqGdsRDoBXGvdt/b/pZT9AA4B2Ere4wCAAwBw0UUXUW85Y/y7COss27llqLybJWKy04JgZmnBit3CQNxXVHbj4S1oWymHbMcOW67xtujOqh5IJ4BLKXNrkR06P981399R1t4tGOEWMQi7Cc6Ct3iMEEXeh21YfoN3njIKtT4KZgsZB9dZHQQ2Q8Zy6CvZBG3G1dSII3vZZ1jCtGYTtLM6DqUSwJ3rwyKAUaHSRrpBTNL/N4qMsyu2Avd180E294GMu5r7isSgmBQZ3Q36Wr8zisGMGmIWmboA7iy8y53V9yYA1w58txPAQeA5a/DwddePcX+qM4lYjpDF2IeMg4b3/fvqxyzELNFC38QS4dvq7QpguYa16EZYpVvQEB5ExgFe6P5/EcChtdBoWLX4HkJ3qEUpZamUsgerrg1ba627Bm61q/vuCICFWuuGIdBYMhTKMBlnexEW4xZmrRkP84iwMs86Ev1+9LE+sWSMspORiLGxjyHE+tJGI+MAHwawb8R3u4b+PoJVgTvqPksD9zg8Ks0omMLMWBkjyDij9PYz7ksjFzyqG6J1Mh7I5O3bmtHa2gIZj70eh6m7QETQwrK19/29O5KMsYktnSEbwzkj2eqrEK2TYbCPJuOGswgXiIyCOaO/d4Y20bwAZn2AM8bQzRijN+PSE/u7LWI2Y1mIyaNy80MTqtkko2tdxIENLYRczGYEtF4zaZoXwAD3ojOGbmmBCBeIjKI842ZBIYTwtrZmPOChFd9+b4tuRouxB70QwN7+fN6Wyr4K7IydoXcEhQjB3Er98CTjO1JZzy4ZN7fK3zYHGd0DWWZ1z0TzAriUkq7wI/xnWyBbeDnLNdniBrdCtsNtImilrPtYdiwZyzrjEb8ZLbotuDT0ZeWyeQEM+HeI2XyAMwbijiBjB83S10gW2QaybPkR46M25EML/WsEGWPeywViNL0QwLNaOFYiTgfL+E4zLp9lsxgLIQTQRjSiCLJZW4F8gnZWxy0J4CkQEa6LRaeDjYe3yM5WV1sh4+Sor/RlcM2O2oQPEe4D7DXeqxyz2kZ7IYC9aSGsGYsEcw4iNvJlxHvwzrg6EEHGvoalhd+g/jIP2dwNrNdkun8WmhfAEZvgxMa0ILqAfMIrY93OOHjLf3E8MtY/72e0MonM1jdFEFEOGd0HMuZpFmleAGck4+wt40lwGTtckYMW/Bcz1u+MgjajKBc+ZBwbWTJGUFCbGI0E8BAtVBTLwKoZog8ZLcAZLbTetODSENGuM8arzdjXZGxDmoBN//5Zn5GxDWWof70QwBkLnyGisrOVMSLmLkuGBrVZIgaxjO0hW9m1Ul+zhWi0PIOlBYHdVxccFrkbzC4ZfncvBLA3GQpykIyDUkZaGGRa8V/Mhlx2xqOFfiNjuWU8RS3iPWUUmxkFbbZ2l7ENjYME8BDZKpaFFn6DhYzLtxk7hhZiPrNktGxmzFME2SzAEWTchOn9njIKwVYEcDZm1WAgASyaIUJQsGT0AW7ByuJNtvwAypMXGcVmxtWpjIJWiM3QCwGsRjV5Mr7TjINMCxujMrpMZBTxGSwa0agfGA921aWPdQnI2a4z1vEWyPBeeyGAhZgU3gNfBN6HmEQc7Znt/lHP6CMZJ5HegraFfkPtQWRHAlikpY9WkIjf7D24ZgyN1cJgbJlYZBRSLBndDZ555hn6GoaMG84y0kLEDzE9JIBFCC342wL+lp+Moj/Ctzrb727Fip1xw5m3KM/Y5jLmiSXj4Ud93UgqJoMEsDDRwqDRwm+w0ELMZ28sg14LS8QtLNWzZPQBzkjGiXDGyD8Z27UYTS8EMFOBW6i8GUVUC1bBjGS0aLQQmzhjdI2MZPwNLYjTjO06IxkNDBG0UNYZyqIXApghQ6EMk1HQsmQM7B6RpxbIaAFpoRwyRqbI+F6VpxxE9K8tTM5llZ4deiGAZ72zyiiAM1qAZ72cLbTSsbVQdi2URQu/QXVpdslYdn08OKgvNC+Aa60pGxWDBPB4ZFzKFD5kDI3VglUm46Rz1vtvMT4RIReFWKN5AQzk60D7OMhk3C0eQUaRw9LCb8jYhpSn2aWFNtHCb8gomDNOhPu4P2EceiGAGTJu1so4KLUgTjOSsePJmKcWViwytuu+vieWjG1CbEzEuNWCVTpCxGdoQ70QwExhttChS5z6kaHRZqeFASBiA0+2fsPyjBb6mogDRtRv9IeMIjujVToDEsCbSGu9poVBQ1YZH1oQjxnLoYXTwSJgf0PGFTOWiJjPGdE45ENE/+dddi3U73GQAN5E2jWydSQZLUUZRZEF78MOWnlP2chobc0omDO6QLBk9KnMiPeKRbZx0UofRXZG32oPwgRwKWUBwIdqrVdtkO4GAG+rta6M+G4RwFEA1wO4GsByrfWgQ3ZTkXHwjsC7kUccd9tHAZyx/kUM3hmXJTO6W2WsHywttFNvMrqWZDQORaD6OpoQAVxK2Q5gGcDCGMl3Atg5VGB7a637uv/fD2ARwLW11gMb3YwNgxYxa21B0GZsUBFH0WY77jajKMqYJ5YIF4iMbShbOVhoYdKZ8XCEjJPIiPG6jyKb/c2zevpniACutR4GNn6ppZR5ADvW0nef7R4QujfWWvP1VmJDMorNFo67zTh4Rxzj2sIkMmOeMpKtzUU8I2O7Zon4DRGCWe20XbL5AC/XWpfW/iil7MSqu8NxlFIWaq1Hxr3prPsitdChR7gbRFhnW3CB6KNFg6Wvm2EztiFvMvYDFrK1a0uf792G+uqWkZEM/V8qATzo91tKmQOwZcgXeL5zp7ip8we+bpQQLqXsBrAbAM477zzPLIeQbcAA+I4ko79tRJ6y3R9oYwNjxtjbbIeeceCLEBTZ2lBGEZ+xzUUIO+/6F9HmMobIy9jXZCCVAB7iXQCuHfxgwA8YpZT9AA4B2Dp8YecycQAALr744nQlL2urzzXZ/HOzkvF3e4tydlCKWCqNcOPItlsciJmoMmQUwBFkFPHe7TqiDWWcKLRQXz3ILIC311r3Dn5QSplbswjXWpc6n+ENYTrQjEGsM1o0MgrgjI0846aOjLSwbJjRBYLNUwuTyIx5ykjGvsl7HMoYti/jnomMOsiDlAJ4IGrE8GeLANYNo7ZZMi7VZ1zaj6i8GS0/3px66qlU+r4u6WVEG/NykHHiHNFfev8Otm+ykDFsH0tGC3C2CCFZmLoA7qy4y0O+vgsAVoaS3oQBl4hug9yGMYBLKVRhZuwMWxDAGTeBZCSjf25GS1EfyegC0Vcy9n8sLeSJbRMWEZ/RAuz9jAijRwbRHBkHeKH7/0UAhwZCnS1i1Zd3MKbvCoAbB+9Ra10ppSyVUvZ032+tte4a8/lMXsdOa70mYgMZC9sxZLSyZCSjoG0Bvaf+kLGs+9j/ZfwNGTeQsXmynAjJjtfeRokMYtZCZBzgwwD2jfjuBBF7sgMuuogPY4c/W4MRkBnj1bYQccFCC8swLfwGlohQQyx9Hbz7+Lszrpj1sRwiiOhfvZ9hKYds/eWsxleeuguEN6UUdwHMXsPO3jJuOMtIRmsAi37D7JJxoGTpY9ll/M0Z88SS8TdkbEMZI1OwVumImM8eNC+AAa5wLAXpbaHNaG1owbIpn0ofIso6YyhB74gLGQfKFjZ3Wci28bkVMv7ujBNb72dEWIAz7BWRAN5EWus1few8MzrJZ3xPLBknFhnJGEow2/0tz8j4nvqK3uvGZGxDGckQoiyC5gUw6wLRggCOQMIrBy24TEQQ0ea8Q9hZLCYZ+6aM/Zk3LfT5YjwyiuwWQjR60LwABvwtwOrcxCzRl85tkIxttAWxmTFPGcno7y36Q0aXsQz0QgBnC4PmXRkzCpyMeWqBjDEkM/oAe/vnWp7RAhknCiwRp4Nl3BjVR1p5rxnzxJLhNzQvgEvJdxAGS19FTh99gCPKIWNZe9OCUBN5aKENZayvGfv8FvKklaDRNC+AM5LxZJmMefKmhTxlFMAZY0JKAPvQQhuK8K3O+J76SMY4wBYy5mkWkQAeooWKklEAW2hBFHmLzYiyzhCuJpo+ilkLGS1qEUQcFZuNjP0lS1/rqxhN8wK41jrzFbgFd4OMZHyvEQLYcvQmQ0afyojB23sjSMbB2zI5amEHu/fBAmxEETEerYxz2VwgLGQoi+YFMNDP5X2WFn5DRrwHb28xC+QUIC2Q0S2jj2XXgmAGJJrFyZFlfTS9EMB9XMIVOWDrHps+ow9wRjJaQMR4ZKyvGVcsvI+vndVQV9lpoX+dVXohgFXBpk9GARJRL7wnX60I4GyhrjLW11bI1h9ny48V7+NrM65YRNBK/RAn0rwArrWGLBPPOn2c3Vs6Nm9BKwvweGQcKNmyy9jmWlgta6F+W/Cuf94W5qxk9KMXk6F5AQy008GNi0UceDfCjJ1hRMfjLWj7GnKshd/Q14GvBYHq/Rss9Y/tY70Fc0R4uQgyttOMbWIWaV4A11pTVmCGjKcIeS+3ZSVbDN2MERf66DIB5LRKe5NxIM4osDO6W3nHMs5oiImgr/VvFmleAAOz36gyhnyKIGO4pGyCNmPn2YIFWIxHX11wWogs5D0uRpwOlnFynnEcytiGMtALAcwUfkbrVcblswgyDjLZBG1fO7aMFuA+kvE9ZRQgLZDRxz1jObTgxpZRB3kgAbyJtFYyWtQy7gj2JqMAzrha0cfJkRiPjP1lRlpo132NAsGScQLWQhvyoHkBzJ4E10JH1VcydjwZ61MfiZh0apAZD72nHGTcBNeCu5/GofHIoGuaF8BAvlisGTeQZbQGtBB+poVoBRlpYSOfGA+VxXi00F9mPE6cpQVx2pcNj70QwAwtLOnNamWMJqPvX8aIHxl9/zL60YvxyOa/mDEeeMZxqC+iaBbINimc1ag5vRDA2TpcMR7yc9qYDJ1Ii/SxLonZJlt/mXEjVSv9Zbb+aVZdxnohgBkyFMpmyWg9yPiMjB10tvsLIWaDFiLOtLCy08c+WRZg4UZGsZkxT330UZvVjmezZBz4hA8tGCVYMvavLfTHGQVwxjz1BQngIWbVlD9IxrAqrfh2sb/j1FNPpdJL2AkhLLSwCY4lo8iOIFu/ny0/49ILAcwUToRFLZtgBvKdcGYhQjy2sEQnhBAZXcwy9pfem4Az9vl9sUr3QgAz9NUC3MKuZpa+ug/0kYwDa1/J2BeIyZPR6GHBO2pOhCHGm2z5GRcJ4AnQggWYpYXObVYb7WaRGBQilhb6/Agy9vl97C9b+A3j0AsBzCxJRPgUSTALcTx9HGSED+prxMno66pfxt+Qwbc6TACXUhYAbO/+fD2At9VaV06Sdh7ATgBHACwAOLCWdr3vTnIvqvAthZJxc4AQQgghhDcW3ZRBlIcI4FLKHIBttdZ93d87AXwcwFUnuWR/rXVHl3YJwCKAa8b47mTP3+xPWBfvM9IzVJTNknHpqYX3mtGiYZkQZrAGDJLxvWakBWtrRgNGC+8147HDrbRR77EuYiNfhrKIsgBvA7AXwIHu78MAbiilzA1bbzsL73PUWpdKKVcDuGa97yaV0RY2wVl+g7eIj6CPArivZBz4MuaJJVtfBuTMk9iYjPU7Qqi1EAWCZVZ/Q4gArrUeLqXsGvhovvt8ZUTyBQAnfN6J35N+V2tdOtnzs4VB875/xtlYxjxZyNa5ZXxHFrxDDWVEws6HjJNzlowrZhFGkhaMGH0drzM+YyPCfIBrrUcG/nwLgH0nSbplxGfLAOY2+O45Sim7Aezu/nzyve99761MXkUTvADA/dPOhAhFZd5PVO79Q2XeP75j0jcMjwLR+QMvrPnxelBrPYDO3aKUclOtdZvXs0ROVO79Q2XeT1Tu/UNl3j9KKTdN+p7TWEdcBLBrne9PsOhi1fK7ssF3QgghhBBCbEioBbiUsgfA3lrryqgNcB1HMMLVodvwhpN9N+m8CiGEEEKINgmzAHehzw4OiN7tA9/Nd64RJ4jZbvPb9Rt9tw4HNvhetInKvX+ozPuJyr1/qMz7x8TLvASdpT0P4OjQx0u11q3d9zcAONT57p5w2MVa/OCNvhNCCCGEEGIjQgSwEEIIIYQQWZj9YJpCCCGEEEIQhIdBi2LYVQLAgZNsuhMzRillAcd8yF8P4G1rZbteuatOtEEpZX+t9ZqBv1XmjdLtHdkCYAlYPVSp+1xl3ihd+W3HatSneYxZtir32aIbxz9Ua71q6HNTGZvKv9ba5D+s+hSv/f88gP3TzpP+TaRc5wDsHvh7J4Cbxyl31YnZ/9d1bHXoM5V5g/+6tr1noOzUznvwb63MB/5eVLm39Q+rE5wT+vLNlLGl/Jt0gehmAs9RV6NHXD2l7IjJsg3A3oG/DwNYKKXMrVfuqhPNMI+BuN8q86ZZrN0m51rrUu0sRSrz5hk+JGsOULm3RK31cD3+dGAA9jK2ln+TAhirM4uV4Q+HX5KYPerqEujgQSrz3ecrWL/cVSdmnFLKzlrrwaGPVeYN0i2Prgz8/yAq88YppRzqjBrbAdzQfaxybx9rGZvKv1UBfMJhGRh9ipyYQYZmjm8BsBYKb71yV52YYbqObNSBNyrzNpkHsNz5AC+VUvZ0YghQmbfOLqyW4x1YDXV6uPtc5d4+1jI2lX+zm+BE+3SHpyzUWoeXzER7LIyw/op22QJg+1rbLqUcwKogumCquRIRrLm5zQPYX0pBVbx/4UCrFuBRyn8LRpjIxUyziOPdIdYrd9WJGaWz/B0+ydcq8zZZxupubgDPuTit+fmrzBtlbTm78xE9AGArgHd1xg6Ve/tYy9hU/q1agI9ghEm8Dh2lLGaXUsoeAHtrrSullLlugDxpuZdScLLvvPMqJsLVXRkCq0JoN1ZFscq8TU7YIINjA53KvF2Om+x2ZXoAq2Wqcm8fUxlby79JATzwQgA8N6u8fno5EpOk8ws8WI/F+Nve/X3ScledmF0GfAABPBcH+MDA34PfqcwboCu75bXJ7ZoPeNfmV1TmzXIYx2K5PseakFG5t411DLeWf7NHIQ8HRZYPURt05Xp06OOlWuvWge9HlrvqxGzTLYPuxqrryz6sxnlcUpm3SVfe78Jqe98K4Np6kqD3KvN26Awca5tetwA4PCCAVe4N0Lm1LeBYX36onuSQm3HL2FL+zQpgIYQQQgghRtHqJjghhBBCCCFGIgEshBBCCCF6hQSwEEIIIYToFRLAQgghhBCiV0gACyGEOI5SykIpZbGUskBet7O7bs4pa0IIMRGajAMshBBiU2zDQNixcam1HhwISk9dK4QQkcgCLIQQQggheoUEsBBCCCGE6BUSwEIIIYQQolfIB1gIIcagO75zP1aP7wSAHbXWXd3n8wCWAby+1rq3S78AYDtWj+YEgJtG+dQO3PcggEMA5gbvM3CvbVg9HnYe3fGw3Wazq7vP5wDMjzoCtDte9l0AlmutO0opewBcA2Cx1npgg9+9AOBDXR7Xnr+1+3uu+40rG91HCCEyIQEshBBjUGs9XEo5DOCqWus1pZSbuvPn99ZadwBAKWVLKWU3gOsBfKjWelX3+SJWheMJInHgvkdrrYe79Cil3NAJ7HmsCtUda9eUUm4GcBWA3VgVw0e6z3eeJO8Hu2d8vBPNS7XWrWP+7iMDv/vAwPNvrLUeBHCklPLgqN8mhBBZkQAWQojxWQFwFHhOGO4BsDIULuyq7r83rX0waM1d577LA+kPllJu6P7ciVXL8CDLneX4IIBDpZQlAIdGWX8H7rlSStkL4OMAdm2Qn2EeQPe7156PVWuwEELMJPIBFkIIjuWB/78Qq9bUI92/A7XWazb7gKE4uhd2/0bmpbPk7gVw4YBoPhlL3b+RlmIhhOgLEsBCCGHnOgDHHRbRWWYPY9Vnd/jz9dgy8P/bccylYP/wMwBs6dwl3gWsWqM3sjJ3onqh1roLwI4x8iOEEM0iFwghhBiDgU1t86WUlVrr4c4NYrFzhTiC1U1hh9fcDYY/3+ARWztROofVTXDXAEC32W2x8y1ewqoYXnNheGDA73cFq2J5VN53Y9VKvCaSDwG4oZRy7XpuEwO/+y1Ydbu4Cau+zNsAXNO5VFwNYK6UsmejewkhRBZKrXXaeRBCiF7TbZJb21Q2ddY28rEnwXXX7gRwpNYqH2EhRFrkAiGEENNnbtoZEEKIPiEXCCGEmCJDrhWynAohRABygRBCCHEcA36/h9ZiE4953W6sHpJxrcV9QgghopAAFkIIIYQQvUI+wEIIIYQQoldIAAshhBBCiF4hASyEEEIIIXqFBLAQQgghhOgVEsBCCCGEEKJXSAALIYQQQohe8f8DmemWxCKmCSQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 1, sharey=True)\n",
    "axs.imshow(d.plot(), clim=o.plot.clim(d[:]),\n",
    "           extent=[d.ax_info[1].o, d.ax_info[1].last, d.ax_info[0].last, d.ax_info[0].o])\n",
    "axs.set_xlabel(d.ax_info[1].l)\n",
    "axs.set_ylabel(d.ax_info[0].l)\n",
    "fig.suptitle(r\"Nonlinear data: $\\mathbf{d}$\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a9821c52",
   "metadata": {},
   "source": [
    "## Instantiate the Born operator..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c0c85a04",
   "metadata": {},
   "outputs": [],
   "source": [
    "B = b.BornSingleSource(velocity=model_smooth, src_pos=src[0], rec_pos=rec, args=args)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "145b636a",
   "metadata": {},
   "source": [
    "### ... and cast it to an Autograd function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6a9970ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "B_ = o.torch.AutogradFunction(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5616b875",
   "metadata": {},
   "source": [
    "## Learning paradigm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c546446b",
   "metadata": {},
   "source": [
    "instantiate the model as a \"learnable\" vector, i.e., a new kind of `VectorTorch` that encapsulate a tensor attached to the AD graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9a46d326",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = o.torch.VectorAD(torch.zeros(model.shape), ax_info=model.ax_info, device=0)\n",
    "\n",
    "m.setDevice(None)\n",
    "\n",
    "assert m[:].requires_grad is True"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0770bb2e",
   "metadata": {},
   "source": [
    "instantiate the optimizer that will work on the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "2b57fde9",
   "metadata": {},
   "outputs": [],
   "source": [
    "opt = torch.optim.SGD([m.getNdArray()], lr=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bded8e00",
   "metadata": {},
   "source": [
    "define the objective function terms: the mean squared error..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "2f8fbadc",
   "metadata": {},
   "outputs": [],
   "source": [
    "mse = torch.nn.MSELoss().to(m.device)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f1cdd2c",
   "metadata": {},
   "source": [
    "and more complex terms (the TV, for example) that can be based on OccamyPy operators!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "10550b86",
   "metadata": {},
   "outputs": [],
   "source": [
    "G = o.torch.AutogradFunction(o.FirstDerivative(model, axis=0))\n",
    "l1 = torch.nn.L1Loss().to(m.device)\n",
    "tv = lambda x: l1(G(x), torch.zeros_like(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58dbb563",
   "metadata": {},
   "source": [
    "### A note on the computation device\n",
    "\n",
    "The system is designed to move from CPU and GPU seamlessly. You can instantiate `m` on a device (let's say a GPU), and then wrap an occamypy Operator that is written for CPU (in this case, Devito's Born). When calling the forward and adjoint functions, the computation steps are:\n",
    "1. move the input vector to the same of the operator (i.e., GPU -> CPU)\n",
    "2. apply the operator\n",
    "3. move the output vector to the input vector device (i.e., CPU -> GPU)\n",
    "\n",
    "A possible use case would be a Deep Prior scheme of imaging, like $J(\\theta) = \\Vert \\mathbf{B} f_\\theta(\\mathbf{z})  - \\mathbf{d} \\Vert_2^2$, in which:\n",
    "* $\\theta$ are the weights of the CNN generator $f$\n",
    "* $\\mathbf{z}$ is a input noise tensor\n",
    "* $\\mathbf{B}$ is the Born modeling operator\n",
    "* $\\mathbf{d}$ is the CSG data vector\n",
    "\n",
    "After the optimization of $J(\\theta)$, the solution to the imaging problem will be $\\hat{\\mathbf{m}} = f_{\\hat{\\theta}}(\\mathbf{z})$\n",
    "\n",
    "This said, we have to move the the vector `d` to the host of `m`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "460c5109",
   "metadata": {},
   "outputs": [],
   "source": [
    "d.setDevice(m.device)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8354954",
   "metadata": {},
   "source": [
    "## Optimization loop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "a8fd5ebc",
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs = 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "1295038b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch 000, obj = 1.44e+00\n",
      "Epoch 001, obj = 1.02e+00\n",
      "Epoch 002, obj = 8.32e-01\n",
      "Epoch 003, obj = 7.27e-01\n",
      "Epoch 004, obj = 6.63e-01\n",
      "Epoch 005, obj = 6.20e-01\n",
      "Epoch 006, obj = 5.89e-01\n",
      "Epoch 007, obj = 5.66e-01\n",
      "Epoch 008, obj = 5.47e-01\n",
      "Epoch 009, obj = 5.32e-01\n"
     ]
    }
   ],
   "source": [
    "for epoch in range(epochs):\n",
    "    \n",
    "    # modeling\n",
    "    d_ = B_(m)\n",
    "    \n",
    "    # note: things can be a lot complicated! Deep Priors, CNN-based denoisers, ecc ecc, possibly on GPU!\n",
    "    \n",
    "    # objective function\n",
    "    loss = mse(d_.float(), d[:]) + 0.1 * tv(m.getNdArray())\n",
    "    \n",
    "    # model update\n",
    "    opt.zero_grad()\n",
    "    loss.backward()\n",
    "    opt.step()\n",
    "    \n",
    "    # log\n",
    "    print(\"Epoch %s, obj = %.2e\" % (str(epoch).zfill(3), loss))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "f3484b7b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsAAAAGmCAYAAAByajI1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABPXklEQVR4nO3dzY9dWVrn+99KA1X55ox0FvWGqlQViSgVQqgVmTkCBjR2IzFB6rYz74Rh2c2AAdKVs+v+A8ipRgz72vwFmc5uwfDKrp4xQKTjXqErSiAcyaVEoYIuO/KlKpME17qDs4/z+Hg/K86z4tkv56zvRwrZsfY5e6/9vmKvZz8r5ZwFAAAAtOKJqSsAAAAAjIkGMAAAAJpCAxgAAABNoQEMAACAptAABgAAQFN+auoKAECflNJFSeckvSTpZs759sRVAgDsiEQaNKxKKV2V9ELO+fWp67KrUkr7Oeej7v97Oefjias0GWtbpJQOJJ3LOd9OKe1Jejfn/Px0NQUA7BJCIPBQSumypENJb3YNYQzj2sr/v51S2peklNL5lNKtlNL9lFJOKd3tfl/+3E0p3ewah1sjpbSfUrrTrdPNtcm920KLJ79XJGmlUbyvkXX75Hr3c3Hs5WOeNj0uUkpXV36udedC77W1m35t5fOXl+Urn9nvPrO8PqzO+9ryO1NKKR1016s7ld+7O1TdgFWEQGDVWyuNjaOJ6zI73Y1oP+d8aYj5d138t1NK1yVdlnQ95/zGyvL3JL0r6WJK6aWc8+EQ9YiWcz5KKV2StPGNbbktpIfrreWT4pHdkvS8pJeXBSmlyznnG0MvuFvvVyVdyjlf6Jl+VdJx9+ve6rFSmGfvMTzWOlnLGvrcGsBjx8W67jy+vnqe9vwBuOzt+GNJ31r77H73+Yd/+HXnwOvdd47W93nXIL85xHbc9NjJOR92n73umX/t97xOOq/QDp4A46HVrviWu+ULbkl6c6qFd/tkeaO5MlU9Kt07xXevSfqNqIpsaqWRcZxzvr0Sgzz4TbNb9quS9rR4Gr4+/aok5ZxvdI2Pw67BdRLrGB6zIdC3rEnPLY/CcbHu1Z4/Ur/V87nvSHp9/bNdY9fbiLwiaW+gJ8GeY6f2fD/NdeJEJ51XaAtPgIENzeQlrL2pKzCmrnv5ek9DYizHq790DYvBQzG69T0sdK9/W9LXVz5/O6V0Syf8YdR3DI+1TqVlzeTc8jje4DN7qzHu0uKP2JTSXyx/7/5oObLWv9uv3kbhTS3+aAx9oj/1sRNhg/MKDaEB3KDur+Bl19rrko4kvSJpGVt6RYuG1mvdZw4kXRnjJtV1T13u6nRBK42flNJ5LS7syycjFyQp5/x6971rku5okTXgrhbxzFdWu+xW5iEtnsbsq8s0kHO+0l3Q72mx7n+wsuyD7nv7OecXC/Pa10AvEa503R0vlzvXfdltx0tdfY4r53Fe0mEXQvHwqVtgHfdkH2sHWmy7Zczmcj0urJRpLURl+bn9rq5vl47ZU9R7X4uQh+OeaeetfVs4hsPWaYPz97Fl9dWrZ/nqpr+xMq8Tz73uODzS4hzYl3S8SahHYbmPHRc557eN2RxKupVSeuR8WwtbOK8u1Kfg2gnT170l6XpK6cD6w7FrAC7ne0WfXj+OtThv72nxZHqv+31Zj42Pne4zBzrFdXGI8wd4KOfMT4M/WjSEshbdWue7sr2u7I6kyyufvS/p7kj1Wl5kl7/f1eJmv/z9Yle/8906XOvKb0q6uLIeZn27795drvfKcq6tL6dnm93dcF4HheXftNa3K7u+sh+udT83u/1wc3V7zHFfdvXNa/vxcld2c+2zvduiW6f7qz8THGvW/r7TM6+Hx1/3+63lMWAdsxvWse84PC8p93z2/modjPkNuk4bbFNrWX31url2Xu1LurU2L/Pc6+p4ee37l9eXbax3abmP1dWYz35Xn7x+fq58Jku6WnHs3tLijwtrej5pXbtzcnW9rurxa+DqvvQcO8vzd+PrYmEZ1edPYTmPnVf8tPfDE+B2HXf/7ufu6URedM9Ji4vU6lOSI336RLHXhvGH0uKiU3oCs6/FhW75maPu9+VTluOufsunJssnHBe1eIKxXI/99e7HFfe0st4ry1l1qM269qx57a/Urdab+dGnccunvd9JKV1aWbfj7t+QfXka3RPAq+qe3K5Mcj1xzosnV0OnPTvpWNtI90T2Yn705aCbWjxZuyL7mI12T0FxjadYp6hteqBF4+nh8vOiJ+DcylPuTc69Symlt/IiXvcopfROwHI30h3/L648KT2vxRPhS9l+ajymt/R4jPHD66gWaQhP8+LpXtB18Vg9x1rgPQeNogGMvouR+6KXFy9fnNryxtM1pJbhCes39b76HWotPvaEi/f6tGM5shRsMK/wFyzy4i3pK1o8/bkl6cW1j4Tsy6579LUNPvpmz418+Vb87LOIbHisbeK8pOOuobP0oh79A2qM7RF5zFWtU+A2fblv/vo0tOL2yu+rjpfLy4twjSuS7qeUDrX2B+Upl7ux/GhGk2taZHxYnjeHevw8Xnb9L8MtzmkRjnVilo/uu3srdS7V6zildJhSutidx8fd960HB16R18W+Y23bXgTGzNAAxqx0T2C+rUUD7y31X8SPe8qua5FH9ltadO3t5MUxL16KkRZxd1E3qvVlvC3nE7shpIFHgtvwWDtpHvta/OG1/iLTel2P62rZy6rnXmHaxk6zTt5tWjiG9zavsS3nfGH5VFfSlZSSTmhIhiy3a4SeX/8DMS/ipK+mTwd9ua3FU9f1ei9TIr6gxZP4jRq/neUfocWn3Z03Jb3W/YHwjhbX0YsppYeNdo+hrkmKPX8ASTSAESSiO6q7aXxH0m/kT1+c2ev+3eTC+i113a0DXYQnlx4dDGLQlEGVljfdvdPMpGu0HHdP8fa0yH8cFhIRcKwtHWjxFO/bfcvIA6QT7Lrkj/vqGfRHQtU6VW7TA/U3km/3LV+Lp8obpUtLXX7a/Omb/ze6+pUak6de7opX1P+H5MOXObsG8cWVp7Dr9uRv/F3RIq3aJt97W4u47Vs55xsppWMt/ng5rgwbsPZnOEIgcFo0gBEiqDtq+Xb7alf+ssvspAvrS91FrvbJ5V7l94ae17rlm9s3hmhcnVbXrXpD0uW1Ro837dByJLjby3jm4KdLtcfaMo5x+cfIMkvFOz2NmFd1+lRUVpfxH2gl1rZ7Wl67rKh12mSbPrasvhl14T63V+Nuuz+KZDQUl/ZW/58eH3TjpLCA2uX2uZxSurX6R0kX2rA+nwuSbqaUjtc+uydnyECXjWFv0yfG3X5+GIbR/S5tfg3baH92Np3nRk55zyEHMMgC0eKPFjejW1q8KXxXi0bV3krZfXVv5+vT7APLt5j3Bq7bMuvB+e5nv6vDxe73W119rurRN5QvrtT9rtayH6yt+02tvH3dzet+N+/za59Zbpu+75XmdUfGG/myMx+srt9y39xa+Xlsvea6L7vtcEuLLtXr3e/L+lw/aVuszWtPw2WBsI619WPgoOd7l3vmd7n7/jIjiXnMnlC3/e7zd1bqcHHtM1eXy9IGb8b3Ha/R61TaptayNqjX5e7nWmldtHburXxvuY2WjcNNj42TlvvIcdFzzF7utsHV1Z8Njsery2Uvt/facbHMsnJ3Zb7LEdRqMkpcXTu+HztOPcfOJvtmg+NzmVWk6vw5zXnFTzs/qTsogK3VPSlZDiV6vPICzjUtGliz6v5KK0OVdi/FXM87GrJxkk22ReoZUhYAgNN4YuoKAAEua9FAOpYW3fBdY+l1LV6gwpZK048EBwDYQTSAsQsO9eloRauW3WfYQunTkeAOU0oHK+mdAAA4FV6Cw9bLXWqw9eFLVR6mFDO2MuiHupdylHMeemAMAEAjiAEGRraazWDAvJlbgW0BAJgCDWAAAAA0hRhgAAAANIUGMAAAAJpCAxgAAABNoQEMAACAptAABgAAQFNoAAMAAKApNIABAADQFBrAAAAAaAoNYAAAADSFBjAAAACaQgMYAAAATaEBDAAAgKbQAAYAAEBTaAADAACgKTSAAQAA0BQawAAAAGgKDWAAAAA0hQYwAAAAmkIDGAAAAE2hAQwAAICm/NTUFTiNlNK+pIuSDiUdSLqRcz6etFIAAACYtZRznroO1VJKt3LOF7r/70t6Ped8ZeJqAQAAYMa2NgSia/A+lHM+kvTqRNUBAADAltjaBrAWIQ/H64XrDWMAAABg1TY3gM/1lN2TtDdyPQAAALBFtvoluE2klC5LuixJZ86ceens2bOPfebBgwe93/3JT35izteKnU4pTVIuSU88Mc3fM6U6eb/jnVfN9vAuwxsnX1Mnb3nN9os81iI+XzMva3v8zM/8TG/5T//0T7uXbe3vw8ND97zwqYODg97ymuPmX//1X3vLP/nkk95y61oe+Q6Md17W573lNd+xtoe3/KR69Ym6/kbWqWabR3w+Us32GLq85l545syZ3vL79+//r5zzz5ozrLDNDeC+p73ntBYWkXO+IemGJJ07dy6fP3/+sRl98MEHvQv48Y9/bC7cajT/1E/1b1LrBm3tbKvcmo8kffazn3XVyeI9iUsNb+/6WXW1yksNnM985jOueVmsfe1driQ99dRTveVPPvlkb7m1T631trZr6Tvecmv7RTWkS8t45plnesu/8pWv9JZ/+ctfdi/73/7t33rLrXNvm18kHoK1v//8z/+8t9x7PkrS97///d7y733ve73lH374YW+5ta9LvI0AaxlWI95bLtnXKOs7H3/8cW/5Rx991Fteuhf+y7/8izmtT+ka1cfafqXlWuttzcsqt7ZrzUMzi/f6aNXV2qeS/Yehd/2s+Vh1Ku1r61747LPP9pbfvHnz/zNnVmmbG8CH6gmD6F6G65VS6r3YWjdVq/GxnJfnO9bN09uYiGzwef9q8z59k+z6WnW15mVt19I+evrpp13Ltk5WaztZ62Y1ZiVpb2+vt9w6Bq2LgbUO1nwk/42nVdb5Yt3cIp98e9U8eRlaTSPAy/rDpuYPnm1R+kPcauBbjUTroY81n+PjY3PZVqPZaoRax6a1ftY6/OhHPzLrZDUGrXKrYWct29sDIfl7IaztFPkHgfePNmv9rO1aejjgfdg1hK2NAV5v6HYvv701UXUAAACwJbb5CbAkXUopXVU3EAY5gAEAAHCSrW4Ad0+B3+h+vX3S5z/zmc/o53/+5/vm0/v5Ute+9Z3nnnuut9wKH7C68K3ujFIXt8WKtbG6M6w6WV041jpLdjfO888/b36nz7lzfUk/yvHHX/jCF1x1KoUuAJahXwKqWUZkaAQxzvNQCmEqXYP7fP7znz9tdcJZoRTW8feDH/zAnJd1r7p3756rTvfv33fV6b333jPnZd2rrPABq/1Rise2WKEtViidVSdru1rrXbreWO2PMcO3tjYEAgAAAKhBAxgAAABNoQEMAACApmx1DLDXM888o1/7tV97rPyHP/xh7+dffPFFc17/+I//2Fv+zW9+s7fcilf62te+1lv+d3/3d73lfTHMS1ZuzF/4hV/oLX/33Xd7y7/xjW+EfF6S/vmf/7m3/Gd/tj+ftRVjNNUgH8CYho63nWPaNEDyv4Px9a9/3b0M6zvW/cV7//rrv/5r97Kt71if/5u/+Zve8lIKwL/927/tLfe2P6z3ar773e/2ln/pS18y63T37t3e8hdeeMH8TjRaFQAAAGgKDWAAAAA0JbWU4uarX/1qvnr16mPlVrqRUpeM1WVojfTlHfrXSnVSk5rNWg8r1Zo1wpk1Oow14lqJlbbHCoE4e/asexlWN5XFSjFnbafSkMdox1tv9Y+/89prr41ck3G9+eabveWvvvrqyDXBHHmvm1aqLosVnlDy/vvv95ZbIRCltGZ9SqPTWfd6a5Q2aztZ6eJqUo5ZKVK9Q0NbIwWW2pfWelh1/b3f+707OeeXzRlW4AkwAAAAmkIDGAAAAE1pKgvEgwcPeh/VW2EFVneJZIc0WF0gVliB1aVgdR9ZIQKS3QVilVvzsupUk6HB2k7erA7WfEqjI1kj2ln726qT1SVjdROVuqK8dbK6zWre5i8dO32s9bPKrWO/1C1odZ1Zx7/VXWiNLGhlQHnppZfMOn3ve9/rLf/t3/7t3nLv6FK7wgrx+N3f/d3e8j/90z/tLf/KV77SW37nzh1z2dab8NbIXdb11+pmtkLZJDvkyyq3zmGr3FK6ZnqvB9Y5bF3rrJHBSsu21s+61lnXp5o6Wfcwq9wSdY+U7P3nvUda27sUbmDNyzr+rTpZ118rZLO0vb37ewg8AQYAAEBTaAADAACgKU2FQOScex/hW4/iS10KVteS1XXrfQPU6ja2Pl9i1dXqnrC6IKxll7aTt7vGWkbNPvJm3rC6Ba23Va1un5p9ZHVReUNkrHLJrq+1fh988EFvufVmdCnUAe2wQkL6BiHaJVYIhJXx5tlnn+0tt7L2WG/sS3ZXtreL2xsmJdnhJVZ9rfXzXpdL3evWPc+6zlrr4A0rLIWgWcuwQkKsfVSTuctatnXMejNBWW2GmjaUN0zlNHgCDAAAgKbQAAYAAEBTaAADAACgKU3FAEv9MSlWnEpNCg9vqjArlYs1Kk4p9YvFG2tjLcObmqq0DOs71jKsOlkxSZI/NZBVJysW1kq5ZMXalpZRE3sHYD5q0gDOjTdGVrJjfa3UhFbss7WMmpSPVjyxlYbSipW26mTF55a2k7UMq06RMcDWMqz7jvfz1n24JkXqmKMT8wQYAAAATaEBDAAAgKY0FwLRpyalibcrwDsijzdEoCQqDZo3NKJmXlZKLqsb0fq8ZK+f1cVi1dUarQwAdol1Xytd461pVoiYNcqelQ7Mun+VRgG1Uq1Zab+sz3tH7rTCHCR/OyMyRMCal7f9YS3bOm5K+6jUvhoLT4ABAADQFBrAAAAAaEpzIRB9j91LbypaakYm86gJy/C+LWt1T1jl1naqqZM3/MLqZiuFQFjdON5lAwBiWFl1vCOilTIueLNZWMv23lNrsh5Y5dYyarJiDB1u4F03qa49EY0nwAAAAGgKDWAAAAA0pbkQiL7H7t4uf8kOgbC6163H+lbCbm9XTalO1jIs3pCQmi4Lq65W2IKVNaK0bt4sENZxUBoQBfiVX/mV3vKXXnqpt/y73/2uOa9bt26F1MnrwoUL5rRvfvObveV37tzpLf+zP/uzkDphN3nvt1Z56dpvTfMORmGVW1kjrPLStNI9vU/NoBNWe8LaTt57urVPS/vIm51qCDwBBgAAQFNoAAMAAKApNIABAADQlKZigFNKvbEqVnxMKaWZFcPqjZ2xPm/FHtXUyYrDsWKGrDgpbxxRaVopPUofa71r4p688dsffvhhb3kpBRva8corr/SW//7v/35v+Ve/+lX3Mqxj0Bpt6/nnn+8tf+aZZ9zL/vu///ve8j/6oz/qLScGGJId82odg97YYOseWVrGs88+6yp/6qmnesutdbM+L9n1tdbbO9ps6V5oLdu6R3rfA6pJVUcMMAAAADAyGsAAAABoSnMhEH1d8tZj+pr0XtaIZRZvKECpTlbXgRWGYJVb87HStZTSuFjhFNZ3vHUtLdtKnWaxwiyee+653vIf/vCHveVWtzS22xe+8IXe8m984xu95TWhDharS7cmpMHLWg9rva3t9IMf/CCsTpgPK9zmhRde6C237gne+21kCIT1eW8IRKlOVjvDagNYoQ5Weakt4Q1r9LLmb5VL9v3WWr8h8AQYAAAATaEBDAAAgKY0FQIh+UaC846gJtndONbjfuvz3u4SyQ5dqHlDs4/VNVHq5vCOaGeVe0fwkeLCUaztevbs2d5yq5tNko6Pj3vLvXW1jidrJL3Sd/Aoq0v3t37rt3rLf/VXf3XI6syWtd7WdvqTP/mT3nJChjZTChGICuOzrrN7e3vmdz73uc/1lluhY1EjbpbuO1ZIw9NPP91bboU6WPcX6/OlOnnbE9a13HufL9XL+o63vVITHlmTzSIaT4ABAADQFBrAAAAAaEpzIRB9IpM4W6wuBW8XRKkbxdul4E2CbSXmLnVzWNOsDA3eUAfrbVypHA7Qx9tFde7cud5yKzRCsrt7vRkrajKUWMeHd2AQa9nWcVPartZ3vMeg9Xnr+Pj85z9v1umXf/mXe8t/8zd/s7f8l37pl8x57TJrvX/nd36nt9za5n/5l3/ZW/5P//RP5rKtQWi8x4f1eatc8p9H3mt/zXlk8YY9WddZKyxIsrM9lK7Nfbxv/5fuz97sDd6sDjVZD7zHoLW/vZ+X7PuwdzAK6/PWepfqNIeQPJ4AAwAAoCk0gAEAANAUGsAAAABoSnMxwH3xM96UH5I/btgb0+uN25XsmNeoVGvWOtTEAHvTnVkxWlasl+TfHt5R6Kz5l9KgWfHBP/7xj3vLvfGDpVQ43vQy3ljEmlR13tg7byy4lYrp537u58w6/eIv/mJv+Te/+U3zO/jUr//6r/eWf/GLX+wt/6u/+qve8n/4h38wl/Hee+/1lluxwVaMvVVuzaf0HSuVoTdm2FL6vBWz6Y3ltK6n1vsOkp0izbq/eEc4s66BpXthVOyut7zmvSHvCK4128MbT2ztC2sdrPtOqQ1lrYf3vDgNngADAACgKTSAAQAA0JSmQiByzr3d1lZXdqlLwZu+yWJ1TVjLLtXJG35Rs4w+pfQ13lCHTz75pLfc2kelEdS8Iwt5R62zumpKadCsEYo++OCD3nJr23rrWprmPQ6sfWd1n1rrLPlH07M+b5VbXbdf/vKXzTp50zdhM1YIyde+9rXe8u9///vmvO7du9db/v7777vKrfPO+rwkffjhh73lVhiTdY2yrmlWuXVtLE2z5mXdd6zzrhQCYY2uZt0LvXW11Ix8FjUyqVVek6rOWm9vWEFNuklr2d4UZd7wjtKyvcfBafAEGAAAAE2hAQwAAICmEAIhu0um5q1KL293Rk1mCm8GCqt7x5p/qZvDmpfVzeF9Q7gUfmHVyxsSYq2D1a1UCsuwwgSst9qtY9Pad6Xu+6hRkKzPW6EOpZAQK0uD9Wa59XlrGaWQEMyDdTy9+OKL5ne+8pWv9JZboQvW+XV8fOz6fGkZVmiENzOF9flSZgprmnV9tM4L7/ko2eEGVoiYNwTCm3lA8mf0iQpFLNUpamTNmhAI657nHX3PYq1baf5RoTCnwRNgAAAANIUGMAAAAJrSXAhE3yN5q8u61KVgde94uy2s+XgHQJD8WSC8g3l4k2mX5mV1jVhdKdYyapLDW6x9532juLRcKwTCepPa6ia1jqdSCIQVomCVW3X1ZnsoDQzizfZQGvgE7bC68D/3uc/1lnvPu1LYjpU5wpsdwiq35mOVS3YIhHXtt8KbrBCIUiYX6z7izdzjHQChZtAfq9yal3egqlIWKO+AKNY13ltes2xvu8T6fCmcwTtwzBB4AgwAAICm0AAGAABAU2gAAwAAoCmziAFOKR1IOt/9+oqkb+Wcj7tp+5IuSjqUdCDpxibT+uSce+Mza0beqUnN4plPTZqVUgxQxOe98VOSP1baOwKN9/Ol71jbwzuaXWm7WrGIVkyeFSdobVerTpI/RteKj7TWwYo/LsXtWt+JSjMISP60gaXrinXtsq4T1rltnaveVF2SvR5WXa1z0op9Ll1XrG1lxXh631Xx3qdqvuO9n0fdv0rL9paXlm3F1XrTkVqs915KbShr2d53d05j8rtMSmlP0ss55ze63y9K+o6kl7qPXM85X+imHUm6JunKBtMAAACAx8whBOJlSa+v/H5b0kFKaa97wvtQzvlI0qvSw6e/vdMAAAAAy+RPgHPOt1NKl1aK9rvy45TSeUnH69/pGr8H1rSuMdyrr+vYeuReSuERFepgdTVYXQelNCveriXv6C3Wsktdc0OHitR0j3lFjRIk+UfZs+ZlHR81afK8Xa7WMqx9UerS8qbCKR3/nmWXugutadayre1qdRtbXeJzZaXis/addd2MSvck+a8f1rKt46N0zHqvaVGjj5XCELyjenpTH5a2d9T+jmSdw94QBW+4QSl0wJv601vXUnvFul9406N517tUJ2/K0yHM4Qmwcs6HK7++JumN7v/nej5+T9LeCdMeSildTim9k1J6x4qnBAAAQDsmfwK8qosHPljG9UbIOd+QdEOSvvjFL/oj1AEAALBTZtUA1uIlttVwiMee6Grx5Pf4hGmmvkf43u6x0jTvSHAWq+ugFG7g7R72dnPUdF156xQ5+o23vt5RbqxjoNRdaO0/q1s8agQfyT8yk7eb2aqrNX8pbuTEqJEFS7z7zsqiYZVLdhd0KazGwzoGSr1jP/rRj1zlVsiE983y0rUjauQu7zkv+bt7o0JISl3i3mPTKrfmU3MvjAx56ROZAchbXjPymTcMwXvvLJ1f3owL3rZBTSiRdxlDmE0DOKV0VdLrXezvXpfO7FA9oQ4556PuZOmdNnRdAQAAsL1mEQPcpT57eyWH73np8cZs9/LbWydNAwAAACyTPwHuGq43u/8vi48kvd39/1L3dPhQi/jg1Ty/pWkbq+kGiwof8Hbplro5vG/ne7txaraTt0vS+9ZtTQiEtxvH4u12l4ZPAl/qgvMmPbe67KzuU2+WiRrerraaEAjvQAfeUBhrPifVK0JN163VhW+FQHz44Ye95dbx5B0sR7K3rbe8RlQIhDc0osQ7KI+1zb1hGaXveK+n3vtXTSYX733be38uDfxgfSdq8KeaDBTe6413O5XmP2aog2XyBnD3JNe8y3fTl1khbm86DQAAAOgzixAIAAAAYCyTPwEeU845rIvR+/je+0anpeatSm/3nzcxd6mb3vqOt7vLUhMCEZkUv0+pC9PqqrS2h7cLv6YLzuqKtTIDPPnkk73l1rrVZC6JysgRme3BO5hN5CAcUWoGCPAOmOMNBfBmh5Ds48Y6j6KyjUj+bCcfffRRb7l3AInSeWSFH1n7ztoX1udL17TSNaePN0zFe404aZrn81HHvlQX6uOZj/c+FbnsyOtWcwNhAAAAAGOhAQwAAICm0AAGAABAU5qKAZ5S1AhxNTGv3jifyNFvalKnRfHGO3pTz9XEXFspxKLSGJW2qxWrZ8WOWTGH3vi3Utovb3x6VHqjmvMoatSkbYsBjorH9sZTlurkvdZ5r7OleErrPPKmO7PW2xsLK/mvzdY6WLG+1jpI/lFRreuKN/a+ZmRSL+8xW7r+eo9Zq9w7wl5pmnWszSFF2Rh4AgwAAICm0AAGAABAUwiBUF3ogPc73tQv3i6Lmjp508vUjE4XNeqapWakGW+dvCOolY4ba15WuRUaYakZ8Wrorv1Sd5oVHlGT+mgq3vUuHbPefeHtJq1JF1ezX6fi7R6uSSfoDXXwjqwZOZpdVKhDaXt407ZZKdus62lNCERU2tHIlGPW/vOOlFnTNrCus957oTfEqKYNNSaeAAMAAKApNIABAADQFEIgVNfl5H1T1/q81R1U8zZz1Nuy1udrwhmiulJqupOjuqa9b/yWQkK8XbTebjPreJLscApvyIv3Lf+a0Qu9YTtepXMiajSxyGPW2+XqHYWupk5RYWBjjC4VlcmiNM27j6xz1TpPS+e2d5RCbwhE6Zpm7W/vCH/W9osMgYg6V73hK6VlR92fS6EU3nPPEhXiedK0sfAEGAAAAE2hAQwAAICmNBUCkVLqfexe83amNS0qNCIyQbU3DMHi7SovfcfbRVXzxmhUF6r3TfsSb9iEte+83aeS/1iz3vyuScTuVfOms2c+pWM/6ty21GSB8IbhRGaBsHivad633WvCMixRWWck+5jyhrNZb+Y/+eSTveWl7efd3zUDOXhF3V9qBrOJymDjvdbVhGhFDRZVk3HBG87jPbdL6+AdCGYIPAEGAABAU2gAAwAAoCk0gAEAANAUYoBlx6/UxADXpC+L4o3v86aEsZQ+761TVAqqk6Z5Pu8tj4w39O47K66wNM0q945sZSnFpnnTnXljg2uOG28MsPc4q4mL826nqJR+JVGp+yyl64o35tW7/UrHrPc8smKDP/vZz/aWe7eTZKc1i4ojLW2PyHMvSlTau5p0ZxZr2d6RAr3v1dTwxjjXxAB7r7ND4AkwAAAAmkIDGAAAAE1pKgRC6n+EXzOKibd7p2bEIS9v6hdvN4c37VFNnaLqetI0D29dS12YQ3dflfZFKTyij5VSzepu9Y78JMWFOnjDEyK7dC0153xUeIm3u7VmtD5vnbypHUvL9Xaje88777ki2ee9dR55l2GlJazhPV8iRyCNCoEoHR9D33dqrhHeUR5r0r95RYVN1rQNoq6zp8ETYAAAADSFBjAAAACa0lQIhJUFIrKrxhvSEPmmp3deUeEG1lvOJdZ28m7zmq7sqO7TqM+XWOtXswxrP3m7bq2312tCI7zH4NCZPUqiRl0rhTl493fUCE+lOnm7Sb0hEzX7KOq4sY790jXNe754R4izzpdSCIT3flGTacISdU5aaq513i58b2iVd59KcSO+WWpCQrxZYbztm5rQGbJAAAAAAAOhAQwAAICmNBUCIfV3UdR0h0Z1SXq7d0pdV0O/dRuZ2NybLaNGTQaAufF2B5WOJ2uaNzTCOpa9WSMku1s3KuwkcoAH67jxhhXUdAtGDWZT073uzRzh3Rc1oRRR182aAWW8WR2sc9jaF96BEUrG7E4eyhj3BO/gLd7jSYoL+7PUDLDjDW+y1Bybc7g/b09LAAAAAAhAAxgAAABNoQEMAACApjQXA9wXd+KNoytNGzplVik+xhuXFFXXUgyTVaeoVGE1qV+iYqVrYqK9y/aOyhc5OpI3RZQVA1lKKfXxxx/3lntTqkXFyJ40zWOM9IpRMcClNGjWNKvcWoZ339XsI2+6M2+qv9J3vCnsrO1XExPtTXfmrWvNvoiKI41Ktyf5Rywb+p5ampd3GaUUZd60iFHXwNJ8hn5naRM8AQYAAEBTaAADAACgKU2FQHhHgiuJ6l6P7PaxpnlHRPPWqWb7ebvwa7aHt9vY2k7eUIdSmIo3pMFbXuLtBrO2rbcLuGZ7WNvWCpmw0nh5Ry4q1cl7DNZ0J3tHbIrqJi2FQFjb1hsCEZVySfKPSmmFNDz55JO95dYxLtmhPtZ6WOE83vOxxHuNj7zOWtOizqOa48Y7Ml5USEhk2I4lch9FhUFGhn4wEhwAAAAwEBrAAAAAaEpTIRBSXBaIodW86RmVxWCMEAiLd16lfeQNgfCGNNSEQFjdO1HhF6XuI+/b6FbXt7UMq2u41J3sDeXwdpN6u+NLvG+Q1yzb2hdRYUze7ST5QyCiQh1Kx7I3E4kV6mCVlzKXRI2y5x3xraZr2BuGUHNeeLvFvcuOzFQQmVnBW6eoUIea0Jmhs1bVzJ8sEAAAAMDIaAADAACgKU2FQFhZICyRb3R6y2vejPYa4w3QqKwYkQMaRL0B7e2WLk3zlntDJiR/CERUF1WpO9k7QIG3+887iINUDgfoY23zmmVb22PKEAhvtgfvtcs74IpkhzpY2R6s48ladmkdrO1hZSixyr0hEJGZX6KuQ6VpUdfZMUIdojK51GwnizcjTCkEImr9IrNAWMgCAQAAAAyEBjAAAACaQgMYAAAATWkqBljyxVHVxLV447ci0355Y4ai4gojY7SiRuqpEZXOpyZWLyoFW4k3LtQawcq7Pay4TMmO87SOTSuW0zrGa+Jwa+Jk+1j7zoo7LU3zHgfefVRat6Fjfa14Xqtcso8Dq9w6nmpGxvPG+lpp0MaIeY26btbcd6LuF5HXU+810FvXmhhgb50ijxuL9/oRqWZ/Vy9rtCUBAAAAM0ADGAAAAE1pLgQiapSRmjQoHjXpaKK6LcbsgthUzahrUbxdVDWjsY0RZuGtkzc0wpveSLLDI7yhANZ8ItN+ectr0ll5Q4CiQiBK+8g7AqQV1uJNeWeN0ibZ+9vaTt7jwApnKE2zjoMxuo0tQ6e/kuLCB8a4pg09uusYaTkj2x5R2yOyTowEBwAAAAyEBjAAAACa0lQIhDUS3BgZBqKyHtSI6nqJHPnMu9412RCGDkex1HSDRWUPmTIkpKY72ZqXd+Qu69i0uspruvytcqsL2Hr73wohmStv9gZr3z311FO95U8//XRveSl7iLW/vaMdWvuitI+s4zxyBDcvb6iUN2xhjMwUXjXXWa/IDBelkdo88xojrGCMYzlqpNHT4AkwAAAAmkIDGAAAAE1pLgSir3vOeqw/5du7Y/C+3T3GW+1jhEB4u1ii6lpatndeVhdwKQPF0BlQaroFre5kb12t9bYyEpQyDERl6vCGRpSmebusre1nbafSoBPe8BIrpOGZZ57pLbf2hbXvSqzjydqu3nAGyX9ODn3eSf7jIzK0ygqR8d5fLDXXlaGzPdRkU/GuxxihAGOEYHqXPWa7a3ZPgFNK19d+308pXU0pne/+3dtkGgAAANBnVk+AU0oHki5LurJSfD3nfKGbfiTp2sr00jQAAADgMbNqAEval3S8/CWltL86Med8lFJ6VdKV0jRr5k888URv91yp29ji7QKJerPc+yapFPdmb02XlrerLaq8NG3osIzS8eQNXbDKvdkQStOiuiQju668oRFRoQCSncXAe05a61DKMODNYuBd75ouT29WB6vcmo+1L2q6k61t7j2easIvasK0+ngHpiktw3uu1nTHW9OGvkfWiAqLGyMEwlITGuG9vwx9jyypaePUmk0IRErpYs757bXiA600iFc+u3/CNAAAAKDXLBrAXaP1qGfSuZ6ye5L2Tpi2Ou/LKaV3UkrvfPjhh6esKQAAALbdLBrAkg5yzodDzDjnfCPn/HLO+WXrLWQAAAC0Y/IY4JTSeUm3jcmPPdHV4snv8QnTep05c0bPPvvsY+WlFEAWbwyalYbHiu2zPl+K/4lKKTVGjGeUmpRjXpHxlN644THij71xcVFpwiR/vJd1vljlNbwxed6R0kop2Lxx9t4RwGrqZH3HWm+LdX0aM+ZvU6XzKOq9hshR17zpzrxx15GjW0YpzT/qnRvvMVsTAxwVl1w6H73XKCsG3nstqLkHl9JERpu8Adx5dWVD7aWULmvRKD5UT6hD98KbrGkD1hMAAABbbvIGcM75kae/KaXrOecbK7+vTtuX9Fb3vSNrGgAAAGCZvAG81A1icbn7/zUtcvweSbqUUrqqxdPgg5zzapqz0rTHnDlzpnc0Iqv7r9QNZnV1WF2xVuqjjz76qLfc6roqdfVa9fV270R2cXtT+njTe5XCDaLCCiK7+KwuIWt7RKXnkYYfqcpSM4LV0OUlUd2nllKKLWsUNWvUtaj0bzVpv6xrmnVdiUyVFJVO0CovHbNTnS+l4yzqvKg5xr1hAt5l19QpKtQh8trvDZ3x3gtLIRDeETGtciuVoTX/0rlt7QurTTSE2TSAc87Hkt7oflbLj1bKbm86DQAAAOgzlywQAAAAwChm8wR4DCml3kf7Vnq0UpeCd8Qmq7vQ6jqwykujSFnL9o6OFNWlJfnDJrzrUNpHVr28oRGRXbpRXbQ1y45629hS0y1oHQdRx2xN5gFvJgZveWkfWcdzTYhCn5owph//+MeueXnLazKaWNspqrwmLGPo82iMrAc113jvsr3HYE3GhaiQBu85XJOFJ+qYLV0jrNCFoUMgarJWjYknwAAAAGgKDWAAAAA0pbkQiL5H9dZj/ZoQCO+bnhZrPqVuDiuBtLd72ApD8L5BW/rOlMnvo7puxxgYJGoQk5rvjJHgPmqbRyWyLy3D4n2LuyZhvTeDhzeTQOkc9l4/vMdHZAiEN2Ro6LCgmmXULHuM68TQIutUOnYi5lOTicGbgcV7/SgN6GVNs9o+VtYZ78AZpePS2h6lMM9oPAEGAABAU2gAAwAAoCk0gAEAANCU5mKA++JOalLeDB03WRMX543DsZZhxRh5092UplmxiJFxYN7RsKLityO3R9SoSTXfiTqWS9svKoWY91guxbx6U7N5U/qVRnMcepQxS+Soht74SO++Ky3by3vsnzStz9DHeM13Ikfc9N4/o+L753iviNxH3vOo1DaI2kdRaeRKdYqK394ET4ABAADQFBrAAAAAaEpTIRBS/yP5mnRWVneolYrs448/7i23Un5Y8yl1n0alJbJ407hIcSmAakbV8nYhedPRRIbORKXxqknvFXV8jDEynjdMxVq3UgiE9xz+6KOPesutEdSs8tKyvee2NzVbKYXSU0895Sq3RpHyplaqCYGI6l6PHPnMK/JaF9WNXhNy4g0l8o4CGnmtGyM1ZtQxaG2n0vni3YZRYT6lOnlH3BwCT4ABAADQFLN5nlI6K+lcxTzv5Zzfr68SAAAAMJxSCMS3Jf2FJO/z6Je7785Ozrm3K8DqHig9ivd2k3q7Q70hE5I/Y0DUCFbWaDKSv3vTG55QCr/wjnLjLbfmX6qT9w3XMUZpGzrjQmmda946j1A6j6yQhvff7/+7/r333ustv3//fm/58fGxuWxrGVZ9rXPe2q7WMXv27FmzTnt7e73lzz//fG/5c88951qGFTJRuq5E8XZLl6Z5s4REhiQNPdpcaXtYXfLWMWvd27zlpXPYm4HFG5bhXefSMrzhTTXZE7z3wqiRWmuyqZTC06KVGsDv5Jz/h3eGKaWap8YAAADAKMzHLDnn/77JDFJK//va9/74tJUCAAAAhuLKApFS+paka5KWz66TpOck/dfgeg3iJz/5iX70ox9t/PlSN6zVBeJ9I9yqjzc0QrLDMqJCIKyufasbpTTN2yVT85apN7m5d/2efvrp3nLr7XjMR6l73Xtslo5BL+uaY10PvCEQ1rFZCoGwQh1eeOGF3vJz5/o7AUvXianUZC4phTi1yDo2reuvN7uBN+OHZHeje8MyosI4JH/4hTcEopTJJSrUITJ0xmqXeNpop+W9cu/lnB+5unWNYgAAAGAreN80ud1T9mZERQAAAIAxeJ8A308p/UdJRytlVyT9blyVAAAAgOF4G8D/WdJ5SfdWyl7SljSAHzx40JuyyBtHV/pOVIyRdz6SHQPsjfPxprkqpX6x6mvFWXrjL60USqV6RY3K5913kj/l05ij4rTEitWXpA8++KC33EpfZqU7s8qttGmS9OGHH/aWW/X1jgRnHX+luDtvCkdrHax0as8++2xveencxulYx4117SqdL973WKzPW+XeERhL06Luw957heRPk2ex2iWl+Xjfh/Gmp7SWXdOGsq4fQ/A2gO/mnP/LakFK6T8F1gcAAAAYlDcG+Kin7G5ERQAAAIAxeJ8Afz2l9N8k3el+T5IuSvrN0FoN5MGDB73dm9aj+FK6G6sLyUp1YnWNeFOjlLo5vCMUebtPvSPZSP7wASuVize9nDR8mIVVHpkWzjoGrXCUUvqmmpRPHtbxVDpmvSMwWceT9/iwwhxK06zQBW95qYvPqq83xaE3VVIpdd8zzzzTW26N+OYtt0IgrHLJrq91Tlrr7T2/JHvbRoUreUeOK03znl81oXfeEIWhwxNK06zzyFtecy/0jvznPZ7GuM562waldbCWXbqnR/Pe+f6zpPckPd/97EnqTwYJAAAAzJD3CfC3cs7/92pBSqkvNRoAAAAwS2YDOKX073LO/89q2Xrjt68spfTvc87/M6yGgayR4KwuiJoRf6xuCG9Xire7pIb37VMrxKKUZcLbVWRtJ+/odKXvWN2hQ5dLdn2jQh1qwhmsY8ob0uDtZpP8ITJW96n1BrkVbhAZAmF93iovZVzwdgN7ecN/JHvEQ2/oQmQIhBWWYdXVm32l5hz2nqtW93BNKIU3bMJ7rtacw0OXl8INvKGF3nLrnue9p5aUMij0KbUNvMdBKZtFn9Lx4f3OmCEQpSfA/1tKqe+lt5Ik6T9ImmUDGAAAACg1gN+U9JpzflmMDAcAAIAZMxvAXWjDYyEP2yzn3Pto3+qiKnWxeBNCe8MHvN3SkbzLKIVAeLeH1S3iTcxd+o5VbnVhRn2+9B1vSIPVTVpzfHiPTW83ac1AKVGDL1jl77//vlknK3TB+o532aUuPmtbRZ331nFjhQJIdsYFKwzBKj979mxvuRXqYH2+ZtnWOljlpZAQb8YWq9x7fSrxXg+iQiYkf/ahqM9H3ne82yPyPhyVPaRUJ2taVOYob3uotAxv+MVpxOQ/AgAAALYEDWAAAAA0hQYwAAAAmuLNA7zVcs69cSfeUZYkf8xLVCzWGMaIP7Zij7ypgWrip6LibWtGhIqK97LUjBblTRnkTVdkpS4rTfPG+kamQfN+pxQfOTfWOVwz0te9e/d6y604dyvW14rbjUyD5o0NttKmlaZ50yJ6UzvWjPIYpSa+dMp7ofceNsa7NZH3MK+omG/rvh0ZAzzm9ZQnwAAAAGiKqwGcUvp3K///ekrpP6WU/n14rQAAAICBeEMg/mdK6f/MOf8fOed3Jb2bUvq/tCUDX1ghEDXdR95ui6G7iUp18tbVm1KnpntsjC6nKcNILFGjrnnT/JSmedOaRY3eVprmHfHNCpmwPl8zGhs2Y3Vh3r9/v7fcGwYj2WEqVkiDd4S4mhAI76hyUWnTStOiwr3G6Kb3GjrsQ4oNyYu6D1vlpfud9/7iDVH0pnOVtjME4pKkJ1JKf5FSWgZo3QiuEwAAADAY7xPg53LO/yWldF7SYUrpdUn9f9YDAAAAM+RtAF+Q9D9yzrdTSq9IekvS85JeCa/ZQPoe+Ze6jS3e7omoEIga3iwGVp0iRx/zhifUhFJ4t21UGEJNnbwjFHlH8ClNs8qtkIao0AjJDjeIKrfqWnPOYxje0BzJ3t/ea5R17JeWbR3P3lAHK2TCyg5hlZemeUMjakbc9N5fosIvaurkDUPwjspXEwJRk03Iqyaksk9keKT3njcEbwjE7ZTSf0wpnc05H+ec/4MIgQAAAMAWcT0Bzjn/956yP46rDgAAADCspgbCkOJCIKLeQB0j1GHoEIhSV421fkMPJFKaFpVZoaYLx5t4PCpReek73m5gbwhEqU7WfrWOKaur1+pOto4zq1ta8g8M4t0XNW9rD610Dnu7y71d+N5sCKVp1jK8XbelY9Y6zi3e7uGaEAjvvrD2qffzpWlWnaIyVpREDWpkGWNgpjEyckSFRlhK89nGEAgAAABgq9EABgAAQFNoAAMAAKApTcUA55x740u8Ma/LeXm/4+GNnZXiYnetGC0rZqcmRisq/coY8VDeuMyaNGhR6dFq4ku9x3Jk/KB1rFkppax425pYTktUmrxWjRHTaB033thxb5yqFJfeK/I4847QZam5F1qGHgV0jFHXova15D9mI9PCefdF1D0v8l44BJ4AAwAAoCk0gAEAANCUpkIgpP5H8jWP3KO6+cZIexSVss1S6maLWnZN+MUY3Tjez9eEtkSJSodnpaCq2X5Dd4fWdLtHpWnyhhiV5jW0mnRFUSMneudfmjZl93pUN3pNOrChU47VpEGbqlzyn3tTntuR6d+8okZ9jUxTOmYqSJ4AAwAAoCk0gAEAANAUQiBU9+ar9R2r22KqEZ6k2LdJPfOR4jJTWOU1b0ZH1anmuIkKQxgjC8Qcu5O9b/PXjKrlHbHMu4zIblKvmuPGClGIGlnQO5/SNG9ddyVsZ+iwgtIx6w0rGGOEuKnCL2r20ZQZjqI+XxPOMIcQiNk0gFNKFyWdk3QkSTnn2135vqSLkg4lHUi6kXM+PmkaAAAA0GcWDeCu8bufc36ja9TelPRSN/l6zvlC97kjSdckXdlgGgAAAPCYWTSAJV3LOb8oSTnnI3WN364x/FDO+Sil9KqkK6VppQX1PV6veYPR22U9xhudUXXyDoxQE4YwRpe/9zveLrsx6uTtJhqj+yhqYJXSNG+5NXCGFbZgfb7mO9bnvSETkr8b2OI9PkqZXLwhDVb5J598EvL5mu9Y6+ctL03bpnNy17NADB3SUBM+GBXqMFWmGMl/P4oMyRvC5C/BpZQOJB2v/H/Vw2lr39k/YRoAAADQa/IGsKR9Sfe6MIijlNLVlNL5btq5ns/fk7R3wjQAAACg1xxCIM5JOr8Sy3tD0ruSno+YeUrpsqTLknT27FlXCETNG4zern1LZDeHt06Rbyd7t623a6m0j6K6KnchE0ONqCTwpRAIbwYFq9wKT/CW13xnjBCIqDfC5xgC4S2v+Y4300QpA4U3bKKme3hoc8xAERlaFRW6ELmdvKYcCMPiPZanHMxmE3N4AnxPiywOkqQui8NeF8rQ90T3nBahD6VpD+Wcb+ScX845v/zUU08FVhsAAADbaA5PgA97yo71acP4sVCH7oU3WdOC6wcAAIAdMvkT4K7Bei+ltCc9fIntKOd8vN6Y7aa9tfK93mkAAACAZQ5PgCXpkqRvp5TuSnpR0m+sTkspXVU32EXO+cqG03pFxV1FjfjmTS02RmywV00aNIt3u5a2d9S8IuOeLLuQKikyBtgbV+str0mDFlVeigGe40hwUTHA3n1RkwYtqjwyBngXUhzWfGfouNrIlGNRsbtTpigriYr1jbpHlqaNGRs/iwZwF/f7ujHtSNIb3a+3N50GAAAA9JnnnysAAADAQGbxBHhq3nRgUlw3VWTalKmWUeritni3U023SFSYSmSXZFQ3Zs0+jeqStLrwrfnXpP2yyqPCDazP18zLO3pbad8Nfa5Gprmy1nuMrvqo1FvWPq1JC2ddo6zPR3YnW6a83kSFG0SGLXhDFMa4P3uNcX5FhUDUpJIdE0+AAQAA0BQawAAAAGhKcyEQfY/krW6OmlFMvNkbvCOilboNot5AnWNYRk2dxgizmEpk93pUd7I3FKD0HavcG35Rc054z+3IEcCGfku9pqty6OwGltK+qwm38cyn5pi1tpNVp8jtN8fR5rzGCE/Y5ZCGyBCZqMxHc82ItMQTYAAAADSFBjAAAACaQghEQeQAD3PMGjFGGEKUyGVb23zKJOZjvBntDZuICo2o6cr27osx3py3uri961ATpjJlCIR32w5dXqqTJSojjOTfR1GhDjXbIyrjzZTmGLZQwxvuGLkvhs7qMOUxexo8AQYAAEBTaAADAACgKTSAAQAA0JTmYoCjRMXzeOdT4k3BZpkyzngM21SnqPLStKjYYO9yS6JiemvOr6H3xRyPv5Kp4gcj45LHeNci6rwYY2TSOcRfbrMp43MjU9vtcuz4JngCDAAAgKbQAAYAAEBTCIFQXTept5sqKjSiZNu6Voc2x+0xdLqzmq79qM9bao5lK+VYzQiJnvlEf2eOy+gzRhfmGN2qc+zS9V7755huco5d3HOs05RqQoaizovIdJNj4gkwAAAAmkIDGAAAAE1pKgQi59z72L2m22fobA9zzA5hqXk7OaqbbxfCHErfmXL9orqoIt9aBqa2TW/IR41YWnONH9oYdYrc13PM9uD9/JRhTEPgCTAAAACaQgMYAAAATWkqBGJKQ4dGlL4ztMhQkaguu5p5tSqqW2uOb5CP0QU3h7eZ52SMbCNzPLc5bsY1xmAUYyw7ypRZU7YVT4ABAADQFBrAAAAAaAoNYAAAADSFGOBgVuqSM2fOuOYTmY7GO3qWNz1aZDxeZNq0OY605DXl9phjXNw2jSbmnf9cDZ2ycIwY4F0456ee11S2aXvUpEHzpjubMrXYrqW05AkwAAAAmkIDGAAAAE0hBKJgjBRbu54GzbJNXZKWOa5DZEq6KU01clHksqM+P7WocIPIcB5vOMUuh0NNbZuuK17blCZsjqOuzXE7reIJMAAAAJpCAxgAAABNaSoEIqU0uzCBMTIoWKKyQ0yZkWCM/Rm5jCm7XOfeHbVqjmEI27T9pjTGuc2+GEZkOMOUoRFThSWV5hOV7WGO2y/SmPdIngADAACgKTSAAQAA0JSmQiAiDd0lvwtvLUeaY7L8XQmNGNo2JbIv8e6jXT+H57geQw/2UprXHE15zA4dCrMrGW/G0Hq2BwtPgAEAANAUGsAAAABoCg1gAAAANIUY4C0zx5HgprTrscF41DbF4XIcDIdtO4xtivUlTd6jWlzn0+IJMAAAAJpCAxgAAABNIQRiYlOmndl11noTAvGoobvOptwehDqMj67Y0+F42kyrx9mupJWcA54AAwAAoCk0gAEAANAUQiAa4O3mmLLbOKp8jDp55zOGyJGttql7bIx9tE3HwRi8x8fQn6/9zlQYSfJRc7xPjXE8WcvYpmN5W/EEGAAAAE2hAQwAAICmEAJRMMeu7F3o6hrLHMMshhaZHH4X1jsy1GEXjo8xeLt0GehgM3M9zqaqV4vHQK2pttXcz2GeAAMAAKApNIABAADQFBrAAAAAaAoxwNjYT37yk97yM2fOjFyTeZhjTF4prmqb0p1FpRyridv1LvuJJ/qfI4yRai1KZMox6zrhPf6mTEE1pW06bkrmuG2HZh37mCeeAAMAAKApNIABAADQFEIgsDGrq3euSFu1e6ZMbWeF+oyRgm1okSEQ1jo8ePDA9flILXbHbxvvcTDHfVq6RxIeMT/b1aIBAAAATmkWT4BTSvuSzku6J2lf0o2c8/HKtIuSDiUdbDoNAAAA6DOLBrCkiznnN5a/pJSuSXq9+/V6zvlCV34k6ZqkKxtM6xXV3UY3+vwNPfLZHLvgSuZY36HPo5r5e0N9diEEIpJ1nFnb1Zs1orQMrzlmRml1ZLwps4HgdGrOozlc6+YSAnFh7fc96eET3odyzkeSXj1pGgAAAGCZSwNYKaVbKaW9lNJ5STe74gNJxz2f3T9hGgAAANBrLiEQlyR9R9K7kv5gJRziXM9n72nxhLg07aGU0mVJlyXp7NmzMbUdwRjdA3PogpiTOXa1zbGLdo626XyZ43kX2b0+xjHb4nlRs25z3E67vI9qcL5MZy4N4Je1iPndl3Q9paTVmODTyDnfkHRDkr70pS+1vbcBAAAwfQN4Gc6wbPCmlG5LupNSuqGeJ7paPPk9PmEaAAAA0GvyBrAW6c9uL3/JOR91jd9zWqQ3eyzUofuMrGlRFdumbtXSvObY5eoV2VVT85b1VHahi6qme33ojBw1+3rKfRF1bEYOeDH0W/ul+YyxjKGXPYY51jXquBnj+JjSGFmJhg6NmHvmkjm8BHdbi0bwI3LOR+uN2e5p8VvL6dY0AAAAwDL5E+Duae5RSumqpCMtnupeX/nIpW7aoRahElc2nAYAAAA8ZvIGsCTlnN8uTDuStHwh7vam0wAAAIA+s2gAj2noOE9rtCPvcqeM541ahxremNBtiw2eY5yxZZvqGskamcw7QpylZnSkoePlxogNtrZrq+YQA7mpMa6zY8xn6G1eumZ6R0KMXLbFe8+LXIc53F/mEAMMAAAAjIYGMAAAAJrSXAhEFG+YgLdLoSYEYqqwiSlDI2q+M0Z6GcscRxOL6navmc8cu/ytOnm7+WpCJqbqFqzZ3t7tEZk2bejUbFN2r+9KuEGUOYaKTDl6W6SokIahwzuGwhNgAAAANIUGMAAAAJrSXAiEp1uypgti6FCHmm7VoTNT1Gi1Sy1qm085ml1knYbO+hG5nbzzevDgQdiypzR0+MAYoRFRnx9jGbtSpxbVnNtnzpzpLfdmUykte6qQBmvdpHkcUzwBBgAAQFNoAAMAAKApTYVApJRcXRSRGRcisz1Yhl7GHEMg5tCN4hHV5b9NoQ5zzAJREpX9osY2ZYEYel7blolhjqEfrYY6zDHbTpSa8ATrO0MfN3M/nngCDAAAgKbQAAYAAEBTaAADAACgKU3FAEv9MT1jjLrmTV9Wk7psm2J9vbGfkXFxU45O443FmjK+dMr4Y6/IWMcpRwps0a5vvzHikoceGW+O18zI9Ipj3COHfqeiZnvUpFTzzD/ymB0CT4ABAADQFBrAAAAAaEpzIRB93SlThkBEfb5kF7qsa+YzZbedxapTzQh/Q5vjcTOGXe+Sx7jmmOrPsivXzKHToM0xNK3mO0OPuEkIBAAAADAjNIABAADQlOZCIDxZIGq6pc+cObPxckvLiOx6GfrtU96ob8uUWUW2ya4fy7u8XyMzDOz6cbALxsgOMcbImhbvMrzZIR48eOD6fGkZY+IJMAAAAJpCAxgAAABNaSoEIqXU+0jeCkOoecvUOy9v18uUb5/W8HbXRNbJ2uZz6HoZSmTmkjkOrDKlKY/lXcD2GM4uh1lMOShU5PynvJ56wxe9oRHe5ZaMeX/mCTAAAACaQgMYAAAATaEBDAAAgKY0FQMs+eJnSp+NivWNHPFtF0TGskWmt4sSFe81RlzcLsen1xxnUXXatnjNqa5Fc3xHQZpnurOo9x2sNJ5TrtsYx8EY18Ap7+neYzZqxNKaeN4xtxNPgAEAANAUGsAAAABoSnMhEH2P/CNHV5tjd7LXrnR3eVO8jGHK42aqZc8xnGeOdaoxx2N8m0SmlJojq8u6ZlTPuZkyDVqkXajTGPftIfAEGAAAAE2hAQwAAICmNBcCEWWbuoosUetQms+UXW1zfHs9qmuppotvqtEI5/pmNGAZIzuEt3yMEbJqMlxwL/xUZEjNNl0bt/UY4AkwAAAAmkIDGAAAAE1pKgQi59zbjRTZ7eNNID1l10FU19wYIRBz7GKJfNN4FzIuzLFOUw56guEMHQ4wx+wQUw4MMsftMcdll47LOV4fo8JwItsGY4T6LHF3AAAAQFNoAAMAAKApNIABAADQlKZigKX++JI5xuDMsU6RcT5Rn48Utc1L67BNqW2isD2wbeb4zsGUxtgeUy7D+x7QGCNxjiEqBtgqf/Dggbls77yGwBNgAAAANIUGMAAAAJrSVAiElQYtkjX/M2fO9JZP2S0SFZ5Q03U1ZjfHprwjpW1T+Io0fKq1MdK/eeczx+MMw5nymjZ0Wsk51qmmXnNIf7UpbxrFMUYBjeTd31ZIQ80+teZFCAQAAAAwEBrAAAAAaEpTIRBS/2N37+gmkt01YpV7R4ibo8guuClZXUtWl0xkV5Q3zCKy22zoZYwRGjHV5zG+Ma43U4YVRM1/jGVEdktvU2Yga73HuKZNydsmqsno4A2nGML2tL4AAACAADSAAQAA0JSmQiCsLBDet/xL3/GGRsyx6zayC2KbBrYYOhRAiutSGyMsw7vsOYY6DD0fDGfK69AYoRGWbQs3GDr0I+rzkSKvH3O8p0eFOtSE7RACAQAAAAyEBjAAAACaQgMYAAAATWkqBljqj1WJjF30xngSizgfUSnptikWdoxlb1u83BzVvKeAT21TPPE21bX2O33mOBJcq+Y4euEQeAIMAACApozWAE4pHaSU7vSU76eUrqaUznf/7p12GgAAAGAZJQQipXRe0j1JBz2Tr+ecL3SfO5J0TdKVU04z9T1eH+OR+xgjfeF0puyCm3K/ckwB8zFlWEury951Y4xGuI1GeQKcc76dcz5cL08p7a997kjSq6eZBgAAAJRMHQN8IOl4vbBr4NZOAwAAAExTZ4E411N2T9LeKaY9IqV0WdLl7td/+cM//MP/t6Ke2F6fk/S/pq4ERsU+bw/7vD3s87Z8I3qGUzeAB5dzviHphiSllN7JOb88cZUwIvZ5e9jn7WGft4d93paU0jvR85w6BKLvqe05LcIbaqcBAAAApqmfAB+qJ5wh53zUvZnunjZAHQEAALBDJn0CvN5g7V5ie+s0005w4zT1xVZin7eHfd4e9nl72OdtCd/faaQcuOe1yNxwTdIbkm7lnG930/YlXdTiafBBzvmNle9VTQMAAAAsozSAAQAAgLmY+iU4AAAAYFRTvwQ3ivVwCUk3cs7Hk1YKp5ZSOpB0vvv1FUnfWu7X0j7neNgNKaXrOecrK7+zz3dUSumiFi8+H0mL0UW7cvb5jun223ktsj3ta8N9yv7eLt39+49zzi+tlVft46r9n3Pe+R8tYo6X/9+XdH3qOvFz6n26J+nyyu8XJd3ZZJ9zPGz/T3eBy2tl7PMd/OnO7asr+47zfId/lvt65fdr7O/d+tHiD5zHruGn2cc1+3/nQyDWh0fOiwwSr05UHcR5WdLrK7/flnSQUtor7XOOh52xr5W83+zznXYtdy8555yPcvfEiH2+sy6s/b4nsb93Sc75ds75cL28dh/X7v+dbwBr8VfG8Xrh+gbDdsmLLtBLK0X7Xfmxyvuc42HLpZQu5pzfXitmn++grpv0eOX/q9jnOyqldKt7mHFe0s2umP29+2r3cdX+b6EB/NiAGeofSQ5bZu0vyNe0SLEnlfc5x8MW6y5ofQPesM93076ke10M8FFK6WrXKJLY57vqkhb7710tUpze7srZ37uvdh9X7f8mXoLDbksp7WlxoVzvOsPuOeh5+ovddU7S+eW5nVK6oUXD6PlJa4UhLcPb9iVdTykpk+cfA2jhCXDfXwHn1PO4HFvrmh4Nhyjtc46HLdU9+bttTGaf76Z7WrzVLelhiNMyzp99vmOW3dldjOgNSS9K+nb3kIP9vftq93HV/m/hCfCheh6P57XhlLGdUkpXJb2ecz5OKe11N0hzn6eUZE0buq4I8Wq3D6VFQ+iyFo1i9vlueuxFGX16w2Of755H/sjt9uUNLfYl+3v3Ve3j2v2/8w3glY0j6eFfmG9NVyNE6eIC386f5vo73/1u7nOOh+21Egso6WEe4Bsrv69OY5/vgG7f3Vv+cbuMAe/O+WP2+c65rU9zuT60bMiwv3db7b27dv83MRTyeoJk4om2X7dP764VH+WcX1yZ3rvPOR62W9cdelmL0Jc3tMj3eMQ+303d/v62Fuf7i5L+IBvJ79nn2697sLF82fWcpNsrDWD29w7owtkO9Ok1/FY2BrfZdB/X7P8mGsAAAADAUgsvwQEAAAAP0QAGAABAU2gAAwAAoCk0gAEAANAUGsAAsMVSSgcppWsppQPn9y5239sbqGoAMFs7nwcYAHbcy1pJDbapnPPbKwnkXd8FgG3HE2AAAAA0hQYwAAAAmkIDGAAAAE2hAQwAM9G9mHYnpXSr+/1qSuluSunyht8/6L5/OaV0vvv3Wkppv5t2ddN5AcAu4yU4AJiJ7sW025K+02VnOMo5v+j4/mH3/ZdyzjckKaV0R9Jf5JzflnSYUrov6cYA1QeArcETYACYkS6bw+uSviPpsGIWP5R0Z+X3e5KOTl8zANgdNIABYH6Oup+LU1cEAHYRDWAAmJEu9OEg53xJ0oWU0vmJqwQAO4cGMADMRPeC2mr4wi1JN1NKVzf8/oGk1yRd6l56u6jFQBlXUkp73fz3Np0fAOyqlHOeug4AgEpdo/Yt70hw3XcvSjrMORMjDKApPAEGAABAU2gAAwAAoCk0gAFgu70j6dvel+W60IlXtEiTBgBNIQYYAAAATeEJMAAAAJpCAxgAAABNoQEMAACAptAABgAAQFNoAAMAAKApNIABAADQlP8f95zeXS2GdW8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 1, sharey=True)\n",
    "axs.imshow(unpad(m.plot()).T, clim=o.plot.clim(m.plot()), extent=model_extent)\n",
    "axs.set_xlabel(m.ax_info[0].l)\n",
    "axs.set_ylabel(m.ax_info[1].l)\n",
    "fig.suptitle(r\"$\\hat{\\mathbf{m}} = \\arg \\min \\Vert \\mathbf{B}\\mathbf{m} - \\mathbf{d} \\Vert_2^2 \\quad$\" + f\"after {epochs} iterations of SGD with lr={opt.param_groups[0]['lr']}\")\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d3ce6d6",
   "metadata": {},
   "source": [
    "## Shot-wise Dask distribution: WIP"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "836f8234",
   "metadata": {},
   "source": [
    "In principle, an `AutogradFunction` can be casted also for Dask-distributed operators.\n",
    "However, distributing a shot to each workers (as done in the [previous tutorial](./2D%20LS-RTM%20with%20devito,%20dask,%20and%20regularizers.ipynb)) will require the objective function to be computed on different workers, too.\n",
    "   \n",
    "Therefore, we need to merge the distribution paradigms of Dask and PyTorch. Any help will be appreciated!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "827a2fec",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}