{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f389f465-1007-443d-873b-d27bfa74cbc7",
   "metadata": {},
   "source": [
    "# Basic model builder tutorial\n",
    "\n",
    "In this tutorial you will learn the basics of contructing a model_builder interface.\n",
    "\n",
    "In our example, we construct a model_builder for a ROM of the advection diffusion equation. \n",
    "\n",
    "We will wrap around the adr_1d.py code.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "358fb709-9c81-45a3-a871-2549de121580",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "module 'mpi4py' is not installed\n"
     ]
    }
   ],
   "source": [
    "#First, let's import the relavant modules:\n",
    "import romtools\n",
    "import os\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from romtools.workflows import sampling\n",
    "from romtools.workflows import ParameterSpace\n",
    "import sys\n",
    "sys.path.append('adr_1d/')\n",
    "from adr_1d import advectionDiffusionProblem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b1ffa4be-246c-44a9-b486-f517b9f7acb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Here, we will interface around a ROM model for solving the 1D advection diffusion problem.\n",
    "The model code is given in adr_1d/adr_1d_rom.py, and it solves a ROM of the 1d advection diffusion equation:\n",
    "\n",
    "c u_x - nu * u_xx = 1\n",
    "\n",
    "\n",
    "To run the code, we require:\n",
    "    A file \"params.dat\" is required in the run directory. This .dat file contains the parameters c,nu.\n",
    "    A \"rom_data.npz\" file in an \"offline_data_directory\" containing the basis, Phi, and reduced operators for the diffusion and reaction term\n",
    "    The code is run as \"python adr_1d_rom.py --offline_data_dir = path_to_offline_data_directory\n",
    " \n",
    "Running the code outputs a solution.npz file with the keys x (grid) and u (solution)\n",
    "\n",
    "'''\n",
    "\n",
    "# We will start with defining a ROM class that meets the romtools model API\n",
    "class adrRom:\n",
    "    def __init__(self,base_dir,offline_directory,):\n",
    "        self.offline_directory_ = offline_directory\n",
    "        self.exec_dir_ = base_dir + '/adr_1d/'\n",
    "\n",
    "    def populate_run_directory(self, run_directory: str, parameter_sample: dict):\n",
    "        # Here, we need to make a params.dat file in the run_directory\n",
    "        c = parameter_sample['c']\n",
    "        nu = parameter_sample['nu']\n",
    "        np.savetxt(run_directory + '/params.dat',np.array([c,nu]))\n",
    "\n",
    "    def run_model(self, run_directory: str, parameter_sample: dict):\n",
    "        os.chdir(run_directory)\n",
    "        os.system('python ' + self.exec_dir_ + '/adr_1d_rom.py -offline_data_dir ' + self.offline_directory_)\n",
    "        return 0\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b1115b8d-0859-4a6a-ae0f-8cf8c888650c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now, we can construct the model builder, which should return an adrRom. \n",
    "#The main logic that we need to add is how to construct the basis given training directories\n",
    "\n",
    "class AdrRomModelBuilder:\n",
    "    def __init__(self,myFom,adrRomFunctor):\n",
    "        self.myFom_ = myFom\n",
    "        self.adrRomFunctor_ = adrRomFunctor\n",
    "        self.base_dir_ = os.getcwd()\n",
    "        \n",
    "    def build_from_training_dirs(self,offline_data_dir,training_dirs):\n",
    "        # The offline_data_dir is a location where we should store any data required for running our ROMs\n",
    "        # (e.g., input yamls, precomputed operators). For workflows leveraging the model builders, romtools will\n",
    "        # create these directories\n",
    "        \n",
    "        #Similar, training_dirs is a list of directories where our FOMs have been run and training data is stored.\n",
    "\n",
    "        # The advection diffusion FOM saves a solution.npz file with a key word 'u' containing the solution for each solution\n",
    "        # We will loop through these solutions to construct a snapshot matrix\n",
    "        for (i,training_dir) in enumerate(training_dirs):\n",
    "            u = np.load(training_dir + '/solution.npz')['u']\n",
    "            if i == 0:\n",
    "                snapshots = u[:,None]\n",
    "            else:\n",
    "                snapshots = np.append(snapshots,u[:,None],axis=1)\n",
    "\n",
    "        # We now have our snapshot matrix. We can now assemble our basis. Here we will just do a reduced basis\n",
    "        # Note that ROM tools requires the snapshots to be in tensor form (n_vars, n_dofs , n_snaps)\n",
    "        my_trial_space = romtools.vector_space.DictionaryVectorSpace(snapshots[None])\n",
    "\n",
    "        # We can now assemble our ROM. First, grab the basis for the state variable \n",
    "        basis = my_trial_space.get_basis()[0]\n",
    "\n",
    "        # We need to instatiate a realization of the FOM to get access to the FOM operators\n",
    "        #Let's assemble the ROM operators\n",
    "        Adr = basis.transpose() @ self.myFom_.Ad_ @ basis\n",
    "        Acr = basis.transpose() @ self.myFom_.Ac_ @ basis\n",
    "        fr = basis.transpose() @ self.myFom_.f_\n",
    "\n",
    "        #We will save this to the offline_data_dir\n",
    "        # Note we will also save the FOM operators so that we can evaluate the FOM\n",
    "        # residual we if want\n",
    "        np.savez(offline_data_dir + '/rom_data',Adr=Adr,Acr=Acr,fr=fr,basis=basis,\n",
    "                Ac=self.myFom_.Ac_,Ad=self.myFom_.Ad_,f=self.myFom_.f_)\n",
    "\n",
    "        #Now, we can instatiate the ROM:\n",
    "        myRom = self.adrRomFunctor_(self.base_dir_,offline_data_dir)\n",
    "        return myRom\n",
    "\n",
    "#That's it! This model builder can be used in, e.g., a greedy workflow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b05c0ec9-c35a-4968-a7cf-52e62426d791",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=======  Sample 0 ============\n",
      "Running\n",
      "Sample complete, run time = 0.14898300170898438\n",
      " \n",
      "=======  Sample 1 ============\n",
      "Running\n",
      "Sample complete, run time = 0.09169697761535645\n",
      " \n",
      "=======  Sample 2 ============\n",
      "Running\n",
      "Sample complete, run time = 0.08579683303833008\n",
      " \n",
      "=======  Sample 3 ============\n",
      "Running\n",
      "Sample complete, run time = 0.09924101829528809\n",
      " \n",
      "=======  Sample 4 ============\n",
      "Running\n",
      "Sample complete, run time = 0.08455801010131836\n",
      " \n"
     ]
    }
   ],
   "source": [
    "#Let's give this model builder a try.\n",
    "if __name__ == \"__main__\":\n",
    "    #First, let's create some training data following the basic_sampling tutorial\n",
    "    from ipynb.fs.full.external_model import adrExternalRomToolsModel\n",
    "    myModel = adrExternalRomToolsModel()\n",
    "    \n",
    "    from ipynb.fs.full.parameter_space import BasicParameterSpace\n",
    "    myParameterSpace = BasicParameterSpace()\n",
    "    \n",
    "    #The sampling algorithm requires a directory argument of where to put all the generated samples, files, etc.\n",
    "    work_directory = os.getcwd() + '/model_builder_tutorial/'\n",
    "\n",
    "    #Now we can run the sampling algorithm.\n",
    "    sample_directories = sampling.run_sampling(myModel,myParameterSpace,work_directory,number_of_samples=5)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "39df9105-d79a-44a6-8399-04cee6914135",
   "metadata": {},
   "outputs": [],
   "source": [
    "    #Great! Now, let's create our modelBuilder\n",
    "    #To build an intrusive ROM, we need access to the operators of the FOM\n",
    "    myIntrusiveFom = advectionDiffusionProblem(nx=33)\n",
    "    myModelBuilder = AdrRomModelBuilder(myIntrusiveFom,adrRom)\n",
    "\n",
    "    #Let's use this to build a ROM model using the data we just generated in our FOM sampling\n",
    "    #We will create a directory for our offline data (e.g., rom bases, etc.)\n",
    "    offline_data_dir =  os.getcwd() + '/model_builder_example_offline_data'\n",
    "    if os.path.isdir(offline_data_dir):\n",
    "        pass\n",
    "    else:\n",
    "        os.mkdir(offline_data_dir)\n",
    "    myRomModel = myModelBuilder.build_from_training_dirs(offline_data_dir,sample_directories)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "id": "55ac49c2-2b1a-4df7-bedb-ea897585e656",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ROM-FOM error =  4.2469193422409324e-05\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAABL4UlEQVR4nO3deZyN5eP/8dc5s8+YzTa2yVb27BESZYsilVJki0op+5qypWyRVEiRQkiy9EkyKbsUWRIhZMm+zL7PuX9/nE9+Xx+UGTNznXPm/Xw8zqOH233PvOduzP2e67ru+9gsy7IQERERMcRuOoCIiIjkbSojIiIiYpTKiIiIiBilMiIiIiJGqYyIiIiIUSojIiIiYpTKiIiIiBilMiIiIiJGeZsOcDMcDgenTp0iODgYm81mOo6IiIjcBMuyiIuLo1ixYtjtNx7/cIsycurUKSIjI03HEBERkSw4ceIEJUqUuOHfu0UZCQ4OBpxfTEhIiOE0IiIicjNiY2OJjIy8ch2/EbcoI39PzYSEhKiMiIiIuJl/W2KhBawiIiJilMqIiIiIGKUyIiIiIkapjIiIiIhRKiMiIiJilMqIiIiIGKUyIiIiIkapjIiIiIhRKiMiIiJilMqIiIiIGKUyIiIiIkapjIiIiIhRbvFGeSIicrWY+FQO/XWRI2cucOz8Bf66fIHTMRc4H3+RxLQEbi9Qhtqly9OkWgUqlyqE3f7Pb1QmYpLKiIiIi7kQncT0VRvY8McOLiZdIDr1AvEZF0jkAqleF0n3uwB+sdc/2Ab4wrY4WLAH2AMkhxGUVIEIrwqUDStPjRIVaFC+PI2rlSUkyDcXvzKR67NZlmVl5oANGzYwadIkduzYwenTp1m2bBlt27b9x2PWr19P//79+e233yhWrBiDBw+mZ8+eN/05Y2NjCQ0NJSYmhpCQkMzEFRFxeQ6Hxcqt+5m9/lu2nP2WS8HrwSf5Jg60Y08pgE9aQfwdBQmyFSTUpyA+dl9OJx/msvfvpAcdA9sNfsw7vPCJL0N4RgUiA8tzZ0Qlhj7ahvKRBbL3C5Q862av35keGUlISKBatWp069aNxx577F/3P3r0KK1ateLZZ59l/vz5bN68mRdffJFChQrd1PEiIp7o8F+XeW/VWr4+8C1HbN+Ske+E8y/yO//jlVCcklZjCgcUp2BgQSKCC1AsrCC3FSxI6YiClC1WkMiCYXjZ/3npX3R8Et/vPsSWgwfYdfJ3jsQc4KzjdxIDDoBvPGkhhzjHIc7xFTuiYe4H/tye3IFRLV+iY5MaOXsSRP4r0yMjVx1ss/3ryMiQIUNYuXIl+/fvv7KtZ8+e7N69m61bt97U59HIiIi4u5TUDD5du50FP37LjpjVxIduA7vj/++Q7kfB+EbcU7QF3Ru3oFXtSjm6zsPhsNh99BTrfj3AT0d/5/fzBziQvJ6k0N1X9sl3uT5dKrzMuM6PEhyo6RzJvBwbGcmsrVu30rx586u2tWjRgtmzZ5OWloaPj881x6SkpJCSknLlz7GxN5gbFRFxYalpGQyft5JFvy7iL7/vsAIuOe9hDHf+vV9sRSr7t6Bd9RY8/8C95A8OzLVsdruNGmWLU6NsceB+wFlQZq/ZypvfvcefgUuID9/C+2e3MGNUERoEPM/bHZ+jVrliuZZR8o4cLyNnzpwhIiLiqm0RERGkp6dz4cIFihYtes0x48aNY/To0TkdTUQkR5y7nMhLH81l2ZkppIcchjDndltKKMVTmtKkVAtebN6COuVvM5rzf9ntNp59oD7PPlCfPUcm03f+LNYnzMQRdIaNjKb2/DeIjG/HkMYv8cJD9XWHjmSbXHnOiM129Tfs3zND/7v9b8OGDSMmJubK68SJEzmeUUTkVu05fJZ7RrxGkQm3sSSxF+khh7Elh3N3+mBm1d1M0ugLnJj8BXNfftblisj/qlqmKN+PGEns6GP0Lb6IkOh7wCudE6GLeGnnPQQNqEmXqbM5H51oOqp4gBwfGSlSpAhnzpy5atu5c+fw9vamQIHrr9j28/PDz88vp6OJiGSLlVv3M2jpFA4GzAPvFAgA77gyPFq0H+/16UahsCDTEbMsyN+Xt3u0523a8/mGXYz4z/sc8F1ActguPo3pwbzxg6nl1Z2Pnx1IlVKFTccVN5XjIyP16tUjKirqqm1r1qyhdu3a110vIiLiDhwOi8lL11G470M8vKYSB4M/Au8UgqPvZnDJL0gYd5DFA15y6yLyv564tzq/T/yQw71O8pDfJLzjSmEFXGK77ySqzajKhCXfmY4obirTZSQ+Pp5du3axa9cuwHnr7q5duzh+/DjgnGLp3Lnzlf179uzJsWPH6N+/P/v372fOnDnMnj2bgQMHZs9XICKSixKT0+g1cyH5BtRm4N77OB/+NVg2isU8wgd1NhP79lYmdH0MXx8v01FzTJmi+flq6EASxv3BiNu/wi+mMo7Aswz9rTn1Xn2FxOQ00xHFzWT61t5169Zx3333XbO9S5cuzJ07l65du/Lnn3+ybt26K3+3fv16+vXrd+WhZ0OGDNFDz0TErfx1IY4XP/yIry9OJSPY+csXaf5USe/G1Pb9aFLjDrMBDboUm0TDN/qzL3AmAPku12NV989oeGcps8HEuJu9ft/Sc0Zyi8qIiJiSkWHRc8Y8Zp8YhBV4DgBbUiHuD3qZGc+8wB3FCxpO6DoGzvmCyYd6gH8MJIfS//bZTO6uh1vmZSojIiK36PP1e+ixrBdx4ZsA8Im7naduG8w7zzxNWL4Aw+lc06a9f9Jy9lPEh/0IQMWE59k4/G0KhOp85UU3e/3OlVt7RUTcybGzMVQf2pf239d0FpHUQFr5jOfy2N/4pPezKiL/4J4qpTg7bgP10oeCZWN/0AcUH12HFVt+Mx1NXJjKiIjIfzkcFj2nz6fMlPLsDngH7BmUiG3Htk6/8/UrQwjy1yPRb0agvw9bXh/H+CrfYk+MICV0L21X3UXnqR/hcLj8YLwYoDIiIgJ8uWkv4f0b88H5TjgCz+IbW47xlb/lxOQl1KkQaTqeWxrSrhm7X9hNgejm4JPEvJhnKTnwSY6djTEdTVyMyoiI5Gknz8dSa9gAHouqTmz4BkgLoIX3m1x4fQ9D2jX/9w8g/6hKqQjOvPUNrXwmQIY3J0M/5/ZJNZi9epvpaOJCVEZEJE9yOCxe+mAhJSdV4Bf/KWDPoHjMo2zt+Durhw8jOFBPgc4u3l52vn5lMLPv2YR3XCnSg4/SY8s9tHpjIukZjn//AOLxVEZEJM9ZuXUfBfo34f0zHXAEncYn9nbGVviGk1OWcndF137PGHf2TPO6/DFoJ5Gxj4NXOt+kD6HsoKdJS1chyetURkQkz4hPSqXeq8N4+JtqRIf/AGkBNLWP5fyYXxne/gHT8fKEkhFh/DlpMZ3CZkGGD8dDF1JzeG8tbM3jVEZEJE9Yu/MPIl5pwI8+48ErnaIxbdn45D6iXhtOaJC/6Xh5it1u49M+z/Jy5Cdg2dgb+D5NXn/ddCwxSGVERDxer5mf0XRJTRLDtmNLys+QUl9yasoy7qlSynS0PG3as0/xeL5pAKxjJE9OnmE4kZiiMiIiHuvMpXhuH9iN6Wc7gl8coZfvZVu33Yzv8ojpaPJfnw98iXutEQAsjutF3w8/N5xITFAZERGPtHj9Lm4bW5vDwXPBYacxozg36XvuKl/CdDT5Hz+MGEXlxJ5gs3jnxNOMXxJlOpLkMpUREfEoDodFu4nv8eR3dUkLPYA9vjjv1PyeH0aOxNfHy3Q8uQ673cYvY98jMuYJ8Epj2K5HmLvmZ9OxJBepjIiIxzh08iLFBjzC0qSXwTuViOjW/N5nN70fbmQ6mvwLXx8v9r7+KfkvNwXfBJ75viWrfvrddCzJJSojIuIRpq3YQMV3qnM2bAWk+9IucBqnJq/gjhIFTEeTmxQS5Mevr31JUPRdWAEXabOkOT/9ftJ0LMkFKiMi4tZSUjNoPGo0fX65j4x8J/GJLceipttYMuhl7Hab6XiSScUKBLO9/9f4xpYnI98JGs5qzqGTF03HkhymMiIibuun309SePD9rLeNAruD2+O6cnz4Dto3qm46mtyCCpGFWNd9Dfb44qSG7qfG5Ac5dznBdCzJQSojIuKWXp33FXfPreZ8c7vUfLwQMZ9Db31Mkfz5TEeTbFCv0m0se2wNtuRwEsK2UWnMY8QnpZqOJTlEZURE3IrDYdF0zBu8caQNVsAlAqNr8V27nUzv2dF0NMlmbe6uxKxGqyA1kIth31Ll1a56Yz0PpTIiIm4jNiGF2wd1Zq31KgBVk17m7JtbaFLjdsPJJKf0eOBuXr9zKWR4cyxkIbWG99X72HgglRERcQv7j5+nxPAmHA2ZDw4vngqewe7x08gX4Gs6muSwV598gF7FPwFgT8C7NB/7huFEkt1URkTE5a38cR9Vp9UlLnwzJIcy4c7VfNa/p+lYkovee74DjwU638dmrfUanad+aDiRZCeVERFxaW8s/paHV9YjPfgo3nFl+PqRHxncrqnpWGLAF4NepqHjNQDmXejNNz8fMJxIsovKiIi4rCcnT+fVfQ+CXywhlxuyt882WtWpYDqWGLRu5GgKXG4OPsm0X/AMqWkZpiNJNlAZERGXk5yaTrWhvVkc3wvsGZSN68KJN6IoH1nQdDQxzG638dVzH0JKMHHhW3hi8jTTkSQbqIyIiEs5cS6WyMFt2BPwLgAtvMZxcOLHhAT5GU4mrqJepdvoVHgyACviXyFqxyHDieRWqYyIiMvYsOcod4yvz4XwbyAtgEG3LWX1q0P1WHe5xtzePZxvqueTTLt5z+j5I25OZUREXMIHq7bQeEFdUkJ/w55QlHn3bWRit0dNxxIXZbfbWPnsR5Caj9jwTbSf/J7pSHILVEZExLheMz+j59b7sQLPExBdgx+7/8TTTWqZjiUurkHlkjxVYBIAX8YOZe3OPwwnkqxSGRERYxwOi0ajRjL9bEfwTqFodFv+HLGRu8qXMB1N3MSnfZ4j/PL94JPEY59013SNm1IZEREjUtMyqDL0BTbYxgBQJ20wx99aSuHwIMPJxJ14e9lZ/sxHkBpETPgGnpoy3XQkyQKVERHJdYnJaZQb2on9QR+AZaNT6Cy2jZ2At5d+JEnm3Vu1NO3zTwTgi+ghrNt9xHAiySz9yxeRXHU5Lokywx7lWMhCyPCmb4lFfNr3WdOxxM3N79uTsMuNwTeRRz7WdI27URkRkVxz6mIcZV57kLNh/4E0f0ZVWMHbPZ4wHUs8gLeXnWXdZkNqINHh63h66kzTkSQTVEZEJFccPnWJcmObEh3+A6QEM/Wu1Yzs0Mp0LPEgjauVoV3YBAAWXxrMhj1HDSeSm6UyIiI5bs+RM1R+qxEJYT9hS8rP3PvW0ufhRqZjiQda2P9FQi/fC74JPDxH0zXuQmVERHLU5t+OUfv9hqSE7sWeUJSlrdfTpdldpmOJh/L2svNl1zmQFkB0+A90fmeW6UhyE1RGRCTHfPPzARrNvYe0kD/wjitFVIeNPNKgiulY4uHur16WR0PGA7Dw4iA2/3bMcCL5NyojIpIjFq/fxYNfNCQj30l8Yyqy9blN3F+9rOlYkkcsHvASIZcbgm88bT7sgcNhmY4k/0BlRESy3QertvDk6sb/fbx7TXb1XU/tcsVNx5I8xNvLzhedZkNaAJfCv6PrtI9MR5J/oDIiItlqwhdR9NzcDPxjCLl8D78P+56KtxUyHUvyoGa17uDhfG8AMO/cALbuO244kdyIyoiIZJuhnyxj6O6HwDeRgpcf4PDob7mtcKjpWJKHfT6gN8GX64NfHK1nPavpGhelMiIi2eL56Z8y4cjj4J1K8ZjHOPLmCgqGBpqOJXmcr48XizvOgTR/Loavodu7s01HkutQGRGRW/bMu3OYdb4L2DO4Pa4rf4xfRHCgr+lYIgC0vKs8DwWNBWDemcGcPB9rOJH8L5UREbklPWfM4+OLPQC4M+kl9k+Yjb+vt+FUIldbOrAvvjEVsPwv03XGNNNx5H+ojIhIlvX5cDEfnOkKNovKiS+w681peuddcUm+Pl48V34EAN8nTeb4uRjDieT/0k8NEcmSIXO/ZNrJjmB3UC6+B7vefA+73WY6lsgNTX7mCXxjKmL5R9NtpkZHXInKiIhk2mvzvmLikfZgz6BMXGf2jvtAIyLi8nx9vOhZcSQAPyRN4djZaLOB5Ar99BCRTBm7aDVjD7YDr3Rui3mK/ePm4OOtHyXiHiZ1bYdfTCUs/2i6znzHdBz5L/0EEZGbNnHpd7z2W9srt+/+Pu5TfH28TMcSuWm+Pl68UMk5OrIu+W2NjrgIlRERuSnvrFjPkJ1twDuFiOg2HHhjIQF+umtG3M+kbu3wi6kC/jF0nvG26TiCyoiI3IQZX2+m708Pgk8SBS+35ODYzwkK8DEdSyRLvL3s9KrsHB3ZkDqVo6cvG04kKiMi8o8+XvMTL25uCb4J5L/cjIOvf0lIkJ/pWCK3ZELXR/GPrgp+sXSeOcV0nDxPZUREbmjB97/wzLrm4BdH2OXGHBi9nPBgf9OxRG6Zt5edl6s6R0c2pb3D4VOXDCfK21RGROS6lmzcTac1zcAvhuDLDdg/4iu914x4lDc7t8U/uhr4xdFFoyNGqYyIyDVWbv2N9l83xQq4RFB0XfYNX0WR/PlMxxLJVt5ednpXc46ObM54h0MnLxpOlHepjIjIVb75+QCPLG+CFXCBwOha/DZ0NSUKhZiOJZIjxnVuS0B0dfCNp8sHk03HybNURkTkinW7j/DQkvtxBJ7FP7oaewatoWREmOlYIjnGbrfRt8YoALY63uXAiQtmA+VRKiMiAsDeo2dp9mlzHEGn8IupzM5+31G2WH7TsURy3Nin2xAQXfO/oyNvmY6TJ6mMiAgnz8dSZ1pL0kMO4x1Xmh97RVHhtoKmY4nkCrvdxoCaowDYZr3H/uPnjebJi7JURqZPn07p0qXx9/enVq1abNy48R/3X7BgAdWqVSMwMJCiRYvSrVs3Ll7UQiERVxCTkMydb7YlKWwntsTCrHrqW6qXLWo6lkiuGt3xIQKja4NvAl0+mGQ6Tp6T6TKyePFi+vbty/Dhw9m5cycNGzakZcuWHD9+/Lr7b9q0ic6dO9O9e3d+++03lixZws8//0yPHj1uObyI3JrUtAwqv/Y00WE/QEow81p8Q7Nad5iOJZLr7HYbg2qPAuBn3ue3P8+ZDZTHZLqMTJkyhe7du9OjRw8qVqzI1KlTiYyMZMaMGdfd/8cff6RUqVL07t2b0qVLc8899/D888+zffv2Ww4vIlnncFjUePUl/gpdCum+TKq1nI731zQdS8SYEU+1Iij6LvBNpOuHGh3JTZkqI6mpqezYsYPmzZtftb158+Zs2bLlusfUr1+fkydPsmrVKizL4uzZs3zxxRc8+OCDN/w8KSkpxMbGXvUSkex1/+uj2Rc4Eywb/UstYOBj95uOJGKU3W5j0F2jANhue5+9R8+aDZSHZKqMXLhwgYyMDCIiIq7aHhERwZkzZ657TP369VmwYAHt27fH19eXIkWKEBYWxrvvvnvDzzNu3DhCQ0OvvCIjIzMTU0T+xZOTp7Oe0QC0D36fyd3bGU4k4hpee7IlQdF1wSeJrh9ONB0nz8jSAlabzXbVny3Lumbb3/bt20fv3r0ZMWIEO3bsYPXq1Rw9epSePXve8OMPGzaMmJiYK68TJ05kJaaIXEe/jz5ncdxLANxrjWTRgBcMJxJxHXa7jaF1RwGwwz6DPUeu/4u2ZK9MlZGCBQvi5eV1zSjIuXPnrhkt+du4ceNo0KABgwYNomrVqrRo0YLp06czZ84cTp8+fd1j/Pz8CAkJueolIrfuraVrmXr8abBZVEroyQ8jRpqOJOJyXnmiBfmi7/7v6MgE03HyhEyVEV9fX2rVqkVUVNRV26Oioqhfv/51j0lMTMRuv/rTeHl5Ac4RFRHJHQu+38GgX9qCVxrFY9qx8433sNuvP6IpkpfZ7TZeqeecxtzpNZNdh6//i7Nkn0xP0/Tv35+PPvqIOXPmsH//fvr168fx48evTLsMGzaMzp07X9m/devWfPnll8yYMYMjR46wefNmevfuTZ06dShWrFj2fSUickPf/XKITt+2BN94wi7fx29j5uPr42U6lojLGtKuGcGX64NPMl0/Gm86jsfzzuwB7du35+LFi4wZM4bTp09TpUoVVq1aRcmSJQE4ffr0Vc8c6dq1K3Fxcbz33nsMGDCAsLAw7r//fiZM0NCXSG7Ydfg0LT9rjhV8noDoGvw6fDmh+fxMxxJxaXa7jVcajGLYvubs9v6A7QcHU7tccdOxPJbNcoO5ktjYWEJDQ4mJidH6EZFMOHY2mgrjG5Ectgfv2LLsfGkzVUpff32XiFzN4bAI69+QuPDNVE16md3jp5mO5HZu9vqt96YR8VCXYpOoOr4NyWF7sCcUIarzGhURkUyw22281tC5dmSPzywOn7pkOJHnUhkR8UDJqelUHvUUsWEbISWEzx5cTeNqZUzHEnE7Ax65H//oquCdwuD5803H8VgqIyIexuGwqP1ab86EroB0P96us5L2jaqZjiXilux2Gw8VexaAr898iMPh8isb3JLKiIiHaTtxCr8FzgDLxsDSC+jbtpHpSCJubeLTHSHNn5TQvXwc9ZPpOB5JZUTEgwyZ+yVfJQ8CoLXfW0x65jHDiUTcX+mi4ZRJfhyAid99aDiNZ1IZEfEQH6/5iYl/OJ+uWjnxBZYP6Wc6kojH6NfIOVVz0HcRpy7GGU7jeVRGRDzApr1/0n1ta/BJotDlVmx/fZqeriqSjV588B58Y8uDbwKD5y00HcfjqIyIuLljZ6NpMqcVVuA5/KOrseu1Rfj7Zvp5hiLyD+x2G80K9gBg+bGPDKfxPCojIm4sITmVGuMfIzV0P/aEYmzo+R+KFQg2HUvEI03q2AUyfEgI+5nPN+w2HcejqIyIuCmHw6LGay9wOex7SA1iQauvuat8CdOxRDxWxdsKUSK+LQBjV2kha3ZSGRFxUw+8MY5D+eaAw87ISp/zZOPqpiOJeLwX6zmnan61zedCTKLhNJ5DZUTEDfWetYgox3AAngh+l1EdWxlOJJI3DHq0Kd5xpcA/hlcXLDUdx2OojIi4mRlfb+bdE10BqJHcj8UDXzQbSCQP8fay0yikOwALD2qqJruojIi4kbU7/6DXxofBO4Ui0W35ccwk05FE8pzxT3YDh53Y8I2s+ul303E8gsqIiJs4dPIiLee3wgq4SGB0bXaPmI+vj5fpWCJ5Tu1yxYmIfRCAESt0m292UBkRcQOxCSnUfusR0kIO4RVXkq0vf0Xh8CDTsUTyrO41nQtZf8n4hLjEVMNp3J/KiIiLczgsqo54htjwjZASwhePfk3VMkVMxxLJ015r3wp7QjGsgAuMXLjCdBy3pzIi4uLuGzOKYyGfQYY342sspW39yqYjieR5/r7e1PPvBsAnv2oh661SGRFxYT2nz2ODbQwAnQvMZMjjTQ0nEpG/vdnOeVfNpfAoNuw5ajiNe1MZEXFRH63+kQ9OO+el704fyid9uhtOJCL/171VS5M/2vkLwitfzDGcxr2pjIi4oJ8PnOT5H9qCdypFo9uycdQbpiOJyHV0rvwsAFuT55Ccmm44jftSGRFxMRdiEmk8sy2OwLP4x9zJL6/Nw9tL/1RFXNHopx7GllQQR9Ap3lj8jek4bks/4URciMNhUWtMdxLDdmBLKsiabispkj+f6VgicgMhQX7U9OoCwIe/aCFrVqmMiLiQFm+8yfGQRZDhzdR6S2l4ZynTkUTkX4x52Lm262zI12w/+JfhNO5JZUTERQz7ZDnfOV4FoGP+9+n98L2GE4nIzWhVpwIh0feA3cHQRXNNx3FLKiMiLmDppl8Zf/BpAO5Meon5fZ8znEhEMqP97c6FrOtjZ5Oe4TCcxv2ojIgYtv/4eZ5c3gZ8Ewi/3IRto982HUlEMunNp9tBcijpwUd568u1puO4HZUREYPik1KpN6Ud6cF/4h1blp8Hf06An7fpWCKSSQVDA7nTco5uvr9VC1kzS2VExBCHw+KukS8TE74BUoJZ+vhKyhbLbzqWiGTR8JbOhawn8y1n//HzhtO4F5UREUOenDKd34NmgWVjZOVFtLm7kulIInIL2jeqTmB0bfBKY/CCT03HcSsqIyIGTFq6liVxfQBo5TOBUR1bGU4kItnhkZLOhaxrLnyIw2EZTuM+VEZEctnanX8wZPvjYM+gTFwnvho20HQkEckmEzs9BalBpIYcYPrXm0zHcRsqIyK56OT5WFrNa4Plf5mgy3XZMWoWdrvNdCwRySbFCgRTLvVJAN5e/5HhNO5DZUQkl6SmZVDrzQ6khu7HnlCMTS8vIyyfv+lYIpLNBjd1LmQ94r+EP89Emw3jJlRGRHJJw9GvcC7sa0jz5+PmK6hetqjpSCKSA7o1q4tfTBXwSWLw/AWm47gFlRGRXPDizAX85DMRgJci59C5aW3DiUQkp9jtNh4s4lzI+tVfWsh6M1RGRHLYwnU7mXHSOWxbL30Y7z7/lOFEIpLTJj79NKT7khy2m69+3Gc6jstTGRHJQYdOXqTzfx4Fn2QKRbdk3YjXTUcSkVxQtlh+CsU3BeC9tcvNhnEDKiMiOSQ1LYO733rqyqPetw1egK+Pl+lYIpJLWpZ6BIAtl5abDeIGVEZEckijMa9yKTwKUgNZ9MgyShcNNx1JRHLRoIdbg2UjMWw72/afMB3HpamMiOSAQXOW8qP3eABeLjmbx+6503AiEcltVUpFEBLdAIC3/rPCcBrXpjIiks1W/riPtw53BaBW6gCmPfek2UAiYsy9EW0BWPvXMrNBXJzKiEg2On4uhnZL2oJvPGGX72PTyPGmI4mIQf1atgXgcuh6Dp+6ZDaMC1MZEckm6RkO6ozrTFrIIbziI9nSfzH+vt6mY4mIQfdXL4t/zJ1gz2Disq9Nx3FZKiMi2aT52Dc4G7YS0v2Y0+JLKt5WyHQkEXEBdUKcd9V89Yemam5EZUQkG4xasIofrJEAdCs8Q09YFZErejZuC8DpoNVciEk0G8ZFqYyI3KK1O/9g9N4OYLOolNiTOS93Mx1JRFxI+3ur4xVXEnySmLIiynQcl6QyInILzlyK58F5j4B/DPku12PbqHdMRxIRF2O327jTpy0An+9ZbjSLq1IZEckih8OiztgepITuxZ5QhPW9viBfgK/pWCLigrrUbQvAEe+vSE5NNxvGBamMiGTRwxOmcCJ0MWR48+69S6h5RzHTkUTERfVsdQ+2pAJYAReZuWqT6TguR2VEJAveWvo9/0keDMDjIVN58aF7DCcSEVfm7+tN2YzWAMzdprtq/pfKiEgmbfntOIN/bg92B2XjOrOo/4umI4mIG3iiqvMW373py3E4LMNpXIvKiEgmXIpNotlHj2IFXCAguiY/jZiJ3W4zHUtE3MCAh5tBaiAZ+Y6zaP0u03FcisqISCbc/fpLJIbtwJZUgKgeX5I/JMB0JBFxE/lDAiiW9AAAM9dpqub/UhkRuUndps3hUL454LAzvvYiGlQuaTqSiLiZ1re3BeCn2OVGc7galRGRm7Bo3S7mnusFQFPvMQxu19RwIhFxR4MfeRAcXqSE/cranYdNx3EZKiMi/+LY2Wg6fdUOfJIpFN2Kb14ZZjqSiLipMkXzEx7TGICpq5ebjOJSVEZE/oHDYXH3+G6khxzGK74kWwfNw9tL/2xEJOualmgLwPqzWjfyN/1UFfkHbcZP5kzYckj3ZU6LLyhbLL/pSCLi5gY+9DAAcWFb2Hv0rOE0riFLZWT69OmULl0af39/atWqxcaNG/9x/5SUFIYPH07JkiXx8/OjbNmyzJkzJ0uBRXLLuys38HXKUACeDJ+qd+IVkWxRp0IkgdG1wWYxccVK03FcQqbLyOLFi+nbty/Dhw9n586dNGzYkJYtW3L8+PEbHvPEE0+wdu1aZs+ezYEDB1i4cCEVKlS4peAiOWnPkTP03dQe7BmUiuvIgr49TUcSEQ/SIL/zAWirjy03G8RF2CzLytRj4OrWrUvNmjWZMWPGlW0VK1akbdu2jBs37pr9V69ezZNPPsmRI0fInz9rQ9yxsbGEhoYSExNDSEhIlj6GyM1KTk2nyOBmxISvwy+mEsdf+4nC4UGmY4mIB1n54z4e/rYypPtyovd5ShTyzGvbzV6/MzUykpqayo4dO2jevPlV25s3b86WLVuue8zKlSupXbs2EydOpHjx4pQrV46BAweSlJR0w8+TkpJCbGzsVS+R3HLfmBHEhK+D1CCWPbVURUREst1DdSriE1sOvFOZvGK16TjGZaqMXLhwgYyMDCIiIq7aHhERwZkzZ657zJEjR9i0aRN79+5l2bJlTJ06lS+++IJevXrd8POMGzeO0NDQK6/IyMjMxBTJslfnfcWPPs4Rvj6lZ9PyLk0nikj2s9tt1AhoC8CX+5YbzeIKsrSA1Wa7+r04LMu6ZtvfHA4HNpuNBQsWUKdOHVq1asWUKVOYO3fuDUdHhg0bRkxMzJXXiRMnshJTJFPW7T7CG/s6A1At+WWm9mhvOJGIeLJn73GuGznu/zVxiamG05iVqTJSsGBBvLy8rhkFOXfu3DWjJX8rWrQoxYsXJzQ09Mq2ihUrYlkWJ0+evO4xfn5+hISEXPUSyUnR8cm0mvs4+EeTL7oum0a8ZTqSiHi4rs3qYE8oAn6xTFv5g+k4RmWqjPj6+lKrVi2ioqKu2h4VFUX9+vWve0yDBg04deoU8fHxV7YdPHgQu91OiRIlshBZJPvVG9OHpLBfnG+A99zn5AvwNR1JRDyct5edCjbnM0cW/JK3H4CW6Wma/v3789FHHzFnzhz2799Pv379OH78OD17Om99HDZsGJ07d76yf4cOHShQoADdunVj3759bNiwgUGDBvHMM88QEKB3PBXznnv/U34PmgWWjbE1F3B3xdtMRxKRPKJDTedUzQFWkJ7hMJzGnEyXkfbt2zN16lTGjBlD9erV2bBhA6tWraJkSec7mJ4+ffqqZ47ky5ePqKgooqOjqV27Nh07dqR169ZMmzYt+74KkSxauvFXPjztLNKNbCN45YkWhhOJSF7Sp819kBKCI+gMc9ZsMx3HmEw/Z8QEPWdEcsKJc7GUnVCbtJBDFLjcnFOTVuHr42U6lojkMaUGdOBYyELqpA1m29gJpuNkqxx5zoiIp3A4LOqN705ayCG84kuwecACFRERMeKRim0B2Jm0DIfD5ccHcoTKiORJj02cxl+hX0CGNzObLKF8ZEHTkUQkjxr0SEtI9yUt5BD/+Wm/6ThGqIxInjN79U8sTxwEwGPBk+nxwN2GE4lIXlasQDCF4psC8G7UcrNhDFEZkTzl6OnLPP9de/BKo3jMY3w+4GXTkUREaFnKeVfNlkt58xZflRHJMxwOi/oTnyEj+E+840qzZehs7PbrPzlYRCQ3DWrTGiwbiWHb2bY/7z11XGVE8ox2k97lTNhyyPBhTsvPua1w6L8eIyKSG6qUjiAkxvnw0Lf+s8JwmtynMiJ5widR21mWMBCAx4LfolOT2oYTiYhc7d7CzqmatX8tNxvEAJUR8XjHzkbT49snwCuNYjGPaJ2IiLikfi3bAnA5dB2HT10yGyaXqYyIR3M4LOqP70F68FG840qxZcgcrRMREZd0f/Wy+MfcCfYMJi3/2nScXKUyIh6t/eTpnApbChk+fNhiMSUjwkxHEhG5obtC2gKw8lDeuqtGZUQ81vy1v/BFbH8A2gZNpGuzOoYTiYj8sxcaO9eNnA5azYWYRMNpco/KiHikE+di6bbqCfBOpUjMwywd1Md0JBGRf9X+3up4xUeCTxIffrvJdJxcozIiHse5TuRZ0kMO4xV/G5sHaZ2IiLgHu91GKet+AFb++oPhNLlHZUQ8TocpMzkZ+jlkePNB08WUKZrfdCQRkZt2X+n7APg1XmVExC0t/GEXi2P6AdAmcDzdW+h9Z0TEvTxzn7OMJIRs5+T5WMNpcofKiHiMk+dj6fKfJ8A7hYiYh1g2uL/pSCIimVav0m14x5YFewYfrdloOk6uUBkRj+BwWNQf9zxpIYfwio9k04C5WiciIm6rrN05OvKffXljqkZlRDxCp6kfciJ0ETi8mH7/Im4vXsB0JBGRLGtS1llG9iWqjIi4hcXrd/PZ5d4APOg/juda1jecSETk1jzbzFlGkkJ3cuTUZcNpcp7KiLi1Uxfj6LTSuU6kcHQrVgwZYDqSiMgtq162KL4xFcBm8WHUetNxcpzKiLgth8Oi/psvkBZyEK+E4mwc8Aledn1Li4hnuMPHOTryze+eP1Wjn9zitnq8/wnHQhaAw4t3Gy2iXImCpiOJiGSb5uWcZeRAisqIiEtate0AH5/pBUAzn9G88OA9hhOJiGSvZ5s1BiA59Ff2Hz9vNkwOUxkRtxMTn8Jji54E30TCo+/j66FDTUcSEcl2FW8rhH/MnQB8GLXObJgcpjIibufesYNJDtuFLakga3vNx8fby3QkEZEcUd7POVWz5qBnT9WojIhbGf7pV+wJmAbAqGpzqXF7McOJRERyTssKzjJyKE1lRMQl/Pz7X4zb1w2Amql9GfHUg4YTiYjkrOeaNwLLRmro7+z847TpODlGZUTcQkpqBs1mdsQKuEhgdE3WDR9vOpKISI4rXTScgJgaAHz0neeOjqiMiFt44M03iQlfD6n5WNFpEcGBfqYjiYjkikqBzqmatYdVRkSMeXfFJtZZowDoWWI6TWveYTaQiEgueqiSs4wcdqiMiBjxx1+X6LepA9gdlInrxIwXOpmOJCKSq3o0bwgOL9JDDrN133HTcXKEyoi4LIfDouFb3cnIdwKfuNvZ9Mr7piOJiOS6EoVCCIqtDcCcHzxzdERlRFzWk1NmcCZsOWT4MPehRRTNH2w6koiIEXfmc07V/HBUZUQk13y+fg9LYvoD0DbfBDo0rmU4kYiIOW3udJaRP20/4HBYhtNkP5URcTnnLifQacWT4J1C4ZhWLB3Y13QkERGjujdrABk+ZOQ7zro9R0zHyXYqI+JyGr7Rl9TQ/dgTirK+31zsdpvpSCIiRhUODyIkti4An6z3vKkalRFxKb1nLeZg8Edg2ZhUbz4VIguZjiQi4hKqhTqnatYfVxkRyTHrdh/l3T+fA+Ae6xX6P3K/4UQiIq6jbXVnGTnh5XnrRlRGxCUkJKXx0MdPgV8swdH1iRo+ynQkERGX8kyzepDuhyPoNKu3HzAdJ1upjIhLaDJ2BAnh27Alh7H62c/w9/U2HUlExKWE5fMnLK4+APM2etZUjcqIGDdxyVq2+UwAYFD5j6hfqaThRCIirqlmfudUzaa/VEZEss2BExcZ9nNnsFlUTHyWCV0fMx1JRMRlPVbTWUb+8llHeobDcJrsozIixjgcFo0n98ARdArf2PKsH/a26UgiIi6tc5M6kBqIFXieFVt/Mx0n26iMiDEd357FmfDlkOHDp20WUigsyHQkERGXli/AlwIJ9wDw2RbPmapRGREjVm7dz6LL/QB4OGg87RvVMJxIRMQ91C7onKrZelplRCTLYuJTaP/5U+CTRIHo5iwd1Nd0JBERt/H4Xc4ycsZ/PalpGYbTZA+VEcl1973xCslhu7ElFeT7l+fiZde3oYjIzep4Xy1ICcbyv8ySjbtNx8kWugpIrnpj0Rp2+k8BYETVj6lapqjhRCIi7sXf15vCSfcCsGibZ0zVqIxIrtn75zlG7OwMQNXkXozq8JDhRCIi7qlOYedUzU/nVEZEblpGhsX9U7vjCDyLX0xl1r0yyXQkERG39WRdZxk5F7CB5NR0w2luncqI5Ir2k6dzPvw/kO7HZ499RnhwgOlIIiJu6/GG1bAlh4FfHAt+2GE6zi1TGZEct3TTXpbGDQCgXehEHm1Q1XAiERH35uvjRZHkxgAs+dn9p2pURiRHXYpNouPSp8A7hcIxLVnc/2XTkUREPEK9os6pmu0XVEZE/lHjN4aQErYXe2IEP/SZi91uMx1JRMQjdKjvLCMXgzYRn5RqOM2tURmRHDNi3tf8GvguAK/XmkulkoUNJxIR8RwP16uMLakg+Cby6dqfTMe5JSojkiN2/XGGsXu7AVAztS+vPPGA4UQiIp7F28tO8VTn6MjSX9x7qkZlRLJdeoaDpu92wwo8j39MVb5/ZZzpSCIiHume4s4y8ssllRGRqzw2cRoX86+GNH+WtF9IaJC/6UgiIh6pU0NnGYkO3kJ0fLLhNFmnMiLZatEPu1mZOASADgWn8FDdSoYTiYh4rgdql8eeUAS8U5gTtdV0nCxTGZFscz46kS5fPQXeqRSNacO83j1NRxIR8Wh2u43IjPsBWL7LfadqslRGpk+fTunSpfH396dWrVps3Ljxpo7bvHkz3t7eVK9ePSufVlzcfW8OJjV0P/aEoqwfMFu38YqI5IJGtzmnanbH5KEysnjxYvr27cvw4cPZuXMnDRs2pGXLlhw/fvwfj4uJiaFz5840adIky2HFdY2cv4rfgt4H4M06c7mjeEHDiURE8oYujZxlJDZkG+cuJxhOkzWZLiNTpkyhe/fu9OjRg4oVKzJ16lQiIyOZMWPGPx73/PPP06FDB+rVq5flsOKa9v55jtd/fQaAmql9GNKuueFEIiJ5R+OqZfCKjwSvNGZHbTYdJ0syVUZSU1PZsWMHzZtffbFp3rw5W7ZsueFxH3/8MYcPH2bkyJE39XlSUlKIjY296iWuyeGwaDK1B1bgWfxjqvD9K+NNRxIRyVPsdhulLOe6kZW/uudUTabKyIULF8jIyCAiIuKq7REREZw5c+a6xxw6dIihQ4eyYMECvL29b+rzjBs3jtDQ0CuvyMjIzMSUXPT02x9yLvwrSPdlwWMLdBuviIgB95V2TtX8Gp8HysjfbLarFyZalnXNNoCMjAw6dOjA6NGjKVeu3E1//GHDhhETE3PldeLEiazElBz2zU8HWXipHwBtg8fp3XhFRAx55j5nGUkI2c7J8+43m5CpMlKwYEG8vLyuGQU5d+7cNaMlAHFxcWzfvp2XXnoJb29vvL29GTNmDLt378bb25vvv//+up/Hz8+PkJCQq17iWhKS0mj3WUfwTSR/dBOWDOhrOpKISJ5Vr9JteMeVAXsGc9feeNmEq8pUGfH19aVWrVpERUVdtT0qKor69etfs39ISAi//voru3btuvLq2bMn5cuXZ9euXdStW/fW0osxzd4YTWL4dmzJ4ax5cS7eXnpkjYiISSUs5w0iPxz82XCSzLu5RRz/R//+/enUqRO1a9emXr16zJo1i+PHj9Ozp/MBV8OGDeOvv/7i008/xW63U6VKlauOL1y4MP7+/tdsF/fx7opNbPVyvt/MgHKzqHVHCcOJRESkRsRd/Jm0gL2X3e8dfDNdRtq3b8/FixcZM2YMp0+fpkqVKqxatYqSJUsCcPr06X995oi4r2NnY+i38WkIdnBHQhcmdWtnOpKIiAAP3FmHZT/BBd+fcTgst3rwpM2yLMt0iH8TGxtLaGgoMTExWj9iWOkBnfkzZB7ecaU5NmwXxQro/4eIiCu4FJtEgbdCwCudLY8fo16l20xHuunrtyb65aa9/MFi/gyZBw4705vOVxEREXEh+UMCCIi7E4Avf3SvdSMqI3JTtu47wXt/OtcFNbIP59kHrl2wLCIiZpX0uQuADUfca92Iyoj8q7R0By1ndQH/aPJF1+Gboa+ZjiQiItdRt0QdAA7Ga2REPMzD4ycTE/4DpAaxvMt8Avx8TEcSEZHreKiGc2QkOnA7aekOw2lunsqI/KPPvt/FN6nDAehaZCpNqt9hOJGIiNzIQ3UrQWog+MWxevsB03FumsqI3NDFmCS6/acDeKVRLLYts3t1Nx1JRET+gb+vNyEJtQBYucN91o2ojMgNNX5zMKmh+7EnFuGHfh+61T3rIiJ51R2BzqmaH4+rjIibG/3ZN+wNfA+AsbXnUq5EQcOJRETkZjQo7VzEejTVfRaxqozINfYdO8+YXc8AUD3lZYY93sJwIhERuVmP1HGOjCQE7yI2IcVwmpujMiJXcTgsmkx9DkfQGfxiK/H9sAmmI4mISCbce2dpbEkFwCuNZVv2mI5zU1RG5Co93vuEM2HLIcOHTx6eT3hwgOlIIiKSCXa7jQIpztGRr3e7x7oRlRG5Yt3uo3x8pjcALf3H0P7eGoYTiYhIVlQKda4b+eWMe6wbURkRAFJSM3j44y7gF0dIdAOWDxpkOpKIiGRRo9udIyMnHBoZETfy8ITJxIZvhNR8/Kf7p/j6eJmOJCIiWfREA2cZSQ35nZPnYw2n+XcqI8KiH3bzbdqrAHQv9g4Nq5QxnEhERG5FldIReMWVBJvF55t2mI7zr1RG8rjLccl0/arjlaesznqhm+lIIiKSDYo4nKMjUftcf6pGZSSPa/LmcFJCf8OeWJjv+83SU1ZFRDxEtYLORay/XnT9RawqI3nYxC++Z6f/FABGVJ9N+RKFDCcSEZHs0qySc2TkjJdGRsRFHT0dzSvbugJQKek5Rj71kNlAIiKSrZ5oWAssGxn5TrDnyBnTcf6Rykgedd+kl8nIdwKfuLJ8P2Sy6TgiIpLNihUIxi+2IgBLNrv2VI3KSB7UZ9bnHAudDw47M5rNIyI8n+lIIiKSAyLtznUj6/9QGREX8vPvf/Hu0Z4A3Gt/he4t6hlOJCIiOaVWUee6kX2xrr1uRGUkD0nPcPDAzGew/C8TGFOLb4aMMB1JRERyUKtqzpGRS34/43BYhtPcmMpIHtL+relcCl8Daf580WE+gf4+piOJiEgOerRBVUj3xQq4xLo9R0zHuSGVkTziq62/82Wc8/1m2uefRMvaFQwnEhGRnJYvwJeg+OoALNvmulM1KiN5QEJSGu0XPw0+yRSMbs6Cvi+ajiQiIrmkjJ9zqmbLMdddxKoykge0ePN1ksJ3YEsO59tec/Cy63+7iEhecXekcxHroUSNjIghH3z9I5vtbwDQ746Z1Ly9uOFEIiKSm9re5RwZicv3C8mp6YbTXJ/KiAc7fTGel77vBHYHZeI7MvmZJ0xHEhGRXNa8VjlICQGfJFb++JvpONelMuLB7h8/iPSQP/CKj+T7ge+ZjiMiIgZ4e9kJT6wNwNc7XXPdiMqIhxq14Bt+zzcTgEkN5lIyIsxsIBERMaZcPue6kZ/+cs11IyojHujQyUu8vrs7ADVS+tCv7f2GE4mIiEn3lnWuGzmWppERyQWWBfdP7oUj6DS+sRVYO2yc6UgiImJYu3rOMpIU8isXYhINp7mWyoiHefmDRZwMWwQOLz5s9SnhwQGmI4mIiGG1yxXHnlAE7Bks3bzLdJxrqIx4kO0HTjH9mPOBZo3tr9K5yV2GE4mIiCuw220UTnOOjqz+1fXWjaiMeIiMDIsHpne/8iZ4q4YONx1JRERcSOVw5y+oO8+pjEgOefrtWVzMvxrS/Vjy5DwC/PQmeCIi8v81qeAcGTmF6y1iVRnxAFE7/mDR5f4APBoynlZ1KhpOJCIiruaJe5zPGkkL+YPDpy4ZTnM1lRE3l5KawaPzOoNvIuHRjVncv7fpSCIi4oLKFsuPT+ztAHy+abvhNFdTGXFzrSdMIj58K6QE8/Vzc/H20v9SERG5vmI4142s/d211o3oyuXGFq/bTVTaCACeLTGNehVLGk4kIiKurEZh57qR3y671roRlRE3FROfQpcVncArjWIxbZnZs4vpSCIi4uJaVHGWkXM+P+FwWIbT/H8qI26q2ZsjSQn7FVtSIaL6fIDdbjMdSUREXFy7e6qDwwtH0Bl2HPrLdJwrVEbc0LsrNvGz70QAhlX+kEolCxtOJCIi7qBgaCD+sXcCsGSL66wbURlxM6cuxNN/QxewWdyR0JU3Oj1sOpKIiLiRUj7ORawbj6iMSBY1mTCQ9JAjeMXfxtpBU03HERERN1OnuHPdyIE411nEqjLiRkbN/4bf830AwFv3zCWyUKjhRCIi4m4erOEcGbkcuJ30DIfhNE4qI27iwImLvL6nOwC1UvvS9+H7DCcSERF31ObuypAWAH6xrNlx0HQcQGXELVgWNJ3SC0fQaXxjKxA19E3TkURExE35+3oTHF8TgOU/u8a6EZURN/DSzEWcDFsMDi/mPDSP8OAA05FERMSN3RHoXDfy4wnXWDeiMuLidhw8zYzjLwJwn9erdLyvtuFEIiLi7uqXdJaRIykaGZF/4XBYPDC9B5b/ZQKja/H1kOGmI4mIiAd4pK5zEWtCvl3EJ6UaTqMy4tK6TpvDhfBVkO7Hoic/IcDPx3QkERHxAI2rlsGWlB+8U/ly8x7TcVRGXNWGPX8y73xfAFoHjaV13cpmA4mIiMew223kT3GOjnyz2/y6EZURF5Se4aDN7GfAN56Q6AYsHdjPdCQREfEwlUKc60a2nza/bkRlxAU9Puk9YvL/AKmBLO86Fx9vL9ORRETEwzS63TkycsKhMiL/Y9W2AyyPHwJAh4JvcV+12w0nEhERT/R4A2cZSQnZz6mLcUazqIy4kOTUdJ5Y2BV8kikQ3Yx5fXqajiQiIh6qapkieMVHgs3i8407jGZRGXEhrcdPIiH8R0gJ4Zues7HbbaYjiYiIByuS4Vw3ErXP7CLWLJWR6dOnU7p0afz9/alVqxYbN2684b5ffvklzZo1o1ChQoSEhFCvXj2+/fbbLAf2VIvX7eG79JEAPBc5jbvKRxpOJCIinq5qAWcZ2X3B7LqRTJeRxYsX07dvX4YPH87OnTtp2LAhLVu25Pjx49fdf8OGDTRr1oxVq1axY8cO7rvvPlq3bs3OnTtvObyniEtMpeuKzuCVRtGYh5nxfGfTkUREJA9oWsm5buSM3ezIiM2yLCszB9StW5eaNWsyY8aMK9sqVqxI27ZtGTdu3E19jMqVK9O+fXtGjBhxU/vHxsYSGhpKTEwMISEhmYnrFhq89hpbvMdiSyrAnp6/UaVUhOlIIiKSB5w8H0vk+2Fgs9jb5SyVSxXO1o9/s9fvTI2MpKamsmPHDpo3b37V9ubNm7Nly5ab+hgOh4O4uDjy589/w31SUlKIjY296uWpZq/+iS12Z4kbUH6mioiIiOSaEoVC8I2tAMCSzeZGRzJVRi5cuEBGRgYREVdfMCMiIjhz5sxNfYzJkyeTkJDAE088ccN9xo0bR2ho6JVXZKRnrp+4EJ3EC2s6gz2DUnFPMalbO9ORREQkj4m018E7tiyxyYnGMnhn5SCb7eq7PCzLumbb9SxcuJBRo0axYsUKChe+8VDQsGHD6N+//5U/x8bGemQhaTZ+OGmhB7AnFGXtgPdMxxERkTxo75sf4e+bpTqQbTL12QsWLIiXl9c1oyDnzp27ZrTkfy1evJju3buzZMkSmjZt+o/7+vn54efnl5lobuftL9ezy38qAKNqfkSZojeethIREckpposIZHKaxtfXl1q1ahEVFXXV9qioKOrXr3/D4xYuXEjXrl357LPPePDBB7OW1IOcPB/HoC1dwWZRIbEHrz3ZynQkERERYzJdh/r370+nTp2oXbs29erVY9asWRw/fpyePZ1PCx02bBh//fUXn376KeAsIp07d+add97h7rvvvjKqEhAQQGhoaDZ+Ke6j6YSBZAT/iXd8SdYOnWw6joiIiFGZLiPt27fn4sWLjBkzhtOnT1OlShVWrVpFyZIlATh9+vRVzxz54IMPSE9Pp1evXvTq1evK9i5dujB37txb/wrczOgFqzkQPAuAt+6ZS7ECnnersoiISGZk+jkjJnjKc0YO/3WZcu9UwRF0ihopffjlzammI4mIiOSYHHnOiNyaJpN74wg6hW9sOaKGvmk6joiIiEtQGcklg+cs51jofHDYmdHiEwqEBJqOJCIi4hJURnLB/mMXeOvA8wDUswbxTPO7DScSERFxHSojOcyyoNnUF7ECz+EXU5nVQ0abjiQiIuJSVEZyWJ9Zn/NX2BJweDGnzSeEBHn2w9xEREQyS2UkB+3+4yzvHX0RgEb24XRoXMtwIhEREdejMpJDHA6LFu89jxVwkYCY6nw9ZLjpSCIiIi5JZSSHPD99PmfDV0CGD/Me+4Qgf1/TkURERFySykgO+Pn3v/jor94ANPcdyWMNqhpOJCIi4rpURrKZw2HRauaz4B9NUPRdrBg8xHQkERERl6Yyks26vjOHC+HfQLofnz/1iUu8NbOIiIgrUxnJRpt+Pca88/0AaB00llZ1KhpOJCIi4vpURrJJRoZF69ndwS+OkOj6LB3Yz3QkERERt6Aykk2emjKT6PC1kBbAl53n4uPtZTqSiIiIW1AZyQZrdx5hSfQgANqFjadJjTsMJxIREXEfKiO3KC3dwaOfdAPfBMKiG7Ow30umI4mIiLgVlZFb9Pikd4kN3wCp+fiq+xy8vXRKRUREMkNXzlvwzU8HWZEwDICnC7/FPVVKG04kIiLiflRGsiglNYPHP+sKPkkUiG7GJy8/ZzqSiIiIW1IZyaK2EyeTEL4VUkJY9fxs7Hab6UgiIiJuSWUkC5Zv3sfqlNcA6F58KnUqRBpOJCIi4r5URjIpKSWdjku6gHcqETEPMeuFrqYjiYiIuDWVkUxqPWEiieHbsSWHsbrXB5qeERERuUUqI5nwxYZfWZs+CoCeJd+letliZgOJiIh4AJWRm5SQlEbnZV3BK40iMW1477mOpiOJiIh4BJWRm9R6wgSSwn7BlpSfNS9rekZERCS7qIzchMXrdvODYwwAvcq8y52lixhOJCIi4jlURv5FQlIaXVd0Ba80isU8wjs9njIdSURExKOojPyLVuPfJDlsF7akAnzbe4amZ0RERLKZysg/WPjDLjZYYwHoc/v7VCkVYTiRiIiI51EZuYG4xFSeWdkFvNIpHvMYU7o/YTqSiIiIR1IZuYFW48eSHLYHW1JB1vSZjs2m6RkREZGcoDJyHfPX/sIm25sA9Cs3nUolCxtOJCIi4rlURv5HbEIKPf7TBewZRMY+weRnHjcdSURExKOpjPyPluNfJyVsL7bEQqzp+57pOCIiIh5PZeT/mLvmZ7bYxwMwsOIMKkQWMpxIRETE86mM/Fd0fDLPr+4K9gxKxj7JxK6PmY4kIiKSJ6iM/NcD40eTGroPe2IEUf00PSMiIpJbVEaA2au3sc1rIgBDKs/kjhIFDCcSERHJO/J8Gbkcl8yLa7qC3UGpuI682bmt6UgiIiJ5Sp4vIy3GjyA19HfsCUX4rv8003FERETynDxdRmat2srP3pMBeKXqB5Qtlt9wIhERkbwnz5YRh8OiT1RPsDsoE9eJ159uYzqSiIhInpRny4jdbuPz9gspEtOG7wa8YzqOiIhInuVtOoBJre+uxOm7V5iOISIikqfl2ZERERERcQ0qIyIiImKUyoiIiIgYpTIiIiIiRqmMiIiIiFEqIyIiImKUyoiIiIgYpTIiIiIiRqmMiIiIiFEqIyIiImKUyoiIiIgYpTIiIiIiRqmMiIiIiFFu8a69lmUBEBsbaziJiIiI3Ky/r9t/X8dvxC3KSFxcHACRkZGGk4iIiEhmxcXFERoaesO/t1n/VldcgMPh4NSpUwQHB2Oz2bLt48bGxhIZGcmJEycICQnJto8r19K5zh06z7lD5zl36Dznjpw8z5ZlERcXR7FixbDbb7wyxC1GRux2OyVKlMixjx8SEqJv9Fyic507dJ5zh85z7tB5zh05dZ7/aUTkb1rAKiIiIkapjIiIiIhRebqM+Pn5MXLkSPz8/ExH8Xg617lD5zl36DznDp3n3OEK59ktFrCKiIiI58rTIyMiIiJinsqIiIiIGKUyIiIiIkapjIiIiIhRHl9Gpk+fTunSpfH396dWrVps3LjxH/dfv349tWrVwt/fnzJlyjBz5sxcSureMnOev/zyS5o1a0ahQoUICQmhXr16fPvtt7mY1r1l9nv6b5s3b8bb25vq1avnbEAPkdnznJKSwvDhwylZsiR+fn6ULVuWOXPm5FJa95XZ87xgwQKqVatGYGAgRYsWpVu3bly8eDGX0rqnDRs20Lp1a4oVK4bNZmP58uX/ekyuXwstD7Zo0SLLx8fH+vDDD619+/ZZffr0sYKCgqxjx45dd/8jR45YgYGBVp8+fax9+/ZZH374oeXj42N98cUXuZzcvWT2PPfp08eaMGGC9dNPP1kHDx60hg0bZvn4+Fi//PJLLid3P5k913+Ljo62ypQpYzVv3tyqVq1a7oR1Y1k5z23atLHq1q1rRUVFWUePHrW2bdtmbd68ORdTu5/MnueNGzdadrvdeuedd6wjR45YGzdutCpXrmy1bds2l5O7l1WrVlnDhw+3li5dagHWsmXL/nF/E9dCjy4jderUsXr27HnVtgoVKlhDhw697v6DBw+2KlSocNW2559/3rr77rtzLKMnyOx5vp5KlSpZo0ePzu5oHier57p9+/bWq6++ao0cOVJl5CZk9jx/8803VmhoqHXx4sXciOcxMnueJ02aZJUpU+aqbdOmTbNKlCiRYxk9zc2UERPXQo+dpklNTWXHjh00b978qu3Nmzdny5Yt1z1m69at1+zfokULtm/fTlpaWo5ldWdZOc//y+FwEBcXR/78+XMiosfI6rn++OOPOXz4MCNHjszpiB4hK+d55cqV1K5dm4kTJ1K8eHHKlSvHwIEDSUpKyo3Ibikr57l+/fqcPHmSVatWYVkWZ8+e5YsvvuDBBx/Mjch5holroVu8UV5WXLhwgYyMDCIiIq7aHhERwZkzZ657zJkzZ667f3p6OhcuXKBo0aI5ltddZeU8/6/JkyeTkJDAE088kRMRPUZWzvWhQ4cYOnQoGzduxNvbY/+5Z6usnOcjR46wadMm/P39WbZsGRcuXODFF1/k0qVLWjdyA1k5z/Xr12fBggW0b9+e5ORk0tPTadOmDe+++25uRM4zTFwLPXZk5G82m+2qP1uWdc22f9v/etvlapk9z39buHAho0aNYvHixRQuXDin4nmUmz3XGRkZdOjQgdGjR1OuXLnciucxMvM97XA4sNlsLFiwgDp16tCqVSumTJnC3LlzNTryLzJznvft20fv3r0ZMWIEO3bsYPXq1Rw9epSePXvmRtQ8JbevhR77q1LBggXx8vK6pmGfO3fumsb3tyJFilx3f29vbwoUKJBjWd1ZVs7z3xYvXkz37t1ZsmQJTZs2zcmYHiGz5zouLo7t27ezc+dOXnrpJcB50bQsC29vb9asWcP999+fK9ndSVa+p4sWLUrx4sWveqv0ihUrYlkWJ0+e5I477sjRzO4oK+d53LhxNGjQgEGDBgFQtWpVgoKCaNiwIWPHjtXodTYxcS302JERX19fatWqRVRU1FXbo6KiqF+//nWPqVev3jX7r1mzhtq1a+Pj45NjWd1ZVs4zOEdEunbtymeffab53puU2XMdEhLCr7/+yq5du668evbsSfny5dm1axd169bNrehuJSvf0w0aNODUqVPEx8df2Xbw4EHsdjslSpTI0bzuKivnOTExEbv96suWl5cX8P9/c5dbZ+RamGNLY13A37eNzZ4929q3b5/Vt29fKygoyPrzzz8ty7KsoUOHWp06dbqy/9+3M/Xr18/at2+fNXv2bN3aexMye54/++wzy9vb23r//fet06dPX3lFR0eb+hLcRmbP9f/S3TQ3J7PnOS4uzipRooTVrl0767fffrPWr19v3XHHHVaPHj1MfQluIbPn+eOPP7a8vb2t6dOnW4cPH7Y2bdpk1a5d26pTp46pL8EtxMXFWTt37rR27txpAdaUKVOsnTt3XrmF2hWuhR5dRizLst5//32rZMmSlq+vr1WzZk1r/fr1V/6uS5cuVqNGja7af926dVaNGjUsX19fq1SpUtaMGTNyObF7ysx5btSokQVc8+rSpUvuB3dDmf2e/r9URm5eZs/z/v37raZNm1oBAQFWiRIlrP79+1uJiYm5nNr9ZPY8T5s2zapUqZIVEBBgFS1a1OrYsaN18uTJXE7tXn744Yd//JnrCtdCm2VpbEtERETM8dg1IyIiIuIeVEZERETEKJURERERMUplRERERIxSGRERERGjVEZERETEKJURERERMUplRERERIxSGRERERGjVEZERETEKJURERERMUplRERERIz6f6OeV4k/sxn2AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "    # Let's test running our model and comparing it to the FOM!\n",
    "    parameter_sample = {}\n",
    "    parameter_sample['c'] = 0.5\n",
    "    parameter_sample['nu'] = 1e-1\n",
    "    run_directory = os.getcwd()\n",
    "\n",
    "    #Let's run the ROM\n",
    "    myRomModel.populate_run_directory(run_directory,parameter_sample)\n",
    "    myRomModel.run_model(run_directory,parameter_sample)\n",
    "    rom_solution = np.load(run_directory + '/solution.npz')\n",
    "    plt.plot(rom_solution['x'],rom_solution['u'],color='blue',label='ROM')\n",
    "    u_rom = rom_solution['u']\n",
    "\n",
    "    #Now let's run the FOM\n",
    "    myModel.populate_run_directory(run_directory,parameter_sample)\n",
    "    myModel.run_model(run_directory,parameter_sample)\n",
    "    fom_solution = np.load(run_directory + '/solution.npz')    \n",
    "    u_fom = fom_solution['u']\n",
    "\n",
    "    print('ROM-FOM error = ' , np.linalg.norm(u_rom - u_fom)/np.linalg.norm(u_fom))\n",
    "    plt.plot(fom_solution['x'],fom_solution['u'],color='green',label='FOM')\n",
    "    plt.show()\n",
    "    "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}