{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "yNiFcBU9f34Q"
      },
      "source": [
        "### Libraries"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "j91ZRxA3f34Y",
        "outputId": "5af71c05-b88f-418f-dfa0-091b9794434c"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Mounted at /content/drive\n"
          ]
        }
      ],
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "import os\n",
        "import torch\n",
        "import torch.nn as nn\n",
        "import glob\n",
        "import matplotlib.tri as tri\n",
        "import pandas as pd\n",
        "from PIL import Image\n",
        "from sklearn.preprocessing import MinMaxScaler\n",
        "import imageio\n",
        "from typing import Tuple, Union\n",
        "from torchsummary import summary\n",
        "from tqdm import tqdm as tqdm\n",
        "from sklearn.model_selection import train_test_split\n",
        "import torch.optim as optim\n",
        "import torch.nn.functional as F\n",
        "from torch.utils.data import DataLoader, TensorDataset\n",
        "import math\n",
        "from scipy.interpolate import interp1d\n",
        "from torch.utils.data import Dataset, DataLoader\n",
        "import copy\n",
        "from torchvision import transforms\n",
        "import torchvision\n",
        "import torch.nn.init as init\n",
        "from scipy.linalg import norm\n",
        "from skimage.transform import resize\n",
        "import shutil\n",
        "import cv2\n",
        "from scipy.linalg import svd\n",
        "import gzip\n",
        "\n",
        "import shutil\n",
        "from google.colab import drive\n",
        "drive.mount('/content/drive', force_remount=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "z6k9P75cf34a"
      },
      "source": [
        "### Load CFD Data"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "metadata": {
        "id": "c-D8zzrRg9fy"
      },
      "outputs": [],
      "source": [
        "model = 'Sten' # or Sten\n",
        "bc_name = 'Mid' # or Wall\n",
        "state_lstm = 'load'  # train | load"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {
        "id": "9QyXL92-f34c"
      },
      "outputs": [],
      "source": [
        "data = np.load(f'/content/drive/MyDrive/Paper/RoSo/{model}/{bc_name}/rylScty_patientSpec_StenosisModel_miduvw.npz')\n",
        "\n",
        "x = data['x']\n",
        "y = data['y']\n",
        "z = data['z']\n",
        "v_x = data['v_x']\n",
        "v_y = data['v_y']\n",
        "v_z = data['v_z']\n",
        "del data"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "metadata": {
        "id": "0dqGc54nf34d"
      },
      "outputs": [],
      "source": [
        "variable = v_y.T # (loc, time)\n",
        "variable_name = 'v_y'"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ob-RwiNyf34d"
      },
      "source": [
        "### General Model Info"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {
        "id": "zj4OeyScf34d"
      },
      "outputs": [],
      "source": [
        "n_snapshots = 999"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qyURirjMf34e"
      },
      "source": [
        "### Read BCs"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {
        "id": "O6H0DvJZf34e"
      },
      "outputs": [],
      "source": [
        "df = pd.read_csv('/content/drive/MyDrive/Paper/RoSo/BC.csv')\n",
        "\n",
        "#Velocity Function Interpolation\n",
        "tv_values = df['t1'].tolist()\n",
        "v_values = df['v'].tolist()\n",
        "#Create interpolation function\n",
        "interp_func = interp1d(tv_values, v_values)\n",
        "# Generate 2000 uniform values\n",
        "uniform_tv_values = np.linspace(min(tv_values), max(tv_values), n_snapshots)\n",
        "interpolated_v_values = interp_func(uniform_tv_values)\n",
        "\n",
        "#Pressure Function Interpolation\n",
        "tp_values = df['t2'].tolist()\n",
        "tp_values = [x for x in tp_values if not math.isnan(x)]\n",
        "P_values = df['P'].tolist()\n",
        "P_values = [x for x in P_values if not math.isnan(x)]\n",
        "#Create interpolation function\n",
        "interp_func = interp1d(tp_values, P_values)\n",
        "#Generate 2000 uniform values\n",
        "uniform_tp_values = np.linspace(min(tp_values), max(tp_values), n_snapshots)\n",
        "interpolated_P_values = interp_func(uniform_tp_values)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "q7hVSxV9f34f"
      },
      "source": [
        "### POD"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 54
        },
        "id": "kv-AckWaf34g",
        "outputId": "b8206207-47b6-4e8d-b6a2-08fca4764e4b"
      },
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA20AAAJFCAYAAABOeh8lAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy81sbWrAAAACXBIWXMAAA9hAAAPYQGoP6dpAACK5UlEQVR4nOzdd3hUVeLG8XfSJoUUICFAgACREpqggIAICIIIoiDqoqgoKoq9I7i7oqs/7F10RcUGskpTBJFIVaSrNCkBQgk9tPR+f3+wmc1kZpJJZpKZJN/P88zjzLnnnnsmOQ7z5t57jskwDEMAAAAAAK/k4+kOAAAAAAAcI7QBAAAAgBcjtAEAAACAFyO0AQAAAIAXI7QBAAAAgBcjtAEAAACAFyO0AQAAAIAXI7QBAAAAgBcjtAEAAACAFyO0AQBqhMmTJ8tkMjl89OvXz9NdrBauuOKKUn+OkydP9nQXAaDWIbQB8Frfffed1ZfFiIgIZWZmerpb8FL16tVTXFyc4uLiPN2Vai0mJkZxcXFq3LhxpR2jtFBY0Ufz5s0rrb8A4Gl+nu4AADgybdo0q9fnzp3Tt99+qzFjxnioR/BmDz30kB566CFJ50OBN7B3Vqp58+a6/fbbq7wvzvr8888lSStWrNDll19eKcfw9fW1PC8oKKiUYwBATUJoA+CVkpOTtXjxYpvyadOmEdpQbTz33HM2ZX379vXq0FYV8vPzLc+bN2+uAwcO2NSZMGGCBg8eXGZbixcv1ssvv+zW/gGAtyG0AfBKn376qd2/wK9evVp//fWX2rVr54FeAagqbdu2deo+xP3791d6XwDA07inDYDXKSws1Keffupwe8nLJgEAAGoyQhsAr7NkyRK7l0sV+fLLL5WTk1OFPQLgraKionTJJZeoS5cunu4KAFQaQhsAr1N0Js3RJZCnTp3SnDlzqrJLALzU0KFDtXbtWs2bN8/TXQGASkNoA+BVjh8/rgULFkiS3nnnHcXHx9utV5FLJA3D0OLFi3X//ferS5cuio6OVkBAgIKDg9WoUSP17NlTY8eO1bRp07R3716n2kxKStLzzz+vgQMHqmnTpgoJCVFAQIDq1aun9u3ba/jw4XrxxRe1atUqq8kXJOn2228vdQrzkpNV3HXXXaXWX7FihU3/ypomvWifjRs3asyYMWrRooUCAwPVsGFDXXrppfrggw+UnZ1t1WZycrKeeeYZdezYUeHh4QoJCVF8fLweffRRHT582O7Paf/+/WX2pfi9ScnJyVW65tqpU6f01Vdf6f7771efPn3UtGlT1alTR35+fgoLC1OTJk3Uv39/Pf7441qzZk2pbRVfL86elStXlvm7Lmn9+vV64okn1KNHDzVq1Ehms1lhYWGKi4vTjTfeqOnTp5fr7HN2dramTp2qwYMHKyYmRmazWdHR0erfv79mzJjhdDuesGLFilJ/vkWcHfubN2/Wvffeq9atWys0NFQhISFq166dnnjiCR0/frxcfXP191TWZ0LRuDcMQ999952GDBmiZs2ayc/Pz1Lns88+s2k3MTFRTz75pLp06aL69esrICBADRs21KBBgzRt2jQVFBQ49Xnk7PILJXXo0MGpzzgApTAAwItMmTLFkGS0bNnSKCwsNF577TVDkt3H7t27nW533bp1RocOHey2ExkZaTRp0sTw9fW1Ku/YsaNx8uRJu+2lpqYat99+u80+kozg4GCjWbNmRkhIiFV5/fr1jRkzZljaeOqpp4y4uDgjNjbWbr/GjBljdcyJEycacXFxRt26de3WX758uU0/4+Liytzn7bffNvz8/Bz+nC+99FLjzJkzhmEYxqJFixy2JcmoW7eusW7dOpt+JCcnW/oSGhpqd9+kpCRL/aNHj1rq26vbt2/fUn/f5dlnz549Nr9HX19fIyYmxqhfv77dtvr06WMkJibabe/tt98ute+BgYGW7UWPp556ym5bBw8eNK644gqbNho1amQEBgZalTVr1sxYvHhxqT8XwzCM3377zWjatKnD36EkY9SoUcbSpUvtbnv22WfLPEZ5OBr/06dPt1t/+fLlljqlcWbsT5061fD393f4c4iKijL++OOPMt+Du35PRZ8JpY37vLw847bbbnPY55I/t8mTJ9t9j/Xq1TOio6MNSUbv3r2N66+/3m57sbGxljEaFxdnNGvWzG49Hx8fS99Luuqqq4y4uDgjPDzc8jlR2rgHYIvQBsBrFBYWWr6svPDCC4ZhGMbx48cdfqly9h/8r7/+2m4bo0aNsgoKZ86cMSZNmmSYTCa7QaLI0aNHjdatW9u016JFC2Px4sVGfn6+YRiGkZ+fb/zwww9G48aNS/3Cm5SUZPf9lQxtRZ599lmHX0IdcbTPQw89ZEgyzGaz0axZM8PHx8duvbvuusv45ZdfjICAAEOS0bRpUyMoKMhu3ebNmxtZWVkO+zJmzBi7+9n7WRtG+QJYRfbZsWOH1RfZmTNnGnl5eZbtu3btMoYPH27TXuPGjY39+/e7rR8lbd682YiMjLTaNzg42FixYoVhGIaRlpZm0y8/Pz+rPwyUtHr1aoe/t6CgIKN58+aW3/E111xjt151CW1FHI39Rx55xDCZTEZAQIARGxvr8HOmdevWRk5OjsP2K+P3ZBiOx87EiRMtr2NiYoyIiAiHP7eHH37Ypo2QkBBjwYIFljpr1641GjdubPW5V9r/l7m5uVafacUfW7ZsKfU9XXTRRYYkY+HChaXWA2CL0AbAaxT9Zd/X19c4fPiwpfy6666z+wWhQYMGRm5ubqltbtiwwTCbzTb7Dhw40CgsLLS7z6OPPurwC0teXp7RvXt3m/bMZrOxZ88eu+39/vvvljM53hbaio6TmppqGMb5kNypUyebOj4+Pka9evWMCy64wNi2bZthGOdD6R133GG3za+++sphX7w5tDk6A5KXl2d06dLFps2RI0e6rR/FnT592m6YefLJJ63qHT9+3OYsYXBwsLFjxw6bNjMyMowWLVrY7dODDz5oZGRkGIZhGGfPnjX+9re/ORwvNSW0STJuvvlmy1nkjIwM49prr7Vb7z//+Y/dtivj91TEXj9iY2ONgIAA49JLLzV27dplqfuf//zH8geXop/bokWL7Lbxr3/9y+ZY8+bNc/gzsvf/5fPPP2+37t133+3w/WzatMmQzv/Bp6CgwGE9APZxTxsAr1F0n9qQIUPUuHFjS/mdd95pt/6JEyf03XffldrmY489ZvcekgkTJji8J2bSpEny9fW1u+3zzz/X+vXrbcpvuOEGxcXF2d2nS5cuGjp0aKn99JT4+Hh98sknCg0NlSQ1aNBADz30kE29wsJCnT59WjNmzFD79u0lSb6+vvrnP/9pt92EhITK67Sb+fv7KyYmRpdddpmuvPJKu3X8/Px0xx132JTPnTtXJ06ccHufXn31VbszqA4ZMsTqdYMGDXTxxRdblWVmZurvf/+7zb5ffPGFkpKSbMrbtGmjt956S8HBwZKk8PBwff7552rUqJErb8Fld9xxh937oC6//HK3tB8bG6vPPvtMERERkqTg4GC9+eabdusuXrzYbnll/J5Kc+DAATVp0kQ//vijWrdubSm/8cYbNWzYMKu69hZ2l6TrrrvOpmzYsGEKCwtzuh/jxo1TQECATfmMGTN05swZu/sUfb7fcccd8vHh6ydQXvxfA8ArnDp1yjL721133WW17corr1RMTIzd/UqbkOSPP/7QL7/8YlNuNpvVu3dvh/tFRkbqwgsvtLvtnXfesVt+xRVXOGzPme2ecsstt9gE1LZt29qt27p1a3Xv3t2qrHnz5na/vO3YscN9naxkcXFxSk5O1qpVq0qtV/xLchHDMLRy5Uq39qegoEAffPCB3W32JuZp06aNTdl3332no0ePWpV98cUXdtu89tprbb5Em81mjRw50tkuV4oGDRooLi7O5lH8DzquuPHGG+Xv729V1qJFC8sfMIr766+/bMoq6/dUln/84x92+3jLLbfozjvvVOvWrbV7926tW7fOpo6fn5/dfvj6+qpjx45O9yE6Olo33nijTXlmZqY++eQTu+UzZ86Uj4+Pxo4d6/RxAPwPoQ2AV/jiiy+Uk5OjRo0a2ZyV8vX1dTjLWEJCgtWsgyW32RMbGyuz2Vxqf5YsWaJDhw6padOmlrLjx49ry5Ytduvb+0Jf3F133aVDhw7pscceK7VeVevQoYNNWb169ezWLTrDVlLdunVtys6ePetSv7yRozHjaMbMitqwYYPdn5/JZFJUVJRNub2y/Px8qzCZnZ2tTZs22T2eo99r586dnetwJXn55Ze1Z88em4e7Zrfs1KmT3fLw8HCbMntnjyrj91QWk8mk4cOH2912/fXX6+OPP1avXr3066+/2q3ToEEDh1cRlPfM6oMPPmi3fOrUqSosLLQq++abb5SamqpBgwYpNja2XMcBcB6hDYBX+PjjjyWdn/La3peKsWPH2r2c0TAMy74lbd261W55w4YNy+xP/fr11aRJE6u+bNu2zWH9stoMCgpSkyZNynUJUlWw9xd7R+HEXl1Jds+0VcfFzw8dOqTXXntNV199tS644AKFh4dbTaXu6LK8c+fOubUf27dvt1tuNpvtXlYWFBRkt37xy3gPHDig3Nxcu/UcfVl35v+T6sxeiJJk9/PH3niujN9TWZo3b265nLM0iYmJdstDQkIc7lN0eayzunfvbnPmXTq/DMoPP/xgVVZ0RcTdd99drmMA+B8/T3cAAFavXq2//vpLJpPJ4f1rLVu2VL9+/bR8+XKbbdOnT9dzzz1n82Xr1KlTdtty9OWpLI7ac6VNT7MXhB3db+LoHsDqfn+KYRj65z//qVdeecUm2Pj4+Khp06YKCAhQVlaWjhw5Ynd/d3I0znJycnTBBRfYlJ8+fdpu/eKX3Tm6z0hy/EW+vF/iq5uyzraXpTJ+T2Vp0KCBU/Uc/b5LXg5aXFnr3tnz4IMP6tZbb7Upf/fdd3XNNddIOn9p6W+//abo6Gib++4AOI/QBsDjPvroI0nnv/za+7JTliNHjmjhwoWWLwlAeUyYMEGvvvqqVZnJZNJLL72khx56SIGBgZLOL+rsrkkwSuMoBBqG4fSi75Lzl6g6+iJfkS/xVaFfv35uD8oVUdW/J0mWsVhReXl5DrdV5Gd644032l2E/Oeff9aOHTsUHx9vOct2++23lxoaAZSuev95FEC1d/bsWX377bcut1MU/IqrX7++3brZ2dkVOoaj9lxpszSOvkRVxrG8gb33W9nv9fDhw3r99ddtyseOHaunnnrK5S/JFeFonAUGBso4v1SPU48ff/zRsq+9+w6LOPoi7w3ByJtVxu/JXRz9vjMyMhzuk5WVVe7jBAQEaNy4cXa3vffee8rJydGXX34pk8lkM8EUgPIhtAHwqBkzZigrK0tRUVHKyckp8wvOo48+aredxYsXKzk52arM3iQbUvkuR3KmPVfalOzfEybJ4T1IJ0+erPCxvEF53m9lv9clS5bYTJogSYMGDarU45bG0eyd2dnZys/Pr1Cbjmb5lKT09HS75ZmZmRU6Vm1RGb8nd3E0MdKJEydUUFBgd1tFP8PGjx9v9wzaF198oenTp+vUqVPq169fha6iAPA/hDYAHlV06cyYMWMcfqksztE9bwUFBfr000+tyhx98d6/f3+pE2WcPHlSnTt3VufOnTVlyhRLeXR0tMMZ53bt2lVqv8ePH6/OnTvr0ksvtQkJjiYWcHTPjKMZLKsLR2cB7L3fyn6vjr6o2rufy90TjjjSvXt3h/eT2Vtnrbj169friSee0BNPPGE1cY7ZbFbXrl3t7uPoZ3Ds2DEne+xZx44dU79+/dSvX78q7XNl/J7cxdGSJvn5+dq5c6dNeWFhYYX70ahRI7vLQ6Snp1v+yMYEJIDrCG0APGbDhg3avHmzJDm8xKak9u3bq0ePHna3ffLJJ1aBqEuXLurTp49NvdzcXLvrtxVZunSpNm/erM2bN9tMh/7II4/Y3efnn3922F5eXp5mzJihzZs3Kzo62mbijuDgYDVv3txmP3trne3atUu///67w2NVB+3atbNbbm8trJkzZ1ZqXxx96bY386i9SXAq0n7JszDvv/++evToob/97W+Szp+JdLTERVnTw7/zzjt6/fXX9c4779jM/jhmzBi7+ziaZfXPP/8s9VjeIjs7WytXrtTKlSur9NLhyvo9uUOrVq3Us2dPu9vmzJljU7Zo0SKXlulwNP1/dna26tevb3dBbwDlQ2gD4DFFZ9n69eunVq1aOb2fo7NtBw8e1E8//WRV9vrrr9u9L+nVV1+1e89OTk6OXnjhBUnnLzG6+uqrrbbfdtttdkPj7NmztW/fPrv9evnll5WWliZJDtdpGzx4sE1ZcnKy3n77bUs/9+zZo1GjRtm9nK86GThwoN1p1d966y3LJa75+fl6++239fXXX1dqXxwtov72229bQothGPr888/1/vvvl7t9e2G85BnFVatWad26dVYzU06aNMnuennvvvuuw3vQ1q5dq9mzZ0s6vy5gZGSk1fZbb71VcXFxNvvNnz/f5pK5nJwcu1/uYa0yfk/uMnnyZLvlr7zyitV9dL///rvuu+8+l47Vq1cvXXTRRXa33XrrrS7P1AmA0AbAQ9LT0y1fyJ09y1Zk1KhRqlOnjt1tRUGwSNeuXfXFF1/YXHq5ZMkS3X777VZflPfs2aORI0dq+/btCggI0GeffWZzVszX11fz58+3uZ8lJydHV155pdVf2DMzM/X666/r+eefl3T+r9GOLlt6+OGH7X6xeeSRRxQWFqbY2Fi1bt1aBw8edHjGpLpo3LixRo8ebVO+fft2NWvWTM2aNVO9evX0yCOPOLyH0V369Olj916b48ePq0uXLmrcuLFCQ0N1++232z1rW5aBAwfalO3Zs8eyxtemTZv0/fffS5L69+9vqRMTE6NvvvnGZimJLVu26Nprr7W6HDcjI0Mffvihhg0bppycHMXHx+ull16yOW5QUJBmzJhhc/Zv7969euCBByz3tp05c0a33XabS/dpusNff/2ln3/+uczH6tWrPdbHyvg9ucugQYP0+OOP25RnZGRoyJAhio6OVtOmTXXxxRfr7NmzLs+M6uhsG5dGAm5iAEAVWrt2rREXF2fExMQYkgxJRmxsrBEXF2fExcWVum9ycrKlntlstuxf/OHr62ups3btWsu+GzduNDp37mxT32QyGY0aNTKio6MtZQ0bNjQSEhJK7UtaWppx9913G35+fjZthoWFGc2aNTP8/f0NSYaPj4/xxBNPGAUFBaW2+fnnnxu+vr5235cko0GDBsYvv/xiPPvss3a3N27c2IiLizNuvvlmS5tFP4u6deuWus/atWstv5vY2Fi7dUNDQ424uDjjqaeesmrb3s/Az8/Psn3OnDk27/XMmTNGt27dHL5XScZtt91m5OXl2d0WGBho83t+++23LWXO7mMYhrF+/XojIiLCYT/8/f2NF154wVi+fLnd7XXr1nU4fg8fPmw0btzY7jht2LCh5XWLFi2MlJQUm/23bNlidOnSxe5xIyMjrdqQZFx99dXGyZMnSx1na9asMZo1a2bTntlsNpo2bWr4+/sbZrPZePzxx8v9fp1V2u/dlUdSUpLlGM6O/aLx6cp4dtfvqTxjuG/fvk79rCdPnmz5LHLU5ty5c40xY8aU+TMtTVZWlhEZGWm1b69evZzaF0DZCG0AqpSjL75Fj9IkJSWV6wvc8uXLbdpYtmyZ8fDDDxvdunUzGjZsaJjNZsNsNhuNGzc2Bg8ebLz33ntGamqq0+/n0KFDxpQpU4whQ4YYsbGxRmhoqOHr62vUrVvX6Nq1q/HYY48Z27Ztc7q9P//807jllluMpk2bGgEBAUbdunWNTp06GX//+9+No0ePGoZhOAxtRY/iX+bK87Mq63dT9BgzZky52p4+fbrd95qdnW28+eabRo8ePYywsDBLaLj22muN7777zun3UPR7LuvnUtrYOHDggPHggw8arVu3NgIDA42QkBCjTZs2xv3332/5/Tnz87Hn5MmTxj//+U+jV69eRlRUlBEQEGAZI927dzf++c9/GqdPny51XCxZssS4//77jS5duhhRUVGWYNWgQQPj0ksvNR5//HFj3bp1pbZRXFZWlvH+++8bgwYNMho2bGj4+/sbYWFhRrt27Yz777/f2Lt3b4XfrzPK8/9xeR7FA0Z5x6er49kwXP89lWcMx8bGOv3zTkxMNJ544gnjwgsvNCIiIgyz2Wy0aNHCuOOOO4wtW7YYhmEYN998s93jFH3uOGPixIlO/6wAlI/JMFiIBQAAoDa75pprtGDBAqsyHx8f5ebm2r0H1Z5Dhw6pRYsWKigoUFhYmI4ePepwsh8A5cM9bQAAADXQsmXLSp0ptzh7Eyl16dLF6cAmSU2bNrVMqnTzzTcT2AA38vN0BwAAAOB+gwYNUmFhoQ4fPqxGjRo5rHf8+HG7S4zcdNNNdut/9dVXCg4OtpnKPy8vT3/88YckJiAB3I0zbQAAADWUYRh6+eWXS63zyiuv2Cwl0qZNGz3wwAN2699666266aabbNYcnD59uo4cOaKuXbs6XAIAQMUQ2gAAAGqwt99+W/fcc48SExOtypOSkvTggw/qjTfesCqPjY3Vjz/+WOr6arm5uXr++ect60guW7ZMTz31lCTpmWeecfM7AMBEJAAAADVQUFCQsrOzrcoCAwPVoEEDpaWl6cyZM1bbfH19dfvtt+vll19W/fr1HbZrMpksz0NDQ2U2m5WSkiJJuu222/T555+78V0AkAhtAAAANVJaWpp++eUXrVmzRn/++aeSkpJ05MgRpaenyzAMhYWFKSoqSp06dVKvXr00atQoNWzYsMx233vvPa1YsUJ//vmnjh49qoCAALVq1Up333237rzzTvn4cCEX4G6Etip08cUX6+jRo4qOjtbKlSs93R0AAAAAHmIYhtLS0tS4ceMy/9hBaKtCjRs31tGjRz3dDQAAAABe4tChQ2rSpEmpdZjyvwoVXQPeqFEj7dy506N9ycvL05IlSzRo0CD5+/t7tC+ofhg/cAXjB65g/KCiGDtwRWWMn9TUVDVt2lShoaFl1iW0VaGi0Obj46OwsDCP9iUvL0/BwcEKCwvjgwvlxviBKxg/cAXjBxXF2IErKnP8FJ/cx5Fqf6doYWGh3n//fYWFhclkMmn//v1ua/vIkSN6+OGHFRcXp8DAQEVHR+vqq6/WTz/95LZjAAAAAEBpqnVo2759u3r37q0HHnhAaWlpbm177dq16tChg6ZNm6Z7771Xq1at0tSpU3Xo0CENHjxYkyZNcuvxAAAAAMCeahvann32WV100UXy9fXV008/7da2T548qWHDhunMmTOaOXOmnnzySXXv3l0jR47UqlWr1LRpU02ZMoV1SAAAAABUumob2t566y29+eabWrVqldq0aePWtp9//nmlpKTokksu0fDhw622hYeHa+LEiZKkCRMmKCsry63HBgAAAIDiqm1o++uvv3Tfffc5deNeeeTm5urLL7+UJI0cOdJunaLy48eP64cffnDr8QEAAACguGob2mJiYiql3dWrV+vcuXOSpG7dutmt06BBAzVr1kyStHDhwkrpBwAAAABI1Ti0VZYtW7ZYnjdv3txhvaJtxesDAAAAgLuxTlsJBw8etDyPiopyWK9o26FDhxzWycnJUU5OjuW1YRiW53l5ea5002VFx/d0P1A9MX7gCsYPXMH4QUUxduCKyhg/5WmL0FZC8aUDAgMDHdYr2paamuqwzpQpU/Tcc8/ZlGdnZ2vRokUu9NJ9EhISPN0FVGOMH7iC8QNXMH5QUYwduMKd4yczM9PpuoS2SjRx4kQ99thjltfx8fE6cuSIAgMDNWTIEA/27HyyT0hI0MCBA92+qjtqPsYPXMH4gSsYP6goxg5cURnjp7STPyUR2koIDQ21PM/OzlZISIjdetnZ2ZKksLAwh22ZzWaZzWbL6+IzXXrLh4W/v7/X9AXVD+MHrmD8wBWMH1QUYweucOf4KU87TERSQtGskNL5RbYdKdrWtGnTSu8TAAAAgNqL0FZCp06dLM/379/vsF7RtuL1AQAAAMDdCG0l9OrVS+Hh4ZKkjRs32q1z4sQJyyyTQ4cOrbK+AQAAALXdrIRduubx7zQrYVeV7utJ3NNWgtls1q233qr33ntPc+bM0RNPPGFTZ+7cuZKk6OhoXX311eU/iGGoIOOcq111SUF+vnzzs1WQmSofP4YByofxA1cwfuAKxg8qirFTeeat3KM5y/ZoZP8LNKLvBZW637yVezR/2R4Fm6T5P/0pv7z0Ktm3MsZPQYbzE5GYjOKLh1VTn332me644w5JUlJSUqmLYkvSRx99pCeffFIdOnTQwoULFRERYbX95MmTateunVJSUvTdd9/pmmuusWxLTU1Vp06ddODAAX322WcaM2aM0/1s0qSJDh8+rIZhgVr9lGdnjwQAAADgOWnZeer8wgKdO3eu1MkNpWp8eeSJEye0bds2bdu2TYcPH7aU796921KekZFhd993331Xqamp+u2337Rs2TKb7VFRUVqwYIHq1q2rm266Sa+99po2bNigefPmqU+fPjpw4IAmTpxYrsAGAAAAABVRbUPb1KlT1bFjR3Xs2FF///vfLeVXXnmlpXzDhg12933ggQcUFhamnj17qn///nbr9OjRQ9u2bdOdd96pDz74QJdddpnuueceNWnSRIsXL9b//d//Vcr7AgAAAIDiasTlkdUFl0cCAAAAkMp3eSR3YXqAb516in3kU4/2IS8/Xz///LOuuOIK+XMzLsqJ8QNXMH7gCsYPKjp5xexleyyvr3dy35L7Fd//2j5xysktUFZuvrJz8pWdW6CsnPPPf/3ziNbvOG6zX9PoOoqMCFJOTtF+BcrJ+9/+hZxKcasAPx8ZkvLyC8usW9aYqIzPntTUVOmFBU7V5dPOE0wm+YaEe7QLhXl5KvALlG9wmHzdtKo7ag/GD1zB+IErGD81x6yEXZq5eKduHtxWowa2cXqfGUuTJQXq86XJyvevU+a+xfcp8vnSZJ3ODdDlFzdRZna+MrPzlZWTr6zsPGXmnH+9OfGkEg+dtdqv+P6fL00uo7e2++08li8dSytRapLk/99HzeNjkgL8fRXg76v8gkJlZueXuU9883rq1CpS5v/uF+Dnoz92ndSabUcd7jO8b5xGXt5KAf4+8vfzkZ+vj0wmk6T/joHFOx3uO3pwW11fxjiqjM8e3wKT03UJbQAAAKiQigSvov2KvkQX/dep8FXii/eMxTt18FiaLmoTpfSsfGVm5ykjK08Z//1v0pFUHT+dabe9Bb/s04Jf9jnd55rG19ekgoKyT+21ia2rDi3ryxzgJ7O/jyVIbdp5Qqu3HHG43w0DWummQW3l52uyhCfJuQBlbyxc1auFw30d7VOkaFtF9vUWhDYAAIBarsJnvcoZvAoLDX314w59uyzRqnzG4p3akXRKrZrWVUZWntKLwtd/HyfPZCkzx/4Zml/+PKxf/jxsd1tN1LB+sFo0Dpc5wFeBAX4KDPDVnuSz2rb3lMN9ru7dQiP6XiBzgO/5h7+vTCZThQOUJA28JLZCIcqVAGVvX2dDlyv7egNCGwAAQA1RFeGr4L/Ba7ad4LVtb4paxkScD15ZuUrPzFN61vlHRmauMkq5NO73XSf1+66TTvXZG/n5+sjHx6TcvIIy63Zt20CXXthYQWZ/BZrPh6/gQD/9vP6gvi/l7F9pIaOqA5Sj/Stzv+L7VuQMryv7ehqhDQAAwItU1SWHhYWGvnQQvv7cfVLNGoYqIzNPaZm5SsvKU3pmrtIyz5/5cmRzYoo2J6Y43WdvYPb3Vd0ws4LN/goKPB+egsx+OpKSoT2Hzjrcb+ilLTSi3wUKMp+v7+93fiUtV85e3T28o0JDAqr0LJSrZ6AqGoRcDV8VDVyu7OtJhDYAAAAvUZFLDkvuV2TG4p3afeC0LmhaV+fSc3QuI1ep6bk6l5Fj+a+jhZ+27zul7fscX27njaLrBatlTLhCAv0VEuSvkEA/7Tp4Rpt2nnC4T1nhpDqdvSq+b1UGqKL9KxKEqmuA8gRCGwAAQCUo7xkzR8ErIytPl3WO0bn0HJ1Ny9HZEv89cCxV59Jz7ba5YccJbdjhOLR4owuahKvTBVHng9d/H5t2HteKTY5na3T3ZYNFqtvZq6J9CVA1D6ENAADAgcq4VDG/oFDn0nN0JjVHZ9KydSYtR6v+SHZ4WeH8lXs1f+VeF99J1WrdLEIXtYlWnWB/1flv8KoT5K9fNx/RwtVJDvdzFGz6XdREMVF1qnzyiuL7V5ezV6iZCG0AAAB2lOdSRcMwlJmdr9Op2ZqzPFFLNxyy2j5j8U798Os+GYaUmmH/rJg3at+ynrq3a6Q6wf4KDfZXneAAra5g8JKkDnGRigg1e2TmwIKCAs1KSNSoga24/A/VDqENAADUeO66VHFf8jm1ahahM2k5On0uW6dTs3UmLVunU3PKnDXQ0SWMVeGCJuHq0qaBwkLMCq8ToPAQs9ZuO6of1+x3uI+jQNTRheAleW7mwBsHtFKdnEQNGdDK6X0Ab0FoAwAA1cY3SxM1a5WUbk7U6MHtnNrH0RmzrJx8nTqXpdOp2Tp9Llun/hvC/kw8qYPH0uy2tWbbUa3ZdtQ9b6acgsx+ahQZoohQsyLq/PcRatbO/af121bHfXIUiC5q20D1wgOr1SWHRfty9gq1DaENAABUC7MSdmlWQuJ/nyfK19fX7pf3/ILC80EsNVvfr9pns/DyjMU79Z+E3covKKySfjvi42NSRB2z6oaZlZWdryMpGQ7r3jyojW66sq3D7RWdbMNTZ72K9id8Ac4htAEAgCrl6gLQRWYs3qmte1LUOKqO5YzZqXPZOpfueCr7Ip4MbEMvbaGbBrVRaHCAfHxMlnJPzHJYfF/OegHei9AGAACqTGmTexQUGjp9LlsnzmTq5NksnTyTqRNnsrQ58aSOOjgLtWVPirbsqfrFnFvGhOnittGqGxqoemHnH3XDzPp5w0H9J2G3w/3Ku74XlxwCkAhtAACggsp7xmzmTzv19ZJdVmUzFu9UwroDksmkU2ezVFBYximyShAS6Kd64YGqHxakeuGBOnYqQ38lnXZYv7QgdcvgePn5+rh8xoxLDgEUR2gDAKAWc+c6ZDcOaK0zadk6dipTx08XPTJ0/HSm9h0+p8zsfLttnTiT5fobqaC/XdFat1wVb1PuqUsVi/YneAEojtAGAEAtVZ51yIqkZ+Xpi4V/2UwVP2Px+bNohR44UyZJHePqq+9FTVU/PFD1w89frrjotyTN/GmXw3289VJFACiJ0AYAQA3grnXICgoMXd61yf/Olp3K+O/z8/9Nz8pz2Ka7AluAn4+i6gYpqm6woiKCdPx0Zqn3rTkKUzcNaiuTycSligCqPUIbAADVXHnPmH2+8C/NXpbosK1ZCY7PTlWmwT1idfPgtoqoY5bJZLLa5okp7Yv2J3gB8DRCGwAAXsRdZ8wys/LUrX1DHUvJ0NH/ni07eipDB46mKi+/8qe7j4mqo+j6wYquF6yG9YK1J/mczXppxXnzemIA4GmENgAAvER5zpgVFBTq0wXb9f0v++xun7dyr+at3Fs5HS2Do4WgmzUMrfClitL5n0VBQYFmJSRq1MBWTGkPoNYgtAEA4AUcnTFLy8hVlzYNdDQlQ0dS0nU0JeP841RGmQtIu8Lfz0fNGoaqYb2Q82fL6gcrun6INu047jAoSpU7uYck3TiglerkJGrIgFZO7wMA1R2hDQAANyvPJY6FhYY+XbBd362yf1bs+1/2lRqSKsMN/VvptqHt7G67qE0DhYYEeGxyDwCojQhtAAC4kaNLHDOy8nT4ZLqST6Tr8Ml0Hf7vfw8dT6vUBaXD6wSobWw9NYoMUcP6IWpUP0QNI4O1YlOyzULXEuuQAYA3IrQBAOBAeScF+XrJTpt1wWYs3qk5yxKVnVtQWd106OYr2+imQbb3lp3f1lY+PibWIQOAaoDQBgCAHWVNCpKbV6BDx9OUdCRVSUfPad22Yzp+OtNuW+4IbPXDA9UxLlKNI0PU6L+PxlF1tHB1EuuQAUANR2gDAKAER5OCbN93SmEhAdp/NFXJJ9Ldtph0WSpzcg+CFwB4P0IbAKDGc/Yyx4KCQk37bpsWrk6yu/3P3Sdd7kujyBDFRNVRkwZ1FBNVRzEN6qhJVB0tXrvf5tJKick9AACENgBADefoMsfsnHztP5aqfYfPWT0qc1KQUQPbaPRg+/eY3TSorUwm1+4xI6wBQM1EaAMA1FiOLnNc8Ms+pWfmyh35rGH9YLVoHK70zFxt3XvKYT3OmAEAKorQBgCoFpy9xDE7N197k8/pm5936fdd9i9nTM3Idakv3dtF6/r+rRXbKFTBgf5WfazopCBFOGMGACiJ0AYA8HqOLnEsLDR0+GS6dh04o90Hz2jXwTPafzS1UicIqcxJQQAAsIfQBgDwao4ucVy6/qDSMnOVkZ3vUvtNo+uoZeMItYwJV8uYMLVoHK4f1+z32DT6AACURGgDAFSpb5YmatYqKd2cqNGD2zmsV1Bo6KN5W7Tot/12tx9zsCZaefztita65ap4m3Km0QcAeBNCGwCgysxK2KVZCYn/fZ4oX19fS7jJzM7T7oNntGP/Ge1IOqVt+04pL7+wQsfxMUlhIWadTc9xWKesEMYZMwCAtyC0AQCqhKPLHDf8dUz5BYb2HzlX4dkc64cHqnWzumrTrK5ax9bVBU0iFGT2c3liEM6YAQC8AaENAFDpvl6y0+7C0ZK0++BZl9oe0S9OY4d1sLuNiUEAADUBoQ0AUG7OTL9/Ji1bWxJTNG/FHu09fK5S+sHEIACA2oDQBgAoF0fT72dm52n7vlP6M/GktiSmaP/R1HK37etjUlyTcMU3r6/4FvUU37yelqw7wNpnAIBajdAGAHCao/vSflp7QGdSs1XgwvpoV3Rvpnuv6ySzv69VOZc4AgBqO0IbAMAppd2XlnI2y6W2mckRAADHCG0AUIuVdW9aakauNiee1JxliRW6Ly3I7KsOcZG6sFWUjqZkaOHqJJs6zOQIAEDpCG0AUEvZuzft+v6ttOvAGf2+64T+2HVCe5LPyijnFY/tW9ZX59ZRuvCCKLVqFiE/Xx/LtohQM5c5AgBQToQ2AKiFHN2b9p+EXcovqPh9aX+7orVuuSre4fZRA9uooKBAsxISNWpgKwIbAABOILQBQC3z1Y879J+fd9vd5kpgc/as2Y0DWqlOTqKGDGhV4WMBAFCbENoAoAYo6960oykZ2rjjuBb8uk9HUzLK1XaAv686xtVXlzYNdOxUhn74teL3pQEAgPIjtAFANWfv3rSRl1+gbXtPaePO49q047gOnyxfUCtqo0ubBmrXop78/f43DX94He5LAwCgKhHaAKAaq6x700oLYky/DwBA1SK0AUA1Vdq6aZV9bxrT7wMAUHUIbQBQjRQUFGrbvlP6YtFf2n3wbLn2jahj1kVtG6hrfLSSDp/Tt8sSbepwqSMAAN6H0AYAXsLRZCI5eQX6Y9cJrdl6VBv+Oqa0zLxytTt6cFtd3LaB4mIi5ONjkiRd1jlGAQG+3JsGAEA1QGgDAC9QcjKR3NwCNYkO1dptR/X7rhPKyS2oULvcmwYAQPVHaAMAD7M3mYi9SxfLi3vTAACoGQhtAOBBn/2wXXOW7yn3fnFNwtWzQyOdScvRwtWsmwYAQE1GaAOAKpaemavfth7Vt0t369ipTKf28TFJ7VtGqkeHhurRoZEa1Au2bIsIZd00AABqMkIbALiZvQlFMrPztG77Ma3647D+3H2i3FPyfzF5sMLrmO1u4940AABqNkIbALhRyQlFDhxNVUGhoY07jisvv7BCbY4e3NZhYCvCvWkAANRchDYAcBN7E4r8uvmIS21yqSMAACC0AYCLDMPQ1NmbtXjtAaf38ffzUdf4aF12YYwOHEvVf37ebVOHwAYAACRCGwBUWMrZLC3fdEjzV+5VakZumfVNJunittG6rHOMenRoqOBAf0nSZYqRn58Pk4kAAAC7CG0A4IC9CUWyc/O1dutRLd14SJsTT8oox3wihiE9e1cPu9uYTAQAADhCaAMAO0pOKHL8dKZkSKu3HFZWTkGF2hw9uG2p25lMBAAA2ENoA4AS7E0o8vP6gy61yeWOAACgoghtAFDMVz/usDspiCNBZj/1vrCxBnRrpi17TmrmT7ts6hDYAACAKwhtAGo9wzCUeOisPpizWXuSz5VZ3ySpS5sG6t+1qS7p0FCBAec/Stu3rC+TycSEIgAAwK0IbQBqrdSMXK3YdEgJ6w9q/9FUp/czJD03rqfdbUwoAgAA3I3QBqBGKzkDZGGhoc2JJ7Vk3QGt3XZM+QWF5W6TCUUAAEBVIrQBqLFKzgC5JfGkjp/O1IkzWRVuk8sdAQBAVSO0AaiR7M0AuXXvqVL3ia4XrIHdmykjK0/zVu612U5gAwAAnkBoA1DjfLHwL327LNGpuv5+PurZsZEGdY9Vxwsi5eNjkiQFB/kzoQgAAPAKhDYANUbyiTS9MfN3JR46W2bdiFCz/nZFa/W9qIlCgwNstjOhCAAA8BaENgDVmmEY+nP3SX3/yz5t3HHc6f3OpeXo6t4tS63DhCIAAMAbENoAVAslZ4HMySvQik3J+v6XvTp4LK3c7d1cxgyQAAAA3sLH0x1wRU5Ojl5++WV16dJFoaGhioiIUM+ePfXhhx+qsLD803gXt2DBAo0YMUIxMTEym82qX7+++vTpo2nTprncNoDyKZpUxND5WSAnvv+rxv5rid779s8KBTbuTwMAANVJtQ1tKSkp6tatm55++ml1795dP/74o+bOnatGjRpp/PjxGjhwoLKzs8vdbm5urq6//npdc801SkxM1Kuvvqo1a9boyy+/VGhoqMaNG6f+/fsrIyOjEt4VgJLszQK5bd8ppWbk2q1fJ8hf1/dvpRF94+xuJ7ABAIDqptpeHnnDDTdo69atevjhh/XWW29Zyi+//HKNGDFC3333ncaPH6/p06eXq90HHnhAc+bMUcOGDbV69WqFh4dbtl111VW69tprtWDBAo0bN04zZsxw19sBYIe9wOZIkwZ1dM1lLXX5xU0VaD7/0cYMkAAAoCaolmfa5syZoxUrVigwMFCTJ0+22mYymTRlyhRJ0ueff65NmzY53W5SUpI++eQTSdIjjzxiFdiK2i463syZM7Vhw4aKvwkApXr/2z+dCmyN6odo8t099P6T/XVVrxaWwCadn0hk9OC2MonABgAAqq9qGdo+/vhjSVL//v0VERFhsz0+Pl7x8fEyDEOffvqp0+0uXbrUcr9a9+7d7dbp0qWLQkJCJEmfffZZ+ToOoEx7k8/qX5+s0+K1B5yqf+xUhi5uG21ZX62kUQPb6PvXryWwAQCAaqvahbbc3FwtXbpUktStWzeH9Yq2LVy40Om2T5w4YXneuHFju3VMJpPq168vSVq5cqXTbQMo3f6jqfq/z9brkTdXav1fx5zej1kgAQBATVft7mnbsWOH8vLyJEnNmzd3WK9o24EDB3Tu3DmbSx3tCQ0NtTxPSUlRmzb2/zJ/6tQpSdKuXbtUUFAgX19fJ3sP1G4lp+2XpEPH0zTzp536dfORcrfHJY8AAKA2qHah7eDBg5bnUVFRDusV35acnOxUaOvatavl+fbt23XppZfa1Dl06JBl5sj8/HydO3dO9erVc6rvQG1WfFKRGYt3KjUjV2kZuVr5R7IMw/4+YSEBiosJ1x+7T9psI7ABAIDaotqFtrS0/63JFBgY6LBe8W2pqalOtd2jRw917txZf/75p9566y3deeedNmfR3nzzTavXWVlZDtvLyclRTk6O5bVR7Jtp0dlCTyk6vqf7geqpvOPnm6WJmpWQaFW24Jd9DuvXCfLXtX1aaEiv5goy+9nsP2pgK43s15LxW03x+QNXMH5QUYwduKIyxk952qp2oa0ymUwmzZgxQ3379tWOHTs0cuRITZkyRW3atFFKSoqmTZumd999V23bttXOnefPGNSpU8dhe1OmTNFzzz1nU56dna1FixZV2vsoj4SEBE93AdWYM+PnjwPSJufmFFGAr9SxidQ+Jk8BWbu1fOluSVIdSRfHnm/n4lipTk6iFi1KLL0xeD0+f+AKxg8qirEDV7hz/GRmZjpdt9qFtuL3nZW2eHbxbWFhYU63365dO/3+++969tlnNWfOHLVr106S5OPjo0suuUSLFy/W4sWLtXPnTvn6+lr1p6SJEyfqscces7yOj4/XkSNHFBgYqCFDhjjdp8qQl5enhIQEDRw4UP7+/h7tC6ofZ8fPN0sTtelA2eHKz9dHI/q21LDeLVQn2H57nv0/Bu7E5w9cwfhBRTF24IrKGD/OXg0oVcPQ1qxZM8vzkydt73Oxt61JkyblOkbTpk316aefatq0aTpy5Iiys7MVHR1tCX8zZ86UdD6E+fg4noDTbDbLbDZbXptM/5uS3Fs+LPz9/b2mL6h+yho/JS+JdCS/oFC3DW3vrm6hmuDzB65g/KCiGDtwhTvHT3naqXahLT4+Xv7+/srLy9P+/fsd1ivaFhsb69QkJPb4+vqqadOmNuWJiee/iPbo0aNC7QI13fHTmfp0wTan649m2n4AAACHql1oCwgI0IABA7R48WJt3LjRYb0NGzZIkoYOHerW42dlZVnaHj16tFvbBqq7nLwCzV2WqNnLEpWbX+jUPswCCQAAULpqt7i2JN11112SpKVLl+rcuXM223fu3KkdO3bIZDJp7Nix5Wr7nXfeUa9evVRQUGB3+9dff63s7Gz17NlTffv2LX/ngRrIMAz9tuWI7nt5qWYu2UVgAwAAcKNqGdpGjhypvn37Kjs722Z2RsMwNGnSJEnSmDFjdPHFF1ttX7BggaKiotShQwe7l1ceOXJEa9as0Zw5c2y2HThwQJMmTVJYWJimTZtmdY8aUJt8szRRH686/98Dx1L1j3//pimfb9CJM7ZLYIQGB+j+6y/UTYOswxmBDQAAwDnV7vLIIrNnz1b//v315ptvKisrS7fccotyc3P1/vvva968eerfv78++OADm/0++ugjpaSkKCUlRXPnzrWa3bG4O++8U3v27FG/fv1kGIZ+/fVXvfLKKwoMDNRPP/2k9u2ZNAG106yEXZYJRmYlJOo/PyfaXRzbxyQNubSFRl/ZVnWCA86X+Zg0c/FO3UxgAwAAcFq1DW2RkZHasGGD3nrrLX399df68ssv5evrq/j4eE2dOlX33HOP3Zkdx40bpzVr1ig6OlrXXXedzfa//e1vysvL06+//qp3331XkydPVlhYmFq1aqWnnnpK9913X6nT/AM12ayEXZqxeKdVmb3A1iGuvu4Z0UnNG1kvtzFqYBvCGgAAQDlV29AmnZ9Sf8KECZowYYLT+wwbNkwpKSkOt3fp0kVdunRxR/eAGsVeYCspMjxQY6/poN4XNubyYQAAADep1qENQNWYsXiHZiXsLrPegG7NdFnnmCroEQAAQO1RLSciAVB1tu5JcSqwSdI3PztXDwAAAM7jTBsAu9IyczV9wXYlrD/o9D43s0g2AACA2xHaAFgxDEMr/zisj7/bqnPpuU7vxxT+AAAAlYPQBsDi2KkMfTBni37fdcLu9vrhgTp1LtumnMAGAABQeQhtQC02K2GXZi7eqVGD2igwwE8zl+xUTm6BTb0gs5/GDInX4F4t9O3S3VazSBLYAAAAKhehDailik/h//WSXQ7rXdK+oe69rpMiI4IknV9rraCgQLMSEjVqYCsCGwAAQCUjtAG1kDNrrtULC9S913VUz46NbbbdOKCV6uQkasiAVpXVRQAAAPwXoQ2oZZwJbEN6NddtQ9opJMi/inoFAAAAR1inDahFnAlsklQ3LJDABgAA4CUIbUAtUVhoOBXYJGmmk/UAAABQ+bg8EqgFMrPz9ObXvztdn0WyAQAAvAehDajhjqSk68Xp63XwWJpT9ZnCHwAAwLsQ2oAa7PddJ/TKlxuVkZVnVW4ySYZhW5/ABgAA4H0IbUANZBiG5q/cq89+2K7CEuEsNNhfE27tph0HTrNINgAAQDVAaANqmJy8Ar337Z9asSnZZlvzRmF65o7ualg/RBe2jpJ0ftKRmwlsAAAAXovQBlRzsxJ2WYLXgK7N9H+frdOe5HM29Xp1aqRHRl2kIPP//rcfNbANYQ0AAMDLEdqAaqz4umszFu/UnGWJys4tsKl3y+C2uvGK1jKZTFXdRQAAALiI0AZUU/YWyi4Z2ILMfnr85ot0SYdGVdk1AAAAuBGhDaiG7AW2khpHhujvYy9R0+jQKuoVAAAAKoOPpzsAoHycCWySdOmFjQlsAAAANQChDahmZjoR2CRp9tLESu4JAAAAqgKhDahmhlzawql6Nw9uW8k9AQAAQFUgtAHVyObEk1q+6VCZ9VgoGwAAoOZgIhKgmli+6ZDe+c8fyi8wSq1HYAMAAKhZCG2AlzMMQ98uTdSXP+6w2RYZEaSUs1mW1wQ2AACAmofLIwEvVlBQqPdnb7Yb2Hp1aqQPnx6g0YPbyiQCGwAAQE3FmTbAS2Xl5OuVLzdq447jNtuu7ROnscPay8fHpFED2xDWAAAAajBCG+CFzqRm6/lP1mpP8jmrcpNJuvOaDrq2T5yHegYAAICqRmgDvMSshF2auXinhvZuofV/HdeJ05lW2wP8fPTY6It1aafGHuohAAAAPIHQBniBWQm7NOO/i2b/8GuSzfbQYH/9Y2wPxbeoV9VdAwAAgIcR2gAPKx7Y7ImuF6zJd/dQkwahVdgrAAAAeAtCG+BBZQW2emGBevWhy1Q3NLAKewUAAABvwpT/gIeUFdgk6XRqtn5ae6CKegQAAABvRGgDPGRmGYGtvPUAAABQMxHaAA8wDENtmtd1qu7Ng9tWcm8AAADgzbinDahihYWGPpq/VTv3nymz7ujBbVk4GwAAoJYjtAFVqKDQ0NTZm7VkXdn3qRHYAAAAIHF5JFBlCgoK9das320Cm49J6tWxkVUZgQ0AAABFCG1AFcgvKNSrMzZpxaZkq3JfH5OevLWrJt7eXaMHt5VJBDYAAABY4/JIoJLl5Rfo5S82at32Y1blfr4+evq2rrqkw/mzbKMGtiGsAQAAwAahDahE2bn5mvLZBv2+64RVeYCfjybd0V0Xt432UM8AAABQXRDagEqSlZOvf32yTlv3pliVmwN89Y+xl+jCVlEe6hkAAACqE0Ib4GazEnZpxuKdiowIUsrZLKttQWY/Tb67h9q1qO+h3gEAAKC6IbQBblQU2CTZBLaQIH89P66nWjdzblFtAAAAQCK0AW5TPLCVFBocoBfu7aWWMeFV3CsAAABUd0z5D7hBaYFNkvp0iSGwAQAAoEIIbYCLygpskrRwdZJmJeyqoh4BAACgJiG0AS6aWUZgK289AAAAoDhCG+CiUYOcWxD75sFtK7knAAAAqIkIbYCLzP6+ZdYZPbitRg10LtwBAAAAxRHaABecOJ2pr8u4V43ABgAAAFcQ2gAXfDR/q3JyCxxuJ7ABAADAVYQ2oILWbjuqdduPWZVd1bO5Rg9uK5MIbAAAAHAPFtcGKiA7J18fzd9qVRZRx6zbhsSrTnAAYQ0AAABuw5k2oAJmJezSyTNZVmVjr2mvOsEBHuoRAAAAaipCG1BOB46mav7KvVZlnS6IVL+LmnioRwAAAKjJCG1AORQWGpo6Z7MKCg1LmZ+vSfde10kmk8mDPQMAAEBNRWgDymHZxoP6K+m0Vdl1l7dS0+hQD/UIAAAANR2hDXBSakauPl3wl1VZdL1g3XhFaw/1CAAAALUBoQ1w0mc/bFdaZq5V2b3XdZLZ39dDPQIAAEBtQGgDnPBX0iklrD9oVdarUyN1jY/2UI8AAABQWxDagDLkFxRq6uzNVmVBZl/dfW1HD/UIAAAAtQmhDSjD96v26cCxNKuym6+MV2REkId6BAAAgNrE5dDWsmVLvfTSS+7oC+B1TpzJ1NdLdlqVtWgcpmG9W3ioRwAAAKhtXA5t+/fv1/Hjx93RF8DrTJu/Vdm5BZbXJpN03/UXyteXk9QAAACoGm755vnOO++oT58++vzzz5WZmemOJgGPmpWwS8Me/05rtx2zKh90SazaxtbzUK8AAABQG7kltMXFxWnbtm2644471KhRI91zzz1av369O5oGqtyshF2asXinTXl4nQCNGdrOAz0CAABAbeaW0DZ8+HAdPXpUX3zxhbp06aJp06apZ8+e6tChg9566y2lpKS44zBApXMU2CRp7LD2Cg0OqOIeAQAAoLZzObT961//0lVXXSWz2axbbrlFK1asUGJiop566imdOXNGjz32mJo0aaIbb7xRixcvlmEY7ug34HalBTZJOn6aS38BAABQ9VwObc8884wuv/xyq7K4uDhNmTJFBw8e1Lx58zRo0CDNmzdPQ4cOVWxsrJ599lnt37/f1UMDblNWYJOkmT/t0qyEXVXUIwAAAOC8Sp0Cz9fXV40aNVJ4eLh8fX1lGIaSk5P1wgsv6IILLtAVV1yhWbNmKScnpzK7AZRpZhmBrbz1AAAAAHeplNCWm5urzz//XN27d1fPnj01c+ZM5eXlyWQyyWQyyTAMFRYWavny5Ro9erQaNWqkf/zjHzp37lxldAco082D27q1HgAAAOAuLoc2X19fTZgwQZJ08OBBTZw4UU2aNNHYsWO1ceNGGYbh8D62xo0b6/7771d4eLhefPFFde3aVQcPHnS1S0C5jRrYRu1blj6V/+jBbTVqYJsq6hEAAABwnsuhzTAMbdmyRcOHD1dcXJxeeeUVpaSkWIJa0dm1orqSNGDAAM2ePVv79+/XO++8oz179ui1117Tvn37NHHiRFe7BJTbkZR07Tpw1uF2AhsAAAA8xS2XRy5ZskQLFixQQUGBDMOwBLWiSyENw1B4eLgefvhh7dixQwkJCbruuuvk6+sr6fzZuscee0zDhg3TsmXLnD5uTk6OXn75ZXXp0kWhoaGKiIhQz5499eGHH6qwsNCl9/TDDz9o+PDhiomJUUBAgEJCQtS+fXs99NBD2rdvn0ttw/t8/N025RfYHzMENgAAAHiS2+5pKx7Wil4bhqGLLrpIH3/8sQ4fPqw333xTrVu3dthGdHS0UlNTnTpeSkqKunXrpqefflrdu3fXjz/+qLlz56pRo0YaP368Bg4cqOzs7Aq9jzvvvFPDhg3T8uXL9eijj2r58uX69ttvddFFF+ndd99Vhw4d9P3335e7bXinjTuOa8Nfx63K2javK5MIbAAAAPA8P3c1VDysBQYG6sYbb9R9992n7t27l7lvfn6+1q5dqzlz5qhx48ZOHe+GG27Q1q1b9fDDD+utt96ylF9++eUaMWKEvvvuO40fP17Tp08v1/v4/PPP9emnn8pkMunHH39Ur169LNuGDBmi4OBgffTRR7r11lu1d+9eRUZGlqt9eJe8/AJ9NH+rVVlEqFnP3d1TwYH+HuoVAAAA8D9uPdPWsmVLvfrqqzp8+LA+++wzpwKbJN16663q27evTp8+rW7dupVZf86cOVqxYoUCAwM1efJkq20mk0lTpkyRdD6Abdq0qVzv48svv5Qkde3a1SqwFXn44YclSampqVq4cGG52ob3mb9yr46mZFiV3T60HYENAAAAXsMtZ9o6d+6s//u//9OVV15Zof2vv/56tW17fir1K664osz6H3/8sSSpf//+ioiIsNkeHx+v+Ph47dixQ59++qkuvvhip/ty+PBhSVKLFi3sbm/evLnl+bFjx5xuF97n1LksffPzbquyNrF1dfnFTT3UIwAAAMCWW0LbgAEDKhzYJGnkyJEaOXKkU3Vzc3O1dOlSSSr1rFy3bt20Y8cOLVy4UO+//77TfYmNjdWuXbscBrLi5RdccIHT7cL7fLpgu7JzCyyvTSbp3hGd5ONj8mCvAAAAAGsuXx7Zt29ftWrVyh19ccqOHTuUl5cnyfqsV0lF2w4cOFCuRbtvu+02SdK6devszhL59ddfSzof2IYOHep0u/Au2/amaNUfh63KBl0SqwuaRnimQwAAAIADLp9pW7lypVP3oblL8cW3o6KiHNYrvi05OVnh4eFOtT969Gj99ddfeumllzRs2DC9++676tGjh1JTU/X111/rhRdeUPfu3fXVV18pMDCw1LZycnKUk5NjeV18kfGi4OkpRcf3dD88oaCgUB/O3WJVFhLop1FXXFArfx4VUZvHD1zH+IErGD+oKMYOXFEZ46c8bbnl8sjt27friy++KNc+fn5+ioqKUteuXVW3bl2n90tLS7M8Ly00Fd/m7DICRV588UVdf/31evzxxzVgwABLeUBAgB566CE98cQTio6OLrOdKVOm6LnnnrMpz87O1qJFi8rVp8qSkJDg6S5Uub+OSAdKXP3aKSZfq1ct9UyHqrHaOH7gPowfuILxg4pi7MAV7hw/mZmZTtd1S2hbvHixFi9eXKF9fX19dd111+nNN99Uo0aN3NEdl+Tm5mry5Ml67bXX1KRJE/373/9Whw4dlJqaqmXLlumtt97S1KlT9corr+j+++8vta2JEyfqscces7yOj4/XkSNHFBgYqCFDhlT2WylVXl6eEhISNHDgQPn7156ZElMzcjXrtZWS/veXjdiGoXps7KXy9XXbZKo1Xm0dP3APxg9cwfhBRTF24IrKGD/lObHktnXail/6Vx75+fn69ttvtW7dOq1atUpNm5Y+c19oaKjleWmLZxffFhYW5nR/brjhBn3//fdq2bKltmzZopCQEMu2wYMH6/LLL9eQIUP0wAMPyNfXV/fee6/Dtsxms8xms+V10Vp2krzmw8Lf399r+lIVvk7YrvQs61PR91zXSYGBZgd7oDS1bfzAvRg/cAXjBxXF2IEr3Dl+ytOO204tmEymCj8Mw9CBAwd0xx13lHmcZs2aWZ6fPHnSYb3i25o0aeLUe/jtt9/0/fffS5L+/ve/WwW2IldddZUuu+wySdILL7zgVLvwDnsOndWSdQesyvp0jlHHOBZIBwAAgPdyS2grOstmGEaZD3v1iixfvly//PJLqceKj4+3pNL9+/c7rFe0LTY21ulJSH777TfL806dOjmsd+GFF0o6v6bbiRMnnGobnlVYaOjf87ao+Alhc4Cv7hjW3nOdAgAAAJzg8uWR06dP15dffqnly5frsssuU69evdSkSRMFBwfLx8dHhYWFyszMVHJysn777Tf9+uuv6t69u8aPH6/s7GwdPXrUKqz95z//sZzJsicgIEADBgzQ4sWLtXHjRof1NmzYIEnlmpa/Ipd4+vm57QpTVKIVvx/SzgNnrMr+dkVrRUYEeahHAAAAgHNcThzNmjXTunXrtHLlSvXu3bvM+qtWrdJVV12l+vXr6+qrr5YkTZ48WR999JHuvfderV+/vsw27rrrLi1evFhLly7VuXPnbM6k7dy5Uzt27JDJZNLYsWOdfi8dOnSwPN+yZYsuvvhiu/U2b94s6fxll/Xq1XO6fXhGZnaepv/wl1VZo8gQDe8b56EeAQAAAM5z+fLIt956SyNGjHAqsElSnz59NHLkSE2dOtWqfNy4cWrfvr2SkpLKbGPkyJHq27evsrOzbabUNwxDkyZNkiSNGTPGJngtWLBAUVFR6tChg83llVdccYXatGkj6fy0/xkZGTbH/vHHHy1nBR944IEy+wrP+3rJLp1Ny7Equ/vaDvL38/VQjwAAAADnuRza1q1bp6Cg8l1iFhQUpE2bNtmU9+jRw2odttLMnj1bHTt21Jtvvqnx48dr9erVWr58uW644QbNmzdP/fv31wcffGCz30cffaSUlBRt375dc+fOtdrm7++vefPmqVmzZtq7d686duyoadOm6bffftNPP/2kCRMm6Nprr5V0PhA+8cQT5XrfqHr/nrdF81futSrrGh+tbu0aeqhHAAAAQPm4HNrOnDmj+fPn6/Dhw07VP3jwoObNm6ezZ8/abAsODraaIr80kZGR2rBhg1566SWtWbNGV155pYYPH67k5GRNnTpVCQkJdhffHjdunOrXr6927drpuuuus9keHx+v7du364033lBsbKwmTZqkvn37asSIEZo9e7auv/56LVmyRJ999pl8fTlT482+XrJTP/xqfebWz9dHdw/v4GAPAAAAwPu4fE9beHi4UlJS1KFDB91www26+OKLFRMTo9DQUPn5+Sk/P19paWlKTk7Whg0bNGfOHKWmpioy0naa9SNHjqh+/fpOH9tsNmvChAmaMGGC0/sMGzZMKSkppdapU6eOHn30UT366KNOtwvvMithl2b+tMumvE1shBpH1vFAjwAAAICKcTm0de7cWT///LPOnTunTz75RJ988kmp9Q3DkMlkUpcuXazK09PTtWzZMl1wwQWudgm13KyEXZqxeKfdbdv3ndashF0aNbBNFfcKAAAAqBiXL48cPXq0JFkWyXZmnbbi+0nS6dOnddNNN+ns2bOlro8GlKW0wFZkxuKdmpVgexYOAAAA8EYuh7bbbrtNl112meUMWlkPSbrssst02223STo/MUjDhg21aNEiSdKll17qapdQi80sI7CVtx4AAADgaS6HNpPJpO+//16DBg2yOZtWXNG2K6+8Ut99952lPCgoSJdccol69eqlXr16qW/fvq52CbXYzYPburUeAAAA4Gkuhzbp/GQkixcv1ty5czVkyBCFh4dbXRIZERGhq6++WvPnz9ePP/5otRj2rbfeql9++cXyaNGihTu6hFrqmstaKsC/9GE9enBb7mkDAABAteHyRCTFDR8+XMOHD5ckpaamKi0tTaGhoQoLC3PnYQCHvlu1T7l5hQ63E9gAAABQ3bh8ps3X19fq8eyzz0qSwsLCFBMTQ2BDlUnLzNX8lXscbiewAQAAoDpyObQVvwyyc+fO6tCBhYvhGfNW7FFmdr7ltckkDb20hUwisAEAAKD6cvnySH9/f+Xn5+vJJ5/USy+95I4+AeV2Ni1HC37ZZ1V2WecY3XtdJ917HctIAAAAoPpy+UxbTEyMJOnGG290uTNARc1Znqjs3ALLax+TdNMgzqwBAACg+nM5tA0cOFCSlJKS4vQ+CQkJev755109NCBJOnUuS4tWJ1mVXd61qZo0CPVQjwAAAAD3cTm0TZgwQREREXr++eeVk5Pj1D5LlizRc8895+qhAUnSt0sTlZv/vxkjfX1M3L8GAACAGsPl0NayZUv99NNPOnXqlC6++GJ9/fXXOnPmjDv6BpTpxOlM/bR2v1XZwEti1bB+iGc6BAAAALiZyxORtGzZUpKUm5urI0eO6JZbbpEk1a1bV6GhoTKZTDb7EOrgLv/5ebfyCwzLa38/H/3titYe7BEAAADgXi6Htv3791sFM8M4/wX69OnTOn36tMP97IU5oDyOpKTr5w0HrcoG92yuyIggD/UIAAAAcD+XQ1sRwzCcDmJFwQ5wxawlu1RY+L+xFODvqxv6t/JgjwAAAAD3c/metiKcOUNVOngsVSt+T7YqG9a7heqGBXqoRwAAAEDlcNuZtkaNGsnf39+pumfOnFFaWpq7Do1aaOaSXSp+wjbI7KsR/S7wXIcAAACASuK20LZs2TK1bu3cBBBPPvmk3njjDXcdGrVM0pFzWr35iFXZNX3iFF7H7KEeAQAAAJXHbZdHlodhGNzXhgqbsXin1euQIH8N78tZNgAAANRMLoe2pKQk7du3zzL1vzNee+01FRYWll0RKGH3wTNat/2YVdmIfnGqE+TcpbkAAABAdePy5ZGxsbHu6AfglJJn2cJCAjSst/N/MAAAAACqG7fd01Zk3bp1Wr58ubZt26azZ8/qhx9+kCT9/PPPuuSSSxQaGuruQ6KW2L7vlH7fdcKqbOTlrRQcyFk2AAAA1FxuC21btmzRXXfdpU2bNkmyXbftySef1P79+/Xcc8/poYcectdhUUsYhqGvFu+wKqsbataQS5t7pkMAAABAFXFLaFu1apUGDx6snJwchxOMhIWF6dy5c3r00Ud14sQJvfDCC+44NGqJLYkp2rb3lFXZDQNaKzDA7SeLAQAAAK/i8kQk6enpuuGGG5SdnS3p/CLb9hbaXrlypRISEtS0aVO99NJL2rZtm6uHRi3x9ZKd+vu/f7Mqi4wI0uCe3E8JAACAms/l0PbRRx/p5MmTMplMlqn8HZ1tGzBggBISEuTr66vp06e7emjUArMSdmnmT7tsykcNbC1/P18P9AgAAACoWi5fW7Zo0SJJko+Pj3r06KH27durXr16WrlypdauXWtTv1WrVho4cKBWrVrl6qFRw81K2GUzW6Qk1Qny14BuzTzQIwAAAKDquRzatm/frsjISK1atUpt27a1lD/55JN2Q5sktWzZUr/99pvdbYDkOLBJUnpWnmYvS9SogW2quFcAAABA1XP58sjTp09r5MiRVoGtLCdOnFB6erqrh0YNVVpgKzJj8U7NSrC9bBIAAACoaVwObUFBQTpx4kTZFf8rOTlZixYtUlBQkKuHRg01s4zAVt56AAAAQHXmcmhr3Lix5s+fr2eeeUb79+93WO/MmTP66quvdNlllyk9PV3R0dGuHho11M2DnTtr62w9AAAAoDpz+Z62nj17aufOnXrppZf00ksvKTg4WPXr11dqaqokKS4uTqmpqTp9+rSk/y263b17d1cPjRpq1MA2MgoNzVzi+PLH0YPbck8bAAAAagWXz7TdeuutlueGYSgjI0MHDx7U2bNnZRiGkpKSdOrUKZulAG688UZXD40arF3L+g63EdgAAABQm7gc2vr166dBgwZZzqAVf0iyKTOZTOratauuueYalzuPmuvnDQftlhPYAAAAUNu4HNokacaMGerYsaPNotpFwa2IYRhq1qyZvv32W3ccFjVUZnaeftty1KacwAYAAIDayC2hrX79+vrtt9/06KOPKjAw0HIpZPGHn5+f7rzzTm3cuFHNmrEwMhz75c8jys0rsLz29THpy8mDCWwAAAColVyeiKRISEiIXn/9df3rX//SL7/8osTERKWmpiosLEwtWrTQZZddprCwMHcdDjXY0hKXRnaNj1ZEqNlDvQEAAAA8y22hrUhwcLCuvPJKXXnlle5uGrXA4ZPp2rH/tFXZFd05MwsAAIDayy2XR5bXhx9+qP79+3vi0PByJc+yhdcJUNd41vQDAABA7eWR0LZ3716tXLnSE4eGFysoNLRs4yGrsn4XNZWfr0eGKQAAAOAV3HZ5ZHp6uhYsWKAtW7bo9OnTysvLc1h3w4YN7josapDNu0/q1Llsq7IB3Zp6qDcAAACAd3BLaPvyyy/14IMPKi0tzan6RWu6AcWVXJstrkm4WjQO91BvAAAAAO/gcmj7+eefdfvtt9us0QaUR3pmrtZus16b7YpuTEACAAAAuBzaXn311XKfOSPgoaRVfx5WXn6h5bWfr4/6dGniwR4BAAAA3sHl0LZhwwaZTCaCGFzy83rrSyMvad9QYSEBHuoNAAAA4D1cDm2ZmZmSpEceeUQjRoxQo0aN5O/vX+o+//d//6ePP/7Y1UOjhjhwLFWJh85alTEBCQAAAHCey6Gtfv36uuiii/TGG284vU9sbKyaNeN+JZy3dIP1NP/1wsy6qE0DD/UGAAAA8C4uL4DVo0cPxcbGlmufSZMmKSkpydVDowbILyjU8k3Woe3yi5vKl7XZAAAAAEluCG333XeffvzxR2VnZ5dd+b8+/PBDDRgwwNVDowb4fdcJnU3LsSobwKyRAAAAgIXLoW3AgAG6+eab1bNnT82aNUv79+9Xbm5uqfvs3btXK1ascPXQqAFKTkDSplldNY0O9VBvAAAAAO/j8j1tvr6+ks5P4z969GiXO4Ta41x6jjb8dcyqjAlIAAAAAGsuh7bia7SVZ9r/8qzrhppp5R/Jyi/435gJ8PPRZazNBgAAAFhxObQVYXFtlNfS9dYTkPTo2Eh1gkpfLgIAAACobdwW2ghiKI99h89p35FzVmVMQAIAAADYcktoa9KkiVq2bOl0/b179+rw4cPuODSqqaUbrCcgiQwP1IWtojzUGwAAAMB7uSW0/e1vf9Mrr7zidP0nn3yyXItxo2bJyy/Uit+Trcou79pUvj7c5wgAAACU5PKU/7GxsapXr1659qlbt66aNeNSuNpq445jSs2wXhbiCi6NBAAAAOxy+UxbUlJSufeZNGmSJk2a5OqhUU39XGICknYt6qlxVB0P9QYAAADwbk6Htueff97qdXBwsJ544gmn9s3IyNCwYcOsykwmk5YuXers4VFDnEnL1sadx63KmIAEAAAAcMzp0DZ58mSZTCbLLJGRkZFOh7b8/HytWLHCqox12mqnFZuSVVj4v5lGzQG+6n1hYw/2CAAAAPBu5bo8cvDgwbrxxhslSYGBgU7vFxYWZnUZ5YMPPqiFCxeW59CoAQzD0M8lZo3s1bGRggNZmw0AAABwpFyhrX379hozZoxV2apVqxzW79Onj6TzZ9ViY2Mt5SEhIeU5LGqIPclndfBYmlXZFd25NBIAAAAojcsTkfTr18/hpY4FBQWuNo8a5IM5W6xeN6gXrA4tIz3UGwAAAKB6cDm03XbbbZbQ9vnnn6t9+/bq2rWryx1DzTJj8Q4lHjprVTaga1P5sDYbAAAAUCqXQ9tnn31mef7555/rqquuKtdC26j5ZiXs0qyE3Tblmdl5HugNAAAAUL24vLg2UJpZCbs0Y/FOu9u+W7VPsxJ2VXGPAAAAgOqF0IZKU1pgKzJj8U6CGwAAAFAKQhsqzcwyAlt56wEAAAC1UbnuaUtNTdXBgwftbitadLu0OkUyMjLKc1hUUzcPblvmmbaiegAAAADsK1domzZtmqZNm+ZyHdQOowa2kaRSg9vowW0t9QAAAADYKvflkYZhOHyUtb1kXdR8owa2Ub0ws91tBDYAAACgbOWe8t/RQtrS+cBW2vbi9VA75BcUKi3Tdmp/AhsAAADgHLdOROJMYEPtknwiXXn5hVZl1/dvRWADAAAAnFSuM22cIUN57Tt8zup1VN0gjRnazkO9AQAAAKqfcoW2cePGaeLEiS4f9MEHH9TChQtdbicnJ0dvvfWWZs2apT179sjX11fx8fEaM2aMxo0bJx+f8p9ILO/ZQoJs6ZKOWIe2lo3DPdQTAAAAoHoqV2gLCwtTbGysywcNCQlxuY2UlBT1799fW7du1bhx4/Tuu+8qNzdX7733nsaPH69vv/1WCxcuVGBgYLnbNpvN8vNz/KMpLCxUVlaWWrRo4cpbqBVKnmlrQWgDAAAAysXpU1HXXnutOnTo4JaDduvWTddcc41Lbdxwww3aunWrHn74Yf373/9W79691b9/f82ZM0fXXnutli1bpvHjx1eo7Q8//FDp6ekOH++8844k6f7773fpPdR0hmHYnmmLCfNQbwAAAIDqyenQNm/ePN12221uOehjjz2mefPmVXj/OXPmaMWKFQoMDNTkyZOttplMJk2ZMkWS9Pnnn2vTpk2udNWud999V8HBwbrzzjvd3nZNknI222bmSM60AQAAAOXj1tkjq8rHH38sSerfv78iIiJstsfHxys+Pl6GYejTTz8tV9tbt27V8OHDHW5ftWqVtmzZoltvvdXusfE/Jc+yhQT6KbpesId6AwAAAFRP1S605ebmaunSpZLOX2bpSNG28k540qFDh1LD2Lvvvivp/GQqKN2+EqGtRUw4y0IAAAAA5VTtQtuOHTuUl3f+krvmzZs7rFe07cCBAzp37pzDeuVx6NAhzZ8/X/3791f79u3d0mZNVnISEmaOBAAAAMqvXLNHeoODBw9ankdFRTmsV3xbcnKywsNdDwwffPCB8vPznT7LlpOTo5ycHMvr4ssDFAVPTyk6fmX2Y9/hs1avm0XX8fj7hntUxfhBzcX4gSsYP6goxg5cURnjpzxtVbvQlpaWZnle2nT+xbelpqa6fNzs7GxNmzZNsbGxGjZsmFP7TJkyRc8995zdthYtWuRyn9whISGhUtrNyZeOn7YuO7xvixad2FIpx4NnVNb4Qe3A+IErGD+oKMYOXOHO8ZOZmel03WoX2jxl1qxZSklJ0VNPPSVfX1+n9pk4caIee+wxy+v4+HgdOXJEgYGBGjJkSGV11Sl5eXlKSEjQwIED5e/v7/b2t+87Lf221vLaz9ekUdddKX+/andFLuyo7PGDmo3xA1cwflBRjB24ojLGT3lOLFW70BYaGmp5np2d7bBe8W1hYa6vDVaRaf7NZrPMZrPldfFJOLzlw8Lf379S+nLwRLrV66bRoQoOMjuojeqqssYPagfGD1zB+EFFMXbgCneOn/K0U+1OezRr1szy/OTJkw7rFd/WpEkTl47522+/6ffff9fo0aNVr149l9qqLZIOW//lgPXZAAAAgIqpdqEtPj7ekkr379/vsF7RttjYWJcnIXnnnXckMc1/eZSc7j8uhtAGAAAAVES1C20BAQEaMGCAJGnjxo0O623YsEGSNHToUJeOd+TIEc2dO1d9+/ZVx44dXWqrtsjLL9TBY2lWZS0IbQAAAECFVLvQJkl33XWXJGnp0qV212DbuXOnduzYIZPJpLFjx7p0rA8//FB5eXl66KGHXGqnNkk+kab8gkKrMi6PBAAAACqmWoa2kSNHqm/fvsrOzraZUt8wDE2aNEmSNGbMGF188cVW2xcsWKCoqCh16NCh1MsrJSk3N1cfffSRmjVrpmuvvdat76EmSypxaWSDesGqE8QNvwAAAEBFVMvQJkmzZ89Wx44d9eabb2r8+PFavXq1li9frhtuuEHz5s1T//799cEHH9js99FHHyklJUXbt2/X3LlzSz3GN998o+PHj+u+++5zepp/SHsPW4e2lo1dn70TAAAAqK2qbWiLjIzUhg0b9NJLL2nNmjW68sorNXz4cCUnJ2vq1KlKSEiwu/j2uHHjVL9+fbVr107XXXddqcd49913FRgYaLkcE84pOXNkSy6NBAAAACqs2q3TVpzZbNaECRM0YcIEp/cZNmyYUlJSnKq7bt26inat1jIMw2bmSCYhAQAAACqu2p5pg3c6eSZLGVl5VmWcaQMAAAAqjtAGtyp5lq1OkL+i6gZ5qDcAAABA9Udog1sllZyEJCZcJpPJQ70BAAAAqj9CG9zK5n42Lo0EAAAAXEJog1vtO1Ji5sgYpvsHAAAAXEFog9ukZ+bqxOlMqzLOtAEAAACuIbTBbZJKnGXz8/VRkwahHuoNAAAAUDMQ2uA2Je9na9YwVP5+DDEAAADAFXyjhtvsKzlzJJdGAgAAAC4jtMFtko7YTvcPAAAAwDWENrhFXn6hDh1PsyojtAEAAACuI7TBLQ4dT1N+gWFV1rwR0/0DAAAAriK0wS1K3s/WsH6wQoL8PdQbAAAAoOYgtMEtSt7PxvpsAAAAgHsQ2uAWe0vOHMn9bAAAAIBbENrgMsMwbGeO5EwbAAAA4BaENrjs+OlMZWbnW5VxeSQAAADgHoQ2uKzkWbbQ4ABFRgR6qDcAAABAzUJog8v2HU61et0yJkwmk8lDvQEAAABqFkIbXMbMkQAAAEDlIbTBZftKTkLCzJEAAACA2xDa4JK0zFydPJNlVcbMkQAAAID7ENrgkpKXRvr7+SimQR0P9QYAAACoeQhtcMm+EotqxzYMlZ8vwwoAAABwF75dwyUlQxuTkAAAAADuRWiDS5KOWE/3H8ckJAAAAIBbEdpQYbl5BTp0PM2qrAWhDQAAAHArQhsq7ODxNBUUGlZlzRuFeag3AAAAQM1EaEOFJZW4n61RZIiCA/091BsAAACgZiK0ocJsFtVmEhIAAADA7QhtqLCSk5C0iOHSSAAAAMDdCG2okMJCw2a6f860AQAAAO5HaEOFHD+dqaycfKuylswcCQAAALgdoQ0VUvJ+tvA6AaoXFuih3gAAAAA1F6ENFVJy5sgWjcNlMpk81BsAAACg5iK0oUKYORIAAACoGoQ2VIjNmTbuZwMAAAAqBaEN5XYuPUcp57Ktylo2Zrp/AAAAoDIQ2lBu+0uszxbg56OYqDoe6g0AAABQsxHaUG4l72eLbRQmX1+GEgAAAFAZ+KaNcrNZVJv72QAAAIBKQ2hDudnMHEloAwAAACoNoQ3lkpNXoOQT6VZlTPcPAAAAVB5CG8rl4LFUFRYaltcm0/l72gAAAABUDkIbymXfYeuZIxtHhijI7Oeh3gAAAAA1H6EN5ZJU4n62FlwaCQAAAFQqQhvKZd32Y1av0zNzPdQTAAAAoHYgtMFpX/+0Uylns6zK/kxM0ayEXR7qEQAAAFDzEdrglFkJuzRzif1wNmPxToIbAAAAUEkIbSjTrIRdmrF4Z6l1CG4AAABA5SC0oUwzywhs5a0HAAAAwHmENpTp5sFt3VoPAAAAgPMIbSjTqIFtNLqMQDZ6cFuNGtiminoEAAAA1B6ENjhl5OWtZDLZ30ZgAwAAACoPoQ1OOZqSLsOwLSewAQAAAJWL0AanHDqeblNGYAMAAAAqH6ENTjl4PM3qdefWUQQ2AAAAoAoQ2uCUg8dSrV43iw71UE8AAACA2oXQBqccKnGmrSmhDQAAAKgShDaUqaCgUIdPWt/TRmgDAAAAqgahDWU6eipD+QXWU0cS2gAAAICqQWhDmUpeGhkRalZYSICHegMAAADULoQ2lKnkzJFMQgIAAABUHUIbynToGPezAQAAAJ5CaEOZSl4e2awhoQ0AAACoKoQ2lKqg0FDyCab7BwAAADyF0IZSHT+dodz8Qqsy7mkDAAAAqg6hDaU6dMz6LFtYSIDC65g91BsAAACg9iG0oVQlZ47k0kgAAACgahHaUCqbSUgIbQAAAECVIrShVCVDG2faAAAAgKpFaINDhYWGDp2wXqONM20AAABA1SK0waGTZ7OUk1tgVdaUNdoAAACAKkVog0MHj6Vava4T5K+6ocwcCQAAAFSlah3acnJy9PLLL6tLly4KDQ1VRESEevbsqQ8//FCFhYVlN1CGP//8U+PHj1fr1q1Vp04dhYaGqlWrVho+fLhef/11paenl91INWbvfjaTyeSh3gAAAAC1U7UNbSkpKerWrZuefvppde/eXT/++KPmzp2rRo0aafz48Ro4cKCys7Mr3P4//vEPde3aVSdPntTLL7+slStXas6cOerTp4++++47PfHEE0pOTnbjO/I+Jaf7b8alkQAAAECV8/N0Byrqhhtu0NatW/Xwww/rrbfespRffvnlGjFihL777juNHz9e06dPL3fbkydP1gsvvKA333xTjzzyiNW2QYMGyd/fX//+979dfAfej5kjAQAAAM+rlmfa5syZoxUrVigwMFCTJ0+22mYymTRlyhRJ0ueff65NmzaVq+0tW7boxRdfVM+ePW0CW5Gnn35aH3zwgRo2bFiR7lcLhmEQ2gAAAAAvUC1D28cffyxJ6t+/vyIiImy2x8fHKz4+XoZh6NNPPy1X26+++qry8/N1++23O6zTvHlz3XvvvXaPXVOknM1WVo71zJFM9w8AAABUvWoX2nJzc7V06VJJUrdu3RzWK9q2cOFCp9vOycnR3LlzJUk9evRwoZfVX8mzbEFmP9UPD/RQbwAAAIDaq9qFth07digvL0/S+TNejhRtO3DggM6dO+dU21u2bFFmZqYkKTY2Vt98840GDRqkqKgohYSEKC4uTnfeeae2bdvm0nuoDmwmIWHmSAAAAMAjql1oO3jwoOV5VFSUw3rFtzk7y+Nff/1leX733XfrrrvuUv/+/fXDDz/op59+0nXXXacvvvhCF110UYUmOKlOSq7RxsyRAAAAgGdUu9kj09L+dwYoMNDx5XrFt6WmpjqsV9zp06ctz2fPnq2VK1fqsssus5T17t1bcXFxGj9+vMaNG6f27dure/fuDtvLyclRTk6O5bVhGJbnRWcLPaXo+I76UTK0NY4M9nif4T3KGj9AaRg/cAXjBxXF2IErKmP8lKetahfaKlNGRobl+cCBA60CW5F77rlHL730kg4cOKAXX3xR3333ncP2pkyZoueee86mPDs7W4sWLXJPp12UkJBgU2YYUtJh67JjB3dq0aKdVdQrVBf2xg/gLMYPXMH4QUUxduAKd46fotuynFHtQlto6P8u0ytt8ezi28LCwpxqOygoyPK8T58+duuYTCb17dtXX3zxhZYuXaqCggL5+vrarTtx4kQ99thjltfx8fE6cuSIAgMDNWTIEKf6VFny8vKUkJCggQMHyt/f32rb6dRsffLLMquyEUMvV4O6QQKk0scPUBbGD1zB+EFFMXbgisoYP85eDShVw9DWrFkzy/OTJ086rFd8W5MmTZxqu169epbn0dHRDuvFxMRIOn9m7vTp0w7vrTObzTKbzZbXxSfy8JYPC39/f5u+HD11xup1YICvGkWGyseHiUhgzd74AZzF+IErGD+oKMYOXOHO8VOedqrdRCTx8fGWN7h//36H9Yq2xcbGKjw83Km2O3ToYHleUFDgsF7xe9NqopIzRzaJJrABAAAAnlLtQltAQIAGDBggSdq4caPDehs2bJAkDR061Om2L7zwQsuC2cVnqSzp8OHzN3yFhYWpfv36TrdfXRw6nm71mkW1AQAAAM+pdqFNku666y5J0tKlS+2uwbZz507t2LFDJpNJY8eOdbrdgIAA3XTTTZa27TEMQytXrpQkDRkyRD4+1fJHWKqSC2s3JbQBAAAAHlMtE8fIkSPVt29fZWdn28zOaBiGJk2aJEkaM2aMLr74YqvtCxYsUFRUlDp06GD38spnn31WERERWrdunRYsWGCz/d///rcOHjyo4OBgPfvss+57U17CMAzWaAMAAAC8SLWbiKTI7Nmz1b9/f7355pvKysrSLbfcotzcXL3//vuaN2+e+vfvrw8++MBmv48++kgpKSlKSUnR3LlzrWZ3lM5PQPLDDz9o6NChGjVqlJ5++mkNGjRI+fn5mj9/vt566y2FhYVp1qxZatu2bVW93SpzNj1HaZnWa0ZweSQAAADgOdXyTJskRUZGasOGDXrppZe0Zs0aXXnllRo+fLiSk5M1depUJSQk2F18e9y4capfv77atWun6667zm7bl156qXbs2KF77rlHX331lS6//HJdeeWVWrx4sR555BH99ddfuuqqqyr7LXpEyUsjA/x9FVU32EO9AQAAAFBtz7RJ56fUnzBhgiZMmOD0PsOGDVNKSkqZ9Ro1aqQ33nhDb7zxhitdrHYOHSsxc2SDOvJl5kgAAADAY6rtmTZUjpLT/XNpJAAAAOBZhDZYKTndPzNHAgAAAJ5FaIMVpvsHAAAAvAuhDRbn0nN0Nj3Hqozp/gEAAADPIrTBouRZNj9fHzWsx8yRAAAAgCcR2mBRMrQ1aVBHvr4MEQAAAMCT+EYOC2aOBAAAALwPoQ0WNpOQcD8bAAAA4HGENlgwcyQAAADgfQhtkCSlZ+bqdGqJmSMJbQAAAIDHEdogyXZRbV8fkxpFhnioNwAAAACKENogSTp4PNXqdeOoOvJj5kgAAADA4/hWDknMHAkAAAB4K0IbJEmHjpUIbcwcCQAAAHgFQhskMXMkAAAA4K0IbVBmdp5SzmVblXF5JAAAAOAdCG2wOcvm42NS4yhmjgQAAAC8AaENNqGtUf0Q+fv5eqg3AAAAAIojtEEHS6zRxiQkAAAAgPcgtEEHj1mv0cYkJAAAAID3ILSBmSMBAAAAL0Zoq+WycvJ14kyWVRkzRwIAAADeg9BWyx0+aX0/m49JimlQx0O9AQAAAFASoa2WO1RiEpLo+iEy+zNzJAAAAOAtCG21XPKJEjNHcmkkAAAA4FUIbbXcoRKhjUlIAAAAAO9CaKvlSp5pI7QBAAAA3oXQVovlF0jHT2dalXF5JAAAAOBdCG212NksyTCsy5owcyQAAADgVQhttdjZDOvXDeoFK9Ds55nOAAAAALCL0FaLnbG+MpJLIwEAAAAvRGirxc4S2gAAAACvR2irxUqGNmaOBAAAALwPoa2WyssvUGqWdVmzhoQ2AAAAwNsQ2mqpIykZKjFxJDNHAgAAAF6I0FZLHTpuvah2ZESQggP9PdQbAAAAAI4Q2mqpkqGNSUgAAAAA70Roq6UOnbAObUxCAgAAAHgnQlstlUxoAwAAAKoFQlstlJdfqKMpGVZlscwcCQAAAHglQlstdDQlXQWF1nNHNuFMGwAAAOCVCG210KyE3Vavg8x+qhPEzJEAAACANyK01TKzEnbplz8PW5Vl5eRrVsIuD/UIAAAAQGkIbbXIrIRdmrF4p91tMxbvJLgBAAAAXojQVkuUFtiKENwAAAAA70NoqyVmlhHYylsPAAAAQNUgtNUSNw9u69Z6AAAAAKoGoa2WGDWwjUaXEchGD26rUQPbVFGPAAAAADiD0FaLlBbcCGwAAACAdyK01TL2ghuBDQAAAPBehLZaaNTANho1sNV/n7cisAEAAABezM/THYBn3DiglerkJGrIgFae7goAAACAUnCmDQAAAAC8GKENAAAAALwYoQ0AAAAAvBihDQAAAAC8GKENAAAAALwYoQ0AAAAAvBihDQAAAAC8GKENAAAAALwYoQ0AAAAAvBihDQAAAAC8GKENAAAAALwYoQ0AAAAAvBihDQAAAAC8mMkwDMPTnagtAgIClJeXJx8fHzVq1MjT3VF2drYCAwM93Q1UU4wfuILxA1cwflBRjB24wt3jxzAMHTlyRJ07d9Yff/xRal0/tx0VZSooKJAkFRYW6vDhwx7uDQAAAABPO3HiRJl1CG1VKDAwUNnZ2fL19VWDBg082peiZN+4cWOZTCaP9gXVD+MHrmD8wBWMH1QUYweuqKzxYxiGU1fgcXlkLZWamqrw8HCdO3dOYWFhnu4OqhnGD1zB+IErGD+oKMYOXOHp8cNEJAAAAADgxQhtAAAAAODFCG21lNls1rPPPiuz2ezprqAaYvzAFYwfuILxg4pi7MAVnh4/3NMGAAAAAF6MM20AAAAA4MUIbQAAAADgxQhtAAAAAODFCG21TE5Ojl5++WV16dJFoaGhioiIUM+ePfXhhx+qsLDQ092DhxUWFur9999XWFiYTCaT9u/f7/S+R44c0cMPP6y4uDgFBgYqOjpaV199tX766afK6zC8Ql5enubMmaPbbrtNbdu2VUhIiAIDA9WsWTONHDlSCxYsKLONc+fO6ZlnnlF8fLyCg4MVGRmp/v37a9asWVXwDuBJOTk5+uGHH/TII4+oZ8+eql+/vvz8/BQaGqpOnTrpkUce0d69e0ttg/GD4q6//nqZTCan/h1j7NReRWOktMcDDzzgcP8qHzsGao2TJ08aHTt2NCQZ48aNM3755Rdj6dKlxogRIwxJRv/+/Y2srCxPdxMesm3bNqNnz56GJMsjKSnJqX3XrFlj1K1b1wgKCjJeeeUVY926dcbs2bONTp06GZKMiRMnVm7n4TGHDh0yYmJiDElGs2bNjHfffddYuXKlsXbtWuO1114z6tevb0gyhg8fbmRnZ9ttIzEx0WjatKnh4+NjTJw40VizZo2xcOFCo2/fvoYkY/To0UZBQUEVvzNUlTvuuMOQZISFhRnPPPOMsWTJEmP9+vXG7NmzjUGDBhmSDLPZbMybN8/u/owfFPfNN984/e8YY6d2k2QEBgYaISEhDh+PPfaY3X09MXYIbbVIv379DEnGww8/bFVeWFhoXHvttYYk4/bbb/dM5+BR//znP42AgACjd+/extNPP12u0HbixAkjMjLSkGTzpers2bNG06ZNDUnGZ599Vjmdh0dt3brVkGQ0adLEOHXqlM32zZs3G35+foYk4/7777fZnp2dbbRu3dqQZLz55ptW23JycoyuXbsakozJkydX1luAh40ZM8aQZKxcudLu9quvvtqQZERERBiZmZlW2xg/KO7kyZNGgwYNjDp16pT57xhjB5KM5cuXl3s/T40dQlstMXv2bMtfFM6cOWOz/a+//jIkGSaTydi4cWPVdxAeFRYWZrz//vtGYWGhMX369HKFtgceeMCQZFxyySV2t0+dOtWQZERHR9t84UL1VxTaXn/9dYd1br31VsvZkrS0NKttr732miHJaNy4sZGfn2+z76JFiwxJRlBQkHH48GG39x+e98wzzxjDhg1zuH3WrFmWz6T169dbbWP8oLhRo0YZ4eHhxosvvljmv2OMHVQ0tHlq7HBPWy3x8ccfS5L69++viIgIm+3x8fGKj4+XYRj69NNPq7h38LS//vpL9913n0wmU7n2y83N1ZdffilJGjlypN06ReXHjx/XDz/84FpH4XUiIyP1+OOP69prr3VY58ILL5R0/t6lXbt2WW0r+mwaPny4fH19bfYdNGiQQkNDlZWVpRkzZrix5/AWL7zwgr7//nuH24svZBsaGmq1jfGDIvPnz9esWbP0+uuvq3HjxmXWZ+ygojw1dghttUBubq6WLl0qSerWrZvDekXbFi5cWCX9gveIiYmp0H6rV6/WuXPnJDkeWw0aNFCzZs0kMbZqooYNG+q1115TXFycwzrF/1GrU6eO5XlSUpJ27twpyfH48fX1VZcuXSQxfmqrr7/+WpJ06aWXqm3btpZyxg+KnDlzRuPHj9fAgQN15513llmfsYOK8uTYIbTVAjt27FBeXp4kqXnz5g7rFW07cOCA5Ys4UJotW7ZYnjsztorXR+2RmJgo6XzAu+CCCyzljB84kp6ertWrV+tvf/ubvvnmG40YMULz5s2zqsP4QZGHH35Y6enpmjZtmlP1GTsosmbNGt10001q3bq16tSpo6ioKPXu3VuvvPKK3e/Cnhw7hLZa4ODBg5bnUVFRDusV35acnFypfULNUN6xdejQoUrvE7xLfn6+5syZI0l6/PHHrc66lXf8nDlzRhkZGZXUU3iDvXv3ytfXV6Ghoerdu7d+//13zZ49W3PnzrUZI4wfSOfPZHz55ZeaMmWKYmNjndqHsYMizz77rOrXr6/3339fq1at0r///W8FBwdrwoQJ6tChg/744w+r+p4cO4S2WiAtLc3yPDAw0GG94ttSU1MrtU+oGco7thhXtc8nn3yi48ePq3v37nr44YettvHZhJKaNm2qzZs3a/369fryyy/VqFEjXX/99Ro4cKDVlyWJ8YPz62Tdc889uuyyy3T//fc7vR9jB5LUr18//fjjj3rvvfc0cOBAXXTRRbruuuv0008/6ZZbblFycrKuuuoqnTx50rKPJ8cOoQ0AUCl2796tJ598Ug0aNNCsWbPk7+/v6S7BywUEBKhDhw7q1q2bbrnlFq1cuVJ33nmnfv75Z/Xu3dvqyxPw+OOP69SpU/r444/LPZEWsHz5cg0YMMCm3GQy6c0331RAQICOHz+u119/3QO9s0VoqwWKz7aVnZ3tsF7xbWFhYZXaJ9QM5R1bjKva4/jx4xo6dKj8/Py0ZMkStWjRwqYOn00oi8lk0htvvKGQkBAdOnRIL7zwgmUb46d2W7JkiT755BM9//zzat26dbn2ZeygLJGRkerataskWc187cmxQ2irBYpm7pNU6l8pi29r0qRJpfYJNUN5x1bTpk0rvU/wvGPHjql///46deqUfvrpJ8uU/yWVd/zUrVtXISEh7u0svF5YWJh69OghSVZLAzB+aq+0tDTdfffd6tatmx577LFy78/YgTOKxklSUpJNmVT1Y8fPLa3Aq8XHx8vf3195eXnav3+/w3pF22JjYxUeHl41nUO11qlTJ8vz/fv3O5xJqWhsFa+Pmik5OVkDBgzQmTNntHz5coeBTbIdP44wfhAdHS1JOnz4sKWM8VN7bdq0SQcPHlRycrLVOn5FDMOwPC8+Y+2YMWP0ySefMHbglOLjqIgnxw5n2mqBgIAAyzW7GzdudFhvw4YNkqShQ4dWSb9Q/fXq1csS8B2NrRMnTlgmEGBs1Wz79+9Xnz59lJaWphUrVtgEtv379ys9Pd3yukWLFpZ1txyNn4KCAsvsXYyfmufw4cNq27atfv3111LrFU29XfwPioyf2qtbt27aunWrNm/erD///NPm8fzzz1vqLlq0yKacsYNx48Zp+vTppdYp+u5S/A/Snhw7hLZa4q677pIkLV261O66Ezt37tSOHTtkMpk0duzYqu4eqimz2axbb71VkizTupc0d+5cSef/Un711VdXWd9QtRITE9WnTx/l5+dr1apVateunU2dFi1aaPbs2VZlRZ9N8+fPV2Fhoc0+CQkJSktLU2BgoG6++ebK6Tw8Ji8vT7t27dLatWsd1snKytKaNWskST179rTaxvipnUJCQtShQweHj5iYGEvd1q1b2y1n7NRuS5Yscfi9RTr/B+eiUFYyeHls7BioNfr27WtIMh599FGr8sLCQmPEiBGGJOP222/3UO/gLaZPn25IMiQZSUlJZdY/ceKEERkZaUgyvvvuO6tt586dM2JjYw1JxmeffVZJPYanbd++3WjUqJHRsmVLY//+/Q7rSTKmT59uVZadnW20bt3akGS8/fbbVttyc3ONbt26GZKMyZMnV0bX4WFJSUmGJKNJkybGsWPH7NZ5+OGHDUmGyWQyVq1aZbWN8QN7nPl3jLFTu8XGxhq+vr7G6tWrbbYVFhYaN910kyHJiIyMtPls8tTYIbTVIidPnjQ6duxoSDLuvfde49dffzWWLVtmjBw50pBk9O/f38jKyvJ0N+EBx48fN7Zu3Wps3brVeOGFFyz/2P3000+W8vT0dIf7r1mzxqhbt64RHBxsvPrqq8b69euNuXPnGhdeeKEhyZg4cWIVvhtUpT179hhRUVGGJCMgIMAICQlx+LAX2gzDMBITE42mTZsavr6+xjPPPGOsWbPGWLRokdGvXz9DkjF69GijoKCg6t8cKt3hw4cNs9lsSDLq1atnPP/888aPP/5o/P7778bcuXONq666ypBkmM1m4+OPP7bbBuMHhmEY6enpZf47VhJjp/Yq+j4cFBRkTJgwwVi0aJGxadMmY/bs2Ub//v0NSUbjxo2NdevW2d3fE2OH0FbLZGdnGy+99JJx4YUXGiEhIUZYWJhxySWXGFOnTuWDqRZ79tlnLf/AOXosX7681DYOHz5sPPjgg0bLli0Ns9lsREVFGUOHDjUWL15cNW8CHjFv3rwyx07xh73QZhiGcfbsWWPSpElG27ZtjcDAQKNevXrG5Zdfbnz99ddV+4ZQ5VJSUowPP/zQGDlypNGqVSsjODjY8PX1NSIiIoyuXbsaTz75pLFnz55S22D8YPny5WV+/tjD2KmdcnJyjLlz5xp33nmn0bFjRyM0NNTw8/Mz6tWrZ1x66aXGSy+9ZJw+fbrUNqp67JgMw87UKAAAAAAAr8BEJAAAAADgxQhtAAAAAODFCG0AAAAA4MUIbQAAAADgxQhtAAAAAODFCG0AAAAA4MUIbQAAAADgxQhtAAAAAODFCG0AAAAA4MUIbQAAAADgxQhtAACv0bt3b5lMJplMJv3666+W8o0bN1rKe/fu7dIxitpx9Jg0aZLTbc2ZM6fM9lasWOFSfyti27ZtDvtz++23V3l/AACu8fN0BwAAkKT8/Hz9/vvvkiQ/Pz9dfPHFlm1r1661PO/Ro4dLx3n88cclSdu3b9fixYtttk+bNk3//Oc/FRgYWGZb7777rt3y9u3ba/DgwZKkpk2butDbiomMjNTjjz+u1NRUTZs2rcqPDwBwL0IbAMArbNmyRVlZWZKkTp06KSgoyLKteGi75JJLXDrOa6+9Jkn67LPP7Ia2lJQUzZw5U2PHji21na1bt2rlypV2t3Xt2tVyHE9o2LChXnvtNe3fv5/QBgA1AJdHAgC8QmnBzJ2hraRWrVrZlDk6g2avTuvWrd3aHwAASiK0AQC8wrp16yzPi18CmZKSor1790qSGjVqpGbNmrn1uCNHjlRUVJRV2Z9//qlVq1Y53OfMmTOaMWOGOnfurJ49e7q1PwAAlERoAwB4heKhrfjZtMo8yyZJZrNZd911l035O++843CfTz75RJmZmbr//vvd3h8AAEoitAEAqtzkyZNtZjXctWuXZXvbtm0t5cOGDbOUz58/32qfzz77zC39ueeee+TjY/1P4vz583Xo0CGbuoWFhZo6daoiIiJ08803V/iYSUlJeuaZZ9S7d29FR0crICBAYWFhuuCCC3TTTTfpP//5j/Lz88tsZ+XKlRo1apSaNm0qs9msunXrqmfPnnr//fed2r+k5cuXa+zYsYqPj1d4eLgCAgLUsGFDDRgwQK+//rrS0tJK3X/nzp168skndckllygyMlIBAQEKCQlR8+bN1b17d91666164403tHHjxnL3DQBqKyYiAQDUerGxsbr66qv1/fffW8oKCgr0/vvv66WXXrKq+8MPPygpKUmPPPKIgoODy32svLw8/f3vf9frr7+ugoICSVJcXJyuv/56HT9+XN9//7327t2rWbNmqXXr1po5c6bVTJpFDMPQo48+qrfffttSFhAQoEGDBqlu3bp68cUX9dVXXzndr5SUFN1yyy366aefLGXXXnutoqOjNWfOHC1btkzLli3Tyy+/rG+//VZ9+/a1aWPq1Kl66KGHLO+rTp06GjZsmBo3bqzTp09rzZo12rBhg7766ivFxMQoOTnZ6f4BQG1GaAMAVLng4GDVr1/f8jotLU25ubmSpNDQUAUEBEg6H5zOnj0rSfL19VVERIRVO85My++s+++/3yq0SdLHH3+sZ5991momy3fffVcmk0n33XdfuY9hGIZuueUWffPNN5ayzp0767fffrMcY/bs2brhhhskSbt371afPn20cuVKde3a1aqtKVOmWAU2Sfrmm2907bXXSpL+9a9/qVOnTk71Kz09XQMGDNCWLVssZRMmTLAE1kceeUQdO3ZUQUGBTp48qSuvvFJr165V586dLfW3b99uFdhCQ0O1bds2q3sQCwoK9Mwzz+jll192ql8AgPO4PBIAUOWeeuoppaSkWB7FZ2BcvXq1pfw///mPpbx///5W+6SkpGjUqFFu69PAgQNtZpI8deqUZsyYYXm9Y8cO/fzzz3brOuPTTz+1CmyS9OCDD1qFwuuvv17Nmze3vM7MzNSoUaOUl5dnKTt+/Liee+45q3bi4uIsgU2SoqKidPfddzvVr2effdYqsEnSuHHjLM/j4+OtJofJycmxuQ9w9uzZlsAmnZ9Vs+SkMb6+vpoyZYrV+wMAlI3QBgDwqLNnz2r79u2SpIiICHXo0MGybfXq1ZbnvXr1qtR+mEwmjR8/3qa8+PT/7733niRVeAKSKVOm2JTZm1ylW7duVq/37t1rFfa+/PJLy5nJIvYuoXTmZ5aRkaEPPvjAqiw4OFgtW7a0Kmvfvr3V602bNmnNmjWW10ePHrXa/ueff2r27NkyDMOq3GQy6a233tI//vGPMvsGADiP0AYA8KjVq1dbvtj36tVLJpPJaluRyg5tknTHHXfY3Ke2ZcsWrVixQqmpqfriiy8s97+V1+7duy1LFxTXpEkTm7KYmBibsh9//NHy/Ndff3VqH3ttl/TLL79YFjUvUvzS1SL16tWzKVu+fLnlecmzagUFBbrhhhsUGxursWPH6osvvtDx48clnb9X7p577imzbwCA87inDQDgUcWDWe/evS3PCwoKLMsA+Pj4WF2eV1kiIiJ000036ZNPPrEqf+edd9SnTx+lp6frmWeesZlp0hm7d++2Wx4aGmpTVqdOnVL3P3DggFPtODNRSmJiok1ZamqqnnjiCauy3377zabetm3bLM9HjBihyZMnW13GKUmHDh3S9OnTNX36dJlMJvXt21dPPPGEhg4dWmbfAADnEdoAAFVm//79atGihcPtkyZN0qRJk2zKCwsLFR4ebnndt29frVixojK6qPvvv98mtH3//fdav369zGaz7rzzzgq1e+7cOZsyHx8fuwHQ39/fpqxoQhZJdqfd9/Oz/SfdmXCZmppqU3bu3Dm9/vrrZe575swZy/P4+Hi99957euCBB2yCWxHDMLRixQqtWLFCDz30kM1EKgAA+7g8EgCAYrp06aKePXtalRUUFOjw4cO68cYbFRUVVaF2i4fOIoWFhSosLLQptxd6is+cae+smr012ey1XVJYWJhNWfPmzWUYRpmP4pdsSucnL9mxY4cee+wxXXDBBaUe95133rG6vBIA4BihDQBQZXx8fBQeHm55hISEWLb5+flZbSt+5qhOnTpW2+xdPuhOjqbzr+gEJJLUpk0bu+X2zprZKyu+v73ZF+3tk5mZWWa/7J35tHf2zVlxcXF6/fXXlZiYqOTkZH311Ve67rrr7J71Kxn6AAD2EdoAAFWmWbNmOnv2rOXxr3/9y7LtgQcesNrWuHFjy7Y9e/ZYbfvhhx8qtZ833HCDzRm1iy++2O5Mj85q1aqV3WUCDh06ZFNmb9HpIUOGWJ4Xv/evtH2cWbz6sssus7kc8/Tp0zp16pTd+oZhqG/fvurcubOWLFliKf/3v/+tHj16WB0zJiZGo0eP1pw5czRr1iybtlwJhwBQmxDaAAAeU/y+tMsvv9zyfO/evTp48KAkqV27doqOjq7SfpnNZpt1yCqymHZJ9u7XW79+vU3Zhg0brF63atXKsuC2JN1yyy0ym81WdTZu3GjTjr3JQ0oKDw/XbbfdZlO+cOFCu/WXLVumVatWKSkpySrEHjp0SOvWrdPcuXPt7nfVVVfZlLFeGwA4h9AGAPCIwsJCrVq1StL5yyb79Olj2Vb8Xqf+/ftXed8k6V//+pfS0tIsjzvuuMPlNseMGaObbrrJquzdd9+1mnJ/1qxZlsAqSSEhIZo1a5bV5aLR0dF69tlnrdpJSkrSnDlzLK9PnjypadOmOdWvl156SbGxsVZl//jHP2zOAiYmJlomYpk8ebLd+/QmTZpk9161kmfaAgMD9be//c2p/gFAbcfskQAAj/jjjz8sMyJ26dLFaqKNZcuWWZ4XPwPnDkVT2Rct6L1kyRKlp6dLkl577TVLPV9fX4f3zp05c0YvvviiJPtnuDZu3Gg5zvjx4xUXFyfp/MLSRWu9vfrqqyooKNCff/6pTp06adCgQTpx4oTmz59vaadNmzb6+uuv1aVLF5tjPP300zp27JjeeecdS9lNN92kESNGqF69evr+++/VvHlzHTt2zG7f4uLiLIuJR0ZGasWKFbrhhhss7+fgwYNq166dhgwZogYNGmjfvn1KSEhQfn6+nn76aT366KN2fzYZGRnq37+/evToofbt28vPz09bt261OusXHByszz77rNSZRAEA/2MyilY0BQCgCr322mt68sknJUmPP/64VWBq3Lixjh49KpPJpJSUFLsLO1dU8cW7S3L2n8Syli4obvny5erXr5/dNqZNm6aVK1dq9+7dOnv2rMxms6Kjo9WtWzeNGDFCI0eOlK+vb6ntr1y5UlOnTtXq1at14sQJBQUFqW3btrr55ps1ZMgQtW7d2u5+9pZNKCz8//bt0EZAAIii4CoUEoFH4OgAQydQHhJPAwhCDyQkGAQGd7kCzpxcMVPAF+ue2J9YliXmeY5t2+K+7/i+L8qyjKZpou/7mKYpuq77s3eeZ6zrGvu+x3EccV1XPM8T7/tGURRRVVW0bRvDMMQ4jlHX9b/uB4BoAwAASM1PGwAAQGKiDQAAIDHRBgAAkJhoAwAASEy0AQAAJCbaAAAAEhNtAAAAiYk2AACAxEQbAABAYqINAAAgMdEGAACQmGgDAABITLQBAAAkJtoAAAAS+wWFNF490/lJHAAAAABJRU5ErkJggg==",
            "text/plain": [
              "<Figure size 1000x600 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "U, S, VT = svd(variable, full_matrices=False)\n",
        "\n",
        "energy_retained = np.cumsum(S) / np.sum(S)\n",
        "\n",
        "\n",
        "plt.rc('font', family='DeJavu Serif', size=18)\n",
        "# plt.rcParams['font.family'] = 'Times New Roman'\n",
        "colors = sns.color_palette('deep', n_colors=100)\n",
        "\n",
        "plt.figure(figsize=(10, 6))\n",
        "plt.plot(energy_retained[:50], color=colors[0], linewidth=3, marker='D', markersize=5, markevery=1)\n",
        "\n",
        "ax = plt.gca()\n",
        "ax.spines['bottom'].set_linewidth(2)  # X-axis line width\n",
        "ax.spines['left'].set_linewidth(2)    # Y-axis line width\n",
        "\n",
        "plt.xlabel('# Modes', fontsize=20, weight='bold')\n",
        "plt.ylabel('Energy', fontsize=20, weight='bold')\n",
        "\n",
        "#plt.legend(frameon=True, edgecolor='black')\n",
        "# handles, labels = plt.gca().get_legend_handles_labels()\n",
        "# plt.legend(handles=handles, labels=labels, frameon=True, edgecolor='black', loc='upper center', bbox_to_anchor=(0.5, -0.15), fancybox=True, shadow=True, ncol=2)\n",
        "\n",
        "plt.grid(True)\n",
        "plt.axhline(y=0.99, linewidth=3, color=colors[1])\n",
        "plt.title('Accumulated Energy', fontsize=24, weight='bold')\n",
        "# plt.savefig('fig.jpg', dpi=400, bbox_inches='tight')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gaDApmlSf34l"
      },
      "source": [
        "### Reduction"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ctXSPLzRf34l",
        "outputId": "9ba826ba-a0ac-4fb5-d228-de0debb1ff32"
      },
      "outputs": [],
      "source": [
        "num_modes = 45  # correpsond to 99% energy\n",
        "\n",
        "# Energy content calculation\n",
        "total_energy = np.sum(S)\n",
        "retained_energy = np.sum(S[:num_modes])\n",
        "energy_content = retained_energy / total_energy * 100\n",
        "print(\"With\", num_modes, \"number of modes, the energy content is:\", energy_content)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 28,
      "metadata": {
        "id": "Z6VOqaEjf34l"
      },
      "outputs": [],
      "source": [
        "q = U[:,:num_modes].T @ variable\n",
        "\n",
        "variable_reconst = U[:,:num_modes] @ q\n",
        "\n",
        "RelErr1 = np.zeros((q.shape[1]))\n",
        "for t in range(q.shape[1]):\n",
        "        # print(t)\n",
        "        RelErr1[t] = (norm(variable[:,t] - variable_reconst[:,t])/norm(variable[:,t]) ) * 100"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "metadata": {
        "id": "2oQSlh77f34m"
      },
      "outputs": [],
      "source": [
        "np.savez(f'/content/drive/MyDrive/Paper/RoSo/{model}/{bc_name}/POD_data_{model}_{bc_name}_{variable_name}.npz', S=S, Reconst_error=RelErr1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 54
        },
        "id": "ZwCHnXRBf34m",
        "outputId": "e4515278-cb0d-4660-dc60-a4c56be2c3d4"
      },
      "outputs": [],
      "source": [
        "plt.rc('font', family='DeJavu Serif', size=18)\n",
        "# plt.rcParams['font.family'] = 'Times New Roman'\n",
        "\n",
        "plt.figure(figsize=(10, 6))\n",
        "plt.plot(RelErr1, color=colors[3], linewidth=3)\n",
        "\n",
        "ax = plt.gca()\n",
        "ax.spines['bottom'].set_linewidth(2)  # X-axis line width\n",
        "ax.spines['left'].set_linewidth(2)    # Y-axis line width\n",
        "\n",
        "plt.xlabel('Snapshot', fontsize=20, weight='bold')\n",
        "plt.ylabel('L2 Norm Difference %', fontsize=20, weight='bold')\n",
        "\n",
        "plt.yscale('log')  # Set y-axis to logarithmic scale\n",
        "plt.xlim(0, 1000)  # X-axis range from 0 to 1\n",
        "plt.ylim(0.01, 100)  # Y-axis range from 0.01 to 100\n",
        "\n",
        "plt.yticks([0.01, 0.1, 1, 10, 100], fontsize=16, weight='normal')\n",
        "\n",
        "\n",
        "#plt.legend(frameon=True, edgecolor='black')\n",
        "# handles, labels = plt.gca().get_legend_handles_labels()\n",
        "# plt.legend(handles=handles, labels=labels, frameon=True, edgecolor='black', loc='upper center', bbox_to_anchor=(0.5, -0.15), fancybox=True, shadow=True, ncol=2)\n",
        "\n",
        "plt.grid(True)\n",
        "\n",
        "\n",
        "plt.title('Reconstruction Error', fontsize=24, weight='bold')\n",
        "\n",
        "# plt.savefig('fig.jpg', dpi=400, bbox_inches='tight')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "cktgziBff34m"
      },
      "source": [
        "# LSTM"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-8HDoX2bf34n"
      },
      "source": [
        "### Generating Input/Output Matrices"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 31,
      "metadata": {
        "id": "j60BSVHGf34n"
      },
      "outputs": [],
      "source": [
        "t = np.linspace(1, n_snapshots, num=n_snapshots)\n",
        "inputt = np.concatenate((t.reshape(-1, 1), interpolated_v_values.reshape(-1, 1), interpolated_P_values.reshape(-1, 1), q.T), axis=1)\n",
        "outputt = np.copy(q).T"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mwk8XwbIf34n"
      },
      "source": [
        "### Scaling Input/Output Matrices"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 32,
      "metadata": {
        "id": "kZrcP3cif34n"
      },
      "outputs": [],
      "source": [
        "min_vals_input = np.min(inputt, axis=0)\n",
        "max_vals_input = np.max(inputt, axis=0)\n",
        "input_scaled = inputt - min_vals_input[np.newaxis, :]\n",
        "input_scaled /= (max_vals_input - min_vals_input)[np.newaxis, :]\n",
        "\n",
        "min_vals_output = np.min(outputt, axis=0)\n",
        "max_vals_output = np.max(outputt, axis=0)\n",
        "output_scaled = outputt - min_vals_output[np.newaxis, :]\n",
        "output_scaled /= (max_vals_output - min_vals_output)[np.newaxis, :]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "SKPpgg0Of34o"
      },
      "source": [
        "### Time Windowing (Sequencing)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 33,
      "metadata": {
        "id": "hPKPvo_Kf34o"
      },
      "outputs": [],
      "source": [
        "time_window = 100\n",
        "\n",
        "x_train = []\n",
        "y_train = []\n",
        "\n",
        "for i in range(0,len(input_scaled) - time_window -1):\n",
        "    x_train.append( input_scaled[i : (i+time_window) , :] )\n",
        "    y_train.append( output_scaled[i+time_window,:])\n",
        "\n",
        "total_x = np.array(x_train)\n",
        "total_y = np.array(y_train)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "pJiPrTGPf34o"
      },
      "source": [
        "### Train/Test Split"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 34,
      "metadata": {
        "id": "shK2UYaof34o"
      },
      "outputs": [],
      "source": [
        "Test_split = 0.8\n",
        "\n",
        "x_train = total_x[:int(Test_split * total_x.shape[0]),:,:]\n",
        "y_train = total_y[:int(Test_split * total_y.shape[0])]\n",
        "\n",
        "x_test = total_x[int(Test_split * total_x.shape[0]):, :,:]\n",
        "y_test = total_y[int(Test_split * total_y.shape[0]):]\n",
        "\n",
        "# Convert data to PyTorch tensors\n",
        "x_train = torch.tensor(x_train, dtype=torch.float32)\n",
        "y_train = torch.tensor(y_train, dtype=torch.float32)\n",
        "x_test = torch.tensor(x_test, dtype=torch.float32)\n",
        "y_test = torch.tensor(y_test, dtype=torch.float32)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "d27e7H1Af34o"
      },
      "source": [
        "### LSTM Network"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 35,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "82CyZz2mf34p",
        "outputId": "ac10ce22-838d-43a4-ce19-fd3008813ec2"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Using device: cpu\n"
          ]
        }
      ],
      "source": [
        "# Device setup\n",
        "device = torch.device('cpu') if state_lstm=='load' else torch.device('cuda' if torch.cuda.is_available() else 'cpu')\n",
        "print(f'Using device: {device}')\n",
        "\n",
        "# Hyperparameters\n",
        "input_size = num_modes + 3  # 3 = (t & v & P)\n",
        "hidden_size = 512\n",
        "output_size = num_modes\n",
        "sequence_length = time_window\n",
        "learning_rate = 0.00001\n",
        "num_epochs = 200000\n",
        "batch_size = 20\n",
        "early_stop_patience = 100\n",
        "dropout_prob = 0.2\n",
        "\n",
        "# LSTMNet definition\n",
        "class LSTMNet(nn.Module):\n",
        "    def __init__(self, input_size, hidden_size, output_size, dropout_prob):\n",
        "        super(LSTMNet, self).__init__()\n",
        "        self.hidden_size = hidden_size\n",
        "        self.lstm1 = nn.LSTM(input_size, hidden_size, num_layers=3, batch_first=True, dropout=dropout_prob)\n",
        "        self.fc1 = nn.Linear(hidden_size, output_size)\n",
        "\n",
        "        # Initialize LSTM weights\n",
        "        for name, param in self.lstm1.named_parameters():\n",
        "            if 'weight' in name:\n",
        "                init.xavier_normal_(param)\n",
        "            elif 'bias' in name:\n",
        "                init.constant_(param, 0.0)\n",
        "\n",
        "        # Initialize linear layer weights\n",
        "        init.xavier_normal_(self.fc1.weight)\n",
        "        init.constant_(self.fc1.bias, 0.0)\n",
        "\n",
        "    def forward(self, x):\n",
        "        # Initialize hidden and cell states on the correct device\n",
        "        h01 = torch.zeros(3, x.size(0), self.hidden_size, device=device).requires_grad_()\n",
        "        c01 = torch.zeros(3, x.size(0), self.hidden_size, device=device).requires_grad_()\n",
        "        out1, (h01, c01) = self.lstm1(x, (h01.detach(), c01.detach()))\n",
        "        out = nn.functional.relu(self.fc1(out1[:, -1, :]))\n",
        "        return out\n",
        "\n",
        "model_lstm = LSTMNet(input_size, hidden_size, output_size, dropout_prob).to(device)\n",
        "\n",
        "criterion = nn.MSELoss()\n",
        "optimizer = torch.optim.Adam(model_lstm.parameters(), lr=learning_rate)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "03St5AYaf34p"
      },
      "source": [
        "### Train"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 36,
      "metadata": {
        "id": "2ffGyBro0-IJ"
      },
      "outputs": [],
      "source": [
        "# Move the Data to GPU\n",
        "x_train = x_train.to(device)\n",
        "y_train = y_train.to(device)\n",
        "x_test = x_test.to(device)\n",
        "y_test = y_test.to(device)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 37,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "WcuV_oGVf34p",
        "outputId": "1d06779c-1091-40cc-fde5-883f6f461185"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Model is loaded completely!\n"
          ]
        },
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "<ipython-input-37-2ed09f8d722d>:53: FutureWarning: You are using `torch.load` with `weights_only=False` (the current default value), which uses the default pickle module implicitly. It is possible to construct malicious pickle data which will execute arbitrary code during unpickling (See https://github.com/pytorch/pytorch/blob/main/SECURITY.md#untrusted-models for more details). In a future release, the default value for `weights_only` will be flipped to `True`. This limits the functions that could be executed during unpickling. Arbitrary objects will no longer be allowed to be loaded via this mode unless they are explicitly allowlisted by the user via `torch.serialization.add_safe_globals`. We recommend you start setting `weights_only=True` for any use case where you don't have full control of the loaded file. Please open an issue on GitHub for any issues related to this experimental feature.\n",
            "  state_dict = torch.load(f'/content/drive/MyDrive/Paper/RoSo/{model}/{bc_name}/{name}', map_location=torch.device('cpu'))\n"
          ]
        }
      ],
      "source": [
        "name = f'LSTMNet_{model}_{bc_name}_{variable_name}.pt'\n",
        "\n",
        "best_loss = float('inf')\n",
        "early_stop_count = 0\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "if state_lstm == 'train':\n",
        "    for epoch in range(num_epochs):\n",
        "        model_lstm.train()\n",
        "        for i in range(0, x_train.shape[0], batch_size):\n",
        "            # Move the batch data to GPU\n",
        "            batch_x = x_train[i:i+batch_size].to(device)\n",
        "            batch_y = y_train[i:i+batch_size].to(device)\n",
        "\n",
        "            # Forward pass\n",
        "            outputs = model_lstm(batch_x)\n",
        "\n",
        "            # Compute loss and backpropagation\n",
        "            loss = criterion(outputs, batch_y)\n",
        "            optimizer.zero_grad()\n",
        "            loss.backward()\n",
        "            optimizer.step()\n",
        "\n",
        "        # Evaluate the model on training and test data\n",
        "        with torch.no_grad():\n",
        "            model_lstm.eval()\n",
        "            train_outputs = model_lstm(x_train)\n",
        "            train_loss = criterion(train_outputs, y_train)\n",
        "            test_outputs = model_lstm(x_test)\n",
        "            test_loss = criterion(test_outputs, y_test)\n",
        "\n",
        "        # Print training and test loss every epoch\n",
        "        print('Epoch [{}/{}], Train Loss: {:.6f}, Test Loss: {:.6f}'.format(epoch+1, num_epochs, train_loss.item(), test_loss.item()))\n",
        "\n",
        "        # Check if the current test loss is the best so far\n",
        "        if test_loss < best_loss:\n",
        "            best_loss = test_loss\n",
        "            early_stop_count = 0\n",
        "\n",
        "            torch.save(model_lstm.state_dict(), name)\n",
        "            shutil.copy(name, f'/content/drive/MyDrive/Paper/RoSo/{model}/{bc_name}/')\n",
        "\n",
        "        else:\n",
        "            early_stop_count += 1\n",
        "            if early_stop_count == early_stop_patience:\n",
        "                print('Early stopping after {} epochs'.format(epoch+1))\n",
        "                break\n",
        "\n",
        "\n",
        "elif state_lstm == 'load':\n",
        "    state_dict = torch.load(f'/content/drive/MyDrive/Paper/RoSo/{model}/{bc_name}/{name}', map_location=torch.device('cpu'))\n",
        "    model_lstm.load_state_dict(state_dict)\n",
        "    model_lstm = model_lstm.to(torch.device('cpu'))\n",
        "\n",
        "    # Recheck the model is loaded completely!\n",
        "    missing_keys, unexpected_keys = model_lstm.load_state_dict(state_dict, strict=False)\n",
        "    if missing_keys or unexpected_keys:\n",
        "        print(\"Missing keys:\", missing_keys)\n",
        "        print(\"Unexpected keys:\", unexpected_keys)\n",
        "    else:\n",
        "        print(\"Model is loaded completely!\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "S8MKhkfaf34p"
      },
      "source": [
        "### Inference"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 38,
      "metadata": {
        "id": "o4t4mH-2f34p"
      },
      "outputs": [],
      "source": [
        "model_lstm.eval()\n",
        "with torch.no_grad():\n",
        "    x = x_test[0:1,:,:].cpu()\n",
        "    predict = model_lstm(x)\n",
        "    pred_final = []\n",
        "\n",
        "    for i in range (1,x_test.shape[0]):\n",
        "        pred_final.append(predict)\n",
        "        x = x_test[i:i+1, :, :]\n",
        "        x = x.to(device)\n",
        "        predict = model_lstm(x)\n",
        "\n",
        "pred_final = torch.cat(pred_final, dim=0).cpu().detach().numpy()"
      ]
    }
  ],
  "metadata": {
    "accelerator": "GPU",
    "colab": {
      "collapsed_sections": [
        "yNiFcBU9f34Q",
        "z6k9P75cf34a",
        "ob-RwiNyf34d",
        "qyURirjMf34e",
        "q7hVSxV9f34f",
        "gaDApmlSf34l",
        "-8HDoX2bf34n",
        "mwk8XwbIf34n",
        "SKPpgg0Of34o",
        "pJiPrTGPf34o",
        "d27e7H1Af34o",
        "03St5AYaf34p"
      ],
      "gpuType": "A100",
      "machine_shape": "hm",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "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.11.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}