{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adaptive scaling of total TEC for sub-orbital TEC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Go to directory: /Users/yunjunz/Papers/2021_Geolocation/figs_src/topTEC\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import os\n",
    "import string\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt, ticker, patches\n",
    "from mintpy.objects import timeseries, sensor\n",
    "from mintpy.utils import readfile, utils as ut\n",
    "from mintpy import add\n",
    "from ipynb.fs.full import utils\n",
    "plt.rcParams.update({'font.size': 12})\n",
    "\n",
    "work_dir = os.path.expanduser('~/Papers/2022_Geolocation/figs_src/topTEC')\n",
    "os.chdir(work_dir)\n",
    "print('Go to directory:', work_dir)\n",
    "\n",
    "proj_dirs = [os.path.expanduser('~/data/geolocation/ChileSenAT149/mintpy_offset'),\n",
    "             os.path.expanduser('~/data/geolocation/ChileSenDT156/mintpy_offset'),\n",
    "             os.path.expanduser('~/data/geolocation/KyushuAlos2DT23/mintpy_offset')]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for proj_dir in proj_dirs:\n",
    "    os.chdir(proj_dir)\n",
    "    #!iono_tec.py timeseriesRg.h5 -g inputs/geometryRadar.h5 -s jpl\n",
    "    #!iono_tec.py timeseriesRg.h5 -g inputs/geometryRadar.h5 -s jpl --ratio 0.69\n",
    "    #!iono_tec.py timeseriesRg.h5 -g inputs/geometryRadar.h5 -s jpl --ratio adaptive\n",
    "    #!iono_tec.py timeseriesRg.h5 -g inputs/geometryRadar.h5 -s cod\n",
    "    #!iono_tec.py timeseriesRg.h5 -g inputs/geometryRadar.h5 -s cod --ratio 0.69\n",
    "    #!iono_tec.py timeseriesRg.h5 -g inputs/geometryRadar.h5 -s cod --ratio adaptive\n",
    "os.chdir(work_dir)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "for proj_dir in proj_dirs:\n",
    "    suffix = '' if 'Alos2' in proj_dir else '_S1Bias'\n",
    "    infile = os.path.join(proj_dir, f'timeseriesRg{suffix}_SET_ERA5.h5')\n",
    "    fbase = os.path.splitext(infile)[0]\n",
    "\n",
    "    tbases = [\n",
    "        'TECclr', 'TECclrR69', 'TECclrRA',\n",
    "        'TECjlr', 'TECjlrR69', 'TECjlrRA',\n",
    "        'TECjhr', 'TECjhrR69', 'TECjhrRA',\n",
    "    ]\n",
    "\n",
    "    for tbase in tbases:\n",
    "        tec_file = os.path.join(proj_dir, 'inputs', f'{tbase}.h5')\n",
    "        outfile = f'{fbase}_{tbase}.h5'\n",
    "        if not os.path.isfile(outfile):\n",
    "            cmd = f'{infile} {tec_file} -o {outfile}'\n",
    "            print(f'add.py {cmd}')\n",
    "            add.main(cmd.split())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ChileSenAT149: RMSE\n",
      "    SAR - S1Bias - SET - ERA5               :   18.1 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECclr      :   21.2 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECclrR69   :   19.9 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECclrRA    :   18.5 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECjlr      :   19.2 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECjlrR69   :   18.4 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECjlrRA    :   16.8 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECjhr      :    6.4 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECjhrR69   :    8.7 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECjhrRA    :    7.5 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECsub      :    4.9 cm\n",
      "ChileSenDT156: RMSE\n",
      "    SAR - S1Bias - SET - ERA5               :    5.7 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECclr      :    5.1 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECclrR69   :    5.1 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECclrRA    :    5.1 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECjlr      :    5.1 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECjlrR69   :    4.9 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECjlrRA    :    5.0 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECjhr      :    5.6 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECjhrR69   :    5.4 cm\n",
      "    SAR - S1Bias - SET - ERA5 - TECjhrRA    :    5.4 cm\n",
      "KyushuAlos2DT23: RMSE\n",
      "    SAR - SET - ERA5                        :  271.2 cm\n",
      "    SAR - SET - ERA5 - TECclr               :  138.1 cm\n",
      "    SAR - SET - ERA5 - TECclrR69            :  166.6 cm\n",
      "    SAR - SET - ERA5 - TECclrRA             :  168.1 cm\n",
      "    SAR - SET - ERA5 - TECjlr               :  131.3 cm\n",
      "    SAR - SET - ERA5 - TECjlrR69            :  164.8 cm\n",
      "    SAR - SET - ERA5 - TECjlrRA             :  162.6 cm\n",
      "    SAR - SET - ERA5 - TECjhr               :   55.7 cm\n",
      "    SAR - SET - ERA5 - TECjhrR69            :   88.8 cm\n",
      "    SAR - SET - ERA5 - TECjhrRA             :   92.7 cm\n"
     ]
    }
   ],
   "source": [
    "proj_names = []\n",
    "rDicts = []\n",
    "for proj_dir in proj_dirs:\n",
    "    suffix = '' if 'Alos2' in proj_dir else '_S1Bias'\n",
    "    infile = os.path.join(proj_dir, f'timeseriesRg{suffix}_SET_ERA5.h5')\n",
    "    fbase = os.path.splitext(infile)[0]\n",
    "    fnames = [infile] + [f'{fbase}_{x}.h5' for x in tbases + ['TECsub']]\n",
    "    proj_name, tsDict, rDict = utils.read_ts_files(fnames, print_msg=True, print_max=False)\n",
    "    proj_names.append(proj_name)\n",
    "    rDicts.append(rDict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "save figure to file /Users/yunjunz/Papers/2021_Geolocation/figs_src/topTEC/topTEC_scale_rmse.pdf\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAGACAYAAAAArQ1sAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABqOElEQVR4nO3dd3hUZfbA8e8JPfReJSggIEhoil1QWVYUZMH1p6KACIFFFBSscQEFRFEQKxJEEYOususiESm6S7GACC4KCAhID0gooSUCSc7vj5nESZ8kM3NnkvN5nvuQufVMcpgz9973vq+oKsYYY4yTwpwOwBhjjLFiZIwxxnFWjIwxxjjOipExxhjHWTEyxhjjOCtGxhhjHGfFyBhjjONKVDESkZoiMkVEtonI7yJyWERWiUh/ESntXmeOiHyZxz5WiIi6p/MisltEXhORal7G8GcR+VxEEkTkrIjsEpE4EektImHudZq493+Nx3bpx+yZwz4XuJe9XeBfisnGMwdEZLyI7Mhn3fS/TaqI7BeRuSLS0MtjXSEi80Uk3p0P+0TkSxG5V0TKeqynInKPx+vd7nkP5rDP6e5lueaxKThv88L9t3na47V9ZnihxBQjEWkE/AD0BZ4FOgBXA7OBMUCbAuzuA6A+cCEwDOgDvOlFDGOBRcAe4A6gJXAXEAeMAxrks4u9wJAs+6wP9AD2FSB+41tf4cqHxsDdQHtgfn4bich9wNfulwOAS4DewHu4/s6X5bOLnPKhPHAvrhwzwcM+M/JR2ukAAmgGUA5oq6onPOZvF5EPgLI5b5ajZFU95P55v4j8AxiY1wYi0gl4BnhMVV/0WLQLWAPEiIjkc9x3gKdFpKGqHnDPux/Xh2GpAsRvfOucRz4cEJEY4FURqaKqJ3PawH3mNAOYqaoPZFm8Hnjfi3z4B/CAiHRW1e/c824HjgPfAF6dnZmAsM+MfJSIMyMRqYHrm8DrWQoRAKp6XlXPFHLfzdz7PpfPqvcAZ4Dpua2g+ffNtBNYCdznPnYYrsSa5WW4xs9EpAGugpDqnnLzV1xfjibltoIX+XAKV0Hy/OYbBbwNWD9fQco+M3JWIooR0AzXe/3ZR/sbICKnReR3YDuuU+dcP1TcLgZ2qur59Bkicqt7P+lTPy+OHQPc706qPwGVgU8K9zaMj3Rx//2SgAPAdcD0fL7gXAycVNX49BkicmmWfHjKi2PHAHeKSGURaQlcAbxbhPdivHdRlr/XaRE5jetybVb2mZGPknKZLv1U1lffFv8NPAVUBB4AapL/9d+cTqeXA+3cP28Fynh57NeAm3B9C35PVc/lf7Zu/Og7XPd8yuO6rt8N+Hs+2+T0B9vGH/nwH7y4dKyqa0VkO677CC2BOFX9zfIhIPYBN+Ywf0UO8+wzIx8lpRhtB9KA1rj+MEV1UlV3AIjIUOBb4GlcDSNysw24TkTKquo5APc35/T9eHVgVT0vIu8B0cCVQNvCvgnjM8np+QBsEpGLgTeAQXlssw2o4nkt350X6flwPo9ts5oF/A24APDmm7LxjfMef/cMIpKSw7r2mZGPEnGZTlWPAYuBESJSNetyESkjIhULuW/F1arlSXeLvdzMA8KBRwpznCxigGuBNaq61Qf7M741HtdlmU55rPNP4Cz5n0F5IxZoDpwGvvDB/owf2WdGzkrKmRHAcFwtjNa7m0tuwHUD8QrgUVyXWTa4160kIu2ybP97bn9EVV0mIttwJdiQXNb5XkSeBSaJyIW4bjzvAqoCf8b1xSCvG96e+9ohIrWA371Z3xRZ2RzyIU1Vf8ppZVXdKiKfAZNxXbLLaZ39IjICmOn+W8bg+sYbjutDoy7e58NJd+u8NFVN82Yb4yz7zMiuxBQjVd0rIh2AJ3B9c20MnAS2AC8CmzxW7wz8L8sutuG6Jp+bF4H3ROQlVd2WSwzjROQ74EFcz6FUw9UMdx2uljP/KMD7OebtuqbAwgDPSy0XkD0fzuK6R5SbKcDXInKjqv4npxVU9W0R+RkYDczFdR/hNPAj8BiuZ+C8klMrUeNzWfOiqOwzw4PYSK/GZCYiy4ADqnqf07GY4GF54V8l4p6RMd4QkVoichtwPXbvxbhZXgRGiblMZ4wX5uNqCDAN+MjhWEzwsLwIALtMZ4wxxnF2mc4YY4zjrBgZY4xxXFDdM6pVq5Y2adKkUNsmJCRQu3Zt3wZURMEYExQ+rvXr1x9RVcfekOWH/xUlJifzo0qVKlquXDmqVatG1arZnmvPU3H7O/hLYWI6ceIEiYmJHDly5KSq5v2HUdWgmTp27KiFVZRt/SUYY1ItfFzAOrX88JniFpOT+WG54X/+zg27TGeMMcZxVoyMMcY4rtgUo6ioKKdDyCYYY4LgjcufgvE9W0zBIRjfc0mMKaieM+rUqZOuW7fO6TBMLkRkvarm1RO1X1l+BDcn88NyI7h5kxtBdWZ04sQJoqKiiIuLczoU4yEuLi79W1HBmin5mOVHcAqW/DChzc6MjNfszMjkxc6MTG5C7szIGGNMyRSyxejJJ59k+vTpXq3bp08flixZ4t+ATNDwzI0VK1bQqFHug2k+8sgjvPXWWwGKzDjN288N+8xwQH4PIgVy8vahqsOHD2uDBg00KSnJq/W/++477dChg1frmtwRAg+9Zs2N5cuXa8OGDXNdPz4+Xhs1aqRnz571/hdhcuRkfhQmN/Jinxm+5U1uhOSZ0Zw5c+jRowcVKlTwav3LL7+ckydPYteUi7+C5kb9+vVp2bIlCxcu9HNkxmkFyQ37zAi8kCxGixcv5vrrr894ffz4cW699VZq165N9erVufXWW9m/f3+mbbp06cKiRYsCHaoJsKy5ke65556jVq1aNGnShHnz5mVaZrkR+rxpaZlTbnz66ae0a9eOKlWq0LRp00yX5iwviq4gLS1Dshht3LiRFi1aZLxOS0vjvvvuY8+ePezdu5cKFSowYsSITNu0atWKH3/8MdChmgDLmhsAhw4d4siRIxw4cID33nuPqKgotm3blrHcciP0Va1alZiYGHr27JnrOllzY+3atfTv358XX3yRxMREVq1ahWdHvJYXRdezZ09iYmIATuS3bkgWo8TERCpXrpzxumbNmvTt25fw8HAqV65MdHQ0K1euzLRN5cqVSUxMDHCkJtCy5ka6CRMmUK5cOa6//npuueUWPv7444xllhslQ9bcmD17NoMGDaJbt26EhYXRsGFDWrZsmbHc8iKwQrIYVa9enVOnTmW8TkpKYujQoURERFClShWuu+46EhMTSU1NzVjn1KlTVKtWzYFoTSBlzY30eRUrVsx4HRERQXx8fMZry42SIWtu7Nu3j6ZNm+a6vuVFYIVkMWrbti2//PJLxuupU6eybds2vvvuO06ePMmqVasAV0vBdFu2bCEyMjLgsZrAypob4LqneObMmYzXe/fupUGDBhmvLTdKhqy5ccEFF7Bz585c17e8CKyQLEY9evTIdBnu1KlTVKhQgWrVqnHs2DGeeeaZbNusXLmSm2++OZBhGgdkzY1048aN49y5c3z11Vd89tln/PWvf81YZrlRMmTNjfvvv593332X//znP6SlpXHgwAG2bt2asdyJvJg3bx5NmjQhLCwsx8Y2xVp+bb8DOTVr1kyHDBmiCxcuzLPNekJCgjZs2DDjeYEDBw7o9ddfrxUrVtTmzZvrW2+9pYCeP39eVVXXrl2r7dq1K1T7eKO6cOFCHTJkiALbNcjzI2tupD9nNHHiRK1Zs6ZecMEFOnfu3Iz14+PjtWHDhvacUREEQ35485xR1txQVf3kk0/00ksv1UqVKmnTpk11yZIlqurMZ0ZsbKyGh4crkDGFh4drbGxsQOPwB7x4zsjxAuQ5FWQkwSeffFJffvllr9bt06ePLlq0yOt9m5x5k1D+nLzNj4LkxiOPPKJvvPGGV+uavDmZH77ODSc+MyIiIjIVovQpIiIioHH4gze5YR2lGq9ZR6kmL9ZRatGEhYWR0+exiJCWluZARL5jHaUaY0yIqFOnTo7zy5Qpw969ewMcTeBZMTLGGIedOXPGdalKJNP8smXLIiK0b9+exYsXOxRdYFgxMsYYhz322GMkJCQQHR1NREQEIkJERATvvPMOP/30E40aNaJHjx48/fTTpKSkOB2uX4R8MSrRTSGNMUBojwL8xRdf8Oabb/Lwww8zYcIEdu/eTVpaGrt376Zfv35cfPHFrFmzhvvvv59JkybRrVs3Dh065HTYXinQKMD5tXAI5FSQ1nSqxbspZDAiRFrTGWc4mR+hmhvHjx/XRo0aaatWrbwa2mLOnDlaoUIFrVevni5fvtz/AfqIN7kR0mdG0dHRJCUlZZqXlJREdHS0QxEZY4z3Ro4cycGDB5k7d65XQ1sMGDCA7777jqpVq3LjjTfy3HPPhXxLu3QhXYxya2FSElqeGGNC24IFC5g7dy7R0dF06uR9i/hLL72U77//njvuuIPo6Gh69uzJ0aNH/RhpYIR0MWrcuHGO81WVP/3pTyxevLjYfGswxhQfCQkJREVF0b59+0JdyalcuTIffPABb7zxBl9++SUdOnTgu+++80OkgRPSxWjSpEmEh4dnmlehQgX++te/smnTJnr06EGbNm2YOXNmtst5xhjjBFVl2LBhnDhxgrlz51K2bNlC7UdEGD58ON988w1hYWFce+21vPLKKzk+OBsS8rupFMjJ277pPMXGxmpERISKiEZERGQ0Xjh79qzOnTtX27dvr4DWrFlTn3rqKT1w4ID3d92MqgZH32NayPww/hcM+RFKDRjef/99BXTKlCk+2+exY8e0V69eCmjfvn01MTHRZ/v2BYpz33TeSktL0xUrVuhtt92mIqJlypTRe++9V3/44QefH6u48yah/DmF0gdOSeRkfoRKbuzbt0+rVq2qV199taakpPh032lpaTplyhQtVaqUNmvWTP/3v//5dP9F4U1uhPRlOm+ICNdffz0LFizgl19+YdiwYXzyySd06NCBLl268Omnn2YahM8YE1xEpImIfC4ix0XkkIi8LiKlnY6roFSVwYMHc/78eebMmUOpUqV8un8R4dFHH2XFihUkJSVxxRVX8Pbbb7vOOkJAsS9Gnpo1a8arr77K/v37efHFF9m1axe9e/emRYsWvPbaa5w+fdrpEI0x2b0JHAbqA+2A64HhTgZUGDNnzmTp0qW89NJLNGvWzG/Hueaaa/jf//7Htddey5AhQxg4cGCmwSWDVYkqRumqVavGmDFj2LlzJx999BG1a9fmoYceolGjRjz66KPWNNyY4HIh8LGq/q6qh4AlQGuHYyqQnTt3MmbMGLp168awYcP8frw6deqwZMkSxo8fz/vvv0/nzp3ZsmWL349bFCWyGKUrXbo0d9xxB6tXr2b16tV0796dl19+mYsuuoj/+7//Y82aNU6HaIyBV4A7RSRcRBoCN+MqSCEhNTWVgQMHUrp0ad55551snaH6S6lSpRg3bhxLly7l8OHDXHbZZXz44YcBOXZhlOhi5OmKK67go48+4tdff+Xhhx9m6dKlXHnllVx55ZXMnz+/2HZOaEwIWInrTOgksB9YByzwXCEhIYFOnTplTDExMYGPMhcvv/wyX3/9Na+99hqNGjUK+PG7devG//73P9q1a8fdd9/N8OHD+f333/1+3JiYmIy/B1Ar3w3ya+EQyCmYWsScPHlSX331VW3atGnGaIsvvfRS0DWZDCSsNZ3Jgz/yA9cX5r1ANFAOqAl8CkzREMiNTZs2admyZbV3796alpbmaCznzp3TRx99VAHt0KGD7ty5M2DH9iY3fHJmJCLlRGS2iOwRkVMi8j8Rudlj+Y0islVEkkRkuYhE+OK4/lS5cmUefPBBtm3bxoIFC2jSpAljxoyhUaNGjBw5kp07dzodYsgojvlhAqYGcAHwuqqeVdWjwLtAD2fDyt/58+fp378/VatWZebMmQG7PJebMmXKMGXKFD799FN+/fVXOnTowKeffupoTJ58dZmuNLAPVyuXqsDfgY/dTTJrAZ+459XAdYr9kY+O63elSpXitttuY8WKFaxfv56//OUvzJgxg+bNm/OXv/yFVatWhUzTSQcV2/ww/qWqR4BdwN9EpLSIVAMGAD86GpgXJk2axA8//MBbb72V6yiuTujVqxc//PADzZo1o3fv3owZM4bz5887HZb/LtMBPwF9gSjgW4/5FYFkoGXWbYL1VDurAwcOaHR0tNasWTPjlPf999/Xs2fPOh2aX+HDyzDFOT9KKl/mh2bOlXbACuA4cASYD9TRIM6NdevWaalSpfSee+5xOpRc/f777zp8+HAF9Oqrr9Z9+/b57Vje5IZfGjCISF3gYmAzrhuPGd9iVPUMsJMQa5rpqUGDBkycOJG9e/dm9Ht37733cuGFF/Lcc88xc+ZMG/AvD8U9P4xvqeoGVe2iqtVVtZaq/lVVDzsdV25+//13+vfvT7169XjttdecDidX5cqV44033uDDDz/kxx9/pH379ixbtsy5gPKrVgWdgDLAl8BM9+vZwPNZ1vkGGJh128aNG2vHjh0zppkzZ/qnTPtYamqqfv7559qtW7dMA/1RDAb8mzlzZsbfA9itlh/Gg6/zo7BTMJ0ZjR49WgFdunSp06F4bevWrdqmTRsVER07dqzPuyoi0H3T4boH9Q/gc6CMe94rwJtZ1tsI9M26fTAlVGHVr18/x4LUuHFjp0MrMm8SKq/J8qN4K2p+FGUKltxYuXKliogOGzbM6VAK7MyZMzpgwAAF9KabbtLffvvNZ/v2Jjd8dplOXE1FZgN13R8k6XfENgORHutVBJq65xc7uY1Nv3fvXl599VVOnjwZ4IiCg+WHKe5Onz7NwIEDufDCC3nxxRedDqfAwsPDmTNnDrNnz+brr7+mXbt2fPXVVwE7vi/vGc0AWgE9VTXZY/6/gTYi0ldEygNjgZ9UdasPjx00chvwr2zZsowcOZKGDRvy0EMP8csvvwQ4MsdZfphibcyYMezevZv33nuPSpUqOR1OoQ0aNIg1a9ZQqVIlunbtyp133klERIT/74Hnd+rkzQRE4Loc9Ttw2mPq515+E7AVVyupFUCTnPYTLKfaRREbG6vh4eE53jNau3at3nvvvVqmTBkF9Oabb9bFixdramqq02F7hUJehrH8KBkKmx++mJzOjcWLFyugjz76qKNx+NKJEyf08ssv98k9cG9yw5HEyW1yOqF8JbcB/9IdPHhQn3nmGa1Xr54CevHFF+urr76qJ06ccChi7zj5YaPFKD+Kq5JajI4dO6YNGjTQ1q1ba3JysmNx+EPjxo1zvAceERFRoP14kxvWN50f9OvXj927d5OWlsbu3bvp169fpuX16tVj7Nix7Nmzhw8++IAaNWpk9Bo+cuRItm/f7lDkxoSmEydOEBUVRVxcXMCP/eCDD3L48GHmzp1L+fLlA358f9q3b1+O870d2SAuLo6oqChwPeyet/yqVSCnkvzN97vvvtN77rkn4xJejx49dMmSJUF1CY8QPDPK7yzV+I6T+eHUZ8c///lPBfSZZ55x5Pj+FhEREbAzI8cLkOdUkotRuoMHD+r48eO1bt26CmiLFi309ddf15MnTzodWsgVo7zu3xnfK2nF6NChQ1qrVi3t2LGjnjt3LuDHDwRf/R+yYhTCzp49q7GxsXrZZZcpoFWqVNGRI0fq9u3bHYsp1IpRbt/qisMzX8GoJBWjtLQ0ve2227RcuXK6efPmgB470HxxdcGKUTGxZs0avfvuu7VMmTIqInrLLbfo0qVLA94lfagVIxHJsRgB2q1bN33uuef022+/LbbfagOtJBWjOXPmKKBTp04N6HFDlTe5YQ0YQkDnzp2ZN28ee/bsYezYsXz//fd0796dSy65hDfffJPTp087HWJQyu2Zr8qVK3Po0CGeeuoprrrqKqpXr86f//xnXnjhBdauXWsDKZo87du3j4ceeohrr72WkSNHOh1OsRFUxcjJFjGhoH79+owfP569e/fy/vvvU6lSJR544AEaNWrEI4884rcxlgrUIsaPCpofkyZNIjw8PNO88PBwZsyYwU8//cThw4eZP38+AwcOZN++fTzxxBN07tyZGjVqcMstt/DSSy+xfv16UlNT/fF2io1gyY9ASEtLY9CgQaSmpjJnzhxKlSrldEjFR36nToGc7DJdwaSlpenq1av1rrvu0tKlS6uIaM+ePfWLL77wyyU8QuwynWrBrncfOnRI//GPf+iwYcO0RYsWGZf0qlatqj179tRp06bpDz/8EFQtHIOJk/kRqM+O119/XQHrpLeAvMkNca0XHDp16qTr1q1zOoyQFB8fz1tvvcVbb71FQkICrVq14qGHHuLee++lYsWKPjmGiKxX1U4+2VkhBDo/4uPjWblyJcuXL2f58uXs2LEDgOrVq3P99dfTpUsXunbtSps2bQgLC6qLDI5wMj8CkRvbt2+nXbt2XHfddXz++eeOj9waSrzKjfyqVSAnOzMquuTkZH3vvfe0Q4cOGd/qH3nkEd25c2eRW8UQgmdGvrRv3z59//33ddCgQXrRRRdlnDnVrFlT+/btq6+99ppu2rQp4A1LgoWT+eHv3EhJSdGrrrpKq1Wrpvv37/frsYojb3LD8QLkOTn9YVOcpKWl6TfffKP/93//p6VLl1ZAS5UqVaTnBUp6Mcpq9+7dOmfOHB0wYECmblNq166tf/3rX/XNN9/ULVu2ZCtOxfVB3OJcjJ5//nkFdN68eX49TnHlTW7YZboS4MCBA7Ru3ZoTJ05kWxYREcHu3bu92k9Ju0xXULt27WL58uWsWLGC5cuXs3//fsDV/VP6Jb1Tp04xduxYkpKSMrYLDw8nJiYmW7dRoaa4XqbbuHEjnTp1omfPnsyfP98uzxWCN7lhxaiECAsLI6e/tYiQlpbm1T6sGHlPVdm5c2dGYVq+fDkHDx7Mdf2CfCkIVk7mR/PmzbVr16707NmTnj17+my/586do3PnzsTHx7Np0yZq167ts32XBHFxccTFxTFr1qwdqto8z5XzO3UK5BRsl2GKE1/0MYVdpiu0tLQ03bZtW64P4YqI0yEWmZP54a/cePrppxXQTz/91C/7Lym8yQ1rAlRC5PbMzaRJkxyKqGQRES6++GIiIiJyXF6uXDnWrFkT4KhMXtauXcvkyZMZMGAAvXr1cjqcYs+KUQnRr18/YmJiiIiIQESIiIgoFvcpQk1OXwrKlClD6dKlufLKK7nlllsIlUuRxVlycjL9+/enQYMGvPLKK06HUyJYMSpB8htnyfhfTl8K3n33XQ4ePMjkyZNZs2YNl112GbfddhsbNmxwOtwS66mnnmLbtm288847VK1a7DuWCApWjIwJsJy+FFSqVIknnniCXbt2MWHCBFatWkX79u25/fbb2bRpk9MhlygrVqxg+vTpPPDAA9x0001Oh1NiBFUxsr7pglOw9D1WEvKjSpUqPP300+zatYuxY8eybNky2rZty1133cXWrVudDi9HwZIfvnDq1CkGDhxIs2bNeOGFF5wOp2TJr4VDIKdQbi1VEmCt6QLu6NGj+tRTT2nFihU1LCxM7733XkfHtMqLk/nhq9wYPHiwhoWF6bfffuuT/RkXb3IjqM6MjDGZ1ahRg0mTJrFr1y5Gjx7NP//5T1q2bMmgQYPYtWuX0+EVK4sWLeLtt9/mscce48orr3Q6nBLHipExIaB27dpMmTKFX3/9lQcffJAPPviAiy++mKFDh7J3716nwwt5R48eZfDgwVx66aWMHz/e6XBKJCtGxoSQevXq8fLLL7Nz506GDh3KnDlzaN68OSNGjODAgQNOhxeyRowYwdGjR5k7dy7lypVzOpwSyYqRMSGoYcOGvP7662zfvp377ruPmTNn0rRpU0aNGsWhQ4ecDi8kzJs3jyZNmhAWFsY//vEPbrvtNtq1a+d0WCWWFSNjQljjxo156623+OWXX+jXrx+vv/46F110EY8++igJCQlOhxcwBW1pOW/ePKKiotizZw+u++vw+eefM2/ePH+GWeIUpKWldZRqvGYdpQa/HTt2MGHCBGJjY6lQoQIPPfQQo0ePpmbNmn4/dij12t2kSRP27NmTbX5x6LA2GHmTG8XmzGj+/Pm0b9+edu3a0bJlS+6+++4i7S8xMZEpU6Zkmjd48GC++uqrIu23S5cufPbZZzkuW7ZsGZ06daJcuXKMGTOmSMcxJVOzZs1477332Lx5M7169eL555/nwgsvZOzYsSQmJjodXtDIrdGHNQZxTrEoRgcPHmT48OEsXLiQDRs2sGXLFh577LEi7TOnYvT2229z7bXXFmm/ebnooouYNWsWjz76qN+OYUqGli1b8sEHH/DTTz/RvXt3JkyYQJMmTZgwYQInT550OjzHNW7cuEDzjf8Vi2J06NAhypQpk3EpQkQy3Yj87rvv6Nq1Kx07dqRjx44sWrQIgN27d1OrVi2io6Np3749LVq04OuvvwbggQceIDExkXbt2nHVVVcBmc9qBg4cyLBhw7jhhhto3rw5/fv3z7j2fPLkSQYPHszll19O27ZtGTlyJKmpqfm+j2bNmtG+fXtKly6d53opKSl0796dTp060bp1a+677z7OnTsHwLfffkuHDh1o164drVu35sMPPwRc19QHDRrEpZdeSmRkJCNGjPD212tCWJs2bZg/fz4bNmygS5cujB07lgsvvJDnn3+e06dPOx2eY6wX+yCU31OxgZwK+xR1amqq3nbbbVqzZk3t27evvvzyy3rkyBFVVT1+/Li2a9dO4+PjVVU1Pj5eGzZsqMePH9ddu3YpoHFxcarqGg76qquuUlXVXbt2ac2aNTMd5/rrr89Yd8CAAXr11VdrcnKynj17Vi+55BJdtmyZqqref//9Onfu3IzY7rzzTo2Jicm2j9yMGzdOR48enevytLS0jPeXlpam9957r86YMUNVVXv16pVx7LS0ND1+/Liqqg4cOFBHjBihqampqqqakJCQZww5wXpgCHnr1q3TW265RQGtVauWvvjii3rmzBmfDIXuZH4UJjeK6/Dvwcib3Mj7K3iApbeIKehojWFhYSxYsIBNmzaxcuVKFixYwIsvvsjGjRtZs2YNu3bt4uabb85YX0TYsWMHtWrVolKlStx6660AXHHFFYwePdrr4/bu3Zvy5csD0KFDB3bu3Em3bt1YuHAha9euZerUqQAkJSXRqFEjr/ebn7S0NF566SUWL15Mamoqx48fz/iW17VrVyZPnsyePXvo1q0bnTt3BuCzzz5j/fr1hIW5ToZr1arl9fHSR2skSPqm8/VoniVJx44d+eyzz1izZg3jxo3j0UcfZcKECSQnJ3P+/HkA9uzZk94Cyque3YMlPwqqX79+1nN9EAmqYlS1alViYmIKvX2bNm1o06YNDzzwAJdccgkrVqygXLlytG3bllWrVmVbf/fu3ZkecCtVqhQpKSleHy+9EGXdVlVZsGABF110Ua7bHj16lBtvvBGAFi1a8NFHH3l93A8++ICvv/6ar776isqVK/Pcc8/xyy+/ADBq1Ch69uzJl19+yYMPPsif/vQnJk6c6PW+c5L+4T9r1qwTRdpRERU1P8wfrrjiCpYuXcrXX3/NTTfdlFGI0iUlJREdHe3Vh3Ww5IcJbcXintGBAwdYvXp1xuv9+/eTkJDAhRdeyFVXXcX27dtZvnx5xvLvv/8+4/5ObqpUqUJSUlKBilO69FZM6feJjhw5kq0fsZo1a7JhwwY2bNhQoEIErsYVtWrVonLlypw4cYIPPvggY9kvv/xC06ZNGTp0KCNHjmTt2rUA3Hrrrbz44osZ7/vIkSMFfl+m+Lnmmmsy7jdmFUwty0TkThHZIiJnRGSniPivJZFxRLEoRikpKYwbN44WLVrQrl07evTowcSJE2nfvj3Vq1dn4cKFPPPMM0RGRtKqVSvGjx+fbzGqUaMG/fr149JLL81owOCt6dOnU6pUKSIjI7n00kv585//7FVXLV9//TWNGjVi2rRpzJw5k0aNGrF06dJs6/Xv359Tp07RunVr/vrXv2Zq4ffqq6/SunVr2rdvz2uvvZZxQ/bll1/m1KlTtGnThsjISJ599lkAFi5cyODBgwv0/kzxEuwty0SkG/ACcB9QGbgO+NXRoIzP2UOvxmv20GvxlN4bQVJSUsa88PDwAg9L76/8EJFvgdmqOju3dSw3gluJeujVGFM4OQ2FXtBC5C8iUgroBNQWkR0isl9EXheRCk7HZnwrqBowGGOcEcQty+oCZYDbgWuB88CnwNNAdPpKCQkJdOr0xxfvqKiojBaBxhkxMTGeDY7ybb5rxcgYE8yS3f++pqoHAURkGlmKUe3atbHLdMHF8wuBiOTbYspnl+lEZISIrBORsyIyJ8uycBF5U0SOiMgJEcneztoUW5YbprBU9TiwHwiem9vGL3x5ZhQPTAS6A1mv58a4j9UKOAa08+FxTfCz3DBF8S7woIgswXWZbhSQc2/DJmT5rBip6icAItIJyOhuQERaAL2ARqqa3kPjel8d1wQ/yw1TRBNw3XP4Bfgd+BiwTuSKmUC0pusM7AGecV+K2SgifQNwXBP8LDdMvlT1vKoOV9VqqlpPVR9S1d+djsv4ViCKUSOgDXACaACMAN4TkVZZV0xvEZM+WdcvzouJicn4e+BFi5gC8jo3wPIjGPk5P0wJEojWdMm4rvNOVNUUYKWILAf+BGzxXNFaxASfgraIKSCvcwMsP4KRn/PDlCCBODP6KQDHMKHJcsMYA/i2aXdpESkPlAJKiUh5ESkNrAL2Ak+617ka6AJk73StCILxkk0wxgSBj8vp3IDg/FtYTMEhGN9ziYwpvwGPvJ2A8bieBfCcxruXtQZWA2eAn4G/5LSPogyeFowDrwVjTKqFj4tCDp7mi9xQy4+AKEpMhc0PX0yWG/7n79zwZdPu8e4PnZyWbQau9NWxTGix3DDG5Ceoeu0WkQRcTX0LoxYQbDdQgzEmKHxcEapa29fBeMvyIyCKEpNj+SEiJ4BzQCKu1pkFUdz+Dv5SmJiqAtWAsqqa50jAQVWMjDHGlEw2hIQxxhjHWTEyxhjjOCtGxhhjHGfFyBhjjOOsGBljjHGcFSNjjDGOs2JkjDHGcVaMjDHGOM6KkTHGGMdZMTLGGOM4K0bGGGMcZ8XIGGOM46wYGWOMcZzPxjPyhSpVqmi5cuWoVq0aVavm2dt4NgkJCdSu7djoBjkKxpig4HGdOHGCxMREjhw5cjK/buD9yfLD/woTUzDkh+WG//k9N/IbfS+Qk43WGBiBHunVV5Plh//ZSK/BobjF5E1u2GU6Y4wxjrNiZIwxxnHFphhFRUU5HUI2wRgTBG9c/hSM79liCg7B+J5LYkxBNex4p06ddN26dU6HYXIhIutVtZNTx7f8CG5O5oflRnDzJjeKzZmRMcaY0GXFyBhjjOOsGBljjHGcFSNjjDGOC9li9OSTTzJ9+vRcl4sIO3bsAOCRRx7hrbfeClBkxmn55Yanyy+/nM2bN/s3IBM0subGjBkzqFu3LpUqVeK3336jZcuWHD582LkAS7L8nooN5OTtE76HDx/WBg0aaFJSUl5P/Or27dtVVTU+Pl4bNWqkZ8+e9Wr/JmeEQA8M3uSGp48++kj79Onj1bomb07mR2Fy49y5c1q+fHndsGFDxjovvPCCPvLII0X4LZiceJMbQXVmdOLECaKiooiLi8tzvTlz5tCjRw8qVKjg1X7r169Py5YtWbhwoS/CLHHi4uLSnzFwrF868C4/CpobvXr1Yvny5Rw8eNBXYZY4wZAfhcmN3377jd9//53WrVtnrHP33Xfz3nvvcfbsWb/HXBIUKDfyq1aBnLw9M+ratau+//77meZNmTJF69Wrp/Xr19fZs2dnOjNSVZ04caIOHDjQq/2bnBECZ0Y55UZcXJxGRkZq1apV9corr9Qff/wx0/KbbrpJ58yZU7BfhsnGyfwoaG5s27ZNw8PDFdCKFStq165dM9Zr1qyZrlixogi/CZOVN7kRVGdG3tq4cSMtWrTIeL1kyRJeeuklvvjiC7Zv386XX36ZbZtWrVrx448/BjJM44CsufHDDz8waNAgZs6cydGjRxk6dCi9evXK9M3XcqNk8MyNiy++OONeYWJiIv/9738z1rN8cEZIFqPExEQqV66c8frjjz/mvvvuo02bNlSsWJHx48dn26Zy5cokJiYGLkjjiKy5MWvWLIYOHUrnzp0pVaoUAwYMoFy5cqxZsyZjHcuNkiFrbuTG8sEZIVmMqlevzqlTpzJex8fHc8EFF2S8joiIyLbNqVOnqFatWiDCMw7Kmht79uxh6tSpVKtWLWPat28f8fHxGetYbpQMWXMjN5YPzgjJYtS2bVt++eWXjNf169dn3759Ga/37t2bbZstW7YQGRkZkPiMc7LmxgUXXEB0dDSJiYkZU1JSEnfddVfGOpYbJUPW3MiN5YMz8ixGInKRl1OTAMULQI8ePVi5cmXG6zvuuIM5c+bw888/k5SUxDPPPJNtm5UrV3LzzTcHMkzjgKy5MWTIEN566y2+++47VJUzZ86waNGijG/IZ8+eZf369XTr1s2pkE2AZM2NnBw4cIBjx45xxRVXBCgqky6/Ycd3AApIPuslAxV9EpEX+vfvT7t27UhOTqZChQrcfPPNjBo1ihtuuIGwsDAmTpzIvHnzMtY/ePAgP//8M7179w5UiMYhWXOjU6dOzJo1ixEjRrB9+3YqVKjANddcw3XXXQfAwoUL6dKlCw0aNHA4cuNvWXMjJx988EHGfUUTWHkOISEip1Q13zt+InJcVasXNZiCdAP/1FNPUadOHUaNGpXvuqNHj6Zp06YMHz68iBGWbKEyhERBcqNz587Mnj2bNm3a+CDCki0UhpDIKzfOnj1LZGQkq1atok6dOn6IsuTyJjfyK0aDVPUdLw40UFXneLHencA4oDFwCBioql+lL7cxSYKbPz9s8ssNsPwIdk7mh+VGcPMmN/K8TOdNIXKvN8eLYLoBLwD/B6wF6nuzb1P8WW6YvFh+lAz53TPKxN1QoS1QyXO+qn7gxebPAM+qavoDHgcKcmxTrFlumLxYfpQAXjftFpEngS3AWOBvHtMwL7YtBXQCaovIDhHZLyKvi0imu4gJCQl06tQpY4qJiSnIezF+EBMTk/H3AGr5ev/e5gZYfgSjYMkPy43gU9DcyPOeUaYVRY4A16nqzwUNSkQa4Po2sx7oCZwHPgVWqGp0+np23Te4+eOegLe5AZYfwc7J/LDcCG7e5EZBHno9CuwuZCzJ7n9fU9WDqnoEmAb0KOT+TPFhuWHyYvlRQhTkntEoIEZEpgOZRp9S1exdHmReflxE9uN6ZsmYDJYbJi+WHyVHQc6MygJ/wtWaZbfHtMvL7d8FHhSROiJSHVdx+6wAxzfFl+WGyYvlRwlQkDOjN4GngH/wx6lzQUzAdRPrF+B34GNgUiH2Y4ofyw2TF8uPEqAgxag08K6qphbmQKp6HhjunozJYLlh8mL5UTIU5DLdS8ATIpJfP3XGGGNMgRTkzOghoB7wlIgc9Vygqo19GpUxxpgSpSDF6B6/RWGMMaZE87oYqWreA4EYY4wxhVSQ7oA+EZFrs8y7VkT+6fuwjDHGlCQFacBwPfBtlnmrga6+CubEiRNERUURFxeX77rz58+nffv2tGvXjpYtW3L33XcX6diJiYlMmTIl07zBgwfz1Vdf5bKFd7p06cJnn+X8SMSyZcvo1KkT5cqVY8yYMV7vc/fu3dSq5fNuwHIVFxdHVFQUQNWAHTQHBckPEzjBkB+WG8GpQLmhql5NuPqHqpJlXjXgkLf7yG/q2LGjeiM+Pl5r1aqle/fuVVXVtLQ0/d///ufVtrnZtWuX1qxZs0j7yMn111+vcXFxOS7bvn27/vDDDxodHa2jR4/2ep/+ijU/wDr10d+6MJO3+WGc4WR+WG4EN29yoyBnRkuBmSJSBcD97+vAkgLswycOHTpEmTJlqFmzJu5YaNeuXcby7777jq5du9KxY0c6duzIokWLgD/OKKKjo2nfvj0tWrTg66+/BuCBBx4gMTGRdu3acdVVVwGZz2oGDhzIsGHDuOGGG2jevDn9+/dPL8icPHmSwYMHc/nll9O2bVtGjhxJamr+j2M1a9aM9u3bU7p0/rfu3njjDZo1a8a1117L7NmzMy37/PPPufrqq+nYsSNXXnkla9a4etrftm0bV155JZGRkbRp04aXXnoJgHPnzjFmzBjatGlDZGQkf/nLX/I9vjHFzbx582jSpAlhYWE0adKEefPmOR1SyZZftUqfgOrAIiAFV990KUAcUM3bfeQ3efvtJjU1VW+77TatWbOm9u3bV19++WU9cuSIqqoeP35c27Vrp/Hx8arqOotq2LChHj9+XHft2qVAxplKbGysXnXVVaqa89mG51nNgAED9Oqrr9bk5GQ9e/asXnLJJbps2TJVVb3//vt17ty5GbHdeeedGhMTk20fuRk3blyeZ0Y//vij1q9fXw8dOqSqqn/7298yYt2xY4deccUVeuLECVVV3bRpk15wwQWqqvrQQw/ps88+m7GfY8eOqarq+PHj9S9/+YuePXtWVVUTEhLyjC8ddmZk8uBkfhQ0N2JjYzU8PFxx9XmngIaHh2tsbGyh37/JnTe5UZDWdMeBW0SkHnABsE9VD/muLHovLCyMBQsWsGnTJlauXMmCBQt48cUX2bhxI2vWrGHXrl3cfPPNGeuLCDt27KBWrVpUqlSJW2+9FYArrriC0aNHe33c3r17U758eQA6dOjAzp076datGwsXLmTt2rVMnToVgKSkJBo1auSz97tixQpuueUW6tatC0BUVBQff/wxAEuXLmXnzp1cd911GeunpKTw22+/cd111zFmzBjOnTtH165d6drVdXvvs88+Y+rUqZQtWxYgoPefjAkG0dHRJCUlZZqXlJREdHQ0/fr1cyiqkq1AI70CuAuQI0UoqzZt2tCmTRseeOABLrnkElasWEG5cuVo27Ytq1atyrb+7t27KVeuXMbrUqVKkZKS4vXx0gtR1m1VlQULFnDRRRfluu3Ro0e58cYbAWjRogUfffSR18dVzb3DYlXlz3/+M3Pnzs22rG/fvlx55ZUsW7aM559/nnfeeYfY2Ng892dMSbB3b84DDezZs4effvqJSy+9FOtsJrDyvGckIr96sxMR2e6bcLxz4MABVq9enfF6//79JCQkcOGFF3LVVVexfft2li9fnrH8+++/z/cDuEqVKiQlJRWoOKXr1asXzz//fMZ9oiNHjrBrV+bOzGvWrMmGDRvYsGFDgQoRQNeuXfn88885fNg1cofnPaM//elPLFmyhM2bN2fM+/777wHYsWMH9erVY+DAgYwbN461a9cC0LNnT6ZPn865c+cy4jWmJGncOPdOYyIjI2nYsCEDBw7kww8/tP8fAZLfmVFDEXnWi/3U9UUw3kpJSWHcuHHs2bOHChUqkJaWxsSJE2nfvj0ACxcu5NFHH2XUqFGcO3eOiy66KN8mnzVq1KBfv35ceumlVK9enW+/zdqKPXfTp0/nscceIzIyEhGhXLlyTJ8+nQsvvDDP7b7++mvuvPNOTp48iaryj3/8g9mzZ9O9e/dM67Vt25annnqKq6++mnr16nHLLbdkLGvevDmxsbHcf//9JCcnc+7cOa6++mouu+wyPv74Y+bNm0fZsmUREV555RUAnnjiCZ588knatWtH2bJladasGf/85z9Zt24dY8eO5fPPP/f6vRsTiiZNmkRUVFSmS3Xh4eG88MILVKxYkaVLlxIXF8d7772HiNCxY0e6d+9O9+7dueKKKyhTpoyD0RdPeQ47LiLvermf86oaVdRgbOjg4OaPYaULwvIjuDmZH4XJjXnz5hEdHc3evXtp3LgxkyZNynS/KDU1lfXr17N06VKWLl3KmjVrSE1NpUqVKtxwww0ZxSm/L53Gu9zIsxgFmn3YBDcrRiYvoVaMCioxMZH//ve/LFmyhKVLl2bcd2revHlGYerSpQuVKlXyaxyhyJvcKHADBmOMKYmqVatGnz596NOnD6rKtm3bMs6aZs+ezeuvv06ZMmW45pprMopT+qV7k7+CPPRqjDEG1+MiLVu2ZOTIkXz++eccO3aML774gpEjR3LkyBGeeOIJ2rdvT/369enfvz/z5s0jISHB6bCDmp0ZGWNMEZUvX56bbrqJm266iRdffJH4+HiWLVvG0qVL+fzzz3n//fcB1/OJ6WdNV111lTWE8BBUZ0bW2WFwCoaOMMHyI1gFQ34EW240aNAgo2n4b7/9xtq1a5kwYQLh4eFMmTKFLl26UKNGDXr37s2MGTP49VfXUzTB2EVRUWLyaUepwKtZXt+f5fW/8tuHt5N19xLcsO6ATB6czI9Qyo3ExET95JNPdOjQodqkSZOM7ojq1KmjpUuXztRFUYUKFXTmzJl68uRJTUpK0nPnzmlaWlrAYvVVt0ne5Ea+relE5KSqVvF4fUxVa+S2vCistVRws9Z0Ji/FvTWdP6gqv/zyC0uXLuXxxx/n999/92q7UqVKUbp06TynMmXK5LtOfuvGxsZy6tSpbMePiIhg9+7dXr9PX7Wmy9oUxJqGGGOMD4gILVq0oEWLFowaNSrX9V588UVSUlLynM6fP5/vOunTmTNnvFovp0IEuXenVBTeFKOsp07B82CSMcYUE40bN2bPnj3Z5kdERBRo8E1fatKkSY4x5dWdUmF504ChtIh0FZEbROSGHF6X8nlUxhhTwkyaNInw8PBM88LDw5k0aZJDEQU2Jm/OjA4D73i8Pprl9WGfRmSMMSVQeldEeXVRVJxjyrcYqWoTnx/VGGNMNv369Qu68ZQCFVOhnjMSkRYi8hcRifB1QMYYY0qefIuRiEwVkXs8XvcHNgMxwFYRuTnXjY0xxhgveHNm1BvwHDb1OeAhVa0NDAPG+SEuY4wxJYg3xai2qu4FEJE2QE0gfajRWOBiP8VmjDGmhPCmGJ0QkfSRXK/F1a3DWffrMvjwIdhg61/KuARD32Ng+RGsgiE/LDeCU0Fyw5vugKYCHYB/A6OB51V1hnvZ9cBUX3UBEqpdepQU1h2QyYt1B2Ry401ueHNm9ASwAuiGq9HCTI9l7dzzjDHGmELLtxip6nlVfUZVe6rqJFVN81j2iqo6WoyCscv1YGW/K2NMsMr3oVd3U+48qepc34RTMPPmzSMqKoqkpCQA9uzZk359MugeHHOa/a6MMcHMm+6A5gA7gEPk3FhBAUeKUXR0dMaHa7qkpCSGDRvG+vXrqVy5csZUpUqVTK89p/DwcJ+OUz9v3jy/dp+hqvz++++cOnUq1+n06dOZXs+dOzfH31V0dHSxLkb+/lsYY3zDm2L0KnA7cApX0Vng0ZrOUbl1Y3769GlmzZrF6dOnvdpPWFhYroUqryKW07LPPvuMYcOGZTsDSUlJ4dZbb82zgORXUDyn1NRUr95buXLlqFy5crZClN/vsDiws0FjQoc3fdONEpHRwJ+B/sB0EfkMeE9Vv/Z3gHnJq8v13bt3k5aWxpkzZzh16hQnT57M9YM9t2W//fZbptfnz58vVJxJSUkMHDjQq3XLlCmTrcBVrVqVRo0aUblyZSpVqpRrccw6VapUiTJlygCB7Qo+WOR25vzII49w0003UadOHZ+eERtjCs+bMyNUNRVYBCwSkSrA08AKEemmqsu92YeIrACuAFLcsw6oaouCh/yHSZMmZfrmC5m7N/c842nQoEFRDgXA2bNn8y1gDz/8cK7bv/LKKzkWDM/X5cqVK3KcOcnvd+Ukf+QG5H7Wd/jwYerVq0fNmjVp3bo1rVu3pk2bNhk/16pVq6iHDjnBfDnTX/lhgkx+45KnT7geWhoKfAtsx9UNUPUCbL8CGJzXOoUZxz42NlYjIiJURDQiIqLAY7P7WkRERKbx4tOniIgIR+NSLfrvCi/GsS/M5E1uaCHyI7e/RZ06dfSVV17RIUOG6FVXXaVVq1bNtvyGG27QBx98UN966y396quv9NixYwU6diiJjY3V8PDwTL+D8PDwkMqPwnx2mMDxJje8SYRbgflAPPAWcHV+2+SynxKRUL76jx2MQq0Yefu3SEtL0/379+uSJUt06tSpOmjQIL388su1UqVKmbZt0KCBduvWTUeNGqVvv/22rl69Wk+cOFHA32JwSE1N1SNHjui2bdu0bt26PvkCZcXI5Mab3PCmB4Y0YBvwGZCcy9nV2Dx3QsapdmtcLfK2AdGqusJzneLyFHUwX/IoCn89Ye9NbkDh8qMof4u0tDT27dvH5s2b2bRpE5s3b2bz5s38/PPPJCf/8V/hggsuyHSZr3Xr1lxyySVUrFjR5zHlJDk5maNHjxZoOn78OF783yctLS3PdbKs71h+FJfPjuLKm9zwphjNwfVNKVeqep8XwXQGfgbOAXcCrwPtVHVn+joRERFau3btjG2ioqIyWj8ZZ8TExBAT43quef369XvUD4MtepMbEDz5kZqayu7duzOKU3qh2rp1K2fP/tHQ9MILL8x2T2rDhg2MGDEi2727mJgY7rrrLhITEwtcWDwLY1YVK1akZs2aGVONGjUyva5ZsyajR48mISEh27bpDYHyEiz5ESy5Yf5Q0NzItxj5i4gsARap6mvp8+zbTXALVN9jOeUGBH9+pKSk8Ouvv2Y6i9q0aRO//PJLvi0xw8JcnaHkdiYSFhaWYyHJa6pRowbly5fPN+6sTeDhjwJZkDM2J/Mj2HOjpPMmN7xqTZfHAdoCf1fVvxZic8WHPX6bYiUkc6N06dJcfPHFXHzxxfTp0ydj/vnz59m+fTubN2/mjjvuyHHbtLQ0nn766VwLS9WqVTMKlq+lF5wQurQckvlh8uZNd0DhwJO4OkXdDowHagFTcXWe+p4X+6gGdAZW4mqe+X/AdcCowgRtio+SkBtlypThkksu4ZJLLiEiIiLXZ+MmTJjgQHQu/fr1C8riUxLyw7h481XrDaAnrmu2NwH/wpUYm4EmqvqAF/soA0wEEoAjwINAb1XdVpigTbFSonJj0qRJhIeHZ5oXLM97BakSlR8lmTeX6brjull4WEReA/YC16vqV94eRFUTgMsKGaMpxkpaboTgJTFHlbT8KMm8KUaVVPUwgKruF5HTBSlExpjMgvWSmDFO8qYYlRaRrnjcMMz6WlX/64fYjDHGlBDeFKPDwDser49mea3ARb4MyhhjTMniTa/dTQIQhzHGmBLMPw8uGGOMMQVgxcgYY4zjgqoYnThxgqioKOLi4pwOxXiIi4tL7+erqpNxWH4Ep2DID8uN4FSQ3HCsb7qcWP9SwS1QfY/lxvIjuDmZH5Ybwc2b3AiqMyNjjDElkxUjY4wxjrNiZIwxxnFWjIwxxjjOipExxhjHWTEyxhjjOCtGxhhjHFdsilFMTIzTIWQTjDFB8MblT8H4ni2m4BCM77kkxmTFyI+CMSYI3rj8KRjfs8UUHILxPZfEmIpNMTLGGBO6gqo7IBE5AZwDEoETBdy8FnDE1zEVUTDGBAWPqypQDSirqo71P2b5ERCFicnx/LDcCAi/5kZQFSNjjDElk12mM8YY4zgrRsYYYxxnxcgYY4zjrBgZY4xxnBUjY4wxjrNiZIwxxnFWjIwxxjjOipExxhjHWTEyxhjjOCtGxhhjHGfFyBhjjOOsGBljjHGcFSNjjDGOK+10AJ5q1aqlTZo0KdS2CQkJ1K5d27cBFVEwxgSFj2v9+vVHVNWxN2T54X9FiSm3/BCREcBA4FLgQ1Ud6LHsRuANoDHwHTBQVfe4lwnwPDDYvfps4HHNYagByw3/80duZKKqQTN17NhRC6so2/pLMMakWvi4gHVq+eEzxS2m3PID6AP0BmYAczzm18I19tBfgfLAi8Aaj+VDgW1AI6Ah8DMwLKdjWG74nz9yw3Oyy3TGGL9S1U9UdQFwNMuiPsBmVZ2vqr8D44FIEWnpXj4AmKqq+1X1ADAV1xmWKYasGBljnNIa+DH9haqeAXa652db7v65NaZYCqp7RkURFRXldAjZBGNMELxx+VMwvmeLiUpAQpZ5J4DKHstPZFlWSUTEfeknQ0JCAp06dcp4HRUV5fV7sb+DdwoaU0xMDDExMekva+W3flANO96pUyddt26d02GYXIjIelXtlP+a/mH5Edzyyw8RmQg0UncDBhF5BSijqsM91tkIjFfVf4nICaCbqq51L+sIrFDVyln3bbkR3Lz57LDLdMYYp2wGItNfiEhFoKl7frbl7p83Y4olK0bGGL8SkdIiUh4oBZQSkfIiUhr4N9BGRPq6l48FflLVre5N5wKPiEhDEWkAjAbmOPAWTABYMTLG+NvTQDLwBHCP++enVTUB6AtMAo4DnYE7PbabCcQBG4FNwCL3PFMMBVUxOnHiBFFRUcTFxTkdSrE0b948mjRpQlhYGE2aNGHevHlebRcXF5d+87KqXwPMh+VHcMovP1R1vKpKlmm8e9mXqtpSVSuoahdV3e2xnarqY6pawz09lrXhgilG8nsQKZBTQR6qeuKJJ/Tll1/2at2HH35YZ8yY4fW+i6PY2FgNDw9XIGMKDw/X2NhYr/dBiDz06pkby5cv14YNG+a6ruWG7ziZH4XJjVWrVunFF1+cseyyyy7TTZs2Fe7Nmzx5kxuOFyDPyduEOnz4sDZo0ECTkpK8Wj8+Pl4bNWqkZ8+e9Wr9UJaWlqanTp3S3bt36w8//KBffPGFfvTRR1qjRo1MhSh9ioiI8HrfoVCMsuZGfsWoJOWGvwV7Mcrvc+Ojjz7SPn36FPyNm3x5kxsh+ZzRnDlz6NGjBxUqVPBq/fr169OyZUsWLlzI7bff7ufoXJfDoqOj2bt3L40bN2bSpEn069evQPtQVc6cOcOxY8c4duwYR48e9frn8+fPe32cvXv3FvTtBbVgzw3jnPxyo1evXgwbNoyDBw9Sv379AEdnQrIYLV68mEGDBmWaN2XKFF5++WVEhGeffZYhQ4awfft2mjVrBkCXLl1YtGiR3z9w5s2bR1RUFElJSQDs2bOHIUOGcPToUa6//vp8C4nnvHPnzuV6nPDwcGrUqEGNGjWoWbMmrVq1yvjZc376v927d+fAgQPZ9tO4cWO//S6ckFNuAEydOpUXXniBUqVK8dxzz3HfffdlLAtUbhhnZc2NFStWcM8997B//34AypcvT8eOHVm2bBkDBgxwKsyg44sv194IyWK0ceNGWrRokfF6yZIlTJs2jf/85z9ceOGFDB06NNs2rVq14l//+pffY4uOjs4oROmSk5MZOXJkjutXqFAhU/Fo2bJlnkWlRo0aVK9e3etv/uleeOGFTEUSXAVt0qRJBX+TQSxrbgAcOnSIEydOcODAAb744gtuv/12evfuTfXq1YHA5YZxVk65kVWrVq348ccf81ynJMnpy3V6Twy+LkghWYwSExOpXPmPh7A//vhj7rvvPlq3dnVbNW7cOGJjYzNtU7lyZRITE/0eW16Xvf71r39lKzAFLSqFlZ44gfiG46SsuQFQpkwZxo4dS+nSpenRoweVKlVi27ZtXHHFFUDgciNdoL5pmsxyyo2sKleuzMGDBwMUUfDL6ct1UlIS0dHRVowAqlevzqlTpzJex8fHZ+qX6oILLsi2zalTp6hWrZrfY7vgggtyLEgRERH06dPH78fPS79+/Yr9h17W3ACoWbMmpUv/kerh4eGcPn0643WgcgMC+03TZJZTbmQVyFwINufPn2fbtm389NNP/Pjjj/z444/s2bMnx3X9ca85JItR27Zt+eWXX7jssssA103o9Ou+APv27cu2zZYtW4iMjMw239ciIyOz/aGK4+WwYJU1N7wRqNyA3L9pPvjgg1SoUIGGDRvSsGFD6tatS5kyZQISU0nhTW5s2bKFe+65J4BROePo0aMZBSe9+GzevDnjPnWZMmW45JJLqFixImfOnMm2vT/uNYdkMerRowcrV67M+CZ5xx13MGjQIO69914iIiJ49tlns22zcuVKBg8enG2+Ly1ZsoS4uDi6du3Kr7/+apdhHJA1N7wRiNxIl9s3yuPHj9O3b9+M1yJC3bp1adCgQUaBSv/Z898aNWrgGhDV5Ce/3Dh79izr16/nvffeC3Bk/pOSksL27dszCk968fFszFS3bl0iIyN56KGHiIyMJDIykpYtW1KmTJlsZ/Lgvy/XIVmM+vfvT7t27UhOTqZChQrcfPPNPPTQQ3Tt2pWwsDD+/ve/8/7771OuXDkADh48yM8//0zv3r39FtOBAwe49957adu2LYsWLQrYvSCTWdbcyE8gcsNTvXr1crwn0ahRIxYuXMiBAweIj4/P9O+ePXtYvXo1R44cybZd+fLladCgQY6FyrOI5fe7KAn3sfLLjYULF9KlSxcaNGjgQHRFd/z48RzPdn7//XcASpcuTatWrejatWtG0Wnbti1169bNdZ+BvNccskNIPPXUU9SpU4dRo0ZlW7ZlyxbatGnD2bNnKV26NKNHj6Zp06YMHz48+458ICUlhRtuuIEffviB9evX59tiJ1SFyhASeeVGVv7ODU+pqak0b96cXbt2ZZofHh5OTExMvv/Bz549y8GDBzMVqpyKV9bLgOC6X5LbWdZPP/3E5MmTSU5OLnBMnpzMD1/kRufOnZk9ezZt2rTxQ4Te8eZLQWpqKjt27Mh2tuN5e6J27doZxSa98LRq1YqyZcsG+i0B3uVGyBajrP79739zyy23cObMGQYMGEBYWBgLFizwbYC5ePrpp5k0aRLvv/9+sb7eHCrFKFhNnz6dhx9+mBEjRhAXF+eXb5qqysmTJ3MtVOk/Hzp0iNTU1Dz3FRERwe7du70+tpP50bx5c+3atSs9e/akZ8+eToRQZDldEqtQoQKjR4+mbt26GYVn06ZNGV8cSpUqRcuWLTOd6URGRlKvXr2guHwbFxdHXFwcs2bN2qGqzfNcOb8uGgI5FaRvuqy6d++uVapU0erVq2vv3r01Pj6+0PsqiGXLlqmI6KBBgwJyPCcRAt0BBavdu3drxYoV9ZZbbtG0tDSnw9GUlBSNj4/X77//XkUkx66i3AOqes3J/Ajl3EgXERGR498hfapZs6becMMNOmrUKH333Xf1hx9+0N9//93psL3iTW54/ccGRgDrgLPAnCzLwoE3gSO4hgZe5bFMgBeAo+5pCu4zsqxTqCXUgQMHtHbt2tq6dWs9c+aM0+H4XV4JZfmRu7S0NL355pu1YsWKumfPHqfDySa3D8GC9FuoasWoqPL6UrB///6g+BJTWN7kRkGGkIgHJgLv5LAsBqgBtHL/+7DHsiigN65RGtsCtwLZu0gIMampqfTr148zZ84wf/58wsPDnQ7JaZYfufjoo49YvHgxkyZNCsrulyZNmpQtf+1xhMDL6flIcDWjbtiwYVBcdvOr/KpV1gnXB84cj9ctgJNAlVzW/xaI8nh9P7Amp3VD6dvN2LFjFdA5c+Y4HUrA4M2ptuVHJkePHtU6deroZZddpikpKU6Hk6vY2FiNiIhQEdGIiIgCDS2Szpv88NcUirmR1fDhw7OdFRV0mJdg5dVnR34rZNsg+4dNf1wjMb6M6zLMRqCvx/ITQGeP152AUzntu3HjxtqxY8eMaebMmQH4NRXcl19+qSKiAwYMcDoUv5s5c2bG3wPYrZYfBTJo0CAtVaqUbtiwwelQ/KKg+eGvKdSLUUJCgtaqVUubNm2qjRs3LtKXgmAUqGL0lLuKjwfKAtcDp4FW7uWpQEuP9Zu71892XyAUEurgwYNat25dbdWqlZ4+fdrpcAKqkGdGJSo/PP33v/9VQJ944gmnQwkIOzMqvHvvvVdLly6tGzdudDoUv/AmN3wx7HgycB6YqKrnVHUlsBz4k3v5aaCKx/pVgNPuAENK+n2ikydPMn/+fCpWrOh0SKGgxOSHp+TkZIYOHUrTpk0ZO3as0+GYILZ06VLef/99nnzySUefcXKaL3pg+Cmf5Ztx3Zxe634d6Z4XciZNmsR///tfZs+endFDuMlXickPTxMnTmT79u18+eWX1huHydXp06cZOnQoLVu2JDo62ulwHOX1mZGIlBaR8kApoJSIlBeR0sAqYC/wpHudq4EuwFL3pnOBR0SkoYg0AEYDc3z4HgJixYoVPPPMM9xzzz2ZBmYzLiU9Pzxt3LiRKVOmMGDAAG688UanwzFB7O9//zt79uzh7bffzui+rMTK7zpe+oTrmn/W1h7j3ctaA6uBM8DPwF88thNcz44cc08h9xzJb7/9pvXr19cWLVroqVOnnA7HMeT9nFGJzQ9PKSkp2rlzZ61Vq5YeOXLE6XACKq/88PcUCrmR1Zo1a1REdPjw4U6H4nfe5IbXl+lUdbz7AyenZZuBK3NZpsBj7inkpKWlcc8993D8+HGWLFlCpUqVnA4pKJXU/MhqxowZfPfdd8TGxlKzZk2nwzFB6ty5cwwZMoSGDRsyefJkp8MJCiHZa3cgTZ48mS+++IKYmBjatm3rdDgmiO3bt48nn3yS7t27c/fddzsdjgliU6ZMYePGjcTFxVGlSpX8NygBfNGarthatWoVY8eO5a677grYeDcmNKkqDzzwAGlpacyYMaP4Py1vCm3r1q1MmDCB//u//+PWW291OpygYcUoFwkJCdx11100bdqUmTNn2oeLydO//vUv4uLiePbZZ7nwwgudDiekiEgTEflcRI6LyCERed3d+AURuVFEtopIkogsF5EIp+MtirS0NIYMGULFihV55ZVXnA4nqFgxykFaWhr9+/fn6NGjfPzxx1SuXNnpkEwQS0xM5MEHH6RDhw6MHDnS6XBC0ZvAYaA+0A7Xg9HDRaQW8Anwd1x9Gq4DPnIoRp+IiYnh66+/Ztq0aXkOalcS2T2jHLz44ossWbKEGTNm0K5dO6fDMUHu8ccf5/DhwyxatIjSpe2/VCFcCLyuqr8Dh0RkCa4WmH2Azao6H0BExgNHRKSlqm51LNpC2r9/P4899hg33XQTAwYMcDqcoGNnRll8/fXXREdHc8cddzB0aLHqPNr4wVdffUVMTAwPP/wwHTp0cDqcUPUKcKeIhItIQ+BmIL0g/Zi+kqqeAXa654cUVWX48OGkpKTYZf9cWDHycPToUe666y6aNGnCrFmzLGFMns6ePUtUVBRNmjThmWeecTqcULYSV4E5CezHdTluAVAJV0e6nk4A2a6bJyQk0KlTp4wpJibGvxEX0D//+U/i4uKYMGECF110kdPhBERMTEzG3wOold/6dk3BLS0tjQEDBnD48GFWr15tzS1NviZPnszWrVtZsmSJ9VNYSCIShqs3jpnAVbgK0Du4BlzM2m8h7tensu6ndu3aBOuQ9MeOHWPEiBF07NixRN1TjIqKIioqCgAROZLf+nZm5DZt2jQWLVrEtGnT7HKLydfPP//Mc889R79+/ejevbvT4YSyGsAFuO4ZnVXVo8C7QA/+6LcQABGpCDQlxPouHDNmDEePHuXtt9+2e4p5sGIErF69mieffJLbb7+d4cOHOx2OCXJpaWlERUVRuXJlpk2b5nQ4IU1VjwC7gL+5+y6sBgzAda/o30AbEenr7vdwLPBTKDVe+PLLL3n33Xd57LHHrDFUPkp8MTp27Bh33nknF1xwAW+//bbdJzL5iomJ4ZtvvmHatGnUqVPH6XCKgz7An4EEYAeQAjysqglAX2AScBzoDNzpVJAFlZSUxNChQ2nevDl///vfnQ4n6AXVOeOJEyeIioqiZ8+e9OzZ0+/HU1Xuu+8+Dh48yLfffkvVqlX9fsxQFBcXR1xcHICjv6BA50dO4uPjefzxx7nxxhvp37+/IzEEm6Lmh6puwNWTe07LvgRaFjY2J40fP55ff/2VFStW2DAi3sivJ9VAToHueXfatGkK6PTp0wN63FCFg70ya5D0zNynTx8tX768bt++3elQgo6T+REMueFp3bp1GhYWpkOGDHE6lKDgTW4E1ZlRIK1du5bHH3+c3r1789BDDzkdjgkBCxYs4JNPPmHy5Mk0a9bM6XBMkDp//jyDBw+mbt26TJkyxelwQkaJLEbHjx/njjvuoGHDhrzzzjt2n8jk6+TJk4wYMYK2bdsyevRop8MxQWzatGls2LCBTz75hGrVqjkdTsgoccVIVRk0aBAHDhzg66+/pnr16k6HZELAU089RXx8PJ988gllypRxOhwTpLZv38748ePp06cPf/nLX5wOJ6SUuGL02muvsWDBAqZOnUrnzp2dDseEgNWrV/Pmm2/y0EMPcfnllzsdjglSqkpUVBTlypXjtddeczqckFOiitG6desYM2YMPXv25OGHH3Y6HBMC0kfkbNSoERMmTHA6HBPEZs+ezYoVK4iJiaFBgwZOhxNySkwxSkxM5I477qB+/frMmTPH7hMZr0yZMoXNmzfz2Wef2VAiJlcHDx5kzJgxdOnSxQbiLKQSUYxUlcGDB7Nv3z5WrVpFjRo1nA7JhIBt27YxYcIE7rjjDm655RanwzFB7MEHH+T3338nJibGvugWUokoRm+++Sb/+te/mDJlCldeeaXT4ZgQkJaWxtChQwkPD7cROU2e/v3vf/Ovf/2LyZMn07x5c6fDCVnFvhj98MMPPPLII/To0cOa5Bqvvfvuu6xcuZJZs2ZRr149p8MxQSoxMZEHHniAdu3a2edLERXrYnTy5EnuuOMO6tSpw3vvvUdYWInvis944dChQ4wZM4brr7+e+++/3+lwTBB7/PHH+e2334iLi7Mm/0VUbD+d05tZ7t69m3/84x/UqpXv2E7GADBq1CiSk5NtRM4Qkt5vobuPvIBYuXIlMTExPPLII3Ts2DFgxw0lcXFx6WMa5d9vYX79BQVy8mX/UjNmzFBAJ0+e7LN9lnSUgL7pPvvsMwV0woQJfj9WceNkfgS6b7qkpCRt3ry5XnTRRXrmzJmAHjsUeZMbXp8ZicgIEVknImdFZE4u64wTERWRmzzmiYi8ICJH3dMU8fPXzQ0bNjBq1Cj+/Oc/89hjj/nzUMYtlPIjN6dPn+Zvf/sbrVu3trwxeZowYQLbt29n5syZhIeHOx1OsVCQe0bxwESgO5CtP3QRaQrcDhzMsigK6I1rxEYFvgB+Bd4qeLj5O3XqFHfccQc1a9Zk7ty5dp8ocEIiP/Ly9NNPs3//fr755hvKli0b6MObEPHjjz8yZcoUBg4cyE033ZT/BsYrXn9Sq+onqroAOJrLKq8DjwPnsswfAExV1f2qegCYCgwseKhexciwYcPYuXMnH374IbVr1/bHYUwOQiE/8rJ27VpeffVV/va3v1nzf5OrlJQUBg8eTM2aNZk6darT4RQrPmlNJyJ/Bc6p6uc5XGFpjWsI4XQ/uudlk5CQQKdOnTJeR0VFpd/88srs2bP54IMPmDhxItddd53X25ncxcTEEBMTk/6yUK1AgiU/cnP+/HmGDBlCgwYNmDx5cpH3V5L4Ij9CySuvvMK6dev46KOP7OF5X8vvplLWCdelmDkerysB24EL3a93Azd5LE8FWnq8bo7rcoxk3XdhbkLGxsZqRESEiogC2qZNG01NTS3wfkz+8OImZLDlhzeef/55BfTf//63X/ZfUniTH/6aAtGAYefOnVqhQgXt2bOnpqWl+f14xYk3ueGLGyrPAO+r6q5clp8Gqni8rgKcdgdYJPPmzSMqKoo9e/akf5BlXKIzQcOx/PDGjh07Mrr87927dyAOaUKQqjJ06FBKly7Nm2++aU3+/cAXxehG4CEROSQih4ALgI9F5HH38s24bk6ni3TPK7Lo6GiSkpIyzUtOTiY6OtoXuze+4Vh+5Efd9xjLli1rXf6bPM2dO5cvv/ySF154gUaNGjkdTrHk9T0jESntXr8UUEpEygMpuD5sPB89/h54BFjsfj0XeEREPsd1+WU04JP/+Xv37i3QfOM/wZgf+Zk7dy7/+c9/mDFjhnX5b3L122+/8fDDD3P11VczdOhQp8MptgpyZvQ0kAw8Adzj/vlpVT2qqofSJ1z3AI6r6mn3djOBOGAjsAlY5J5XZI0bNy7QfONXQZcfeUlISOCRRx7h6quv9kkjCFN8jRw5kjNnzjBr1ix7VMSf8rupFMipoDchY2NjNTw8XHF9o1ZAw8PDNTY2tkD7Md6hGPXA0K9fPy1Tpoxu3rzZZ/ss6ZzMD381YFi4cKEC+uyzz/pl/yWFN7kR0mW+X79+xMTEEBERgYgQERFBTEwM/fr1czo0E8SWLl3KvHnzePLJJ7nkkkucDscEqZMnTzJ8+HDatGnD448/nv8GpkhCvtfufv36WfExXjtz5gzDhg2jRYsWPPXUU06HY4LYk08+yYEDB/jnP/9pPXIEQEifGRlTUOPHj2f37t3MmjWLcuXKOR2OcRORO0Vki4icEZGdInKte/6NIrJVRJJEZLmIRAQinm+++YY333yThx56iM6dOwfikCWeFSNTYvzwww9MmzaNqKgorr32WqfDMW4i0g14AbgPqAxcB/wqIrWAT4C/AzWAdcBH/o7n7NmzDB48mMaNGzNx4kR/H864hfxlOmO8kZKSwpAhQ6hTpw4vvPCC0+GYzJ4BnlXVNe7XBwBEJArYrKrz3a/HA0dEpKWqbvVXMM899xxbt25l8eLFVKpUyV+HMVnYmZEpEV555RV++OEHXnvtNapVq+Z0OMZNREoBnYDaIrJDRPaLyOsiUoEs/Raq6hlgJ7n0XegLmzZtYvLkydxzzz38+c9/9tdhTA7szMgUa/PmzeOxxx4jPj6eChUqcPbsWadDMpnVxfVQ9O3AtcB54FNcz61VAhKyrH8C16W8THzRiW5qaiqDBw+matWqvPzyywXa1mRX0E50rRiZYiu978L0LqOSk5MzPqCsBWbQSHb/+5qqHgQQkWm4itEqMvdbiPv1qaw7qV27NuvWrStSIG+88QbfffcdsbGx1KpV7Dsg9zvPLwQiciS/9Yv1Zbrjx49Tvnx5Ro0alWn+nDlzuP3223PcZtasWbRu3ZpWrVpx8cUXM2nSJNLS0jKWv/nmm1x66aVERkbSsmVLxowZk+N+fvnlF7p27UrLli1p06YN9913H8nJyTmum9XAgQN5/fXXvXuTJlc59V2YlJRkfRcGEVU9DuzH9dB6Vpn6LRSRikBT/NB34Z49e3jqqae4+eabufvuu329e+OFYl2M5s2bx5VXXsmHH37IuXNZx3TL7v3332f69OksXryYLVu28N1337Fs2TKee+45AL7//ntefvllvvrqK3788Uc2b95M//79c9xX2bJlmTZtGlu3buWnn34iKSmJl156yafvz+TN+i4MGe8CD4pIHRGpDowCPgP+DbQRkb7uvg7HAj/5uvGCujvMBZgxY4b1yO2QYl2M3nnnHZ5++mkuvfRSFi5cmO/648aNY+rUqRl921WvXp233nqL5557jrNnz7J//36qVq2a0cKmVKlStG3bNsd9NWnShPbt2wMQFhbG5Zdfzp49e3Jc98CBA9x4441ERkbSu3dvjhz544z25MmTDB48mMsvv5y2bdsycuRIUlNTAXjmmWdo2bIl7dq1o3379iQmJgKwevVqrrnmGiIjI4mMjGTZsmXe/cKKGeu7MGRMwNWB7i/AFuB/wCRVTQD6ApOA40Bn4E5fHXTevHk0adKEsLAwlixZQp8+fYiICMhjTCYn+fUXFMipWbNmOmTIEF24cGGR+0LasGGDRkREaFpamsbGxurNN9+csezdd9/Vvn37Zlr/5MmTCujx48ez7atq1aq6ceNGPX36tHbu3FkbNGigd911l86cOVPPnDmTbyxJSUl6ySWX6Keffprj8j59+uj48eNV1TWAV6VKlfS1115TVdX7779f586dq6qqqampeuedd2pMTIweO3ZMK1WqpElJSRnxnz9/Xo8ePap169bVb775RlVVU1JS9NixY/nGmJeFCxfqkCFDFNiuIZQf1ndhYARDfli/lsENbwbmzG+FQE6+7OzwwQcf1L///e+q6ioG1atX1/3796tqwYtRlSpVdOPGjaqqmpaWpmvXrtWpU6fqZZddpm3bttWzZ8/mGsf58+e1V69eOmLEiFzX8YxNVfW2227LKEa1a9fWSy+9VCMjIzUyMlKbN2+uY8aM0ZSUFO3YsaP26dNHY2JiND4+XlVVP/vsM+3atasXv6GC8yah/DkVdSTgiIgI+7DxIyfzo6C5ERERkakQpU8RERGFffsmD97kRrFsTXfu3Dk++OADypcvz9y5cwE4f/487733Xq79kVWuXJkLL7yQNWvWZHq+YMuWLZw/f55mzZoBICJcdtllXHbZZYwYMYI6deqwadMmFi9ezPz58wF4+eWX6dq1K6mpqfTr14/q1avz6quvFuq9qCoLFizgoosuyrZszZo1fPPNN/z3v/+lY8eOLFmyxPUNw2SwvgtNTux+YvAplveMFixYQMuWLdm/fz+7d+9m9+7dLFu2jHfffTfP7caNG8eYMWPYt28f4GqN97e//Y0nnniC8uXLs3XrVjZt2pSx/rZt2zh37hyNGjUiOjqaDRs2sGHDBrp27UpaWhoDBw6kVKlSzJ49O8+bojfccENGbLt27eI///lPxrJevXrx/PPPZ9wnOnLkCLt27eLUqVMkJCRw/fXX88wzz9CmTRs2bdrEVVddxc8//8zq1asB17MTx48fL9wv0phiyu4nBp9ieWb07rvvZvs2fOWVV5KWlsaqVasA+PzzzzMNH3zfffcxYcIEkpOT6d69O6pKSkoK/fv3z2gKnJSUxKhRozh8+DDly5enVKlSxMbGUqdOnWwxLF68mNjYWNq0aUPHjh0BuPrqq3njjTeyrfvKK6/Qv39/5s+fT4sWLejWrVvGsunTp/PYY48RGRmJiFCuXDmmT59OmTJl6Nu3L8nJyaSlpdGhQwf69OlD+fLl+eSTT3jkkUc4c+YMYWFhvPTSS9x0000MHjyYXr160atXr6L/ko0JYZMmTcr0DBpAeHg4kyZNcjCqkk2C6bJOp06dtKgPrhn/EZH1qtop/zX9w/IjuDmZH4XJjXnz5hEdHc3evXtp3LgxkyZNsku6fuJNbhTLMyNjjMmP3U8MLsXynpExxpjQYsXIGGOM46wYGWOMcZwVI2OMMY6zYmSMMcZxVoyMMSHvxIkTREVFERcX53QoxkNcXFz6mEZV81vXmnYbY0Je1apVPUcVNUGiZ8+e9OzZk1mzZp3Ib12vz4xEZISIrBORsyIyx2P+FSLyhYgcE5EEEZkvIvU9louIvCAiR93TFPHDgCHBmIjBGBP4Jy7Lj4KzmIJDML7nkhhTQS7TxQMTgXeyzK8OxABNgAhcQwJ7dgIXBfTGNWJjW+BWYGihos1DSfzjFZaf4rL8KCCLKTgE43suiTF5fZlOVT8BEJFOQCOP+Ys91xOR14GVHrMGAFNVdb97+VRgCPBW4cM2wcbywxhTFP64Z3Qdmceobw386PH6R/e8bNavX39aRDzP1hKAIzmtm4NaIuLtuoESjDFBweKqBdR2/9zCB8e2/PhDcYjJ1/lRKOvXrz8iIjkPpZy/4vB3CISixJTvELo+LUYi0hbXOPW3ecyuBHjevDoBVBIR0Sy9tKpqZV/GY4KL5YfxF1Wtnf9aJpj5rGm3iDQDFgMjVfUrj0WngSoer6sAp7N+0JjizfLDGJMXnxQjEYkAvgQmqOr7WRZvxnVzOl0kmS/TmGLO8sMYk5+CNO0uLSLlgVJAKREp757XEPgv8Iaq5nTTeS7wiIg0FJEGwGhgjg9iD1oiclpEso8THvg45ojIxAAdy/LDSyUxP4x3SnJuFOTM6GkgGXgCuMf989PAYOAiYJz7F3laRE57bDcTiAM2ApuARe55BSYiu0XkNxGp6DFvsIisKMz+fMEdU3KW936xqv7q5+OuEJHjIlLOn8cpAMuP3GOy/HCY5Uam4wZlbnhdjFR1vKpKlmm8qj7j/rmS5+SxnarqY6pawz09VsT7AaWBkUXY3h96Znn/8f48mIg0Aa4FFAiKMcQtP/JU4vMjSFhuBHFuhGLfdC8CY0SkWtYFInKViHwvIifc/14V+PBARFREmolIWRHZICIPuueXEpFvRGRsEQ/RH1iD63LWAI/jtheRH0TklIh8BJT3WFZdRD4TVy8Ix90/N/JYvkJEJovIWvfv71MRqVHEOJ1g+WH5kRvLjSDOjVAsRuuAFcAYz5nuN78IeBWoCUwDFolIzUAHmE5Vz+G6ZPWsiLTCdQmrFDCpiLvuD8xzT91FpK6IlAUWAO8DNYD5QF+PbcJw9XwQATTGdRnt9Rz2OwhoAKTg+l2GGssPy4/cWG4Ec26oashMwG7gJqANrudRauO6J7ECuBdYm2X91cDAAMR0Gkh0TwtwnQI381hnNLAVOA40L+LxrgHOA7Xcr7cCD+N6mDQeEI91vwUm5rKfdsBxj9crgOc9Xl8CnANKOf13t/yw/LDcKP65EYpnRqjqJuAzXN8W0jUAsj6BvQdoGICQeqtqNffUO4fl7+Hqm+1zVd1exGMNAJapavqT0B+45zUADqg7G9wyfh8iEi4iM0Vkj4icBFYB1USklMf6+7JsWwbXE/YhxfLD8iM3lhvBmxshWYzcxuHqwyw9YeLJ3uVEY+BAIIPKxZu4/gN0F5FrCrsTEakA3AFcLyKHROQQrm82kcBBoKFIph6vG3v8PBpXdy2dVbUKrm9DAJ7rX5Bl2/N4391OsLH8sPzIjeVGEOZGyBYjVd0BfAQ85J71OXCxiNwtrudb/g/X6eJnTsUIICL3Ah2BgbhifU9EKuW5Ue56A6m43lc799QK+Mq9LAV4yP3++wCXe2xbGde13kT3NfJxOez/HhG5RETCgWeBf6pqaiFjdZTlh+VHbiw3gjQ3An3t1gfXWG/yeH0B8DuwwuOa6Hpc14TXA9cEOib3PAWa4fqGcBS42mPZR8CsQh5rCa4errPOvwM4BHQC/odrmIaP3NNE9zoNcF3bPQ38gmuYBgVK6x/XfScDa4GTuJ79qeVxjNPAtU7ngOWH5YflRvHMDXGvbEo4cT38F6uqbzsdiwk+lh8mN77KjZC9TGeMMab4sGJkjDHGcXaZzhhjjOPszMgYY4zjrBgZY4xxnBUjY4wxjrNiZEwORGSgiHzt8TooBj0rKPHhIGkiMl5EYn2xL2OysmJkQp6IXCMi37q7rz/m7mr/Ml8eQ11jzfh10LOiylpAjQklpZ0OwJiiEJEquLpt+RvwMVAW1+BhZ52MyxhTMHZmZELdxQCq+qGqpqpqsqouU9Wf0lcQkSEissU9cNjPItLBPf8JEdnpMf8vuR1E3IOeuX+eIyJviMgi97bfiUhTj3X/JCLb3Gdqb4rIShEZnMt+x4vIfBGJde9ro4hcLCJPishhEdknIn/yWL+qiMwWkYMickBEJopr4LVWwFvAle5Liokeh6meR6y5DionIhe6Yz8lIl8QYj10m9BixciEul+AVBF5T0RuFpHqngtF5K/AeFyDf1XBNdTyUffinbjOoqoCzwCxIlLfy+Pe5d6mOrAD96BnIlIL+CfwJK6B2rYB+Y0a2hPXwGbVcfUPthTX/82GuDqdnOmx7nu4OrVsBrQH/gQMVtUtwDBgtfuSYjUvYs1vULkPcPXTVguYgMfIoMb4mhUjE9JU9SSuTi4VmAUkiMhCEanrXmUwMEVVv1eXHaq6x73tfFWNV9U0Vf0I2E7m3orz8omqrlXVFFyjZrZzz+8BbFbVT9zLXsXVEWVevlLVpe715+Ma+O15VT0P/ANoIiLV3O/pZmCUqp5R1cPAy8CdhYz1FmC7qr6vqimq+iGuAdd6ikhj4DLg76p6VlVX4eoA0xi/sGJkQp6qblHVgaraCNdIng2A6e7FF+A6A8pGRPqLyAYRSXRf1mqD95eiPAtMEpDetX8DPAYaU1cXJ/vz2ddvHj8nA0f0j+73k93/VsI15k4Z4KBHzDOBOkWINbdB5RrgGs3zTJZlxviFFSNTrKjqVmAOrsICrsLQNOt6IhKB60xqBFDTfVlrE5kHDCuMg0Ajj+OI5+si2oerYUYt/WN00Cqq2tq9vKB9e+U1qNxBXPeaKmZZZoxfWDEyIU1EWorIaBFp5H59Aa57JGvcq7wNjBGRjuLSzF2IKuL68E5wb3cffxSwolgEXCoivUWkNPAAUM8H+0VVDwLLgKkiUkVEwkSkqYhc717lN6CRiJT1cpe5DirnvpS5DnhGRMqKa5TRnr54H8bkxIqRCXWngM7AdyJyBlcR2oRrqGRUdT6uG/YfuNddANRQ1Z+BqcBqXB/ilwLfFDUYVT0C/BWYgquhxCW4PtR91dS8P67m6z8Dx3E1lkhvdPFfYDNwSETyHfJZVY8Ct+L6XR0FHgNudb8HgLtx/W6P4Rrdc66P3oMx2Viv3cb4kYiE4bpn1E9VlzsdjzHBys6MjPExEenubv1WDngK132oNflsZkyJZsXIGN+7ElcLviO47rP0VtXkvDcxpmSzy3TGGGMcZ2dGxhhjHGfFyBhjjOOsGBljjHGcFSNjjDGOs2JkjDHGcVaMjDHGOO7/AcVa9NFIcVe0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x396 with 9 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = [0,1,2]\n",
    "fig, axs = plt.subplots(nrows=3, ncols=3, figsize=[6, 5.5], sharex=True)\n",
    "for i in range(3):\n",
    "    tecs = list(rDicts[i].values())\n",
    "    clr, jlr, jhr = tecs[1:4], tecs[4:7], tecs[7:10]\n",
    "    for j, tec in enumerate([tecs[1:4], tecs[4:7], tecs[7:10]]):\n",
    "        axs[i, j].plot(x, tec, 'o-', color='k')\n",
    "# axis format\n",
    "nums = ['(a)', '(b)', '(c)',\n",
    "        '(d)', '(e)', '(f)',\n",
    "        '(g)', '(h)', '(i)']\n",
    "gnames = ['CLR GIM', 'JLR GIM', 'JHR GIM']\n",
    "for i, ax in enumerate(axs.flatten()):\n",
    "    ax.tick_params(which='both', direction='in', top=True, bottom=True, left=True, right=True)\n",
    "    ax.set_xlim(-0.4, 2.4)\n",
    "    ax.annotate(f'{nums[i]}', xy=(0.08, 0.8), xycoords='axes fraction', ha='left')\n",
    "for ax in axs[0,:2]:  ax.set_ylim(15.5, 23.5); ax.yaxis.set_major_locator(ticker.MultipleLocator(4));  ax.yaxis.set_minor_locator(ticker.MultipleLocator(1))\n",
    "for ax in axs[0,2:]:  ax.set_ylim(5.5 , 9.5 ); ax.yaxis.set_major_locator(ticker.MultipleLocator(2));  ax.yaxis.set_minor_locator(ticker.MultipleLocator(0.5))\n",
    "for ax in axs[1,: ]:  ax.set_ylim(4.5,  6.5 ); ax.yaxis.set_major_locator(ticker.MultipleLocator(1));  ax.yaxis.set_minor_locator(ticker.MultipleLocator(0.25))\n",
    "for ax in axs[2,:2]:  ax.set_ylim(120 , 175 ); ax.yaxis.set_major_locator(ticker.MultipleLocator(20)); ax.yaxis.set_minor_locator(ticker.MultipleLocator(10))\n",
    "for ax in axs[2,2:]:  ax.set_ylim(47  , 102 ); ax.yaxis.set_major_locator(ticker.MultipleLocator(20)); ax.yaxis.set_minor_locator(ticker.MultipleLocator(10))\n",
    "for i in range(3):    axs[0,i].set_title(gnames[i])\n",
    "axs[1,0].set_ylabel('RMSE [cm]')\n",
    "#plt.xticks(x, ['no\\nscale', 'scale\\nw/\\nfixed\\nratio', 'scale\\nw/\\nadap.\\nratio'])\n",
    "plt.xticks(x, ['No', 'Fix', 'Adap.'])\n",
    "axs[2,1].set_xlabel('Scaling method')\n",
    "fig.tight_layout()\n",
    "fig.subplots_adjust(wspace=0.3)\n",
    "axs[0,0].annotate('Sentinel-1 asc.',  fontsize=11, xy=(0.52, 0.15), xycoords='axes fraction', ha='center')\n",
    "axs[1,0].annotate('Sentinel-1 desc.', fontsize=11, xy=(0.52, 0.55), xycoords='axes fraction', ha='center')\n",
    "axs[2,0].annotate('ALOS-2 desc.',     fontsize=11, xy=(0.52, 0.13), xycoords='axes fraction', ha='center')\n",
    "\n",
    "# outputs\n",
    "out_fig = os.path.join(work_dir, 'topTEC_scale_rmse.pdf')\n",
    "print('save figure to file', out_fig)\n",
    "plt.savefig(out_fig, bbox_inches='tight', transparent=True, dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take home message:\n",
    "\n",
    "1. For Sentinel-1, scaling is better than no scaling, except for JHR GIM at Sentinel-1 ascending orbit.\n",
    "2. For Sentinel-1, adaptive scaling is better than fixed scaling\n",
    "3. For ALOS-2, different altitude, different local time, applying ratio get worse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}