{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "7a04d7d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.insert(0, \"../\")\n",
    "import time\n",
    "import numpy as np\n",
    "from pinn import *\n",
    "from grad_stats import *\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "from torch.autograd import grad\n",
    "from torch.optim import Adam\n",
    "from torch.optim.lr_scheduler import ExponentialLR, MultiStepLR\n",
    "\n",
    "from tqdm import tqdm_notebook as tqdm \n",
    "device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a351054a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "seed 1\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfcAAAGiCAYAAAD6No9jAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABmNElEQVR4nO29f5RdVXn//z7nTjJJIBkKMTOTj0kIa4UfEquYUCAgCSLRgFbERUEUoT+sESmELIpGtAwUErWWlWURFGoJFKmsLqDFQpW0mqDfoEIkSgONWiNJadIIxQyRZCZz7/7+MeTOPs8+59l333PvzD133q+su1bO3c/e+zn73pk953nvZ+/IGGNACCGEkLYhHmsHCCGEENJYOLkTQgghbQYnd0IIIaTN4OROCCGEtBmc3AkhhJA2g5M7IYQQ0mZwcieEEELaDE7uhBBCSJvByZ0QQghpMzi5E0IIIW0GJ3dCCCEkgCeeeALvfe97MXPmTERRhH/6p39S7Tds2IAoipzXf/7nfzbNx46mtUwIIYS0Ib/97W/xlre8BX/4h3+ID3zgAzXX27ZtG6ZNm1a9fsMb3tAM9wBwcieEEEKCWLZsGZYtWxZcb8aMGTjiiCMa71AKLTe5VyoV/M///A+mTp2KKIrG2h1CCCGBGGPw6quvYubMmYjj5qm/Bw4cwODgYO52jDHOfNPZ2YnOzs7cbducdNJJOHDgAN70pjfhM5/5DM4666yGtm/TcpP7//zP/2DWrFlj7QYhhJCc7Ny5E2984xub0vaBAwcwe/Zs/PrXv87d1uGHH459+/Yl3rvhhhvQ19eXu20A6O3txZ133okFCxZgYGAAf//3f4+zzz4bGzZswJlnntmQPiQtN7lPnToVAHDC77wVpagEAIisE+cjsQYwMtl/FfpsI0TiOrb+X3Jas4lt20hvN7baku1KH+NIlkfW/7PLZF3ZLqLkR23XjWWZEzER/cQdVomwjZT7keMUSR+jzDLXJWXMnbqaj3p0KLKfOjyRJL0fDwFRqtGKaBlj/EYjxkpZpfZ+PH2ait2WkYWZ7RpT9tjK4rJVJm3LwnbI8ijZUMUqk21VjCiDaFeUl81Qpm1Z2Fas8goOJsuE/7atgWgnyrYFABNZ92P9v2zKeP6VLdXf581gcHAQv/71r/HDjT/A4YcfXnc7+/btwymLT8XOnTsTengjn9qPO+44HHfccdXr0047DTt37sQXv/jF8TO5H/rFVYpKKMWHJnfrl74zQQdM7s519uQeB03unkm3UZO7U5Zd17nXXJO76GdMJne9n7ond9+Ezcm9VmOlrIGTuzWBpMzIme36J/dkW/bk6Z3crZ9pZ3IXP+923Yr4XSYnbFlu/x40ngeVivUzK0e/4thGmbbyuyb7sT+PKLb8r6TXbwaHH34YpuaY3A/9kTht2rTE5N5sTj31VNx3331Na7/lJvd2R/7wy6++ET9ebgShXuSPrf0LKaxP+5ed/IPDsbXuN5K/JIVt4g+SBv5OkL+4o8T9Zv+BkdJQ7bZAcgLxTfR22552tUm3JdepyHvXJnvpf8gfGJoLzqQkfRKTv/2HpvyOSP+tP3YhntSdP2Ct755s1/2eij+UrV/XFTFMsRPFSO8TAEyUtI21IXZ/SBNU7K9tJftnvZkYY8L+EE2pPxY888wz6O3tbVr7nNwJIYQUGAP3r5DQ+mHs27cPv/jFL6rX27dvx5YtW3DkkUdi9uzZWLVqFV588UXce++9AIC1a9fi6KOPxoknnojBwUHcd999ePDBB/Hggw/m8FuHkzshhJDiYipe6cdbP5Cnn346sdJ95cqVAIDLLrsM69atw65du7Bjx45q+eDgIK699lq8+OKLmDx5Mk488UQ8+uijOPfcc+v320NkxiomkUF/fz+6urow/8gFtWnuRoYhs8Npjsau6FZu3YCFbiK0mNTcpYYVoqOLuooG72r32Xq3679nrYLtk0dH17VwZfxj8Xen1P6UMHdD9XmlnSC93kcLau42DdPfIRfFeepq2rgj82TXdTR1+QtdKXc0dkWDlwvoZEi8ktDc9QV07sI3a0GdybOgTvhktVuRC+qcBXRyvYG1oM7S3MuVITz7ymbs3bu3aTr2obni2ad+jKmH179w79V9r+LNJ7+tqb6OBXxyJ4QQUnBa6hm1JeDkTgghpLgYk2/hZWsFrxsGJ3dCCCGFpair5ZtNm0/ueh6Hm5YWKWUyD9YukxpcUmu1N4Jw8+dFCo4gkWoitXznOznyRln4FEPmuVdS/z987VtvYNWV6xaUHPkgfV5oir6Ngmw9XNPjXzeo2TZxPz5b5/ulGdeux8sxHa1fRWE6e7aO7mjhSl1fnwm926ujW+3Kn1+5uYwsT/gkN6IRuewmexMbTUeX+rbU0aXebdd1NrFR2qpI3V/R1RP7CMBNk5Wb2ph4pNz6byMzWf3wyT2VNp/cCSGEtDPGlN1NigLrtyOc3AkhhBSY0c9zLwKc3AkhhBQWU6no6ZU11G9HCje5azq5z1b+hSblYTv/Mxa6s9S/bO1calpSn4eViy+3jXTz3rPvLzJin2rnfrL3A5AaXQRlH3pHRx/KLHc0a23vAE/+fLJPfX97n36f5cPwG1q/Ae34CNlbvkYf8uLo0jVX9OwPrzz9ePvUtHEnp9y+9u0Pb+vmek68o6vb+dvOXvLZee7O3vLqIS2eQ1mcupZPUmN3tpi1bWXYWfgYZ/sfi31ESuJjjsv2eRwjtkOjqWOPwSY2RaBwkzshhBAyAsPyaXByJ4QQUlhMpZw4oree+u0IJ3dCCCHFhalwqRR+cpfHFyYv3QNVs6+SGrzcf9k9Lzm7TOrDJqGb+9YMZOeNR45uru2zX/u++u6Zzfr+/cmjWfW6arshdb06tL1Xvu/c9YAcc8Unt93WO241JFdd6rbJdurX3F1tPDuXWtPNnWuvrcksc3V05Tx3p67sJ1uzlrp6QstXdHLftVsmtPJEzrm+P0aHtRBIrjUqibolpbxk/VwNBfzc5GV4E5scC+o4uRNCCCEtRs48d7lhVrvAyZ0QQkhx4Wr5VAoxuSe2eo1kCoswtncLlRlpasuycvZRssP9GqtEhpfllyVSbGX4SsgBzpG2dl1tq1fZj9KOkr4m23XrekLgik9y+PWQfsCRqCEhwRYMpXtpVKjds0pYD9MraWdOu7WH9H0hcCj96KH12kP4gAx7y7JsW1+43L6UPzsyXC6/xbFy9LUMn9tpaW5oXVxH6elsANAhJKwO8fMywSrvsDw+6Bx92zyMybmgjk/uhBBCSGthTCWn5s4nd0IIIaS14Gr5VDi5E0IIKSymMuSc8hdavx1p68ldT5MD9CNg9eNiba3WtZTacu26ua4t16+j+9qq3QdfudKv1Nhz6N3BW8FmMUp/tHuPPbVtg54kwlLUdD9CdPTsdh1bqVErbbl1hWWkaO7yB9y6Vope70amio78X+rQkjiho9euhct25XVJ6N12udS+Na1c2nYgW0eXfTp1leuOeOT/A5XRS4Xjk3s6bT25E0IIaW+Gd6jL8+TOBXWEEEJISzG8iU39T9/cxIYQQghpNZjnnkrxJ3df3nvCVr4hNfmQ42M122xTqe35/2asXxtPdtQC+rZEuXlfnyEadh6C+nEXdQS0Ywu3OXzw6dBau442nl0mqyZztOXaELWbRHmQvu3o5LXr286WrLIfu66ifQNJ3dnVzbN9kjnkUp93dPUou66qhcsycbNJ3TxZJv2X5ROsU6hLVtmB8ihuP8sFdakUf3InhBAyfuGCulQ4uRNCCCksplKBqeTYxCZH3VamcJO7928sO/3FqevZbjPKusiHHQIcq78RGxla10K9Desjz0h5/AtqWWtLz9ryNKPJLbW74G6S7NviV/HA2YZ4BDesnd2vz1Yr94XA7Sst5C3LQ8LlstxJZ3NC4CM+l8TNaiF8O3Usra4WTpf+l2Q43bq2Q+dptol2S8kvTIe4Lom2Yuu6w5pN9g8Z4OcYFXieezqFm9wJIYSQKlxQlwond0IIIYVl+OCYHAvqeHAMIYQQ0mJwQV0qxZjc69R4c31kDdSVg1pq0vcsz8mmYUNRu5asIY+/DOrVs7Voze146/p06Ox2tH58Pmi6uZPipRz969PCk+lguv+2Buymmel6t6aja3p3kG7u08IV/dvxSbRlXzvtCH3bLu4I0MKBpB4utXGnHztFTfYjrzuUsglijJ261hhPGHn/tYOjmArHU+FSKcbkTgghhKRRLg+/8tRvQzi5E0IIKSx8ck+HkzshhJDiwtXyqbTf5B6yW2hAs43SnUMtGyX968pxmCBfuwqtl4a0E6JZh7Sl6eReW0XPlva+4z3tpnyadVynrSx3bbM1a01jl235NGrps13u5pBn9yvvzdHVlW1VY4++bbft6NCy3xAt3Lp29OsAbVzTvmW5UyZ+68dWeSxsS2Jg4gnZ15E1MNHA6E0tzHNPp/0md0IIIeMHPrmnwsmdEEJIYTHlCkyORXGmzMm9NfCEqdVQr3/v2praqa2ttFbTyvKkadVuHVK3kWlbIT40KgTuhH0DQuCaT1oKV2q/SlhbC4k7J5Sp4XI9BK6F6X0paslwOVTbRBhbCY8DenqY7EemiyVsfaH12A6XZ7eTWm5dy340W7csOyTus5VpaJHtkyd8HtmhdhlKLymhdbFXbSyuI+F0NKHDsh35f+XAIEYN5rmnUrzJnRBCCDkEw/KpcHInhBBSWEylnC8szwV1hBBCSIvBJ/dUCj+5O7qn0UpD6nps69TK3Vq1a+NjtYWppn/n8UnTnX0+aSleIba+dDBbnZS6f55tVbXUMp9mHSfKsttJL89uN0QL144YdTRqeT/Ktqo+LTyKs8s0bdyXdhaJjkO08bhDKZPXlo4eS53csRXauGXv1O2QWvjIta2Fp9t2WGXCdsJEcT0hszy2/j+4/wBGDWruqRR+cieEEDJ+YZ57OpzcCSGEFJeKGX7lqd+GcHInhBBSWEylkvPJnZr72GH9YeXXyRXdNiBH3qep16uj+3TnkPxtTZ93tzfN7sebrz0GW6X6tmvVdWh9TEO2VS0p/nvzuRO53/pnV2+euLOewKONh2yrmtgqNa5/W1VvnnhAXXVbVUWfd3LIRd1I9mvniTsau9TCre+TU5Z97WjfUkcXWrmtnWv55sPXE1L/P9yP1M0t245OUZbU3OOJk5LXln00YeT/E3/7GkYNLqhLpRiTOyGEEJIGF9SlwsmdEEJIYTHlnHnuPM+dEEIIaTH45J5KW0/uIRr78LWmQ2fXlnuBhxz92UgdXcv11rTlPMeGanq27Ncty64b4oOsG6Kja9q97Dckh1y2rWnfTj8+W6vcd5Sp1Mrtck03B8L2Vw/JP5f7omv95MkpjwJ088jxKc4uc/Zbz9bRtRxzRyd3dHRFG5f55zI/PaGFC918gtDVOyYqZcnr0gRFc7f+f7Djtxg1cj65g0/uhBBCSKuR88nddxpZQeHkTgghpLgwLJ9KG0zuevhZK2vUlqy+UHsyjI3MstRyK97pq5vYVjVgu1NfWD4k1O6mzWWH2nOlnSltueFxZNo2K7Qu7bXQurzWji6VtiHbtcpyLQwvyxu5XauTdmaVu+Hz7NC7G0qXYXllu1bn2NPstDRv2pkdane2dk2G1mP7uqRv9RqLa1ih99iXspaRopZWN7bD8rJMhOFlmD4uTbD+P9LOBJO0aybDee71p7PVU/eJJ57AX/3VX2Hz5s3YtWsXHn74YZx//vlqnY0bN2LlypXYunUrZs6cieuuuw7Lly+v02s/UmIkhBBCikOlkv8VyG9/+1u85S1vwW233VaT/fbt23Huuefi7W9/O5555hl8+tOfxlVXXYUHH3wwuO9aaYMnd0IIIeMVYwxMjtD6obr9/f2J9zs7O9HZmR6BWLZsGZYtW1ZzH1/5ylcwe/ZsrF27FgBwwgkn4Omnn8YXv/hFfOADH6jPcQ98cieEEFJcKibnk/vw5D5r1ix0dXVVX2vWrGmYi08++SSWLl2aeO9d73oXnn76aRw8eLBh/dgU7sndm96mHtuqKfJhW72qR4wG6PO+rVJte20bWMc2pF2vbo7McnfLVW0dQLa/sm7I1q6yrl9HV/xVtHGnT28amuKD3M7V1sJ9+nxAilpIXTXtLEeKmtT9NR3dewxqQgv3pLfZx5467UqNXTsGVdpKzd2yDdnqtZR9fOrwdbZWrmns8trR5zXNvSR8ED7atgAQRda9x/aagCGMGg06OGbnzp2YNm1a9e2sp/Z62L17N7q7uxPvdXd3Y2hoCC+99BJ6e3sb1tchCje5E0IIIYdo1IK6adOmJSb3RiP3JTkkB8j3GwUnd0IIIcWlAKlwPT092L17d+K9PXv2oKOjA0cddVRT+uTkTgghpLCYioHJEZbPU7dWTjvtNHzzm99MvPf4449j4cKFmCCknEbR1pO7L9jhlmfXcDX47Nx1TZ8P0diBpL6taeyAvn2rpsFLjT1ER9dy1Z1+mpSrLstl3Q6Zk21dOpq1onf7cso1Xd3RwrV+Ao5MDclVHy6Psm21/HPxmyIWN2Dnrmua+nB5tq4epLl7ctdtrVzbBna4rqa5Z2vsw7b2kalyW9hsXd2nubu57Nk6utTCs7aFTbUtKT7J60jce2JBRdb/m0yd6WyJ+oHs27cPv/jFL6rX27dvx5YtW3DkkUdi9uzZWLVqFV588UXce++9AIDly5fjtttuw8qVK/HRj34UTz75JL72ta/hH/7hH+r320NbT+6EEELaG4OcqXB1bD/79NNP46yzzqper1y5EgBw2WWXYd26ddi1axd27NhRLZ87dy4ee+wxXHPNNfjyl7+MmTNn4ktf+lLT0uCAwD+vhoaG8JnPfAZz587F5MmTccwxx+Cmm25CxfrLxxiDvr4+zJw5E5MnT8aSJUuwdevWhjtOCCGEVFfL53kFsmTJkmp+vf1at24dAGDdunXYsGFDos7ixYvx4x//GAMDA9i+fXtTd6cDAp/cP//5z+MrX/kK7rnnHpx44ol4+umn8Yd/+Ifo6urC1VdfDQD4whe+gFtvvRXr1q3Dsccei5tvvhnnnHMOtm3bhqlTpzb8BpqzztAfak+W6Wlzel3PdWSX6SH9pFSg22rpbO4pd0niKLtM2zY2ZPvZkDC8W7f+9DYnjG2H+wPC8LKtZqW3hYThnX5828Ta/ShheCAZiveG4Z22oprKALktbHYY3rmOw1Lh7EGXYXi5bSzs0LX8sGTd2G7XEwKX5XZd8eFFjm7SkW0rPzx7i2tpK34buCu7x36rlCJo7mNB0OT+5JNP4n3vex/OO+88AMDRRx+Nf/iHf8DTTz8NYPipfe3atbj++utxwQUXAADuuecedHd34/7778fHPvaxBrtPCCFkXFOA1fJjQdCfXWeccQb+/d//HT/72c8AAD/5yU/w/e9/H+eeey6A4UUFu3fvTuzE09nZicWLF2PTpk2pbQ4MDKC/vz/xIoQQQmrBlCu5X+1I0JP7Jz/5SezduxfHH388SqUSyuUybrnlFnzwgx8EgGoeX9pOPC+88EJqm2vWrMGNN95Yj++EEELGOXxwTydocn/ggQdw33334f7778eJJ56ILVu2YMWKFZg5cyYuu+yyql3aTjxZu/CsWrWqutIQGN68f9asWSFu1U2z9Hq9z/p71bTvRiI/K59+PxbIe4+i7LI85Nk8qkkbTxUOqWlKXT0Iuy0Zdwz4Le2srhbXtoeu/+JJz1TS/w+4i7Vie/Fx0jYSPphKWfQ7oofLum5b2T7J1eGRde3s9FaS9yo1eLt8jPR3zu6pBE3uf/7nf45PfepTuPjiiwEAb37zm/HCCy9gzZo1uOyyy9DT0wNg+Ane3it3z549ztP8IbSTdwghhBANUzYw5RwL6nLUbWWC/tR67bXXEItVp6VSqZoKN3fuXPT09GD9+vXV8sHBQWzcuBGLFi1qgLuEEEKIxaEn9zyvNiToyf29730vbrnlFsyePRsnnnginnnmGdx66634oz/6IwDD4dwVK1Zg9erVmDdvHubNm4fVq1djypQpuOSSS5pyA4QQQsYvppIvnU0qKe1C0OT+N3/zN/jsZz+LK664Anv27MHMmTPxsY99DH/xF39Rtbnuuuuwf/9+XHHFFXjllVdwyimn4PHHH29Kjnte5NdBUwJDbLW6rt6lt2T/UWmEqfw+y7zrpA9JY/v7LKs5cqQwyFYnU/qxKksN2vUpe0tfOf56JE3/tOxjgeXPtTwyWN2ZUsS9IjFQFUUflv3Yv2Acn3Jp98pYOO2KPH3LoDwkjqgVlSvKLl/uMczJO4ytwakoZbLcH3a0jh2VPywezT3ZqdSdlV+bsh2n34mZtpWKopsLP6Q+796Ppe0L7T4W17C3qi0JW5PcqtbR4DOOfHW/xc3j0AYyeeq3I0GT+9SpU7F27VqsXbs20yaKIvT19aGvry+na4QQQogHg3x/S7Tn3M695QkhhBSXRp3n3m4UfnLPEy7XWjOeLWVNwlaWZQeV3ZB3QHhTXDtfSbuq44LMHbPj/Z6OpM/W/93QaEA/akduOqWNzKayx9Ety5YKxA6mThTV3jnV+Zw9kdFSxfJJ9iO2qi0ltnoV7YoPWtuqVmZpOVvMWjcoI6wVeXpex4itzyd7q1oR2UVUEileop+KpbHEQ+Jz70jWjcoj5RVxs3GHCGtbgxp1iHCzc6Kc2Hb14EhIP3JOlBsStgetsuwT4wAgGhocuRDbzcbiVDgcPJC4NHb4XJRVlNPo3BPksq/lFrjyBDnnhDnLPopH/l8+uB+jBTPh0in85E4IIWQcw9k9FU7uhBBCCospD7/y1G9HOLkTQggpLFwtn05bT+6BUnJyy0nHNttalkjd3E4pcjR1uW2klvckd7KU5VZV+YWVR6YaY/ukdpNy55Zu66TCJa/Ltr4tU/lEy7ZtHGWP4XBbQre17KXmXpHb6Vq2JTmmom65kq3lO0ezKuVu3eQbHQnNXdy71NWHlDLlKFnHJ8dW9HtQs82+Lot2IqeudgytLJM6evbRslGHrGulVsrjYX3H0tZ5tKzvKNnY1u+FPl8pSW0/WbdiafKOti+Pi7VtnTJFN5dlQmN3NflO6/8jfQ699luMFjzyNZ22ntwJIYS0OQb50tnac27n5E4IIaS4DO9Ql69+O8LJnRBCSHHJGZZ38l/bhMJN7lq6tix3yjxbvyYPL8zeGnW4bu356Fqf0qdY8dlEcitOWdfSIz3bg9rl7jgkr5NZvUmdvSxsY/EBxZbeLbeMjZWjZaWmXhIfptzq1ZZQ5b07+n1k95P0SfZrf85SN3c0dumzprkr11KP17R+TecHgJLU0aNsW3kdKfq8065V7mjsPm0/oaN71hskdPTsMulHrOjxw7bZGrzU8lVbR5+PFVtdn3euLXunTObT25q8o+XLXHxFn3dsszX42Co78Brz3Meawk3uhBBCyAh5T3Zrz9mdkzshhJDCUql4DniqoX47wsmdEEJIcamYfLo5NfdWRexbrR1p6XyG2Rq8c+ymErpxbbN7ke24Gnx2WxWxlZKsa2u+skyuP0jo2zI33fEhW1eXenCIraPP2+sApL4t7yfKruv044xTbWWy3NHjHQ0+e28B+Z3Q9Hu3XcV/T669M04BmntiHYMv9z6gXTf3PnsPe3cdgLVHQUmWyX5sLd/nk/hsldx7px9Lk3f0eCWfXraj6fPDPsUBtoo+H6Dly1z8SOr3ltZv6/MD+5N73zcTPrmn0waTOyGEkHEL89xT4eROCCGksHD72XSKMbkntlXNLgOAyDJwjl71hOnttDM3nS07fO6G2rORZbKum3Jn/19Px7PTw9x+5PXIOzLVzQmtO/1a/Thb4maHrmX6mvQpEcb22LqSRLYPmk++cL995abJBdyPUubtxwn319YO4Em586UBJmx1+UJL+5Mygx7+l7bZcoBvO+BICffLkL7sR08DlP0o4X9FOnBSHsVvYy1dT6b2uVKBnXInJQcpHWSn/cXerXfTU/sGBwYxWnATm3SKMbkTQgghKTDPPR1O7oQQQgoLn9zT4eROCCGksPBUuHTabnJPbj+rpMkBujju2bpWa0ZqsdpyTJ8Gb1to6XiyLZ8e7PpolSlpc24/HlurLfU4W/jWDHi0ZSWNTvc/WaZp2JrO7+vHTbmT6wKseh4dvaSlPCq2jk9K2p8s1/z12ippf65P4vOIFf8Djtx1vi8NTQMMSeWz+vBttRuwHbCr9deXBhiypS+Q1P5t298OHsRowbB8Om03uRNCCBk/MM89HU7uhBBCik2bPn3ngZM7IYSQwsKwfDrFm9yFeKZ+MI7E6wjpmRifPqw0I7X+ENxe7dx73aeseuntamW137uvLOGzo+UrtiHtyrKAfnwjGgXo2yGau9tPbe3Icu/2v4qO3sh1DFpOf8i6hoZuHWz1o60RkLayXFsjIOuWpE4ubQOO681z1K+6pa9T1/od490PIPva1vb3D42e5s6wfDrFm9wJIYSQ1+GTezqc3AkhhBQW5rmnU/zJXYl2eg6Bg1Hq+lZoqFWVQhkyHi1CAvp5ZIWG+RBkHdZPkOhgMkuC2vWFwJP95Ehj9F4H1K1TktDC+8Pt1l5XSz8M2/63dgliuJ/skP5YnEjo1JUpg4p0oJ0qKOu6qYjZ7QLZUsGBcv0/v8HkfHJv18V4xZ/cCSGEjFsYlk+HkzshhJDCwgV16XByJ4QQUlx4nnsqxZjcdXF8TKj7+zBGtxLkr8fH0flZCOylUYsZAprJ81Hq7tbesncdQI5hiZR1AEHtKO2mWyh1E//3rGMIsM21xbL1/1jkkjVybUIiNTFArw9JRfQda1xr3cEKt58da4oxuRNCCCEpMCyfDid3QgghhYVP7ulwcieEEFJY+OSeTuEn9yAlMEC7D9EYo6A1AaMjuufRSPP4mK9fhRZcd9ESNHDjBHmsscn4f9Y7tbabj+y2jDIWBuI3uGfcEqW+MVb2Qsi1vkCZdJq190HI+gJZbuv8Q2Yo059GY5Dzyb1hnrQWhZ/cCSGEjF8Ylk+HkzshhJDCUqkAlRyBPYblWwTvZ2jsEJNno1E11NvIEH7tPmnhZ+/9eJKB6m9XbsBZW7u+tvV+tT5T2lJO8QuSWAL7rbcfl/r7rZdGSiiNDb3XivitHOCCCXhcc0L6TrkiB6j9hLWrj7G0zW5ba9fbZ6SVj/y/rHjaaPjknk7hJndCCCHkEJzc0xn9xwVCCCGkQRxaLZ/nVQ+333475s6di0mTJmHBggX43ve+l2m7YcMGRFHkvP7zP/+zzrv2wyd3QgghhWa0H74feOABrFixArfffjtOP/10fPWrX8WyZcvw3HPPYfbs2Zn1tm3bhmnTplWv3/CGNzTNx+JP7kamalj6tqNf59HRlW0ZFR9ytYswvTuZ0uJubpnZimfLSc0HXz/6WAhbVTevv11nH0/NNqBc89dfV+vX51PA91jtR7YbYOu59zCKHkBszIqsRq4DSLZb/5oBeW9aW7ZPQ5VBAD+oxb3cjMWCultvvRV//Md/jD/5kz8BAKxduxbf/va3cccdd2DNmjWZ9WbMmIEjjjiiTk/DKPpPFSGEkHHMIc09zwsA+vv7E6+BgYHU/gYHB7F582YsXbo08f7SpUuxadMm1deTTjoJvb29OPvss/Hd7363IfefBSd3QgghhaVRk/usWbPQ1dVVfWU9gb/00ksol8vo7u5OvN/d3Y3du3en1unt7cWdd96JBx98EA899BCOO+44nH322XjiiScaOhY2xQ/LE0IIGbdUzPArT30A2LlzZ0IP7+zsVOtJacoYkylXHXfccTjuuOOq16eddhp27tyJL37xizjzzDPr9FynEJN7Yrg8+nZSZw/NwVY0a03b9+joMUqZZa6+na1hh2jjIVq4T/vW9G6fvp30yaMlR5pPsq7ik3PsZohunt2vT792xrzeutJ/ra5P93faUuyVfrw6v1bXp88HrAvw3W+y2QYFJhu6vqAFyJH7ZWoUqA+WDwD/V3c3QTQqFW7atGmJyT2L6dOno1QqOU/pe/bscZ7mNU499VTcd999Qb6GwLA8IYSQwnLoyT3PK4SJEydiwYIFWL9+feL99evXY9GiRTW388wzz6C3tzes8wAK8eROCCGEpNGosHwIK1euxKWXXoqFCxfitNNOw5133okdO3Zg+fLlAIBVq1bhxRdfxL333gtgeDX90UcfjRNPPBGDg4O477778OCDD+LBBx+s33EPnNwJIYQUlrHYoe6iiy7Cyy+/jJtuugm7du3C/Pnz8dhjj2HOnDkAgF27dmHHjh1V+8HBQVx77bV48cUXMXnyZJx44ol49NFHce6559bvuIfIhCRXjgL9/f3o6urC/CMXoBQPa9UJHd2ruYfohtl1fVq4riVr+nwpWZZDR3faUvX5pI+JdQBejT1bt3XLkn8vhmjhSZ913dm9LmWXhWjhciyi7DF1tHG1H72u6oPUju3x97Wr5fg7dZW2Yp+2b68Z8Gnsnn6z2pVIn9T1BIG6eUi/9bbjrVpcrX9waD8e2LQce/furUnHrodDc8X1x5yHSaUJdbdzoHwQt/zy0ab6OhbwyZ0QQkhhqSBnWL5hnrQWnNwJIYQUFh4ck07hJ3ft2NaQMLy0923BqoXlnbp2yNjbrgy1j9jHUXYY3vHfYxsnfPKEvB0pQbl32a8SAldtlRD38KXWj6duIoRfe1hbhsfdcHlI+Lz2ELha19OuUzeuvW7i2hMCV8Pymg+A/LGs3Sdf2Nr+OLwpgwEhfkeFq1MOCI26h/hUbzuhZLQVHRy9qcXkXFDHyZ0QQghpMfjkng4nd0IIIYWlbIZfeeq3I5zcCSGEFBaDfEe+tuncXrzJPeS4y9C66pasquaerbHL8ljo11Isk+WxptcrurpXn7dS1mS7sUxnc/RVq24OHT2KxdfP6icWmrqmmzs+StuSbEvRwh0tX9F4S3paY1Kz9mwDG6KFl+LMMr+2r+jomq7u+Ju8DFpfoGnWTl3FNo/2rbUr2/am9tVn6011C7q/2vpMLw9J7avRbHCo9jZzYowJOi43rX47UrjJnRBCCDkEw/LpcHInhBBSWBiWT4eTOyGEkMIyFnvLF4E2n9xr19j99jI/XdkqVT3ytXaNXbYVkrsudXOnXzXPXfgUK1vKKmVOudD2NF1d1dQBVxu3NG1f/nmiPERHV7aBTfVR0cY1Hd2rWSfaTRa52rHwWdGSXf9DcuI1H4SPzjgG1E3o80qZKHe/P7XXDdKspa2W0++T55vlk0TdCTlgvYHd5IBcW9Q8KsagkkM3z1O3lWnzyZ0QQkg7w7B8OsHnub/44ov48Ic/jKOOOgpTpkzBW9/6VmzevLlaboxBX18fZs6cicmTJ2PJkiXYunVrQ50mhBBCgOGwejnHi2F5AK+88gpOP/10nHXWWfjXf/1XzJgxA//1X/+FI444omrzhS98AbfeeivWrVuHY489FjfffDPOOeccbNu2DVOnTm20/6g5NwNwTpTTW/Wc+GWXeU5V09rxhb30bW61a8+WuNBC4J66cXYaXUh6m+w3EZZ3QvaiHzV8HhDSl1u9aiewOSFuTwjc7rckpQ+lrm97Vk1WcBQJxadGhsvjEKlAKffce6T5HxLG9qXGJVITZV3FJ5/Sp4bwlXYlIf146qq/ggLaTe6gHPzcWDfm9X956rcjQZP75z//ecyaNQt333139b2jjz66+n9jDNauXYvrr78eF1xwAQDgnnvuQXd3N+6//3587GMfa4zXhBBCCLigLougP68eeeQRLFy4EBdeeCFmzJiBk046CXfddVe1fPv27di9ezeWLl1afa+zsxOLFy/Gpk2bUtscGBhAf39/4kUIIYTUwqEFdXle7UjQ5P7LX/4Sd9xxB+bNm4dvf/vbWL58Oa666irce++9AIDdu3cDALq7uxP1uru7q2WSNWvWoKurq/qaNWtWPfdBCCFkHFJpwKsdCQrLVyoVLFy4EKtXrwYAnHTSSdi6dSvuuOMOfOQjH6naSa3IGJOpH61atQorV66sXvf393OCbyOkft+8jtR8HqXIm49UXxng6schdWu1DdmiVNqH3HuerV5zaONqCpt3S9kA21J2uZqSFtqPjU9jz9FPvTq698fVp/UfYvQkd2ruGQR9BL29vXjTm96UeO+EE07Ajh07AAA9PT0A4Dyl79mzx3maP0RnZyemTZuWeBFCCCG1kGelfN6ta1uZoMn99NNPx7Zt2xLv/exnP8OcOXMAAHPnzkVPTw/Wr19fLR8cHMTGjRuxaNGiBrhLCCGEjGAa8K8dCQrLX3PNNVi0aBFWr16NP/iDP8CPfvQj3HnnnbjzzjsBDIeVVqxYgdWrV2PevHmYN28eVq9ejSlTpuCSSy5pyg0QQggZv5SNQTnHorg8dVuZoMn95JNPxsMPP4xVq1bhpptuwty5c7F27Vp86EMfqtpcd9112L9/P6644gq88sorOOWUU/D44483Kcd99JB/3Wl570Ys0YicY14tW/HFklpZsi09n962jZxlIiXFVq6RED5IIc5YbWtlQHK1ipO7Lu7druvsSSD7SdY1lbLlkrjXivg8rKZMGckykY8Oq66Rx7aWxb3Kuna5XIeStBR6qyeglhg3JTcdgInl91bJR4+lV/ZASR+0uvKz8tVFJs7XwGrL0ZXVrWqFE0q7stw47WZfOz5p1+LmjPRR0+Q995Mo9erzI9byx7fePPfRXKXGHerSCd5+9j3veQ/e8573ZJZHUYS+vj709fXl8YsQQgjxY3Keyd6mszv3lieEEFJYyjCIc8zQ5Tad3dtgctdifjJ05cSbM9tyQ+syBG61Jf9qdEJkdrmIAwvkbkmJ09tEP5VoKGlrRj7OCkSZ9En5SzcS4yai2oiikTfkKXYmShrb5UbE/CJ5bW8/65NBZPzQjgmW5Och61o+R1JG0LaqFbZyu1Z3oEb+q4T7AcAoW9Wqp7fJULpn+9bE56OdGAckxkY9MU766AuBa+lt2ha4oq4bLtdS7AK2tZU+BaSoGZ9trP1+ym7X8aPWlDRZL8XWGcdEVVk3+/eG/eu0MjB6cXmG5dNpg8mdEELIeKVsDGIuqHPg5E4IIaSwcBObdDi5E0IIKSwMy6dTuMk9JCXNV1f/WGWZliMidHST1KErkZWmJfyVcpfjoy3tSz1b9mMZy/Q7N+XO9kN+DWJxJbXykXIjdP8oSrZl7HuXR8uKPDR7bLSjYwEgqoi0NEsbj8pSJ8y2ddqV6W6Jo0zFugtNCwcSer2RR8s6Wqato+vHwyba0toBVP1e9UFcS/9VPdir5StjEbJ1rVcLr337XOc3gapv59DntXY86wA0bTysX09du0/fr9eMo3Erg6OnuTMsn07hJndCCCHkEAzLp8PJnRBCSGFhWD4dTu6EEEIKC8Py6RR+cpdbNkaWSOTq8wJRN/kZS81ICl5lpcxx0vJJ6rZJH2JIHd3eklXJtUdSo9P07OHrkX4qjm2yn4rQ5GJbVxc+RTgofBoplznxctwSee5SU5f+C00+US41drlFrnU/sbMlrtKv066u7SeP0lS0fFnXp+VH2T7JPHetrpuTreXT167Pu5p0QF2JujYhQHMPaRcIOho3aM1AgA/uOoCAtkL6rbUswAdzcFBvp4EwLJ9O4Sd3Qggh4xcDk1hMXE/9doSTOyGEkMJi4G4SGlq/HeHkTgghpLBUcj6556nbyrTh5G5/UFLD8uSu23qTk1SanefuavvavvQyZ1z6KPRve295536EbaItj2Zt3Z/UZaXmLn1MHhuq103sja9o7NJnx1+pd5c1/VvvJ7mHvScnXrOVaFq/o2cr/mu6vyh3NHWPT9o6AK/WX2u7IZq6bEvq6NJUzXOvX0t29yiovS39eFWlHc+9qser+sa0Vh9qKbfx+fw6Zui12tvMCSf3dNpwcieEEDJe4IK6dDi5E0IIKSzMc0+nEJN7ItCupq8hEYl3TyfUw/R2qdzq1a2bMBaW8rqSuNJtRbjTOtrUFxZOlutpW4mtXp10Nu+ek1ZdT0jfDtd6UgYTPklbT/i/Vh9kuSsNSGqXCjQ//GHU2n1qmFQgTT0SRa3t+CQI1zwg/G8XBd57zX0GtPN6YwFN+b5v9bXrVq2/brKhAH8tyuUDjem/BiomZ1ieee6EEEJIa8En93Q4uRNCCCksXFCXDid3QgghhYUL6tIp/uSuafCO7JS9Va1bql0JvU/0Yzy6etK4dlu3LMDWyPLE/prZ/qW0pWudcs2AVc+XjqRp7o7/AT4o/bhlni1mFduQfoLayuVD7ZppnntNthOq0wa0HTCmiXp1asfpPjSuraB+G6WjjwIHK6O3/SxJp/iTOyGEkHELw/LpcHInhBBSWLigLh1O7oQQQgpNu07QeSjc5K6ktb/+xoiFL33ROS5W02IdfV5r3HN0Y7ZpSnGIzhai1yuYEK04wDagXe/A5PAjWc93XG/Sul7q9a/ZbSUZGy3Zpki6chpmzHKm5b4c9dGohWlDZqgR7tTcLxfUuRRucieEEEIOYaKg5xG3fuNcaSk4uRNCCCksfHJPp/CTux7yDv3Qsu3z/GWodzn2YUh329Gx8WO0CJMHmoOUhIpG034hNnBccrVUrx8B1XzfwjAXNEkuoJUcg2b3U8boheVJOmMvshFCCCF1Yhrwqofbb78dc+fOxaRJk7BgwQJ873vfU+03btyIBQsWYNKkSTjmmGPwla98pc6ea4OTOyGEkMJiMKK71/Wqo88HHngAK1aswPXXX49nnnkGb3/727Fs2TLs2LEj1X779u0499xz8fa3vx3PPPMMPv3pT+Oqq67Cgw8+mOveNSIzdss7U+nv70dXVxfmH7kApbg0/OYohVHHJFjbAiHi5q2+bk0Yls8Pw/L5O23rsLwZwjO/2Yy9e/di2rRp9TeqcGiuePORC1CK61eYy5UhPPt/m7Fz586Er52dnejs7Eytc8opp+Btb3sb7rjjjup7J5xwAs4//3ysWbPGsf/kJz+JRx55BM8//3z1veXLl+MnP/kJnnzyybp91yiG5j5KvwhH488cLXUvd9tN878xk2Ejp9SxmBvz+O/bDrj+dptpr1nLLYkb06c8fjikLX1bYU+/uerW3q5qG7Q9s8+n2m1j29Z3uq1yP3bZwUrxNPdZs2Ylrm+44Qb09fU5doODg9i8eTM+9alPJd5funQpNm3alNr2k08+iaVLlybee9e73oWvfe1rOHjwICZMmJDP+RSKMbkTQgghaeTOhRuum/bknsZLL72EcrmM7u7uxPvd3d3YvXt3ap3du3en2g8NDeGll15Cb29v/f5nwMmdEEJIYYle/5enPgBMmzYtSEKQkRZjjH74Uop92vuNgpM7IYSQwhIhyrWOJvQPg+nTp6NUKjlP6Xv27HGezg/R09OTat/R0YGjjjoqzOEaKf7krmiXwXpknn1ia7RspFZcr/7Y2H7qr90onXbYvt5+fDpnfbYhZb62nYN+1aNZ8/hYuz4cS1vFOEh3FqVxjna1sZDtNureHf+VduTn6PM/Vvpx2s6ol1a3XlvNp4HKQdW/xhIh74qYECZOnIgFCxZg/fr1eP/73199f/369Xjf+96XWue0007DN7/5zcR7jz/+OBYuXNgUvR1gKhwhhJACE5ko9yuUlStX4m//9m/xd3/3d3j++edxzTXXYMeOHVi+fDkAYNWqVfjIRz5StV++fDleeOEFrFy5Es8//zz+7u/+Dl/72tdw7bXXNmwcJMV/cieEEDJuiRDDewCUp34oF110EV5++WXcdNNN2LVrF+bPn4/HHnsMc+bMAQDs2rUrkfM+d+5cPPbYY7jmmmvw5S9/GTNnzsSXvvQlfOADH6jbbx/FyHO38XibCKN67ywgXJgjfzXo5DTPtVYaEhYOCgPnsK3XJ28Y2wlp1uejE2bU2skRyvWFm+tNcwoOo9Z5f/LXnxZSDrnXNB9rrauFiPPYAsmxcO9d1rXb9YTao+wy6ZPmo1Om+OSUKbYlcbOarSy3/3+gfBCf3fboqOS5n3TE6ShFOfLczRCe+c3/11RfxwI+uRNCCCks8ev/6sW0qTrNyZ0QQkiBGd0FdUWBkzshhJACky/PnZN7i+JoQCa7VLf1tNsw3dyX/pLdb4i+3SjtO7U8imq3DdC3tZSiPGlCWt2gNCHF31SfErZN8ilUc09oyQHrAKQeLNrVteQAn0I0Xo//ie9EgO7s+CR1aM0n6YOom/xOSNvkLyTXZ5NZpmnwTjuKru7cq2fcbHu7bP+QAbZhVBiLBXVFoPCTOyGEkPFLo3aoazc4uRNCCCksfHJPh5M7IYSQwjL83J5ncueT+9hhSVF+3VzRqD256moOsHpdu47u6rL1a9ZB+rYQ2kL0eU1rdv1vjBZe8um0AXq91lZITnZpjPRtt1+7ndp9kNc+bVzVrGXdWCnz6NCJz8Nna5U746Jo1D4tebQ06yyN2tdueF1rDYfcLiSoXfH9kv1Y17FlO/FgBaMFw/LpFGNyJ4QQQlKIUEKElA3Paq4/en+IjCac3AkhhBQWau7pcHInhBBSWKIoynUmOsPyBSGh3zVUY9d0W73doFxvpV9XX82u62i8DdLNpX2IrS9PvJTQnXX/G6WN+7R9uzxECx/20f7sdFu7LbcfbX2B6NOTk11S9FVNd9bytWV5JH7wQvRun+YOzac8+dshdRulbyvtpJXbbftsba3cuddSdl1XY/dd2+2OXJjB0XsajlFCnCMsHzMsTwghhLQaMfKdXs6wPCGEENJSxCghjnI8uRs+ubcoteslPkstRUoLp4eE5d0ooww3a3UDQsiesLyWphWSSuaWZdfV/JV1faF1rTykri8ErqWO+VKxtJC+lqIWkjrmC03nCZ+HhMsjLdyfJ3zepHC5L+0sK8ULAOQ8kpBJPKljdlsyPO6G8GW/Vl1PP/aH5w2tJ9qNM8vSym2f7bD8wYH6J9tQmAqXThtM7oQQQsYr3MQmHU7uhBBCCksc5QzL56jbynByJ4QQUmC4oC6Nwk3uvvQ2/RhXXbOOaiyT195UMq3Mk/KV0JI9qWR2XZ8Wrtsm0bXx+rX9EC1/rHR0LZ0tRBsvKXq2rOtLuatXN/f14+rD2e2qW6MGaOxuXfFZKallPs09qW9nl8l2ZVua9i1tQ1LJHN1ffEl0bVzvx/7Cue1kfyCu5q5/mFGin5H/d0yYgNEiimJEzqKDsPrtSOEmd0IIIeQQUc489zxb17YynNwJIYQUFj65p8PJnRBCSGGJohKiHIvi8tRtZQo/uWtJDOEJDtk52a5en51TrtX1byErtOTElqzJdjUN25fnbtv6tHxZHiv9aFq5bwtZdR2AR4cuKdq4pqvLz8ORI23/RJmmhfvqam25mnS2rq5tGZvqk6Kj58k/V/X5oGNDff3Up6Nr7Qy3FaCjO7ZRTWVeW2dL2fp19IR2rpXJuh7N3clzt66j0sgHom3t3Wh4cEw6hZ/cCSGEjF+GD47JE5Zvzzz3XH+yrFmzBlEUYcWKFdX3jDHo6+vDzJkzMXnyZCxZsgRbt27N6ychhBDicCgsn+fVjtT95P7UU0/hzjvvxO/+7u8m3v/CF76AW2+9FevWrcOxxx6Lm2++Geeccw62bduGqVOn5nZ4rAj5286XRle/be0hfdmOFqp2UqKcMH12P74QuG6rpKh5wvDauOVJb3PGzY5uekLgWl1fSD8ZAq89vS30tDPNp5AQfp50Ns0nJ9ws78dqywmtayejOT4IHwuw1WuiXJY1KtTuC8OXSpnldlgeHWWMFlxQl05dd7Vv3z586EMfwl133YXf+Z3fqb5vjMHatWtx/fXX44ILLsD8+fNxzz334LXXXsP999/fMKcJIYQQANW95fP8a0fqmtw/8YlP4LzzzsM73/nOxPvbt2/H7t27sXTp0up7nZ2dWLx4MTZt2pTa1sDAAPr7+xMvQgghpBaiuCP3qx0JvqtvfOMb+PGPf4ynnnrKKdu9ezcAoLu7O/F+d3c3XnjhhdT21qxZgxtvvDHUDUIIIeT1BXX1P32364K6oMl9586duPrqq/H4449j0qRJmXZysIwxmQO4atUqrFy5snrd39+PWbNmhbjV8tgKajO/RsbqKSTUZMSWvW36Xa9i6vxEKmKcpG4e5EOOMbfrNvKzc9pqUD+yLoxSLoxDfvG6/Vg/DxV9gYSRle0jvqULRq5BCUC5V4hjxY3cS9sql9sOy7bs+4lkPzlsURFOWp+PsX9A5A9LE4miDkRR/U/fUTR66wNGk6AR2bx5M/bs2YMFCxZU3yuXy3jiiSdw2223Ydu2bQCGn+B7e3urNnv27HGe5g/R2dmJzs7OenwnhBAyzuGTezpBmvvZZ5+NZ599Flu2bKm+Fi5ciA996EPYsmULjjnmGPT09GD9+vXVOoODg9i4cSMWLVrUcOcJIYSMc6JoOIWh7ld7Tu5BT+5Tp07F/PnzE+8ddthhOOqoo6rvr1ixAqtXr8a8efMwb948rF69GlOmTMEll1zSOK8JIYQQAFFcQuQcHxhQnwfH1MZ1112H/fv344orrsArr7yCU045BY8//njTctylstOov8E8MmGiHyNKjbIdrdT2jFyfIPqpJHR03aeKpQXKP0alBGaXV0RLTj9SdkvckJ6XXE5UFlu9Ct2wbJKjmkTWrb0fSTJcpfdjLIlR5sRLjVTGwey6jpQpY2YJW3fNSqIbu67UaUU/sSKh+nR0u1xrB0hKsb48fXdbVbsh0Y9oK7LHyTkyNXkdW+NYrsh2pI8yHz3b1sldL9eeEx+pW8iKa/lDa9mbcnKgtJx4E5AT79h68t5NqZxaZg4OYLRgWD6d3JP7hg0bEtdRFKGvrw99fX15myaEEEJ0DoXX89RvQ9ozwY8QQsi4IIpyhuUNw/KFIzxknx3adUPv2S27If3a23UjvVZo0YmjStsRyp6bt0O98gQ2J1XJ8cku0yUJu20nNK246PSpSQNIbilbkWFsZytbRb6Isj9Z3wlyIlKqb10rbK3IrhvGrsgxVWwjPSSubfGrbT/rs40DbCMZatfux7k/a+vgkgyXy36tcscnGRJXwvTeLX3tbW71lLs4rmTaOnKFs/VundvcNvC0OfnhJrefHfl/+cAgRo0oyrcojmF5QgghpLXIvaCOT+6EEEJIi8En91Q4uRNCCCksXC2fTuEmdyN1KKnjRkpZAxPn7PQxJavJ6SfyCNpOuSXw+VLU7OCS1F5lepW9Pa07Ksl3pCYfWW2VnLQ/sYbAZGt/Q852rpFVpqfNRUq53HpXk0F9x7gm9Hlxb7HU2KW+am8XKn2SUqZdTymTdR19W35Wiv7tnFyqaNZuP8JHdR2AuFbKvcfQ2p9zyDG0jg/is5T6d0A/yaNxPZq7nWLn6Ney3ezP0tXrpY/Z6wA0bd9n6/icsQ7g4MAoau6l2DmKNgyulieEEEJaCj65p8PJnRBCSHGJYzfEE4Kzo1R7wMmdEEJIceEmNqm0weReu47uHkIYkmktsfKqfdu3JvR5PSfe1dUrmWUynFSxLGJHYxcaY8J/aZvE1b8tbVz6JGxtyU76JMfCvh1Hy5dan+g3qY2H2CZx5UhrfYFSNuyj7Me2lbn32ePk8ylpW7uWL33U2h22jZQypV2pmwsfVB1d0c0Bz9oExSe/Pp+de+/zCYl7z27HqauMN5DU591+ZF1xbefeyzLRbnKc9O90resADgwexKgRR+46ghDkQq42oQ0md0IIIeMWhuVT4eROCCGkwEQ5c9X55E4IIYS0FtzEJpW2m9wTipfUizx7m+t7y4u2kJ1XLfVVTcl3baWOrvRjsvvx6c6xnT/v5Np76ipHy7o+Zpc5mq/drsdW81HT8mVdJ4cf2bYhPkh7f+59fesAnL3kPXVjZb2BO06aD9ntevV5p9/sulpbvnbVPHfZj+OT3Y6uozfOJ8/ZEQH3Y7fl0+eh+B+2rmHk4rUhuRqneUSl2N0DP6Q+w/KEEEJIi8En91Ta808WQggh44M4yv9qIq+88gouvfRSdHV1oaurC5deeil+85vfqHUuv/zy6uY8h16nnnpqUL/FeHK3xl5uueqG3kcM5PGjIWF6uVWt+/GPvOOz9QX4k1f1h9qT5R6fjMkuU9vVfZJo4WaZNldvaN3tU7fVwuWNDMsn+gy4HzdEXPv45/KxUVJBgC0gQuCiri4VeD67KLvMJzMk6+pjnFXP265X6qi9bpB8oUgSIeF+t5+Rsv3OmdNNJO8E3eTJ/ZJLLsF///d/41vf+hYA4E//9E9x6aWX4pvf/KZa793vfjfuvvvu6vXEiROD+i3G5E4IIYSk0cJh+eeffx7f+ta38IMf/ACnnHIKAOCuu+7Caaedhm3btuG4447LrNvZ2Ymenp66+2ZYnhBCSGGJ4ghRHOd4DU/u/f39idfAwEBu35588kl0dXVVJ3YAOPXUU9HV1YVNmzapdTds2IAZM2bg2GOPxUc/+lHs2bMnqG9O7oQQQopLgzT3WbNmVXXxrq4urFmzJrdru3fvxowZM5z3Z8yYgd27d2fWW7ZsGb7+9a/jO9/5Dv76r/8aTz31FN7xjncE/cFR+LC8o+zYGpAQ6KWtUcIxTon3+NhsdI1ab0fT631qd622Po1atQ9Mo6vb1jPcmubr2Cpb1brtpvfhswVcnT3pg6xb/xqCRJ8t2K42DoCubzt+JOrV/13zyawhawi0fupd45C7rlJZW9fg3mvAz4d1MeCee908IuiDVUt9ADt37sS0adOqb3d2dmZW6evrw4033qg2+9RTTw03nzKGxhj1NLqLLrqo+v/58+dj4cKFmDNnDh599FFccMEFar+HKPzkTgghZByTd/vZ1+tOmzYtMblrXHnllbj44otVm6OPPho//elP8b//+79O2a9//Wt0d3fX7GJvby/mzJmDn//85zXX4eROCCGkwORcUFfHY//06dMxffp0r91pp52GvXv34kc/+hF+7/d+DwDwwx/+EHv37sWiRYtq7u/ll1/Gzp070dvbW3Mdau6EEEKKSynK/2oSJ5xwAt797nfjox/9KH7wgx/gBz/4AT760Y/iPe95T2Kl/PHHH4+HH34YALBv3z5ce+21ePLJJ/GrX/0KGzZswHvf+15Mnz4d73//+2vuu3hP7h7JOpGr7skjVfVu2a76+df+5ZB57E55w9JDQ/R54UMuAUu21RxbzdrbTmJL3JAewzR31Vb5Qo2ZTzn6CelT0xob2k+e74hFyJqBRvoQsgaiUf161zzUeK8HK6O3/Wwrp8IBwNe//nVcddVVWLp0KQDg93//93HbbbclbLZt24a9e/cCAEqlEp599lnce++9+M1vfoPe3l6cddZZeOCBBzB16tSa+y3e5E4IIYQcokEL6prFkUceifvuu0+1Mdbi78mTJ+Pb3/527n45uRNCCCkupXj4VS+V9lSniz+5B/zVFRLx1sPwnpbzhNab/Fdkbd2M4taRGTROnsjLKH0gLehB0z6DUfpsmzVuTui6Zb6rIzTq3uttZ8iMZlgeLf3kPlYUf3InhBAyfmnxveXHCk7uhBBCCsuhU9Py1G9HOLkTQggpLnxyT6UYk3uYAJ5JQz/ChvnUnC+WlmpVR2sNbKv5NGtMmzkOTfO5od+DbJo35q3VZzqt4scIYz42o6m5c3JPpRiTOyGEEJIGF9SlwsmdEEJIcWnxTWzGCk7uhBBCikuccwvZMif3lkT9WDx6o7pFo1erHC3Nvc6tOgO0Vv8WpvqGlXnaTtpqm0nk+CwDNxvNbMd7/GVj7tXfTu2bbuhbverthKwiznc/tfsU8pSl+9SYMQxtN4pCNkwJ8FH7/gf0GTIuw22n93uwMghgY1BbdcOwfCqFn9wJIYSMY7igLhVO7oQQQooLn9xTKdzk7nwO6ulaybI8ofaQ0HRIuNzxsUH9uOG1kLqyrPZQtXs/IeHa2m3Vdh2fYnEdMG4BYe08PiVtZVnAuET1+6R/zrF6rfkET1jYbssbwrdtfeFm6/58tmro3RnTANuAdn2SQxQr9yC/Xw3q1ytJ2OVWuwfL+4Ff61UbBp/cUync5E4IIYQcIoojRDkm6Dx1WxlO7oQQQooLw/KpcHInhBBSXDi5p9J2k7uth7kae0gqk0ejttr26+aabe3aeEiKVIgW7tO+1boB+rZPz076pGuXeTTrhE8ejTrZr647u22l65GhtpJEXUdr9Wi+mr6t1fXotLn0YPV+pK3iU+ypq/ikXTvtyDBuwidPu1oIOOjefT4pPmj9Ou0G+Gj9v3LwNWCr3m3DiJFTc2+YJy1F203uhBBCxhFcUJcKJ3dCCCHFJUK+p+/2nNs5uRNCCCkuPM89nUJM7omh9+WJJ8p9GruiewbloyfLYpScnmprZ7h2otz64sl2tbo+bV9vt3HaeBzZXzGPlh/V3q5sK45H+vFryQG6czwyNiF6ttO2r67dj9SvVZ0zwFb049WoA7TkRF0Z5vTcT8K+kbqz0m6Itqz6IOuG6PMePdu9rtEHWe786Gifnc8nWZ5uGw2O4oQZwfUrtH4bUojJnRBCCEmFmnsqnNwJIYQUFz65p8LJnRBCSGGJYm/2qLd+O1L4yT0kl927x7jVVhyUf65r4cm8ap+t9LFUs20cKbZOv9n6dixtlbxxx1bRyqOoQ5Qp+edxR2bZ6wbisnZtPDEWjdTC41Jmuap9A4nQoGub3a9jW/Lo2/Z1ybeGINvW1V5t3Vbmnyu2EOMaokM7GrvmU6Dmrur14rJUu23iWlsjkNZPwielXemH87VU+vH5oPRrj0Nl4CBGDYblUyn85E4IIWQcw7B8KpzcCSGEFBaG5dMp3OTu34I121YLwwMyFC/LstPb3BC+DHMrKV6Oj9lhbjcEnh1u9obaE1vI6iHw2AmnB4TaIy1cnh3+lyFuNbQ+XCGzrhMGLtk+6eHyZJpQdih9uF8lfO6EqpXQuyfUHqkh8IBUrJIecw0Kl2thbE8/ibpaeFzUdU7zUkPrnnaVttQwtqefyLl3pU9deUree/avo9dNs8fUzSqNMss0H+RlwvaAdLCJRFHKTQXWb0MKN7kTQgghh+Dcng4nd0IIIcUlRr7tZxmWJ4QQQloLbj+bTptP7vqHpun32rGtbnntx5P60uakNh4jW0fX0tvcdLZsXV3T1NPqxkoqmexX08JjRxu31wEI/xupoyf0bTGmWipZHGAr7X36vGIbopv70vWgpW1pGrZX21fKHC1f8ylZpOrSzr3n0LcDtHFXh66xHemjT992tonVfMj20aej2z57bZX7SZTJn40mEsUpaxAC67cjbT65E0IIaWuYCpcKJ3dCCCHFhZN7KpzcCSGEFJYojlxJJrB+O9IGk3vAB+NsVZvdlnPsqWNpa016Pn1yq9psPd53rWnsgNh+NiB3Xcs3l+0CSGwN69yP3DY2ztbCVV3d0cIV3VzYu3q2qFtSBElNV/dtCxuSj67lfntsI0Vzd7eqTV4m7INsA3R0j+7stKVo7mouu7OeQKnr8SEkp1zT9oO08ABbQHxcnjUDmm2s3Lt352NZN8u2MnoTJlPh0mmDyZ0QQsi4hWH5VILWCa5ZswYnn3wypk6dihkzZuD888/Htm3bEjbGGPT19WHmzJmYPHkylixZgq1btzbUaUIIIQQYCcvnebUjQU/uGzduxCc+8QmcfPLJGBoawvXXX4+lS5fiueeew2GHHQYA+MIXvoBbb70V69atw7HHHoubb74Z55xzDrZt24apU6c25SbqJc9Wtrqtvh1twtYT0lf7ccLcWj+eE8C0dpTT23y5MpGWv6NcyzC8L0cn0rZv1ULKmr+y3JMOpp7m5j2VzLLVwvCyPCS0Lu29ttlhed0n0a5vO1ptq1clPcyXdpawVdLXUtsK8CkkvS0zdSzVNttHX91Y8UnLwpQNa7skS3O7z5IczybCsHw6QZP7t771rcT13XffjRkzZmDz5s0488wzYYzB2rVrcf311+OCCy4AANxzzz3o7u7G/fffj4997GON85wQQsi4hwfHpJPrtvbu3QsAOPLIIwEA27dvx+7du7F06dKqTWdnJxYvXoxNmzaltjEwMID+/v7EixBCCKmFQ0/ueV7tSN2TuzEGK1euxBlnnIH58+cDAHbv3g0A6O7uTth2d3dXyyRr1qxBV1dX9TVr1qx6XSKEEDLeiBrwakPqXi1/5ZVX4qc//Sm+//3vO2VSozLGZGq8q1atwsqVK6vX/f39ozbBG5jEdYjeHUbF+r++T6LrU+0YM9KPpr8P24704/3L1VSSl4kjbH0+1dePEak0jgYvfbLSHL17RVs+eZ2ybY0Rhe73PFGauA7pR3Yj6tq3HgljzVbae21t/5QyQPjv0cIris9Ou9nXzng7Rzob67/ZKVyA83VKuOT8TMr70z5nx//s8ZdfL/nRQvuKiDcSP3fO91RxSbRsPD5m9Sk/m2YSxZGT3hdavx2pa3L/sz/7MzzyyCN44okn8MY3vrH6fk9PD4DhJ/je3t7q+3v27HGe5g/R2dmJzs7OetwghBBC2vbpOw9BYXljDK688ko89NBD+M53voO5c+cmyufOnYuenh6sX7+++t7g4CA2btyIRYsWNcZjQggh5HWouacT9OT+iU98Avfffz/++Z//GVOnTq3q6F1dXZg8eTKiKMKKFSuwevVqzJs3D/PmzcPq1asxZcoUXHLJJU25AUIIIeOXOHZT9ELrtyNBk/sdd9wBAFiyZEni/bvvvhuXX345AOC6667D/v37ccUVV+CVV17BKaecgscff7yJOe4+MdAukmJTtq1RBcikjiX1JRNVhO2IXiz1O0c3dCS7hMCaKJPSZTK1tRJgK9qtDCVtxZayjkCZqCyu4xrLgMTn4ZNenW1iFdtIfj72lriyH0eQVNYXSFvhU0L39NgmPhDxvTROrvfIfyNZqG1rK9v25Z8rtkbLiRc/Z24OudT6La3Wd9xtbI2pWJdh5M93IicewhbCVlxn/9g5y2YSwyTakR9PYhg9ee3O0oVytk9OP6Xsz0N+9SpKTn8U67+fkqcnjxQeTP4KaS5MdE8laHKvZZFEFEXo6+tDX19fvT4RQgghNcEn93S4tzwhhJDCwgf3dNp8cvekLjlpZ1Z40KkpQrvWtQyBGxnmxkg8TYv+pXVsh8lijwIhQ+9J02TPdr9GhNljcaJcpSJkhsSJcnJbW9GWnd4mbY1s1w6Bi8+mUk5cy5PeEul5JZm6J0ONVnlFi5sCsFPwykpMMqWufaKc9EELy7un2mXbmljJ4QLUML3bD8R1wFa1iX1IPSF8bXtgaarUdeQK5d6dELdysttw29Ybzk7I0v9sH7Rtb33b57o+Wr9zRLtlLZzuaTd5KpzJLAupO3hAke4aDCf3dNp8cieEENLOcPvZdNr0tgghhIwHoijK/Womt9xyCxYtWoQpU6bgiCOOqKlOI05X5eROCCGksLR6nvvg4CAuvPBCfPzjH6+5zqHTVW+77TY89dRT6OnpwTnnnINXX3215jYKF5bPs2WsU1fqbIktG2XdbA1JauxSGzfI1i6ltuzuFmrphsI2VsZC+is198hY6wBEnpBBWdgmyyMrJ8dpV7RVsWJesh+pdyfKHS1fHvEqfbK05LIUSZXjbmU7MsUutnVzj8YuU+ES6wA82n7IkbWabUm5V2FvHJ+Sl7pPoqrqk0dzty+FrbtWQUnlU46H9Y+pqJv4jiSLHJ8SaWeyHXGprC/w+RQpPqlav+OTXBuSbetbm5BM7Ru5GE3NvdVXy994440AgHXr1tVk36jTVfnkTgghpLA06sldnk46MDAwJvdTz+mqaXByJ4QQUlji1w+OyfMCgFmzZiVOKF2zZs2Y3E89p6umwcmdEEJIYYmQ88n99XZ27tyJvXv3Vl+rVq3K7LOvr8+7SO/pp5/Od19yV0jldNU0Cqe5S+SWk/al95RK2Zido+0k0cpc9uwjFd2Na0dsK87WrVJ3ztbKnX5MUhuPYW+rKnRzqY1bbVVEIq/c0tTV1Ufqypx4OW52P45uruj1bv58drsAEFlb5DrrMJS2NO1e1nX71LV9LafcyTGPsjVrN6/a6lfmHctfCIrm7ua5K4KqR9u3+9F0/rS6iWH11s32Se3X2RbWs95A07e1up51AEZbB+C4Lz5Luy2PFp787ESZkl/vjqFsV1zH6eN0cOAgRotG5blPmzYN06ZNq6nOlVdeiYsvvli1Ofroo+vyp57TVdMo/OROCCFk/BJF7t86ofVDmT59OqZPn15/pwr26aonnXQSgJHTVT//+c/X3A7D8oQQQgpLq6fC7dixA1u2bMGOHTtQLpexZcsWbNmyBfv27avaHH/88Xj44Ydfv5+R01Uffvhh/Md//Acuv/zy4NNV+eROCCGksJTiCKUcj+556tbCX/zFX+Cee+6pXh96Gv/ud7+LJUuWAAC2bduGvXv3Vm0acbrquJrc3Rx5WW4XCo1dbE6d1F+FNqbkn0ud2dHy1X70HOaK3Y+ifcu2ItGnzEd3c+SV+5E6aGLNgO5/wtbZE9KjuUe19wOlH1XbD1wHkNDrfTny2joAtV19fYGai+/VrENss8fJty5A1fa1fuUe9iHtausLZD++dQBxtk9qXd/e+I2qq+nxoi11DH1tWe1UDo5eGlmjNPdmsW7dOm+Ou3v0d/7TVcfV5E4IIaS9aPXJfazg5E4IIaSw2Lnq9dZvRwoxudsBC3kkoVHy3SLlCFTZbkpPol0Zerctha0SWpfHw7podZ19MDNt3Va10LRoR56P6bSVnR6mpcJp7Ug/tNQ9b1tKuF9e+/JG80gFiTJvGmC2/KK1FSIN+H0KkBXcxrLb9UgH6merteWRIJLN+GSR2sctLISvtBMiV6Rd1+qTTyoI6bPGtipD+7P7aDB8ck+nEJM7IYQQkkYc5UuFa9MHd07uhBBCikurHxwzVnByJ4QQUlii1//lqd+OFH9yVzR4Z2tazwa0yeNVJdkCviP7K/q8m0si6nr0e9W2xjLpR4jGrrUT6odPr0/YerVxq9zIshD9XvEhRPeXZZ51AI3qxycgaqmJblO1r5dI1vM9CtW/VqFWH4br1v45q/147idkzYDWrvdntl4ffd8JLbUyxAfL/4OVA3o7DSSO3R2NQ+u3I8Wf3AkhhIxbuKAuHU7uhBBCCkscRYhzzNB56rYynNwJIYQUljjKuaCuPef24k3untR1Jw8+WVeUCa1ZavRqP4o+5mjwHk07Wdez3WODCFtEMkqLVbRxcnT00VpA06B+PF8tn34sresl5N7Dxqn+364hZ1TX+7mHje/o3HvovYTeQ7Ve0Pjm+Byt+xkyo3fkK5/c0ync5E4IIYQcIkK+Z6D2nNo5uRNCCCkwpZyr5fPUbWXabnKvPbAOJ40urJ+Auu36p2ELwSEmuQiQzoifMoZGra/h1fI5pKo2/ejbbnInhBAyfuCTezqc3AkhhBQWLqhLh5M7IYSQQtOe03M+ijG516uN1y+p51t9maPffGlOddbL5W/95FvhOvo/zo3sMV9bY/AdqbvHvGmLxep37O41T93G+zw0imsYSnGEUo5k9Tx1W5liTO6EEEJICjzyNR1O7oQQQgpLKYpQyqGb56nbynByJ4QQUliiKMqZCsfJvTXwbuOplHm1Ze2Iyzzt1tPjofJ6t9vMU15/qV6Wp93m1M0zTrl8co6ArfO7l8OH2FdXO8k015hmW/hCpM0aJ63fZo2x93RVbZz0qp7jegP6qdPHwcpobj/LsHwaxZvcCSGEkNfhgrp0OLkTQggpLNxbPp3CT+5OiEkNkctQaEjd7JbCTtrSr7UavrohIctIadetmz1uIXUdH9TQoeezCrhXGRq1uw0JuUrbkJBlnlCobywSdZV7lW2FhHZHzf8m1fV9T9W6okwbN3e8ZT8BP6OetpLtZtuGhP+ddgJ8sG0PVADs0fttFNzEJp3CT+6EEELGL3GULyzPyZ0QQghpMRiWT4eTOyGEkMLCBXXptPnkXr/Gnk/PzrYI1c01zTdERw+y9aRphfhka8Ca/utrN0SLDambR5+X2qs7bg3yqUG20t5ra1269yrr2u16vqdq3ewyWdc7pglbjw8hPjWoH58P2pi7dU2Arbi28ou9ttq19f/9Q6O3pzU193TafHInhBDSzkSv/8tTvx3h5E4IIaSwMCyfDid3QgghhSWGP63TV78dKcbkbsk3ft3cytkM3Ko2bLvH+vRtR9P1+JAhaaW3lchLrr2fUkButOtTskwewlCv7hzqU5ywlWXZ/frGqRRlf84h+rYzLoq+6sstTvjksdU0bKk15tGs7bohGrXjk/hB0+r61gHEcQ4tuU6f5EC5PgX4oNT19aNp7vJ+kLAV31PfvWfcz2sHR+9pmHvLp1OMyZ0QQghJgafCpcPJnRBCSHHJ+eTu3cKvoHByJ4QQUlj45J5OG0zuumaqlYXtma5ovCG53r5caTUvOc4sk/349hhXfVJspb2msTs+OVq4sFXWDPh80nT0kqbt+3Rz69K1FT4pGrarm4u2NP+l7mn3oZRJH6S9V0dX9GHVf7EQxqdZJ8bJpyUHaNaa/z6f7MHx+q9q4eI7UdJsdR+TPulfRnucvLb2+AtbZ4yVH4KE7eDo5blHUb6H7zad29thcieEEDJeYZ57OpzcCSGEFBbmuadTuMnddyyrvqWs/iFqAWYtfO5Lb0uGwLPLAD2NywnXOqHp7H7c8Hm2/94QeINC7VpI3xfuD6mrhZC9skKU/v+0ay2c7rSrhNOdMHxAGFu263weiXBt9palof0kQ7vZZb663rSzOkPtTrhZ2spx08L/Sug6KiGzTNYNaXe47WydRwu9y1C6lu4WaRpQmo/2/Vj/Lw+MXvZ4hHxP3+05tRdwcieEEEIOQc09HU7uhBBCCkuMyIkuhtZvRzi5E0IIKSzcoS6dtp7cfR+ZW56tWefaFlYpy5Pe5urdthau+1TStHBvP9llmq7uSzuzy73rAFSfkrba1q/SJ3VLXKnLKlq+bFvTs4GwtLOkbf1pZyH6dkg6mJbC5WvLKVO0ck0nl+369WzZVra+7fab+OIK2yjbVtPU0Tht3N9u9sKSSAyqmwoXpdp2TBg9zZ1HvqbT1pM7IYSQ9iZCvkVx7Tm1c3InhBBSYJjnng4nd0IIIYUlilyZK7R+O1L4yb1Zn4tfr8/WrDW93tG7FFtAzyl39HlFs3Z9tLTkAI1dlmv588Pl2T5p6wJCctWlvWe3zaRPnmNPVY1atqto8j7NuqTaZuvqPt3f1bBr90nT0fPo845mbd281zahudeuo4ds1+rYB+jo2vassm6QFi7tnXYVbdzbbmyV6R+eo/XHdl3rZxCjCQPzaRR+cieEEDJ+YSpcOk37A+v222/H3LlzMWnSJCxYsADf+973mtUVIYQQQiya8uT+wAMPYMWKFbj99ttx+umn46tf/SqWLVuG5557DrNnz25Gl3UR8vfaaAV9wmz1kLiNL1St9uMJvWvt6rbZIXGfDub2k40WuvaFtZOntdXerlM3yLb29LbQ087qresLl4fYxk44WrFV0uqals4GJB571HQ2+MLl2bZBYXhZ7oTLpW323scytJ6wlWWeL0lWSD8qiQ+jiTAon05TntxvvfVW/PEf/zH+5E/+BCeccALWrl2LWbNm4Y477nBsBwYG0N/fn3gRQgghtRA14F8zueWWW7Bo0SJMmTIFRxxxRE11Lr/88urmPIdep556alC/DZ/cBwcHsXnzZixdujTx/tKlS7Fp0ybHfs2aNejq6qq+Zs2a1WiXCCGEtClRA17NZHBwEBdeeCE+/vGPB9V797vfjV27dlVfjz32WFD9hoflX3rpJZTLZXR3dyfe7+7uxu7dux37VatWYeXKldXrvXv3Yvbs2SibMlB5/U0rSikPfXPCs0YpE+9UIMle7amtPndPTcuuq62G9/ZjdNuKcipcJcqu69v5ruxELK26Shkgd8JDZtnw9cj/3dPxRD+isUQ/Tll2XemDdFJZNO3eu8kul2XqbnZGhOWVuj4fnNuzvvTOCXJBO8dBvFG7bVxpTFjeCY+LH+hcYfmETwHh80aG5bWT7JoVlpdh9zrD8r89MAgAMMZzjGcD2LdvX0Pqy6hxZ2cnOjs7c7UNADfeeCMAYN26dUH1Ojs70dPTU3e/TVstL78UxpjUPXzlAB4a4Odf2dIs1wghhIwCr776Krq6uprS9sSJE9HT04NTF5+Su63DDz/ciRrfcMMN6Ovry912vWzYsAEzZszAEUccgcWLF+OWW27BjBkzaq7f8Ml9+vTpKJVKzlP6nj17nKf5NGbOnImdO3fCGIPZs2dj586dmDZtWqPdbBv6+/sxa9YsjpMHjlNtcJxqg+OkY4zBq6++ipkzZzatj0mTJmH79u0YHBzM3Vbaw2cjntrrZdmyZbjwwgsxZ84cbN++HZ/97Gfxjne8A5s3b67Zr4ZP7hMnTsSCBQuwfv16vP/976++v379erzvfe/z1o/jGG984xurT/DTpk3jD08NcJxqg+NUGxyn2uA4ZdOsJ3abSZMmYdKkSU3vR9LX11cNt2fx1FNPYeHChXW1f9FFF1X/P3/+fCxcuBBz5szBo48+igsuuKCmNpoSll+5ciUuvfRSLFy4EKeddhruvPNO7NixA8uXL29Gd4QQQsioceWVV+Liiy9WbY4++uiG9dfb24s5c+bg5z//ec11mjK5X3TRRXj55Zdx0003YdeuXZg/fz4ee+wxzJkzpxndEUIIIaPG9OnTMX369FHr7+WXX8bOnTvR29tbc52m7VB3xRVX4Fe/+hUGBgawefNmnHnmmUH1Ozs7ccMNN4yp7lEEOE61wXGqDY5TbXCcSK3s2LEDW7ZswY4dO1Aul7FlyxZs2bIlscr/+OOPx8MPPwxgePX+tddeiyeffBK/+tWvsGHDBrz3ve/F9OnTE1K3j8iMRq4CIYQQMg65/PLLcc899zjvf/e738WSJUsADGeX3X333bj88suxf/9+nH/++XjmmWfwm9/8Br29vTjrrLPwl3/5l0H7wHByJ4QQQtqM0T2ZjxBCCCFNh5M7IYQQ0mZwcieEEELaDE7uhBBCSJvRspP77bffjrlz52LSpElYsGABvve97421S2PGmjVrcPLJJ2Pq1KmYMWMGzj//fGzbti1hY4xBX18fZs6cicmTJ2PJkiXYunXrGHncGqxZswZRFGHFihXV9zhOw7z44ov48Ic/jKOOOgpTpkzBW9/6VmzevLlaznEChoaG8JnPfAZz587F5MmTccwxx+Cmm25CpTJyQg3HibQspgX5xje+YSZMmGDuuusu89xzz5mrr77aHHbYYeaFF14Ya9fGhHe9613m7rvvNv/xH/9htmzZYs477zwze/Zss2/fvqrN5z73OTN16lTz4IMPmmeffdZcdNFFpre31/T394+h52PHj370I3P00Ueb3/3d3zVXX3119X2OkzH/93//Z+bMmWMuv/xy88Mf/tBs377d/Nu//Zv5xS9+UbXhOBlz8803m6OOOsr8y7/8i9m+fbv5x3/8R3P44YebtWvXVm04TqRVacnJ/fd+7/fM8uXLE+8df/zx5lOf+tQYedRa7NmzxwAwGzduNMYYU6lUTE9Pj/nc5z5XtTlw4IDp6uoyX/nKV8bKzTHj1VdfNfPmzTPr1683ixcvrk7uHKdhPvnJT5ozzjgjs5zjNMx5551n/uiP/ijx3gUXXGA+/OEPG2M4TqS1abmw/ODgIDZv3oylS5cm3l+6dCk2bdo0Rl61Fnv37gUAHHnkkQCA7du3Y/fu3Ykx6+zsxOLFi8flmH3iE5/Aeeedh3e+852J9zlOwzzyyCNYuHAhLrzwQsyYMQMnnXQS7rrrrmo5x2mYM844A//+7/+On/3sZwCAn/zkJ/j+97+Pc889FwDHibQ2TTvPvV5eeukllMtl53jY7u5u5xjZ8YgxBitXrsQZZ5yB+fPnA0B1XNLG7IUXXhh1H8eSb3zjG/jxj3+Mp556yinjOA3zy1/+EnfccQdWrlyJT3/60/jRj36Eq666Cp2dnfjIRz7CcXqdT37yk9i7dy+OP/54lEollMtl3HLLLfjgBz8IgN8n0tq03OR+CHm2rkk5b3c8cuWVV+KnP/0pvv/97ztl433Mdu7ciauvvhqPP/64egzkeB+nSqWChQsXYvXq1QCAk046CVu3bsUdd9yBj3zkI1W78T5ODzzwAO677z7cf//9OPHEE7FlyxasWLECM2fOxGWXXVa1G+/jRFqTlgvLT58+HaVSyXlK37Nnj/MX8njjz/7sz/DII4/gu9/9Lt74xjdW3+/p6QGAcT9mmzdvxp49e7BgwQJ0dHSgo6MDGzduxJe+9CV0dHRUx2K8j1Nvby/e9KY3Jd474YQTsGPHDgD8Ph3iz//8z/GpT30KF198Md785jfj0ksvxTXXXIM1a9YA4DiR1qblJveJEydiwYIFWL9+feL99evXY9GiRWPk1dhijMGVV16Jhx56CN/5zncwd+7cRPncuXPR09OTGLPBwUFs3LhxXI3Z2WefjWeffbZ66tKWLVuwcOFCfOhDH8KWLVtwzDHHcJwAnH766U4q5c9+9rPqkcz8Pg3z2muvIY6TvyJLpVI1FY7jRFqaMVzMl8mhVLivfe1r5rnnnjMrVqwwhx12mPnVr3411q6NCR//+MdNV1eX2bBhg9m1a1f19dprr1VtPve5z5muri7z0EMPmWeffdZ88IMfZEqOMYnV8sZwnIwZThPs6Ogwt9xyi/n5z39uvv71r5spU6aY++67r2rDcTLmsssuM//v//2/aircQw89ZKZPn26uu+66qg3HibQqLTm5G2PMl7/8ZTNnzhwzceJE87a3va2a9jUeAZD6uvvuu6s2lUrF3HDDDaanp8d0dnaaM8880zz77LNj53SLICd3jtMw3/zmN838+fNNZ2enOf74482dd96ZKOc4GdPf32+uvvpqM3v2bDNp0iRzzDHHmOuvv94MDAxUbThOpFXhka+EEEJIm9FymjshhBBC8sHJnRBCCGkzOLkTQgghbQYnd0IIIaTN4OROCCGEtBmc3AkhhJA2g5M7IYQQ0mZwcieEEELaDE7uhBBCSJvByZ0QQghpMzi5E0IIIW3G/w/hvKH5kCVkkAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# experiment setup\n",
    "lx=0\n",
    "ly=0\n",
    "rx=1\n",
    "ry=1\n",
    "\n",
    "seed   = 1\n",
    "#omega = 6*np.pi\n",
    "\n",
    "print(\"seed\", seed)\n",
    "\n",
    "def kg_equation(x,y):\n",
    "    return x*np.cos(5*np.pi*y) + (x*y)**3\n",
    "\n",
    "x = np.linspace(lx,rx, 100)\n",
    "y = np.linspace(ly, ry, 100)\n",
    "\n",
    "\n",
    "xx,yy = np.meshgrid(x,y)\n",
    "u_sol = kg_equation(xx,yy)\n",
    "\n",
    "X = np.vstack([xx.ravel(), yy.ravel()]).T\n",
    "plt.imshow(u_sol, cmap=\"twilight\", origin=\"lower\",vmin=-1.5,vmax=1.5)\n",
    "plt.colorbar()\n",
    "\n",
    "\n",
    "def u(x):\n",
    "    return x[:, 0:1] * np.cos(5 * np.pi * x[:, 1:2]) + (x[:, 1:2] * x[:, 0:1])**3\n",
    "\n",
    "def u_tt(x):\n",
    "    return - 25 * np.pi**2 * x[:, 0:1] * np.cos(5 * np.pi * x[:, 1:2]) + 6 * x[:,1:2] * x[:,0:1]**3\n",
    "\n",
    "def u_xx(x):\n",
    "    return np.zeros((x.shape[0], 1)) +  6 * x[:,0:1] * x[:,1:2]**3\n",
    "\n",
    "def f(x, alpha=-1.0, beta=0.0, gamma=1.0, k=3.0):\n",
    "    return u_tt(x) + alpha * u_xx(x) + beta * u(x) + gamma * u(x)**k\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "02500524",
   "metadata": {},
   "outputs": [],
   "source": [
    "def sampler(num_r=2500, num_b=100,lx=0,rx=1,ly=0,ry=1,seed=1):\n",
    "    # generate training data\n",
    "    x = np.linspace(lx, rx, 100)\n",
    "    y = np.linspace(ly, ry, 100)\n",
    "    xb = np.linspace(lx,rx,num_b)\n",
    "    yb = np.linspace(ly,ry,num_b)\n",
    "    \n",
    "    xx,yy = np.meshgrid(x,y)\n",
    "    \n",
    "    X = np.vstack([xx.ravel(), yy.ravel()]).T\n",
    "\n",
    "\n",
    "    # X boundaries\n",
    "    lb   = lx*np.ones((yb.shape))\n",
    "    rb   = rx*np.ones((yb.shape))\n",
    "    Xlb  = np.vstack((lb,yb)).T\n",
    "    Xrb  = np.vstack((rb,yb)).T\n",
    "    UXlb = kg_equation(Xlb[:,0:1],Xlb[:,1:2])\n",
    "    UXrb = kg_equation(Xrb[:,0:1],Xrb[:,1:2])\n",
    "\n",
    "    # Y boundaries\n",
    "    lb   = ly*np.ones((xb.shape))\n",
    "    rb   = ry*np.ones((xb.shape))\n",
    "    Ylb  = np.vstack((xb,lb)).T\n",
    "    Yrb  = np.vstack((xb,rb)).T\n",
    "    UYlb = kg_equation(Ylb[:,0:1],Ylb[:,1:2])\n",
    "    UYrb = kg_equation(Yrb[:,0:1],Yrb[:,1:2])\n",
    "\n",
    "    seedc = seed\n",
    "    np.random.seed(seedc)\n",
    "    torch.manual_seed(seedc)\n",
    "\n",
    "    # training tensors\n",
    "    idxs = np.random.choice(xx.size, num_r, replace=False)\n",
    "    X_train = torch.tensor(X[idxs], dtype=torch.float32, requires_grad=True,device=device)\n",
    "    X_rb = torch.tensor(Xrb, dtype=torch.float32, device=device)\n",
    "    X_lb = torch.tensor(Xlb, dtype=torch.float32, device=device)\n",
    "    Y_rb = torch.tensor(Yrb, dtype=torch.float32, device=device)\n",
    "    Y_lb = torch.tensor(Ylb, dtype=torch.float32, requires_grad=True,device=device)\n",
    "    # compute mean and std of training data\n",
    "    X_mean = torch.tensor(np.mean(np.concatenate([X[idxs], Xrb, Xlb, Yrb, Ylb], 0), axis=0, keepdims=True), dtype=torch.float32, device=device)\n",
    "    X_std  = torch.tensor(np.std(np.concatenate([X[idxs], Xrb, Xlb, Yrb, Ylb], 0), axis=0, keepdims=True), dtype=torch.float32, device=device)\n",
    "    U_Train= torch.tensor(f(X[idxs]), dtype=torch.float32, requires_grad=True,device=device)\n",
    "    U_X_rb = torch.tensor(UXrb, dtype=torch.float32, device=device).reshape(num_b,1)\n",
    "    U_X_lb = torch.tensor(UXlb, dtype=torch.float32, device=device).reshape(num_b,1)\n",
    "    U_Y_rb = torch.tensor(UYrb, dtype=torch.float32, device=device).reshape(num_b,1)\n",
    "    U_Y_lb = torch.tensor(UYlb, dtype=torch.float32, requires_grad=True, device=device).reshape(num_b,1)\n",
    "    \n",
    "\n",
    "    return X_train, X_lb, X_rb, Y_lb, Y_rb, U_Train,U_X_lb, U_X_rb, U_Y_lb, U_Y_rb, X_mean, X_std\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "437c82ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "def KG_res(uhat, data):\n",
    "    poly = torch.ones_like(uhat)\n",
    "    \n",
    "    du = grad(outputs=uhat, inputs=data, \n",
    "              grad_outputs=torch.ones_like(uhat), create_graph=True)[0]\n",
    "    \n",
    "    dudx = du[:,0:1]\n",
    "    dudy = du[:,1:2]\n",
    "    \n",
    "    dudxx = grad(outputs=dudx, inputs=data, \n",
    "              grad_outputs=torch.ones_like(uhat), create_graph=True)[0][:,0:1]\n",
    "    dudyy = grad(outputs=dudy, inputs=data, \n",
    "              grad_outputs=torch.ones_like(uhat), create_graph=True)[0][:,1:2]\n",
    "    \n",
    "    xin = data[:,0:1]\n",
    "    yin = data[:,1:2]\n",
    "    \n",
    "    \n",
    "    \n",
    "    residual = dudyy - dudxx + uhat**3\n",
    "                \n",
    "    return residual\n",
    "def u_t(uhat, data):\n",
    "    poly = torch.ones_like(uhat)\n",
    "    \n",
    "    du = grad(outputs=uhat, inputs=data, \n",
    "              grad_outputs=torch.ones_like(uhat), create_graph=True)[0]\n",
    "    \n",
    "    \n",
    "    dudy = du[:,1:2]\n",
    "    return dudy - 0*uhat\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d378fed2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_function(lx,rx,ly,ry,u_sol,out,method,extras=None):\n",
    "    methods=[\"W1 (uniform)\",\"W2 (max_by_mean)\",\"W3 (std)\",\"W4 (mean+std)\",\"W5 (mean*std)\",\"W6 (kurtosis)\"]\n",
    "    plt.subplot(1,3,1)\n",
    "    plt.imshow(u_sol, cmap=\"twilight\", origin=\"lower\",vmin=-1.0,vmax=1.5)\n",
    "    plt.xticks(np.arange(0,101,50),np.linspace(lx,rx,3),fontsize=12)\n",
    "    plt.yticks(np.arange(0,101,50),np.linspace(ly,ry,3),fontsize=12)\n",
    "    plt.xlabel(r\"$x$\",fontsize=15)\n",
    "    plt.ylabel(r\"$y$\",fontsize=15)\n",
    "    plt.title(\"Ground Truth\",fontsize=18)\n",
    "    plt.colorbar(fraction=0.046, pad=0.04)\n",
    "\n",
    "    plt.subplot(1,3,2)\n",
    "    plt.imshow(out, cmap=\"twilight\", origin=\"lower\",vmin=-1.0,vmax=1.5)\n",
    "    plt.xticks(np.arange(0,101,50),np.linspace(lx,rx,3),fontsize=12)\n",
    "    plt.yticks(np.arange(0,101,50),np.linspace(ly,ry,3),fontsize=12)\n",
    "    plt.xlabel(r\"$x$\",fontsize=15)\n",
    "    plt.ylabel(r\"$y$\",fontsize=15)\n",
    "    plt.title(\"Prediction\",fontsize=18)\n",
    "    plt.colorbar(fraction=0.046, pad=0.04)\n",
    "\n",
    "    plt.subplot(1,3,3)\n",
    "    plt.imshow(np.abs(out-u_sol)/np.max(np.abs(u_sol)), cmap=\"nipy_spectral\", origin=\"lower\",vmin=0,vmax=0.2)\n",
    "    plt.xticks(np.arange(0,101,50),np.linspace(lx,rx,3),fontsize=12)\n",
    "    plt.yticks(np.arange(0,101,50),np.linspace(ly,ry,3),fontsize=12)\n",
    "    plt.xlabel(r\"$x$\",fontsize=15)\n",
    "    plt.ylabel(r\"$y$\",fontsize=15)\n",
    "    plt.title(\"Point-wise Error\",fontsize=18)\n",
    "    plt.colorbar(fraction=0.046, pad=0.04)\n",
    "    \n",
    "    \n",
    "    \n",
    "    plt.gcf().set_size_inches(15,5)\n",
    "    plt.tight_layout()\n",
    "    plt.suptitle(\"Klein-Gordon Equation using PINN_{}\".format(methods[method]),fontsize=25)\n",
    "    plt.savefig(extras+\"KGEqn_{}_based_Tanh\".format(methods[method]),dpi=800)\n",
    "    return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "684e9bbc",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "#######Training with#####\n",
      " 1000+100\n",
      "epoch 755/40001, loss=8999.0664062500, lambda=1.0000, lr=0.00001\t\t\t\t\r"
     ]
    }
   ],
   "source": [
    "losses_boundary_global=[]\n",
    "losses_residual_global=[]\n",
    "all_losses=[]\n",
    "lambdas_global=[]\n",
    "list_of_l2_Errors=[]\n",
    "for i in range(6):\n",
    "    for j in range(5):\n",
    "        #sets = [[300,25],[500,25],[600,30],[750,50],[1000,100]]\n",
    "        mm         = 5\n",
    "        alpha_ann  = 0.5\n",
    "        n_epochs   = 40001\n",
    "        method = i\n",
    "        num_r= 1000\n",
    "        num_b= 100\n",
    "        extras=str(num_r)+ \"+\"+ str(num_b)\n",
    "        print(\"#######Training with#####\\n\",extras)\n",
    "        #print(extras)\n",
    "        X_train, X_lb, X_rb, Y_lb, Y_rb, U_Train,U_X_lb, U_X_rb, U_Y_lb, U_Y_rb, X_mean, X_std= sampler(num_r,num_b,lx,rx,ly,ry,seed=j)\n",
    "        net = PINN(sizes=[2,50,50,50,50,50,1], activation=torch.nn.Tanh()).to(device)\n",
    "        lambd       = 1\n",
    "        lambds      = [];\n",
    "        losses_boundary  = [];\n",
    "        losses_residual = [];\n",
    "        params = [{'params': net.parameters(), 'lr': 1e-5}]\n",
    "        #milestones = [[10000,20000,30000]]\n",
    "        optimizer = Adam(params)\n",
    "        #scheduler = MultiStepLR(optimizer, milestones[0], gamma=0.1)\n",
    "        #print(\"training with shape\", X_train.size())\n",
    "        start_time = time.time()\n",
    "        for epoch in range(n_epochs): \n",
    "            \n",
    "            uhat  = net(X_train)\n",
    "            res   = KG_res(uhat, X_train)\n",
    "            l_reg = torch.mean((res-U_Train)**2)\n",
    "            predl = net(X_lb)\n",
    "            predr = net(X_rb)\n",
    "            l_bc  = torch.mean((predl - U_X_lb)**2, dim=0)\n",
    "            l_bc += torch.mean((predr - U_X_rb)**2, dim=0)\n",
    "            predl = net(Y_lb)\n",
    "            #predr = net(Y_rb)\n",
    "            l_bc += torch.mean((predl - U_Y_lb)**2, dim=0)\n",
    "            #l_bc += torch.mean((predr - U_Y_rb)**2, dim=0)    \n",
    "            gpreds=u_t(predl,Y_lb)\n",
    "            l_bc += torch.mean((gpreds)**2)\n",
    "            with torch.no_grad():\n",
    "                if epoch % mm == 0:\n",
    "                    stdr,kurtr=loss_grad_stats(l_reg, net)\n",
    "                    stdb,kurtb=loss_grad_stats(l_bc, net)\n",
    "                    maxr,meanr=loss_grad_max_mean(l_reg, net)\n",
    "                    maxb,meanb=loss_grad_max_mean(l_bc, net,lambg=lambd)\n",
    "                    if method == 2:\n",
    "                        # inverse dirichlet\n",
    "                        lamb_hat = stdr/stdb\n",
    "                        lambd     = (1-alpha_ann)*lambd + alpha_ann*lamb_hat\n",
    "                        if lambd < 1:\n",
    "                            lambd = torch.tensor(1.0, dtype=torch.float32, device=device)\n",
    "                       \n",
    "                    elif method == 1:\n",
    "                        # max/avg\n",
    "                        lamb_hat = maxr/meanb\n",
    "                        lambd     = (1-alpha_ann)*lambd + alpha_ann*lamb_hat \n",
    "                        if lambd < 1:\n",
    "                            lambd = torch.tensor(1.0, dtype=torch.float32, device=device)\n",
    "            \n",
    "                    elif method==3:\n",
    "                        # mean + std weighing\n",
    "                        covr= stdr + maxr\n",
    "                        covb= stdb + meanb\n",
    "                        lamb_hat = covr/covb\n",
    "                        lambd     = (1-alpha_ann)*lambd + alpha_ann*lamb_hat\n",
    "                        if lambd < 1:\n",
    "                            lambd = torch.tensor(1.0, dtype=torch.float32, device=device)\n",
    "            \n",
    "                    elif method == 5:\n",
    "                        # kurtosis based weighing\n",
    "                        covr= stdr/kurtr\n",
    "                        covb= stdb/kurtb\n",
    "                        lamb_hat = covr/covb\n",
    "                        lambd     = (1-alpha_ann)*lambd + alpha_ann*lamb_hat\n",
    "                        if lambd < 1:\n",
    "                            lambd = torch.tensor(1.0, dtype=torch.float32, device=device)\n",
    "            \n",
    "                    elif method == 4:\n",
    "                        # mean * std weighing\n",
    "                        covr= stdr * meanr\n",
    "                        covb= stdb * meanb\n",
    "                        lamb_hat = covr/covb\n",
    "                        lambd     = (1-alpha_ann)*lambd + alpha_ann*lamb_hat\n",
    "                        if lambd < 1:\n",
    "                            lambd = torch.tensor(1.0, dtype=torch.float32, device=device)\n",
    "            \n",
    "                    else:\n",
    "                        # uniform weighing \n",
    "                        lambd = 1;\n",
    "            if(method == 0):\n",
    "                loss = l_reg + l_bc\n",
    "            elif(method == 1 or method == 2 or method==3 or method==4 or method == 5):\n",
    "                loss = l_reg + lambd*l_bc\n",
    "            if epoch%250==0:\n",
    "                all_losses.append(loss.item())\n",
    "                inp = torch.tensor(X, dtype=torch.float32, device=device)\n",
    "                out = net(inp).cpu().data.numpy().reshape(u_sol.shape)\n",
    "                list_of_l2_Errors.append(np.linalg.norm(out.reshape(-1)-u_sol.reshape(-1))/np.linalg.norm(out.reshape(-1)))\n",
    "                if method !=0:\n",
    "                    lambds.append(lambd.item())\n",
    "            optimizer.zero_grad()\n",
    "            loss.backward()\n",
    "            optimizer.step()\n",
    "            #scheduler.step()\n",
    "            print(\"epoch {}/{}, loss={:.10f}, lambda={:.4f}, lr={:,.5f}\\t\\t\\t\"\n",
    "                  .format(epoch+1, n_epochs, loss.item(), lambd, optimizer.param_groups[0]['lr']), end=\"\\r\")\n",
    "        elapsed_time = time.time() - start_time\n",
    "        #print('CPU time = ',elapsed_time)\n",
    "        inp = torch.tensor(X, dtype=torch.float32, device=device)\n",
    "        out = net(inp).cpu().data.numpy().reshape(u_sol.shape)\n",
    "        print(\"\\n.....\\n\")\n",
    "        print(\"Method:\",method)\n",
    "        print(\"pred rel. l2-error = {:e}\".format(np.linalg.norm(out.reshape(-1)-u_sol.reshape(-1))/np.linalg.norm(out.reshape(-1))))\n",
    "        print(\"\\n.....\\n\")\n",
    "        #if j==0:\n",
    "            #plot_function(lx,rx,ly,ry,u_sol,out,method,extras=extras)\n",
    "        #list_of_l2_Errors.append(np.linalg.norm(out.reshape(-1)-u_sol.reshape(-1))/np.linalg.norm(out.reshape(-1)))\n",
    "        #losses_boundary_global.append(losses_boundary)\n",
    "        #losses_residual_global.append(losses_residual)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc0ee37f",
   "metadata": {},
   "outputs": [],
   "source": [
    "list_of_l2_Errors\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa195472",
   "metadata": {},
   "outputs": [],
   "source": [
    "def split(list_a, chunk_size):\n",
    "\n",
    "  for i in range(0, len(list_a), chunk_size):\n",
    "    yield list_a[i:i + chunk_size]\n",
    "List = list(split(list_of_l2_Errors,5))\n",
    "for i in List:\n",
    "    \n",
    "    arr = np.array(i)\n",
    "    #print(\"\\n\")\n",
    "    print(np.mean(arr),\"+_\",np.std(arr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f5c8f8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "List = list(split(list_of_l2_Errors,161))\n",
    "errs=[]\n",
    "stds= []\n",
    "\n",
    "All_list = list(split(List,5))\n",
    "for i in All_list:\n",
    "    errs.append(np.mean(i,axis=0))\n",
    "    stds.append(np.std(i,axis=0))\n",
    "fig,ax = plt.subplots(1,1,figsize=(10,14))\n",
    "colors = [\"red\", \"blue\",\"yellow\" ,\"orange\", \"purple\", \"green\"]\n",
    "labels = [\"uniform PINN\", \"mean based\",\"mean*std\" ,\"mean+std\" ,\"std based\",\"kurt based\"]\n",
    "markers = [\"o\", \"d\", \"s\", \"p\", \"v\",\"x\"]\n",
    "linestyle= [\"dotted\", \"dashed\",\"dashdot\", \"\" ]\n",
    "for i in range(6):\n",
    "\n",
    "    ax.plot(range(0,40001,250), errs[i], colors[i], label=labels[i],linestyle=\"dashdot\")\n",
    "    if i ==0 :\n",
    "        ax.fill_between(range(0,40001,250), errs[i]-(stds[i]), errs[i]+(stds[i]), color=colors[i], alpha=0.2)\n",
    "    else:\n",
    "        ax.fill_between(range(0,40001,250), errs[i]-(.5*stds[i]), errs[i]+(.5*stds[i]), color=colors[i], alpha=0.2)\n",
    "    \n",
    "ax.set_xlabel(\"Epochs\", fontsize=14)\n",
    "ax.set_ylabel(\"Rel. L2 error\", fontsize=14)\n",
    "ax.legend(fontsize=11)\n",
    "ax.set_yscale(\"log\")\n",
    "#ax.set_xscale(\"log\")\n",
    "\n",
    "plt.gcf().set_size_inches(12,12)\n",
    "plt.savefig(\"kg_erros_rms_1e5.eps\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.17 ('newtorch')",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.17"
  },
  "vscode": {
   "interpreter": {
    "hash": "813e009bea45919a754463d31638190cf3c033fdb5fd5c2966a42653051a1217"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}