1_Y2O3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BBO-Rietveld Y2O3\n",
    "\n",
    "In this notebook, you can see how to automate Rietveld refinement with black-box optimization.  \n",
    "Run the following cells to run the optimization.\n",
    "\n",
    "## Documentations\n",
    "- [JupyterLab](https://jupyterlab.readthedocs.io/en/stable/)\n",
    "- [GSASIIscriptable: Scripting Interface](https://gsas-ii.readthedocs.io/en/latest/GSASIIscriptable.html)\n",
    "- [Optuna](https://optuna.readthedocs.io/en/stable/)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import packages\n",
    "%matplotlib inline\n",
    "\n",
    "import os\n",
    "import sys\n",
    "from multiprocessing import Process, Queue\n",
    "import pandas as pd\n",
    "import optuna\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.gridspec import GridSpec\n",
    "from dataclasses import dataclass\n",
    "sys.path.append('/root/g2full/GSASII/')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configurations\n",
    "RANDOM_SEED = 1024\n",
    "STUDY_NAME = 'Y2O3'\n",
    "DATA_DIR = 'data/' + STUDY_NAME\n",
    "# all output files include GSAS project file (*.gpx) will be saved in WORK_DIR\n",
    "WORK_DIR = 'work/' + STUDY_NAME"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make directories\n",
    "! rm -f $WORK_DIR/$STUDY_NAME*\n",
    "! mkdir -p $WORK_DIR"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "@dataclass\n",
    "class ProjectConfig:\n",
    "    work_dir: str\n",
    "    random_seed: int\n",
    "    data_dir: str\n",
    "    cif_file: str\n",
    "    csv_file: str\n",
    "    prm_file: str"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "class Project:\n",
    "    def __init__(self, config: ProjectConfig, trial_number: int):\n",
    "        import GSASIIscriptable as G2sc\n",
    "        self.gpx = G2sc.G2Project(newgpx=os.path.join(\n",
    "            config.work_dir, f'project_seed{config.random_seed}_trial_{trial_number}.gpx'))\n",
    "\n",
    "        self.hist1 = self.gpx.add_powder_histogram(\n",
    "            os.path.join(config.data_dir, config.csv_file), os.path.join(config.data_dir, config.prm_file))\n",
    "        \n",
    "        self.phase0 = self.gpx.add_phase(os.path.join(config.data_dir, config.cif_file),\n",
    "                                         phasename=config.cif_file.split(\".cif\")[0],\n",
    "                                         histograms=[self.hist1])\n",
    "        \n",
    "        self.hist1.data['Instrument Parameters'][0]['I(L2)/I(L1)'] = [0.5, 0.5, 0]\n",
    "\n",
    "        # Set to use iso\n",
    "        for val in self.phase0.data['Atoms']:\n",
    "            val[9] = 'I'\n",
    "\n",
    "    def refine_and_calc_Rwp(self, param_dict):\n",
    "        self.gpx.do_refinements([param_dict])\n",
    "        for hist in self.gpx.histograms():\n",
    "            _, Rwp = hist.name, hist.get_wR()\n",
    "        return Rwp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example usage:\n",
    "config = ProjectConfig(\n",
    "    work_dir=WORK_DIR,\n",
    "    random_seed=RANDOM_SEED,\n",
    "    data_dir=DATA_DIR,\n",
    "    cif_file='Y2O3.cif',\n",
    "    csv_file='Y2O3.csv',\n",
    "    prm_file='INST_XRY.PRM'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the objective function, first, configurations for the following parameters are sampled.\n",
    "Then, the crystal structure is refined by GSAS-II.\n",
    "Finally, the calculated Rwp is returned.\n",
    "\n",
    "| Parameter                                 | Parameter type   | Search space                                                                                                      |\n",
    "|-------------------------------------------|------------------|-------------------------------------------------------------------------------------------------------------------|\n",
    "| initial 2theta bound refinement           | categorical      | {True, False}                                                                                                      |\n",
    "| initial 2theta lower bound (deg.)         | real             | [15, 130]                                                                                                          |\n",
    "| initial 2theta upper bound (deg.)         | real/conditional | [l + 20, 150] (l is the initial 2theta lower bound)                                                                |\n",
    "| background function                       | categorical      | {Chebyshev, cosine, Q^2 power series, Q^-2 power series, linear interpolate, inverse interpolate, log interpolate} |\n",
    "| degree of background function             | initegral        | [1, 15]                                                                                                            |\n",
    "| instrument zero position refinement       | categorical      | {True, False}                                                                                                      |\n",
    "| sample shift X refinement                 | categorical      | {True, False}                                                                                                      |\n",
    "| sample shift Y refinement                 | categorical      | {True, False}                                                                                                      |\n",
    "| scale factor refinement                   | categorical      | {True, False}                                                                                                      |\n",
    "| peak shape refinement (Gaussian: U, V, W) | categorical      | {True, False}                                                                                                      |\n",
    "| peak shape refinement (Lorentzian: X, Y)  | categorical      | {True, False}                                                                                                      |\n",
    "| peak shape refinement (SH/L)              | categorical      | {True, False}                                                                                                      |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def objective(trial, config):\n",
    "    \"\"\"\n",
    "    objective function for Optuna.\n",
    "    \n",
    "    You can apply black-box optimization for other data analysis \n",
    "    by modifying this function for that task.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    trial : optuna.trial object\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    Rwp : float\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    ### define search space ###\n",
    "    # Limits (acute angle)\n",
    "    limits_lb = trial.suggest_float('Limits lower bound', 15, 130)\n",
    "    limits_ub = trial.suggest_float('Limits upper bound', limits_lb + 20, 150)\n",
    "    limits_refine = trial.suggest_categorical('limits refine', [True, False])\n",
    "    refdict0 = {'set': {'Limits': [limits_lb, limits_ub]}, 'refine': limits_refine}\n",
    "\n",
    "    # Background\n",
    "    background_type = trial.suggest_categorical(\n",
    "        'Background type', ['chebyschev',\n",
    "                            'cosine',\n",
    "                            'Q^2 power series',\n",
    "                            'Q^-2 power series',\n",
    "                            'lin interpolate',\n",
    "                            'inv interpolate',\n",
    "                            'log interpolate'])\n",
    "    no_coeffs = trial.suggest_int('Number of coefficients', 1, 15)  # [1, 15]\n",
    "    background_refine = trial.suggest_categorical('Background refine', [True, False])\n",
    "    refdict0bg_h = {\n",
    "        'set': {\n",
    "            'Background': {\n",
    "                'type': background_type,\n",
    "                'no. coeffs': no_coeffs,\n",
    "                'refine': background_refine\n",
    "            }\n",
    "        }\n",
    "    }\n",
    "\n",
    "    # Instrument parameters\n",
    "    instrument_parameters1_refine = []\n",
    "    for p in ['Zero']:\n",
    "        if trial.suggest_categorical('Instrument_parameters refine %s' % (p), [True, False]):\n",
    "            instrument_parameters1_refine.append(p)\n",
    "    refdict1_h = {'set': {'Cell': True, 'Instrument Parameters': instrument_parameters1_refine}}\n",
    "\n",
    "    sample_parameters1_refine =[]\n",
    "    for p in ['DisplaceX', 'DisplaceY', 'Scale']:\n",
    "        if trial.suggest_categorical('Sample_parameters refine %s' % (p), [True, False]):\n",
    "            sample_parameters1_refine.append(p)\n",
    "    refdict1_h2 = {\"set\": {'Sample Parameters':sample_parameters1_refine }}\n",
    "\n",
    "    instrument_parameters2_refine = []\n",
    "    for p in ['U', 'V', 'W', 'X', 'Y', 'SH/L']:\n",
    "        if trial.suggest_categorical('Peakshape_parameters refine %s' % (p), [True, False]):\n",
    "            instrument_parameters2_refine.append(p)\n",
    "    refdict2_h = {'set': {'Instrument Parameters': instrument_parameters2_refine}}\n",
    "\n",
    "    refdict3_h = {'set': {'Atoms': {'all': 'XU'}}}\n",
    "\n",
    "    # Limits (wide angle)\n",
    "    refdict_fin_h = {'set': {'Limits': [15, 150]}, 'refine': True}\n",
    "\n",
    "    # Evaluate\n",
    "    refine_params_list = [refdict0,\n",
    "                          refdict0bg_h,\n",
    "                          refdict1_h,\n",
    "                          refdict1_h2,\n",
    "                          refdict2_h,\n",
    "                          refdict3_h,\n",
    "                          refdict_fin_h]\n",
    "    \n",
    "    def evaluate(config, trial_number, refine_params_list):\n",
    "        ERROR_PENALTY = 1e9      \n",
    "        try:\n",
    "            print(config)\n",
    "            print(trial_number)\n",
    "            print(refine_params_list)\n",
    "            # todo: ここを実装（trial_number）を渡す\n",
    "            project = Project(config, trial_number)\n",
    "            for params in refine_params_list:\n",
    "                Rwp = project.refine_and_calc_Rwp(params)                    \n",
    "            # validate Uiso >= 0\n",
    "            phase_Y2O3 = project.gpx.phases()[0]\n",
    "            u_iso_list = [atom.uiso for atom in phase_Y2O3.atoms()]\n",
    "            if min(u_iso_list) < 0:\n",
    "                # Uiso < 0\n",
    "                Rwp = ERROR_PENALTY\n",
    "            return Rwp\n",
    "            \n",
    "        except Exception as e:\n",
    "            # Refinement failed\n",
    "            print(e, file=sys.stderr)\n",
    "            return ERROR_PENALTY\n",
    "            \n",
    "    Rwp = evaluate(config, trial.number, refine_params_list)\n",
    "    \n",
    "    return Rwp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[I 2024-04-03 12:39:27,443] Using an existing study with name 'Y2O3_seed1024' instead of creating a new one.\n"
     ]
    }
   ],
   "source": [
    "# Create Optuna study\n",
    "\n",
    "study = optuna.create_study(study_name=STUDY_NAME + '_seed%s' % (RANDOM_SEED),\n",
    "                            storage=f\"sqlite:///{config.work_dir}/history_sqlite.db\", \n",
    "                            load_if_exists=True,\n",
    "                            sampler=optuna.samplers.TPESampler(n_startup_trials=20, seed=RANDOM_SEED))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run 200 refinements to find the best configuration. It may take abount an hour to complete."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Optimize\n",
    "study.optimize(func=lambda trial: objective(trial, config), n_trials=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Results\n",
    "df = study.trials_dataframe()\n",
    "df.columns = [''.join(col).replace('params', '').strip() for col in df.columns.values]\n",
    "df.rename(columns={'value':'Rwp', 'number':'trial'}, inplace=True)\n",
    "df.sort_values('Rwp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'Limits lower bound': 91.81741956486644,\n",
       " 'Limits upper bound': 124.16688310979242,\n",
       " 'limits refine': False,\n",
       " 'Background type': 'log interpolate',\n",
       " 'Number of coefficients': 14,\n",
       " 'Background refine': True,\n",
       " 'Instrument_parameters refine Zero': True,\n",
       " 'Sample_parameters refine DisplaceX': True,\n",
       " 'Sample_parameters refine DisplaceY': True,\n",
       " 'Sample_parameters refine Scale': False,\n",
       " 'Peakshape_parameters refine U': False,\n",
       " 'Peakshape_parameters refine V': False,\n",
       " 'Peakshape_parameters refine W': True,\n",
       " 'Peakshape_parameters refine X': True,\n",
       " 'Peakshape_parameters refine Y': True,\n",
       " 'Peakshape_parameters refine SH/L': False}"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Best configuration\n",
    "study.best_params"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7.946615605079828"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Best Rwp\n",
    "study.best_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAG2CAYAAACDLKdOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAt1ElEQVR4nO3deXxU9b3/8fdkkkxCQnZNiJAQEWULiwKKYAVFES1IFxWLNMjvoVVRpCgqvUUULrJUuSClWr1VWqu4XIWL3oJaFCgu7CFR2Q2CQKAISUgCIZn5/v7QjMaENbN+83o+HvOQOXPmO5/hg+TN93zPOQ5jjBEAAIAlIoJdAAAAgC8RbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVUIq3KxYsUKDBg1SZmamHA6HFi5cWG+fTZs2afDgwUpMTFRcXJx69OihXbt2Bb5YAAAQkkIq3FRUVKhLly6aO3dug6/v2LFDffr0Ubt27bRs2TIVFBRowoQJiomJCXClAAAgVDlC9caZDodDCxYs0JAhQ7zbhg4dqqioKL300kvBKwwAAIS0yGAXcLo8Ho/+7//+Tw899JAGDBigDRs2KCcnR+PHj68TgH6sqqpKVVVVdcY5dOiQUlNT5XA4AlA5AABoLGOMjhw5oszMTEVEnOLAkwlRksyCBQu8z/ft22ckmWbNmpmZM2eaDRs2mKlTpxqHw2GWLVt2wnEmTpxoJPHgwYMHDx48LHjs3r37lBkibA5L7d27V+edd55uvfVWvfLKK979Bg8erLi4OM2fP7/BcX48c1NaWqqsrCzt3LlTCQkJPq3Z7Xbr888/V8eOHeV0On06NnyDHoU+ehTa6E/os7VHZWVlat26tUpKSpSYmHjSfcPmsFRaWpoiIyPVoUOHOtvbt2+vlStXnvB9LpdLLper3vbk5GS/hJv4+HglJydb9QfKJvQo9NGj0EZ/Qp+tPar9LqezpCSkzpY6mejoaPXo0UNbtmyps33r1q3Kzs4OUlUAACDUhNTMTXl5ubZv3+59XlRUpPz8fKWkpCgrK0vjxo3TLbfcop/85Cfq16+flixZorffflvLli0LXtEAACCkhFS4Wbt2rfr16+d9PnbsWElSXl6e5s2bp5/97Gd69tlnNXXqVI0ePVoXXXSR3nzzTfXp0ydYJQMAgBATUuGmb9++OtX65pEjR2rkyJEBqggAAISbsFlzAwAAcDoINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWCblws2LFCg0aNEiZmZlyOBxauHDhCfe966675HA4NGvWrIDVBwAAQlvIhZuKigp16dJFc+fOPel+CxYs0KeffqrMzMwAVQYAAMJBZLAL+LGBAwdq4MCBJ91nz549uu+++/Tuu+/qhhtuCFBlAAAgHIRcuDkVj8ej4cOHa9y4cerYseMp96+qqlJVVZX3eVlZmSTJ7XbL7Xb7tLba8Xw9LnyHHoU+ehTa6E/os7VHZ/J9wi7cTJ8+XZGRkRo9evRp7T916lQ9/vjj9bYXFBQoPj7e1+VJkgoLC/0yLnyHHoU+ehTa6E/os61H5eXlp71vWIWbdevWafbs2Vq/fr0cDsdpvWf8+PEaO3as93lZWZlatWqlzp07KyEhwaf1ud1uFRYWKjc3V06n06djwzfoUeijR6GN/oQ+W3tUe+TldIRVuPnXv/6lAwcOKCsry7vN7XbrgQce0KxZs7Rz585673G5XHK5XPW2O51OvzXdn2PDN+hR6KNHoY3+hD7benQm3yWsws3w4cPVv3//OtsGDBig4cOH6/bbbw9SVQAAIJSEXLgpLy/X9u3bvc+LioqUn5+vlJQUZWVlKTU1tc7+UVFRysjI0EUXXRToUgEAQAgKuXCzdu1a9evXz/u8dr1MXl6e5s2bF6SqAABAuAi5cNO3b18ZY057/4bW2QAAgKYr5K5QDAAA0BiEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrhFy4WbFihQYNGqTMzEw5HA4tXLjQ+1p1dbUefvhh5ebmKi4uTpmZmfr1r3+tvXv3Bq9gAAAQUkIu3FRUVKhLly6aO3duvdcqKyu1fv16TZgwQevXr9dbb72lLVu2aPDgwUGoFAAAhKLIYBfwYwMHDtTAgQMbfC0xMVHvv/9+nW1//OMf1bNnT+3atUtZWVmBKBEAAISwkAs3Z6q0tFQOh0NJSUkNvl5VVaWqqirv87KyMkmS2+2W2+32aS214/l6XPgOPQp99Ci00Z/QZ2uPzuT7hHW4OXbsmB5++GHdeuutSkhIaHCfqVOn6vHHH6+3vaCgQPHx8X6pq7Cw0C/jwnfoUeijR6GN/oQ+23pUXl5+2vs6jDHGj7U0isPh0IIFCzRkyJB6r1VXV+sXv/iFvv76ay1btuyE4aahmZtWrVrp0KFDJ3zP2XK73SosLFRubq6cTqdPx4Zv0KPQR49CG/0Jfbb2qKysTCkpKSotLT3lz++wnLmprq7WzTffrK+++koffPDBSb+ky+WSy+Wqt93pdPqt6f4cG75Bj0IfPQpt9Cf02dajM/kuYRduaoPNtm3b9OGHHyo1NTXYJQEAgBAScuGmvLxc27dv9z4vKipSfn6+UlJS1KJFC/3yl7/U+vXr9c4778jtdqu4uFiSlJKSoujo6GCVDQAAQkTIhZu1a9eqX79+3udjx46VJOXl5emxxx7TokWLJEldu3at874PP/xQffv2DVSZAAAgRIVcuOnbt69OtsY5hNc/AwCAEBByVygGAABoDMINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKtEBrsAW+w+VKnfvVWgsiNHlLBhjeRwBLukM9a1VZIeuPaiYJcBAECjEG585Gi1W//a/s23T/Z/E9xiztK/th3ULy9pqezUuGCXAgDAWSPc+Eh6QoyeuqmzvvrqK2VnZysiIrxmbp56b6u+PnxUXx6sINwAAMIa4cZHEmOjNKRrpvJ1QF27ZsrpdAa7pDOy5LNifX34qHYerJA4MgUACGMsKIYkqXXat7M1Ow9WBLkSAAAah3ADSVLOd4eiir6pDHIlAAA0TqMPS61fv15PPvmkSkpKlJubq7vuuks5OTm+qA0BxMwNAMAWjZ65ueWWW/TTn/5UU6ZM0YUXXqif//zneu+993xRGwIo57tw8/XhSh2v8QS5GgAAzl6jZ26Sk5P1q1/9SpLUrVs3DRkyRFdddZWuvfbaRheHwDm3uUvNop2qPO7W7sOVanNOfLBLAgDgrDR65qZNmzZ66qmnZIyRJCUlJTV2SASBw+HwngLOoSkAQDhrdLipqqrSM888o6ysLF133XXq1KmT+vfvrz179viiPgRQTlozSVIR4QYAEMYafVjqrbfekiRVVFSosLBQBQUFKigo0K233qq9e/dq+/btjS4SgdG6dubmG8INACB8NTrcXHHFFbrzzjt1880367LLLtNll13mi7oQBN+fMcXp4ACA8NXow1Kvvfaadu/erT59+mjs2LHaunWrL+pCENSeMcVhKQBAOGt0uMnMzNT48eO1ZMkSdejQwXu2FMJP7WGpvaVHdazaHeRqAAA4O40+LJWQkKB27drpggsuUHx8vK6//nolJCT4ojYEWFp8tOJdkSqvqtHuQ5Vqm9482CUBAHDGGh1uXn75Zb3wwgtyOBwaPny4rrjiCl/UhSBwOBxqndZMn+0pU9HBCsINACAsNfqw1KBBg7RgwQJNmTJFixcv1tVXX61nnnnGF7UhCGqvdfMV95gCAISpRs/cXHfddaqoqFBFRYXcbrc8Ho/mzJmju+++2xf1IcC+v4Emi4oBAOGp0eFm7ty5SkxMVGJioqKionxRE4KIG2gCAMLdGR2W+uqrr/Tee++puLjYu61NmzZKS0tTVFSU9u7d6/MCEVi1Vykm3AAAwtVph5v58+frggsu0HXXXafzzz9fL730kiRp165dmjZtmnr27KmsrCy/FYrA+P508GOcDg4ACEunHW4mT56s++67T4WFhbrmmmt09913a8KECWrTpo3mzZunHj166I033vBnrQiAlLhoNY/59mgli4oBAOHotNfc7NixQ/fff7+ys7M1d+5cZWVl6aOPPlJBQYHat2/vzxoRQA6HQzlpcSr4ulRFByt0UQangwMAwstpz9xUV1crNjZWktSyZUvFxMToySefJNhYiBtoAgDC2RktKH7llVe0efNmSZLT6VRycrJfikJwccYUACCcnXa4ueKKKzRx4kR17NhRaWlpOnbsmGbPnq3XX39dX3zxhWpqavxZJwKo9owpbqAJAAhHp73mZvny5ZKkbdu2ad26dVq/fr3Wr1+vv/3tbyopKVF0dLQuvPBCFRQU+K1YBAaHpQAA4eyML+LXtm1btW3bVkOHDvVuKyoq0tq1a7VhwwafFofgyPnusNT+sipVHq9Rs+hGX+sRAICA8clPrZycHOXk5Oimm27yxXAIsqRm0UpqFqWSymrtPFipDpnc5R0AED4afeNM2IlDUwCAcEW4QYNqD02xqBgAEG4IN2iQd+aGcAMACDOEGzSode0NNDksBQAIM4QbNOj7w1LcXwoAEF4IN2hQ7VWKD5ZX6cix6iBXAwDA6eMCJmhQQkyUUuOi9U3FcT3/ryKd29x1VuPknpeoLq2SfFscAAAnQbjBCZ1/Tpy+qTiup5duO+sxYqIitO731yjOxR81AEBg8BMHJzRuQDv99eOdqvF4zur9H2w+oGPVHv37SBXhBgAQMPzEwQn1zElRz5yUs35/72kfaE/JUR2uPK7WivNhZQAAnBgLiuE3Sc2iJEklR1mQDAAIHMIN/Ca5WbQkqaTyeJArAQA0JYQb+E3idzM3hyuYuQEABA7hBn6TzGEpAEAQEG7gN0mxHJYCAAQe4QZ+411QXMnMDQAgcAg38JvaBcWHmbkBAAQQ4QZ+UztzU8qaGwBAABFu4DdJzNwAAIIg5MLNihUrNGjQIGVmZsrhcGjhwoV1XjfG6NFHH1WLFi0UGxur/v37a9u2s7/3EfzHu+aGU8EBAAEUcuGmoqJCXbp00dy5cxt8fcaMGXr66af17LPPatWqVYqLi9OAAQN07NixAFeKU6ldc3OkqkbV7rO7PxUAAGcq5O4tNXDgQA0cOLDB14wxmjVrln7/+9/rxhtvlCT97W9/U3p6uhYuXKihQ4cGslScQkLM93+8So9WKy3eFcRqAABNRciFm5MpKipScXGx+vfv792WmJioSy+9VJ988kmD4aaqqkpVVVXe52VlZZIkt9stt9vt0/pqx/P1uOHKoW8DTtmxGh0qP6bk2OD/caNHoY8ehTb6E/ps7dGZfJ/g/7Q5A8XFxZKk9PT0OtvT09O9r/3Y1KlT9fjjj9fbXlBQoPj4eN8XKamwsNAv44ajZk6jMklrNn6hI2nRwS7Hix6FPnoU2uhP6LOtR+Xl5ae9b1iFm7Mxfvx4jR071vu8rKxMrVq1UufOnZWQkODTz3K73SosLFRubq6cTqdPxw5X5378iYorSpV2Xmt1bX9usMuhR2GAHoU2+hP6bO1R7ZGX0xFW4SYjI0OStH//frVo0cK7ff/+/eratWuD73G5XHK56q/1cDqdfmu6P8cON8lx387WlB2rCanfE3oU+uhRaKM/oc+2Hp3Jdwm5s6VOJicnRxkZGVq6dKl3W1lZmVatWqVevXoFsTKcCLdgAAAEWsjN3JSXl2v79u3e50VFRcrPz1dKSoqysrI0ZswY/ed//qfatm2rnJwcTZgwQZmZmRoyZEjwisYJ1Z4OXnKUC/kBAAIj5MLN2rVr1a9fP+/z2vUyeXl5mjdvnh566CFVVFTozjvvVElJifr06aMlS5YoJiYmWCXjJGpnbg4zcwMACJCQCzd9+/aVMeaErzscDk2aNEmTJk0KYFU4W0mx391finADAAiQsFpzg/BTu6CY+0sBAAKFcAO/SozlsBQAILAIN/Cr2gXFpczcAAAChHADv2JBMQAg0Ag38Kuk72Zujla7dazarvucAABCE+EGfpUQEylnhEPSt3cGBwDA3wg38CuHw/GDRcWsuwEA+B/hBn7HLRgAAIFEuIHf1V7Ir4SZGwBAABBu4Hfe+0sxcwMACADCDfyu9owpTgcHAAQC4QZ+9/2aGw5LAQD8j3ADv0tmQTEAIIAIN/C7xGbcPBMAEDiEG/idd+aGi/gBAAKAcAO/S4qtPVuKmRsAgP8RbuB3XMQPABBIhBv4XXLc99e5McYEuRoAgO0IN/C72isUH3d7VHmcO4MDAPyLcAO/axbtVLTz2z9qLCoGAPgb4QZ+53A4lPjdupvDFSwqBgD4F+EGAVF7OngpMzcAAD8j3CAgkriQHwAgQAg3CIjaRcXcPBMA4G+EGwRE8nczN6XM3AAA/Ixwg4CovZAfMzcAAH8j3CAgatfccJViAIC/EW4QEN6bZ3JYCgDgZ4QbBEQSdwYHAAQI4QYBwangAIBAIdwgILgzOAAgUAg3CIhk74Li4/J4uDM4AMB/CDcIiMTvLuLnMdKRqpogVwMAsBnhBgERE+VUbJRTklTKoSkAgB8RbhAwyd4L+bGoGADgP4QbBEwiZ0wBAAKAcIOAqZ25KeVaNwAAPyLcIGC895eqYOYGAOA/hBsEjPf+UszcAAD8iHCDgEnmQn4AgAAg3CBgkmK/v5AfAAD+QrhBwHjX3DBzAwDwI8INAiapGTM3AAD/I9wgYLxrblhQDADwI8INAoZTwQEAgUC4QcDUHpYqO1YjN3cGBwD4CeEGAZP03Z3BJa5SDADwH8INAibSGaHmrkhJ3F8KAOA/kcEuAE1LUlyUjlTVaMOuEh097g7457vdbn15uFqRe0rldDr9+lmJsVE6LylWEREOv34OAKAuwg0CKrlZtHYfOqoH39gY3EL++UlAPiYmKkLnp8XrgnO/fWQmxerHWcfhkGKjIhXncqpZdKTiXZFqFu2Us4FQlJEQQ1gCgFMg3CCgftUzS9+Ub5fHBGdBsZFUffy4oqKj5c+IYIz0TUWVjlV79MW+Mn2xr8wn43bMTNArd1ymxB+sXwIA1EW4QUAN7ZmloT2zgvb5brdb+fn56tq1q98PS9W4Pdp9+Ki2HyjX9gPl2vHvcu0vO1ZvP48xOnrcrcrjblUcr1Fl1bf//fEJZTVujz7fW6YHXs/Xc8O7M4MDACdAuAH8JNIZoZy0OOWkxemaDumNHq/g6xL98tlP9M9NB/TM8h0a1e8CH1QJAPbhbCkgTHRumaTJN3aUJD353hat2PrvIFcEAKGJcAOEkVt6ZGloj1YyRhr96gbtPlQZ7JIAIOQQboAw89jgjurcMlElldW65+X1OlYd+FPqASCUseYGCDMxUU79adjFGjRnpQr3lGrs6/ka2KmFmsdEKiE2SgkxUYp3RdY75fyEHFKEwyGHvvuvQ4p3RSrSyb99AIQnwg0QhlomN9PTt3ZT3gur9Y/CYv2jsNin46cnuPT8r7urc8skn44LAIHAP82AMHVF23P0p2GX6JoO6bo0J0UdWiSoZXKsEmIi5YxwKMKh03o4Gpjh2V9Wpdv+e5UKvi4J+PcCgMZi5gYIY9d1ytB1nTJ8MpYxRh4jlVfV6P/NW6O1Xx3WsP9epb//v0vVpVWSTz4DAAKBmRsAkiSHwyFnhEOJsVGaN7KnerRO1pFjNbrtL6u0cXdJsMsDgNNGuAFQT7wrUi/eXjfg5BNwAIQJwg2ABsW7IjXv9p7q2TpFR47VaPh/r9LHOw4GuywAOKWwCzdut1sTJkxQTk6OYmNj1aZNG02ePFkmSDdiBGwW54rUi7f3UM+cFB2pqtHwv6zWCyuL+P8NQEgLu3Azffp0PfPMM/rjH/+oTZs2afr06ZoxY4bmzJkT7NIAK8W5IvW3kT31s27nye0xmvTOFxr7+kYuHgggZIXd2VIff/yxbrzxRt1www2SpNatW2v+/PlavXp1kCsD7BUT5dTMm7so97xETfnHJi3YsEdb9x/Rn4dfopbJzYJdHgDUEXbh5vLLL9dzzz2nrVu36sILL9TGjRu1cuVKzZw5s8H9q6qqVFVV5X1eVlYm6dvDW263b//lWTuer8eF79CjxsnrlaUL0+M0en6+Pt9bpkFzVuqa9umKjXaqWbRTsVFONXM5FRVx9pPCHuPR3j2VKji6UxGO78dJjY/WBefGKzulmaIjw27S2Rr8PxT6bO3RmXwfhwmzg+cej0e/+93vNGPGDDmdTrndbk2ZMkXjx49vcP/HHntMjz/+eL3ty5YtU3x8vL/LBaz070q3Znx0WF+W1AT8syMcUka8Uy2bRyqtmVOREVJkhMP7X6dDDV+ZMMRESIpySi6nQ9E/eIRB6cApxUQ6dFFqtE/HLC8vV9++fVVaWqqEhIST7ht24ebVV1/VuHHj9Ic//EEdO3ZUfn6+xowZo5kzZyovL6/e/g3N3LRq1UqHDh065W/OmXK73SosLFRubq6cTqdPx4Zv0CPfOVbt1jsF+3TgSJUqj7tVedyto8fdqqyukdtz9n+tGGNUWlqqxMREOb77Se8x0r6SY9rx73JVHLfrX6OAjS5Mj9fi0X18OmZZWZlSUlJOK9yE3WGpcePG6ZFHHtHQoUMlSbm5ufrqq680derUBsONy+WSy+Wqt93pdPrth5s/x4Zv0KPGi3M6dUvPbJ+P63a7lZ+fr65du9brkTFGxWXHtP1AubbtL9fB8irVeIyq3R5Vuz2qcRtVu8Pj32seY1RV820gPFbt0dFqt6pqPGFxJtrRo0cVGxsb7DJwEsHuUevUOJ//HXsm44VduKmsrFTEj47nO51OeTyeIFUEIFAcDodaJMaqRWKsrmh7TrDLaZJOFj4RGuhRGIabQYMGacqUKcrKylLHjh21YcMGzZw5UyNHjgx2aQAAIASEXbiZM2eOJkyYoHvuuUcHDhxQZmamfvOb3+jRRx8NdmkAACAEhF24ad68uWbNmqVZs2YFuxQAABCCuFgEAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABgFcINAACwCuEGAABYhXADAACsQrgBAABWIdwAAACrEG4AAIBVCDcAAMAqhBsAAGAVwg0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCphGW727Nmj2267TampqYqNjVVubq7Wrl0b7LIAAEAIiAx2AWfq8OHD6t27t/r166fFixfrnHPO0bZt25ScnBzs0gAAQAgIu3Azffp0tWrVSi+++KJ3W05OThArAgAAoSTsws2iRYs0YMAA3XTTTVq+fLnOO+883XPPPbrjjjsa3L+qqkpVVVXe56WlpZK+nQFyu90+rc3tdqu8vFyHDx+W0+n06djwDXoU+uhRaKM/oc/WHpWVlUmSjDGn3tmEGZfLZVwulxk/frxZv369+fOf/2xiYmLMvHnzGtx/4sSJRhIPHjx48ODBw4LH7t27T5kVHMacTgQKHdHR0erevbs+/vhj77bRo0drzZo1+uSTT+rt/+OZG4/Ho0OHDik1NVUOh8OntZWVlalVq1bavXu3EhISfDo2fIMehT56FNroT+iztUfGGB05ckSZmZmKiDj5+VBhd1iqRYsW6tChQ51t7du315tvvtng/i6XSy6Xq862pKQkf5UnSUpISLDqD5SN6FHoo0ehjf6EPht7lJiYeFr7hd2p4L1799aWLVvqbNu6dauys7ODVBEAAAglYRdufvvb3+rTTz/VE088oe3bt+uVV17Rc889p1GjRgW7NAAAEALCLtz06NFDCxYs0Pz589WpUydNnjxZs2bN0rBhw4JdmlwulyZOnFjvMBhCBz0KffQotNGf0EePpLBbUAwAAHAyYTdzAwAAcDKEGwAAYBXCDQAAsArhBgAAWIVw4yNz585V69atFRMTo0svvVSrV68OdklN1tSpU9WjRw81b95c5557roYMGVLv2kjHjh3TqFGjlJqaqvj4eP3iF7/Q/v37g1Qxpk2bJofDoTFjxni30aPg27Nnj2677TalpqYqNjZWubm5Wrt2rfd1Y4weffRRtWjRQrGxserfv7+2bdsWxIqbFrfbrQkTJignJ0exsbFq06aNJk+eXOfeS022R76651NT9uqrr5ro6GjzwgsvmM8//9zccccdJikpyezfvz/YpTVJAwYMMC+++KL57LPPTH5+vrn++utNVlaWKS8v9+5z1113mVatWpmlS5eatWvXmssuu8xcfvnlQay66Vq9erVp3bq16dy5s7n//vu92+lRcB06dMhkZ2ebESNGmFWrVpkvv/zSvPvuu2b79u3efaZNm2YSExPNwoULzcaNG83gwYNNTk6OOXr0aBArbzqmTJliUlNTzTvvvGOKiorMG2+8YeLj483s2bO9+zTVHhFufKBnz55m1KhR3udut9tkZmaaqVOnBrEq1Dpw4ICRZJYvX26MMaakpMRERUWZN954w7vPpk2bjCTzySefBKvMJunIkSOmbdu25v333zdXXnmlN9zQo+B7+OGHTZ8+fU74usfjMRkZGeYPf/iDd1tJSYlxuVxm/vz5gSixybvhhhvMyJEj62z7+c9/boYNG2aMado94rBUIx0/flzr1q1T//79vdsiIiLUv3//Bm/kicArLS2VJKWkpEiS1q1bp+rq6jo9a9eunbKysuhZgI0aNUo33HBDnV5I9CgULFq0SN27d9dNN92kc889V926ddPzzz/vfb2oqEjFxcV1epSYmKhLL72UHgXI5ZdfrqVLl2rr1q2SpI0bN2rlypUaOHCgpKbdo7C7cWaoOXjwoNxut9LT0+tsT09P1+bNm4NUFWp5PB6NGTNGvXv3VqdOnSRJxcXFio6OrncD1fT0dBUXFwehyqbp1Vdf1fr167VmzZp6r9Gj4Pvyyy/1zDPPaOzYsfrd736nNWvWaPTo0YqOjlZeXp63Dw393UePAuORRx5RWVmZ2rVrJ6fTKbfbrSlTpniv2N+Ue0S4gdVGjRqlzz77TCtXrgx2KfiB3bt36/7779f777+vmJiYYJeDBng8HnXv3l1PPPGEJKlbt2767LPP9OyzzyovLy/I1UGSXn/9db388st65ZVX1LFjR+Xn52vMmDHKzMxs8j3isFQjpaWlyel01juLY//+/crIyAhSVZCke++9V++8844+/PBDtWzZ0rs9IyNDx48fV0lJSZ396VngrFu3TgcOHNDFF1+syMhIRUZGavny5Xr66acVGRmp9PR0ehRkLVq0UIcOHepsa9++vXbt2iVJ3j7wd1/wjBs3To888oiGDh2q3NxcDR8+XL/97W81depUSU27R4SbRoqOjtYll1yipUuXerd5PB4tXbpUvXr1CmJlTZcxRvfee68WLFigDz74QDk5OXVev+SSSxQVFVWnZ1u2bNGuXbvoWYBcffXVKiwsVH5+vvfRvXt3DRs2zPtrehRcvXv3rncJha1btyo7O1uSlJOTo4yMjDo9Kisr06pVq+hRgFRWVioiou6PcafTKY/HI6mJ9yjYK5pt8OqrrxqXy2XmzZtnvvjiC3PnnXeapKQkU1xcHOzSmqS7777bJCYmmmXLlpl9+/Z5H5WVld597rrrLpOVlWU++OADs3btWtOrVy/Tq1evIFaNH54tZQw9CrbVq1ebyMhIM2XKFLNt2zbz8ssvm2bNmpm///3v3n2mTZtmkpKSzP/+7/+agoICc+ONNzaJ04xDRV5enjnvvPO8p4K/9dZbJi0tzTz00EPefZpqjwg3PjJnzhyTlZVloqOjTc+ePc2nn34a7JKaLEkNPl588UXvPkePHjX33HOPSU5ONs2aNTM/+9nPzL59+4JXNOqFG3oUfG+//bbp1KmTcblcpl27dua5556r87rH4zETJkww6enpxuVymauvvtps2bIlSNU2PWVlZeb+++83WVlZJiYmxpx//vnmP/7jP0xVVZV3n6baI4cxP7iUIQAAQJhjzQ0AALAK4QYAAFiFcAMAAKxCuAEAAFYh3AAAAKsQbgAAgFUINwAAwCqEGwABt3PnTjkcDuXn5we7FK/NmzfrsssuU0xMjLp27dro8R577LEzHsfhcGjhwoWN/mygqSPcAE3QiBEj5HA4NG3atDrbFy5cKIfDEaSqgmvixImKi4vTli1b6tyL54f69u2rMWPGnNZ4Dz744AnHAeBfhBugiYqJidH06dN1+PDhYJfiM8ePHz/r9+7YsUN9+vRRdna2UlNTz3ocY4xqamoUHx/fqHEAnD3CDdBE9e/fXxkZGZo6deoJ92no0MqsWbPUunVr7/MRI0ZoyJAheuKJJ5Senq6kpCRNmjRJNTU1GjdunFJSUtSyZUu9+OKL9cbfvHmzLr/8csXExKhTp05avnx5ndc/++wzDRw4UPHx8UpPT9fw4cN18OBB7+t9+/bVvffeqzFjxigtLU0DBgxo8Ht4PB5NmjRJLVu2lMvlUteuXbVkyRLv6w6HQ+vWrdOkSZPkcDj02GOP1RtjxIgRWr58uWbPni2HwyGHw6GdO3dq2bJlcjgcWrx4sS655BK5XC6tXLmy3u/dmjVrdM011ygtLU2JiYm68sortX79+hP+3h8/flz33nuvWrRooZiYGGVnZ5+0VwC+R7gBmiin06knnnhCc+bM0ddff92osT744APt3btXK1as0MyZMzVx4kT99Kc/VXJyslatWqW77rpLv/nNb+p9zrhx4/TAAw9ow4YN6tWrlwYNGqRvvvlGklRSUqKrrrpK3bp109q1a7VkyRLt379fN998c50x/vrXvyo6OlofffSRnn322Qbrmz17tp566ik9+eSTKigo0IABAzR48GBt27ZNkrRv3z517NhRDzzwgPbt26cHH3ywwTF69eqlO+64Q/v27dO+ffvUqlUr7+uPPPKIpk2bpk2bNqlz58713n/kyBHl5eVp5cqV+vTTT9W2bVtdf/31OnLkSIM1P/3001q0aJFef/11bdmyRS+//HKdUAngJIJ8404AQZCXl2duvPFGY4wxl112mRk5cqQxxpgFCxaYH/61MHHiRNOlS5c67/2v//ovk52dXWes7Oxs43a7vdsuuugic8UVV3if19TUmLi4ODN//nxjjDFFRUVGkpk2bZp3n+rqatOyZUszffp0Y4wxkydPNtdee22dz969e7eR5L2r8ZVXXmm6det2yu+bmZlppkyZUmdbjx49zD333ON93qVLFzNx4sSTjvPjO5cbY8yHH35oJJmFCxfW2d7Q790Pud1u07x5c/P22297t0kyCxYsMMYYc99995mrrrrKeDyek9YEoD5mboAmbvr06frrX/+qTZs2nfUYHTt2VETE93+dpKenKzc31/vc6XQqNTVVBw4cqPO+Xr16eX8dGRmp7t27e+vYuHGjPvzwQ8XHx3sf7dq1k/Tt+phal1xyyUlrKysr0969e9W7d+8623v37t2o7/xj3bt3P+nr+/fv1x133KG2bdsqMTFRCQkJKi8v165duxrcf8SIEcrPz9dFF12k0aNH67333vNZrYDtIoNdAIDg+slPfqIBAwZo/PjxGjFiRJ3XIiIiZIyps626urreGFFRUXWeOxyOBrd5PJ7Trqu8vFyDBg3S9OnT673WokUL76/j4uJOe0x/OlUdeXl5+uabbzR79mxlZ2fL5XKpV69eJ1wEffHFF6uoqEiLFy/WP//5T918883q37+//ud//scf5QNWYeYGgKZNm6a3335bn3zySZ3t55xzjoqLi+sEHF9em+bTTz/1/rqmpkbr1q1T+/btJX37w/3zzz9X69atdcEFF9R5nEmgSUhIUGZmpj766KM62z/66CN16NDhjOqNjo6W2+0+o/f88PNGjx6t66+/Xh07dpTL5aqzOLohCQkJuuWWW/T888/rtdde05tvvqlDhw6d1ecDTQnhBoByc3M1bNgwPf3003W29+3bV//+9781Y8YM7dixQ3PnztXixYt99rlz587VggULtHnzZo0aNUqHDx/WyJEjJUmjRo3SoUOHdOutt2rNmjXasWOH3n33Xd1+++1nHDDGjRun6dOn67XXXtOWLVv0yCOPKD8/X/fff/8ZjdO6dWutWrVKO3fu1MGDB89oJqpt27Z66aWXtGnTJq1atUrDhg1TbGzsCfefOXOm5s+fr82bN2vr1q164403lJGRoaSkpDOqGWiKCDcAJEmTJk2q98O6ffv2+tOf/qS5c+eqS5cuWr16dYNnEp2tadOmadq0aerSpYtWrlypRYsWKS0tTZK8sy1ut1vXXnutcnNzNWbMGCUlJdVZ33M6Ro8erbFjx+qBBx5Qbm6ulixZokWLFqlt27ZnNM6DDz4op9OpDh066JxzzjnhepmG/OUvf9Hhw4d18cUXa/jw4Ro9erTOPffcE+7fvHlzzZgxQ927d1ePHj20c+dO/eMf/zjj7w40RQ7z4wPqAAAAYYx/AgAAAKsQbgAAgFUINwAAwCqEGwAAYBXCDQAAsArhBgAAWIVwAwAArEK4AQAAViHcAAAAqxBuAACAVQg3AADAKoQbAABglf8Poh/vg6wGJjsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Rwp plot\n",
    "def rwp_plot():\n",
    "    minvalues = [df.iloc[0]['Rwp']]\n",
    "    for i in range(1, df.shape[0]):\n",
    "        minvalues.append(min(minvalues[-1], df.iloc[i]['Rwp']))\n",
    "    minvalues = pd.DataFrame(minvalues)\n",
    "    \n",
    "    minvalues.plot(legend=None)\n",
    "    plt.ylim([6, 16])\n",
    "    plt.grid(color='#cccccc')\n",
    "    plt.ylabel('$R_{wp}$')\n",
    "    plt.xlabel('Number of trials')\n",
    "    plt.show()\n",
    "    \n",
    "rwp_plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk8AAAG0CAYAAAAresMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAACXw0lEQVR4nOzdeXgT1dcH8O8kTdJ9X9JCN2SHQoFCKQiCIBUQZBEBEQqCKFIUqiCoLKICorIoCD9X9BUUUEEERBEQQcpWWik7aqEsTcvShe5tct8/hkwzbbamaZO25/M8PKQzd+7cOZlkTu7cmeEYYwyEEEIIIcQsEls3gBBCCCGkPqHkiRBCCCGkGih5IoQQQgipBkqeCCGEEEKqgZInQgghhJBqoOSJEEIIIaQaKHkihBBCCKkGB1s3oCHSaDS4efMm3NzcwHGcrZtDCCGEEDMwxnDv3j0EBQVBIjHcv0TJUy24efMmgoODbd0MQgghhFjg2rVraNq0qcH5lDzVAjc3NwB88N3d3S2uR61W4/Tp0+jQoQOkUqm1mtfgUJxMoxiZh+JkGsXIPBQn0+wxRnl5eQgODhaO44ZQ8lQLtKfq3N3da5w8ubq6wt3d3W52LHtEcTKNYmQeipNpFCPzUJxMs+cYmRpyQwPGCSGEEEKqgZInQgghhJBqoOSJEEIIIaQaaMwTIYSQaispKTF6KXdjp1arAQDFxcV2N57HXtgiRjKZzCrrouSJEEKI2UpLS8FxHK5cuUL3sTOCMQYHBwdcvXqV4mSArWLk6ekJpVJZo3VS8kQIIcQsjDFkZmbCyckJISEh1KNiBGMMRUVFcHJyouTJgLqOEWMMhYWFyMrKAgAEBgZaXBclT4QQQsxSXl6OwsJC+Pj4wNnZmZICIxhj0Gg0cHR0pDgZYIsYOTk5AQCysrLg7+9v8Q8AOmFNCCHELNoxKjKZzMYtIcRyzs7OAICysjKL66DkiRBCCCGNhjV6uSh5IiYlp2fjwMUsWzeDEEIIsQs05omYNPzjIwCAv+Y+jCaeTjZuDSGEEGJb1PNEzKbKLbZ1EwghpN7r3bs3Nm3aVO3ltLeHSElJsWi9t2/fhr+/P65fv27R8qQCJU/EbHTBCCGkvlGr1ejRowdGjBghmp6bm4vg4GC8/vrrACoSE+0/b29vPPTQQzh06JBouUWLFgllHBwc4Ovri969e2PVqlUoKSkx2Z4dO3YgMzMTY8aMEU0/cuQIBg0aBC8vLzg6OiIiIgIrVqwQBulbg6+vLyZMmICFCxdarc7GipInYhRjTHgtoeyJEGIl2gSktkmlUmzYsAF79uzBxo0bhekzZsyAt7d3lUTi999/R0ZGBv78808EBQXhscceQ2ZmpqhMu3btkJGRgfT0dBw4cACjRo3C0qVL0aNHD9y7d89oez788ENMmjRJdHf2bdu24aGHHkLTpk1x4MABXLhwAS+99BLefvttjBkzRvQ9XFOTJk3Cxo0bcffuXavV2RhR8kSM0v3MSil5IoTUQy1btsSyZcswY8YMZGRk4KeffsJ3332Hr7/+GnK5XFTWx8cHSqUS7du3x2uvvYa8vDwcO3ZMVMbBwQFKpRJBQUGIiIjAjBkzcPDgQZw5cwbvvvuuwXbcunUL+/fvx5AhQ4RpBQUFePbZZzF06FB88skniIyMRFhYGKZMmYKvvvoK33//PbZs2SKq58KFC+jRowccHR3Rvn17HDx4UJiXnZ2NcePGwc/PD05OTmjRogW+/PJLYX67du0QFBSEbdu2WRRLwqPkiRil1smeKHcihNRXM2bMQMeOHTF+/HhMnToVCxYsQMeOHQ2WLyoqwtdffw0AVRIsfVq3bo2BAwfixx9/NFjm8OHDcHZ2Rps2bYRpv/32G+7cuYNXXnmlSvkhQ4agZcuW+Pbbb0XTZ8+ejZdffhnJycmIiYnBkCFDcOfOHQDA/Pnzce7cOfzyyy84f/481q1bB19fX9Hy3bp1q3I6klQPXW1HjNLQaTtCiJXoO02nO82ap6f0rWfdunVo06YNIiIiMHfuXL3levToAYlEgsLCQjDG0KVLF/Tr18+sdbRu3Rq//fabwflXr15FQECA6JTdpUuXAECUUFWuU1tGKz4+HiNHjgQArFu3Dnv27MHnn3+OOXPmID09HZ06dUJUVBQAICwsrEqdQUFBSE5ONmubiH7U80SM0mgqXtMD1Akh9dkXX3wBZ2dnpKWlGbzibPPmzUhOTsYPP/yA5s2bY8OGDWbfUZ0xZnQcV1FRERwdHQ0ua66YmBjhtYODA6KionD+/HkAwLRp0/Ddd98hMjISc+bMwZEjR6os7+TkhMLCQrPXR6qiwyHRSzuY08XVVZjWISKizgZ5EkIaHsaY8M/YtNpw5MgRrFy5Ejt37kS3bt0wefJkvesMDg5GixYtMHz4cCxZsgTDhw836yo6ADh//jzCw8MNzvf19UV2drZoWsuWLYVlDdWpLWOOgQMH4urVq5g1axZu3ryJfv36VTklePfuXfj5+ZldJ6mKkidiHFexizCmMVKQEELsU2FhISZOnIhp06ahb9+++Pzzz3H8+HGsX7/e6HJPPPEEHBwc8PHHH5tcx4ULF7Bnzx7hdJo+nTp1gkqlEiVQAwYMgLe3Nz744IMq5Xfs2IHLly9j7NixoulHjx4VXpeXlyMpKUl02s/Pzw9xcXH45ptvsGrVKnzyySei5c+cOYNOnTqZ3CZiGCVPRC/tL8Hs7Bxh2oXz5+vkFyIhhFjTvHnzwBjDsmXLAPDjgN5//33MmTMHV65cMbgcx3F48cUXsWzZMtFprvLycqhUKty8eROpqan46KOP8NBDDyEyMhKzZ882WF+nTp3g6+uLv/76S5jm4uKC//3vf/jpp58wdepUnD59GleuXMHnn3+OiRMn4oknnsCTTz4pqmft2rXYtm0bLly4gOnTpyM7OxvPPPMMAGDBggX46aef8M8//+Ds2bPYuXOnKLEqLCxEUlISBgwYUK0YEjG7SZ5u3LiBp59+Gj4+PnByckJERAROnjwpzGeMYcGCBQgMDISTkxP69++Py5cvi+q4e/cuxo0bB3d3d3h6emLy5MnIz88XlTl9+jR69eoFR0dHBAcHY/ny5VXasnXrVrRu3Vq4Udnu3btrZ6PrATUNGCeE1IK6+iF28OBBrF27Fl9++SWcnZ2F6c899xx69Ohh8PSdVlxcHMrKyrBmzRph2tmzZxEYGIiQkBD06dMHW7Zswbx583Do0CG46gx1qEwqlQr3WdL1xBNP4MCBA0hPT0evXr3QqlUrrFy5Eq+//jq+++67KkMlli1bhmXLlqFjx444fPgwduzYIVxRJ5fLMW/ePHTo0AG9e/eGVCrFd999Jyz7008/ISQkBL169TIvgEQ/Zgfu3r3LQkND2cSJE9mxY8fYf//9x3799Vf2zz//CGWWLVvGPDw82Pbt29nff//Nhg4dysLDw1lRUZFQ5tFHH2UdO3ZkR48eZYcOHWLNmzdnY8eOFebn5uaygIAANm7cOHbmzBn27bffMicnJ/a///1PKPPXX38xqVTKli9fzs6dO8feeOMNJpPJWGpqqtnbk5ubywCw3NzcGsWlvLycnTx5kpWXl9eonpq4da+Yhb66k4W+upNduZ1vs3YYYw9xsncUI/NQnIwrKipiZ8+eZXfu3GEajcbWzbFrGo2G5efnV4lTRkYG8/b2ZleuXLFJu6Kjo9nGjRttsu7KDMWothUVFbFz586J8gctc4/fdnGrgnfffRfBwcGiG3npDrpjjGHVqlV444038PjjjwMAvv76awQEBGD79u0YM2YMzp8/jz179uDEiRPCJZofffQRBg0ahPfffx9BQUHYuHEjSktL8cUXX0Aul6Ndu3ZISUnBihUrMHXqVADA6tWr8eijjwpdr2+99Rb27t2LNWvWGDw/XlJSIhpQmJeXB4B/LEBNbq2vXdaat+evrrKycuG1Wq2xaVsMsYc42TuKkXkoTsbpxoXR6XujtPGpHKeAgAB89tlnuHr1KkJCQuq0Tbdv38bw4cOtftdySxmKUV2slzGm9xht7mffLpKnHTt2IDY2FqNGjcLBgwfRpEkTvPDCC3j22WcBAGlpaVCpVOjfv7+wjIeHB6Kjo5GYmIgxY8YgMTERnp6eQuIEAP3794dEIsGxY8cwfPhwJCYmonfv3qIbnsXGxuLdd99FdnY2vLy8kJiYiISEBFH7YmNjsX37doPtX7p0Kd58880q00+fPm20C9dcqampNa7DUneKKnakc+fPIdvVLnYZvWwZp/qCYmQeipNhDg78d0BRUZGNW1I/6IuTdrxRXd8uwNnZGfHx8Xb33tV1e0pKSlBWVoYLFy5UmVd5qI8hdnEk/O+//7Bu3TokJCTgtddew4kTJ/Diiy9CLpcjLi4OKpUKAJ+x6woICBDmqVQq+Pv7i+Y7ODjA29tbVKbyZaTaOlUqFby8vKBSqYyuR5958+aJEq68vDwEBwejQ4cOcHd3r04oRNRqNVJTUxEREQGpVGpxPTVxM6cI2Mnf+r9NmzYI83GxSTuMsYc42TuKkXkoTsYVFxfj6tWrAPh7BdFtSwxjjKGoqIjiZIStYiSRSCCTydC8efMq993SnjkyxS6SJ41Gg6ioKCxZsgQAf0XCmTNnsH79esTFxdm4daYpFAooFIoq06VSqVW+gK1VT3Vod2Spuz+aTvsCANC6VWuU52QAsM8ue1vEqb6hGJmH4qSfbkzonm/moTiZVtcx0q5P3+fc3M+9XVxtFxgYiLZt24qmtWnTBunp6QAApVIJAFWebJ2ZmSnMUyqVyMrKEs0vLy/H3bt3RWX01aG7DkNltPMbG9EOTV8AhBBCiH0kTz179sTFixdF0y5duoTQ0FAA/OBxpVKJffv2CfO1T7rW3qY+JiYGOTk5SEpKEsrs378fGo0G0dHRQpk///wTZWVlQpm9e/eiVatW8PLyEsrorkdbRvd2+I2BdkDdpX/+FaZdunSJ7vNECCGk0bOL5GnWrFk4evQolixZgn/++QebNm3CJ598gunTpwPgez9mzpyJt99+Gzt27EBqaiomTJiAoKAgDBs2DADfU/Xoo4/i2WefxfHjx/HXX38hPj4eY8aMQVBQEADgqaeeglwux+TJk3H27Fls3rwZq1evFo1Xeumll7Bnzx588MEHuHDhAhYtWoSTJ08iPj6+zuNiDzSVHqNACCGENHZ2Meapa9eu2LZtG+bNm4fFixcjPDwcq1atwrhx44Qyc+bMQUFBAaZOnYqcnBw8+OCD2LNnj2iw18aNGxEfH49+/fpBIpFg5MiR+PDDD4X5Hh4e+O233zB9+nR06dIFvr6+WLBggXCbAoB/ovamTZvwxhtv4LXXXkOLFi2wfft2tG/fvm6CYWc0GkqYCCGEEF12kTwBwGOPPYbHHnvM4HyO47B48WIsXrzYYBlvb29s2rTJ6Ho6dOiAQ4cOGS0zatQojBo1yniDGwnKnQghhBAxuzhtR+yXWid7ojyKEEJqrnfv3iZ/6Fti4sSJwlAWU27fvg1/f39cv37d6u1oDCh5IkZpaJwTIaQeU6vV6NGjB0aMGCGanpubi+DgYLz++usAgCtXrgiXsHMcB29vbzz00ENVzlQsWrRIKOPg4ABfX1/07t0bq1atEj1pwpAdO3YgMzMTY8aMEU1PTk7GqFGjEBAQAEdHR7Ro0QLPPvssLl26VMMI6Ofr64sJEyZg4cKFtVJ/Q0fJEzFKPGDchg0hhBALSKVSbNiwAXv27BE9kHfGjBnw9vaukjz8/vvvyMjIwJ9//omgoCA89thjVW5f065dO2RkZCA9PR0HDhzAqFGjsHTpUvTo0QP37t0z2p4PP/wQkyZNgkRScfjduXMnunfvjpKSEmzcuBHnz5/HN998Aw8PD8yfP98KUdBP+5Diu3fv1to6GipKnohRNOaJEGIIYwyFpeU2+Vedq39btmyJZcuWYcaMGcjIyMBPP/2E7777Dl9//bXocV0A4OPjA6VSifbt2+O1114Tboujy8HBAUqlEkFBQYiIiMCMGTNw8OBBnDlzBu+++67Bdty6dQv79+/HkCFDhGmFhYWYNGkSBg0ahB07dqB///4IDw9HdHQ03n//ffzvf/8DwPegTZ48GeHh4XByckKrVq2wevVqo9ut0WiwfPlyNG/eHAqFAiEhIXjnnXeE+e3atUNQUBC2bdtmdiwJz24GjBP7JD5tR5kUIaRCUZkabRf8apN1n1scC2e5+YewGTNmYNu2bRg/fjxSU1OxYMECdOzY0WD5oqIifP311wBQJcHSp3Xr1hg4cCB+/PFHvP3223rLHD58GM7OzmjTpo0w7ddff8Xt27cxZ84cvct4enoC4BOhpk2bYuvWrfDx8cGRI0cwdepUBAYG4sknn9S77Lx58/Dpp59i5cqVePDBB5GRkVHleW7dunXDoUOHMHnyZJPbSCpQ8kSMolsVEEIaAo7jsG7dOrRp0wYRERGYO3eu3nI9evSARCJBYWEhGGPo0qUL+vXrZ9Y6Wrdujd9++83g/KtXryIgIEB0yu7y5cvCssbIZDLRA+jDw8ORmJiILVu26E2e7t27h9WrV2PNmjXCY84eeOABPPjgg6JyQUFBSE5ONr1xRISSJ2IU5U6EEEOcZFKcWxxrs3VX1xdffAFnZ2ekpaXh+vXrCAsLq1Jm8+bNaN26Nc6cOYM5c+Zgw4YNkMlkZtXPGDP6jLaioqIqD6KtzunHtWvX4osvvkB6ejqKiopQWlqKyMhIvWXPnz+PkpISk4mfk5MTCgsLzW4D4VHyRIwS3aqAEilCiA6O46p16syWjhw5gpUrV+K3337D22+/jcmTJ+P333+vkuwEBwejRYsWaNGiBcrLyzF8+HCcOXNG78PfKzt//jzCw8MNzvf19UV2drZoWsuWLQEAFy5cMPoYsO+++w6vvPIKPvjgA8TExMDNzQ3vvfdelfFYWk5OTibbCwB3796Fn5+fWWVJBRowToyiR7IQQuq7wsJCTJw4EdOmTUPfvn3x+eef4/jx41i/fr3R5Z544gk4ODjg448/NrmOCxcuYM+ePRg5cqTBMp06dYJKpRIlUAMGDICvry+WL1+ud5mcnBwAwF9//YUePXrghRdeQKdOndC8eXP8+++/epcBgBYtWsDJyanKs1orO3PmDDp16mS0DKmKkidilJrRTTIJIfXbvHnzwBjDsmXLAABhYWF4//33MWfOHFy5csXgchzH4cUXX8SyZctEp7bKy8uhUqlw8+ZNpKam4qOPPsJDDz2EyMhIzJ4922B9nTp1gq+vL/766y9hmouLCz777DPs2rULQ4cOxe+//44rV67g5MmTmDNnDp5//nkAfDJ08uRJ/Prrr7h06RLmz5+PEydOGFyXo6MjXn31VcyZMwdff/01/v33Xxw9ehSff/65UKawsBBJSUkYMGCAyRgSMUqeiFE05okQUp8dPHgQa9euxZdffglnZ2dh+nPPPYcePXpg8uTJRnvY4+LiUFZWhjVr1gjTzp49i8DAQISEhKBPnz7YsmUL5s2bh0OHDsHV1dVgXVKpVLi3kq7HH38cR44cgUwmw1NPPYXWrVtj7NixyM3NFa7ce+655zBixAiMHj0a0dHRuHPnDl544QWj2z5//ny8/PLLWLBgAdq0aYPRo0cjKytLmP/TTz8hJCQEvXr1MloPqYpjdF7G6vLy8uDh4YHc3Fy4u7tbXI9arUZKSgoiIyMhlVZ/cKQ1HLiQhUkb+F83v83qjZYBbjZphzH2ECd7RzEyD8XJuOLiYvz3339QKpXw8vIyOji6sWOMobCwEM7OzqI4qVQqtGvXDqdOnUJoaKgNWwh0794dL774Ip566imbrN9QjGpbcXEx0tLSEB4eXmUAv7nHb+p5IkbRHcYJIcR6lEolPv/8c6Snp9u0Hbdv38aIESMwduxYm7ajvqofl0kQm1HTeTtCCLEqcx/eW5t8fX0N3piTmEY9T8Qo3dyJ0ZBxQgghhJInYhwNiSOEEELEKHkiRqlpzBMhhBAiQskTMYoSJkIIIUSMkidCCCGEkGqg5IkYpdvxRL1QhBBCCCVPhBBCiFWEhYVh1apVtm6GVW3YsAGenp4my33++ee18piXP/74AxzHCc/4M2XMmDH44IMPrN6Oyih5IkYx0bPtqOuJEFL/TJw4ERzHCf98fHzw6KOP4vTp07ZuWoNQXFyM+fPnY+HChaLpeXl5eP3119G6dWs4OjpCqVSif//++PHHH2vtSu433ngD77zzDnJzc2ulfi1KngghhDR4jz76KDIyMpCRkYF9+/bBwcEBjz32mK2bZVJpaamtm2DS999/D3d3d/Ts2VOYlpOTgx49euDrr7/GvHnzcOrUKfz5558YPXo05syZU2vJTfv27fHAAw/gm2++qZX6tSh5IoQQ0uApFAoolUoolUpERkZi7ty5uHbtGm7duiWUefXVV9GyZUs4OzujWbNmmD9/PsrKykT1/Pzzz+jatSscHR3h6+uL4cOHG1znhg0b4OXlhX379gEA7t27h3HjxsHFxQWBgYFYuXIl+vTpg5kzZwrLhIWF4a233sKECRPg7u6OqVOnAgB++OEHtGvXDgqFAmFhYVVOTXEch+3bt4umeXp6YsOGDQCAK1eugOM4/Pjjj+jbty+cnZ3RsWNHJCYmVmlzSEgInJ2dMXz4cNy5c8dkbL/77jsMGTJENO21117DlStXcOzYMcTFxaFt27Zo2bIlnn32WaSkpAgPUP7222/RtWtXuLm5QalU4qmnnhI9vFifv/76C3369IGzszO8vLwQGxuL7OxsYf6QIUPw3XffmWx3TVDyRMxGA8YJISKMAQUFtvlXgy+k/Px8fPPNN2jevDl8fHyE6W5ubtiwYQPOnTuH1atX49NPP8XKlSuF+bt27cLw4cMxaNAgJCcnY9++fejWrZvedSxfvhwLFizAr7/+in79+gEAEhIS8Ndff2HHjh3Yu3cvDh06hFOnTlVZ9v3330fHjh2RnJyM+fPnIykpCU8++STGjBmD1NRULFq0CPPnzxcSo+p4/fXX8corryAlJQUtW7bE2LFjUV5eDgA4duwYJk+ejPj4eKSkpKBv3754++23TdZ5+PBhREVFCX9rNBp89913GDduHIKCgqqUd3V1hYMD/3S4srIyLF68GH///Te2b9+OK1euYOLEiQbXlZKSgn79+qFt27ZITEzE4cOHMWTIEKjVaqFMt27dcPz4cZSUlJgblupjxOpyc3MZAJabm1ujesrLy9nJkydZeXm5lVpWfdtOXWehr+5koa/uZKnXc2zWDmPsIU72jmJkHoqTcUVFRezs2bPszp07TKPRMJafzxifxtT9v/x8s9sdFxfHpFIpc3FxYS4uLgwACwwMZElJSUaXe++991iXLl2Ev2NiYti4ceMMlg8NDWUrV65kc+bMYYGBgez48eN8nBhjeXl5TCaTsa1btwrlc3JymLOzM3vppZdEdQwbNkxU71NPPcUeeeQR0bTZs2eztm3bCn8DYNu2bROV8fDwYF9++SVjjLG0tDQGgH322WfC/LNnzzIA7Pz584wxxsaOHcsGDRokqmP06NHMw8PD4DZnZ2czAOzPP/8UpmVmZjIAbMWKFQaXY4wxjUbD8vPzhRgxxtiJEycYAHbv3j3GGGMHDhxgAFh2drbQxp49exqt9++//2YA2JUrV/TOLyoqYufOnWNFRUVV5pl7/KaeJ0IIIQ1e3759kZKSgpSUFBw/fhyxsbEYOHAgrl69KpTZvHkzevbsCaVSCVdXV7zxxhtIT08X5mt7PYz54IMP8Omnn+LQoUNo27atMP2///5DWVmZqKfKw8MDrVq1qlKHbi8OAJw/f140nggAevbsicuXL4t6XMzRoUMH4XVgYCAACKfJzp8/j+joaFH5mJgYo/UVFRUBABwdHYVprBq9gsnJyRg6dChCQkLg5uaGhx56CABEcddlznvg5OQEACgsLDS7HdVFyRMxiq6wI4QY5OwM5Ofb5p+zc7Wa6uLigubNm6N58+bo2rUrPvvsMxQUFODTTz8FACQmJmLcuHEYNGgQdu7cieTkZLz++uuiAdvag7IxvXr1glqtxpYtW6oXy0ptrS6O46okLZXHawGATCYTLQPwp9ks5ePjA47jRGOO/Pz84OnpiQsXLhhdtqCgAI8//jjc3d2xceNGnDhxAtu2bQNgeKC8Oe/B3bt3hXbUFkqeiNlozBMhRITjABcX2/y7f+C3vOkcJBKJ0HNy5MgRhIaG4vXXX0dUVBRatGgh6pUC+F4b7eBvQ7p164ZffvkFS5cuxerVq4XpzZo1g0wmw4kTJ4Rpubm5uHTpksm2tmnTBn/99Zdo2l9//YWWLVtCKpUC4BOFjIwMYf7ly5er3fPSpk0bHDt2TDTt6NGjRpeRy+Vo27Ytzp07J0yTSCQYM2YMNm7ciJs3b1ZZJj8/H+Xl5bhw4QLu3r2LpUuXolevXmjdurXJweLmvAdnzpxB06ZN4evra7RcTVDyRIyihIkQ0hCUlJRApVJBpVLh/PnzmDFjBvLz84WrxFq0aIH09HR89913+Pfff/Hhhx8KvSBaCxcuxLfffouFCxfi/PnzSE1NxbvvvltlXT169MCuXbuwdOlS4aaZbm5uiIuLw+zZs3HgwAGcPXsWkydPhkQiEXqADHn55Zexb98+vPXWW7h06RK++uorrFmzBq+88opQ5uGHH8aaNWuQnJyMkydP4vnnnxf1MpnjxRdfxJ49e/D+++/j8uXLWLNmDfbs2WNyudjYWBw+fFg07Z133kFwcDCio6Px9ddf49y5c7h8+TK++OILdOrUCfn5+QgJCYFcLsdHH32E//77Dzt27MBbb71ldF3z5s3DiRMn8MILL+D06dO4cOEC1q1bh9u3bwtlDh06VCs37BQxOiKKWKQhDRj/IemaMGD89DUaMF5fUYzMQ3EyrsqA8XoiLi6OgX/aFAPA3NzcWNeuXdn3338vKjd79mzm4+PDXF1d2ejRo9nKlSurDJb+4YcfWGRkJJPL5czX15eNGDFCmKcdMM4YPxh6z549zMXFhX344YeMMX7Q+FNPPcWcnZ2ZUqlkK1asYN26dWNz587VW4eu77//nrVt25bJZDIWEhLC3nvvPdH8GzdusAEDBjAXFxfWokULtnv3br0DxpOTk4VltIO9Dxw4IEz7/PPPWdOmTZmTkxMbMmQIe//9940OGGeMH3ju5OTEcnLEx4icnBw2d+5c1qJFCyaXy1lAQADr378/27ZtG9NoNEyj0bAvvviChYWFMYVCwWJiYtiOHTtE7aw8YJwxxv744w/Wo0cPplAomKenJ4uNjRXmFxUVMQ8PD5aYmGiwvdYYMG4XydPChQtFOzYA1qpVK2F+UVERe+GFF5i3tzdzcXFhI0aMYCqVSlTH1atX2aBBg5iTkxPz8/Njr7zyCisrKxOVOXDgAOvUqROTy+XsgQceEHYqXWvWrGGhoaFMoVCwbt26sWPHjlV7expS8vT9yYrk6e9r2TZrhzH2ECd7RzEyD8XJuPqaPNmCvivJKsvPz2ceHh6iK+DqqyeeeIItWbKkWsuYE6Pq+vjjj6tcmVhZg7rarl27dsLdXzMyMkRdgLNmzcLPP/+MrVu34uDBg7h58yZGjBghzFer1Rg8eDBKS0tx5MgRfPXVV9iwYQMWLFgglElLS8PgwYOFKy5mzpyJKVOm4NdffxXKbN68GQkJCVi4cCFOnTqFjh07IjY21uQ5WEIIIcSU5ORkfPvtt/j3339x6tQpjBs3DgDw+OOP27hlNffee+8JN760JZlMho8++qj2V1SjFM9KFi5cyDp27Kh3Xk5OTpV7Y5w/f54BELrldu/ezSQSiag3at26dczd3Z2VlJQwxhibM2cOa9eunaju0aNHs9jYWOHvbt26senTpwt/q9VqFhQUxJYuXVqt7WmoPU8p6dk2a4cx9hAne0cxMg/FyTjqeTKfvl6VU6dOsc6dOzMXFxfm5eXF+vfvz06fPm3DVtpWbfQ8mcMaPU8OtZ+emefy5csICgqCo6MjYmJisHTpUoSEhCApKQllZWXo37+/ULZ169YICQlBYmIiunfvjsTERERERCAgIEAoExsbi2nTpuHs2bPo1KkTEhMTRXVoy2hvi19aWoqkpCTMmzdPmC+RSNC/f/8qt6+vrKSkRHQn07y8PAB8j1h178GhS7tsTeqoKbXOJawajcambTHEHuJk7yhG5qE4GacbF0ZXkxiljY9unCIjI3Hy5EmDZRsbfTGqq/UyxvQeo8397NtF8hQdHY0NGzagVatWyMjIwJtvvolevXrhzJkzUKlUkMvl8PT0FC0TEBAAlUoFAFCpVKLESTtfO89Ymby8PBQVFSE7OxtqtVpvGVP3qli6dCnefPPNKtNPnz5tlW7M1NTUGtdhqfT0IuH1pUuXoLldvas36pIt41RfUIzMQ3EyTPtYDe0l/sQ4ipNpdR2jkpISlJWV6T225+fnm1WHXSRPAwcOFF536NAB0dHRCA0NxZYtW8y6IZatzZs3DwkJCcLfeXl5CA4ORocOHeDu7m5xvWq1GqmpqYiIiBDu5VHX/tFcB07wT79u0aIFOgZ72qQdxthDnOwdxcg8FCfjiouLhXsfOTk5mbzEvjFjjKGoqIjiZIStYsRxHGQyGZo3by66MzpQcebIFLtInirz9PREy5Yt8c8//+CRRx5BaWkpcnJyRL1PmZmZUCqVAAClUonjx4+L6sjMzBTmaf/XTtMt4+7uDicnJ0ilUkilUr1ltHUYolAooFAoqkzX1llT1qrHEhKu4poCiQ3bYQ5bxqm+oBiZh+Kkn1wuB8DfuZrjOEoKzEBxMq2uY1RUVASO4+Do6Fjlc27u594uk6f8/Hz8+++/GD9+PLp06QKZTIZ9+/Zh5MiRAICLFy8iPT1deOZOTEwM3nnnHWRlZcHf3x8AsHfvXri7uwvPFoqJicHu3btF69m7d69Qh1wuR5cuXbBv3z4MGzYMAD/GZ9++fYiPj6+LzbZLumeiG+t5eUIIz8HBAc7Ozrh79y5cXV0pwTSCMYaSkhKzboLZWNV1jBhjKCwsRFZWFjw9PWu0/9pF8vTKK69gyJAhCA0Nxc2bN7Fw4UJIpVKMHTsWHh4emDx5MhISEuDt7Q13d3fMmDEDMTEx6N69OwBgwIABaNu2LcaPH4/ly5dDpVLhjTfewPTp04Ueoeeffx5r1qzBnDlz8Mwzz2D//v3YsmULdu3aJbQjISEBcXFxiIqKQrdu3bBq1SoUFBRg0qRJNokLIYTYE47jEBAQgPPnz+Pq1auUFBjBGENZWRlkMhnFyQBbxcjT09PkGSVT7CJ5un79OsaOHYs7d+7Az88PDz74II4ePSo81G/lypWQSCQYOXIkSkpKEBsbi48//lhYXiqVYufOnZg2bRpiYmLg4uKCuLg4LF68WCgTHh6OXbt2YdasWVi9ejWaNm2Kzz77DLGxsUKZ0aNH49atW1iwYAFUKhUiIyOxZ8+eKoPIGxWm9yUhpJGSy+VgjCEsLAwSid3cKtDuqNVqXLhwAc2bN6ceOgNsESOZTGaVddlF8vTdd98Zne/o6Ii1a9di7dq1BsuEhoZWOS1XWZ8+fZCcnGy0THx8fKM+TUcIIeZQKBSUFBihveRd37gawqvPMaKfDcQoRv1NhBBCiAglT8RsNF6cEEIIoeSJmEAJEyGEECJGyROpBsqkCCGEEEqeCCGEEEKqgZInYhT1NRFCCCFilDwRs9H4J0IIIYSSJ2ICJUyEEEKIGCVPxGyURxFCCCGUPBET6CaZhBBCiBglT8RsdAqPEEIIoeSJmEAJEyGEECJGyRMhhBBCSDVQ8kSM0u14YtQNRQghhFDyRAghhBBSHZQ8EeN0epuo34kQQgih5IkQQgghpFooeSKEEEIIqQZKnohR4gHjNmsGIYQQYjcoeSKEEEIIqQZKnohRur1N9KgWQgghhJInQgghhJBqoeSJGMXEXU+EEEJIo0fJEyGEEEJINVDyRIyiziZCCCFEjJInYjZKpAghhBBKnogJdG8nQgghRIySJ2I2SqQIIYQQSp6ICZQvEUIIIWKUPBFCCCGEVAMlT8Qo3fs80R3GCSGEEEqeCCGEEEKqxS6Tp2XLloHjOMycOVOYVlxcjOnTp8PHxweurq4YOXIkMjMzRculp6dj8ODBcHZ2hr+/P2bPno3y8nJRmT/++AOdO3eGQqFA8+bNsWHDhirrX7t2LcLCwuDo6Ijo6GgcP368Njaz3qEB44QQQogdJk8nTpzA//73P3To0EE0fdasWfj555+xdetWHDx4EDdv3sSIESOE+Wq1GoMHD0ZpaSmOHDmCr776Chs2bMCCBQuEMmlpaRg8eDD69u2LlJQUzJw5E1OmTMGvv/4qlNm8eTMSEhKwcOFCnDp1Ch07dkRsbCyysrJqf+MJIYQQYvccbN0AXfn5+Rg3bhw+/fRTvP3228L03NxcfP7559i0aRMefvhhAMCXX36JNm3a4OjRo+jevTt+++03nDt3Dr///jsCAgIQGRmJt956C6+++ioWLVoEuVyO9evXIzw8HB988AEAoE2bNjh8+DBWrlyJ2NhYAMCKFSvw7LPPYtKkSQCA9evXY9euXfjiiy8wd+5cve0uKSlBSUmJ8HdeXh4APqFTq9UWx0O7bE3qqCmNRiO8Vmtqtj21xR7iZO8oRuahOJlGMTIPxck0e4yRuW2xq+Rp+vTpGDx4MPr37y9KnpKSklBWVob+/fsL01q3bo2QkBAkJiaie/fuSExMREREBAICAoQysbGxmDZtGs6ePYtOnTohMTFRVIe2jPb0YGlpKZKSkjBv3jxhvkQiQf/+/ZGYmGiw3UuXLsWbb75ZZfrp06fh6upa7ThUlpqaWuM6LHXjRoHw+t9//4N7wQ2btcUUW8apvqAYmYfiZBrFyDwUJ9PsKUb5+flmlbOb5Om7777DqVOncOLEiSrzVCoV5HI5PD09RdMDAgKgUqmEMrqJk3a+dp6xMnl5eSgqKkJ2djbUarXeMhcuXDDY9nnz5iEhIUH4Oy8vD8HBwejQoQPc3d1NbLlharUaqampiIiIgFQqtbiemjiZnwacvggAeOCBZohs4WeTdhhjD3GydxQj81CcTKMYmYfiZJo9xkh75sgUu0ierl27hpdeegl79+6Fo6OjrZtTbQqFAgqFosp0qVRqlR3CWvVYguM4ndcSu9nB9bFlnOoLipF5KE6mUYzMQ3EyzZ5iZG477GLAeFJSErKystC5c2c4ODjAwcEBBw8exIcffggHBwcEBASgtLQUOTk5ouUyMzOhVCoBAEqlssrVd9q/TZVxd3eHk5MTfH19IZVK9ZbR1kEIIYSQxs0ukqd+/fohNTUVKSkpwr+oqCiMGzdOeC2TybBv3z5hmYsXLyI9PR0xMTEAgJiYGKSmpoquitu7dy/c3d3Rtm1boYxuHdoy2jrkcjm6dOkiKqPRaLBv3z6hTGOje3sCulMBIYQQYien7dzc3NC+fXvRNBcXF/j4+AjTJ0+ejISEBHh7e8Pd3R0zZsxATEwMunfvDgAYMGAA2rZti/Hjx2P58uVQqVR44403MH36dOGU2vPPP481a9Zgzpw5eOaZZ7B//35s2bIFu3btEtabkJCAuLg4REVFoVu3bli1ahUKCgqEq+8IIYQQ0rjZRfJkjpUrV0IikWDkyJEoKSlBbGwsPv74Y2G+VCrFzp07MW3aNMTExMDFxQVxcXFYvHixUCY8PBy7du3CrFmzsHr1ajRt2hSfffaZcJsCABg9ejRu3bqFBQsWQKVSITIyEnv27KkyiLyxoN4mQgghRMxuk6c//vhD9LejoyPWrl2LtWvXGlwmNDQUu3fvNlpvnz59kJycbLRMfHw84uPjzW5ro0GZFCGEEGIfY56I/aJHshBCCCFilDwRszHqeiKEEEJqnjzFxcXhzz//tEZbCCGEEELsXo2Tp9zcXPTv3x8tWrTAkiVLcOOG/T6+g1Sfbm8TncIjhBBCrJA8bd++HTdu3MC0adOwefNmhIWFYeDAgfj+++9RVlZmjTYSQgghhNgNq4x58vPzQ0JCAv7++28cO3YMzZs3x/jx4xEUFIRZs2bh8uXL1lgNsQHqbSKEEELErDpgPCMjA3v37sXevXshlUoxaNAgpKamom3btli5cqU1V0VsgBIpQgghxArJU1lZGX744Qc89thjCA0NxdatWzFz5kzcvHkTX331FX7//Xds2bJFdLNKQgghhJD6qsY3yQwMDIRGo8HYsWNx/PhxREZGVinTt29feHp61nRVxMao44kQQgixQvK0cuVKjBo1Co6OjgbLeHp6Ii0traarIjbA6FwdIYQQIlLj03YHDhzQe1VdQUEBnnnmmZpWTwghhBBiV2qcPH311VcoKiqqMr2oqAhff/11TasnNqbb8US9UIQQQkgNTtvl5eWBMQbGGO7duyc6badWq7F79274+/tbpZGEEEIIIfbC4uTJ09MTHMeB4zi0bNmyynyO4/Dmm2/WqHHE9piB14QQQkhjZXHydODAATDG8PDDD+OHH36At7e3ME8ulyM0NBRBQUFWaSQhhBBCiL2wOHl66KGHAABpaWkICQkBx3FWaxSxTzTkiRBCCLEweTp9+jTat28PiUSC3NxcpKamGizboUMHixtHbI8SJkIIIUTMouQpMjISKpUK/v7+iIyMBMdxeq/E4jgOarW6xo0khBBCCLEXFiVPaWlp8PPzE16ThovRkHFCCCFExKLkKTQ0VO9rQgghhJCGzio3ydy1a5fw95w5c+Dp6YkePXrg6tWrNa2e2Jj4Jpm2awchhBBiL2qcPC1ZsgROTk4AgMTERKxZswbLly+Hr68vZs2aVeMGEkIIIYTYkxo/GPjatWto3rw5AGD79u144oknMHXqVPTs2RN9+vSpafXExqiziRBCCBGrcc+Tq6sr7ty5AwD47bff8MgjjwAAHB0d9T7zjtRflEgRQgghVuh5euSRRzBlyhR06tQJly5dwqBBgwAAZ8+eRVhYWE2rJ7ZGA50IIYQQkRr3PK1duxYxMTG4desWfvjhB/j4+AAAkpKSMHbs2Bo3kNgPyqMIIYQQK/Q8eXp6Ys2aNVWm00OBGwbKlwghhBCxGidPAJCTk4Pjx48jKysLGo1GmM5xHMaPH2+NVRA7wCiVIoQQQmqePP38888YN24c8vPz4e7uLnpAMCVPhBBCCGloajzm6eWXX8YzzzyD/Px85OTkIDs7W/h39+5da7SR2BCNczLuyD+3cfS/O7ZuBiGEkDpU456nGzdu4MUXX4Szs7M12kPsGCVSYnnFZXjqs2MAgItvPwqFg9TGLSKEEFIXatzzFBsbi5MnT9aojnXr1qFDhw5wd3eHu7s7YmJi8Msvvwjzi4uLMX36dPj4+MDV1RUjR45EZmamqI709HQMHjwYzs7O8Pf3x+zZs1FeXi4q88cff6Bz585QKBRo3rw5NmzYUKUta9euRVhYGBwdHREdHY3jx4/XaNvqOxrnZFheUZnwurRcY6QkIYSQhqTGPU+DBw/G7Nmzce7cOUREREAmk4nmDx061GQdTZs2xbJly9CiRQswxvDVV1/h8ccfR3JyMtq1a4dZs2Zh165d2Lp1Kzw8PBAfH48RI0bgr7/+AgCo1WoMHjwYSqUSR44cQUZGBiZMmACZTIYlS5YAANLS0jB48GA8//zz2LhxI/bt24cpU6YgMDAQsbGxAIDNmzcjISEB69evR3R0NFatWoXY2FhcvHgR/v7+NQ1VvUdplJju+D5CCCGNCKshjuMM/pNIJBbX6+XlxT777DOWk5PDZDIZ27p1qzDv/PnzDABLTExkjDG2e/duJpFImEqlEsqsW7eOubu7s5KSEsYYY3PmzGHt2rUTrWP06NEsNjZW+Ltbt25s+vTpwt9qtZoFBQWxpUuXVqvtubm5DADLzc2t1nKVlZeXs5MnT7Ly8vIa1VMT7/5ynoW+upOFvrqT/ZRyw2btMMZWcbqRXSjEJreotE7XXV32sC/VBxQn0yhG5qE4mWaPMTL3+F3jnifdWxNYg1qtxtatW1FQUICYmBgkJSWhrKwM/fv3F8q0bt0aISEhSExMRPfu3ZGYmIiIiAgEBAQIZWJjYzFt2jScPXsWnTp1QmJioqgObZmZM2cCAEpLS5GUlIR58+YJ8yUSCfr374/ExESjbS4pKUFJSYnwd15enrAtarW6RrHQ/d8WmM5AJ41GY9O2GGKrOOnu++XlaqjVNT4LXmvsYV+qDyhOplGMzENxMs0eY2RuW6xynyet4uJiODo6WrRsamoqYmJiUFxcDFdXV2zbtg1t27ZFSkoK5HI5PD09ReUDAgKgUqkAACqVSpQ4aedr5xkrk5eXh6KiImRnZ0OtVustc+HCBaNtX7p0qd6bgp4+fRqurq6mN96E1NTUGtdhKVXmPeH1lStXkKLJNFLatuo6TneKKj5kqampcJXbb/KkZct9qT6hOJlGMTIPxck0e4pRfn6+WeVqnDyp1WosWbIE69evR2ZmJi5duoRmzZph/vz5CAsLw+TJk82qp1WrVkhJSUFubi6+//57xMXF4eDBgzVtXp2YN28eEhIShL/z8vIQHBwsDIK3lFqtRmpqKiIiIiCV2uZKrr2ZF4ELaQCA0NBQRHYMskk7jLFVnFS5xcDOPwAA7du3h6ezvM7WXV32sC/VBxQn0yhG5qE4mWaPMdKeOTKlxsnTO++8g6+++grLly/Hs88+K0xv3749Vq1aZXbyJJfL0bx5cwBAly5dcOLECaxevRqjR49GaWkpcnJyRL1PmZmZUCqVAAClUlnlqjjt1Xi6ZSpfoZeZmQl3d3c4OTlBKpVCKpXqLaOtwxCFQgGFQlFlurbOmrJWPRbhKnpTJBKJ3ezg+tRVnLQDxaWu3mg6/WsAgL+/PzTF/C8WZsf3dLDpvlSPUJxMoxiZh+Jkmj3FyNx21Pg8w9dff41PPvkE48aNE620Y8eOJk93GaPRaFBSUoIuXbpAJpNh3759wryLFy8iPT0dMTExAICYmBikpqYiKytLKLN37164u7ujbdu2QhndOrRltHXI5XJ06dJFVEaj0WDfvn1CGUJEREkSXXlHCCGNhVVukqntMdKl0WhQVlamZ4mq5s2bh4EDByIkJAT37t3Dpk2b8Mcff+DXX3+Fh4cHJk+ejISEBHh7e8Pd3R0zZsxATEwMunfvDgAYMGAA2rZti/Hjx2P58uVQqVR44403MH36dKFH6Pnnn8eaNWswZ84cPPPMM9i/fz+2bNmCXbt2Ce1ISEhAXFwcoqKi0K1bN6xatQoFBQWYNGlSTcNUb+ne58mOO1TqlLZnKeteMbq9wyfbWbduwce1au8jIYSQhqfGyVPbtm1x6NAhhIaGiqZ///336NSpk1l1ZGVlYcKECcjIyICHhwc6dOiAX3/9FY888ggAYOXKlZBIJBg5ciRKSkoQGxuLjz/+WFheKpVi586dmDZtGmJiYuDi4oK4uDgsXrxYKBMeHo5du3Zh1qxZWL16NZo2bYrPPvtMuMcTAIwePRq3bt3CggULoFKpEBkZiT179lQZRE5IZZRXEkJI41Hj5GnBggWIi4vDjRs3oNFo8OOPP+LixYv4+uuvsXPnTrPq+Pzzz43Od3R0xNq1a7F27VqDZUJDQ7F7926j9fTp0wfJyclGy8THxyM+Pt5omUaFsgLDdGKjoW45QghpNGo85unxxx/Hzz//jN9//x0uLi5YsGABzp8/j59//lnoOSINAz2qRUw3GpQ7EUJI42GV+zz16tULe/futUZVhNQbugkTJU+EENJ41LjnqVmzZrhz506V6Tk5OWjWrFlNqyc2Rr0r5qHTdoQQ0njUOHm6cuWK3tuZl5SU4MaNGzWtnhC7JboS0YbtIIQQUrcsPm23Y8cO4bX2lgJaarUa+/btQ1hYWI0aR2zPnm/4aGu6odFoKE6EENJYWJw8DRs2DAB/t+W4uDjRPJlMhrCwMHzwwQc1ahyxL5RHidEpTUIIaZwsTp60T5QPDw/HiRMn4Ovra7VGEftBSYFhur1ydCUiIYQ0HjW+2i4tLc0a7SD1AKUHhtFZO0IIaTyscquCffv2Yd++fcjKyhJ6pLS++OILa6yC2AjlBIaJb1VAkSKEkMaixsnTm2++icWLFyMqKgqBgYHCE+dJw0MJgmHU80QIIY1HjZOn9evXY8OGDRg/frw12kPsDOVLhlHPEyGENE41vs9TaWkpevToYY22EFKv0H2eCCGkcapx8jRlyhRs2rTJGm0hdogSBMNE93minidCCGk0anzarri4GJ988gl+//13dOjQATKZTDR/xYoVNV0FIXaPcidCCGk8apw8nT59GpGRkQCAM2fO1LQ6Ys8oQRDRDQf1PBFCSONR4+TpwIED1mgHsVOUExgmukkmxYkQQhoNi5OnESNGmCzDcRx++OEHS1dBiF2jx7MQQkjjZHHypPsgYNI40CNIxGjAOCGENE4WJ09ffvmlNdtBSL1GqRMhhDQeNb5VAWnYaFyPMRUBoZ4nQghpPCh5IsRC4juM264dhBBC6hYlT8QoZuA1qTxgnKJDCCGNBSVPhFhIPGDcdu04eOkWhq39C5cz79muEYQQ0ohQ8kSMog4Vw0SPrrFhoOK+OI6Uazl47pskm7WBEEIaE0qeiNkokTLMlj1PWtkFpbZuAiGENAqUPBGj6N5OhokGjFOcCCGk0aDkiZiNEgQxW15tx3Gc8E/r9p07VaYRQgixPkqeCLEQo/s8EUJIo0TJEzGKcgLDbNnzxBgT/ml5e/tUmUYIIcT6KHkiZqNjsmHU80QIIY0HJU/EKEoJzGMPcaIeJ0IIqRuUPBGz1eTQ3BAHMotP21HiQgghjYVdJE9Lly5F165d4ebmBn9/fwwbNgwXL14UlSkuLsb06dPh4+MDV1dXjBw5EpmZmaIy6enpGDx4MJydneHv74/Zs2ejvLxcVOaPP/5A586doVAo0Lx5c2zYsKFKe9auXYuwsDA4OjoiOjoax48ft/o21xeUExgmvkmmDRtCCCGkTtlF8nTw4EFMnz4dR48exd69e1FWVoYBAwagoKBAKDNr1iz8/PPP2Lp1Kw4ePIibN29ixIgRwny1Wo3BgwejtLQUR44cwVdffYUNGzZgwYIFQpm0tDQMHjwYffv2RUpKCmbOnIkpU6bg119/Fcps3rwZCQkJWLhwIU6dOoWOHTsiNjYWWVlZdRMMe0YZgoi9PJ6FEEJI3XKwdQMAYM+ePaK/N2zYAH9/fyQlJaF3797Izc3F559/jk2bNuHhhx8GAHz55Zdo06YNjh49iu7du+O3337DuXPn8PvvvyMgIACRkZF466238Oqrr2LRokWQy+VYv349wsPD8cEHHwAA2rRpg8OHD2PlypWIjY0FAKxYsQLPPvssJk2aBABYv349du3ahS+++AJz587V2/6SkhKUlJQIf+fl5QHgEzq1Wm1xXLTL1qSOmtJoNDqvWbXa4uHhIbx2cXEBALi6ugrTcnNzrdBC28VJd33lNXyvrYHBcAzsYV+qDyhOplGMzENxMs0eY2RuW+wieapMe1D19vYGACQlJaGsrAz9+/cXyrRu3RohISFITExE9+7dkZiYiIiICAQEBAhlYmNjMW3aNJw9exadOnVCYmKiqA5tmZkzZwIASktLkZSUhHnz5gnzJRIJ+vfvj8TERIPtXbp0Kd58880q00+fPi1KFiyVmppa4zosdedORYJz/fp1pKTcNXvZgwcPGp2fkpJiabP0qus4Xb5b8TiUtLQ0pJRm1On6K1OXq03G1Jb7Un1CcTKNYmQeipNp9hSj/Px8s8rZXfKk0Wgwc+ZM9OzZE+3btwcAqFQqyOVyeHp6isoGBARApVIJZXQTJ+187TxjZfLy8lBUVITs7Gyo1Wq9ZS5cuGCwzfPmzUNCQoLwd15eHoKDg9GhQwe4u7tXY+vF1Go1UlNTERERAalUanE9NeGddgZIuw4AaNK0CSIjQy2qR9sLZa3eJl02i9O1HGDfUQBAaGgYItsr627durbyPbcSqRSRkZF6i9jDvlQfUJxMoxiZh+Jkmj3GSHvmyBS7S56mT5+OM2fO4PDhw7ZuitkUCgUUCkWV6VKp1Co7hLXqsYRE5wo5iURicTu049dqczvqOk6cRCJ6bQ8fflNtsOW+VJ9QnEyjGJmH4mSaPcXI3HbYxYBxrfj4eOzcuRMHDhxA06ZNhelKpRKlpaXIyckRlc/MzIRSqRTKVL76Tvu3qTLu7u5wcnKCr68vpFKp3jLaOhobW95F296JB4xTcAghpLGwi+SJMYb4+Hhs27YN+/fvR3h4uGh+ly5dIJPJsG/fPmHaxYsXkZ6ejpiYGABATEwMUlNTRVfF7d27F+7u7mjbtq1QRrcObRltHXK5HF26dBGV0Wg02Ldvn1CGWKZhPjZE99l2NmwGIYSQOmUXp+2mT5+OTZs24aeffoKbm5swRsnDwwNOTk7w8PDA5MmTkZCQAG9vb7i7u2PGjBmIiYlB9+7dAQADBgxA27ZtMX78eCxfvhwqlQpvvPEGpk+fLpxSe/7557FmzRrMmTMHzzzzDPbv348tW7Zg165dQlsSEhIQFxeHqKgodOvWDatWrUJBQYFw9R2xTJlaAw6Ag9Qu8nWroJtkEkJI42QXydO6desAAH369BFN//LLLzFx4kQAwMqVKyGRSDBy5EiUlJQgNjYWH3/8sVBWKpVi586dmDZtGmJiYuDi4oK4uDgsXrxYKBMeHo5du3Zh1qxZWL16NZo2bYrPPvtMuE0BAIwePRq3bt3CggULoFKpEBkZiT179lQZRN5YiG8EaVmCUK7WIGbpPjjKpDg0p2+Du9M4YB+nNCmBI4SQumEXyZM5X/qOjo5Yu3Yt1q5da7BMaGgodu/ebbSePn36IDk52WiZ+Ph4xMfHm2wTMc+NnCLczucv6y8p18BRZh8DA2uKiV7bPnGxfQsIIaRxaDjnUEitEJ2asrAODg2vpwmoNGBcY7gcIYSQhoWSJ0IspNtjag9X29lBEwghpFGg5IkYJTo1ZYWDsz0kGdbCDLwmhBDSsFHyROpUQ72knwZrE0JI40HJEzHKGjmB7sV1DarnSXSTTNu1gxBCSN2i5ImYzRr5AWtAA6sZ7GzME508JISQOkHJEzHK2gdke0gyrEZ0tZ3tt6shhZYQQuwZJU/EbNYY19OQkifdLVHbQfJECCGkblDyRIyzQk7QUMcG6W6X2g62yw6aQAghjQIlT6TW6fY2NdSr0uzhtB0hhJC6QckTMcoaKYFGdDNJK1RoJ3THg6kbaFJICCGkKkqeSK3TiE7b1f8kg+M4cByHAQMqHij96tx5wnRCCCENGyVPxGyW5z32dUm/1egkShxnBx+lBhRaQgixZ3bwjU/smXWusNOtr8bV2RxjDIwx/PLLL8K0RYsXC9Nt1i7KngghpE5Q8kTMZunBWWNnD9C1Ft0toQHjhBDSeFDyRIyyyoBxnbuKN9Qcwx4GjNtBEwghpFGg5InUOnt7jInV6N7nqQE9doYQQohxlDwRo5gVxiuJ62g4yVODTQoJIYQYRckTqXUN9j5Pop6nBrRhhBBCjKLkiRjFDLyujoZ2tZ2WvSVPtm8BIYQ0DpQ8kVrHGujVdroa6nYRQgipipInYpS17/PUkJIM0a0K7GC7GtJ4MkIIsWeUPBGzWT5gXPfBwFZqjB3Q3S57uNquAYWWEELsGiVPxCjrPBhY93XDOcTTTTIJIaRxouSJmM3SO4yzxnC1XQNKCgkhhBhHyVM9sPdcJj479J+tm2GxhtrzJHrgcUPKCgkhhBjlYOsGENOe35gMAOga5o2OwZ51u3Kr3CST6X3dkNhDz5MdNIEQQhoF6nmqR27dK7F1Eywi7nmyXTuszd7u80QIIaRuUPJEjLJ0nJOhOhrS6S3dLdm8ZavN2kEIIaRuUfJEal1j6HniOPooEUJIY0Hf+HaK4zh4eHiIpg19fCg4jgPHcXXWDmuMo9E00DFPol45CX2UCCGksbCbb/w///wTQ4YMQVBQEDiOw/bt20XzGWNYsGABAgMD4eTkhP79++Py5cuiMnfv3sW4cePg7u4OT09PTJ48Gfn5+aIyp0+fRq9eveDo6Ijg4GAsX768Slu2bt2K1q1bw9HREREREdi9e7fVt7c+sjTxaYi3KuA4DqNHj9aZIBES27pMbgkhhNQ9u0meCgoK0LFjR6xdu1bv/OXLl+PDDz/E+vXrcezYMbi4uCA2NhbFxcVCmXHjxuHs2bPYu3cvdu7ciT///BNTp04V5ufl5WHAgAEIDQ1FUlIS3nvvPSxatAiffPKJUObIkSMYO3YsJk+ejOTkZAwbNgzDhg3DmTNnam/j9WCMITc3VzRtx087wBir094ba6yKNdhbFVQkSRz1PBFCSKNhN7cqGDhwIAYOHKh3HmMMq1atwhtvvIHHH38cAPD1118jICAA27dvx5gxY3D+/Hns2bMHJ06cQFRUFADgo48+wqBBg/D+++8jKCgIGzduRGlpKb744gvI5XK0a9cOKSkpWLFihZBkrV69Go8++ihmz54NAHjrrbewd+9erFmzBuvXr6+DSBhWX9OOhnifJ8YYdvx9Ey9+y99GApykQZ2SJIQQYpjdJE/GpKWlQaVSoX///sI0Dw8PREdHIzExEWPGjEFiYiI8PT2FxAkA+vfvD4lEgmPHjmH48OFITExE7969IZfLhTKxsbF49913kZ2dDS8vLyQmJiIhIUG0/tjY2CqnEXWVlJSgpKTiNgJ5eXkAALVaDbVabfF2q9VqUbKh0WhqVJ8ldNdfrrZs/eU6y1hahzHa+uo6NrrrkzrI6nz92jFxvvGbhWmurq7Ca92eS1vFqL6hOJlGMTIPxck0e4yRuW2pF8mTSqUCAAQEBIimBwQECPNUKhX8/f1F8x0cHODt7S0qEx4eXqUO7TwvLy+oVCqj69Fn6dKlePPNN6tMP336tOhgZgndzoy0tP+QUnKzRvVVV25ujvD6ZkYGUlLyDRc24L/rFadW//33X3gV3bBG06pITU2tlXoNuZpeJLzu1LkzUlJS6nT9Bw8eBACM3KqqMg2A3vbUdYzqK4qTaRQj81CcTLOnGFUeJ21IvUie7N28efNEvVV5eXkIDg5Ghw4d4O7ubnG9arUap/4+LfwdHt4MkW38jSxhfe5nTgE3sgAAvv4BiIxsWe06VA4qIDEFABBWC9ugVquRmpqKiIgISKVSq9ZtzFXcBI7x74+TsysiIyPrbN0iW/cILw21wVYxqm8oTqZRjMxDcTLNHmOkPXNkSr1InpRKJQAgMzMTgYGBwvTMzEzhYKFUKpGVlSVarry8HHfv3hWWVyqVyMzMFJXR/m2qjHa+PgqFAgqFosp0qVRa4x1CI7oaXlLnO5julWMMnGXr170HEmdhHWawRryrQ6IzSFzDmF18+NOzixHu62Jwfl3HqL6iOJlGMTIPxck0e4qRue2oF5cIhYeHQ6lUYt++fcK0vLw8HDt2DDExMQCAmJgY5OTkICkpSSizf/9+aDQaREdHC2X+/PNPlJWVCWX27t2LVq1awcvLSyijux5tGe166pqtByHrjvWy9O7gjeE+T2o72awrdwps3QRCCGnw7CZ5ys/PR0pKijBOIy0tDSkpKUhPTwfHcZg5cybefvtt7NixA6mpqZgwYQKCgoIwbNgwAECbNm3w6KOP4tlnn8Xx48fx119/IT4+HmPGjEFQUBAA4KmnnoJcLsfkyZNx9uxZbN68GatXrxadcnvppZewZ88efPDBB7hw4QIWLVqEkydPIj4+vq5DAsC+7otk6fPbNA3wPk9ApVsw2MmGeTrJbN0EQghp8OzmtN3JkyfRt29f4W9tQhMXF4cNGzZgzpw5KCgowNSpU5GTk4MHH3wQe/bsgaOjo7DMxo0bER8fj379+kEikWDkyJH48MMPhfkeHh747bffMH36dHTp0gW+vr5YsGCB6F5QPXr0wKZNm/DGG2/gtddeQ4sWLbB9+3a0b9++DqJQlT111Kit0JiGcqsCwD4fDNyQ4ksIIfbKbpKnPn36GD2lw3EcFi9ejMWLFxss4+3tjU2bNhldT4cOHXDo0CGjZUaNGoVRo0YZb3Ad0dhgnbrjnPxGvCG8Xv3hR1g09H8Aqnf6rcH2POm8tmXS4uemwK17/K0y1LbYYQghpJGxm9N2RD/dZMPW44U4iWUD+jQ6B3Rbb0NtsWXPkz32gBFCSENGyZOd00026uqwqH0EDL/uil6oqc89Z9HjYcQ9T/ZxcLfGM+h042CNU5o1aInwyl7iSwghDRklT3aO2VHPk6W9GqLTWw3otJLudtnirdEmgJmZFbfo6P/IAHo4MSGE1DJKnuycRnRKxgYN0DkIW7p+Zoc9T1ZhL6fLdBMlekAxIYTUOvqmtXO6+YotEo8hjz1m8fq1PSBTnq24mvGZyVNs1jOiXa/uuvVNM5foPk82SJ60p1B9fH2FaTt37rLo1CohhBDzUfJk52zda6O7xnKLEwSdxIRrOLuc6D5PNkxWaMC4Zej0JiHEUg3nSNZAaezkAA3U4EaQotNKtrsFv+5AeM5BXmVaTdj2arsGelqUEELsFCVPdk6UPNlgzJPoijIrJE+W3u7Amo78exshL/8I95gna1SPvdznSXfNdJ8nQgipfZQ82TnRgdEGB2jdZ7ZVd/3aHp01a9YK01asXGnzMTlvbDsDAPDqPaFG9djN6bIavEeNjfZUnTygGQInfQSnZlE1GvdGCGmcKHmyc7a+SWZZeUVXhjUeDGz5uCnrkUqsc5C09YDxinZUsJdn7Nk7/+GvQ+4fDv9Ri2zdFEJIPUTJk50Tnzar+/WX6azU0sRHY0YPTV3+8rda8iQaj2aVKi2iscap1UZC2+vpFdC0yjS6QpEQYi5KnuycrQeM6yZPlq5frTNYq1xt+wOUg9RaPU8VapK01DRxZHTajhBC6hQlT3bO1slTqbrmvRplojpsP6LZoRZuJGnLpEX39CGdtjMPRYkQUhOUPNk5W49n0e15sjx50n/qzxo3reQ4Dh4eHtVqj4OVTtvpdvnU6L2ROkDq7meNZlDPEyGE1AEHWzeAGCd+qG7dr98ap+10T9XpO7hzDnK4dxuBwsvHUHYrzaJ1VEetnLaz4O7rWsoJK6AIbAlFk9YovXmRr7sa9dk6wSaEkMaGep7slLZHRfcYOjMhoc4vqda92s7iniedU3VqnURKO0h3xS+p8Oz1NIKe+ahWB+9qY3dg//4q0yx6PAsTv7a0zYrAlgAA1/b9LVrebp6xRwghjQQlT3ZO91hoi/vQWGPMk27Pk74r9s7cyKtWfboJj8TRDZyjGwDAw8PDvBhp1NVanyGVk6XqxEdfkvjs1KkWJY6iWybcf1lQUg63yEchdfGsVl2EEEJMo+TJTjHGkJubK0qelr37Xp1fUi0a82TxaTvjvVdyB0uTQg5NnvsUPlM+Q3G56YHo2tg9NmigMK2kTG21mNZ0vJE5yZe+XjLxXej5P978+Sx8Hp0B/1Fv1qhNxH7QjTwJsR+UPNk5Wz8CJDu3olfo1KkUi+oo0xjvearu1W/aZKegpBQSR1cAwD93y5Cbm2tWEqR7ACqtwc2zKm9LTS8ktDg5LS+vaMP9On5JVQEAFMrmNWtUHaLkgBBSX1DyZOd0EyZb3MSPk+pcU2DhJf7inqeqGUZ1B3BrD7JuHl7CtJv31GafttMtUmpGj5Wh9b/86jzRdBc3txod/K3x4OU5r84Fx3HIuXtbZ3bFFYleIa3AOcgsbiOxMY6+sgmxB/RJtHPiu3PX7boZY+CkOgdaC7+4xfd5qjpfZmFSptu2wjLzg6N7eszc5Elfr4jEQVGpUM0+TobuH2rqlg6c7nrvx5Kpy6rUczqzBN7jVyNg9Ns1amdDoPtDxFjSauveMN332TVyIIJf+g6OTdvUqF223ibSsDWW/YuSJzsnfgRI3fQ8CYOxK/VQcBIHi65OK6tGz5M54360p+3S0q8L00rUzOzTdrt+2SO8tqTnSbv+hDlzRdPv3Llbo95Bc3qeHEM6wKvfVHAOcgAAVymB0yZSTF1xKk/mGwIA+D2tiK8juL1F7auNL0XR/sRJAAv3MXPXo8/3Sdf1Tre3A4BP7HRIFM7wGfhSra+rsRwACbEUJU92zpZjnkSn7ABI5E4W1WPoajvtF/S6tWuEaTInF7O/uHXHKxWXmx8bTlKxXaVq8668k7r5QOYbKppWXCZeVnfMkiUHny3ff693uu6VeQFjl8A9aijco0fy65FV6v3S0/Pk1mUoAEAh1d9zZQtV180h8JmPEPTMGlEPXm0la/fu3ROmxS/7xGAsFE3bQeLkbtU2VIe+qzLbtGlTo4scFMHt4RjWqdrLUUJFSAW6Saad002YyurouXDaL+XsglJ0emuvMF3i6AqNRlPtL9ByjYl7RbGK+ZxMAVZaZFa9ur1GpdW5+4BOUlhipOdJdztDX90JAHDw8Ic67xYAwGfQLLhG9BPKWDJmiYkSLvN/y8iVLcAYw/XsQjz47oGqdeg+LLgwBwAgk1Ys7xjWCcVXku8vw1VpS2U5haUInPghCs7/aXYbzaFd5538EnR5+3cAgNTVB+V5WeAkUsh8Q8FxErD7+4g5bbWgEVUmeXh4wPOhOHh0HwWmUfPjxO735plad620UYfCoWpyaW6bwEkQOmcHAEDq7AFNUZ7e5eXKFlAEteJP3dcgYartWNgD3fg05O00Rt8+Up24eHh4oKCgwKyy9oJ6nuyc7rG9qLTccMFakF1YKvqbc5AZTTYMKTPQ86T99Tx9RsVpiH+vXDf7V7Vu8mTqVgWiMUI6PU+do7pV6+CgvaElAHAyuWiebs+Tomk7OHgqTdYnSoglUsMFK9H2AhZVzho5vg5NcUXPiqY4HwAg1+l58o6dbtZ6tDELeyQO8oBm8OozsVZ6rbILK3rKtHH1engKgp75CG5RQ81qo6Hppto5YvhwYZ/TfdSPa4cBfD0SKbwfnmL2tkDqAEiM/y41p12GyiSdOGZ+WyrXKXMUXkudxY810h1DFxi3Et6PPA/n8E7376fmCs+H4qr0vprTXlPlLN2XrLkP2lcvbN0ub606tBRNWsOl3cNWqcueUfJkp7RXSBXpnI5a+8lndfohv1PAJ09l2Rlg928smVdUdSCyKXt+q+i9+mnHz1Xmf/y/T4TXhXq6kAxts+5pu5Jq9Dxx8ooDiHbskD7aA6pu0vr1xu+ExK7yeKPgkFBwHAeZTzCU495Fk+c+g7EDOMdxcHR2Ef6WyBSVBoJXLFflgCNzBMdxaB/ZRVypRALGGAY9UvHlpR1YLzPwTD+3zo/BMaRDlbaJ1ifVnwxUZ3/U3SaJkzsUTVoLf7fvHF2xCQr+1K37/aTJu9+zQjmZTzCc2/Q2mlzra5NjWCTcug7TW/6HH/SfLtXtEXXr/JhZ21iu1qDJs58gaMrHtfeonPKKHzVOLWPg4N1U+NvQ+6Hdl29m3ako6yDX+0NFO0YOAKTu/gAA74enwKP7KAROXC0kmOY8U1Ie2BJuUUOrHQt7PUVoql21nczVdVzM/QGie3pZ+fT78H0sAfKg1nXUStug5KkWVfeBtfroXkXGyZ1rXJ85tB+WhwbwBwxNYS40JXyXatNmLY0tqr8+nR4VTs+VdaKr5qrRu6bb81SiZ8yT7ode98MtUVTE8bff95vs5SooqWjTxGefq0hoKg/Wvp+Iyf10fp2bOBWnm5SYGlvD6RlzptuTAFScttONjeR+Gd3jV2nGJQD8+BfvR55HwNgler8otTF7c+ECYZqiabsa31g0cNJHUD79PhxDO/JtdK7YdlH8Kgmasg5+Q+fA+YEovT0YurRtLC4rR8Dot+H98BSh7aKrSHWSpNzcXOF/VuniBiFprrQe3bg5eivh4OEPmVcQZC4eBg88jmGd4OARYH5Pjc7FG4zxp86dwiLhP/x1NHl2vdkH1QKdzxenEH+faON18GiSMM130EuQBzwAxf2LDAwl0VXae/9f4IQV8O43FW7tHjL4fsn8m8E7Nh5Z94pN1i1ekQQOXkE2SbRE73l4ZwTP2grnVj1F8wwtU5kiJEL0Q8Labawt+n7c6X7Xy/3DTS6jPUZyChf4P7kYLu0f1vs9VJfbYi5KnmqZoZ4EU7R3GC8qqzhADRwyrE7vMC69f0BTF+ZCU8SfBpK4eFZ7RxMdqPScmtL9QtbX8wROojcJESVP1RgPJtFJQo1dbafdTmVwmDBN91RH5dN2UhcvMMaw/ccfKtbl7G78PdPZ9iqPUtGJFWMM1zJ1ew1kfBJQecC4VAaO4/D7gT902qlAbm6u6FE7Wg4eAfrbpYPjOMx/reKeVv6jFln0ZaObwDq4+QAAnJrzPU66cTXnajKFzq9a7wEvwOv+aTVtmyRO7kIbPcI7CmWV496tkoxol6u6PdX/nOm+H9obuFZehyKoFQJGv4Umz39eZd263xNe/abC86G4KnU5hXWCc+teQkKjy8G7KZxbP2jwwFOo00XrGvGI3h7OfgPFvWweD44TJ033tzE3N1e0vMw/HL5D50Dq7lelXdreLE7mCPfuo+DgFSSsL2jSh3CLfBTtJiw2mRDrrs9n4Aw0mfoJnFrGVJlnjL7vZImjK1wiHoFE4Vzt7+uAJxdDIneCz6MzxHEz8ONZt26piyeUY5dC+fT7wnecofdO33SO4+D96AwEPrMGErmT3jabSuYM1WuovL7pWhInt4rXlZJzY7F06jQETuGd4Ts4wWj9+tph6n0y9V5ampRR8mQD+j68+t5gDw8PFOokT7/t/8OsD4Y5H35zPggO3k0AAOV5t1CWfRMAINM5RWCsPtH6dXue7h/cRf90kqtHhzxeZb5y/AcIenY9OAe5aPrgocOE5YrLNcJNMoUvGwdFlZ4ZjUbc8/T4iJF6t0e0bTo9PlJnz4r4uHgBgHBKU+rqDY7jMOyJJyvKu3jrj4lOPCrKegljZaTufgh+cRO8B7wgLNOsdYRQVuLoyn/pV0qeHFy9q9SrPaDrJk+czJFP6HSe86c7Rsu5TW/xQVha8R7qXnXp4N0EUjefah1sdEnun0LV1+tWlp0hvK481kabTDh4N4Vbp0Fw7zpMqMPjwXEIfnETHJtF8U2vnJRyEqEnFRD3IIoOeMz0xQQSmSNc2vcT1n006XRFG53c9CbOcmWLivX1HKu3fkWT1nCPGgqP7qPAKVyqJGLOLaL19gI1eXY9/B6fC8dmXfR+Jrv2eFD427X9w6K2CNtUuUdKXSZal8cQPpGunBz4PT4XLm16o+m0L/l46OyD2qs/PXqMgddDcQicsKJK27U9FZzcCYETV8Oz7zNCu/TtV64RjwAAPHuMqTKvuvuiz+AE+A56Cd79nzdYF8dxUDRpwyeI95N/XRJHV8j8qva2GLt5r8wnWHgtvf99AgAyvzDRqXSO4wCpA9y7jRBOpd6fA7eOsZD7hQk9X0L5anBu0xtBz30GecADojrkAQ8gOOF7uHcbIZpuqH7d70eZdxN49ZsKRUiE3mOULombTsJtZLyguUmSvoSzchlD9Zl7xoiSp1o0IbAVOnk1gcP9XxQO3k2E0xTG6L65RbqnX+73mFj6i8HUcpWXVdz/Yi3N/Bdlt9MBAPKAZkbXqa8eqWtFAuHg4Y/KdA/0kkrJjsTRFYrAFpB5BcGxaVvxcjq9BznF4gMdJ3VA4KSPEDjxQ3AOFQmbg5NLpXJ6krlK26DbUyW5fyDmHOSQ3t+W4vTU+9vJf6FKFBXrkHnrP60gfHil4h4Q7ReoW+RASBQucOs0qGLdOr0zEoXr/XXysdVejab91S9KnhwU8PDwECVPwq9UneTLd+gcfp6jK/yGzoHf43OFU8USWUXCVJr5L78uNz8EPbMWynHLARjepyr/032PXTsMADhJlcHLlTm3jBG1VTseSvfL3m/4awAAz/sJiU/sC3zbKyVmUhdP0Xuk7+HJHh4eeu6hxUHRtB1cIwcKf3v2Hg/fwbPgN4xPKB58+BGhvP+oReAk0qrvv6bi1Jnng+P09jgomlTs6549x0JauRdL5ii68CHgqWWi+W4dHxVeO3gqhQNS5c+XommbKtuu3be0WFmJKBGSBVVdBuAPllqBkz4SJQNevSdA5hMsJAQSR1c4t+ktXt4rCJA6wLl5NOQBD8Cj2wi4doyt2GYDB055wANwat6Nf61sXiWJqPK5ljrAo+dTon3H+f7yulfP6qN8+j24tOkN737P8RMqHeiDnvkIEp0EQpfUzQ/+o94Ujb3TjVnT6V/Bq9+zAIDACSsQMHYJ5EGthPmevcbDq+8zUD79njBN9zvB6YGu8HxoYsXwgaBWoh9Exr7r/YbOgcxTCf8nF4um+zz2MiQyR3j1feb+flqxnGOzKFEMAfF3vWuHAXCPGgrl2KXgZApI7p+ac9X5ThO2Qydhd2nTC8Ezt8Cjx5j77eTg0v5hOD3QVbSMqQTZVAIpOm7d/26pTtJJyVMt+ijjIk5l38A9iRTJviH4umk7vNyyBx5+oBtcKh00Dckv1TngKepmzBPA/yJS3E9WSq6lojid/0Xt3CJGdOAxReLoKpyiAQAHd/8qp5qk7r7C68o9RdqEBADkAeLntOn+Gs4u1gC6vWaegZB5B0HmHSQ6xVOlK7nyaS9926CzjPaAIFe2AMdJoC66h5JrZwAAjk3aVCmvaCI+0Di36S36QpM6Vzqw348V0719w/146/6qk8gdAalMiEnRP8f57fYKBDiJKLHUbqPuqU1tb5rugVJ7JaGDzhe69tl4uoPshW1p1QOc1AEOHgGQKx+oMt8QSeVtdvWu2vMklYl6uGR+YVV6X6Su3pD7Vvxyr3zzTwd3f0DqAGmlunV/7QPQe5oJEO+L5blZADgox70Ln9jpcO8+CgDgGsknKY4hfK+gKClzctc77qNKMqdzwNHXJt1eNaEOZw8hkQf4bdftDXFuGQPH0I5wahaFJs99xieSEmmVcU4O7v5wbtUTTaZ9CYWwDVW/Zyonb/puOaIuzK2o180HisAWovnej8ZDo3MbEr+hc6BoIh5U7Df8dVFPr8+jM4SrHnVVvtDDf+QCQOKAwLhV8Bs2Dx4Pjquyv2h59hoPzwefQuDE1XBu81CV+T6DZgISKZyaRcGr//N87CslSc6tesC146OQunpVWd65RTQgdYDEM1DU6+j9yHNwataFH3zf8yk4+DQVDc4HAPeoxyF19RG2Tzl2GZRPvw+pm6+Q4Dm4+cDr4SlwbtUTDjr7iUvrB+HR/Qm4xzwJpxYxCBz/AQLjVkERHAEHz0C9sahM6uyBwMkf30/GOch0vquaTv8KAeOWgXOQw7FZFwSMWgT/JxfzP3CD2/PlDYxXlCtbwCWiP5zCO8NnwAtw8AjgY+DP/xiX6pzR8H3sZUgUzvDs9bQQa9/BCfB/YiEcvILgGN5ZlEgpmrSB77B5wqngyjz7PoOgZ/9Xkbg7uYvOhjg174bgmVvgO2ye3uUNofs81aLjPiGIvncLbqVFiLydjsj7vTcAoOYk+M9DidOuXrgulaGU41AKDqUcUAYOJRzAFM6Q/7oLQ+GAUqkDyqQOyG71IAryslCqLkMpOGg4DhLGwIFBIpGC46SQyBWQQAKoyyDhOHASCSScBByn/Z8T/uYkEqglDlBLpFBLpdBIHAAPfzg0bQd2Ox0F18/CMUcFdW4WfK6ehsRLCf+RC5B38TBKcjJRVpIPaWArOHcehKKsq7j759eARAJO4gCJsztcWvWETF2Gsrs3wTnI4ODuB/++k1F0ORGstBic1AGuXk2Acr5b371pO7Dr58DKSwFOAufQDlDcv7rIM6wj1BcO8QcBqQPcm7YX5gGAe1hncNcuAozBKSwSirISAIBv24eQm5sFiaMLHEMi4FhWMTDV54GukNy8CFZSeP+Ljt0f6nI/nnInuD8QDadSfhkn/2bgujwOl/Z9ISstRsGFQ8Dlo3CKfgJO4Z0h7zMJcmVLKO6XlzWLgiZ5N5hGzX8ZPzQR6uJ7uPXjO2DlJXBtESPUDQA+YZEoKMqDl28IHO9PD+w6DAUXD8O9aTtR2ZBHZ8CxWWdISotx79ROyP3DIXXyQOjQOXBy8wPul5X7N4M8qC1UWflwuj+mTOHqjSLfEHj4Bgt1MnUpXB3d4NqkrTAtsOdY3L6dDk9PpbBNju5+yHf1hk/TdkIblVGP4/aulUZOdXH8b1ZOAufAlnDWOYj6BEfA1SMACp1pvqEd4eogA3d/mtw3BPB/QLxc07Zw9g2Bk840N0c3UZmgqMch8w4STVN2GgSn0iIwdTmfgMudUOgXDsY0cHJ2xu1befDwC4cr0wD3l2MyBfxa9RTqcY55Ejh7AK6MCWWa9BoPMCZal3dQaxSVFEAR1AZlt9NRnpcFT4+AqmX+Oym0h5PK4OUTIoqHsvMQOOr87eShBHP1gbSksGJbOw+Bi87fTXtPgIO7HyQlhXBp9SDcHd2gzrslKiNr0gby+8mJx4j5KPznODgHGZx0yjgGtRatBwBS/8uGm9wJklI133vrIIMbx4HTrbvrMMh1/nbxC0dZdgZkOtOa9HgKjrplmrZDMTjRNFnHWGRdSuTjI5GAkyng2LS9aDsAILjX08I0l6jHUeQbiju/rgXAAI4Dx0kAjkNgxCNCO11ipyPX3Q+uOnW5tYyBS94tuHYYAKlMgabt+qC8IBv3jmyB+/3bfnCMweOhOJS0iIZjUcXD0znG4Nw8GmgWBcembfHSxssI6zMFeef/grdXEKQFOQAA304DUd68K5haA3l+ts5WMLg9PBmO+XfBaRMvdz/4954AOQDJvdvgGKBs9SBYi+4oufI3nO7fdw4AODAENOuCcu+mcMzN4uMyeBagLkPO4U0ov3sd0GjASTj+bv4SKaQuXvDKUVWsTyoD5+qNJjFPwvXe7YrpAODsiZCYJyHzCYH87g1wjKHF0+8DAIqv/g0HNx/I7lwTDRfkwBDc7mFInFwhv3UVANBi5AJwYNCERuDQb6fRtqQQuHVFFEcAaBr1OBTKFpBn/QcAaN57vNBbd8/dH+W5mfDoORZSqQzlD09B4cW/4Ny2N8qy0lD03ylIXTzhERoJlJWgvOswFHs3gWvEIyjPv4v8lF8AMLh2GADZrSuApxL3Oj8GJnfG4aNbYArH6ssdqerY2rVr8d5770GlUqFjx4746KOP0K1bN7OWzcvLq7iKAEALVx90cfNBZwc5OqnL0D4nEwH3b1xICCGEEPuQB8AD/AUR7u6Gr4Cm5EmPzZs3Y8KECVi/fj2io6OxatUqbN26FRcvXoS/f9UxO5XpJk+GBDl5IMrdF51ljvABg5xp/2kgE16r0TzIH9fPnIK0vBxOLp5QyBSQg4Nco4ZcUw6OMTCOA99nwv+vAZ/4V0zH/en3/+b4myrzbzyDlGngoNFAwjSQMgaJphzS8jJINeWQAvDz9EReTg6kAKQAHMDBwYIrkQghhFSl21/LOE7o/Wbg7g8z4u5PYoB29v3xR/z0imELTDs2iWn4urRlta/vHxeE9eH+bTjU5WDaW8lwHD/GiTEwdTmYRMqXvH/GgnH88zP51/cPKBIJf1xRl91ft4Rvl0QqjJlj6jK+XomUv1hF6gCZXIHSkiJoykr4+iVS/tmcMjk4iQMYY9CUl/JjnBxkwi1E+GEbHBjT3L+oQca3mzGh/Qz3219exl9cw0nANPfbrT11xzTQlBbxPahSGe4xDVqXFFDyZIno6Gh07doVa9bwz1zTaDQIDg7GjBkzMHfu3CrlS0pKUFJSIvydl5eH4OBgODs7V/uqB13Ozs745ZdfMHDgQBQWFppeoJZo731TOSHkGLufTEH4vzrM3fFMlXNycsIPP/yAkSNHorDIvEe7WGvdAJCRkYHAQPPGFNQWU21wdnbGtm3bMHz48NrZlyQO/NgqiVS4lxdjmvvfzPcPDdpbFWjKRTd6tJj2y7lSO3QHZFeXs7Mztm3fjuHDhpmOk/YAYKqNHFdRTu8y/CBmTiLlr9xUl8OcPc9WX9zOzs7Yvn07hg0ficLiYrPbW2MSB349pmJuocoJRU1fOzs7Y3fl7+8aHA8aIns5xulijKGwsJCSp+oqLS2Fs7Mzvv/+ewwbNkyYHhcXh5ycHPz0009Vllm0aBHefPPNKtP/+OMPuLrqH7RICCGEEPuSn5+PPn36mEyeaMB4Jbdv34ZarUZAgPjmgQEBAbhw4YLeZebNm4eEhIqbe2l7njp06GA0+FqGTvFps/KIiAhIpVKr3LHcXNreJn0saYeh3itrsOTXS25ubrXbUnmZyjGqybZZ0h5T7dCtU3df8vauenVXfWJprPS9f5XrsdYvYXPaaEkZcz5HtbVuLVMxMifO5myXuW2svJwt6GuHufuS7ufXnO0wFCtzlzdWr77lzWmf7rLGvhcrr0M3RhkZGQaXq2267TK3P4mSJytQKBRQKKpe8i6VSiHVubmgIfn5+Xqnq9VqpKSkCPUYKlfXatKO2tgGbZwyMjLMindN2mJsmZpum7Vio1uP9rXuvmQv+1FNWLoNlZer/Lel+5I567JmGVPL1ea6zYmRqTibO83c99le9mnddliyL1VnOyzZL2q6fku+/4y9p8ZiVNfvqXZ95oxZBug+T1X4+vpCKpUiMzNTND0zMxNKpdLAUoQQQghpLCh5qkQul6NLly7Yt2+fME2j0WDfvn2IiYmxYcsIIYQQYg/otJ0eCQkJiIuLQ1RUFLp164ZVq1ahoKAAkyZNsnXTCCGEEGJjlDzpMXr0aNy6dQsLFiyASqVCZGQk9uzZU2UQOSGEEEIaH0qeDIiPj0d8fLytm0EIIYQQO0NjngghhBBCqoGSJ0IIIYSQaqDkiRBCCCGkGih5IoQQQgipBkqeCCGEEEKqga62qwVqNf/U7+vXr5v1bDtj9WRmZuL69es1flREQ0ZxMo1iZB6Kk2kUI/NQnEyzxxjl5eUBqDiOG0LJUy34559/AADt2rWzcUsIIYQQUl3//PMPunbtanA+x8x9hDAxW3Z2Nry9vXHt2rUa9TwRQgghpO7k5eUhODgYd+/ehZeXl8Fy1PNUC7Tdj+7u7pQ8EUIIIfWMqdOINGCcEEIIIaQaKHkihBBCCKmGepM8LVq0CBzHif61bt1amF9cXIzp06fDx8cHrq6uGDlyJDIzM0V1pKenY/DgwXB2doa/vz9mz56N8vJyUZk//vgDnTt3hkKhQPPmzbFhw4a62DxCCCGE1BP1JnkC+KvXMjIyhH+HDx8W5s2aNQs///wztm7dioMHD+LmzZsYMWKEMF+tVmPw4MEoLS3FkSNH8NVXX2HDhg1YsGCBUCYtLQ2DBw9G3759kZKSgpkzZ2LKlCn49ddf63Q77QljwJkzQFGRrVtCCCGE2Id6c7XdokWLsH37dqSkpFSZl5ubCz8/P2zatAlPPPEEAODChQto06YNEhMT0b17d/zyyy947LHHcPPmTQQEBAAA1q9fj1dffRW3bt2CXC7Hq6++il27duHMmTNC3WPGjEFOTg727Nljdlvz8vLg4eGB3Nzcej9g/OefgaFDgchIIDnZ1q0hhBBCao+5x+961fN0+fJlBAUFoVmzZhg3bhzS09MBAElJSSgrK0P//v2Fsq1bt0ZISAgSExMBAImJiYiIiBASJwCIjY1FXl4ezp49K5TRrUNbRluHISUlJcjLyxP9ayg+/5z/X0/OSgghhDRK9SZ5io6OxoYNG7Bnzx6sW7cOaWlp6NWrF+7duweVSgW5XA5PT0/RMgEBAVCpVAAAlUolSpy087XzjJXJy8tDkZHzVkuXLoWHh4fwLzg4uKabazc0Glu3gBBCCLEv9eY+TwMHDhRed+jQAdHR0QgNDcWWLVvg5ORkw5YB8+bNQ0JCgvC39iZbDYGJO9QTQgghjU696XmqzNPTEy1btsQ///wDpVKJ0tJS5OTkiMpkZmZCqVQCAJRKZZWr77R/myrj7u5uNEFTKBTCDTEb2o0xKXkihBBCxMzuedqxY4fZlQ4dOtSixlRHfn4+/v33X4wfPx5dunSBTCbDvn37MHLkSADAxYsXkZ6ejpiYGABATEwM3nnnHWRlZcHf3x8AsHfvXri7u6Nt27ZCmd27d4vWs3fvXqGOxohO2xFCCCFiZidPw4YNM6scx3Emn0ZsiVdeeQVDhgxBaGgobt68iYULF0IqlWLs2LHw8PDA5MmTkZCQAG9vb7i7u2PGjBmIiYlB9+7dAQADBgxA27ZtMX78eCxfvhwqlQpvvPEGpk+fDoVCAQB4/vnnsWbNGsyZMwfPPPMM9u/fjy1btmDXrl1W3576gnqeCCGEEDGzkyeNjbsgrl+/jrFjx+LOnTvw8/PDgw8+iKNHj8LPzw8AsHLlSkgkEowcORIlJSWIjY3Fxx9/LCwvlUqxc+dOTJs2DTExMXBxcUFcXBwWL14slAkPD8euXbswa9YsrF69Gk2bNsVnn32G2NjYOt9ee0HJEyGEECJWb+7zVJ80pPs89eoFaO9FSnsKIYSQhszc47fFV9sVFBTg4MGDSE9PR2lpqWjeiy++aGm1xM5QzxMhhBAiZlHylJycjEGDBqGwsBAFBQXw9vbG7du3hWfGUfLUcFDyRAghhIhZdKuCWbNmYciQIcjOzoaTkxOOHj2Kq1evokuXLnj//fet3UZiQ3S1HSGEECJmUfKUkpKCl19+GRKJBFKpFCUlJQgODsby5cvx2muvWbuNxIao54kQQggRsyh5kslkkEj4Rf39/YVnzHl4eODatWvWax2xOUqeCCGEEDGLxjx16tQJJ06cQIsWLfDQQw9hwYIFuH37Nv7v//4P7du3t3YbiQ3RaTtCCCFEzKKepyVLliAwMBAA8M4778DLywvTpk3DrVu38Mknn1i1gcS2qOeJEEIIEbOo5ykqKkp47e/vjz179litQcS+UPJECCGEiNXbBwM3JmlpQEmJbdZNp+0IIYQQMYt6nsLDw8FxnMH5//33n8UNImJ//QU8+CDQtStw/Hjdr596ngghhBAxi5KnmTNniv4uKytDcnIy9uzZg9mzZ1ujXeS+L77g/z9xwjbrLyureR0aDTBwIODsDGzbVvP6CCGEEFuyKHl66aWX9E5fu3YtTp48WaMGEbHyctuuv7i45nVcugT89hv/uqgIcHKqeZ2EEEKIrVh1zNPAgQPxww8/WLPKRq8hJE/37lm3PkIIIcSWrJo8ff/99/D29rZmlY2eNU6b1URRUc3ryM2teN0QkyfGbN0CQgghdcnim2TqDhhnjEGlUuHWrVv4+OOPrdY4YvueJ2sMGL9zp+K1NZIxe3L1KtCzJ/DCCwA9mYgQQhoHi5KnYcOGif6WSCTw8/NDnz590Lp1a2u0i9xn6+TJGm7frnhtD8lTURHw7bf8IPb793q12NtvAzduAK+/btvkiTE+0XWw6BNNCCHWcfEiEBQEuLnZuiW1y6Kv2oULF1q7HcQAW5+2swbd5MkeTtstXQq89RbQqhVw4ULN6pLYyZ3Sxo0D/vgDOHcO8PS0dWsIIY3R6dNAx45AdDRw9Kj5yx0/zl9Z/vbbgK9v7bXPmsz+6s/LyzP7X0Owdu1ahIWFwdHREdHR0Thui5sswbY9T9a6x5O9nbb76Sf+/4sXa16Xi0vN66iO/Hzg77+rTv/2WyAjA/i//6vb9hDS2F25Ahw5YutW1I7kZGDZMvOPBZs28f8fO1a940d0NPC//wELFlS/jbZidvLk6ekJLy8vs/7Vd5s3b0ZCQgIWLlyIU6dOoWPHjoiNjUVWVladt8WWydOuXdapx9RpuwsXgAceAN5/3zrrM8Xd3Tr1fPMNsHJlxd91MXB88GAgMhLYv1//fN3B+TWRny++SrIubN7Mx7Q2nD8P0FOkSG0ID+fHPaammi5b+UkRt28Dzz9f9QfR7dv8kyV0McYnarrHBMaALVv428FUlpICFBaaswXAkiVAQgJQWir+HuvcGZg3D6j8yNrycmD37qrf57o98Tt3Au+9V9GGu3eBH380fkw7e5Y/HqxaVVHu6lW+bceO6V9GpeLLGFJQYPwJHUeO8Mee0lLDZfThGDPvK//gwYPC6ytXrmDu3LmYOHEiYmJiAACJiYn46quvsHTpUsTFxVWvFXYmOjoaXbt2xZo1awAAGo0GwcHBmDFjBubOnWty+by8PHh4eOCrr3JRXu4ODw/+FNGmTUDz5kBODv+Gl5YCISHAmTN8d+fo0UBoKP9hLCrip0+fXlHvgQOAXM6fnsnKAsLCgFOn+J6dQYP4HbewEEhK4pf/7TdgxAhg8mQgMxOQSvnTZteuAX5+/IezXTu+HR06AI6OfM/F778DM2YA338P/Pxzxfo1Gn58z9mz/DozM/n/U1L4L46IiIo29OwJNGnCn/seMqSijg8/5JOXdu34np/Ll/meoJQUfv7cufx9oFxd+TiFhfFfGAUF/Pr69uXbUV4O/PMP8OWXwOOP8x/y7t2B9HTg4EHA25tvZ/fufOxv3+bbxnHAhAkV7dm2jZ/n5sb33DAGNG3KrzM0lJ9+8yYfy8BA/k7vGg3w33/AokXi933dOn6bWrYEYmL49+bMGT7uXbrw9WZk8Nsql/M3Db1yhY+/qysfN6WSj2NWFn/Av3aNfx0ayr9HuuOq5s8H1q8HfHzEpx8nTODjtmMHP2/YML6H7N9/+f1tyxagXz++nbdvAx4e/Ot//+W31ccHGD+er4vjgB49+O0ZNQpo356v99df+fcpJ4d/b9q35/ebZs2AoUMBhYLfritX+P2grIyP28GD/Psgk/GxGT2aX8d//wEffcSvs317YPhwvh1btgC3bgFz5vDtzMwEnniCX//p03yZ7GwgL4+fplAAcXFA69b8jWVv3+bbqa176lTg0CH+M7Z6NX/6tqSE/1w2a8Z/wbu48MmjQgE8+ij/WTl9mk++0tP5eGRm8vtxcTH/3hQU8O/vX3/xya2HBx9PiaQiFjdu8LHw9+eXa9aMj3dWFn+6VS7n3//cXKBNG35fcXTkp5eX8/tNWhq/bN++/HsWHs634+BB/te7VMrH4s4d/r1xc+Njf+8e/3r7dr6u558HgoP5WHAcv16Aryc0FLh+nT9YPfEEH4sffuDj8uijfNmcHL68kxMfMzc3vt4//wQ6deI/1z168HU1a8bvnzdv8u/5xx/zp3aeeorfJ5o35+tOTubb7u/Pv1+envxBrbiY3y6Nhm+rqyu/7o8/5g+aY8fyn5mrV/kLNwoLgYkTgV69+LYUFfExDQrit+XwYf7A/vTT/I+20lK+3f/+CwwYAERF8b0mu3fz+9CiRRVJy+nT/P4wcCD//XPoEPDuuxWfvdat+TZv3crH78MP+dNQoaHAO+8AK1bw33HBwfx6Z80Sf4c8+STwyCP8Z1ul4rctKor/Hn/uOf57yM2Nnx8VBaxdy783ANC2Lb/vLVrEt+GTT/gygwfzn/cHH+Q/D8XF/Lp9fPh969IlPk4A4OXFf56eeIJ//xISKtr24IP8ftmpE/+5/OMPvu6oKP6z37078PnnFfuGruHD+ff3yhX+fX/sMWDvXqBPH36/mDev6jK9e/Pfd9p9EACWL+ePS4MH8/OPHuU/w1Ip/2PfwYGPW/v2/Pfrc8/x0/39+ePnqVP8vtC7N/+5f+st/vsMAGbO5L/j1q7Nw6ZNHsjNzYW7kV/aZidPuvr164cpU6Zg7NixoumbNm3CJ598gj/++KO6VdqN0tJSODs74/vvvxcNjI+Li0NOTg5+0p7z0VFSUoISndQ2Ly8PwcHBAHIBWKmbww506cInZoQQQuoXDw/r9Uw3bHkATCdPFg13TUxMRFRUVJXpUVFRNhsbZC23b9+GWq1GQECAaHpAQABUKpXeZZYuXQoPDw/hH5848b+onJwMX3Xw8MN8r4mWi0tF74iPj+E2Ojvz/zdpwv+q0/L25n99DhnC9zY1a8b/cpVK+cy7aVP+A6TL1Cmsbt0qXusmTu3b870GzZrx6+zWje9x0KdVK+Pr6NIF0Ia78mBBhYKPh1ZAQMWvJGOaNKl47ePDL9etG/9LS5f2LLOnJx8bhYL/herszMesY0e+t0nL2dn0FXoymfjvLl34fcHPj//l1rkz/wtRJtN/dZy2F6AyI4+TNMrPj++hqCw6mt9erbAwfh+qqVat+J4q7YW3oaF8T2SXLlXLSqX8L8vKOnasOq1pU/HfI0eabkufPlXXZ0jLlvy+EhHB91oY4+paMd6t8j6lS6kEWrTg329vb34/qEz3PfD25vexiAhxW9u1E+/ThigU/Gc+IgJ46CH+fVAq+fdWX0y9vflyup8p7feLsW3q2pWvu00bfpqnZ8VnydA4wKCgis+Ot7f4u0VX5fdZKziYX4eLi/h7z9GR38f1jRjx9uZ7l7TfL7rff66uFdta+Tu6JletmnOVmVJp+P3s0KHqd4guHx9+//Tw4LddH33fUdVJnPR91+j7Xqocp65dDbdJ9z3TJZfzvePabQ4NrZjXogXf46UVFGSwyWZp0oTfV3RV3lZjx14RZoGWLVuy2bNnV5k+e/Zs1rJlS0uqtBs3btxgANiRI0dE02fPns26deumd5ni4mKWm5sr/Lt27RoDwHJzc5lGw5e5c4exu3cZu3mTCdP0KStjrLy84u+iIsb4k0mMtWtXMb2w0Hg9uvSVy8sTTy8s5P/pLpOdzf/Trh9g7LnnxOVMUav5/7/7jrHISMYSE/ltKiur2i7t3xoNY6Wl+uvRKi3l43TyJGOzZjF261ZFubIyw9utVjOWlMTY++8zlpnJlzEnjjk5jH39tXjbX3utIi6jRjGWn18x7/BhxnbsMK/u997j63j3XfH07GzGrlyp+LusjLFp0xjr0IExlYqfdveu+P3Zv5/fRt1pWo8+WjFNdzt27ODfF13XrlXEkTHGDhyoWDY7m5+m0TB24gTfhurQaBj75x9+We37evUqYy1bMtanD2PFxfy0vn0r1qlty/r1jCmVfHu0yw0dytjvv1fUff06Y2vX8vsZY4wlJ1fUo41b27YV07TLVhYYWFHmr7/4aZX3Qy21mrHcXL6db7zBWK9e/P6lT1kZY++8w1h8fEXcy8oYu3fPdNzKy/l/ZWX8OouK+H1TO99ceXmM3b6tfx2M8TGbNIlvY24uP+3ePf691t0vGOM/i7qfXe3r4mLGdu1ibPZs8WdDV3ExP2/37optqtweY59RbfnKZbRx0qUb35wcfr5GU7G/aevT7iOM8ftSTg5jWVnicrqxKC3l98fK32naNuTk8N9PpaX892BysrhdRUV83ffuib/3btzg36fs7Ir1aLe1srIyPo7Hj1d8PlUqxk6d4j/rKhU/78YNxvbsYezCBX7bsrMZ+/dffrvXr6/aNo2G37aCgort1u6HWVkV+15eXsUyGRmMzZ/Pr7Pyd/nNm4ydPcu/Linh69Eum5+vf59kjP+cnzsnjtnt2xXv+zff8McC7XtcVsbHVLvukycZ++03fp0lJfxnMyuror7duxnbskX7Oc4Vjt/GWHTabvfu3Rg5ciSaN2+O6OhoAMDx48dx+fJl/PDDDxg0aFB1q7Qblpy2q0w75slUt585dO/dM2iQ9QZxm+vevYreqY4dK8YmEX7MjPYZ2QkJwAcfWF5XWRn/PlvSu6S7jPbT3K8fP6g8NJQfZwDw430+/VRcrjoOHOB/CRrqGbC27Gx+zN2YMcZ7d8zx33/8L0ptz8P16/z4G19ffsyGoWU2bgReftl0bwwhpGEw9/ht0Wm7QYMG4dKlSxgyZAju3r2Lu3fvYsiQIbh06VK9TpwAQC6Xo0uXLti3b58wTaPRYN++fcLg+Lqk23VvTrd9ba6fHugrpnsKoKYHd5nM8tNy+u7r9H//x19ssHt3xbQ33+QHdX7xhWXr6du37hIngD8NM2FCzWML8KcMdN+vpk35wc+GEiftMvPnU+JECKnK4jO7wcHBWLJkiTXbYjcSEhIQFxeHqKgodOvWDatWrUJBQQEmTZpk03bV5YFLS/fSU92xGUQ8XsyWsQkLq9ojGBQE3L9YVBAYCCQm1lWrCCGk4TI7eTp9+jTat28PiUSC06dPGy3boUOHGjfMlkaPHo1bt25hwYIFUKlUiIyMxJ49e6oMIq9rtuh5ouTJMN3kyRq9I5bq0IFOpxJCSF0yO3mKjIyESqWCv78/IiMjwXEc9A2X4jgOamvdmtqG4uPjER8fb+tmAODHmVy+zN/boq7pnraj5ElM96oaWyZPS5bwY+GefNJ2bSCEkMbE7OQpLS0Nfvev8UurfOtTUquSk/mb8Pn71/26dXueDF2C2ljZy2m7Jk0qboJKCCGk9pmdPIXq3HxB9zWpfS4udf8MNS3dQczU8yRmL6ftAEqcCCGkLll0td1XX32FXTrXzM+ZMweenp7o0aMHrhp7yAyp1yh5EvPxqTh199BDtm0LIYSQumNR8rRkyRI43b9uPTExEWvWrMHy5cvh6+uLWZUf1kMajJAQW7fAvjg68s/E+vvvirtpE0IIafgsulXBtWvX0Pz+cwa2b9+OJ554AlOnTkXPnj3Rp/LzEEi916MH/9DMF16wdUvsT2SkrVtACCGkrlnU8+Tq6oo7d+4AAH777Tc88sgjAABHR0cUFRVZr3XELmifOF75uXOEEEJIY2RRz9MjjzyCKVOmoFOnTqK7ip89exZhYWHWbB+xA1Kp7QasE0IIIfbGop6ntWvXIiYmBrdu3cIPP/wAn/uPIU5KSsLYsWOt2kBCCCGEEHti0YOBiXHWfDAwIYQQQupGrT4YGAAOHTqEp59+Gj169MCNGzcAAP/3f/+Hw4cPW1olIYQQQojdsyh5+uGHHxAbGwsnJyecOnUKJSUlAIDc3NwG+7BgQgghhBDAwuTp7bffxvr16/Hpp59CJpMJ03v27IlTp05ZrXGEEEIIIfbGouTp4sWL6N27d5XpHh4eyMnJqWmbCCGEEELslkXJk1KpxD///FNl+uHDh9GsWbMaN4oQQgghxF5ZlDw9++yzeOmll3Ds2DFwHIebN29i48aNePnllzFt2jRrt5EQQgghxG5YdJPMuXPnQqPRoF+/figsLETv3r2hUCgwe/ZsTJkyxdptJIQQQgixGxb1PHEch9dffx13797FmTNncPToUdy6dQseHh4IDw+3dhsJIYQQQuxGtZKnkpISzJs3D1FRUejZsyd2796Ntm3b4uzZs2jVqhVWr16NWbNm1UpDw8LCwHGc6N+yZctEZU6fPo1evXrB0dERwcHBWL58eZV6tm7ditatW8PR0RERERHYvXu3aD5jDAsWLEBgYCCcnJzQv39/XL58uVa2iRBCCCH1T7WSpwULFmDdunUICwtDWloaRo0ahalTp2LlypX44IMPkJaWhldffbW22orFixcjIyND+DdjxgxhXl5eHgYMGIDQ0FAkJSXhvffew6JFi/DJJ58IZY4cOYKxY8di8uTJSE5OxrBhwzBs2DCcOXNGKLN8+XJ8+OGHWL9+PY4dOwYXFxfExsaiuLi41raLEEIIIfUIq4bw8HD2008/McYYS01NZRzHsUmTJjGNRlOdaiwSGhrKVq5caXD+xx9/zLy8vFhJSYkw7dVXX2WtWrUS/n7yySfZ4MGDRctFR0ez5557jjHGmEajYUqlkr333nvC/JycHKZQKNi3335rdltzc3MZAJabm2v2MoQQQgixLXOP39Xqebp+/Tq6dOkCAGjfvj0UCgVmzZoFjuNqIa2ratmyZfDx8UGnTp3w3nvvoby8XJiXmJiI3r17Qy6XC9NiY2Nx8eJFZGdnC2X69+8vqjM2NhaJiYkAgLS0NKhUKlEZDw8PREdHC2X0KSkpQV5enugfIYQQQhqmal1tp1arRcmJg4MDXF1drd4ofV588UV07twZ3t7eOHLkCObNm4eMjAysWLECAKBSqaoMVg8ICBDmeXl5QaVSCdN0y6hUKqGc7nL6yuizdOlSvPnmmzXbQEIIIYTUC9VKnhhjmDhxIhQKBQCguLgYzz//PFxcXETlfvzxR7Pqmzt3Lt59912jZc6fP4/WrVsjISFBmNahQwfI5XI899xzWLp0qdAeW5k3b56ofXl5eQgODrZhiwghhBBSW6qVPMXFxYn+fvrpp2u08pdffhkTJ040WsbQHcujo6NRXl6OK1euoFWrVlAqlcjMzBSV0f6tVCqF//WV0Z2vnRYYGCgqExkZabCNCoXC5gkcIYQQQupGtZKnL7/80qor9/Pzg5+fn0XLpqSkQCKRwN/fHwAQExOD119/HWVlZcLDivfu3YtWrVrBy8tLKLNv3z7MnDlTqGfv3r2IiYkBAISHh0OpVGLfvn1CspSXl4djx47RndMJIYQQAsDCm2TWtcTERKxatQp///03/vvvP2zcuBGzZs3C008/LSRGTz31FORyOSZPnoyzZ89i8+bNWL16teh02ksvvYQ9e/bggw8+wIULF7Bo0SKcPHkS8fHxAPibf86cORNvv/02duzYgdTUVEyYMAFBQUEYNmyYLTadEEIIIfambi7+q5mkpCQWHR3NPDw8mKOjI2vTpg1bsmQJKy4uFpX7+++/2YMPPsgUCgVr0qQJW7ZsWZW6tmzZwlq2bMnkcjlr164d27Vrl2i+RqNh8+fPZwEBAUyhULB+/fqxixcvVqu9dKsCQgghpP4x9/jNMcaYrRO4hiYvLw8eHh7Izc2Fu7u7rZtDCCGEEDOYe/yuF6ftCCGEEELsRbUGjBPzaDvz6GaZhBBCSP2hPW6bOilHyVMtuHfvHgDQvZ4IIYSQeujevXvw8PAwOJ/GPNUCjUaDmzdvws3NrUaPrtHebPPatWs0dsoIipNpFCPzUJxMoxiZh+Jkmj3GiDGGe/fuISgoCBKJ4ZFN1PNUCyQSCZo2bWq1+tzd3e1mx7JnFCfTKEbmoTiZRjEyD8XJNHuLkbEeJy0aME4IIYQQUg2UPBFCCCGEVAMlT3ZMoVBg4cKF9Nw8EyhOplGMzENxMo1iZB6Kk2n1OUY0YJwQQgghpBqo54kQQgghpBooeSKEEEIIqQZKngghhBBCqoGSJ0IIIYSQaqDkyQ4sXboUXbt2hZubG/z9/TFs2DBcvHhRVKa4uBjTp0+Hj48PXF1dMXLkSGRmZtqoxba3bNkycByHmTNnCtMoRsCNGzfw9NNPw8fHB05OToiIiMDJkyeF+YwxLFiwAIGBgXByckL//v1x+fJlG7a47qnVasyfPx/h4eFwcnLCAw88gLfeekv0LKvGGKc///wTQ4YMQVBQEDiOw/bt20XzzYnJ3bt3MW7cOLi7u8PT0xOTJ09Gfn5+HW5F7TIWo7KyMrz66quIiIiAi4sLgoKCMGHCBNy8eVNUR0OPEWB6X9L1/PPPg+M4rFq1SjTd3uNEyZMdOHjwIKZPn46jR49i7969KCsrw4ABA1BQUCCUmTVrFn7++Wds3boVBw8exM2bNzFixAgbttp2Tpw4gf/973/o0KGDaHpjj1F2djZ69uwJmUyGX375BefOncMHH3wALy8voczy5cvx4YcfYv369Th27BhcXFwQGxuL4uJiG7a8br377rtYt24d1qxZg/Pnz+Pdd9/F8uXL8dFHHwllGmOcCgoK0LFjR6xdu1bvfHNiMm7cOJw9exZ79+7Fzp078eeff2Lq1Kl1tQm1zliMCgsLcerUKcyfPx+nTp3Cjz/+iIsXL2Lo0KGicg09RoDpfUlr27ZtOHr0KIKCgqrMs/s4MWJ3srKyGAB28OBBxhhjOTk5TCaTsa1btwplzp8/zwCwxMREWzXTJu7du8datGjB9u7dyx566CH20ksvMcYoRowx9uqrr7IHH3zQ4HyNRsOUSiV77733hGk5OTlMoVCwb7/9ti6aaBcGDx7MnnnmGdG0ESNGsHHjxjHGKE6MMQaAbdu2TfjbnJicO3eOAWAnTpwQyvzyyy+M4zh248aNOmt7XakcI32OHz/OALCrV68yxhpfjBgzHKfr16+zJk2asDNnzrDQ0FC2cuVKYV59iBP1PNmh3NxcAIC3tzcAICkpCWVlZejfv79QpnXr1ggJCUFiYqJN2mgr06dPx+DBg0WxAChGALBjxw5ERUVh1KhR8Pf3R6dOnfDpp58K89PS0qBSqUQx8vDwQHR0dKOJEQD06NED+/btw6VLlwAAf//9Nw4fPoyBAwcCoDjpY05MEhMT4enpiaioKKFM//79IZFIcOzYsTpvsz3Izc0Fx3Hw9PQEQDHS0mg0GD9+PGbPno127dpVmV8f4kQPBrYzGo0GM2fORM+ePdG+fXsAgEqlglwuFz6AWgEBAVCpVDZopW189913OHXqFE6cOFFlHsUI+O+//7Bu3TokJCTgtddew4kTJ/Diiy9CLpcjLi5OiENAQIBoucYUIwCYO3cu8vLy0Lp1a0ilUqjVarzzzjsYN24cAFCc9DAnJiqVCv7+/qL5Dg4O8Pb2bpRxKy4uxquvvoqxY8cKD72lGPHeffddODg44MUXX9Q7vz7EiZInOzN9+nScOXMGhw8ftnVT7Mq1a9fw0ksvYe/evXB0dLR1c+ySRqNBVFQUlixZAgDo1KkTzpw5g/Xr1yMuLs7GrbMfW7ZswcaNG7Fp0ya0a9cOKSkpmDlzJoKCgihOxCrKysrw5JNPgjGGdevW2bo5diUpKQmrV6/GqVOnwHGcrZtjMTptZ0fi4+Oxc+dOHDhwAE2bNhWmK5VKlJaWIicnR1Q+MzMTSqWyjltpG0lJScjKykLnzp3h4OAABwcHHDx4EB9++CEcHBwQEBDQ6GMUGBiItm3biqa1adMG6enpACDEofIViI0pRgAwe/ZszJ07F2PGjEFERATGjx+PWbNmYenSpQAoTvqYExOlUomsrCzR/PLycty9e7dRxU2bOF29ehV79+4Vep0AihEAHDp0CFlZWQgJCRG+y69evYqXX34ZYWFhAOpHnCh5sgOMMcTHx2Pbtm3Yv38/wsPDRfO7dOkCmUyGffv2CdMuXryI9PR0xMTE1HVzbaJfv35ITU1FSkqK8C8qKgrjxo0TXjf2GPXs2bPKLS4uXbqE0NBQAEB4eDiUSqUoRnl5eTh27FijiRHAXxUlkYi/+qRSKTQaDQCKkz7mxCQmJgY5OTlISkoSyuzfvx8ajQbR0dF13mZb0CZOly9fxu+//w4fHx/RfIoRMH78eJw+fVr0XR4UFITZs2fj119/BVBP4mTrEeuEsWnTpjEPDw/2xx9/sIyMDOFfYWGhUOb5559nISEhbP/+/ezkyZMsJiaGxcTE2LDVtqd7tR1jFKPjx48zBwcH9s4777DLly+zjRs3MmdnZ/bNN98IZZYtW8Y8PT3ZTz/9xE6fPs0ef/xxFh4ezoqKimzY8roVFxfHmjRpwnbu3MnS0tLYjz/+yHx9fdmcOXOEMo0xTvfu3WPJycksOTmZAWArVqxgycnJwpVi5sTk0UcfZZ06dWLHjh1jhw8fZi1atGBjx4611SZZnbEYlZaWsqFDh7KmTZuylJQU0Xd5SUmJUEdDjxFjpvelyipfbceY/ceJkic7AEDvvy+//FIoU1RUxF544QXm5eXFnJ2d2fDhw1lGRobtGm0HKidPFCPGfv75Z9a+fXumUChY69at2SeffCKar9Fo2Pz581lAQABTKBSsX79+7OLFizZqrW3k5eWxl156iYWEhDBHR0fWrFkz9vrrr4sOcI0xTgcOHND7PRQXF8cYMy8md+7cYWPHjmWurq7M3d2dTZo0id27d88GW1M7jMUoLS3N4Hf5gQMHhDoaeowYM70vVaYvebL3OHGM6dxWlxBCCCGEGEVjngghhBBCqoGSJ0IIIYSQaqDkiRBCCCGkGih5IoQQQgipBkqeCCGEEEKqgZInQgghhJBqoOSJEEIIIaQaKHkihBBCCKkGSp4IIYQQQqqBkidCCCGEkGqg5IkQ0qAxxrBixQqEh4fD2dkZw4YNQ25urtFl7ty5A39/f1y5csVouT59+mDmzJnWa6weY8aMwQcffFCr6yCEVA8lT4SQemvp0qXo2rUr3Nzc4O/vj2HDhuHixYuiMrNnz8a6devw1Vdf4dChQ0hKSsKiRYuM1vvOO+/g8ccfR1hYWO013kxvvPEG3nnnHZMJHyGk7lDyRAiptw4ePIjp06fj6NGj2Lt3L8rKyjBgwAAUFBQAAI4dO4YVK1Zg8+bN6N27N7p06YJnn30Wu3fvNlhnYWEhPv/8c0yePLmuNsOo9u3b44EHHsA333xj66YQQu6j5IkQUm/t2bMHEydORLt27dCxY0ds2LAB6enpSEpKAgC8//776NevHzp37iwsExAQgNu3bxusc/fu3VAoFOjevbtoekFBASZMmABXV1cEBgZWOZWm0WiwdOlShIeHw8np/9u5n5Co1jCO41/lTiY2SdNk+CeiEkJLwoQJGaN/9odMItwEBS6qlS6ymiAigqIMIlsUTC2ijJKKYLKajf0hm5HKIFtMUZqBQeAiGk5azoygdyH3cM+tqc5VC+P3gbOY57znOS/v6pnnfWfSWbRoEdevX7eM6evrY8uWLWRkZJCdnc3Jkyd/auuvsrKSK1eu/MySiMgvoOJJRP4Y/2xtuVwu4vE4wWCQTZs2WcbEYjEyMzOT5giFQpSUlHwV9/l8tLa20tzcTEtLCw8ePODZs2fm/fr6ei5evMiZM2d48eIFdXV1bN26ldbWVnPMrl27aGtr4+bNm9y5c4dQKGTJkYzH46G9vZ14PP7DsSIy/v763RMQERkLQ0ND7Ny5E6/Xy8KFC3n06BEDAwPs3r2bvXv3muMGBwdZsWJF0jw9PT3k5ORYYv39/Zw7d45Lly6xatUqABobG8nLywMgHo9z9OhR7t69S2lpKQBz584lHA5z9uxZli1bRl9fH42NjTQ1NZk5zp8//9W7viUnJ4dEIkFvby+zZ8+2tzAiMuZUPInIH6GmpoZIJEI4HAags7OTjIwMnj9/bhlXUVGB1+tNmmdgYIDJkydbYt3d3SQSCZYsWWLGXC4X8+fPB+DNmzd8+fKF1atXW55LJBIUFxcD8PbtWwYHB/F4POb9zMxMM8f3pKenAyPnsUTk91PxJCITXm1tLbdv3+bhw4dmN+jTp0+43W7y8/PNcT09PXR1dVFVVZU0l9vtJhqN2np/f38/AMFgkNzcXMu9tLQ0W7m+5ePHjwDMmDFj1LlEZPR05klEJqzh4WFqa2sJBALcv3+fOXPmmPfcbjeGYTA8PGzGjhw5wvr16yksLEyas7i4mJcvX1pi8+bNw+Fw8OTJEzMWjUbp7OwEoLCwkLS0NN69e0d+fr7lmjVrFjCyjedwOHj69KmZwzAMM8f3RCIR8vLycLvdPxwrIuNPnScRmbBqampoamqiubkZp9NJb28vMLIdtnLlSmKxGMeOHWPz5s1cvnyZW7du0d7e/t2ca9euZd++fUSjUaZNmwbAlClT2LZtGz6fj+nTp5OVlcX+/ftJTR35/ul0OtmzZw91dXUMDQ1RVlaGYRi0tbUxdepUqqurcTqdVFdX4/P5cLlcZGVlcfDgQVJTU0lJSTHff/r0aQKBAPfu3TNjoVCINWvWjPXyicj/pM6TiExYfr8fwzBYvnw52dnZ5nX16lVmzpzJhQsX8Pv9LFiwgMePHxMOh81OUDJFRUUsXryYa9euWeLHjx9n6dKlVFZWUl5eTllZmeVXeYcPH+bAgQPU19dTUFDAunXrCAaDlm5YQ0MDpaWlbNiwgfLycrxeLwUFBZYzVh8+fKC7u9v8HIvFuHHjBjt27BjtconIGEkZ/ndPW0RECAaD+Hw+IpGI2V0aD58/fyY3N5cTJ04k/VNOv99PIBCgpaVl3OYhIvZo205E5D8qKiro6uri/fv3P+xU2dHR0cGrV6/weDwYhsGhQ4cA2LhxY9JnHA4Hp06dGrM5iMjoqfMkIvKLdHR0sH37dl6/fs2kSZMoKSmhoaGBoqKi3z01EbFBxZOIiIiIDTowLiIiImKDiicRERERG1Q8iYiIiNig4klERETEBhVPIiIiIjaoeBIRERGxQcWTiIiIiA0qnkRERERsUPEkIiIiYoOKJxEREREb/gZtXT8GTl0DUQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Rietveld plot\n",
    "def rietveld_plot():\n",
    "    import GSASIIscriptable as G2sc\n",
    "\n",
    "    gpx = G2sc.G2Project(\n",
    "        f'{WORK_DIR}/project_seed{RANDOM_SEED}_trial_{study.best_trial.number}.gpx')\n",
    "\n",
    "    hist1 = gpx.histograms()[0]\n",
    "    phase0 = gpx.phases()[0]\n",
    "\n",
    "    hist = hist1\n",
    "    i = 5\n",
    "    two_theta = hist.getdata(\"X\")[::i]\n",
    "    Yobs = hist.getdata(\"Yobs\")[::i]\n",
    "    Ycalc = hist.getdata(\"Ycalc\")[::i]\n",
    "    bg = hist.getdata(\"Background\")[::i]\n",
    "    residual = hist.getdata(\"Residual\")[::i]\n",
    "\n",
    "    fig = plt.figure()\n",
    "    gs = GridSpec(5, 1, figure=fig)\n",
    "    ax1 = fig.add_subplot(gs[:4, :])\n",
    "    ax2 = fig.add_subplot(gs[4, :])\n",
    "    fig.subplots_adjust(hspace=0)\n",
    "    ax1.grid(color='#cccccc')\n",
    "\n",
    "    ax1.scatter(two_theta, Yobs, marker='P', lw=0.0001, c='Black', label='XRD (Obs)')\n",
    "    ax1.plot(two_theta, Ycalc, label='XRD (Calc)')\n",
    "    ax1.plot(two_theta, bg, color='red', label='Background (Calc)')\n",
    "    ax1.set_ylabel('Intensity')\n",
    "    ax1.legend()\n",
    "    ax2.plot(two_theta, residual, color='blue')\n",
    "    plt.setp(ax1.get_xticklabels(), visible=False);\n",
    "    ax2.set_ylim(-6600, 6600)\n",
    "    plt.xlabel(r'$2\\theta$ (deg.)')\n",
    "    ax2.set_ylabel('Residual')\n",
    "    ax1.set_xlim(15, 150)\n",
    "    ax2.set_xlim(15, 150)\n",
    "    plt.show()\n",
    "    \n",
    "rietveld_plot()"
   ]
  },
  {
   "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.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}