{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Welcome to Caustic!\n",
    "\n",
    "In need of a differentiable strong gravitational lensing simulation package? Look no further! We have all your lensing simulator needs. In this tutorial we will cover the basics of caustic and how to get going making your own lensing configurations. Caustic is easy to use and very powerful, you will get to see some of that power here, but there will be more notebooks which demo specific use cases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import torch\n",
    "from torch.nn.functional import avg_pool2d\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.io import fits\n",
    "import numpy as np\n",
    "\n",
    "import caustic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "FlatLambdaCDM(\n",
       "    name='cosmo',\n",
       "    static=[h0, critical_density_0, Om0],\n",
       "    dynamic=[],\n",
       "    x keys=[]\n",
       ")"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Specify the image/cosmology parameters\n",
    "n_pix = 100\n",
    "res = 0.05\n",
    "upsample_factor = 2\n",
    "fov = res * n_pix\n",
    "thx, thy = caustic.get_meshgrid(res/upsample_factor, upsample_factor*n_pix, upsample_factor*n_pix, dtype=torch.float32)\n",
    "z_l = torch.tensor(0.5, dtype=torch.float32)\n",
    "z_s = torch.tensor(1.5, dtype=torch.float32)\n",
    "cosmology = caustic.FlatLambdaCDM(name = \"cosmo\")\n",
    "cosmology.to(dtype=torch.float32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulating an SIE lens\n",
    "\n",
    "Here we will demo the very basics of lensing with a classic `SIE` lens model. We will see what it takes to make an `SIE` model, lens a backgorund `Sersic` source, and sample some examples in a simulator. Caustic simulators can generalize to very complex scenarios. In these cases there can be a lot of parameters moving through the simulator, and the order/number of parameters may change depending on what lens or source is being used. To streamline this process, caustic impliments a class called `Parametrized` which has some knowledge of the parameters moving through it, this way it can keep track of everything for you. For this to work, you must put the parameters into a `Packed` object which it can recognize, each sub function can then unpack the parameters it needs. Below we will show some examples of what this looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# demo simulator with sersic source, SIE lens. then sample some examples. demo the model graph\n",
    "\n",
    "class Simple_Sim(caustic.Simulator):\n",
    "    def __init__(\n",
    "        self,\n",
    "        lens,\n",
    "        src,\n",
    "        z_s=None,\n",
    "        name: str = \"sim\",\n",
    "    ):\n",
    "        super().__init__(name) # need this so `Parametrized` can do its magic\n",
    "        \n",
    "        # These are the lens and source objects to keep track of\n",
    "        self.lens = lens\n",
    "        self.src = src\n",
    "        \n",
    "        # Here we can add a parameter to the simulator, in this case it is `z_s` which we will need later\n",
    "        self.add_param(\"z_s\", z_s)\n",
    "\n",
    "    def forward(self, params):# define the forward model\n",
    "        # Here the simulator unpacks the parameter it needs\n",
    "        z_s = self.unpack(params)\n",
    "\n",
    "        # Note this is very similar to before, except the packed up `x` is all the raytrace function needs to work\n",
    "        bx, by = self.lens.raytrace(thx, thy, z_s, params)\n",
    "        mu_fine = self.src.brightness(bx, by, params)\n",
    "        \n",
    "        # We return the sampled brightness at each pixel location\n",
    "        return avg_pool2d(mu_fine.squeeze()[None, None], upsample_factor)[0, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\n",
       "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\n",
       " \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\n",
       "<!-- Generated by graphviz version 2.43.0 (0)\n",
       " -->\n",
       "<!-- Title: %3 Pages: 1 -->\n",
       "<svg width=\"926pt\" height=\"188pt\"\n",
       " viewBox=\"0.00 0.00 926.00 188.00\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 184)\">\n",
       "<title>%3</title>\n",
       "<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-184 922,-184 922,4 -4,4\"/>\n",
       "<!-- sim -->\n",
       "<g id=\"node1\" class=\"node\">\n",
       "<title>sim</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"243\" cy=\"-162\" rx=\"93.68\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"243\" y=\"-158.3\" font-family=\"Times,serif\" font-size=\"14.00\">Simple_Sim(&#39;sim&#39;)</text>\n",
       "</g>\n",
       "<!-- sim/z_s -->\n",
       "<g id=\"node2\" class=\"node\">\n",
       "<title>sim/z_s</title>\n",
       "<polygon fill=\"lightgrey\" stroke=\"lightgrey\" points=\"173,-108 119,-108 119,-72 173,-72 173,-108\"/>\n",
       "<text text-anchor=\"middle\" x=\"146\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">z_s</text>\n",
       "</g>\n",
       "<!-- sim&#45;&gt;sim/z_s -->\n",
       "<g id=\"edge1\" class=\"edge\">\n",
       "<title>sim&#45;&gt;sim/z_s</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M220.01,-144.41C207.49,-135.37 191.78,-124.04 178.07,-114.14\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"179.99,-111.21 169.84,-108.2 175.9,-116.89 179.99,-111.21\"/>\n",
       "</g>\n",
       "<!-- sie -->\n",
       "<g id=\"node3\" class=\"node\">\n",
       "<title>sie</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"243\" cy=\"-90\" rx=\"51.99\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"243\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">SIE(&#39;sie&#39;)</text>\n",
       "</g>\n",
       "<!-- sim&#45;&gt;sie -->\n",
       "<g id=\"edge2\" class=\"edge\">\n",
       "<title>sim&#45;&gt;sie</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M243,-143.7C243,-135.98 243,-126.71 243,-118.11\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"246.5,-118.1 243,-108.1 239.5,-118.1 246.5,-118.1\"/>\n",
       "</g>\n",
       "<!-- src -->\n",
       "<g id=\"node10\" class=\"node\">\n",
       "<title>src</title>\n",
       "<ellipse fill=\"none\" stroke=\"black\" cx=\"639\" cy=\"-90\" rx=\"65.79\" ry=\"18\"/>\n",
       "<text text-anchor=\"middle\" x=\"639\" y=\"-86.3\" font-family=\"Times,serif\" font-size=\"14.00\">Sersic(&#39;src&#39;)</text>\n",
       "</g>\n",
       "<!-- sim&#45;&gt;src -->\n",
       "<g id=\"edge9\" class=\"edge\">\n",
       "<title>sim&#45;&gt;src</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M309.28,-149.28C383.41,-136.18 502.61,-115.11 575.56,-102.21\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"576.3,-105.64 585.54,-100.45 575.08,-98.74 576.3,-105.64\"/>\n",
       "</g>\n",
       "<!-- sie/z_l -->\n",
       "<g id=\"node4\" class=\"node\">\n",
       "<title>sie/z_l</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"54,-36 0,-36 0,0 54,0 54,-36\"/>\n",
       "<text text-anchor=\"middle\" x=\"27\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">z_l</text>\n",
       "</g>\n",
       "<!-- sie&#45;&gt;sie/z_l -->\n",
       "<g id=\"edge3\" class=\"edge\">\n",
       "<title>sie&#45;&gt;sie/z_l</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M203.92,-77.84C160.62,-65.43 93.75,-46.11 63.85,-35.88\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"64.92,-32.55 54.33,-32.44 62.54,-39.13 64.92,-32.55\"/>\n",
       "</g>\n",
       "<!-- sie/x0 -->\n",
       "<g id=\"node5\" class=\"node\">\n",
       "<title>sie/x0</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"126,-36 72,-36 72,0 126,0 126,-36\"/>\n",
       "<text text-anchor=\"middle\" x=\"99\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x0</text>\n",
       "</g>\n",
       "<!-- sie&#45;&gt;sie/x0 -->\n",
       "<g id=\"edge4\" class=\"edge\">\n",
       "<title>sie&#45;&gt;sie/x0</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M214.21,-75C191.6,-64.01 159.97,-48.64 135.58,-36.78\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"136.8,-33.49 126.28,-32.26 133.74,-39.78 136.8,-33.49\"/>\n",
       "</g>\n",
       "<!-- sie/y0 -->\n",
       "<g id=\"node6\" class=\"node\">\n",
       "<title>sie/y0</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"198,-36 144,-36 144,0 198,0 198,-36\"/>\n",
       "<text text-anchor=\"middle\" x=\"171\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">y0</text>\n",
       "</g>\n",
       "<!-- sie&#45;&gt;sie/y0 -->\n",
       "<g id=\"edge5\" class=\"edge\">\n",
       "<title>sie&#45;&gt;sie/y0</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M226.3,-72.76C217.3,-64.01 206.01,-53.03 195.98,-43.28\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"198.22,-40.59 188.61,-36.12 193.34,-45.6 198.22,-40.59\"/>\n",
       "</g>\n",
       "<!-- sie/q -->\n",
       "<g id=\"node7\" class=\"node\">\n",
       "<title>sie/q</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"270,-36 216,-36 216,0 270,0 270,-36\"/>\n",
       "<text text-anchor=\"middle\" x=\"243\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">q</text>\n",
       "</g>\n",
       "<!-- sie&#45;&gt;sie/q -->\n",
       "<g id=\"edge6\" class=\"edge\">\n",
       "<title>sie&#45;&gt;sie/q</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M243,-71.7C243,-63.98 243,-54.71 243,-46.11\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"246.5,-46.1 243,-36.1 239.5,-46.1 246.5,-46.1\"/>\n",
       "</g>\n",
       "<!-- sie/phi -->\n",
       "<g id=\"node8\" class=\"node\">\n",
       "<title>sie/phi</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"342,-36 288,-36 288,0 342,0 342,-36\"/>\n",
       "<text text-anchor=\"middle\" x=\"315\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">phi</text>\n",
       "</g>\n",
       "<!-- sie&#45;&gt;sie/phi -->\n",
       "<g id=\"edge7\" class=\"edge\">\n",
       "<title>sie&#45;&gt;sie/phi</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M259.7,-72.76C268.7,-64.01 279.99,-53.03 290.02,-43.28\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"292.66,-45.6 297.39,-36.12 287.78,-40.59 292.66,-45.6\"/>\n",
       "</g>\n",
       "<!-- sie/b -->\n",
       "<g id=\"node9\" class=\"node\">\n",
       "<title>sie/b</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"414,-36 360,-36 360,0 414,0 414,-36\"/>\n",
       "<text text-anchor=\"middle\" x=\"387\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">b</text>\n",
       "</g>\n",
       "<!-- sie&#45;&gt;sie/b -->\n",
       "<g id=\"edge8\" class=\"edge\">\n",
       "<title>sie&#45;&gt;sie/b</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M271.79,-75C294.4,-64.01 326.03,-48.64 350.42,-36.78\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"352.26,-39.78 359.72,-32.26 349.2,-33.49 352.26,-39.78\"/>\n",
       "</g>\n",
       "<!-- src/x0 -->\n",
       "<g id=\"node11\" class=\"node\">\n",
       "<title>src/x0</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"486,-36 432,-36 432,0 486,0 486,-36\"/>\n",
       "<text text-anchor=\"middle\" x=\"459\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">x0</text>\n",
       "</g>\n",
       "<!-- src&#45;&gt;src/x0 -->\n",
       "<g id=\"edge10\" class=\"edge\">\n",
       "<title>src&#45;&gt;src/x0</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M599.52,-75.53C572.48,-66.07 535.68,-52.73 495.47,-36.03\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"496.59,-32.71 486.01,-32.07 493.88,-39.16 496.59,-32.71\"/>\n",
       "</g>\n",
       "<!-- src/y0 -->\n",
       "<g id=\"node12\" class=\"node\">\n",
       "<title>src/y0</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"558,-36 504,-36 504,0 558,0 558,-36\"/>\n",
       "<text text-anchor=\"middle\" x=\"531\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">y0</text>\n",
       "</g>\n",
       "<!-- src&#45;&gt;src/y0 -->\n",
       "<g id=\"edge11\" class=\"edge\">\n",
       "<title>src&#45;&gt;src/y0</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M614.49,-73.12C600.11,-63.79 581.69,-51.85 565.81,-41.56\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"567.59,-38.55 557.3,-36.04 563.78,-44.42 567.59,-38.55\"/>\n",
       "</g>\n",
       "<!-- src/q -->\n",
       "<g id=\"node13\" class=\"node\">\n",
       "<title>src/q</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"630,-36 576,-36 576,0 630,0 630,-36\"/>\n",
       "<text text-anchor=\"middle\" x=\"603\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">q</text>\n",
       "</g>\n",
       "<!-- src&#45;&gt;src/q -->\n",
       "<g id=\"edge12\" class=\"edge\">\n",
       "<title>src&#45;&gt;src/q</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M630.29,-72.05C626.15,-64.01 621.11,-54.2 616.49,-45.23\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"619.57,-43.57 611.88,-36.28 613.35,-46.77 619.57,-43.57\"/>\n",
       "</g>\n",
       "<!-- src/phi -->\n",
       "<g id=\"node14\" class=\"node\">\n",
       "<title>src/phi</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"702,-36 648,-36 648,0 702,0 702,-36\"/>\n",
       "<text text-anchor=\"middle\" x=\"675\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">phi</text>\n",
       "</g>\n",
       "<!-- src&#45;&gt;src/phi -->\n",
       "<g id=\"edge13\" class=\"edge\">\n",
       "<title>src&#45;&gt;src/phi</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M647.71,-72.05C651.85,-64.01 656.89,-54.2 661.51,-45.23\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"664.65,-46.77 666.12,-36.28 658.43,-43.57 664.65,-46.77\"/>\n",
       "</g>\n",
       "<!-- src/n -->\n",
       "<g id=\"node15\" class=\"node\">\n",
       "<title>src/n</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"774,-36 720,-36 720,0 774,0 774,-36\"/>\n",
       "<text text-anchor=\"middle\" x=\"747\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">n</text>\n",
       "</g>\n",
       "<!-- src&#45;&gt;src/n -->\n",
       "<g id=\"edge14\" class=\"edge\">\n",
       "<title>src&#45;&gt;src/n</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M663.51,-73.12C677.89,-63.79 696.31,-51.85 712.19,-41.56\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"714.22,-44.42 720.7,-36.04 710.41,-38.55 714.22,-44.42\"/>\n",
       "</g>\n",
       "<!-- src/Re -->\n",
       "<g id=\"node16\" class=\"node\">\n",
       "<title>src/Re</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"846,-36 792,-36 792,0 846,0 846,-36\"/>\n",
       "<text text-anchor=\"middle\" x=\"819\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">Re</text>\n",
       "</g>\n",
       "<!-- src&#45;&gt;src/Re -->\n",
       "<g id=\"edge15\" class=\"edge\">\n",
       "<title>src&#45;&gt;src/Re</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M678.48,-75.53C705.52,-66.07 742.32,-52.73 782.53,-36.03\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"784.12,-39.16 791.99,-32.07 781.41,-32.71 784.12,-39.16\"/>\n",
       "</g>\n",
       "<!-- src/Ie -->\n",
       "<g id=\"node17\" class=\"node\">\n",
       "<title>src/Ie</title>\n",
       "<polygon fill=\"none\" stroke=\"black\" points=\"918,-36 864,-36 864,0 918,0 918,-36\"/>\n",
       "<text text-anchor=\"middle\" x=\"891\" y=\"-14.3\" font-family=\"Times,serif\" font-size=\"14.00\">Ie</text>\n",
       "</g>\n",
       "<!-- src&#45;&gt;src/Ie -->\n",
       "<g id=\"edge16\" class=\"edge\">\n",
       "<title>src&#45;&gt;src/Ie</title>\n",
       "<path fill=\"none\" stroke=\"black\" d=\"M691.8,-79.13C733.82,-70.59 794.1,-56.79 854.1,-35.98\"/>\n",
       "<polygon fill=\"black\" stroke=\"black\" points=\"855.46,-39.21 863.72,-32.57 853.13,-32.61 855.46,-39.21\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>\n"
      ],
      "text/plain": [
       "<graphviz.graphs.Digraph at 0x7fca9dbb87f0>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sie = caustic.lenses.SIE(cosmology, name = \"sie\")\n",
    "src = caustic.sources.Sersic(name = \"src\")\n",
    "\n",
    "sim = Simple_Sim(sie, src, torch.tensor(0.8))\n",
    "\n",
    "sim.get_graph(True, True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simple_Sim(\n",
      "    name='sim',\n",
      "    static=[z_s],\n",
      "    dynamic=[],\n",
      "    x keys=[('sie': ['z_l', 'x0', 'y0', 'q', 'phi', 'b']), ('src': ['x0', 'y0', 'q', 'phi', 'n', 'Re', 'Ie'])]\n",
      ")\n",
      "SIE(\n",
      "    name='sie',\n",
      "    static=[],\n",
      "    dynamic=[z_l, x0, y0, q, phi, b],\n",
      "    x keys=[('sie': ['z_l', 'x0', 'y0', 'q', 'phi', 'b'])]\n",
      ")\n"
     ]
    }
   ],
   "source": [
    "print(sim)\n",
    "print(sie)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGdCAYAAABU0qcqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABTvUlEQVR4nO29e4xd1Xn3/z1nbh7fxsGEsR1smEZIJkAKwVwMqM2vWLVS1EJx0yI5lUOi0iQ2YCyF4DYmSgOYpGpDSRooqKXkVwgNaiEB/QriNa3zopqbKTSUxtAXFHgDM05E7DHGnss5+/fHwJm1nufsZ511zj6zztjfj2Tp7Flrr7325Zzl/XyfSynLsgyEEELIDFNOPQFCCCFHJ1yACCGEJIELECGEkCRwASKEEJIELkCEEEKSwAWIEEJIErgAEUIISQIXIEIIIUnoTj0BSbVaxZtvvokFCxagVCqlng4hhJBIsizDgQMHsGzZMpTL+e85HbcAvfnmm1i+fHnqaRBCCGmRN954A8cff3xue8ctQAsWLAAAfHzJ5egu90790X0Tkm9FZbFt9RXbmRqrnNs35jh63JxjhPoCgNssjyO6unNU40qMcUNkTb6YlmKSPskMUYF9S27/qt9ZTdftK8etNngMOU69bWesqH2rEX1FW9RxQlm4qkZfa9/Q/CGbjXNXnY0bZB1Hzd8YRxKYk5nNrBpxnOD9iBgrYlj/GC1kZnOu6WQ2gf9d+WHt9zyPjluA3je7dZd70V3ue/+Pbgd/h1YWILlvRy5A1nFEXy5AU8gFSO5r/Qgb5xa9ALm3LmZfOeGYBaiFxUvT5AIkV/GYBSh0o5tdgNQPa8SPeeDBzaw5l2IWjdBD3uQCFBrXO0YrqUH1/EIySsctQDVKpegfxrpjOKgFp5XxS8aPvbGIRC04al+/KesSO0csKlnMAhThqiLPz/3hDT7a3vMrTzb0w+T0L/t91bpgLVbyXN35q9+w/HN97w/5c6ggHzmHmP9Eq3sl5hjzVuOOJecgnxl3X2nzD/zP3R1Jz19QdTrIxcick2irBg7kjq32Fc+M+x1V1zDiWljzrzdWg+PKRcB8Ywucq4n35WnsR4NecIQQQpLABYgQQkgSOtcE5+K+QobMaEW5bseY6yLMaFF9AWRdhvlOTskw12lNyzlGaFxBUF9y+4YGcw9bNl73Q8KTZ1YLmDJizHWVadNGSc6hZO/rzkOa5zJpa69kuX1NTShG84Fv4oqSFLoC19QayzIdAZ75KPS0eCa6kBnNmpQyc8prbpj6DDNV0NwVY0YLmeSscTvBJNcAfAMihBCSBC5AhBBCksAFiBBCSBJmhwY0U8RoG1ZfKy4oGBxrHEe2WTpV2e5rjSux4o2KxNOLArZmLYs4eovUktRFzdeLZBxN5mgBpjs3gFIlX7sJunC7faX93nL3DsUMGRqR0qGsmKKQG3yzLtuAr1+kctm2xgppTWVnbMNFWxLtsu2OFbwfEdfUciNX4xarCfENiBBCSBK4ABFCCEkCFyBCCCFJOLo1oJiYoQjNR2Es8zHpdUKphNyYoZBe5B3Xaqt3HOt0W5GHXHNyl9Q9RFeVhNNpk6bzGL1Imsurht1dpjwz0hBZ+pAautTldxUBO17fiHRAaueKrQ2YMUOWrhPUpYyDdkrMkJnGx4gZKgeCq1qJGXK1nE6IEQL8a9OEHsQ3IEIIIUngAkQIISQJXIAIIYQk4ejSgIrSfICArmNoNcE6Q8Y8LM0HQudR2pKh46iyFXIOYixTD8tvikLKCPIwlsyQ2fqRG9+i9CJxP/xYGHFyUtcxtA5LHwKERiT1IfGwWXFAUWUfpN7YUsyQcRyJVW5ipmKGLB1HTUpgjSXHsfLIBWKGTE0olDfOHyh/nHpjxczJG9fVSxv7IeAbECGEkCRwASKEEJKEo8sEF8JIMROTeifK7TqmbHjANOa5bFvmOTkPaTVQpj7jODGlv2MIVSeW7a4VRLlSy5IE+aYCmYrHLcGgXanlHAwTXYQLd0m4RysziGsyCbhhZ/JARomIVso+eCa62DLhxr7tctlOlsanyUqrgDB/tavSakTaHjWnJuAbECGEkCRwASKEEJIELkCEEEKScERrQCp1TZE069IdcLvWpQ+cNqnjyG23fLfqK45jumwbfeUcQ/MvqkS6Sr0jD+y2yX3zdR6tF4ld3ewnlos2AhpRwIXbcw2XrryGy3ApVBbcSOsTLBHh5/wRfdEw8tnrBJfttqXxUXOQOo8xlkzjY2hCbSv1HaMthebUAHwDIoQQkgQuQIQQQpLABYgQQkgSjmgNaKYIxsIUNLaKzzH0F1Pzga/zhPUiMSl3TiqGSPZFLla2DqXjKLEpv7/SGFRcihsLI8aR5nI3lY0q1y00E0MjCsYQuccVWoyepGuzb7wsOOCn9QnFAWXOuZekNlBUzJBon7GYobal8THKOgB2+QKp/RmlHUI6WlTMkFXqu8m0PY3+IvINiBBCSBK4ABFCCEkCTXB5tGJWM0xlodQ7duVVsW25ZUe4ViszmkrjI/d1PsvjGG7mwQS5nt0j0NVMxWO7bLvmC2VyU9mxjTaZadow0anbLPb1zBehKq2V/DZlFjTS+mTSnCLdsHEEuWxLuuQXQl7zJtP4KDfrJrNqA3Zm7Yg0PoVWWo00ZYbgGxAhhJAkcAEihBCSBC5AhBBCktC5GlCWTdsfm9RjpBYgtY2OxCi5YLldq76hlD+ujhPQfKox+pE1/9B/d6zbY2g+gNCEDM1Hbiu9SGobhl6kPJNNF26/Td4f32VbjCN1HfdkjZQ+gJ3WpyQfNumG7Zxg0L3b1eAC42q9Jd+1un0u2wGRsck0Ptpl29BxYso6yLGsSquAf5+LrLTaQimHevANiBBCSBK4ABFCCEkCFyBCCCFJ6FwNyMWzO84CHccipONAtlulHIyYm1AqnnJ+X6X5GBqR7KvKexvHscuci21l3pf2f7dN9DVKLkjNpyz1F0+bsbUyta8VQzQp5uSVuLBjPLw4ICulj+wrxpb6iooZch5Oq7Q3AJQc+78q4R4oG+5NuhNjhmLS0YjtoCZkTcJM42NoS6pvgaW+Gy3lEBR839uloV6EEEJIwXABIoQQkoTZYYJzkeaIdrlWK5fOmTH9xWTW1m7ZRioeI2VOyFynTXL1P9cdq8s4jpU5O1iy0qiIGkqvY7g8y/NxzW7a7Vq6Leeb89QcpPXUzbIyGci67d7mGJdtwHPbVuYiqxKraFOVWJt02Z6aRzm3b5TLtjBhJUnjE3DZ9iutBlLvWCa5UKXVmUjbM7Vz/fk1CN+ACCGEJIELECGEkCRwASKEEJKE2acBSWK0mkS6TktY5QysVDzS1TIiRU5Y15n+XJUu2t35fUNpe0zPzaBbtrMR1ICc0gdSQxHbZUfPyFRfofmo85net1yxtRnPHTxQusF3w45w2ZZjR7hs6/oXxbhsA+J5a8VlW9yAzEjF07Y0PgHNxKy0KrFKO1gu2nLfdqXtAYxSDo39tvINiBBCSBK4ABFCCEkCFyBCCCFJmP0akMTLCx+p8bi23dlQusFC6UPCtmvpRYFy3q7uIzWfqqEByTazRETkf43iNCDns9R8RIqcrJKvF6l4HZmap5yvoah4HXeOoYxLznFiYoYA/xzaFTMUU+Zhak75eZSU7laOiBmq5AfzdEIan2Com5nGpwPS9gAsx0AIIWR2wgWIEEJIErgAEUIISULnakDVDDU7p8y9ROriaSihMg9WmYRQfI5Xzlv0lZpQl9Em9SInbkaVdQicj1+SW7TJGBxHV1Caj8oFl99XxlOUJqV93BlX2NJVXrlKvnYptRk3pqgqNR6rzIMcu8iYIQelQxllHgARNyTjgMSBPQkiEF+ELlefsPUis/S3oSUB8J/VkAwSoZnYeeTS542bOkx8/jcXvgERQghJAhcgQgghSehcE5yL+ZoXkXpH+f02nppHpgTxvL1VWnhrTvlN0UR4ipsmrIgUP4BvVpOlGlTpBucJq/b4J1/tEX29Mg/SRBKYo1WOQVU5deYgTFblCWEqc9sDpsmyrDrr7FtWZVrFppH2Rpo93H2Va7h8wCpy34KqC8s5uS7bskyFunDi/IxyDHLfkv/Fa7qv5aItybrEuJY5L8Y13CrrAMSl8UmQtkfBiqiEEEJmC1yACCGEJIELECGEkCTMDg3IIqJEt7JNW+l2WindENKeiqIVPSlCE5Klpr32gC5S9TQg0WZpQt3yXkWcbMgN29F5SpNCt5Hak9NXaz7isLLdy05vaTEyk73QocQJVb02UT49lMbH9bgV40aV/jbTz4ht6bItPYadAymXben2Wy3IZVtMySy5IMtHtJLGx3INl0Sk8UmStgf+uTfjks03IEIIIUngAkQIISQJXIAIIYQkYfZrQJJWyjF4NszidBuvuq/RpucAZKEcNAWgyi8EYm7MND4qDmj6fGTqnWqvOFd3u0doAd1i29CEMmHvzyZFqWZXA5oQbbLEuJmySGg11nVTtzFf55E6juzr60VCS5JlwpUg5qTxkUexUvHElP6WcVhiWKU1OWO1LWZI6kNmX9nfnpObxieYtsclpqwDYKbxSZK2R/T3ypywJDchhJBOJmoBqlQq2LZtG4aGhtDf348Pf/jD+NrXvuZ5P2RZhuuvvx5Lly5Ff38/1qxZg1deeaXwiRNCCJndRJngvv71r+O2227D3XffjVNOOQXPPvssLr/8cgwMDOCqq64CAHzjG9/ArbfeirvvvhtDQ0PYtm0b1q5di5deeglz5sxp/GBZNv1qWDXSV4TGcJFmBKtdmhikrcy1k1htAHyrQUTaHoHK8tHwnu0jylwnXZyFGzZ6p+9z1xw/R05PT0Vs++1e4mYxqclJ/5mZGJ9+7CtjflvWLU1ybvlR25W6qNQ2ltv1VLvbFnLZzs+sLV241fl4LuqNn6t8JoKVV61hY1y25Txca3yoSqvKlGRUXrUqrXYFTIiWia6gSqvADKXtAbRJLpKoBejf//3fcfHFF+Oiiy4CAJx44on43ve+h6effnpqLlmGW265BV/+8pdx8cUXAwC++93vYnBwEA8++CAuu+yyliZLCCHkyCHKBHfeeedhx44dePnllwEAL7zwAp544gl84hOfAAC89tprGB4expo1a2r7DAwM4JxzzsGuXbvqjjk2NobR0VHvHyGEkCOfqDeg6667DqOjo1i5ciW6urpQqVRw4403Yv369QCA4eFhAMDg4KC33+DgYK1Nsn37dnz1q19tZu6EEEJmMVEL0Pe//33cc889uPfee3HKKafg+eefx+bNm7Fs2TJs2LChqQls3boVW7ZsqW2Pjo5i+fLl+TsE0z0UVJ4h1oXbwnLvVilAhL3cqRKq9CJpt3bdIAMSlle+IFZMiknj42xLF22U/QOXHZ2nt9fXeOb3j3nb83rHve2+run+FWH0PjTh5wA6ONZb+/xud6/XNtHlfyUy5ytSUQYDoyZE3Xajq1fRNaDruF6zYpiQJuTOqSralCbkpkMyyzoIbUbnuRH7GuUlDM1nqt1w2S5LvcXpq3Qc0Vdps1bKH4G7b8AN29WIVFkHWSG1IE2o0LQ9kpKrkbpaUmO/n1EL0Be/+EVcd911NS3ntNNOw09/+lNs374dGzZswJIlSwAAIyMjWLp0aW2/kZERnH766XXH7OvrQ19fX8w0CCGEHAFEaUDvvvsuymKl7erqQvW91XZoaAhLlizBjh07au2jo6N46qmnsHr16gKmSwgh5Egh6g3ot3/7t3HjjTdixYoVOOWUU/Af//Ef+Mu//Et85jOfATDlbrl582bccMMNOOmkk2pu2MuWLcMll1zSjvkTQgiZpUQtQN/61rewbds2fOELX8DevXuxbNky/PEf/zGuv/76Wp9rr70WBw8exBVXXIF9+/bhggsuwCOPPBIXA9QKloHcsCcDceUZSl7wLXLbAKHdSLuvyt9iIPYtyRLQno5gD2WGYlj6RANjNzywNOE7cUJ9Is5nYM5hb3uw/4C3fUzvwdzDvj0+z9v++aH5021dc722UfjP6YRjy1bJWlTtCUOvkOmBjMrNujKA0HWcDlaM0FS7pQmJUg7iDEuGrlMWpczNtEOhND4OuiRE8zFDVpsVMzTV3mzKH3ncfJ1HlXWQd8/ShFKl7VHNrUUjlrJmiji0kdHRUQwMDGDN4B+hu/yeQOw+lFJsU4Km21d+a+x9/aBDq7iL364WLrGvG9wo84eZgY8Asm5HNJZ5ysS+Vbdvj/hxkdvOfz0qvbKvt6naK331PwNAVWxPzpl+vCpzRT6ufn+7e95E7fP8ef6Cc+x8f4EpbAE6JBagg2IBOjR9MbLDvhdF+XBZbPvXqcvxkyiPi4VA/ICXJ+p/BoCS7Dvp5h4TbbJuj2rP31f+h8xzapG/Q5NyBW1snLrtrkOAcqyx+srjGNtWW72xjOOYDgxy0bAcDQLzV2Nl+edj/o/GGgeijo8aN+CU4PZ3+k5m49ix7//F/v37sXDhQuTBXHCEEEKS0LHZsLMsq63M3v8d5WrertQ8RVY1dcYqCXOKOoz6n4xjBgmk4il5x5HjykmV8tsC29Fu23moqqDT2z3d/n/PF/T4b0Qf6t/nbZ/c/2bt89yy77L9f8cXe9t7uv04NZfJiv9/snec7UnhipypbX+sqmN2U+7Fcrsrv01efy9aIPAfYWtfTcitfBr3bRuQ/3u3x7Errxbnsu1dHGX+Nd54EHLDzq+8qsxolpE0VGnVMsnFZtIuCmWLdedQzmnIh29AhBBCksAFiBBCSBK4ABFCCElCx2pADWMatQO6TYzOY5VnkOUXrPIMKit8vuYj56iqo1r7BnQby8NJ72voUnJf6fHkeUdJ/St/W7aVxQks6PI1oZW9b9U+nzvH1wWfPPy2t33YcfPbP9Hvtb0z4afmGZtwUvGo6qkx1VSFxiC9IY3Kq6ogp7OvHDek3/nXNbAvjL4Ko28o/METdjrAZRvwUszIcS037AzS5dkuA+HPKUITiknbU2TpBqv8gvtM67K+9XdpqBchhBBSMFyACCGEJIELECGEkCTMfg3IQqW9Ee1WqnclmuTrLyr1jmUuV5HEol3a9N1YEhm7YMQU6ehzOQ93HHtca6zQcbxtOQcRR+PGzch4nMMVPz3DgYqfsWBfdTqjQSU75LUtKPulG+Z2TccJ9Xf5aQd6yv4ky841l9dfPiMqA72h66jtUn6bKm9g9ZXXWD62bn85YfnsWQNZmpAVj1NvLHdT6kMpYobUdvOlv5Um5MaKyfLdEkMTkt9RhXmfiyrdAJjlvBuAb0CEEEKSwAWIEEJIErgAEUIISULnakDVDO/bmT3bp+rXptxwVnbZ+jPJxSrdoI30otnLTyfs4TIrrxsfEpHPLUYvku1BDcixeZcnA7nUJqdvtBt/AwCjY77m89ZhP8PuC90rpvtW93ptE5k/1i8np7NjHxLa0qSwcUeZtUNaR7PIR6SU26R1KGPb0pbUttKSDE1I6ZriPkv5xdV9QjFDJeMEVF+3a0TMECDOXWom+XnkQnFAbrvMRh6jCckYtJJVClzXSPdptnQDkF/OO3TM9w/dUC9CCCGkYLgAEUIISULnmuCapbDUPIG0N57ZwEi9I9tlanf5XwAzvY5097Zcw/OnAAhTWcjtWpoKqvltel83bY8wXUwIN+zx6YsxPu4/mgcO+5Xu3uoeQB5vjfttXcIs8vbEtAlu35ifimds0j9uxXEHl1VN9XUzzKlFlbCIxXjktbnOKmcQGNbZN5QeSLlEuyY66douzHfe86bc4OV3y00NIyYhn2llFnTGReME3bCddtNFG7BNclJ6sExe0lxqpSGSLtoVOali4RsQIYSQJHABIoQQkgQuQIQQQpJw5GlAFhGpeZT7Z5ehCalUNvmpeXRbRGoepR9Jd1BnXMt2Dl970q7UYlxx3LLrdirnNCnGcrziSxOyTezrlDuoHPYfzYNlv0yCdKOdqEwf6O3eeV6bLOVwaHLa9frAmK8tHRoXbtmT0+O6buKAr28BsEtRBHSQkqUXBfQ8k6K0p0C0gOfeLTpnsrfSPd3PgQO5Jy+vv5Ejx0rpA9T5fri/BfKCqxRA+Sm8yuL74Llhx6TtAXwtR6bhstyyrTLagP/cWqEpCLhluy7ZStyuD9+ACCGEJIELECGEkCRwASKEEJKEztWAsup0+VfHtqjL1Apc/3jp0x5KD+GJNYE4IGsspes4fQNxP1aJBTNGCLbmoPSvLF/vkpqPWY5B2KnLYjtzdSmZimdS6GHjTiyJsGlPln1t5h3/MJ4GdGDM14tkGETFOT+Z8md8zD9OZcwRsURJbojzkamG3Gsj9aJWSlxY8UXB+C/rGWkFV55QTXYqHl9PlW1yu1Tv49Rx5PfDvebBtENGWp9AzJAXFyc0n2q3/8y4mlBM2p567R4iFY+rAav9rHLeEaUbpvq39hDxDYgQQkgSuAARQghJQuea4NpByOXZdQsOZb9226V5S7psm3OwKzOWjGzYVmoe5UYu3bKd9hiTGyBNS9IVPN8tuyQSlZeFW7Y7x6wszRE+E/IyOia4sS7/sVbWVOd8K6LyalWa2RwTnEwdJOevXHkdk1woZREiTHD+vcsfB4Bt0lLmufbkCwoV73RNdKZ5DjDdymVqKj/UwDDPoY5p333GAy7b3q7iF9UyyVku2kAdk5xbuVRNKn9bnZv1O9JgFut68/CPw2zYhBBCOhguQIQQQpLABYgQQkgSZocGlDm2UVmx0nLLjq2WatjAtaZi7Gek5lHe3EbaHgB22QeVPsRNkSPapM7j2riV/VukExHtrlaj3K6FXbvcld+mKmWW823RFWkPl5fCSZNTEa6vpke9SK8jtQFX9ymPCw1Iul3Lc3f1L3mvDPd1q/yFbA/rRfnbSvOxXLhbSQcUwNSIxM3zCqYofUt8l7xh7DIoSiNy+0tdTYxVRb42U5a6juuGLVNRGZoPAC91VSb6mopLqHSDuy1fSWTaIcMt2ysiy4qohBBCOhkuQIQQQpLABYgQQkgSOlcDyrJpm26kb7o5poVrKrXKLwDwc4/IoBTDYB6I+zFLdgdS8XglucU4suyDnfqlhVQ8okSEm35Ha0ly3/zSzNLKXZVaYI+jf8mUPyq2xBlVpvQX226sj9J8RBxQeSJfE5JamaUJaX3IuHehkuhWDFFM2Yf2hAi9d6DmdtPaUf5ASiuW6adknJwbQ6TifsQz7nSW96oq/o/va0JWXQQdr+PqPqVMiq8yhY5RukGeqyXCFfXbmwPfgAghhCSBCxAhhJAkdK4JrlEM92gzUzZgu2WHKgM6xzWrpYrtkNu1WTHVqpYKeOY7q1rq1B+MNDERbtnqdV6a1ZyM17KvNMn5N8ywmwHKVOll3Zb3wzLBKRdnw11dulnL1Dyy3dpXVo71zHX2vYtzwzbMd0VWXo0glJrH7JtvWTKR56LNUKK/cw902hux6dw79RWVv0/upFVmfMOMBvHbkAVctt05R6TtUT9IKpW84ZbtzoFu2IQQQjoZLkCEEEKSwAWIEEJIEjpXA8qyaX3HWyalkbuFNdR0yw65YRttVnXFgNu1Trdj2KJNd285pXy37LCOkO+WLfUKSxOSuof6749pNhaN3fm2dZXyx9KAZFdLb1EVUAP7Tjbe172OUWl7DI2n3jaMR1GnkEIhqMqlskKq+znwTLhjyWdN4fZVkxJd5bm7FUUDFVF9t2z5PRPlPpwvppRw1YUSVYHd+2Gl6VFjyetkuGWbLtlA4W7ZfAMihBCSBC5AhBBCksAFiBBCSBI6VwNqFjc+xyrVEDEOgDp6hWM3lb7+lvHcKNVQp1mUcpD6kRjbNSqrlD5GbImMcYqIC1Lp5410O2oOMrWNESKhkHZs51pk0rhuaUCyq1W2Quo44vpbOo/ScSbz74el+Uwd19CLpAakzs8VEuy+HnZmJL9rMAwrXzNV8TpGHJDSh+SzZ/TV556vc5bURcw/eZWmx9B5VMkUpePk7qqlb/l9cK9pTByQWbsEdlyQ9+VnHBAhhJAOhgsQIYSQJHABIoQQkoTZpwHpmtX+ZkxckMwNp5KTeQMH5mG0eeV9A4Z3Ky5I6Ucy5sBpk8cx4oKsUg31t0v5bSpvnGtfhuicLywoeUtsK/u4q4XIUt9WtvnQ4+TFPIk2I+4HkCW5G8/vpo8TU0ojcC+tZ0QSEwfk6S0BEUhVITDyKsopuc+trA4tUzu6OQtD4S3Gc1CWuodVKiSQX6+c/3X288TVO44TF6S0MhkX5JVxsXPMebpUqFRDKOYx9HcB34AIIYQkgQsQIYSQJHSuCS6rovYe7Fa/DKXfaBeWW7ZRqgEozi1bmmLkuKVK1WmT9ol8t2yrVEO99rJnx7H39dywG7cEqLZQJYrMMb8o91Y0jk7LX/9zvW3LdBYy33npdQyTm+obTMUjL1R+36Kqnkq35UxVPVG5bKb7hubgpeIRbYZJTqXtCbllV41n3DLHS3duw1Ss3Mal+UuaI11TpTwd5ZY9/QezVMPUH+p/bmTbPa4s1dAAfAMihBCSBC5AhBBCksAFiBBCSBI6VwNqFMstu+qvr8HUPFXP9zJwYKOMrcQTciJKKgC+dqMMyoZbdnBcd0r2nJRW4JStVml6lF6Rf+7S49zVlqri7igTt3JDdT7HlHkIuM16Ls8Rmo9sD6bMcdPrGJqP6ivb5DNh7BvUfJp2w7Ya9cCZlYpH7lmq/xnQbthZdym3LSQMeqXZlTCYPyl5/WWkR8n5xZVapXRBz/QXZHocJZLmhx6E0vZEYaXqsbSkHPgGRAghJAlcgAghhCSBCxAhhJAkzA4NKMvXdToiLsgo1QD4dncVExRIf+6XYxC7WnFBocAZV5+QcRuBktyuTVyn4hHbXlr4LLdtatuZrugrNSElDjg2fjknM6W/ofnI7VDaIRVPZewrS5l78Uah9DpWHFBAE/LmpASWfNEnWCbBbVPpXOxaDu59V6U0LLFGxtEIDaXa43yWv3RBDShnglNH8je9+2zH0Ln6pC4fka+NTU0jP2WO/B30rrn8/qr8QPnCmrqXMtanZHyZGiD6DehnP/sZPvWpT2Hx4sXo7+/HaaedhmeffXZ6ClmG66+/HkuXLkV/fz/WrFmDV155JXpihBBCjmyiFqBf/vKXOP/889HT04N/+Zd/wUsvvYS/+Iu/wAc+8IFan2984xu49dZbcfvtt+Opp57CvHnzsHbtWhw+fLjwyRNCCJm9RJngvv71r2P58uW46667an8bGhqqfc6yDLfccgu+/OUv4+KLLwYAfPe738Xg4CAefPBBXHbZZY0frJpNv/fHmNm8131hy4hxy56pTNmhKqcl41Va2kWMSqXaLTvfjVyZ3FRqIcespkxuYttzLxau1cbrvPJAFa/30iTnJQ2PccOWXU037Nis4e6+jfdV5jnLZdtys66zr3uvdQZoaf5FLoYnsjbXBbJJe5mdAi7C7r2VrtWuyU1uSxOcekbkfVfpg9y++WYpmXZIbVesNvs6uSZG9R213LJleizDtBedisf97fDm0AY37B/+8IdYtWoVPvnJT+K4447DGWecgTvvvLPW/tprr2F4eBhr1qyp/W1gYADnnHMOdu3aFXMoQgghRzhRC9Crr76K2267DSeddBIeffRRfP7zn8dVV12Fu+++GwAwPDwMABgcHPT2GxwcrLVJxsbGMDo66v0jhBBy5BNlgqtWq1i1ahVuuukmAMAZZ5yBF198Ebfffjs2bNjQ1AS2b9+Or371q03tSwghZPYStQAtXboUH/nIR7y/nXzyyfinf/onAMCSJUsAACMjI1i6dGmtz8jICE4//fS6Y27duhVbtmypbY+OjmL58uX5k5A+wgW6ZbuVQZWGYmGVapgabPpjTKmGqUlN72uUapDtssqp6ZatdJx8l+2p9vxxreqpMqWMkrQ813CYSE3IO11po495JIy0/EGXc0sDUtqf6JsZbTFu2KFyDNb5mNplfhNgSop6KKvch9RqpOTgujF3+41SA6r0Om294nkJaEDlifz0QPK6VR3tpjxpH8fTwwIpfqRbeckox6BCP9xxZFsLqXnMiqnBWhqaKBPc+eefjz179nh/e/nll3HCCScAmHJIWLJkCXbs2FFrHx0dxVNPPYXVq1fXHbOvrw8LFy70/hFCCDnyiXoDuuaaa3Deeefhpptuwu///u/j6aefxh133IE77rgDwNRbw+bNm3HDDTfgpJNOwtDQELZt24Zly5bhkksuacf8CSGEzFKiFqCzzjoLDzzwALZu3Yo/+7M/w9DQEG655RasX7++1ufaa6/FwYMHccUVV2Dfvn244IIL8Mgjj2DOnDmFT54QQsjspZQpsSAto6OjGBgYwIULP4XuUq/uoGIKIjQg5dOev69Zthbw44JkW5cMMDJ88mVfGW/k9u+ScUyir2MzVm3dYl+3rxxX2Nal3bfq7FvtyR93aixnP2mzl6lT3FT1si0UX2GkOFG2dUMUUuUMjHQ6Uv+yym5bqXdku2qblHPKT7mkj5OvPSn9KKQbWrjXXz3j8j6L7d7pnSfn+G2Tc/ybOTnX/ez3nZjrbaLS75T3ED8lqqS1vI4T05+7Dvudu0RMfbez3TXmX8OuMYjt6fauCVG+Y1xsT1TFdpbbVpoUfd1t0VaSJbrddhn/WJGCl9jXbXeen8nqGP7XT/8a+/fvN2UVJiMlhBCSBC5AhBBCktCx2bCzLEP2nr3AdIm23LJjXbJdN8dQ9dQGx5mahztQvmslEHDLNjJly2bVpkozuiZE2+1azdkw46gLZVVPFcctO32rgczZUSiX1cazPnummYA7rnatzoy2/H1DLtuWGU09e5ZbuTyONN/FGOfdaxEwgVrpdpS5V/xCVXqn2yeFrOya3ABgcq7zfe61T6Y0qey0tY9l6RouUwA529LybbqRW1nagTqSgZv2xt43piKql1orNhVPjP99HfgGRAghJAlcgAghhCSBCxAhhJAkdKwGlIvSVwxtIJQiR+WRiViPXU3FLNUAeAbaYJ4So13pLdL+b6T5kJUaPXtywP1WlU1wPgb0IrN6quqbP26oHIMp0oVs7W6TqbfYfW1dR7YZeoulJYl5hOag9nW3Q5qPqWFJ7dLREdQwUnPw2123bJVOp09s909/djUeAJicJ56Ruc5D1GM/49m4/x0uTU5vyxCAsgoFccYxNB/ZV+k2If3IO45RUgHwfxeDFVGdzzGu9/A1O0/ba1Cz5RsQIYSQJHABIoQQkgQuQIQQQpIwKzSgqDIJblxQjKYjkfE5RcUFGaUagDpxQW4ak0A5b78cA3LbVLtVrrvecdxLHCj9HVO+29OLDH2o7r4GKu1Ks+UZDI2n3rat6+TvG4ov8oiI+wHgl22X41qpeQLlukvOgavd9vdOpubxYnv6RSqeeSLdzjynbYHQfBZMets9c6fz6ZSEzlmZ9IWditBJvFifUMxNOb/Nis8JxwGhYazyDPK7pPSjnP0AaI1dCVOtZXLjGxAhhJAkcAEihBCSBC5AhBBCkjArNCCTdsUFNRsTBNhxQdK+H1GSW49lCAuq3nV+XJAqwR2IC/J0HdFXpuYrefEI+VrSexPJn5OKezDa1bj+porNMIjKhzZTZPk6jllWW/YPlGPwSjkE9Ecf/wK75RYAYGKevz2+YPoGjS/0b9b4gH+c8UXOzR2Y8NrmL/TrJMzrG58+ZsU/5sFDfoBRpcvXhKycbVYFa13qw8iVGOorMct5R46Vs6/SSyP2tcrO5ME3IEIIIUngAkQIISQJnWuCq1anbTSOSUsWcI0q1RBjVrPKPMB/RVczkCY58XrvDxQwmRglIiyziOl2LdpVm0jRotw2jYqcllu2Gkel8WnBDdszSUhXUtHXnXOEq2tHmuNCKJfu/GdEV081THByXOcaTyz08+kcXuz/zBxeJExwA86+wuQ2MeA/JD2LpkuMfvADB7y2wbn+dtm5Yb84NN+fw7jI+dNCtQ9rX8tcp/vaZnMzFY9h6lNmNev7IAmVY6AbNiGEkNkIFyBCCCFJ4AJECCEkCZ2rAbnElD7w0t4EDLARpbNbwnSlFttdlo01YH91zydQjiHGZbtU8SeZuUZkq8wD4JX31rqN2DbKdyvtT+X0d+cnxjU8YdWdacGb1STgyutVNpZpn6RWZhEzX5XGx0jbU7EnMXLOB2qfDx/rn1xVlMOeFKWzK/Onxy4v8F2rjxk46G0PLXq79nnlghGvbaD7XW/7p4eOrX3+5dhcr02WiLBSJQW1vyafkZA+FNR5XOTPonW7LF0nVJrFwC3NoK5vDnwDIoQQkgQuQIQQQpLABYgQQkgSZocGZNC2uCCl2+THBQVLNVgaVoRfvSrVYMUFWaWYRXM4ZqiU267S/Svtxi3Zmz/dqT/kjxuKC/KqAQfKS3itEZqP0m1kX0vnUWZ3mSLfjYGyx/VTp8jrhMKw0vpU5/hxNGWnEsLhQV+AKB8z7m3P6fe3B+Yeqn1eOm/Uaztl4Vve9oUL/qv2+aTud7y2/334Q972a4c+WPs8URHlF0RqHld/BOB/3UOlNCLwdR1DnKyHl+omMAlP17FLpjQ8Tr2xvOCk+AvDNyBCCCFJ4AJECCEkCZ1rgssy1F5RvUzHMZmnIzJly/6hvkURSnFiVESV2aQ9c4yVKRu+a2+m8umI1EGGW3YmzZiq6mlmtEH8wR1Hunf7XbOyvBaueUL0NUxlyixozSlkchPN7jOkzJyWG7Y0z6nUKY65Tjyn2spsuOMrc4rY7HbNzMKENc83wfW8O71z/5t+33e7/b7ZHN8Et6B3Or3OCXPf9tpO6f+/3vavOGa3n072e20vH17qbQ8fWlD7/M5Yr9c2OSGecRlOUK3/GUChJrmiiHLZNiqkagmjSXMds2ETQgjpZLgAEUIISQIXIEIIIUnoXA2oSVy3bKs0AACgXJBbtlGqARBSQUjDikl3bvWNqWCp0umI1Dtijp5bdsDl2T2uTM9humUb5SMAofkAfvr5wL6uABOTLj+k+aihjPQ6Sm8x3FkzWc3Dnb+chKxQK0+wnLuhdDWvrcfvOzknv8SIyIiDnl9KTchPi/Mz57hl8VAcqvrazX++u2K6reJrSz87vMjb3vvutAZ0SGhA1QnxTE/KUiGOLiK/H1b4QIvlCfwD5TeFKpf6IQAR+pAaR/42WOIl3bAJIYTMErgAEUIISQIXIEIIIUmYHRqQl/tF6jr5mkpUmh7AD6IQuk5HxAUFSnKXvJIERqkGMVYwFY8q0e3GG4k5SR3BK8eQ3zbV7m6IcWVckFGuwUodFIupEUWFTNid/dtcym0DgLJRTl0+JFWxt/sd0NqY2HSe+WqPmFOf2O5yNRN/nK7Dft/KO74mdLBnOp7nDX9XvDPe523P7ZmOIZL61jsTvs4zemhO7fP4mPipGw9oQO45BOKAvO1AjFArGpGXbqrpUergXoqI8gtFH5oQQgiZMbgAEUIIScLsMMG1i5hUPWbW5/xM2YDvgmtmygYCFV8jXLRDbtjuYWOqp4qxdGVSwy1bmfry3VutTNlAnWzZXjVG0TfCI920bQQrWOZva+9V/w9l5wSVyc0yyRlesQBQki70XpbnQNZwt/Btr9822SdMfcYviazoWh4XYx2aNsm92+Wb3GTW6r4ev91lfNI37Y2PTbtpV8f9ttKkNMHlz1maFNV2B6TiiXluTRdulf06kCrM/c20Hvgc+AZECCEkCVyACCGEJIELECGEkCR0rgZUrU4bW107dsBubVUfbVv11NC+MdVTLQLVR91roaqnWobqiOqpgO+Wrd29xdiuDdkq1SDa1a0xKq0CIn2+6iuOa6QHUljpdIy+Uzs4n+Wc5HVy753sK+zwZVfvks+TrKShygy4Y9n5XFzX6mq3cKUWmlDVyYqjUgcJlIYy4RxnzN95TMzJK6Mgnp+qrHrquloLt2sYqXfkHEPlGPxUPIG+M4VVebWwvmgtxgF8AyKEEJIILkCEEEKSwAWIEEJIEjpXA2oUSxOKirFBXElur69oa9AHvv64MSXH5XEj4oK8cQKBM1ZcUCWgq7ljB+zHrh6jYncCsRieFqI0LRHfggiM9EBBTcg4ppy/FzqmHgExf1cGkbpaMGbFKEVhnJ/WgPy+bpxQ1a+SUOf7IefozEnE51T96t3IuvO/W5nQceCUXJCpdsoy7kdeJ3esYHoddxK50ysW+ftUSSU2tQbfgAghhCSBCxAhhJAkcAEihBCShNmhARmxPa3Qtrigosp3A/65dwUCLAydR8UFdbk2bjtZmhUXpNuschKhcgxOfI663HbOuZI3p0A5hghzuVcpW2kkdhyNFWJklnWW90rpE25fo4yAGFceN6QBVZ3HLRO/FLI8g6sJybxwctvKmafuTUXeS+N7JzUgR/exSm5Pbfu7enFAKg+hOK7xPHVMXJCL8V0KluBWD1RO/jfmgiOEENLJcAEihBCShNlhgnMJuSk3WT1VHydQqiFF9dRQ2gv39LrkfCNqEsS4Zat0OoZZ0yrVIOYkTUtyX2X+8kxY4jDSzGmkTrGusDymtAYp65CbQSqU8cdz5bXNala1zlCV0xgTnGumrQrrr3S1drer0kW7R5p/A8+X11lsu8+FNIVJc6RjZot1V3fHbikVj0GRZRyU6QzG72ALdkB1HKbiIYQQMhvhAkQIISQJXIAIIYQkoWM1oCwDsvdslaZ7dKHHdDQIK6UMEKfzeHbSxst3A8IsH9KwjPLXEtfVV2kXEW7ZVqkDAMjcOUsxRpX7dUsSBFyRlb3f6Rs8n4KeJ6WZiG13HlIvsvSXkERi6F1KP4pxOTc0LXluUhPKHA2o0ivunXTDlr867nFDcoWllSnX6pLRhty+gK2z6fuT5bYF72W7MFI7KWLKMRQM34AIIYQkgQsQIYSQJHABIoQQkoSO1YAaxtJFYsp3i31bStPTSlyQindpcA6Af76qNHZ+X1XaQMX9GIZsq1SDaFfZi8xxRZM4d3V/rDggQxeRaW/s/Dl2V6WhdDXWJseW45q6QQuaQ0w5BqVvdeXrPNXukAYUiAuycCdtxO4AQucJ9TV0kqhUPKk0n6JQ39GYfY2HuMHDRXHzzTejVCph8+bNtb8dPnwYGzduxOLFizF//nysW7cOIyMjrRyGEELIEUjTC9AzzzyDv/mbv8FHP/pR7+/XXHMNHnroIdx///3YuXMn3nzzTVx66aUtT5QQQsiRRVMmuHfeeQfr16/HnXfeiRtuuKH29/379+Nv//Zvce+99+I3fuM3AAB33XUXTj75ZDz55JM499xzm5qk6R4dN5C/bZnkZqp6qrIFRGTZjqr4aqTmCVQqtbatTNmyPZPztdyypZu1cMuWbrRumpxgVVB3O2AyaaXapWnCkhmju6xnROxrZEJqBeW+bmUClyY4L3N2wMSmjhNhJ/TMp7Ybtmf+DbjxWyY5bdLN3w7eD+95isgDlYqYbNhNpOVp6g1o48aNuOiii7BmzRrv77t378bExIT395UrV2LFihXYtWtXM4cihBByhBL9BnTffffhueeewzPPPKPahoeH0dvbi0WLFnl/HxwcxPDwcN3xxsbGMDY2VtseHR2NnRIhhJBZSNQb0BtvvIGrr74a99xzD+bMmVPIBLZv346BgYHav+XLlxcyLiGEkM4m6g1o9+7d2Lt3Lz72sY/V/lapVPCjH/0I3/72t/Hoo49ifHwc+/bt896CRkZGsGTJkrpjbt26FVu2bKltj46OmotQ0D16tlVPDe3ruHS3VD1VpRJyjynOLcYt2yrVINuVruN3tccVwwrdzXOnDrnYWjqClZZf2bgb1yPlvasabszhMg/5tvaoLEMhN2x3Wx5TzdG9qGJcY76qf8it3CjHoLSaitHWgk5opwcK9I2hTZqQWbohhHwWK3V7NUzUAnThhRfixz/+sfe3yy+/HCtXrsSXvvQlLF++HD09PdixYwfWrVsHANizZw9ef/11rF69uu6YfX196Ovra3L6hBBCZitRC9CCBQtw6qmnen+bN28eFi9eXPv7Zz/7WWzZsgXHHHMMFi5ciCuvvBKrV69u2gOOEELIkUnhmRC++c1volwuY926dRgbG8PatWvxne98p+jDEEIImeWUMiluJGZ0dBQDAwP4f/p+H92lnmD/qLigkCZkjWXsq+ZgxQWp2Bepmcjtcm5fdVx32yrVILe77L6ZbHfHluOKvpnbt1u0deUfJ1N9xXa3v2/VLR/dE+jbM71d6ZFt3qbXLktNV3rFvqrdnYNsExqQW9JaxdFAbOfHDGVSb7G+HqGvTsnQddS20VeidJKIMtuTpbqfAaA8Ifvmt5Un/MOUJ5HbrvpO+CdQnjT6Thp9VZvQYsWcSk65+3JF9hVjVaq5bVIPLk86F13E6pUqcl+j3WmbrIxhx56/xP79+7Fw4ULkwWSkhBBCksAFiBBCSBI6Nxt2NZt+rTdMWlHu0SGs9DQxFFY9FfBsEh1QPRWo45btNcos1VYqHt+2VPJctgM5coxs2SrDdUyaFVkZs8ttk33t7bKzXZGXX5iWvOc4lA7bS/Ej2yLS3sSY1WSTFYVgmdiAAqucijnF3OeIbTP7tdyOcdmOxLsd1fzv2WyCb0CEEEKSwAWIEEJIErgAEUIISULnakAuEaUPzNINUeUL5ByOouqpXdJmHyjPYMzJTMUjjOuuy7Z0/1S6k7T/O+1WqQbZLh8Bqau5rrFyHOk2K++dO5Z61ISLbWZ4tmtX6/y2IEZ6HXXcqHKqRlvAtdormyBLKigNCLl9tSbXuF5kVkxV6YDyt81yHnI7lHaoFV3HrPDaOXoR34AIIYQkgQsQIYSQJHABIoQQkoTZoQG1C0sT6oTy3YD4L8IMle+OiAtSpQJaKd/t6D5S8ylVAufu7CtLNci0Ja4EEcqM5N73bFK2+ZsyDU7ZfWasUgcAuoyaBJnaefpjVcXYyDQ+Rtln2WScu8LSLwyNp96+nlYzGdB1Jg1dR6bmce5XTMwWIOOAZMyN3zcuDsj9PuCoh29AhBBCksAFiBBCSBK4ABFCCEnC7NOAIvKsFZonLgLzuGr+Leg6UhOaifLdYh5R5buVhiWO6+4r080rDS4/wKIs9KKquKZe1YpgPrTpjyGNRIfGuHO0dRDvdMRAUr7zyoTLmCfxja6KUhRu7I/Sh0LXwsJNMac0INFXxoNFxPa4Wo2l+Uy119+v3rjWdihmKEYvsuKAwjnnOlw/8jTPxh4evgERQghJAhcgQgghSZh9JjhJkSY5y7XacssOuDib6YHUJI1UPTFpemIIpeaQpoGIcgz+tZGmMJGKx2kviYPKy1ISY3mmPyNNj5xFVdwPs8hs0O4hxnLNUqpMhdFXlYiQ245ruKyeKt2YhZXTfYakCU6bGHM+1yPCtGSZ6ELnXrbS68iqpp5pzx5XzzE/tVNUKh7DrKZS4gRT8yAfs60T7XVT8A2IEEJIErgAEUIISQIXIEIIIUnoXA0oq6JmHLZck9tFKHWNVb47Ju1NTPnukHt0u8p3Swybsirf7Q4VKOvgle820vQAQCY0Is9Ob6TpmWp3Pk7645TN/5PJexOwrbuiitIGpFZmaBviW+pqG1VROiMTmk9ZbLu6jywvodIqxWhA3kBimIB7se/GHCqxkD+OcrU2UvHovllue7h8d35Ygn7G3TYxjtVXEpJ1Olj3ceEbECGEkCRwASKEEJIELkCEEEKS0LkakIsbBBLSg5os3w3MTKqecCySkaondO7e2Hbpb68EtBwnpAl519je143X0WW18+sGqzifQBofL9ZnMn/cKcp1Pk1RReOaUEz5ZakNVFUsjKPNGHE/gF8WQmo8VfGNlrqOq/sozUf2bVYDkoTiWaw4IEtPko+plU5HteVrPqF9ZXmPkjV/Kw4oIvUOAO8Z0XpRhOYjf2PcfUNfHQt3nAbnwzcgQgghSeACRAghJAmzwwTXLDEuzgikzGlT9dQoM6CVpgeIS9VjVDUNjmK9sluZs4XJLesyTHsqq7ZMmWOl8RHI58CbtD8HyyQnzXHKymBsK/dimTXZaVcmOOEi7Ga8lm7Y5Ql/36pKxeOMo8xzoq+ZCRyNEzLBRaXxaawNgHebLbMZYLt7h92w3c+BbNiey7adiidk4vWwTJeteGS32Z2bb0CEEEKSwAWIEEJIErgAEUIIScLs04BUXv5Ea2ibSjfo4zTuVm66RysjccR1i0otJLabLt0ghpG7ytIOzoGtNgBAt1eQQYwsNCHPs1RUWhXCSNnQL6R7tKUjyHQ6JZlux9EnyuKgKhWPofPEpOJpSQMKYUh/dtqegIZiVI4N6TquZhSTiieUcsnbN+iebpxfqkw7rbhp14FvQIQQQpLABYgQQkgSuAARQghJwuzTgCQxmlCR5btjsEo3xBzXStMD+OfertINgK1/SazSDVZckNxP9LU0IVXO29KEQrqaJ37k60NAnWoThl5R7ZZxQc5+Uq+QJRXcy6TabF3H14CM0hnAzGlA3iTEYazYKiu8C0IDCsYXie+dGQeUr+uoeCOlPUXoRUa5Bl3OO387WObBGrfN8A2IEEJIErgAEUIIScLsN8G1QlEmuVjXam8O9r52eqCZz5wNIK6aqukaLsZ10wOpVMeyymm+SS7khm2m7ZGZtB2zYEmYezNhRssyaf6a3lYu2+IwrilNuWFXjL4yU3bArOaZ4ErFuWGr9iYx080Azbtsy+9vK+l1rIqolskNEJm/87NdAzBTO5mpd1TfgLnOItS3RZMd34AIIYQkgQsQIYSQJHABIoQQkoSO1YCyalaroFmKKTMQUz1VEpP2xhynoNINYt8ZK90Q2DeqmqqnlYm+RjXVYEVUQxMKpu1x/Zyl8CHS3ri29awrIHzI8hKOVqNctqVW45VjEOOKvq4Lt6Xb1Gv3KqKWbL0oTgMynqcWSjeoobyqm7Itf6yQG7b8Xtou3JbLdr7mI9vDFXWl+7SXF8rv3FKF1Ma7Fu2mzTcgQgghSeACRAghJAlcgAghhCShYzUgl8y1mzarB03tXNCMAvE5kmZLN4T2VccxNCxL11G2/8avsdKldIf8nS1NqCrjfCI0oUDaHlfQkONYxylJIURdY3/TK7MtU+TIEgtuKh5VYsGIIZJ9ZWyPkcYnFMvj9w2kkHIFjHal6VEH9TftFDNx+5opc2Ssj1Gi3kqvU1Jlwu1tGLqU0s4Mrcn8Trai8bjzlXPPgW9AhBBCksAFiBBCSBK4ABFCCEnCrNCAXDJpJ22XJpSqdEMEcaUbjDnFlG4AhH7Ugv5lHcaIEQICmlDo+jtzVuPI47h6UVnGE9lp7j3tRgksIt7ImbOlD03NI1/rU3E/Qq9wNSMVuyOn6OacCz3SVsyQ0TcaQ1ow88gFYmwsnSfY1yuzLWOEDA0olPvNLP0dEfcj21QeOSO+KBQjxFxwhBBCZiNcgAghhCRh1pngkjFTpRsst+yZKt3QrmqqMWUrItL2AMJUJlPiyOMaNiJ53bzSDdKelTVuvpPPizQde27ahnkOECY6lU5H9JWn6qbikY+etMi5pyuPEzDf+XPKb2snvmlMtAXLDBj7WpVLQ6UbrFRCFWmSM7YNt2s1pyJLKjAVDyGEkCMBLkCEEEKSwAWIEEJIEma9BjRjbtmSmSjdAPg21yJLN3jzD1wHUxMqsHRDK5qQQ6kiz0dqKPlpe2Rfb0uW3Jaaj1WuQbZVxbbrGm7pQ7BT/JQCqXlccSBKL5LDWJpPRN8ZI5S2R/U3UtmYeos4juqbn7ZHl4TId9OWz7iZbsfSoSQxfQH/Wlju3DnwDYgQQkgSuAARQghJAhcgQgghSZj1GpCk6dINUztPf+6E0g1A8+W8jRih4DxitDFl683XhMwYIQCouPlPpEYSClpx9Rb73rn2c6mDmIKF1FuMmKGp/q6OIHUdI62PoQ+pbdnX0HGm2l2tKUIvMtoklpYUIpjGJwJT5wmUCzB1HTPtjehrxPYEyySosdzSIKF93eNExBdJYvo2QdSv7Pbt23HWWWdhwYIFOO6443DJJZdgz549Xp/Dhw9j48aNWLx4MebPn49169ZhZGSk0EkTQgiZ/UQtQDt37sTGjRvx5JNP4rHHHsPExAR+8zd/EwcPHqz1ueaaa/DQQw/h/vvvx86dO/Hmm2/i0ksvLXzihBBCZjelTNpmIvj5z3+O4447Djt37sSv/dqvYf/+/fjgBz+Ie++9F7/3e78HAPjJT36Ck08+Gbt27cK5554bHHN0dBQDAwP4eOkSdJd63ptlMeawaJOcv3PjfSOOE5U5O5RN2jQf5e+r5hCav3stgn3dtDHSRNV4SpygSc47TuB8jL4qjY81bug4rlu86utvumNJ85xZsVbOP2D+8uYRcsPOO2adOanz8/aN6BvC2rUF61Ao47XZ18qGbZnKQqa9mPQ6xpxKKhWV0TfWXOe6gzttk5Ux7Pg/f4X9+/dj4cKFyKOlX/b9+/cDAI455hgAwO7duzExMYE1a9bU+qxcuRIrVqzArl27WjkUIYSQI4ymnRCq1So2b96M888/H6eeeioAYHh4GL29vVi0aJHXd3BwEMPDw3XHGRsbw9jYWG17dHS02SkRQgiZRTT9BrRx40a8+OKLuO+++1qawPbt2zEwMFD7t3z58pbGI4QQMjto6g1o06ZNePjhh/GjH/0Ixx9/fO3vS5Yswfj4OPbt2+e9BY2MjGDJkiV1x9q6dSu2bNlS2x4dHZ39i1CqaqpFpe1RaW+kET+/HIPu6+WqF3PK31fNSUwpRRqfoMu2FFHckACpLUm3bK/Cgij7oFyt3RQ5MhWP31Xu61VelcNaepEUM9ScLNFEpjcKlftIgFX5MyCTx5U+cPeL0Hzk2AVWOTX7BubYKlG3PssybNq0CQ888AAef/xxDA0Nee1nnnkmenp6sGPHjtrf9uzZg9dffx2rV6+uO2ZfXx8WLlzo/SOEEHLkE/UGtHHjRtx77734wQ9+gAULFtR0nYGBAfT392NgYACf/exnsWXLFhxzzDFYuHAhrrzySqxevbohDzhCCCFHD1EL0G233QYA+PjHP+79/a677sKnP/1pAMA3v/lNlMtlrFu3DmNjY1i7di2+853vFDJZQgghRw4txQG1g7pxQC4FpsiZsbggl8hjmppQKC7IH6jpcaLihKzrEtJQIvY144SsGCFJCzFDKgYnZt+IeCMzhig0p5h9JYZeZOqaEeW6G5pH7nEa7xssv2AeqAVdR+LqghHpdKb2NY5jbYdKLJh6UUyZcBEH9Oqt7Y0DIoQQQpqFCxAhhJAkzL5s2K1UMVVDFVRNtdlKqkDQJGdm0o5xRbYyaQfGMd20LRdtwL82ZmVVwPSFNSqtAsLKI89HUpDLtkS6r8s5egQqsXpN8hpXGzfXmW7ZQXOd8+zJSRmu4dLNOmgqi8qW7R6nfdhVQ/PbYiqtBk17LblL5+zXyHFnEL4BEUIISQIXIEIIIUngAkQIISQJs08DkrRJE2paD4qdQ5Fpe4yKqHUGcwfKH6fOWJ4uFZBbUDa0MksTUmlu7DQ+nt5S5HWqGNcp4FbuVUFtxWVbVtU0SiqoczVS/kBURDV1HqVRGdVUQ5Vv1f0wrrHctZWaC0URo6EE0uC4RGk1raTXscaKcbsW2yXPxbyxa8Q3IEIIIUngAkQIISQJXIAIIYQkYfZrQG2isBihqZ0LmFEkRcUIBcYK61JOezlwXby+eZOtdRab+TskixkydLaSqoUQU/bBte8H+loakZy/Wd7bjuGKGdfUj0IaTyf8lzmkezqYWkiMNhM4blR6HUnV6BvYLrVYnqETbichhJCjEC5AhBBCksAFiBBCSBKOPA2o2RxtwWFnKG9cu8p5t0sTCuWN87QBeW6N536zS30DTeeRayVmSNnHrRIRYk7GNbb0IUBoRKF7FdKIXEIxRda4jnajcsqF9KNG5wf4t7mV8vVF0kpckEtAWyostsfq20reuJj5vQffgAghhCSBCxAhhJAkHHkmOJdOcIeWxM7Jc01uoXSDGrfJtD1Tg9cfp85YdjkJw9wYSL2jaCGNj9fVMg9VKqJRujwbc7TS6QD+/ANmNSvFj3JxDrl0e52bN9f5uwWePaM+Q0yVU5lKqCMoykUbiEvFo+bRpLku0Fe5XeftSxMcIYSQToYLECGEkCRwASKEEJKEI1sDknRC6QY9kL9tzakTXLSnBncHbnisuLQ9EaW+1b5ywpLmS397U1D7GuO2VPZB9LWuv9SLDO1J6UMhl26vr9EW0oBa0Y8anUOnUJSLNmC7UxeU8ieYWqfg8t58AyKEEJIELkCEEEKSwAWIEEJIEjpXA8oy1GIa2pVyoyBNqKU0PXowf7tNmpB3iMiS3GJgOVj+WIG0Pd4wMaW+JYlihsy0Pq2UfTBjiALPgFHKQaX8kTSrF4WeCXk+Vl+LVr5nnUBsKYNmdZ7Y2B5r30aPwzggQgghnQwXIEIIIUnoXBOci+V2Wuhxismk3TaTXIFpe7xDxGSEloQyRDdphopy2QbsND4z5LItcU10wbvR7DVW10XNQowVk0m7OXOdmkHoWbTMd1Zfy5Rn7TeTtOKmXJALd7Rrdcwc8rJuN2he5BsQIYSQJHABIoQQkgQuQIQQQpIwOzQgl5h0IS0d5+hJ26MOU2QaH8uGXFSlVcC2ORflsq13tvd1xs4C7tLm3SlKgwPslD8F6UUx6YDqHtdtCVbCtcZtcL9OoYU5RrlSx7S3kh6oAfgGRAghJAlcgAghhCSBCxAhhJAkzD4NSDLLNKGOSdvjElHqGyiwtENBpb4l7YsZitAjpnY2mloo++CN00KKnxDm0AXpRfXa3aagfmTMyerZAWl8gvE5FkWWeYjp26jWxFQ8hBBCOhkuQIQQQpIw+01wqZjtaXtcIl22XfNXsIKl4VpdZ+DpzzHpgMTY6Vy2G3dNjsq6HWHuUudaqUB0yD1OVMofSbPmOnkcSch8Z/U1UKa9GDrB3btdJrii0vbQBEcIIaST4QJECCEkCVyACCGEJOHI04BmqnSDd8z2pO2ZGqrJc2hlTi1UWo1y0Za4GkRMWn459ky5bEssjaiVsg8Rx7HcuYEWXbpdWtGLWnGttr4PBelDQTohrU+Ruk6z4zbTX8A3IEIIIUngAkQIISQJXIAIIYQk4cjTgFxmKk2POu4sL+UgaUET8g7ZrpihqcHrj1MPI2ZIUrKGUvEsbSr7EIohijhO0yl/1HFa0ItawSwnYdF8eZKOocXSBzXaFkPkPBMsyU0IIaST4QJECCEkCVyACCGEJOHI1oAknaAJdUKM0NRg7kBx+1olCqxDtqusw9Tg+W0RMUN6WCPvXcjO3WzZBzWOcYwYfQhoPuecPKx9FHHMRHqRRSgX3Ez9NsRQVPxRK+PE3MsG4BsQIYSQJHABIoQQkoSjywQnSWGSa2PaHn/YiHMp0kVbkiKNjz6Qv90ul+2iyj4UZZ6rSzHHDaX88YZtvGtrJp4izXcpUnrFUGQ6oFaued48WI6BEEJIJ8MFiBBCSBK4ABFCCEnC0a0BSWZ5KQd/2BZctkOuvM1qRBGlvtUh25nGxz9Q/riSIss+eOMWU5phaqyIktYFlYSQtE0vkhTsIlxjplzDY2jXubagLbnPf0gvfZ8OvLKEEEKOBrgAEUIISQIXIEIIIUmgBpRHJ6TtATpDE9KDuQM1vl9EWQd1yE6IGWrhOFExRK2k+JFEaDVx/x1NrxdJ2vYNbZfe0k4KihNqVMtplra9Af31X/81TjzxRMyZMwfnnHMOnn766XYdihBCyCykLQvQP/7jP2LLli34yle+gueeew6/+qu/irVr12Lv3r3tOBwhhJBZSClrwzvWOeecg7POOgvf/va3AQDVahXLly/HlVdeieuuu87cd3R0FAMDA/g4LkZ3qafoqRVHivQcBZnj9LAFnktRc2xhTkGXbfO4MT7DEcdpwZU36nxirlvMvWrlGYmZfyvPT5uqnLb0PHUAbTWj5ZiLJ6vj2PGLv8X+/fuxcOHC3N0L/0UbHx/H7t27sWbNmumDlMtYs2YNdu3apfqPjY1hdHTU+0cIIeTIp/AF6Be/+AUqlQoGBwe9vw8ODmJ4eFj13759OwYGBmr/li9fXvSUCCGEdCDJveC2bt2KLVu21Lb379+PFStWYBITkPWxOosUr+VtMsFlRZ5LQXNsYU4tnU2UG1ZMdokWTHBRx4k5+xgPxllggmvTd5ImOAPDBNfIsQtfgI499lh0dXVhZGTE+/vIyAiWLFmi+vf19aGvr6+2/b4J7gn8f0VPrVhSLI7tOuYs9DIlhHQ+Bw4cwMDAQG574QtQb28vzjzzTOzYsQOXXHIJgCknhB07dmDTpk3B/ZctW4Y33ngDWZZhxYoVeOONN0wR62hndHQUy5cv53UKwOvUGLxOjcHrZJNlGQ4cOIBly5aZ/dpigtuyZQs2bNiAVatW4eyzz8Ytt9yCgwcP4vLLLw/uWy6Xcfzxx9fehBYuXMgb3AC8To3B69QYvE6NweuUj/Xm8z5tWYD+4A/+AD//+c9x/fXXY3h4GKeffjoeeeQR5ZhACCHk6KVtTgibNm1qyORGCCHk6KRjk5H29fXhK1/5iuegQDS8To3B69QYvE6NwetUDG3JhEAIIYSE6Ng3IEIIIUc2XIAIIYQkgQsQIYSQJHABIoQQkoSOXYBY0G6a7du346yzzsKCBQtw3HHH4ZJLLsGePXu8PocPH8bGjRuxePFizJ8/H+vWrVPpkI42br75ZpRKJWzevLn2N16nKX72s5/hU5/6FBYvXoz+/n6cdtppePbZZ2vtWZbh+uuvx9KlS9Hf3481a9bglVdeSTjjmadSqWDbtm0YGhpCf38/PvzhD+NrX/ual9+M16lFsg7kvvvuy3p7e7O/+7u/y/7rv/4r+6M/+qNs0aJF2cjISOqpJWHt2rXZXXfdlb344ovZ888/n/3Wb/1WtmLFiuydd96p9fnc5z6XLV++PNuxY0f27LPPZueee2523nnnJZx1Wp5++unsxBNPzD760Y9mV199de3vvE5Z9vbbb2cnnHBC9ulPfzp76qmnsldffTV79NFHs//5n/+p9bn55puzgYGB7MEHH8xeeOGF7Hd+53eyoaGh7NChQwlnPrPceOON2eLFi7OHH344e+2117L7778/mz9/fvZXf/VXtT68Tq3RkQvQ2WefnW3cuLG2XalUsmXLlmXbt29POKvOYe/evRmAbOfOnVmWZdm+ffuynp6e7P7776/1+e///u8MQLZr165U00zGgQMHspNOOil77LHHsl//9V+vLUC8TlN86Utfyi644ILc9mq1mi1ZsiT78z//89rf9u3bl/X19WXf+973ZmKKHcFFF12UfeYzn/H+dumll2br16/PsozXqQg6zgQXW9DuaGT//v0AgGOOOQYAsHv3bkxMTHjXbOXKlVixYsVRec02btyIiy66yLseAK/T+/zwhz/EqlWr8MlPfhLHHXcczjjjDNx555219tdeew3Dw8PedRoYGMA555xzVF2n8847Dzt27MDLL78MAHjhhRfwxBNP4BOf+AQAXqciSF4PSGIVtPvJT36SaFadQ7VaxebNm3H++efj1FNPBQAMDw+jt7cXixYt8vrmFQE8krnvvvvw3HPP4ZlnnlFtvE5TvPrqq7jtttuwZcsW/Mmf/AmeeeYZXHXVVejt7cWGDRtq16LRopJHKtdddx1GR0excuVKdHV1oVKp4MYbb8T69esBgNepADpuASI2GzduxIsvvognnngi9VQ6jjfeeANXX301HnvsMcyZMyf1dDqWarWKVatW4aabbgIAnHHGGXjxxRdx++23Y8OGDYln1zl8//vfxz333IN7770Xp5xyCp5//nls3rwZy5Yt43UqiI4zwcUWtDua2LRpEx5++GH867/+K44//vja35csWYLx8XHs27fP63+0XbPdu3dj7969+NjHPobu7m50d3dj586duPXWW9Hd3Y3BwUFeJwBLly7FRz7yEe9vJ598Ml5//XUAqF2Lo/07+MUvfhHXXXcdLrvsMpx22mn4wz/8Q1xzzTXYvn07AF6nIui4BcgtaPc+7xe0W716dcKZpSPLMmzatAkPPPAAHn/8cQwNDXntZ555Jnp6erxrtmfPHrz++utH1TW78MIL8eMf/xjPP/987d+qVauwfv362mdeJ+D8889Xbvwvv/wyTjjhBADA0NAQlixZ4l2n0dFRPPXUU0fVdXr33XdRLvs/kV1dXahWp0oI8zoVQGoviHrcd999WV9fX/b3f//32UsvvZRdccUV2aJFi7Lh4eHUU0vC5z//+WxgYCD7t3/7t+ytt96q/Xv33XdrfT73uc9lK1asyB5//PHs2WefzVavXp2tXr064aw7A9cLLst4nbJsykW9u7s7u/HGG7NXXnklu+eee7K5c+dm//AP/1Drc/PNN2eLFi3KfvCDH2T/+Z//mV188cVHnXvxhg0bsg996EM1N+x//ud/zo499tjs2muvrfXhdWqNjlyAsizLvvWtb2UrVqzIent7s7PPPjt78sknU08pGQDq/rvrrrtqfQ4dOpR94QtfyD7wgQ9kc+fOzX73d383e+utt9JNukOQCxCv0xQPPfRQduqpp2Z9fX3ZypUrszvuuMNrr1ar2bZt27LBwcGsr68vu/DCC7M9e/Ykmm0aRkdHs6uvvjpbsWJFNmfOnOxXfuVXsj/90z/NxsbGan14nVqD5RgIIYQkoeM0IEIIIUcHXIAIIYQkgQsQIYSQJHABIoQQkgQuQIQQQpLABYgQQkgSuAARQghJAhcgQgghSeACRAghJAlcgAghhCSBCxAhhJAkcAEihBCShP8f21qsljUV7uYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Reading the x_keys above we can input the parameters that we would like the simulator to evaluate\n",
    "x = torch.tensor([\n",
    "    z_l.item(), # sie z_l\n",
    "    0.7,     # sie x0\n",
    "    0.13,    # sie y0\n",
    "    0.4,     # sie q\n",
    "    np.pi/5, # sie phi\n",
    "    1.,      # sie b\n",
    "    0.2,     # src x0\n",
    "    0.5,      # src y0\n",
    "    0.5,      # src q\n",
    "    -np.pi/4, # src phi\n",
    "    1.5,       # src n\n",
    "    2.5,      # src Re\n",
    "    1.,       # src Ie\n",
    "])\n",
    "plt.imshow(sim(x), origin=\"lower\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Where to go next?\n",
    "\n",
    "The caustic tutorials are generally short and to the point, that way you can idenfity what you want and jump right to some useful code that demo's the particular problem you face. Below is a list of caustic tutorials and a quick description of what you will learn in each one::\n",
    "\n",
    "- `LensZoo`: here you can see all the built-in lens mass distributions in `caustic` and how they distort the same background Seric source.\n",
    "- `Playground`: here we demo the main visualizations of a lensing system (deflection angles, convergence, potential, time delay, magnification) in an interactive display so you can change the parameters by hand and see how the visuals change!\n",
    "- `VisualizeCaustics`: here you can see how to find and display caustics, a must when using `caustic`!\n",
    "- `Simulators`: here we describe the powerful simulator framework and how it can be used to quickly swap models, parameters, and other features and turn a complex forward model into a simple function.\n",
    "- `InvertLensEquation`: here we demo forward ray tracing in `caustic` the process of mapping from the source plane to the image plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "PY39",
   "language": "python",
   "name": "py39"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}