{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Meep Adjoint Solver - Introduction\n",
    "\n",
    "This tutorial serves as a basic introduction to MEEP's adjoint solver. The adjoint solver provides a simple and flexible interface to calculating \"sensitivities\" or gradients of user specified objective functions with respect to **any number** of design variables -- **all with the cost of just two simulations.**\n",
    "\n",
    "Practical electromagnetic design problems are often constrained by complicated design requirements that aren't easily satisfied with physical intuition. We often formulate our problem as functions of Fourier transformed fields, (e.g. S-parameters, Poynting Fluxes, mode overlap coefficients, etc.). We can \"inverse design\" a structure that maximizes (or minimizes) this cost function with various optimization algorithms. The adjoint method provides an extremely cheap gradient, accelerating many out of the box optimization algorithms.\n",
    "\n",
    "In this tutorial, we'll demonstrate how you can quickly begin leveraging this interface. As a toy example, we'll examine a 2D integrated photonic waveguide bend. We'll code up a cost function that will tell our optimizer to maximize power transport around the bend. We'll discretize the bend into several individual pixels and calculate the gradient of this cost function w.r.t. all these pixels. Finally, we'll compare our gradients with finite difference approximates, which are much more expensive to calculate. Hopefully by the end of the tutorial, you appreciate the simplicity and flexibility this particular package offers. \n",
    "\n",
    "First, we'll import `meep` and `autograd`, a widely used automatic differentiation package. `autograd` wraps around `numpy` and allows us to quickly differentiate functions composed of standard `numpy` routines. This will be especially useful when we want to formulate our objective function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using MPI version 3.1, 1 processes\n"
     ]
    }
   ],
   "source": [
    "import meep as mp\n",
    "import meep.adjoint as mpa\n",
    "import numpy as np\n",
    "from autograd import numpy as npa\n",
    "from matplotlib import pyplot as plt\n",
    "mp.quiet(quietval=True)\n",
    "seed = 240\n",
    "np.random.seed(seed)\n",
    "Si = mp.Medium(index=3.4)\n",
    "SiO2 = mp.Medium(index=1.44)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll begin to specify the domain of our waveguide bend simulation, just as we would with any other simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "resolution = 30\n",
    "\n",
    "Sx = 6\n",
    "Sy = 5\n",
    "cell_size = mp.Vector3(Sx,Sy)\n",
    "\n",
    "pml_layers = [mp.PML(1.0)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll then define our sources. We'll use a narrowband Gaussian pulse, even though our objective function will only operate over a single wavelength (for this example). While we could use the CW solver, it's often faster (and more dependable) to use the timestepping routines and a narrowband pulse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "fcen = 1/1.55\n",
    "width = 0.2\n",
    "fwidth = width * fcen\n",
    "source_center  = [-1,0,0]\n",
    "source_size    = mp.Vector3(0,2,0)\n",
    "kpoint = mp.Vector3(1,0,0)\n",
    "src = mp.GaussianSource(frequency=fcen,fwidth=fwidth)\n",
    "source = [mp.EigenModeSource(src,\n",
    "                    eig_band = 1,\n",
    "                    direction=mp.NO_DIRECTION,\n",
    "                    eig_kpoint=kpoint,\n",
    "                    size = source_size,\n",
    "                    center=source_center)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we'll define our waveguide geometry and \"design region\". The design region takes a 10x10 grid of randomly generated design variables and maps them to a point within the specified volume. Since meep operates under the \"illusion of continuitiy\", it's important we provide an interpolator to perform this mapping.\n",
    "\n",
    "In this case, we'll use a builtin bilinear interpolator to take care of this for us. You can use whatever interpolation function you want, provided it can return either a medium or permittivity (as described in the manual) and you can calculate the gradient of the interpolator with respect to the inputs (often just a matrix multiplication). The built-in interpolator takes care of all of this for us."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "design_region_resolution = 10\n",
    "Nx = design_region_resolution\n",
    "Ny = design_region_resolution\n",
    "\n",
    "design_variables = mp.MaterialGrid(mp.Vector3(Nx,Ny),SiO2,Si,grid_type='U_MEAN')\n",
    "design_region = mpa.DesignRegion(design_variables,volume=mp.Volume(center=mp.Vector3(), size=mp.Vector3(1, 1, 0)))\n",
    "\n",
    "\n",
    "geometry = [\n",
    "    mp.Block(center=mp.Vector3(x=-Sx/4), material=Si, size=mp.Vector3(Sx/2, 0.5, 0)), # horizontal waveguide\n",
    "    mp.Block(center=mp.Vector3(y=Sy/4), material=Si, size=mp.Vector3(0.5, Sy/2, 0)),  # vertical waveguide\n",
    "    mp.Block(center=design_region.center, size=design_region.size, material=design_variables), # design region\n",
    "    mp.Block(center=design_region.center, size=design_region.size, material=design_variables,\n",
    "             e1=mp.Vector3(x=-1).rotate(mp.Vector3(z=1), np.pi/2), e2=mp.Vector3(y=1).rotate(mp.Vector3(z=1), np.pi/2))\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now formulate or simulation object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sim = mp.Simulation(cell_size=cell_size,\n",
    "                    boundary_layers=pml_layers,\n",
    "                    geometry=geometry,\n",
    "                    sources=source,\n",
    "                    eps_averaging=False,\n",
    "                    resolution=resolution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll start defining our objective function. Objective functions must be composed of \"field functions\" that transform the recorded fields. Right now, only the Eigenmode Decomposition monitors are readily accessible from the adjoint API. That being said, it is easy to extend the API to other field functions, like Poynting fluxes.\n",
    "\n",
    "In our case, we just want to maximize the power in the fundamental mode at the top of the waveguide bend. We'll define a new object that specifies these details. We'll also create a list of our objective \"quantities\" (just one in our case)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "TE0 = mpa.EigenmodeCoefficient(sim,mp.Volume(center=mp.Vector3(0,1,0),size=mp.Vector3(x=2)),mode=1)\n",
    "ob_list = [TE0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our objective function will take the magnitude squared of our predefined objective quantity. We'll define the function as a normal python function. We'll use dummy parameters that map sequentially to the list of objective quantities we defined above. We'll also use autograd's version of numpy so that the adjoint solver can easily differentiate our objective function with respect to each of these quantities. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def J(alpha):\n",
    "    return npa.abs(alpha) ** 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now define an `OptimizationProblem` using our simulation object, objective function, and objective quantities (or arguments). We'll also tell the solver to examine the `Ez` component of the Fourier transformed fields. The solver will stop the simulation after these components have stopped changing by a certain relative threshold (default is 1e-6)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "opt = mpa.OptimizationProblem(\n",
    "    simulation=sim,\n",
    "    objective_functions=J,\n",
    "    objective_arguments=ob_list,\n",
    "    design_regions=[design_region],\n",
    "    fcen=fcen,\n",
    "    df = 0,\n",
    "    nf = 1,\n",
    "    decay_fields=[mp.Ez]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now initialize our object with a set of random design parameters. We'll use `numpy`'s `random` library to generate a uniform random variable between 1 and 12, to correspond to the refractive index of air and silicon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "x0 = np.random.rand(Nx*Ny)\n",
    "opt.update_design([x0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we can visualize our final simulation domain with any extra monitors as defined by our objective function. This `plot2D` function is just like `Simulation`'s `plot2D`, only it takes an additional first argument. We'll set it to `True` to tell the solver to initialize the solver and clear any stored fields."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAToAAAEGCAYAAAD1+lmKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO2de5RcxZ3fP9U93dM9LzECCfQAjwwSwg/APLxr72Zj1sHBvMRDYlkjnrvGnF1vvHFsJyecBJG1s9jOMZtkWWdJwIZdbK8FOMZYlkFGHOI9OEGsjUEWYPQwkkZCQpoZzaPfXfmjuy49o3l0T9/urqr+fc6po+mZe39dunV/31u36lf1U1prBEEQfCbS6goIgiA0GhE6QRC8R4ROEATvEaETBMF7ROgEQfCejlZXoBZUt9IsaHUtBEGwhgO8rbVeNOdxWmtnCkvQYXLv8/dqtUHpe5+/N1S7WmvNBvTW3Vv1SV85SW/dvTV0+1prr+wPDg7qK6+8UgMa0Eqp4OdaSiQSCX7+nRt+Ry+8Z6EX16dZ9tkQno810r+MfWCbrkY7qjnIlhKm0DW6EdiAd07QSPtTha5SsOYldAPo+J1x/diLjzWl/r7YD0vomiFyaoMSoavmIjX0SVPu0TUCH50sVKEbQPP5Uo/uwIEDTam/L/bDELpmiZz06Kq8SI0geNKE+ApQia9OFpbQqRVK8/lSj27NmjWhC52v199Q733bTJHTWovQVXORwmbSk6YBQtdqJ2ik/VCErtyTY6D0OWyh8/n6G+q5b5stclqL0FV1kcLkuCdNyEJngxM00n7dQlcWObXinUmMMIWu1denWfbne9+2QuS0FqGr6iKFxbRPmhCFzhYnaKT9uoSuoidXeV5YQmfD9WmW/fnct60SOa1F6Kq6SGEw45MmJKGzyQkaaX/eQjfldTVsobPl+jTLfq33bStFTmsRuqouUr3M+qQJQehsc4JG2p+X0E0RubCFzqbr0yz7tdy3rRY5rUXoWt4I9QqdjU7QSPs1C900Ihem0Nl2fZplv9r7ttX+ZWhrobOhEeoROludoJH2axK6GUQuLKGz8fo0y341960N/mVoW6GzpRHmK3Q2O0Ej7Ve9BGwWkQtD6Gy9Ps2yP9d9a4t/GdpS6GxqhPkIne1O0Ej7VfXo5hC5eoXO5uvTLPuz3bc2+Zeh7YTOtkaoVehccIJG2p9T6KoQuXqEzvbr0yz7M923tvmXoa2EzsZGqEXoXHGCRtqfVeiqFLn5Cp0L16dZ9qe7b230L0PbCJ2tjVCt0LnkBI20P+MYXQ0iNx+hc+X6NMv+1PvWVv8ytIXQ2dwI1Qida07QSPvT9uhqFLlahc6l69Ms+5X3rc3+ZfBe6GxvhLmEzkUnaKT943p0FbuQVCtytQida9enWfbNfWu7fxmsFzrgVGAr8CtgO/CZOc9Z4k4jzCZ0rjpBI+1PErp59ORqEToXr0+z7LMBJ/zL4ILQLQHOK//cC7wOvGfWc5a40wgzCZ3LTtBI+4HQlUWucheSMIXO1evTLPtswAn/MlgvdMdVBL4PXDzrMUvcaYTphM51J2ik/cHBQf3hT3x42l1IwhI6l69Ps+ybHl0jaEQnxSmhAwaAN4G+af52O7AN2MYCdxphqtD54ASNtP/otkd1/M548LoaRnKcSqFz/fo0y36jdsZu1JuYM0IH9AAvAtfMeWzIWcAMDXnSVNwwvjhBI+0vvGdhqUc3jWDVK3Q+XJ9m2W+E0DVyuMkJoQNiwI+Bz1Z1fAOErmFPmvIN45MTNNL+o9seDT3d4Zo1a/RjLz7mxfVplv2wha7RY+rWCx2ggIeBv6r6HMnrOgmf7Iee7hDJ6zof+5LXNXyh+93yDflL4Bflcums50he1wDfnCz0BNYDktd1PvYlr6sFRfK6lvDRySSvqx32Ja+rBUXyuvrrZJLX1Q77ktfVgiJ5Xf11Msnraod9yetqQZG8rv46meR1tcO+5HW1oEhe162h2bTNvuR1tcO+5HW1oEhe1/Cxxb7kdbXDvuR1taBIXtdwscm+5HW1w77kdbWgSF7X8LDNfs1xdNOIXJhCZ9v1aZZ9yetqQZG8ruFgo33J62qHfcnrakGRvK71Y6v9qoVuFpELQ+hsvT7Nsi95XS0okte1Pmy2L3ld7bAveV0tKJLXdf7Ybl/yutphX/K6WlAkr+v8cMG+5HW1w77kdbWgSF7X2nHF/oxCV4PIzUfoXLk+zbIveV0tKJLXtTZcsj9teEmNIler0Ll0fZplX/K6WlAkr2v1uGZ/qtDNN+Wh5HWtD8nrakGRvK7V4aJ9yetqh31f87qq0rFuoJYqzadaXQtBEKxhAy9qrS+Y87hq1NCWwhJ3tidv1MabvtKInBGN2HjTd1zbGZsqe3SR0BW2wWxct5GPDHwkdLvP7nmWdRvXNcy+UBvaoTcN3/DRv5wTOh8bQTgepVSrqyCEhA3+5ZzQhY0NjSAcj/To/MAW/2probOlEYTjkR6d+zTav57d82zVx7at0InICULjaIbIrdu4rurj21LobHrSCNMjr67u0iyR27huY9XntJ3Q2fakEaZHXl3dpJkiV4v9thI6G580guALtooctJHQ2dwIguA6tvtXWwid7Y0gCC7jgn95L3QuNIIguIor/uW10LnSCILgIi75l7dC51IjCMcj4SV245p/eSl0rjWCcDwSXmIvLvqXd0LnYiMIgiu46l9eCZ2rjSAILuCyf7VU6JRSDyqlDimlXqnXlsuNIByPjNHZhev+1eoe3TeBS+o14nojCMcjY3T24IN/tVTotNbPAUfrseFDIwiCrXjjX9Xst97IAgwAr1R17JS8rjZnW3IpZwRIaXRxhcr71mb/MuBKusO5hA64HdgGbGOBO40gQifFZaGz3b8M3gjdpGOXuNMILgmdDd4oWcDswGQBs92/DNUKXasnI2rGmzEDYVZK97DQCnz0r1aHl3wbeB44Uym1Tyn1R3Od42MjCMcjs66tw0f/6mjqt01Ba/2HtZ7jYyMIxyM9utbho3859+rqYyMIxyM9On+wwb+cE7qwsaERBMFXbPGvthY6WxpBEHzEpmx7bSt0InKC0Dhsy7bXlkJn05NGEHzDxmx7bSd0tj1pBMEnbI1zbSuhs/FJIwi+YKvIQRsJnc2NIAiuY7t/tYXQ2d4IguAyLviX90LnQiMIxyMrI9zAFf/yWuhcaQTheGRlhP245F/eCp1LjSCUenBhi5tSSnqGDcI1//JS6FxrBKFxoiQ9w/Bx0b9auntJI3CxEdqVYFPESIRCoTBJ6IrF4rxsVp5XKBSCz5W/n0n8RBTnxlX/8kroXG2EdkUpFYhLd3c3sVgMgHg8TjweR2s9p+BFIpFJJRqNkslkGB0dpauri3g8Hhwn1IfT/lXNNsS2lKnJcSqxbftn2Up95vwJNpd2Z6b71jb/MuDrVurT4fSTRhAsxwf/cv7V1YdG8J2XX34FgGw2G5ShoSFeffVVXnrpJX7+85+zZ88eTjzxRFasWMHixYtZvHgxXV1dvPnmu/nhD2/m4x//BsuW/Rp4J8YuGo2SSCTo7u6mp6eHrq4uCoUCx44dY3h4mMHBQQ4dOgRAIpEgkUjQ2dkZlHg8zqFD7+Gxx/6AVKqnNRfHcnzxL6eFzpdG8J1oNBqEj5gSiURIJBL09fWxePFiUqkUkUiEXC5HLpcjk8kwOLiKp566iY997H+yaNFrZLOTw1Dy+TzApM/FYpGxsTHGxsYA6Ovro6Oj4ziBi8ViDA6u4nvfu47rrvsuDz10S0uujc345F/OCp1PjeA7uVyOYrFIOp0mlUoxMTHByMgI4+PjZLNZAGKxWDBhkM/n2bv3dH760z/mox/9W04++XVmmpPIZrNorclkMsRiMfL5PKlUilQqRS6XmzzGC8EEx549A3z/++u4/vrHWLp0Z1Oug0v45l9OCp1vjeA72WyWQqHA+Pg4o6OjjI2NcfToUUZGRkin02it6ejoCHp7hw+/l1/+8s/4vd/7b5x44uvkcpFJPUGY3IszNgqFQtCrKxQKFAoFoDTjagSuWCzym9+sYNOmT3DVVY+wbNk+Uql0y66NjfjoX84JnY+N4DvZbJZ8Ps/4+DjHjh1jaGiIoaEhRkdHyWQyKKWIx+MopRgZ+QBvvHEn55//l/T1vU4mE6Ojo4NoNEo0GqVYLBKNRoFS7yyXyzExMRGUYrFIJBIJzjHhJaZXt3//SjZvvoXLL3+IJUt2kcvFg16lUMJH/3JO6HxsBN+ZmJggl8sxOjrKyMgIQ0NDHD16lLGxMfL5fCB0o6Pns3Pnl3jf++7mhBNeJZcrxdVV9uQikUjw2fTiUqkUIyMjHDt2jHw+T2dnJ93d3XR1dQXnK6U4eHA1W7bcxmWXfZPly3eitaJQKJDL5Vp2bWzER/9yTuh8bASfqJwsMOzbt49iscjQ0BDDw8NBby6VSlEoFIhEIqTTv83OnX/B+953NyefvINYLE5HR0cQPGx+Nj07I3yFQiHovRlM4HAsFgsmII4ePZstWz7FFVc8zIoVe1GqMzjW9PqEEj76l3NC52Mj+M727duJxWKMjo4yPj5OOp0mnU4H42mjo+fz2mt38oEP/GdOOWUnHR1dweyomSE1/0aj0UnCppQin8+Tz+fJ5XLk83mSySQ9PT1Br+7Ikffzk598kmuu+TYrVx4iGn0nlESWfTUeG/zLOaELGxsawXdefvllurq6yOVywfpT0/M7evRsXnrpTs4//x6WLHmDeLxrUgjI1B5dR0cHsVhp3A5KoSvm9dMInomtSyaTDA2dw9NP/xFr136X1asP09nZG7z2wuTXYiF8bPGvthY6WxrBd7Zv387ChQsBglfKWCzG+PiFbN9+Jxde+BWWLdtJPF5am2oCe43YmZ5c5efKWVoThJzP58lms4HQHT16Nj/+8a1cd92jrF59iESiJKJThU56dY3Bpmx7zgndww8/HIqdHakd3HfoPv508Z/y5nNv8vBz4ditJKy6Npqbyv/WW18jHgDJZBK4DoDx8XE6OjqC8TDzuvrGG3eVX1ffoLOzi2QyOWkFQ6WYQSmUxIy/QUmkTE8vkUgEi/g7Ozs5cuT9bN5cErlVqwbp6IgHW0FV1lN6dI3Btmx7qrLRbUctVbrjcP3aXHxXkeI1RSKPR4j8JvybvPiuIsUbi3R80Y3nSK68wiDWEcK1LYd3JBIJxsZGAXj/+88GSr05pRQTEx9k586/5Oyzv8ipp+6kq6srGE/r6ioJXkkoIZPJkE6nyWQyFAoFOjs7SSQSJJNJYrEYuVwuiM9LpVIUi0UOH34vmzffxrp1G1m5cn8wU1sZi2dKR0cH6XSaL3zh80BpaX87o+5W6LvquwjNDAG7aMVFL2qtL5jrHDc8sQKz7GfeDADXABuhuKdIkfntezanfUKoa5MJq75mGZZhyZIlwSqFoaFz2LVrA+997wZOPvl14vGeYGbU9Mp6enpIJpPk83kymQwTExMMDw+TTqeD11IzBmdi6UyIiomTW7v2u6xcOUgsFqNQKJDNZoNxPCNwZgxQ4ujCw9Y4V+eErqOOXseknty+SOj/+0r7xRuLddW1qZQFrt76mhlQpRSLFy/m4MHS7xcvXkw2m2VwcBU7dnyOc875IgsX7iAeTwaTDfF4nGQyGSzQTyaTZLNZRkdHSafTDA0NMTIyEqyPNeNylZMS+/evZNOm9eWe3GAwHmfWzk5MTJBOp1FKBeJqvkeoH1tFDhwUunn3OgZoTk+ubB/at0c3NZauv7+fnTtP5YUXPsWHPvQ1+vpeQ+t3BM705sysamX4SC6XI51OB4HGiUQCeGfDTSNmg4OreOKJG1m79h84/fR9RKMxlFLBcjAjdBMTE8HvlVLEYrF572YsvIPNIgcOCt28GADWARuBPQ7ad5zBwVU888zNXHTR39Df/zrpdOm2qwwV6ejooFgskslkKBaLjI+PMzExwdGjR4PxNxN3B0wad9u37wyefHI9V131CMuX7yGbjQSvvZWbCZg1t2YWNxaLkUgkgjWxwvywXeSgHYRuABG5FlA5ybVp0y1cfPH9nHTSa+RyJYExr4+mR2fWsY6OjpLP54NNAI4cOcLIyAjZbBalVBAwbFY9HDhwJk8+uZ4rr/x7li3bSSZTnLSAv1gsBls/FQqFYKa2s7MzGA8U5o8LIge+C90AInIWcMklD7Jo0Wvk8yroSZnJg8o4uVwux9jYGCMjI4yNjQW9uomJiWDsz/QA4/E4Bw+u5okn1pd3IdkZ7F5SucFnZdIdI5QmHKW7u5ve3l4JL5knrogctFjolFKXAP8ViAL/S2t9T2jGBxCRs4SlS18nm9VBoK6Jqatc2hWLxYJxtKGhIQ4fPszY2FiwisKcZ4Ty4MHVPP30DVx77XdYvnw32ew7vbd0Os3ExESwJ13ld5l1spXrYE1SHqF6XBI5mEXolFKbgD/RWu+p+1umtx8F7gMuBvYBLyilntBa/6pu4wOIyFmEGVerjGMz43JTx+lMUK/piVXGu5nX3EOH3sPPfvbHXHrpN1i0aDeZTOl7zHFKqSCkxGwDZdbCFgoF4vE4mUyGTCbj3ISRDbgmcjB7j+4bwFNKqYeAr2itw97L5oPAG1rrXQBKqe8Aa4D6hG4AETlLmbrdUqXgGREzC/JTqRTxeDzofZmxuaNHz+aFF/6cD3/4a3R2vsbwcOK41RS5XC4414zTmVffyo0+Te+wcucTYXZcFDmYRei01huVUj8C/gOwTSn1d/BOTIbW+mt1fvcyYG/F533Ab009SCl1O3A7AAvmsDiAiJylmHAT8/paKXLm1dXsI5dOp4PYt8q/v/XWWWzb9idceOFXSCRe5O23CySTSfr7+4MQFSDYjNP0Dk2gsOlZmq2hjNDJGF11uCpyMPcYXRYYBzqBXgg7+GxutNb3A/dDaQkYIzMcOICInKVEIpFgxhOYtAzL7C9nBM306MyOJGZs7eDB1Tz33M1cdNF9dHa+wpEjpQQ4vb29JBIJ+vv7g96Z6dmZYGIzbmfCTfL5fCByZmdjYXZcFjmYfYzuEuBrwBPAeVrriZC/ez9wasXn5eXf1c4AInIWY2LkzNibGX8zM6O5XC4QJbPwP5lMBhMVBw6cyaZNpdnVBQv2cOhQdFJvrbKYV1UTRmLG/8yx2Wy2vNFnOvjuAwfObPEVshvXRQ5m79HdCazTWm9vyDfDC8BKpdQKSgJ3PfCJmq0MICJnIZW79nZ1daGUIpPJkMvlSKVSAMH612KxOCmbl9kYIBKJlIOBb+Dqq7/FqafuYny8g+7ubgqFAolEItgowKytjUajwQYA8Xic3t7eQOTM95teX1dXF8PD5/LUUze05Bq5gA8iB7OP0f2zhn1ryX5eKfVp4MeUwkserFlUBxCRs5TKAf6enp5g1UMmkwlyO5hdhs2rpQka1loTjUZ588138+STN3LllX/P8uW7KRRK64d7enpIJBJBL64y8Q4w6RXWCK75DtN7TCaTHDt2Hs88s55bb93Effetbcl1shlfRA5aHEentd4EbJrXyQOIyFlMZ2dn8HN3dzfZbJaxsbEgKHh8fDxYlmV6W5XBwwcOnMkPfnBDecXDr8nlSkG/ZlWDedUdGxvjyJEjQXYxrTV9fX309/fT29tLLBZDa006nWZsbIx0Ok1nZydjYxewbdvt3HHHFs47L81997XkMlmLTyIHrq6MGEBEznIqe3SVqQorA3rNWJqZ9TR70A0OruKHP1zP1Vd/i+XLd5PL6SDezfTUkslkcN7IyAiZTIaRkRGKxWLQizOvp2btrNmpJJP5EDt2fI7rr3+Uc85RJJOyDKwS30QOHBS6yOmRpm21VI/9PPm23KbJrESo3BHELKrP5/PHhZaYeDazPOvgwdVs3nwza9f+AytW7COdJth9xMyYGkzMnNaaeDxOd3c3xWKRrq6uIK7OiOKCBQvI5/OkUr/Fjh13c/HF9zMwMEY2u4zx8fG6/t++4ZvIgYNCV7ym2LStluq171rUfZgbbx44cCD4bNauZspLGEyPy6w3NTsMl7Y/v421a7/LqlWDFIuTt2oy+89lMhnGx8eJxWJks1mKxSLJZJJFixYRjUZZsGABvb29wSuu6dUdPvxeXn31dn7/97/O6acPks/3MjIyU7xS++KbyIGDQhd9PEpkfwRCXp5YPK1I4ZpCaPZz5NxZQ1lO4Fxvfc1Sq0gkQm9vL8PDpd8fOXIkGI8DJgmQ2Sl4qsjF4/FAeCv3kzM9PzOWZ7Zr6unpoa+vLwhNMfkjzHH79p3B1q1rufbab3HaacNASfyOHTsmPbop+CZy4KDQfeOub4Ruc0dqB3/91l/z6ZM/zVl3nRWKzZt23cQDDzwQiq2Gc1MpPU4Y9TWbbiaTSdaVc5e89dZbwfpSsz1TJBIJBGl4+NyKHA+DwQac5vXX7C1nEl7DO4lxTK/Q5HE1G3hWrp/dtes0Nm68nFtv3cQZZ4yTzfaTSqWChf9DQ0Mz/j+E+mm1yIGDQnfjjTeGau/ZPc/y2Y2f5Ykbnwi1EW66+6bQ69owykLXqPpms9lgBULl7iGmJ/f007dx3XWPcsYZ+4lGo2itg+2WTBCwCTY2uVuNCJoxQbP1UuU+eNFolN2738XDD1/GzTf/kHe/ey+Fgg5smH3vZOPNxmGDyIGDQhcmtjSCT0zXE5puqZVZoL9lyydZt24jZ555EKU6grg604tLp9PB2lSTKtEEFptlYmY2N5PJ0NPTQ09PD729vezdezqPPHIZ69c/wSmn/Jq3306RTqcnJdDu6CgFIE9FenP1I3ldLUBErnmYGdDKkJO33jqLZ565nWuu+TarV78dbGluJh5SqVQgdMViMdgw06ygMAHHZquldDrN+Pg4fX19FAqF8s7Da7j11h+xdOlujh49xpEjR4IVFGZ34d7eXtm9pAHYlte1LYXOpieNj0zt1U0Vuv37V/LMM5/kqqseYfXqQyQSpZSHlTkeTHCvea00Y3ZmwsQsJas8Ph6PlzOGncMLL6znxhv/N6tWHWF4OM/o6CiHDx/m7bffJhaLcdJJJ5FMJoPwEyE8mprX9QsXVXVO2wmdbU8a35julc9MDkQiEfbuPZ2nn76Nyy77JqeffpBYrCeYODCzrOZVtHL7dNOjM8Jm7JnXUPPv8PC57Nr1b7j00gcYGBinUOgKxvvMtuxmOyczTli5LleoD1uDjdtK6Jr6pHmouidNO2DWnO7adRqbN9/EZZd9k9NO24VS3ZO2bjKTFZV5Jcw+dWYp19QVF1rrIDB4fPxCtm+/i4su+joDA2+Ry50QhLWYEBQTc9ff3x/E24nQhYOtIgdtJHQ2N4LvxONx9u9fyZNP/gGXX/4wy5fvCv5m8kGY193KJV7FYjFYxmVCUhKJRNCbSyQSQZhKae3qn3HFFQ9x2mmHiMXiwW4l2WyWzs5OFi5cSG9vL93d3Zx44omceOKJ9PX1uRPvaDG2+1dbCJ3tjeA7Bw+u5vHHr+faa7/DsmVvksupSXlZDWZpWGdnZxA+YiYNurq6AnEzy8a6u7vp6Ojg2LHz+Md/vJVbbnmSM8/Mk04vJpUqzbCardPNduuxWIzu7m76+/s54YQTWLBggTtL9SzFBf/yvoVdaATfefzx67nhhu+zdOlestnIpNdPE/dm/jWvriZjlwkGNrOjZhNNs761tDb2Ju644ydceKEGTgkW95vE12YVhinGnundyazr/HHFv7wWOlcawXeuuOLvOPnk/YyOTgS7CgPBOlXT4yoUCqRSqaA3VrmvnJmoMBMLqVSKnTtPZcuW9XziE4+xatUEkUh/sMnm2NgYw8PDjI6O0tnZOSkWzwQJV+5GLNSOU/5VOUZie2EJulq27t6qT/rKSXrr7q1Vn1MLc9lnQ/V1bTlQKg0y61Jpd6q9b1vtXwZgm65CO7xMf+TUk0YQHMNF//Lu1dXFRvCVimWngie46l9eCZ2rjeAzOgS1mzqGNp3NmY7R06y9lTG5+eGyf3kjdC43gs80QlSqsVkZiCzUj+v+5cUYneuNIAg244N/OS90PjSCINiKL/7ltND50giCYCM++ZezQudTIwiCbfjmX04KnW+NIAg24aN/OSd0PjaCINiEj/7lnND52AiCYBM++pdzQudjIwiCTfjoX84JnY+NIAg+Y4N/OSd0YWNDIwiCr9jiX20tdLY0giD4iE3Z9tpW6ETkBKFx2JZtry2FzqYnjSD4RjNDwKql7YTOtieNIPiErXGuLRE6pdQ6pdR2pVRRKXVBs77XxieNIPiCrSIHrevRvQJcAzzXrC+0uREEwXVs96+WbLyptd4BzdsU0fZGEASXccG/vB+jc6ERBMFVXPGvhvXolFJbgFOm+dOdWuvv12DnduB2ABbUVgdXGkEQXMQp/6omJ2KjCvAscEHVx0teV0FoKJLX1SGcetIIgmM46V/VqGHYBbga2AdkgLeAH1d1XhU9OmueNNKjExxkrvvWFv8yUGWPrqWvrrWWuYTOpkYQoRNcZLb71ib/MlQrdN68ujrZnRYER3Ddv7wQOtcbQRBsxgf/cl7ofGgEQbAVX/zLaaHzpREEwUZ88i9nhc6nRhAE2/DNv5wUOt8aQRBswkf/ck7ofGwEQbAJH/3LOaHzsREEwSZ89C/nhM7HRhAEm/DRv5wTOh8bQRB8xgb/ck7owsaGRhAEX7HFv9pa6GxpBEHwEZuy7bWt0InICULjsC3bXlsKnU1PGkHwDRuz7bWd0Nn2pBEEn7A1zrWthM7GJ40g+IKtIgegSnvXuYFaqjSfanUtBEGwhg28qLW+YM7jqtmd05ZSS3KcSu59/l6tNih97/P3zuv8+dgPc4dhG3d2Fft+2p/PfdsK/zLQjlup13qRwmAm+2EJnU1OIPb9t1/rfdtKkdNahK6qi1Qvsz5pQhA625xA7Ptvv5b7ttUip7UIXcsboV6hs9EJxL7/9qu9b1vtX4a2FjobGqEeobPVCcS+//aruW9t8C9D2wqdLY0wX6Gz2QnEvv/257pvbfEvQ1sKnU2NMB+hs90JxL7/9me7b23yL0PbCZ1tjVCr0LngBGLff/sz3be2+ZehrYTOxkaoRehccQKx77/96e5bG/3L0DZCZ2sjVCt0Lgj9xAwAAAclSURBVDmB2Pff/tT71lb/MrSF0NncCNUInWtOIPb9t19539rsXwbvhc72RphL6Fx0ArHvv31z39ruXwavhc6FRphN6Fx1ArHvv3024IR/GbwVOlcaYSahc9kJxL7/9tmAE/5l8FboXGmE6YTOdScQ+/7bNz26RtCIToq3QudKI0wVOh+cQOz7bz/M7cUqadSbmLdC1wga8qSpuGF8cQKx77/9RghdI4ebROiqpGFPmvIN45MTiH3/7YctdI0eU7da6ICvAq8CvwS+B5xQ1XkhC11DnzQb8M4JxL7/9sMUumZMHNoudB8DOso/fxn4clXnhSh0DX/SbMA7JxD7/tsPS+iaFR1htdBNqgBcDTxS1bEhCV1TnjTlHl0j8NXJxH7r7YchdM0MAXNJ6H4ArJ/l77cD24BtLHCnERo1e+Wzk4n91tuv975tdpxry4UO2AK8Mk1ZU3HMneUxOlWVzTp7dE190jRA6FrtBGLff/v13LetCOZvudDN+cVwC/A80FX1OXUIXdOfNCELnQ1OIPb9tz/f+7ZVK5asFjrgEuBXwKKazpO8rqHZFPtifzokr2u4QvcGsBf4Rbn8j6rOk7yuoSP2xX4lktfVgiJ5XcNF7Iv9qUheVwuK5HUND7Ev9qdD8rpaUCSvaziIfbE/E5LX1YIieV3rR+yL/dmQvK4WFMnrWh9iX+zPheR1taBIXtf5I/bFfjVIXlcLiuR1nR9iX+xXi+R1taBIXtfaEftivxYkr6sFRfK61obYF/u1InldLSiS17V6xL7Ynw+S19WCInldq0Psi/35InldLSiS13VuxL7YrwfJ62pBkbyusyP2xX699iWvqwVF8rrOjNgX+2HYl7yuFhTJ6zo9Yl/sh2Vf8rpaUCSv6/GIfbEfpn3J62pBkbyukxH7Yj9s+77mdVVaa1xBKTUKvNbqelTJScDbra5EDbhUX5fqCm7V16W6Apypte6d66COZtQkRF7TWl/Q6kpUg1Jqmyt1Bbfq61Jdwa36ulRXKNW3muMija6IIAhCqxGhEwTBe1wTuvtbXYEacKmu4FZ9XaoruFVfl+oKVdbXqckIQRCE+eBaj04QBKFmROgEQfAe54ROKfUXSqlfKqV+oZR6Sim1tNV1mgml1FeVUq+W6/s9pdQJra7TbCil1imltiulikopK0MMlFKXKKVeU0q9oZT6d62uz2wopR5USh1SSr3S6rrMhVLqVKXUVqXUr8r3wGdaXaeZUEollFL/Tyn1Urmud895jmtjdEqpPq31sfLP/wp4j9b6jhZXa1qUUh8DntFa55VSXwbQWv/bFldrRpRSZwFF4G+Bz2mtq4pRahZKqSjwOnAxsA94AfhDrfWvWlqxGVBK/R4wBjystX5fq+szG0qpJcASrfU/KaV6gReBq2y8tkopBXRrrceUUjHgp8BntNY/m+kc53p0RuTKdAPWKrXW+imtdb788WfA8lbWZy601ju01javPPkg8IbWepfWOgt8B1jT4jrNiNb6OeBoq+tRDVrrA1rrfyr/PArsAJa1tlbTU14BNlb+GCuXWXXAOaEDUEp9SSm1F7gB+I+trk+V3Ab8qNWVcJxlwN6Kz/uw1BldRik1AHwA+L+trcnMKKWiSqlfAIeAp7XWs9bVSqFTSm1RSr0yTVkDoLW+U2t9KvAI8Gmb61o+5k4gT6m+LaWa+grti1KqB3gM+PMpb09WobUuaK3PpfSW9EGl1KxDA1auddVa/4sqD30E2ATc1cDqzMpcdVVK3QJcDnxUWzAgWsO1tZH9wKkVn5eXfyeEQHm86zHgEa31462uTzVorYeVUluBS4AZJ32s7NHNhlJqZcXHNcCrrarLXCilLgG+AFyptZ5odX084AVgpVJqhVIqDlwPPNHiOnlBeYD/AWCH1vprra7PbCilFpkIBqVUktLk1Kw64OKs62PAmZRmB38D3KG1tvKprpR6A+gEjpR/9TNbZ4gBlFJXA/8dWAQMA7/QWv/L1tZqMkqpS4G/AqLAg1rrL7W4SjOilPo28BFKWx+9BdyltX6gpZWaAaXU7wL/B3iZkm8B/Hut9abW1Wp6lFJnAw9RugciwHe11v9p1nNcEzpBEIRace7VVRAEoVZE6ARB8B4ROkEQvEeEThAE7xGhEwTBe0ToBCco766xWym1sPy5v/x5oLU1E1xAhE5wAq31XuDrwD3lX90D3K+13tOySgnOIHF0gjOUlyi9CDwIfBI4V2uda22tBBewcq2rIEyH1jqnlPo8sBn4mIicUC3y6iq4xseBA4DVG1kKdiFCJziDUupcSgu4fxv41+VdcQVhTkToBCco767xdUr7pL0JfBX4L62tleAKInSCK3wSeFNr/XT5898AZyml/nkL6yQ4gsy6CoLgPdKjEwTBe0ToBEHwHhE6QRC8R4ROEATvEaETBMF7ROgEQfAeETpBELzn/wNAVD00iU3IQgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "opt.plot2D(True,frequency=1/1.55)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see an expected monitor on the top of the waveguide, but we also see an additional monitor spanning our design region. This monitor records the Fourier transformed fields needed to calulcate the gradient.\n",
    "\n",
    "Since everything looks good, we can go now calculate the gradient and the cost function evaluation. We do so by calling our solver object directly. The object returns the objective function evaluation, `f0`, and the gradient, `dJ_du`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting forward run...\n",
      "Starting adjoint run...\n",
      "Calculating gradient...\n"
     ]
    }
   ],
   "source": [
    "f0, dJ_du = opt()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visualize these gradients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f18a6d40be0>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPUAAAD4CAYAAAA0L6C7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAMFklEQVR4nO3dXWjddx3H8c+nSZo2ydas6MXazqUMV6mTMQk6NxTcvPAJdzOx4iYKUhAf5gPI9MZ7EdELUcrUG4e76AbKGD7gAyJIMe2KtamT0dU0sdZ27rFbk5zk60Ui1HbJ+ef09/OffHm/YNDkpN99Sc+7/5OTk18dEQKQx6a2FwBQFlEDyRA1kAxRA8kQNZBMf42hfUPDMTC6vcboKjZ1Ksycr/NdhU2dOnM7W+r8/d63fb7K3MEKf2gvzw0WnylJI5tni898+czLuvj8Rb/WbVWiHhjdrrFPfan43Kj0uGLobPlQhs8uFJ8pSYPP1onkuT1bqszdtm+mytw3jDxXfOYfp8aKz5Skd954svjMn3788RVv4+E3kAxRA8kQNZAMUQPJEDWQDFEDyTSK2vZ7bT9l+2nbD9ZeCkDvukZtu0/SdyW9T9JeSR+1vbf2YgB60+RK/TZJT0fEyYiYk/SIpHvqrgWgV02i3inp9CVvTy+/73/Y3m97wvbEwisXSu0HYI2KPVEWEQciYjwixvuGhkuNBbBGTaKekXTDJW/vWn4fgHWoSdR/kvRG27ttb5a0T9LP6q4FoFddf0orIjq2PyvpF5L6JP0wIo5X3wxATxr96GVEPCHpicq7ACiAV5QByRA1kAxRA8kQNZAMUQPJVDl4UKpzSODcdYvlh0ra/FL5Za+ZqPP6nM50nbkLt95RZe43bjpYZe7hi2PFZ/5+5s3FZ0rSsZHri898dWFgxdu4UgPJEDWQDFEDyRA1kAxRA8kQNZAMUQPJEDWQDFEDyRA1kAxRA8kQNZAMUQPJEDWQDFEDyRA1kAxRA8kQNZAMUQPJEDWQDFEDyVQ5TXRTRxo6G8Xn1jj1U5JU4ZDSzo7t5YdK6o/yn1dJmru2ylgNuM4JsJOv7Cg+c2Sqzv3r7GD5+0Jntm/F27hSA8kQNZAMUQPJEDWQDFEDyRA1kAxRA8l0jdr2DbZ/a3vS9nHbD/w/FgPQmyYvPulI+nJEHLF9jaTDtn8VEZOVdwPQg65X6og4ExFHln/9kqQTknbWXgxAb9b0NbXtMUm3STr0Grfttz1he6Lz6oUy2wFYs8ZR2x6R9KikL0TEi5ffHhEHImI8Isb7tw6X3BHAGjSK2vaAloJ+OCIeq7sSgKvR5NlvS/qBpBMR8a36KwG4Gk2u1HdKul/SXbaPLv/3/sp7AehR129pRcQfJPn/sAuAAnhFGZAMUQPJEDWQDFEDydQ5eHA+NHx2ofjcayZmis+U6hwS+OxbRorPlKR/f6zOC3tuuuV0lbmTs9dXmfvrqZuLzxw6X+dQx9nR8pm5s/Jz11ypgWSIGkiGqIFkiBpIhqiBZIgaSIaogWSIGkiGqIFkiBpIhqiBZIgaSIaogWSIGkiGqIFkiBpIhqiBZIgaSIaogWSIGkiGqIFk6pwm2gkNPjtffG5nus5pov1R/hTJWqd+nvzw96vMfWHx1SpzP/K3e6vM7RwdLT5zsb/OaaKLAxXmrvIPYXGlBpIhaiAZogaSIWogGaIGkiFqIBmiBpJpHLXtPttP2n685kIArs5artQPSDpRaxEAZTSK2vYuSR+Q9FDddQBcraZX6m9L+oqkxZU+wPZ+2xO2J+bmLxRZDsDadY3a9gcl/SsiDq/2cRFxICLGI2J880Cd1z0D6K7JlfpOSR+yfUrSI5Lusv3jqlsB6FnXqCPiqxGxKyLGJO2T9JuIuK/6ZgB6wvepgWTW9PPUEfE7Sb+rsgmAIrhSA8kQNZAMUQPJEDWQDFEDyVQ5TbSzZZOe27Ol+NyFW+8oPlOS5q4tP/OmW06XH6p6p35++u/vrzL37KM3Vpk7PFf+hM7zb18oPlOSBrbNlh86uOIrtrlSA9kQNZAMUQPJEDWQDFEDyRA1kAxRA8kQNZAMUQPJEDWQDFEDyRA1kAxRA8kQNZAMUQPJEDWQDFEDyRA1kAxRA8kQNZAMUQPJVDlNtG/7vLbtmyk+9xs3HSw+U5IGvPLJjL2anL2++ExJ+sjf7q0yt9apn6/7c53TT6ffvbX4zD03l7/PStKp89urzF0JV2ogGaIGkiFqIBmiBpIhaiAZogaSIWogmUZR2x61fdD2X22fsP2O2osB6E3TF598R9LPI+Je25slDVXcCcBV6Bq17W2S3iXpE5IUEXOS5uquBaBXTR5+75Z0TtKPbD9p+yHbw5d/kO39tidsT8w//0rxRQE00yTqfklvlfS9iLhN0gVJD17+QRFxICLGI2J8YJRH50BbmkQ9LWk6Ig4tv31QS5EDWIe6Rh0R/5R02vae5XfdLWmy6lYAetb02e/PSXp4+Znvk5I+WW8lAFejUdQRcVTSeOVdABTAK8qAZIgaSIaogWSIGkiGqIFkqpwmOripozeMPFd87uGLY8VnStLkKzuKz/z11M3FZ0pS5+holbnDc1Flbo1TPyXp4s754jOfOlnnBNitpzaXH/rqytdjrtRAMkQNJEPUQDJEDSRD1EAyRA0kQ9RAMkQNJEPUQDJEDSRD1EAyRA0kQ9RAMkQNJEPUQDJEDSRD1EAyRA0kQ9RAMkQNJFPl4MGX5wb1x6mx4nN/P/Pm4jMlaWSq/N9tQ+frHOS32F9n7vm3L1SZu+fmmSpzaxwSeO2xCgcEShp8ofyfWd/cyrdxpQaSIWogGaIGkiFqIBmiBpIhaiAZogaSaRS17S/aPm77L7Z/YntL7cUA9KZr1LZ3Svq8pPGIuEVSn6R9tRcD0JumD7/7JW213S9pSNI/6q0E4Gp0jToiZiR9U9KUpDOSXoiIX17+cbb3256wPbHw4oXymwJopMnD7+sk3SNpt6QdkoZt33f5x0XEgYgYj4jxvmuHy28KoJEmD7/fI+mZiDgXEfOSHpN0R921APSqSdRTkm63PWTbku6WdKLuWgB61eRr6kOSDko6IunY8u85UHkvAD1q9PPUEfF1SV+vvAuAAnhFGZAMUQPJEDWQDFEDyRA1kEyV00RHNs/qnTeeLD732Ej5EyQl6ezg9uIzZ0erfGq1OFDnNNGBbbNV5p46X/5zK0lbT5U/+bPGqZ+SNHudi89cXOXuxZUaSIaogWSIGkiGqIFkiBpIhqiBZIgaSIaogWSIGkiGqIFkiBpIhqiBZIgaSIaogWSIGkiGqIFkiBpIhqiBZIgaSIaogWSIGkjGEeVPULR9TtLfG3zo6ySdL75APRtp3420q7Sx9l0Pu94YEa9/rRuqRN2U7YmIGG9tgTXaSPtupF2ljbXvet+Vh99AMkQNJNN21BvtH6/fSPtupF2ljbXvut611a+pAZTX9pUaQGFEDSTTWtS232v7KdtP236wrT26sX2D7d/anrR93PYDbe/UhO0+20/afrztXVZje9T2Qdt/tX3C9jva3mk1tr+4fD/4i+2f2N7S9k6XayVq232SvivpfZL2Svqo7b1t7NJAR9KXI2KvpNslfWYd73qpBySdaHuJBr4j6ecR8SZJt2od72x7p6TPSxqPiFsk9Una1+5WV2rrSv02SU9HxMmImJP0iKR7WtplVRFxJiKOLP/6JS3d6Xa2u9XqbO+S9AFJD7W9y2psb5P0Lkk/kKSImIuI59vdqqt+SVtt90sakvSPlve5QltR75R0+pK3p7XOQ5Ek22OSbpN0qN1Nuvq2pK9IWmx7kS52Szon6UfLXyo8ZHu47aVWEhEzkr4paUrSGUkvRMQv293qSjxR1pDtEUmPSvpCRLzY9j4rsf1BSf+KiMNt79JAv6S3SvpeRNwm6YKk9fz8ynVaekS5W9IOScO272t3qyu1FfWMpBsueXvX8vvWJdsDWgr64Yh4rO19urhT0odsn9LSlzV32f5xuyutaFrSdET895HPQS1Fvl69R9IzEXEuIuYlPSbpjpZ3ukJbUf9J0htt77a9WUtPNvyspV1WZdta+prvRER8q+19uomIr0bErogY09Ln9TcRse6uJpIUEf+UdNr2nuV33S1pssWVupmSdLvtoeX7xd1ah0/s9bfxP42Iju3PSvqFlp5B/GFEHG9jlwbulHS/pGO2jy6/72sR8USLO2XyOUkPL//lflLSJ1veZ0URccj2QUlHtPRdkSe1Dl8yystEgWR4ogxIhqiBZIgaSIaogWSIGkiGqIFkiBpI5j8Rnr0td0OtzAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.imshow(np.rot90(dJ_du.reshape(Nx,Ny)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now have a gradient with respect to our design parameters. To verify the accuracy of our method, we'll perform a finite difference approximation. \n",
    "\n",
    "Luckily, our solver has a built finite difference method (`calculate_fd_gradient`). Since the finite difference approximates require several expensive simulations, we'll only estimate 20 of them (randomly sampled by our function)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "db = 1e-3\n",
    "choose = 10\n",
    "g_discrete, idx = opt.calculate_fd_gradient(num_gradients=choose,db=db)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compare the results by fitting a line to the two gradients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "(m, b) = np.polyfit(np.squeeze(g_discrete), dJ_du[idx], 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and subsequently plot the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0wAAAGeCAYAAACn9IG6AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdeXhU1f3H8feXCCRABEXLWgUEcWFJMCxSQdAKVkXRWqt1CxasVTFKi4pWRWuVFjdQW3dRq4hlCSr8RK2mSKmCZAEUcEGsBHAtyBIU4/f3x9zEJCQhQCZ3MvN5Pc88zJx7Z+73TPIw88m551xzd0RERERERGRnDcIuQEREREREJFYpMImIiIiIiFRBgUlERERERKQKCkwiIiIiIiJVUGASERERERGpggKTiIiIiIhIFRSYREQqYWZTzOzWvXj+O2Y2qBZLihtmNsjM1oZdh4iISE0oMIlIzDOzNWZWZGZbzGxDEGaahV1XicrClbsf6e45dXDsv5jZJ2b2tZl9bGbXVdieZmZLzGxb8G9aNa91pJm9bGZfmdnGYP+Tot2H6pjZj8xsqpmtM7NNZvZvM+tbxb6PmZmbWecybfub2Swz2xq8P7+q5ljjg+efVaZtn6CtQy30ZbCZvR70Y00l2zsE27eZ2Uoz+2k1rzUlqKtPmbbOZlYrF1c0s7PMbGFQS04l22v8eyUiUt8pMIlIfTHM3ZsBaUA6MC7kemLFo8Bh7r4v0B8418zOADCzRsBs4O/AfsATwOygvTIvAK8ArYEfAVcAX0e3/F1qBiwGjgL2J9KHORUDs5kdAxxSyfPvB74FWgHnAn8zsyOrOd5XwM1mllQLtVe0FXgMGFvF9qlAHtASuB6YbmYHVvN6XwF7PAq6C18B9wATKm7Yg98rEZF6TYFJROoVd98AzCMSnAAws8ZmdoeZ/dfMPjWzB8wsJdh2gJm9GIyYfGVmb5hZg2Db4WaWE2x7x8xOreyYZpZpZgsqtHnwF/2LiXwRvzoYAXsh2L6mZIQgqO+eYJRkXXC/cbBtkJmtNbPfmdlnZrbezEbsxvuxyt23lmn6HigZYRkE7APc4+7fuPtkwIDjKunjAUBH4GF3/za4/dvdF5TZ5xQzyw/er4Vm1qPMtrZmNsPMPjezj8zsijLbUoIRkf+Z2btA793o32p3v8vd17t7sbs/BDQCupZ5/X2Ae4HRFfrUFPg5cIO7bwn68jxwfjWHfIlIwDqvkveoUdD/0cHjpGDE68Ya9mWRuz8FrK7ktQ8FegE3uXuRu88AlgX1V+UJoIeZHVvJ6/3CzJZUaBtjZrNrWOur7v4csK6SzYOo4e+ViEg8UGASkXrFzNoDPwM+KNM8ATiUSIjqDLQDSr7E/g5YCxxIZJThOsDNrCGREZWXiYymjAaeNrOu7IbgC/zTwF/cvZm7D6tkt+uBfkF9PYE+wB/KbG8NNA/q/jVwv5ntF/T3V2a2tLoazOxaM9sS9LMp8Eyw6UhgqbuXPU1radBe0ZdE3tO/m9lwM2tV4RjpREZHfkNkBORB4PkgDDYg8l4WBH04HrjSzIYGT7+JyOjPIcBQ4MIKr/1XM/trdX0ss28akcBU9ud/FTDf3Su+T4cC37n7e2XaCqrofwkHbgBuCn5HftjgXhKkbjGzw4FrgSTgT0Ftx5jZxpr0oxJHAqvdffNu1LoNuK3k+BU8D3QM6ixxPvBkUOsuf692UWtNf69EROo9BSYRqS+yzWwz8AnwGZEv4ZiZARcDV7n7V8EXztuAs4Pn7QDaAAe7+w53fyP4otePyOleE4LRlNeAF4FzolD7ucAt7v6Zu38O3Ez5UY4dwfYd7j4X2EIwguLuz7h7j51esQx3nwCkEhmheArYFGxqVuZ+iU3BvhVfw4HBwBrgTmC9mc03sy7BLhcDD7r7W8FIzxPAN0Tex97Age5+S/BergYe5oefwVnAn4KfzyfA5ArHvtTdL62ujwBmtm/Qv5vdfVPQ9mMiIa6yUZ5m7HxKYaX9r1DP88DnwMhKti0nchpcNvB74Hx3Lw62LXD3FrvqRxVq/LOq4EHgIDP7WYU6vwGmEYyUBachdiDyO16j36so1CoiUi8pMIlIfTHc3VOJnA50GHBA0H4g0ARYEpwqtpHIaVUlcz8mEhmNeNnMVpvZtUF7W+ATd/++zDE+JjJCUtvaBq9d9jhtyzz+0t2/K/N4G5EvpTXmEXlAEZFABpHgtW+FXfcFNlMJd1/r7pe7+yHAwUTm3DwZbD4Y+F3Jexy8zz8O+nEw0LbCtuuIjOgR7PNJmUOVfS9qJDjF8gXgTXe/vcyme4iEzYpf4GE3+1/BH4iMDCZXsu0JIn2e6+7v1+C1amKPag2C0R+DW0VPAL8K/qhwPvBcsH8otYqI1FcKTCJSr7j7v4ApwB1B0xdEQsKR7t4iuDUPFojA3Te7++/cvRNwKjDGzI4nMjfjxyXzmQIHAYWVHHYrkVAGgJm1rljWLspeR+QLdtnjVDY3pDbsww+LH7xDZI6LldneI2ivVjASdD/QLWj6hMgoUYsytybuPjXY9lGFbanuXrLC3noi4arEQbvToWC+VzaRUw5/U2Hz8cBEi6yeuCFo+49FVsN7D9inzCgZRE6JrEn/XyEStCsb+forkZGaoRZZbKI2vAN0MrOyozQ1qhV4HGgBnFG20d3fJDIfawDwKyKjc7Vhj3+vRETqIwUmEamP7gFOMLOewQjRw8DdZvYjADNrVzJ/JliooHPw5W4TUExkYYS3iIzkXG1mDS1yzaRhwLOVHK8AONIiSyknA+MrbP8U6FRNvVOBP5jZgcHiCjcSWWFsr5hZAzP7jZntZxF9gMuAfwa75BDp7xXBXKPLg/bXKnmt/czs5uC9ahDUeRHwZrDLw8AlZtY3OFZTMzs5+IK/CNhsZtcECzwkmVk3MytZ3OE5YFxwjPZUWJxhF31sCEwnEoovrDAiCJF5Sj2JzA8rWQhkGDArWAxjJpE5R03N7CfAadQ8OFwPXF2hnvOJrNiXSWQVwSeshkvcB+9rMtAw8tCSLVhZLphnlU9k7lSymZ1OJITM2NXrBqOTNwHXVLL5SeA+YEfZBTxqUGtSUOs+QIOgppI5XTnU8PdKRCQeKDCJSL0TzAN6kh/mrVxDZDTgTTP7GniVH1ZR6xI83gL8B/iru78eTOAfRmQBiS+IjBpc4O4rKznee8Atweu8D1T84vkocERwOlp2JSXfCrxNZGL8MiCXGi4HbWbnmll1f7k/HfiQyOlQfyeyWty9Qd3fAsOBC4CNRALQ8KC9om+JzHF5lci8n+VE5ihlBq/1NjCKyJfv/xF5v0u2FQOnEAksHxF5Px8hspAFRE4R/DjY9jIVAotFVjV8oIr+9Q9eewiw0SIrEW4xswHBsT9z9w0lt+A5X7h7UXD/UiCFyLy3qcBv3b1GIyHu/m8iYbCkzoOIhPULglX3niHyc7072D7AIotvVGUgkeA3l8goW1HwfpQ4G8gg8v5OAM4MftdrYiqRkbyKniIySlguoNfg9+r8oL6/ERmhKiISmnf390pEpN6z8ovciIiISLwI5n59BvSqxflWIiIJRSNMIiIi8eu3wGKFJRGRPbdP2AWIiIhI7TOzNUQuKDs85FJEROo1nZInIiIiIiJSBZ2SJyIiIiIiUgUFJhERERERkSooMImIiIiIiFRBgUlERERERKQKCkwiIiIiIiJVUGASERERERGpggKTiIiIiIhIFRSYREREREREqqDAJCIiIiIiUgUFJhERERERkSrsE3YBte2AAw7wDh06hFrD1q1badq0aag1RFM89099q7/iuX+x1LclS5Z84e4Hhl1HLImFz52yYun3JdoSqa+QWP1VX+NXrPa3us+3uAtMHTp04O233w61hpycHAYNGhRqDdEUz/1T3+qveO5fLPXNzD4Ou4ZYEwufO2XF0u9LtCVSXyGx+qu+xq9Y7W91n286JU9ERERERKQKCkwiIiIiIiJVUGASERERERGpQtzNYarMjh07WLt2Ldu3b6+T4zVv3pwVK1bUybHCUN/7l5ycTPv27WnYsGHYpYiIiIhIjEuIwLR27VpSU1Pp0KEDZhb1423evJnU1NSoHycs9bl/7s6XX37J2rVr6dixY9jliIiIiEiMS4hT8rZv307Lli3rJCxJbDMzWrZsWWejjSIiIiJSvyVEYAIUlqSUfhdEREREpKYSJjBJ3Vq3bh1nnnlm2GWIiIiIiOwVBSapdd999x1t27Zl+vTpYZciIiIiIrJXFJjqwI033sg999xT+vj6669n0qRJu3zepk2b6Nq1K6tWrQLgnHPO4eGHH95pv8WLF9O/f3969uxJnz592Lx5M9u3b2fEiBF0796d9PR0Xn/9dQCmTJnC8OHDOeGEE+jQoQP33Xcfd911F+np6fTr14+vvvoKgEGDBpGVlUVaWhrdunVj0aJFACxatIjjjz+e9PR0+vfvX1rblClTOPXUUznuuOM4/vjjWbNmDd26dQPgnXfeoU+fPqSlpdGjRw/ef/99AO666y66detGt27dSt+fNWvWcPjhhzNq1CiOPPJIhgwZQlFR0R697yIiIiIieyshVskr6+YX3uHddV/X6mse0XZfbhp2ZJXbL7roIs444wyuvPJKvv/+e5599lkWLVrE5s2bGTBgQKXPeeaZZzjiiCO47777yMzMJCsri//973+MGjWq3H7ffvstv/zlL5k2bRq9e/fm66+/JiUlhUmTJmFmLFu2jJUrVzJkyBDee+89AJYvX05eXh7bt2+nc+fO/PnPfyYvL4+rrrqKJ598kiuvvBKAbdu2kZ+fz/z587noootYvnw5hx12GPPmzWO//fbj1Vdf5brrrmPGjBkA5ObmsnTpUvbff3/WrFlTWuMDDzxAVlYW5557Lt9++y3FxcUsWbKExx9/nLfeegt3p2/fvhx77LHst99+vP/++0ydOpWHH36Ys846ixkzZnDeeeftzY9IRERERGSPJFxgCkOHDh1o2bIleXl5fPrpp6Snp9OyZUsA8vPzq33uCSecwD/+8Q8uu+wyCgoKdtq+atUq2rRpQ+/evQHYd999AViwYAGjR48G4LDDDuPggw8uDUyDBw8mNTWV1NRUmjdvzrBhwwDo3r07S5cuLX3tc845B4CBAwfy9ddfs3HjRjZv3syll17KRx99hJmxY8eOcrXuv//+O9V49NFH86c//Ym1a9dyxhln0KVLFxYsWMDpp59O06ZNATjjjDN44403OPXUU+nYsSNpaWkAHHXUUeXCl4jsmTmr5zApdxIbtm6gddPWZPXK4uROJ4ddlojsQnZeIRPnrWLdxiL2TzZuaF7I8PR2YZclklASLjBVNxIUTSNHjmTKlCls2LCBiy66CKBGI0zff/89K1asoEmTJvzvf/+jffv2e11L48aNS+83aNCg9HGDBg347rvvSrdVXE3OzLjhhhsYMGAAL7zwAmvWrGHQoEGl20vCT0W/+tWv6Nu3L3PmzOGkk07iwQcfrHF9SUlJOiVPZC/NWT2H8QvHs704spz++q3rGb9wPIBCk0gMy84rZNzMZRTtKAbgy+3OuJnLABSaROqQ5jDVkdNPP52XXnqJxYsXM3ToUABSU1PJz8+v9HbEEUcAcPfdd3P44YfzzDPPMGLEiHIjOgBdu3Zl/fr1LF68GIiEsO+++44BAwbw9NNPA/Dee+/x3//+l65du+5WzdOmTQMio1XNmzenefPmbNq0ibZt2wKReUs1sXr1ajp16sQVV1zBaaedxtKlSxkwYADZ2dls27aNrVu3MmvWrCrDo4jsnUm5k0rDUontxduZlLvruZQiEp6J81aVhqUSRTuKmThvVUgViSSmhBthCkujRo0YPHgwLVq0ICkpqUbPWbVqFY888giLFi0iNTWVgQMHcuutt3LzzTeXe91p06YxevRoioqKSElJ4dVXX+XSSy/lt7/9Ld27d2efffZhypQp5UZuaiI5OZn09HR27NjBY489BsDVV1/N+eefz5133snJJ9fsL9PPPfccTz31FA0bNqR169Zcd9117L///mRmZtKnTx8gMgKXnp6u0+9EomDD1g271S4isWHdxsrPsKiqXUSiQ4Gpjnz//fe8+eab/OMf/6jxc7p27cqKFStKH991112V7te7d2/efPPNndoff/zxndoyMzPJzMwsfVw2oFTcdt5555Vb3Q8i85Hy8vJITU0F4NZbb630uR06dGD58uUAXHvttVx77bU71TJmzBjGjBlTrq3s8wB+//vf7/Q8Edk9LS2FL3zbTu2tm7YOoRoRqam2LVIorCQctW2REkI1IolLp+TVgXfffZfOnTtz/PHH06VLl7DLEZEEsuwff+L3n/6XRt+Xn5OYnJRMVq+skKoSkZoYO7QrKQ3Ln5WS0jCJsUN37xR7Edk7GmGqA0cccQSrV68Ou4zdkpOTE3YJIrKX5r65lGOW38fWlKO4/uhf88Dyv2mVPJF6pGRhh3Kr5J3WXQs+iNQxBSYRkTiUnVfImNmfMPzHk7l1xCn0S0nhjMNOD7ssEdlNw9PblQaknJwcBiksidQ5BSYRkTizbPptrMr/mL4dR3Br5lCaNNJ/9SIiIntKn6IiInFk2fTb6L78z2xtNoDRF/ZSWBIREdlLWvRBRCROlISlN5MH0CPrHzRp3CjskkREROo9BaY6kpSURFpaGj179qRXr14sXLhwj17nnnvuYdu2nZcHBnjjjTc48sgjSUtLo7CwkDPPPBOA/Px85s6du8e1i0js2ykspWjZYRERkdqgwFRHUlJSyM/Pp6CggNtvv51x48bt0etUF5iefvppxo0bR35+Pu3atWP69OmAApNIvMvOK+SZvC8VlkRERKJAgSkEX3/9Nfvtt1/p44kTJ9K7d2969OjBTTfdBMDWrVs5+eST6dmzJ926dWPatGlMnjyZdevWMXjwYAYPHlzuNR955BGee+45brjhBs4991zWrFlDt27d+Pbbb7nxxhuZNm0aaWlpTJs2rU77KiLR9dLCJYx5Lp81B/+CHlfNUlgSERGpZYk5G/jxSq49cuRw6DMKvt0GT/9i5+1pv4L0c2Hrl/DcBeW3jZizy0MWFRWRlpbG9u3bWb9+Pa+99hoAL7/8Mu+//z6LFi3C3Tn11FOZP38+n3/+OW3btmXOnMhrb9q0iebNm3PXXXfx+uuvc8ABB5R7/ZEjR7JgwQJOOeUUzjzzTNasWQNAo0aNuOWWW3j77be57777dv3eiEi9sWz6bRy77G5+2f5ubtBqeCIiIlGhEaY6UnJK3sqVK3nppZe44IILcHdefvllXn75ZdLT0+nVqxcrV67k/fffp3v37rzyyitcc801vPHGGzRv3jzsLohIDCmZs1SQ0pcbLvq5wpKIiEiUJOYnbHUjQo2aVL+9acsajShV5+ijj+aLL77g888/x90ZN24cv/nNb3baLzc3l7lz5/KHP/yB448/nhtvvHGvjisi8UELPIiIiNQdjTCFYOXKlRQXF9OyZUuGDh3KY489xpYtWwAoLCzks88+Y926dTRp0oTzzjuPsWPHkpubC0BqaiqbN2/erePtyXNEJDYtnDdNYUlERKQOhTrCZGaPAacAn7l7t0q2DwJmAx8FTTPd/Za6q7D2lMxhAnB3nnjiCZKSkhgyZAgrVqzg6KOPBqBZs2b8/e9/54MPPmDs2LE0aNCAhg0b8re//Q2Aiy++mBNPPJG2bdvy+uuv1+jYgwcPZsKECaSlpTFu3Dh++ctfRqeTIhJV2XmF/D6nCeN+NJpzLrleYUlERKQOhH1K3hTgPuDJavZ5w91PqZtyoqe4uLjKbVlZWWRlZZVrO+SQQxg6dOhO+44ePZrRo0dX+jpTpkwpvd+hQweWL18OwP7778/ixYv3oGoRiRV52ZO5463m9O7YlXMyx2vOUgwws2HAsM6dO4ddioiIRFGop+S5+3zgqzBrEBGJdTveySY9/wau3/81Hs3MUFiKEe7+grtfrEV5RETiW32Yw3S0mRWY2f+Z2ZFhFyMiUpeWTb+dEz5/nLeSj+HYyx9QWBIREaljsf7Jmwsc7O5bzOwkIBvoUnEnM7sYuBigVatW5OTklNvevHnzOl30oLi4OK4XWYiH/m3fvn2n3xOALVu2VNoeD+K5bxCf/dvxTjYnfP44bzToTVGvK1n01lthlyQiIpJwYjowufvXZe7PNbO/mtkB7v5Fhf0eAh4CyMjI8EGDBpV7nRUrVpCamloHFUds3ry5To9X1+Khf8nJyaSnp+/UnpOTQ8Xfn3gRz32D+Ovf7CUf0enT13mryTEU9bqSIUNOCLskERGRhBTTp+SZWWszs+B+HyL1fhluVSIi0ZWd+wlXTX+Xe9reSfes6TRq1DDskkRERBJW2MuKTwUGAQeY2VrgJqAhgLs/AJwJ/NbMvgOKgLPd3UMqV0Qk6pZOn0Dzgv/jJx3+yL0j+mvOkoiISMjCXiXvHHdv4+4N3b29uz/q7g8EYQl3v8/dj3T3nu7ez90Xhlnv3srOzsbMWLlyZZX7ZGZmMn36dABGjhzJu+++W+1rnnTSSWzcuLHafaZMmcK6det2v+AQ9e/fP+wSROrc0ukT6LH8dpo0bcaDF/RWWBIREYkBMX1KXljmrJ7DkOlD6PFED4ZMH8Kc1XNq5XWnTp3KMcccw9SpU2u0/yOPPMIRRxxR7T5z586lRYsW1e5TnwLTd999B8DChfU6G4vstpKw9FbyMXTPmq6L0oqIiMQIBaYK5qyew/iF41m/dT2Os37resYvHL/XoWnLli0sWLCARx99lGeffba03d25/PLL6dq1Kz/96U/57LPPSrcNGjSIt99+G4iEre7du9OtWzeuueaa0n06dOjAF198wZo1azj88MMZNWoURx55JEOGDKGoqIjp06fz9ttvc+6555KWlkZRUVG5uj744AN++tOf0rNnT3r16sWHH36IuzN27Fi6detG9+7dmTZtGhCZVH/sscdy9tln06lTJ6699lqefvpp+vTpQ/fu3fnwww+ByCjZJZdcQkZGBoceeigvvvgiAGvWrGHAgAH06tWLXr16lYainJwcBgwYwKmnnloaEJs1awbA+vXrGThwIGlpaXTr1o033nij2vejWbNmXH/99fTs2ZN+/frx6aef7tXPTaQu5M+8U2FJREQkRikwVTApdxLbi7eXa9tevJ1JuZP26nVnz57NiSeeyKGHHkrLli1ZsmQJALNmzWLVqlW8++67PPnkk5WOrKxbt45rrrmG1157jfz8fBYvXkx2dvZO+73//vtcdtllvPPOO7Ro0YIZM2Zw5plnkpGRwdNPP01+fj4pFb6InXvuuVx22WUUFBSwcOFC2rRpw8yZM8nPz6egoIBXX32VsWPHsn79egAKCgq45557WLFiBU899RTvvfceixYtYuTIkdx7772lr7tmzRoWLVrEnDlzuOSSS9i+fTs/+tGPeOWVV8jNzWXatGlcccUVpfvn5uYyadIk3nvvvXL1PfPMMwwdOrS0nrS0tGrfj61bt9KvXz8KCgoYOHAgDz/88B7+xETqRnZeIX94O5nXmpyosCQiIhKDFJgq2LB1w26119TUqVM5++yzATj77LNLT8ubP38+55xzDklJSbRt25bjjjtup+cuXryYQYMGceCBB7LPPvtw7rnnMn/+/J3269ixI2lpaQAcddRRrFmzptqaNm/eTGFhIaeffjoQWWq7SZMmLFiwoLSmVq1aceyxx7J48WIAevfuTevWrWncuDGHHHIIQ4YMAaB79+7ljnfWWWfRoEEDunTpQqdOnVi5ciU7duxg1KhRdO/enV/84hfl5mf16dOHjh077lRj7969efzxxxk/fjzLli0jNTW12vejUaNGnHLKKTV+D0TC9K9/zmXMc/mkdjiKflc+rbAkIiISgzSjuILWTVuzfuv6Stv31FdffcVrr73GsmXLMDOKi4sxMyZOnLg3pe6kcePGpfeTkpJ2Ov2uto/RoEGD0scNGjQonX8EEKwGX+7x3XffTatWrSgoKOD7778nOTm5dHvTpk0rPd7AgQOZP38+c+bMITMzkzFjxtC8efMq62vYsGHpsZOSksrVJBK2OavnMCl3Ehu2bqClpfD7T//LFa2v4uLMMVrgQUREJEZphKmCrF5ZJCcll2tLTkomq1fWHr/m9OnTOf/88/n4449Zs2YNn3zyCR07duSNN95g4MCBTJs2jeLiYtavX8/rr7++0/P79OnDv/71L7744guKi4uZOnUqxx57bI2Pn5qayubNmyttb9++fenpbN988w3btm1jwIABpTV9/vnnzJ8/nz59+uxWn//xj3/w/fff8+GHH7J69Wq6du3Kpk2baNOmDQ0aNOCpp56iuLh4l6/z8ccf06pVK0aNGsXIkSPJzc3d6/dDJAwV50d+4du48YADOWhoF4UlERGRGKbAVMHJnU5mfP/xtGnaBsNo07QN4/uP5+ROJ+/xa06dOrX0tLcSP//5z0vbu3TpwhFHHMEFF1zA0UcfXW4/M6NNmzZMmDCBwYMH07NnT4466ihOO+20Gh+/ZBGGyhZ9eOqpp5g8eTI9evSgf//+bNiwgdNPP50ePXrQs2dPjjvuOP7yl7/QuvXujbAddNBB9OnTh5/97Gc88MADJCcnc+mll/LEE0/Qs2dPVq5cWeWoUlk5OTn07NmT9PR0pk2bRlZW1l6/HyJhqGx+5LcNnAeW/y2kikRERKQmLN6uA5uRkeElK8uVWLFiBYcffnid1bB582ZSU1P3+nW6d+/O888/X+ncnjDtqn+ZmZmccsopnHnmmXVY1e6p6nciJyeHQYMG1X1BdSCe+wax378eT/TA2fn/W8NYeuHSap8bS30zsyXunhF2HbGkss+dMMXS70u0JVJfIbH6q77Gr1jtb3WfbxphilEnnHAC3bt3j7mwJCJ7Zt+GB1bavjfzI0VERCT6dOJ8jHrllVfCLmGPTZkyJewSRGLK0hl/Zr+PD+HrNhtx+7a0fW/nR4qIiEj0aYRJRCSKlk6fQI9lt3Ftg08Z3/+mWp0fKSIiItGXMCNM7r7TUteSmOJt3p7ErqXTJ9Bj+e0sSv4J6Vc8S5OUFM449NSwyxIREZHdkBAjTMnJyd7EoOEAACAASURBVHz55Zf6oiy4O19++WW5a0CJREPZsNQta4YuSisiIlJPJcQIU/v27Vm7di2ff/55nRxv+/btcf2FvL73Lzk5mfbt24ddhsSx7Ny1bCtYxPamCksiIiL1XUIEpoYNG9bpanM5OTmkp6fX2fHqWrz3T2RvvLD4PcbMfJ9+HX7PIxekKyyJiIjUcwlxSp6ISF1YOn0CvV78GSce5Dwyoo/CkoiISBxQYBIRqQUlc5YKU7pyR+bxNGmUEAP4IiIicU+BSURkL2mBBxERkfilwCQishcWv/CwwpKIiEgcU2ASEdlD2XmFjFzYnFnNzlFYEhERiVM6yV5EZA+8NWcK4xY0J63jQQzNvE9zlkREROKURphERHbT0ukT6Ls4i/EHvM6jmRkKSyIiInFMgUlEZDeUXeBh2KUTFJZERETinAKTiEgNlV8NbzpNUpqEXZKIiIhEmQKTiEgNzHnrXdouu19hSUREJMEoMImI7EJ2XiGjsz/ij60mKSyJiIgkGJ18LyJSjYIZf6Ewbzl9O/yW20f01pwlERGRBKNPfhGRKhTM+As9l/2Jb5r9hBEXpissiYiIJCCdkiciUomSsFQ6Zyk5OeySREREJAQKTCIiFewUljRnSUREJGEpMImIlJGdV8hjuZt4M3mAwpKIiIgoMImIlJj370WMeS6fzw4eRo+rZiksiYiIiAKTiAhETsMb9PLPOLfdZzyamUGTxg3DLklERERigAKTiCS8kjlLS1P6MO7Xv9RqeCIiIlJKgUlEElpJWFqc3J8js2boNDwREREpR4FJRBLWG69kKyyJiIhItRSYRCQhZecVcuFr+/BgizEKSyIiIlIlnagvIgknL3sykxc1oW/Hbpyfeb3mLImIiEiVNMIkIgmlYMZfSM+/geta/DOyGp7CkoiIiFRDgUlEEsYPCzwcTf/RjyosiYiIyC4pMIlIQigblo7Mmqk5SyIiIlIjCkwiEvdmL/mYbwumKyyJiIjIbtP5KCIS17Jz/8uY6cs59uAJ3H9hP4UlERER2S0aYRKRuFUwYyKtsn/JwIObcv9FAxWWREREZLcpMIlIXCqYMZGey25lnybN+euFfbTAg4iIiOwRBSYRiTslYUlzlkRERGRv6U+uIhJX8rInk66wJCJRkJ1XyMR5q1i3sYi2LVIYO7Qrw9PbhV2WiERZqCNMZvaYmX1mZsur2G5mNtnMPjCzpWbWq65rFJH6IzuvkLGLknmtyYkKSyJSq7LzChk3cxmFG4twoHBjEeNmLiM7rzDs0kQkysI+JW8KcGI1238GdAluFwN/q4OaRKQe+u+HyxnzXB4HduhBvyufVlgSkVo1cd4qinYUl2sr2lHMxHmrQqpIROpKqIHJ3ecDX1Wzy2nAkx7xJtDCzNrUTXUiUl8UzJjIBZ9cz9hWuTyamaEFHkSk1q3bWLRb7SISP8IeYdqVdsAnZR6vDdpERIAfFnj4T9JRXPib3yssiUhUtG2RslvtIhI/4uKbhZldTOSUPVq1akVOTk6o9WzZsiX0GqIpnvunvtUv37z7AkM/e4T/JB3F+m5X8M1bi8IuKSri8WcnUt+MHdqVcTOXlTstL6VhEmOHdg2xKhGpC7EemAqBH5d53D5oK8fdHwIeAsjIyPBBgwbVSXFVycnJIewaoime+6e+1R8vL3iL4z59nLdT+tEzaxbfvLUorvpXVrz97ETqo5LV8LRKnkjiifXA9DxwuZk9C/QFNrn7+pBrEpGQZecVMmbOF1zQ9s9cPfJcLfAgInVieHo7BSSRBBRqYDKzqcAg4AAzWwvcBDQEcPcHgLnAScAHwDZgRDiVikisyJ95Jy8v+Zq+HX/G1ZlDNWdJREREoirUbxrufs4utjtwWR2VIyIxrmDGRNKW3cqFqcfQ/cI/KCyJiIhI1MX6KnkiIsAPq+G9ndyP7ln/oEnjhmGXJHHKzDqZ2aNmNj3sWkREJHwKTCIS88qGpSOyZmnOklTJzB4zs8/MbHmF9hPNbJWZfWBm11b3Gu6+2t1/Hd1KRUSkvtD5LCIS07LzCvk0L58dzRSWpEamAPcBT5Y0mFkScD9wApHr+S02s+eBJOD2Cs+/yN0/q5tSRUTqr+y8woRZNVKBSURi1ouLVjJm1of07XA551+YTpPk5LBLkhjn7vPNrEOF5j7AB+6+GiBYefU0d78dOKVuKxQRqf+y8wrLXZescGMR42YuA4jL0KTAJCIxqWDmHfQpuJdTfjyJCSN6a4EH2RvtgE/KPF5L5FIVlTKzlsCfgHQzGxcEq4r7xNQF08tKpAsdJ1JfIbH6q77Gtj/mbKNoh5drK9pRzB9nF9Bi0/vVPrc+9lffQEQk5hTMvIOeS//I2yn9mDBiiMKS1Cl3/xK4ZBf7xNQF08tKpAsdJ1JfIbH6q77Gtq9emlN5+3bfZV/qY3+16IOIxJTSsKQFHqT2FAI/LvO4fdAmIiJ7oG2LlN1qr+8UmEQkZrw15wmFJYmGxUAXM+toZo2As4HnQ65JRKTeGju0KykNk8q1pTRMYuzQriFVFF0KTCISE7LzCrloQTOmNztPYUn2mJlNBf4DdDWztWb2a3f/DrgcmAesAJ5z93fCrFNEpD4bnt6O28/oTrsWKRjQrkUKt5/RPS4XfADNYRKRGLBozmPcsGBfenRsz0mZkzRnSfaYu59TRftcYG4dlyMiEreGp7eL24BUkUaYRCRUBTPvoM/iq7il5as8mpmhsCQiIiIxRYFJREJTssDDksb9GHrZXQpLIiIiEnMUmEQkFGXD0uFXas6SiIiIxCYFJhGpcy8ueo8DC+5XWBIREZGYp8AkInUqO6+QK2a9z22tJiksSb1mZsPM7KFNmzaFXYqIiESRApOI1Jn8mXeyccZV9O2wP38ZeZLCktRr7v6Cu1/cvHnzsEsREZEo0gxrEakT+TPvJG3pLRQ37cdZF6RpgQcRERGpFzTCJCJRVxKWljTux+FZM2mSkhJ2SSIiIiI1osAkIlGVP+uu8mGpSdOwSxIRERGpMQUmEYma7LxC7nt7C28lH6OwJCIiIvWSApOIRMXLC/7DmOfy2XrwELqPma2wJCIiIvWSApOI1Lr8mXdy3CsnMaLtJzyamaEFHkRERKTe0rcYEalVpQs8JPfjd78+X2FJRGJOdl4hE+etYt3GItq2SGHs0K4MT28XdlkiEqP0TUZEas1Oq+HpNDwRiTHZeYWMm7mMoh3FABRuLGLczGUACk0iUimdkicitSLntf9TWBKRmDdx3qrSsFSiaEcxE+etCqkiEYl1Ckwistey8wq56JVi7m9xtcKSiMS0dRuLdqtdRESn5InIXsnNnswDixrRt2MaIzKv0ZwlSRhmNgwY1rlz57BLkd3QtkUKhZWEo7YtdEFtEamcRphEZI/lz7yTXvk3MK75q1oNTxKOu7/g7hc3b9487FJkN4wd2pWUhknl2lIaJjF2aNeQKhKRWKdvNyKyR35Y4KEvvUc/obAkIvVCycIOWiVPRGpK33BEZLeVDUuHZ83SnCURqVeGp7dTQBKRGtMpeSKyW7JzP+Hr/FkKSyIiIpIQNMIkIjU2e8nHjJm+nGM6/IkHzs9QWBIREZG4pxEmEamR/Fl3cfDs0xl0cDIPjPiJwpKIiIgkBI0wiUil5qyew6TcSWzYuoH9rQljP/2Y1indue+CflrgQURERBKGRphEZCdzVs9h/MLxrN+6Hsf50rdy4wEH8PHw32hkSaLKzI42s/vNbKmZfW5m/zWzuWZ2mZlp/W4REalzCkwispNJuZPYXry9XNu3DeCBdx4MqSJJBGb2f8BIYB5wItAGOAL4A5AMzDazU8OrUEREEtEuz6sxs6fc/fxdtYlI/NiwdcNutYvUkvPd/YsKbVuA3OB2p5kdUPdliYhIIqvJCNORZR+YWRJwVHTKEZFYsF/SvpW2t27auo4rkURSSVjCzI43s2Fm1rCqfURERKKpysBkZuPMbDPQw8y+Dm6bgc+A2XVWoYjUqfxZd3H1+g/Zx5PKtScnJZPVKyukqiQRmdmdwE+AnuhzR0REQlJlYHL32909FZjo7vsGt1R3b+nu4+qwRhGpI/mz7iKt4GbafNeNG/reQJumbTCMNk3bML7/eE7udHLYJUocM7M7zaxFmaaDgD8Cfwrux5Rg5OuhTZs2hV2KiIhE0S7nMLn7ODNrBxxcdn93nx/NwkSkbpWEpdzGfTgsK5teTZpyxuE/D7ssSSwzgWfNbC5wP/Ak8DqRBR8eDrOwyrj7C8ALGRkZo8KuRUREoqcmiz5MAM4G3gWKg2YHFJhE4sRLC3MZlH8bucmRsKSlwyUM7v5v4EQzO4/ISnmT3X1QuFWJiEiiq8nVJ08Hurr7N9EuRkTqXnZeIWNeWM+v2t3BdRf9QmFJQmNm+wBDicyVHQ5cZWYjgRvcvSDU4kREJGHVJDCtBhoCCkwicSZv1t0seHsDfTueznWZQ2nSqCb/JYhETTbwH6AJcK67X2hmbYFbzMzdXae+iYhInavJt6NtQL6Z/ZMyocndr4haVSISdXmz7ia9YDzfN+vH4RfeqrAkseBgdz/FzBoBbwK4+zpgpJmlhVuaiIgkqpp8Q3o+uIlInCgJS3mN+3B41kyaNG4YdkkiAA+Z2X+C+3eV3eDu+SHUIyIiUqNV8p4wsxTgIHdfVZsHN7MTgUlAEvCIu0+osD0TmAgUBk33ufsjtVmDSKIpG5a6aoEHiSHufi9wb9h1iIiIlFXldZhKmNkwIB94KXicZmZ7PeJkZklElo39GXAEcI6ZHVHJrtPcPS24KSyJ7IXsvEJeX7JcYUlikpn9wcz2q2b7cWZ2Sl3WJCIiUpNT8sYDfYAciJwWYWadauHYfYAP3H01gJk9C5xGZPlyEalluZ9s5L55+fTtOIpLLkijSXJy2CWJVLQMeNHMtgO5wOdErsHUBUgDXgVuC688ERFJRLscYQJ2uHvFy5h/XwvHbgd8Uubx2qCtop+b2VIzm25mP66F44oknLxZd5P5wWWc3n4rj2ZmKCxJTHL32e7+E+AS4B0ip2t/Dfwd6OPuV7n752HWKCIiiacmI0zvmNmvgCQz6wJcASyMblmlXgCmuvs3ZvYb4AnguIo7mdnFwMUArVq1Iicnp47Kq9yWLVtCryGa4rl/8di3be/O5aTPHuStpDSGdmnBooULwi4pKuLxZ1cinvtWGXd/H3g/7DpERESgZoFpNHA9kSXFpxK5+vofa+HYhUDZEaP2/LC4AwDu/mWZh48Af6nshdz9IeAhgIyMDB80aFAtlLfncnJyCLuGaIrn/sVb3/Jm3c2gzx4kr3EfNh01liFDhoRdUtTE28+urHjum8S+7LxCJs5bxbqNRbRtkcLYoV0Znl7ZCSEiIvFpl6fkufs2d7/e3Xu7e0Zwf3stHHsx0MXMOgbX3DibCsuXm1mbMg9PBVbUwnFFEsJ/Xnqm3Gp4jRo1CrskEalnsvMKGTdzGYUbi3CgcGMR42YuIzuvcJfPFRGJF1UGJjO7J/j3BTN7vuJtbw/s7t8BlxMZsVoBPOfu75jZLWZ2arDbFWb2jpkVEDkVMHNvjyuSCLLzChnxryZMS71Qq+FJvWJmSWZ2Vdh1SMTEeaso2lFcrq1oRzET59XqVUZERGJadafkPRX8e0e0Du7uc4G5FdpuLHN/HDAuWscXiUeLXnyUm//djPSOBzMs8y6aNKrJmbciscHdi83sHODusGvZleCyG8M6d+4cdilRs25j0W61i4jEoyq/Sbn7kuDff9VdOSKyN/Jm3U2fgvHcuv8ZDM58WGFJ6qt/m9l9wDRga0mju+eGV9LO3P0F4IWMjIxRYdcSLW1bpFBYSThq2yIlhGpERMJR5bcpM1sGeFXb3b1HVCoSkT2SN+vu0jlLgy+7V2FJ6rO04N9byrQ5laySKtE1dmhXxs1cVu60vJSGSYwd2jXEqkRE6lZ136hKrqZ+WfBvySl651FNkBKRulc2LHXNmkWTJs3CLklkj7n74LBrkIiS1fC0Sp6IJLLqTsn7GMDMTnD39DKbrjGzXODaaBcnIrv2/Nsf0iP/b+QlKyxJfDCz5sBNwMCg6V/ALZVcRF3qwPD0dgpIIpLQdrmsOGBm9pMyD/rX8HkiEmXZuWu5csZK/tL6boUliSePAZuBs4Lb18DjoVYkIiIJqyaTHH4NPBb8xc+A/wEXRbUqEdmlvOxJ7FjyT/p1uIY7RvTVnCWJJ4e4+8/LPL7ZzPJDq0ZERBLaLr9hBavl9QwCEzolQiR8edmTSM+/EZr04ZHz0xSWJN4Umdkx7r4AIDjLQetYi4hIKGr0LcvMTgaOBJLNDAB3v6XaJ4lIVJSEpR8WeNBFaSXuXAI8WfKHOiJnNlwYYj0iIpLAdhmYzOwBoAkwGHgEOBNYFOW6RKQSedmTK4QlzVmS+GJmDYCu7t7TzPYFcPevQy5LREQSWE0Wb+jv7hcA/3P3m4GjgUOjW5aIVJSdV8gdi7bzVvIxCksSt9z9e+Dq4P7XCksiIhK2mgSm7cG/28ysLbADaBO9kkSkolfnL2DMc/l832Eg3cfMVliSePeqmf3ezH5sZvuX3MIuSkREElNN5jC9YGYtgIlALpGL1j4c1apEpFRe9iQG593Eb9rcwujMoVrgQRLBL4N/LyvT5kCnEGoREZEEV+03r+Bc8n+6+0Zghpm9CCRrpTyRulG6wENyb0aPvEhhSeJe8Llznrv/O+xaREREYBen5AXnkt9f5vE3CksideOH1fB60zUrW6fhSUIIPnfuC7sOERGREjWZw/RPM/u5lawnLiJR91rOawpLksj0uSMiIjGjJuf3/AYYA3xnZtsBA9zd941qZSIJKjuvkDHzirii9fVcPOpShSVJRPrcEZHdlp1XyMR5q1i3sYi2LVIYO7Qrw9PbhV2WxIFdBiZ3T62LQkQEcmffx5S3iunbsS8XZ47RnCVJSPrcEZHdlZ1XyLiZyyjaUQxA4cYixs1cBqDQJHutylPyzCzJzJqVedzPzAYGN32YidSyvOzJ9Mq7nqubv8ajmRkKS5JwzOy8Mvd/UmHb5XVfUfXMbJiZPbRpk6b2ioRt4rxVpWGpRNGOYibOWxVSRRJPqpvD9Gfg0jKPpwJjgRuAP0SzKJFEk5c9mfT8G8hr3Ju00U8rLEmiGlPm/r0Vtl1Ul4XUhLu/4O4XN2/ePOxSRBLeuo1Fu9Uusjuq+1Z2PNC7zOON7j4smIT7RnTLEkkcZcOSFniQBGdV3K/ssYhIqbYtUiisJBy1bZESQjUSb6obYWrg7t+VeXwNRGbdAvpGJ1ILsnPX8tmS5xWWRCK8ivuVPRYRKTV2aFdSGiaVa0tpmMTYoV1DqkjiSXUjTI3MLNXdNwO4+8sAZtYcSK6L4kTi2ewlHzFm+rv073ATD52frrAkAoeZ2VIio0mHBPcJHncKrywRiXUlCztolTyJhuoC08PANDO7xN3/C2BmBwN/Ax6pi+JE4lVe9mQOzX2Ynx50B/eMOFpzlkQiDg+7ABGpv4ant1NAkqio8luau99lZtuABWbWNGjeAkxw97/VSXUicahkzlJ+Sgb3XPgThSWRgLt/HHYNIiIiFVX7Tc3dHwAeKFlGvOT0PBHZM6VhqXEGh2bN1ml4IiIiIjGuRn/aVlAS2XuLX3yE3gpLIiIiIvVKdavkiUgtyc4rZPS/k3mlySkKSyI1YGYpZqblrUREJHS7HGEys8bu/s2u2kSkcgtensHY1xuS0bEzP8l8QnOWRHbBzIYBdwCNgI5mlgbc4u6nhluZiIjEkuy8wjpZGbEmI0z/qWGbiFSQlz2ZYxZexE0H/otHMzMUlkRqZjzQB9gI4O75QMcwCxIRkdiSnVfIuJnLKNxYhAOFG4sYN3MZ2XmFtX6sKr+9mVlroB2QYmbp/HCV9X2BJrVeiUicKVngoaDxUZxxyc0KSyI1t8PdN5lZ2TZduFZEREpNnLeKoh3F5dqKdhQzcd6qWh9lqu4b3FAgE2gP3FWmfTNwXa1WIRJnyoalLlnPa86SyO55x8x+BSSZWRfgCmBhyDWJiEgMWbexaLfa90Z112F6AnjCzH7u7jNq/cgicWrum8sYkHcbBckKSyJ7aDRwPfAN8AwwD7g11IpERCSmtG2RQmEl4ahti5RaP1ZNzhF6MfhLX4ey+7v7LbVejUg9l51XyJjZ/+WX7e/ghhHDFZZE9sxh7n49kdAkIiKyk7FDuzJu5rJyp+WlNExi7NDaX2C1JoFpNrAJWELkr30iUonc7MnkLf6Qvh3P4YbMoZqzJLLn7gzm0U4Hprn78rALEhGR2FIyT6kuVsmryTe69u5+Yq0fWSSO5GZPJi3vRho0O4pDL+ylsCSyF9x9cBCYzgIeNLN9iQQnnZYnIiKlhqe3i0pAqqgmy4ovNLPuUa9EpJ4qCUvLkntFLkrbuFHYJYnUe+6+wd0nA5cA+cCNIZckIiIJqiZ/Bj8GyDSzj4ickmeAu3uPqFYmUg+UDUta4EGkdpjZ4cAvgZ8DXwLTgN+FWpSIiCSsmgSmn0W9CpF6KDuvkGVvv0eDpkdFRpYUlkRqy2NEQtJQd18XdjEiIpLYqrtw7b7u/jWR6y6JSBlz31zOmNkf07fjefzugnSaJDcOuySRuOHuR4ddg4iISInqRpieAU4hsjqeEzkVr4QDnaJYl0jMys2ezIC82/hF+zu5SavhidQaM3vO3c8ys2VEPmdKNxGDp4Kb2TBgWOfOncMuRUREoqi6C9eeEvzbse7KEYltZecs3TTiNIUlkdqVFfx7SqhV1JC7vwC8kJGRMSrsWkREJHpqskoeZnaqmd0R3OrFB5lIbdMCDyLR5e7rg7uXuvvHZW/ApWHWJiIiiWuXgcnMJhD5q9+7wS3LzG6LdmEisWTBy9MVlkTqzgmVtGkBIhERCUVNzic6CUhz9+8BzOwJIA+4LpqFicSK7LxCxr7eiJsPGMHw3/5RYUkkSszst0RGkjqZ2dIym1KBf4dTlYiIJLqaTsBoAXwV3G8epVpEYs6iFx/ltn8nk9HxEIZnTtScJZHoegb4P+B24Noy7Zvd/avKnyIiIhJdNfn2dzuQZ2avE1mpaCDlP8hE4lJu9mQy8m7k1v1O4ZjMKQpLIlHm7puATcA5AGb2IyAZaGZmzdz9v2HWJyIiiWmXc5jcfSrQD5gJzACOdvdptXFwMzvRzFaZ2QdmtlMIM7PGZjYt2P6WmXWojeOK7ErZBR6OufwBhSWROmRmw8zsfeAj4F/AGiIjTyIiInWuysBkZocF//YC2gBrg1tbM0s3s4P35sBmlgTcT2Qi7xHAOWZ2RIXdfg38z907A3cDf96bY4rUhFbDEwndrUT+UPdecGmL44E3wy1JREQSVXV/Nv8dMAq4s4rtLc2swN3P38Nj9wE+cPfVAGb2LHAakZX4SpwGjA/uTwfuMzNz97IXNBTZK3NWz2FS7iQ2bN1ACs0Z+en/SFJYEgnTDnf/0swamFkDd3/dzO4JuygREUlM1V24dlTw7+Cq9jGzl/fi2O2AT8o8Xgv0rWofd//OzDYBLYEv9uK4IqXmrJ7D+IXj2V68HYBtbOS+HzXipr6X0VNhSSQsG82sGTAfeNrMPgO2hlyTiIgkqCoDk5mdUd0T3X2muw+p/ZJ2n5ldDFwM0KpVK3JyckKtZ8uWLaHXEE3x1L8/r/1zaVgq8b3t4J4lk9j/05YhVRUd8fRzq0w89y+e+1aF04DtwFXAuURWZ70l1IpERCRhVXdK3rDg3x8B/YHXgseDgYVEFoHYG4XAj8s8bh+0VbbPWjPbh8iH5pcVX8jdHwIeAsjIyPBBgwbtZWl7Jycnh7BriKZ46t/GJzZW3l68MW76WCKefm6Vief+xXPfKuPuZUeTngitEBEREao/JW8ElJ52d4S7rw8etwGm1MKxFwNdzKwjkWB0NvCrCvs8D1wI/Ac4E3hN85ekNu1vTfnSt+zU3rpp6xCqEUlsZrYZKPt/vAWPDXB33zeUwkREJKHtcllx4MclYSnwKXDQ3h7Y3b8DLgfmASuA59z9HTO7xcxODXZ7lMjiEh8AY9D1n6QWLZn9V3634WMafV++PTkpmaxeWeEUJZLA3D3V3fctc0st+2/Y9YmISGKqycVl/mlm84CpweOzgVdr4+DuPheYW6HtxjL3twO/qI1jiZSVnVfIk299y9X79uL6fpk88M5DbNi6gRZJLbim/zWc3OnksEsUSWhmdgzQxd0fN7MDgFR3/yjsukREJPHsMjC5++VmdjowMGh60N1nRbcskej5579yGPPSVvp27EePzMvp12gfzjj8TCCYK9JpUKj1iSQ6M7sJyAC6Ao8DjYC/Az8Jsy4REUlMNTklD3ef5e5XuftVwBdmdn+U6xKJitzsexn82nCuaP0Oj2Zm0KRRTQZZRaSOnQ6cSrCUuLuvA1JDrUhERBJWjb4tmlk6cA5wFvARe79Cnkidy82+l7S8G1ienM7FI3+rsCQSu751dzczBzCzpmEXJCIiiau66zAdSiQknUPkQrHTAKvuQrYisapsWOp8xfM0aao/VovEsOfM7EGghZmNAi4CHgm5JhERSVD/3969x1lV1/sff31EGZghIbwAmgmoebS8gCRmHcPyknkjj5aeLmAXslIpzNL04aHLScvTOdFJJUR/Wr8u1tGGm4WXnJ8d0WQAFW+kkaUjZHYEgx+a4vf8sdfYZpw9zDAze+3L6/l4rMesvfaatT/f2XvmO++91/e7unqL/VHg18AJKaXHASLic2WpSupDt/76Lt5lWJKqRkrp3yLiaOB5CuOYLkkp3ZpzWZKkOtXVGKZTgDXAHRFxdUS8m8K1MKSq0byijU/evI7/GH6JYUmqIimlW1NK56eUPk9httYP5l2TJKk+lQxMKaXmlNLpwD8AdwCfBXaNiKsi4phyFShtq2XzruQnP/sxE8fsxKc/da5hSapwEbFjRFwYEd+NiGOixldNUQAAGyxJREFU4GxgNYUxtJIkld1WZ8lLKW1MKf0opXQi8AZgBfDFfq9M6oXl877LuOVf4rwdb3c2PKl6/IDCKXgrgY9TeLPuNGBySunkPAuTJNWvHv0XmVJ6DpiTLVJFWj7vuxy8/GIeHDSON5/zU8OSVD3GppQOAIiIuRROC39jdhFzSZJy0a3rMEnVojgsOWZJqjovta+klDYDTxmWJEl586131Yzm5U+xufWXbN9oWJKq1EER8Xy2HsDg7HYAKaW0Y36lSZLqlYFJNWF+6++ZcePDHDb6QuZ++CAaG4fkXZKkHkopDci7BkmSOvKUPFW9ZfOu4M3zj+PYN8LcMycaliSVRUScGBFz1q9fn3cpkqR+ZGBSVVs27wrGLb+I/z9oJN/6yBFO8CCpbFJKC1JK04YOHZp3KZKkfmRgUtVqD0sPNRzMXtMXOGZJkiRJfc7ApKr0m0XXG5YkSZLU7wxMqjrNK9r49H8P5JYhJxmWJEmS1K8c8KGqctfin/HFloGMHzOGI6Ze65glSZIk9Ss/YVLVWDbvCt625BN8defbuGbqBMOSJEmS+p2BSVWheIKHE876umFJkiRJZeF/napIi1YvYtbyWazduJbh0cR5a//AQCd4kCRJUpn5CZMqzqLVi5i5ZCZrNq4hkfhL2sDMXXbisfedbViSJElSWRmYVHFmLZ/FC5tf2GLb37aD2Q/NyakiSZIk1SsDkyrO2o1re7RdkiRJ6i8GJlWc4dHU6faRTSPLXIkkSZLqnYFJFWXZvCsLEzy8suX2QQMGMX389HyKkiRJUt0yMKliLJt3JeOWf4nRL7+Jiw77F0Y1jSIIRjWNYubhMzl+7PF5lyhJkqQ647TiqgjNK9q4597VbD9kPPucO48Dml7HKfudmndZkiRJqnMGJuWi+DpLQ7ffibV/eDcT9jyVfaZ8ncaGgXmXJ0mSJAEGJuXga/d8jRtW3fDq7XUvP0vTbjdy2uEHGpYkSZJUURzDpLJatHrRFmGp3eZ4mdkPfDeHiiRJkqTSDEwqq1nLZ5W8z+ssSZIkqdIYmFRWXYUir7MkSZKkSuMYJvW74gkeIIDU6X5eZ0mSJEmVxsCkfrVo9SJmLpnJC5tfyLakQl6KLff7wL4f8DpLkiRJqjgGJvWrWctnFYWlTMB2sR0pJUY2jWT6+OmGJUmSJFUkA5P6VakxSyklHpjyQJmrkSRJknrGSR/Ur4Zv19Tpdid4kCRJUjXwEyb1m+blf+T9azZyzS7wt6JoPmjAICd4kCRVlOYVbVy+eBVPr9vEbsMGc/6x+zJ53O55lyWpAhiY1C+alz/FjJ+t5F17fouLDlvP7IfmsHbjWscsSZIqzsXNK/nhPX98dQ7XtnWbuPCmlQCGJkkGJvW9ZfOvYselN/CO0V/hO2ceTuPA7Tllv1PzLkuSpNdoXtG2RVhqt+mlzVy+eJWBSZJjmNS3ls2/inHLLmTXwYnZHzqYxoFmcklS5bp88aoSVweEp9dtKmstkiqTgUl9pj0sPdxwEGOnL6Cxace8S5IkqUtdhaLdhg0uYyWSKpWBSX1i6YI5hiVJUtUpFYoCOP/YfctbjKSKZGBSrzWvaONf7n6Zewe/w7AkSaoq5x+7L4N3GLDFtgA+eNgbHb8kCXDSB/XSr1puZ8biF5g4ZjwHTp3mmCVJUlVpD0VOKS6plFz+u42I4cANwGjgCeD9KaXnOtlvM7Ayu/nHlNJJ5apRW9c6/yomLbuQ80ecx5SpFxqWJElVafK43Q1IkkrK65S8C4DbU0r7ALdntzuzKaV0cLYYlipI6/yrGJ+NWZryiemGJUmSJNWkvALTycD12fr1wOSc6tA2WP/o7a+GJccsSZIkqZblFZhGpJTWZOtrgREl9hsUEa0RcU9EGKoqwOK77uW9a64wLEmSJKku9Nt5VBFxGzCyk7suKr6RUkoRUeqacXumlNoiYizwq4hYmVL6XSePNQ2YBjBixAhaWlp6V3wvbdiwIfca+sOSp1/m6gde5IzXnceRhxzCs0uX511Sn6vV5w5qu21Q2+2r5bZJklTp+i0wpZSOKnVfRPwpIkallNZExCjgmRLHaMu+ro6IFmAc8JrAlFKaA8wBmDBhQpo0aVLvG9ALLS0t5F1DX2udP5vHH3qWw8YexTvHHs7RRx2Zd0n9ohafu3a13Dao7fbVctuqWUScCJy49957512KJKkf5XVK3nxgSrY+BZjXcYeIeH1ENGTrOwNvBx4uW4V6Vev8qxi37ALOGfL/uGbKITRsH3mXJEm5SyktSClNGzp0aN6lSJL6UV6B6TLg6Ih4DDgqu01ETIiIudk++wGtEXE/cAdwWUrJwFRmhbB0IY80HMS+595IY8MOeZckSZIklU0uc0GnlP4CvLuT7a3Ax7P1JcABZS5NRYrDkhM8SJIkqR558Rx1qnlFG88vvZ1Bgw1LkiRJql8GJr3G/NbfMePGRzls9AxO/fCBNDYOybskSZIkKRd5jWFShWqdfxXjFhzL8Xu8zNwzDzUsSZIkqa4ZmPSq9jFLzzfsxjc+MonGgX4AKUmSpPpmYBLw97D0aMOBjDl3AY1DnCZXkiRJMjCJu3/xI8OSJEmS1AkDU51rXtHGWXfuwOIhkw1LkiRJUgcOUqljdy3+KV9q2Z6DxuzBpKlzHbMkSZIkdeAnTHWqdf5sDlsyjUt3+iXXTJ1gWJIkSZI6YWCqQ63zZzNu2QU82nAgR591uWFJkiRJKsHAVGeKw5JjliRJkqSuGZjqyIKlj7HHsksNS5IkSVI3GZjqRPOKNqbf9FsuG/Etw5IkSZLUTQ5eqQOt82fTdu8SJo7+FP965lsdsyRJkiR1k/8517j2MUuDGw/kzA8faFiSJEmSesBT8mrYayZ4aGzKuyRJkiSpqhiYalTrgu85G54kSZLUS56fVYOaV7Rx+z1PMnDIOPY+p9mwJEmSJG0jP2GqMb+4+z5m/PQ+nt3zBPaecYthSZIkSeoFA1MNaV3wPY785VFM3a2Na6ZOoLFhh7xLkiRJkqqagalGtC74HuNav8jvGvbj82d+wNnwJEmSpD5gYKoB7WHp0YYDGHPuQk/DkyRJkvqIganKtfxqsWFJkiRJ6icGpirWvKKNj976EtcP+4xhSZIkSeoHDnSpUvcuupb/vOsVJo55C6dP/bJjliRJkqR+4H/ZVWjpgjkc0voFvjL0nYybOtWwJEmSJPUTT8mrMksXzGF86xd4tOEAxn3m+4YlSZIkqR/533aFW7R6EbOWz2LtxrUM324I5615wgkeJEmSpDLxE6YKtmj1ImYumcmajWtIJP7yyl+ZuctOrJp8rmFJkiRJKgMDUwWbtXwWL2x+YYttf9sOZj98dU4VSZIkSfXFwFTB1m5c26PtkiRJkvqWgamCDd9uSKfbRzaNLHMlkiRJUn0yMFWopQvmcN6aJxj4ypbbBw0YxPTx03OpSZIkSao3BqYKdO+iaxnf+gX2enkvLpp4CaOaRhEEo5pGMfPwmRw/9vi8S5QkSZLqgtOKV5jmFW38512b+eqOR3Dw2T9g/yFDOWX/0/IuS5IkSapLBqYKcscdtzLjlheZOOYADp56phellSRJknLmKXkVYumCqzmi5TQu3vVurpk6wbAkSZIkVQADUwVYuuBqxreez28b3sLpn/iCYUmSJEmqEAamnBWHpT3PXUTjkKF5lyRJkiQpY2DK0S/uvo8DWi80LEmSJEkVysCUk+YVbXxmfhuX7/w1w5IkSZJUoRwsk4OlC67m1nv+yMQxx3Pe1GMdsyRJkiRVKP9TL7P2MUuDmg5irymXGJYkSZKkCuYpeWVUPMHDXufOp7Fhh7xLkiRJktQFA1OZOBueJEmSVH08H6wMmle0seY3dzKk0bAkSZIkVRMDUz9bsPQxZtz0WyaO/jRTPnwgjY1NeZckSZIkqZtyOSUvIk6LiIci4pWImNDFfu+JiFUR8XhEXFDOGvvC0gVX89aFR3PyHi9wzZlvNSxJkiRJVSavMUwPAqcAd5baISIGAFcAxwH7A2dExP7lKa/32scsPdfwBv71I+92NjxJkiSpCuXyX3xK6RGAiOhqt0OBx1NKq7N9fwKcDDzc7wX20rpHW/jHNd/OJnhYSOOQYXmXJEmSJGkbVPIsebsDTxbdfirbVtHuWvwzTjQsSZIkSTWh3z5hiojbgJGd3HVRSmleHz/WNGAawIgRI2hpaenLw3fbkqdf5gcPDODLTe9l+CEf4k+t9+VSR3/bsGFDbj/j/mbbqlctt6+W2yZJUqXrt8CUUjqql4doA/Youv2GbFtnjzUHmAMwYcKENGnSpF4+dM8tWXwDP1m5HQeP3YOhYz/BUUcdWfYayqWlpYU8fsblYNuqVy23r5bbJklSpavkU/KWAvtExJiIGAicDszPuaZOLV04l0OXnMU3hi/kmqkTaNi+y7FZkiRJkqpEXtOKvy8ingLeBiyKiMXZ9t0i4maAlNLLwNnAYuAR4KcppYfyqLcrSxfOZdzS83ms4c2881OznA1PkiRJqiF5zZL3c+DnnWx/Gnhv0e2bgZvLWFqPFIclJ3iQJEmSak8ln5JX0ea3rmbk0ssMS5IkSVINMzBtg+YVbXz2xkf45oh/MyxJkiRJNcwBNz20dOFc1t1zCxNHf45vnPlWxyxJkiRJNcz/9nugfczSkMb9ef+H3mJYkiRJkmqc//F3098neNifPc9dRGPT6/IuSZLUDyJiMnA8sCNwTUrplpxLkiTlyDFM3XBvx7DkmCVJqkgRcW1EPBMRD3bY/p6IWBURj0fEBV0dI6XUnFL6BHAW8IH+rFeSVPkMTFvRvKKNOXev4ZGGgwxLklT5rgPeU7whIgYAVwDHAfsDZ0TE/hFxQEQs7LDsWvStF2ffJ0mqY56S14VfLlnGjAVrmTjmaPaacgGNDTvkXZIkqQsppTsjYnSHzYcCj6eUVgNExE+Ak1NKlwIndDxGRARwGfCLlNLy/q1YklTpDEwl3LtwLkcu/SLTRs3k3KnHOsGDJFWv3YEni24/BUzsYv9zgKOAoRGxd0ppdscdImIaMA1gxIgRtLS09F21vbRhw4aKqqc/1VNbob7aa1trVzW21xSQWbR6EbOWz2LtxrUM324IM9Y8wY4N+3LuRz9sWJKkOpJS+g7wna3sMweYAzBhwoQ0adKkbXqs5hVtXL54FU+v28RuwwZz/rH7Mnnc7tt0rHYtLS1saz3Vpp7aCvXVXttau6qxvY5hohCWZi6ZyZqNa0gk/vLKX/nyLjvx8OTPOmZJkqpfG7BH0e03ZNty1byijQtvWknbuk0koG3dJi68aSXNK3IvTZJUxMAEzFo+ixc2v7DFtr9tB7MfnptTRZKkPrQU2CcixkTEQOB0YH7ONXH54lVsemnzFts2vbSZyxevyqkiSVJnDEzA2o1re7RdklSZIuLHwN3AvhHxVER8LKX0MnA2sBh4BPhpSumhPOsEeHrdph5tlyTlw8E5wMimkazZuKbT7ZKk6pFSOqPE9puBm8tcTpd2GzaYtk7C0W7DBudQjSSpFD9hAqaPn86gAYO22DZowCCmj5+eU0WSpFp3/rH7MniHAVtsG7zDAM4/dt+cKpIkdcZPmIDjxx4P8OoseSObRjJ9/PRXt0uS1NfaZ8Pr61nyJEl9y8CUOX7s8QYkSVJZTR63uwFJkiqcp+RJkiRJUgkGJkmSJEkqwcAkSdI2iIgTI2LO+vXr8y5FktSPDEySJG2DlNKClNK0oUOH5l2KJKkfGZgkSZIkqQQDkyRJkiSVYGCSJEmSpBIMTJIkSZJUgoFJkiRJkkowMEmSJElSCQYmSZIkSSohUkp519CnIuLPwB9yLmNn4Nmca+hPtdw+21a9arl9ldS2PVNKu+RdRCWpkH6nWCW9XvpbPbUV6qu9trV2VWp7S/ZvNReYKkFEtKaUJuRdR3+p5fbZtupVy+2r5bap79XT66We2gr11V7bWruqsb2ekidJkiRJJRiYJEmSJKkEA1P/mJN3Af2slttn26pXLbevltumvldPr5d6aivUV3tta+2quvY6hkmSJEmSSvATJkmSJEkqwcDUByLitIh4KCJeiYiSs35ExHsiYlVEPB4RF5Szxt6IiOERcWtEPJZ9fX2J/TZHxH3ZMr/cdfbE1p6LiGiIiBuy+38TEaPLX+W26UbbpkbEn4ueq4/nUee2iIhrI+KZiHiwxP0REd/J2v5ARIwvd43bqhttmxQR64uet0vKXaMqU633QcVqsT/qqJb7p87Ucp/VUS33YR3VWp9mYOobDwKnAHeW2iEiBgBXAMcB+wNnRMT+5Smv1y4Abk8p7QPcnt3uzKaU0sHZclL5yuuZbj4XHwOeSyntDfwH8I3yVrltevA6u6HouZpb1iJ75zrgPV3cfxywT7ZMA64qQ0195Tq6bhvAr4uet6+UoSZVh1rvg4rVVH/UUS33T52pgz6ro+uo3T6so+uooT7NwNQHUkqPpJRWbWW3Q4HHU0qrU0p/A34CnNz/1fWJk4Hrs/Xrgck51tIXuvNcFLf5v4B3R0SUscZtVc2vs61KKd0J/E8Xu5wMfD8V3AMMi4hR5amud7rRNqlTddAHFau1/qijWu6fOlMrr8tuqeU+rKNa69MMTOWzO/Bk0e2nsm3VYERKaU22vhYYUWK/QRHRGhH3REQld2LdeS5e3Sel9DKwHtipLNX1TndfZ/+Ufdz/XxGxR3lKK4tq/j3rjrdFxP0R8YuIeHPexaiq1MrvRq31Rx3Vcv/UmXrvszqqld/T7qqaPm37vAuoFhFxGzCyk7suSinNK3c9fa2r9hXfSCmliCg1teKeKaW2iBgL/CoiVqaUftfXtarXFgA/Tim9GBGfpPBO5btyrklbt5zC79iGiHgv0EzhtA3VgVrvg4rZH6kD+6zaVFV9moGpm1JKR/XyEG1A8bsib8i2VYSu2hcRf4qIUSmlNdlHw8+UOEZb9nV1RLQA44BK7KC681y07/NURGwPDAX+Up7yemWrbUspFbdjLvDNMtRVLhX9e9YbKaXni9ZvjogrI2LnlNKzedal8qj1PqhYnfVHHdVy/9SZeu+zOqqa39PeqrY+zVPyymcpsE9EjImIgcDpQLXM3DMfmJKtTwFe825mRLw+Ihqy9Z2BtwMPl63CnunOc1Hc5lOBX6XquGjZVtvW4Xzok4BHylhff5sPfCSbaegwYH3R6TtVLSJGto9TiIhDKfz9rtZ/klR+1dwHFau1/qijWu6fOlPvfVZHNduHdVR1fVpKyaWXC/A+CueZvgj8CVicbd8NuLlov/cCv6XwLtdFedfdg/btRGE2oseA24Dh2fYJwNxs/XBgJXB/9vVjede9lTa95rkAvgKclK0PAn4GPA7cC4zNu+Y+bNulwEPZc3UH8A9519yDtv0YWAO8lP3OfQw4Czgruz8ozLj0u+x1OCHvmvuwbWcXPW/3AIfnXbNLZSy13gd1aGvN9UedtLFm+6dtbG/V9lmdtLVm+7BtaGtV9WmRFS1JkiRJ6sBT8iRJkiSpBAOTJEmSJJVgYJIkSZKkEgxMkiRJklSCgUmSJEmSSjAwKXcRsTki7itaRkfEhIj4Tje+d0n2dXRE/HMvHvuhiLg/Is6LiO2y+16tISIaIuK2bN8PRMQ/Zt9zX0QM7unjlktEzIiIRyNiZda+f4+IHXpxvNER8WC23q3nqItjfWlbv1eSys2+qv/YV6nSOa24chcRG1JKQ3p5jEnA51NKJ2zrY0fErsCPgLtSSv/SYb/DgK+l7Ar0ETEb+O+U0v/t5uMEhd+3V3pSX29ExFnAZOD0lNK67KKAM4ArU9EVtrN9B6SUNnfjmKOBhSmlt/RBfb1+3iWpXOyr+od9lapC3heCcnEBNnSybRKFP3YAM4FrgRZgNXBux++lcNGz9cB9wOeAAcDlFK4i/gDwye48NjCWwpWmo70GYFcKFwhsP/4ngf8Bfg/8MPu+84se68vZttHAKuD7FC7OtidwDHA3sJzChQeHZPs+AXw5276S7MJ8wBDg/2TbHgD+Kdve6XE6tOVJYExXP3fgWxQuGvcO4JKsDQ8Cc/j7GyqHZPvcn/1MH+zkOWrKnqN7gRXAydn2qcBNwC8pXGjym9n2y4DN2c/zh3m/Bl1cXFy2tthX2Vfl/Rp0yW/JvQAXl6I/RvcBP8+2deyElgANwM5ZJ7FDdt+Gjvtnt6cBF2frDUBrZ3+QO3ZC2bZ1wIgONXQ8/nXAqdn6Me1/tCmc5roQOCLrhF4BDsv22xm4E2jKbn8RuCRbfwI4J1v/NH+/Yv03gG8XPe7ruzpO0X47As9t5eeegPcX3R5etP4D4MRs/QHgiGy9VCf0deBD2fowCldtb8o6odXAUApXp/8DsEepn72Li4tLpS72VfZVLvW7bI+Uv00ppYO3ss+ilNKLwIsR8QyFTuKpLvY/BjgwIk7Nbg8F9qHwTltfOyZbVmS3h2SP9UfgDymle7LthwH7A3cVznpgIIV33trdlH1dBpySrR8FnN6+Q0rpuYg4YSvHeY2IOJZChzYM+OeU0hIKnf+NRbsdGRFfABqB4cBDEfFrYFhK6c5snx8Ax5X4GZwUEZ/Pbg8C3pit355SWp/V8TCFdy+f7KpeSapA9lUF9lWqOwYmVYsXi9Y3s/XXblB4F2xxTx4kIsZmx38G2K+73wZcmlL6XodjjQY2dtjv1pTSGSWO097GrbVva8chpfR8RGyIiDEppd9nP4fFEbGQQqcF8ELKzgWPiEHAlcCElNKTETGTQkfSXUHhFIxVW2yMmEjPnztJqlb2Vd0/jn2Vqoaz5KlW/BV4XdHtxcCn2mfZiYg3RURTVweIiF2A2cB3U0o9mQ1lMfDRiGgfkLt7Nii3o3uAt0fE3tl+TRHxpq0c+1bgM0U1vr4Hx7kUuCoihmX7BaU7lvbtz2btOBUgpbQOWBcR78ju/2CJ718MnJM9BhExbivtAnipN7MgSVIVsq96LfsqVTzTs2rFA8DmiLifwjnbsyicl708+8P4Zwqz8HQ0OCLuA3YAXqbwMf6/9+SBU0q3RMR+wN3Z3+ANwIcovENVvN+fI2Iq8OOIaMg2X0zhHOpSvgZckU2PupnCIN2bunmcqyicm/2biHgxq+su/n46RnFt6yLiagqDaNdSGFDb7kzg2ohIwC0l6vwq8G3ggShMdft7YGuzQM3J9l+eUirVuUlSLbGvsq9SFXJacUmSJEkqwVPyJEmSJKkEA5MkSZIklWBgkiRJkqQSDEySJEmSVIKBSZIkSZJKMDBJkiRJUgkGJkmSJEkqwcAkSZIkSSX8L9eGc3ioWXlOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "min_g = np.min(g_discrete)\n",
    "max_g = np.max(g_discrete)\n",
    "\n",
    "fig = plt.figure(figsize=(12,6))\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "plt.plot([min_g, max_g],[min_g, max_g],label='y=x comparison')\n",
    "plt.plot([min_g, max_g],[m*min_g+b, m*max_g+b],'--',label='Best fit')\n",
    "plt.plot(g_discrete,dJ_du[idx],'o',label='Adjoint comparison')\n",
    "plt.xlabel('Finite Difference Gradient')\n",
    "plt.ylabel('Adjoint Gradient')\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.axis(\"square\")\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "rel_err = np.abs(np.squeeze(g_discrete) - np.squeeze(dJ_du[idx])) / np.abs(np.squeeze(g_discrete)) * 100\n",
    "plt.semilogy(g_discrete,rel_err,'o')\n",
    "plt.grid(True)\n",
    "plt.xlabel('Finite Difference Gradient')\n",
    "plt.ylabel('Relative Error (%)')\n",
    "\n",
    "fig.tight_layout(rect=[0, 0.03, 1, 0.95])\n",
    "fig.suptitle('Resolution: {} Seed: {} Nx: {} Ny: {}'.format(resolution,seed,Nx,Ny))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We notice strong agreement between the adjoint gradients and the finite difference gradients. Let's bump up the resolution to see if the results are consistent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting forward run...\n",
      "Starting adjoint run...\n",
      "Calculating gradient...\n"
     ]
    }
   ],
   "source": [
    "resolution = 20\n",
    "opt.sim.resolution = resolution\n",
    "f0, dJ_du = opt()\n",
    "g_discrete, idx = opt.calculate_fd_gradient(num_gradients=choose,db=db)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0kAAAGeCAYAAABB3Ur+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdeXhU5d3/8feXAIKAoKhsLqAoyiIJBBQskAASqoK4i1GKVqyPFrH2oYpaRatPbWmttCooiogi4oKoxV9RH02BagVkEdnE8mDNJmHJJIRBknD//pgz6SQkIfuZJJ/XdeVK5j4z53zPYXTyyb0cc84hIiIiIiIiIU38LkBERERERCSaKCSJiIiIiIhEUEgSERERERGJoJAkIiIiIiISQSFJREREREQkgkKSiIiIiIhIBIUkEZEIZjbPzB6txus3mVlCDZbUYJhZgpml+l2HiIjI0SgkiUjUMrOdZhY0s/1mlukFmNZ+1xVWWqByzvVyzqXUwbH/YGbbzSzXzLaa2YQS22PN7AszO+B9jy1nX73M7AMz22tm2d7zL67tcyiPmZ1sZgvNLN3MAmb2DzM7v4znzjUzZ2bdI9pOMLO3zSzPzL41s+vLOdZ07/XXRLQ19dq61sC5JJrZJ9557Cxle1dv+wHv33JkOfua59U1MKKtu5nVyE0PzewaM/vUqyWllO0Vfl+JiNRnCkkiEu3GOOdaA7FAHDDN53qiRR4wBmgL/ASYaWaDAcysOfAO8ApwPPAS8I7XXpr3gA+BjsDJwJ1ATq1Wf3StgdVAf+AEQuewtGRINrMfAWeW8vqngUNAByAZmGVmvco53l7gYTOLqYHaS8oD5gJTy9i+EFgHtAfuB940s5PK2d9eoMq9nUexF3gSeLzkhiq8r0RE6i2FJBGpF5xzmcAyQmEJADM7xutR+beZfW9ms82spbftRDP7q9czstfMVphZE2/buWaW4m3bZGZjSzummU00s5Ul2pz3l/tbCf3y/Suvp+s9b/vOcE+AV9+TXm9IuvfzMd62BDNLNbNfmtkuM8sws5sqcT0ecs5tdc4dds59DqwABnmbE4CmwJPOuR+cc38GDBheyjmeCHQD5jjnDnlf/3DOrYx4zqVmtt67Xp+a2XkR2zqb2VtmlmVm/2dmd0Zsa+n1fOwzs83AgEqc3w7n3BPOuQznXKFz7jmgOdAjYv9Ngb8Ak0ucUyvgSuDXzrn93rm8C9xYziH/RihU3VDKNWrunf9k73GM17P1YAXPZZVz7mVgRyn7PhvoBzzknAs6594CNnr1l+Ul4DwzG1bK/q42sy9KtN1tZu9UsNaPnHOvA+mlbE6ggu8rEZH6TiFJROoFMzsF+DHwTUTz48DZhIJTd6ALEP7F9ZdAKnASod6E+wBnZs0I9Zx8QKjXZDKwwMx6UAneL+0LgN8751o758aU8rT7gQu8+voCA4EHIrZ3JNQT1AX4KfC0mR3vne/1ZvZlRWrxguEAYJPX1Av40jkXOQTrS6+9pD2ErukrZjbOzDqU2HccoV6QnxHq6XgWeNcLgE0IXcsN3jmMAO4ysyTv5Q8R6uU5E0gi1OMVue9nzOyZCp5jLKGQFPnv/wtguXOu5HU6Gyhwzn0d0bahjPMPc8CvgYe898h/NjgXDk+PmNm5wL1ADPCYV9uPzCy7IudRil7ADudcbiVqPQD8T/j4JbwLdPPqDLsRmO/VWuH3VRm1VvR9JSJSrykkiUi0W2JmucB3wC5Cv3hjZgbcCvzCObfX+yXzf4DrvNflA52A051z+c65Fd4vdxcQGsr1uNdr8jHwV2B8LdSeDDzinNvlnMsCHqZ4b0a+tz3fOfc+sB+vp8Q596pz7rwj9li62YR+sV7mPW4NBEo8JwC0KflC75okAjuBPwIZZrbczM7ynnIr8Kxz7nOvR+cl4AdC13EAcJJz7hHvWu4A5vCff4NrgMe8f5/vgD+XOPbtzrnbj3ZyZnYc8DLwsHMu4LWdSii4ldab05ojhwuWev4l6nkXyAJuKWXbV4SGuC0B/hu40TlX6G1b6Zxrd7TzKEOF/61KeBY4zcx+XKLOH4BFeD1i3hDDroTe45V9X9VUrSIi9Y5CkohEu3HOuTaEhvqcA5zotZ8EHAt84Q0DyyY0ZCo8l2MGoV6HD8xsh5nd67V3Br5zzh2OOMa3hHpCalpnb9+Rx+kc8XiPc64g4vEBQr+IVpiZzQB6A9dE/IV/P3BciaceB+RSCudcqnPu5865M4HTCc2hme9tPh34Zfgae9f5VO88Tgc6l9h2H6GeO7znfBdxqMhrUdHza0mot+qfzrnfRmx6klDALPlLO1Ty/Et4gFAPYItStr1E6Jzfd85tr8C+KqJKtXph6DfeV0kvAdd7f0i4EXjde74vtYqI1EcKSSJSLzjn/g7MA/7gNe0GgkAv51w776utt8gDzrlc59wvnXNnAGOBu81sBKG5FqeG5yd5TgPSSjlsHqEgBoCZdSxZ1lHKTif0S3XkcUqb61ElZvYwoSGIo5xzkT0nmwjNWbGItvP4z3C8Mnk9Pk8TCl4QCjmPRVzjds65Y51zC71t/1diWxvnXHhlvAxCgSrstEqe3zGEem5SCfUaRRoBzLDQqoeZXttnFlrF7mugaURvGISGO1bk/D8kFK5L6+F6hlCPTJKFFoyoCZuAM8wssjemQrUCLwLtgCsiG51z/yQ0v2oIcD2hXriaUOX3lYhIfaOQJCL1yZPARWbW1+sJmgP8ycxOBjCzLuH5MN5iA929X+gCQCFwGPicUI/Nr8ysmYXuaTQGeK2U420Aello2eMWwPQS278Hziin3oXAA2Z2krdAwoOEVgarNjObRugX4JHOuT0lNqcQOt87vblDP/faPy5lP8eb2cPetWri1Xkz8E/vKXOA28zsfAtpZWaXeL/UrwJyzeweb5GGGDPrbWbhBRpeB6Z5xziFEgssHOX8mgFvEgrCPynR8weheUd9Cc33Ci/mMQZ42zmXBywmNIeolZldCFxGxcPC/cCvStRzI6GV9iYSWv3vJavgcvTedW0BNAs9tBbmrQjnzZtaT2guVAszu5xQ8HjraPv1eiEfAu4pZfN84CkgP3IRjgrUGuPV2hRo4tUUnqOVQgXfVyIi9Z1CkojUG968nvn8Zx7KPYT+6v9PM8sBPuI/q5+d5T3eD3wGPOOc+8SbhD+GUA/MbkK9AxOcc1tLOd7XwCPefrYDJX/ZfAHo6Q01W1JKyY8CawhNbt8IrKWCSzebWbKZlfcX+v8h1DPzjYVW19tvZvd5dR8CxgETgGxCoWec117SIUJzVj4iNI/nK0JzjiZ6+1oDTCL0C/c+Qtc7vK0QuJRQSPk/QtfzeUKLUUBoDta33rYPKBFSLLQa4ewyzm+wt+9RQHbEOQ7xjr3LOZcZ/vJes9s5F/R+vh1oSWge20Lgv5xzFerxcM79g1AADNd5GqGAPsFbLe9VQv+uf/K2DzGz/eXsciihsPc+oX+zoHc9wq4D4gld38eBq7z3ekUsJNRjV9LLhHoDi4XyCryvbvTqm0WoJypIKChX9n0lIlKvWfFFakRERKS+8+Zy7QL61eD8KRGRRkM9SSIiIg3PfwGrFZBERKqmqd8FiIiISM0xs52EbvI6zudSRETqLQ23ExERERERiaDhdiIiIiIiIhEUkkRERERERCIoJImIiIiIiERQSBIREREREYmgkCQiIiIiIhJBIUlERERERCSCQpKIiIiIiEgEhSQREREREZEICkkiIiIiIiIRFJJEREREREQiNPW7gMo48cQTXdeuXf0uwzd5eXm0atXK7zKilq5P+XR9yqfrA1988cVu59xJftcRjRr7509Z9N9N1ei6VZ2uXdXoupWuvM+9ehWSunbtypo1a/wuwzcpKSkkJCT4XUbU0vUpn65P+XR9wMy+9buGaNXYP3/Kov9uqkbXrep07apG16105X3uabidiIiIiIhIBIUkERERERGRCApJIiIiZTCzMWb2XCAQ8LsUERGpQ/VqTlJp8vPzSU1N5eDBg36XUuvatm3Lli1b/C4jarVu3Zr8/HyaNWvmdyki0kA4594D3ouPj5/kdy0iIlJ36n1ISk1NpU2bNnTt2hUz87ucWpWbm0ubNm38LiMqOedITU0lNTWVbt26+V2OiIiIiNRj9X643cGDB2nfvn2DD0hSPjOjbdu2jaJHUURERERqV70PSYACkgB6H4iIiIhIzWgQIUnqVnp6OldddZXfZYiIiIiI1AqFJKmUgoICOnfuzJtvvul3KSIiIiIitUIhqZoefPBBnnzyyaLH999/PzNnzjzq6wKBAD169GDbtm0AjB8/njlz5hzxvNWrVzN48GD69u1LQkICubm5HDx4kJtuuok+ffoQFxfHJ598AsC8efMYN24cF110EV27duWpp57iiSeeIC4ujgsuuIC9e/cCkJCQwJQpU4iNjaV3796sWrUKgFWrVjFo0CDi4uIYPHhwUW3z5s1j7NixDB8+nBEjRrBz50569+4NwKZNmxg4cCCxsbGcd955bN++HYAnnniC3r1707t376Lrs3PnTs4991wmTZpEr169GDVqFMFgsErXXURERESkttT71e0iPfzeJjan59ToPnt2Po6HxvQqc/vNN9/MFVdcwV133cXhw4d57bXXWLVqFbm5uQwZMqTU17z66qv07NmTp556iokTJzJlyhT27dvHpEnFV5g9dOgQ1157LYsWLWLAgAGkpaXRsmVLZs6ciZmxceNGtm7dyqhRo/j6668B+Oqrr1i3bh0HDx6ke/fu/O53v2PdunX84he/YP78+dx1110AHDhwgPXr17N8+XJuvvlmvvrqK8455xxWrFhB06ZN+eijj7jvvvt46623AFi7di1ffvklJ5xwAjt37iyqcfbs2UyZMoXk5GQOHTpEYWEhX3zxBS+++CKff/45zjnOP/98hg0bxvHHH8/27dtZuHAhc+bM4ZprruGtt97ihhtuqM4/kYiIiIhIjWpQIckPXbt2pX379qxbt47vv/+euLg42rdvD8D69evLfe1FF13EG2+8wR133MGGDRuO2L5t2zY6derEgAEDADjuuONo2rQpK1euZPLkyQCcc845nH766UUhKTExkTZt2tCmTRvatm3LmDFjAOjTpw9ffvll0b7Hjx8PwNChQ8nJySE7O5vc3Fx+8pOfsH37dsyM/Pz8YrWecMIJR9Q4aNAgHnvsMVJTU7niiis466yzWLlyJZdffjmtWrUC4IorrmDFihWMHTuWbt26ERsbC0D//v2LBS4RqZqlO5Yyc+1MMvMy6diqI1P6TeGSMy7xuyzxLFmXxoxl20jPDtK5XUumJvVgXFwXv8sSEZFyNKiQVF6PT2265ZZbmDdvHpmZmdx8880AFepJOnz4MFu2bOHYY49l3759nHLKKdWu5Zhjjin6uUmTJkWPmzRpQkFBQdG2kivBmRm//vWvSUxM5O2332bnzp0kJCQUbQ8HnpKuv/56zj//fJYuXcrFF1/Ms88+W+H6YmJiNNxOpJqW7ljK9E+nc7AwtPx9Rl4G0z+dDqCgFAWWrEtj2uKNBPMLAUjLDjJt8UYABSURkSimOUk14PLLL+dvf/sbq1evJikpCYA2bdqwfv36Ur969uwJwJ/+9CfOPfdcXn31VW666aZiPTcAPXr0ICMjg9WrVwOh4FVQUMCQIUNYsGABAF9//TX//ve/6dGjR6VqXrRoEQArV66kbdu2tG3blkAgQJcuoQ/tefPmVWg/O3bs4IwzzuDOO+/ksssu48svv2TIkCEsWbKEAwcOkJeXx9tvv11mYBSR6pm5dmZRQAo7WHiQmWuPPjdSjs7MxpjZc4FAoEqvn7FsW1FACgvmFzJj2baaKE9ERGpJg+pJ8kvz5s1JTEykXbt2xMTEVOg127Zt4/nnn2fVqlW0adOGoUOH8uijj/Lwww8X2++iRYuYPHkywWCQ5s2b88knn3D77bfzX//1X/Tp04emTZsyb968Yj00FdGiRQvi4uLIz89n7ty5APzqV7/iJz/5CY8++iiXXFKxv0C//vrrvPzyyzRr1oyOHTty3333ccIJJzBx4kQGDhwIhHra4uLiNLROpBZk5mVWql0qxzn3HvBefHz8pKM+uRTp2aX3lpfVLiIi0cGcc37XUGHx8fFuzZo1xdq2bNnCueee61NFIYcPH6Zfv3688cYbnHXWWbV2nNzcXNq0aVPt/SQkJPCHP/yB+Pj4GqgqeuTm5pKamur7+yFapaSkFBtCKcXVx+uTk/Udl7/7Y3Y1PfJGyp1adeKDqz6o1P7M7AvnXMP6H0MNKe3zpyIufPxj0koJRF3ateQf9w6vidJ8VR//u4kGum5Vp2tXNbpupSvvc0/D7app8+bNdO/enREjRtRqQBIRiZST9R25s5O4fU8Oza1ZsW0tYlowpd8UnyqTSFOTetCyWfERBi2bxTA1qXJDpEVEpG5puF019ezZkx07dvhdRqWkpKT4XYKIVEMgmM8Hcx/hxwV7OG3oXB45Fa1uF6XCizNodTsRkfpFIUlEpB4JBPOZ8MLnbM25jC5jJjL4ggsBrWQXzcbFdVEoEhGpZzTcTkSknsjJ+o5tT1zM3oz/4+nkAUUBSURERGqWepJEROqBnKzvyJk9ml4FWTyRdBIDenbwuyQREZEGSz1JIiJRLhyQji/IYvPwFxkwdLTfJYmIiDRoCkk1ICYmhtjYWPr27Uu/fv349NNPq7SfJ598kgMHDpS6bcWKFQwcOJDY2FjS0tK46qqrAFi/fj3vv/9+lWsXkeh2REAaprlHIiIitU0hqQa0bNmS9evXs2HDBn77298ybdq0Ku2nvJC0YMEC7r77btavX0+XLl148803AYUkkYYsEMzn9oVfkpnfSgFJRESkDikk1bCcnByOP/74osczZsxgwIABnHfeeTz00EMA5OXlcckll9C3b1969+7NokWL+POf/0x6ejqJiYkkJiYW2+fzzz/P66+/zmOPPUZycjI7d+6kd+/eHDp0iAcffJBFixYRGxvLokWL6vRcRaT25OxO5+bnV/L590b2te8qIImIiNShhrdww4ul/CLRaxwMnASHDsCCq4/cHns9xCVD3h54fULxbTctPeohg8EgsbGxHDx4kIyMDD7++GMAPvjgA7Zv386qVatwzjF27FiWL19OVlYWnTt3ZunS0L4DgQBt27bliSee4JNPPuHEE08stv9bbrmFlStXMmLECG688UZ27twJQPPmzXnkkUdYs2YNTz311NGvjYjUC+Ehdjfln0qL5JcYqUUaRERE6pR6kmpAeLjd1q1b+dvf/saECRNwzvHBBx/wwQcfEBcXR79+/di6dSvbt2+nT58+fPjhh9xzzz2sWLGCtm3b+n0KIhIlIucgdRgxWQFJRETEBw2vJ6m8np/mx5a/vVX7CvUclWfQoEHs3r2brKwsnHNMmzaNn/3sZ0c8b+3atbz//vs88MADjBgxggcffLBaxxWR+k+LNEQfMxsDjOnevbvfpYiISB1ST1IN27p1K4WFhbRv356kpCTmzp3L/v37AUhLS2PXrl2kp6dz7LHHcsMNNzB16lTWrl0LQJs2bcjNza3U8aryGhGJPoEDh0h/9ioFpCjjnHvPOXerevxFRBqXhteT5IPwnCQA5xwvvfQSMTExjBo1ii1btjBo0CAAWrduzSuvvMI333zD1KlTadKkCc2aNWPWrFkA3HrrrYwePZrOnTvzySefVOjYiYmJPP7448TGxjJt2jSuvfba2jlJEak1gWA+E+auwoI3cv+obgpIIiIiPlNIqgGFhYVlbpsyZQpTpkwp1nbmmWeSlJR0xHMnT57M5MmTS93PvHnzinqMunbtyldffQXACSecwOrVq6tauoj4LCfrO16dN4vN2T9iVvLVDNAcJBEREd8pJImI+CQn6zsCs0czoSCL3pdfwxAFJBERkaigOUkiIj4IB6QTCrLYMnwuQ+Jj/S5JREREPApJIiJ1rGRAih92qd8liYiISASFJBGROhQI5jNr/iscX7BHAUlERCRKaU6SiEgdCRz4gQlzV7N5b18uuPIThvU71++SREREpBTqSRIRqQM5Wans/uMgjs9cyazk/gpIIiIiUUwhqYYsWbIEM2Pr1q1lPmfixIm8+eabANxyyy1s3ry53H1efPHFZGdnl/ucefPmkZ6eXvmCfTR48GC/SxCpUzlZqQRmJ9GxIJWfX9SbkVrFTkREJKo1upC0dMdSRr05ivNeOo9Rb45i6Y6lNbLfhQsX8qMf/YiFCxdW6PnPP/88PXv2LPc577//Pu3atSv3OfUpJBUUFADw6aef+lyJSN0JByQt0iAiIlJ/NKqQtHTHUqZ/Op2MvAwcjoy8DKZ/Or3aQWn//v2sXLmSF154gddee62o3TnHz3/+c3r06MHIkSPZtWtX0baEhATWrFkDhAJWnz596N27N/fcc0/Rc7p27cru3bvZuXMn5557LpMnT6ZXr16MGjWKYDDIm2++yZo1a0hOTiY2NpZgMFisrm+++YaRI0fSt29f+vXrx7/+9S+cc0ydOpXevXvTp08fFi1aBEBKSgrDhg3jsssu44wzzuDee+9lwYIFDBw4kD59+vCvf/0LCPWG3XbbbcTHx3P22Wfz17/+FYCdO3cyZMgQ+vXrR79+/YqCUEpKCkOGDGHs2LFFobB169YAZGRkMHToUGJjY+nduzcrVqwo93q0bt2a+++/n759+3LBBRfw/fffV+vfTaS25ezdpYAkIiJSDzWqkDRz7UwOFh4s1naw8CAz186s1n7feecdRo8ezdlnn0379u354osvAHj77bfZtm0bmzdvZv78+aX2oKSnp3PPPffw8ccfs379elavXs2SJUuOeN727duZNGkSmzZtol27drz11ltcddVVxMfHs2DBAtavX0/Lli2LvSY5OZk77riDDRs28Omnn9KpUycWL17M+vXr2bBhAx999BFTp04lIyMDgA0bNjB79my2bNnCyy+/zNdff82qVau45ZZb+Mtf/lK03507d7Jq1SqWLl3KbbfdxsGDBzn55JP58MMPWbt2LYsWLeLOO+8sev7atWuZOXMmX3/9dbH6Xn31VZKSkorqiY2NLfd65OXlccEFF7BhwwaGDh3KnDlzqvgvJlL7AsF8JizYyieHeikgiYiI1DONKiRl5mVWqr2iFi5cyHXXXQfAddddVzTkbvny5YwfP56YmBg6d+7M8OHDj3jt6tWrSUhI4KSTTqJp06YkJyezfPnyI57XrVs3zjvvPAD69+/Pzp07y60pNzeXtLQ0Lr/8cgBatGjBsccey8qVK4tq6tChA8OGDWP16tUADBgwgE6dOnHMMcdw5plnMmrUKAD69OlT7HjXXHMNTZo04ayzzuKMM85g69at5OfnM2nSJPr06cPVV19dbL7VwIED6dat2xE1DhgwgBdffJHp06ezceNG2rRpU+71aN68OZdeemmFr4GIX3KyUvnlc++xKXM/ncf/RQFJRESknmlUS4B3bNWRjLyMUturau/evXz88cds3LgRM6OwsBAzY8aMGdUp9QjHHHNM0c8xMTFHDK2r6WM0adKk6HGTJk2K5hMBmFmx15kZf/rTn+jQoQMbNmzg8OHDtGjRomh7q1atSj3e0KFDWb58OUuXLmXixIncfffdtG3btsz6mjVrVnTsmJiYYjWJRIvwHKSpBcZ14z/QIg0iIiL1UKPqSZrSbwotYloUa2sR04Ip/aZUeZ9vvvkmN954I99++y07d+7ku+++o1u3bqxYsYKhQ4eyaNEiCgsLycjI4JNPPjni9QMHDuTvf/87u3fvprCwkIULFzJs2LAKH79Nmzbk5uaW2n7KKacUDVX74YcfOHDgAEOGDCmqKSsri+XLlzNw4MBKnfMbb7zB4cOH+de//sWOHTvo0aMHgUCATp060aRJE15++WUKCwuPup9vv/2WDh06MGnSJG655RbWrl1b7esh4qfIRRpyhz/OyN6d/S5JREREqqBR9SRdcsYlQGhuUmZeJh1bdWRKvylF7VWxcOHCYosLAFx55ZUsXLiQZ555ho8//piePXty2mmnMWjQoGLPMzM6derE448/TmJiIs45LrnkEi677LIKHz+8kELLli357LPPis1Levnll/nZz37Ggw8+SLNmzXjjjTe4/PLL+eyzz+jbty9mxu9//3s6duxY7tLlJZ122mkMHDiQnJwcZs+eTYsWLbj99tu58sormT9/PqNHjy6z9yhSSkoKM2bMoFmzZrRu3Zr58+dX+3qI+EWr2DVMZjYGGNO9e3e/SxERkTpkzjm/a6iw+Ph4F14RLmzLli2ce279uyljnz59ePfdd0udq1OW3Nxc2rRpU4tVHd3EiRO59NJLueqqq3ytozS5ubmkpqbWy/dDXUhJSSEhIcHvMqJWda5PIJjP+ievIv7gZ/U6IJnZF865eL/riEalff6I/r9SVbpuVadrVzW6bqUr73OvUfUkRYuLLrqIPn36VCogiUh0CgTzmfDC56TuT2bW6DsYOGSU3yWJiIhINSkk+eDDDz/0u4Qqmzdvnt8liESNnKxU/vHCr/gm9xpmJg9joBZpEBERaRAa1cINIiI1JTwHKSH4EXMvPk6r2ImIiDQgDSIk1ad5VVJ79D6QuhIOSO0Lstg8fC7nX3jkPdBERESk/qr3IalFixbs2bNHvyA3cs45AoFAsfszidSGkgGpvi7SICIiImWr93OSTjnlFFJTU8nKyvK7lFp38OBBhYBy5OXl0bdvX7/LkAYsEMzngQUp3FNwUAFJRESkAav3IalZs2aNZpW4lJQU4uLi/C4jaqWkpNCsWTO/y5AGKhDYx4T5X7E56wTGXZfCiD6n+l2SiIiI1JJ6P9xORKS25WSlkvPnISTseolZyf0VkERERBq4et+TJCJSm0JzkEbTvmAXQ0aOI16r2ImIiDR46kkSESlDZEDSHCQREZHGQyFJRKQUgf0H2Dv7YgUkERGRRkjD7URESggE85kwby1n/jCG6y8apIAkIiLSyCgkiYhEyMlK5Y8vv83mPWcyOflOzUESERFphHwLSWZ2KjAf6AA44Dnn3Ey/6hERCc9BurtgL4nXrCBRAUlERKRR8qh8z5cAACAASURBVLMnqQD4pXNurZm1Ab4wsw+dc5t9rElEGqmCvL3/WaQh8QUS+57pd0kiIiLiE98WbnDOZTjn1no/5wJbgC5+1SMijVdOVird1zxQFJDiE8b4XZKIiIj4KCpWtzOzrkAc8Lm/lYhIYxMI5vPOi7/j5MO7FZBEREQEiIKFG8ysNfAWcJdzLqeU7bcCtwJ06NCBlJSUui0wiuzfv79Rn//R6PqUT9fnSHn5jj+sOci/c5I4cEYfetBG10hERET8DUlm1oxQQFrgnFtc2nOcc88BzwHEx8e7hISEuiswyqSkpNCYz/9odH3Kp+tTXM7uNL559gaaHJjAszdeTNNdW3R9RBqxJevSmLFsG+nZQTq3a8nUpB6Mi9MsAJHGyrfhdmZmwAvAFufcE37VISKNT87uNLJnjeacQ5v43aiTGKlV7KQMZjbGzJ4LBAJ+lyK1aMm6NKYt3khadhAHpGUHmbZ4I0vWpfldmoj4xM85SRcCNwLDzWy993Wxj/WISCMQDkgnFnwfmoOkG8VKOZxz7znnbm3btq3fpUgtmrFsG8H8wmJtwfxCZizb5lNFIuI334bbOedWAubX8UWk8cnZnV48IGmRBhEB0rODlWoXkYYvKla3ExGpbYFgPj97dSPf5R+ngCQixXRu17JS7SLS8CkkiUiDl7M7nVue/ztrvi8keO1bCkgiUszUpB60bBZTrK1lsximJvXwqSIR8ZvvS4CLiNSm8Byk2/JPxCW/pkUaROQI4VXstLqdiIQpJIlIgxW5SEPW8Mfor4AkImUYF9dFoUhEimi4nYg0SJEBaUvi8/RPGOt3SSIiIlJPKCSJSIMTCOaz89nrFJBERESkSjTcTkQalEAwnwkvfM6hAzfym4s6aZEGEalVS9alaS6TSAOkkCQiDUbO7jTefPEJNmcPZ1byZcRrDpKI1KIl69KYtnhj0Y1o07KDTFu8EUBBSaSe03A7EWkQwnOQrt8/n3ljT9IqdiJS62Ys21YUkMKC+YXMWLbNp4pEpKaoJ0lE6r1wQDqpIJPNiS9w4fkD/S5JRCI01CFp6dnBSrWLSP2hniQRqddKBiQt0iASXcJD0tKygzj+MyRtybo0v0urts7tWlaqXUTqD4UkEam3AsF8npy/iOMLshSQRKJUQx6SNjWpBy2bxRRra9kshqlJPXyqSERqiobbiUi9FMg7yIQX17B5zzkMvervJMTplxKRaNSQh6SFhww2xKGEIo2dQpKI1Ds5u9PZM+sSuvwwjsnJt5GgRRpEolbndi1JKyUQNZQhaePiuigUiTRAGm4nIvVKzu50smcl0akglZ+OjNMqdiJRTkPSRKQ+Uk+SiNQb4YCkRRpE6g8NSROR+kghSUTqhUD2XgIKSCL1koakiUh9o5AkIlEvEMxnwsubuPhQLPEjrlZAEhERkVqlkCQiUS1ndzr3LljB5l1tmJz8B/prDpKIiIjUMoUkEYla4TlI9xb8wFXjP2aEApKIiIjUAa1uJyJRKXKRht2JMxjR+xS/SxIREZFGQiFJRKJO8VXsnqd/wmV+lyQiIiKNiEKSiESVQDCff865SwFJREREfKM5SSISNQLBfCa88Dk791/HnNG3MHDIKL9LEhERkUZIIUlEokLO7nQ+n3MXO/dfxx+TL2SgFmkQERERnygkiYhvlu5Yysy1M8nMy+TEAsfkmOxQD5ICkoiIiPhIc5JExBdLdyxl+qfTycjLwOHIagqPdjyJrC75fpcmIiIijZxCkoj4YubamRwsPFis7ZDLZ+bamT5VJCIS/ZasS+PCxz+m271LufDxj1myLs3vkkQaJA23ExFfZOZlVqpdxA9mNgYY0717d79LEWHJujSmLd5IML8QgLTsINMWbwRgXFwXP0sTaXDUkyQidS6QvZeYw8eXuq1jq451XI1I2Zxz7znnbm3btq3fpYgwY9m2ooAUFswvZMaybT5VJNJwKSSJSJ3K2Z3Ovr8kMHxXK5o1OabYthYxLZjSb4pPlYmIRLf07GCl2kWk6hSSRKTO5OxOZ++s0XQoyOD683/Gby58mE6tOmEYnVp1Yvrg6VxyxiV+lykiEpU6t2tZqXYRqTrNSRKROhEZkDYnPk//hMsAFIpERCpoalKPYnOSAFo2i2FqUg8fqxJpmBSSRKTWBfIOsnvWpXQuEZBERKTiwoszzFi2jfTsIJ3btWRqUg8t2iBSCxSSRKRWBYL5THhxDZ1+GMekkX0VkEREqmFcXBeFIpE6oJAkIrUmZ3c6M+e/xuY95zA5+Wf079nB75JEREREjkohSURqRXgO0t0Fuxh2VQrDFJBERESkntDqdiJS4yIXadiWMJthcef4XZKIiIhIhSkkiUiNKraKXcIc+ieO87skERERkUpRSBKRGhMI5rPoxT8pIImIiEi9pjlJIlIjAsF8JrzwOZuzE+g59mouPH+g3yWJiIiIVIl6kkSk2nL2ZLDziREcytjErOR4BSQRERGp19STJCLVkrMng73PJHF2QQaPXtRRy3yLiIhIvaeeJBGpsnBACs1Bek43ihUREZEGQSFJRKokZ09m8YCUeLnfJYmIiIjUCIUkEam0QDCfn776FTvy2ysgiYiISIOjOUkiUik5ezK4bcGXrP++kMPJr2kOkoiIiDQ4CkkiUmHhOUiT81tz4PrFjFRAEhERkQZIIUlEKiRykYY9ic8xuFdHv0sSERERqRWakyQiR3XEKnaagyQiIiINmEKSiJQrEMzn69k3KCCJiIhIo6HhdiJSpkAwnwkvfM7+Azfw+4va0z9hrN8liYiIiNQ6hSQRKVXOngzeeeG3bA4kMSv5Yq1iJyIiIo2GQpKIHCE0B2k01xSkcdal1zJIAUlEREQaEYUkESkmHJA6FqSxKWEOgwZd6HdJIiIiIgAsWZfGjGXbSM8O0rldS6Ym9WBcXJcaP45CkogUKRmQtEiDiIiIRIsl69KYtngjwfxCANKyg0xbvBGgxoOSVrcTESC0SMPvX1rM8QVZCkgiIiISdWYs21YUkMKC+YXMWLatxo+lniQRIZAXZMKLX7B5TzdGXL2cxNjufpckIiIiUkx6drBS7dWhniSRRi5nTyZ7nhhM98z3mZXcXwFJREREolLndi0r1V4dCkkijVjOnkz2PpNE54LvuH7EAEZqFTsRERGJUlOTetCyWUyxtpbNYpia1KPGj6XhdiKNVDggaZEGqWtmNgi4ARgCdAKCwFfAUuAV51zAx/JERCRKhRdn0Op2IlIrAoFs9ikgiQ/M7P8B6cA7wGPALqAFcDaQCLxjZk845971r0oREYlW4+K61EooKumoIcnMXnbO3Xi0NhGpHwLBfCa8/BUJhwYwdPh0BSSpazc653aXaNsPrPW+/mhmJ9Z9WSIiIv9RkTlJvSIfmFkM0L8mDm5mc81sl5l9VRP7E5Hy5ezJ5IFnF7E5I4c+4x9TQJI6V0pAwsxGmNkYM2tW1nNERETqUpkhycymmVkucJ6Z5XhfuYSGRrxTQ8efB4yuoX2JSDkOHQiw95kk7tv3a2Zf10eLNEhUMLM/AhcCfam5zxYREZFqKTMkOed+65xrA8xwzh3nfbVxzrV3zk2riYM755YDe2tiXyJStpw9mZy5+gE6FqSRkfBHRvQ51e+SpJEysz+aWbuIptOA3xCan3SaP1WJiIgUd9Q5Sc65aWbWBTg98vlewBGRKBdexa7L4Qw2JzxHv8Qr/C5JGrfFwGtm9j7wNDAf+ITQ4g1z/CxMREQkrCILNzwOXAdsBgq9ZgfUSUgys1uBWwE6dOhASkpKXRw2Ku3fv79Rn//R6PocKS/fcfCfs7mkII23O/03newEXaMy6P1TN5xz/wBGm9kNwDLgz865BH+rEhERKa4iS4BfDvRwzv1Q28WUxjn3HPAcQHx8vEtISPCjjKiQkpJCYz7/o9H1KS4QzGfCC5/zzcHxnDr6p3QqbK7rUw69f+qGmTUFkgjNbx0H/MLMbgF+7Zzb4GtxIiIinoqsbrcDaFbbhYhIzcnZk8maJ68lNSOdmcmDGDhklN8liYQtAWKBYcDTzrnfALcBk81Mw+1ERCQqVKQn6QCw3sz+FyjqTXLO3Vndg5vZQiABONHMUoGHnHMvVHe/Io1Zzp5M9jwzmgsLUpk96qcM0Cp2El1Od85dambNgX8COOfSgVvMLNbf0kREREIqEpLe9b5qnHNufG3sV6SxCgekTgWpbB72HAOGXeJ3SSIlPWdmn3k/PxG5wTm33od6REREjlCR1e1eMrOWwGnOuW11UJOIVMLSHUuZuXYmmXmZnFjgmHzMPvZd+Bz9hmsVO4k+zrm/AH/xuw4REZHyHHVOkpmNAdYDf/Mex5pZrfQsiUjlPPrPR7l3xb1k5GXgcGQ1hUc7nkxG12P8Lk2kVGb2gJkdX8724WZ2aV3WJCIiUlJFhttNBwYCKRAaDmFmZ9RiTSJSAUt3LGXRtkVHtB9y+cxcO5NLztBQO4lKG4G/mtlBYC2QRegeSWcRWtDhI+B//CtPRESkYiEp3zkXMLPItsO1VI+IVNDMtTPL3JaZl1mHlYhUnHPuHeAdMzsLuBDoBOQArwC3OueCftYnIiICFQtJm8zseiDG+1C7E/i0dssSkaMpLwh1bNWxDisRqTzn3HZgu1/H90ZE3A+0dc5d5VcdIiISnSpyn6TJQC9Cy38vJPQXv7tqsygRKV/OntAiDWWZ0m9KHVYjUrfMbK6Z7TKzr0q0jzazbWb2jZndW94+nHM7nHM/rd1KRUSkvqrI6nYHCP217f7aL0dEjiZw4BCZs8Yyufk+Hu14ModcfrHt1/a4VvORpKGbBzwFzA83mFkM8DRwEZAKrPYWGYoBflvi9Tc753bVTakiIlIflRmSzOxJ59xdZvYecMSfrJ1zY2u1MhE5QiCYz4S5qzj+4BXcOfhsHul6TNHy3x1bdWRKvykKSBL1vEBzp3PuT1V5vXNuuZl1LdE8EPjGObfDO8ZrwGXOud8CWi1PREQqpbyepJe973+oi0JEpHw5ezKZNe8lNu/tw6zkm+jXswOAQpHUO865QjMbD1QpJJWhC/BdxONU4Pyynmxm7YHHgDgzm+aFqZLPuRW4FaBDhw6kpKTUYLkNw/79+3VdqkDXrep07apG163yygxJzrkvvO9/r7tyRKQ0OXsy2fPMaO4sSGfwFR8z1AtIIvXYP8zsKWARkBdudM6trYuDO+f2ALcd5TnPAc8BxMfHu4SEhDqorH5JSUlB16XydN2qTteuanTdKq+84XYbKWWYXZhz7rxaqUhEigkHpE4FqWwe9ixD+/f2uySRmhDrfX8kos0Bw6u4vzTg1IjHp3htIiIilVbecLvwGO47vO/h4Xc3UE54EpGaUzIg9Rt+pd8lidQI51xiDe9yNXCWmXUjFI6uA66v4WOIiEgjUd5wu28BzOwi51xcxKZ7zGwtUO7yqiJSNUt3LC1ajKFVYQt+dcxu9l2ogCQNi5m1BR4ChnpNfwcecc4FKvDahUACcKKZpQIPOedeMLOfA8sIrWg31zm3qVaKFxGRBq8iN5M1M7vQOfcP78FgKnZ/JRGppKU7ljL90+kcLDwIwP6YII927MAjXVv4XJlIjZsLfAVc4z2+EXgRuOJoL3TOjS+j/X3g/ZoqUEREGq+KhKSfAnO9v/oZsA+4uVarEmmkZq6dWRSQwg65Q8xcO1Or2ElDc6ZzLrJ79GEzW+9bNSIiIhEqcjPZL4C+XkiiIkMhRKRqMvMyK9UuUo8FzexHzrmVAGZ2IRD0uSYRERGgYj1JmNklQC+ghZkB4Jx7pNwXiUil5Oz9nhMLHFml/FfZsVXHui9IpHbdBswP/wGO0CiFn/hYT6nMbAwwpnv37n6XIiIideioc4vMbDZwLTCZ0HC7q4HTa7kukUYlZ+8udj+dxOQ9+2huzYptaxHTgin9pvhUmUjNM7MmQA/nXF/gPOA851ycc+5Ln0s7gnPuPefcrW3btj36k0VEpMGoyAIMg51zE4B9zrmHgUHA2bVblkjjEQjmc9OCTWzJ78Dp8TN55Ee/oVOrThhGp1admD54uuYjSYPinDsM/Mr7Occ5l+NzSSIiIsVUZLhdeBb5ATPrDOwBOtVeSSKNR87e77njlTV8+b1xTPJ8+vXsAKBQJI3BR2b238AiIC/c6Jzb619JIiIiIRUJSe+ZWTtgBrCW0I1k59RqVSKNQM7e79n9dBL/nR/D7uuXMtILSCKNxLXe9zsi2hxwhg+1iIiIFFNuSPLGjf+vcy4beMvM/gq00Ap3ItUTDkhdClLJTniWkb20MIM0Ht5nyw3h+++JiIhEm3LnJHnjxp+OePyDApJI9UQGpE3DnqXf8CuP/iKRBsT7bHnK7zpERETKUpGFG/7XzK608NrfIlJlgWA+X82+SQFJRJ8tIiISxSoyJ+lnwN1AgZkdJLQMuHPOHVerlYk0MIFgPhNe+Jy9eeN5cuSt9E8c53dJIn7SZ4uIiESto4Yk51ybuihEpCHL2fs978+ZztacS3g6eRT9tUiDNHL15bNFN5MVEWmcyhxuZ2YxZtY64vEFZjbU+6oXH24i0SA0B2k0Vxx4g5cubqVV7KRRM7MbIn6+sMS2n9d9ReXTzWRFRBqn8uYk/Q64PeLxQmAq8GvggdosSqShCAekLgXfsWnYs1xwYaLfJYn47e6In/9SYtvNdVmIiIhIWcobbjcCGBDxONs5N8abZLuidssSqf9KBiQt0iAChOYelfZzaY9FRER8UV5PUhPnXEHE43sgNKsWaF36S0QEQos0PPbSOxxfsEsBSaQ4V8bPpT0WERHxRXk9Sc3NrI1zLhfAOfcBgJm1BVrURXEi9VFgfx4T5q1j8+5TGH3NchL7nul3SSLR5Bwz+5JQr9GZ3s94j8/wrywREZH/KC8kzQEWmdltzrl/A5jZ6cAs4Pm6KE6kvgkNsfsx5x1KYHLyPSRqkQaRks71uwAREZGjKTMkOeeeMLMDwEoza+U17wced87NqpPqROqR8BykUwq+4/LEwfRTQBI5gnPuW79rEBEROZpy75PknJsNzA4v+R0eeicixRVfpGE2/YZf5XdJIiIiIlJFR72ZLCgciZQnkLufPQpIIiIijd6SdWnMWLaN9Owgndu1ZGpSD8bFdfG7LKmCCoUkESldIJjPhJfWM/DQhfw4MUEBSaQSzKwlcJpzbpvftZTFzMYAY7p37+53KSIS5ZasS2Pa4o0E8wsBSMsOMm3xRgAFpXqovCXAATCzYyrSJtLY5Oz9noeffYXNGTmcP/7XCkgileCFj/XA37zHsWb2rr9VHck5955z7ta2bdv6XYqIRLkZy7YVBaSwYH4hM5ZF7d+BpBxHDUnAZxVsE2k0wnOQ7tv3IM9eey4jtUiDSGVNBwYC2QDOufVANz8LEhGpjvTsYKXaJbqVOdzOzDoCXYCWZhbHf+6EfhxwbB3UJhKVIhdp2Dx0FsPP0+91IlWQ75wLmFlkm24mKyL1Vud2LUkrJRB1btfSh2qkusqbk5QETAROAZ6IaM8F7qvFmkSiVsmAFDfiar9LEqmvNpnZ9UCMmZ0F3Al86nNNIiJVNjWpR7E5SQAtm8UwNamHj1VJVZV3n6SXgJfM7Ern3Ft1WJNIVAoE8/lwzv2MUUASqQmTgfuBH4BXgWXAo75WJCJSDeHFGbS6XcNQkdXt/ur9ta9r5POdc4/UVlEi0SYQzGfCC5/zdc5YTv3xdZz/o5F+lyRS353jnLufUFASEWkQxsV1UShqICoSkt4BAsAXhP7iJ9Ko5OzdxfpnJ5G5/zr+kjyc87VIg0hN+KM39/VNYJFz7iu/CxIREQmrSEg6xTk3utYrEYlCOXt3kfX0aC4o+DfPXHQz/RWQRGqEcy7RC0nXAM+a2XGEwpKG3ImIiO8qsgT4p2bWp9YrEYky4YB0SsG/2Tx0Fv0TLvO7JJEGxTmX6Zz7M3AboXsmPehzSSIiIkDFepJ+BEw0s/8jNNzOAOecO69WKxPxUcmApEUaRGqWmZ0LXAtcCewBFgG/9LUoERERT0VC0o9rvQqRKBII5nPnK//knoICBSSR2jOXUDBKcs6l+11MWcxsDDCme/fufpciIiJ1qLybyR7nnMshdF8kkUYhsG83E1/ZxFffNyV9/N8Y2buz3yWJNEjOuUF+11ARzrn3gPfi4+Mn+V2LiIjUnfJ6kl4FLiW0qp0jNMwuzAFn1GJdInUuZ+8udj89mhvzO3Nc8guM1CINIjXOzF53zl1jZhsJfZYUbUJDuUVEJEqUdzPZS73v3equHJG6tXTHUmaunUlmXiYnFjjuPGYf3QZPI04BSaS2TPG+X+prFSIiIuWoyOp2mNlYM/uD96UPNmkQlu5YyvRPp5ORl4HDkdUUftPxZNK7Het3aSINlnMuw/vxdufct5FfwO1+1iYiIhJ21JBkZo8T+svfZu9ripn9T20XJlLbZq6dycHCg8XaDrl8Zq6d6VNFIo3KRaW0aaEgERGJChVZ3e5iINY5dxjAzF4C1gH31WZhIrUtMy+zUu0iUn1m9l+EeozOMLMvIza1Af7hT1UiIiLFVSQkAbQD9no/t62lWkTqTM7eXRxb2JK8mANHbOvYqqMPFYk0Gq8C/w/4LXBvRHuuc25v6S8RERGpWxWZk/RbYJ2ZzfN6kb4AHqvdskRqT/hGsffuTqO5NS+2rUVMC6b0m1LGK0WkupxzAefcTufceG8eUpDQKnetzew0n8sTEREBKtCT5JxbaGYpwACv6R7nnMYjSb0UDkinFPybnMGzeKTbsUWr23Vs1ZEp/aZwyRmX+F2mSIPn3aT1CaAzsAs4HdgC9PKzLhERESj/ZrLnOOe2mlk/rynV+97ZzDoBe72/AorUC5EBafPQWcSNuJo4UCgS8cejwAXAR865ODNLBG7wuSYRERGg/J6kXwKTgD+Wsb29mW1wzt1Y82WJ1KxAMJ85Lz7H5IiAJCK+ynfO7TGzJmbWxDn3iZk96XdRIiIiUP7NZCd53xPLeo6ZfVAbRYnUpMCBQ0yYu4rNe/tx/uUfMSQ+1u+SRASyzaw1sBxYYGa7gDyfaxIREQHKH253RXkvdM4tds6NqvmSRGpOzt5dfPfM5Rx78EpmJSczpGcHv0sSkZDLgIPAL4BkQiunPuJrRSIiIp7yhtuN8b6fDAwGPvYeJwKfAotrsS6RasvZu4tdT4/m7IJv+VXiKcQpIIlEDedcZK/RS74VchTeAhNjunfv7ncpIiJSh8obbncTFA2p6+mcy/AedwLm1Ul1IlUUDkinFXzLpqGziBtxjd8liQhgZrmElvwuavIeG+Ccc8f5UlgZnHPvAe/Fx8dP8rsWERGpOxW5meyp4YDk+R7QvSwkagWy95ClgCQSlZxzbfyuQURqz5J1acxYto307CCd27VkalIPxsV18bsskUqryM1k/9fMlpnZRDObCLwPfFQTBzez0Wa2zcy+MbN7j/4KkfIFgvlMfPkrNuSfooAkEuXM7EdmFh61cKKZdfO7JhGpuiXr0pi2eCNp2UEckJYdZNrijSxZl+Z3aSKVdtSQ5Jz7OTAb6Ot9Peucm1zdA5tZDPA08GOgJzDezHpWd7/SeB06kMOUOf+PrzLzaHvdHAUkkShmZg8B9wDTvKbmwCv+VSQi1TVj2TaC+YXF2oL5hcxYts2nikSqriLD7XDOvQ28DWBmQ8zsaefcHdU89kDgG+fcDm+/rxFa7WhzNfcrjVDO3l10Xf0g9x4uIHX8B4zUIg0i0e5yIA5YC+CcSzczDcUTqcfSs4OVaheJZhUZboeZxZnZ781sJ6ElWrfWwLG7AN9FPE712kQqJbxIw+mHvyM49NeM7N3Z75JE5OgOOecc3iIOZtbK53pEpJo6t2tZqXaRaFbefZLOBsZ7X7uBRYCVd3PZ2mBmtwK3AnTo0IGUlJS6PHxU2b9/f6M+/9IcOpBD19UPcvrh73in4y84OaZxv0fKo/dP+XR96tzrZvYs0M7MJgE3A8/7XJOIVMPUpB5MW7yx2JC7ls1imJrUw8eqRKqmvOF2W4EVwKXOuW8AzOwXNXjsNODUiMeneG3FOOeeA54DiI+PdwkJCTVYQv2SkpJCYz7/kgLBfNY+eTVdD3/HpqHPcHJMB12fcuj9Uz5dn7rlnPuDmV0E5AA9gAedcx/6XJaIVEN4FTutbicNQXkh6QrgOuATM/sb8Bqh+1jUlNXAWd5qRmnesa6vwf1LAxYI5jPhhc/J3H8tz4ycSP/EceoFEKlnvFD0IYCZNTGzZOfcAp/LEpFqGBfXRaFIGoTybia7BFjijRO/DLgLONnMZgFvO+c+qM6BnXMFZvZzYBkQA8x1zm2qzj6lcQjsy+LD56axPWcsf04eTn8t0iBSb5jZccAdhOagvksoJN0B/DewAVBIEmlkdG8liUZHXd3OOZcHvAq8ambHA1cTWra1WiHJ2/f7hO67JFIhgX1ZZD2VxNiCbzlt9DUMVEASqW9eBvYBnwG3APcRGqUwzjm33s/CRKTuhe+tFJ7HFL63EvD/27v3OLvK+t7jnx8BSQiaQDEXwyXEIKdBCCGpXJMmMIRLuIkogdjU0jZYEOOp0mr1Zcdz+hLP4YiMJoBBIuDhaG+Axek5UawpVjSi4ZoEBAKUYIZcZ4eESdgkz/lj1uDOkJnMnszM2pfP+/Var6z9zNpr/+aZy5PvrPU826CkXPVoCfAOKaXNtM8PWtQ/5Uhd6whIR775Eium3cIHps7MuyRJ5RuXUjoeICK+BawFjkwpbc+3LEl56O69lQxJylOPlgCX8tY5IE066/K8S5LUO8WOnZTSTmCNAUmqX763kipVWVeSpIHSvLqZpuVNtGxrYcRBo/i9dZNY+OZ6A5JU/SZGxJZsP4Ah2eMAUkrpXfmVJmmgvWf4EF7ZQyDyvZWUN68kqeI0r26m8eFG1m5bSyLx6utrWXHQj1g8478ZkKQql1IalFJ6V7a9lAy1XwAAFfRJREFUM6W0f8m+AUmqM9efcyxDDhi0W5vvraRKYEhSxWla3sT2nbvffRP7FXlww3dzqkiSJPWHSyaN4YZLj2fM8CEEMGb4EG649HjnIyl33m6nitOyraWsdknqLxFxIXDh+PHj8y5Fqlm+t5IqkVeSVHFGDHn3HttHDR01wJVIqncppQdSSvOGDRuWdymSpAFkSFJFKWxrY/aaTQzetWu39sGDBjP/pPk5VSVJkqR6YkhSxSi0FZn77V/zwqYz+dihH2T00NEEweiho2k8rZFZ42blXaIkSZLqgHOSVBEKm9fz5TvvY+WGw7luzl/RMGEk1+ZdlCRJkuqSIUm5K2xez/qF5/I3xbWc++GHmDFhZN4lSZIkqY55u51y1RGQjiy+yOqpNzHjRFeQkiRJUr4MScpNaUBaMXUhkxpm512SJEmSZEhSPgptRR64/UsGJEmSJFUc5yRpwBXaisy9YxlPF87lmPMv5eTTz8y7JEmSJOkthiQNiObVzTQtb6JlWwuHvLkfbVvPY+GcT3GyizRIkiSpwni7nfpd8+pmGh9uZO22tSQSm/bfyc73LGHH4F/lXZokSZL0NoYk9bum5U1s37l9t7Y3UpGm5U05VSRJkiR1zZCkfteyraWsdkmSJClPhiT1q0Jbkf13Ddvjx0YNHTXA1UiSJEl7Z0hSvyls3sBV3/oPtr16Dgfsd+BuHxs8aDDzT5qfU2WSJElS11zdTv2isHk96xaex9VvDGe/2fewY/Ckt1a3GzV0FPNPms+scbPyLlOSJEl6G0OS+lxHQDqquJpt027hxAkjgVmGIkmSJFUFb7dTnyoNSCun3sKJDbPzLkmSJEkqiyFJfabQVuT5Wz5iQJIkSVJV83Y79YlCW5G5dyyD1y+jccbVTDIgSZIkqUoZkrTPCps3cPcdTazcfCq3zrmMSRNG5l2SJEmS1GuGJO2TwuYNrFt4LlcXV3PSRRdwugFJUg2JiAuBC8ePH593KZKkAeScJPVaR0DqmIN0+skn512SJPWplNIDKaV5w4bt+U2xJUm1yZCkXukckFykQZIkSbXCkKSyFdqKLFj8bY4ovmhAkiRJUs1xTpLKUnj9DeYu/iUrNx3H1Ev/jWmTT8i7JEmSJKlPeSVJPVbYvIE1N03n0JafcuucyQYkSZIk1SRDknqkYw7SMcWnmT99LA2uYidJkqQaZUjSXu2+SMNCTmy4Iu+SJEmSpH5jSFK3CoXNBiRJkiTVFUOSulRoK/LHdz/JL3YcbUCSJElS3XB1O+1RYfMGPn33Q6xYN5hPzFnAic5BkiRJUp0wJOltCps38OrC8/h8cQtXXPFjzjIgSZIkqY54u5120xGQxhafpzC1kbPef3jeJUmSJEkDypCkt5QGJOcgSZIkqV4ZkgS0L9Lw829ea0CSJElS3XNOkii0FZl7xzJefu0ybm+4gskzLsm7JEmSJCk3hqQ6V9i8gZ988zM8/9ol3DxnGpNdpEGSJEl1zpBUxzrmIJ1ffJ4xMz/MHxiQJEmSJOck1avOizT8wR/OyrskSZIkqSIYkuqQq9hJkiRJXTMk1ZlCW5HP3/VD3lVcb0CSJEmS9sA5SXWksGULc+96nJXrh3Pp5f/OmScclXdJkiRJUsUxJNWJjlvsGna8n+vmfIUzXaRBkiRJ2iNvt6sDpXOQpk1roMGAJEmSJHXJkFTjSgPSqjO+wcSz5+RdkiRJklTRvN2uxjSvbqZpeRMt21oYcdBILn95I39cfNGAJEmSJPWQIamGNK9upvHhRrbv3A7Aq6+3sOCQQRTf+3GuMSBJkiRJPeLtdjWkaXnTWwGpw679dnJ/8Vc5VSRJkiRVH0NSDWnZ1lJWuyRJkqS3MyTVkBFDRuyxfdTQUQNciSRJklS9DEk1otC6kdlrNjJ4167d2gcPGsz8k+bnVJUkVbeIuDAiFhUKhbxLkSQNIENSDSi0Ffmn27/M3MJL/MmhFzN66GiCYPTQ0TSe1siscbPyLlGSqlJK6YGU0rxhw4blXYokaQC5ul2VK7QVmXvHMla2nslxsy7kmtOmc03eRUmSJElVzCtJVazQupEVX7uIrWt/w61zpnDKadPzLkmSJEmqeoakKlVo3cirC85lyo5l/K8ZB9EwYWTeJUmSJEk1wZBUhToC0tji86w64xtMapidd0mSJElSzcglJEXEhyNiRUTsiogpedRQrToHpIlnz8m7JEmSJKmm5HUl6SngUuChnF6/Km0rJq6++9dsemN/A5IkSZLUT3JZ3S6ltAogIvJ4+apUaN3IgkdaeW7bgWy98j5OOc43iJUkSZL6g0uAV4FC60ZaFpzH5954B29c+U80GJAkSZKkftNvISkiHgT29L/5z6eUvl/GeeYB8wBGjhzJ0qVL+6bAKrGjbStH/PJvOWbXCzw56pMctv5pli59Ou+yKtLWrVvr7vujHPZP9+wfSZLUod9CUkqpoY/OswhYBDBlypQ0ffr0vjhtVei4gjRu1wusPOMbHHbAGOrp8y/X0qVL7Z9u2D/ds38kSVIHlwCvUIW2IqtuuZJxxedY6SINkiRJ0oDJawnwD0bEGuBUoDkiluRRR6UqtBWZe8cyGrddxoqpCw1IkiRJ0gDKa3W7+4D78njtSldo3cj3br+Rla1TuXXORZw4YWTeJUmSJEl1xdXtKkihdSNrF5zPVcVnOeGCWZxqQJIkSZIGnHOSKkRHQHpv8VlWnvENTj11at4lSZIkSXXJkFQBOgck5yBJkiRJ+TEk5azQVuSmxfdwVHG1AUmSJEmqAM5JylHh9R3MXfwIKzeO56wPLWXaScflXZIkSZJU97ySlJNC60bW3DSDw1se5NY5kw1IkiRJUoUwJOWgYw7S+4pPM+8Px9PgKnaSJElSxTAkDbDdF2n4OhPP/mjeJUmSJEkqYUgaQIUtBQOSJEmSVOFcuGGAFNqKzL3rCWbtOIYd064zIEmSJEkVypA0AAqtm/jru37EynXv5Lo5X2Wic5AkSZKkimVI6meF1k2sXXAeXyyu5yOzl3KmAUmSJEmqaM5J6kcdAem9xWdZf/qXOPP4I/MuSZIkSdJeGJL6SWlAWnl6ExNn/lHeJUmSJEnqAUNSPyi0Ffn32z5lQJIkSZKqkHOS+lihrcjcO5bxwmuXMObsDzF5+sV5lyRJkiSpDIakPlRo3cRPb/skL732QW6acwaTXaRBkiRJqjqGpD7SMQfpnOKzjG64zIAkSZIkVSnnJPWBzos0TJ5xSd4lSZIkSeolQ9I+chU7SZIkqbYYkvZBoa3I9Xf9mKHFjQYkSZIkqUY4J6mXClsKzL3rCVauG8rls3/CWccfkXdJkiRJkvqAIakXOm6xu3jHOK6bczNnuUiDJEmSVDMMSWUqnYO0Y+q1TDQgSZIkSTXFOUll2H2RhpuZOHNu3iVJkiRJ6mOGpB4qvP4G/7ngIgOSJEmSVOO83a4HCm1F5i7+JaO2z+SaafMMSJIkSVINMyTtRaF1E19bfDcrN76P6+Zc7RwkSZIkqcYZkrpRaN3Ebxecz2eLz3HmpT9hmgFJkiRJqnnOSepCR0AaX/wNz5x+E9MmH593SZIkSZIGgCFpD0oDkos0SJIkSfXFkNRJoa3IPbd/1YAkSZIk1SnnJJUotBWZe8cyVraexqQLzuPUU8/IuyRJkiRJA8wrSZlC6yZWfe0CimtXcOucKQYkSaphEXFJRNweEX8fETPzrkeSVFkMSfxuDtLkHY9ww/ShNLiKnSRVrIhYHBHrIuKpTu3nRsQzEfFcRHy2u3OklO5PKf058HHg8v6sV5JUfer+drvSRRpWnX4zE8/+aN4lSZK6dyewALi7oyEiBgELgbOBNcAjEfEvwCDghk7PvyqltC7b/0L2PEmS3lLXIalzQDrBRRokqeKllB6KiLGdmj8APJdSWg0QEd8DLk4p3QBc0PkcERHAV4D/m1Ja3r8VS5KqTd2GpEJbkT/7znKufeNAimcYkCSpyo0BXi55vAY4uZvjrwMagGERMT6ldFvnAyJiHjAPYOTIkSxdurTvqq0RW7dutV96wX7rPfuud+y38tVlSCq0buLjdz/CY6/u5M0r/54TjhuVd0mSpAGUUvo68PW9HLMIWAQwZcqUNH369AGorLosXboU+6V89lvv2Xe9Y7+Vry5CUvPqZpqWN9GyrYURQ0bwkTWb+XTrIFqv/D4NBiRJqgWvAEeUPD48a5MkqWw1H5KaVzfT+HAj23duB+DVtle5/ZBdXDXuQv7iuNE5VydJ6iOPAMdExNG0h6PZwJX5liSpO/c/+go3LnmG37a28Z7hQ7j+nGO5ZNKYvMuSgDpYArxpedNbAanD9v324743nacrSdUoIr4L/Bw4NiLWRMSfppTeBD4BLAFWAf+QUlqRZ52Sunb/o6/wuXuf5JXWNhLwSmsbn7v3Se5/1AvAqgw1fyWpZVtLWe2SpMqWUrqii/Z/Bf51gMuR1As3LnmGtuLO3draiju5cckzXk1SRaj5K0mjhu55zlFX7ZIkSepfv21tK6tdGmg1H5LmnzSfwYMG79Y2eNBg5p80P6eKJEmS6tt7hg8pq10aaDUfkmaNm0XjaY2MHjqaIBg9dDSNpzUya9ysvEuTJFW4iLgwIhYVCoW8S5FqyvXnHMuQAwbt1jbkgEFcf86xOVUk7a7m5yRBe1AyFEmSypVSegB4YMqUKX+edy1SLemYd+TqdqpUdRGSJEmSVFkumTTGUKSKVfO320mSJElSOQxJkiRJklTCkCRJkiRJJQxJkiRJklTCkCRJkiRJJQxJkiRJklTCkCRJUhd8M1lJqk+GJEmSupBSeiClNG/YsGF5lyJJGkCGJEmSJEkqYUiSJEmSpBKGJEmSJEkqESmlvGvosYhYD7yUdx05OgzYkHcRFcz+6Z790z37B45KKb077yIqkeNPl/y56R37rffsu96x3/asy3GvqkJSvYuIX6WUpuRdR6Wyf7pn/3TP/pHK589N79hvvWff9Y79Vj5vt5MkSZKkEoYkSZIkSSphSKoui/IuoMLZP92zf7pn/0jl8+emd+y33rPvesd+K5NzkiRJkiSphFeSJEmSJKmEIanKRMSHI2JFROyKCFcpASLi3Ih4JiKei4jP5l1PpYmIxRGxLiKeyruWShMRR0TETyJiZfZzNT/vmqRq47hUHses8jmO9Y5j3L4xJFWfp4BLgYfyLqQSRMQgYCFwHjABuCIiJuRbVcW5Ezg37yIq1JvAp1NKE4BTgGv9/pHK5rjUQ45ZvXYnjmO94Ri3DwxJVSaltCql9EzedVSQDwDPpZRWp5TeAL4HXJxzTRUlpfQQsCnvOipRSmltSml5tv8asAoYk29VUnVxXCqLY1YvOI71jmPcvjEkqdqNAV4uebwGfwGoFyJiLDAJWJZvJZJqmGOWcuEYV7798y5AbxcRDwKj9vChz6eUvj/Q9Ui1LiIOBv4Z+FRKaUve9UiVxnFJql6Ocb1jSKpAKaWGvGuoIq8AR5Q8Pjxrk3okIg6gffC4J6V0b971SJXIcanPOGZpQDnG9Z6326naPQIcExFHR8Q7gNnAv+Rck6pERARwB7AqpXRT3vVIqnmOWRowjnH7xpBUZSLigxGxBjgVaI6IJXnXlKeU0pvAJ4AltE9I/IeU0op8q6osEfFd4OfAsRGxJiL+NO+aKsjpwB8BZ0bEY9l2ft5FSdXEcannHLN6x3Gs1xzj9kGklPKuQZIkSZIqhleSJEmSJKmEIUmSJEmSShiSJEmSJKmEIUmSJEmSShiSJEmSJKmEIUn9JiJ2liw5+VhEjI2IKRHx9R489+Hs37ERceU+vPaKiHg8Ij4dEftlH3urhog4MCIezI69PCKmZs95LCKGlPu6AyUi/jIino6IJ7PP76bsDeN6e76xEfFUtt+jr1E35/qb3j5XkiqB41f/cfxStXAJcPWbiNiaUjp4H88xHfhMSumC3r52RIwA/g/ws5TS33Y67hTg7zreTT4ibgP+I6X0v3v4OkH7z9GucurbFxHxceASYHZKqTV7Q8K/BG5JKW3pdOyglNLOHpxzLPCDlNL7+6C+ff66S1KeHL/6h+OXqkpKyc2tXzZg6x7aptP+ywygEVgMLAVWA5/s/FzgF0ABeAz4r8Ag4Eba37X8CeDqnrw2MA7YCERHDcAI4LmS818NbAJeAO7Jnnd9yWt9KWsbCzwD3A2sAI4CZtL+RnfLgX8EDs6OfRH4Utb+JPBfsvaDgW9nbU8AH8ra93ieTp/Ly8DR3fU78FXgceAM4IvZ5/AUsIjf/XFkcnbM41mfPrWHr9HQ7Gv0S+BR4OKs/WPAvcD/A54F/mfW/hVgZ9af9+T9Pejm5ubWm83xy/Er7+9Bt/y33Atwq92t5JfNY8B9WVvnQeZh4EDgsGwQOCD72NbOx2eP5wFfyPYPBH61p1+4nQeZrK0VGNmphs7nvxO4LNuf2fFLmfZbU38ATMsGmV3AKdlxhwEPAUOzx38NfDHbfxG4Ltu/BvhWtv8/gJtLXveQ7s5Tcty7gM176fcEfKTk8aEl+98BLsz2nwCmZftdDTJfBj6a7Q8HfpMNPB+j/T8Gw4DBwEvAEV31vZubm1s1bY5fjl9ubvsj9Z+2lNKJezmmOaW0A9gREetoHwTWdHP8TOCEiLgsezwMOIb2v571tZnZ9mj2+ODstf4TeCml9Ius/RRgAvCz9rsXeAftf03rcG/276+BS7P9BmB2xwEppc0RccFezvM2EXEO7QPWcODKlNLDtA/u/1xy2IyI+CvgIOBQYEVE/BQYnlJ6KDvmO8B5XfTBRRHxmezxYODIbP/HKaVCVsdK2v8i+XJ39UpSlXD8auf4pbplSFLedpTs72Tv35NB+1+2lpTzIhExLjv/OuD3e/o04IaU0jc7nWsssK3TcT9KKV3RxXk6Pse9fX57Ow8ppS0RsTUijk4pvZD1w5KI+AHtgxLA9pTdxx0Rg4FbgCkppZcjopH2gaKngvZbKZ7ZrTHiZMr/2klSLXH86vl5HL9UdVzdTpXuNeCdJY+XAH/RsRJORLwvIoZ2d4KIeDdwG7AgpVTOSiVLgKsiomMC7ZhsEm1nvwBOj4jx2XFDI+J9ezn3j4BrS2o8pIzz3ADcGhHDs+OCrgeOjvYN2edxGUBKqRVojYgzso/P6eL5S4DrstcgIibt5fMCKO7LSkWSVCMcv97O8UtVw+SsSvcEsDMiHqf9fusm2u+pXp794ltP+0o5nQ2JiMeAA4A3ab8cf1M5L5xS+mFE/D7w8+x37Fbgo7T/1an0uPUR8THguxFxYNb8Bdrvf+7K3wELs2VLd9I+qfbeHp7nVtrvq14WETuyun7G726rKK2tNSJup33SawvtE2A7/AmwOCIS8MMu6vzvwM3AE9G+BO0LwN5WalqUHb88pdTV4CVJtc7xy/FLVcwlwCVJkiSphLfbSZIkSVIJQ5IkSZIklTAkSZIkSVIJQ5IkSZIklTAkSZIkSVIJQ5IkSZIklTAkSZIkSVIJQ5IkSZIklfj/HU9H+hRq8fYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "min_g = np.min(g_discrete)\n",
    "max_g = np.max(g_discrete)\n",
    "\n",
    "fig = plt.figure(figsize=(12,6))\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "plt.plot([min_g, max_g],[min_g, max_g],label='y=x comparison')\n",
    "plt.plot([min_g, max_g],[m*min_g+b, m*max_g+b],'--',label='Best fit')\n",
    "plt.plot(g_discrete,dJ_du[idx],'o',label='Adjoint comparison')\n",
    "plt.xlabel('Finite Difference Gradient')\n",
    "plt.ylabel('Adjoint Gradient')\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.axis(\"square\")\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "rel_err = np.abs(np.squeeze(g_discrete) - np.squeeze(dJ_du[idx])) / np.abs(np.squeeze(g_discrete)) * 100\n",
    "plt.semilogy(g_discrete,rel_err,'o')\n",
    "plt.grid(True)\n",
    "plt.xlabel('Finite Difference Gradient')\n",
    "plt.ylabel('Relative Error (%)')\n",
    "\n",
    "fig.tight_layout(rect=[0, 0.03, 1, 0.95])\n",
    "fig.suptitle('Resolution: {} Seed: {} Nx: {} Ny: {}'.format(resolution,seed,Nx,Ny))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}