hurricane_prediction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Time</th>\n",
       "      <th>Longitude</th>\n",
       "      <th>Latitude</th>\n",
       "      <th>Velocity component u</th>\n",
       "      <th>Velocity component v</th>\n",
       "      <th>Geopotential</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1800</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>-1.156552</td>\n",
       "      <td>0.424945</td>\n",
       "      <td>14506.583520</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1770</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.016667</td>\n",
       "      <td>-0.824795</td>\n",
       "      <td>0.036197</td>\n",
       "      <td>14507.684626</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1740</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.033333</td>\n",
       "      <td>-0.958053</td>\n",
       "      <td>-0.408673</td>\n",
       "      <td>14517.080724</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1710</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.050000</td>\n",
       "      <td>-1.131566</td>\n",
       "      <td>-0.248520</td>\n",
       "      <td>14533.670709</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1680</th>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.066667</td>\n",
       "      <td>-0.588818</td>\n",
       "      <td>0.034828</td>\n",
       "      <td>14539.690085</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>159359</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.933333</td>\n",
       "      <td>2.721812</td>\n",
       "      <td>11.960073</td>\n",
       "      <td>14993.198635</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>159329</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.950000</td>\n",
       "      <td>3.345071</td>\n",
       "      <td>13.250880</td>\n",
       "      <td>14964.349677</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>159299</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.966667</td>\n",
       "      <td>4.263996</td>\n",
       "      <td>14.024268</td>\n",
       "      <td>14942.841421</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>159269</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.983333</td>\n",
       "      <td>5.542441</td>\n",
       "      <td>13.880541</td>\n",
       "      <td>14926.838691</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>159239</th>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>7.033265</td>\n",
       "      <td>14.014687</td>\n",
       "      <td>14889.474520</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>161040 rows × 6 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "        Time  Longitude  Latitude  Velocity component u  Velocity component v  \\\n",
       "1800     0.0        0.0  0.000000             -1.156552              0.424945   \n",
       "1770     0.0        0.0  0.016667             -0.824795              0.036197   \n",
       "1740     0.0        0.0  0.033333             -0.958053             -0.408673   \n",
       "1710     0.0        0.0  0.050000             -1.131566             -0.248520   \n",
       "1680     0.0        0.0  0.066667             -0.588818              0.034828   \n",
       "...      ...        ...       ...                   ...                   ...   \n",
       "159359   1.0        1.0  0.933333              2.721812             11.960073   \n",
       "159329   1.0        1.0  0.950000              3.345071             13.250880   \n",
       "159299   1.0        1.0  0.966667              4.263996             14.024268   \n",
       "159269   1.0        1.0  0.983333              5.542441             13.880541   \n",
       "159239   1.0        1.0  1.000000              7.033265             14.014687   \n",
       "\n",
       "        Geopotential  \n",
       "1800    14506.583520  \n",
       "1770    14507.684626  \n",
       "1740    14517.080724  \n",
       "1710    14533.670709  \n",
       "1680    14539.690085  \n",
       "...              ...  \n",
       "159359  14993.198635  \n",
       "159329  14964.349677  \n",
       "159299  14942.841421  \n",
       "159269  14926.838691  \n",
       "159239  14889.474520  \n",
       "\n",
       "[161040 rows x 6 columns]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "from pennylane import numpy as np\n",
    "data = pd.read_csv('dataFranklin.csv')              #Read the file\n",
    "data.drop(columns='Unnamed: 0', inplace = True)     # Drup unwanted column\n",
    "data.sort_values(by = ['Time', 'Longitude', 'Latitude'], ascending=[True, True, True], inplace=True)    #Reorder the rows\n",
    "\n",
    "radius_earth = 6371000 # radius earth in meters\n",
    "\n",
    "# Transform longitud and laititude to meters (x and y)\n",
    "data['Longitude'] = radius_earth*data['Longitude']\n",
    "data['Latitude'] = radius_earth*data['Latitude']\n",
    "\n",
    "\n",
    "# Function to minmax scale and save linear transformation coefficients\n",
    "def minmaxscale(data, min_val, max_val):\n",
    "    min_dat = data.min()\n",
    "    max_dat = data.max()\n",
    "    \n",
    "    alpha = (min_val - max_val)/(min_dat-max_dat)\n",
    "    beta = min_val - alpha*min_dat\n",
    "    rescaled_data = alpha*data + beta\n",
    "    return rescaled_data, alpha, beta\n",
    "\n",
    "# We now which time in the future we want to predict so we can scale this from the begginning\n",
    "data['Time'], alpha_t, beta_t = minmaxscale(data['Time'], 0, 1)\n",
    "data['Longitude'], alpha_x, beta_x = minmaxscale(data['Longitude'], 0, 1)\n",
    "data['Latitude'], alpha_y, beta_y = minmaxscale(data['Latitude'], 0, 1)\n",
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.00000000e+00,  0.00000000e+00,  0.00000000e+00,\n",
       "        -1.15655232e+00,  4.24944533e-01,  1.45065835e+04],\n",
       "       [ 0.00000000e+00,  0.00000000e+00,  1.66666667e-02,\n",
       "        -8.24795318e-01,  3.61968769e-02,  1.45076846e+04],\n",
       "       [ 0.00000000e+00,  0.00000000e+00,  3.33333333e-02,\n",
       "        -9.58053360e-01, -4.08672800e-01,  1.45170807e+04],\n",
       "       ...,\n",
       "       [ 1.00000000e+00,  1.00000000e+00,  9.66666667e-01,\n",
       "         4.26399619e+00,  1.40242684e+01,  1.49428414e+04],\n",
       "       [ 1.00000000e+00,  1.00000000e+00,  9.83333333e-01,\n",
       "         5.54244054e+00,  1.38805412e+01,  1.49268387e+04],\n",
       "       [ 1.00000000e+00,  1.00000000e+00,  1.00000000e+00,\n",
       "         7.03326489e+00,  1.40146866e+01,  1.48894745e+04]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import torch\n",
    "\n",
    "# feature lengths \n",
    "t_len = 88\n",
    "x_len = 61\n",
    "y_len = 30\n",
    "\n",
    "# Agnostic code in case you have a GPU\n",
    "torch_dev = 'cuda' if torch.cuda.is_available() else 'cpu'\n",
    "\n",
    "# Transform dataframe to numpy\n",
    "np_dataset = data.to_numpy()\n",
    "np_dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "#Fit Transforming the X (time, latitude, longitude)\n",
    "X = np_dataset[:10*1830, :3]\n",
    "y = np_dataset[:10*1830, 3:]\n",
    "\n",
    "\n",
    "# Dividing in Train and Test\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.5, shuffle=False)\n",
    "\n",
    "# Fit transform y_train data\n",
    "y_train[:, 0], alpha_u, beta_u = minmaxscale(y_train[:, 0], -1, 1)\n",
    "y_train[:, 1], alpha_v, beta_v = minmaxscale(y_train[:, 1], -1, 1)\n",
    "y_train[:, 2], alpha_G, beta_G = minmaxscale(y_train[:, 2], -1, 1)\n",
    "\n",
    "# Transform test data based on previous transfromation to avoid data leakage\n",
    "y_test[:, 0] = alpha_u*y_test[:,0] + beta_u\n",
    "y_test[:, 1] = alpha_v*y_test[:,1] + beta_v\n",
    "y_test[:, 2] = alpha_G*y_test[:,2] + beta_G"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x1d25bd0c890>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg4AAAGdCAYAAAB6oftJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAABIZklEQVR4nO3dfXRU1bk/8O8585r3EAgZUgIGsSBVULHGFNuLyBVYXS4srFatt0WvC29tsJVo1bhUBF9i7bpqvSvCrZdCu26V1v4utraVVqngrxpQsPn50pYrmDaxkKDYZPI2r2f//pgydXjz7GcmZw7D9+Oatcxk9pw9Z84MO3s/z34MpZQCERERkQ1mvjtAREREJw8OHIiIiMg2DhyIiIjINg4ciIiIyDYOHIiIiMg2DhyIiIjINg4ciIiIyDYOHIiIiMg2b747cCTLsrB//36UlZXBMIx8d4eIiFxMKYWBgQHU1tbCNEfvb+FIJIJYLJb18/j9fgSDwRz0KH9cN3DYv38/6urq8t0NIiI6iXR3d2PixImj8tyRSAT1k0vRczCZ9XOFQiF0dnae1IMH1w0cysrKAADz/88y+Er8ttslLP2R5tjAkHYbADiv/C/abQatgHabvniJdpty74h2m5CvX7sNAJwX7NJuU2Va2m0648XabdpHTtduAwC/+2CqdpvBmP57Oyao/z7NrPirdpszg/u12wDABG+fdptiU/+vsX5L/8uzN1Gp3QYA3h75hHab/rh+/z6MlWq3OTii/1kHgL5h/c/GyIhPu00ypv9PhYrpfycbCf1ZZisSwV/vuD/9b8doiMVi6DmYROfuySgvk89qhAcs1M/+C2Kx2Kk1cPjrX/+K2267Dc899xyGh4cxdepUbNiwAeeffz6A1LTRqlWr8MQTT6Cvrw9z5szB2rVrccYZZ9h6/sPLE74Sv9bAwRAMHPzBuHYbACgq1f8QJSz9D2sgrt8m6E1otynyycaPpUX651zymSuJ6zcKevTPHQB4R/QHAV6ffhtfUP8vl0Cp/msqLvJotwGAEq/+OS8RTBMnLP3+FSdkrykguCb8cfvfQYf5ovptvKb+NQQAHui3Mw39/imPYODgcWbgkG7rwNJ2eZmZ1cChUGidgb/97W+YM2cOfD4fnnvuOfzhD3/Av//7v2PMmDHpxzz00EN47LHHsG7dOuzcuRMlJSVYsGABIpFIzjtPRETklKSysr4VAq1h5Le//W3U1dVhw4YN6fvq6+vT/6+UwqOPPoo777wTixcvBgD88Ic/RE1NDZ555hlceeWVOeo2ERGRsywoWJAXlM6mrZtozTj8/Oc/x/nnn48vfvGLGD9+PM4991w88cQT6d93dnaip6cH8+fPT99XUVGBhoYGtLe3H/M5o9EowuFwxo2IiMhtrBz8Vwi0ZhzeffddrF27Fs3Nzbjjjjvw2muv4Rvf+Ab8fj+WLVuGnp4eAEBNTU1Gu5qamvTvjtTa2orVq1cfdf9Q3AevxvqiaeiP5AYTsnXFt4b0g6wCpn7swbTiY5+zE5lb/L/abab7ZIFZFvSDe9oFK1avR07TbvOXkbH6BwJwetkH2m0ayt7VbnN24D3tNlWmfkxOmSmLB5D8ZRRX+m0GLEGbpCyozFL6a+AjSf14gL6Yfv/6hou02wDASEQ/bsMSxAwhKohXEARHmnFBnIKgb5QdrTNuWRbOO+88PPDAAzj33HNx/fXXY/ny5Vi3bp24Ay0tLejv70/furu7xc9FREQ0WpJKZX3T9dJLL+Gyyy5DbW0tDMPAM888k/H7a665BoZhZNwWLlyYo1d8bFoDhwkTJmDGjBkZ95155pno6kql5YVCIQBAb29vxmN6e3vTvztSIBBAeXl5xo2IiMhtDsc4ZHPTNTQ0hFmzZqGtre24j1m4cCEOHDiQvj311FPZvMyPpbVUMWfOHOzZsyfjvv/93//F5MmTAaQCJUOhELZu3YpzzjkHABAOh7Fz507ccMMNuekxERHRKWLRokVYtGjRCR8TCASO+8f5aNAaOKxcuRKf+cxn8MADD+BLX/oSXn31VXzve9/D9773PQCpPNqbbroJ9913H8444wzU19fjrrvuQm1tLS6//PLR6D8REZEjLCgkc5BVcWQSQCAQQCAgi7kDgG3btmH8+PEYM2YM5s2bh/vuuw9jx8rivOzQGjh8+tOfxubNm9HS0oI1a9agvr4ejz76KK6++ur0Y2699VYMDQ3h+uuvR19fHy666CJs2bLlpN4li4iIKFfpmEeWVVi1ahXuuece0XMuXLgQS5YsQX19Pfbt24c77rgDixYtQnt7OzweWXD0xzGUEkRrjKJwOIyKigrM/j8r4S2xPwIrC+iH61cF9Lf9BYDZFfpbTp9frB95P8Wrn5oq2b1PEg0PAN2CrJTXRuo//kFHGBBsSzzF/752GwA4w9/78Q86giTbISjY5S5o6H8JBAzZrqCStLG9cf02zw/N+PgHHeHdkWrtNgAwINg++oOofsZR76D+1scDQ7I/rBJR/WvCium3MSL63yseQRszJtty+t3Vd6C/v3/UYuQO/7u0708hlGWxc+TAgIXTp/egu7s7o692ZxwMw8DmzZtPOIP/7rvv4vTTT8cLL7yASy65RNzXE2EeCxERkQ25yqo4MiEgm2WKI02ZMgXjxo3D3r17c/acR3JdkSsiIiI3sv5+y6b9aHvvvfdw6NAhTJgwYdSOwYEDERGRSw0ODmbMHnR2dqKjowNVVVWoqqrC6tWrsXTpUoRCIezbtw+33norpk6digULFoxanzhwICIisiGZZVaFpO2uXbtw8cUXp39ubm4GACxbtgxr167FG2+8gR/84Afo6+tDbW0tLr30Utx77705Xf44EgcORERENiRV6pZNe11z587FiXIYfv3rX8s7JOTagUNZIAJvwP5ZltQYmFasH0EPABP9h7Tb9CWLtdv8P0t//3pTsIoWUfr73QNAxNJvV+Pr125ztke/rkOxGdNuAwCVDtWDGGPqXw9xpV/v5JAlK2e/I1Lz8Q86wvbwdO027w6O024Ts2QpZglLPxb80LB+VsXwiH59C0l2BACohKCGhCDbwUjqZzsYCf02go8flKCN1MkQ4+AEZlUQERGRba6dcSAiInITCwaSEFTw/Ej7QsCBAxERkQ2WSt2yaV8IuFRBREREtnHGgYiIyIZklksV2bR1E9cOHE4vOwR/qf2o/THe4VHsTab3E/r7oQcNQbS+R7+WhsfBC9NvJLXbjPUMarepNPXPQ7Gp3zcAKBPU+pDUkPhzYkC7TUdUv2zunqh+pgMA/GlQf9e5fWH9anzRuP5XkMeUzff2D+vXg4hG9DOHlCADQcWFxYjikswFQVaFIHNB8PUA0RYJDk7/c+CQwqUKIiIiss21Mw5ERERuYikDlsoiqyKLtm7CgQMREZENXKpI4VIFERER2cYZByIiIhuSMJHM4u9tWci2+7h24PBOeBy8SfvVvfwe/bek2Cvb5Nxv6tcM8Br6u5QHPPrHiQv28a8O6Gc6AEDI3ydqp8vjlYRN62diAECf4Pz9Oa6fTbA3ql8L4lC8TLvNOwPV2m0A4C/9Y7TbjET1MxC8HkFtFUEtCABIRvTfW1EtCEFBAiMm+8dIdCz9rxUYljNT7EpwGiRtpFSWMQ6KMQ5ERESnDsY4pDDGgYiIiGzjjAMREZENSWUimcXaSLJAalVw4EBERGSDBQNWFhP1lpPbXI4iLlUQERGRbZxxICIisoHBkSmuHTgcGi6BB/bTMQ3B++ExBblMAMoDEe02pqE/RSVJF63066chekXVaIBiM6bdZsiy/54e9nb0E9pthi1Zyt7BmH4Bs/djpdptwvEi7Tb7B/X7dmigRLsNAMRHnPlqiMT0UyTNIVlBKI+gIJRoZtnBfxsMyWuSHEf2VSk4kENthLKPceBSBREREZ1iXDvjQERE5Cap4MgsilxxqYKIiOjUYWW55TSzKoiIiOiUwxkHIiIiGxgcmeLagYNShlZBkCKffgZCTcmAdhsAqA7qF4Wq9OlnO5R79NtI1tAiln5xIgDYMxzSbtMXL9ZuM5jQz8Toiwa12wDAgKBdwtL/IonG9T96Q2HBaxIUdgIACKL1PRH98+Ab1j+OYxH+gChi38k6RoJkLccyRZTk0nPq9QhZMLkBFFw8cCAiInKTpDKQzGJkmE1bN2GMAxEREdnGGQciIiIbkllmVSS5VEFERHTqsJQJK4vgSKtAgiO5VEFERES2uXbGoap4CN6ShO3HS+o6+E1ZjQZJgMtgQr92QjiuH0U/IMhACMdkGQj9gnaSbIJIXD/rQ3IcALCS+u+tldAffycj+v0zRvTD1D1RWTCWd0iSVaF/HEmGhPQPPklJFlEZF0kGgvA1KcFlLjmWJTiOU+dbyb7GRbhUkeLagQMREZGbWMguM8LJTOLRxKUKIiIiso0zDkRERDZkvwFUYfytzoEDERGRDdlvOV0YA4fCeBVERETkCNfOOEQSPngT9qPpR+L6WQsDXvtZGx/1/kiJdhtDsKm8Tq2Ow+KWfuT9UEz/3AHA8Ih+u0RcsIG9IBBZkukAABgRfCQS+u+TKXhNnhH940iyI1LH0m8jqZsgyaqQ1qow9ROvRCSxc8JyMaKMAsmxPIKvSsn1IDl3KqbfRsqCIaoH9NH2hcC1AwciIiI34VJFCgcORERENmS/j0NhDBwK41UQERGRIzjjQEREZIOlDFjZbABVIGW1OXAgIiKywcpyqYL7OIyygUgQHtN+3QVJ0bFEsli/EQCPqR/WLcqQEGQgWHH9C1NJMxAkJFHTgtdkjshekymo7WBYgqwKQSS4V1ILQphJ4FS2g6iNcLt/SVyaKJbNwT8qJefPI7j2RNkOgn9dJBkfVqHs43wcL730Er7zne9g9+7dOHDgADZv3ozLL788/XulFFatWoUnnngCfX19mDNnDtauXYszzjhj1Pqk9bG45557YBhGxm369Onp30ciETQ1NWHs2LEoLS3F0qVL0dvbm/NOExEROe1wWe1sbrqGhoYwa9YstLW1HfP3Dz30EB577DGsW7cOO3fuRElJCRYsWIBIRPCXhk3aY8JPfepTeOGFF/7xBN5/PMXKlSvxy1/+Ek8//TQqKiqwYsUKLFmyBC+//HJuektERJQnSRhIZjGlJGm7aNEiLFq06Ji/U0rh0UcfxZ133onFixcDAH74wx+ipqYGzzzzDK688kpxX09Ee+Dg9XoRCoWOur+/vx/r16/Hk08+iXnz5gEANmzYgDPPPBM7duzAhRdemH1viYiITnLhcDjj50AggEDA/tL8YZ2dnejp6cH8+fPT91VUVKChoQHt7e2jNnDQnjd55513UFtbiylTpuDqq69GV1cXAGD37t2Ix+MZL2D69OmYNGkS2tvbj/t80WgU4XA440ZEROQ2uVqqqKurQ0VFRfrW2toq6k9PTw8AoKamJuP+mpqa9O9Gg9aMQ0NDAzZu3Ihp06bhwIEDWL16NT772c/irbfeQk9PD/x+PyorKzPafNwLaG1txerVq0WdJyIickoSsuWGj7YHgO7ubpSXl6fvl8w25JPWwOGj6ywzZ85EQ0MDJk+ejJ/85CcoKioSdaClpQXNzc3pn8PhMOrq6hAZ8sFU9mshWAlBDQRpNK7kuonpB8UYEUEbQd0EFRCeCMF5kPTPG9Fv4xG0AWT1DJxqI7penYzwl9SqEGZISChJmRRBVoWgXIyob4Dw/Dl0zi1Jzp7kej0Jt0YoLy/PGDhIHQ4b6O3txYQJE9L39/b24pxzzsn6+Y8nqzy8yspKfPKTn8TevXsRCoUQi8XQ19eX8Zje3t5jxkQcFggE0icxVyeTiIgo1/KRVXEi9fX1CIVC2Lp1a/q+cDiMnTt3orGxMafH+qisXsXg4CD27duHCRMmYPbs2fD5fBkvYM+ePejq6hrVF0BEROSEw0WusrnpGhwcREdHBzo6OgCkAiI7OjrQ1dUFwzBw00034b777sPPf/5zvPnmm/jqV7+K2trajL0eck1rMumWW27BZZddhsmTJ2P//v1YtWoVPB4PrrrqKlRUVOC6665Dc3MzqqqqUF5ejhtvvBGNjY3MqCAiopOeyrKsthK03bVrFy6++OL0z4eX9pctW4aNGzfi1ltvxdDQEK6//nr09fXhoosuwpYtWxAMBsX9/DhaA4f33nsPV111FQ4dOoTq6mpcdNFF2LFjB6qrqwEAjzzyCEzTxNKlSxGNRrFgwQI8/vjjo9JxIiKiQjd37lyoE2yNbBgG1qxZgzVr1jjWJ62Bw6ZNm074+2AwiLa2tuPucEVERHSyki43fLR9IXBtrQr1twBURCNFxSMJ6dZvAgBGXFDPIKnfxjeg3yZpPxHlI2QXsyRzwTuifxwzKmgjrdEguIxEx5Ls/S+JvJdmDjn0/SY5d5L6DIAsyl90zgXXkKDcSepQknoQgtck+Vw4lcXi5L/FrI6ZUhjDHyIiInKEa2cciIiI3CSZZVntbNq6CQcORERENnCpIqUwhj9ERETkCM44EBER2WDBhJXF39vZtHUT1w4czKgB09CY1tF5bJa8w/rHkkSPS2a1vAn9Np4+2bnzCLIdJG0MwWsSk9TfEET5JyU1bQTHkZxvAM7NRUqSoYRZFZKXZAn650l+/GOOZAraAIAleFGmJFPEoSwgyyc4joPfD0llIJnFckM2bd2kMIY/RERE5AjXzjgQERG5CYMjUzhwICIiskFlWeFScedIIiKiU0cSBpJZFLnKpq2bFMbwh4iIiBzh2hkH75ABT0Inq0L/GNJZIzMmaCOsnaBLVAtC2DdPTD/k3CM4d5JlQUldAgBQpv7BRJHggmvPkETeC69xST0DCVE2gSDTAQDMuH5DQ+c7KAuimhgQZkg4RHINSWpvSGrZSFkquzgFSZaOG7l24EBEROQmVpYxDtm0dZPCeBVERETkCM44EBER2WDBgJVFgGM2bd2EAwciIiIbuHNkCpcqiIiIyDbXzjj4+wGPxn7+koGcNJJZQhITEwjrt/FEJJv/6zcBZDUDJJHWhoP1DCQHMyxB7RLB/vqSayjp128jJriOkoLrwSONTE/qd9CxWDbhZ1DSP0nGkegPZYfaOBlvyODIFNcOHIiIiNzEQpZbThdIjENhDH+IiIjIEZxxICIiskFlmVWhCmTGgQMHIiIiG1gdM4UDByIiIhsYHJni2oGDJw54dEpVSKPoBZKC2gSS/nmHBbUgBPvxWzon+iMkWSmiDAlBPQPpuF5yGZmS4gmirAr9VyXNqhBlHEmyHQRtxHVIHMowkdTfENUhEbaTnIekRoZbNiTfD5aD3/2U4tqBAxERkZtwqSKFAwciIiIbuOV0SmEsuBAREZEjOONARERkA5cqUjhwICIisoEDhxTXDhySfgA6tSoEiy7STAxPVL+NPyzIdhC8O5IMCUkkMwAYgswAKEktCMFhhItwZkJQz0ASCS7JSLEEWTYxYcaM4Pw5VvtFeL06VtfBoWwjwLkaPZL+ia4hN193lObagQMREZGbcMYhhQMHIiIiGzhwSGFWBREREdnGGQciIiIbFLLbi0EYyuI6HDgQERHZwKWKFNcOHCwPYGhEyyaDgoMI30NJNoF3RD81IFGsv5IkipgW1LcA5JHg+gdy6DhCogwOh16TkZS9SaaoloZ+G1GdFCevB8GxJPVBpLUqJJzKdpD8eS3JoHLy3HHgkMIYByIiIrLNtTMOREREbsIZhxQOHIiIiGzgwCGFSxVERERkGwcORERENihlZH3Tcc8998AwjIzb9OnTR+nV2efapQpvBPDoRNhKZoCEWQGeuH4bSbSwJ6LfwWRQEqUurGcgCpsW1IKQ7McvHRI7le0gqZPiYBK4KLpdVPtFUsNF9iY5lQUkifI3hZkBovcpIjiQQ+dOlBUmqB0kZcHIah8HSdtPfepTeOGFF9I/e735/2c7/z0gIiKiY/J6vQiFQvnuRgYuVRAREdlwODgymxsAhMPhjFs0evxpk3feeQe1tbWYMmUKrr76anR1dTn1co+LAwciIiIbchXjUFdXh4qKivSttbX1mMdraGjAxo0bsWXLFqxduxadnZ347Gc/i4GBASdf9lG4VEFEROSg7u5ulJeXp38OBALHfNyiRYvS/z9z5kw0NDRg8uTJ+MlPfoLrrrtu1Pt5PBw4EBER2ZCrfRzKy8szBg52VVZW4pOf/CT27t0r7kMuuHbg4BtS8GjUUDCTzkTrA7JaFZJIcG9EEjKtv/qU9Auj1AURwk7teS+uQyI45ZK6DoVI8t7Kzp00xF//opBcD5YgM8ASfhOLshAEnw1TkkkmeZskX3kOZhtJUiqPbJ+NwcFB7Nu3D1/5yleyep5sMcaBiIjIBpVlYKTuwOGWW27B9u3b8ec//xmvvPIKvvCFL8Dj8eCqq64apVdoT1YDhwcffBCGYeCmm25K3xeJRNDU1ISxY8eitLQUS5cuRW9vb7b9JCIiOqW89957uOqqqzBt2jR86UtfwtixY7Fjxw5UV1fntV/ipYrXXnsN//mf/4mZM2dm3L9y5Ur88pe/xNNPP42KigqsWLECS5Yswcsvv5x1Z4mIiPJFAVBZLI3oNt20aZP8YKNINOMwODiIq6++Gk888QTGjBmTvr+/vx/r16/Hww8/jHnz5mH27NnYsGEDXnnlFezYsSNnnSYiInLa4Z0js7kVAtHAoampCZ///Ocxf/78jPt3796NeDyecf/06dMxadIktLe3H/O5otHoUZthEBERkTtpL1Vs2rQJr7/+Ol577bWjftfT0wO/34/KysqM+2tqatDT03PM52ttbcXq1auPut+MK3g0Jna8Uef2vJdEjyeKJBHd+m08McE8mnDqTYkWuiSFJ/Q7KKubIIvQltQmEJH0zZK9uZLzJ/lcWJKMHmlkuuA6kqTemaL3Sb8N4FxGgeUTtJFkb7g8QynfWRVuofVR7+7uxje/+U386Ec/QjAYzEkHWlpa0N/fn751d3fn5HmJiIhyKVdbTp/stAYOu3fvxsGDB3HeeefB6/XC6/Vi+/bteOyxx+D1elFTU4NYLIa+vr6Mdr29vcct0hEIBNKbYUg3xSAiIiJnaE02X3LJJXjzzTcz7rv22msxffp03Hbbbairq4PP58PWrVuxdOlSAMCePXvQ1dWFxsbG3PWaiIjIYUplmVXh4GZVo0lr4FBWVoazzjor476SkhKMHTs2ff91112H5uZmVFVVoby8HDfeeCMaGxtx4YUX5q7XREREDmOMQ0rOt5x+5JFHYJomli5dimg0igULFuDxxx/P9WGIiIgoD7IeOGzbti3j52AwiLa2NrS1tWX71Fr8/fqh7dJI5kSxfvh4XNAmGRDUgpBEJQunzyQ1O0xLeNI1iaL1hSQ1RUxB9ouZ0D+O5ZVtDis5f5Icdck1ZAjqMwCyz7tOvZw0SV0HZz4WAGTZLyqq38apbCPJZ0mKMw4pri1yRURE5CaWMmDkoDrmyY4DByIiIhsYHJnC6phERERkG2cciIiIbEjNOGQT45DDzuQRBw5EREQ2MDgyhUsVREREZJt7ZxxMA8pjf3SWDOiPgTyCNDoA8A7p5xmZCWcKNSWC+ufBkl4FgsGzERMeS5MnIpsTlLxPEvFS/ZxCZeq3MZKy1yNKQxRcD5KU42RA/zgARGnHpiCl0KkCYYCs3pdX8NnwCIoISq49Zeq/oETcuVxWBXH2erp9IXDvwIGIiMhFuFSRwqUKIiIiso0zDkRERHZwrQIABw5ERET2ZLlUIQpKcSEOHIiIiGzgzpEprh04WB69YjaSwlPSSGZJxLmkEIth6bfxCbICkoJMDACiT4EkalpyHiyvbGQfrdDPXDAF0eO+IWciwQ1hlohORlP6WIKXJDl3ngH94wBwbprYwT8qJd9hkuJTkj+UDUNwDUmSeQrkH+OTiWsHDkRERG7CrIoUDhyIiIjsUEZ2cQoFMnBgOiYRERHZxhkHIiIiGxgcmcKBAxERkR3cxwGAiwcOZlJvn3hLPxgeEGYTJIoEbYKSsGT9JqYg4yPQL4vw9wgyRZTgfZLU35AcB5Dt4+9U7RLLJ9jHv1h2IiyffhtJfRB/v/65k9aYkWTnSNpIMlKSfuGqsSBzQUJy7TmVXSJ5jyg7rh04EBERuQmzKlI4cCAiIrKLExzMqiAiIiL7OONARERkA5cqUjhwICIisoNZFQDcPHDQfIOSAf2RnFeYTSCqTTCsfxzJ4FRSs2N4vCzy3kzoH6voA/0oet+gYHN94cBeEhEvOZakloYkYyZ4KKbdBpCdB0kdEjOu/xkUvUcALK/+9epUxL643oJDdVxEmSyCt0l0DQnrscgYyC5dpDBmHBjjQERERLa5d8aBiIjITbhUAYADByIiIns4cADApQoiIiLSwBkHIiIiO1hWG4CLBw7h0wx4NDIlTEHweGSsbMKlcq9+lH9SUG8heCih3UYacS7hHdaPtPaHBa9JsB9/MiB7byXnT5Lt4P9bRLuNkdA/35Zf9hG3JLUTPPptLJ+gjfAaV4JsAkm2gyHIFJHWnDAEGV7eiCBLSXAeROdbkCViCj4XUqyOmcKlCiIiIrLNtTMORERErsLgSACccSAiIrLncIxDNjeBtrY2nHbaaQgGg2hoaMCrr76a4xemhwMHIiIil/rxj3+M5uZmrFq1Cq+//jpmzZqFBQsW4ODBg3nrEwcORERENhgq+5uuhx9+GMuXL8e1116LGTNmYN26dSguLsb3v//93L9Am1wb45D0AfDbf7wSDIEsn34bABiu0a/tkAjoH2ewVr+D3hH94wT6hQtvgkjwkWqNN/XwYQSR49L9632CrA/vcFz/QIJzF68I6h9GGMYtOeeSrA9ZtUBhxozgVFg+/f5ZRfrfD0nBcQBACcrMeEfcmw0lvV4dk6MYh3A4nHF3IBBAIHD0PxKxWAy7d+9GS0tL+j7TNDF//ny0t7dn0ZHscMaBiIjIjhzFONTV1aGioiJ9a21tPebhPvjgAySTSdTU1GTcX1NTg56enlF/ucfj2hkHIiKiQtTd3Y3y8vL0z8eabXAzDhyIiIjsyNFSRXl5ecbA4XjGjRsHj8eD3t7ejPt7e3sRCoWy6Eh2uFRBRERkh8rBTYPf78fs2bOxdevW9H2WZWHr1q1obGzM8sXIccaBiIjIpZqbm7Fs2TKcf/75uOCCC/Doo49iaGgI1157bd76VDADB0lWRTIgm3MamiCpZ6B/nEC/fpuiD/QjpmNlsomnRFD/PPiGBdH6gsjx4gP6tSAAWZ2BZKkwPUeTLxzVbmMOys6DiKDOgKS+hRLW30iU668jJ0r039ukoJaGpEbD31tqt0gUCz7vhv558A0KMpQGBEWHEoI2UnnYOfKKK67A+++/j7vvvhs9PT0455xzsGXLlqMCJp1UMAMHIiKiUZWn6pgrVqzAihUr5MfNMcY4EBERkW2ccSAiIrJBuvvjR9sXAg4ciIiI7GB1TABcqiAiIiINrp1x+MTLUXi99gNJkj79MZAVEO55L4hv8Q0ltdt4B/VTMeJl+tHPyiPLCvBE9YfPkj3vPRH9c2f5BakYACBtp8n7N/1sB7NvQP9AEf1MDABAQv+cwxR8MCT1DIr0a3YAgNfSz5gxBDVPjKT+50n6XSQ4lKg+iCTrQ1Krwozofz+YScG1SlnRulrXrl2LmTNnpne9amxsxHPPPZf+fSQSQVNTE8aOHYvS0lIsXbr0qB2viIiITkYGsqyOme8XkCNaA4eJEyfiwQcfxO7du7Fr1y7MmzcPixcvxttvvw0AWLlyJZ599lk8/fTT2L59O/bv348lS5aMSseJiIgclaMiVyc7raWKyy67LOPn+++/H2vXrsWOHTswceJErF+/Hk8++STmzZsHANiwYQPOPPNM7NixAxdeeGHuek1ERER5IQ6OTCaT2LRpE4aGhtDY2Ijdu3cjHo9j/vz56cdMnz4dkyZNymvdcCIiopxwuFaFW2kHR7755ptobGxEJBJBaWkpNm/ejBkzZqCjowN+vx+VlZUZj/+4uuHRaBTR6D8CuMLhsG6XiIiIRh/TMQEIBg7Tpk1DR0cH+vv78dOf/hTLli3D9u3bxR1obW3F6tWrj7rfdygCr8f+WZbkBVjFsmwCJdiL3hTUQEBSkLXQpx9FL2kDAOaIoACHQHxMkXabpDBKXZLB4RNkSBiSDImk/jWk4rL3SA2N6DcS1J0wS0v0j+MgpZHZlW4juPQsYTKP5FhQggwJQcaM5dNvE6/S/6wnEoURN3Ay0b7s/H4/pk6ditmzZ6O1tRWzZs3Cd7/7XYRCIcRiMfT19WU8/uPqhre0tKC/vz996+7u1n4RREREoy2rjIosd510k6w3gLIsC9FoFLNnz4bP58uoG75nzx50dXWdsG54IBBIp3cevhEREbkOYxwAaC5VtLS0YNGiRZg0aRIGBgbw5JNPYtu2bfj1r3+NiooKXHfddWhubkZVVRXKy8tx4403orGxkRkVREREBUJr4HDw4EF89atfxYEDB1BRUYGZM2fi17/+Nf75n/8ZAPDII4/ANE0sXboU0WgUCxYswOOPPz4qHSciInIUgyMBaA4c1q9ff8LfB4NBtLW1oa2tLatOERERuQ2rY6a4tlZFstQHw+u3/XjPiKAGwqBwH3/BHuySfdslDEHENOLO7fWerBRkSBTpX6aGICMFEGa/OHj+dBle2UdcCTIkrOFh7TZmkf71AENYYyagn0UlufYSxfopEsmA7PtB0s4QXOISytR/n5I+wXsUA/B/tZtRFlw7cCAiInKVbLeNPhW3nCYiIjplMcYBAAcOREREtjDGISXrfRyIiIjo1MEZByIiIju4VAHAxQMH5TWhvBoTIpZ+qLAREdZaEEScSzIxIKiJIeKTbZSfLLaf9fKPNvqXnKRegCfqXKaDEpwHCSOunzkkuu4AmCX633BmaLx2G6tSv1aFJFofSGVqabcR1DxJBgU1Goqdy6qwBCV6JDUxEoKEmViFfhtLv1SMXLbbRhfIwIFLFURERGSba2cciIiIXIVLFQA4cCAiIrKHAwcAXKogIiIiDZxxICIisoH7OKS4duDgGUnA49WIIpfUjwg49/Il++Rbfv1sB8PSvzIldQkAIBkQ7MnvE9RAkCV9yAiuI1EtDUEmi6QOiRELareRktQHsfz6505yvgEgXq7fLl6if70mAtpNkCiSZVUkBG9vUpDtkCgWtCkVZLpV6dcPsoadTKsggEsVREREpMG1Mw5ERESuwuBIABw4EBER2cIYhxQOHIiIiOwqkH/8s8EYByIiIrLNtTMOZiQB02O/loQSRKkny4QR55IMDkHdCUunVsdhgiaWsCaG5Rfs4y/Y+98SXKVmUnZpG35BJLjgL5BkUP88xIv123iisj+PJK9JktEjufYsQe0SQJa5EC9xphZEUljuRHIsyedJeQTvbZF+vZhgUL8eS9IS1HCRYowDABcPHIiIiNyEMQ4pXKogIiIi2zjjQEREZAeXKgBw4EBERGQLlypSuFRBREREtrl3xiGZBKARlSvIQFCSrAXIajRI6gwoU5C9IQg4V4L6EQBg+SQR8c4cJ1ohK3BhxvTPhSSbIFauf5xohaCOhjBxSBKtbwiC272CMgOGfrA+AOm1p99GCT5OSliPRfIXrKR/EuaI/ouK+vXTS6wR/UwoMS5VAHDzwIGIiMhNOHAAwKUKIiIi0sCBAxERkQ2HgyOzuY2m0047DYZhZNwefPDBnB+HSxVERER2nARLFWvWrMHy5cvTP5eVleX8GBw4EBER2XESDBzKysoQCoVG9RjuHTh4PKmbXZL6EYKsBUBWZ0ASyewVRAt7YvptLGFQsuWRZCBI0j70m0j3/k8EBP0TXHuSDIm44A+HWIXsm0p59dt5hgW1SwRZH9KsAEnmgqhmhyDrw5AmBggyWUTHEnxuzZj+YZLDgn+SRtz7z9jxhMPhjJ8DgQACgUBOnvvBBx/Evffei0mTJuHLX/4yVq5cCa83t+fo5DvjREREeZCrDaDq6uoy7l+1ahXuuece+RP/3Te+8Q2cd955qKqqwiuvvIKWlhYcOHAADz/8cNbP/VEcOBAREdmRo6WK7u5ulJeXp+8+0WzD7bffjm9/+9snfNo//vGPmD59Opqbm9P3zZw5E36/H//2b/+G1tbWnM1oABw4EBEROaq8vDxj4HAiN998M6655poTPmbKlCnHvL+hoQGJRAJ//vOfMW3aNN1uHhcHDkRERDbko1ZFdXU1qqurRcfr6OiAaZoYP368qP3xcOBARERkh4uzKtrb27Fz505cfPHFKCsrQ3t7O1auXIl/+Zd/wZgxY3J6LNcOHAamlMHrsx9yHfwwrn0MSc0JAEgUObRvliBa30gKMh0Ssqs5Uax/rHiRoL6FoFaFdO9/p2oTxO3NUma2KdF/n6xiYbi+IDNA+QT9ExRXSZQ4V5vASAiyCeL6bTxR7SapYzmU9WEKsjcswefCjAgaRbmPIZCKk9i0aRPuueceRKNR1NfXY+XKlRlxD7ni2oEDERGRq7h4xuG8887Djh07Ru8AH8GBAxERkQ3G32/ZtC8EnOMhIiIi2zjjQEREZIeLlyqcxIEDERGRDflIx3Qj1w4cDlwMmEUaDSTh8JK6CQACHzqzwlP+rv5VJtmHXnLqACBeLMh2EFxxkroTiWL9NgBEi5CSKHWn/vIwo7JrXJDsIDoPyaD+iVB+2cmTROxLvugldT5UTLj6LcmqEHxHSNpIsksMSe2NiIORA5xxAMAYByIiItLg2hkHIiIi1ymQWYNscOBARERkA2McUrhUQURERLZxxoGIiMgOBkcC4MCBiIjIFi5VpLh24GBGTJgaKylm7bD2MRKH7BfRymhX5EwK2cAkZ1aSTEk6IQBLUkhK8JIkqYGSwlOALO1MUlArXiq4hgRFpKTpmJJWovfJI8l31G8CAKYgbU+Ujil4TdJ/UKTXuS5Jqq3ksyR6b2OCNpQVrcuutbUVn/70p1FWVobx48fj8ssvx549ezIeE4lE0NTUhLFjx6K0tBRLly5Fb29vTjtNRETkOJWDWwHQGjhs374dTU1N2LFjB55//nnE43FceumlGBoaSj9m5cqVePbZZ/H0009j+/bt2L9/P5YsWZLzjhMRETnp8FJFNrdCoLVUsWXLloyfN27ciPHjx2P37t343Oc+h/7+fqxfvx5PPvkk5s2bBwDYsGEDzjzzTOzYsQMXXnhh7npOREREjstqhay/vx8AUFVVBQDYvXs34vE45s+fn37M9OnTMWnSJLS3tx/zOaLRKMLhcMaNiIjIdbhUASCLgYNlWbjpppswZ84cnHXWWQCAnp4e+P1+VFZWZjy2pqYGPT09x3ye1tZWVFRUpG91dXXSLhEREY0eDhwAZJFV0dTUhLfeegu/+93vsupAS0sLmpub0z+Hw2HU1dXBjBkwTftR0J49JdrHVoLsCAAI/E0/OjteKigIJYjWj47VD2X2Dgsj7xP67bwj+seJlwuyCQQFdgBhEa5iSWUxwbkb1B/nS6PuJZkBScHnyRRcQ6bwevU4FX0vKFhlCoo7Ac5lAYn6JynAJcnwigraCDEdM0U0cFixYgV+8Ytf4KWXXsLEiRPT94dCIcRiMfT19WXMOvT29iIUCh3zuQKBAAKBgKQbRERE5DCtv0eUUlixYgU2b96M3/72t6ivr8/4/ezZs+Hz+bB169b0fXv27EFXVxcaGxtz02MiIqJ84FIFAM0Zh6amJjz55JP42c9+hrKysnTcQkVFBYqKilBRUYHrrrsOzc3NqKqqQnl5OW688UY0NjYyo4KIiE5qhlIwlPxf/2zauonWwGHt2rUAgLlz52bcv2HDBlxzzTUAgEceeQSmaWLp0qWIRqNYsGABHn/88Zx0loiIiPJLa+CgbIyWgsEg2tra0NbWJu4UERGR67DIFQAX16qwfArQ2pvfmX3oASAyXrAXvSAqOT5WP8TY96H+WyraUx6yCGhJlL9vQBZFL5EUlC/xDDtTMMDy6l93VlB2kStT0E7QPyOiH+LvGRFmVQgyekQJ65JTLq1VIciQkHzejbhDx5FkiQj6JsWsihSHSqQQERFRIXDtjAMREZGrcKkCAAcOREREtnCpIoVLFURERGQbZxyIiIjs4FIFABcPHNS4KFSx/ehpSxAFrj5wbqtra6wkrUJ/QkgSeS9lJAV78kv2opdEZwtPgyhTRBDZLskukdTREGVHQFZDQgmuB6fqnaQOJmznBOHcryhzQXCNm05lLrj5PQKXKg5z7cCBiIjIVTjjAIAxDkRERKSBMw5EREQ2FcpyQzY4cCAiIrJDqdQtm/YFgEsVREREZJtrZxz8RXF4iuyHq5umfnjxcLkgHB4ARgTtovptzIhgXCcZCgprVSitWiIpScne+oJofSnJNKQSdC9ZrH/SlYMZM1aRfv/MYf03V3kE9S2E14Mkk8WpKH/JNQRA1j/BZST5XEjaWE7VBhFiVkWKawcORERErsKsCgBcqiAiIiINnHEgIiKywbBkm259tH0h4MCBiIjIDi5VAOBSBREREWlw7YxDtKcEZlHQ9uNF0aolksIJgPLpzzcZMf0xmqTOgCQTwxSU0QAAy6/fP0tSb0GQvSHZjx+QpVmbcUGNBktQC0LQOWkGgjnszN8Ukv7Fy2XH8g7ptzEkNRpc/ueYpOZJoki/jSem30byuXUyU4FZFSkuv8SJiIhc4vAGUNncRtH999+Pz3zmMyguLkZlZeUxH9PV1YXPf/7zKC4uxvjx4/Gtb30LiYTeX4+unXEgIiJyE7fPOMRiMXzxi19EY2Mj1q9ff9Tvk8kkPv/5zyMUCuGVV17BgQMH8NWvfhU+nw8PPPCA7eNwxoGIiKgArF69GitXrsTZZ599zN//5je/wR/+8Af893//N8455xwsWrQI9957L9ra2hCL2V9b4sCBiIjIDpWDG4BwOJxxi0ajjnS/vb0dZ599NmpqatL3LViwAOFwGG+//bbt5+HAgYiIyIbDSxXZ3ACgrq4OFRUV6Vtra6sj/e/p6ckYNABI/9zT02P7eQomxsEKCjIdBNHwgCwS3EjotzGjguNI1tCE++RL+peQZLIINvL3D8nGxIli/ROY8AuuPUFWhXdAUO9EmDETLxfUqhB8niQ1MVRMdsF6hwXtJE0E51xcEkPSUJDtIKnzIcnEkBwn6cwf6znV3d2N8vJ/pAcFAoHjPvb222/Ht7/97RM+3x//+EdMnz49Z/37OAUzcCAiIhpVOSqrXV5enjFwOJGbb74Z11xzzQkfM2XKFFvPFQqF8Oqrr2bc19vbm/6dXRw4EBER2ZCPrIrq6mpUV1fLD/oRjY2NuP/++3Hw4EGMHz8eAPD888+jvLwcM2bMsP08HDgQEREVgK6uLnz44Yfo6upCMplER0cHAGDq1KkoLS3FpZdeihkzZuArX/kKHnroIfT09ODOO+9EU1PTCZdLjsSBAxERkR0ur1Vx99134wc/+EH653PPPRcA8OKLL2Lu3LnweDz4xS9+gRtuuAGNjY0oKSnBsmXLsGbNGq3jcOBARERkg9s3gNq4cSM2btx4wsdMnjwZv/rVr7I6jmsHDrpvkOj9EETrA7LocV+/oFaFoHtWQP9MKFN6HvTbeAQ1EJSgJoYkOwKAKErdN6D/mpKC90lSG8QjzEAofk/wmor1j5P0C7KNhJkikoh9yRe9qN6CfpMUQf9ENSQE5aCjFfptEiX6bZTweiA51w4ciIiIXMVSqVs27QsABw5ERER2uDzGwSkcOBAREdlgIMsYh5z1JL+45TQRERHZxhkHIiIiO3K0c+TJzrUDB1WWgCqyHy5rDApeSoksHNcY8mm3sQTdk0SPS+pHSIkiziVR4GH9iTFJdgkAmJIsBMGhJMfxRPSPI54blbQTRN77hvTbSLIjAFkNCVOQISHK+hC+Jklmk6SNN6J/kUcrBLV2JBkpgutOyu3pmE7hUgURERHZ5toZByIiIldhVgUADhyIiIhsMZSCkUWcQjZt3YRLFURERGQbZxyIiIjssCAKAs5oXwBcO3Aww16YMfvdMwRR6iqunx0BAIagVoXy6k9RJSXvjiWIZBZezEmfJK1C0ERwGP/fZOkEXkHmQrxUv03SfgXbNOXRb+Md1m8DyK4JSdaH5DjS61Vy/iTHUoLPrSSbAJD1L+l35jj+sH6buOR6iOq3keJSRQqXKoiIiMg21844EBERuQqzKgBw4EBERGQPd44EwIEDERGRLdw5MoUxDkRERGSba2cclKG3J71XUKNBCfZsB2QZEpYggcMT0X9NyaB+3yR71wPCWhUJh7I+hDUaJHUQfAPOtJFkb0gyCQBhzQBJzQ7BcSRZAQBkGT2CiH1RLQ3h9ZqUfK/E9NsoU1BbJaZ/QXg/0G6CpOA4YlyqAODigQMREZGbGFZ2RbWcLMg1mrTHxi+99BIuu+wy1NbWwjAMPPPMMxm/V0rh7rvvxoQJE1BUVIT58+fjnXfeyVV/iYiIKI+0Bw5DQ0OYNWsW2trajvn7hx56CI899hjWrVuHnTt3oqSkBAsWLEAkIqkJTERE5BKHlyqyuRUA7aWKRYsWYdGiRcf8nVIKjz76KO68804sXrwYAPDDH/4QNTU1eOaZZ3DllVdm11siIqJ84T4OAHKcVdHZ2Ymenh7Mnz8/fV9FRQUaGhrQ3t5+zDbRaBThcDjjRkRERO6U0+DInp4eAEBNTU3G/TU1NenfHam1tRWrV68+umMDJjxx++MaSSSzJAMBgGy4JQg5l2RiKI8gq2JYNn60BFePpJ6BJDMgUazfBgBMQcS5J6HfRlKrwjuk30YYrA/LoboOkuwNUdYCIPrcis6DfhMYgmsIAExJO0EHJd9FUPoH8g/qf38p4bmTYK2KlLzv49DS0oL+/v70rbu7O99dIiIiOhpjHADkeOAQCoUAAL29vRn39/b2pn93pEAggPLy8owbERERuVNOBw719fUIhULYunVr+r5wOIydO3eisbExl4ciIiJylgJgZXErjAkH/RiHwcFB7N27N/1zZ2cnOjo6UFVVhUmTJuGmm27CfffdhzPOOAP19fW46667UFtbi8svvzyX/SYiInIUYxxStAcOu3btwsUXX5z+ubm5GQCwbNkybNy4EbfeeiuGhoZw/fXXo6+vDxdddBG2bNmCYDCYu14TERE5TSHLLadz1pO80h44zJ07F+oEJ84wDKxZswZr1qzJqmPKoxnRbDr3jliCWhVm3JkaDaL6Fn7ZuXPsNY3ot5Fu7Rqr0G8jiWyXZG9IMjEkWSyALMrfI6h5EivTbyOtrSIhqqUheJ+kJOdCkiEhucYlGTOJoOD7S1BHg7LDWhVERER2sMgVAA4ciIiI7LEg3xzlcPsCkPd9HIiIiOjkwRkHIiIiG5hVkcKBAxERkR2McQDg4oGD5VeATrS/YN1JGnkvibR2Kg1HsD08jKRs0c7ROgOaPFFZu5L9+m0SRYIDCU655DVJsjcA2WfD8gsOJPlcCNeYJZkiTl2v0q8HyWfQELwmyXHMuP6rilYKsipc+69Y4eIpJyIisoMzDgA4cCAiIrKHAwcAzKogIiIiDZxxICIisoP7OADgwIGIiMgWpmOmuHbgYPktIGB/eOaJ6K+6iCOmLf0hpyfqzH7qSlB3QlILInUwYTtdklMnfG8twSdCku0gqRcgqQXhiTn3RWUIDpUoktRW0T8OAPgl17ngNUkyUqTfRUqnns/fOZUNlQw4U8tGmh0n4vIYh/vvvx+//OUv0dHRAb/fj76+vqMeYxhHvy9PPfUUrrzyStvHce3AgYiIiOyLxWL44he/iMbGRqxfv/64j9uwYQMWLlyY/rmyslLrOBw4EBER2WEp2dTaR9uPotWrVwMANm7ceMLHVVZWIhQKiY/DrAoiIiI7Di9VZHNzgaamJowbNw4XXHABvv/970Np9oszDkRERA4Kh8MZPwcCAQQCAUeOvWbNGsybNw/FxcX4zW9+g69//esYHBzEN77xDdvPwRkHIiIiW7KdbUj9ZV9XV4eKior0rbW19bhHvP3222EYxglvf/rTn2y/grvuugtz5szBueeei9tuuw233norvvOd72idhYKZcRBFJUsTHSRRvJLobNE+9PovSlLfAhAmOwgyAyT1FuLF+m0AwCuIvJdEtkt4h/UvIlFdFQBxSbaD4A8myXsbPCSb7k0G9V+TpA6J5BqSvLeA7HtP1kb/3PlGBNfr3/TbJOIOplXkKKuiu7sb5eXl6btPNNtw880345prrjnh006ZMkXcpYaGBtx7772IRqO2Zz0KZuBARER0MigvL88YOJxIdXU1qqurR60vHR0dGDNmjNZSCQcOREREdlj/WG6Qtx89XV1d+PDDD9HV1YVkMomOjg4AwNSpU1FaWopnn30Wvb29uPDCCxEMBvH888/jgQcewC233KJ1HA4ciIiI7FBW6pZN+1F099134wc/+EH653PPPRcA8OKLL2Lu3Lnw+Xxoa2vDypUroZTC1KlT8fDDD2P58uVax+HAgYiIqABs3LjxhHs4LFy4MGPjJykOHIiIiOxw+ZbTTnHvwEEZeuH+pv4bYsRk6QSSKHoz4d7jSPd6l2ygZgj6ZwlqE4hfkyALISHI4PAN6bcxE5ITLrvGvRFBzRNBhoTkfZLW35Acy/I6U2MmIcj4AIBAWFLcQb+JqFaFoB6LU3U0xFwe4+AU9w4ciIiI3IQzDgC4ARQRERFp4IwDERGRHQpZzjjkrCd5xYEDERGRHVyqAMClCiIiItLAGQciIiI7LAuyYkUfbX/yc+3AwbT0CvRICjWZSWk6pv50k6R/otckSHcUF/uSpE5JUrQEhbE8gjaALAVW1EZwzpN+/UaWT/jmCmZUJSmckpReyXkA5AW/tAn+bVDCb+Johf6ksW9Iktar3yQhKJQmuV6TMQcnzrlUAYBLFURERKTBtTMORERErsIZBwAcOBAREdnDnSMBcKmCiIiINHDGgYiIyAalLKgsSmNn09ZNXDtwUNCMPJdkLQiyIwBh5oJgbkd5Bf2L658IaUS3KShqZAmO5ZEUCNNvkmonOReiiHPBYSzn3tukoLBYIqrfP++I/nGMpOzdlRSsEl3jgnMn5R3WPxfxUv3zEC/RbgLfoH4bUWaOsOiZiFLZLTcwxoGIiOgUorKMcSiQgQNjHIiIiMg2zjgQERHZYVmAkUWcAmMciIiITiFcqgDApQoiIiLSUDAzDkZCUthBdiynZpssQV0HjyBKXVI3AZBFj0ui1CW1IJJB/TaAMGNG8EdEolS/jWO1FoQkmRiiWV/hBSs5lhJ8R0g+t9I/Yodr9M+Fb0j/OIbg2jMT+i/KEnzWJW2klGVBZbFUwXRMIiKiUwmXKgBwqYKIiIg0cMaBiIjIDkvJasEfViAzDhw4EBER2aEUgGzSMQtj4MClCiIiIrLNtTMOBvRmhER1HZLC6GzJoQLOjDSTRc6NaJMe/fMnCkiWZDoIM0USAf02kqhuM67fRpIpIsmyAWT9SwjqGSQF14N0pliU0SP400qUgSA434DsXIjqTggyMeLF+h/CRLH+cZKCGilSylJQWSxVKM44nFhbWxtOO+00BINBNDQ04NVXXx2tQxEREY0+ZWV/KwCjMnD48Y9/jObmZqxatQqvv/46Zs2ahQULFuDgwYOjcTgiIqJRpyyV9a0QjMrA4eGHH8by5ctx7bXXYsaMGVi3bh2Ki4vx/e9/fzQOR0RERA7JeYxDLBbD7t270dLSkr7PNE3Mnz8f7e3tRz0+Go0iGo2mf+7v7wcAWJGI1nElm8kZljDGQbLmLiDZMRHCtVIRwUlP6r2tAGRrxtIYB9FOgYL3SQmuISU5D9GPf8ixSGJRkpJvE0H/pEvMSvDZkFwPEpLrAUhlB+qShHaJ4kMEx0lKdomNpr5UnIgfSKhoVssNCUe/oEdPzgcOH3zwAZLJJGpqajLur6mpwZ/+9KejHt/a2orVq1cfdX/3qvty3TUiIipQhw4dQkVFxag8t9/vRygUwu96fpX1c4VCIfj9gv3ZXSTvWRUtLS1obm5O/9zX14fJkyejq6tr1C6Ck0E4HEZdXR26u7tRXl6e7+7kFc9FCs9DCs9DCs9DSn9/PyZNmoSqqqpRO0YwGERnZydiMcHUyxH8fj+CQWExHZfI+cBh3Lhx8Hg86O3tzbi/t7cXoVDoqMcHAgEEAkfnwFVUVJzSH4bDysvLeR7+jucihechhechhechxTRHd10pGAye9P/g50rOz7Tf78fs2bOxdevW9H2WZWHr1q1obGzM9eGIiIjIQaOyVNHc3Ixly5bh/PPPxwUXXIBHH30UQ0NDuPbaa0fjcEREROSQURk4XHHFFXj//fdx9913o6enB+eccw62bNlyVMDksQQCAaxateqYyxenEp6Hf+C5SOF5SOF5SOF5SOF5cJ6hCmUPTCIiIhp1LHJFREREtnHgQERERLZx4EBERES2ceBAREREtrlu4HCql+O+5557YBhGxm369On57taoe+mll3DZZZehtrYWhmHgmWeeyfi9Ugp33303JkyYgKKiIsyfPx/vvPNOfjo7ij7uPFxzzTVHXR8LFy7MT2dHUWtrKz796U+jrKwM48ePx+WXX449e/ZkPCYSiaCpqQljx45FaWkpli5detTGcyc7O+dh7ty5R10TX/va1/LU49Gxdu1azJw5M73ZVWNjI5577rn070+Fa8FNXDVwYDnulE996lM4cOBA+va73/0u310adUNDQ5g1axba2tqO+fuHHnoIjz32GNatW4edO3eipKQECxYsQESzGJrbfdx5AICFCxdmXB9PPfWUgz10xvbt29HU1IQdO3bg+eefRzwex6WXXoqhoaH0Y1auXIlnn30WTz/9NLZv3479+/djyZIleex17tk5DwCwfPnyjGvioYceylOPR8fEiRPx4IMPYvfu3di1axfmzZuHxYsX4+233wZwalwLrqJc5IILLlBNTU3pn5PJpKqtrVWtra157JWzVq1apWbNmpXvbuQVALV58+b0z5ZlqVAopL7zne+k7+vr61OBQEA99dRTeeihM448D0optWzZMrV48eK89CefDh48qACo7du3K6VS77/P51NPP/10+jF//OMfFQDV3t6er26OuiPPg1JK/dM//ZP65je/mb9O5cmYMWPUf/3Xf52y10I+uWbG4XA57vnz56fvO1E57kL2zjvvoLa2FlOmTMHVV1+Nrq6ufHcprzo7O9HT05NxbVRUVKChoeGUuzYAYNu2bRg/fjymTZuGG264AYcOHcp3l0Zdf38/AKQLGe3evRvxeDzjmpg+fTomTZpU0NfEkefhsB/96EcYN24czjrrLLS0tGB4eDgf3XNEMpnEpk2bMDQ0hMbGxlP2WsinvFfHPEy3HHehamhowMaNGzFt2jQcOHAAq1evxmc/+1m89dZbKCsry3f38qKnpwcAjnltHP7dqWLhwoVYsmQJ6uvrsW/fPtxxxx1YtGgR2tvb4fF48t29UWFZFm666SbMmTMHZ511FoDUNeH3+1FZWZnx2EK+Jo51HgDgy1/+MiZPnoza2lq88cYbuO2227Bnzx78z//8Tx57m3tvvvkmGhsbEYlEUFpais2bN2PGjBno6Og45a6FfHPNwIFSFi1alP7/mTNnoqGhAZMnT8ZPfvITXHfddXnsGbnBlVdemf7/s88+GzNnzsTpp5+Obdu24ZJLLsljz0ZPU1MT3nrrrVMi1udEjncerr/++vT/n3322ZgwYQIuueQS7Nu3D6effrrT3Rw106ZNQ0dHB/r7+/HTn/4Uy5Ytw/bt2/PdrVOSa5YqdMtxnyoqKyvxyU9+Env37s13V/Lm8PvPa+NoU6ZMwbhx4wr2+lixYgV+8Ytf4MUXX8TEiRPT94dCIcRiMfT19WU8vlCvieOdh2NpaGgAgIK7Jvx+P6ZOnYrZs2ejtbUVs2bNwne/+91T7lpwA9cMHFiO+9gGBwexb98+TJgwId9dyZv6+nqEQqGMayMcDmPnzp2n9LUBAO+99x4OHTpUcNeHUgorVqzA5s2b8dvf/hb19fUZv589ezZ8Pl/GNbFnzx50dXUV1DXxcefhWDo6OgCg4K6JI1mWhWg0espcC66S7+jMj9q0aZMKBAJq48aN6g9/+IO6/vrrVWVlperp6cl31xxz8803q23btqnOzk718ssvq/nz56tx48apgwcP5rtro2pgYED9/ve/V7///e8VAPXwww+r3//+9+ovf/mLUkqpBx98UFVWVqqf/exn6o033lCLFy9W9fX1amRkJM89z60TnYeBgQF1yy23qPb2dtXZ2aleeOEFdd5556kzzjhDRSKRfHc9p2644QZVUVGhtm3bpg4cOJC+DQ8Ppx/zta99TU2aNEn99re/Vbt27VKNjY2qsbExj73OvY87D3v37lVr1qxRu3btUp2dnepnP/uZmjJlivrc5z6X557n1u233662b9+uOjs71RtvvKFuv/12ZRiG+s1vfqOUOjWuBTdx1cBBKaX+4z/+Q02aNEn5/X51wQUXqB07duS7S4664oor1IQJE5Tf71ef+MQn1BVXXKH27t2b726NuhdffFEBOOq2bNkypVQqJfOuu+5SNTU1KhAIqEsuuUTt2bMnv50eBSc6D8PDw+rSSy9V1dXVyufzqcmTJ6vly5cX5MD6WOcAgNqwYUP6MSMjI+rrX/+6GjNmjCouLlZf+MIX1IEDB/LX6VHwceehq6tLfe5zn1NVVVUqEAioqVOnqm9961uqv78/vx3PsX/9139VkydPVn6/X1VXV6tLLrkkPWhQ6tS4FtyEZbWJiIjINtfEOBAREZH7ceBAREREtnHgQERERLZx4EBERES2ceBAREREtnHgQERERLZx4EBERES2ceBAREREtnHgQERERLZx4EBERES2ceBAREREtnHgQERERLb9fzzApT7fBnhMAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "# Lets see if this is well structured\n",
    "wind_u_real = np_dataset[:61*30, 3].reshape(y_len, x_len)\n",
    "plt.pcolor(wind_u_real.T)\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a torch customized dataset\n",
    "\n",
    "from torch.utils.data import Dataset, DataLoader\n",
    "class HurricaneDataset(Dataset):\n",
    "    def __init__(self, X, y):\n",
    "        self.X = X\n",
    "        self.y = y\n",
    "    \n",
    "    def __len__(self):\n",
    "        return len(self.X)\n",
    "    \n",
    "    def __getitem__(self, idx):\n",
    "        return torch.tensor(self.X[idx], dtype = torch.float32, device = torch_dev), torch.tensor(self.y[idx], dtype = torch.float32, device = torch_dev)\n",
    "        #return list(self.X[idx]), self.y[idx]\n",
    "train_dataset = HurricaneDataset(X_train, y_train)\n",
    "test_dataset = HurricaneDataset(X_test, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Featrue map function "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pennylane as qml\n",
    "from pennylane import numpy as np\n",
    "import torch\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def feature_map(inputs, \n",
    "                circ_wires,\n",
    "                input_weights,\n",
    "                weights, \n",
    "                rot_seq = [qml.RX]):\n",
    "    \"\"\" \n",
    "    **Args**:\n",
    "        * inputs: *Sequence*[float] features from the training dataset, should be in the form (t, x, y)\n",
    "        * input_weights: *Sequence*[float] trainable parameters to linear transform the inputs (array should be twice as long as the inputs array)\n",
    "        * train_weights: *Sequence*[float] trainable parameters of basic entangling layers\n",
    "        * rotation_sequence: *Single or Sequence*[pennylane.ops.Operation]  if Single, use the same roation operation for each basic entangle embedding, if Sequence should be # of inputs -1 long.\n",
    "    \n",
    "    **Returns**:\n",
    "        * Trainable feature map quantum circuit\n",
    "    \"\"\"\n",
    "\n",
    "    n_qubits = len(circ_wires)\n",
    "\n",
    "    #################################### Quantum Circuit ####################################\n",
    "    # First layer of rotations (time)\n",
    "    k = 0   # This will keep track on the input_weights indices\n",
    "    for wire in circ_wires:\n",
    "        qml.RY(input_weights[k]*inputs[0]+input_weights[k+1], wire) # time\n",
    "        k+=2\n",
    "    \n",
    "    # If len(rot_seq) == 1 we want to multiply this so it has the disired length\n",
    "    if len(rot_seq) == 1:\n",
    "        rot_seq = rot_seq*3\n",
    "    kk = 0\n",
    "\n",
    "    # Layers of BasicEnganledLayers and rotations (t, x, y)\n",
    "    for ii in range(1, 3):\n",
    "        qml.BasicEntanglerLayers(weights=weights[kk:kk+n_qubits].unsqueeze(dim = 0), wires = circ_wires, rotation = rot_seq[ii-1])\n",
    "        for wire in circ_wires:\n",
    "            qml.RY(input_weights[k]*inputs[ii] + input_weights[k+1], wire)\n",
    "            k+=2\n",
    "    \n",
    "    # Add a z constant feature for the derivatives\n",
    "    qml.BasicEntanglerLayers(weights=weights[kk:kk+n_qubits].unsqueeze(dim = 0), wires = circ_wires, rotation = rot_seq[-1])\n",
    "    for wire in circ_wires:\n",
    "        qml.RY(input_weights[k]*0.5 + input_weights[k+1], wire)\n",
    "        k+=2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create full circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def full_circuit(qnode, \n",
    "                 inputs,\n",
    "                FM_input_weights,\n",
    "                FM_train_weights,\n",
    "                NN_weights,\n",
    "                n_layers,\n",
    "                rot_FM_seq = [qml.RX],\n",
    "                rot_Train_seq = [qml.RX]):\n",
    "    \"\"\" \n",
    "    **Args**:\n",
    "        * qnode: *Single*[pennylane.workflow.qnode.QNode] Qnode where circuit runs\n",
    "        * inputs: *Sequence*[float] features from the training dataset, should be in the form (t, x, y)\n",
    "        * FM_input_weights: *Sequence*[float] trainable parameters to linear transform the inputs (array should be twice as long as the inputs array)\n",
    "        * FM_train_weights: *Sequence*[float] trainable parameters of basic entangling layers on the feature map\n",
    "        * FM_train_weights: *Sequence*[float] trainable parameters of basic entangling layers on the neural network\n",
    "        * n_layers: *Sequence*[int] number of layers of the neural network\n",
    "        * rot_FM_seq: *Single or Sequence*[pennylane.ops.Operation]  Rotations used in BasicEntangledLayers on the featuremap if Single, use the same roation operation for each basic entangle embedding, if Sequence should be # of inputs -1 long.\n",
    "        * rot_NN_seq: *Single or Sequence*[pennylane.ops.Operation]  Rotations used in BasicEntangledLayers on the Neural network if Single, use the same roation operation for each basic entangle embedding, if Sequence should be # of inputs -1 long.\n",
    "    \n",
    "    **Returns**:\n",
    "        * Expected value: Z(0)@Z(1)@...@Z(n_qubits)\n",
    "    \"\"\"\n",
    "\n",
    "    n_qubits = len(qnode.device.wires) # Numbers of qubits\n",
    "\n",
    "    # Apply featuremap\n",
    "    feature_map(inputs, \n",
    "                qnode.device.wires,\n",
    "                FM_input_weights,\n",
    "                FM_train_weights, \n",
    "                rot_seq = rot_FM_seq)\n",
    "\n",
    "    if len(rot_Train_seq) == 1:\n",
    "        rot_Train_seq = rot_Train_seq*n_layers\n",
    "\n",
    "    # Create NN layers\n",
    "    for ii in range(n_layers):\n",
    "        qml.BasicEntanglerLayers(weights=NN_weights[ii*n_qubits:(ii+1)*n_qubits].unsqueeze(0),\n",
    "                                 wires = qnode.device.wires,\n",
    "                                 rotation = rot_Train_seq[ii])\n",
    "\n",
    "    ####################################### Cost Operators ##########################################\n",
    "    CostOperator = qml.Z(0)\n",
    "    for ii in range(1, n_qubits):\n",
    "        CostOperator = CostOperator@qml.Z(ii)\n",
    "\n",
    "    return qml.expval(CostOperator)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lets draw the circuit to be certain of how is implemented"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(<Figure size 1300x500 with 1 Axes>, <Axes: >)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABSgAAAIHCAYAAACczLLlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAACK+ElEQVR4nO3de3hU5b3+/3uSwCSoAUQSBGxCQiqCVoFK6imAnOqJogJatydqA+G0gYqN+BW2VHcFUWFvsBKgGrC1lqhI1VoEBcFuDwXqKYhNIkGpQhA5BxKSzO8PfgmsTBZMJplZz5p5v67La3etTCafzJ17HvaTySyPz+fzCQAAAAAAAAAcEOP0AAAAAAAAAACiFxuUAAAAAAAAABzDBiUAAAAAAAAAx7BBCQAAAAAAAMAxbFACAAAAAAAAcAwblAAAAAAAAAAcwwYlAAAAAAAAAMewQQkAAAAAAADAMWxQAgAAAAAAAHAMG5QAAAAAAAAAHMMGJQAAAAAAAADHsEEJAAAAAAAAwDFsUAIAAAAAAABwDBuUAAAAAAAAABzDBiUAAAAAAAAAx7BBCQAAAAAAAMAxbFACAAAAAAAAcAwblAAAAAAAAAAcwwYlAAAAAAAAAMewQQkAAAAAAADAMWxQAgAAAAAAAHAMG5QAAAAAAAAAHMMGJQAAAAAAAADHsEEJAAAAAAAAwDFsUAIAAAAAAABwDBuUAAAAAAAAABzDBiUAAAAAAAAAx7BBCQAAAAAAAMAxbFACAAAAAAAAcAwblAAAAAAAAAAcwwYlAAAAAAAAAMewQQkAAAAAAADAMWxQAgAAAAAAAHAMG5QAAAAAAAAAHMMGJQAAAAAAAADHsEEJAAAAAAAAwDFsUAIAAAAAAABwDBuUAAAAAAAAABzDBiUAAAAAAAAAx7BBCQAAAAAAAMAxbFACAAAAAAAAcAwblAAAAAAAAAAcwwYlAAAAAAAAAMewQQkAAAAAAADAMWxQAgAAAAAAAHAMG5QAAAAAAAAAHMMGJQAAAAAAAADHsEEJAAAAAAAAwDFsUAIAAAAAAABwDBuUAAAAAAAAABzDBiUAAAAAAAAAx7BBCQAAAAAAAMAxbFACAAAAAAAAcAwblAAAAAAAAAAcwwYlAAAAAAAAAMewQQkAAAAAAADAMWxQAgAAAAAAAHAMG5QAAAAAAAAAHMMGJQAAAAAAAADHsEEJAAAAAAAAwDFsUAIAAAAAAABwDBuUAAAAAAAAABzDBiUAAAAAAAAAx7BBCQAAAAAAAMAxbFACAAAAAAAAcAwblAAAAAAAAAAcwwYlAAAAAAAAAMewQQkAAAAAAADAMXFODwCEQ0VFhb788ksVFRWpuLhY3333nY4ePaqKigqnR7Pwer2Kj4/XOeeco65duyojI0NpaWnyer1Oj9asyMMs5GEW8jALeZiFPMxCHmYhD7OQh1nIwyzkgYawQYmIVFVVpXfeeUcFBQVatWqVtm/fLp/P5/RYQfF4PEpJSdGQIUM0YsQI9e3bV3Fx7qoueZiFPMxCHmYhD7OQh1nIwyzkYRbyMAt5mIU8EBAfEEGKiop8OTk5vvbt2/skReR/7du39+Xk5PiKi4udfrhPizzMQh5mIQ+zkIdZyMMs5GEW8jALeZiFPMxCHmgMNigREQ4fPuybPn26r2XLlo4/QYXrP6/X65s+fbrv8OHDTj/8fsjDLORhFvIwC3mYhTzMQh5mIQ+zkIdZyMMs5IFgeHw+l76uFvj/vfnmmxo9erS2b99+2tsmJyera9euOu+889SqVSu1bNlSMTFmXCuqpqZGlZWVKi8v19dff63i4mLt2rXrtJ+XmpqqvLw8DR48OAxTnh55kEcokAd5hAJ5kEcokAd5hAJ5kEcokAd5hAJ5mJWHqzi9Qwo0xcsvv+yLi4uz/S1G7969fbNmzfJt3LjRt3//fqfHbbT9+/f7Nm7c6Js1a5avV69ett9nXFyc7+WXX3Z6XPIgj7AiD7OQh1nIwyzkYRbyMAt5mIU8zEIeZnFbHm7DBiVcy+7JLzY21jd58mRfSUmJ0yM2u+LiYt/kyZN9sbGxxj0Jkgd5OI08zEIeZiEPs5CHWcjDLORhFvIwC3mYxeQ83IgNSrjSqlWrGnzyy8rK8n366adOjxdyn376qS8rK6vBJ8FVq1aFfR7yIA+TkIdZyMMs5GEW8jALeZiFPMxCHmYhj/Bbv3697cdMy8Ot2KCE6xw+fNiXkpLiV/7s7GxfdXW10+OFTXV1tS87O9vvcUhNTfWVl5eHbQ7yOI48zEIeZiEPs5CHWcjDLORhFvIwC3mYhTzCo7q62pebm+v75S9/edrbmZCHm7FBCdeZPn16xD75NZbdk+CMGTPCNgN5nEAeZiEPs5CHWcjDLORhFvIwC3mYhTzMQh6hVV5e7hs+fLhPku/3v//9aW9vQh5uxgYlXKWoqMjn9XotZc/KyoqIJ79gVVdX+6666irLY+L1en3FxcUh/9rk4Y88zEIeZiEPs5CHWcjDLORhFvIwC3mYhTxCY+fOnb4+ffrUfV9btmwJ6POczMPt2KCEq+Tk5FiKHhsbGzHvadEUn3zyid8b8+bk5IT865JHw8jDLORhFvIwC3mYhTzMQh5mIQ+zkIdZyKN5ffbZZ5Y/W2/Tpk2jNl2dysPt2KCEaxw7dszXvn17S8knT57s9FjGmDx5suWxSUpK8h07dixkX488To08zEIeZiEPs5CHWcjDLORhFvIwC3mYhTyax+rVq32JiYmW72vIkCGNvp9w5xEJ2KCEa6xZs8ZScEm+kpISp8cyRnFxsd/j89Zbb4Xs65HHqZGHWcjDLORhFvIwC3mYhTzMQh5mIQ+zkEfTLV68uMGrkT/00EONvq9w5xEJYgS4REFBgeW4d+/eSktLc2ga86Snp6tXr16Wc/Ufs+ZEHqdGHmYhD7OQh1nIwyzkYRbyMAt5mIU8zEIewaupqVFubq6ys7NVVVXl9/HLLrus0fcZ7jwiARuUcI1Vq1ZZjkeMGOHQJOaq/5jUf8yaE3mcHnmYhTzMQh5mIQ+zkIdZyMMs5GEW8jALeTTekSNHNHLkSD322GMNftzj8SgzMzOo+w5nHpGADUq4QkVFhbZv3245N3DgQIemMdegQYMsx9u3b1dFRUWzfx3yCAx5mIU8zEIeZiEPs5CHWcjDLORhFvIwC3k0zq5du9SvXz+99NJLtrfp3r27WrduHdT9hyuPSMEGJVzhyy+/lM/ns5z74Q9/6NA05srIyLAc19TUaNu2bc3+dcgjMORhFvIwC3mYhTzMQh5mIQ+zkIdZyMMs5BG4wsJCZWZm6sMPPzzl7YL58+5a4cojUrBBCVcoKiqyHCcnJ+uss85yaBpzJSYmKikpyXKu/mPXHMgjMORhFvIwC3mYhTzMQh5mIQ+zkIdZyMMs5BGY1atX6/LLL/d7FWjr1q0VE2PdJmvKBmW48ogUbFDCFYqLiy3HXbt2dWgS89X/LU0ongDJI3DkYRbyMAt5mIU8zEIeZiEPs5CHWcjDLORxaosXL9Y111yjAwcOWM6npqZq+fLlqqmpsZxvygalFJ48IgUblHCF7777znJ83nnnOTSJ+Tp37mw53rNnT7N/DfIIHHmYhTzMQh5mIQ+zkIdZyMMs5GEW8jALeTSs9krdo0ePVnV1teVjmZmZev/997V//37L+TZt2uj8889v0tcNRx6RIs7pAYBAHD161HLcqlUrhyYxX/3Hpv5j1xzII3DkYRbyMAt5mIU8zEIeZiEPs5CHWcjDLOThr7y8XHfeeWeDF8MZMWKEli5dqoSEBL333nuWj2VmZvr9yXdjhSOPSMEGJVyh/pWuWrZs6dAk5vN6vZbjUDwBkkfgyMMs5GEW8jALeZiFPMxCHmYhD7OQh1nIw+r777/XNddcY3sxnA4dOighIUGS/DYom/rn3VJ48ogU/Ik3XKmpv8WIZE48NuRhjzzMQh5mIQ+zkIdZyMMs5GEW8jALeZiFPKzatGmjnJwctW/fvsGPz58/Xx6PR6tXr9bmzZstH2uODUqTHxvT8EgBAAAAAAAg4sTExGjUqFH64osvNH78eNsNw8GDB6uysrLu2OPxKDMzM1xjQmxQAgAAAAAAIIK1bdtWCxYs0MaNGwO6fffu3dW6desQT4WTsUEJAAAAAACAiFf/St12muPPu9E4bFACAAAAAAAgovl8PvXv3z+g27JBGX5cxRsAAAAAAACuc/ToUW3ZskWFhYU6ePBg3VWy4+PjddZZZ6lHjx7q0aOHvF6vJk2a1OB9bN68WePHj7dcxZsNyvBjgxIAAAAAAADGO3jwoAoKCvTWW2/p448/1tatW1VdXX3Kz4mNjVVGRoa2bt3q97Hi4mKlp6fr3Xff1dKlS5Wbm6tjx47p/PPPD9W3ABtsUAIAAAAAAMBIPp9P69ev17PPPquCggKVl5c36vOrq6sb3Jy8+OKLlZaWJunE1b6HDRumV155xfZq3wgdHnEAAAAAAAAYZ82aNerRo4f69eunpUuXNnpz8lQ+/vhj9ejRQ2vWrKk717ZtW40aNarZvgYCFxUblP/4xz907bXXqk2bNjrjjDP0k5/8RMuXL3d6LLhAaWmpPB6P5b8WLVqoU6dOGjlypDZu3ChJmjdvnjwezymfyNatW6eYmBhdeumlqqqqCte3AIQM/QDs0Q/AHv0AAJzON998o1tvvVWDBg3S559/ftrbp6enKysrS4MGDdKgQYOUlZWl9PT0037e559/rkGDBunnP/+5vv322+YYHUGK+D/xXrt2rYYMGaL4+HjdeuutOuuss/TSSy/plltu0ddff617773X6RHhAunp6br99tslSYcPH9amTZtUUFCgV155RWvWrNGkSZO0cuVK5efn66abbtINN9xg+fxDhw5p1KhR8nq9WrZsmeLiIr56iCL0A7BHPwB79AMAUJ/P59OiRYt033336eDBgw3eJjY2Vtddd51++tOf6uKLL9ZFF12ks846q8Hb5uTkKC8v77Rf94UXXtDrr7+uOXPmaMyYMU36HhCciF7Fq6qqlJ2drZiYGK1fv16XXHKJJGnGjBnq06ePHnjgAQ0fPlwpKSnODgrjde3aVQ899JDl3KxZszRt2jRNnz5d77zzjvLz8/WjH/1I2dnZKiwsVLt27epue++996q0tFRz587VBRdcEObpgdCiH4A9+gHYox8AgJNVVVVp8uTJeuqppxr8eLdu3XTPPffo9ttvV4cOHU57f/v27WtwczI3N1crV670e1/KgwcPKicnR5999pnmzp3LL77CLKL/xPvtt99WSUmJbrvttrrNSUlq3bq1HnjgAVVWVmrp0qXODQhXu+eeeyRJmzZtkiSlpKRo3rx52rVrl8aOHVt3u1WrVmnRokXq37+/Jk2a5MisQLjRD8Ae/QDs0Q8AiE6VlZUaOXJkg5uTbdu2VV5engoLCzV16tSANicl6dxzz/U7l5WVpVmzZqmwsFALFy5UmzZt/G6zYMECjRw5UpWVlY3+PhC8iN6gXLdunSRp8ODBfh8bMmSIJOmdd94J50iIQCf/VmXUqFEaOnSoCgoK9Kc//Un79u3TL3/5SyUmJurZZ5+Vx+NxcFIg/OgHYI9+APboBwBEj2PHjmn48OFasWKF38fuuusubd26VaNHj27UlbXXrVuno0eP+p2vvSBOTEyMxowZoy+++EJ33XWX3+1WrFih4cOH69ixY434TtAUEb1BWVRUJEnKyMjw+1iHDh105pln1t0GaKwlS5ZIkq688krL+UWLFumcc87R+PHjdffdd2vHjh2aN28ebyWAqEI/AHv0A7BHPwAg+kyfPl2vvvqq5VzLli31/PPPKz8/X0lJSY26P5/Pp/79+/udX7JkiVq0aGE5l5SUpPz8fP3xj39Uy5YtLR979dVXNWPGjEZ9bQQvov+gfv/+/ZKO/0l3QxITE+tu4xY+n0/l5eVOjxF2Tv/Wori4uO49kmrfxH3t2rVKTk7WnDlzLLdNTk5WXl6ebr75Zq1cuVJDhw495dUnQ+3YsWM6fPhws98nghOJedAP//tEcCIxD/rhf58ITiTmQT/87xPBIQ+zkIdZTM9j9erVmj17tuXcGWecoddee039+vUL6j7t3vqj9m1EGnLbbbepY8eOuv766y2P1+zZszVgwAANHDgwqFnQCL4INmjQIJ8kX1FRUYMf79ixoy8xMTHMUzXNoUOHfJKi/r9x48aF5fHetm2b7QwdOnSw/dny+Xy+Pn36+CT5tmzZEpZZa40bNy5i83CjSM6DfpiVhxtFch70w6w83CiS86AfZuXhRuRhFvIwi5vy2LVrly85OdlyXy1atPBt2LAh6O9/7969Dc5YXFwc0Odv2LDB16JFC7+1adeuXUHNUz8PfnbtRfSfeNe+ctLuVZIHDhywfXUlcLIhQ4bI5/PJ5/OprKxMc+bMUVlZmYYOHapDhw41+DkJCQmW/wtEKvoB2KMfgD36AQDRLScnR7t27bKce+yxx/ze5qMx7C6Mk56eHtDnX3nllX6v6Ny5c6flQm0IjYj+E+/a954sKipS7969LR/buXOnDh06pD59+jgxWtBatWpl+w+2SDZlyhQtXrzY6TEkSe3bt9fUqVO1f/9+PfLII3rwwQc1b948p8eylZ2drblz5zbrfZqUh9tEeh70w6w83CbS86AfZuXhNpGeB/0wKw+3IQ+zkIdZTM2jsLDQ76I41157re2fZwfidBfGCdSkSZO0evVqvfHGG3XnXn75ZW3ZskXdu3cPej6cWkRvUPbt21ePPvqo3nzzTd16662Wj61ataruNm7i8Xh0xhlnOD1G2NV/I1sTPPDAA3rmmWf0u9/9TpMnT1ZqaqrTIzWoRYsWzf4zY2IebhEtedAPBCNa8qAfCEa05EE/EAzyMAt5mMXUPJ588knL8dlnn638/Hx5PJ6g7s/XiAvjnE5MTIzy8/PVrVs37d271zJz7cXc0Pwi+k+8BwwYoLS0ND3//PP66KOP6s7v379fv/3tb9WyZUvdeeedzg0IV0tISFBubq6OHTumhx9+2OlxAKPQD8Ae/QDs0Q8AiHw7d+7UH/7wB8u5cePGqX379kHfZzAXxjmVpKQkjRs3znLuueee086dO4O6P5xeRG9QxsXFacmSJaqpqVFWVpZGjx6te++9VxdffLH+9a9/6be//a2xv5WFO4wePVodO3bUsmXLVFJS4vQ4gFHoB2CPfgD26AcARLa8vDxVVlbWHbds2VITJkwI+v727dun+fPn+50vLi4O+j4lacKECWrZsmXdcWVlpfLy8pp0n7AX0RuUktS/f3+9++67uuKKK/TnP/9ZTz/9tJKTk/XCCy/o3nvvdXo8uFx8fLymTZumqqoqzZw50+lxAKPQD8Ae/QDs0Q8AiGxvv/225fiOO+5QcnJy0PfX1Avj2OnQoYNuv/12y7m1a9c26T5hL6Lfg7JWnz59LG9uCgQqNTVVPp/vlLeZMGFCg7/tWbduXYimAsxAPwB79AOwRz8AIHpVVVVp48aNlnPDhg0L+v6a68I4doYNG6Znnnmm7njjxo2qqqpSXFxUbKeFVcS/ghIAAAAAAADO++yzz1ReXm45l5mZGdR9NeeFcezUn+3w4cMqLCxslvuGFRuUAAAAAAAACLn333/fcpyenh70xXGa+8I4DUlKSlJaWprlXP3vAc2DDUoAAAAAAACE3FdffWU57tmzZ1D3E6oL4zSk/oz1vwc0DzYoAQAAAAAAEHJHjhyxHLdu3Tqo+wnVhXEaUn/G+t8Dmgfv6gkAAAAAAICQu/fee3XrrbfqyJEjOnLkiDp27Njo+wj1hXHqmzhxooYPH66EhAQlJCSoU6dOIfk60Y4NSgAAAAAAAIRc586d1blz56A/PxwXxqnvkksu0SWXXBKS+8YJ/Ik3AAAAAAAAjBeOC+PAGWxQAgAAAAAAwGjhvDAOwo8NSgAAAAAAABgtnBfGQfixQQkAAAAAAABjhfvCOAg/NijhSjU1NU6PYCweG/AzYI/HBvwM2OOxAT8D9nhsAMBeqJ8jnbgwTnNh/QgcG5RwBa/XazmurKx0aBLzVVRUWI7j4+MdmgThQj8CRz+iD/0IHP2IPvQjcPQDAE4I9/rh5gvjsH4Ejg1KuEL9EpeXlzs0ifnqPzY8AUY++hE4+hF96Efg6Ef0oR+Box8AcEI41w+3XxiH9SNwbFDCFc455xzL8ddff+3QJObbsWOH5bhdu3YOTYJwoR+Box/Rh34Ejn5EH/oROPoBACeEc/1w+4VxWD8CxwYlXKFr166WY7f8tsQJRUVFluOMjAyHJkG40I/A0Y/oQz8CRz+iD/0IHP0AgBPCtX5s2bLF9RfGYf0IHBuUcIX6Jd61a5cOHDjg0DTmOnDggMrKyizneAKMfPQjMPQjOtGPwNCP6EQ/AkM/AMAqHOuHz+fThAkT/M674cI4tVg/GocNSrhCWlqaPB6P5Vz930TA/zGJiYlRly5dHJoG4UI/AkM/ohP9CAz9iE70IzD0AwCswrF+VFZWqnv37oqJObFtdcUVV7jiwji1WD8ahw1KuILX61VKSorlnJte1h0uq1evthynpKT4XWENkYd+BIZ+RCf6ERj6EZ3oR2DoBwBYhWP98Hq9WrBggTZu3KjLLrtMiYmJevHFF5v1a4Qa60fjsEEJ1xgyZIjluKCgwKFJzFX/Man/mCFy0Y/Tox/Ri36cHv2IXvTj9OgHAPgL1/rRs2dPvfvuu3rvvffUoUOHkHyNUGH9aBw2KOEaI0aMsBxv2rRJX375pUPTmKekpESbN2+2nKv/mCFy0Y9Tox/RjX6cGv2IbvTj1OgHADQsnOtHTEyMunfvHpL7DhXWj8ZjgxKu0bdvX7Vv395ybv78+Q5NY54FCxZYjpOSkpSVleXQNAg3+nFq9CO60Y9Tox/RjX6cGv0AgIaxfpwa60fjsUEJ14iLi9PNN99sOTd//nx99tlnDk1kjk8//dRvMbjpppsUFxfn0EQIN/phj36AftijH6Af9ugHANhj/bDH+hEcNijhKlOnTrW8qWx1dbXGjx+vmpoaB6dyVk1NjcaPH6/q6uq6c16vV1OnTnVwKjiBfvijH6hFP/zRD9SiH/7oBwCcHuuHP9aP4LFBCVdJT0/Xr3/9a8u59evXKycnJyqfBGtqapSTk6MNGzZYzufm5io9Pd2hqeAU+mFFP3Ay+mFFP3Ay+mFFPwAgMKwfVqwfTcMGJVzn/vvvV0pKiuXc4sWLo+5JsPbJb/HixZbzqampuv/++x2aCk6jH8fRDzSEfhxHP9AQ+nEc/QCAxmH9OI71o+nYoITrtGrVSosWLfJ7/4bFixerX79++vTTTx2aLHw+/fRT9evXz+/JLy4uTnl5eUpISHBoMjiNftAP2KMf9AP26Af9AIBgsH6wfjQXNijhSoMHD9by5cv9ngQ3bNignj17asqUKSopKXFoutApKSnRlClT1LNnT7+XjcfFxWn58uUaPHiwQ9PBFPSDfsAe/aAfsEc/6AcABIP1g/WjObBBCde68cYbG3wSrK6u1rx589S1a1f17t1bs2bN0qZNm3TgwAGHJg3egQMHtGnTJs2aNUu9e/dW165dNW/ePMsb7konnvxuvPFGhyaFaejHCfQD9dGPE+gH6qMfJ9APAAgc68cJrB/B4RrncLUbb7xRr7/+usaMGaPS0lK/j2/evFmbN2/WtGnTJElJSUnKyMhQ586d1apVK3m9XsXEmLFPX1NTo4qKCpWXl2vHjh0qKipSWVnZaT8vNTVVeXl5/GYGfugH/YA9+kE/YI9+0A8ACAbrB+tHU7BBCdcbPHiwtmzZolmzZmn27NmqqKiwvW1ZWVlATypu4PV6lZubq/vvv5/3tIAt+kE/YI9+0A/Yox/0AwCCwfrB+hEsM7amgSZKSEjQzJkzVVhYqJycHLVv397pkUImKSlJOTk5Kiws1MyZM3nyw2nRD8Ae/QDs0Q8AQDBYPxAMNigRUdLT0/X000/rm2++0VtvvaWcnBx16dJFHo/H6dGC5vF41KVLF+Xk5Oitt97Sv//9bz399NNKT093ejS4DP0A7NEPwB79AAAEg/UDjcGfeCMixcXF6eqrr9bVV18tSaqoqNC2bdtUVFSkoqIi7dmzR0ePHtXRo0cdntQqPj5e8fHxateunTIyMpSRkaEuXbrI6/U6PRoiCP0A7NEPwB79AAAEg/UDgWCDElHB6/WqW7du6tatm9OjAMahH4A9+gHYox8AgGCwfqAh/Ik3AAAAAAAAAMewQQkAAAAAAADAMWxQAgAAAAAAAHAMG5QAAAAAAAAAHMMGJQAAAAAAAADHsEEJAAAAAAAAwDFsUAIAAAAAAABwDBuUAAAAAAAAABwT5/QAABAKpaWlWr16teLj43XjjTfqzDPPdHokwBj0A7BHPwB79MMs5GEW8gCahldQAnC13/72t+rSpYv27t1bd27dunW68MILlZOTo7vvvlu9evXS999/7+CUgDPoB2CPfgD26IdZyMMs5AGEBhuUAFztlVdeUWpqqtq2bVt3Ljc3VzU1NZo5c6bGjh2r4uJizZs3z7khAYfQD8Ae/QDs0Q+zkIdZyAMIDTYoAbhaaWmpunfvXnf8zTff6B//+IfGjx+vBx98UAsWLNCAAQO0YsUKB6cEnEE/AHv0A7BHP8xCHmYhDyA02KAE4GoHDhxQmzZt6o7Xr18vj8ejG264oe5cr1699NVXXzkwHeAs+gHYox+APfphFvIwC3kAocEGJQBXS05O1vbt2+uOV69eLa/Xq8zMzLpzR48elcfjcWI8wFH0A7BHPwB79MMs5GEW8gBCg6t4A3C1Sy+9VCtXrtRrr72m+Ph4/fnPf1b//v3l9XrrbrNt2zZ17NjRwSkBZ9APwB79AOzRD7OQh1nIAwgNXkEJwNUeeOABVVVV6Wc/+5mGDBmio0eP6oEHHqj7eEVFhdavX2/5jSYQLegHYI9+APboh1nIwyzkAYQGr6AE4Gq9evXS+++/r+eee06SNHLkSPXp06fu4//85z/Vv39/3XbbbU6NCDiGfgD26Adgj36YhTzMQh5AaLBBCcDVvvrqKyUnJ+vxxx9v8OM/+clPuIIeohb9AOzRD8Ae/TALeZiFPIDQ4E+8Abhaly5d9P/+3/9zegzASPQDsEc/AHv0wyzkYRbyAEKDV1AiKlRUVOjLL79UUVGRiouL9d133+no0aOqqKhwejQLr9er+Ph4nXPOOeratasyMjKUlpZmecNlWLVt21Znn32202O4Gv2IXPSj6ehH5KIfTUc/Ihf9MAt5mIU8mo71Aw1hgxIRqaqqSu+8844KCgq0atUqbd++XT6fz+mxguLxeJSSkqIhQ4ZoxIgR6tu3r+LiqG6tq666Sh988IHTY7gK/Yge9KPx6Ef0oB+NRz+iB/0wC3mYhTwaj/UDgeBPvBFRiouLNXbsWHXs2FEDBw5UXl6eSktLXfvkJ0k+n0+lpaXKy8vTwIED1bFjR40dO1YlJSVOj2aERx99VJ988ol+85vfqKqqyulxjEY/og/9CBz9iD70I3D0I/rQD7OQh1nII3CsH2gMtnkREcrLyzVr1izNnj1blZWVTo8TUrt379bChQv17LPP6te//rXuv/9+tWrVyumxHPPYY4/poosu0syZM5WXl6eLL75YycnJ8ng8ltt5PB79/ve/d2hKZ9EP+kE/7NEP+kE/7NEP+kE/zEAeZiGP02P9iN71oynYoITrvfnmmxo9erS2b99+2tsmJyera9euOu+889SqVSu1bNlSMTFmvJC4pqZGlZWVKi8v19dff63i4mLt2rXL9vYVFRV6+OGH9dxzzykvL0+DBw8O47TmyM/Pr/vf3377rb799tsGbxet/0CgH/SjFv3wRz/oRy364Y9+0I9a9MN55GEW8jg11o/oXj+agg1KuNqKFSs0cuRI25fW9+7dWyNGjNDAgQOVkZGhxMTEME/YNAcOHFBRUZHWrFmj5cuXa/PmzX63KS0t1XXXXafly5frxhtvdGBKZ23bts3pEYxFP+gH/bBHP+gH/bBHP+gH/TALeZiFPOyxfrB+NAUblHAtuye/2NhYTZw4URMnTlRaWppD0zWPxMRE9e7dW71791Zubq5KSkq0YMECzZ8/X9XV1XW3q6qq0siRI6PySTAlJcXpEYxEP+iHRD/s0A/6IdEPO/SDfkj0wzTkYRbyaBjrB+tHU5nx2lmgkd58880Gn/yysrL00Ucfae7cua5/8mtIenq65s6dq48++khZWVmWj9U+Cb755psOTWeG77//Xl9//bXTYziKftAPO/SDftAPe/SDftAPe/TDLORhFvKIjvVjw4YNfudYP5oXG5RwnfLyco0ePdrvyS87O1tr167VhRde6NBk4XPhhRdq7dq1ys7OtpyvqqrSmDFjdOTIEYcmc8b+/fs1adIkJScnq3379urSpUvdxz744ANde+212rRpk4MThg/9oB/10Y8T6Af9qI9+nEA/6Ed99MMs5GEW8jgh0tePmpoa3X///Vq2bJntbVg/mgcblHCdWbNm+b3hbnZ2thYuXGjMG+qGQ0xMjBYuXOj3JFhaWqpZs2Y5NFX4ff/998rMzNT8+fN13nnn6YILLpDP56v7+I9+9CP9/e9/1x//+EcHpwwf+nEc/TiOfljRj+Pox3H0w4p+HEc/jqMfZiEPs5CHVSSvH0eOHNEtt9yi2bNn67LLLjvlbVk/ms7dPy2IOsXFxXrssccs57KysiLiyS8YtU+CV111leX87NmzVVJS4tBU4fXQQw/pX//6l1544QVt3LhRI0aMsHw8ISFBffv21dtvv+3QhOFDP6zoB/04Gf2woh/042T0w4p+0A/TkIdZyOOESF4/du3apX79+unFF1+UpNNuUEqsH03l7p8YRJ0nnnhCFRUVdcexsbF66qmnXP/k1xQxMTF66qmnFBsbW3euoqJCjz/+uINThc9f/vIXXX/99Ro5cqTtbVJTU7Vjx44wTuUM+uGPftCPWvTDH/2gH7Xohz/6QT9MQh5mIY8TInX9KCwsVGZmpj788ENJUps2bXT++ecH9LnRvn40hbt/ahBVqqqq9NJLL1nOTZw40fXvadEcLrroIk2cONFy7uWXX/Z7H5BI9O2336p79+6nvI3X69Xhw4fDNJEz6Ic9+kE/6Ic9+kE/6Ic9+kE/TEEeZiGP4yJ1/VizZo0uv/xyy5+tZ2ZmNmrTNZrXj6ZggxKu8c4772j37t2Wc/VLH80mTJhgOS4rK9P69esdmiZ82rVrd9qr5m3dulXnnntumCZyBv04Nfphj36AftijH6Af9qKhH6YgD7OQx3GRuH4sWbJE11xzjQ4cOGA5H8ifd9cXretHU7BBCdcoKCiwHPfu3VtpaWkOTWOe9PR09erVy3Ku/mMWibKysrRy5UrbP6HYsmWL/va3v2ngwIFhniy86Mep0Q/6cTL6YUU/6MfJ6IcV/YjufpiCPMxCHsdF0vpRU1Oj3NxcZWdnN/gqx2A2KKN1/WgKNijhGqtWrbIc138zYvg/JvUfs0j0//7f/1N1dbWuuOIK/fGPf9R3330nSfr888/1+9//XldffbW8Xq/uu+8+hycNLfpxevSDftSiH/7oB/2oRT/80Y/o7YcpyMMs5HFcpKwfR44c0ciRI/0u9lPL4/EoMzMzqPuOxvWjKeKcHgAIREVFheU9ICRF/G+kgjFo0CBNmzat7nj79u2qqKiQ1+t1cKrQuuiii/TnP/9Zd9xxh+68805Jks/n04UXXiifz6ezzjpLy5cvV0ZGhsOThg79CAz9oB+16Ic/+kE/atEPf/QjOvthEvIwC3lEzvqxa9cuDR06tO5iOA3p3r27WrduHdT9R+P60RRsUMIVvvzyS/l8Psu5H/7whw5NY676i2BNTY22bdumbt26OTRReAwdOlTbtm3T0qVL9cEHH+j7779XYmKiMjMzNWrUKJ1zzjlOjxhS9CMw9IN+1KIf/ugH/ahFP/zRj+jsh2nIwyzRnkckrB+FhYW67rrr/DZa6wvmz7trRev6ESw2KOEKRUVFluPk5GSdddZZDk1jrsTERCUlJamsrKzuXFFRUVQ8AZ599tmaMmWK02M4gn4Ehn7QD4l+2KEf9EOiH3boR3T2w0TkYZZozsPt68fq1as1fPhwv4vhtG7dWgcPHlRNTU3duaZsUEbz+hEM3oMSrlBcXGw57tq1q0OTmK/+b2nqLx6R5s9//rOOHTvm9BiOoh+Box/Rh34Ejn5EH/oROPoBJ5GHWcjD3evH4sWLG7xSd2pqqpYvX27ZnJSatkEpRd/60RRsUMIVat94uNZ5553n0CTm69y5s+V4z549Dk0SHj//+c/VqVMnTZ06VVu3bnV6HEfQj8DRj+hDPwJHP6IP/Qgc/YCTyMMs5OHO9aP2St2jR49WdXW15WOZmZl6//33tX//fsv5Nm3a6Pzzz2/S14229aMp2KCEKxw9etRy3KpVK4cmMV/9x6b+YxdpHnzwQcXHx+vJJ59Ujx49lJWVpeeeey7iv++T0Y/A0Q/6QT/s0Q/6QT/s0Y/o64dJyMMs5OG+9aO8vNz2St0jRozQ2rVrlZycrPfee8/ysczMTMXENG3bLNrWj6ZggxKuUFFRYTlu2bKlQ5OYr/4VwSL9CfA3v/mNSktL9eqrr2ro0KF6//33dffdd+vcc8/VxIkT9fHHHzs9YsjRj8DRD/pBP+zRD/pBP+zRj+jrh0nIwyzk4a714/vvv1f//v310ksvNfjxDh06KCEhQZL8Niib+ufdUvStH03BBiVcqam/xYhk0fjYxMTE6LrrrtOKFSu0Y8cO/fa3v1X79u311FNPqVevXurTp4+WLFmiQ4cOOT1qWETjz0CgovGxoR9W0fgzEKhofGzoh1U0/gwEKhofG/phFvIwC3lYmfwc2aZNG+Xk5Kh9+/YNfnz+/PnyeDxavXq1Nm/ebPlYc2xQmvzYmIZHCkBESUpKUm5urv71r39p1apVOvfcc7Vp0yaNGTNGHTt21Lhx47R9+3anxwQcQT8Ae/QDsEc/zEIeZiEPs8XExGjUqFH64osvNH78eNsNw8GDB6uysrLu2OPxKDMzM1xjQmxQAohAW7Zs0ZQpU3Tbbbfpm2++UatWrfQf//EfSk1N1cKFC9W9e3e98cYbTo8JOIJ+APboB2CPfpiFPMxCHuZr27atFixYoI0bNwZ0++7du6t169YhngonY4MSQEQ4cuSInn32WV1++eW66KKL9D//8z/q1KmTnnrqKX3zzTdatmyZPvnkE73++us688wzlZub6/TIQNjQD8Ae/QDs0Q+zkIdZyMOd6l+p205z/Hk3GifO6QEAoCk2btyoJUuW6IUXXtDBgwcVHx+vO++8Uzk5OQ2+JP+aa67RPffco8cff9yBaYHwoh+APfoB2KMfZiEPs5CHe/l8PvXv3z+g27JBGX5sUAJwtT59+kg6/hL8MWPG6M477zztS/F/8IMfqFOnTuEYD3AU/QDs0Q/AHv0wC3mYhTzMcvToUW3ZskWFhYU6ePBg3VWy4+PjddZZZ6lHjx7q0aOHvF6vJk2a1OB9bN68WePHj7dcxZsNyvBjgxKAq/3Hf/yHxowZoyuvvDLgz8nJyVFOTk4IpwLMQD8Ae/QDsEc/zEIeZiEPZx08eFAFBQV666239PHHH2vr1q2qrq4+5efExsYqIyNDW7du9ftYcXGx0tPT9e6772rp0qXKzc3VsWPHdP7554fqW4ANNigBuNpzzz3n9AiAsegHYI9+APboh1nIwyzkEX4+n0/r16/Xs88+q4KCApWXlzfq86urqxvcnLz44ouVlpYm6cTVvocNG6ZXXnnF9mrfCB0ecQAAAAAAABhnzZo16tGjh/r166elS5c2enPyVD7++GP16NFDa9asqTvXtm1bjRo1qtm+BgIX8RuUf/jDHzRmzBj9+Mc/ltfrlcfjUX5+vtNjwSVKS0vl8Xgs/7Vo0UKdOnXSyJEjtXHjRknSvHnz5PF4TvlEtm7dOsXExOjSSy9VVVVVuL6FqHDw4EE9+uijGjBggC644AKlpaX5/Zeenu70mBGHfrgD/XAG/XAH+uEM+uEO9MMs5GEW8gi9b775RrfeeqsGDRqkzz///LS3T09PV1ZWlgYNGqRBgwYpKysroAw+//xzDRo0SD//+c/17bffNsfoCFLE/4n3gw8+qO3bt+ucc87Rueeeq+3btzs9ElwoPT1dt99+uyTp8OHD2rRpkwoKCvTKK69ozZo1mjRpklauXKn8/HzddNNNuuGGGyyff+jQIY0aNUper1fLli1TXFzEVy9sdu/ercsvv1wlJSVKTEzUgQMH1Lp1a1VWVurIkSOSpI4dO6pFixYOTxq56Ie56Ifz6Ie56Ifz6Ie56IdZyMMs5BFaPp9PixYt0n333aeDBw82eJvY2Fhdd911+ulPf6qLL75YF110kc4666wGb5uTk6O8vLzTft0XXnhBr7/+uubMmaMxY8Y06XtAcCL+FZRLlixRaWmpdu/ezZvSImhdu3bVQw89pIceekhz5szR22+/rUcffVTHjh3T9OnT616Zm5iYqOzsbO3Zs8fy+ffee69KS0v16KOP6oILLnDou4hMDz30kEpKSrRs2TLt3btXkjRlyhQdPnxYH3zwgfr06aPU1FQVFhY6PGnkoh/moh/Oox/moh/Oox/moh9mIQ+zkEfoVFVVaeLEicrJyWlwc7Jbt26aM2eOduzYoZUrV2rs2LG6/PLLbTcn9+3b1+DmZG5urrp16+Z3/uDBg8rJydHEiRN5Vb4DIn6DcuDAgUpJSXF6DESge+65R5K0adMmSVJKSormzZunXbt2aezYsXW3W7VqlRYtWqT+/ftr0qRJjswayf76179qwIABuv322+XxeCwfu/TSS/XGG2+otLRUM2fOdGjC6EQ/zEA/zEQ/zEA/zEQ/zEA/zEIeZiGP0KisrNTIkSP11FNP+X2sbdu2ysvLU2FhoaZOnaoOHToEdJ/nnnuu37msrCzNmjVLhYWFWrhwodq0aeN3mwULFmjkyJGqrKxs9PeB4EX8BiUQaif/OdGoUaM0dOhQFRQU6E9/+pP27dunX/7yl0pMTNSzzz7rt4Ch6b799lv17Nmz7jg2NrbuTyuk44vZNddco+XLlzsxXtSjH86iH2ajH86iH2ajH86iH2YhD7OQR/M7duyYhg8frhUrVvh97K677tLWrVs1evToRl1Ze926dTp69Kjf+doL4sTExGjMmDH64osvdNddd/ndbsWKFRo+fLiOHTvWiO8ETcEGJRCkJUuWSJKuvPJKy/lFixbpnHPO0fjx43X33Xdrx44dmjdvHq/kDZHWrVtbFo22bdtqx44dltskJiZq165d4R4tqtEPM9APM9EPM9APM9EPM9APs5CHWcij+U2fPl2vvvqq5VzLli31/PPPKz8/X0lJSY26P5/Pp/79+/udX7Jkid97gyYlJSk/P19//OMf1bJlS8vHXn31Vc2YMaNRXxvB452kXcbn86m8vNzpMcLO6d9aFBcX66GHHpJ04k3c165dq+TkZM2ZM8dy2+TkZOXl5enmm2/WypUrNXTo0FNefTLUjh07psOHDzf7fZoiLS1NpaWldcc9e/bU6tWrtWfPHrVr105HjhzRq6++qh/84AfODXmSSMyDfvjfpynoh/N50A//+zQF/XA+D/rhf5+moB/k0RTkQR6NsXr1as2ePdty7owzztBrr72mfv36BXWfdm/9Ufs2Ig257bbb1LFjR11//fWWx2v27NkaMGCABg4cGNQsaARfFHn00Ud9knzPPvus06ME7dChQz5JUf/fuHHjwvJ4b9u2zXaGDh06+IqKimw/t0+fPj5Jvi1btoRl1lrjxo2L2DwaMmPGDN9ZZ53lO3z4sM/n8/leeukln8fj8XXq1Mk3fPhwX1pami8mJsb36KOPOjJfJOdBP8zKoyH0g340RiTn0RD6QT8aI5LzaAj9II/GIA/yCDaPXbt2+ZKTky331aJFC9+GDRuC/v737t3b4IzFxcUBff6GDRt8LVq08Fubdu3aFdQ89fNw8mfXdPyJNxCAIUOGyOfzyefzqaysTHPmzFFZWZmGDh2qQ4cONfg5CQkJlv+L0MjJydHixYvrXll80003ac6cOTp8+LBeeukl7dy5U7/61a80depUhyeNXPTDXPTDefTDXPTDefTDXPTDLORhFvJoPjk5OX5/Cv/YY4/5vc1HY9hdGCc9PT2gz7/yyiv9XtG5c+dOy4XaEBr8ibfLtGrVyvYfbJFsypQpWrx4sdNjSJLat2+vqVOnav/+/XrkkUf04IMPat68eU6PZSs7O1tz585t1vs0KY9zzz1Xt9xyi+Xcvffeq8mTJ+u7775TUlKSUW+eH+l50A+z8qAfZuVBP8zKg36YlQf9MCsP+kEeTUEe5BGIwsJCv4viXHvttbZ/nh2I010YJ1CTJk3S6tWr9cYbb9Sde/nll7VlyxZ179496PlwamxQuozH49EZZ5zh9BhhV/+NbE3wwAMP6JlnntHvfvc7TZ48WampqU6P1KAWLVo0+8+MiXnUFxsbq+TkZEnS73//e/3973/XM8884/BU0ZMH/TAb/XAW/TAb/XAW/TAb/TALeZiFPBrnySeftByfffbZys/PD3pz19eIC+OcTkxMjPLz89WtWzft3bvXMnPtxdzQ/PgTbyBICQkJys3N1bFjx/Twww87PQ5O4d1339XSpUudHiOq0A/3oB/hRz/cg36EH/1wD/phFvIwC3mc2s6dO/WHP/zBcm7cuHFq37590PcZzIVxTiUpKUnjxo2znHvuuee0c+fOoO4PpxfxG5RLlizR3XffrbvvvlsFBQV+59j9RlOMHj1aHTt21LJly1RSUuL0OIBR6Adgj34A9ugHAES2vLw8VVZW1h23bNlSEyZMCPr+9u3bp/nz5/udLy4uDvo+JWnChAlq2bJl3XFlZaXy8vKadJ+wF/EblLW/uVi6dKk2b94sSfr73/9ed+7dd991eEK4WXx8vKZNm6aqqirNnDnT6XEAo9APwB79AOzRDwCIbG+//bbl+I477qj78/hgNPXCOHY6dOig22+/3XJu7dq1TbpP2Iv496DMz89Xfn6+02PApVJTU+Xz+U55mwkTJjT4255169aFaCrADPQDsEc/AHv0AwCiV1VVlTZu3Gg5N2zYsKDvr7kujGNn2LBhlvcS3bhxo6qqqhQXF/HbaWEX8a+gBAAAAAAAgPM+++wzlZeXW85lZmYGdV/NeWEcO/VnO3z4sAoLC5vlvmHFBiUAAAAAAABC7v3337ccp6enB31xnOa+ME5DkpKSlJaWZjlX/3tA8+A1qQBcp/7V1E7nvffeC9EkgHnoB2CPfgD26IdZyMMs5NF8vvrqK8txz549g7qfUF0YpyE9e/bUl19+WXdc/3tA82CDEoDrLFy4sNGf4/F4QjAJYB76AdijH4A9+mEW8jALeTSfI0eOWI5bt24d1P2E6sI4Dak/Y/3vAc2DDUoArsOV0wB79AOwRz8Ae/TDLORhFvJoPvfee69uvfVWHTlyREeOHFHHjh0bfR+hvjBOfRMnTtTw4cOVkJCghIQEderUKSRfJ9qxQQnAdfr27ev0CICx6Adgj34A9uiHWcjDLOTRfDp37qzOnTsH/fnhuDBOfZdccokuueSSkNw3TuAiOQBc7eqrr9b06dOdHgMwEv0A7NEPwB79MAt5mIU8nBWOC+PAGWxQAnC1Dz74QNXV1U6PARiJfgD26Adgj36YhTzMQh7OCeeFcRB+bFACcLVu3bpp+/btTo8BGIl+APboB2CPfpiFPMxCHs4J54VxEH5sUAJwtYkTJ2rlypXasmWL06MAxqEfgD36AdijH2YhD7OQhzPCfWEchB8XyYEr1dTUOD2CsaLtsUlLS1O/fv30k5/8RGPGjNGll16q5ORkeTwev9tmZWU5MGH4RdvPQGNE22NDP/xF289AY0TbY0M//EXbz0BjRNtjQz/MQh5mIQ9/oX6OdOLCOM0l2taPpmCDEq7g9Xotx5WVlQ5NYr6KigrLcXx8vEOThEe/fv3k8Xjk8/n0xBNPNPgPg1qR+l4x9CNw9IN+0A979IN+0A979CP6+mES8jALeYR//XDzhXGibf1oCjYo4Qr1S1xeXu7QJOar/9hE+hPgjBkzTvmPgmhAPwJHP6IP/Qgc/Yg+9CNw9ANOIg+zkEd41w+3Xxgn2taPpmCDEq5wzjnnWI6//vprhyYx344dOyzH7dq1c2iS8HjooYecHsFx9CNw9CP60I/A0Y/oQz8CRz/gJPIwC3mEd/1w+4Vxom39aAoukgNX6Nq1q+XYLb8tcUJRUZHlOCMjw6FJEC70I3D0I/rQj8DRj+hDPwJHPwDghHCtH1u2bHH9hXFYPwLHKyjhCvVLvGvXLh04cECJiYkOTWSmAwcOqKyszHIuWp4A//nPf+pPf/qTtm7dqvLy8rpFa/v27frggw80cOBAnX322Q5PGRr0IzD0g35I9MMO/aAfEv2wQz+isx8mIg+zRHMe4Vg/fD6fJkyY4HfeDRfGqRXN60cw2KCEK6SlpdW9EXGtoqIi9e7d28GpzFP/tzMxMTHq0qWLQ9OEz69//Ws98cQTdT8fJ78njM/n02233aYnnnjC9s2V3Y5+BIZ+0I9a9MMf/aAfteiHP/oRnf0wDXmYJdrzCMf6UVlZqe7du+udd96puxL2FVdc4YoL49SK1vUjWPyJN1zB6/UqJSXFcs5NL+sOl9WrV1uOU1JS/K6wFmmeffZZPf7447r++uv1ySefaNq0aZaPp6amqk+fPvrLX/7i0IShRz8CQz/oRy364Y9+0I9a9MMf/YjOfpiEPMxCHuFZP7xerxYsWKCNGzfqsssuU2Jiol588cVm/RqhFo3rR1OwQQnXGDJkiOW4oKDAoUnMVf8xqf+YRaLf/e53uuCCC/TSSy/pwgsvVMuWLf1u061bN7/fXkUa+nF69IN+1KIf/ugH/ahFP/zRj+jthynIwyzkcVy41o+ePXvq3Xff1XvvvacOHTqE5GuESjSuH03BBiVcY8SIEZbjTZs26csvv3RoGvOUlJRo8+bNlnP1H7NItGXLFg0aNEhxcfbvWJGcnOz33h+Rhn6cGv2gHyejH1b0g36cjH5Y0Y/o7ocpyMMs5HFcONePmJgYde/ePST3HSrRun40BRuUcI2+ffuqffv2lnPz5893aBrzLFiwwHKclJSkrKwsh6YJn7i4OFVWVp7yNt98843OPPPMME3kDPpxavTDHv0A/bBHP0A/7EVDP0xBHmYhj+NYP04tWtePpmCDEq4RFxenm2++2XJu/vz5+uyzzxyayByffvqp32Jw0003nfK3epHioosu0ttvv63q6uoGP157Rb1If8N/+mGPftAP+mGPftAP+mGPftAPU5CHWcjjONYPe9G8fjQFG5RwlalTp1reVLa6ulrjx4+vu6pXNKqpqdH48eMtC6TX69XUqVMdnCp8fvGLX+hf//qXcnJyVFFRYfnYgQMHdPfdd2vnzp3Kzs52aMLwoR/+6Af9qEU//NEP+lGLfvijH/TDJORhFvI4gfXDX7SvH03BBiVcJT09Xb/+9a8t59avX6+cnJyofBKsqalRTk6ONmzYYDmfm5ur9PR0h6YKr1/84he69dZb9fvf/17t27fX73//e0lSnz591KlTJ7344ou66667NHz4cIcnDT36YUU/6MfJ6IcV/aAfJ6MfVvSDfpiGPMxCHiewflixfjQNG5Rwnfvvv18pKSmWc4sXL466J8HaJ7/Fixdbzqempur+++93aCpnPP/888rLy1OXLl3073//Wz6fTxs3btQPfvADPf3003rmmWecHjFs6Mdx9OME+nEC/TiOfpxAP06gH8fRjxPoh1nIwyzkcQLrx3GsH03HBiVcp1WrVlq0aJHf+zcsXrxY/fr106effurQZOHz6aefql+/fn5PfnFxccrLy1NCQoJDkzknOztbH3/8sQ4dOqQdO3bowIEDKiws1JgxY5weLazoB/1oCP04jn7Qj4bQj+PoB/1oCP0wC3mYhTyOY/1g/WgubFDClQYPHqzly5f7PQlu2LBBPXv21JQpU1RSUuLQdKFTUlKiKVOmqGfPnn4vG4+Li9Py5cs1ePBgh6YzQ0JCgjp27BjxV807FfpBP+zQD/pBP+zRD/pBP+zRD7OQh1nIg/WD9aN5sEEJ17rxxhsbfBKsrq7WvHnz1LVrV/Xu3VuzZs3Spk2bdODAAYcmDd6BAwe0adMmzZo1S71791bXrl01b948vyvG1T753XjjjQ5NCtPQjxPoB+qjHyfQD9RHP06gHwAQONaPE1g/gsM1zuFqN954o15//XWNGTNGpaWlfh/fvHmzNm/erGnTpkmSkpKSlJGRoc6dO6tVq1byer2KiTFjn76mpkYVFRUqLy/Xjh07VFRUpLKystN+XmpqqvLy8qLmNzNpaWlBfZ7H44nI39qdCv2gH4GiH6V+H6cfkYd+BI5+0I9ARWM/woE8zEIegWP9iL71ozmxQQnXGzx4sLZs2aJZs2Zp9uzZqqiosL1tWVlZQE8qbuD1epWbm6v7778/qt7ToqamRh6Pp9Gf5/P5QjCN+egH/QgE/aAf0YB+NA79oB+BiNZ+hBp5mIU8Gof1I7rWj+bEBiUiQkJCgmbOnKk777xTjz/+uF566SXt3r3b6bFCIikpSTfddJOmTp2q9PR0p8cJu4Z+E4dTox/Rg340Hv2IHvSj8ehH9KAfZiEPs5BH47F+IBhmvHYWaCbp6el6+umn9c033+itt95STk6OunTpEtRvvEzh8XjUpUsX5eTk6K233tK///1vPf300zz5odHoB2CPfgD26AcAIBisH2gMXkGJiBQXF6err75aV199tSSpoqJC27ZtU1FRkYqKirRnzx4dPXpUR48edXhSq/j4eMXHx6tdu3bKyMhQRkaGunTpIq/X6/RoiCD0A7BHPwB79AMAEAzWDwSCDUpEBa/Xq27duqlbt25Oj4Jm9pvf/Oa0t4mJiVFiYqLOP/989evXjwWlHvoRuehH09GPyEU/mo5+RC76YRbyMAt5NB3rBxrCBiUAV3vooYcsfyJw8ptR1z/v8XjUtm1bPfnkk7rzzjvDOifgBPoB2KMfgD36YRbyMAt5AKHBe1ACcLW1a9fq+uuvl9frVXZ2tpYuXaq//e1vWrp0qbKzs+X1enXDDTfoxRdf1LRp03Ts2DH94he/0Jo1a5weHQg5+gHYox+APfphFvIwC3kAocErKAG4WlFRkd555x1t3rzZ708E7rjjDk2ePFmZmZkaOnSoHnnkEd12223q1auXnnjiCQ0cONChqYHwoB+APfoB2KMfZiEPs5AHEBq8ghKAq/3P//yPbrnlFtv3L+nWrZtuueUWzZ07V5LUvXt33XDDDfrwww/DOSbgCPoB2KMfgD36YRbyMAt5AKHBBiUAVysuLtbZZ599ytu0a9dOJSUldcfp6ek6dOhQqEcDHEc/AHv0A7BHP8xCHmYhDyA02KAE4Grt27fXG2+8YXlz6pP5fD698cYbateuXd25vXv3qnXr1uEaEXAM/QDs0Q/AHv0wC3mYhTyA0GCDEoCr3Xrrrfrkk080dOhQffLJJ5aPffLJJ/rZz36mTz/9VD//+c/rzn/44Ye64IILwj0qEHb0A7BHPwB79MMs5GEW8gBCg4vkAHC1mTNnauPGjXr99df117/+VWeccYbat2+v3bt36/Dhw/L5fMrKytLMmTMlSTt37lRqaqpGjhzp8ORA6NEPwB79AOzRD7OQh1nIAwgNXkEJwNUSEhK0Zs0aLV68WFlZWWrRooW++uortWjRQn379tXixYv19ttvKyEhQZLUoUMHrVixwvIbTSBS0Q/AHv0A7NEPs5CHWcgDCA1eQQnA9WJiYnTPPffonnvucXoUwDj0A7BHPwB79MMs5GEW8gCaH6+gBAAAAAAAAOAYXkEJICJUVVXpiy++0L59+1RdXd3gbbKyssI8FWAG+gHYox+APfphFvIwC3kAzYsNSgCu5vP5NGPGDM2fP18HDx485W3t/uEARCr6AdijH4A9+mEW8jALeQChwQYlAFd7+OGH9d///d9q06aN7rzzTnXu3FlxcTy1ARL9AE6FfgD26IdZyMMs5AGEBi0C4GrPPPOMUlJStHHjRrVr187pcQCj0A/AHv0A7NEPs5CHWcgDCA0ukgPA1Xbu3Klhw4bxjwOgAfQDsEc/AHv0wyzkYRbyAEKDDUoArtalSxcdOHDA6TEAI9EPwB79AOzRD7OQh1nIAwgNNigBuNrYsWP12muvqayszOlRAOPQD8Ae/QDs0Q+zkIdZyAMIDd6DEoCr/exnP9OGDRt0+eWXa8aMGerVq5cSExMbvO0PfvCDME8HOIt+APboB2CPfpiFPMxCHkBosEGJqFBRUaEvv/xSRUVFKi4u1nfffaejR4+qoqLC6dEsvF6v4uPjdc4556hr167KyMhQWlqavF6v06MZq0uXLvJ4PPL5fBo1apTt7Twej6qqqsI4mXvQj8hFP5qOfkQu+tF09CNy0Q+zkIdZyKPpWD/QEDYoEZGqqqr0zjvvqKCgQKtWrdL27dvl8/mcHisoHo9HKSkpGjJkiEaMGKG+ffsqLo7q1rrzzjvl8XicHsNV6Ef0oB+NRz+iB/1oPPoRPeiHWcjDLOTReKwfCASPIiJKcXGxnnjiCb300kvavXu30+M0C5/Pp9LSUuXl5SkvL0/t27fXzTffrKlTpyo9Pd3p8RyXn5/v9AiuQT+iD/0IHP2IPvQjcPQj+tAPs5CHWcgjcKwfaAwukoOIUF5erhkzZqhHjx5auHBhxDz5NWT37t1auHChevTooRkzZqi8vNzpkWA4+gHYox+APfoBAAgG6weCwSso4XpvvvmmRo8ere3bt5/2tsnJyeratavOO+88tWrVSi1btlRMjBn79DU1NaqsrFR5ebm+/vprFRcXa9euXba3r6io0MMPP6znnntOeXl5Gjx4cBinhVvQD/oBe/SDfsAe/aAfABAM1g/Wj2CxQQlXW7FihUaOHGn75sO9e/fWiBEjNHDgQGVkZNheXc1UBw4cUFFRkdasWaPly5dr8+bNfrcpLS3Vddddp+XLl+vGG290YErnHTx4UAsWLNCaNWv0zTffNPjmyh6PRyUlJQ5M5xz6QT8k+mGHftAPiX7YoR/0Q6IfpiEPs5BHw1g/WD+agg1KuJbdk19sbKwmTpyoiRMnKi0tzaHpmkdiYqJ69+6t3r17Kzc3VyUlJVqwYIHmz5+v6urquttVVVVp5MiRUfkkuHv3bl1++eUqKSlRYmKiDhw4oNatW6uyslJHjhyRJHXs2FEtWrRweNLwoh/0Q6IfdugH/ZDohx36QT8k+mEa8jALeTSM9YP1o6nMeO0s0Ehvvvlmg09+WVlZ+uijjzR37lzXP/k1JD09XXPnztVHH32krKwsy8dqnwTffPNNh6ZzxkMPPaSSkhItW7ZMe/fulSRNmTJFhw8f1gcffKA+ffooNTVVhYWFDk8aPvSDftSiH/7oB/2oRT/80Q/6UYt+mIU8zEIe/qJh/diwYYPfOdaP5sUGJVynvLxco0eP9nvyy87O1tq1a3XhhRc6NFn4XHjhhVq7dq2ys7Mt56uqqjRmzJi639xFg7/+9a8aMGCAbr/9dnk8HsvHLr30Ur3xxhsqLS3VzJkzHZowvOgH/TgZ/bCiH/TjZPTDin7Qj5PRD7OQh1nIwyrS14+amhrdf//9WrZsme1tWD+aBxuUcJ1Zs2b5veFudna2Fi5caMwb6oZDTEyMFi5c6PckWFpaqlmzZjk0Vfh9++236tmzZ91xbGysZQFo27atrrnmGi1fvtyJ8cKOfhxHP46jH1b04zj6cRz9sKIfx9GP4+iHWcjDLORhFcnrx5EjR3TLLbdo9uzZuuyyy055W9aPpnP3TwuiTnFxsR577DHLuaysrIh48gtG7ZPgVVddZTk/e/bsqHlD5tatW+vYsWN1x23bttWOHTsst0lMTDzlFdciBf2woh/042T0w4p+0I+T0Q8r+kE/TEMeZiGPEyJ5/di1a5f69eunF198UZJOu0EpsX40lbt/YhB1nnjiCcsV0mJjY/XUU0+5/smvKWJiYvTUU08pNja27lxFRYUef/xxB6cKn7S0NJWWltYd9+zZU6tXr9aePXskHf+t16uvvqof/OAHDk0YPvTDH/2gH7Xohz/6QT9q0Q9/9IN+mIQ8zEIeJ0Tq+lFYWKjMzEx9+OGHkqQ2bdro/PPPD+hzo339aAp3/9QgqlRVVemll16ynJs4caLr39OiOVx00UWaOHGi5dzLL7/s9z4gkWjw4MF66623VF5eLkkaM2aMysrKdPHFF2vEiBG68MILVVJSorvvvtvZQUOMftijH/SDftijH/SDftijH/TDFORhFvI4LlLXjzVr1ujyyy+3/Nl6ZmZmozZdo3n9aAo2KOEa77zzjnbv3m05V7/00WzChAmW47KyMq1fv96hacInJydHixcvrvsHwk033aQ5c+bo8OHDeumll7Rz50796le/0n333efwpKFFP06NftAP+mGPftAP+mGPfkR3P0xBHmYhj+Micf1YsmSJrrnmGh04cMByPpA/764vWtePpmCDEq5RUFBgOe7du7fS0tIcmsY86enp6tWrl+Vc/ccsEp177rm65ZZbdM4559Sdu/fee/Xdd9/p22+/1aFDhzRnzhzLS+wjEf04NfpBP05GP6zoB/04Gf2woh/R3Q9TkIdZyOO4SFo/ampqlJubq+zs7AZf5RjMBmW0rh9NwQYlXGPVqlWW4xEjRjg0ibnqPyb1H7NItH79en311Vd+52NjY5WcnCyPx6Ovv/464n9bRT9Oj36cQD/oR3304wT6QT/qox8nRFs/TEEeZiGP4yJl/Thy5IhGjhzpd7GfWh6PR5mZmUHddzSuH03BBiVcoaKiwvIeEJI0cOBAh6Yx16BBgyzH27dvt7xpcSTq37+/8vPzT3mbZcuWqX///uEZyAH0IzD0o2H0AxL9sEM/INEPO5HeD5OQh1nII3LWj9orddd/L82Tde/eXa1btw7q/qNx/WgKNijhCl9++aV8Pp/l3A9/+EOHpjFXRkaG5bimpkbbtm1zaJrwqP9z0ZCamhp5PJ4wTOMM+hEY+tEw+gGJftihH5Doh51I74dJyMMs5BEZ60f9K3XbCebPu2tF4/rRFGxQwhWKioosx8nJyTrrrLMcmsZciYmJSkpKspyr/9hFo6KioqB/6+UG9CMw9KNh9AMS/bBDPyDRDzuR3g+3IQ+zRHoebl8/Vq9e7Xelbklq3bq139W6m7JByfrROHFODwAEori42HLctWtXhyYxX0ZGhsrKyuqOI/EJ8Be/+IXl+JVXXlFpaanf7aqrq+ve/+Waa64J03ThRz8CRz9OoB+oj36cQD9QH/04IVr64TTyMAt5WLl5/Vi8eLHGjh2r6upqy/nU1FTl5eVpyJAhlvNN2aCUomP9aC5sUMIVvvvuO8vxeeed59Ak5uvcubPleM+ePQ5NEjonv+eLx+PRRx99pI8++qjB23o8Hl166aWaO3dueIZzAP0IHP2woh84Gf2woh84Gf2wioZ+OI08zEIeVm5cP2pqajRt2rQGL4aTmZmplStX+l3YqE2bNjr//POb9HWjYf1oLmxQwhWOHj1qOW7VqpVDk5iv/mNT/7GLBLXv2+Hz+ZSWlqbJkydr0qRJfreLjY1V27ZtdcYZZ4R7xLCiH4GjHyfQD9RHP06gH6iPfpwQLf1wGnmYhTys3LZ+lJeX684772zwYjgjRozQ0qVLlZCQoPfee8/ysczMTL8/+W6saFg/mgsblHCF+le6atmypUOTmM/r9VqOI/EJMCUlpe5/P/vss+rZs6flXLShH4GjH9GHfgSOfkQf+hE4+oFwIw+zkIeVm9aP77//Xtdcc43txXA6dOighIQESfLboGzqn3dL0bF+NBc2KOFKTf0tRiSLtsfmrrvucnoE40Tbz0BjRNtjQz/8RdvPQGNE22NDP/xF289AY0TbY0M/zEIeZiEPfyY/R7Zp00Y5OTnatm2bdu/e7ffx+fPna/78+XrzzTe1efNmy8eaY4PS5MfGNGxQAogIH374of7xj39o3759fm94LB1/L5jp06c7MBngPPoB2KMfgD36YRbyMAt5uENMTIxGjRqlYcOGafr06Xr66adVU1Pjd7vBgwdbjj0ejzIzM8M1JsQGJQCX+/777zVs2DD9/e9/l8/ns70d/0BANKIfgD36AdijH2YhD7OQhzu1bdtWCxYs0D333KNevXqd9vbdu3dX69atwzAZarFBCcDVfvWrX+ndd99Vv379dNddd6lz586Ki+OpDZDoB3Aq9AOwRz/MQh5mIQ93279/f0C3a44/70bj0CIArvbaa6+pT58+euutt+TxeJweBzAK/QDs0Q/AHv0wC3mYhTzcy+fzqX///gHdlg3K8GODEoCrHTlyRFlZWfzjAGgA/QDs0Q/AHv0wC3mYhTzMcvToUW3ZskWFhYU6ePBg3VWy4+PjddZZZ6lHjx7q0aOHvF6vJk2a1OB9bN68WePHj7dcxZsNyvBjgxKAq11yySUqLS11egzASPQDsEc/AHv0wyzkYRbycNbBgwdVUFCgt956Sx9//LG2bt3a4EWKThYbG6uMjAxt3brV72PFxcVKT0/Xu+++q6VLlyo3N1fHjh3T+eefH6pvATa43jkAV/uv//ov/eUvf9H777/v9CiAcegHYI9+APboh1nIwyzkEX4+n0/vvPOO7r77bnXo0EH33HOPnn/+eRUWFp52c1KSqqurG9ycvPjii5WWlibpxNW+v/jiCz355JOKiWG7LNx4BSUAV9u5c6euu+469e3bV//xH/+hXr16KTExscHb3nnnnWGeDnAW/QDs0Q/AHv0wC3mYhTzCa82aNfrP//xPff75581+3x9//LF69Oih//3f/9XAgQMlHb/a96hRo5r9a+H0InqD8t///rcKCgr017/+VVu3btXOnTt19tln64orrtCvf/1rZWZmOj0iDFdaWqouXbpYzsXFxSkpKanu5+jHP/6x5s2bpylTpujuu+/Ws88+2+B9rVu3TldffbV69+6t9957jyu9NZO7775bHo9HPp9P+fn5ys/P93s/GJ/PJ4/Hwz8Qmhn9MB/9cA79MB/9cA79MB/9MAt5mIU8wuObb77Rr371K/35z38O6Pbp6enq1KmTvF6vJKmiokL//ve/VVJScsrP+/zzzzVo0CDdeuutevLJJ3Xuuec2eXYEJ6JX8Pnz52v27NlKT0/X4MGD1b59exUVFemVV17RK6+8oueff1633HKL02PCBdLT03X77bdLkg4fPqxNmzapoKBAr7zyitasWaNJkyZp5cqVys/P10033aQbbrjB8vmHDh3SqFGj5PV6tWzZMv7x3Izs/h8WhA/9MBf9cB79MBf9cB79MBf9MAt5mIU8Qsvn82nRokW67777dPDgwQZvExsbq+uuu04//elPdfHFF+uiiy7SWWed1eBtc3JylJeXd9qv+8ILL+j111/XnDlzNGbMmCZ9DwhORK/iffr00bp169S3b1/L+Q0bNmjAgAEaO3ashg0bVrfDDtjp2rWrHnroIcu5WbNmadq0aZo+fbreeecd5efn60c/+pGys7NVWFiodu3a1d323nvvVWlpqebOnasLLrggzNNHtrvuusvpEaIe/TAX/XAe/TAX/XAe/TAX/TALeZiFPEKnqqpKkydP1lNPPdXgx7t166Z77rlHt99+uzp06HDa+9u3b1+Dm5O5ublauXKl3/tSHjx4UDk5Ofrss880d+5cfvEVZhH9rp833XST3+akJF111VXq37+/9u7dq08//dSByRAJ7rnnHknSpk2bJEkpKSmaN2+edu3apbFjx9bdbtWqVVq0aJH69++vSZMmOTIrEG70A7BHPwB79AMAolNlZaVGjhzZ4OZk27ZtlZeXp8LCQk2dOjWgzUlJDf65dlZWlmbNmqXCwkItXLhQbdq08bvNggULNHLkSFVWVjb6+0DwInqD8lRatGghSeyIo8lO/hkaNWqUhg4dqoKCAv3pT3/Svn379Mtf/lKJiYl69tln/d6bBMH5zW9+o/Xr11vOlZWV6ZNPPmnw9n/+85910003hWM01EM/wo9+uAf9CD/64R70I/zoh1nIwyzkEVrHjh3T8OHDtWLFCr+P3XXXXdq6datGjx7dqCtrr1u3TkePHvU7v2bNGknHr9o9ZswYffHFFw2+KnbFihUaPny4jh071ojvBE0RlRuUX331ldasWaNzzz1XF110kdPjwKWWLFkiSbryyist5xctWqRzzjlH48eP1913360dO3Zo3rx5SklJcWLMiPTQQw9p3bp1lnNPP/20evbs2eDtt27dqpUrV4ZhMtSiH86hH+ajH86hH+ajH86hH2YhD7OQR2hNnz5dr776quVcy5Yt9fzzzys/P19JSUmNuj+fz6f+/fv7nV+yZEndi9VqJSUlKT8/X3/84x/VsmVLy8deffVVzZgxo1FfG8GLupcPHjt2THfccYcqKio0e/ZsxcbGOj1So/h8PpWXlzs9Rtg5/VuL4uLiuvdIqn0T97Vr1yo5OVlz5syx3DY5OVl5eXm6+eabtXLlSg0dOlSjRo1yYOrjjh07psOHDzf7fSI4kZgH/fC/TwQnEvOgH/73ieBEYh70w/8+ERzyMAt5mMX0PFavXq3Zs2dbzp1xxhl67bXX1K9fv6Du0+6tP2rfRqQht912mzp27Kjrr7/e8njNnj1bAwYM0MCBA4OaBYGLqg3Kmpoa3X333Vq/fr2ys7N1xx13OD1So5WXl+vMM890eoyoU1JSopkzZ1rOdejQQRs2bFDXrl39bn/TTTepT58++vDDDzVr1qxwjdmgxYsXa/HixY7OgBMiMQ/6geYSiXnQDzSXSMyDfqC5kIdZyMMsJudRVlbmty/TokUL/e1vf/N7JX2g9u3bp/nz5/udLy4uPu3n9uvXT3/729909dVX123C+nw+3XHHHfr4448b/UpONE7U/Il3TU2NfvGLX+j555/X7bffroULFzo9ElxkyJAh8vl88vl8Kisr05w5c1RWVqahQ4fq0KFDDX5OQkKC5f8CkYp+APboB2CPfgBAdMvJydGuXbss5x577LGgNycl+wvjpKenB/T5V155pd8rOnfu3Gm5UBtCIypeQVlTU6NRo0Zp2bJl+vnPf678/PxGvbmqSVq1amX7D7ZINmXKFGN+69O+fXtNnTpV+/fv1yOPPKIHH3xQ8+bNc3osW9nZ2Zo7d26z3qdJebhNpOdBP8zKw20iPQ/6YVYebhPpedAPs/JwG/IwC3mYxdQ8CgsL/S6Kc+2119r+eXYgTndhnEBNmjRJq1ev1htvvFF37uWXX9aWLVvUvXv3oOfDqUX8BuXJm5O33HKLnnvuOde97+TJPB6PzjjjDKfHCLv6b2RrggceeEDPPPOMfve732ny5MlKTU11eqQGtWjRotl/ZkzMwy2iJQ/6gWBESx70A8GIljzoB4JBHmYhD7OYmseTTz5pOT777LOVn58vj8cT1P015sI4pxMTE6P8/Hx169ZNe/futcxcezE3NL+I3qCs/bPuZcuWacSIEfrDH/7g6s1JmCUhIUG5ubmaNGmSHn74Yf3+9793eqSo8tlnn2n58uWWY0kqKCiQz+fzuy3Ci344i36YjX44i36YjX44i36YhTzMQh7NZ+fOnfrDH/5gOTdu3Di1b98+6PsM5sI4p5KUlKRx48bpv//7v+vOPffcc3rkkUfUoUOHoO4TpxbRG5S/+c1vtHTpUp155pn64Q9/qEceecTvNsOGDdMll1wS/uEQEUaPHq3Zs2dr2bJleuCBBwJ+Xws03UsvvaSXXnqp7rj2HwW33nqr3219Pl/Qv4lD8OiHc+iH+eiHc+iH+eiHc+iHWcjDLOTRfPLy8lRZWVl33LJlS02YMCHo+2vKhXFOZcKECZozZ07drJWVlcrLy9N//dd/Nel+0bCI3qAsLS2VJB06dMiy632y1NRUNigRtPj4eE2bNk0TJ07UzJkztWzZMqdHigosCO5AP5xBP9yBfjiDfrgD/XAG/TALeZiFPJrX22+/bTm+4447lJycHPT9NfXCOHY6dOig22+/Xc8880zdubVr1/LzECIRvUGZn5+v/Px8p8eAi6Wmpvq9XL++CRMmNPjbnnXr1oVoKrAgmIF+mIl+mIF+mIl+mIF+mIl+mIU8zEIezaeqqkobN260nBs2bFjQ99dcF8axM2zYMMsG5caNG1VVVaW4uIjeTnOEOy9lDQAAAAAAAFf57LPPVF5ebjmXmZkZ1H0154Vx7NSf7fDhwyosLGyW+4YVG5QAXG3Lli363//9X+3evbvBj5eVlel///d/9fnnn4d5MsB59AOwRz8Ae/TDLORhFvJomvfff99ynJ6eHvTFcZr7wjgNSUpKUlpamuVc/e8BzYMNSgCuNmvWLM2ePVvt2rVr8OPt2rXTnDlz9Nhjj4V5MsB59AOwRz8Ae/TDLORhFvJomq+++spy3LNnz6DuJ1QXxmlI/Rnrfw9oHmxQAnC1DRs2aMCAAYqJafjpLDY2VgMGDND69evDPBngPPoB2KMfgD36YRbyMAt5NM2RI0csx61btw7qfkJ1YZyG1J+x/veA5sG7egJwtZ07d+q888475W06deqkb7/9NkwTAeagH4A9+gHYox9mIQ+zkEfT3Hvvvbr11lt15MgRHTlyRB07dmz0fYT6wjj1TZw4UcOHD1dCQoISEhLUqVOnkHydaMcGJQBXO+OMM1RWVnbK25SVlSk+Pj5MEwHmoB+APfoB2KMfZiEPs5BH03Tu3FmdO3cO+vPDcWGc+i655BJdcsklIblvnMCfeANwtV69eumVV17Rvn37Gvz43r17tWLFCvXq1Su8gwEGoB+APfoB2KMfZiEPs5CHs8JxYRw4gw1KAK42fvx47dmzR/379/d7n5d33nlH/fv31969ezVhwgSHJgScQz8Ae/QDsEc/zEIeZiEP54TzwjgIP/7EG4Cr/exnP9OUKVM0d+5c9e/fX16vVx06dNDOnTtVUVEhn8+n++67T8OGDXN6VCDs6Adgj34A9uiHWcjDLOThnHBeGAfhxysoAbjeE088ob/85S8aMmSIzjjjDO3YsUNnnnmmrrnmGr3++uuaPXu20yMCjqEfgD36AdijH2YhD7OQR/iF+8I4CD9eQQlXqqmpcXoEY0XrY3P99dfr+uuvd3oMI0Trz0AgovWxoR8nROvPQCCi9bGhHydE689AIKL1saEfZiEPs5DHCaF+jnTiwjjNJVrXj2DwCkq4gtfrtRxXVlY6NIn5KioqLMdcPS7y0Y/A0Y/oQz8CRz+iD/0IHP0AgBPCvX64+cI4rB+BY4MSrlC/xOXl5Q5NYr76jw1PgJGPfgSOfkQf+hE4+hF96Efg6AcAnBDO9cPtF8Zh/QgcG5RwhXPOOcdy/PXXXzs0ifl27NhhOW7Xrp1Dk4RGTEyM4uLi9K9//avuODY29rT/xcVF7jta0I/A0Q/6QT/s0Q/6QT/s0Y/o64eTyMMs5OEvnOuH2y+ME+nrR3OK3MYgonTt2tVy7JbfljihqKjIcpyRkeHQJKGRlZUlj8ejVq1aWY6jGf0IHP2IPvQjcPQj+tCPwNEPhBN5mIU8/IVr/diyZYvrL4wT6etHc2KDEq5Qv8S7du3SgQMHlJiY6NBEZjpw4IDKysos5yLtCXDdunWnPI5G9CMw9CM60Y/A0I/oRD8CQz8QbuRhFvLwF471w+fzacKECX7n3XBhnFrRsH40J/7EG66Qlpbm91uq+r+JgP9jEhMToy5dujg0DcKFfgSGfkQn+hEY+hGd6Edg6AcAWIVj/aisrFT37t0VE3Ni2+qKK65wxYVxarF+NA4blHAFr9erlJQUyzk3vaw7XFavXm05TklJ8bvCWrQoLS3V4sWL9dxzz+nQoUNOjxNS9CMw9OME+kE/6qMfJ9AP+lEf/TghmvrhBuRhlmjKIxzrh9fr1YIFC7Rx40ZddtllSkxM1IsvvtisXyPUWD8ahw1KuMaQIUMsxwUFBQ5NYq76j0n9xywS/fa3v1WXLl20d+/eunPr1q3ThRdeqJycHN19993q1auXvv/+ewenDD36cXr04zj6QT8aQj+Oox/0oyH047ho7IcpyMMs5HFcuNaPnj176t1339V7772nDh06hORrhEo0rh9NwQYlXGPEiBGW402bNunLL790aBrzlJSUaPPmzZZz9R+zSPTKK68oNTVVbdu2rTuXm5urmpoazZw5U2PHjlVxcbHmzZvn3JBhQD9OjX7Qj5PRDyv6QT9ORj+s6Ed098MU5GEW8jgunOtHTEyMunfvHpL7DpVoXT+agg1KuEbfvn3Vvn17y7n58+c7NI15FixYYDlOSkpSVlaWQ9OET2lpqWWx+uabb/SPf/xD48eP14MPPqgFCxZowIABWrFihYNThh79ODX6cRz9OIF+nEA/jqMfJ9CPE+jHcdHaD1OQh1nI4zjWj1OL1vWjKdighGvExcXp5ptvtpybP3++PvvsM4cmMsenn37qtxjcdNNNiouLc2ii8Dlw4IDatGlTd7x+/Xp5PB7dcMMNded69eqlr776yoHpwod+2KMfbeqO6ccJ9OM4+tGm7ph+nEA/jqMfbeqOo7UfpiAPs5DHcawf9qJ5/WgKNijhKlOnTrW8qWx1dbXGjx+vmpoaB6dyVk1NjcaPH6/q6uq6c16vV1OnTnVwqvBJTk7W9u3b645Xr14tr9erzMzMunNHjx71u8pcJKIf/ugH/ahFP/zRD/pRi374ox/0wyTkYRbyOIH1w1+0rx9NwQYlXCU9PV2//vWvLefWr1+vnJycqHwSrKmpUU5OjjZs2GA5n5ubq/T0dIemCq9LL71UK1eu1GuvvaY1a9boz3/+s/r3729ZKLdt26aOHTs6OGV40A8r+kE/TkY/rOgH/TgZ/bCiH/TDNORhFvI4gfXDivWjadighOvcf//9SklJsZxbvHhx1D0J1j75LV682HI+NTVV999/v0NThd8DDzygqqoq/exnP9OQIUN09OhRPfDAA3Ufr6io0Pr16y2/0Yxk9OM4+nEc/bCiH8fRj+PohxX9OI5+HEc/zEIeZiEPK9aP41g/mo4NSrhOq1attGjRIr/3b1i8eLH69eunTz/91KHJwufTTz9Vv379/J784uLilJeXp4SEBIcmC79evXrp/fff15QpUzRlyhT93//9n6644oq6j//zn/9U//79ddtttzk4ZfjQD/pxMvphRT/ox8nohxX9oB8nox9mIQ+zkIcV6wfrR3PhHTrhSoMHD9by5cs1cuRIVVVV1Z3fsGGDevbsqYkTJ2rChAkR9zLqkpISLViwQPPnz7e8p4V0/Mlv+fLlGjx4sEPTOefiiy/WxRdf3ODHfvKTn0T8FfTqox/042T0w4p+0I+T0Q8r+kE/TkY/zEIeZiEPK9YP1o/mwCso4Vo33nijli9f7vebmurqas2bN09du3ZV7969NWvWLG3atEkHDhxwaNLgHThwQJs2bdKsWbPUu3dvde3aVfPmzbN98rvxxhsdmhSmoR8n0A/URz9OoB+oj36cQD8AIHCsHyewfgSHV1DC1W688Ua9/vrrGjNmjEpLS/0+vnnzZm3evFnTpk2TJCUlJSkjI0OdO3dWq1at5PV6FRNjxj59TU2NKioqVF5erh07dqioqEhlZWWn/bzU1FTl5eVF/W9mduzYobVr1+qbb75RRUWF38c9Ho+mT5/uwGTOoR/0oxb98Ec/6Ect+uGPftCPWvTDLORhFvLwx/rB+tEUbFDC9QYPHqwtW7Zo1qxZmj17doOLQ62ysrKAnlTcwOv1Kjc3V/fff3/Uv6fFfffdp//5n/+x/ObK5/PJ4/FY/ne0/QNBoh/0g36cCv2gH/TDHv2gH/TDLORhFvKwx/rB+hEsM7amgSZKSEjQzJkzVVhYqJycHLVv397pkUImKSlJOTk5Kiws1MyZM6P+yW/x4sV64okn1L9/f7344ovy+Xy666679Kc//Uk5OTmKi4vTiBEj9Pbbbzs9qmPoR/SiH6dHP6IX/Tg9+hG96IdZyMMs5HF6rB8IBq+gRERJT0/X008/rfnz52v9+vUqKCjQqlWrVFpaKp/P5/R4QfF4PEpNTdWQIUM0YsQIZWVl+b2vRzRbtGiRUlNT9cYbb9T9OUBqaqpuueUW3XLLLRo5cqQGDRqkESNGODyp8+hH9KEfgaMf0Yd+BI5+RB/6YRbyMAt5BI71A43Bo4iIFBcXp6uvvlpXX321JKmiokLbtm1TUVGRioqKtGfPHh09elRHjx51eFKr+Ph4xcfHq127dsrIyFBGRoa6dOkir9fr9GjG2rp1q+644w7Le5WcfOW4vn376rrrrtPjjz+u4cOHOzGicehH9KAfjUc/ogf9aDz6ET3oh1nIwyzk0XisHwgEG5SICl6vV926dVO3bt2cHgUh0KZNm7r/fcYZZ2jPnj2Wj59//vlas2ZNmKdyD/oR2ehH09CPyEY/moZ+RDb6YRbyMAt5NA3rBxrCe1ACcLVOnTppx44ddcfp6en64IMPLLf57LPPdMYZZ4R7NMBx9AOwRz8Ae/TDLORhFvIAQoMNSgCudsUVV+j999+vO/7Zz36mf/7znxozZoxef/11TZs2TW+88YaysrIcnBJwBv0A7NEPwB79MAt5mIU8gNDgT7wBuNodd9yhb775Rtu3b1dKSoruu+8+vfbaa1q8eLGWLFkin8+n1NRUzZkzx+lRgbCjH4A9+gHYox9mIQ+zkAcQGmxQAnC1fv36qV+/fnXHZ555pt5//3298sor+vLLL5WSkqIbbriBP7FAVKIfgD36AdijH2YhD7OQBxAabFACiDgtWrTQiBEjnB4DMBL9AOzRD8Ae/TALeZiFPICm4z0oAbiSz+fTunXr9OKLL6q4uLju/D/+8Q8NGTJEZ599ttq0aaOf/exn+vzzzx2cFAg/+gHYox+APfphFvIwC3kAocUrKAG4zuHDhzV48OC6N6f2eDyaM2eOBgwYoP79+6u8vLzutq+++qr+7//+T//85z/VuXNnp0YGwoZ+APboB2CPfpiFPMxCHkDo8QpKAK7zxBNP6L333tMll1yiKVOm6OKLL9aDDz6o6dOnKykpSatXr9aBAwe0fft2/ed//qf27Nmj2bNnOz02EBb0A7BHPwB79MMs5GEW8gBCj1dQAnCdl156SV26dNH777+vFi1a6NixY7rgggv02muvaeXKlRowYICk429YPW/ePL377rtatWqVw1MD4UE/AHv0A7BHP8xCHmYhDyD0eAUlANcpKSnRNddcoxYtWkg6/qbUP/3pTyVJV1xxhd/tr7jiCu3YsSOsMwJOoR+APfoB2KMfZiEPs5AHEHpsUAJwnfLycrVv395y7pxzzpEktW3b1u/2Z599tioqKsIyG+A0+gHYox+APfphFvIwC3kAoccGJQBX8ng8pzwGohn9AOzRD8Ae/TALeZiFPIDQYoMSAAAAAAAAgGO4SA4AV/rDH/6g999/v+64uLhYknTttdf63bb2Y0C0oB+APfoB2KMfZiEPs5AHEFpsUAJwpeLi4gYX/r/97W8N3p4/wUA0oR+APfoB2KMfZiEPs5AHEFpsUAJwnW3btjk9AmAs+gHYox+APfphFvIwC3kAoccGJQDXSUlJcXoEwFj0A7BHPwB79MMs5GEW8gBCj4vkAAAAAAAAAHAMG5QAAAAAAAAAHMOfeCMqVFRU6Msvv1RRUZGKi4v13Xff6ejRo6qoqHB6NAuv16v4+Hidc8456tq1qzIyMpSWliav1+v0aIhg9AOwRz8Ae/QDABAM1g80hA1KRKSqqiq98847Kigo0KpVq7R9+3b5fD6nxwqKx+NRSkqKhgwZohEjRqhv376Ki6O6CB79AOzRD8Ae/QAABIP1A4HgT7wRUYqLizV27Fh17NhRAwcOVF5enkpLS1375CdJPp9PpaWlysvL08CBA9WxY0eNHTtWJSUlTo8Gl6EfgD36AdijHwCAYLB+oDHYoEREKC8v14wZM9SjRw8tXLhQu3fvdnqkkNm9e7cWLlyoHj16aMaMGSovL3d6JBiOfgD26Adgj34AAILB+oFg8DpUuN6bb76p0aNHa/v27ae9bXJysrp27arzzjtPrVq1UsuWLRUTY8Y+fU1NjSorK1VeXq6vv/5axcXF2rVrl+3tKyoq9PDDD+u5555TXl6eBg8eHMZp4Rb0g37AHv2gH7BHP+gHAASD9YP1I1hsUMLVVqxYoZEjR6qqqqrBj/fu3VsjRozQwIEDlZGRocTExDBP2DQHDhxQUVGR1qxZo+XLl2vz5s1+tyktLdV1112n5cuX68Ybb3RgSpiKftAP2KMf9AP26Af9AIBgsH6wfjSFGVvTQBDsnvxiY2M1efJklZSUaOPGjcrNzVXv3r1d9+QnSYmJierdu7dyc3O1adMmFRcXa/LkyYqNjbXcrqqqSiNHjtSKFSscmhSmoR8n0A/URz9OoB+oj36cQD8AIHCsHyewfgSHDUq40ptvvtngk19WVpY++ugjzZ07V2lpaQ5NFzrp6emaO3euPvroI2VlZVk+Vvsk+Oabbzo0HUxBP+gH7NEP+gF79IN+AEAwomH92LBhg9851o/mxQYlXKe8vFyjR4/2e/LLzs7W2rVrdeGFFzo0WfhceOGFWrt2rbKzsy3nq6qqNGbMGB05csShyeA0+kE/YI9+0A/Yox/0AwCCEenrR01Nje6//34tW7bM9jasH82DDUq4zqxZs/zecDc7O1sLFy405g11wyEmJkYLFy70exIsLS3VrFmzHJoKTqMfx9EPNIR+HEc/0BD6cRz9AIDGieT148iRI7rllls0e/ZsXXbZZae8LetH07n7pwVRp7i4WI899pjlXFZWVkQ8+QWj9knwqquuspyfPXu2SkpKHJoKTqEfVvQDJ6MfVvQDJ6MfVvQDAAITyevHrl271K9fP7344ouSdNoNSon1o6nc/RODqPPEE0+ooqKi7jg2NlZPPfWU65/8miImJkZPPfWU5Y15Kyoq9Pjjjzs4FZxAP/zRD9SiH/7oB2rRD3/0AwBOL1LXj8LCQmVmZurDDz+UJLVp00bnn39+QJ/L+hE8d//UIKpUVVXppZdespybOHGi69/TojlcdNFFmjhxouXcyy+/7Pc+IIhc9MMe/QD9sEc/QD/s0Q8AsBep68eaNWt0+eWXW/5sPTMzs1GbrqwfwWGDEq7xzjvvaPfu3ZZz9UsfzSZMmGA5Lisr0/r16x2aBuFGP06NfkQ3+nFq9CO60Y9Tox8A0LBIXD+WLFmia665RgcOHLCcD+TPu+tj/Wg8NijhGgUFBZbj3r17Ky0tzaFpzJOenq5evXpZztV/zBC56Mep0Y/oRj9OjX5EN/pxavQDABoWSetHTU2NcnNzlZ2d3eCrHIPZoGT9aDw2KOEaq1atshyPGDHCoUnMVf8xqf+YIXLRj9OjH9GLfpwe/Yhe9OP06AcA+IuU9ePIkSMaOXKk38V+ank8HmVmZgZ136wfjcMGJVyhoqLC8h4QkjRw4ECHpjHXoEGDLMfbt2+3vGkxIhP9CAz9iE70IzD0IzrRj8DQDwCwipT1o/ZK3fXfS/Nk3bt3V+vWrYO6f9aPxmGDEq7w5ZdfyufzWc798Ic/dGgac2VkZFiOa2pqtG3bNoemQbjQj8DQj+hEPwJDP6IT/QgM/QAAq0hYP+pfqdtOMH/eXYv1o3HYoIQrFBUVWY6Tk5N11llnOTSNuRITE5WUlGQ5V/+xQ+ShH4GhH9GJfgSGfkQn+hEY+gEAVm5fP1avXu13pW5Jat26td/VupuyQcn60ThsUMIViouLLcddu3Z1aBLz1f8tDU+AkY9+BI5+RB/6ETj6EX3oR+DoBwCc4Ob1Y/HixQ1eqTs1NVXLly9XTU2N5XxTNigl1o/GYIMSrvDdd99Zjs877zyHJjFf586dLcd79uxxaBKEC/0IHP2IPvQjcPQj+tCPwNEPADjBjetH7ZW6R48ererqasvHMjMz9f7772v//v2W823atNH555/fpK/L+hG4OKcHAAJx9OhRy3GrVq0cmsR89R+b+o8dIg/9CBz9iD70I3D0I/rQj8DRDwA4wW3rR3l5ue68884GL4YzYsQILV26VAkJCXrvvfcsH8vMzPT7k+/GYv0IHBuUcIX6V7pq2bKlQ5OYz+v1Wo55Aox89CNw9CP60I/A0Y/oQz8CRz8A4AQ3rR/ff/+9rrnmGtuL4XTo0EEJCQmS5LdB2dQ/75ZYPxqDP/GGKzX1txiRjMcG/AzY47EBPwP2eGzAz4A9HhsAsGfyc2SbNm2Uk5Oj9u3bN/jx+fPny+PxaPXq1dq8ebPlY82xQWnyY2MaHikAAAAAAABEnJiYGI0aNUpffPGFxo8fb7thOHjwYFVWVtYdezweZWZmhmtMiA1KAAAAAAAARLC2bdtqwYIF2rhxY0C37969u1q3bh3iqXAyNigBAAAAAAAQ8epfqdtOc/x5NxqHDUoAAAAAAABENJ/Pp/79+wd0WzYow4+reAMAAAAAAMB1jh49qi1btqiwsFAHDx6su0p2fHy8zjrrLPXo0UM9evSQ1+vVpEmTGryPzZs3a/z48ZareLNBGX5sUAIAAAAAAMB4Bw8eVEFBgd566y19/PHH2rp1q6qrq0/5ObGxscrIyNDWrVv9PlZcXKz09HS9++67Wrp0qXJzc3Xs2DGdf/75ofoWYIMNSgAAAAAAABjJ5/Np/fr1evbZZ1VQUKDy8vJGfX51dXWDm5MXX3yx0tLSJJ242vewYcP0yiuv2F7tG6HDIw4AAAAAAADjrFmzRj169FC/fv20dOnSRm9OnsrHH3+sHj16aM2aNXXn2rZtq1GjRjXb10DgInqD8ujRo/rVr36lrKwsdezYUfHx8erQoYOuuOIKPfvsszp27JjTI8JwpaWl8ng8lv9atGihTp06aeTIkdq4caMkad68efJ4PKd8Ilu3bp1iYmJ06aWXqqqqKlzfAhAy9AOwRz8Ae/QDAHA633zzjW699VYNGjRIn3/++Wlvn56erqysLA0aNEiDBg1SVlaW0tPTT/t5n3/+uQYNGqSf//zn+vbbb5tjdAQpov/E+9ChQ3r66afVp08fXXfddWrfvr327t2rN954Q7/4xS/0wgsv6I033uCluzit9PR03X777ZKkw4cPa9OmTSooKNArr7yiNWvWaNKkSVq5cqXy8/N100036YYbbrB8/qFDhzRq1Ch5vV4tW7ZMcXERXT1EGfoB2KMfgD36AQCoz+fzadGiRbrvvvt08ODBBm8TGxur6667Tj/96U918cUX66KLLtJZZ53V4G1zcnKUl5d32q/7wgsv6PXXX9ecOXM0ZsyYJn0PCE5Er+Jnn3229u/fr5YtW1rOV1VVadCgQXrzzTf1xhtv6LrrrnNoQrhF165d9dBDD1nOzZo1S9OmTdP06dP1zjvvKD8/Xz/60Y+UnZ2twsJCtWvXru629957r0pLSzV37lxdcMEFYZ4eCC36AdijH4A9+gEAOFlVVZUmT56sp556qsGPd+vWTffcc49uv/12dejQ4bT3t2/fvgY3J3Nzc7Vy5Uq/96U8ePCgcnJy9Nlnn2nu3Ln84ivMIvqlgzExMX6bk5IUFxenG2+8UdLxKzYBwbjnnnskSZs2bZIkpaSkaN68edq1a5fGjh1bd7tVq1Zp0aJF6t+/vyZNmuTIrEC40Q/AHv0A7NEPAIhOlZWVGjlyZIObk23btlVeXp4KCws1derUgDYnJencc8/1O5eVlaVZs2apsLBQCxcuVJs2bfxus2DBAo0cOVKVlZWN/j4QvIjeoLRTU1Ojv/3tb5KkCy+80OFp4HYn/1Zl1KhRGjp0qAoKCvSnP/1J+/bt0y9/+UslJibq2WeflcfjcXBSIPzoB2CPfgD26AcARI9jx45p+PDhWrFihd/H7rrrLm3dulWjR49u1NvzrVu3TkePHvU7X3tBnJiYGI0ZM0ZffPGF7rrrLr/brVixQsOHD+faJWEUFa9Xrays1G9/+1v5fD7t2bNHb731lrZu3apRo0ZpwIABTo8Hl1qyZIkk6corr7ScX7Rokf7v//5P48ePV1ZWlnbs2KFnnnlGKSkpTowJOIJ+APboB2CPfgBA9Jk+fbpeffVVy7mWLVsqPz9fP//5zxt9fz6fT/379/c7v2TJErVo0cJyLikpSfn5+Ro8eLBGjRpledXkq6++qhkzZujRRx9t9AxovKjZoJw5c2bdscfj0dSpU135Q+bz+VReXu70GGHn9G8tiouL694jqfZN3NeuXavk5GTNmTPHctvk5GTl5eXp5ptv1sqVKzV06NBTXn0y1I4dO6bDhw83+30iOJGYB/3wv08EJxLzoB/+94ngRGIe9MP/PhEc8jALeZjF9DxWr16t2bNnW86dccYZeu2119SvX7+g7tPurT9q30akIbfddps6duyo66+/3vJ4zZ49WwMGDNDAgQODmgWN4Isi1dXVvq+//tr3u9/9ztemTRvfFVdc4du/f7/TYzXKoUOHfJKi/r9x48aF5fHetm2b7QwdOnTwFRUV2X5unz59fJJ8W7ZsCcustcaNGxexebhRJOdBP8zKw40iOQ/6YVYebhTJedAPs/JwI/IwC3mYxU157Nq1y5ecnGy5rxYtWvg2bNgQ9Pe/d+/eBmcsLi4O6PM3bNjga9Gihd/atGvXrqDmqZ8HP7v2ouo9KGNiYtS5c2eNHTtWixYt0t///nf993//t9NjwQWGDBkin88nn8+nsrIyzZkzR2VlZRo6dKgOHTrU4OckJCRY/i8QqegHYI9+APboBwBEt5ycHO3atcty7rHHHvN7m4/GsLswTnp6ekCff+WVV/q9onPnzp2WC7UhNKLiT7wbMnjwYEnH3zjVTVq1amX7D7ZINmXKFC1evNjpMSRJ7du319SpU7V//3498sgjevDBBzVv3jynx7KVnZ2tuXPnNut9mpSH20R6HvTDrDzcJtLzoB9m5eE2kZ4H/TArD7chD7OQh1lMzaOwsNDvojjXXnut7Z9nB+J0F8YJ1KRJk7R69Wq98cYbdedefvllbdmyRd27dw96Ppxa1G5QfvPNN5Lk9wappvN4PDrjjDOcHiPsTMzpgQce0DPPPKPf/e53mjx5slJTU50eqUEtWrRo9p8ZE/Nwi2jJg34gGNGSB/1AMKIlD/qBYJCHWcjDLKbm8eSTT1qOzz77bOXn58vj8QR1f75GXBjndGJiYpSfn69u3bpp7969lplrL+aG5hfRf+K9ZcuWBi8oU15erl/96leSju/QA8FISEhQbm6ujh07pocfftjpcQCj0A/AHv0A7NEPAIh8O3fu1B/+8AfLuXHjxql9+/ZB32cwF8Y5laSkJI0bN85y7rnnntPOnTuDuj+cXkRvUC5fvlwdOnTQtddeq3Hjxun+++/XHXfcoR/84Af629/+pquuukpTpkxxeky42OjRo9WxY0ctW7ZMJSUlTo8DGIV+APboB2CPfgBAZMvLy1NlZWXdccuWLTVhwoSg72/fvn2aP3++3/ni4uKg71OSJkyYoJYtW9YdV1ZWKi8vr0n3CXsRvUF5/fXX69Zbb9VXX32lP/3pT3riiSf0xhtv6Ec/+pHy8vL09ttv8wbbaJL4+HhNmzZNVVVVmjlzptPjAEahH4A9+gHYox8AENnefvtty/Edd9yh5OTkoO+vqRfGsdOhQwfdfvvtlnNr165t0n3CXkS/B+WPf/xj/fjHP3Z6DLhYamqqfD7fKW8zYcKEBn/b47YLMAGNRT8Ae/QDsEc/ACB6VVVVaePGjZZzw4YNC/r+muvCOHaGDRumZ555pu5448aNqqqqUlxcRG+nOSKiX0EJAAAAAAAAM3z22Wd+1wrJzMwM6r6a88I4durPdvjwYRUWFjbLfcOKDUoAAAAAAACE3Pvvv285Tk9PD/riOM19YZyGJCUlKS0tzXKu/veA5sEGJQAAAAAAAELuq6++shz37NkzqPsJ1YVxGlJ/xvrfA5oHG5QAAAAAAAAIuSNHjliOW7duHdT9hOrCOA2pP2P97wHNg3f1BAAAAAAAQMjde++9uvXWW3XkyBEdOXJEHTt2bPR9hPrCOPVNnDhRw4cPV0JCghISEtSpU6eQfJ1oxwYlAAAAAAAAQq5z587q3Llz0J8fjgvj1HfJJZfokksuCcl94wT+xBsAAAAAAADGC8eFceAMNigBAAAAAABgtHBeGAfhxwYlAAAAAAAAjBbOC+Mg/NigBAAAAAAAgLHCfWEchB8blHClmpoap0cwlhOPDXnYIw+zkIdZyMMs5GEW8jALeZiFPMxCHmaJxDycuDBOc+FnNXBsUMIVvF6v5biystKhScxXUVFhOY6Pj2/2r0EegSMPs5CHWcjDLORhFvIwC3mYhTzMQh5micQ83HxhnHDkESnYoIQr1C9xeXm5Q5OYr/5jE4onQPIIHHmYhTzMQh5mIQ+zkIdZyMMs5GEW8jBLpOXh9gvjhCOPSMEGJVzhnHPOsRx//fXXDk1ivh07dliO27Vr1+xfgzwCRx5mIQ+zkIdZyMMs5GEW8jALeZiFPMwSaXm4/cI44cgjUrBBCVfo2rWr5dgtvy1xQlFRkeU4IyOj2b8GeQSOPMxCHmYhD7OQh1nIwyzkYRbyMAt5mCWS8tiyZYvrL4wTjjwiBRuUcIX6Jd61a5cOHDjg0DTmOnDggMrKyiznQvEESB6BIQ+zkIdZyMMs5GEW8jALeZiFPMxCHmaJpDx8Pp8mTJjgd94NF8apFa48IgUblHCFtLQ0eTwey7n6v4mA/2MSExOjLl26NPvXIY/AkIdZyMMs5GEW8jALeZiFPMxCHmYhD7NEUh6VlZXq3r27YmJObFtdccUVrrgwTq1w5REp2KCEK3i9XqWkpFjOuell3eGyevVqy3FKSorfFdaaA3kEhjzMQh5mIQ+zkIdZyMMs5GEW8jALeZglkvLwer1asGCBNm7cqMsuu0yJiYl68cUXm/VrhFq48ogUbFDCNYYMGWI5LigocGgSc9V/TOo/Zs2JPE6PPMxCHmYhD7OQh1nIwyzkYRbyMAt5mCUS8+jZs6feffddvffee+rQoUNIvkaohDOPSMAGJVxjxIgRluNNmzbpyy+/dGga85SUlGjz5s2Wc/Ufs+ZEHqdGHmYhD7OQh1nIwyzkYRbyMAt5mIU8zBLJecTExKh79+4hue9QCXcekYANSrhG37591b59e8u5+fPnOzSNeRYsWGA5TkpKUlZWVsi+HnmcGnmYhTzMQh5mIQ+zkIdZyMMs5GEW8jALeZgl3HlEBB/gIjk5OT5Jdf/Fxsb6Pv30U6fHctwnn3zii42NtTw2OTk5If+65NEw8jALeZiFPMxCHmYhD7OQh1nIwyzkYRbyMItTebgdG5RwleLiYp/X67UUPSsry1ddXe30aI6prq72XXXVVZbHxOv1+oqLi0P+tcnDH3mYhTzMQh5mIQ+zkIdZyMMs5GEW8jALeZjFyTzcjg1KuM706dMtZZfky87Ojsonwerqal92drbf4zFjxoywzUAeJ5CHWcjDLORhFvIwC3mYhTzMQh5mIQ+zkIdZTMjDzdighOscPnzYl5KSEvVPgnZPfqmpqb7y8vKwzUEex5GHWcjDLORhFvIwC3mYhTzMQh5mIQ+zkIdZTMnDzdighCutWrXKFxcX51f+q666yvfJJ584PV7IffLJJ34vG5fki4uL861atSrs85AHeZiEPMxCHmYhD7OQh1nIwyzkYRbyMAt5mMW0PNyKDUq41ssvv9zgk2BsbKxv8uTJEfkeD8XFxb7Jkyf7veFu7ZPfyy+/7Nhs5EEeTiMPs5CHWcjDLORhFvIwC3mYhTzMQh5mMTkPN2KDEq5m9yRY+1+vXr18jz76qG/jxo2+/fv3Oz1uo+3fv9+3ceNG36OPPurr1auX7fdpypMfeZBHOJGHWcjDLORhFvIwC3mYhTzMQh5mIQ+zuC0Pt/H4fD6fABd78803NWbMGJWWlp72tklJScrIyFDnzp3VqlUreb1excTEhH7IANTU1KiiokLl5eXasWOHioqKVFZWdtrPS01NVV5engYPHhyGKU+PPMgjFMiDPEKBPMgjFMiDPEKBPMgjFMiDPEKBPMzKw1Wc3iEFmkN5eblvxowZPq/Xa/tbjEj7z+v1+mbMmGHkG+6Sh1nIwyzkYRbyMAt5mIU8zEIeZiEPs5CHWcgDwWCDEhGluLjYl5OT42vfvr3jT1Ch+i8pKcmXk5PjivfwIA+zkIdZyMMs5GEW8jALeZiFPMxCHmYhD7OQBxqDP/FGRKqqqtL69etVUFCgVatWqbS0VG79Ufd4PEpNTdWQIUM0YsQIZWVlKS4uzumxGoU8zEIeZiEPs5CHWcjDLORhFvIwC3mYhTzMQh4IBBuUiAoVFRXatm2bioqKVFRUpD179ujo0aM6evSo06NZxMfHKz4+Xu3atVNGRoYyMjLUpUsXeb1ep0drVuRhFvIwC3mYhTzMQh5mIQ+zkIdZyMMs5GEW8kBD2KAEAAAAAAAA4BgzLo8EAAAAAAAAICqxQQkAAAAAAADAMWxQAgAAAAAAAHAMG5QAAAAAAAAAHMMGJQAAAAAAAADHsEEJAAAAAAAAwDFsUAIAAAAAAABwDBuUAAAAAAAAABzDBiUAAAAAAAAAx7BBCQAAAAAAAMAxbFACAAAAAAAAcAwblAAAAAAAAAAcwwYlAAAAAAAAAMewQQkAAAAAAADAMWxQAgAAAAAAAHAMG5QAAAAAAAAAHMMGJQAAAAAAAADHsEEJAAAAAAAAwDFsUAIAAAAAAABwDBuUAAAAAAAAABzDBiUAAAAAAAAAx7BBCQAAAAAAAMAxbFACAAAAAAAAcAwblAAAAAAAAAAcwwYlAAAAAAAAAMewQQkAAAAAAADAMWxQAgAAAAAAAHAMG5QAAAAAAAAAHMMGJQAAAAAAAADHsEEJAAAAAAAAwDFsUAIAAAAAAABwDBuUAAAAAAAAABzDBiUAAAAAAAAAx7BBCQAAAAAAAMAxbFACAAAAAAAAcAwblAAAAAAAAAAcwwYlAAAAAAAAAMewQQkAAAAAAADAMWxQAgAAAAAAAHAMG5QAAAAAAAAAHMMGJQAAAAAAAADHsEEJAAAAAAAAwDFsUAIAAAAAAABwDBuUAAAAAAAAABzDBiUAAAAAAAAAx7BBCQAAAAAAAMAxbFACAAAAAAAAcAwblAAAAAAAAAAcwwYlAAAAAAAAAMewQQkAAAAAAADAMWxQAgAAAAAAAHAMG5QAAAAAAAAAHMMGJQAAAAAAAADHsEEJAAAAAAAAwDFsUAIAAAAAAABwDBuUAAAAAAAAABzDBiUAAAAAAAAAx7BBCQAAAAAAAMAxbFACAAAAAAAAcAwblAAAAAAAAAAcwwYlAAAAAAAAAMewQQkAAAAAAADAMWxQAgAAAAAAAHAMG5QAAAAAAAAAHMMGJQAAAAAAAADHsEEJAAAAAAAAwDFsUAIAAAAAAABwDBuUAAAAAAAAABzDBiUAAAAAAAAAx/x/q954a1INX58AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1300x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_qubits = 4\n",
    "n_layers = 3\n",
    "\n",
    "full_dev = qml.device('default.qubit', wires = 4)\n",
    "inp_params = torch.rand(2*4*n_qubits)\n",
    "weight_params = torch.rand(3*n_qubits)\n",
    "nn_params = torch.rand(n_layers*n_qubits)\n",
    "\n",
    "full_qnode = qml.QNode(full_circuit, full_dev)\n",
    "\n",
    "qml.draw_mpl(full_circuit, range(4))(full_qnode, torch.tensor([0, 1, 2, 3]), inp_params, weight_params, nn_params, n_layers)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate circuits expectation values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will use a circuit for each variable this will help us to evaluate everything\n",
    "def get_predicts(qnodes, X_batch, FM_input_weights, \n",
    "                FM_train_weights, NN_train_weights, TL_train_weights,\n",
    "                n_layers, FMrot_seq, NNrot_seq, calc_features):\n",
    "    \"\"\" \n",
    "    **Args**:\n",
    "        * qnode: *Single*[pennylane.workflow.qnode.QNode] Qnode where circuit runs\n",
    "        * X_batch: *Sequence*[float] features from the training dataset, should be in the form (t, x, y)\n",
    "        * FM_input_weights: *Sequence*[float] trainable parameters to linear transform the inputs (array should be twice as long as the inputs array)\n",
    "        * FM_train_weights: *Sequence*[float] trainable parameters of basic entangling layers on the feature map\n",
    "        * FM_train_weights: *Sequence*[float] trainable parameters of basic entangling layers on the neural network\n",
    "        * TL_train_weights: *Sequence*[float] trainable parameters for linear transformation on output\n",
    "        * n_layers: *Sequence*[int] number of layers of the neural network\n",
    "        * FMrot_seq: *Single or Sequence*[pennylane.ops.Operation]  Rotations used in BasicEntangledLayers on the featuremap if Single, use the same roation operation for each basic entangle embedding, if Sequence should be # of inputs -1 long.\n",
    "        * NNrot_seq: *Single or Sequence*[pennylane.ops.Operation]  Rotations used in BasicEntangledLayers on the Neural network if Single, use the same roation operation for each basic entangle embedding, if Sequence should be # of inputs -1 long.\n",
    "        * calc_features: *Array*[int] which variables to calculate 0 = wind_speed_u, 1 = wind_speed_v, 2 = wind_speed_w, 3 = wind_speed_geopotential\n",
    "    **Returns**:\n",
    "        * Expected value: Z(0)@Z(1)@...@Z(n_qubits)\n",
    "    \"\"\"\n",
    "\n",
    "    # Lets run each circuit with their own parameters\n",
    "    predictions = []\n",
    "    for ii in calc_features:\n",
    "        res = TL_train_weights[ii][0]*qnodes[ii](qnodes[ii], X_batch, FM_input_weights[ii], \n",
    "                                              FM_train_weights[ii], NN_train_weights[ii], \n",
    "                                              n_layers, FMrot_seq, NNrot_seq) + TL_train_weights[ii][1]\n",
    "        predictions.append(res)\n",
    "    return predictions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# PDE residual"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate f0 and beta from the differential equations\n",
    "\n",
    "phi0 = 11 # hurricanes initial latitud\n",
    "phi0_rad = np.pi*phi0/180 # transform to radians\n",
    "\n",
    "Omega = 7.2921e-5\n",
    "f0 = torch.tensor(2*Omega*np.sin(phi0_rad), dtype = torch.float32).to(torch_dev)\n",
    "beta = torch.tensor(2*Omega*np.cos(phi0_rad)/radius_earth, dtype = torch.float32).to(torch_dev)\n",
    "\n",
    "def pde_loss(qnodes, X_batch, FM_input_weights, FM_train_weights, NN_train_weights, TL_train_weights, n_layers, FMrot_seq, NNrot_seq):\n",
    "        \"\"\" \n",
    "        **Args**:\n",
    "                * qnode: *Single*[pennylane.workflow.qnode.QNode] Qnode where circuit runs\n",
    "                * X_batch: *Sequence*[float] features from the training dataset, should be in the form (t, x, y)\n",
    "                * FM_input_weights: *Sequence*[float] trainable parameters to linear transform the inputs (array should be twice as long as the inputs array)\n",
    "                * FM_train_weights: *Sequence*[float] trainable parameters of basic entangling layers on the feature map\n",
    "                * FM_train_weights: *Sequence*[float] trainable parameters of basic entangling layers on the neural network\n",
    "                * TL_train_weights: *Sequence*[float] trainable parameters for linear transformation on output\n",
    "                * n_layers: *Sequence*[int] number of layers of the neural network\n",
    "                * FMrot_seq: *Single or Sequence*[pennylane.ops.Operation]  Rotations used in BasicEntangledLayers on the featuremap if Single, use the same roation operation for each basic entangle embedding, if Sequence should be # of inputs -1 long.\n",
    "                * NNrot_seq: *Single or Sequence*[pennylane.ops.Operation]  Rotations used in BasicEntangledLayers on the Neural network if Single, use the same roation operation for each basic entangle embedding, if Sequence should be # of inputs -1 long.\n",
    "        \n",
    "        **Returns**:\n",
    "                * Navier stokes residual: *Sequence[float]*\n",
    "        \"\"\"\n",
    "        # Num qubitsF\n",
    "        n_qubits = len(qnodes[0].device.wires)\n",
    "        \n",
    "        # U and V components of wind speed val\n",
    "        u, v, w, _ = get_predicts(qnodes, X_batch, FM_input_weights, \n",
    "                FM_train_weights, NN_train_weights, TL_train_weights,\n",
    "                n_layers, FMrot_seq, NNrot_seq)\n",
    "\n",
    "        # partials of t\n",
    "        par_t_u = 0\n",
    "        par_t_v = 0\n",
    "        for ii in range(n_qubits):\n",
    "                # Zeros vec\n",
    "                zeros_vec = torch.zeros(FM_input_weights.shape, dtype = torch.float32, device = torch_dev)\n",
    "                zeros_vec[:, 2*ii+1] = torch.pi/2       #element to modify\n",
    "\n",
    "                # Plus and minus shifts\n",
    "                plus_train_FM = FM_input_weights+zeros_vec\n",
    "                minus_train_FM = FM_input_weights-zeros_vec\n",
    "\n",
    "                # calculate\n",
    "                pos_term_u, pos_term_v, _, _ = get_predicts(qnodes, X_batch, plus_train_FM, \n",
    "                                                        FM_train_weights, NN_train_weights, TL_train_weights,\n",
    "                                                        n_layers, FMrot_seq, NNrot_seq, [0, 1])\n",
    "                neg_term_u, neg_term_v, _, _ = get_predicts(qnodes, X_batch, minus_train_FM, \n",
    "                                                        FM_train_weights, NN_train_weights, TL_train_weights,\n",
    "                                                        n_layers, FMrot_seq, NNrot_seq, [0, 1]) \n",
    "\n",
    "                # Partial accumulations\n",
    "                par_t_u += (FM_input_weights[0][2*ii]/2)*(pos_term_u - neg_term_u)\n",
    "                par_t_v += (FM_input_weights[1][2*ii]/2)*(pos_term_v - neg_term_v)\n",
    "\n",
    "        # partials of x \n",
    "        par_x_u = 0\n",
    "        par_x_v = 0\n",
    "        par_x_G = 0\n",
    "        for ii in range(n_qubits):\n",
    "                # Zeros vec\n",
    "                zeros_vec = torch.zeros(FM_input_weights.shape, dtype = torch.float32, device = torch_dev)\n",
    "                zeros_vec[:, 2*n_qubits + 2*ii+1] = torch.pi/2       #element to modify\n",
    "\n",
    "                # Plus and minus shifts\n",
    "                plus_train_FM = FM_input_weights+zeros_vec\n",
    "                minus_train_FM = FM_input_weights-zeros_vec\n",
    "\n",
    "                # calculate\n",
    "                pos_term_u, pos_term_v, _, pos_term_G = get_predicts(qnodes, X_batch, plus_train_FM, \n",
    "                                                        FM_train_weights, NN_train_weights, TL_train_weights,\n",
    "                                                        n_layers, FMrot_seq, NNrot_seq, [0, 1, 3])\n",
    "                neg_term_u, neg_term_v, _, neg_term_G = get_predicts(qnodes, X_batch, minus_train_FM, \n",
    "                                                        FM_train_weights, NN_train_weights, TL_train_weights,\n",
    "                                                        n_layers, FMrot_seq, NNrot_seq, [0, 1, 3]) \n",
    "\n",
    "                par_x_u += (FM_input_weights[0][2*n_qubits + 2*ii]/2)*(pos_term_u - neg_term_u)\n",
    "                par_x_v += (FM_input_weights[1][2*n_qubits + 2*ii]/2)*(pos_term_v - neg_term_v)\n",
    "                par_x_G += (FM_input_weights[3][2*n_qubits + 2*ii]/2)*(pos_term_G - neg_term_G)\n",
    "\n",
    "        # partials of y \n",
    "        par_y_u = 0\n",
    "        par_y_v = 0\n",
    "        par_y_G = 0\n",
    "        for ii in range(n_qubits):\n",
    "                # Zeros vec\n",
    "                zeros_vec = torch.zeros(FM_input_weights.shape, dtype = torch.float32, device = torch_dev)\n",
    "                zeros_vec[:, 4*n_qubits + 2*ii+1] = torch.pi/2       #element to modify\n",
    "\n",
    "                # Plus and minus shifts\n",
    "                plus_train_FM = FM_input_weights+zeros_vec\n",
    "                minus_train_FM = FM_input_weights-zeros_vec\n",
    "\n",
    "                # calculate\n",
    "                pos_term_u, pos_term_v, _, pos_term_G = get_predicts(qnodes, X_batch, plus_train_FM, \n",
    "                                                        FM_train_weights, NN_train_weights, TL_train_weights,\n",
    "                                                        n_layers, FMrot_seq, NNrot_seq, [0, 1, 3])\n",
    "                neg_term_u, neg_term_v, _, neg_term_G = get_predicts(qnodes, X_batch, minus_train_FM, \n",
    "                                                        FM_train_weights, NN_train_weights, TL_train_weights,\n",
    "                                                        n_layers, FMrot_seq, NNrot_seq, [0, 1, 3]) \n",
    "                \n",
    "                par_y_u += (FM_input_weights[0][4*n_qubits + 2*ii]/2)*(pos_term_u - neg_term_u)\n",
    "                par_y_v += (FM_input_weights[1][4*n_qubits + 2*ii]/2)*(pos_term_v - neg_term_v)\n",
    "                par_y_G += (FM_input_weights[3][4*n_qubits + 2*ii]/2)*(pos_term_G - neg_term_G)\n",
    "\n",
    "        # Partial of p\n",
    "        par_p_u = 0\n",
    "        par_p_v = 0\n",
    "        par_p_w = 0\n",
    "\n",
    "        for ii in range(n_qubits):\n",
    "                # Zeros vec\n",
    "                zeros_vec = torch.zeros(FM_input_weights.shape, dtype = torch.float32, device = torch_dev)\n",
    "                zeros_vec[:, 6*n_qubits + 2*ii+1] = torch.pi/2       #element to modify\n",
    "\n",
    "                # Plus and minus shifts\n",
    "                plus_train_FM = FM_input_weights+zeros_vec\n",
    "                minus_train_FM = FM_input_weights-zeros_vec\n",
    "\n",
    "                # calculate\n",
    "                pos_term_u, pos_term_v, pos_term_w, _ = get_predicts(qnodes, X_batch, plus_train_FM, \n",
    "                                                        FM_train_weights, NN_train_weights, TL_train_weights,\n",
    "                                                        n_layers, FMrot_seq, NNrot_seq, [0, 1, 2])\n",
    "                neg_term_u, neg_term_v, neg_term_w, _= get_predicts(qnodes, X_batch, minus_train_FM, \n",
    "                                                        FM_train_weights, NN_train_weights, TL_train_weights,\n",
    "                                                        n_layers, FMrot_seq, NNrot_seq, [0, 1, 2]) \n",
    "                \n",
    "                par_p_u += (FM_input_weights[0][6*n_qubits + 2*ii]/2)*(pos_term_u - neg_term_u)\n",
    "                par_p_v += (FM_input_weights[1][6*n_qubits + 2*ii]/2)*(pos_term_v - neg_term_v)\n",
    "                par_p_w += (FM_input_weights[2][6*n_qubits + 2*ii]/2)*(pos_term_w - neg_term_w)\n",
    "\n",
    "        # calculate the square of each equation\n",
    "        sq1 = torch.abs((par_t_u*alpha_t + \n",
    "                                    ((u - beta_u)/alpha_u)*par_x_u + \n",
    "                                    ((v - beta_v)/alpha_v)*par_y_u + \n",
    "                                    w*par_p_u)/alpha_u - \n",
    "                                    (f0 + beta*((X_batch[2]-beta_y)/alpha_y))*((v-beta_v)/alpha_v) + \n",
    "                                    par_x_G*alpha_x/alpha_G)**2\n",
    "        sq2 = torch.abs((par_t_v*alpha_t + \n",
    "                                    ((u - beta_u)/alpha_u)*par_x_v + \n",
    "                                    ((v - beta_v)/alpha_v)*par_y_v + \n",
    "                                    w*par_p_v)/alpha_v +\n",
    "                                    (f0 + beta*((X_batch[2]-beta_y)/alpha_y))*((u-beta_u)/alpha_u)+ \n",
    "                                    par_y_G*alpha_y/alpha_G)**2\n",
    "        sq3 = torch.abs(par_x_u*alpha_x/alpha_u + par_y_v*alpha_y/alpha_v + par_p_w)**2\n",
    "        \n",
    "        \n",
    "        return torch.mean(sq1 + sq2 + sq3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cost function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\gilbe\\anaconda3\\envs\\QuantumAI\\Lib\\site-packages\\tqdm\\auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n",
      "  from .autonotebook import tqdm as notebook_tqdm\n"
     ]
    }
   ],
   "source": [
    "from tqdm.auto import tqdm\n",
    "import time\n",
    "\n",
    "# Create cost function \n",
    "def cost_fn(qnodes, X_batch, y_batch, n_layers, weights, gamma, FMrot_seq, NNrot_seq, epoch):\n",
    "    n_qubits = len(qnodes[0].device.wires)\n",
    "\n",
    "    # Indx lengths\n",
    "    idx0 = 8*n_qubits\n",
    "    idx1 = 3*n_qubits\n",
    "    idx2 = n_layers*n_qubits\n",
    "    idx3 = 2\n",
    "    \n",
    "    # Weights split u\n",
    "    FM_input_weights = weights[:4*idx0].reshape(4, idx0)\n",
    "    FM_train_weights = weights[4*idx0:4*(idx0 + idx1)].reshape(4, idx1)\n",
    "    NN_train_weights = weights[4*(idx0 + idx1):4*(idx0 + idx1 + idx2)].reshape(4, idx2)\n",
    "    TL_train_weights = weights[4*(idx0 + idx1 + idx2):4*(idx0 + idx1 + idx2 + idx3)].reshape(4, idx3)\n",
    "\n",
    "    # Transpose X and y batches\n",
    "    X_batch = torch.t(X_batch)\n",
    "    y_batch = torch.t(y_batch)\n",
    "\n",
    "    global cte_data\n",
    "\n",
    "    # If gamma = 1 we don't need to calculate the data loss\n",
    "    if gamma != 1:\n",
    "        predict_u, predict_v, _, predict_G = get_predicts(qnodes, X_batch, FM_input_weights, \n",
    "                    FM_train_weights, NN_train_weights, TL_train_weights,\n",
    "                    n_layers, FMrot_seq, NNrot_seq)\n",
    "        loss_u = ((predict_u - y_batch[0])*alpha_u)**2\n",
    "        loss_v = ((predict_v - y_batch[1])*alpha_v)**2\n",
    "        loss_G = ((predict_G - y_batch[2])*alpha_G)**2\n",
    "        static_loss = torch.mean(torch.sqrt(y_batch[0]**2 + y_batch[1]**2)*(loss_u + loss_v + loss_G))\n",
    "        if epoch == 0:\n",
    "            cte_data = static_loss.item()\n",
    "    else:\n",
    "        static_loss = 1\n",
    "        cte_data = 1\n",
    "\n",
    "    global cte_phys\n",
    "    \n",
    "    # If gamma = 0 we don't need to calculate the physics loss\n",
    "    if gamma != 0:\n",
    "        physics_loss = pde_loss(qnodes, X_batch, FM_input_weights, FM_train_weights, NN_train_weights, TL_train_weights, n_layers,  FMrot_seq, NNrot_seq)\n",
    "        if epoch == 0:\n",
    "          cte_phys = physics_loss.item()\n",
    "    else:\n",
    "      physics_loss = torch.tensor(0)\n",
    "      cte_phys = 1\n",
    "        \n",
    "    # Calculate the physics and data loss (return them individually)\n",
    "    print(f'Epoch: {epoch} | Physics Loss: {physics_loss.item()/cte_phys} | Static Loss: {static_loss.item()/cte_data} | Cte Phys: {cte_phys} | Cte Data: {cte_data}')\n",
    "\n",
    "    # Total cost\n",
    "    cost =  (1-gamma)*static_loss/cte_data +  gamma*physics_loss/cte_phys\n",
    "    return cost"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dataloader"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create dataloaders \n",
    "BATCH_SIZE = 61*30*5\n",
    "train_dataloader = DataLoader(train_dataset, BATCH_SIZE, shuffle = False)\n",
    "test_dataloader = DataLoader(test_dataset, BATCH_SIZE, shuffle = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " # Training NN (gamma = 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tqdm.auto import tqdm\n",
    "import torch\n",
    "import torch.optim as optim\n",
    "import random\n",
    "\n",
    "# Create quantum circuit qnode and parameters\n",
    "n_layers = 4\n",
    "n_qubits = 4\n",
    "\n",
    "# Training lengths\n",
    "len0 = 8*n_qubits\n",
    "len1 = 3*n_qubits\n",
    "len2 = n_layers*n_qubits\n",
    "len3 = 2\n",
    "\n",
    "# Initialize random parameters\n",
    "params = torch.rand(4*(len0+len1+len2+len3), dtype = torch.float32).to(torch_dev)\n",
    "params.requires_grad = True\n",
    "\n",
    "print(params)\n",
    "\n",
    "# Create qnodes\n",
    "device = qml.device('default.qubit.torch', \n",
    "                    wires = n_qubits,\n",
    "                    torch_device = torch_dev)\n",
    "qnode_u = qml.QNode(full_circuit, device = device, interface = 'torch', diff_method = 'best')\n",
    "qnode_v = qml.QNode(full_circuit, device = device, interface = 'torch', diff_method = 'best')\n",
    "qnode_w = qml.QNode(full_circuit, device = device, interface = 'torch', diff_method = 'best')\n",
    "qnode_G = qml.QNode(full_circuit, device = device, interface = 'torch', diff_method = 'best')\n",
    "\n",
    "qnodes = [qnode_u, qnode_v, qnode_w, qnode_G]\n",
    "\n",
    "def torch_optimize(qnodes, X_train, X_test, params, epochs, FMrot_seq, NNrot_seq, gamma, verbose = False):\n",
    "    train_cost_vec = []\n",
    "    test_cost_vec = []\n",
    "    best_cost = 1e10\n",
    "    # Optimizer\n",
    "    opt = torch.optim.LBFGS([params], history_size = 200)\n",
    "    for i in range(epochs):\n",
    "        avg_train_cost = 0\n",
    "        # Calculate loss and optimize from train batch\n",
    "        for X_batch, y_batch in tqdm(X_train):\n",
    "            def closure():\n",
    "                opt.zero_grad()\n",
    "                cost = cost_fn(qnodes, X_batch, y_batch, n_layers, params, gamma, FMrot_seq, NNrot_seq, i)\n",
    "                if cost.requires_grad:\n",
    "                    cost.backward()\n",
    "                return cost\n",
    "            avg_train_cost += opt.step(closure)\n",
    "\n",
    "        avg_train_cost /= len(X_train)\n",
    "\n",
    "        # Calculate loss from test batch\n",
    "        avg_test_cost = 0\n",
    "        with torch.no_grad():   \n",
    "            for X_batch, y_batch in X_test:\n",
    "                avg_test_cost += cost_fn(qnodes, X_batch, y_batch, n_layers, params, gamma, FMrot_seq, NNrot_seq, i+1)\n",
    "        avg_test_cost /= len(X_test)\n",
    "        test_cost_vec.append(avg_test_cost.item())\n",
    "        \n",
    "        # Print status\n",
    "        if verbose and i%1 == 0:\n",
    "            print(f\"Epoch: {i} | Train loss: {avg_train_cost} | Test loss: {avg_test_cost}\")\n",
    "        train_cost_vec.append(avg_train_cost.item())\n",
    "\n",
    "        # Save parameters, train and test cost\n",
    "        if train_cost_vec[-1] < best_cost:\n",
    "            torch.save(params, 'best_params_0.pth')\n",
    "            best_cost = train_cost_vec[-1]\n",
    "        torch.save(torch.tensor(train_cost_vec), 'train_cost_0_00.pth')\n",
    "        torch.save(torch.tensor(test_cost_vec), 'test_cost_0_00.pth')\n",
    "\n",
    "    return train_cost_vec, test_cost_vec, params\n",
    "\n",
    "n_epochs = 41\n",
    "gamma = 0\n",
    "\n",
    "train_cost, test_cost, new_params= torch_optimize(qnodes, train_dataloader, test_dataloader, params, n_epochs, [qml.RX], [qml.RX], gamma, True)\n",
    "print(new_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train NN (gamma = 0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tqdm.auto import tqdm\n",
    "import torch\n",
    "import torch.optim as optim\n",
    "import random\n",
    "\n",
    "# Create quantum circuit qnode and parameters\n",
    "n_layers = 4\n",
    "n_qubits = 4\n",
    "\n",
    "# Training lengths\n",
    "len0 = 8*n_qubits\n",
    "len1 = 3*n_qubits\n",
    "len2 = n_layers*n_qubits\n",
    "len3 = 2\n",
    "\n",
    "# Initialize random parameters\n",
    "params = torch.rand(4*(len0+len1+len2+len3), dtype = torch.float32).to(torch_dev)\n",
    "params.requires_grad = True\n",
    "\n",
    "print(params)\n",
    "\n",
    "# Create qnodes\n",
    "device = qml.device('default.qubit.torch', \n",
    "                    wires = n_qubits,\n",
    "                    torch_device = torch_dev)\n",
    "qnode_u = qml.QNode(full_circuit, device = device, interface = 'torch', diff_method = 'best')\n",
    "qnode_v = qml.QNode(full_circuit, device = device, interface = 'torch', diff_method = 'best')\n",
    "qnode_w = qml.QNode(full_circuit, device = device, interface = 'torch', diff_method = 'best')\n",
    "qnode_G = qml.QNode(full_circuit, device = device, interface = 'torch', diff_method = 'best')\n",
    "\n",
    "qnodes = [qnode_u, qnode_v, qnode_w, qnode_G]\n",
    "\n",
    "def torch_optimize(qnodes, X_train, X_test, params, epochs, FMrot_seq, NNrot_seq, gamma, verbose = False):\n",
    "    train_cost_vec = []\n",
    "    test_cost_vec = []\n",
    "    best_cost = 1e10\n",
    "    # Optimizer\n",
    "    opt = torch.optim.LBFGS([params], history_size = 200)\n",
    "    for i in range(epochs):\n",
    "        avg_train_cost = 0\n",
    "        # Calculate loss and optimize from train batch\n",
    "        for X_batch, y_batch in tqdm(X_train):\n",
    "            def closure():\n",
    "                opt.zero_grad()\n",
    "                cost = cost_fn(qnodes, X_batch, y_batch, n_layers, params, gamma, FMrot_seq, NNrot_seq, i)\n",
    "                if cost.requires_grad:\n",
    "                    cost.backward()\n",
    "                return cost\n",
    "            avg_train_cost += opt.step(closure)\n",
    "\n",
    "        avg_train_cost /= len(X_train)\n",
    "\n",
    "        # Calculate loss from test batch\n",
    "        avg_test_cost = 0\n",
    "        with torch.no_grad():   \n",
    "            for X_batch, y_batch in X_test:\n",
    "                avg_test_cost += cost_fn(qnodes, X_batch, y_batch, n_layers, params, gamma, FMrot_seq, NNrot_seq, i+1)\n",
    "        avg_test_cost /= len(X_test)\n",
    "        test_cost_vec.append(avg_test_cost.item())\n",
    "        \n",
    "        # Print status\n",
    "        if verbose and i%1 == 0:\n",
    "            print(f\"Epoch: {i} | Train loss: {avg_train_cost} | Test loss: {avg_test_cost}\")\n",
    "        train_cost_vec.append(avg_train_cost.item())\n",
    "\n",
    "        # Save parameters, train and test cost\n",
    "        if train_cost_vec[-1] < best_cost:\n",
    "            torch.save(params, 'best_params_0_05.pth')\n",
    "            best_cost = train_cost_vec[-1]\n",
    "        torch.save(torch.tensor(train_cost_vec), 'train_cost_0_05.pth')\n",
    "        torch.save(torch.tensor(test_cost_vec), 'test_cost_0_05.pth')\n",
    "\n",
    "    return train_cost_vec, test_cost_vec, params\n",
    "\n",
    "n_epochs = 41\n",
    "gamma = 0.05\n",
    "\n",
    "train_cost, test_cost, new_params= torch_optimize(qnodes, train_dataloader, test_dataloader, params, n_epochs, [qml.RX], [qml.RX], gamma, True)\n",
    "print(new_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train NN (gamma = 0.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tqdm.auto import tqdm\n",
    "import torch\n",
    "import torch.optim as optim\n",
    "import random\n",
    "\n",
    "# Create quantum circuit qnode and parameters\n",
    "n_layers = 4\n",
    "n_qubits = 4\n",
    "\n",
    "# Training lengths\n",
    "len0 = 8*n_qubits\n",
    "len1 = 3*n_qubits\n",
    "len2 = n_layers*n_qubits\n",
    "len3 = 2\n",
    "\n",
    "# Initialize random parameters\n",
    "params = torch.rand(4*(len0+len1+len2+len3), dtype = torch.float32).to(torch_dev)\n",
    "params.requires_grad = True\n",
    "\n",
    "print(params)\n",
    "\n",
    "# Create qnodes\n",
    "device = qml.device('default.qubit.torch', \n",
    "                    wires = n_qubits,\n",
    "                    torch_device = torch_dev)\n",
    "qnode_u = qml.QNode(full_circuit, device = device, interface = 'torch', diff_method = 'best')\n",
    "qnode_v = qml.QNode(full_circuit, device = device, interface = 'torch', diff_method = 'best')\n",
    "qnode_w = qml.QNode(full_circuit, device = device, interface = 'torch', diff_method = 'best')\n",
    "qnode_G = qml.QNode(full_circuit, device = device, interface = 'torch', diff_method = 'best')\n",
    "\n",
    "qnodes = [qnode_u, qnode_v, qnode_w, qnode_G]\n",
    "\n",
    "def torch_optimize(qnodes, X_train, X_test, params, epochs, FMrot_seq, NNrot_seq, gamma, verbose = False):\n",
    "    train_cost_vec = []\n",
    "    test_cost_vec = []\n",
    "    best_cost = 1e10\n",
    "    # Optimizer\n",
    "    opt = torch.optim.LBFGS([params], history_size = 200)\n",
    "    for i in range(epochs):\n",
    "        avg_train_cost = 0\n",
    "        # Calculate loss and optimize from train batch\n",
    "        for X_batch, y_batch in tqdm(X_train):\n",
    "            def closure():\n",
    "                opt.zero_grad()\n",
    "                cost = cost_fn(qnodes, X_batch, y_batch, n_layers, params, gamma, FMrot_seq, NNrot_seq, i)\n",
    "                if cost.requires_grad:\n",
    "                    cost.backward()\n",
    "                return cost\n",
    "            avg_train_cost += opt.step(closure)\n",
    "\n",
    "        avg_train_cost /= len(X_train)\n",
    "\n",
    "        # Calculate loss from test batch\n",
    "        avg_test_cost = 0\n",
    "        with torch.no_grad():   \n",
    "            for X_batch, y_batch in X_test:\n",
    "                avg_test_cost += cost_fn(qnodes, X_batch, y_batch, n_layers, params, gamma, FMrot_seq, NNrot_seq, i+1)\n",
    "        avg_test_cost /= len(X_test)\n",
    "        test_cost_vec.append(avg_test_cost.item())\n",
    "        \n",
    "        # Print status\n",
    "        if verbose and i%1 == 0:\n",
    "            print(f\"Epoch: {i} | Train loss: {avg_train_cost} | Test loss: {avg_test_cost}\")\n",
    "        train_cost_vec.append(avg_train_cost.item())\n",
    "\n",
    "        # Save parameters, train and test cost\n",
    "        if train_cost_vec[-1] < best_cost:\n",
    "            torch.save(params, 'best_params_0_5.pth')\n",
    "            best_cost = train_cost_vec[-1]\n",
    "        torch.save(torch.tensor(train_cost_vec), 'train_cost_0_5.pth')\n",
    "        torch.save(torch.tensor(test_cost_vec), 'test_cost_0_5.pth')\n",
    "\n",
    "    return train_cost_vec, test_cost_vec, params\n",
    "\n",
    "n_epochs = 41\n",
    "gamma = 0.5\n",
    "\n",
    "train_cost, test_cost, new_params= torch_optimize(qnodes, train_dataloader, test_dataloader, params, n_epochs, [qml.RX], [qml.RX], gamma, True)\n",
    "print(new_params)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "QuantumAI",
   "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.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}