python/examples/adjoint_optimization/01-Introduction.ipynb

{
 "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.verbosity(0)\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 = 20\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.1\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",
    "    # \n",
    "    # The commented lines above impose symmetry by overlapping design region with the same design variable. However,\n",
    "    # currently there is an issue of doing that; We give an alternative approach to impose symmetry in later tutorials.\n",
    "    # See https://github.com/NanoComp/meep/issues/1984 and https://github.com/NanoComp/meep/issues/2093\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",
    ")"
   ]
  },
  {
   "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": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python3.9/site-packages/meep/visualization.py:357: UserWarning: The frequency parameter of plot2D has been deprecated. Use the frequency key of the eps_parameters dictionary instead.\n",
      "  warnings.warn(\"The frequency parameter of plot2D has been deprecated. \"\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAToAAAEGCAYAAAD1+lmKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAsW0lEQVR4nO2de3Rc1X3vP1uaGWlmNHrLRsSAKCUhiUMNsXmkPMszBEgCYV3yapJmhWblUnKbNOkjq417e9MbkqzQNKu3vWnDbZuSkBIgpcE0MWAnIdggG0gwOMbY2BhkI2ust2ZG89j3j9E+jGQ9RtIZzd5nfp+19rJHmvnNT2ef3/fs909prREEQQgyddV2QBAEodKI0AmCEHhE6ARBCDwidIIgBB4ROkEQAk+o2g4sBhVXmpZqeyEIgjUcZkBr3bXg+7TWzhS60X5yx7Y7tNqo9B3b7vDVrtZasxG95aUtuvMrnXrLS1t8t6+1dtJ+oVDQk5OT0+zf+rVbdVtbmwY0oMPhsFZKea8XU5RSOhwOF1/3oNUfK33b12+b5sPk5KQuFArL/ltcvP4L2WejfzFWyfgy9oEduhztKOdNthQ/ha7SlcBGAhcEfmCErtT+N77xDf+FrgfN59CJtyX0N7/5zWk++CF0rl7/hez7JXQrIXJqoxKhK+ciVfRJM9WiqwQuB1mhUNCb926eZt9voas/rV7zuWKLrr293Xehc/n6L2TfD6FbKZGTFl2ZF6kSeE8aH7sApbgeZI/uf1R33j7dvp9CZ1py9BRf+y10rl//hewv975dSZHTWovQlXOR/Gbak6YCQlftIPDL/ua9m6f93DehmxK5+tPqvZ/5KXRBuf7z2V/OfbvSIqe1CF1ZF8lPjnvS+Cx0NgSBH/Yf3f+oNxlh8EXoSlpy3mSEj0IXlOu/kP2l3rfVEDmtRejKukh+MeuTxkehsyUI/LBfOutqWLbQlYjctFlXn4QuSNd/IZZy31ZL5LQWoSvrIvnBnE8an4TOpiDww77vQjdjTM5voQva9V+Ixd631RQ5rUXoyrpIy2XeJ40PQmdbEPhh31ehmyFyfgtdEK//Qizmvq22yGktQlf1Sliu0NkYBH7Y903oZhE5P4UuqNd/Icq9b6sdX4aaFjobKmE5QmdrEPhh3xehm0Pk/BK6IF//hSjnvrUhvgw1K3S2VMJShc7mIPDD/rKFbh6R80Poqn19qm1/ofvWlvgy1KTQ2VQJSxE624PAD/vLEroFRG65QmfD9am2/fnuW5viy1BzQmdbJSxW6FwIAj/sL1noyhC55QidLden2vbnum9tiy9DTQmdjZWwGKFzJQj8sL8koStT5JYqdDZdn2rbn+2+tTG+DDUjdLZWQrlC51IQ+GF/0UK3CJFbitDZdn2qbX/mfWtrfBlqQuhsroRyhM61IPDD/qKEbpEit1ihs/H6VNt+6X1rc3wZAi90tlfCQkLnYhD4Yb9coVuKyC1G6B7d/6iV16fa9s19a3t8GawXOuAkYAvwPPAc8OkFP9PtTiXMJ3SuBoEf9ssRutLz5BYjcuUK3czz8PzG5uu/EGzEifgyuCB03cDZU/9PAC8Ab5n3M93uVMJcQudyEPhhf0GhW2JLrlyh2/LSFt15e/EUlUpg+/VfCDbiRHwZrBe64xyB/wCumPc93e5UwmxC53oQ+GF/XqGbErn60+r9yRkxQ+iM/5v3bvYlZ8RMXLj+C9k3LbpKUIlGilNCB/QALwPNs/zuFmAHsIMWdyphptAFIQj8sD+X0DWtbZp2npzfQlfqv1/JcUpx5fovZL9SJ2NXqifmjNABTcBO4IYF3+tzFjBDRZ40JTdMUILAD/uzCd2tX7tVq88rr7vqt9Dd9vXbpvnvt9C5dP0Xsl8JoavkcJMTQgeEgR8Dnynr/RUQuoo9aaZumCAFgR/MFLotL23R8b+MF1t0+C90ibcldNNfNk3z30+hc+36L2Tfb6Gr9Ji69UIHKOBfgb8p+zOS13UaLtqXvK5225e8rv4L3QVTN+evgGemyjXzfkbyunq4GmSS19Vu+5LX1YIieV2LuBxkktfVbvuS19WCInld3Q8yyetqt33J62pBkbyuwQgyyetqr33J62pBkbyu7geZ5HW1277kdbWgSF7XLb7ZrJZ9yetqt33J62pBkbyu/rPS9iWvq932Ja+rBUXyuvpLNexLXle77UteVwuK5HX1j2rZl7yudtuXvK4WFMnr6g/VtC95Xe22L3ldLSiS13X5VNu+5HW1277kdbWgSF7X5WGDfcnrard9yetqQZG8rkvHFvuS19Vu+5LX1YIieV2Xhk32Ja+r3fYlr6sFRfK6Lh7b7EteV7vtS15XC4rkdV0cNtqXvK5225e8rhYUyetaPrbal7yudtuXvK4WFMnrWh4225e8rnbbD2peV1V8rxuoE5Xm96vthSAI1rCRnVrr9Qu+rxw1tKXQ7c7x5JU6eNN1KnJ6iQ/LS4Qirp2MTZktujrfFbbC3HPTPVzSc4nvdrce2MpN99xUMfuC4ApBjC/nhC6IlSAIQcaG+HJO6PzGhkoQhKBiS3zVtNDZUgmCEEQqHV9bD2wt+701K3QicoJQOVZC5G6656ay31+TQmfTk0YQgsZKidw9N91T9mdqTuhse9IIQpBYSZFbjP2aEjobnzSCEBRsFTmoIaGzuRIEwXVsj6+aEDrbK0EQXMaF+Aq80LlQCYLgKq7EV6CFzpVKEAQXcSm+Ait0LlWCILiGa/EVSKFzrRIEwSVcjK/ACZ2LlSAIruBqfAVK6FytBEFwAZfjq6pCp5S6UynVr5TatVxbLleCINiO6/FV7RbdPwNXL9eI65UgCDYThPiqqtBprX8GHFuOjSBUgiDYSmDiq5zz1itZgB5gV1nvnZHX1eZsSy7ljAAplS6uUHrf2hxfBlxJd7iQ0AG3ADuAHbS4UwkidFJcFjrb48sQGKGb9t5udyrBJaFbyWiULGB2Y7KA2R5fhnKFrtqTEYsmMGMGgmApQYyvai8v+R6wDXiTUuoVpdTHF/pMECtBEGwiiPEVWtFvm4HW+v2L/UwQK0EQbCKI8eVc1zWIlSAIQcaG+HJO6PzGhkoQhKBiS3zVtNDZUgmCEERsyrZXs0InIicIlcO2bHs1KXQ2PWkEIWjYmG2v5oTOtieNIAQJW9e51pTQ2fikEYSgYKvIQQ0Jnc2VIAiuY3t81YTQ2V4JguAyLsRX4IXOhUoQBFdxJb4CLXSuVEKtUzyE4nWUUtTX1y/JVn19PUop73Vd3fG3eKFQmO80HaFMXIqvqu51rSQuVUItMZugzCZ0SqlZRcq8d7bPlH525s9nfn+hUPDsG1szPyvMjWvxFUihc60Saom6urrjBCwSiUx7nclkFm3XiFU2m53288nJyWnfp5SisbFx0faF13ExvgLXdXWxEmqNmd3SmULnJ5lMZkndYOnGzo6r8RWoFp2rlRBE5hr3KhQK5PN5IpEIdXV1pFIpjh2bnh+pubkZpdSsYlMoFMjlct44GxRbiWZcz4zRjYyMoLWmo6ODdDpNKpUiGo0yNjZGJpMhHA4f13Wtr6+noaHhuHE+oYjT8VXOMcS2lJnJcUqx7fjnWj9KPZ1OVz1Pw1JKX1/ftL+j1o5cn+u+tS2+DAT1KPXZcPpJE1BmjpW5wsDAwLTXWrqwgYgv57uuQaiEoPPcc8+jlCKXy6GUIhQKMTIywnPPPcfjjz/O9u3bOXToEC0tLaxbt47169dzxhln0NbWxje+cTo/+tHJXHPNAT7+8ac5cuQI+/btY+/evRw9epTGxkZOOukkenp6WLVqFZFIhFQqxfj4OJOTk4yPj3Ps2DGGhobIZDIUCgVCoRChUIhYLMaOHR9lz57LgGJXNZfLeX6b1kAtE5T4clroglIJQcdMBsxc39bY2EhTU5M3Czo5OUk6nWZsbIzBwUH+8R9/i82bT+ayy/bw/vfvYGQkx9jYGKlUyhuny+VyZDIZxsbGiEQiRCIRb0wuk8mQy+WIRCK0tLSQzWbJ5/PU19cTDod58snfZc+ei3nzm7ewe/fvABAOhz0fa325SZDiy1mhC1IlBJFCoeD9P5PJoJQik8mQz+fRWjM2Nsb4+LjXyjMTA5lMhqGhIb797bN48snTufDCZ7nxxl8wOdlAoVBAKUVDQwNNTU3kcjnq6urQWjM6OkoulyMUCqG1Jp/PA8VWmZlkgKLA1tfX09v7EZ577iI2bOjl8ssfnlXoapmgxZeTQhe0SggCWutprZ9SoUun02itmZiY8LqUIyMjHDlyhOHhYbLZLKFQCKUUqVSKzZvfzUsvreess7Zx+eUPk0o1EQ6HCYVCRKNRWlpayOfzxONxMpkMWmuGhoYYGBjw1snF43Gi0SjhcJhcLsfk5CS5XI5cLscTT3yYPXsu4pxzdvCBD2wjn+/0fK3kUhdXCGJ8OSd0j+x7hJvvu5m7b7ibC0+60Hty+8XWA1t9s++3b5XCrDJbjr9mLCsUKt5SpWNdpos5Pj7O4OAgw8PDDA0N8dprrzE8PEw+n/eWk7zwwm0MDb2TNWse4Dd/8y76+tpZvXo10WiUSCRCQ0MDzc3NhEIhmpubGR4eZmBggKNHj3oTCR0dHXR3d3tLRUwLb3Jykt7ej7B//2WcddY2PvShnXR0dE0T6JlCV4td16CJHDgodEaELj7lYl8HipVS00SuXPsmQOda91X6s9neU/qz2WzNZXc+XxZDaSAv53rObNGViqYZKxsfH2doaIhkMumViYkJbzvW0ND/Ip3+AC0td9HW9hUOH24BIJFIeF3cSCRCPB4nFosxOTmJUopkMsnw8DCvvvoqWhe3d7W2tlIoFKivryefz1NXV8eOHR9l//6rOOOMR7n++p/T1HQKra2t07qr0nUNZkpR54SukpVw8303+2p/qRvTq4VpjflBacvo0KFDAJ7IDQ4OMjo6SjqdJp/PEwqFGBv7Mun079La+l1OOeVr1NeHve6q6daaRcHhcNj7f0NDA42NjTQ2NtLQUBzHM4uBzULicDjM9u0fYt++i3jTmx7h6qsfJBYrtvjC4fA0cZttf22tETSRAweFLoiVEETS6bT3/6effppYLMbExASpVIp0Os3k5CRQbEEdPvwFxsbeR2fnv7N27T8Rj6/xuqhtbW20trYSi8UIh8Ne669QKHii1tzcTHd3t9d6a2trI5FIeBMVDz/8HnbtOo+zz97Otdf+nHi8m5aWFsLhMPl8ftp4Yq0vJ6kENsSXc0L37LPPzvrzpXTbAHqP9vJHT/wRXzv3a3SMdsxpfyn4aauSvG3q3+X4awSoqamJfD7Pjh07gJsB2LlzJ7FYjFAo5M16mpbTr399K/3913LiiT/krLO+Q2vrG2htbaW5uZmmpibi8TiJRILm5mYaGhrIZrNorb1ua0NDA+3t7dTV1dHZ2UkulyMcDnut6Yceehe/+tV5nH/+03zkI78iHj/d235mhHPm2jnBP2wQOXBQ6P70T//UN1vJRJKnTnuKs/edzX2/uI/7uM8322zw19dK8qOpf5fjrxE6Ix7F/atFoXv22WeJxWI0NTWRSCRIJBLEYjF27fokL798DT09mzj77H8jkeikvb2drq4u2traaGoqzraa7qnZg2pac4VCwft5PB4nm816JZ1O89BD72LXrgvYsKGXj33sl3R3F0W0rq6OTCbD5OQk+Xx+2niiCJ1/2JRtzzmhe/DBB/0x1APcBNwF2w9s98dmKRt89HWFqJS/g4ODjI+Pk81mvVbY3r2f5tChositX//PhMMNRKNRmpqaaGlpob29nWg0ChTFJ5vNMjk56S0qNguCoTgWmkgkUEqRz+dJpVJ873u/za5dG1i37nFuvPExGhtPJBqNEo1GvRahmYkt7WaXdmOFpWNbtj3nhM4XeiiK3D3AgQrZFzwSiQT19fVEo1EaGhrYv/8z9PVdR0/PJs4669vk88VJALOY2Ew0KKVIp9OMj497LTAjeOl02lsgHIvFiMfjNDY2Eg6HufvuC+jtfSvnn/80N9ywjcbGOACjo6Nks1kKhQKpVMoTzlQq5fla2o0VlsZKrsO79POXlvWZ2hO6HiovcjWY1tUsLTFdv0Qiweho8XcnnHACSilisRivvfbn9PW9x2vJ5fN4W7PMGJzZDwswMTFBf38/yWSSsbExstmst+0rn88Ti8Xo6urytpTdddf5/PSnb+SKK/byiU/sA07xFgwPDAx4i4bNMU/5fJ6JiQnv7zCTJIaZy2aE+bF1sbFzQrecJRCFUwoUbihQd18dda/U+f7Xl9ovfLjg63KNijLVilmOv0YMzKklDQ0NntC1t7ejtebQoT+mr+89/MZv/BfnnnsXWofI5/PkcjmvO6m1JhwOE4vFvIMAjh07xssvv0wymfRaZOl0mkKhQHt7O6FQiNbWVr773XfwyCNv5NprX+azn+0jHu8B8BYnm/V2xkczezs+Pu79HTOFTigfW0UOHBS6JXcteoAbgHugcKBAAZ/HYmbYB/e6QX76WzrWFY1G2bPnD0pacv/K6/sx8Bb5mlL6emJigsHBQY4cOcLRo0e9hcOl43OpVIr777+MJ598E+9610E+85mXiMebaG5uBvAOARgdHfVshEIhb/3c2NhYRa5BLWGzyIGDQrckeliZ7mql7DtI6cJbI3InnvhDzjjjW6RS4WkzqGZBsJkNTSaT3phdf38/w8PDpFIpryVmlo+Ycb/iKSTrueiiXXz4w88zMtLojfeZfbBDQ0MMDg4yNDSE1ppoNOotRjYb/qE2t3wtF9tFDmpB6HoQkasCpS2jvr73sGrVDzj55K8zOlrviVokEiEcDnuTCMUdEmMcPHiQw4cPk8vlGBoa4tixY0xOThIKhYhEIiQSCaLRKI2NjTz//KfYu/cy1q17nCuu2MIrrxQnPKLRKPF4HK01yWSSI0eOkEwmGR0d9Q4HMKegxGIxz1fZArY4XBA5CLrQ9SAiVyVKu64dHd/nhBP+NxMTr2/wLz2JpKmpiWg0Sj6fZ3BwcNoRTrlczlsGYpamtLW10dXVRW/vR9i797c588zHuPDCezhy5PW1fLFYjGg0SqFQYGRkhP7+fo4ePcr4+DjxeJy2tjZisRjt7e3T/JbTS8rHFZGDKgudUupq4BsUB2z+SWv9Zd+M9yAiZwknnXQ72WzeEywzswrQ1NREKBSisbHRWx93+PBhBgYGyGQyNDQ0EIlEqK+v91p+0WiUHTs+yi9/+Q7OPXcnF1/8I4aGMoyPj5PP5733mL2vY2NjDAwMMDg4SCqV8iY5jM3SSZiZ+5OlKzs7LokczCN0SqlNwKe01geW/S2z268H/g64AngF6FVKPaC1fn7ZxnsQkasypYLR3NzM5OSkN4FgZlnNaSSmm1napYxEImitiUQi3vYx87sdOz7Kiy++gw0bennnOx9ibOz1lpzZyG8mLbLZLKlUyjsLb2Jigvr6esbGxrxS6uvMo6pkecnxuCZyMH+L7v8BP1FK/QvwFa2139lOzgFe1FrvB1BK3Q28G1ie0PUgImcBpUmiTz75ZE9khoeHgaKAxGIxEokELS0ttLW1eefSRSIRurq6vJOJzWkldXV17N793zlw4EpOP30za9feTV9f2GvBtbQUj3UyS1BmHnVlRK++vp5kMukd9VTazS5daiIcj4siB/MIndb6HqXUQ8CfAzuUUt+B19dkaK2/vszvfgNwqOT1K8C5M9+klLoFuAWAlgUs9iAiVyVmHqpgtm8BnHrqqQwPDzMyMuKNgWUyGZqbm2lvb2fVqlV0dRUPwGxtbWX16tVMTEx4+1bNWN3Wre/jwIHLOOWUB+nu/iq7d0+SSCRYs2YNXV1drF69mlAoxOjoKCMjI97uB621t2jY/My09o4ePTptZ8SoWfyHJMeZiasiBwuP0U0C40ADkAC/F58tjNb6W8C3ANSJSjM8xxt7EJGziFKhW7NmDbFYjEgkQi6X87qL8Xjc29fa1dVFOBymtbXVO6hzcnLSO9rpBz+4hD17NnDGGY/S0/N3HDjQz9DQEJ2dnXR2dhKPx71ThQcGBsjn895ylFwu5230n5iY8HZEpNNpGhoaGBkZ8Xwt3SUBssnf4LLIwfxjdFcDXwceAM7WWk/M9d4l8ipwUsnrNVM/Wzw9iMhVmZktutJlGm1tbd62LTOhYETIiI5JdGOWnZhSV1fHD35wCb29a9mwoZf16+/l5ZfrKRQK3l7V0laf+V6zPq6xsZFIJOKNs5nxwVJfR0b+eiUukbO4LnIwf4vuC8BNWuvnKvLN0AucrpQ6laLA3Qx8YNFWehCRs5xQKOQN+JsdD9lslvHxcQYGBgiHw4yNjR2XqUtrzXe/+w5+8Yu1nHfeU1x11SYGB8PE43G6urq8iYpCoUB/f7+XGMfsZTXbyNLpNMlkkoaGBu9U4VWrVtHd3c3Bg58jnb7W87WpqWma77U+EREEkYP5x+gurNi3Fu3nlFK3Aj+muLzkzkWLag8icpZSKhBmJ4Np9eVyOW+ng1k7Z45CNyUSibBly408/fRvcfbZ27nyyv+iUNDeQZuRSMTrhtbV1XHw4EEOHjxINBqltbWVE044gdWrV3sZw0yy68bGRpqbmzn55JPp7/8LXn75Is477ym2b387UEysU/o31LLQBUXkoMrr6LTWm4BNS/pwDyJyFlO6Ni0Wi01br5bP58lkMkxMTHhjYlprT4Sampp44YXb2LfvfN7ylq1cfPF/ks/HvNnVSCRCR0cH+Xye0dFR+vv76evro7+/n4aGBk4//XQ6OztpaWkhFouRzWZJJBLeTox4PM4rr/wJL7xwEdddd4g/+7Ms559f9LWtrW3a31GrQhckkQNXd0b0ICLnEGZMrHTjvjlXzpxaks1maWxsJJPJsGfPH9Dff9XUycN3MzER9wSnrq7OG38zn83lcgwODtLX10djYyNdXV3kcjkaGhpIJBLeNi/TohsY2MjRo5dz1VX7uP32LK2tp3i+lk6iQG0KXdBEDhwUurrT6lbsqKXl2M+Rq5ljmswaNyNgRZOv73U1Ox2Gh4eZmJjwfmdySJjFwQ0NDRw9+kWSyfexZs0DnHnmnaTTEU/QxsfHp01YAN6Gf3O8k5nAKJ2sMAuPE4kE+/d/hmTyJtate5yPfewVUqnTp8YPu4/zu1YJmsiBg0JXuKGwYkctLde+a0GzHH9Nxi1D6dFHe/fu9TbnDw8Pk8lkPLFqbm4mFovR0NDA889/imTyOk4++UesW3cnUFxvl0qlvAXAZqfCzPG+hoYGVq9e7R3Fbva4aq2ZmJggHA5z8ODnSCav4rTTfszllz9CMnmqd5gAvBEoLi8xxztBbe6MCJrIgYNCJ91VO5mZa8HkcgXYvXs3Y2NjpFIpb8N+Pp8nHo97C4SffvrjHDp0EWec8Shvf/v3yWYbvWPUx8fHGRkZIZPJHJeaMBKJeBMZzc3NXvKdQqHA0NCQd+z6I4+8l337zuUtb9nKpZf+iLq6OMeOHSOdTk+1Qq8Hjt8ZUYtCFzSRAxeF7kAFbPYgIuczpa3DQ4cOeSKVyWTIZrNeiy4ej/PLX37Cy9Z1wQUPMjzc4OVbnZycZHh42DvVxOxFNUtIEokEHR0d3n5ZM+lgkuQA3HvvpezceRZvf/sTXHnlZqDYYkun094iYoOcMOw/1RY5cFDorr/+el/tDTQNsOPUHax/aT2dZ3bCmf7YfYAHfPe1YjzwALC8a2u6rWb9WyqV4ic/Kf7OnChiWkfmwMumpiZ27foke/dewPnnP8311z9KMlnvfcaImjlfrnTiIpPJoLWetovCLBw2S05CoRD3338Z27b9Fhs29PI7v/ND0ulJb/+sWd9XevBm6YGhwvKxQeTAQaH70pe+5Jut3qO9fHb7Z/mH8/6BDV0bfLML8MC9D/jqa0WZErrl+GtELJFIkM/n2blzpyd05uw3M1lhdkHs2fMHHDhwBWedtY1rrnmUVCrrjcmZLF/hcJj29nba29u9CYnh4WGGhoZIpVJEIhEvq5fB7L546KF38YtfrOXCC5/lggt+yOHDrzE2NubNzLa1tdHa2jptbLH0MAKozVlXv5C8rstg7dq1vtjZemArn3/o89z3/vsq86S51z9fVwo//S2dwTWTBIB3kOa2bR/kwIFLOfPMx7jqqk2k0/VkMhnvlJPR0VHy+by3ANh8fmhoiL6+Pm+5imnVmS5uNpulubmZhx9+D729a7nuukO8733P8MwzaW+tXVtbG4lEgng8zhve8IZpLbqZy0uEpSF5XS3ApidNUCltJXV2dnpCpbVm06Zr2L37PM4+eztXXfUQoKYdymladYA3i9rd3U04HKa/v59UKsXQ0BD19fXe+01CaqUUTz31exw6tIHrrjvEF794lP37o2SzWfr7+3n11VfJZrOceuqp3lay0lnWUtEDadEtBcnr6gPLXbLx04M/5eb7bubuG+7mgjUX+L4ExNgHd5aXmJtgOf4WCgVvXAyY1pVsbW31Tvu9//7LeOqpDaxf/yRXX70JrYtjYkakzMGZxpY5Z66jo8PrVg4PD3stPrOGzhzrNDLy14yPX8c55+zgU59KMjnZ4glhJpPxxvQA72TjeDzu+TrzhGFhcdi62Ng5oVvOE3brga2eyFWqEoz9y//tcudaA8vxd+ZnSwf1zZjZ979/Ib29b2ft2p+xfv3dJJPK261gJgVisRgtLS3U19cTCoXo6OigpaXFS4iTy+Xo7u72ZlzT6bR35NK+fX9IMvnfOO20H3PFFY/R1/cbHDt2jNdee807kqm5uZnm5mai0ShKvd6SNMixTEvHVpEDB4VuqU9cI0KVrISZ9l1rHSzHX9MKM5SO0TU2NnLPPRfzxBNn8ta3/pQzzvg7Dh4cJxqN0tHR4a2FM7sa6urqaGpqmjYRkUgkaGxs9FIkxuNxxsbGvNba1q3vI5m8lDe+8WEuvvg+Jiaa2L9/v5czIp1O09TURF1dHZ2dnV7mL7Pn1jDzKHWhPGwWOXBQ6JaC7ZUQBGa26EpF84c/vJzHHnszGzb08ta3/gsHDxbTDpoJAbPw1ySVNuvgIpEIbW1tdHR00NbWRjQa9fatrlq1inQ6TSaT4TvfOZfdu9/Geec9xdVXP874eCtjY2PegmBz3l08Hicej9PR0UEsFqO+vn7aAZ1wfIuuFhcMLxYX4ivwQudCJQSRUnHYsuXNXHLJ81x00YPs3j1OMplkaGiIbDZLa2sr2Wx22to3M55nFgGb1IVmq1goFCKRSJDL5bjjjtP4+c97uPzyF/jgB59leLjVm8w4cuQIIyMjnk2TLtEk0TFr9UrHJqXrujhcia9AC50rlRBESrdqXXLJ89x44xb27x9lcHCQ/v5+BgcHyWazXostFotN2w1hMMcsmZ0QZklJoVDgb//2jWzadBJXXvkiH/rQk6TTWS8vxNDQEP39/QwPD3uZxszkQ2NjI+Pj4163t7TLPXMrmzA3TsWX2SjtQqEbXS5bXtqiO7/Sqbe8tKXszyyGheyzsXxfqw4UyzIoFArTXh88eNAz61JJJpPz/l1Bp9z7ttrxZQB26DK0I5D7XZx60gSURCJRbReWhOyMWBgX4ytwXVcXKyEIzBSElpYWMplizlRzpFKhUPDGxcy6O5MforT7qPXr+VjNbO5sv5/tX1NKu7kz/Zz5vcb34u9kZ8R8uBpfgRI6VyshiJQejikEA5fjKzBdV5crQRBsx/X4CkSLzvVKCCozu42lXUvDQpm2Sn83M3fsfN+50Ptm+14Zj5udIMSX80IXhEoIKjOFqfQI9NKflWOnnPfrksW95t/5vkuEbWGCEl9OC11QKiHIzCYmlRKYlfyuWiBI8eXsGF2QKkEQbCNo8eWk0AWtEgTBJoIYX84JXRArQRBsIojx5ZzQBbESBMEmghhfzgldECtBEGwiiPHlnNAFsRIEIcjYEF/OCZ3f2FAJghBUbImvmhY6WypBEIKITdn2alboROQEoXLYlte1JoXOpieNIASNlVwCVi41J3S2PWkEIUjYus61KkKnlLpJKfWcUqqglFq/Ut9r45NGEIKCrSIH1WvR7QJuAH62Ul9ocyUIguvYHl9VOb1Ea70bVu5kCdsrQRBcxoX4CvwYnQuVIAiu4kp8VaxFp5R6GDhhll99QWv9H4uwcwtwCwAti/PBlUoQBBdxKr7KyYlYqQJsBdaX/X7J6yoIFUXyujqEU08aQXAMJ+OrHDX0uwDvBV4BMsBrwI/L+lwZLTprnjTSohMcZKH71pb4MlBmi66qXdfFloWEzqZKEKETXGS++9am+DKUK3SB6bo62ZwWBEdwPb4CIXSuV4Ig2EwQ4st5oQtCJQiCrQQlvpwWuqBUgiDYSJDiy1mhC1IlCIJtBC2+nBS6oFWCINhEEOPLOaELYiUIgk0EMb6cE7ogVoIg2EQQ48s5oQtiJQiCTQQxvpwTuiBWgiAEGRviyzmh8xsbKkEQgoot8VXTQmdLJQhCELEp217NCp2InCBUDtuy7dWk0Nn0pBGEoGFjtr2aEzrbnjSCECRsXedaU0Jn45NGEIKCrSIHoIpn17mBOlFpfr/aXgiCYA0b2am1Xr/g+8o5ndOWspjkOKXcse0OrTYqfce2O5b0+aXY9/OEYRtPdhX7wbS/lPu2GvFloBaPUl/sRfKDuez7JXQ2BYHYD779xd631RQ5rUXoyrpIy2XeJ40PQmdbEIj94NtfzH1bbZHTWoSu6pWwXKGzMQjEfvDtl3vfVju+DDUtdDZUwnKEztYgEPvBt1/OfWtDfBlqVuhsqYSlCp3NQSD2g29/ofvWlvgy1KTQ2VQJSxE624NA7Aff/nz3rU3xZag5obOtEhYrdC4EgdgPvv257lvb4stQU0JnYyUsRuhcCQKxH3z7s923NsaXoWaEztZKKFfoXAoCsR98+zPvW1vjy1ATQmdzJZQjdK4FgdgPvv3S+9bm+DIEXuhsr4SFhM7FIBD7wbdv7lvb48sQaKFzoRLmEzpXg0DsB98+G3EivgyBFTpXKmEuoXM5CMR+8O2zESfiyxBYoXOlEmYTOteDQOwH375p0VWCSjRSAit0rlTCTKELQhCI/eDb9/N4sVIq1RMLrNBVgoo8aUpumKAEgdgPvv1KCF0lh5tE6MqkYk+aqRsmSEEg9oNv32+hq/SYutVCB3wV+DXwK+B+oLWsz/ksdBV90mwkcEEg9oNv30+hW4mJQ9uF7kogNPX/24Hby/qcj0JX8SfNRgIXBGI/+Pb9ErqVWh1htdBNcwDeC9xV1nt9EroVedJMtegqQVCDTOxX374fQreSS8BcErr/BD40z+9vAXYAO2hxpxIqNXsV5CAT+9W3v9z7dqXXuVZd6ICHgV2zlHeXvOcLU2N0qiyby2zRreiTpgJCV+0gEPvBt7+c+7Yai/mrLnQLfjF8FNgGxMr+zDKEbsWfND4LnQ1BIPaDb3+p9221dixZLXTA1cDzQNeiPid5XX2zKfbF/mxIXld/he5F4BDwzFT5h7I+J3ldfUfsi/1SJK+rBUXyuvqL2Bf7M5G8rhYUyevqH2Jf7M+G5HW1oEheV38Q+2J/LiSvqwVF8rouH7Ev9udD8rpaUCSv6/IQ+2J/ISSvqwVF8rouHbEv9stB8rpaUCSv69IQ+2K/XCSvqwVF8rouHrEv9heD5HW1oEhe18Uh9sX+YpG8rhYUyetaPmJf7C8FyetqQZG8ruUh9sX+UpG8rhYUyeu6MGJf7C8HyetqQZG8rvMj9sX+cu1LXlcLiuR1nRuxL/b9sC95XS0oktd1dsS+2PfLvuR1taBIXtfjEfti30/7ktfVgiJ5Xacj9sW+3/aDmtdVaa1xBaXUKLCn2n6USScwUG0nFoFL/rrkK7jlr0u+ArxJa51Y6E2hlfDER/ZorddX24lyUErtcMVXcMtfl3wFt/x1yVco+lvO++oq7YggCEK1EaETBCHwuCZ036q2A4vAJV/BLX9d8hXc8tclX6FMf52ajBAEQVgKrrXoBEEQFo0InSAIgcc5oVNK/ZVS6ldKqWeUUj9RSp1YbZ/mQin1VaXUr6f8vV8p1Vptn+ZDKXWTUuo5pVRBKWXlEgOl1NVKqT1KqReVUn9SbX/mQyl1p1KqXym1q9q+LIRS6iSl1Bal1PNT98Cnq+3TXCilGpVSTyqlfjnl618u+BnXxuiUUs1a65Gp/98GvEVr/ckquzUrSqkrgUe11jml1O0AWus/rrJbc6KUejNQAP4v8Eda67LWKK0USql64AXgCuAVoBd4v9b6+ao6NgdKqYuAMeBftdZrq+3PfCiluoFurfVTSqkEsBN4j43XVimlgLjWekwpFQYeAz6ttd4+12eca9EZkZsiDlir1Frrn2itc1MvtwNrqunPQmitd2utbd55cg7wotZ6v9Z6ErgbeHeVfZoTrfXPgGPV9qMctNaHtdZPTf1/FNgNvKG6Xs3O1A6wsamX4akyrw44J3QASqkvKaUOAR8E/qLa/pTJ7wEPVdsJx3kDcKjk9StYGowuo5TqAc4CnqiyK3OilKpXSj0D9AObtdbz+mql0CmlHlZK7ZqlvBtAa/0FrfVJwF3ArTb7OvWeLwA5iv5WlXL8FWoXpVQTcC/wP2b0nqxCa53XWq+j2Es6Ryk179CAlXtdtdaXl/nWu4BNwBcr6M68LOSrUuqjwLXAZdqCAdFFXFsbeRU4qeT1mqmfCT4wNd51L3CX1vq+avtTDlrrIaXUFuBqYM5JHytbdPOhlDq95OW7gV9Xy5eFUEpdDXweuF5rPVFtfwJAL3C6UupUpVQEuBl4oMo+BYKpAf5vA7u11l+vtj/zoZTqMisYlFJRipNT8+qAi7Ou9wJvojg7eBD4pNbayqe6UupFoAFITv1ou60zxABKqfcC3wS6gCHgGa31VVV1agZKqWuAvwHqgTu11l+qrkdzo5T6HnAJxaOPXgO+qLX+dlWdmgOl1AXAz4FnKcYWwJ9prTdVz6vZUUqdCfwLxXugDvh3rfX/nPczrgmdIAjCYnGu6yoIgrBYROgEQQg8InSCIAQeETpBEAKPCJ0gCIFHhE5wgqnTNV5SSrVPvW6bet1TZdcEBxChE5xAa30I+Hvgy1M/+jLwLa31gao5JTiDrKMTnGFqi9JO4E7gE8A6rXW2ul4JLmDlXldBmA2tdVYp9Tngv4ArReSEcpGuq+Aa7wQOA1YfZCnYhQid4AxKqXUUN3CfB/zh1Km4grAgInSCE0ydrvH3FM9Jexn4KvC16noluIIIneAKnwBe1lpvnnr9f4A3K6UurqJPgiPIrKsgCIFHWnSCIAQeETpBEAKPCJ0gCIFHhE4QhMAjQicIQuARoRMEIfCI0AmCEHj+P6PNGScq9M6XAAAAAElFTkSuQmCC\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 0x7f96eff75ee0>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPUAAAD4CAYAAAA0L6C7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAMHklEQVR4nO3dX2zddRnH8c+npxtbOxibA5F1ssU/mElUsBqEqASIUSFwoTGYYCJe7ELRQUgMemPiNSESNSYLwo2LXEwuDCGiCRpF42LZprANIk5kGxPGv/1pYVvbx4vWZG6059fT79df++T9Ski2nvLwpJx3f6enp986IgQgj762FwBQFlEDyRA1kAxRA8kQNZBMf42hncHB6F+9uvhcTxYfOTV3os7cGlzpmxWTVe4JkpfX+eB2+irdGSpYtfTN4jNfO/imRl8/6be7rcr/yv7VqzW0+c7icztvFR8pSVpy9G0/NvPSV+kTRd+pOnPfWlNn7pKPvF5l7nnLThSfWeubu18Y2lV85n1f+vOMt/HwG0iGqIFkiBpIhqiBZIgaSIaogWQaRW37s7aftf2c7btrLwWgd12jtt2R9GNJn5O0UdKXbW+svRiA3jS5Un9c0nMRsS8iTkp6SNLNddcC0KsmUa+VtP+0vx+Yftv/sL3J9ojtkcnjo6X2AzBHxZ4oi4gtETEcEcN9KwZLjQUwR02iPihp3Wl/H5p+G4AFqEnUf5H0PtsbbC+VdIukX9ZdC0Cvuv6UVkSM275d0mOSOpIeiIjd1TcD0JNGP3oZEY9KerTyLgAK4BVlQDJEDSRD1EAyRA0kQ9RAMpXOkJSiU/4Yt+grf0CgJA28XH7XgZfqnBA4sazO5+GxizpV5n79/b+vMvdPR95TfuYTHyw+U5J+9NI1xWe+PDrzd5W5UgPJEDWQDFEDyRA1kAxRA8kQNZAMUQPJEDWQDFEDyRA1kAxRA8kQNZAMUQPJEDWQDFEDyRA1kAxRA8kQNZAMUQPJEDWQDFEDydQ5TdTSZIXJMVD+1E9JOnZJ+c9t4SXFZ0rS8tcmqsxVnQ+tzu+MVZn7znOOFZ/ZOVF8pCRp8liF+8LkzCfrcqUGkiFqIBmiBpIhaiAZogaSIWogGaIGkukate11tn9re4/t3bY3/z8WA9CbJi8RGZd0V0TssH2upCdt/yYi9lTeDUAPul6pI+JQROyY/vMxSXslra29GIDezOlratvrJV0uafvb3LbJ9ojtkYnjo4XWAzBXjaO2vULSLyTdERFHz7w9IrZExHBEDHdWDJbcEcAcNIra9hJNBb01Ih6uuxKA+Wjy7Lcl/VTS3oi4t/5KAOajyZX6aklfkXSt7V3T/3y+8l4AetT1W1oR8YSkmX94E8CCwivKgGSIGkiGqIFkiBpIpsrBg/1j0pqd5Z9bGzx0qvhMSTp6SfmD4V756GTxmZIUK+rMHVhZ51WA+05cWGXuP46vKT5zcmnxkVPOq3C/7cx8UiRXaiAZogaSIWogGaIGkiFqIBmiBpIhaiAZogaSIWogGaIGkiFqIBmiBpIhaiAZogaSIWogGaIGkiFqIBmiBpIhaiAZogaSIWogmSqniXbenNDqv571227nLXbuLj5Tks694WPFZ37oa88UnylJD777D1Xm3vr8NVXmPvDUVVXmTr5a/ujPZUfq/Hap8VWd8kNnOVSWKzWQDFEDyRA1kAxRA8kQNZAMUQPJEDWQTOOobXds77T9SM2FAMzPXK7UmyXtrbUIgDIaRW17SNINku6vuw6A+Wp6pf6BpG9rlhen2d5ke8T2yKnxsRK7AehB16ht3yjp5Yh4crb3i4gtETEcEcNL+geKLQhgbppcqa+WdJPt5yU9JOla2z+ruhWAnnWNOiK+ExFDEbFe0i2SHo+IW6tvBqAnfJ8aSGZOP08dEb+T9LsqmwAogis1kAxRA8kQNZAMUQPJEDWQTJXTRE+d29Gha1YWnzt5fZ2TKceXV5j5xgXlh0rauuodVeb+8W/vrzL3nU/UuW4cHyo/d+ziieIzJalvtMZpojOffMqVGkiGqIFkiBpIhqiBZIgaSIaogWSIGkiGqIFkiBpIhqiBZIgaSIaogWSIGkiGqIFkiBpIhqiBZIgaSIaogWSIGkiGqIFkiBpIpsppoloxofj068XH3vbe7cVnStIP/3xt8Zmjj1xUfKYkff+iL1WZu+KNmU+nnI836xyqqvjYkeIzr7jwpeIzJenpx8uf1OpZDj7lSg0kQ9RAMkQNJEPUQDJEDSRD1EAyRA0k0yhq2+fb3mb7Gdt7bX+i9mIAetP0xSf3SfpVRHzR9lJJAxV3AjAPXaO2vVLSpyR9VZIi4qSkk3XXAtCrJg+/N0g6LOlB2ztt32978Mx3sr3J9ojtkYmjY8UXBdBMk6j7JV0h6ScRcbmkUUl3n/lOEbElIoYjYrhzHo/OgbY0ifqApAMR8d+fptimqcgBLEBdo46If0vab/vS6TddJ2lP1a0A9Kzps9/flLR1+pnvfZJuq7cSgPloFHVE7JI0XHcVACXwijIgGaIGkiFqIBmiBpIhaiCZKqeJ9jm0bMl48bn7Kh1NObBvafGZ73ro2eIzJenUB9ZVmbv/M8urzD25/kSVue8+d7T4zD0v1TkB9vy/TxafefCtmW/jSg0kQ9RAMkQNJEPUQDJEDSRD1EAyRA0kQ9RAMkQNJEPUQDJEDSRD1EAyRA0kQ9RAMkQNJEPUQDJEDSRD1EAyRA0kQ9RAMlUOHpw83q/RP5Y/JPCxlWuKz5SkcybKzzz2yfeWHypp7MI6n4dPrYgqc/1q+UMdJenFF99VfOaSoy4+U5L6xssfPDjbplypgWSIGkiGqIFkiBpIhqiBZIgaSIaogWQaRW37Ttu7bT9t++e2l9VeDEBvukZte62kb0kajojLJHUk3VJ7MQC9afrwu1/Sctv9kgYkvVhvJQDz0TXqiDgo6R5JL0g6JOlIRPz6zPezvcn2iO2R8bHyvzsYQDNNHn6vknSzpA2SLpY0aPvWM98vIrZExHBEDPcPDJbfFEAjTR5+Xy/pnxFxOCJOSXpY0lV11wLQqyZRvyDpStsDti3pOkl7664FoFdNvqbeLmmbpB2Snpr+d7ZU3gtAjxr9PHVEfE/S9yrvAqAAXlEGJEPUQDJEDSRD1EAyRA0kU+U00c4JadWz5Y/oHLugzuegNy4rv+vopRWOKJWkk3VO/ewcr/OxPefVOnOjU37miXeUP/VTkl45r/wppeO/n/k2rtRAMkQNJEPUQDJEDSRD1EAyRA0kQ9RAMkQNJEPUQDJEDSRD1EAyRA0kQ9RAMkQNJEPUQDJEDSRD1EAyRA0kQ9RAMkQNJEPUQDKOKH86pe3Dkv7V4F3XSHql+AL1LKZ9F9Ou0uLadyHseklEXPB2N1SJuinbIxEx3NoCc7SY9l1Mu0qLa9+FvisPv4FkiBpIpu2oF9svr19M+y6mXaXFte+C3rXVr6kBlNf2lRpAYUQNJNNa1LY/a/tZ28/ZvrutPbqxvc72b23vsb3b9ua2d2rCdsf2TtuPtL3LbGyfb3ub7Wds77X9ibZ3mo3tO6fvB0/b/rntZW3vdKZWorbdkfRjSZ+TtFHSl21vbGOXBsYl3RURGyVdKekbC3jX022WtLftJRq4T9KvIuIDkj6sBbyz7bWSviVpOCIuk9SRdEu7W52trSv1xyU9FxH7IuKkpIck3dzSLrOKiEMRsWP6z8c0dadb2+5Ws7M9JOkGSfe3vctsbK+U9ClJP5WkiDgZEW+0ulR3/ZKW2+6XNCDpxZb3OUtbUa+VtP+0vx/QAg9Fkmyvl3S5pO0tr9LNDyR9W1Kd36JezgZJhyU9OP2lwv22B9teaiYRcVDSPZJekHRI0pGI+HW7W52NJ8oasr1C0i8k3RERR9veZya2b5T0ckQ82fYuDfRLukLSTyLickmjkhby8yurNPWIcoOkiyUN2r613a3O1lbUByWtO+3vQ9NvW5BsL9FU0Fsj4uG29+niakk32X5eU1/WXGv7Z+2uNKMDkg5ExH8f+WzTVOQL1fWS/hkRhyPilKSHJV3V8k5naSvqv0h6n+0Ntpdq6smGX7a0y6xsW1Nf8+2NiHvb3qebiPhORAxFxHpNfVwfj4gFdzWRpIj4t6T9ti+dftN1kva0uFI3L0i60vbA9P3iOi3AJ/b62/iPRsS47dslPaapZxAfiIjdbezSwNWSviLpKdu7pt/23Yh4tL2VUvmmpK3Tn9z3Sbqt5X1mFBHbbW+TtENT3xXZqQX4klFeJgokwxNlQDJEDSRD1EAyRA0kQ9RAMkQNJEPUQDL/AUiMqsQgz5/lAAAAAElFTkSuQmCC\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": "iVBORw0KGgoAAAANSUhEUgAAA0kAAAGeCAYAAABB3Ur+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABbLElEQVR4nO3deXhU5dnH8e9NAImRRRBlUwFRtqAEI1aqFtywAoKKiuKCCNbaWuyCFa0WrW+tpbXS1y6CCy6IKCJKsWKtxqVuKCCLglhEQxYQYoCQREK43z9mkndCFhJI5sxkfp/rmouZ58ycc58z0ZNfnuc8x9wdERERERERCWkSdAEiIiIiIiKxRCFJREREREQkgkKSiIiIiIhIBIUkERERERGRCApJIiIiIiIiERSSREREREREIigkiYhEMLOuZuZm1nQ/Pz/WzF6p77oaCzPLMLMJQdchIiJSE4UkEYlpZrbBzIrMrMDMcs1slpkdEnRdUHWgcvfZ7n5OFLb9HTP7l5nlmdnXZvasmXWMWG5mdq+ZbQ0/7jUzq2F9t5rZF+HjvNHM5jb0PuyLmV1tZh+Z2fZwTb+vKrya2bFmVmxmT+7VfrmZfWlmO81sgZm1rWFbbmYrzaxJRNvdZjarnvZlhpmtNbM9ZjauiuU/Df98bzezR8zsoGrWU/Yz99Je7U+a2dRYqVNEJN4pJIlIPBjh7ocA/YE0YEqw5cSEQ4EZQFfgaGAH8GjE8uuAUcAJwPHACOAHVa3IzK4GrgTOCh/ndODfDVR3XRwM3AQcBpwMnAn8oor3/QVYEtlgZn2BBwnt1xFAIfDXfWyvEzDmgCqu3sfADcDSvReY2VDgFkL7dzTQHbhzH+s72cwG1XeR1H+dIiJxSSFJROKGu+cCiwmFJaC8R+UdM8s3s4/NbHDEsnFmtt7MdoR7ScaG25uY2a/CvQybzexxM2td1TbDPVlnRbyeGtFj8Wb43/xwD8wp4W2+HfH+QWa2xMy2hf8dFLEsw8x+Y2b/Cdf4ipkdVstj8U93f9bdt7t7IfAA8N2It1wN/NHdN7p7FvBHYFw1qzsJWOzu/w2vO9fdZ0TU2drMHjazHDPLCvewJEUsH29mn5rZN2a22MyOjlh2tpmtCe//A0C1vVlV7OPf3P0td98V3ofZe+0jZjYGyKdyqBsLLHT3N929ALgduNDMWtawyd8Dd1bTW3Vp+GeoVfj198M9Ku1ruS9/cfd/A8VVLL4aeNjdV7v7N8BvqP67iqz1f6paYGarzGxExOtmZrbFzNICqFNEJC4pJIlI3DCzLsD3gc/DrzsDi4C7gbaEehmeM7P2ZpYC/Bn4vru3BAYBy8OrGhd+DCH01/BDCIWMujo9/G8bdz/E3d/dq9624fr+DLQD7gMWmVm7iLddDlwDHA40J6KnxMxWmNnldahldcTrvoR6Bcp8HG6rynvAVWY22czSIwNQ2CxgN9CDUE/eOcCEcI0jgVuBC4H2wFvAnPCyw4D5wK8I9Qb9l4iQY2ZHhcPtUfuzj+HAchfwsyreW2H/wwFwF3BcDeufD2ynil/83X0u8A7w5/D39zAwwd2/DtfyDzO7pZb7UWOt4edH7PVzsre/AsdFBvgIjwNXRLw+D8hx92UB1CkiEpcUkkQkHiwwsx1AJrAZ+HW4/QrgJXd/yd33uPu/gA8J/VIIsAdINbNkd89x97JfsMcC97n7+nAvwxRgTFU9CAdoGLDO3Z9w993uPgdYQ2joW5lH3f0zdy8CniGil8zdj3f3p/a1ETM7HrgDmBzRfAiwLeL1NuAQs8rXJbn7k8CNwFDgDWCzmf0yvO4jCB3Pm9x9p7tvBv7E/w9Lux64x90/dffdwG+B/uHepPOA1e4+z91LgPuB3IjtfuXubdz9q1rs43hCwwD/ENH8G0I9Gxur+Mje+192DGrqSXJCPU63m1nzKpb/CDgDyCDUS/WPiH0Z7u6/29d+VKOq74p91FpEqCfp7iqWPQmcV9brRWjI4RMB1SkiEpcUkkQkHowK9wYNBnoR6pWA0HURF4d7I/LNLB84Fejo7juBSwn9Ep9jZovMrFf4c52ALyPW/yXQlNC1K/Vp7+2UbatzxOvciOeFhH4RrTUz6wH8E5jk7m9FLCoAWkW8bgUUuLtXtZ7whBNnAW0IHbPfhK9BORpoRugYlh3jBwn1fBFePj1iWR6hIXWdCe1/ZsQ2PPJ1HfZxFHAPoV7BLeG2/sBZhAJbVfbef8Kvd9S0LXd/CdhIFddvuXs+8CyQSmj4Yn2p6ruCfdQKPESoJycydOPu2cB/gIvMrA2h3tfZAdYpIhJ3FJJEJG64+xuEhn6V9SZkAk+EeyPKHillfyl398XufjbQkVAPzszw57IJ/XJf5ihCw8k2VbHZnYQmECjTIbKkfZS893bKtpW1j8/VSri35lXgN+7+xF6LVxOatKHMCVQcjlcldy9x92eBFYTCQCbwLXBYxDFu5e5lQ/cygR/s9R0ku/s7QA5wZES9Fvm6lvt4LqHvbYS7r4xYNJjQpBVfmVkuoWGKF5lZ2YQDFfbfzLoDBwGf1WKztxEaQhj5vZcFs/GEhhP+uS77sQ9VfVeb3H1rTR9y912EJk74DZWv9XqMUE/rxcC74Wu6AqlTRCQeKSSJSLy5HzjbzE4gNKxohJkNNbMkM2thZoPNrIuZHWFmI8PXJn1L6K/ge8LrmAP81My6WWg68d8Cc8PDxfa2nNBQvGZmlg6Mjlj2dXid3aup9SVC141cbmZNzexSoA/wj2reX2vh67FeAx5w979X8ZbHgZ+ZWWcz6wT8nFDArGpd48xsmJm1tNCkFt8ndP3J++6eA7wC/NHMWoWXH2Nm3wt//O/AFAvNJlc2ycPF4WWLgL5mdmF4KONPqBgy97WPZxDqAbnI3T/Ya/EM4BhCwxP7h+tYRGjIIOHPjTCz08I/A3cB8919n70e7p4BrCI0UUFZLS0I/bzdSugass5mdkMd9qV5eB0GNAv/rJadgx8HrjWzPuGen19RzXdVhSeAFsC5e7UvAAYAk8LrD7pOEZG4opAkInElfKH848Ad7p4JlE0c8DWhXo3JhP7f1oTQBf3ZhIaAfQ/4YXg1jxD65fJN4AtCM3ndWM0mbyf0y/g3hP5qX36NUHhWuf8B/hMebvadvWrdCgwnFFC2AjcDw8uGjO2Lma228Ix8VZhAKJxNtdDMegVmVhCx/EFgIbCS0C/8i8JtVdlO6Bh+RWimuN8DP3T3sln6riI0qcQnhI7DPEK9c7j788C9wNNmtj28re+Hl20h1JPxu/D+H0toGFjZ/h0Vrru6iRtuB1oDL0Xs4z/D6y4Mz8KXG571sAAoLptIIXz92fWEwtJmQtfN1DrUEAoAkfdVugfIDM+49y2hXpq7zezY8L7808xurWF9rxC6jmgQoYBXRHjiD3d/mdAxf53Qd/Al/3/dXY3cvZTQ9Wht92ovAp4DuhGakIIg6xQRiTdWzfB0ERERiWNmdgdwnLtfsc83i4hIBfU9k5OIiIgEzELTz19LaGY7ERGpIw23ExERaUTMbCKhoaf/dPc39/V+ERGpTMPtREREREREIqgnSUREREREJIJCkoiIiIiISASFJBERERERkQgKSSIiIiIiIhEUkkRERERERCIoJImIiIiIiERQSBIREREREYmgkCQiIiIiIhJBIUlERERERCSCQpKIiIiIiEiEpkEXUBeHHXaYd+3atd7Xu3PnTlJSUup9vfEikfdf+56Y+w6Jvf/V7ftHH320xd3bB1BSzGuo80+kRP6ZrI6OSWU6JpXpmFSmY1JZVcekpvNeoCHJzNoADwGpgAPj3f3d6t7ftWtXPvzww3qvIyMjg8GDB9f7euNFIu+/9n1w0GUEJpH3v7p9N7Mvo19NfGio80+kRP6ZrI6OSWU6JpXpmFSmY1JZVcekpvNe0D1J04GX3X20mTUHDg64HhERERERSXCBhSQzaw2cDowDcPddwK6g6hEREREREYFgJ27oBnwNPGpmy8zsITPT4EkREREREQlUkMPtmgIDgBvd/X0zmw7cAtwe+SYzuw64DuCII44gIyODvZaTkpJCUlLSfhfSqlUrli1btt+fj3eNZf9LS0vZuXMn7l7rzxQUFFT6mUoUibzvkNj7n8j7LiIiUhtBhqSNwEZ3fz/8eh6hkFSBu88AZgCkp6f73hdcffHFF7Rs2ZJ27dphZvtVyI4dO2jZsuV+fbYxaAz77+5s3bqVHTt20K1bt1p/LpEvbEzkfYfE3v9E3ve6MrMRwIgePXoEXYqIiERRYMPt3D0XyDSznuGmM4FP6rqe4uLiAwpI0jiYGe3ataO4uDjoUkSkEXH3he5+XevWrYMuRUREoijo2e1uBGaHZ7ZbD1yzPytRQBLQz4GIiIiI1I8gJ27A3Ze7e7q7H+/uo9z9myDrkdrJzs5m9OjRQZchIiIiItIgAg1JEn92795Np06dmDdvXtCliIiIiIg0CIWkA3THHXdw//33l7++7bbbmD59+j4/t23bNnr27MnatWsBuOyyy5g5c2al9y1ZsoRBgwZxwgknMHDgQHbs2EFxcTHXXHMN/fr1Iy0tjddffx2AWbNmMWrUKM4++2y6du3KAw88wH333UdaWhrf+c53yMvLA2Dw4MFMmjSJ/v37k5qaWn4X+Q8++IBTTjmFtLQ0Bg0aVF7brFmzOP/88znjjDM488wz2bBhA6mpqQCsXr2agQMH0r9/f44//njWrVsHwH333Udqaiqpqanlx2fDhg307t2biRMn0rdvX8455xyKior246iLiIiIiDScoK9Jqld3LlzNJ9nb6/y50tLSaqcQ79OpFb8e0bfaz44fP54LL7yQm266iT179vD000/zwQcfsGPHDk477bQqP/PUU0/Rp08fHnjgAcaNG8ekSZP45ptvmDhxYoX37dq1i0svvZS5c+dy0kknsX37dpKTk5k+fTpmxsqVK1mzZg3nnHMOn332GQCrVq1i2bJlFBcX06NHD+69916WLVvGT3/6Ux5//HFuuukmAAoLC1m+fDlvvvkm119/PZ988gm9evXirbfeomnTprz66qvceuutPPfccwAsXbqUFStW0LZtWzZs2FBe49///ncmTZrE2LFj2bVrF6WlpXz00Uc8+uijvP/++7g7J598Mt/73vc49NBDWbduHXPmzGHmzJlccsklPPfcc1xxxRW1/apERERERBpcowpJQejatSvt2rVj2bJlbNq0ibS0NNq1awfA8uXLa/zs2WefzbPPPsuPfvQjPv7440rL165dS8eOHTnppJOA0P2MAN5++21uvPFGAHr16sXRRx9dHpKGDBlCy5YtadmyJa1bt2bEiBEA9OvXjxUrVpSv+7LLLgPg9NNPZ8eOHeTn57Njxw6uvvpq1q1bh5lRUlJSoda2bdtWqvGUU07hf/7nf9i4cSMXXnghxx57LG+//TYXXHABKSmhewNfeOGFvPXWW5x//vl069aN/v37A3DiiSdWCFwisn8WrV/E9KXTyd2ZS4eUDkwaMIlh3YcFXZZIVC1YlsW0xWvJzi+iU5tkJg/tyai0zkGXJSJxqlGFpJp6fGpyoPcJmjBhArNmzSI3N5fx48eXr3NfPUl79uzh008/5eCDD+abb76hS5cu+11DmYMOOqj8eZMmTcpfN2nShN27d5cv23smODPj9ttvZ8iQITz//PNs2LChwn1UygLP3i6//HJOPvlkFi1axHnnnceDDz5Y6/qSkpI03E7kAC1av4ip70yluDQ0/X3OzhymvjMVQEFJEsaCZVlMmb+SopJSALLyi5gyfyWAgpKI7Bddk1QPLrjgAl5++WWWLFnC0KFDAWjZsiXLly+v8tGnTx8A/vSnP9G7d2+eeuoprrnmmgo9NwA9e/YkJyeHJUuWAKHgtXv3bk477TRmz54NwGeffcZXX31Fz549qYu5c+cCoV6pVq1a0bp1a7Zt20bnzqGTyaxZs2q1nvXr19O9e3d+8pOfMHLkSFasWMFpp53GggULKCwsZOfOnTz//PPVBkYROTDTl04vD0hlikuLmb5039dGijQW0xavLQ9IZYpKSpm2eG1AFYlIvGtUPUlBad68OUOGDKFNmzbVXtu0t7Vr1/LQQw/xwQcf0LJlS04//XTuvvtu7rzzzgrrnTt3LjfeeCNFRUUkJyfz6quvcsMNN/DDH/6Qfv360bRpU2bNmlWhh6Y2WrRoQVpaGiUlJTzwwAMA3HzzzVx99dXcfffdDBtWu79AP/PMMzzxxBM0a9aMDh06cOutt9K2bVvGjRvHwIEDgVBPW1pamobWiTSA3J051bTnRrkSkeBk51c9KqG6dhGRfVFIqgd79uzhvffe49lnn631Z3r27Mmnn35a/vq+++6r8n0nnXQS7733XqX2Rx99tFLbuHHjGDduXPnryFCy97IrrriifNa5HTt2AKHri8qubQK4++67q/xs165dWbVqFQC33HILt9xyS6Vafvazn/Gzn/2sQlvk5wB+8YtfVPqciNTe10vm06FkNznNKv+vvENKhwAqEglGpzbJZFURiDq1SQ6gGhFpDDTc7gB98skn9OjRgzPPPJNjjz026HJEJEF8vWQ+bRZN5IJvDqJ5k4o9yS2SWjBpwKSAKhOJvslDe5LcrOJIjuRmSUweWreh6CIiZdSTdID69OnD+vXrgy6jTjIyMoIuQUQOQGZeISv/+ThH0pXvXfQ8R7FSs9tJQiubnEGz24lIfVFIEhGJI5lbtjPmoQ8p2vMDnry6H326daEPRykUScIbldZZoUhE6o2G24mIxImvP1zArr98l+TizTw+cRB9uh34bQNERESkMoUkEZE48PWHC2jzj2sp8mb8+apBpHZuHXRJIiIijZZCkohIjCsLSGs5miZXLaBP96OCLklERKRRU0iqB0lJSfTv358TTjiBAQMG8M477+zXeu6//34KCwurXPbWW2/Rt29f+vfvT1ZWFqNHjwZg+fLlvPTSS/tdu4jEts3LX1ZAEhERiTKFpHqQnJzM8uXL+fjjj7nnnnuYMmXKfq2nppA0e/ZspkyZwvLly+ncuTPz5s0DFJJEGrPMvEKu+WcR/+S7CkgiIiJRpJBUz7Zv386hhx5a/nratGmcdNJJHH/88fz6178GYOfOnQwbNowTTjiB1NRU5s6dy5///Geys7MZMmQIQ4YMqbDOhx56iGeeeYbbb7+dsWPHsmHDBlJTU9m1axd33HEHc+fOpX///sydOzeq+yoiDWfTJ29zxYNvsXFXCt0nPqGAJCIiEkWNbwrwR6uYBrfvKBg4EXYVwuyLKy1u2utCOOVa2LkVnrmq4sJrFu1zk0VFRfTv35/i4mJycnJ47bXXAHjllVdYt24dH3zwAe7O+eefz5tvvsnXX39Np06dWLQotO5t27bRunVr7rvvPl5//XUOO+ywCuufMGECb7/9NsOHD2f06NFs2LABgObNm3PXXXfx4Ycf8sADD+z72IhIXPj6oxc4dOF4rmI4J0/8syZpEBERiTL1JNWDsuF2a9as4eWXX+aqq67C3XnllVd45ZVXSEtLY8CAAaxZs4Z169bRr18//vWvf/HLX/6St956i9at9QuQiIR8/dELtF44ns84mlOuuksBSUREJACNryeppp6f5gdXuXz3jh2hJyntatVzVJNTTjmFLVu28PXXX+PuTJkyhR/84AeV3rd06VJeeuklfvWrX3HmmWdyxx13HNB2RST+lQWkdRxFk6uep0/3o4MuSUREJCGpJ6merVmzhtLSUtq1a8fQoUN55JFHKCgoACArK4vNmzeTnZ3NwQcfzBVXXMHkyZNZunQpAC1btmRHWWCrpf35jIjEno2btpC08EbWcRR21QIFJBERkQA1vp6kAJRdkwTg7jz22GMkJSVxzjnn8Omnn3LKKacAcMghh/Dkk0/y+eefM3nyZJo0aUKzZs3429/+BsB1113HueeeS6dOnXj99ddrte0hQ4bwu9/9jv79+zNlyhQuvfTSBtlHEWk4mXmFjHl0BV24jV9fda4CUgwxsxHAiB49egRdioiIRJFCUj0oLS2tdtmkSZOYNGlShbZjjjmGoUOHVnrvjTfeyI033ljlembNmlX+vGvXrqxatQqAtm3bsmTJkv2oWkRiwddLF/LCS/+ioHQ4t08cQx9dgxRT3H0hsDA9PX1i0LWIiEj0aLidiEhAvl66kNYvjuOM3W/y1DX9NUmDiIhIjFBPkohIAMoC0ucciV35An2OOjzokkRERCRMPUkiIlG2d0DqfYyuQRIREYklCkkiIlGUmVfIwy/9JzSLnQKSiIhITFJIEhGJko05uYyZ8R5zSs/Ar/2XApKIiEiMUkgSEYmCr5f+g1YPnkiP4lXMnnAyqUceFnRJIiIiUg2FpHqyYMECzIw1a9ZU+57Bgwfz4YcfAnDeeeeRn59f4zoHDRq0z+3ef//9FBYW1qnWIGVnZzN69OigyxCJqq+X/oNWL44ji/bccuVIzWInIiIN6p3sEr77u9fodssivvu711iwLCvokuJOwoWkResXcc68czj+seM5Z945LFq/qF7WO2fOHE499VTmzJlTq/e/9NJLtGnTpsb3vPPOO/tcTzyFpN27d9OpUyfmzZsXdCkiUVMWkNbTBXQNkoiINLAFy7KYtWoXWflFOJCVX8SU+SsVlOoooULSovWLmPrOVHJ25uA4OTtzmPrOVBZ/tfiA1ltQUMDbb7/Nww8/zNNPP13eXlRUxJgxY+jduzcXXHABRUVF5cu6du3Kli1bALjvvvtITU0lNTWV+++/v/w9hxxyCAAZGRkMHjyY0aNH06tXL8aOHYu78+c//5ns7GyGDBnCkCFDKtW1ZMkSBg0axAknnMDAgQPZsWMHxcXFXHPNNfTr14+0tDRef/11IHSz2lGjRnH22WfTtWtXHnjgAe677z7S0tL4zne+Q15eHhDqDZs0aRL9+/cnNTWVDz74AIAPPviAU045hbS0NAYNGsTatWvL13v++edzxhlncOaZZ7JhwwZSU1MBWL16NQMHDqR///4cf/zxrFu3rtrjsWHDBnr37s3EiRPp27cv55xzToXjKRKLctcuUUASEZGomrZ4Lbv2VGwrKill2uK1wRQUpxIqJE1fOp3i0uIKbcWlxfx99d8PaL0vvPAC5557Lscddxzt2rXjo48+AuBvf/sbBx98MJ9++il33nlneXukjz76iEcffZT333+f9957j5kzZ7Js2bJK71u2bBn3338/n3zyCevXr+c///kPP/nJT+jUqROvv/56edgps2vXLi699FKmT5/Oxx9/zKuvvkpycjJ/+ctfMDNWrlzJnDlzuPrqqykuDh2TVatWMX/+fJYsWcJtt93GwQcfzLJlyzjllFN4/PHHy9ddWFjI8uXL+etf/8r48eMB6NWrF2+99RbLli3jrrvu4tZbby1//9KlS5k3bx5vvPFGhRr//ve/M2nSJJYvX86HH35Ily5dajwe69at40c/+hGrV6+mTZs2PPfcc/vzdYlERWZeIRfPz+chLlRAEhGRqMnOr/qPyNW1S9USKiTl7sytsn1z0eYDWu+cOXMYM2YMAGPGjCkfcvfmm29yxRVXAHD88cdz/PHHV/rs22+/zQUXXEBKSgqHHHIIF154IW+99Val9w0cOJAuXbrQpEkT+vfvz4YNG2qsae3atXTs2JGTTjoJgFatWtG0aVPefvvt8pp69erF0Ucfzeeffw7AkCFDaNmyJe3bt6d169aMGDECgH79+lXY3mWXXQbA6aefzvbt28nPz2fbtm1cfPHFpKam8tOf/pTVq1eXv//ss8+mbdu2lWo85ZRT+O1vf8u9997Ll19+SXJyco3Ho1u3bvTv3x+AE088cZ/HQCQomz9+hZ88+A+273K+N3GaApJIAluwLEvXhkhUdWqTXKd2qVpChaQOKR2qbD88ef/vdJ+Xl8drr73GhAkT6Nq1K9OmTeOZZ57B3fd7nVU56KCDyp8nJSWxe/fuel3/3tto0qRJ+esmTZpU2J6ZVficmXH77bczZMgQVq1axcKFC8t7pwBSUlKq3N7ll1/Oiy++SHJyMueddx6vvfZaretrqGMgcqA2L3uJ1s9fzo+/nRmaxU6TNIgkrAXLspgyf6WuDZGomjy0J833+g0/uVkSk4f2DKagOJVQIWnSgEm0SGpRoa1FUguu73v9fq9z3rx5XHnllXz55Zds2LCBzMxMunXrxltvvcXpp5/OU089BYSGsq1YsaLS50877TQWLFhAYWEhO3fu5Pnnn+e0006r9fZbtmzJjh07KrX37NmTnJwclixZAsCOHTvYvXs3p512GrNnzwbgs88+46uvvuLYY4+t0z7PnTsXCPWCtW7dmtatW7Nt2zY6d+4MhK5Dqo3169fTvXt3fvKTnzBy5EhWrFhxwMdDJEibl71E6xeu4gs60enKmQpIIglu2uK1FJWUVmjTtSHS0EaldWZcanM6t0nGgM5tkrnnwn6MSuscdGlxpWnQBUTTsO7DgNC1Sbk7c+mQ0oFJAyZxevvT93udc+bM4Ze//GWFtosuuog5c+Zw3333cc0119C7d2969+7NiSeeWOF9ZsaAAQMYN24cAwcOBGDChAmkpaXVevvXXXcd5557bvm1SWWaN2/O3LlzufHGGykqKiI5OZlXX32VG264gR/+8If069ePpk2bMmvWrAo9NLXRokUL0tLSKCkp4ZFHHgHg5ptv5uqrr+buu+9m2LBhtVrPM888wxNPPEGzZs3o0KEDt956K23btq3yeGhoncS6yIDkV75I72O6Bl2SiARM14ZIUAZ1asatlw8Ouoy4ZvU9LKwhpaene9l9hsp8+umn9O7d+4DWu2PHDlq2bHlA66iL0tJSDj/8cHJzc2nWrFnUtluduuz/4MGD+cMf/kB6enoDV7V/6vrzUDZzYCJK5H2H+t3/zK07yX9gCM28OC4CUnX7bmYfuXts/scdsKrOP/Ut0f+brEq8H5Pv/u41sqoIRJ3bJPOfW87Yr3XG+zFpCDomlemYVFbVManpvJdQw+1iRd++fZkwYUJMBCQROTCZeYWMmfk+N/gv4yIgiUj0TB7ak+RmSRXadG2ISHxIqOF2sWLNmjVBl7DfMjIygi5BJGZsXv5P1r04nW/33MisiWfRW9cgiUiEsmtApi1eS3Z+EZ3aJDN5aE9dGyISBxSSRET2w+bl/6T1gqvoTEeeuDJVAUlEqjQqrbNCkUgcahTD7eLpuippOPo5kGgpC0gb6Mge3ShWRESk0Yn7kNSiRQu2bt2qX5ATnLuzdetWWrRose83ixyAzctf3isgdQu6JBEREalncT/crkuXLmzcuJGvv/56v9dRXFyc0L9cN5b9b9GiBV26dAm6DGnEMvMKuf2ljdzAsbS8crYCkoiISCMV9yGpWbNmdOt2YL+oZGRk1OneRI1Nou+/SG1kf7GGMXOzKSg5ioMnvETvLm2CLklEREQaSNwPtxMRaWibP15Mu8dO4/vFLzF7wsmkKiCJiIg0anHfkyQi0pA2f7yY1s9fwZd04KIrbtAsdiIiIglAIUlEpBqRAan0ihfo3aN70CWJSIxasCxL90MSaUQUkkREqpC18SsOff4qBSQR2acFy7KYMn8lRSWlAGTlFzFl/koABSWROKVrkkRE9pKZV8glT37ObfxYAUlE9mna4rXlAalMUUkp0xavDagiETlQ6kkSEYmw+eNX+NM/VlCw+3iunfgTXYMkIvuUnV9Up3YRiX3qSRIRCdv88Su0en4s40ueYva1J5GqgCQitdCpTXKd2kUk9ikkiYjw/wEpkyNIumIeqV0ODbokEYkTk4f2JLlZUoW25GZJTB7aM6CKRORAabidiCSkJQVL+O2835K7M5f2zdrw46wv6McR7NY1SCJSR2WTM2h2O5HGQyFJRBLOovWLmJM3hxIvAWBzyTf8z2Gtue67P+O6HscEXJ2IxKNRaZ0VikQakcCH25lZkpktM7N/BF2LiCSG6UunlwekMt82MeZlPx1QRRKrzGyEmc3Ytm1b0KWIiEgUBR6SgEnAp0EXISKJI3dnbp3aJXG5+0J3v651a03iISKSSAINSWbWBRgGPBRkHSKSWNo3a1Nle4eUDtEtRERERGJS0D1J9wM3A3sCrkNEEsTmFa/y46wvOGiPV2hvkdSCSQMmBVSViIiIxJLAJm4ws+HAZnf/yMwG1/C+64DrAI444ggyMjLqvZaCgoIGWW+8SOT9175nBF1GVJVu/pTvfnIH/bw9Zx5xPu82eZtvSr/h0KRDGdFmBClfpZDxVUbQZTa4RPzuRURE6iLI2e2+C5xvZucBLYBWZvaku18R+SZ3nwHMAEhPT/fBgwfXeyEZGRk0xHrjRSLvv/Z9cNBlRE1mXiETHtzBj/kOPa68n+9vzOTewb8NuqxAJNp3LyIiUleBDbdz9ynu3sXduwJjgNf2DkgiIvUhd827XPXgG+TuOohuE5+kt6b5FhERkRroPkki0qhtXvEqredfzvV8j74THyK1s2YpExERkZrFREhy9wwgI+AyRKSR2bziVVrOv5xs2nP82HvorYAkCWrBsiymLV5Ldn4RndokM3loT934VESkBjERkkRE6tvmlf8uD0i7xr5A72N7BF2SSCAWLMtiyvyVFJWUApCVX8SU+SsBFJRERKoR9BTgIiL1LvPrfHbPv55sDlNAkoQ3bfHa8oBUpqiklGmL1wZUkYhI7FNIEpFGJTOvkDEPL+XHfgu7xr6ogCQJLzu/qE7tIiKikCQijcjmlf9m8V9/SkFxCXdNHK2AJAJ0apNcp3YREVFIEpFGYvPKf3PIc5dzxu63mHN1qmaxEwmbPLQnyc2SKrQlN0ti8tCeAVUkIhL7NHGDiMS9TStfo+Vzl5NLO3aNfYE+3XQxukiZsskZNLudiEjtKSSJSFwLBaTLygNSr2OPDbokkZgzKq2zQpGISB1ouJ2IxK3MvEL+8uJbZIWn+VZAEhERkfqgkCQicWljzibGzHiPF3afwrfXZiggiYiISL1RSBKRuLNp5Wu0fHAA/Yo/YvaEk0k98rCgSxIREZFGRCFJROJK2TVIebTiprEjNYudiIiI1DuFJBGJG5tWvU7L5y5jE235duwL9Dr2uKBLEhERkUZIIUlE4kLO+tW0nDdGAUlEREQanEKSiMS8zLxCRs/N5SEuUEASERGRBqf7JIlITNu0+g1+/sJGCkrac8bEe+mla5BERESkgSkkiUjM2rQqg5bzLuXn9CBl4kuapEFERESiQiFJRGJSWUDazKG0vvwR9SCJiIhI1CgkiUjMiQxIxZe/QK/jegZdkoiIiMSABcuymLZ4Ldn5RXRqk8zkoT0Zlda53rejkCQiMSUzr5Cs+VPpoIAkIiIiERYsy2LK/JUUlZQCkJVfxJT5KwHqPShpdjsRiRmZeYWMmfEeP93zU4rHvqiAJCIiIuWmLV5bHpDKFJWUMm3x2nrflkKSiMSETavfYMMD51NaXMDMiYM1zbeIiIhUkJ1fVKf2A6GQJCKB27T6DQ559lKO2rORWWN7aRY7ERERqaRTm+Q6tR8IhSQRCVRZQNpCG4ou141iRUREpGqTh/YkuVlShbbkZklMHlr/w/M1cYOIBGbT6jcrBiRdgyQiIiLVKJucQbPbiUijlZlXyC9f+IKf05WUyx9TQBIREZF9GpXWuUFC0d4UkkQk6rI3rGXMnI0UlHTkoAkv06tLm6BLEhERESmna5JEJKo2rX6L1rO+x+hv5zN7wsmkKiCJiIhIjFFPkohEzabVb3HIsxezhdZ8/7Ib6aVZ7ERERCQGqSdJRKKiLCBtpRVFl71Ar569gi5JREREpEoKSSLS4DbmbKLFs2PYSisKL3tRAUlERERimkKSiDSozLxCLn1sNb/ixwpIIiIiEhd0TZKINJhNn7zNX59/k4Ld6fxg4g26BklERETigkKSiDSITZ+8TcozFzORQxl77XWkKiCJiIhInNBwOxGpd2UBKY+WlFz2HKlHHhZ0SSIiIiK1ppAkIvWqLCB9Q0uKLnuBnj17B12SiIiISJ0oJIlIvcnMK+Sl52bxDS0pVEASERGROKVrkkSkXmRu2cGYh5ZQUHox37nydnofc3TQJYmIiIjsF/UkicgB2/TJf9jzl5NpX7yB2RO/o4AkIiIicU0hSUQOyKZP/sPBz4ymie/m3su+o1nsREREJO5puJ2I7LeygLSNQyi8bAE9e/YJuiQREQlbsCyLaYvXkp1fRKc2yUwe2pNRaZ2DLkskLigkich+yf3sI1LCAWnnZS8oIImIxJAFy7KYMn8lRSWlAGTlFzFl/koABSWRWtBwOxGps8y8QsY+t4nXOUkBSUQkBk1bvLY8IJUpKill2uK1AVUkEl/UkyQidZL72YeMf24TW3Y1o/vEJ+ipa5BERGJOdn5RndpFpCKFJBGptdxP3yFl7kX8lDSOmjhHkzSIiMSoTm2SyaoiEHVqkxxANSLxR8PtRKRWygLSdlI4Zsw0BSQRkRg2eWhPkpslVWhLbpbE5KE9A6pIJL6oJ0lE9in303c4eO5otpNCwZgX6Nmrb9AliRwQM+sO3Aa0dvfRQdcjUt/KJmfQ7HYi+0chSURqlLm1gJJnrqOUgxWQJCaY2SPAcGCzu6dGtJ8LTAeSgIfc/XfVrcPd1wPXmtm8hq5XJCij0jorFInsJ4UkEalWZl4hY2Z+QFufzB/GpCsgSayYBTwAPF7WYGZJwF+As4GNwBIze5FQYLpnr8+Pd/fN0SlVRETikUKSiFQpd827ZDw7g4LSMTw4cZRmsZOY4e5vmlnXvZoHAp+He4gws6eBke5+D6FeJxERkVpTSBKRSnLXvMvBT1/EGRxM+hW301sBSWJfZyAz4vVG4OTq3mxm7YD/AdLMbEo4TO39nuuA6wCOOOIIMjIy6rXgvRUUFDT4NuKNjkllOiaV6ZhUpmNSWV2PiUKSiFRQFpB2cDAFYxbQu0f3oEsSqXfuvhW4fh/vmQHMAEhPT/fBgwc3aE0ZGRk09DbijY5JZTomlemYVKZjUlldj4mmABeRcrlr3uPgp0eXB6SevVL3/SGR2JAFHBnxuku4TUREpM4UkkQECE3S8Mfn3uAbWiogSTxaAhxrZt3MrDkwBngx4JpERCROKSSJCBtzNjFmxnu8sjuNHde+rYAkMc3M5gDvAj3NbKOZXevuu4EfA4uBT4Fn3H11kHWKiEj8CuyaJDM7ktD0rUcADsxw9+lB1SOSqHLXvEfK0xdzMtcxfuJPSNUkDRLj3P2yatpfAl6KcjkiItIIBTlxw27g5+6+1MxaAh+Z2b/c/ZMAaxJJKKFrkC6igBb84FJN8y3RZ2YpQLG7lwZdi4iISJnAhtu5e467Lw0/30FoeIRuCy0SJSVff14ekHZcuoCevfsFXZIkADNrYmaXm9kiM9sMrAFyzOwTM5tmZj2CrlFERCQmrkkK3xQwDXg/4FJEEkLWV+sZuPoOBSQJwuvAMcAUoIO7H+nuhwOnAu8B95rZFUEWKCIiss/hdmZ2r7v/cl9t+8vMDgGeA25y9+1VLG/wm/kl+g23Enn/E3Hfvy7cw+8+KGJM6Ui69TuNwzZtJWdTRtBlRV0ifvdlAt73s9y9ZO9Gd88jdC54zsyaRb+sqpnZCGBEjx7q4BIRSSS1uSbpbGDvQPT9KtrqLHwifA6Y7e7zq3pPNG7ml+g33Erk/U+0fc9d+wEPPbea3daZ9umjGT3yzKBLCkyiffeRgtz3vQOSmbUArgCSgafcfWtVISoo7r4QWJienj4x6FpERCR6qg1JZvZD4Aagu5mtiFjUEvjPgW7YzAx4GPjU3e870PWJSM1y135A8pwLmMLh+ITX2fL58qBLEgGYTuicUgwsAE4LtBoRERFqvibpKWAEoZvxjYh4nOju9TFe/LvAlcAZZrY8/DivHtYrInspC0iFHETTSx8jtUuboEuSBGVmc8zsmIimtsCzhEYVHBpMVSIiIhVV25Pk7tuAbcBlZpZE6H5GTYFDzOwQd//qQDbs7m8DdiDrEJF9iwxI2y9dQM/exwddkiS224C7zSwH+A3wB+B5oAUwNcC6REREytVm4oYfEzpxbQL2hJsd0G9aIjEuM6+QtXPvpK8CksQId18PXG5mpwJzgUXAMN0nSUREYkltJm64Cejp7lsbuBYRqUeZeYWMmfEeJXuu54nLjqFnzz5BlySCmR0KXA6UABcDI4HFZjY9PEmCiIhI4Gpzn6RMQsPuRCRO5K5dQuYDw7HibTwy8XQFJIklC4B8QiMSnnD3Jwhd75pmZgpJIiISE2rTk7QeyDCzRcC3ZY2akU4kNuWuXUKLORfQjeY8fGkPenZuHXRJIpHaAfMITfn9AwB3LwLuMrOOQRYmIiJSpjYh6avwo3n4ISIxqiwgFdGcHZc+T8/e/YIuSWRvvwZeBkqBWyIXuHtOIBWJiIjsZZ8hyd3vBDCzg929sOFLEpH9kfvZh7SYcwHF4YB0XO8Tgi5JpBJ3f47QdN9xwcxGACN69OgRdCkiIhJF+7wmycxOMbNPgDXh1yeY2V8bvDIRqbXMvEJ+PG8d6+nMdgUkiWFmNtPMUqtZlmJm481sbLTrqo67L3T361q31rBVEZFEUpvhdvcDQwndVBZ3/9jMTm/IokSk9rK//IzLnvqKHSXtaD5hMcfpRrES2/4C3GFm/YBVwNeE7pF0LNAKeASYHVx5IiIitQtJuHumWYX7vup+FiIxIPezD0l+6gKu5kxOmTidVE3SIDHO3ZcDl5jZIUA60BEoAj5197VB1iYiIlKmNiEp08wGAW5mzYBJwKcNW5aI7EvuZx9y0FMX8C3N+N4lP+U4BSSJI+5eAGQEXYeIiEhVanOfpOuBHwGdgSygf/i1iASkLCDtoinbL3me4/roGiQRERGR+lKb2e22ADFzEa1IosvcnEfzpy5hF03ZdskCBSQRERGRelZtSDKzm93992b2v4TujF6Bu/+kQSsTkUoy8woZ88hyUvkhP7/kHAUkiUtmlgTc6+6/CLoWERGRqtTUk1R23dGH0ShERGqWu24pDz/zIgW7B3HjxOt0DZLELXcvNbNTg65DRESkOtWGJHdfGP73seiVIyJVyV23lINmj+QHNOPia66nrwKSxL9lZvYi8Cyws6zR3ecHV5KIiEhITcPtFlLFMLsy7n5+g1QkIhWUBaRdNGXHJfPpe3SHoEsSqQ8tgK3AGRFtDigkiYhI4GoabveH8L8XAh2AJ8OvLwM2NWRRIhJSFpBKSGLbJc9zXJ/+QZckUi/c/Zqga6gNMxsBjOjRo0fQpYiISBRVOwW4u7/h7m8A33X3S919YfhxOXBa9EoUSUyZeYU8O/cxSkgi/5IFCkjSqJhZFzN73sw2hx/PmVmXoOvaW/i8d13r1hriKiKSSGpzn6QUM+te9sLMugEpDVeSiGRu2cGYGe/xUOkwtl75hgKSNEaPAi8CncKPheE2ERGRwNUmJP0UyDCzDDN7A3gduKlBqxJJYLnrllL6l5PpUryO2RNOpvcxRwddkkhDaO/uj7r77vBjFtA+6KJERESgdjeTfdnMjgV6hZvWuPu3DVuWSGIquwapCUn85uKBmuZbGrOtZnYFMCf8+jJCEzmIiIgEbp8hKexYoCeh2YhOMDPc/fGGK0sk8eSsW0aL8CQN+Rc/z3F904IuSaQhjQf+F/gToVnt3gHiYjIHERFp/PYZkszs18BgoA/wEvB94G1AIUmknmSv/0QBSRKGmSUBv9WtJEREJFbV5pqk0cCZQG54ytYTAI0BEqknmXmFjJ37FW8yQAFJEoK7lwJHm1nzoGsRERGpSm2G2xW5+x4z221mrYDNwJENXJdIQsj5/GMmPruBvF3J9Jj4mK5BkkSyHviPmb0I7CxrdPf7gitJREQkpDYh6UMzawPMBD4CCoB3G7IokUSQs24ZB80eya0cQ9uJC0hVQJLE8t/wownQMuBaREREKqgxJJmZAfe4ez7wdzN7GWjl7iuiUZxIY1UWkEoxOlz8R/UgSUIJX5N0nLuPDboWERGRqtR4TZK7O6HJGspeb1BAEjkwkQHpm4uf57i+A4IuSSSqdE2SiIjEutoMt1tqZie5+5IGr0akkcvcupP8p35ABwUkEV2TJCIiMas2IelkYKyZfUnoRGaEOpmOb9DKRBqZzLxCxsx8n0N8Ev97cW8FJEl0cXFNkpmNAEb06NEj6FJERCSKahOShjZ4FSKNXM66Zbz99H0Ull7OgxOH6xokSXjufufebWZW2xucR427LwQWpqenTwy6FhERiZ593ifJ3b909y+B3YTuiu5AVkMXJtJY5KxbRvPZozi79A2evqyrZrGThGZmb0c8f2KvxR9EuRwREZEqVftXOzObAjRz97vCTe8C+UBz4DHgngavTiTO5Xy+nOazR+E4eRc/T8+evYMuSSRoKRHPU/daZtEsREREpDo19SRdDPwx4vXW8HVIfYFhDVqVSCOQ8/lymj85sjwgHdf3xKBLEokFXs3zql6LiIgEosbx3+6+M+Ll9HBbqZklN2hVInEuM6+Q38/N4Jc0p/DipxWQRP5fGzO7gNAf6dqY2YXhdgM0FlVERGJCTSHpEDNr5u4lAO4+C8DMDgJaRaE2kbi0MXczY2atomB3H64f/w59j2ofdEkiseQN4PyI5yMilr0Z/XJEREQqqykkzQMeNLMfu3shgJmlAA+El4nIXnI+/5iDnhzFEC5jzMRf0leTNIhU4O7XBF2DiIjIvtR0TdLtwGbgKzP7yMw+AjYAm8LLRCRCzucf0+zJkRilXD36As1iJyIiIhKnqu1JcvdS4BYzuxMou4ve5+5eFJXKROJIWUCCPXwz+jmOTT0p6JJEREREZD/t88Z94VC0Mgq1iMSlrOwsDlJAEhEREWk09nkzWRGpXmZeIZc8vpaHuEABSaQOzOxgM7vdzGaGXx9rZsODrktERARq0ZMkIlXL+e8Kps59j4KSrgyfOJVjdQ2SSF08CnwEnBJ+nQU8C/wjsIpERETC9tmTZGb/rk2bSCLJ+e8Kmj1xPrfvuo/Z40/UJA0idXeMu/8eKLvNRCGheyWJiIgErtqeJDNrARwMHGZmh/L/J69WQOco1CYSMxatX8T0pdPJ3ZlL+4PaMTEzk7PZQ8lFj5N6ZLugyxOJR7vCNyZ3ADM7Bvg22JJERERCahpu9wPgJqAToSERZSFpO6F7JYkkhEXrFzH1nakUlxYDsPnbLfzhsBbsPOlGru03MODqROLWVOBl4Egzmw18FxgXZEEiIiJlqh1u5+7T3b0b8At37+7u3cKPE9xdIUkSxvSl08sDUplvmxhztywKqCKR+OfurwAXEgpGc4B0d88IsqaqmNkIM5uxbdu2oEsREZEoqs0U4P9rZoOArpHvd/fHG7AukZiRuzO3Tu0ism9mthB4CnjR3XcGXU913H0hsDA9PX1i0LWIiEj01GbihieAPwCnAieFH+kNXJdIzGh/UNXXHHVI6RDlSkQalT8ApwGfmNk8MxsdvhZWREQkcLWZAjwd6OPu3tDFiMSanP+uZELmRv542EF82+T/J95qkdSCSQMmBViZSHxz9zeAN8wsCTgDmAg8QmhyIBERkUDV5mayqwD9yVwSTs5/V9L0ifMZurOYH3abSMeUjhhGx5SOTB00lWHdhwVdokhcC89udxFwPaFRCo8FW5GIiEhIbXqSDiM0HOIDIqZndffzG6wqkYDlrF9F0yfOJ4ndbL3oOa7tN5BrUc+RSH0xs2eAgYRmuHsAeMPd9wRblYiISEhtQtLUhi5CJJZk5hXy47lr+DWH0/Ki/+VYTfMt0hAeBi5z99KgCxEREdlbbWa3eyMahYjEguwvP2fsU1+wbVcrmk9YzLFd2gRdkkijYmZnuPtrQAow0swqLHf3+YEUJiIiEqHakGRmb7v7qWa2g/Ad0csWAe7uurhWGpWc9ato+vj53EA6qRNnktq5ddAliTRG3wNeA0ZUscwBhSQREQlctSHJ3U8N/9syeuWIBCNn/SqSHj+fpuzixIt+zrEKSCINwt1/HX56l7t/EbnMzLoFUJKIiEgltZndDjM7wcx+HH4cX18bN7NzzWytmX1uZrfU13pF6iIyIG296DmO7Xdy0CWJJILnqmibF/UqREREqrDPa5LMbBKh+1eUDYGYbWYz3P1/D2TD4Xtj/AU4G9gILDGzF939kwNZr0hdbNlZQrcnLqElu9h64TwFJJEGZma9gL5AazO7MGJRK0A3kxURkZhQm9ntrgVOdvedAGZ2L/AucEAhidDUr5+7+/rwep8GRgIKSRIVmXmF3LOkhH78gNsuPJljj/9O0CWJJIKewHCgDRWvS9pB6A9yIiIigatNSDIgcorW0nDbgeoMZEa83gjoz/gSFTnrV/PEU7MpLv0uv7h+nK5BEokSd38BeMHMTnH3d4OuR0REpCq1CUmPAu+b2fPh16MI3d8iKszsOuA6gCOOOIKMjIx630ZBQUGDrDdeJNr+F+ZlMWDFr7jeS+jQ63i2rFtGxrqgq4q+RPve95bI+x8j+77MzH5EaOhd+TA7dx8fXEkiIiIhtblP0n1mlgGcGm66xt2X1cO2s4AjI153Cbftvf0ZwAyA9PR0Hzx4cD1suqKMjAwaYr3xIpH2P2f9appkjKc5JWy5aB7d84oTZt/3lkjfe1USef9jZN+fANYAQ4G7gLHAp4FWJCIiElbt7HZm1ir8b1tgA/Bk+PGlmR0annjhQCwBjjWzbmbWHBgDvHiA6xSpVs761TR5fATN2cWWC+fpGiSRYPVw99uBne7+GDAMDbkWEZEYUVNP0lOELq79iMo3kwU4xMxmuvut+7Nhd99tZj8GFgNJwCPuvnp/1iWyL5l5hcx+6gl+wC62XvisApJI8ErC/+abWSqQCxweYD0iIiLlarqZ7PDwv1Xe3C/ck7QK2K+QFF73S8BL+/t5kdrI3LKDMQ8toaD0DEZeeT29j+kadEkiAjPM7FDgdkKjCA4B7gi2pMrMbAQwokePHkGXIiIiUVRtSDKzATV90N2XAr3rvSKRepS9/hN2P3EJx/gPuHniVfTWLHYiMcHdHwo/fQPoHmQtNXH3hcDC9PR0TU8uIpJAahpu98fwvy2AdOBjQkPtjgc+BE5p2NJEDkz2F5/S5PERHEoxd1x4Ij0UkEQCZ2Y/q2m5u98XrVpERESqU9NwuyEAZjYfGODuK8OvU4GpUalOZD9lf/EpTR4bTguK2XLhs/Q4flDQJYlISMugCxAREdmX2twnqWdZQAJw91VmpmF2ErOyv/qvApJIjHL3O4OuQUREZF+qnQI8wgoze8jMBocfM4EVDV2YyP7IzCtkzOz/8jb9FZBEYpiZHWdm/zazVeHXx5vZr4KuS0REBGoXkq4BVgOTwo/VwLgGrElkv2R/sYYbHvwn23ZBr4mPKCCJxLaZwBTCU4G7+wpC98sTEREJ3D6H27l7MfCn8AMzOw24D/hRw5YmUnuha5BGcBdtaTbhFVI1SYNIrDvY3T8ws8i23UEVIyIiEqk21yRhZmnAZcAlwBfA/IYsSqQusr9YQ5PHRpBMIS0veIweXdoEXZKI7NsWMzuG8M3KzWw0kBNsSSIiIiE13SfpOELB6DJgCzAXsLJZ70SCtGj9IqYvnU7uzlzal5Ryfcpu0s55lh4nfDfo0kSkdn4EzAB6mVkWoT/AjQ22JBERkZCarklaA5wBDHf3U939f4HS6JQlUr1F6xcx9Z2p5OzMwXE2N2vCvUe0Y23L/KBLE5Facvf17n4W0B7oBXwPODXYqkREREJqCkkXEhr68LqZzTSzMwndTFYkUNOXTqe4tLhC27e+i+lLpwdUkYjUlpm1MrMpZvaAmZ0NFAJXA58TGtItIiISuJpuJrsAWGBmKcBI4CbgcDP7G/C8u78SlQpF9pK7M7dO7SISU54AvgHeBSYCtxH6A9wF7r48wLpERETK1WZ2u53AU8BTZnYocDHwS0AhSaIue8Na2peUsrlZ5U7QDikdAqhIROqou7v3AzCzhwiNWDgqPJOqiIhITKjNfZLKufs37j7D3c9sqIJEqpO9YS02axjXf7OTg6x5hWUtklowacCkgCoTkTooKXvi7qXARgUkERGJNbWaAlwkKFXPYvckd7bcVt7eIaUDkwZMYlj3YUGXKyL7doKZbQ8/NyA5/NoAd/dWwZUmIiISopAkMatsFruySRrKZrG7s+U2hnUfplAkEofcPSnoGkRERPalTsPtRKJJs9iJiIiISBAUkiRmaRY7EREREQmCQpLEpOwNa2m/u+p7F2sWOxERERFpSLomSWJO9oa18NhwfnhwCb87oh3f+q7yZZrFTkSkcViwLItpi9eSnV9EpzbJTB7ak1FpnYMuS0QEUE+SxJiNm7bgj43gEN9J/7Of4M5T76JjSkcMo2NKR6YOmqoJG0RE4tyCZVlMmb+SrPwiHMjKL2LK/JUsWJYVdGkiIoB6kiSGZOYVMubRFZzjoxg7agQ9+p9GD1AoEhFpZKYtXktRScUh1UUlpUxbvFa9SSISExSSJCZkf/kZv539CgUlx3HRxNvo0bl10CWJiEgDyc4vqlO7iEi0KSRJ4LK//AyfNZypXsKWa96jrwKSiEij1qlNMllVBKJObZIDqEZEpDJdkySBKgtILX0HO0Y9Rt+jjwi6JBGRcmY2wsxmbNu2LehSGpXJQ3uS3KzifYWTmyUxeWjPgCoSEalIPUkSVYvWL2L60unk7syl/UHtGLcxh5G+g82j5tKj/+lBlyciUoG7LwQWpqenTwy6lsak7LojzW4nIrFKIUmiZtH6RUx9ZyrFpcUAbP52C9PbNWXXgJ9zrQKSiEhCGZXWWaFIRGKWhttJ1ExfOr08IJX5tokxN++VgCoSEREREalMIUmiJndnbp3aRURERESCoJAkUdP+oMOqbO+Q0iHKlYiIiIiIVE8hSaIi+8t1jNuYzUF7vEJ7i6QWTBowKaCqREREREQqU0iSBpf95ef4rGGMKsjnhqPH0TGlI4bRMaUjUwdNZVj3YUGXKCIiIiJSTrPbSYPKzCvkB0+t4jfeluKRMxmf9j3G84ugyxIRERERqZZCkjSYrMwNXPHkWvJ3HcxBE17mmC5tgi5JRERERGSfFJKkQWR/uQ6fNZzJfixdJ84mtXProEsSEREREakVhSSpd9lfrmPPrOG09m30HvkLjlFAEhEREZE4ookbpF5FBqTNI5/mmLTBQZckIiIiIlInCklSbzK37mTbrDGhgHT+HAUkEREREYlLGm4n9SIzr5AxM9/nKJ/A3ef34pgBQ4IuSURERERkvygkyQHL/vJznn5iBgWlZ3HbxMt1DZKIiIiIxDWFJDkg2V9+Tums4fzQ8xlx+bX0UkASERERkTina5Jkv5UFpEM9n03nz6HXcT2DLklERERE5ICpJ0lqbdH6RUxfOp3cnbm0P+gwrt6Yw4XhgKRrkERERESksVBIklpZtH4RU9+ZSnFpMQCbv/2aP7dLYnfazxivgCQiIiIijYiG20mtTF86vTwglfm2ifH0N68GVJGIiIiISMNQSJJayd2ZW6d2EREREZF4pZAktdL+oMOqbO+Q0iHKlYiIiIiINCyFJNmnrK/Wc9XGHFrs8QrtLZJaMGnApICqEhERERFpGApJUqOsrI2UPnoeowvy+OFRV9IxpSOG0TGlI1MHTWVY92FBlygiIiIiUq80u51UKzOvkMseX8OP/AQGnn8d4wecyXh+GXRZIiIiIiINSiFJqpT11XomzV7Cjl2H0m/iDI7p3DrokkREREREokIhSSrJ+mo9pY8O4x5vRsmEN0hVQBIRERGRBKKQJBWUBaR2nsemEU/Ss8uhQZckIiIiIhJVmrhBymVlflEhIHU/8eygSxIRERERibpAQpKZTTOzNWa2wsyeN7M2QdQh/y8zr5BPHv2xApKIiIiIJLygepL+BaS6+/HAZ8CUgOoQ4OvCPYyZ8R537RlP7qhnFZBEREREJKEFEpLc/RV33x1++R7QJYg6JDTEzj/4K7uKC/nbxLM4pv/pQZckIiIiIhKoWLgmaTzwz6CLSERZmV+w+5FhDPc3mXNhW81iJyIiIiJCA85uZ2avAh2qWHSbu78Qfs9twG5gdg3ruQ64DuCII44gIyOj3mstKChokPXGsh35W+i3/HYO9628cvTNtMnbxcYEOwaQmN99mUTed0js/U/kfa8rMxsBjOjRo0fQpYiISBQ1WEhy97NqWm5m44DhwJnu7jWsZwYwAyA9Pd0HDx5cj1WGZGRk0BDrjVVZmV9Q8sgNHM5Wckc8SZuC5gm1/5ES7buPlMj7Dom9/4m873Xl7guBhenp6RODrkVERKInqNntzgVuBs5398IgakhUmXmFTHnidZr5LnKHP0n39HOCLklEREREJKYEdTPZB4CDgH+ZGcB77n59QLU0aovWL2L60unk7sylfYv2bM8Zyq6S/uSPf5e+R7UPujwRERERkZgTSEhydw3ujoJF6xcx9Z2pFJcWA7C5eDNNW8/mh6d2U0ASEREREalGLMxuJw1k+tLp5QGpzO4me5i/YWZAFYmIiIiIxD6FpEYsd2dundpFREREREQhqVFr3+LwKts7pFQ1M7uIiIiIiIBCUqOVmVfI9pxzSNqTVKG9RVILJg2YFFBVIiIiIiKxTyGpEcra+CW//9sMdm3rz/X9bqNjSkcMo2NKR6YOmsqw7sOCLlFEREREJGYFNQW4NJCsjV+y6+Fh3O15ZF39AX26deH69IuDLktEREREJG6oJ6kRKQtIHXwzecMeoU+3LkGXJCIiIiISdxSSGonIgJQ77HG6nXRu0CWJiIiIiMQlDbdrBDLzCnlp1u+4SgFJRETkgC1YlsW0xWvJzi+iU5tkJg/tyai0zkGXJSJRpJAU5zLzChkz4z127h7JkIsncFzfE4MuSUREJG4tWJbFlPkrKSopBSArv4gp81cCKCiJJBANt4tjWRu/YuMDw2hdnMWTE09RQBIRETlA0xavLQ9IZYpKSpm2eG1AFYlIENSTFKeyNn7Ftw+fR3/fxAPDD6d759ZBlyQiIhL3svOL6tQuIo2TepLiUFlA6uibyB32ON3ThwZdkoiISKPQqU1yndpFpHFSSIozewekbid9P+iSREREGo3JQ3uS3CypQltysyQmD+0ZUEUiEgSFpDiSmVfINU98zBZvRe55jykgiYiI1LNRaZ2558J+dG6TjAGd2yRzz4X9NGmDSILRNUlxIisrk6sfX8XWXQdx8ISX6NalTdAliYiINEqj0jorFIkkOIWkOJC18SuKHx7GHd6ewyY+T6omaRARERERaTAKSTGuLCB18lyanPc7uikgiYiIiIg0KF2TFMPKAlJnzyH3vFl0Gzgs6JJERERERBo9haQYlZlXSPYjY+nsOeSc95gCkoiIiIhIlGi4XQzKzCtkzIz3OGLPOP54Xke6DTwv6JJERERERBKGQlKMycrKZN6jf6Jg91AenDha1yCJiIiIiESZQlIMycrKpOihYfzQsxl2ydUcp4AkIiIiIhJ1uiYpRpQFpC6eTc55j3FcnxOCLklEREREJCEpJMWArKxMCh8aHg5ImsVORERERCRICkkBy8wr5E+z5tDJc8n5/qN0Gzg86JJERERERBKarkkKUObWAsbM/ICC3f259sr36H1Mt6BLEhERERFJeOpJCkhWViY7HvgeJxa/y+wJJysgiYiIiIjECPUkBSArayOFDw2nu2fxs+8fT1fNYiciIiIiEjPUkxRloYA0jCM9i9zvP0rXk0cEXZKIiIiIiERQT1IUbczdTJECkoiIiIhITFNPUpRk5hVy6aMr+Y/3U0ASEREREYlh6klqQIvWL2L60unk7swlaXdLdjcbTvpVf9E1SCIiIiIiMUw9SQ1k0fpFTH1nKjk7c3Cc3U230/yI5/jy27eDLk1ERERERGqgkNRApi+dTnFpcYW2XXu+ZfrS6QFVJCIiIiIitaGQ1EByd+bWqV1ERERERGKDQlIDyMwrpNXupCqXdUjpEOVqRERERESkLhSS6llmXiFjZrxHcd4oDrJmFZa1SGrBpAGTAqpMRERERERqQyGpHmVlZ/HBX66htHgHsy6/iTtP/Q0dUzpiGB1TOjJ10FSGdR8WdJkiIiIiIlIDTQFeT7Kys9gxczjDPZO0ERPp3rk1qQxTKBIRERERiTPqSaoHZQGpm2eSc+7DdD/xrKBLEhERERGR/aSQdID2DkhdvzMy6JJEREREROQAKCQdgMy8Qn4263Va+g4FJBERERGRRkLXJO2njZu2MOaRjykoOZz88e/S96j2QZckIiIiIiL1QCFpP2zMzqZg5jCu9X4MnPi/9O3cOuiSRERERESknigk1VFZQOruX9Fi6B10VUASEREREWlUdE1SHYQC0nC6+1dkD32IrqdcEHRJIiIiIiJSzxSSailzawHfzBxFd/9SAUlEREREpBHTcLtayMwrZMzMDxjoI5g0NFUBSURERESkEVNI2oeN2dnc9+gcCnb349qJP9E1SCIijYCZjQKGAa2Ah939lWArEhGRWKLhdjXYmJ3Njpkj+G3JvTw99hhSFZBERAJnZo+Y2WYzW7VX+7lmttbMPjezW2pah7svcPeJwPXApQ1Zr4iIxB/1JFWjLCAd4xvIOWcGvXscE3RJIiISMgt4AHi8rMHMkoC/AGcDG4ElZvYikATcs9fnx7v75vDzX4U/JyIiUk4hqQplAamHf0H2OTM5etBFQZckIiJh7v6mmXXdq3kg8Lm7rwcws6eBke5+DzB873WYmQG/A/7p7ksbuGQREYkzCkl7ycwr5NlH7uNGBSQRkXjSGciMeL0ROLmG998InAW0NrMe7v73vd9gZtcB1wEcccQRZGRk1F+1VSgoKGjwbcQbHZPKdEwq0zGpTMeksroek0BDkpn9HPgD0N7dtwRZC4RnsZvxHgW7z2bE6Es5NvWkoEsSEZEG4O5/Bv68j/fMAGYApKen++DBgxu0poyMDBp6G/FGx6QyHZPKdEwq0zGprK7HJLCJG8zsSOAc4Kugaoi0MSeHLx84n8OLv2D2xO8oIImIxJcs4MiI113CbSIiInUW5Ox2fwJuBjzAGgDI376D7TOGc/KepfxpaFvNYiciEn+WAMeaWTczaw6MAV4MuCYREYlTgYQkMxsJZLn7x0FsP9LGnBx6LP01x/oXZJ8zQzeKFRGJcWY2B3gX6GlmG83sWnffDfwYWAx8Cjzj7quDrFNEROJXg12TZGavAh2qWHQbcCuhoXa1WU+DXTibv30HPZb+mp6+gTeOuZmmuw7jiwS8yC2RL+7TvmcEXUZgEnn/433f3f2yatpfAl6KcjkiItIINVhIcvezqmo3s35AN+Dj0AysdAGWmtlAd8+tYj0NcuFsZl4h4x7M4C4O4c1jJnPWVTXed7BRS+SL+7Tvg4MuIzCJvP+JvO8iIiK1EfXZ7dx9JXB42Wsz2wCkR3N2u405uYyf9SFbdjWn9YQXKfl8ebQ2LSIiIiIiMS7h7pO0MSeXbTOG81tvRvKEf5LapQ0ZnwddlYiIiIiIxIrAQ5K7d23obSxav4jpS6eTuzOXdrudmw7Oo/+g+zm6S5uG3rSIiMQxMxsBjOjRo8d+r2PBsiymLV5Ldn4RndokM3loT0alda6/IkVEpN4FOQV4VCxav4ip70wlZ2cOjrOlKfzmiCNY1fHgoEsTEZEY5+4L3f261q3379YQC5ZlMWX+SrLyi3AgK7+IKfNXsmCZbuEkIhLLGn1Imr50OsWlxRXavvUSpi+dHlBFIiKSKKYtXktRSWmFtqKSUqYtXhtQRSIiUhuNPiTl7qw0YV6N7SIiIvUlO7+oTu0iIhIbGn1I6pBS1a2aqm8XERGpL53aJNepXUREYkOjD0mTBkyiRVKLCm0tklowacCkgCoSEZFEMXloT5KbJVVoS26WxOShPQOqSEREaiPw2e0a2rDuwwDKZ7frkNKBSQMmlbeLiIg0lLJZ7DS7nYhIfGn0IQlCQUmhSEREgjAqrbNCkYhInGn0w+1ERERERETqQiFJREREREQkgkKSiIhINcxshJnN2LZtW9CliIhIFCkkiYiIVMPdF7r7da1btw66FBERiSKFJBERERERkQgKSSIiIiIiIhEUkkRERERERCIoJImIiIiIiERQSBIREREREYmgkCQiIiIiIhJBIUlERERERCSCuXvQNdSamX0NfNkAqz4M2NIA640Xibz/2vfElcj7X92+H+3u7aNdTDxowPNPpET+mayOjkllOiaV6ZhUpmNSWVXHpNrzXlyFpIZiZh+6e3rQdQQlkfdf+56Y+w6Jvf+JvO+xTN9LZTomlemYVKZjUpmOSWV1PSYabiciIiIiIhJBIUlERERERCSCQlLIjKALCFgi77/2PXEl8v4n8r7HMn0vlemYVKZjUpmOSWU6JpXV6ZjomiQREREREZEI6kkSERERERGJoJC0FzP7uZm5mR0WdC3RYmbTzGyNma0ws+fNrE3QNUWDmZ1rZmvN7HMzuyXoeqLFzI40s9fN7BMzW21mk4KuKdrMLMnMlpnZP4KuJdrMrI2ZzQv/N/+pmZ0SdE1SUSKeh6qTqOenqiTqOas6OpdVL5HPcVXZ3/OeQlIEMzsSOAf4KuhaouxfQKq7Hw98BkwJuJ4GZ2ZJwF+A7wN9gMvMrE+wVUXNbuDn7t4H+A7wowTa9zKTgE+DLiIg04GX3b0XcAKJexxiUgKfh6qTcOenqiT4Oas6OpdVL5HPcVXZr/OeQlJFfwJuBhLqQi13f8Xdd4dfvgd0CbKeKBkIfO7u6919F/A0MDLgmqLC3XPcfWn4+Q5C/7PoHGxV0WNmXYBhwENB1xJtZtYaOB14GMDdd7l7fqBFyd4S8jxUnQQ9P1UlYc9Z1Un0c1l1EvkcV5UDOe8pJIWZ2Uggy90/DrqWgI0H/hl0EVHQGciMeL2RBPyfq5l1BdKA9wMuJZruJ/RL6J6A6whCN+Br4NHwUIyHzCwl6KIkROehfUqU81NVdM6qQYKey6pzP4l7jqvKfp/3mjZsXbHFzF4FOlSx6DbgVkJDHBqlmvbd3V8Iv+c2Qt3Xs6NZmwTDzA4BngNucvftQdcTDWY2HNjs7h+Z2eCAywlCU2AAcKO7v29m04FbgNuDLStxJPJ5qDo6P8mBSMRzWXV0jqvSfp/3EiokuftZVbWbWT9CSfNjM4NQd/5SMxvo7rlRLLHBVLfvZcxsHDAcONMTY174LODIiNddwm0JwcyaETqpzHb3+UHXE0XfBc43s/OAFkArM3vS3a8IuK5o2QhsdPeyv7bOI3SykChJ5PNQdXR+qpWEPmdVJ4HPZdVJ9HNcVfb7vKf7JFXBzDYA6e6+JehaosHMzgXuA77n7l8HXU80mFlTQhcBn0noRLMEuNzdVwdaWBRY6Dewx4A8d78p4HICE/4r2y/cfXjApUSVmb0FTHD3tWY2FUhx98kBlyV7SbTzUHUS8fxUlUQ+Z1VH57KaJeo5rir7e95LqJ4kqdYDwEHAv8J/wXzP3a8PtqSG5e67zezHwGIgCXgkgU423wWuBFaa2fJw263u/lJwJUkU3QjMNrPmwHrgmoDrEalJwp2fqpLg56zq6FwmtbVf5z31JImIiIiIiETQ7HYiIiIiIiIRFJJEREREREQiKCSJiIiIiIhEUEgSERERERGJoJAkIiIiIiISQSFJGpSZlZrZ8ohHVzN7pxafe8jM+oSf37of280ws7VmtsLM1pjZA2bWJmL5OxHPp5nZ6vC/7c3sfTNbZman1XW70WJmPwvv10oz+9jM7gvfVG9/19fVzFaFn6eb2Z8PYF11/r5ERGKJzl0NQ+cuiSeaAlwalJkVuPsh0V6HmWUQuonah+F58e8hdGPG71Xx3m1AW3cvNbMxwFnuPqEO20py99K61HcgzOx6YBQwxt3zw/v3M+Cv7r59f2ozs67AP9w9tR7qO+DvXEQkSDp31T+duyTeqCdJos7MCsL/Dg7/1Wxe+C9Ls8N30C77a1q6mf0OSA7/JW92eNkVZvZBuO1BM0uqaXvuvgu4GTjKzE7Yq4YXgUOAj8zsl8DvgZHhdSeb2Tlm9q6ZLTWzZ83skPDnNpjZvWa2FLh4H++7M9y+0sx6hdsPMbNHw20rzOyicHuV69nLbcAP3T2/bP/c/XdlJxkzKzCzP5rZx8ApZnaHmS0xs1VmNiPiGJ8Y/kvex8CPIr6fwWb2j/DzFDN7JHy8l5nZyHD7ODObb2Yvm9k6M/t9uL3S9yUi0hjo3KVzlyQYd9dDjwZ7AKXA8vDj+XBbQfjfwcA2oAuhwP4ucGp4WQahv56Vvz/8vDewEGgWfv1X4Koqtlv++Yi2BcClVawz8vk44IHw88OAN4GU8OtfAneEn28Abq7l+24MP78BeCj8/F7g/ojtHlrTeiLe1wr4Zh/H3IFLIl63jXj+BDAi/HwFcHr4+TRgVcT38o/w898CV4SftwE+A1LCx2k90BpoAXwJHLn38dRDDz30iMeHzl06d+mhR1NEGlaRu/evYfkH7r4RwMyWA12Bt2t4/5nAicCS8B+VkoHNtazFavm+Mt8B+gD/CW+rOaGTYZm5tXzf/PC/HwEXhp+fBYwpe4O7f2Nmw/exnso7ZDaU0EmrDXC5u79D6OT+XMTbhpjZzcDBQFtgtZm9BbRx9zfD73kC+H4VmzgHON/MfhF+3QI4Kvz83+6+LVzHJ8DRQGZN9YqIxAmdu3TukgSnkCRB+zbieSn7/pk04DF3n1KXjYSHNfQDPq3Lx4B/uftl1SzfWcv3le3jvvZvX+vB3beHhyR0c/cv3H0xsDg8xKB5+G3FHh7LbWYtCP3FMt3dM81sKqGTRW0ZcJG7r63QaHYydf/uREQaC527ar89nbskLumaJIkHJfb/s9/8GxhtZocDmFlbMzu6pg+HP3sPkOnuK+qw3feA75pZj/B6UszsuAN4X6R/UXEs9aF1WM89wN8sPONReJx2dSePsvYt4THiowE8NCY838xODS8fW83nFwM3RowFT9vHfkHF70tEJFHp3FWRzl0SVxSSJB7MAFaY2Wx3/wT4FfCKma0g9D/sjtV8bnb4PasIjUUeWZeNuvvXhMYvzwmv512g1/6+by93A4eGL0j9GBhSh/X8jdAJ9/3w+/4DLAs/9q4tH5hJ6BgsBpZELL4G+Et4qEh1wzl+AzQjdPxXh1/vS/n3VYv3iog0Vjp3VaRzl8QVTQEuIiIiIiISQT1JIiIiIiIiERSSREREREREIigkiYiIiIiIRFBIEhERERERiaCQJCIiIiIiEkEhSUREREREJIJCkoiIiIiISASFJBERERERkQj/Bzd2jUFCHqCAAAAAAElFTkSuQmCC\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 = 30\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+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABYEklEQVR4nO3deXxU9b3/8deHCCZsQVxYVVQQQRCCiKJioVZREEWKIoIbUay1Fnt/V1u0Wtvaa1t7qVi74QYiIoqARbyiVuNGFRSQRUCUohASlgAxCROB5PP7YyZ0AklIIJkzk3k/H495MHPOzDmfcyZ68s53OebuiIiIiIiISFiDoAsQERERERGJJwpJIiIiIiIiURSSREREREREoigkiYiIiIiIRFFIEhERERERiaKQJCIiIiIiEkUhSUQkipl1MDM3syMO8fOjzOz12q6rvjCzLDO7Oeg6REREqqKQJCJxzczWm1nIzArNLNfMJptZ06DrgooDlbtPc/eLY7Dvrmb2sZntiDzeNLOuUevNzH5nZnmRx+/MzKrY3j1m9u/Ied5oZjPq+hgOxsxuMLNPzOybSE2/ryi8mlknMys2s2f3W36tmX1lZkVmNsfMWlaxLzez5WbWIGrZg2Y2uZaOZZKZrTGzUjO7sYL1P4n8fH9jZk+Z2ZGVbKfsZ+7V/ZY/a2YPxEudIiKJTiFJRBLBEHdvCvQEMoDxwZYTFzYBw4GWwDHAP4Dno9aPBYYCPYAzgCHArRVtyMxuAK4Dvhc5z72Bf9ZV4TXQGLiT8PGdDVwI/HcF7/szsCh6gZmdDvyd8HG1AnYBfznI/toC1xxWxZX7FPghsHj/FWY2EPgZ4eM7ETgZ+OVBtne2mZ1b20VS+3WKiCQkhSQRSRjungvMJxyWADCzc8xsgZntNLNPzax/1LobzWydmRVEWklGRZY3MLOfR1oZtpjZM2aWXtE+Iy1Z34t6/UBUi8W7kX93Rlpg+kb2+X7U+881s0Vmlh/599yodVlm9msz+yBS4+tmdkw1z8VOd1/v7g4YUAJ0jHrLDcD/uvtGd88G/he4sZLNnQXMd/cvI9vOdfdJUXWmm9mTZpZjZtmRFpaUqPVjzGxVpEVrvpmdGLXuIjNbHTn+xyK1Vou7/9Xd33P33ZFjmAacF/0eM7sG2MmBoW4UMNfd33X3QuA+YJiZNatil78HfllJa9WIyM9Q88jrSyMtKsdW81j+7O7/BIorWH0D8KS7r3T3HcCvqfy7iq71NxWtMLMVZjYk6nVDM9tmZhkB1CkikpAUkkQkYZhZe+BS4IvI63bAPOBBwi0q/w28ZGbHmlkT4FHgUndvBpwLLI1s6sbIYwDhv4Y3BR47hJIuiPzbwt2buvu/9qu3ZaS+R4GjgQnAPDM7Oupt1wI3AccBjYhqKTGzZWZ2bVUFmNlOwr/Q/gn4n6hVpxNuFSjzaWRZRT4Erjezu8ysd3QAipgM7CUcwjKAi4GbI/u/ArgHGAYcC7wHTI+sOwaYBfyccGvQl0SFHDM7IRJuT6jqGKNcAKyM+nxz4FfAf1Xw3nLHHwmAu4FTq9j+LOAbKvjF391nAAuARyPf35PAze6+NVLLK2b2s2oeR5W1Rp632u/nZH9/AU6NDvBRngFGR70eBOS4+5IA6hQRSUgKSSKSCOaYWQGwAdgC/CKyfDTwqru/6u6l7v4G8DHhXwoBSoFuZpbm7jnuXvYL9ihggruvi7QyjAeuqagF4TANBta6+1R33+vu04HVhLu+lXna3T939xDwAlGtZO5+hrs/V9UO3L0FkA78CFgStaopkB/1Oh9oanbguCR3fxa4AxgIvANsMbOfAphZK8Ln8053L3L3LcAf+U+3tB8AD7n7KnffSzio9Yy0Jg0CVrr7THffAzwC5Ebt92t3b+HuX1d1jJE6xhDuBviHqMW/JtyysbGCj+x//GXnoKqWJCfc4nSfmTWqYP3twHeBLMKtVK9EHctl7v7bgx1HJSr6rjhIrSHCLUkPVrDuWWBQWasX4S6HUwOqU0QkISkkiUgiGBppDeoPnEa4VQLC4yKuirRG7Iy0qpwPtHH3ImAE4V/ic8xsnpmdFvlcW+CrqO1/BRxBeOxKbdp/P2X7ahf1Ojfq+S7Cv4jWSORY/wY8Y2bHRRYXAs2j3tYcKIx0z6toG9Pc/XtAC8Ln7NeRMSgnAg0Jn8Oyc/x3wi1fRNZPjFq3nXCXunaEj39D1D48+nV1mdlQ4CHCrYLbIst6At8jHNgqsv/xE3ldUNW+3P1VYCMVjN9y953Ai0A3wt0Xa0tF3xUcpFbgCcItOdGhG3ffBHwAfN/MWhBufZ0WYJ0iIglHIUlEEoa7v0O461dZa8IGYGqkNaLs0aTsL+XuPt/dLwLaEG7BeTzyuU2Ef7kvcwLh7mSbK9htEeEJBMq0ji7pICXvv5+yfWUf5HOHogHhOssC2ErCkzaU6UFUV7XKuPsed38RWEY4DGwAvgWOiTrHzd29rOveBuDW/b6DNHdfAOQAx5dtO9KKdTw1YGaXEP7ehrj78qhV/YEOwNdmlku4m+L3zaxswoFyx29mJwNHAp9XY7f3Eu5CGP29lwWzMYS7Ez5ak+M4iIq+q83unlfVh9x9N+GJE37NgWO9phBuab0K+FdkTFcgdYqIJCKFJBFJNI8AF5lZD8LdioaY2UAzSzGzVDPrb2btzayVmV0RGZv0LeG/gpdGtjEd+ImZnWTh6cT/B5gR6S62v6WEu+I1NLPehGeUK7M1ss2TK6n1VcLjRq41syPMbATQFXilkvdXW2RChIzIcTcnPN5pB7Aq8pZngP8ys3Zm1hb4f4QDZkXbutHMBptZMwtPanEp4fEnH7l7DvA68L9m1jyy/hQz+07k438Dxlt4NrmySR6uiqybB5xuZsMiXRl/TPmQebBj/C7hFpDvu/vC/VZPAk4h3D2xZ6SOeYS7DBL53BAz6xf5GfgVMMvdD9rq4e5ZwArCExWU1ZJK+OftHsJjyNqZ2Q9rcCyNItswoGHkZ7XsGvwMkGnhad1bEB7DNbmam54KpAKX7Ld8DtALGBfZftB1iogkFIUkEUkokYHyzwD3u/sGoGzigK2EWzXuIvz/tgaEB/RvItwF7DvAbZHNPEX4l8t3gX8Tnvjgjkp2eR/hX8Z3EP6r/b4xQu6+i/C4kA8i3c3O2a/WPOAywgElD7gbuKysy9jBmNlKi8zIV4EWhMNePuEJEU4BLnH3slnJ/g7MBZYT/oV/XmRZRb4hfA6/JjxT3O+B29y9bJa+6wlPKvEZ4fMwk3DrHO4+G/gd8LyZfRPZ16WRddsIt2T8NnL8nQh3Ays7vhMsPCtgZRM33Ed4vNWrkfcVmtn/Rba9KzILX25k1sNCoLhsIoXI+LMfEA5LWwiPm6l2qCEcAKLvq/QQsCEy4963hFtpHjSzTpFj+T8zu6eK7b1OeBzRuYQDXojIxB/u/hrhc/424e/gK/4z7q5K7l4C3L9frUTGuL0EnER4QgqCrFNEJNFYJd3TRUREJIGZ2f3Aqe4++qBvFhGRcmp7JicREREJmIWnn88kPLOdiIjUkLrbiYiI1CNmdgvhrqf/5+7vHuz9IiJyIHW3ExERERERiaKWJBERERERkSgKSSIiIiIiIlEUkkRERERERKIoJImIiIiIiERRSBIREREREYmikCQiIiIiIhJFIUlERERERCSKQpKIiIiIiEgUhSQREREREZEoCkkiIiIiIiJRjgi6gJo45phjvEOHDkGXcUiKiopo0qRJ0GUEJpmPP5mPHXT8iXT8n3zyyTZ3PzboOuJRIl9/KpJIP5exonNSMZ2XA+mcVCwRz0tV172ECkkdOnTg448/DrqMQ5KVlUX//v2DLiMwyXz8yXzsoONPpOM3s6+CriHemNkQYEjHjh0T9vpTkUT6uYwVnZOK6bwcSOekYol4Xqq67qm7nYiISCXcfa67j01PTw+6FBERiSGFJBERERERkSgKSSIiIiIiIlESakxSRfbs2cPGjRspLi4OupQqpaens2rVqqDLCEwsjj81NZX27dvTsGHDOt2PiIiIiNRvCR+SNm7cSLNmzejQoQNmFnQ5lSooKKBZs2ZBlxGYuj5+dycvL4+NGzdy0kkn1dl+RERERKT+S/judsXFxRx99NFxHZCk7pkZRx99dNy3KIqIiIhI/Ev4kAQoIAmgnwMRERERqR31IiRJbG3atInhw4cHXYaIiIiISJ1QSJIa2bt3L23btmXmzJlBlyIiIiIiUicUkg7T/fffzyOPPLLv9b333svEiRMP+rn8/Hw6d+7MmjVrABg5ciSPP/74Ae9btGgR5557Lj169KBPnz4UFBRQXFzMTTfdRPfu3cnIyODtt98GYPLkyQwdOpSLLrqIDh068NhjjzFhwgQyMjI455xz2L59OwD9+/dn3Lhx9OzZk27durFw4UIAFi5cSN++fcnIyODcc8/dV9vkyZO5/PLL+e53v8uFF17I+vXr6datGwArV66kT58+9OzZkzPOOIO1a9cCMGHCBLp160a3bt32nZ/169fTpUsXbrnlFk4//XQuvvhiQqHQIZx1EREREZG6k/Cz20X75dyVfLbpm1rdZte2zfnFkNMrXT9mzBiGDRvGnXfeSWlpKc8//zwLFy6koKCAfv367XtfaWkpDRqEM+lzzz1H165deeyxx7jxxhsZN24cO3bs4JZbbim37d27dzNixAhmzJjBWWedxTfffENaWhoTJ07EzFi+fDmrV6/m4osv5vPPPwdgxYoVLFmyhOLiYjp27Mjvfvc7lixZwk9+8hOeeeYZ7rzzTgB27drF0qVLeffddxkzZgwrVqzgtNNO47333uOII47gzTff5J577uGll14CYPHixSxbtoyWLVuyfv36fTX+7W9/Y9y4cYwaNYrdu3dTUlLCJ598wtNPP81HH32Eu3P22WfTu3dv2rdvz9q1a5k+fTqPP/44V199NS+99BKjR4+uja9KRKTWmdkQYEjHjh2DLkVERGKoXoWkIHTo0IGjjz6aJUuWsHnzZjIyMjj66KMBWLp06b73VTQF9kUXXcSLL77I7bffzqeffnrAttesWUObNm0466yzAGjevDkA77//PnfccQcAp512GieeeOK+kDRgwACaNWtGs2bNSE9PZ8iQIQB0796dZcuW7dv2yJEjAbjgggv45ptv2LlzJwUFBdxwww2sXbsWM2PPnj3lam3ZsuUBNfbt25ff/OY3bNy4kWHDhtGpUyfef/99rrzySpo0aQLAsGHDWLBgAVdffTUnnXQSPXv2BODMM88sF7hEpObmrZvHxMUTyS3KpXWT1ozrNY7BJw8Ouqx6w93nAnN79+59y0HfLCK1Ys6SbB6ev4ZNO0O0bZHGXQM7MzSjXdBlSZKpVyGpqhafunTzzTczefJkcnNzGTNmDEC1WpJKS0tZtWoVjRs3ZseOHbRv3/6waznyyCP3PW/QoMG+1w0aNGDv3r371u0/E5yZcd999zFgwABmz57N+vXr6d+//771ZYFnf9deey1nn3028+bNY9CgQfz973+vdn0pKSnqbidyGOatm8cDCx6guCQ89X1OUQ4PLHgAQEFJRBLSnCXZjJ+1nNCeEgCyd4YYP2s5gIKSxJTGJNWCK6+8ktdee41FixYxcOBAAJo1a8bSpUv3PT744IN9z7t27QrAH//4R7p06cJzzz3HTTfdVK7lBqBz587k5OSwaNEiIBy89u7dS79+/Zg2bRoAn3/+OV9//TWdO3euUc0zZswAwq1S6enppKenk5+fT7t24f8BTZ48uVrbWbduHSeffDI//vGPueKKK1i2bBn9+vVjzpw57Nq1i6KiImbPns25555bo/pE5OAmLp64LyCVKS4pZuLig4+LFBGJRw/PX7MvIJUJ7Snh4flrAqpIklW9akkKSqNGjRgwYAAtWrQgJSWlWp9Zs2YNTzzxBAsXLqRZs2ZccMEFPPjgg/zyl78st90ZM2Zwxx13EAqFSEtL48033+SHP/wht912G927d+eII45g8uTJ5VpoqiM1NZWMjAz27NnDU089BcDdd9/NDTfcwIMPPsjgwdX7K/QLL7zA1KlTadiwIa1bt+aee+6hZcuW3HjjjfTp0wcIt7T16NGDvLy8GtUoIlXLLcqt0XIRkXi3aWfFPUwqWy5SV8zdg66h2nr37u0ff/xxuWWrVq2iS5cuAVUUVlpaSq9evXjxxRfp1KlThe+paExSUPr3788f/vAHevfuHbN9xur44+HnYX9ZWVnlui4mGx1/3R3/96adx+a9B05W06ZJG14f/nqNt2dmn7h77P7HkEAquv4ksmT/77IiOicVi/V5Oe+3b5FdQSBq1yKND3723ZjVURX9rFQsEc9LVdc9dbc7TJ999hkdO3bkwgsvrDQgiYjUth0LpzNu079pVFp+eWpKKuN6jQumKJEamrMkm/N++xYn/Wwe5/32LeYsyQ66JAnYXQM7k9awfK+ctIYp3DWwZsMKRA6Xutsdpq5du7Ju3bqgy6iRrKysoEsQkcOQs3MXa+Y/Q9u9JzDm/B/ycvazmt1OEo4G6EtFyr57zW4nQVNIEhFJIDk7Chj5xMfk7/0RT1/fg9tPac/tXB90WSI1VtUAff1CnNyGZrTTz4AETt3tREQSxI6Fz1PwpwsoLdzGk5nn0vOUw79tgEhQNEBfROKZQpKISALYsfB5mr96GwUlR/Lo9X3pdcJRQZckcljatkir0XIRkVhSSBIRiXNlAWmpdybl+plqQZJ6QQP0RSSeKSTVgpSUFHr27EmPHj3o1asXCxYsOKTtPPLII+zatavCde+99x6nn346PXv2JDs7m+HDhwOwdOlSXn311UOuXUTi2/bFL0cC0qkKSAEwsyFmNik/Pz/oUuqdoRnteGhYd9q1SMMIT/H80LDuGosiInFBEzfUgrS0NJYuXQrA/PnzGT9+PO+8806Nt/PII48wevRoGjdufMC6adOmMX78eEaPHg3AzJkzgXBI+vjjjxk0aNChH4CIxKWc/BA/emMv13s/Trz+LwpIAXD3ucDc3r173xJ0LfWRBuiLSLxSS1It++abbzjqqP+MFXj44Yc566yz6Nu3L7/4xS8AKCoqYvDgwfTo0YNu3boxY8YMHn30UTZt2sSAAQMYMGBAuW0+8cQTvPDCC9x3332MGjWK9evX061bN3bv3s3999/PjBkz6NmzJzNmzIjpsYpI3dm26j1G//0DPt/VlOMzpyggiYiIxFD9a0l6uoL7g5w+FPrcArt3wbSrDlzf81rIGAVFefDCflPp3jTvoLsMhUL07NmT4uJicnJyeOuttwB4/fXXWbt2LQsXLuSbb75h1KhRvPvuu2zdupW2bdsyb1542/n5+aSnpzNhwgTefvttjjnmmHLbv/nmm3n//fe57LLLGD58OOvXrwegUaNG/OpXv+Ljjz/mscceO2idIpIYdix6gaPm3coVfjXnZ/5WkzSIiIjEmFqSakFZd7vVq1fz2muvcf311+PuvP7667z++utkZGTQr18/Vq9ezdq1a+nevTtvvPEGP/3pT3nvvfdIT08P+hBEJE7sWPQCzefdylI/lQuu/7kCkoiISADqX0tSVS0/jRpXvb7J0dVqOapK37592bZtG1u3bsXdGT9+PLfeeisFBQU0a9Zs3/sWL17Mq6++ys9//nMuvPBC7r///sPar4gkvrKA9Kl3ikzScHzQJYmIiCQltSTVstWrV1NSUsLRRx/NwIEDeeqppygsLAQgOzubLVu2sGnTJho3bszo0aO56667WLx4MQDNmjWjoKCgRvs7lM+ISPzJ3ZzLEfPu5FPvRIPrX1JAEhERCVD9a0kKQNmYJAB3Z8qUKaSkpHDxxRezatUq+vbtS2lpKc2bN+fZZ5/liy++4K677qJBgwY0bNiQv/71rwCMHTuWSy65hLZt2/L2229Xa98DBgzgt7/9LT179mT8+PGMGDGirg5TROpITn6Ikc+sorXfy8+uv1wBSUREJGAKSbWgpKSk0nXjxo1j3Lhx5brbnXLKKQwcOPCA995xxx3ccccdFW5n8uTJ+5536NCBFStWANCyZUsWLVp0GNWLSJB2fPwiM15fSN6332NC5rX01BgkERGRwKm7nYhIQHZ8/CLNXrmVC759hyk39dIkDSIiInFCLUkiIgHY8fFMmr1yK8u9Iw2um0nPDscGXZKIiIhEqCVJRCTGwgFpLMu9I3bdTHp2PCHokkRERCSKQpKISAzl5Id49vUFCkgiIiJxTN3tRERiJHfLZkZO+Yy8bwdy3k0/I6PDcUGXJCIiklDmLMnm4flr2LQzRNsWadw1sDNDM9rV+n7UkiQiEgM7Pp5J2l/OpFXhaqZk9qGXAlJCMLMhZjYpPz8/6FJERJLenCXZjJ+1nOydIRzI3hli/KzlzFmSXev7CiwkmVmqmS00s0/NbKWZ/TKoWmrDnDlzMDNWr15d6Xv69+/Pxx9/DMCgQYPYuXNnlds899xzD7rfRx55hF27dtWo1iBt2rSJ4cOHB12GSExt/+Qlmr0ylvXemvHXDdYsdgnE3ee6+9j09PSgSxERSXoPz19DaE/5W++E9pTw8Pw1tb6vIFuSvgW+6+49gJ7AJWZ2Tl3vdN66eVw882LOmHIGF8+8mHnr5tXKdqdPn87555/P9OnTq/X+V199lRYtWlT5ngULFhx0O4kUkvbu3Uvbtm2ZOXNm0KWIxEzqpgU0n3sLK/1k/LpZGoMkIiJyiDbtDNVo+eEILCR5WGHkZcPIw+tyn/PWzeOBBQ+QU5SD4+QU5fDAggcOOygVFhby/vvv8+STT/L888/vWx4Khbjmmmvo0qUL1157LaHQf77ADh06sG3bNgAmTJhAt27d6NatG4888si+9zRt2hSArKws+vfvz/DhwznttNMYNWoU7s6jjz7Kpk2bGDBgAAMGDDigrkWLFnHuuefSo0cP+vTpQ0FBAcXFxdx00010796djIwM3n77bSB8s9qhQ4dy0UUX0aFDBx577DEmTJhARkYG55xzDtu3bwfCrWHjxo2jZ8+edOvWjYULFwKwcOFC+vbtS0ZGBueeey5r1qzZt93LL7+cyy67jAsvvJD169fTrVs3AFauXEmfPn3o2bMnZ5xxBmvXrq30fKxfv54uXbpwyy23cPrpp3PxxReXO58i8Wjbqvfo/fkfFJBERERqQdsWaTVafjgCHZNkZilmthTYArzh7h/V5f4mLp5IcUlxuWXFJcVMXDzxsLb78ssvc8kll3Dqqady9NFH88knnwDw17/+lcaNG7Nq1SruueeefcujffLJJzz99NN89NFHfPjhhzz++OMsWbLkgPctWbKERx55hM8++4x169bxwQcf8OMf/5i2bdvy9ttv7ws7ZXbv3s2IESOYOHEin376KW+++SZpaWn8+c9/xsxYvnw506dP54YbbqC4OHxOVqxYwaxZs1i0aBH33nsvjRs3ZsmSJfTt25dnnnlm37Z37drF0qVL+ctf/sKYMWMAOO2003jvvfdYsmQJv/rVr7jnnnv2vX/x4sU888wzvPPOO+Vq/Nvf/sa4ceNYunQpH3/8Me3bt6/yfKxdu5bbb7+dlStX0qJFC1566aVD+bpEYiInP8SIucX8rfQKBSQREZFacNfAzqQ1TCm3LK1hCncN7Fzr+wp0djt3LwF6mlkLYLaZdXP3FdHvMbOxwFiAVq1akZWVVW4b6enpFBQUVGt/uUW5lS6v7jYqMnXqVG677TYKCgoYOnQoU6ZM4dRTT+Wtt97iBz/4AQUFBXTp0oVu3bpRVFREQUEB7k5hYSFvvvkmgwYNorS0FIDBgwfzxhtv0LFjRwAKCgrYtWsXZ555Junp6RQVFXH66aezatUqevTosW87Rx55ZLmaVq5cyXHHHcdpp51GQUEBZkYoFCIrK4tbb72VgoIC2rVrR/v27VmyZAnFxcWcf/75AKSmptK8eXMGDBhAQUEBnTp1YsWKFRQUFFBSUsIVV1xBQUEBGRkZ5Ofns2HDBgoLC7n77rv58ssvMTP27Nmzr+Wqf//++76nwsJCSktLKSgooGfPnjz44IN8+eWXDBkyhI4dO1Z6PgYNGsSJJ57IKaecQkFBAd26dWPNmjUHfG/FxcUH/IwErbCwMO5qiqVkO/5FhYt4Je8ltnsR3rIFn7e5mG4b15G1cV3QpYmIiCS0slnsYjG7XVxMAe7uO83sbeASYMV+6yYBkwB69+7t/fv3L/fZVatW0axZs2rtp3WT1uQU5VS4vLrb2N/27dt59913WbVqFWZGSUkJZsbEiRM54ogjaNy4Mc2aNaOgoIAGDRrQpEkTmjVrhpnRtGlTUlNTOfLII/ft/8gjjyQ1NXXf62bNmtG4ceN924FwiGnYsGG57exff5MmTUhJSTlgeXRNACkpKTRp0oTU1NRy20lJSeHoo4/et/8GDRrQrFmzfe8ve5+Z0bx5c+677z4uuugi5s6dy/r16+nfvz/NmjUjNTWVFi1a7KuladOm+7aVmZlJ//79mTdvHldffTV///vfKz0fTZs2JS0tbd/yxo0bU1hYeMDxpaamkpGRcUjfZV0p6y6ZrJLp+Oetm8eM96fzLXvAwBru5H2fzcATejH45MFBlyciIpLwhma0q5NQtL8gZ7c7NtKChJmlARcBlU8NVwvG9RpHakpquWWpKamM6zXukLc5c+ZMrrvuOr766ivWr1/Phg0bOOmkk3jvvfe44IILeO655wD47LPPWLZs2QGf79evH3PmzGHXrl0UFRUxe/Zs+vXrV+39lwWw/XXu3JmcnBwWLVoEhFuk9u7dS79+/Zg2bRoAn3/+OV9//TWdO9esiXLGjBkAvP/++6Snp5Oenk5+fj7t2oV/YCdPnlyt7axbt46TTz6ZH//4x1xxxRUsW7bssM+HSJD++OFDfOt7yi3b43sOu0uviIiIxFaQLUltgClmlkI4rL3g7q/U5Q7L/pI7cfFEcotyad2kNeN6jTusv/BOnz6dn/70p+WWff/732f69OlMmDCBm266iS5dutCpUyfOPPPMcu8zM3r16sWNN95Inz59ALj55ptr1BIyduxYLrnkkn1jk8o0atSIGTNmcMcddxAKhUhLS+PNN9/khz/8Ibfddhvdu3fniCOOYPLkyQd01TuYstaaPXv28NRTTwFw9913c8MNN/Dggw8yeHD1zucLL7zA1KlTadiwIa1bt+aee+6hZcuWFZ6P9evX16hGkVjbvng2W3bvBLMD1lXW1VdERETik7nX6YRytap3795edp+hMqtWraJLly4BVVR9BQUF+7qGlZSUcNxxx5Gbm0vDhg0Drqxm+vfvzx/+8Ad69+5do89FH39disefh2TqblaRZDj+nB2FFP3pPG5uW0peBf9Jt2nShteHvx77wmrIzD5x95r9x50kKrr+JLJk+O+ypnROKqbzciCdk4ol4nmp6roX6Ox2yer000/n5ptvTriAJCIHyskPMfKJRdxccg9Xdf1/B3TpbWgND6tLr4iIiMReXEzckGxWr67ToVd1KplmKRM5mO2L57Bi3pPk7/0BT2ZeTK8TjqJD62PLdem9KPUiTdogIiKSYBSSREQOwfbFc2j2j0xaeQeevq47PU84CgiPfYwORfrDgoiISOKpF93tEmlcldQd/RxIrJQFpFXeAR/9km4UKyIiUs8kfEhKTU0lLy9PvyAnOXcnLy+P1NTUg79Z5DBsX/xyuYDUo1OHoEsSERGRWpbw3e3at2/Pxo0b2bp1a9ClVKm4uDipf4GPxfGnpqbSvn37Ot2HJLec/BC/fn0T13kXGo9+VgFJRESknkr4kNSwYUNOOumkoMs4qKysrBrd/6i+Sfbjl8S35evVjJyRQ17oZG4e8w96nNgy6JIkBsxsCDCkY8eOQZciIiIxlPDd7URE6tr2JS/T4qnz6Fv4JlMy+9BLASlpuPtcdx+bnp4edCkiIhJDCd+SJCJSl7YveZlmL49hjZ/INaNvoUdkFjsRERGpv9SSJCJSieiAVDJ6lsYgiYiIJAmFJBGRCmzesJamCkgiIiJJSSFJRGQ/OfkhRjy/kXv8RwpIIiIiSUhjkkREomxf8g9+O/8r8kKncm3mnRqDJCIikoQUkkREIrYv+QdNX76J67wTN4x5lV4KSCIiIklJIUlEhP8EpLV+Ag1HP6/7IImIiCQxjUkSkaQXHZD2jp6tMUgiIiJJTi1JIpLUcvJDLJ43mZP8eAUkERERARSSRCSJ5ewoZOQTi9ix5xamXNeVnh1PDLokERERiQPqbiciSWn70lcIPdqXlMJcns48RwFJRERE9lFIEpGks33pKzSdcwPFpSlMGH2OZrETERGRchSSRCSplAWkL/x49oyeTY9OJwVdkoiIiMQZjUkSkaSRt+INmikgiUgNzFmSzcPz17BpZ4i2LdK4a2Bnhma0C7osEaljakkSkaSQkx9izKu7eM3PVUASkWqZsySb8bOWk70zhAPZO0OMn7WcOUuygy5NROqYQpKI1Hvb1nzAdX9/j3VFqbTPfEYBSUSq5eH5awjtKSm3LLSnhIfnrwmoIhGJFXW3E5F6bfun82g++wZG+SX0yPyTJmkQkWrbtDNUo+UiUn+oJUlE6q28T1+l6ewb+MLb02v0rxWQRKRG2rZIq9FyEak/FJJEpF7K+/RVms2+ni+9ncYgySEzsyFmNik/Pz/oUiQAdw3sTFrDlHLL0hqmcNfAzgFVJCKxopAkIvVOzrY8mPNDvvR27B6lgCSHzt3nuvvY9PT0oEuRAAzNaMdDw7rTrkUaBrRrkcZDw7prdjuRJKAxSSJSr+Tkhxj59DKOKR3PfaMuosepJwddkogksKEZ7RSKRJKQQpKI1Bt5y/6Pl+e+St7uy5iQeTU9NAZJREREDoFCkojUC3nL/o9ms66jv7fl7BvuIUMBSURERA6RxiSJSMIrC0jrvC3fjppNxsltgi5JREREEphCkogktLxl/0fTWdfvC0g9Tj0l6JJEREQkwSkkiUjCyskP8cQr7/Olt1NAEhERkVqjMUkikpByt2xh5JSV5H17ARfdNI5eHY4LuiQREZGkMWdJNg/PX8OmnSHatkhj8Akl9A+6qFqkliQRSTh5y16j8V8y6FC4lCmZfRSQREREYmjOkmzGz1pO9s4QDmTvDDF5xW7mLMkOurRao5AkIgklb9l8ms66jlxvyU9GXU4vzWInIiISUw/PX0NoT0m5ZbtLw8vrC4UkEUkY4YA0mq+8NaFRczQGSUREJACbdoZqtDwRKSSJSELY+uViBSQREZE40LZFWo2WJyKFJBGJezn5Ia5+aQeP+1AFJBERkYDdNbAzaQ1Tyi1r1CC8vL7Q7HYiEtfyVv6TO1/ZxrZdLTg38/f00BgkERGRQA3NaAdwwOx2ZcvrA4UkEYlbectfp8lLo/iBn0F65kxN0iAiIhInhma0KxeKsrKygiumDigkiUhcKgtIX3trjr7275yhgCQiIiIxojFJIhJ3ogNS8bWzOaNzx6BLEhERkSSiliQRiSs5O3exZfYvaOytFJBEREQkEApJIhI3cvJDjHz8I/aU/Dd/vbaHApKIiIgEQt3tRCQu5K14g88fvZKCwkL+lPldBSQREREJjFqSRCRweSveoMnMa2nrrXhqVFdN8y0iIiKBUkuSiASqLCBt8FaErp2tG8WKiIhI4BSSRCQweSveLBeQzujcKeiSRERERBSSRCQYOfkh7n7la5Z5JwUkERERiSsakyQiMbfl6zWMnLGJvFA7bh/zCmec2DLokkRERET2UUuSiMRU3oo3afZUPy4q/AdTMvvQSwFJ4piZDTGzSfn5+UGXIiIiMaSQJCIxk7fiTRrPvJZsP5Yh1/6QXprFTuKcu89197Hp6elBlyIiIjGkkCQiMREdkHZdO0djkERERCRuBRaSzOx4M3vbzD4zs5VmNi6oWkSkbuXmbiR15uhIQNIkDSIiIhLfgmxJ2gv8P3fvCpwD3G5mXQOsR0TqwPbiUq6ZupZ7/PZIQDo16JJEREREqhRYSHL3HHdfHHleAKwC2gVVj4jUvryVb/HeR/8ir3A3N2T+SAFJREREEkJcTAFuZh2ADOCjgEsRkVqSt/ItGr94Dbd4G6676VZN0iAiIiIJI/CQZGZNgZeAO939mwrWjwXGArRq1YqsrKzYFlhLCgsLE7b22pDMx5+Mx25bltPns1+T7cew4JS76bh+BVnrg64qGMn4/YuIiCS6QEOSmTUkHJCmufusit7j7pOASQC9e/f2/v37x67AWpSVlUWi1l4bkvn4k+3Y81a+ReOsX7PJj2HXtXPomLMpqY5/f8n2/YuIiNQHQc5uZ8CTwCp3nxBUHSJSe3LyQ/xz9tNs8mMoGqlJGkRERCQxBTm73XnAdcB3zWxp5DEowHpE5DDk7Chk5KQPeXDPtRSOfpUzTuscdEkiIiIihySw7nbu/j5gQe1fRGpP3sq3+HbmHaSV/pTJmUPooUkaREREJIEF2ZIkIvXAtpVv0/jFkZSWOr8bebZmsRMREZGEp5AkIods28q3afLiNeT40RSOnK0udiIiIlIvKCSJyCHZtmaBApKIiIjUSwpJIlJjOfkhbnh5O2/6WQpIIiIiUu8EfjNZEUksW9cu5MbZ29lU1JB2mVM5Q2OQREREpJ5RSBKRatv2WRZNXxjBLX4eJ2c+qUkaREREpF5SSBKRatn2WRZNXhhBrh/FqSMfUguSiIiI1FsakyQiB7XtsywaRwJSwciXNQZJRERE6jW1JIlIlXK2f0Ppi2MpUEASERGRJKGWJBGpVE5+iJFPfsIPS3+qgCQiIiJJQy1JIlKhbZ+9w/zZ08jbPYwJmcM0BklERESShkKSiBxg22fv0PiFq+nvR9HzunvpqYAkIiIiSUTd7USknLKAtMWPomDkHHp2PDHokkRERERiSi1JIsK8dfOYuHgiuUU5tNpTwqjGLTnz8tl0P+20oEsTERERiTm1JIkkuXnr5vHAggfIKcrBgdyGKfypTWO+bvRl0KWJiIiIBEIhSSTJTVw8keKS4nLLdpfuZuLiiQFVJCIiIhIshSSRJJdblFPJ8twYVyIiIiISHxSSRJLYtlXv0WpPSYXrWjdpHeNqJBmZWRMzSwm6DhERkWgKSSJJatuq90ibcRWjtjuNGjQqty41JZVxvcYFVJnUZ2bWwMyuNbN5ZrYFWA3kmNlnZvawmXUMukYRERGFJJEktOWr1aTNuIptns6ZQ2bxq/N+RZsmbTCMNk3a8MC5DzD45MFBlyn109vAKcB4oLW7H+/uxwHnAx8CvzOz0UEWKCIictApwM3sd+7+04MtE5HEkJMfYuSMTQz1oVx4zTi6d+lCd7ooFEmsfM/d9+y/0N23Ay8BL5lZw9iXJSIi8h/VaUm6qIJll9Z2ISJS97at/oC7/vYSeUV7uCDzIbp36RJ0SZJk9g9IZpZqZjeb2R1mdnRF7xEREYm1SluSzOw24IfAyWa2LGpVM+CDui5MRGrXttXvk/r8VdzpJ9JgzKv0OuGooEsSAZhI+JpSDMwB+gVajYiICFV3t3sO+D/gIeBnUcsLIt0iRCRBlAWk7d6MI695iu4ntgy6JElSZjYd+Lm7l92tuCXwYuT5zyr+lIiISGxVGpLcPR/IB0ZGpmdtFXl/UzNr6u5fx6hGETkM0QEp/5o5dO/SNeiSJLndCzxoZjnAr4E/ALOBVOCBAOsSERHZpzoTN/yI8IVrM1AaWezAGXVXlojUhpz8EOtf/AXtFJAkTrj7OuBaMzsfmAHMAwa7e8U37AqYmQ0BhnTsqJnJRUSSSXUmbrgT6Ozup7t798hDAUkkzuXkhxg56UPu3HsH+de8rIAkccHMjjKz24GuwFXADmB+JIzEHXef6+5j09PTgy5FRERiqDohaQPhbncikiC2rfmAdY8OYVfhN/w18zuaxU7iyRxgJ+EeCVPdfSowBMgws7kB1iUiIrLPQbvbAeuALDObB3xbttDdJ9RZVSJySOatm8cfP/o9W77N49g2cFPnzZrFTuLN0cBMIA24FcDdQ8CvzKxNkIWJiIiUqU5I+jryaBR5iEgcmrduHr/44H6+Ld0NZmxpCE999TgntGunG8VKPPkF8BpQwn6z2bl7TiAViYiI7OegIcndfwlgZo3dfVfdlyQih+KPH/0+HJCiFJcUM3HxRIUkiRvu/hLwUtB1iIiIVOWgY5LMrK+ZfQasjrzuYWZ/qfPKRKTacvJDbP624tuX5RblxrgakcqZ2eNm1q2SdU3MbIyZjYp1XSIiItGq093uEWAg8A8Ad//UzC6oy6JEpPq2bPickdM3QssW0HDnAetbN2kd85pEqvBn4H4z6w6sALYSvkdSJ6A58BQwLbjyat+cJdk8PH8Nm3aGaNsijbsGdmZoRrugy4opnQMRSTTVCUm4+wYzi14Ul/ezEEk229YsoPH073OlX05K/zuY8vnDFJcU71ufmpLKuF7jAqxQpDx3XwpcbWZNgd5AGyAErHL3NUHWVhfmLMlm/KzlhPaEL5vZO0OMn7UcIGlCgs6BiCSi6oSkDWZ2LuBm1hAYB6yq27JE5GC2rVlA6vTvs9Ob8t0Rd9C9azdOOqYJExdPJLcol9ZNWjOu1ziNR5K45O6FQFbQddS1h+ev2RcOyoT2lPDw/DVJExB0DkQkEVUnJP0AmAi0A7KB14Hb67IoEaladEDaMWI23buGh3gMPnmwQpFIHNm0M1Sj5fWRzoGIJKLqzG63DdAgWpE4kbt1G0dOv4ad3pSdI2btC0giEn/atkgju4Iw0LZFWgDVBEPnQEQSUaWz25nZ3ZF//2Rmj+7/iF2JIlImJz/ENZOXc6/fxs4Rs+jWtXvQJYnUmJmlmNkfgq4jFu4a2Jm0hinllqU1TOGugZ0Dqij2dA5EJBFV1ZJUNu7o41gUIiJV27bmX/x11lvkFffm5swf0O2Eo4IuSeSQuHuJmZ0fdB2xUDbmJplndtM5EJFEVGlIcve5kX+nxK4cEanItjX/4sjp3yfTm3HljTeToYAkiW+Jmf0DeBEoKlvo7rOCK6luDM1ol/SBQOdARBJNpSHJzOYCXtl6d7+8TioSkXLKAlK+N6ZgxEwyTmoVdEkitSEVyAO+G7XMgXoXkkREJPFU1d2urL/4MKA18Gzk9Uhgc10WJSJh0QFp54jZGoMk9Ya73xR0DSIiIpWpqrvdOwBm9r/u3jtq1Vwz0zglkTqWkx/itZmTucgbs/NqTdIg9YuZtQf+BJwXWfQeMM7dNwZXlYiISFils9tFaWJmJ5e9MLOTgCZ1V5KI5OwoZOSkD5mw+0ryRr9Bt9PPCLokkdr2NPAPoG3kMTeyTEREJHDVCUk/AbLMLMvM3gHeBu6s06pEktjWzz9i96Nnk164jimZZ9Oj00lBlyRSF45196fdfW/kMRk4NuiiREREoHo3k33NzDoBp0UWrXb3b+u2LJHktPXzjzjyuWHs9jR+c/VZmuZb6rM8MxsNTI+8Hkl4IgcREZHAHTQkRXQCOhOejaiHmeHuz9RdWSLJpywgFXgaO66epS52Ut+NITwm6Y+EZ7VbAGgyBxERiQsHDUlm9gugP9AVeBW4FHgfUEgSqSVb1y1VQJKkYWYpwP/oVhIiIhKvqjMmaThwIZAbmbK1B5Bep1WJJJGc/BCjZ+bytvdSQJKk4O4lwIlm1ijoWkRERCpSne52IXcvNbO9ZtYc2AIcX8d1iSSFrV8uJvOlHDYVHcHxmc9oDJIkk3XAB2b2D6CobKG7TwiuJBERkbDqhKSPzawF8DjwCVAI/KsuixJJBlvXLuTIaVdyh59Bq8zn6KWAJMnly8ijAdAs4FpERETKqTIkmZkBD7n7TuBvZvYa0Nzdl8WiOJH6qiwgFXoqx1/9O7UgSVKJjEk61d1HBV2LiIhIRaock+TuTniyhrLX6xWQRA7P1rULaTRtGIWeyvarZ2sMkiQdjUkSEZF4V53udovN7Cx3X1Tn1YjUczk7iyiafguN/Ui2X6VJGiSpaUySiIjEreqEpLOBUWb2FeELmRFuZDrs3+7M7CngMmCLu3c73O2JxLOc/BAjH19IWsl/8fBVPenWrUfQJYkESWOSREQkblUnJA2sw/1PBh5D91ySem7r2kW8/cJfyNs9gimZl2sMkiQ9d//l/svMrLo3OBcREalTB71Pkrt/5e5fAXsJ3xXdgeza2Lm7vwtsr41ticSrkrwvaTTtSgbsfodp156iWewkqZnZ+1HPp+63emGMyxEREalQpX+1M7PxQEN3/1Vk0b+AnUAjYArwUJ1XJ5Lgtq5dRO/l97PLjyTvqtmc0blT0CWJBK1J1PP9u1lbLAsRERGpTFVdG64C+kW9znP3jMjUre8Qo5BkZmOBsQCtWrUiKysrFrutdYWFhQlbe21IxuMvyfuS3svvp8iP5L2uv6bVth1Jdw4gOb/7aMl+/BXwSp5X9FpERCQQVfb/dveiqJcTI8tKzCytTqsqX8MkYBJA7969vX///rHada3KysoiUWuvDcl2/Dn5IR75yyI6elM+7Ho/I0Yk7+1gku2731+yH38FWpjZlYS7e7cws2GR5QakB1eWiIjIf1QVkpqaWUN33wPg7pMBzOxIoHkMahNJSLlbtzJy8gryis9gxE3v0Wr9Z0GXJBJP3gEuj3o+JGrdu7EvR0RE5EBVhaSZwN/N7EfuvgvAzJoQno1uZm3s3MymA/2BY8xsI/ALd3+yNrYtEoStX3xCo2nD6FU6htGZP6bXCUeRpZAkso+73xR0DSIiIgdTVUi6D/gN8HXkHkkAJwBPRtYdNncfWRvbEYkHW7/4hEbPXkHIGzHmqis0zbeIiIhIgqo0JLl7CfAzM/sl0DGy+At3D8WkMpEEEh2Qtl01i27degZdkoiIiIgcooPeuC8SipbHoBaRhLR501ekKiCJiIiI1BsHvZmsiFQuJz/EiGe/ZJIPVUASqQEza2xm95nZ45HXnczssqDrEhERgWq0JIlIxbZ+uZh7Z35K3q62XJj5a41BEqmZp4FPgL6R19nAi8ArgVUkIiIScdCWJDP7Z3WWiSSTrV8upuHUK7hr1yNMGXMWvRSQRGrqFHf/PVB2m4ldhO+VJCIiErhKW5LMLBVoTHh67qP4z8WrOdAuBrWJxKWygPStH0HpVZPpdWLLoEsSSUS7IzcmdwAzOwX4NtiSREREwqrqbncrcCfQlnCXiLKQ9A3heyWJJJ3ogLTtqlmc3i0j6JJEEtUDwGvA8WY2DTgPuDHIgkRERMpUNQX4RGCimd3h7n+KYU0icSknP8SK5x7gDAUkkcPm7q+b2SfAOYT/CDfO3bcFXJaIiAhQvSnA/2Rm5wIdot/v7s/UYV0icSUnP8TISR9SsPcWpozoQLeu3YMuSSShmdlc4DngH+5eFHQ9IiIi0aozccNU4A/A+cBZkUfvOq5LJG5s/XIx6x8dzO7CHTyeeb4Ckkjt+APQD/jMzGaa2fDIWFgREZHAVWcK8N5AV3f3ui5GJN6UjUE62VOYdNXJmuZbpJa4+zvAO2aWAnwXuAV4ivDkQCIiIoGqTkhaAbQGcuq4FpG4svXLxRwxdSjfegrbhutGsSK1LTK73RBgBNALmBJsRSIiImHVCUnHEO4OsZCo6Vnd/fI6q0okYFu/XMIRU4ey2xuwbfgsTu/eK+iSROoVM3sB6EN4hrvHgHfcvTTYqkRERMKqE5IeqOsiROJJTn6In8xcw93emiOH/00BSaRuPAmMdPeSoAsRERHZX3Vmt3snFoWIxIPNG7/g2ue+Ytuuo2DMa5yuG8WK1Coz+667vwU0Aa4ws3Lr3X1WIIWJiIhEqTQkmdn77n6+mRUQuSN62SrA3V2Da6Ve2frlEhpNvYLR3p+MzIn00iQNInXhO8BbhMci7c8BhSQREQlcVTeTPT/yb7PYlSMSjPAYpCvY4w04Z/g4TldAEqkT7v6LyNNfufu/o9eZ2UkBlCQiInKAg94nCcDMepjZjyKPM+q6KJFYig5IW4e/xOndzwy6JJFk8FIFy2bGvAoREZEKHHRMkpmNI3z/irIuENPMbJK7/6lOKxOJgZy8ndizV+HegK3DZyogidQxMzsNOB1IN7NhUauaA7qZrIiIxIXqzG6XCZzt7kUAZvY74F+AQpIktJz8ECOfWkKn0lu5c/gATu/eO+iSRJJBZ+AyoAXlxyUVEP6DnIiISOCqE5IMiJ6itSSyTCRhbV23lCdnzCGvuC8TMm/WGCSRGHH3l4GXzayvu/8rqDrM7GTgXiDd3YcHVYeIiMSn6oSkp4GPzGx25PVQwve3EElIW9ctJeWZy7nFU7js+pvpqYAkEoQlZnY74a53+7rZufuYg33QzJ4i3Bq1xd27RS2/BJgIpABPuPtvK9uGu68DMs1M46BEROQAB524wd0nADcB2yOPm9z9kTquS6ROlAWkEje2DX+Rnqe0D7okkWQ1FWgNDATeAdoT7nJXHZOBS6IXmFkK8GfgUqArMNLMuppZdzN7Zb/HcbV1ECIiUj9VdZ+k5u7+jZm1BNZHHmXrjgK+0Z3SJZFEB6Qtw1/SGCSRYHV096vM7Ap3n2JmzwHvVeeD7v6umXXYb3Ef4ItICxFm9jxwhbs/RLjVSUREpNqq6m73HOELyycceDNZgKZm9ri731NXxYnUlpz8EC89P4URbmz5/kwFJJHg7Yn8u9PMugG5wOG08LQDNkS93gicXdmbzexo4DdAhpmNj4Sp/d8zFhgL0KpVK7Kysg6jvPhSWFhYr46nNuicVEzn5UA6JxWrb+elqpvJXhb5t8Kb+0W6NqwAFJIkruXsKGTkE4vI+/YS+o26gx6nnhx0SSICkyK9Eu4D/gE0Be6P1c7dPQ/4wUHeMwmYBNC7d2/v379/DCqLjaysLOrT8dQGnZOK6bwcSOekYvXtvFTV3a5XVR9098VAl1qvSKQWbV33Kd9OvZbjSn/EhMxr6KFJGkTigrs/EXn6DlAbf7nIBo6Pet0+skxERKTGqupu97+Rf1OB3sCnhLvanQF8DPSt29JEDs/WdZ+SMvVympY6D3y/N10VkEQCZ2b/VdX6yGRBh2IR0MnMTiIcjq4Brj3EbYmISJKrdHY7dx/g7gOAHKCXu/d29zOBDPTXOYlzW/69jJSpl1Na6mz5/kt0PeOsoEsSkbBmB3kclJlNJ3xT885mttHMMt19L/AjYD6wCnjB3VfWQf0iIpIEqnOfpM7uvrzshbuvMDN1s5O4teXr1RzxzBAFJJE45O6/rIVtjKxk+avAq4e7fRERkYPeJwlYZmZPmFn/yONxYFldFyZyKHLyQ4x6/mveKe2pgCQSx8zsVDP7p5mtiLw+w8x+HnRdIiIiUL2QdBOwEhgXeawEbqzDmkQOyZb1K7j176+TW+ScmDlZAUkkvj0OjCcyFbi7LyM8jkhERCRwB+1u5+7FwB8jD8ysHzABuL1uSxOpvi3/XkbKM5dzd+kJNM58mV6apEEk3jV294VmFr1sb1DFiIiIRKvOmCTMLAMYCVwN/BuYVZdFidREWUCitJSWw/6gWexEEsM2MzuFyM3KzWw44YmCREREAlfVfZJOJRyMRgLbgBmARWa8E4kLW/69fF9A2jxsJl179Am6JBGpntsJ36j1NDPLJvwHuFHBlnQgMxsCDOnYsWPQpYiISAxVNSZpNfBd4DJ3P9/d/wSUxKYskYPL2bmL3Km3KCCJJCB3X+fu3wOOBU4DvgOcH2xVB3L3ue4+Nj09PehSREQkhqoKScMId31428weN7MLCd9MViRwOfkhRj7+Ef9V8mMFJJEEYmbNzWy8mT1mZhcBu4AbgC8Id+kWEREJXKXd7dx9DjDHzJoAVwB3AseZ2V+B2e7+ekwqFNnPln8v571pv2fHnmt5OvNSjUESSSxTgR2EbwZ7C3Av4T/AXenuSwOsS0REZJ/qzG5XBDwHPGdmRwFXAT8FFJIk5rb8ezkNnhnChaUldB1xN90UkEQSzcnu3h3AzJ4g3GPhhMhMqiIiInGhOvdJ2sfdd7j7JHe/sK4KEqlMWUCy0hI2D5tJt67dgy5JRGpuT9kTdy8BNiogiYhIvKnWFOAiQduyfkW5gNS1x9lBlyQih6aHmX0TeW5AWuS1Ae7uzYMrTUREJEwhSeJeTn6I3z7/NneVHsE3w2YoIIkkMHdPCboGERGRg1FIkriWu3UbIycvJy90Kjfe+AEZJ7UKuiQRERERqedqNCZJJJa2rF9Jg7/0oV/ha0zJ7KOAJCIiIiIxoZYkiUtb1q/EplxGSukeRg67UtN8i4iIiEjMqCVJ4k50QModNpOuPc4JuiQRSVJmNsTMJuXn5wddioiIxJBCksSV3M25CkgiEjfcfa67j01PTw+6FBERiSGFJIkbOfkhrnlmFY+XDlVAEhEREZHAaEySxIUtX33GL6e/R16oA5dk3q8xSCIiIiISGIUkCcy8dfOYuHgiuUW5HLu3hOsblDD2pvfopYAkIiIiIgFSdzsJxLx183hgwQPkFOXgOFuOaMCjrZuQU7ow6NJEREREJMkpJEkgJi6eSHFJcbllu303ExdPDKgiEREREZEwhSQJRG5Rbo2Wi4iIiIjEikKSxFxOfggraVHhutZNWse2GBERERGR/SgkSUxt+WoVXz86iLRt36FRgyPLrUtNSWVcr3EBVSYiIiIiEhZoSDKzS8xsjZl9YWY/C7IWqXtbvloFkwdz6t61PHXxIH513i9p06QNhtGmSRseOPcBBp88OOgyRURERCTJBTYFuJmlAH8GLgI2AovM7B/u/llQNUnd+XZnNkweQ8PS3eQOnUHXnn3pCgpFIiIiIhJ3gmxJ6gN84e7r3H038DxwRYD1SB3Z8tVqei79+b6A1CXjvKBLEhGpFjMbYmaT8vPzgy5FRERiKMiQ1A7YEPV6Y2SZ1CM5+SF+8PxKNvoxCkgiknDcfa67j01PTw+6FBERiaHAuttVl5mNBcYCtGrViqysrGALOkSFhYUJW/uh2vXNNh5cciQ79qSyoMv9dM/fw+YkOweQnN99NB1/ch+/iIhIIgoyJGUDx0e9bh9ZVo67TwImAfTu3dv79+8fk+JqW1ZWFola+6HY8tVqfHImd5b2pOPYJ/hm3adJdfzRku2735+OP7mPX0REJBEF2d1uEdDJzE4ys0bANcA/AqxHakk4IA3myNJiug/9Cb1OOCrokkREREREqi2wliR332tmPwLmAynAU+6+Mqh6pHZEB6ScoTPoknF+0CWJiIiIiNRIoGOS3P1V4NUga5Dak7OjkOLJwzmqtJicoc8rIImIiIhIQor7iRskMeTkhxj5xCKOL72Fe4b2oktGv6BLEhERERE5JApJcti2fL2aqc8+Q96332FC5vV00RgkERE5BHOWZPPw/DVs2hmibYs07hrYmaEZujuIiMSeQpIcli1fr6b06csYWxpi4Kib6aGAJCIih2DOkmzGz1pOaE8JANk7Q4yftRxAQUlEYi7I2e0kwZUFpLTSXeQMfZ4ep54cdEkiIpKgHp6/Zl9AKhPaU8LD89cEVJGIJDOFJDkk0QFp09AZGoMkIiKHZdPOUI2Wi4jUJYUkqbGc/BBPPzuV1NIQm65QQBIRkcPXtkVajZaLiNQlhSSpkZwdhYyc9CHPfnsBX416ly69FJBEROTw3TWwM2kNU8otS2uYwl0DOwdUkYgkM03cINW25es1hJ4ezvGlNzMh83pN0iAi9Z6ZDQGGdOzYMehS6r2yyRk0u50cDs2QKLVFIUmqZcvXayh9ejBHl+7init6aZpvEUkK7j4XmNu7d+9bgq4lGQzNaKdfaOWQaYZEqU3qbicHteXrNZSUTdJwxQx1sRMREZG4oxkSpTapJUmqtDl7PaVPX0bj0iIFJBEREYlbmiFRapNakqRSOfkhrp32Be+XnqGAJCIiInFNMyRKbVJIkgpt2fA5t//tVbYUlXBK5hMKSCIiIhLXNEOi1CZ1t5MDbNnwOSVPDeYXpemUjHmdXpqkQUREROKcZkiU2qSQJOWUBaQmpUUcecUTnHZiy6BLEhEREakWzZAotUXd7WSfLRs+Z+9Tl9GktIhNV0zntF7fCbokEREREZGYU0gSIDxJw9rJP6RpaaECkoiIiIgkNXW3E3LyQ4yc9CG+dyx/G9qGLhmapEFEREREkpdakpLclg2fs+BPmeQX7uKPmRcpIImIiIhI0lNLUhKat24eExdPJLcol2P3lnJbowKevfwnnK5Z7ERERERE1JKUbOatm8cDCx4gpygHx9lyhPFQ66NZ33RL0KWJiIiIiMQFhaQkM3HxRIpLisst2+17mLh4YkAViYiIiIjEF3W3SwLR3escr/A9uUW5Ma5KRERERCQ+KSTVc2Xd6/ZvPdpf6yatY1SRiEjiMLMhwJCOHTsGXYqIiMSQutvVcxV1r9tfakoq43qNi1FFIiKJw93nuvvY9PT0oEsREZEYUktSPVdVNzrDaN2kNeN6jWPwyYNjWJWIiIiISPxSSKrnjk09hi3FWw9Y3qZJG14f/noAFYmIiIiIxDd1t6vHcvK2c+OGHI4sLT9Zg7rXiYiIiIhUTi1J9Uz0THZW0oJepWdy6wln8OKON8ktylX3OhERERGRg1BIqkf2n8nOU3awrM0uhp80ktcvvDvg6kREREREEoO629UT89bN45737znwRrGl3+pGsSIiIiIiNaCQVA+UtSCVemmF63WjWBERERGR6lNIqgcOdi8k3ShWRERERKT6FJLqgapaijSTnYiIiIhIzSgkJbic/BBW0qLCdQ2sAQ+c+4BmshMRERERqQGFpAS2OXsdXz06iKbb+tGowZHl1qWmpPI/5/+PApKIiIiISA0pJCWozdnr2PPEILrtXcWTF36PX533S9o0aYNhtGnSRi1IIiIiIiKHSPdJSkBlAalF6U6yhzzLaWcO4DRQKBIRERERqQVqSUow+wekzr2/F3RJIiIiIiL1ikJSAsnJD3HLtOVsLm1B9mUKSCIiIiIidUHd7RLE5k1fc/2zq8ktSsXGvErnE1sGXZKIiIiISL2kkJQANmf/m91PDOInpSfROnMavU44KuiSRERERETqLYWkOFcWkFqWbueUy/5EZwUkEREREZE6pTFJcSw6IG28bBqdz9IYJBGRWDKzIWY2KT8/P+hSREQkhhSS4lTOzl3kPXl1JCA9q4AkIhIAd5/r7mPT09ODLkVERGJI3e3iUE5+iJGPf8RxJWP49WUd6XzWRUGXJCIiIiKSNBSS4szm7PVMn/IX8r69kAmZ12gMkoiIiIhIjCkkxZHN2ev59olLubV0OxeNuIHuCkgiIiIiIjGnMUlxoiwgHV26nY2Dp9K96+lBlyQiIiIikpQUkuJAOCAN4pjSPDYOnkrnPhcHXZKIiIiISNJSSApYTn6IvzzzLC1Lt7Nh8LMKSCIiIiIiAdOYpADl7Cxi5OMLySs+i2GjPqDHqacEXZKIiIiISNJTS1JANmevp/DR8zi1cBFTMvsoIImIiIiIxAm1JAVg86avKH5iMO1Kt/Lfg8/gVM1iJyIiIiISNwJpSTKzq8xspZmVmlnvIGoIyuZNX1H8+CCOLd3KxsHPcGqfS4IuSUREREREogTV3W4FMAx4N6D9ByJ38yYFJBERERGROBdISHL3Ve6+Joh9B2V7cSkjp6zivdLuCkgiIiIiInFMY5JiYPOmr3jmow1sK21N18y/agySiIiIiEgcq7OQZGZvAq0rWHWvu79cg+2MBcYCtGrViqysrNopMEaKvsmj6+L7+L0b7/X8I9+s+5SsdUFXFXuFhYUJ993VlmQ+dtDxJ/vxi4iIJKI6C0nu/r1a2s4kYBJA7969vX///rWx2ZjYvOlrQo/fznG+jbc6/ZzMK5P3RrFZWVkk0ndXm5L52EHHn+zHLyIikojU3a6OhAPSpRxXupWNg6bQNJQWdEkiIiIiIlINQU0BfqWZbQT6AvPMbH4QddSVnPwQnz75o30B6dSzLw26JBERERERqaZAWpLcfTYwO4h917Wc/BAjJ33I7r038tSQcZzW+8KgSxIRERERkRoI6j5J9dLmnK/58E9jKCgs5LHMAQpIIiIiIiIJSGOSasnmnK/ZNWkQl5Ru5rQrbqOLpvkWEREREUlIakmqBWUBqXXpZjYMmkKXXhcEXZKIiIiIiBwihaTDtH9AOvXsQUGXJCIiIiIih0Eh6TDk5If42TNvcWRpSAFJRERERKSe0JikQ5S7bRsjn1pG3q42bL5xARkntQq6JBERERERqQVqSToEm3O+ZtefBzC8aDpTMvsoIImI1FNmNsTMJuXn5wddioiIxJBCUg1tzvmaokmDaVOay8WXDqWXZrETEam33H2uu49NT08PuhQREYkhhaQaiA5IGy6dzKnnDA66JBERERERqWUKSdWUs6OAgseHKCCJiIiIiNRzmrihGnLyQ4x84mP6lFzOzZeeo4AkIiIiIlKPKSQdxOacDTw85SXyQl24JvO/OFVjkERERERE6jWFpCpsztlA4aTBPFC6jRtGL6CHApKIiIiISL2nMUmVKAtI7Uo3sfnSSfTo1CHokkREREREJAYUkioQHZA2XDqZTudcFnRJIiIiIiISIwpJ+8nJDzHnqd8rIImIiIiIJCmNSYqSkx9i5KQPyds9mH7DrqNrjz5BlyQiIiIiIjGmlqSI3NyNrH/0MhoXfs2UzLMVkEREREREkpRakggHpMJJg8koyeaRQcdomm8RERERkSSW9C1JZQGpfUk2Gy55mlPP1o1iRURERESSWVKHpP0DUqe+Q4IuSUREREREApa0ISknP8SYZz5lW0lTNlzylAKSiIiIiIgASTomKTc3m5ueWUZ20RE0GjOXTie2DLokERERERGJE0kXknJzsymcNIh7So6i6ZhZ9FJAEhERERGRKEkRkuatm8fExRPJLcrlmL3OHak76NbvV2pBEhERqUNzlmTz8Pw1bNoZom2LNO4a2JmhGe2CLktE5KDq/Zikeevm8cCCB8gpysFxth4Bv2l9HJ+3Sop8KCIiEog5S7IZP2s52TtDOJC9M8T4WcuZsyQ76NJERA6q3oekiYsnUlxSXG7Zt76HiYsnBlSRiIhI/ffw/DWE9pSUWxbaU8LD89cEVJGISPXV+5CUW5Rbo+UiIiJy+DbtDNVouYhIPKn3Ial1k9Y1Wi4iIiKHr22LtBotFxGJJ/U+JI3rNY7UlNRyy1JTUhnXa1xAFYmIiNR/dw3sTFrDlHLL0hqmcNfAzgFVJCJSffV+9oLBJw8G2De7XesmrRnXa9y+5SIiIlL7ymax0+x2IpKI6n1IgnBQUigSERGJraEZ7RSKRCQh1fvudiIiIiIiIjWhkCQiIiIiIhJFIUlERERERCSKQpKIiIiIiEgUhSQREREREZEoCkkiIiIiIiJRFJJERERERESiKCSJiIiIiIhEUUgSERERERGJopAkIiIiIiISRSFJREREREQkirl70DVUm5ltBb4Kuo5DdAywLegiApTMx5/Mxw46/kQ6/hPd/digi4gnZjYEGAKMANYGXE5tSqSfy1jROamYzsuBdE4qlojnpdLrXkKFpERmZh+7e++g6whKMh9/Mh876PiT/fglPunn8kA6JxXTeTmQzknF6tt5UXc7ERERERGRKApJIiIiIiIiURSSYmdS0AUELJmPP5mPHXT8yX78Ep/0c3kgnZOK6bwcSOekYvXqvGhMkoiIiIiISBS1JImIiIiIiERRSIohM7vKzFaaWamZ1ZvZP6piZpeY2Roz+8LMfhZ0PbFkZk+Z2RYzWxF0LUEws+PN7G0z+yzycz8u6JpixcxSzWyhmX0aOfZfBl2TyP6S8ZpUmWS+VlUm2a9hFUnm61pl6vP1TiEptlYAw4B3gy4kFswsBfgzcCnQFRhpZl2DrSqmJgOXBF1EgPYC/8/duwLnALcn0ff/LfBdd+8B9AQuMbNzgi1J5ABJdU2qjK5VlZpMcl/DKpLM17XK1NvrnUJSDLn7KndfE3QdMdQH+MLd17n7buB54IqAa4oZd38X2B50HUFx9xx3Xxx5XgCsAtoFW1VseFhh5GXDyEMDQCWuJOE1qTJJfa2qTLJfwyqSzNe1ytTn651CktSldsCGqNcbSfL/mSQrM+sAZAAfBVxKzJhZipktBbYAb7h70hy7SILRtUpqLBmva5Wpr9e7I4IuoL4xszeB1hWsutfdX451PSJBM7OmwEvAne7+TdD1xIq7lwA9zawFMNvMurm7+vZLTOmaJFL7kvW6Vpn6er1TSKpl7v69oGuII9nA8VGv20eWSZIws4aELyTT3H1W0PUEwd13mtnbhPv2J/xFQxKLrknVomuVVJuua5Wrb9c7dbeTurQI6GRmJ5lZI+Aa4B8B1yQxYmYGPAmscvcJQdcTS2Z2bOQvaphZGnARsDrQokSkMrpWSbUk83WtMvX5eqeQFENmdqWZbQT6AvPMbH7QNdUld98L/AiYT3hw4wvuvjLYqmLHzKYD/wI6m9lGM8sMuqYYOw+4DviumS2NPAYFXVSMtAHeNrNlhH8Be8PdXwm4JpFyku2aVJlkv1ZVRtewCiXzda0y9fZ6Z+71YgIKERERERGRWqGWJBERERERkSgKSSIiIiIiIlEUkkRERERERKIoJImIiIiIiERRSBIREREREYmikCR1ysxKoqbJXGpmHcxsQTU+94SZdY08v+cQ9ptlZmvMbJmZrTazx8rm8Y+sXxD1/GEzWxn591gz+8jMlphZv5ruN1bM7L8ix7XczD41swmRG9wd6vY6mNmKyPPeZvboYWyrxt+XiEg80bWrbujaJYlEU4BLnTKzQndvGuttmFkW8N/u/nHk5oAPAb3d/TsVvDcfaOnuJWZ2DfA9d7+5BvtKcfeSmtR3OMzsB8BQ4JrI3a0bAf8F/MXdvzmU2sysA/CKu3erhfoO+zsXEQmSrl21T9cuSTRqSZKYM7PCyL/9I381mxn5y9K0yN2sy/6a1tvMfgukRf6SNy2ybrSZLYws+7uZpVS1P3ffDdwNnGBmPfar4R9AU+ATM/sp8Hvgisi208zsYjP7l5ktNrMXzaxp5HPrzex3ZrYYuOog7/tlZPlyMzstsrypmT0dWbbMzL4fWV7hdvZzL3Cbu+8sOz53/23ZRcbMCs3sf83sU6Cvmd1vZovMbIWZTYo6x2dG/pL3KXB71PfT38xeiTxvYmZPRc73EjO7IrL8RjObZWavmdlaM/t9ZPkB35eISH2ga5euXZJk3F0PPersAZQASyOP2ZFlhZF/+wP5QHvCgf1fwPmRdVmE/3q27/2R512AuUDDyOu/ANdXsN99n49aNgcYUcE2o5/fCDwWeX4M8C7QJPL6p8D9kefrgbur+b47Is9/CDwRef474JGo/R5V1Xai3tcc2HGQc+7A1VGvW0Y9nwoMiTxfBlwQef4wsCLqe3kl8vx/gNGR5y2Az4EmkfO0DkgHUoGvgOP3P5966KGHHon40LVL1y499DgCkboVcveeVaxf6O4bAcxsKdABeL+K918InAksivxRKQ3YUs1arJrvK3MO0BX4ILKvRoQvhmVmVPN9syL/fgIMizz/HnBN2RvcfYeZXXaQ7Rx4QGYDCV+0WgDXuvsCwhf3l6LeNsDM7gYaAy2BlWb2HtDC3d+NvGcqcGkFu7gYuNzM/jvyOhU4IfL8n+6eH6njM+BEYENV9YqIJAhdu3TtkiSnkCRB+zbqeQkH/5k0YIq7j6/JTiLdGroDq2ryMeANdx9Zyfqiar6v7BgPdnwH2w7u/k2kS8JJ7v5vd58PzI90MWgUeVuxR/pym1kq4b9Y9nb3DWb2AOGLRXUZ8H13X1NuodnZ1Py7ExGpL3Ttqv7+dO2ShKQxSZII9th/Zr/5JzDczI4DMLOWZnZiVR+OfPYhYIO7L6vBfj8EzjOzjpHtNDGzUw/jfdHeoHxf6qNqsJ2HgL9aZMajSD/tyi4eZcu3RfqIDwfwcJ/wnWZ2fmT9qEo+Px+4I6oveMZBjgvKf18iIslK167ydO2ShKKQJIlgErDMzKa5+2fAz4HXzWwZ4f9ht6nkc9Mi71lBuC/yFTXZqbtvJdx/eXpkO/8CTjvU9+3nQeCoyIDUT4EBNdjOXwlfcD+KvO8DYEnksX9tO4HHCZ+D+cCiqNU3AX+OdBWprDvHr4GGhM//ysjrg9n3fVXjvSIi9ZWuXeXp2iUJRVOAi4iIiIiIRFFLkoiIiIiISBSFJBERERERkSgKSSIiIiIiIlEUkkRERERERKIoJImIiIiIiERRSBIREREREYmikCQiIiIiIhJFIUlERERERCTK/wdVw8AcK6xpUQAAAABJRU5ErkJggg==\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 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}