diffphys-code-control.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "P9P3fJaa30da"
   },
   "source": [
    "# Solving Inverse Problems with NNs\n",
    "\n",
    "Inverse problems encompass a large class of practical applications in science. In general, the goal here is not to directly compute a physical field like the velocity at a future time (this is the typical scenario for a _forward_ solve), but instead more generically compute one or more parameters in the model equations such that certain constraints are fulfilled. A very common objective is to find the optimal setting for a single parameter given some constraints. E.g., this could be the global diffusion constant for an advection-diffusion model such that it fits measured data as accurately as possible. Inverse problems are encountered for any model parameter adjusted via observations, or the reconstruction of initial conditions, e.g., for particle imaging velocimetry (PIV). More complex cases aim for computing boundary geometries w.r.t. optimal conditions, e.g. to obtain a shape with minimal drag in a fluid flow.\n",
    "\n",
    "A key aspect below will be that we're not aiming for solving only a _single instance_ of an inverse problem, but we'd like to use deep learning to solve a _larger collection_ of inverse problems. Thus, unlike the physics-informed example of {doc}`physicalloss-code` or the differentiable physics (DP) optimization of {doc}`diffphys-code-ns`, where we've solved an optimization problem for specific instances of inverse problems, we now aim for training an NN that learns to solve a larger class of inverse problems, i.e., a whole solution manifold. Nonetheless, we of course need to rely on a certain degree of similarity for these problems, otherwise there's nothing to learn (and the implied assumption of continuity in the solution manifold breaks down).\n",
    "\n",
    "Below we will run a very challenging test case as a representative of these inverse problems: we will aim for computing a high dimensional control function that exerts forces over the full course of an incompressible fluid simulation in order to reach a desired goal state for a passively advected marker in the fluid. This means we only have very indirect constraints to be fulfilled (a single state at the end of a sequence), and a large number of degrees of freedom (the control force function is a space-time function with the same degrees of freedom as the flow field itself).\n",
    "\n",
    "The _long-term_ nature of the control is one of the aspects which makes this a tough inverse problem: any changes to the state of the physical system can lead to large change later on in time, and hence a controller needs to anticipate how the system will behave when it is influenced. This means an NN also needs to learn how the underlying physics evolve and change, and this is exactly where the gradients from the DP training come in to guide the learning task towards solution that can reach the goal.\n",
    "[[run in colab]](https://colab.research.google.com/github/tum-pbs/pbdl-book/blob/main/diffphys-code-control.ipynb)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "gaCwYgN_L-We"
   },
   "source": [
    "## Formulation\n",
    "\n",
    "With the notation from {doc}`overview-equations` this gives the minimization problem \n",
    "\n",
    "$$\n",
    "\\text{arg min}_{\\theta} \\sum_m \\sum_i (f(x_{m,i} ; \\theta)-y^*_{m,i})^2 , \n",
    "$$\n",
    "\n",
    "where $y^*_{m,i}$ denotes the samples of the target state of the marker field, \n",
    "and $x_{m,i}$ denotes the simulated state of the marker density.\n",
    "As before, the index $i$ samples our solution at different spatial locations (typically all grid cells), while the index $m$ here indicates a large collection of different target states.\n",
    "\n",
    "\n",
    "Our goal is to train two networks $\\mathrm{OP}$ and $\\mathrm{CFE}$ with weights\n",
    "$\\theta_{\\mathrm{OP}}$ and $\\theta_{\\mathrm{CFE}}$ such that a sequence \n",
    "\n",
    "$$\n",
    "\\newcommand{\\pde}{\\mathcal{P}}\n",
    "\\newcommand{\\net}{\\mathrm{CFE}}\n",
    "\\mathbf{u}_{n},d_{n} = \\pde(\\net(~\\pde(\\net(\\cdots \\pde(\\net( \\mathbf{u}_0,d_0, d_{OP} ))\\cdots)))) = (\\pde~\\net)^n ( \\mathbf{u}_0,d_0, d_{OP} ) .\n",
    "$$\n",
    "\n",
    "minimizes the loss above. The $\\mathrm{OP}$ network is a predictor that determines the state $d_{OP}$ that the action of the $\\mathrm{CFE}$ should aim for, i.e., it does the longer term planning from which to determine the action. Given the target $d^*$, it computes \n",
    "$d_{OP} = \\mathrm{OP}(d,d^*)= f_{\\mathrm{OP}}(d,d^{*};\\theta_{\\mathrm{OP}})$.\n",
    "The $\\mathrm{CFE}$ acts additively on the velocity field by computing\n",
    "$\\mathbf{u} + f_{\\mathrm{CFE}}(\\mathbf{u},d, f_{\\mathrm{OP}}(d,d^{*};\\theta_{\\mathrm{OP}}) ;\\theta_{\\mathrm{CFE}}) $, where we've used $f_{\\mathrm{OP}}$ and $f_{\\mathrm{CFE}}$ to denote the NN representations of $\\mathrm{OP}$ and $\\mathrm{CFE}$, respectively, and $d^{*}$ to denote the target density state. $\\theta_{\\mathrm{OP}}$ and $\\theta_{\\mathrm{CFE}}$ denote the corresponding network weights.\n",
    "\n",
    "For this problem, the model PDE $\\mathcal{P}$ contains a discretized version of the incompressible Navier-Stokes equations in two dimensions for a velocity $\\mathbf{u}$:\n",
    "\n",
    "$$\\begin{aligned}\n",
    "  \\frac{\\partial u_x}{\\partial{t}} + \\mathbf{u} \\cdot \\nabla u_x &= - \\frac{1}{\\rho} \\nabla p \n",
    "  \\\\\n",
    "  \\frac{\\partial u_y}{\\partial{t}} + \\mathbf{u} \\cdot \\nabla u_y &= - \\frac{1}{\\rho} \\nabla p \n",
    "  \\\\\n",
    "  \\text{s.t.} \\quad \\nabla \\cdot \\mathbf{u} &= 0,\n",
    "\\end{aligned}$$\n",
    "\n",
    "without explicit viscosity, and with an additional transport equation for the marker density $d$ given by \n",
    "$\\frac{\\partial d}{\\partial{t}} + \\mathbf{u} \\cdot \\nabla d = 0$.\n",
    "\n",
    "To summarize, we have a predictor $\\mathrm{OP}$ that gives us a direction, an actor $\\mathrm{CFE}$ that exerts a force on a physical model $\\mathcal{P}$. They all need to play hand in hand to reach a given target after $n$ iterations of the simulation. As apparent from this formulation, it's not a simple inverse problem, especially due to the fact that all three functions are non-linear. This is exactly why the gradients from the DP approach are so important. (The viewpoint above also indicates that _reinforcement learning_ is a potential option. In {doc}`reinflearn-code` we'll compare DP with these alternatives.)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "tFb5WYgzL-Wf"
   },
   "source": [
    "\n",
    "---\n",
    "\n",
    "## Control of incompressible fluids \n",
    "\n",
    "The next sections will walk you through all the necessary steps from data generation to network training using [Φ<sub>Flow</sub>](https://github.com/tum-pbs/PhiFlow). Due to the complexity of the control problem,\n",
    "we'll start with a supervised initialization of the networks, before switching to a more accurate end-to-end training with DP.\n",
    "(_Note: this example uses an older version 1.4.1 of Φ<sub>Flow</sub>._)\n",
    "\n",
    "The code below replicates an inverse problem example (the shape transitions experiment) from [Learning to Control PDEs with Differentiable Physics](https://ge.in.tum.de/publications/2020-iclr-holl/) {cite}`holl2019pdecontrol`, further details can be found in section D.2 of the paper's [appendix](https://openreview.net/pdf?id=HyeSin4FPB).\n",
    "\n",
    "First we need to load phiflow and check out the _PDE-Control_ git repository, which also contains some numpy arrays with initial shapes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "pwVPXx_Y30dd"
   },
   "outputs": [],
   "source": [
    "!pip install --quiet phiflow==1.4.1\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from phi.flow import *\n",
    "\n",
    "if not os.path.isdir('PDE-Control'):\n",
    "  print(\"Cloning, PDE-Control repo, this can take a moment\")\n",
    "  os.system(\"git clone --recursive https://github.com/holl-/PDE-Control.git\")\n",
    "    \n",
    "# now we can load the necessary phiflow libraries and helper functions\n",
    "import sys; sys.path.append('PDE-Control/src')\n",
    "from shape_utils import load_shapes, distribute_random_shape\n",
    "from control.pde.incompressible_flow import IncompressibleFluidPDE\n",
    "from control.control_training import ControlTraining\n",
    "from control.sequences import StaggeredSequence, RefinedSequence\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "vQCiicZv30de"
   },
   "source": [
    "## Data generation\n",
    "\n",
    "Before starting the training, we have to generate a data set to train with, i.e., a set of ground truth time sequences $u^*$. Due to the complexity of the training below, we'll use a staged approach that pre-trains a supervised network as a rough initialization, and then refines it to learn control looking further and further ahead into the future. (This will be realized by training specialized NNs that deal with longer and longer sequences.) \n",
    "\n",
    "First, let's set up a domain and basic parameters of the data generation step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "EthZeAV030df"
   },
   "outputs": [],
   "source": [
    "domain = Domain([64, 64])  # 1D Grid resolution and physical size\n",
    "step_count = 16  # how many solver steps to perform\n",
    "dt = 1.0  # Time increment per solver step\n",
    "example_count = 1000\n",
    "batch_size = 100\n",
    "data_path = 'shape-transitions'\n",
    "pretrain_data_path = 'moving-squares'\n",
    "shape_library = load_shapes('PDE-Control/notebooks/shapes')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "uDvxnTavZQk2"
   },
   "source": [
    "The `shape_library` in the last line contains ten different shapes that we'll use to intialize a marker density with at random positions.\n",
    "\n",
    "This is what the shapes look like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 142
    },
    "id": "9m8HFEPzQD70",
    "outputId": "010df794-01e9-49b5-b4a2-492bb4050311"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9YAAABwCAYAAAD7ROjcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXhU933v8fdvRqNZtCIhBAhsMIsRYOwkxHbi7EndOIvj2zrOUjtJnca9t0ma3nSx2+Y+bnqbxk33J86TlMZuvSQhvsRN7MSNa1zvO16wDWIzYINAoBXt0syc3/1DIyEIoGWW3zkzn9fzzAM6M5r5zvfLV5yv5pzfMdZaRERERERERGR2Qq4DEBEREREREQkyDdYiIiIiIiIiWdBgLSIiIiIiIpIFDdYiIiIiIiIiWdBgLSIiIiIiIpIFDdYiIiIiIiIiWSib6gHGmBjwKBDNPH6TtfZGY8xSYCNQDzwPXGOtHT3Tc5WbqI1RkX3UJcZiGaQPDw+w25llDZT/2RnPP4BHOg18Qz1QOCfm3wPsXyr/haUecEs94J56wC31gHu52hcF1WA2irEHTDjEcFOMNbXthDCnfZyHZVtPA7GDw1jPK2CEpzbMAKN25JQBm6muY22MMUCFtbbfGBMBHge+AnwVuNtau9EY8z1gq7X2u2d6rmpTZy8y75/Vmyhl1lrSpHmeR+iju5xZ1kD5n53x/JeZMjbbTS8AKdQDBTM5/0/bzfTT8yzKf0GpB9xSD7inHnBLPeBervZFQTWYjWLsgXB1NS3fWsWOj36HqImc9nGD3ihr7vkSq65vwevrK2CEp/aMfZBe23XKwXrKQ8HtmP7Ml5HMzQLvAzZltt8GXJGDWOUUjDGUmYmDC1SDAjsp/wblv6BOyj8o/wWnHnBLPeCeesAt9YB72hd1Sz0QDNM6x9oYEzbGvAQcBR4AXgN6rLWpzEMOAk35CVFg7DdVA/SCauCEtZan7QMA56P8F9x4/vvpAeXfCfWAW+oB99QDbqkH3NO+qFvqAf+b1mBtrU1bay8AFgEXAqum+wLGmOuMMVuMMVuSjMwyTDHGUEE1zLAGyn9uGGO42PwawMuoBwpuPP+V1IDy74R6wC31gHvqAbfUA+7Ndl80872qQZbUA/43o1XBrbU9wEPA24BaYyaOSVgEtJ7mezZYa9dba9dHiGYVrMy8Bsp/zqVRDzhjxn5kKf9uqQccUg/4gnrAIfWAe5oH3FIP+NeUg7UxpsEYU5v5exz4NaCFsYJemXnYZ4Gf5SvIUjdqR0hmFvhTDQpvcv4ZO7dO+S+gyfm3WFD+C0494JZ6wD31gFvqAfe0L+qWeiAYprzcFrAAuM0YE2ZsEL/LWvtzY8x2YKMx5q+AF4Fb8hhnSRthiG1sGV9m/zlUg4Iazz9jK+ivBr6p/BfO5PxneuAB5b+w1ANuqQfcUw+4pR5wT/uibqkHgmHKy23lkl+Wdw+qMy3vPh3Kf/Y2203PW2vXz/b7VYPsqAfcUw+4pR5wTz3glnrArWzzD6pBtoqlB0ryclsiIiIiIiIicnoarEVERERERESyoMFaREREREREJAsarEVERERERESyoMFaREREREREJAsarEVERERERESyoMFaREREREREJAsarEVERERERESyoMFaREREREREJAsarEVERERERESyUOY6ANdMNEp4bj2UR45vHE2S7ujEjoy4C0xERKQYGEO4vg5TVZmzp7S9faS7usHanD2nSNCFqqoI1dWCZ/E6u/AGB12HJFJSSn6wZu1yWq6tpGbRsYlNx96Yz7m31GJf2u4wMBERkeALV1XRes258O5uQib7QThtDeahJhZ+fyvewEAOIhQpDv0fWM3hj4/gjYZZekcjZQ8+7zokkZJS8oP18PwE177zEb42d8fEthuXrOHRX7ydcodxiYiIFIVolP43D7HrrT8gbLI/Ay1p05zbcx0mFgUN1iITepaF2fS2f6E1XcNfPPrb1LkOSKTElORgbcrK8C5aS9eqOMfOhebYoRPub44d4ofvKaN24duYs2OQ0LPbsKmUo2hFRCSfypoWcuxtixmuzc2yI5FBy5xnj5Desy8nzycyHeGVy+h6awOpuHEdygmq948SfbJFhyXnWOj8ZjovqMWbdCbj4AVD1IWTtKbdxSVSykpzsI7H2fexODdd8QOayrpZW54EYhP3f6TiMEuuupk3UnV87e5Ps2xrVIO1iEiRGl61gLlf3s8Xmh7JyfM93NvM4397EdUarKWAOi+ax3u++hTvqtox9YML6MuP/xbNO+s0WOeSMRx63xyu+8K9LClvByCMZXFZDwvCcV5xHJ5IqSqpwdpEo4Sqq6GuhnTjCB+p6CRqIkD4hMdVhmJcHIN1XjvXz0timuYT7u7F6zmGTY66CT6ojCFcU42JxyE0w0+DPA87MEi6tzc/sYmIAOlYiPVzXufDieEcPWMLD8cvztFziUxPKgHvr9rGpYmk61BI2jSH00N0pSOEytMQ8ten6EETqqggVF11fD8qFGKo0fLRyhbOKpu8KGDcSXwiMmbKwdoYsxi4HWgELLDBWvvPxpg64MfAEmA/cJW1tjt/oWbPnr+SHVcniDf187/OfYSykwbqk0VNGf/zrY9w219fxMjBBlbc2Q9bXi1QtGOG7SDbeI5RhhliAGPMV4KU/1BlJYevXsPAO/oxM/x/1UuHqHkwztw7X3C6QvvkGgBrglaDoAt6DxQD9YBb6gH31AMzsys5yuWPfZnE1jgLXkvjdfVk9Xwl3QOhMD0fO4+jHxqhLHL8GO8PLttCXahwn4+N12CAYxhjthHgeSCISroHAmQ6HZkC/tBa+4Ixpgp43hjzAPA54EFr7U3GmBuAG4Dr8xdq9gbOSvCHH7iPL9YeyGw58yeoYRPi+vrdXP/23fxT9xJ+8sivk9iS/zgnMxhWsI5qM4en7Wb66flikPJvYjF6Lxpixzv/jYg58y8yTtbvDXNez+/TcFe508F6cg02200tQKBqEHRB74FioB5wSz3gnnpgZg6kapnzcIz6W54EwMvy+Uq5B0w4TOd5hsff9W0WlJ18ybrYKb8nL3FkatDCC/TRfTEBngeCqJR7IEimHKyttYeBw5m/9xljWoAm4GPAezIPuw14GB8W0kTKsW9ZRe/SBB3nm4lzUWbqnPKjHLkwRH3sYqr3DWK2bC/IeddREyeaObTHYAB8m38TjeKtb6bv7OM/6JMVhhVNbxBi5oeBRUyYxiVddFy5lrKRzCVaLAXNP5xYA8b2D3xbg2IU9B7Ip8iAR9XzraQOtub1ddQDbgWpB4pVEHqg8lCaP3r1SlY3HOHaxsecHhKexpCDK6tNKJkeCIUx56+i99wqbOazHy9sKFvRRyI0sw8ncm2iBjaY80DQlUwPBNyMjiExxiwB3gQ8AzRmhm6ANsYOFfedUHUlOz+R4P9cdjeLI51cFB1gNr/he3e8k5s//n0OXFHPN39xBSt3VJDuOTb1N+aQRxp8nP/wnFp2Xl3O19+7aWJbuUnz5thBwqZixs8XNRFuWX0HL92wiLQdG8yTtsxZ/jPK8XENil0QeyCf7utcx6FvLSee58H6JOoBh/zeAyXClz1Q8fhuKvbO42jTOXz9q7Vcuu5uV6HkVTH3QKg8wv7La/jSJ+6lKjQEQNhY1pQfotJEHUd3XBDngWJSzD0QdNMerI0xlcBPgD+w1vaaSSfMWmutMaf+3aQx5jrgOoAYieyinQETKSdUWQFz6wgvGOSaqrbM9TNn90lSTSjOpYkkSdvKNxcOQ0M9YRPC6x8oyIJmKZtiiAHwc/7Lyog3DPKZ6o6T7pj5UD1uTXmcNeWdE18nbZpvzE1CuPC/uU3ZFMAy4Bpf1MAYQpWVmKh/r7huR0bx+vvBZv/RRbB7IF9e5rsVKwu2XI3veqDEBKIHipyfeyDd3Q3d3cSPNdHWPycvrzGVfm+Y9nSK/aOLCeXhoLJi64FQIoFJHP8JbuJxhuenuLp6FzWhyT/Z/TNUWywEaB4oNsXWA8VmWoO1MSbCWBP9wFo7/ivQI8aYBdbaw8aYBcDRU32vtXYDsAGg2tTl8MCgM/PWN7P743FC84e5bu1jmaE6exET5gvnPcG/fv0SODSPZZsG4emXc/Lcp+NZj5d5igjlDNvBQOS/2IzXAOjySw+E6+toveZc+s53d/75VKpejNJ05w7SnV1ZPY96wD0/9kApUQ+4px6Y2p+1vZNfbl5P/Ihh4UtdWZ9bPVmx9UAoFqPzqvNpf1cSExoLyYQ9rlr7LAnjz1+Ye9ZjiH4I0DxQTIqtB4rRdFYFN8AtQIu19h8m3XUP8FngpsyfP8tLhLPUuyzOjR/axKeqjsx40aypXF+/m6++ewe39zbx/eevoPrpnD79Cay1bGcLFVSR5IRPxn2d/2IyuQZdHDky6S6nNTDVVYTf28nOt/ywkC87I2+pvRpzTxVkMVirB9zzaw+UCvWAe+qB6bn/tWZWfOcAqYOteDk4UmlcMfaAKS+n/eI0r1x6c+bSr2NCGMI53m/NhfEahAiTssnAzAPFohh7oBhN5xPrS4BrgFeMMS9ltv0ZYwW8yxjzeeB14Kr8hDh9JhqFdSsZWJSgc63hrEhXzofqcRETZnGkk851hrLhi6g4MABbd+X8sPBjdNLGG1RSwyD9ZGrgq/yXLV7EwHkL6GosY3nDvry+VgjDogVddF22knh7isTLB0kdbsvra06uAbDaLzWwIUM4ZPP2bzwXwjlYvSYIPVDs/NoDpUI94F6QesCOjJLaVcXVC9/D+dUH+O2al5kbnv0pWTPmeTk5/WeyoPeAiZTDuhUMnHV8Ve9U3DB3URdRE/H1/+PjxmsQIozf54FiFPQeKBXTWRX8cTjtks7vz2042QnPqaXlMwl+/333c075UdZHB8nnpQguivXyd1fext7L5/HtBz7Iua/XkG6f3arjp1Nr5vIBrgTgGfsgvbbrgkl3+yL/3Zcs4pzf38G75uzinfE9kMdzN8ImxM3n/ogn/2wZPz+6jr6/O4vYz/M7WE+uwWa7abu1dv2ku31Rg2IWhB4oduoBt9QD7gWpB7zubpbf2sbRe8/m39+7mvOuPcgHE/49ZWg6gt4DoepKdn2ymt+77H5CZuwA+YhJc0l8DxEfLUp2JuM1OEX+IQA1CLqg90CpKNyV5fPIlJVholFsTRUVi/r4gzn7M/fk95I3NaE4l1cMQsV+NizqhznVhAYHsSMjBbsUlB+M1Bg+M++JzKU98r8gwrryGOvKW0mERvhu9dkFvIqjiIiIv9lUivSefZg9ULn8bfSkE0D+Buu09ei3Iwx6aVLJ8Ngn1iUuFItB5Pjh3aamGtM0xFfm7DlpzZ9gDNW5ZkIhQpVVrsMILNOfm3WjJPeKYrD2LlrLvo/GSTWO8nsrHnYSw2fPfYbv3fBuIkfnseRng5intjqJQ0RERKRQXhlN8unnv0BqdxXzXrF4vX2uQ3IqXF3N0U+sofNN6eMbox6fW/NkzhbSDbrhRTF2/Gmz6zACa/ib/+U6BDmNohisu1bH+dvfvINL411ETRlQ+B9cf1S3ky9fuo2fDzbwtzs/Td1TBQ9BREREpKC2jS6k6mdVzNn4PDadxvPSU39TETPVVfR9YIBXLtlAaNL+qKv9Uz9aU9vOo5ff7DqMwHrXd3N72qnkTmAHaxONwurljDQm6FkJC8u6SYTcXZ4gbEIkTDlnlXVx7Fyo/OBbibUNYFtew44E+9wmERERCZ54V4pv73svLzbu5fKaF7gklvvBLm0NobTN+eKtQTB5X3Rcb12YpQ0HqAzpRLXTCWGc7rMHXei0S1+Ja4EdrMMNc2n5nUquecfjrIkfZHUkDUSm/L58W1ue5C//x0ZaLlvIHY++g+ab6km1HnIdloiIiJSYxLP7SX19IU81Xsgj1y3nifPv0uHIOTR5X3RcIjTKr1e+SqmePy1SygI7WFMeobKpl6/WbyFCmIRPfjNYGYpxecURLqs4xE+aLoBIcFMsIiIiwZVub8e0t1PZ0MDO31yU0+dO2jSDdpQ+L44p1fXKfmVfdPxTWA3VIqUosFOfPdZL5L4FvHn//2b+qqPc2Xw7SyOVU39jnr2W7Oe3tn2O9l1zmfOqwR7Tp9UiIiJSPNLW42tH38JdT19IrK2Ms3ceI7dXrg4Gv+6LiogbgR2s051dNPzb88yLlNF27QXsXVHD0oj7BTN2J+sZvaeBlbdvxSZTpEvwnCMREREpXh6WTa++idXfaMXr6MQbTboOyQm/7ouKiBuBHawBbHIUmxwl3u7xndb3sbXuNT5QuZ115YU/LPzl0WE296/m0c4VJNo9vMHBgscgIlIKyuY3klw6H1uem3NF+5rKeOjoypw8F8CW7rOJHivVY2PFl1IpQq0x/qpjLctjR/hQ4gBzwompv+8MbDqEHRrCGx7OUZDBdPK+6FO1BybuWxk7zK8n2qgJxR1GKCKFEujBetycJw7Q03kWd81fxubPr+K+c+8reAx//NqVHLt1ERVtSWp2HCBV8AhEREpD1/uWUv35gyyv7sjJ8z3RupSeu5rY3NqYk+cLD3v6f0B8xevrY/kPu3nosUvYuD5CxdW3c0VFv+uwisr4vujm6JKJbXe8vYz6T27g/XF9ii1SCopisE4dbCVysJX6RU3s+ehcBleMEjFhIiac99dO2jRJm2bvkbms/O/9pA63aWdKRCSPBueF+Odz/oMLo7m5EsTvAbt3NhN65MWcPB+g/wfEV2wqhX15B9GXoT5+EW3JGmB2g3XSphmxSaynS/5MNr4vOtmc2os5kKxnJHZ8e4hQQfZP/czDMmJL8/SBXPBKckWDYCiKwXqc7e+n8r/PYk33l2ha0sGtzXewMlKRt9drGR3kd3ZczaH9c6l/Powd0OHfIiIiUpxeGhnhC9uuoeNgLQ1Ph7FDpX0Y+FRqdg/w1z/9Tb5ee/wT64bF3Xx/zR1OTlv0i209Day694uuwwistp5/dh2CnEZRDdbpnmPMu+1FGn9UztGr1vDKnyxgZaQ3b6+3daSJ0bsaad60DTs6SrrEzzMSERGR4vXE0HLKflhH88+3Y0dHS/786qmYF1pY1hLHmOOf7nddvpqn//wc1pWX7lVjYq3DNP/JDtdhBFZPv/rOr4pqsAbGfsgPD5PoSPOjIxfSk36Vt8f30lye3SIdk20bHeLpoaXce/R8Eu1p0r35G96DoLzX8qOOizlUu5O3x/fl9SgBODH/5X1aIEhERGS6Iv0eP227gJCxXBzfO6NPTpM2TNmQLfn9numyqRS2r++EbfGONHcffhOePT5sz48c492xo1kvKBcUNu3p31AWrNW+r18V3WA9rmpLKx3fXMq/zFvBD69u5b+af0rYZL+CbNp6fHHXpxi6cwHxjhSVWw+W/Ll0dU+2srevmVca1/LTz+zlpyvuz9trKf8iIiKzl3jpDQa/dRa3zl3M9z/VyRMXbCz5c34LaTz/d1YumtjWuSbM1z79Y36rqtNhZCKSraIdrFMHW4kebCXR0MCO9y0m1ZwGS1bDddp6pEjzems9zb/YSbqzS0MdkHr9ALHXD1C5qIldlzXAivy9lofljUP1rFL+MZ7Fs2P/Lv1q8m/kRUTEvVTbEcp/eYR4dTUt71iFd4EHTD1Yp62HZ3NzibtSNp7/8knbyobeyp7faCRd2X7CY3PxgZCIFM6Ug7Ux5lbgI8BRa+3azLY64MfAEmA/cJW1tjt/Yc6eHR6m9skozfwuSxd2cPPyjbM6LHzb6BBf2v1JXj9cT91T5QVdsGOb3UIHh0lzfAXFINUg6MbzX050Ypsf8m/7Bkg+vJCVPZ8v5MvOSOULcWzf0ayfRz3gll97oJSoB9wq5R5IW4+/6ljLndsuxLbGWf5Gr5M1iYu5B+IH+rjzP9/N7U0XTWxrqO/j71fdxSUxfwzXxZz/oFAN/G86n1j/O3AzcPukbTcAD1prbzLG3JD5+vrch5c9r6+Pxh+8yvyfxOj48HKevP4cmsvbZvw8jw0uZ/D2haz65R7s8DDeYOFWAF/I2SxmGVt4aPLmwNQg6Mbzv43nJm92nv90RweL/nUEotGpH+zKyAjpk84vmw31gFt+7YFSoh5wq5R7YMSm+Pctb2f1Nzqwx97A63Vz/eti7gG7fQ8r/qYKyo7vlvdfspR7bnwzl8RechjZccWc/6BQDfxvysHaWvuoMWbJSZs/Brwn8/fbgIfxcRG9vj7o6yNxdAn/2bGW2vAg50UPT2uRrZbRQbaPzueX7WtJtKdIt7dP+T25Nsc0MGQHTt7svxqkUgx1JNjYN2diU8SkuSB6iGWRylk95Xj+k3bsMLVRG8Z0RSCdnuI7c8e3+bels4CMb2tQIpR/91QDt4ox/9ZaIj1hftTXxOJIJ+uj/dSE4qd+8EgIr71zbH/KkWKswTibSpHuPvFDxtiRhTx+5Bw2Jl6f2FYbHmR9tIu54fwuEnsqxZz/oFAN/G+251g3WmsPZ/7eBjTmKJ68qnj5EEf//hy+OXcFsSuP8NB5/++MC3aM2CTXtlxDctM8Ep1pql7x1UJZvquB13OM5T9I8k+PfwKbubREsgLmf+J17l358xmfKzQ5/+HRsW3GWpa9NoQ3MJTr8GfKd/kvQaqBW8q/e6qBW4HOvx0aYund/Xzv1d+gu9nwfz/+Q66qPOY6rJkKdA3OpGzPISLfXsI/1n1yYlv/YsPnPn0/f1z3msPITlC0+Q8Q1cBHsl68zFprjTGnPd3GGHMdcB1ADLeXEUgdbCV+sJXK6mpa1q/CO+/MC3YkbZq2ffWsuutVvL4+Pw3VJzhTDQqZf294mNBjL1Lz2PFt4YYGdl6yCG+lncbSKCc6Of+TuTi/63SC1APFyi89MBZL6S3Yph5wz089UIqC2AM2lYJnX6HmWYh9+K3s+ugCOMVg7eHfBTInK7YeSLe3E/1FO5NP9qp52/ls/chi8M9gPSGIPVBsiq0Hgmi2g/URY8wCa+1hY8wC4LSrE1lrNwAbAKpNnS/mITs6St2LYdbN+TzL5nVw09K7T7iO40sjI9yw7zfY115P/YthSCbP8GzOTKsGrvNvh4epei7OutjnOMPP21NKp0OBzz+4r0ER810P2IFBIk+fTTPXzPjf+2wMtSdY8XrhFlI8iXrAPd/1QIkp6h7QvpD/lLX38uzDzaxeefxSXXMqB7lx+b1cmnBSn6LugYAoqR7wu9kO1vcAnwVuyvz5s5xFVADe8DCNP97O/F9U0vmes9h8w2rW1e2duP/+/rUc+/5ilj96ANu7n/Swsx3XMwlEDbz+fhbe0YK5e3a/HVP+5Qx8V4N0Tw+Lbm3B/LhAvw1OtZPu7nF1BIfv8l+CVAO3ijr/2hfyH2//AZb/Yz8mevxiXUOrF3Drje/k0qX/7SKkksq/T6kGPjKdy239iLGT4ucaYw4CNzJWvLuMMZ8HXgeuymeQ+ZDuOQY9x0i0zefxrmWsjR2YuO+JrmVUHB4ldbDVYYTHvWKfoZt2PDwCVwNrxxbk6A7uyv/j+U8yArAuk/Ng5L9IBKYHiuDf+6moB9wLTA8UqWLvgdCIxws9i3mwYjfnlB1jaaSSQa+cWFda+0I+YlO/uohurKaS7e2N/FdjZGLbo32rKBvO7a9clX/3VAP/m86q4J86zV3vz3EsTsR3tNH2nWX8ae3yiW2xbkvtrjd8c071eWbsuobP2AfptV2LJt1VFDXwu/H8A2y2m1621t6S+VL5LxD1gFvqAffUA24Vew+M7wv9UcMKEh9p44G1G12H9CvUA6dx6ChzblnBny78nYlNZYNQ/9wRcnn9FOXfPdXA/7JevCzoUgdbqfpxK1Unb3cSjYiIiEhhje8L1SQS7Fq+juTawl3SUrKT7u4mdu+zxE7e7iQakdJW8oO1iIiI5NHICBUvxbmw9pOEcrBoftozVG2NYodHsn8yOYFNpqjdbnjnlmvpb61m5ZF+X12FQ0TEzzRYi4iISN6k+/pYdPtuzH9UYHMwWRvPYnt3kh4czEF0MplNjtL4k52YzdUw2ku6vcN1SCIigaHBWkREAiU8bHlu6Bxg75SPnY7XeudiRoNxrd5AsnZswaWTFl0Sf0p3dkFnl+swREQCR4O1iIgESuNTPdzCR9gQy8FxxUDiqEfdvn1aW0NERERmTYO1iIgEire1hYatuX1ODdUiIiKSjZDrAERERERERESCTIO1iIiIiIiISBY0WIuIiIiIiIhkQYO1iIiIiIiISBY0WIuIiIiIiIhkQYO1iIiIiIiISBY0WIuIiIiIiIhkQYO1iIiIiIiISBY0WIuIiIiIiIhkQYO1iIiIiIiISBY0WIuIiIiIiIhkIavB2hjzQWPMTmPMHmPMDbkKSqZPNXBL+XdPNXBL+XdPNXBL+XdL+XdPNXBL+fePWQ/Wxpgw8B3gMmA18CljzOpcBSZTUw3cUv7dUw3cUv7dUw3cUv7dUv7dUw3cUv79JZtPrC8E9lhr91prR4GNwMdyE5ZMk2rglvLvnmrglvLvnmrglvLvlvLvnmrglvLvI9kM1k3AgUlfH8xsk8JRDdxS/t1TDdxS/t1TDdxS/t1S/t1TDdxS/n2kLN8vYIy5DrgOIEYi3y8nJ1H+3VMN3FL+3VMN3FL+3VMN3FL+3VMN3FL+CyObwboVWDzp60WZbSew1m4ANgAYY9o3200DQEcWr+sXcyn8+zj7pK+nrMEp8v86bmLPB9c1UA8EsweKJf/gvgbqAfWAa65roB5Q/l0rdA1m/DMIiroG6oHZOgb8LiR+dzoP/hPeOP2drnvgOGvtrG6MDeV7gaVAObAVWDON79sy29f0080P70M1cPs+lH/372M2NfBD3MVSA/WA+/ehHlAPKP+lm38/vJdSr4Hr91Hq+ffbe5n1J9bW2pQx5kvA/UAYuNVau222zyczpxq4pfy7pxq4pfy7pxq4pfy7pfy7pxq4pfz7S1bnWFtr7wPuy1EsMguqgVvKv3uqgVvKv3uqgVvKv1vKv3uqgVvKv39ksyr4bG1w8Jr5EOT3EeTYJwvq+whq3CcL6vsIatynEtT3EtS4TxbU9xHUuE8lqO8lqHGfLKjvI6hxn0pQ30tQ4z5ZUN9HUOM+Fd+8F5M5NuXcSnMAAAK7SURBVF1EREREREREZsHFJ9YiIiIiIiIiRaOgg7Ux5oPGmJ3GmD3GmBsK+drZMMYsNsY8ZIzZbozZZoz5SmZ7nTHmAWPM7syfc1zHeibKv3uqgVvKv1tBzT+oBq4p/+6pBm4p/+6pBm4FIv8FXAo9DLwGnMPx5eBXu14WfZqxLwDenPl7FbALWA18C7ghs/0G4G9cx6r8u49XNXAfr/Lvv1uQ868auL8p/+5vqoHz2JV/9/GrBsr/GW+F/MT6QmCPtXavtXYU2Ah8rICvP2vW2sPW2hcyf+8DWoAmxuK/LfOw24Ar3EQ4Lcq/e6qBW8q/W4HNP6gGrin/7qkGbin/7qkGbgUh/4UcrJuAA5O+PpjZFijGmCXAm4BngEZr7eHMXW1Ao6OwpkP5d081cEv5d6so8g+qgWvKv3uqgVvKv3uqgVt+zb8WL5sBY0wl8BPgD6y1vZPvs2PHH2iJ9TxS/t1TDdxS/t1TDdxS/t1TDdxS/t1TDdzyc/4LOVi3Aosnfb0osy0QjDERxor4A2vt3ZnNR4wxCzL3LwCOuopvGpR/91QDt5R/twKdf1ANXFP+3VMN3FL+3VMN3PJ7/gs5WD8HrDDGLDXGlAOfBO4p4OvPmjHGALcALdbaf5h01z3AZzN//yzws0LHNgPKv3uqgVvKv1uBzT+oBq4p/+6pBm4p/+6pBm4FIv+5WgVtOjfgQ4yt4PYa8OeFfO0s434HY4cVvAy8lLl9CKgHHgR2A5uBOtexKv/u41UN3Mer/PvzFtT8qwbub8q/+5tq4Dxu5d997KqB8n/Gm8kEKiIiIiIiIiKzoMXLRERERERERLKgwVpEREREREQkCxqsRURERERERLKgwVpEREREREQkCxqsRURERERERLKgwVpEREREREQkCxqsRURERERERLKgwVpEREREREQkC/8f+QhFu8t4LFIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1224x360 with 10 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light",
      "tags": []
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pylab\n",
    "pylab.subplots(1, len(shape_library), figsize=(17, 5))\n",
    "for t in range(len(shape_library)):\n",
    "    pylab.subplot(1, len(shape_library), t+1)\n",
    "    pylab.imshow(shape_library[t], origin='lower')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "eQdzVAJH30dg"
   },
   "source": [
    "The following cell uses these shapes to create the dataset we want to train our network with.\n",
    "Each example consists of a start and target (end) frame which are generated by placing a random shape from the `shape_library` somewhere within the domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "MaAzExNW30dh",
    "outputId": "14de4d69-177c-4f0e-c603-5f041ce72cce"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shape-transitions/sim_000000\n",
      "shape-transitions/sim_000100\n",
      "shape-transitions/sim_000200\n",
      "shape-transitions/sim_000300\n",
      "shape-transitions/sim_000400\n",
      "shape-transitions/sim_000500\n",
      "shape-transitions/sim_000600\n",
      "shape-transitions/sim_000700\n",
      "shape-transitions/sim_000800\n",
      "shape-transitions/sim_000900\n"
     ]
    }
   ],
   "source": [
    "for scene in Scene.list(data_path): scene.remove()\n",
    "\n",
    "for _ in range(example_count // batch_size):\n",
    "    scene = Scene.create(data_path, count=batch_size, copy_calling_script=False)\n",
    "    print(scene)\n",
    "    start = distribute_random_shape(domain.resolution, batch_size, shape_library)\n",
    "    end__ = distribute_random_shape(domain.resolution, batch_size, shape_library)\n",
    "    [scene.write_sim_frame([start], ['density'], frame=f) for f in range(step_count)]\n",
    "    scene.write_sim_frame([end__], ['density'], frame=step_count)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AIh2z-be30di"
   },
   "source": [
    "Since this dataset does not contain any intermediate frames, it does not allow for supervised pretraining. This is because to pre-train a CFE network, two consecutive frames are required while to pretrain an $\\mathrm{OP}_n$ network, three frames with a distance of $n/2$ are needed.\n",
    "\n",
    "Instead, we create a second dataset which contains these intermediate frames. This does not need to be very close to the actual dataset since it's only used for network initialization via pretraining. Here, we linearly move a rectangle around the domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "A0QlQ-AL30di",
    "outputId": "5f464216-98ba-4c10-c660-887d3394ba43"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "moving-squares/sim_000000\n",
      "moving-squares/sim_000100\n",
      "moving-squares/sim_000200\n",
      "moving-squares/sim_000300\n",
      "moving-squares/sim_000400\n",
      "moving-squares/sim_000500\n",
      "moving-squares/sim_000600\n",
      "moving-squares/sim_000700\n",
      "moving-squares/sim_000800\n",
      "moving-squares/sim_000900\n"
     ]
    }
   ],
   "source": [
    "for scene in Scene.list(pretrain_data_path): scene.remove()\n",
    "\n",
    "for scene_index in range(example_count // batch_size):\n",
    "    scene = Scene.create(pretrain_data_path, count=batch_size, copy_calling_script=False)\n",
    "    print(scene)\n",
    "    pos0 = np.random.randint(10, 56, (batch_size, 2))  # start position\n",
    "    pose = np.random.randint(10, 56, (batch_size, 2))  # end position\n",
    "    size = np.random.randint(6,  10,  (batch_size, 2))\n",
    "    for frame in range(step_count+1):\n",
    "        time = frame / float(step_count + 1)\n",
    "        pos = np.round(pos0 * (1 - time) + pose * time).astype(np.int)\n",
    "        density = AABox(lower=pos-size//2, upper=pos-size//2+size).value_at(domain.center_points())\n",
    "        scene.write_sim_frame([density], ['density'], frame=frame)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "cwrv1oHs30dj"
   },
   "source": [
    "## Supervised initialization\n",
    "\n",
    "First we define a split of the 1000 data samples into 100 test, 100 validation, and 800 training samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ZcbTw9cF30dk"
   },
   "outputs": [],
   "source": [
    "test_range = range(100)\n",
    "val_range = range(100, 200)\n",
    "train_range = range(200, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "22TTd6xJ30dk"
   },
   "source": [
    "The following cell trains all $\\mathrm{OP}_n \\,\\, \\forall n\\in\\{2,4,8,16\\}$. Here the $n$ indicates the number of time steps for which the network predicts the target. In order to cover longer time horizons, we're using factors of two here to hierarchically divide the time intervals during which the physical system should be controlled.\n",
    "\n",
    "The `ControlTraining` class is used to set up the corresponding optimization problem.\n",
    "The loss for the supervised initialization is defined as the observation loss in terms of velocity at the center frame:\n",
    "\n",
    "$$\n",
    "L_o^\\textrm{sup} = \\left| \\mathrm{OP}\\Big(d_{t_i},d_{t_j}\\Big) - d^*_{(t_i+t_j)/2} \\right|^2 .\n",
    "$$\n",
    "\n",
    "Consequently, no sequence needs to be simulated (`sequence_class=None`) and an observation loss is required at frame $\\frac n 2$ (`obs_loss_frames=[n // 2]`).\n",
    "The pretrained network checkpoints are stored in `supervised_checkpoints`.\n",
    "\n",
    "*Note: The next cell will run for some time. The PDE-Control git repo comes with a set of pre-trained networks. So if you want to focus on the evaluation, you can skip the training and load the pretrained networks instead by commenting out the training cells, and uncommenting the cells for loading below.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "b5ksGZwO30dl"
   },
   "outputs": [],
   "source": [
    "supervised_checkpoints = {}\n",
    "\n",
    "for n in [2, 4, 8, 16]:\n",
    "    app = ControlTraining(n, IncompressibleFluidPDE(domain, dt),\n",
    "                          datapath=pretrain_data_path, val_range=val_range, train_range=train_range, trace_to_channel=lambda _: 'density',\n",
    "                          obs_loss_frames=[n//2], trainable_networks=['OP%d' % n],\n",
    "                          sequence_class=None).prepare()\n",
    "    for i in range(1000):\n",
    "        app.progress()  # Run Optimization for one batch\n",
    "    supervised_checkpoints['OP%d' % n] = app.save_model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "_A9-V1Vg30dl",
    "outputId": "994c7fdc-a5aa-4769-eab8-1f27f29ad082"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'OP16': '/root/phi/model/control-training/sim_000003/checkpoint_00001000',\n",
       " 'OP2': '/root/phi/model/control-training/sim_000000/checkpoint_00001000',\n",
       " 'OP4': '/root/phi/model/control-training/sim_000001/checkpoint_00001000',\n",
       " 'OP8': '/root/phi/model/control-training/sim_000002/checkpoint_00001000'}"
      ]
     },
     "execution_count": 8,
     "metadata": {
      "tags": []
     },
     "output_type": "execute_result"
    }
   ],
   "source": [
    "supervised_checkpoints # this is where the checkpoints end up when re-training:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "jD7nKXCv30dl"
   },
   "outputs": [],
   "source": [
    "# supervised_checkpoints = {'OP%d' % n: 'PDE-Control/networks/shapes/supervised/OP%d_1000' % n for n in [2, 4, 8, 16]}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "imuKMrRG30dm"
   },
   "source": [
    "\n",
    "This concludes the pretraining of the OP networks. These networks make it possible to at least perform a rough planning of the motions, which will be refined via end-to-end training below. However, beforehand we'll initialize the $\\mathrm{CFE}$ networks such that we can perform _actions_, i.e., apply forces to the simulation. This is completely decoupled from the $\\mathrm{OP}$ networks.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "bklXI_o130dm"
   },
   "source": [
    "## CFE pretraining with differentiable physics\n",
    "\n",
    "To pretrain the $\\mathrm{CFE}$ networks, we set up a simulation with a single step of the differentiable solver.\n",
    "\n",
    "The following cell trains the $\\mathrm{CFE}$ network from scratch. If you have a pretrained network at hand, you can skip the training and load the checkpoint by running the cell after."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "kCtOr3ed30dn"
   },
   "outputs": [],
   "source": [
    "app = ControlTraining(1, IncompressibleFluidPDE(domain, dt),\n",
    "                      datapath=pretrain_data_path, val_range=val_range, train_range=train_range, trace_to_channel=lambda _: 'density',\n",
    "                      obs_loss_frames=[1], trainable_networks=['CFE']).prepare()\n",
    "for i in range(1000):\n",
    "    app.progress()  # Run Optimization for one batch\n",
    "supervised_checkpoints['CFE'] = app.save_model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "-KOcgr5M30dn"
   },
   "outputs": [],
   "source": [
    "# supervised_checkpoints['CFE'] = 'PDE-Control/networks/shapes/CFE/CFE_2000'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "34EA3S98L-Wr"
   },
   "source": [
    "Note that we have not actually set up a simulation for the training, as the $\\mathrm{CFE}$ network only infers forces between pairs of states. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "kQZLbOolI4-2"
   },
   "outputs": [],
   "source": [
    "# [TODO, show preview of CFE only?]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "g9HcP2oK30do"
   },
   "source": [
    "## End-to-end training with differentiable physics\n",
    "\n",
    "Now that first versions of both network types exist, we can initiate the most important step of the setup at hand: the coupled end-to-end training of both networks via the differentiable fluid solver. While the pretraining stages relied on supervised training, the next step will yield a significantly improved quality for the control. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "NipTaC2S30do"
   },
   "source": [
    "To initiate an end-to-end training of the $\\mathrm{CFE}$ and all $\\mathrm{OP}_n$ networks with the differentiable physics loss in phiflow, we create a new `ControlTraining` instance with the staggered execution scheme.\n",
    "\n",
    "The following cell builds the computational graph with `step_count` solver steps without initializing the network weights.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "iJjrTrl530dp",
    "outputId": "69c6407e-efb2-41eb-df62-e472e9442834"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "App created. Scene directory is /root/phi/model/control-training/sim_000005 (INFO), 2021-04-09 00:41:17,299n\n",
      "\n",
      "Sequence class: <class 'control.sequences.StaggeredSequence'> (INFO), 2021-04-09 00:41:17,305n\n",
      "\n",
      "Partition length 16 sequence (from 0 to 16) at frame 8\n",
      "Partition length 8 sequence (from 0 to 8) at frame 4\n",
      "Partition length 4 sequence (from 0 to 4) at frame 2\n",
      "Partition length 2 sequence (from 0 to 2) at frame 1\n",
      "Execute -> 1\n",
      "Execute -> 2\n",
      "Partition length 2 sequence (from 2 to 4) at frame 3\n",
      "Execute -> 3\n",
      "Execute -> 4\n",
      "Partition length 4 sequence (from 4 to 8) at frame 6\n",
      "Partition length 2 sequence (from 4 to 6) at frame 5\n",
      "Execute -> 5\n",
      "Execute -> 6\n",
      "Partition length 2 sequence (from 6 to 8) at frame 7\n",
      "Execute -> 7\n",
      "Execute -> 8\n",
      "Partition length 8 sequence (from 8 to 16) at frame 12\n",
      "Partition length 4 sequence (from 8 to 12) at frame 10\n",
      "Partition length 2 sequence (from 8 to 10) at frame 9\n",
      "Execute -> 9\n",
      "Execute -> 10\n",
      "Partition length 2 sequence (from 10 to 12) at frame 11\n",
      "Execute -> 11\n",
      "Execute -> 12\n",
      "Partition length 4 sequence (from 12 to 16) at frame 14\n",
      "Partition length 2 sequence (from 12 to 14) at frame 13\n",
      "Execute -> 13\n",
      "Execute -> 14\n",
      "Partition length 2 sequence (from 14 to 16) at frame 15\n",
      "Execute -> 15\n",
      "Execute -> 16\n",
      "Target loss: Tensor(\"truediv_16:0\", shape=(), dtype=float32) (INFO), 2021-04-09 00:41:44,654n\n",
      "\n",
      "Force loss: Tensor(\"truediv_107:0\", shape=(), dtype=float32) (INFO), 2021-04-09 00:41:51,312n\n",
      "\n",
      "Supervised loss at frame 16: Tensor(\"truediv_108:0\", shape=(), dtype=float32) (INFO), 2021-04-09 00:41:51,332n\n",
      "\n",
      "Setting up loss (INFO), 2021-04-09 00:41:51,338n\n",
      "\n",
      "Preparing data (INFO), 2021-04-09 00:42:32,417n\n",
      "\n",
      "Initializing variables (INFO), 2021-04-09 00:42:32,443n\n",
      "\n",
      "Model variables contain 0 total parameters. (INFO), 2021-04-09 00:42:36,418n\n",
      "\n",
      "Validation (000000): Learning_Rate: 0.0005, GT_obs_16: 399498.75, Loss_reg_unscaled: 1.2424506, Loss_reg_scale: 1.0, Loss: 798997.5 (INFO), 2021-04-09 00:42:59,618n\n",
      "\n"
     ]
    }
   ],
   "source": [
    "staggered_app = ControlTraining(step_count, IncompressibleFluidPDE(domain, dt),\n",
    "                                datapath=data_path, val_range=val_range, train_range=train_range, trace_to_channel=lambda _: 'density',\n",
    "                                obs_loss_frames=[step_count], trainable_networks=['CFE', 'OP2', 'OP4', 'OP8', 'OP16'],\n",
    "                                sequence_class=StaggeredSequence, learning_rate=5e-4).prepare()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "tUzt45Xv30dp"
   },
   "source": [
    "The next cell initializes the networks using the supervised checkpoints and then trains all networks jointly. You can increase the number of optimization steps or execute the next cell multiple times to further increase performance.\n",
    "\n",
    "*Note: The next cell will run for some time. Optionally, you can skip this cell and load the pretrained networks instead with code in the cell below.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "oMTXnLqR30dq"
   },
   "outputs": [],
   "source": [
    "staggered_app.load_checkpoints(supervised_checkpoints)\n",
    "for i in range(1000):\n",
    "    staggered_app.progress()  # run staggered Optimization for one batch\n",
    "staggered_checkpoint = staggered_app.save_model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "xpLvDj5-30dq"
   },
   "outputs": [],
   "source": [
    "# staggered_checkpoint = {net: 'PDE-Control/networks/shapes/staggered/all_53750' for net in ['CFE', 'OP2', 'OP4', 'OP8', 'OP16']}\n",
    "# staggered_app.load_checkpoints(staggered_checkpoint)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "KkQxqh1u30dr"
   },
   "source": [
    "Now that the network is trained, we can infer some trajectories from the test set.\n",
    "(This corresponds to Fig 5b and 18b of the [original paper](https://openreview.net/pdf?id=HyeSin4FPB).)\n",
    "\n",
    "The following cell takes the first one hundred configurations, i.e. our test set as defined by `test_range`, and let's the network infer solutions for the corresponding inverse problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "8AgfnvE830dr"
   },
   "outputs": [],
   "source": [
    "states = staggered_app.infer_all_frames(test_range)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rDEPL8E9fiFt"
   },
   "source": [
    "Via the index list `batches` below, you can choose to display some of the solutions. Each row shows a temporal sequence starting with the initial condition, and evolving the simulation with the NN control forces for 16 time steps. The last step, at $t=16$ should match the target shown in the image on the far right."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 413
    },
    "id": "dJvecqOY30dr",
    "outputId": "27fbfd47-9d1f-4ded-8238-cadd905e3429"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9cAAAGRCAYAAABv6UBxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3xc5ZX/8c+ZolGXLBe5yAWDsQ0B02sg1ACBAMmypEMSEvLbJBtI2Syb3V92syH5bbI1mwYsoWRDQgIbYkgooaxpAUw3GIMb7pYsV8mqU57fH3dkxrIsjzTSPDOa7/v1mpekmbmaZ86ZMzPn3ufea845RERERERERGT4Qr4HICIiIiIiIlLs1FyLiIiIiIiI5EjNtYiIiIiIiEiO1FyLiIiIiIiI5EjNtYiIiIiIiEiO1FyLiIiIiIiI5EjNtYiIiIiIiEiOSqK5NrM1ZnbOMJY7ysxeNLPO9M+jRmN8pWA4OTCzQ81soZm1mtl2M3vIzOaO1hjHsuHWQMbyV5iZM7PPjOS4SkkO70NhM7vezDaZWbuZvWxm9aMxxrEsh/ifZWYvmVmbma02s6tHY3xjUQ4xv8nM3jKzlJl9coDbv2xmzemc3GJmsREZ8Bg0GjkwsyvT34nazGyDmX3fzCIjNugxZLRqION+j6Y/mxX//RjF96HZZvb79OfyVjP7/ogMeAzI9TtnsT52n5JorofDzMqAhcAvgHHA7cDC9PWSH/XAvcBcoBFYTJATySMzGwd8A1jqeywl6lvAKcDJQC3wCaDb64hKhJlFgXuAG4E64EPAv5nZAq8DG/teBT4PvNT/BjM7D7gOOBuYCcwmqBEZWfvNAVAJXAtMAE4kyMXX8je0kjBY/AEws48B0byNqPQM9j5UBjwMPAZMBpoI+gXJkZmFfY8hZ865MX0B/htIAV3AbuDrWS73XmAjYBnXrQPO9/2ciu0y3BwM8H8aAAeM9/2ciumSa/yBGwg+YBYBn/H9fIrxksP70Lj0/Q/2/RyK+ZJD/BvT7zmVGdc9D3zE93Mq9MtIvO8DTwGf7HfdL4HvZvx9NtDs+/kW4mW0cjDAfb4C3Of7+RbaZTTjT7CybzlwUvo9KuL7+RbiZRTfh64GnvT9/ArxMlDMgbuAZmAX8ARweMb9bwN+CtwPdADnAMcALwPt6WV/DVyfscxFwCvATuBPwJEjle+RuIz5LdfOuU8QNMXvd85VO+e+b2Y7B7lcl170cGCJS2crbUn6ehmCHHLQ3+kEX6K25W/0xS+X+JvZCcBxBA22DFMOOTgCSACXpafBLjezL3h7IkVquPF3zrUAvwI+lZ6efzLB1tKn/D2b4jCC7/v9HU6wRanPq0CjmY0f6edQ7EYxB/2djmY27WOU4/9dgoakeVQGP0aMYg5OAtaY2QPpKeGLzOyI0XsmxWOgmAMPAHOASQQzAe7ot9hHge8ANQSzVO8haLobCD6DP9B3RzM7GrgF+BwwnmBm2b1mFtvPY+ddSe6j4ZzLZn/FaoI1LJl2ESRecpRlDvYwsybgxwRryCVH2cQ/PTXnJ8AXnXMpMxv9gZWQLGugiWALxaHAQQQfTo+a2XLn3MOjOb6xbgjvQb8CbgZ+kP77L5xz60dnVGPbUN/396P/Z3Pf7zWAVrwewAjlYA8z+zTBClgdjyMLIxF/MzsOOBW4huAzQoZghGqgCTgTuBh4lCAXC81snnOudwT+/5jinLul73cz+wdgh5nVOef63r8XOueeTt9+FEF/+p/pDZy/NbPFGf/uauBG59xz6b9vN7NvEKzweHyUn0pWxvyW6xzsJti/MVMtwRQFySMzmwj8EfiJc+5XvsdTQj5PMHvjWd8DKWFd6Z//6Jzrcs4tAe4E3udxTCXDzOYRxPsKoIxgq+nXzexCrwMrbf0/m/t+12dznpnZpcD/Ay5wzm31PZ5SYGYhgpXe1zjnEr7HU8K6gKeccw+km+l/IdiKOt/vsApPetbXP5nZKjNrA9akb5qQcbfMFdZTgY39Zg5n3j4T+GrmjANgenq5glAqzXVmgjCz3YNcvpG+21LgSNt7c92RaOrTcA0nB30H0/ojcK9z7jv5HvQYMpz4nw18ID0duZngoFr/amY/yvfgx4jh5GDJAMvu9X8ka8OJ/7uA5c65h5xzKefcW8AfgAvyPfgiNaz3/QNYCmQeUG4B0KLdhfZrNHKAmZ0P/BfB9MvXRnrQY8hIx7+WYKbAr9Ofy8+nr99gZqeN8NjHitGogSX9/6/sJTM2HwUuIdiXug6Ylb7e9nP/zcC0fv3X9Izf1wPfcc7VZ1wqMza+ec9LqUwLbyE4oigAzrnqLJZZBCSBL5nZDcBn09c/NuKjKw1DzoGZ1QIPAU8754a7L5gEhlMDnwTKM/7+LXA38LMRHVnpGHIOnHOrzOxJ4G/N7Evp5T8MfGTURjl2DacGXgbmmNlZwP+ml78I0ClXsjOcmPcdiTdE8OUramblQK9zLgX8HLjNzO4ANgF/R7BvngxsxHOQroc7gA845xYP9n9kZONPsBtE5ha66QT7qB4LtI7QmMea0Xgf+gXB1tNzCD4bvgRsBZaN8NiLVWbMa4Aegt12KgmOFzCYZwj6ry+a2U+BC4ETCPoyCFbq3WNmjxC89iuBM4AnnHPt/R7bDx9HUcv3hWCNyTqCo8p9bQjLHQ28SDD94yXgaN/PpVgvw8kBcCXBGqgOgqmAfZcZvp9PsV2GWwP9/scidLTwvOcAmAY8mH7trwY+5/u5FOMlh/hfDrxOMO14A/A9IOT7+RTDJYeYL0q/92dezsi4/SsEX6DagFuBmO/nWqiX0cgBQTOR6Pe5/IDv51qIl9GqgYz7zUJHC/eSA+CDwMr0+9AiMo6AXeqXfjH/O4LT6LYDawl2s3LAIen73kbGkcDT1x1HcDTw3QRHC/8t8H8zbj+fYNbGToIt3XcBNbnkeyQvlh6IiIiIiIiISMEws+eAG5xzt/oeSzZKZZ9rERERERERKWBm9h4zm2xmETO7kuCYVw/6Hle2SmWfaxERERERESlsc4HfAFUEu8Nd5pzb7HdI2dO0cBEREREREZEcaVq4iIiIiIiISI7UXIuIiIiIiIjkKK/7XJdZzJVTlc+HHHPa2bHVOTdxOMsq/rnLJf6gHIwE1YBfqgH/VAN+qQb8Uw34pRrwTzXg12Dxz6q5NrN64GbgXQTnJvs08Bbwa4Jz7K0BLnfO7Rjs/5RTxYl2dtYDl0Dc9bKMF9lNG0C1mZ2M4p9XfTloZ0e1mS1DNZBXqgH/VAN+qQb8Uw34pRrwTzXgl2qgcDzi7l67v9uynRb+A+BB59w8YAGwDLgOeNQ5Nwd4NP23jILlvMp4JnOKnQfwBop/3vXlAFiKaiDvVAP+qQb8Ug34pxrwSzXgn2rAL9VAcThgc21mdcDpwM8AnHO9zrmdwCXA7em73Q5cOlqDLGUJF2cHrUxlVt9VTvHPr/45UA3kl2rAP9WAX6oB/1QDfqkG/FMN+KUaKB7ZTAs/CGgFbjWzBcCLwDVAY8Y5x5qBxoEWNrOrgasByqnMecClposOyojxBi/Q7nYBzDSzKhT/vMnMAXCYmd2MaiBvVAP+qQb8Ug34pxrwSzXgn2rAL9VA8chmWngEOAb4qXPuaKCDflMOXHCy7AFPmO2cu8k5d5xz7rgosVzHW3IcKdrZSROzOcnOAUih+OdVZg4IpuGoBvJINeCfasAv1YB/qgG/VAP+qQb8Ug0Uj2ya6w3ABufcc+m/7yZotlvMbApA+ueW0RliaYtRSYwK6mx831U7UPzzaoAcqAbySDXgn2rAL9WAf6oBv1QD/qkG/FINFI8DNtfOuWZgvZnNTV91NsEaq3uBK9PXXQksHJURlriYlVNOBR2uve+qWhT/vBogB6qBPFIN+Kca8Es14J9qwC/VgH+qAb9UA8Uj2/Nc/yVwh5mVAauBTxE05r8xs6uAtcDlozNEmcvRvM5inEsBVADfRfHPq74cAIcRzOZQDeSRasA/1YBfqgH/VAN+qQb8K8oaMHvndzfgjOmioRooDuby+EKrtQan86rl5hF394vOueOGs6zin7tc4g/KwUhQDfilGvBPNeCXasC/gq+BzIbOQuBSRd/YZSroGuiLvQ0yOXcM5KPga2CMGyz+2W65FhERERGRwZjt29hZCIwx0dQVrIGa6mALb8Z9Qhk/lQsZHdkc0ExERERERAYzUGO91+2h9H1s//eRocumsd5nmZDyIKNCW65FRERERHJ1oCngZlg4jEs5tOV0hGSu0Biood4rxqk997WQ4VLagi0jT1uuRURERERy0bcVdLBGzTlcIhHcPRyGUFhbT3OxJ+apfRtr5/bNhXOQSu65r4XsnTyIjBBtuRYREZGR1b9h0Jah/BmoWbMQFnrnepdy7zQjyk3u+mKeSgIQmTKZxMxJdE6pILYjTuztVlItraS6u/fczxFO5yS9JTu9rAxR/9ev2YFf087hksk9dWEhwxFWDmREqLmWomKR4CXbt+ZX8qhv6pW+kPmTzZYRGZqBphQqvsOXnvaKhaCvmUul45mO716N3Z7lxt4RlfMu47Wc2UgPeNd0U9dnT04U/+FLHyQrfPAsVn5qMpF5bYyramVjexXRZ6cz9fE6Qm+sJtXZGdzfpYBwusFLBc2dcpC9UJjwuDpaL5lLZ6NR0epoWNZF+OXlpLq6smqwIZWeGp6WTWMu+1j396fQPS1OrDnCzG8+43s43qm5lsIXChOqKCc0rh4Xj0NXN8ndHVrDmE/pL2179lFyKX0IDVdfLMMZX2yTyQPup7dXA6jY5y7jNU04HDSAzt7ZFxL2jrG2xGbFwuFgJWg0ipVFsdoa6O7BdffgenpwySRmQayDBTIPQGRq8oYq83W8vwNpuVT6dZ1eJJQR+/TvfTlR/IfPwmFCtXVsPXUyH77oCb4x4RUAWpM9XFL+aXZtaqBhQ01Gc+1wKYeFUu802NoHOGvhhno6Tj6YsstbuGjKmyzaMoc1L0xlRtk8yl5bi9vdQao3Pvh3RZfxfm8hHUV8qEJhei44hms+tJD/U7+RO9rHc9PiP6P8/hdL+ju6mmspaBaJEKqpITl3OlsWVFG3qpeK1duwrm5cCRduXplhkSiQ0QRCekuTcjBUVlZGqLISK49BWRSSKdzu3aQ6ugZusjO3BLqUVm6MhFA4aALDIQgHv7tEYp/Gr39D4jK3wCr++7BIBCsrwyrKseoqkhNq2X5ELZWtCSrWtmPNrbhd7emGLoyZQSjIAcnkXvEfcOt2f6Uc//T7wp4VGWY45yCZhFQKl0ztd6WdS5Fuyh3mLKPJDmMkB27w9rdfcCnnIIOFw4Tqaug69iAmfHItF9e+TIQIYQvRGA5x6pS3WTR1AvUTxkFzyzsLppI4Z1iYwRtsxX9vZsQPn8H0v1nOTTP+SMwi/N2E13lzTg9XH/0xEj+cTeXqnYSbt5LcsWPw/7Unhqm9G+y9bnvncfe/fIkxIzxxPPff+GMqQ2UAfKxmG5fd+GM+eMIlJDa3lGyDreZaCpZFy0gdP5+151TxZx98kntWHYl7uJayHVVYcwQX7/U9xLGvr7Hr2/KR+UVNjfXQmBGKxbCDZ7LrsHrilSFCSUf59iQVG3cT3tWB6+qGrm5SXd1Bo5G5j2QyHe8DNRyyt/TKIQuHgkYulY5fNBo0dIlEsEU1tfcXpP7Tavvf3ve/S/aLVaa+Ri8Ww8pjbLtwLq3HpzjsiHVc2vACd64+Bnu8gcbnYoRefBPX0wNmuD2zYfbdqvrOkXwHUMo1YEa4rharryM5robehnJCPUmi2zuxLdtwXd17VhQd6MBauOQ7oczYAr4n9n3nZYahneKo1JgRqqkhMXc6az4Q4uYZD3FINEnmV+yKcJxQHKy7Z9/3jT37/7q9ZxWQEWfFfy/u5CNZd0453578KClSpHCEMA6JRrh1/n9z4z+exkO/OYlpj1dgL3cG7zkH/Kf9tmIDe+Ug83rlgO4Lj+cPN/wnlaHyva6PWZRfP/dbzr/2Gqrufs7T6PxScy0FKVRZyfq/PIqaM1r43iF3cUxsCzXhbh6oOZzVM6Yw+38OIrR8Handu/XldhT1TV1+p7FTrIcrVF0NM6ex7qIGOg/rpr6+g3DIsW5rDQ1P1TPh1QjhTQlSzhGqrtozhbZvi94eykFWLFqGlcewWFkwOyCRCBrrcDj4MtvVNegUWOf2s6XonTuM0siLj0WihKoqYOJ4tpzeSOqS7Xzp4Oc4rXI55Zakc1YZv1x3CpXNFTQ0TyGxes2+zR0EDXffS32QKc7BzxKMvxmpUxew/dAKdjcZPeNTuGiK8uYY9SvKqX8thK3bhLMQMMSVn3356D9df7BzNpdiDjKFwoTKYySOn8umEytInbiL6991H8fFdlNpZYQtRNKl6HFxnmudRUVrCnbt5ztL3/6/SQaPfWbBlGj8I7NmsOq8Ks49/0UOiXYTIkoII2whcDAzUsaV4//Es2fOYnNiMtM3TyKxdn12/9y59MqPjK3Y/VcwKQes/dYp/OWf30d1v8a6T3WonCu+fR//cuwlHPQ3pbcPtprrfiIHzWT556YOebmaNTDxhtJ7AY04M0LV1bS/9zDqzmzm7w+5j/dWxoFqLqpZwvzyjXy79yJ2vzSOuk1VpDo6tQU1G6EwobIoVlUJoTCuuxvi8aDhsBAWjew1NZO+Zi7zg0WGzSIRbPoUNp3ZQOM5G3j3xFVMie4E4M0pU7h/8/HUbCinsjUS5CWVntKZcntvMS3RD/KsmBGqqCBUVwuxMlwkDKEQ1hvHdXRhLoVLBY20i2e5ZU8OyCIRQnU1uKkT2XbMOHad3cW1h/yJ91UtozFcRo9zHFLegitzuNABtvT3Px9tVvcrDRYto/fMI3n7gyEaZ7SyoG47k8vbSDlj6c4prJ41id3TG5j2v1FCa5pJtbXhenuHHqv+U2T3XK+Goo/FYsGuPePqaD1tCq2nJFgwbxUfnryY8yo3Um3lQaMHpHC0pxKsf2Myszf04Nrb9/+PM2MPin9/ZkRmNLHmI03MPn0NV0x4inILE86Yqh22EGFCzI4k+eysp7jh7NNZE5vOzN9ESK5el9005YEabFAO0jZedwrnvO9FvlA/+AqLq+s2seK8Z7m//RSavvunPI2uMKi57qdn1nhWXPHTIS936Yrz6LphFAZUKtJTN0N1NSTmNLHpkjg/POR+zqzopu+IpoeXVXB4WSdPT1vOg02nUFtdibWaer8sRKZNIT59PJ1Ty0mFjWhHirJdccK7e7F4EsKGJVJYRxdudweuJz3lPh5Pr0nPYv9H2a9w4yS2L2ggdO42bprzS2pCRtI5uh0cVb6O3884gu6GGBWxsmBfyb4pbINtNZI9W46ssiLYx3d8DR2NlSTLQ7gQhBKOiuZuwtsiWFcPdHcHU+91wKYRY7EYNE5g5+H1bD2jl28ecz/vqVjN1EiMCMEX36gFX2jDcaCrO+tT5QQPYPteV2rMCNVWs+b9YW4+92aOKWsnZEbKObpdik3jy3itqYm7Zx/L+vhBTA4b4beN1Padw999KnOK7D7Xl65QTQ02eSLxqXXsnF1O5PIt3D73fzimrJsKKyNslXvum3QpOl0vz/dMYsKLRtmaVhJ9p+IaTP8me8DbSotFywiPH8fW9zRx9p89z+cmPMEh0QhJFyLEvp+T1aFyPlC9lqPm38GvppzIoy0nM2nLNpLt7dnFMLPBHui2EmSRCKkTDue7V93GxVWdWS3zz5Nf5vyrlvC9RR8ntHhpyZzpR821FIRQLIbNmMbOoyey6/J2nj7+x0wKVxK28D73rQt30VtPcDAoGVx6P7o1n5jBCRe/xq0zniTukvS4ODfvmsd9m49k7ZYGUtti1KwKM255nMqV2wm17YbKCtz2HbiOLjUjObBYjNZzZtJ16S5eOfZXhK0agB4XpzMVJ0UX72raxLq62aRqyglVlJPcnQzWsGtWxn71HewwfvhMdsyroO0QSM3oIhTqIZUMk0oYritCxfpqqjdUUdmaILa9h8iWNlKt23DapSR3ZoQaJ7L1mAa2nJ7g+pN/xyVVG6mwyj1b7hIuyc5kJThIRYBY2dAeQzkiVFGBmz6ZGy64JT3lOEYIA4NaoCGcZHrkbY6auZ7HrprPDRPPY+qTM6l6vYzE+g3Df2DFfi8WiRA/bg6bTy7HjtvFl+c/wAerV1MdihGibM9rHkhPB0/wVjzCVxZ/iENf3EaqdevQHlDxD5gRmjmNzedO5rov/5JzKzZTHYoBEDHbK+6Z6kIVHBWDgyY+Q+OX2/jjklODU6F1dGT3uP1X8JV4PsKNk3jw7tuHvNzZFUnOvvt2LjzlYhJr1o3CyAqPmmspCG7+wbz9gVpOO28J/z7tUapD1fu9b2W4h0SFg4Qaj/0JlZdjTVNY85Ep/POVt7Cg7HEaQmVAGVELE7Uwf1G/gg/VvM7q2ZV8Z+1FLO+ZQc+4KOEFjUx4bTyVjy8Lpt2X6NEeR0KospLWjy5gxhUr+deZ9+xprAFChAilP7THxzpYXWHBPqe9vcHRqxX3AVkkQmh8A73zm9h8cjnu2DZOm/EKZ9e/weFlzZRbkiRGZypCm4vRnYrS6WI8tms+i7fMpGXtJKrWTGXmXZtIbtyc3YFuZF+hMOG5s2k5ZQI7z+rm74/9AxdVbaDCYnu+7CZdiqRzpFywL6SlgN74nqOzS3YsFqNnfAXjwx2UW4RoeqVz0qX2TIMdHwozrsxxUHQpyYuNG2eeTt3j05l8X4LE5mbPz6D4WSxG13sX0Pv5bfz1QfdybuUaGsIxIpTvt7l7obeMb7/9fuZ8v5vUijU6COtwHf8u3r6ghi9/+Hd8oGo7oUFiPpBqi/GFcW9x01+dStNP5xL+35eG9vgl3lQD7PzEydz1nX8G9v/d/EBufOKXXPqPf8X4m8f+LrRqrqUgrP5QLaee8TrXNj5Cdahi0Pu+1DaDyc+mcBs2732gJ3lHKISrjNE1I86Csq1MCVfu82EUsyiTwmFqQr1cO/1hfh47lZauGlatn0TV/+wkqcZ62CwSIVRfR+eJs5n88TVc2/QwTZG9X9cpUqSco8NFWLRiDtOXxQlvaCWhZm9Qofo6euc38fbFZcw/5m0umPQ6x5avYWaki7pQGWELPtaSzpGimyRdpJxjbnQL76sfx0tNs3jpXdN5K3koTY9U4l55w/MzKiKhMJGpk4nPmEDbQRVsP9yoftd2Pj3rFc6qXE217fs+06dic4TaVe2kduxUYz1ErreXaFsv0yNxIrzzPpIZ66DJhmpinFW1jLcPnciD2xZQv3IaETXXOQlVVmKzmtj44TjfO/hhTinfxIRwxZ6VHAN5sjvCdW99kPjCiUxauUSN9XCZsXNeNXZEG++vWk7Uht7c9R3o7JI5r7Fo+sk01NSQGmzfd9nLmu+czEcvepwZkeE31gBNkWo++KXHuPXQs5j99bHdYKu5Fq9C5eXsft8Cjjx1BX/R+BjzorFB7//Hzih/ens2c1btyn5qTwlyySSh7jikYFcqTGPYMdDXgLCFqLZyTixvIz7pWTbGG7g/cgRbjjmI2tVrtZv1MISqqrAZU2k7rIGN5zh+MeMPLCjrJWp7H1Wzbz+x1mQV5W9UULmmldSOnVpLfgBWU03H1DLGz93GZ6Y9yRFlzUwIh6m0fl92+x3suzaUYmpkJ3OiL3Ji1UquPXUSu1fXUPmqTqeVrciURtqOb2LrEWHic7s4btZazh//OieWr2FKuGKfxjqFY7eLc9vqE6lfkSK8fkt2+5zKvpxLrzwafItd2EJMj8Q5vfYtnmg6mN3TaqnP0xDHIotEsJnT2HT2BL541P2cVrGZ8aH9r0QCWNLbzd8s/xgdiyYx/cmtJPVdJSeVWxJs2V7BpHDlge88iNNr3uThqpODc8NL1mLbjdfapsLE3FdEv94+ldj2A5yJYwzIel6FmYXN7GUz+33674PM7DkzW2lmvzazIe5IJUPhnONZ9wjAITAG4m8WrA2ePpXIXzTz3RkLOSEWHfQDK+6S/PPa86l4uZLQ9vyvdXTBF/DDiqEGXE8PbN9JdHuEVfHx7E4NvjW0LlTBMWVbObUiOOJp4hPbCFXm9kE20oqhBiwSgUNm0HLaBJo/2MuPzv05J8TcgKeriFqYJI6N8XFMXBKH5q2kCnyrdSHUgCsvo7famFW3nVPKW2iKxKi22KBbkaBvRVKMGZEKFpS18YW5j9M2PUyoYvCZMoXEaw2EwnTPm8Km042T3vcaN570c/5l+r1cXr2BQ6Pl+8S/72BOz3ZPJPXABOpfaCHR3DJqw8uXvNdA+hziqViEyICrSPdVbVEOi21m3sQWdjcZhLJbrhjkuwZC9XXsOmI8jZeu4y/rVzMpXLXf7yl9p926ofUMehZOYvrvWki+sXwkh1MQ8loDzhF79FUanwjzYFclySzX+Ae7pLxzSeEotzg43jkbSpHKdw1M+dc/senHh/Bg5+Abvw7kj51RVt04j2n/NPaPHD6UQ9FeAyzL+Pt7wL875w4BdgBXjeTAZG/rWEEVNZlXFXX8wxMmED9xHis+08gfD/sth0arDrjMy70p1rzQxOTFXbi2/DfX61gB0JVxVWHnwDmi7UbYUnv27R3MpHAlE8MpjinfwMIjbsXNmEqofOBzGPpQDDUQmj2Tt66q5YzPPceTp/2QCyu799v0JV2KuHNsTdRS+dxqkrvaCn4LakHUQPNWatYnWNc2jkoLE7PBV8oNpNzCTI3uoGN6Cnf4waM00JHnswbCdbWse28Z5737Fb419QHOKI/TFKmmcj9bU7tcLy/3VHHt/Vcw9cFNpNZvGq2h5VXea8C5vmaGFNm9P8QswvRwio9Pfpba01uITJ2891HXi1i+a6DnyFlsfrfjv+bcecD3mR6X4LmeKItvPJopf1hPcsXqkRxKwch3Dbh4L3W/eJafvPs9dLleelyc+CAH++xrphMkSZCkxyXYlerm8y98jHEreoMjhhcxH58DNXc+yw/ffUbWKzcG8oMzz6P+52N7OnifrL6RmFkTcCFwc/pvA84C7k7f5Xbg0tEYoEC362Qrm5nGQUARx9+MyPQm1n7rFNbeMImDv7eMGy676YBbnHpcnMU9cdPHUowAACAASURBVD558zXMuq+L6NJ1JNva8jToQF8OgK1QJDlIJom2w4qeyXRmse902EKMD1UwM1JGXaiMrSeMIzSuMCYUFnINWCRCuL4OO/4I3v5OJf92wR18s/Fppgyyf1LcJblh10xOefQaHvmzY0hu217w+7cXSg0kd+6kaslGUr+ZyKfXXMiq+O6slw1biBBG1MK8p2IbTO2mq7ECi3qf9HBAPmvAIhGSc6dz7Lvf4qoJTww4BTxTZ6qXW3fN5VNPfIq5t+wi1bwlOOdykfNVA6mOTqIbt/Nyb7Bl9EDCFqI6FOPMila+fPAjbD1zRlHN0NiffNeARSLEa8O4mgTTDjAluTPVy70djXzhp59n0t1vkNi4ueBXlg6Hz8+BRMsWPnTqnzP/4f/DN1qOY3HPvk12X2OdIsX2ZA8r4wke6JzA+1+/gvH3VFK+vKWo8+LzcyDRsoWLT72UO9vHDWm5P3SWB0cK37BxNIZVkLJd3f8fwNd556R744Gdzrm+o5JsAKYNtKCZXW1mL5jZC3EKe8pjoVrOq8zhyMyrii7+ofJywnNms+z6SVx2yZN8f8H/cO2kRzkxNvi+SFuTHTzUWcdnX72CGQ+1EV22jtTOXXka9TuKMQepji4mP7WLH758Jvd3HMKOZOega3sh+FIWsygxixC+rJXuw6YVxPTwQo2/RcsIN02l89RDeesz5Vy/YCHvLm+hbpCD8m1O7OaqdWfyH7+/iBm/DePWbyqKD/uCyYFzJLduY8JTzbx111w+tvSTXL91Hk90c8DXNwSv8Qhhyi3C+XOW0XxyGEKFv1XPW/zNsIoKdh5axan1q5gZiQ+6QrTHxXmkq55frz+W6jdihFq2k+ruKYrX+IF4y0Eqievo4kfNZ9OazG7ZqIWptDKOim2i/ZJ2rKJwZiENV77j75JJou1J6AoTttB+t9rFXZLf7G7imy9dTNPDO4NZSAW+snS4vH4OOEdi7Xqa7gnz+9+dzKdf/iRPd0cHzEtLsod/23oaH33l0/z1fR/FbptI/QstpLZtH/rjFhDv8V+zjn//zoe5eMX5WS3y4bfP4lvXfyo4BdcY+AzI1gGbazO7CNjinHtxOA/gnLvJOXecc+64KLnN1y9FrW4TZcSotaGtKepTKPEPjW9gx7ETueHU/+a6Cc9zYWU388sqB9wXtc+WZAePd03hhg1nkHpqHPbGapLbtuf9SLPFmgMX74U3VlL3dDn3bVnAW/EY6xJdB16QoAH5/OxF9IyLQNjv/nqFHH8rj9E9eyKbT4rwf05ZxMVVO5gQ3v8uDp2pXn7d/i6een4+U59IUvXMqqLYoldoOXA9PSRXrWHaw1vpfnASP3vuNP5jw7k80V3G7lT3AaeuBVuwQ5xZt4zEzG6srLC3XPuMv0WihKqr2HWwMTe2iZrQ/mMVd0me6Y7xk/VnsuWVRia+2kuqffeYaDS810BPD8+uncWmZIzOVG9WK5JCGA0huPTgJQX/Gj8QL/F3jnBnglB3+pzt7BvzuEvyYg/815p3U/FsFW7ZqjHbRHivgbTqp1fR9HgXLK7jgbYj9+wu0bfVusfFuW/3fO5efDzhB+uZ9fte6h9ZTmrtRlJd2X0HKkSFEv/6nz/D2t/N5lPrThv0fp/bcDKv3zuPcbeVxlTwTNkcMu9U4GIzex9QDtQCPwDqzSySXlvSBJTO9v482sU2WtnMVnc/qeCNvYYijH9y6nhaTna8tzJO8DIa3O5UN3e1z+P2t0+i7YWJzHq83VsTkpkDYDYwnSLJgevpYfLDzbx6zAzuqT2WhkgHX2t4K6t9VOvDnYS7XXCua48KuQassoLt82LMOGU9V9W/QtQGb6yXxeEHi89m8jNG9RstJLZuy+Noh68ga8A5km8sZ+qGZsYfP4c3T57D355Sxz8eupB5ZTuoC4WJEpzTPZQ+dHjm+ZdDGEkMsyCPFPB+eN5qIBTGymO4hjqSczqZHGknwsBNWo+Ls6w3xTdXXsa2x6cw8+kuQo+/zFg54YDvGnDJJKn1Vbxx1DTqQ2/TEIozLrT39Py+lUqZ16WA1t4aCIWC/a6LtPHzVQOheJJQ3Ei6FN0uQYTwXvHdnerhPzZfxPZnJzPz6bbgYKJjlO8a6JPcuo2ytyI0jJvJg2vn842Jz1CdbhZTpNiZSvCjpWcw+64kZYuXkWpvH2C1SPEppO9Ck//9T7z95vEs+enDHFm273f6Jb3dLPvuEUz73dg/eNlADthcO+f+BvgbADM7A/iac+5jZnYXcBlwJ3AlsHAUx1myDrEjOIQjANjutvAST7QXbfyznHmZdCk+vuoSOv52CuOeeoVxrMjyMC6jIzMHj7i7VwMriikHyZVvM/OeCTz4+il0T3Rc8PHXmB+NHnBf9yPKtgQ587zVqZBrwDU20D0J5tc1D7pFr8fFuWnXofzowfM59KvPAlBMZ/ot5BpItrUReewlZj5TQeiXk/japZ+l/YgeZk7bxrvGbeb0ujeZFG6nMbyb8WHHuPRsmaiF+dZrFzH5nrKCnyroswZC9XXsmlfPsTNX0BBKELaBV47+b1c1X/3ZVcy86U0qdqzz/r4x0nzXgOuNM/m5FN+bdR5nz1rOu2uXc0z5BiotOM1i2IwoRm3GbLAel+Dhzhms/fIhhFqWFm1jDf5qINTWRaSjFoC4S+31PSbpUqxNhHlu6cHMfD6Oe+H1kXzoguO7BjIlt+2gelkN25+dROcxSarTX2fiLsnaRCUVj1VTvnw9iQJeaTpUhfZdKPbA8/z1Eedy17KH95qF2uPiXHfs+6jYtjgfwyhIuZzs7a+BO83seuBl4GcjMyTJUlHF35au4tDbDuK2cydxefUGKvfTiCyPd/CZNz9O9ZejRFa/WehbPYomBxWPv0Hli1VQW81nl3+ZxJ9v4+Bx25hc3sYx1Wv5SM1GYhbdc/8eF+eKa79KzTMrC3mNr//4r1rPuDfreWj+fL7ZuIhYODrg3eY/+jmm3Bdl7pNvF1VTnQX/OQBwjlRXF27dBppub8MqKqAsyurodFbWHkr3pHK2zY+y++AEE2fsoHX9OCJtYaY9nqDyyWUk87yryQga3finkrjOTsq3x9kdjxF3QUPRf2vpf7dP5tsvXMjcm1eQ3L6jqJu4YchLDbh4LzUPvUHNiiZWxA5leXgetxh0NFXQOTFEogJSZdAzPkVoWifx9hiR7RGq1xqTX3qJVPG+xg9kdOPf3MrEJeN595I/50fzfsVaEjSEegkbvNE7ni8+dxXT/hii6o3msfbePhR5/xxw8V7c2g3MuAdOmXsNXznhEaaXbePJtkP5/YMncsh960i2tI72MAqFt8/h1O7d/PmFn+KsOxbzVw2ruGHnNBZ+5HRS29/M1xAK0pCaa+fcImBR+vfVwAkjPyS/Yiu3cNhPPz/k5arXOcYxuufwbLBJ4FgJxRf/VGcn4RXr+debL+PJy5bw+cbHODa2d4P9s12T+ZfXz6XhN5Uk33yhULd8tDvnLoLiykGqowM6O7Gdu5jQ3Uvb7mmsq2zg7ajxZNXxfG88dDfFaZy2g1l121m+bSJTnllDYvtO30PfS6HVQKqzk7rXd9JTN473117B1w5+mIOjrXS7CM3JOkKkuPbpjzD1vgh1z64fE+f5pVBrwDlcIhEceT2DRcuorK2mcs1EeidV0d3QwEHtSSLtPUTXtJDI85kHcpXvGnAdnZRtauPNF2byo5rTeV/9q0wNt7M+UcfNzafz/NLZVK+MMn1pnGRryXyZ9VIDqfZ2Qm+9TYg95xqmfnU9dTVVEI3gomGSVWX0josR7k4R6eggvKOTZHd3PoaXN/msgVRnJzWvtrDz1il8Yu61WEYHHemCpuVxqt7YTGrL1tEaQqHy/jmQ6ukhtH4Ts+6s55YXLyReDeFumPWndpItrcFxZ8aogvku5BypV5dx5w/fy62TjNg2mPRqaU4Fz5TLlusxKbF+A9O/vcH3MMakZFsb029eyp+qj+SFo6ZzxKTNhCzYNp1yIZ5eOofGRWGq73rW80jHKOdwPT0k1m+gcv3er3GLxXBHz2X7/Im8PmkSVZscya2rCnUFR+FwDlavozGeoIVJ/NWpl1Ff30E8GaZzdwwLwcG3poi+voJEgU89HqtcvDdouLdtJ7wUqggabhfvLeUtTVlL9cYJt7Qy9ckGFtoJ/GH64dRWdbOzvYLYS9XMfayN0KqVJHfs8D3UkpDq1ygnmlug+Z2/DfY6VJHewXPjEgkSb6+leu0G6sbV4To6cYkELpncM0ND7yOeOEeqs5Oyh16g8bEyQrXVWDRKornF666EpWjCjaV30LLBqLmWvEru3MWMfwjWavXfxnEoz+d/QAIEBz7j2SU0ZKzX0IdTdlKdnfDWSia8tZIJNw58H33BLSxjeYvGiEslSe7cRcXCxRycsSdf3/FqHXp9SwlIJfeZFSOFY89KVJECkO15rkVERERERERkP9Rci4iIiIiIiORIzbWIiIiIiIhIjtRci4iIiIiIiORIzbWIiIiIiIhIjtRci4iIiIiIiORIzbWIiIiIiIhIjtRci4iIiIiIiORIzbWIiIiIiIhIjtRci4iIiIiIiORIzbWIiIiIiIhIjtRci4iIiIiIiORIzbWIiIiIiIhIjtRci4iIiIiIiOQocqA7mNl04OdAI+CAm5xzPzCzBuDXwCxgDXC5c27H6A21NHW7TpbyPL10AwYwCUDxz5+9c8DhZnaNaiB/VAP+qQb8Ug34pxrwSzXgn2rAL9VA8chmy3UC+Kpz7jDgJOALZnYYcB3wqHNuDvBo+m8ZYYYxhyM52c7jeM4EmKT451dmDoBlqAbySjXgn2rAL9WAf6oBv1QD/qkG/FINFI8DNtfOuc3OuZfSv7cTFNQ04BLg9vTdbgcuHa1BlrKYVVBr4wCIWBSgC8U/rzJzAKRQDeSVasA/1YBfqgH/VAN+qQb8Uw34pRooHgecFp7JzGYBRwPPAY3Ouc3pm5oJpo0PtMzVwNUA5VQOd5wCdLkOgEoUf5/KUA14oxooCKoBj1QDBUE14JFqoCCoBjxSDRS2rA9oZmbVwP8A1zrn2jJvc845gv2x9+Gcu8k5d5xz7rgosZwGW8oSLsESngFYr/j7kXAJgINRDXihGvBPNeCXasA/1YBfqgH/VAN+qQYKX1bNtZlFCRrrO5xzv01f3WJmU9K3TwG2jM4QJeVSLOEZJjMDYGf6asU/j/pyAGxXDeSfasA/1YBfqgH/VAN+qQb8Uw34pRooDgdsrs3MgJ8By5xz/5Zx073AlenfrwQWjvzwxDnHG7xAFTXMtEMzb1L88yQzB0BLxk3KQR6oBvxTDfilGvBPNeCXasA/1YBfqoHikc2W61OBTwBnmdkr6cv7gH8CzjWzFcA56b9lhO1iG82sYwetPOseBjhM8c+vzBwQxF81kEeqAf9UA36pBvxTDfilGvBPNeCXaqB4WDA9Pz9qrcGdaGfn7fHGokfc3S86544bzrKKf+5yiT8oByNBNeCXasA/1YBfqgH/VAN+qQb8Uw34NVj8sz6gmYiIiIiIiIgMTM21iIiIiIiISI7UXIuIiIiIiIjkSM21iIiIiIiISI7UXIuIiIiIiIjkSM21iIiIiIiISI7UXIuIiIiIiIjkSM21iIiIiIiISI7UXIuIiIiIiIjkSM21iIiIiIiISI7UXIuIiIiIiIjkSM21iIiIiIiISI7UXIuIiIiIiIjkSM21iIiIiIiISI7UXIuIiIiIiIjkSM21iIiIiIiISI5yaq7N7Hwze8vMVprZdSM1KMmecuCX4u+fcuCX4u+fcuCX4u+fcuCX4u+fclA4ht1cm1kY+DFwAXAY8BEzO2ykBiYHphz4pfj7pxz4pfj7pxz4pfj7pxz4pfj7pxwUlly2XJ8ArHTOrXbO9QJ3ApeMzLAkS8qBX4q/f8qBX4q/f8qBX4q/f8qBX4q/f8pBAcmluZ4GrM/4e0P6Oskf5cAvxd8/5cAvxd8/5cAvxd8/5cAvxd8/5aCAREb7AczsauDq9J89j7i7Xx/txyxQE4CtI/B/Zg7lzv3iv/sRd/e2ERpHMRqJHAwp/qAayKAa8E814F/ec6Aa2ItqwD/VgF+qAf9UA36Navxzaa43AtMz/m5KX7cX59xNwE0AZvaCc+64HB6zaI3Scz9gDjLjP4rjKAqj8NxVA0OgGvBPNeCfjxyoBt6hGvBPNeCXasA/1YBfo/3cc5kW/jwwx8wOMrMy4MPAvSMzLMmScuCX4u+fcuCX4u+fcuCX4u+fcuCX4u+fclBAhr3l2jmXMLMvAg8BYeAW59zSERuZHJBy4Jfi759y4Jfi759y4Jfi759y4Jfi759yUFhy2ufaOXc/cP8QFrnpwHcZs0bluSsHQzLiz13xHxLVgH+qAf+UA78Uf/+UA78Uf/+UA79G9bmbc240/7+IiIiIiIjImJfLPtciIiIiIiIigpprERERERERkZzlrbk2s/PN7C0zW2lm1+XrcX0xszVm9pqZvWJmL6SvazCzh81sRfrnuDyOR/H3GP/04ysHqoG8Ug78Uvz9Uw78Uvz9Uw78Uvz9y3cO8tJcm1kY+DFwAXAY8BEzOywfj+3Zmc65ozLOpXYd8Khzbg7waPrvUaf4+40/KAe+c1DC8QflwDfF3z/lwC/F3z/lwC/F37+85SBfW65PAFY651Y753qBO4FL8vTYheQS4Pb077cDl+bpcRX/gK/4g3LQRzXgn3Lgl+Lvn3Lgl+Lvn3Lgl+Lv36jlIF/N9TRgfcbfG9LXjWUO+KOZvWhmV6eva3TObU7/3gw05mksin/AV/xBOfCdg1KMPygHvin+/ikHfin+/ikHfin+/uU1Bzmd51oG9W7n3EYzmwQ8bGZvZt7onHNmpvOgjR7F3z/lwD/lwC/F3z/lwC/F3z/lwC/F37+85iBfW643AtMz/m5KXzdmOec2pn9uAe4hmIrRYmZTANI/t+RpOIq/3/iDcuA7ByUXf1AOfFP8/VMO/FL8/VMO/FL8/ct3DvLVXD8PzDGzg8ysDPgwcG+eHjvvzKzKzGr6fgfeC7xO8JyvTN/tSmBhnoak+PuNPygHvnNQUvEH5cA3xd8/5cAvxd8/5cAvxd8/HznIy7Rw51zCzL4IPASEgVucc0vz8dieNAL3mBkEMf6lc+5BM3se+I2ZXQWsBS7Px2AUf7/xB+UA1YAPyoFfir9/yoFfir9/yoFfir9/ec+BOadp/iIiIiIiIiK5yNe0cBEREREREZExS821iIiIiIiISI7UXIuIiIiIiIjkSM21iIiIiIiISI7UXIuIiIiIiIjkSM21iIiIiIiISI7UXIuIiIiIiIjkSM21iIiIiIiISI7UXIuIiIiIiIjkSM21iIiIiIiISI7UXIuIiIiIiIjkSM21iIiIiIiISI7UXIuIiIiIiIjkSM21iIiIiIiISI7UXIuIiIiIiIjkSM21iIiIiIiISI7UXIuIiIiIiIjkSM21iIiIiIiISI7UXIuIiIiIiIjkSM21iIiIiIiISI7UXIuIiIiIiIjkSM21iIiIiIiISI7UXIuIiIiIiIjkSM21iIiIiIiISI7UXIuIiIiIiIjkqCSaazNbY2bnDGO5o8zsRTPrTP88ajTGVwqGkwMzO9TMFppZq5ltN7OHzGzuaI1xLBtuDWQsf4WZOTP7zEiOq5Tk8D4UNrPrzWyTmbWb2ctmVj8aYxzLcoj/WWb2kpm1mdlqM7t6NMY3FuUQ85vM7C0zS5nZJwe4/ctm1pzOyS1mFhuRAY9Bo5EDM7sy/Z2ozcw2mNn3zSwyYoMeQ0arBjLu92j6s1nx349RfB+abWa/T38ubzWz74/IgMeAXL9zFutj9ymJ5no4zKwMWAj8AhgH3A4sTF8v+VEP3AvMBRqBxQQ5kTwys3HAN4ClvsdSor4FnAKcDNQCnwC6vY6oRJhZFLgHuBGoAz4E/JuZLfA6sLHvVeDzwEv9bzCz84DrgLOBmcBsghqRkbXfHACVwLXABOBEglx8LX9DKwmDxR8AM/sYEM3biErPYO9DZcDDwGPAZKCJoF+QHJlZ2PcYcuacG9MX4L+BFNAF7Aa+nuVy7wU2ApZx3TrgfN/Pqdguw83BAP+nAXDAeN/PqZguucYfuIHgA2YR8Bnfz6cYLzm8D41L3/9g38+hmC85xL8x/Z5TmXHd88BHfD+nQr+MxPs+8BTwyX7X/RL4bsbfZwPNvp9vIV5GKwcD3OcrwH2+n2+hXUYz/gQr+5YDJ6XfoyK+n28hXkbxfehq4Enfz68QLwPFHLgLaAZ2AU8Ah2fc/zbgp8D9QAdwDnAM8DLQnl7218D1GctcBLwC7AT+BBw5UvkeicuY33LtnPsEQVP8fudctXPu+2a2c5DLdelFDweWuHS20pakr5chyCEH/Z1O8CVqW/5GX/xyib+ZnQAcR9BgyzDlkIMjgARwWXoa7HIz+4K3J1Kkhht/51wL8CvgU+np+ScTbC19yt+zKQ4j+L7f3+EEW5T6vAo0mtn4kX4OxW4Uc9Df6Whm0z5GOf7fJWhImkdl8GPEKObgJGCNmT2QnhK+yMyOGL1nUjwGijnwADAHmEQwE+COfot9FPgOUEMwS/Uegqa7geAz+AN9dzSzo4FbgM8B4wlmlt1rZrH9PHbeleQ+Gs65bPZXrCZYw5JpF0HiJUdZ5mAPM2sCfkywhlxylE3801NzfgJ80TmXMrPRH1gJybIGmgi2UBwKHETw4fSomS13zj08muMb64bwHvQr4GbgB+m//8I5t350RjW2DfV9fz/6fzb3/V4DaMXrAYxQDvYws08TrIDV8TiyMBLxN7PjgFOBawg+I2QIRqgGmoAzgYuBRwlysdDM5jnnekfg/48pzrlb+n43s38AdphZnXOu7/17oXPu6fTtRxH0p/+Z3sD5WzNbnPHvrgZudM49l/77djP7BsEKj8dH+alkZcxvuc7BboL9GzPVEkxRkDwys4nAH4GfOOd+5Xs8JeTzBLM3nvU9kBLWlf75j865LufcEuBO4H0ex1QyzGweQbyvAMoItpp+3cwu9Dqw0tb/s7nvd30255mZXQr8P+AC59xW3+MpBWYWIljpfY1zLuF7PCWsC3jKOfdAupn+F4KtqPP9DqvwpGd9/ZOZrTKzNmBN+qYJGXfLXGE9FdjYb+Zw5u0zga9mzjgApqeXKwil0lxnJggz2z3I5Rvpuy0FjrS9N9cdiaY+DddwctB3MK0/Avc6576T70GPIcOJ/9nAB9LTkZsJDqr1r2b2o3wPfowYTg6WDLDsXv9Hsjac+L8LWO6ce8g5l3LOvQX8Abgg34MvUsN63z+ApUDmAeUWAC3aXWi/RiMHmNn5wH8RTL98baQHPYaMdPxrCWYK/Dr9ufx8+voNZnbaCI99rBiNGljS///KXjJj81HgEoJ9qeuAWenrbT/33wxM69d/Tc/4fT3wHedcfcalMmPjm/e8lMq08BaCI4oC4JyrzmKZRUAS+JKZ3QB8Nn39YyM+utIw5ByYWS3wEPC0c264+4JJYDg18EmgPOPv3wJ3Az8b0ZGVjiHnwDm3ysyeBP7WzL6UXv7DwEdGbZRj13Bq4GVgjpmdBfxvevmLAJ1yJTvDiXnfkXhDBF++omZWDvQ651LAz4HbzOwOYBPwdwT75snARjwH6Xq4A/iAc27xYP9HRjb+BLtBZG6hm06wj+qxQOsIjXmsGY33oV8QbD09h+Cz4UvAVmDZCI+9WGXGvAboIdhtp5LgeAGDeYag//qimf0UuBA4gaAvg2Cl3j1m9gjBa78SOAN4wjnX3u+x/fBxFLV8XwjWmKwjOKrc14aw3NHAiwTTP14Cjvb9XIr1MpwcAFcSrIHqIJgK2HeZ4fv5FNtluDXQ738sQkcLz3sOgGnAg+nX/mrgc76fSzFecoj/5cDrBNOONwDfA0K+n08xXHKI+aL0e3/m5YyM279C8AWqDbgViPl+roV6GY0cEDQTiX6fyw/4fq6FeBmtGsi43yx0tHAvOQA+CKxMvw8tIuMI2KV+6RfzvyM4jW47sJZgNysHHJK+721kHAk8fd1xBEcD301wtPDfAv834/bzCWZt7CTY0n0XUJNLvkfyYumBiIiIiIiIiBQMM3sOuME5d6vvsWSjVPa5FhERERERkQJmZu8xs8lmFjGzKwmOefWg73Flq1T2uRYREREREZHCNhf4DVBFsDvcZc65zX6HlD1NCxcRERERERHJkaaFi4iIiIiIiOQor9PCyyzmyqnK50OOOe3s2OqcmzicZRX/3OUSf1AORoJqwC/VgH+qAb9UA/6pBvxSDfinGvBrsPhn1VybWT1wM/AugsOnfxp4C/g1wWkA1gCXO+d2DPZ/yqniRDs764FLIO56WcaL7KYNoNrMTkbxz6u+HLSzo9rMlqEayCvVgH+qAb9UA/6pBvxSDfinGvBLNVA4HnF3r93fbdlOC/8B8KBzbh6wgOAk6dcBjzrn5gCPpv+WUbCcVxnPZE6x8wDeQPHPu74cAEtRDeSdasA/1YBfqgH/VAN+qQb8Uw34pRooDgdsrs2sDjgd+BmAc67XObeT4CTdt6fvdjtw6WgNspQlXJwdtDKVWX1XOcU/v/rnQDWQX6oB/1QDfqkG/FMN+KUa8E814JdqoHhkMy38IKAVuNXMFgAvAtcAjRmHRW8GGkdniKWtiw7KiPEGL9DudgHMNLMqFP+8ycwBcJiZ3YxqIG9UA/6pBvxSDfinGvBLNeCfasAv1UDxyGZaeAQ4Bvipc+5ooIN+Uw5ccD6vAc/pZWZXm9kLZvZCnJ5cx1tyHCna2UkTsznJzgFIofjnVWYOCKbhqAbySDXgn2rAL9WAf6oBv1QD/qkG/FINFI9smusNwAbn3HPpv+8maLZbzGwKQPrnloEWds7d5Jw7zjl3XJTYSIy5pMSoJEYFdTa+76odKP55NUAOVAN5pBrwTzXgl2rAP9WAX6oB/1QDuKlbSgAAIABJREFUfqkGiscBm2vnXDOw3szmpq86m2CN1b3AlenrrgQWjsoIS1zMyimngg7X3ndVLYp/Xg2QA9VAHqkG/FMN+KUa8E814JdqwD/VgF+qgeKR7Xmu/xK4w8zKgNXApwga89+Y2VXAWuDy0RmizOVoXmcxzqUAKoDvovjnVV8OgMMIZnOoBvJINeCfasAv1YB/qgG/VAP+qQb8Ug0UBwum5+dHrTU4nVctN4+4u190zh03nGUV/9zlEn9QDkaCasAv1YB/qgG/VAP+qQb8Ug34pxrwa7D4Z3uea5H/z959x8l11ff/f33unbI72qIt0q56sVUsd1u4YmMbcAeT0EnAEPg6BQKEkITw+yZ8kwAhlBBCqKEZAhgwxUAMBhuMbdwkW5abLMmW1duqbG8z957fH3d2NZJX2tXOzlzt7Pupxzy0U/fM+eznznzuOfdcEREREREROQoV1yIiIiIiIiJFUnEtIiIiIiIiUiQV1yIiIiIiIiJFUnEtIiIiIiIiUiQV1yIiIiIiIiJFUnEtIiIiIiIiUiQV1yIiIiIiIiJFUnEtIiIiIiIiUiQV1yIiIiIiIiJFSsTdABERERERkROO2eHXnYunHTJpqLgWEREREREpZAZ2xCRfA1yoIluOSsW1iIiIiMiJ5MjCzoX5/1XUlYUZ5vsj3uVCD1CBLSNTcS0iIiIicqLw/EOFnZeflhxGhZwLAgiDmBo2BeR3ali+313oDu3YyO/siO7zo/sUCzmCimsRERERkRPBUGHtwqh4y9duQ8WeJRMQmAq7iVIwQ2Coj4ccVlhHNxw2m8A8w+ErDnIYFdciIiIicmxDCztpKmzpeH5U4LkwGqEu6GsXEsUgdJhnKuwmwlB/j2C4sH7B33uICw+NbJtnOGfKCxmm4lpEREROXKYvrmU3VEgfsZiTefkRU9CiTiVSWNR5mQyWqY7ikMsRdHRCGODwDxXYKuzGJ/83bqkUlkjgcrn8lPtDo9Xm+y/YyRH9fHiBjXngtJNDImMurs3MB1YDO5xz15vZIuAWoAl4BHizc26wNM0U5xwPcRfAyQDq//Jz0cZ1hZn9XDlQfsqB+CkH4jWlciC/mNBwITe0Qi/EWkhUfA6MNEW2sMD2DPOB0OVrifIW2BWdA8MzA0LM9/EaGsgtnUN/cwoMEr0h1RvbyG3eeqjA9n3M54UFYAlVTA6YR2LBXILmOnJVCfz+HDaQhRAsX2S7dAIvF+K27SLs7j7Ux4UF9lAMyrSzqZJywD95ESRHL0Xdlh2Evb1laNHEOJ6R6/cA64C6/PV/Az7tnLvFzL4IvB34wgS3T/K2spFp1NJN+9BN6v8y28pGgL6CmxSDMlIOxE85EK8pkwNmWCqFl04fGk2C4QWeXBBAEBxj2mbpVHQODB3r61m+eB4aiTs0ImfOIP8Y8xKHHlemqckVnQMFOzG8hfN45j0zueeGT9HopRhwOdZn07x9zVtY8PY6gvaOfIGdzwtzZRs5rZQcsGSC1u/u5/2t32VxMkngHHuCQaoMns7WsyPbwAXVW2gLqnnfh95J0+3rCfYfOPQCQwV2EL3WiKPcJVBJOfCW/72bN9QeHPVxl7zzT8n8+KEytGhieKM/BMxsLnAd8JX8dQOuAG7NP+Rm4FWlaKBAv+tlH7uYwyJA/R+HoRgA+0AxKDflQPyUA0cwOzTSVAZTJQcskcBLp7FEgqC7h7B/AJfNQRAQ9vbiggBLJLB0Gi+VxBLJssWhYnPALOr3VDIqrIMgXySEUdFccHG53HA8AMz38gX5yKcsmkgVnwP5mRleJkO2tZ7lp22j0UuRNJ+Ml+TkZD9/d+oduLmz8KqqDj3Hs0MxKHEuVFIOWCJB20ANbWGGfpcjbQnmJqpp9NNcmO7jD2u2syCRotXvxc86XP/AC1/E5XfwBUG0wrhf2hhUfA5UiDEV18B/AH8LDC2Z1wS0O+dy+evbgTkjPdHMbjKz1Wa2OssIf5gyqg2sZQlnFN6k/i8zxSBe6v+8oYKuzIUdTOEYmEUjeokElkwduqTyl2Qq+lI79MW2RLGZCv1viUQ0KpqXfdnZdL96JQNXnYMtPxkANziIGxjA5aK3bb4XFdhlKCwqLgZDf9u+D74fTfcdy4yAMP+YIMA5F8UgmSh5DCqu/wsVnP4pOONktl+R4aa595C0Q/ngYzQlunFpH7z81/d8zIYLbPMUgzFyg1mW1e5hYaKbjKXwzcPD8PBImk/akgD8rPs0kt3h8DbnhS/kCg5fKW2BXUn9X8lGnRZuZtcDe51zj5jZZcf7C5xzXwa+DFBnjVpx4Ti1uZ2kSFNnDRxwe4/7+er/4hXGgHH04JSMQf54SXy/6GmDyoE8z8errjr0pSoIokJj6AP/iAWICldALfaULVMqB4YWuRn6ouoZZkNTYT0Iw8Mfny8uoh/zby1fnEzUwk9TIgfyOzDwfawqDa0z2HJNAtc4iLc3TdMTDTRuqSXs6or+5kMXPTaZwDwP8z1cEL6w34cUOVWz0nIg2ikUjVhbVRpLJnH9/YQ9fcAYprYOHfNLAIlElCPJgmnihX/zE7DKeEXnQL5/vKo0tnAuOy+aRuslO3hJ1V48qghxBM7R5UJ+cfBMLBscth1yoYuOg/cMwwdnBVP6CygHDudCTqraS4ufPmwnhocRFrzBzf3N+IPh8GyNEYXB8OEShF5JpulXdA5UmLEcc30x8EozuxaoIjrm+jPAdDNL5PeWzAV2lK6ZU1cH+2ljF/vc7YTRcU+1qP+PrWAP8EScB7IwBsBiYB6KwdF5fjRdsyoNyVT0Jaynl3CQccVCOZCfKltbC80NkPAhF2B9A4QdndBXcOib70dfcj0v+tn3oi+7g4Pj7n+YOjlgyVTUZ/6h0byhUaHh4yFdGBV2hYZOyWIFOz6CAMt/mXVBEE2lHWeRPRVywHw/mgkwLUPQ2sTuF9fzkWu/y/LUbr598AJ+mDmPpodaYGPvoSnKQyv6VkXTyI2ory0/ojpsAo4LrpgcGDqefd7s6LMylSTMpMjWpUjv7sbbuTeafp/NP/5YfRYG0UrVFo1c29COqcCLpo/nTxk1pJjjUSsyB6yg71IprHUG269pZvZ1W/j6yd+jzssAEBIy4HI8PdjE7fefzSn7txIUjqLmt0nRtitaPXykcdNivw9VTA4M8X2m+73DI9QAvnkELjyswJ6V6uBxGHmHRQEXBJiXKNjZesSOWPOK2tFakTlQoUYtrp1zfw/8PUB+5Pr9zrk/MrMfAK8hWqHuRuC2ErZzyjrZTudkTgfggNvLo9zTpf4/Bs8nMWcWwcx6guokyb1dsLuNoLtn3B8qhTG40926CdioGIzA8/GmZfCaGyEMcd09hAcPHn0q1RhN9RywZAq/qYGec+ezZ2WSVCfU7Aipe6YdGxjEqqoOfWCbYek0mOEGBnH9/bj+gUMjSeM0FXLAn15PuGQ+uaoEXjYgsacD19sHOcMNZnH9fbhc9oVfjPI7k4YLcs8O/TxUZOdHvl1ufDGo+Bwww6urwc2aya5LG+m7rIsvnvtFTkt1kXWOk6r2UtXSw+DsOhKbfNzQttw5XHYw+ttOpYZHsvH9qNAemjJr+RHtIk5ZVBE5YIbf3MzAGfOp/9BWzqzfwdqOOWzrbCAIB+lc38jc39aRee4gtHdiZuT27D12nzmXz4vw0LT+ZDLaWTI8kyOMpo8z/gK74nLADC+dxmudSf/iGex5UZoVr1jPrfO+xaJEFUmrIXDh8Kj1/QONfPTZa1n+5YPkdu05/PtM/rhf56IZY8MzbWB4hDvq/0A5UMANHHtqtIcBHq+rW8P/Vl9OcgyzOQgMN7Rjm2S0bSpYqK6YFfYrLgcqWDHnuf474BYz+zCwBvjqxDRJxkj9X8ASCbyGBrouWUzdA1uwZzaTcI4wmxv5C/HEUAzyLJ2OzsdZMw3X1Y3r7SPsHyj1CrKV3/9m+E0N9J02ly2vgvNXPEMu9Hhs6zwG6hpo+WUXrid/eor8CJHr7Ts0XXx4YaKSzQCb9DHwp9eTO3URO87L0DMnBINUu8eMNVXUPLkb190T9ePRtiNhQDhw6AtU4UjdkBKuaj3p+x8zvOpq+l50Epv/EK48aw1va76PZckB0pZkgBzT/V7qMv0M1mVIVVcRHDEKPTxK6vuQn3mAc4eds7aEeTBpYuCvWMrW65t41423cUPNeqrMo6fhYbpCj44wzeZlzax52QLu2X0Se55bRMMTHjNv7iDs7z/2C+f7OopBODz7AxgurEu4LZo0/T/M80nMm83uq+dyYGWOc5c/z5+2rOKSqh00+9V42HBhHRLy+KDPvz57LX23tVD3/NoRP1ddEOTPVpcfwfYKCrrS9j9Mxhgch1rPCJM2fC7sYxmKA4mhEWzv0CyaIw+XmDgV3f+T0XEV1865u4G78z9vAs6b+CZVpq0fuoiabY7Grz0w7tdotJngeBbU/0P86fW4hXPoXFrHgRUeqYNQvasFL9sc3b/7ILntEzpDpss5dz1Ufgz8urr8aFsO19c38oeK50fHiWWqoxHTbDYasR7MlqSwnlI5YEZiwTz2XTKHtpcN8Nazfs9p1dvJOp/BMMFzMxbjqlK4g+3Do9ZAvhA8dPqiEnyQV0wO+EsWc3DlTHa/PMfFy59mTnU7fUGSzT1NPHnSHFprZjNtez/Jtm78/QcJDrSP/HftDh1fNzxBoHBBmwmMQUXlgOfjz2ii99wFbH1DwNvPup9rah/nlCQkLU1ISOAcAR6+F5Kr8qJtTV//odHrIUPTxJ0PwdD5sMNS7diYfDng+exb2Yhd0M4b656l3qsBoM6FtPiOAZdlcXI751Rt57xpm/hV06n8KnMqLXfPxduyfUwFNq4gBkeuPzCBhd1kzwFLJug/eSbtL+7ndaet4dXTV7EsmRuOSWFh3RYM8I7H/pTEb+uZfc9+gsLDgAoN7UyywvMwU8qdS5MvB8YgGGGGV0foCBPRrLBRZ+INLW4WhMPTw4Ho0JTQKQemiGJGrmUUflMjXZcuAeBPXnsHn3/wchq/FnOjKohfV0f2tEXsPTdD1zn9vPr0Ndy6eiXZmmlYCMlumPW7HOwY/zSoqcoSCbovX85gjUe6IySzrQtv5z7I5YZHRIFo+mUqiSWTUVHX20c4MKD+LpIlEniLF7D34pkcuKKfT5z/Qy6p3kVX6OhxCRbX7GNd/SLwvOFT4wyNmE7EOgNTgV9Xx/4LW9h3ZT+fPO+HXJPZR0hIV5hjZ5DikZkL+deB66jeliGzp5qa7Q1MW7N19GmyQ5QDR5c/7tef0UzvabPZerXHxy74PpdV76TBqyJpPlkXEDhH1jl6wjS9AymqDUgkhmdpvMBQgVFwXmbFIWLJBD1zjHNbt1Nj6eHbffPwIVod2QXUWMD06l3MSHRycGmGtoWLqN67H0YrrocUzhgYvq285yI/oXk+Xm0Nu89L85rT7uctDQ+wNJkiadXAocJ6yJODTfj31jP7zjaCdRuP/dqF5112NmELKk51g84jTEQ5NCYuisHwdqpwkT+ZElRcl4A3bRp4Hj0Xncy9n/vS8O2f5/IYW1VB8l/MwqXzef6Gas67aB3vnnUnF1T5/NHLHyTAeGJgLp/deBn2m/KerqgSWCKBP6OZJX//NK9tXsUP96/kN78/ndn31pE+mCXZ3o/1DeaPZTSsq4ews4uwr19FXbHMsEQSv7mRTW9o4ZV/cD83Nd7H/EQ1UMWg62PQhdQn+sjV51eMzRfT+tw+Dp5P7vTF5F67n6fP+XZ+QZsUADUezPRDliU3cfFV/8nmXAN3d57Cz587jZnfXED1Lw9Gx/nK+BQca7rvkjnsvWKQH1/2X5yVTgPTgKECI/qXBXYMNtC+v4a67hCXzR7z5VVEjMzLZPD7YU9vHb6NfBZWDyNpPjXAwkQ3L2tax5fnLyXzVAYOHhz7LyuYySGH86rShPNaeflrHuYDzffT4GeO+fhv7LmYxnVZ3JYxzsAbnkGgQYXxGFrQrNDshBEkC45jH01+J8fwae1AOzimGBXXJXDxA/t5X+Pj+PYbIDnq4+U4eD7+4vnsv7CVrhu6+MrZX+LsdA/1XrTX96x0mqwLSNpWljTuY8+8k8gcaCHYd0BfiMfq7FNY+ZVH+ZumR0hbknNm/4onbniQT59zJbu/vZCmtVkS7YOQC3B9+WOrgyC/OnvcjZ+cLJHAqqux1hl0nT6DPa/t546LPs7sRJq0RVMFsy6gyoyshXg4LDRcTfXRR/FkZGYkZrXwzGureM3cx+gKB0n7h2+nffOoJsWiRI4Wfz8LG+/npKq9fPyq61n+uyoCbUvGxdJp/IbpdL9oAdtfl+UVK1bx2oaHOT11qP8LF3Hy8NiYreeOHafQsCrJtAc3ELZ3jLpqr7xQ2N2DOTAb/Qu+b0ajl+ClmQ18+spuBtfNwG/bN+oCUDK68IyT2XBjFT9quZ+Md3hhXThqnXUBW3I5Nt68jNYnt5Lr7T2+X6RCrmiHYhGS7HWH1jcZC+eKXsxVJi8V1xNk4JoX8ZKP3Q/AnzU8QsabFnOLKpNXlabrjJl4b9rLF5f9kHNS/VQXTHGDaHrb4kSWC6Zv4lsLlpF5Sn/mY9X92vNpe3Ufn5v+ENWWwTePJq+a89M9fHThj3nTNW9n88I6qtpqSfQ56jZnqV69ibCjUx8kRQpOW8zOS6eRfvE+vn7qLcxNVA+fe3NoJC9jSdpcjo09M8ls9fH2d0aLO8mYWSJJ2FTHnFP2cHntOjI28g5Q3zyS+CQJqPUC5qX24zUOQCoZHU+tL6/HxautJbtyCRtvSPHGy37PNXVrWZDopdFL4VvqBY/3zdiZG+BfnnsdPXfPZN6du49+zLuMymUHCRNQkxy9QPbw8HA0eh6fPPsH/M1lf8L8/iWw+skytLRy+TNmsOvcGv7r5V8j473wbx6GVqiGjnCQ73ecR/Nj3YQH28vZzCnvyNHrEPAHXXTIm8gYqOoo0tYPXUS2LqR6cSf/NOOp/K0qrEtmyQIOLPN567zHOD+dJW1VIz6s1wV05DKg779j4/n4yxaz91yPd55+D3MT6eGpg7551FgVpyQD3rHs93wn/SL2H6ghHPRJ9CapHjqeSMbFkinc2cvYevU0ms7fzV8s+h3npd1wYV0oS8Cq/vk8tnsONdtDXEen+v54mGG+R1CT5kXNG1iS3E/ajj0t08eoMmO610vttH5oqMe6ujWKdzzMYMEc9pxbxasueZD3Nj1IvRedbuhIfv5csCFwe88p7Fgzi7lrB3HbdqqwLpIFkAuPPrX1UN87fDPSJDgntY/cad10bKqh/rHRV0uWo8uumEvHKQEvq+5ipFmNQ0VdjoAtuWq+t/EcFm1rI9C2pmTacnX0hnuOurMDwIPotHLaoSpjpOJ6HPyWmWSXzwHgC2/5IpdVay5sqVkigTe9nr0rp5M9s5uLMhvzx0mObHfgs767hap2hxsc1EISo/BSSbqXNeAv7uaqaU+PWHAkzefGuqdhIWyY2cqBwQxPPLecJt/Th844eVVV2NxZbL6ylpUvf5q3tdzL+ekekkfsNPLNI3SOrjDg9v1nMLCxjtmbeqPzt6vvj4/nEVQnOH3adlr8xFGPP4VoFMnDI0lIrTfIvOnt9LW0kmzbry+8x8FSKTpXTKf3rD7+fua9NPvH3gGdI2BnboDPr7uUllUhmad3kxvrglpyVKlOx9aO6ewNeph5lBgMFdjggYVkPJ/LF2/kvgVnM30sqyXLUR08uYqmhfuP+d0FYFM2yw8Ovhj/4TqCfRvU56VixvbBBrpdlgxHL66T5uFMh1/J2Km4Pg6WiLpr12tPZs0HPx9za6YQz8ef0cy+ly8i+Zq9fGHpT7i46uhfiAF25urZeHAGTZu6Cdr2a8RjFFZdTZgwwsDjQDjybACABj/DTfUb2JBZz697VlB7wwA7fjgLNF3z+Hk+tng+z7+2mc+8+b85v6qTGkvjH2U2RtJ8tuQyPLRuMXMfCrEH1pa5wRVg6Hy8vtHod1M9wnTkQkNTwwFa/Sx/POsBPrX4TTRvqYX2jnK0uCL4zU3sOR/++PSHRy2ssy4qrD+z7zLmfsLHf3odua6uMrW0ss1Y3cG25ib+c/4F/NOMtaPuWAKPjKV4deNqfj3nDCyTgZ6e8jW4wmTrjHk1R/9bDlxItxvgQ9tfwZO/WMai728np/UdSsarqeHAYEhX6Jh5jLXKqiwRrRY+hvNci4CK6+Pyl888yTmpfVTZPcCxpxJK8bxMBq+xgey8ZjZfPo3Tr3uGv5n9S85KJYgm6hzdf2x9GcH/NuEeeVAje2Pg+vrwBxxBzqMnTANHX5E346U4JZWlxV9LZvqTXHP2X9HQ2U1u1+7yNbgC+KeczJZXNvHzP/k4CxMZ/PypWI7mYNDLOx75U2bdmaBu1Xb0ET8+bnCQqkc2sTPbQI4O/FG2JUnzh1dRvmHaPj48wwjra7C9aU0NHyPX00PYPMjlNeuO+biOsI9P7DuP79x7Ecs//Dy0PUWgnXYTZ+MWWhuq+Pb8C7np+geY42deUGAPHWs6dIquwIVk8XEJB9Nroa0thoZXADNCHxIWErhwxH7vDPs596d/xdxfORb+fj25/QdiauwU4PmwYA4D4U56XGLEmHgYIY7eMMsY1gEUGabiehQH3nYhJ71jPQAvqWqnxnvhMWIy8fzp9ex99QoOnuqYtriDP1r8INfXrWVxgmPubd+V6+bS7/wNs+8LmPXEdnIqrMfEUimy0zymT+/hRekORtt5lMCn0U/j4dH1+k4SffOp/lmbRq/HyMtk2PHyJk66ehPzE9XH/JsG2Jrr5i3P/DFzPp8k/cxmggPHcVocOZxzhB2dfOrnr+Txy9fw1zPvzO/cOHoMomn5AW3BAN0LA4K6NJ6OdR8zNxjtrEtajpF2jAYuZHOul+se+nOq7q1l6QOdhPsPaHsywdzAAIn2Aar2VnEgSNLs5cjmT5k1dMaNIUNFdo6AjQOtJDp9rPs4V6yWQ8wjSENDujc6pr3grgGXZUtukKt+9V4W/yhH1VPbyR1o18BACZnv03NSHRdkHqfRyw3PGDvyNFwhIT/qXozf72CUz2mRISquR1F9MGBdWwtrz/sucPTpsjKxcisWcuDSAV552uNcUfc0l1Tto86rOuoX4L1BD7/unc+HHnkri3/WR+LpLeQ6Osvc6snL5XKkD+Y42J/GG8OxRdGoRrT4ynuX/4b/mveHTEslCfv1ZXgsLJWke2HIO2bfM+LCZUO6w37u76/lK7tvYN9vZjPv4bXkNC2zaC6XY/4dg/ym/2weOHMhV85/hj9qeJBlSf8Fx0Me+WUr0ePh9eW0kNxxcEGA6/dpC+oIXOdh2/HAhfS5QT6996VMu6uGlvv2E254XtMvS8AFAYn2buo21fGxnddwdt02OnLVtA3W0pNL0ZzupiHZy8xkJzMSnbQmOtg4MJv/euxyGp+KTuclRfDAw+VPFxpt97Mu4EAwwKr+Bcy60yf95BZye9tUWJeYJRN0zk9wSvUOpnuHSqEQR0iIl98J2BEOsrZnHv5gqLV7ZMxUXI+i+icPU7N+Cbf8pAGAF1dvY25Co9cl5fnsvnAaN551N+9sXEWDV41/lBV9B1yWdYMhP+48j289dj6n/GMb4Z42wiDEfB+nkY8xCXt7yazfy+D22awdrObi9AunSI3EN4831m7lXxc4Zs2fAxueK0NrK0A6jWsY5Jz0XmDk7cmAy7J2MMWntlzF9rvms+Bn+wn6tKjTREnc9QgnPTuPnhWt3HbBRay7rJU/bHmU5aldtPh9tPgp0hZ9ROYI6Hc5ukKPhqfBb2vXrJjjEQQkDyRY3bOI6zKPQcEUzAGX4+mszy/uPZtlv91L8OxmjViXinOE+w7Q+Gg1j9++nIdnLsXvN1LtRvU+x+ONRrbOMdgYkGzop7Whi+17Gpj7wwS1q7eQ6+6O+x1MWuYZloOObBUd4SBpSxDiOBj283S2ntvazqLh4V3R6ea0bSk936evxbE0uZdqSw2v0t4VDhICGfPx8GgP4cDgNLzBUDtUZcxUXI9BsG4jX1+2AIDP/vJyfn/Gj2JuUWXzqtIsfOUmXle/etTFb37Xl+E9//N/WPQ/u1jy7KO4urroXIT6cDpuua07mL5uDh8761o+u/j7nJQc206kjJfiPdfezqerrmbJX6q4HpUZbkYjXjKk1x19lsB9/VW8/8nX0v9oIy2PZQmeWl/GRk4NuS3bSG/ZxoJfwIDn8/1TLmfXZU30XNLNx8/9Ea/IRKOs+4M+Hh1o5PM7Lqf5V5vI7dkbd9MnFZfLMXN1yHfnreQDl60i68LhmUhrBhP82do/ZtlXDhBu2a7CusTCri548hnmHeuU1Z6PV5XGamtYnusg7Ogkp5kERavZ4VizZR6/mbmAazPb6HEhd/SczM1bL2T/3bNYyI64mzh1BAGZXUanSxMS0OsG2Z6De3uXkvEGODW9k2XJkBbf48y6bezyT9L3ShkzFdfHqf6P2rmu+jpyc5u440ffjLs5lccMfJ8dnXX0umP/eXaH/XzoH/6MhjDk4Hkt1DfUEKx+UhvA8QoDWu/cxf7e+bztzW/m68u/xaz86N1oo9hJy422xpwU8iDoT7AzV8vS5MhTzU5PdVL93em0/u8TuP4BnbK91MKAYN2ztD6/De+2Bj677PX8/SVpchlHssuo3eJo/t12cnt3ahszDrW/eIJpO07mvHXvI31etFBT75MNzFwdMu+ejQQHO1RYnyjCgLC3F/r6ouv6ey+ay+Vo/t12pu2cyScffD3/lobMnpDazf3Ubj/AtLa15Pr6lQNlEvb20vr1x/jwozfSPzONhTBtXRvW0weeh6tO42qqCTJJklv3kWx7Up/BJfCtl5zHtxKjl6LT2h6bVP2v4vo4BfnVG70DBzn3n/4cgJPevIHvL74rzmZVDucgm6WuaoCkHfv4lmpLsf9MY/YKeLG1AAAgAElEQVS9OTLr90F7J4G+BBQl3LmbpvtCOvpmceUl72fZmVs5tX4XDYlerqh5mtNSWaot9YJie033fNJ7jnEuCznEOWzPAaZtaOSLyy7jskW/OezuAZflyUHH6+57F0vWd0UjTVIe+aIi7B8g3dPLwrYWwqoklg3x27vJbdupL7/jFPb24j+zhYVdrfQ/FM2KaW1rx3bvH/5clROMPk8nVLinjarefmZtqQUzrLcf19lF0NurNQZiEPb24j/+LDWZ6LDDoHARRc/Hkgn8VCo6HEK5UBK53XvibkJJqLgep7C3l+YvPQDAurqLeNFLmkglcpoyPgFcLkf3QJqsG3kodFeum9/1zeMzm64gvd+o3tFD8OzzZW5lZQr7+wm3bKN2/0Gq957Mji0L2dSwgDAF/93yYmobelnYcJALGp7n9fWPcFKyhrv7PH711ArmPqGiY6zC9g5mPJZlzbRlXH9JmlPrd5HxBtnc18SG9hns2judOT9J4m3fjno1BmFAcPAgtLdjvg/mkQsCFdZFCto7oL2D5FPRdS0PJFNJ2N8P/f06ndkJJOzpGfnc7WGAGwh0ukUZFxXXE2D2J+6HT0Snj7rnUaiyLPP8AWZp4bNx27+pgVWLFzHbX0+tlyBtSXrdIHuCkFs7zuW/H7qU5e97hrruTYTaozixnCPs6sK7dw2t9x5xnxn9F5zB/1yxmC03NPK25vv4qyffwLzbfKp/8lAszZ2M3MAA6TvXcNLqenK3zeKOCxeRmwYNGwKmr9pBzbZHAVRYx805jSiJiIjImI1aXJvZPOCbQAvggC875z5jZo3A94CFwGbgdc65KX3y1aC9g48sPguAjZ87n01/8KWiX7Pf9fIUqxikHzCAmQCV3P8ul2PJux/i5je9go9fdD1nnbGJN7Q+zFe3XULnV+fS8OvnWNq2qmyjHofHgFPN7D1TNgecwx5Yy7wHYPNHjX+uv4KZHetLOmWqUnPA5XIE+/bDvv20rD50+4lYyikH4lWpOTCZKAfipRyIn3IgXsqByWMsSxDlgL92zq0ALgDeaWYrgA8AdznnlgB35a9L3vJ/fp5rX/56rrn2TRwMesf9OoaxhDO40K7iRVwOMHOq9P/0nz3F8n95noG3VPON615K4k+M6T97imDfvrK2ozAGwDqUAxHnCDo6S34s0lTOgROFciBeyoH4KQfipRyIn3IgXsqByWPUkWvn3C5gV/7nLjNbB8wBbgAuyz/sZuBu4O9K0spJKNizF/bsBc/n0s++n9CHBWuzx/06aasmTTUACUuCo48p0v9hVxecAIs5FcaA6DBB5cCQMkzJn8o5cKJQDsRLORA/5UC8lAPxUw7ESzkweRzXMddmthA4G3gIaMkX3gC7iaaNj/Scm4CbAKrIjLedk1cYMPvj90/IS/W5HoAM6v84pVAOxEY5cEJQDsRIOXBCUA7ESDlwQlAOxEg5cGIb85lpzawG+CHwXudcZ+F9zjkHI5+CzDn3ZefcSufcyiTpoho7leVcjsd5AGCb+j8eOZcDOAnlQCyUA/FTDsRLORA/5UC8lAPxUw7ESzlw4htTcW1mSaLC+tvOuaFzTe0xs1n5+2cBe0vTRAldyOM8QCvzAdrzN6v/y2goBsAB5UD5KQfipxyIl3IgfsqBeCkH4qcciJdyYHIYtbg2MwO+Cqxzzv17wV0/BW7M/3wjcNvEN0+cczzNaqZRywJbWniX+r9MCmMAFJ7xXjEoA+VA/JQD8VIOxE85EC/lQPyUA/FSDkweYxm5vhh4M3CFmT2Wv1wLfAx4uZltBF6Wvy4TrIP97GYrB2njQfdrgBXq//IqjAFR/ysHykg5ED/lQLyUA/FTDsRLORA/5UC8lAOTh7kyrPY7pM4a3fn20rL9vkp0p7v1EefcyvE8V/1fvGL6HxSDiaAciJdyIH7KgXgpB+KnHIiXciB+yoF4Hav/x7ygmYiIiIiIiIiMTMW1iIiIiIiISJFUXIuIiIiIiIgUScW1iIiIiIiISJFUXIuIiIiIiIgUScW1iIiIiIiISJFUXIuIiIiIiIgUScW1iIiIiIiISJFUXIuIiIiIiIgUScW1iIiIiIiISJFUXIuIiIiIiIgUScW1iIiIiIiISJFUXIuIiIiIiIgUScW1iIiIiIiISJFUXIuIiIiIiIgUScW1iIiIiIiISJGKKq7N7GozW29mz5rZByaqUTJ2ikG81P/xUwzipf6Pn2IQL/V//BSDeKn/46cYnDjGXVybmQ98DrgGWAG80cxWTFTDZHSKQbzU//FTDOKl/o+fYhAv9X/8FIN4qf/jpxicWIoZuT4PeNY5t8k5NwjcAtwwMc2SMVIM4qX+j59iEC/1f/wUg3ip/+OnGMRL/R8/xeAEUkxxPQfYVnB9e/62w5jZTWa22sxWZxko4tfJCEaNgfq/pJQD8VMOxEs5ED/lQLyUA/FTDsRLORA/5cAJJFHqX+Cc+zLwZQAz67rT3bq+1L/zBNUM7JuA11lwPA8+ov/b7nS39kxQOyajiYjBcfU/KAcKKAfipxyIX9ljoBw4jHIgfsqBeCkH4qcciFdJ+7+Y4noHMK/g+tz8bcey3jm3sojfOWmZ2eoSvPfjioFzbkaJ2jEplOC9KweOg3IgfsqB+MUdA+WAciBuccdAOaAciFvcMVAOlPa9FzMtfBWwxMwWmVkKeAPw04lployRYhAv9X/8FIN4qf/jpxjES/0fP8UgXur/+CkGJ5Bxj1w753Jm9i7gDsAHvuace2rCWiajUgzipf6Pn2IQL/V//BSDeKn/46cYxEv9Hz/F4MRS1DHXzrnbgduP4ylfLub3TXIlee+KwXGZ8Peu/j8uyoH4KQfipxjES/0fP8UgXur/+CkG8SrpezfnXClfX0RERERERKTiFXPMtYiIiIiIiIig4lpERERERESkaGUrrs3sajNbb2bPmtkHyvV742Jmm83sCTN7zMxW529rNLNfm9nG/P8NZWyP+j/G/s//fsVAOVBWikG81P/xUwzipf6Pn2IQL/V//Modg7IU12bmA58DrgFWAG80sxXl+N0xu9w5d1bBudQ+ANzlnFsC3JW/XnLq/3j7HxSDuGMwhfsfFIO4qf/jpxjES/0fP8UgXur/+JUtBuUauT4PeNY5t8k5NwjcAtxQpt99IrkBuDn/883Aq8r0e9X/kbj6HxSDIcqB+CkG8VL/x08xiJf6P36KQbzU//ErWQzKVVzPAbYVXN+ev62SOeBXZvaImd2Uv63FObcr//NuoKVMbVH/R+Lqf1AM4o7BVOx/UAzipv6Pn2IQL/V//BSDeKn/41fWGBR1nms5phc753aY2Uzg12b2TOGdzjlnZjoPWumo/+OnGMRPMYiX+j9+ikG81P/xUwzipf6PX1ljUK6R6x3AvILrc/O3VSzn3I78/3uBHxNNxdhjZrMA8v/vLVNz1P/x9j8oBnHHYMr1PygGcVP/x08xiJf6P36KQbzU//ErdwzKVVyvApaY2SIzSwFvAH5apt9ddmY2zcxqh34GrgSeJHrPN+YfdiNwW5mapP6Pt/9BMYg7BlOq/0ExiJv6P36KQbzU//FTDOKl/o9fHDEoy7Rw51zOzN4F3AH4wNecc0+V43fHpAX4sZlB1Mffcc790sxWAd83s7cDW4DXlaMx6v94+x8UA5QDcVAM4qX+j59iEC/1f/wUg3ip/+NX9hiYc5rmLyIiIiIiIlKMck0LFxEREREREalYKq5FREREREREiqTiWkRERERERKRIKq5FREREREREiqTiWkRERERERKRIKq5FREREREREiqTiWkRERERERKRIKq5FREREREREiqTiWkRERERERKRIKq5FREREREREiqTiWkRERERERKRIKq5FREREREREiqTiWkRERERERKRIKq5FREREREREiqTiWkRERERERKRIKq5FREREREREiqTiWkRERERERKRIKq5FREREREREiqTiWkRERERERKRIKq5FREREREREiqTiWkRERERERKRIKq5FREREREREiqTiWkRERERERKRIKq5FREREREREiqTiWkRERERERKRIU6K4NrPNZvaycTzvLDN7xMx68/+fVYr2TQXjiYGZLTWz28yszcwOmNkdZrasVG2sZOPNgYLnv8XMnJm9YyLbNZUUsR3yzezDZrbTzLrMbI2ZTS9FGytZEf1/hZk9amadZrbJzG4qRfsqURF9/mUzW29moZm9dYT7/8rMdudj8jUzS09IgytQKWJgZjfmvxN1mtl2M/u4mSUmrNEVpFQ5UPC4u/Kfzer/oyjhdmixmf08/7m8z8w+PiENrgDFfuecrL97yJQorsfDzFLAbcD/AA3AzcBt+dulPKYDPwWWAS3Aw0QxkTIyswbgg8BTcbdlivon4CLgQqAOeDPQH2uLpggzSwI/Br4E1AOvB/7dzM6MtWGVby3wF8CjR95hZlcBHwBeCiwAFhPliEyso8YAyADvBZqB84li8f7yNW1KOFb/A2BmfwQky9aiqedY26EU8GvgN0ArMJeoXpAimZkfdxuK5pyr6AvwLSAE+oBu4G/H+LwrgR2AFdy2Fbg67vc02S7jjcEIr9MIOKAp7vc0mS7F9j/wRaIPmLuBd8T9fibjpYjtUEP+8SfF/R4m86WI/m/Jb3MyBbetAt4Y93s60S8Tsd0H7gPeesRt3wE+WnD9pcDuuN/viXgpVQxGeMz7gJ/F/X5PtEsp+59oZ98G4IL8NioR9/s9ES8l3A7dBNwb9/s7ES8j9TnwA2A30AHcA5xa8PhvAF8Abgd6gJcB5wBrgK78c78HfLjgOdcDjwHtwP3AGRMV74m4VPzItXPuzURF8SucczXOuY+bWfsxLh/IP/VU4HGXj1be4/nb5TgUEYMjXUr0JWp/+Vo/+RXT/2Z2HrCSqMCWcSoiBqcDOeA1+WmwG8zsnbG9kUlqvP3vnNsDfBd4W356/oVEo6X3xfduJocJ3O4f6VSiEaUha4EWM2ua6Pcw2ZUwBke6FM1seoES9/9HiQqS3SVpfIUoYQwuADab2S/yU8LvNrPTS/dOJo+R+hz4BbAEmEk0E+DbRzztTcBHgFqiWao/Jiq6G4k+g/9g6IFmdjbwNeBPgSaimWU/NbP0UX532U3JYzScc2M5XrGGaA9LoQ6iwEuRxhiDYWY2F/gc0R5yKdJY+j8/NefzwLucc6GZlb5hU8gYc2Au0QjFUmAR0YfTXWa2wTn361K2r9Idxzbou8BXgM/kr/+5c25baVpV2Y53u38UR342D/1cC2jH6ygmKAbDzOxPiHbAaj2OMZiI/jezlcDFwHuIPiPkOExQDswFLgdeCdxFFIvbzGy5c25wAl6/ojjnvjb0s5n9P+CgmdU754a237c5536fv/8sovr0P/MDnD8ys4cLXu4m4EvOuYfy1282sw8S7fD4XYnfyphU/Mh1EbqJjm8sVEc0RUHKyMxmAL8CPu+c+27c7ZlC/oJo9saDcTdkCuvL///Pzrk+59zjwC3AtTG2acows+VE/f0WIEU0avq3ZnZdrA2b2o78bB76WZ/NZWZmrwL+FbjGObcv7vZMBWbmEe30fo9zLhd3e6awPuA+59wv8sX0J4lGUU+Jt1knnvysr4+Z2XNm1glszt/VXPCwwh3Ws4EdR8wcLrx/AfDXhTMOgHn5550QpkpxXRggzKz7GJcP5h/2FHCGHT5cdwaa+jRe44nB0GJavwJ+6pz7SLkbXUHG0/8vBf4gPx15N9GiWp8ys/8qd+MrxHhi8PgIzz3sdWTMxtP/pwEbnHN3OOdC59x64H+Ba8rd+ElqXNv9UTwFFC4odyawR4cLHVUpYoCZXQ38N9H0yycmutEVZKL7v45opsD38p/Lq/K3bzezSya47ZWiFDnw+JGvK4cp7Js3ATcQHUtdDyzM325HefwuYM4R9de8gp+3AR9xzk0vuGQKBt9ij8tUmRa+h2hFUQCcczVjeM7dQAC828y+CPyf/O2/mfDWTQ3HHQMzqwPuAH7vnBvvsWASGU8OvBWoKrj+I+BW4KsT2rKp47hj4Jx7zszuBf4/M3t3/vlvAN5YslZWrvHkwBpgiZldAfw2//zrAZ1yZWzG0+dDK/F6RF++kmZWBQw650Lgm8A3zOzbwE7g/xIdmycjm/AY5PPh28AfOOcePtbryMT2P9FhEIUjdPOIjlE9F2iboDZXmlJsh/6HaPT0ZUSfDe8G9gHrJrjtk1Vhn9cCA0SH7WSI1gs4lgeI6q93mdkXgOuA84jqMoh26v3YzO4k+tvPAJcB9zjnuo743fGIYxW1cl+I9phsJVpV7v3H8byzgUeIpn88Cpwd93uZrJfxxAC4kWgPVA/RVMChy/y4389ku4w3B454jbvRauFljwEwB/hl/m9/E/Cncb+XyXgpov9fBzxJNO14O/BvgBf3+5kMlyL6/O78tr/wclnB/e8j+gLVCXwdSMf9Xk/USyliQFRM5I74XP5F3O/1RLyUKgcKHrcQrRYeSwyAPwSezW+H7qZgBeypfjmiz/8v0Wl0u4AtRIdZOeDk/GO/QcFK4PnbVhKtBt5NtFr4j4B/KLj/aqJZG+1EI90/AGqLifdEXizfEBEREREREZEThpk9BHzROff1uNsyFlPlmGsRERERERE5gZnZS8ys1cwSZnYj0ZpXv4y7XWM1VY65FhERERERkRPbMuD7wDSiw+Fe45zbFW+Txk7TwkVERERERESKpGnhIiIiIiIiIkUq67TwlKVdFdPK+SsrThcH9znnZoznuer/4hXT/6AYTATlQLyUA/FTDsRLORA/5UC8lAPxUw7E61j9P6bi2symA18BTiNaPv1PgPXA94hOA7AZeJ1z7uCxXqeKaZxvLx1zwyWSdYOs4xG66QSoMbMLUf+X1VAMujhYY2brUA6UlXIgfsqBeCkH4qcciJdyIH7KgXgpB04cd7pbtxztvrFOC/8M8Evn3HLgTKKTpH8AuMs5twS4K39dSmADa2milYvsKoCnUf+X3VAMgKdQDpSdciB+yoF4KQfipxyIl3IgfsqBeCkHJodRi2szqwcuBb4K4JwbdM61E52k++b8w24GXlWqRk5lOZflIG3MZuHQTU79X15HxkA5UF7KgfgpB+KlHIifciBeyoH4KQfipRyYPMYyLXwR0AZ83czOBB4B3gO0FCyLvhtoKU0Tp7Y+ekiR5mlW0+U6ABaY2TTU/2VTGANghZl9BeVA2SgH4qcciJdyIH7KgXgpB+KnHIiXcmDyGMu08ARwDvAF59zZQA9HTDlw0fm8Rjynl5ndZGarzWx1loFi2zvlOEK6aGcui7nAXgYQov4vq8IYEE3DUQ6UkXIgfsqBeCkH4qcciJdyIH7KgXgpByaPsRTX24HtzrmH8tdvJSq295jZLID8/3tHerJz7svOuZXOuZVJ0hPR5iklTYY01dRb09BNB1H/l9UIMVAOlJFyIH7KgXgpB+KnHIiXciB+yoF4KQcmj1GLa+fcbmCbmS3L3/RSoj1WPwVuzN92I3BbSVo4xaWtiiqq6XFdQzfVof4vqxFioBwoI+VA/JQD8VIOxE85EC/lQPyUA/FSDkweYz3P9V8C3zazFLAJeBtRYf59M3s7sAV4XWmaKMs4myd5GOdCgGrgo6j/y2ooBsAKotkcyoEyUg7ETzkQL+VA/JQD8VIOxE85EC/lwORg0fT88qizRqfzqhXnTnfrI865leN5rvq/eMX0PygGE0E5EC/lQPyUA/FSDsRPORAv5UD8lAPxOlb/j/U81yIiIiIiIiJyFCquRURERERERIqk4lpERERERESkSCquRURERERERIqk4lpERERERESkSCquRURERERERIqk4lpERERERESkSIm4GyAyESyZwptej02rJtixG5cdjLtJIiIiIiIyhai4Pk6WTmOJo3Sbc4S9veVt0FRnhiWS+HNaaT9vNu0neSz8VkBux05wLu7WiYiIiIjIFKHi+jjt/+ECfnvWN0e8b1MO/mbxxRAGZW7V1OWfsoRn3tnAX19+O2+qvYUeF3Jp6/tZ/JMW/N8+GnfzRERERERkitAx12PgL1nMGY8aZzxqfOO0m6nxqka8nJJMctYjAWc8aux+70VxN7tyeT5+cxM7//Yiuj+d5V9editvql1Pg59hlp/hP679JjsuqcKrrY27pSIiIiIiMkVo5HoMXCbNJ1rX5K9VH/VxSfP5t5bHAPhJ6wVlaNnU5C+ez56XtnL6Det496w7OSU1SL2Xie4zjxel9zLQHOA1NRB2dcXc2snNEonoMAjPwwUBhA6Xy2rK/YnADL++DpoboasH19FJODCg2IiIiIjERMX1KBKzWjlwav1xPy/XMoh3xnLCx58pQaumIDO86mq8xgb2X9BC3Wt28p1FvwV8jtzhsTFXg+UMV52OpakVIz9DwNXV4FJJ/PYuXE8vYW8vbnBQRVxczPCnT8fNa6F7cR29TT71zw9StXEPbm8bbmAg7hZWFjOw/CQvHfIjIiIix6DiehTPv2MxT//554//eVd9lb8582weP6cEjZpqzPAyGYIzTmbTdRmWXfI8X1p0K1Az4sNv3vtipj9jBOs2lredFcZvauTgSxbRsdgjW+toWVVP3do9eIkErjcqslVgl5+XTtN+1TJOfs/T/GT+t/jU/tP4yt2XMee3c6i9b4CgrS3uJlYMSyTwMhlIJiB0BO3t+psXkRNLfmHXaHaZdgCKxE3FtZz4nCPs62fjTQnef/7PeX3tMzT7IxfWAPf/8gwWPNqJvgIXwfNZ/+9zecWK1byo5nlaE+3wenjHr97O7N961K/dh7dzj6bdl5nf1MjBq5ay58osn2i5m7QleW/jEzx85kI27D2JujXVoNp6Yng+/Veezc5LE3gDRtOTAfX3PE/Qtl9fYE8EZnjpNFaVJujsVkxkyvGqqrCTFtDxqSxJL2TL1mamP5pi1veeIdh/IO7miUxZKq7lhJVobWHglDnsWVnFSdc/x3/N+TbnpPbRkD+++mgangnxtuxBX7XGJ9Hawq5XLeaDK3/MJZnnaPSgynwAPvnSW/jHma9g/6kzqX92Bs2/30WwdTsul4u51ZXL0mm8ebPZ/MZZNFy8m/NnrOaS2g2clhoAqklbggsbN7F27gIdCjFBvGnTCM5aQtvbe3nN4idIegFPXD6bNVcuZNEP5lO9bjdh2z7C/v64mzolJVpb6Dl3Ptuv8GFWP/O+mSSzdhu5XbvjbppIWXinLWf3ZY3Uv3Inn17yA3wcT8yfy9fmXcyeYBnNX35Qs2xEYjLm4trMfGA1sMM5d72ZLQJuAZqAR4A3O+cGS9NMcc7xEHcBnAxQ6f3vTZtG29WL2bcy5KwzNvKjk3+dv+foI9YDLstf73wxtZt6CNs7JrxNLvqgWmFmP6/YHPB8wpZGBl7eyStqnqPJq8a3QycVeHVNJ12n/pqfzjiTx5fM4eCKWZz8qS6CfftL3rQpkwP5Y3wtmcCf3Urv0hm0nZ3i/Oue4J9m/4JGL0HakiQtWmvANw/PQnBgfaU93noq5IBfV0ewfAHPvbaKj5x2G+dXbaPKoL3e4zfNy/ikdxX1a+czfeMsMhv3ETy3uWxfYqdMDozEDL+2FmbNpHtpI7su9Lnhsod4Sd0zfPDxtzJvdyPsaSv5CPZUyIGhbZCXSmLV0XbGDQwQ9vXFXrBN6RyA4YUs953bQM/FPXxn6Xc4JRUNOCxObKZ9bobPLruaZvPAlSYXKjYHvGgQwW9qxFJJXC5HsGfvmJ5qyRTme4SD2bJsg6Z0DkwCxzNy/R5gHVCXv/5vwKedc7eY2ReBtwNfmOD2Sd5WNjKNWrppH7qpovvf5rSS+8MDfOO073Fp1eiPH3BZ1g7Cvd8+lzlbniXITvx2ZSsbAfoKbqq4GFgywWBThg+e9pMXFNZD3lq3l+un/YRn51ax7JIBXvfLd5JY3U/Y01PStk2VHDDfj0arZzSx/+JZ7Lk04C8v+iXva9zESDuXsi5gbec8qnckCPfuK2nbpkIOuPmz2XVJLd+94TOcm/LxLerzZhewuH4Tf3zVZ/mHMy/lZ4+fSfN9rcw42EFw4GBZio6pkgOH8fyoyKutJbd0DntWZuhaGrDi1M28q/ke5iaqec+yAfrXZEh7hgtL25ypkANedTVefR2uoY5cQwZnkNzfg7d1Z8m386OZkjkwxPPxG6eTXT6PtgsD3rpiFYuTyeG7M16SeckDhNNLO5OsInPADEsm8ObNpueUGQzWeCQGHLXPTIegYKPiWbTzaei2/PUwkyJI+vg9A3gdPYR72ko2s6lic8AMf+lJUZ8WKwgJNjxX/OuM05iKazObC1wHfAR4n5kZcAXwpvxDbgb+H5MxmJNAv+tlH7tYxCnsYRsV3/9mtJ8zgxe1rmVxoptjjVYPubOvlnfd/ccs/cwDBCX4kjsUA2Bf1MTKjIHL5vCCkBl+5zEf1+xPo9kHyFD/4W3s/Nzp1N7yYMnaNWVywPOxVApvRhPbb5hL/bW7+OSiu3jVtHbghTs6AB4eMB68fzkn39EVLTJXIlMhByyZYu/FDSx85SbOSycPuy9pPuCTtiT/OXsV75nxW7569kX8yruYmd9/iqDz2DlTrCmTA3DYCu2JBXPpXTaTtrOSXPGaVfyg9R4yliIkJHApPIzmGV1kaxoo9UERUyEHMCN3zlJ2X5Sh94w+0lVZzBzZDc3Mubue9F2PxXYY0JTKgSOZkZjZzOa3ncSp163nK7PuZmW6m7QdfraUbdlGqrakStaMSs0BSyTxZ7fS8M2DfGr2t1mQSOCb0e9yhM7hmZHExzcjcI60JfDNI3AhIS7/+QB7gx7evfUVtP3jqSTuemTC21nJOeDX1vLj33yXtCVHf/Aons9282eLXhLbWhxjHbn+D+Bvgdr89Sag3Tk3tIXdDsyZ4LZJ3gbWsoQzyJEduqki+3/onMrhYJZdL8/xj42PMss/9vkss9oAABwfSURBVPHVAN/onMm/rr2aBT+2ko0eDcXg4WgqDlRoDAgDknu7edeqN/HYJf9Nxkb/kH5J40a+0rqU6TNmlGyl6imTA8kENncWW29o4e1vvZ0bap5kdiKNf5QPm95wkDff/pcsuj0Lj28oadumQg7YssV0LoY3N49+CsWFiQx/N+Mhrvz7J/nYmjdij2/E/f/t3Xl0XGeZ5/HvW4tKu2zZlrzJuxzbsQ04i21wQkJMNkISmh4gk6bDTKZDQ4cJpzlpMn2aPs2cHhjoPvQwDXTaHaDTh07CkJClyerEJCTYWezYjuN9tyVbtvalSqqqe+87f0gOMpFlyaWqVyr9Puf4WKpF99XznEd1n3vf+94szJg5Y7zUAKEw4YpyvMWzaFxRQseqbq67aAffnvRblhakKQ31NhO+NYSMJWk9WvZVMv9kd9abvvFQA613rCJxcwd/v/zH1EabSdkQx70KDi2uYuNV89m49jIW/qQJTp4m6IrntNEeNzXQT2TaVHqWzODUpTEKPtLMw8u/R004oDQUe+/SIADfBjT6SV5uuoiaFxNZayrytQZMOIQtLuSvZzzE3Ejhe83ygI1evxOrYRMi3O+pilABNUWtHJm6kOHfxPf8xmMNDNcfHbmK1tvKIDjubAwDnwrpxxhzE3DaWntBh2CMMXcZYzYbYzanGXv3X531TDuL/uXLw37f5Vv/E7/97sqMt99oT1BAjHIz8YLeP5biH66uwlu5mIavrORPLnuVDxQ0Dzgtub+uoIcNrYtgXwklb2enkMZTDjAG05kgtqWEpB3aTtOKosN0V1uoqszKkMZT/EM102laXUXtJ/fzh2XvMjNSdM6juIkgxf1ti5i53lK4qy6r97ceLzkw1hLELPMKzn+dXdiEqAgVsbowyYE/L6D7ug8SKhzCNSwXYLzEH/oui5hQTvOyYiLXN3Hvihf4ypQNfCgWUBr6XXzPfDbU+WmqNkP0WHYviRgPOTDRAhqvTPPFi15jTWE7cyLFzItGWRFr46bSvdwz9UWu+eg29t05ieZPLsEsXZizsY2H+EPvSYZweTmRObPwr1rB8dvmcfgOy/Kbd/OtJU9wcbSAieHi95o/6G2st6R8Prvrjzn4H/OJ7NW+0HAFqTQ0tfLXdZ/kpN99/jecw/OJCp7av4zKzSP/9yif4w8QJBKs+l/38ET8/LNVB9OVjuEddddYw9DOXH8EuNkYcyNQSO81198HJhhjIn1HS2YC9QO92Vq7DlgHUG4qx9zShXbLTuZ6i+FPhve+1h2TmffIpoy3304zjZykyT5D0Lv+dRn5Gv9ohERVAekrOvjjCVuoGuSstW8DjngJvt94Na9tWcz0bUHWVortnwNgHlBDvuYAsF1dVG1N0jPEixdro92kJ/h4E4oYgStl3me81ECosJCui6fQ+GGPB2b9BzMj5/6A8W3AUc/jB7/+OIs3H8cb4qIrF2q81IDpShBOVNLolQNDm2IfM1G+d/nPuW/XF5i7YwpBFj7Ux0sNAIQnVxJfVEXLqhTfrl3P2uI6JodLBnxtQMCG+EIq9nURZHlRxfFQA6HKCdTMaGZl8YH3DmSECREL9x7gmxJO86UpL+OtCbNh4iK8ogqmNlTjNZzK+tjGRQ2EwoSnVpOumUzHvCJalhpmXFrPPTPf4MNFh1gQjRE24bPekrY+O1Mef7H/s3Sun8qMDe1Zuw1XXtdA4BO0tLF9/aX82x+c4PaKzcyNDq/JS9o0D51eRWRbKf7+nSM+xHyvAet5VP1oI9u/OItbS3Zd0M/4ZuMS3n1zHvPJ/t+kwZz3zLW19n9Ya2daa+cAnwM2WGtvB34N/GHfy+4AnszaKB0zXsC2ZJJtySTtweBHtN5J9bAtmSTSMzJtxgKzjCvMJ1hjbmQpKwE68zX+titONBGwfNoJpoRj5zxr7duA1qCbf2y6ilf/9TIueqCTsqe3Z21c/XMAHCKfa8Ba/I4uYm8fImF7Y30+VeESwuVpUhWZXyczkHFRA8ZgaqbTsDLMPWvWs7zg3GdAfRtw0k/weMeHWPSNPXgnTmZ9Ma1xUQPGEJxqpOwoPNawYlhvvaG4k+45aVJzJvdeLzzCxkUNABhDauE0Tnw0wvfXPMynS5vO2VgDdAYpHjy6inBzJ0EWZ27AOKiBUJhg+hR8a0gEA1+9HjNRFkYNn5+ykT9ftZ6ua+LEV8zKyfDGQw2ES0uIL5tO3cdK6PhUF1+5+RmeW/Iod1Y0sLjg7LPVcGZfqIf/efwmkj+dSs0v67Hb92Tt8yDfa8CmU8x74Cg/3nglj3Z+YEj7P2ekrU+dl+SNHQuYtqknK9Pyx0MNABxMTOa0f2ELJ/7i4auYf2/mJzYzdd7mehBfp3dxswP0zvn/8cgMafTxd+7l63NX8vW5K7nx3dvP+bpjXhd/sWANX5+7kll/szHbw8q7+PtNzRTvb+GNvfPoDFLn/MN2zEvwf5pX8dZ3LmX640ewuw+6ut9s3uUAgMDHJpP02PD5X9vnU4u3ceKK4dx8YETkR/yNITJjOge/UM3aj2/lrorBr53ekUrz+b238+TffQy/vcP1rXHyIwcA1mI9j+pfHqDt/ll0BUP/m5K2PnPmnKZpWdH5Xzyy8if+QHjyZOrXFHHbDb/hhuLO9zUTv6/RNySer8a2trmsg7zIgQkZeqYWM6koQXEoec7P3+JQAR8o6OK60l3cXLuD+isj793CyJG8iD/GwIxq6taGWX7jHu5f8TP+dMKhQRd3OuknWNd6Cd13T6Li0bfxDh91tYBTfuQA8OrqWfSDDn705tU0B93nbbB9G5AIUryT8vnkW3/KjPWG6JvnX7NjhOVN/AFOre5g9WNfcz2MjAxrb9ha+zLwct/Xh4DLR35Io9uEL/ncUPGfB34yCLBe9oqq0lSB5QDkb/xNooeigwX0rB14R2lbMsnf1t3KnqcXMnvzCfzTTVldRGgAndbamyCPcxCLYWbPIG4jePiEh3AMbnPzLEqPZX9seVkD1hK0tOLPqeIj5fuJmXP/Wf5Ocy0PPL2WGa94THpzb1ZWxh+CvK0B63n4zS1M2AA3/tl/J/XFFr46/yXWFB0fdJp+0no0tJUzsSnL94EiT2vAGEw4TNs18/GWxlldcmDQxtq3AXvSSb7XcC3Vb8QJunN+cDW/asAYTCxGfGqEubE4hcY/5yKKAKWhGDHjc035Ll5ZtoBwaQl+Z2fODnDkYw2YSJS6G6fwqas38YXKjSyIRogOkoPTfpy7D3+axh/MoWz31lzvB0G+1UA/dv9hJv12BZ+ZeTu/WvIIpebcM8k8fJqCFK/El2O2lFN8Ip6V2Uu/Lx9rIFNXfvkuZr22Fzfrg58t56eaxjrv8FHXQ8hrQXsH015Psve/VlAZ6j5rtep96Th/8JuvMOmVGDO3thM0nHbxgZL3jDEEhQU0eBUsibYM6bYIRw5XUbs9e7eBync27Q26X9rkx/lp+3IeeHYtNS+mKNx+rPfeyjLyAh+/uYXS13xaSxbyjeWfo2xJC3cu2Mh/qxj4TNK3Gj+C2V5GxZ42ArczCcYeY3pvP7dgDg2r4foFe6mNNnOuWzAe87rY2F3Dj+vWUP/rGuYeOYSXdnNrqLxiLX7fbPDADt4cRE2YEIbaaDOfmLmTp277KNUPZ/92dHkrFCZUUUbikgTXV7zD3Eh40M9d3wY8E5/L9j2zWPJ6PZ72g0aUTSap2tRMffl0Di2Ei6PBgJcppq1Pe5BiV2oSjxy9hCnb0kSPN+Gn0gP8VBmOWc96zC28i8M3rxvye0r3teFnee2NoVJzLaNK0NV7ve9X3/ks/6X2dRbEGkjbCMfTlTzdsIwZT0YpfW47QSKBdmGzw1pLqCfF94+u5d2qvRSG0jSlS0kEBUSNT4dXSFE4zaxYC9XRdnYmZjBhe5TIrn2j4ojhWGR9n6Ajyon0BAJOn3VrjyY/zouJmfzo1WtY8HQP0XePZm3BGunT12BPeHIHFbvn0HRJJX+/8npaV77K7RWbAUhj6LFhAmt49K1Lmb01DQfdrlA6JpkQoViM9osnMn9pPTdP3Ep1+OxdE98GdNkku1MFPNRyNS8cWkR4axmznmvvPcjk6F6mecNarO9jLMS9AlJDmK0UNiGqwxFuKN/OL2/4AOaZclBzfUFMqHfmwKo5R5gXaSdmBl/M9aDXzT/uv4qJWyN49SdzONLxw9+1j+qJxRxKT2ZxtJ2B5tEEBLQFsKtnBk2HK1m0/zR+Y85nU+alguc3Mz+1Am4+/2u7gh6+fPxaTMLJ5aEDUnMto4u1+K2tTP9UK//6tevpmu0TSoaYuBum/KaByIkdBAmdIc0mm0zi794P376EX1w0k2jCUnG4h2hjnKAkRriuETuhjE2XXUbHXMPMDd1M27Mfv63d9dDHLhtQcjTClvbZJCfsfG9KbNr6vJiYyd9sv4kl3zyKd+q0q6ng41IQj8PWnUx+t4DqZyfzxPVX0/rlYnxCNPSU09xTQtoPM/spKNlyDK+z0/WQx6ZYjKZlhjum7OOiaDNFv9dcNAfdbOyp5hvv3kz0hQpmb40TPrgfv6nJ9ZoD+cP3Kei0tCaL6QwKYQiHSotMAUujHv+8/Gd8s/KPMCcjOb3vdb6wvo9NpVlU2kDhOSYN+DYgwBIQ8HcN1xL7WSUVv3qHQAeWssb4AcdTk0gUNVLK2YvsnrkWu8UvZGfXdApPhjFdCaxm0YwYE1hO+3GqBlnUEmB7qoBTqzuA0XNwT821jFrTvrcJzvwxs4GaihyLvLSFqg19n/TWvrer5QGcbGDCngNMMKHes3yOxpg3rGXmC+1sqbyI/1tWz72TdhE1Yf7q9CU88fRqFvzz8Zzc7kYGZtMpvJOnqHoyyc5ti0lXFBJKB0RPdxKpO4lN1uGpqbgwNsD29BAUwPSCViaEQmftxNZ5Xfxl/Y289dxS5j7UQHDkANZL6/MgCwpbfCbGEkwKJfBt9H1TYX0b0G1TFJkCwqY3TzEiVIcTNKyZyLTu2fj7Djoa/dhlwmFMrIDiUIoe27uGQ/9L4tK29xM2adNsTxVQ//mplB9+m0BnSLMq1OPxSkstt5TtpDDs8/vTJRNBml91rOTVDctY8HgTfmsbDGOFcRlc6JWt3LHoWn6y+3mmDbLmyWik5lpGL2vBqm1zarAdWOVnRJm9h1l4/xReeXwl6yuuIJS2RDuSLGisw8/SPdxlGAKfoK2dkOcRA6wfEKRS2CzfAirvWQtB7w5pWajnrAX9fBvwrVNref3XFzP/yTaCY/Wacpkl1vcpOtTCsY6JHKueyMxII2WhgrOu/e1tpqN4+O+tLRAQUBYyDOPmEjIQYzjQXUVLSQGV4d99rp5prE/63TzZeTHrHvwENY27sZ6u6822UHucd16r5Yela7ir8jVmRc6+G4SP5UhiEuUHgea23rPWOug3oux5rl+/df91JP5yGoZtORrR0Ki5FhEZBYJEguDwUThiiIXDWM/D0jdTQEYF63n4HV29Zye0EzVyrMUvDSgJJQn3W2m3yyZ59u1l1Lzpw8HjOpCRTdbC6SbaNi/hu6Hr+WzNFj5aspfOIEy9N5EwlknhLmoiXVSHI8T6ndVOWEvxqQDTeWH3phXAWvZ3TCE+qQDoBs5MBQ8IEeKEV8QrLbVM3p7Cdnfr708O2NY2Zrw8hV9MX8Hlqw8xOXzqvRkFSZvm+cQsXj88l3l7urHtHVr7IQus73PFz+/lO7f8O58uff+072PtE5jy29HVWIOaaxGR0aXvfssySmkHauQFAXiGTfEFhNjHhHCCQuOzJzWVyW9EKHvrqK5nzwG/vYOZG5Kc6pjKPyxdy9OzltLWU0RzWykmFFBR2kNNeSvLKk4wP3aKqkgnxaEkL3V+kLLDcYL20XPN45jj+xxvrqRtVjFJGydpe8/Ypa1PwqbY1rOQd+qnM78hTqDVqHPCb2sn9vIOJs5ZweO1K5g/7Xku7put3xl4/KZ9IeEjhUT27MXvGT2LaeWVwGf+vZt4cs0H+XTpb8566rGuclrrKpjiaGiDUXMtIiIizljPY9YLPk8c+SgPVV+JV+kRLU/CoRIWbGrCqz/heojjg7WEX36b6S8DxmCBChNiQjSCMQbrB8TTKTbPnMOLH76C1otC9Ez1mfaKoeLQfvwezSy4ENb3sfEE3rESdi+eQVW4k7JQikLTeyDvpcRC1h1YQ+GbpbD3bR3gyyGbSjH1meO8OXcx667p4R+mvwpAwkJNYSt+zOK3ajHXbEsFEdLWf2+xV4Bvf/d2Fj6wyeGozk3NtYiIiDhjPY/Ys28z7YUwhAymb2q4tRY/pWusnTgz7dj62KR/1lpOXl09ZY83Uh4OYa0F38fXbJsLZy1BIsG8x7r5Sfda7p9yNRhLpD1C+QGoer2VqcdPEXQd1iJmuWYtXl0901+ZxnPTlnBf9YtEgbiN8PD+S5i40+hgRw60X5tk6X13s/fOf3I9lCFRcy0iIiJuBT62bydVV5OOfjadwmp28oixnkdk52Hmt1RhC3sXkTNpH9MRJ2hpxdd11u5YS7TLI9QY49GO5USNzw93XsnEx0uo3FindVFyIIjHCad+tx7HZX/1JapfGr2xV3MtIiIiIuKQ39EBHbpufTSKnuqg9Gg1Pzt0GeWFSUqfL6XyjZN4x+pcD23cmLzD4wNv3sb2yx+m6vmjo/pyITXXIiIiIiIiA7DH6pm8o4ITFZPoLLfUPnUAr7lFswlyqOiJNyndW8sjT0wc9Yu+qrkWEREREREZSDhMbE89c+qKscUx/OYWXWvtgL97Pz+9aDZw2vVQBqXmWkREREREZABBPE6QSGCawtjAqrGWQam5FhERERERORdrR/10ZBkdQq4HICIiIiIiIjLWqbkWERERERERydB5p4UbY2qAfwOq6b395Dpr7feNMZXAz4E5wBHgM9ba1uwNdXzqsQl28hYpegADUAWg+OfO2TngYmPMPaqB3FENuKcacEs14J5qwC3VgHuqAbdUA2PHUM5ce8DXrLVLgFXAnxljlgD3AS9Za2uBl/q+lxFmMNSynNXmOi7jaoAqxT+3+ucA2I1qIKdUA+6pBtxSDbinGnBLNeCeasAt1cDYcd7m2lp70lr7dt/XnfQW1AzgFuDBvpc9CNyarUGOZzFTRLmZCEDERAG6Ufxzqn8OgADVQE6pBtxTDbilGnBPNeCWasA91YBbqoGxY1irhRtj5gAfAt4Aqq21J/ueaqB32vhA77kLuAugkOILHacA3TYOUIzi71IBqgFnVAOjgmrAIdXAqKAacEg1MCqoBhxSDYxuQ17QzBhTCjwGfNVa29H/OWutpfd67Pex1q6z1l5qrb00SiyjwY5nnvV4h00AxxV/NzzrAcxHNeCEasA91YBbqgH3VANuqQbcUw24pRoY/YbUXBtjovQ21v9urf1l38OnjDHT+p6fBpzOzhAlsAHvsImpzAJo63tY8c+hMzkAWlQDuacacE814JZqwD3VgFuqAfdUA26pBsaG8zbXxhgD/BjYba39Xr+nngLu6Pv6DuDJkR+eWGvZxWZKKGO2Wdj/KcU/R/rnADjV7ynlIAdUA+6pBtxSDbinGnBLNeCeasAt1cDYMZQz1x8BPg98zBizre/fjcD/Bj5ujNkPrO37XkZYO800cIxWGnndrgdYovjnVv8c0Bt/1UAOqQbcUw24pRpwTzXglmrAPdWAW6qBscP0Ts/PjXJTaVeaa3K2vXz0on10i7X20gt5r+KfuUziD8rBSFANuKUacE814JZqwD3VgFuqAfdUA24NFv8hL2gmIiIiIiIiIgNTcy0iIiIiIiKSITXXIiIiIiIiIhlScy0iIiIiIiKSITXXIiIiIiIiIhlScy0iIiIiIiKSITXXIiIiIiIiIhlScy0iIiIiIiKSITXXIiIiIiIiIhlScy0iIiIiIiKSITXXIiIiIiIiIhlScy0iIiIiIiKSITXXIiIiIiIiIhlScy0iIiIiIiKSITXXIiIiIiIiIhlScy0iIiIiIiKSoYyaa2PM9caYvcaYA8aY+0ZqUDJ0yoFbir97yoFbir97yoFbir97yoFbir97ysHoccHNtTEmDPwQuAFYAtxmjFkyUgOT81MO3FL83VMO3FL83VMO3FL83VMO3FL83VMORpdMzlxfDhyw1h6y1qaAR4BbRmZYMkTKgVuKv3vKgVuKv3vKgVuKv3vKgVuKv3vKwSiSSXM9Azje7/u6vsfOYoy5yxiz2RizOU0yg83JAM6bA8U/q1QD7qkG3FINuKcacEs14J5qwC3VgHuqgVEkku0NWGvXAesAjDGdL9pH92Z7m6PUZKBpBH7O7OG8+Pfi3/iifTQ+QuMYi0YiB8OKP6gG+lENuKcacC/nOVANnEU14J5qwC3VgHuqAbeyGv9Mmut6oKbf9zP7HhvMXmvtpRlsc8wyxmzOwu8+rBxYa6dkaRxjQhZ+d9XAMKgG3FMNuOc6B6oB1YBrrnOgGlANuOY6B6qB7P7umUwLfwuoNcbMNcYUAJ8DnhqZYckQKQduKf7uKQduKf7uKQduKf7uKQduKf7uKQejyAWfubbWesaYu4HngTDwE2vtzhEbmZyXcuCW4u+ecuCW4u+ecuCW4u+ecuCW4u+ecjC6ZHTNtbX2GeCZYbxlXSbbG+Oy8rsrB8My4r+74j8sqgH3VAPuKQduKf7uKQduKf7uKQduZfV3N9babP58ERERERERkbyXyTXXIiIiIiIiIkIOm2tjzPXGmL3GmAPGmPtytV1XjDFHjDE7jDHbjDGb+x6rNMasN8bs7/t/Yg7Ho/g7jH/f9pUD1UBOKQduKf7uKQduKf7uKQduKf7u5ToHOWmujTFh4IfADcAS4DZjzJJcbNuxq621H+y33Pt9wEvW2lrgpb7vs07xdxt/UA5c52Acxx+UA9cUf/eUA7cUf/eUA7cUf/dyloNcnbm+HDhgrT1krU0BjwC35Gjbo8ktwIN9Xz8I3Jqj7Sr+vVzFH5SDM1QD7ikHbin+7ikHbin+7ikHbin+7mUtB7lqrmcAx/t9X9f3WD6zwAvGmC3GmLv6Hqu21p7s+7oBqM7RWBT/Xq7iD8qB6xyMx/iDcuCa4u+ecuCW4u+ecuCW4u9eTnOQ0a24ZFBrrLX1xpgqYL0xZk//J6211hijpdqzR/F3TzlwTzlwS/F3TzlwS/F3TzlwS/F3L6c5yNWZ63qgpt/3M/sey1vW2vq+/08Dj9M7FeOUMWYaQN//p3M0HMXfbfxBOXCdg3EXf1AOXFP83VMO3FL83VMO3FL83ct1DnLVXL8F1Bpj5hpjCoDPAU/laNs5Z4wpMcaUnfkauBZ4l97f+Y6+l90BPJmjISn+buMPyoHrHIyr+INy4Jri755y4Jbi755y4Jbi756LHORkWri11jPG3A08D4SBn1hrd+Zi245UA48bY6A3xg9Za58zxrwF/D9jzJ3AUeAzuRiM4u82/qAcoBpwQTlwS/F3TzlwS/F3TzlwS/F3L+c5MNZqmr+IiIiIiIhIJnI1LVxEREREREQkb6m5FhEREREREcmQmmsRERERERGRDKm5FhEREREREcmQmmsRERERERGRDKm5FhEREREREcmQmmsRERERERGRDKm5FhEREREREcnQ/wc2gptRYeY6cwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1008x432 with 30 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light",
      "tags": []
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "batches = [0,1,2]\n",
    "\n",
    "pylab.subplots(len(batches), 10, sharey='row', sharex='col', figsize=(14, 6))\n",
    "pylab.tight_layout(w_pad=0)\n",
    "\n",
    "# solutions\n",
    "for i, batch in enumerate(batches):\n",
    "    for t in range(9):\n",
    "        pylab.subplot(len(batches), 10, t + 1 + i * 10)\n",
    "        pylab.title('t=%d' % (t * 2))\n",
    "        pylab.imshow(states[t * 2].density.data[batch, ..., 0], origin='lower')\n",
    "\n",
    "# add targets\n",
    "testset = BatchReader(Dataset.load(staggered_app.data_path,test_range), staggered_app._channel_struct)[test_range]\n",
    "for i, batch in enumerate(batches):\n",
    "        pylab.subplot(len(batches), 10, i * 10 + 10)\n",
    "        pylab.title('target')\n",
    "        pylab.imshow(testset[1][i,...,0], origin='lower')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "UCU0AjuUMk-K"
   },
   "source": [
    "As you can see in the two right-most columns, the network does a very good job at solving these inverse problems: the fluid marker is pushed to the right spot and deformed in the right way to match the target.\n",
    "\n",
    "What looks fairly simple here is actually a tricky task for a neural network: \n",
    "it needs to guide a full 2D Navier-Stokes simulation over the course of 16 time integration steps. Hence, if the applied forces are slightly off or incoherent, the fluid can start swirling and moving chaotically. However, the network has learned to keep the motion together, and guide the  marker density to the target location.\n",
    "\n",
    "Next, we quantify the achieved error rate by comparing the mean absolute error in terms of the final density configuration relative to the initial density. With the standard training setup above, the next cell should give a relative residual error of 5-6%. Vice versa, this means that more than ca. 94% of the marker density ends up in the right spot!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "BOgtZlJrM-0N",
    "outputId": "6beb52e6-08b3-4976-96d4-59b660dd91a8"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative MAE: 0.05450168251991272\n"
     ]
    }
   ],
   "source": [
    "errors = []\n",
    "for batch in enumerate(test_range):\n",
    "  initial = np.mean( np.abs( states[0].density.data[batch, ..., 0] - testset[1][batch,...,0] )) \n",
    "  solution = np.mean( np.abs( states[16].density.data[batch, ..., 0] - testset[1][batch,...,0] )) \n",
    "  errors.append( solution/initial )\n",
    "print(\"Relative MAE: \"+format(np.mean(errors)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XKzyhAGjL-Wv"
   },
   "source": [
    "## Next steps\n",
    "\n",
    "For further experiments with this source code, you can, e.g.:\n",
    "\n",
    "- Change the `test_range` indices to look at different examples, or test the generalization of the trained controller networks by using new shapes as targets.\n",
    "- Try using a `RefinedSequence` (instead of a `StaggeredSequence`) to train with the prediction refinement scheme. This will yield a further improved control and reduced density error."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "diffphys-code-control.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}