{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.6.10"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,
"autoclose": false,
"autocomplete": true,
"bibliofile": "biblio.bib",
"cite_by": "apalike",
"current_citInitial": 1,
"eqLabelWithNumbers": true,
"eqNumInitial": 1,
"hotkeys": {
"equation": "Ctrl-E",
"itemize": "Ctrl-I"
},
"labels_anchors": false,
"latex_user_defs": false,
"report_style_numbering": false,
"user_envs_cfg": false
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {
"height": "672px",
"left": "0px",
"right": "1324px",
"top": "107px",
"width": "276px"
},
"toc_section_display": "block",
"toc_window_display": true
},
"varInspector": {
"cols": {
"lenName": 16,
"lenType": 16,
"lenVar": 40
},
"kernels_config": {
"python": {
"delete_cmd_postfix": "",
"delete_cmd_prefix": "del ",
"library": "var_list.py",
"varRefreshCmd": "print(var_dic_list())"
},
"r": {
"delete_cmd_postfix": ") ",
"delete_cmd_prefix": "rm(",
"library": "var_list.r",
"varRefreshCmd": "cat(var_dic_list()) "
}
},
"oldHeight": 383.4,
"position": {
"height": "40px",
"left": "839px",
"right": "20px",
"top": "-2px",
"width": "550px"
},
"types_to_exclude": [
"module",
"function",
"builtin_function_or_method",
"instance",
"_Feature"
],
"varInspector_section_display": "none",
"window_display": true
},
"colab": {
"name": "StructuralModeling.ipynb",
"provenance": [],
"include_colab_link": true
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"![\"Open](\"https://colab.research.google.com/assets/colab-badge.svg\")
"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "vIs1VkuJgaGu",
"colab_type": "text"
},
"source": [
"# Structural Modeling\n",
"\n",
"The following notebook is comprised of 9 primary steps:\n",
"\n",
"1. Initialize required packages, directories and parameters \n",
"2. Load and inspect the surface data\n",
"3. Decluster and normal score transform\n",
"4. Variogram calculation and modeling\n",
"5. Simulation and back-transformation\n",
"6. Simulation checking \n",
"7. Base surface calculation\n",
"8. 3-D grid and keyout calculation\n",
"9. Save project setting and clean the output files"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "oqZj66iGgaGx",
"colab_type": "text"
},
"source": [
"## 1. Initialize required packages and parameters"
]
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:30.989488Z",
"start_time": "2020-06-03T02:32:28.193490Z"
},
"id": "p06M01P7gaGz",
"colab_type": "code",
"colab": {},
"outputId": "d8495efa-0e08-4138-9ee2-a6494d1b1034"
},
"source": [
"import pygeostat as gs\n",
"import numpy as np\n",
"import os\n",
"import pandas as pd\n",
"import matplotlib.pyplot as plt"
],
"execution_count": null,
"outputs": [
{
"output_type": "stream",
"text": [
"Loading default Pygeostat Parameters from C:\\Users\\MostafaHadavand\\.Pygeostat\\Parameters.json\n"
],
"name": "stdout"
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "caLLsYFLgaHG",
"colab_type": "text"
},
"source": [
"### Load project settings and modify where required\n",
"Load the Matplotlib and Pygeostat project defaults."
]
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:31.014489Z",
"start_time": "2020-06-03T02:32:30.991491Z"
},
"id": "O1uEj6VIgaHP",
"colab_type": "code",
"colab": {}
},
"source": [
"gs.Parameters.load('Parameters.json')\n",
"gs.PlotStyle.load('Parameters.json')"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "Tp8pZ_nEgaHp",
"colab_type": "text"
},
"source": [
"Assign the number of realization and griddef to variables for convenience."
]
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:31.020491Z",
"start_time": "2020-06-03T02:32:31.016493Z"
},
"id": "tmmLiRXqgaHq",
"colab_type": "code",
"colab": {}
},
"source": [
"# a path to executables can be defined\n",
"exedir='../pygeostat/executable/'"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:31.033493Z",
"start_time": "2020-06-03T02:32:31.025493Z"
},
"id": "XzOaTn_hgaHw",
"colab_type": "code",
"colab": {},
"outputId": "744b4943-4e7d-4bb3-e2aa-1ef891c69019"
},
"source": [
"# nreal = gs.Parameters['data.nreal']\n",
"nreal=100\n",
"gs.Parameters['data.nreal']=nreal\n",
"print('nreal = {}'.format(nreal))\n",
"griddef = gs.Parameters['data.griddef']\n",
"print(griddef)"
],
"execution_count": null,
"outputs": [
{
"output_type": "stream",
"text": [
"nreal = 100\n",
"120 5.0 10.0 \n",
"110 1205.0 10.0 \n",
"1 0.5 1.0\n"
],
"name": "stdout"
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "6-uJdf89gaH1",
"colab_type": "text"
},
"source": [
"The cmap was set to bwr in the previous section - reset to the \n",
"Matplotlib default"
]
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:31.053491Z",
"start_time": "2020-06-03T02:32:31.038493Z"
},
"id": "NSBmKW_7gaH2",
"colab_type": "code",
"colab": {}
},
"source": [
"gs.Parameters['plotting.cmap'] = 'viridis'"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "LCZpZW0GgaID",
"colab_type": "text"
},
"source": [
"### Directories and files"
]
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:31.065491Z",
"start_time": "2020-06-03T02:32:31.060504Z"
},
"id": "T8Xyb1l4gaIF",
"colab_type": "code",
"colab": {}
},
"source": [
"# create the output directory\n",
"outdir = 'Output/'\n",
"gs.mkdir(outdir)"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"metadata": {
"id": "LZBOUGATgaIU",
"colab_type": "text"
},
"source": [
"## 2. Load and inspect the surface data\n",
"\n",
"### Load the data and note its attributes"
]
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:31.116509Z",
"start_time": "2020-06-03T02:32:31.069491Z"
},
"scrolled": true,
"id": "YDOp27_1gaIV",
"colab_type": "code",
"colab": {},
"outputId": "1f7b998d-13df-4438-dd16-ff23cf1fe660"
},
"source": [
"dat = gs.ExampleData('reservoir_surface')\n",
"print('columns = {}'.format(dat.columns.values))\n",
"print('x column = {}, y column = {}'.format(dat.x, dat.y))"
],
"execution_count": null,
"outputs": [
{
"output_type": "stream",
"text": [
"columns = ['HoleID' 'X' 'Y' 'Top Elevation' 'Thickness' 'Base Elevation']\n",
"x column = X, y column = Y\n"
],
"name": "stdout"
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "OMs1NK-PgaIe",
"colab_type": "text"
},
"source": [
"### Summary statistics\n",
"The DataFile describe excludes special attributes."
]
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:31.158486Z",
"start_time": "2020-06-03T02:32:31.119489Z"
},
"id": "HzAVXjoIgaIf",
"colab_type": "code",
"colab": {},
"outputId": "0c86c59f-7144-491d-cfe2-500a2c89f409"
},
"source": [
"dat.describe()"
],
"execution_count": null,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" Top Elevation | \n",
" Thickness | \n",
" Base Elevation | \n",
"
\n",
" \n",
" \n",
" \n",
" | count | \n",
" 230.000000 | \n",
" 230.000000 | \n",
" 230.000000 | \n",
"
\n",
" \n",
" | mean | \n",
" 379.173739 | \n",
" 50.096391 | \n",
" 329.077348 | \n",
"
\n",
" \n",
" | std | \n",
" 2.607090 | \n",
" 4.376842 | \n",
" 4.571519 | \n",
"
\n",
" \n",
" | min | \n",
" 372.070000 | \n",
" 37.730000 | \n",
" 315.850000 | \n",
"
\n",
" \n",
" | 25% | \n",
" 377.222500 | \n",
" 47.502500 | \n",
" 327.137500 | \n",
"
\n",
" \n",
" | 50% | \n",
" 378.950000 | \n",
" 49.390000 | \n",
" 329.140000 | \n",
"
\n",
" \n",
" | 75% | \n",
" 380.782500 | \n",
" 51.900000 | \n",
" 332.225000 | \n",
"
\n",
" \n",
" | max | \n",
" 386.330000 | \n",
" 62.870000 | \n",
" 340.870000 | \n",
"
\n",
" \n",
"
\n",
"
"
]
},
"metadata": {
"tags": []
}
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "nmHIHQBUgaI9",
"colab_type": "text"
},
"source": [
"## 3. Decluster and Normal Score Transform\n",
"\n",
"Declustering defines the spatially representative distribution of the Top Elevation and Thickness. These distributions are then used for the normal score transformation, generating conditioning data that is used for calculating normal score variograms and conditioning sequential Gaussian simulation.\n",
"\n",
"Note that variograms are often calculated on data that are normal score transformed without declustering weights. Given the controversy of the subject, the more convenient method is selected for this course (since declustering with weights is absolutely required for conditioning).\n",
"\n",
"### Decluster the data\n",
"All variables are homotopically sampled, so only one variable need be considered for declustering.\n"
]
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:32.714487Z",
"start_time": "2020-06-03T02:32:32.394512Z"
},
"id": "l1rStNcagaI-",
"colab_type": "code",
"colab": {}
},
"source": [
"declus = gs.Program(exedir+'declus')"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:32.725489Z",
"start_time": "2020-06-03T02:32:32.716487Z"
},
"id": "3YdcmwvmgaJH",
"colab_type": "code",
"colab": {},
"outputId": "d7a67fea-6fd0-4cdf-a4e4-f453392c3753"
},
"source": [
"os.path.join(outdir, 'declus.out')"
],
"execution_count": null,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
"'Output/declus.out'"
]
},
"metadata": {
"tags": []
},
"execution_count": 11
}
]
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:32.830494Z",
"start_time": "2020-06-03T02:32:32.730506Z"
},
"id": "fsuR5_iKgaJQ",
"colab_type": "code",
"colab": {},
"outputId": "3807cfe8-606a-475a-93fa-35306cbff2fd"
},
"source": [
"# Cell size based on the data spacing study in the intro notebook\n",
"cellsize = 90\n",
"parstr = \"\"\" Parameters for DECLUS\n",
" *********************\n",
"\n",
"START OF PARAMETERS:\n",
"{datafl} -file with data\n",
"{xyzcol} {varcol} - columns for X, Y, Z, and variable\n",
"-1.0e21 1.0e21 - trimming limits\n",
"junk.out -file for summary output\n",
"{outfl} -file for output with data & weights\n",
"1.0 1.0 -Y and Z cell anisotropy (Ysize=size*Yanis)\n",
"0 -0=look for minimum declustered mean (1=max)\n",
"1 {cellsize} {cellsize} -number of cell sizes, min size, max size\n",
"5 -number of origin offsets\n",
"\"\"\"\n",
"pars = dict(datafl=dat.flname, xyzcol=dat.gscol(dat.xyz),\n",
" varcol=dat.gscol('Top Elevation'),\n",
" outfl=os.path.join(outdir, 'declus.out'), cellsize=cellsize)\n",
"declus.run(parstr=parstr.format(**pars))\n",
"gs.rmfile('junk.out')"
],
"execution_count": null,
"outputs": [
{
"output_type": "stream",
"text": [
"Calling: ['../pygeostat/executable/declus', 'temp']\n",
"\n",
"DECLUS Version: 3.000\n",
"\n",
" data file = c:\\users\\mostafahadavand\\anaconda3\\envs\\\n",
" columns = 2 3 0 4\n",
" tmin,tmax = -1.000000E+21 1.000000E+21\n",
" summary file = junk.out \n",
" output file = Output/declus.out \n",
" anisotropy = 1.000000 1.000000\n",
" minmax flag = 0\n",
" ncell min max = 1 90.000000 90.000000\n",
" offsets = 5\n",
"\n",
"\n",
"There are 230 data with:\n",
" mean value = 379.17374\n",
" minimum and maximum = 372.07001 386.32999\n",
" size of data vol in X = 970.03998\n",
" size of data vol in Y = 1080.00000\n",
" size of data vol in Z = 0.00000\n",
"\n",
" declustered mean = 379.24036\n",
" min and max weight = 0.42091 2.19493\n",
" equal weighting = 1.00000\n",
"\n",
"\n",
"DECLUS Version: 3.000 Finished\n",
"\n",
"Stop - Program terminated.\n"
],
"name": "stdout"
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "qeVJ1F_8gaJY",
"colab_type": "text"
},
"source": [
"### Load and inspect the declustering results\n",
"\n",
"In the context of the upcoming modeling steps, Base Elevation is not considered a variable. Use of the notvariables kwarg on initialization excludes it from the variables attribute."
]
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:32.886490Z",
"start_time": "2020-06-03T02:32:32.833490Z"
},
"scrolled": false,
"id": "_n1zQZWVgaJZ",
"colab_type": "code",
"colab": {},
"outputId": "79c09d8f-7781-483c-fde5-4852be91b0cc"
},
"source": [
"dat = gs.DataFile('Output/declus.out', notvariables='Base Elevation')\n",
"print('declustering weight = ', dat.weights)\n",
"print('variables = ', dat.variables)"
],
"execution_count": null,
"outputs": [
{
"output_type": "stream",
"text": [
"declustering weight = Declustering Weight\n",
"variables = ['Top Elevation', 'Thickness']\n"
],
"name": "stdout"
}
]
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:33.384503Z",
"start_time": "2020-06-03T02:32:32.889492Z"
},
"id": "e_XKoPGogaJi",
"colab_type": "code",
"colab": {},
"outputId": "722b6aed-7546-4711-cae4-67edd4246096"
},
"source": [
"gs.location_plot(dat, var=dat.weights)"
],
"execution_count": null,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/plain": [
""
]
},
"metadata": {
"tags": []
},
"execution_count": 14
},
{
"output_type": "display_data",
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAeEAAAHqCAYAAAAzjfbNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi41LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvSM8oowAAIABJREFUeJzs3Xd4VGX68PHvlLRJb6QQWkIJHQIkhBaQjhQVEQUsCBZQUVH5ubqrq64vu4iN1UXXiiy6CisqRUU6EpAmKL0lgZCQhPQ2mXbePwIDkxAJIVOS3B+vc10+Z845zz0wzD1POc9RKYqiIIQQQgiHUzs7ACGEEKKpkiQshBBCOIkkYSGEEMJJJAkLIYQQTiJJWAghhHASScJCCCGEk0gSFqIedOjQgdTUVGeH4XKWLVvG/fff7+wwhHBZkoRFk+bt7U3btm1p06YN/fr1Y8+ePXW6TkVFBSaT6brP+9e//kVRUdF1n/fAAw+wb9++6z7vSvHx8Wzfvt1afvrpp5k4caK1XFZWRkhIyB++r2vFYTQaMRgMV32tru9diMZEkrBo0srKyjh69CgpKSm89NJL3HzzzQ5t0S5YsIDs7OzrPu+DDz4gLi7uhuoeOHAg69evt5Y3bdrE5s2bMZvNAPz888/07NkTrVZrlzjq+t6FaEwkCQtx0fDhw5kxYwbvvPMOANu3b6dPnz60a9eOiRMnWlttGzdupGvXrrRq1Yq5c+faXOPs2bPExMRYy9OnT2fZsmUAzJs3j5iYGDp16kRqaiqTJ0/m3LlzDB06lKeeegqAb775hu7du9OuXTtmzpxpbYU2b96chx56iBYtWnD+/HmGDx/O5s2bAYiKiuL555+nY8eODB06lOLiYgCWL19ObGwsgYGBeHh48PTTT9vEOmLECDZu3AhAYWEhJpOJ2NhYdu3aBVQm5REjRgBw5MgRBg0aRIcOHRg2bBjnzp2z/pldiuNPf/oTMTExeHp6EhAQwPLlywHIyspi3LhxtGzZkgULFgBc9b0L0SQpQjRhgGI0Gq3lb7/9VhkzZoyi1+uVbt26KRkZGYqiKMojjzyizJ8/XykvL1ciIyOVgwcPKhaLRfnxxx8VRVGUVq1aKSdOnFBSUlKU5s2bW683depU5ZNPPlFyc3OVyMhIxWQyKQUFBUpxcbHNeYqiKFlZWUpcXJxSWFiomM1mZdy4ccoXX3xhjfPtt9+2XjcpKUn56aefrK999NFHiqIoyvDhw63/Hx0drWRnZyvl5eVKu3btlIKCApv3XlZWpvj7+ytlZWXK6tWrldmzZysvvPCC8sorryiKoih9+/ZVDhw4oCiKoiQkJCiHDh1SFEVRXnvtNeWhhx6yiePIkSNKYmKiYjKZlCNHjih9+vRRFEVRPvnkEyUqKkpJS0tTUlNTFZ1Op5SXl1d770I0VdISFuIKZrMZDw8Pjh07xokTJxgyZAixsbGsWbOG/Px8Tpw4QUhICJ07d0alUllbitcSEBBAdHQ0TzzxBLm5ufj4+FQ7ZufOnZw8eZL4+Hg6derE/v37uXDhgvX12bNn13j9e+65B4BevXpZW6kajQaTyYTBYLC+ryt5eXnRp08ftm/fzpYtW0hKSmLYsGGsX7+ekpISzp07R7du3SgqKmLv3r3cdtttxMbGsnjxYgoLC22updFosFgsmM1mSkpK0Ol01teSkpJo2bIlrVq1ws/Pj5ycnFr9mQnRFNQ82CNEE7R3715iY2OxWCy0a9eOAwcO2Ly+b98+VCpVjedrNBrrmOqV1Go1mzdv5n//+x8jR47k008/pX///jbHWCwWBg8ezLffflvn+C/VX1FRgb+/P+PHj8doNPLKK6/g6elZ7fhLXdLbtm1j7ty5BAcHc/DgQdatW8eQIUMAUBQFb29vjh49WmO9er2eiooKevbsSVBQEIsWLfrD+IQQlaQlLMRFW7du5eOPP2b69Om0b9+ezMxMfvnlFwCKi4vJycmhQ4cOpKenc+TIERRFYeXKlTbXCA0NpaioiPT0dFJSUqxjrhcuXCArK4s77riDkSNHWmdh63Q6a8uwd+/ebN++ndOnT1vPqevsYbVaTUlJCbm5uZSUlPD++++TkpJS7bgRI0awevVqSkpKCA8Px83NjYSEBF555RVGjhwJgL+/PxEREaxYsQKonAl+qbV9iZeXFydOnMBsNnPu3DkWLlxY46zoS65870I0VZKERZOm0+mIjY0lKiqKF198ke+//5527dqh0+lYsmQJ06dPJyYmhqSkJNLS0vD29ub9999n3LhxtGvXjuTkZAC0Wi0ajQZPT09eeOEFEhMTeeSRRxgyZAharZaysjLGjBlD27ZtOXToENOmTQNg1qxZ3HbbbcydO5eoqChef/11RowYQdu2bRk7dix5eXnW61/ZAtdoNGg0mmqvXdp//PhxfH19OX78OCdOnGDAgAH885//rPb+u3fvTkFBAQMHDrTuGz16NIcPH2b48OHWfZ999hnz588nOjqa3r17s3//fpv6fvrpJ2bOnMnRo0c5evQoGRkZ/PDDDzZxVo37yvcuRFOlUhR5nrAQjc3atWtZuHAhP/30E2q1mhkzZhAbG8u8efPsUt+8efPw8vLipZdeoqSkhMTERD766CPi4+PtUp8QjYW0hIVohIYPH05kZCRt27alffv2mEwmHnvsMbvVN3v2bDZs2EBMTAxxcXFMmzZNErAQtSAtYSGEEMJJpCUshBCiyUpLS2P06NEMHDiQxMTEahMYn3zySQYPHkyXLl1YsmRJvdcvLWEhhBBNVl5eHmq1moCAAJ577jnKy8t58803ra+fOnWKmJgYjh8/TkJCAvn5+fVav9wnLIQQoskKCgqy/r9Wq7VZaAawLkN7tdfqQ4NNwjNnziQqKsrZYQghhHAR6enpfPjhhzb75j8fgr6ieoevV8gzPPvss9ZyXl4eS5cuZcOGDVe99vz583n00UfrN2AacBKOiorir3/9q7PDEEII4SKulhP0FQovPh1cbf/L7+mt/19RUcHEiRN59dVXiY6Ornbse++9R2ZmJosXL67XeKEBJ2EhhBDiWhQULFhqfN1kMnHHHXcwceJEpkyZUu31L774gi+//JI1a9b84WM960qSsBBCiCbrq6++YsOGDRQWFrJixQpCQkKYPHkyRqORyZMnc//999O5c2fGjBkDwMKFC+ndu3e91S9JWAghRKNlURTMSs0t4SlTply1BXxJeXm5PcKykvuEhRBCCCeRlrAQQohGq3JM2HWXw5CWsKgVi6KQqy/DZKm5W0cIIVzNpYlZVf9zFdISFtd0pjifx7eu5kThBcJ1vrwYP5SBkW2cHZYQQjR40hIW1/Tavm38eiGDEqOBk4W5LNy3zdkhCSFErSiAWVGqba5CkrC4pkKD3qZcYLDvbEEhhGgqJAmLa+rTLAoNKmu5Z2ik02IxKxZWpv3K56d3UVLlx4Gj6c1GssuLsbjQr2ohhK1LE7Oqbq5CxoTFNT3aLRFvNzd+z80iXOfLkz0GOCUOs2Lh8V++ZOP5YwB8c2Y/H/Sbhr97/S+qfi1rzv7Om4c3kFtRSif/SN5OmESIp4/D4xBC/DELYHahpFuVJGFxTSqVivs79XF2GCRnn7ImYICDBRl8cXo3D8cmOTQOi6LwztEtnCsrBGBf3hkWHdnEyz3HOTQOIUTDJ0lYNBhXe/S14oRfuEaLmWKjbVd41bIQwjXIfcJC1JP+YW0ZHNbeWu7kH86dbRzfQvfQaOkRdPkxmh5qLf1CYxwehxCi4ZOWsGgwNCo1ixIms/LMfspNBm5p2QM/dy+nxLKwz0TeObKZvIpS+oS05tZWPZwShxDijykKLnVLUlWShEWDolVrmNS6l7PDwFPjxtNdhjs7DCHENSjgQutjVSfd0UIIIYSTSEtYVGM0mVn4vy2czsilWYA3T08aTKCP428DciVmi4UFyds4lJNFsM6b5wYkEeYttyQ5ksWwByp+AW00aq/RDqtXsZShlH4M6MFzHGq3Dg6rW9w4C4rcoiQaltf/t4Wvthywlkv1Rt6aNcGJETnfGzu38+9f91jLeeXlLL3ldidG1LRYytdC0V9BKQDcsJiOovZ90u71KooRpeAhMPxSuaN8LZbAxZKIRb2R7mhRTcr5PNtyVl4NRzYdx3Iv2JRP5uXKSlmOVP7NxQQMYAT99w6pVjH8djkBA1jSQb/KIXWL+lG5dnT1zVVIEhbVhAX62pTDq5SbojAf267nCF9f1CpVDUeL+qe5RtlO1H6AW5WdHo6pWzQJ0h0tqnlqYhKl5QZOn88lLMCXeXcMdnZITvf8gMEU6Ms5knOBUG8dzw8Y7OyQmhbve8F0GCyZgA/o7nJItWq3dli87oLyLwAjuCWg8p7hkLpF/XD12dGShEU1/t6evP6QLMF4JZ2bG++OHu/sMJostUdfLEGfgyEZtJ1Qu3d2XN3+f8aiuwMsBajce6BSuTusbnHjLKgw47q9VpKEhRANglrbHLSTnFO3W/trHyREHUgSFvVCURQUBdRq1/3FKYRoehQFLC40EasqScLihn2zah//+3YvBqOZfgltmTNrKCqZtCSEENckSVjckPRzeXy4ZBslpRUAfLN6HzFtQhk7uruTIxNCiIu3KMmYsGisUtJyrQkYKrt+MrMKnRiREEJc5uoTs+Q+YXFDuneJIqp5oLXs6+NJXPeWToxICCEaDmkJixvi5+fFoH7tWPHtPkwmM2GhvnTp1NzZYQGQkVPI35es51xOIc1D/Xlu+nDCg/2cHZYQwoEUwKJIS1g0UiWlFfyw/hAGgwmLReFkSg5LPt/u7LAAeH3ZJnYeTONsVgE7D6bx+rLNzg5JCCFsSEtY3JDi4nKKistt95VU1HC0Y+Xkl9iULxSU1HCkEKKxUmRMWDRmzUL96BQbaS1769yJ79XGiRFd1qZ5sE25dUSQkyIRQjhL5cQsdbXNVUhLWNwQjUbNK3++lU+Xbae83EDf+BgG9nON1YX+dO8wPNy1ZOQUEhUawBNTkpwdkhBC2JAkLG6Yv78Xj88e5uwwqvH0cOO5+4Y7OwwhhBNVrpgl3dFCCCGEqEJawkIIh1AUhb3nz1FQUcHAFq3w0MjXj7C/Jj0xKy0tjdGjRzNw4EASExNJSUnhs88+Y/DgwfTs2ZPHHnvMeuzy5cuJi4tjwIABLFiwAACj0ciMGTNITExkwIABHDhwwJ7hCiHsRFEUnt/yE3d+8yUz165k2rfLKTEYnB2WaAIsqDAr6mqbq7BrJL6+vnzxxRds27aNIUOGsGjRIgYMGMDmzZvZvXs33333HQcOHKCwsJAXX3yR7du3s2XLFpYtW8bp06f55JNPCA0NZceOHbzxxhvMmTPHnuEKIezk0IVslh85iEmpfJzN7vPn+GD/bidHJYTz2TUJBwUFERAQAIBWq0Wn0xEdHW0tA+h0OpKTk0lISMDLywuNRsOgQYPYuXMnGzZsYPTo0QDEx8dz6NAhe4YrhLCTMqMRo2Kx2Wcwm50UjWhKFFRYUFfbXIVDBmXy8vJYunQpGzZssO5btmwZsbGxtGvXjt27d+Pj42N9LTAwkMLCQvLz823263Q6ysvL8fLyckTYogZlpgr+m7aFCsXIqPBetPEJd3ZIwkGOp2azfsdRPNy1TB3XB093t1qdFxceSf+olmxPPwNAG/8AJnXsbM9QhWgQ7J6EKyoqmDhxIq+++qq1Fbxnzx7+/ve/s27dOgB8fHwoLLz85B2TyURwcHC1/e7u7pKAncxoMfF/Bz7mQMFpALZk/cb87vfTxifMyZEJezuWksWzb3zL+dxiAPYfPcebz05Eq7l2q0KrVvPRmNv4+MBe9GYjt7bvTJuAwGueJ8SNcvVHGdq1TW4ymbjjjjuYOHEiU6ZMAeDw4cNMnz6dr7/+moiICAB69+7N9u3b0ev1mM1mNm3aRFJSEvHx8dZEvX//fjp27GjPcEUtHCg4bU3AABn6PH46v9eJEQlHWZd81JqAAfYcOsOJ1Kxan++p1TK7VwJz4wdIAhbiIru2hL/66is2bNhAYWEhK1asICQkhKKiIsrLy3nggQcAmDZtGjNnzuSRRx5h0KBBeHt78+CDDxIWFsasWbOYNm0agwcPxsPDg8WLF9szXFELflodbioNRuXyeJ6H2t2JEQlHcXfT2JQ93LT4ekvPlHBtFkXlUrOhq7JrEp4yZYq1BXwtc+fOZe7cuTb7/P39WbVqlT1CE3XU3i+Kcc378t25nZgUM3GBbbmzlWsuB1lhMDH/0/WcPJtDkL+OJ+4aTHSV9aSbugMZmby5NZmiigrimkfy3NAk1Kqrd91NHduH/UfT2X/0HG5aDbeN6E5UeICDIxbi+lROzHLd7mi5W15ctyc63MK4yHhKzRV08muJVq259klO8M+vtrI2+XBl4Sz8fcl6/v3cZOcG5UJMFgt//n49R3MuAPB7ZhYh3joeToy/6vE+Og8WPTeJgycy8PfxIrpFiCPDFaJRkiQs6iTGN/LaBzlZZk6RbflCIYqioKqhpdfU5JeVc6ag0Gbfqdy8PzzHTauhZ8cW9gxLiHqlXHyKkqty3ciEuEEtI2wn/7QIC5QEfIVAnRctA227k9sGS3e9EI4kLWHRaD1y+wD0FUaOn80h0FfHk3fVfexarzfy5jvrOJOeR7MQX+bMHkZwkM+1T3RhWrWaV0cN5a1tOyjSVxAXFckDfXs7OyxRS3qzgbeO/ZfUsvMEufvycMyttPSWe/arurRspauSJCwaLa1Ww//dWz+PWFz03np+3FC5YtuRY5kYjGbmvzSxXq7tTN0iI/h48m3ODkPUwQenvmVTzj4AUkph0YnlLOzx2DXOanourZjlqlw3MiFcyLmMAptyRmZBDUcK4RhZFbbj99n6PJSLa3OLhkOSsBC1EBHmZ1MOD/d3UiRCVIryamZTjvQKlTkPV6EoYFZU1TZXId3RQtTCYw8Pw2iycCY9l2YhfsyZVT/d3ELU1YzocRgsRlJKMwl09+XhmFucHZKoA0nCQtSCt7cHf/m/cc4OQwgrN7WWOe3vcHYYLs/i4rcoSRIWQgjRaCmosLjw7GjXjUwI0aSVlRn47XA6+QWlzg7luqRn5vPlN7vZuuO4s0MRDYC0hIUQLuf4qSxeeXMNaWdzCQzQMfu+wYwc4vrPHz50LIOXFq7ifHYRWo2a226O49EZQ5wdVpMmK2YJIcR1+s+KnaSdzQUgv6CML1bucnJEtfPdDwc4n125XKrJbOGnLYfRVxidHJVwZdISFkK4nAqDyaYsiUzUlcXFbkmqSlrCQgiXMyC+LZ4ebgCoVNCnR2vnBlRL40d1w9fHAwC1CoYN6mh9H0JcjbSEhRAuZ9zI7vj5efH74XNEhPlx281xzg6pVrYmn6C4pAIAiwIu/BjbJsPVl62UJCyEcElJie1JSmzv7DCuy4Ej6Tbl3w6n13CkcBTFxR/g4LqRCSFEA+Oj8/jDshBVSUtYuJSiCj2vJf9MXnk5vSIimN6jV4NcD/fsmVw+X7KdigojA5JiuWm4699eI27cjCkDyM0rJT0zn+bhAdw/ZYCzQ2ryFCpXzXJVkoSFS3nihzVsTksFYN2pE1gUmBnXsJ5xW15u4JW/fE3KqWwA9u5OwdvHg4TEtk6OTNhbx/YRfPDG3eQVlBHor8PNTePskISLk+5o4TIqTCYO5mRby2ZFYV9mhhMjqptTJ7OsCRigrLSCX/ekODEi4UharYZmIb6SgF1E5S1K6mqbq3CdSEST567REOzlZbMv0NPTSdHUXXh4AL5+tu8jKNjXSdEI0bRdWjGr6uYqpDu6kTp2Nps3V24jv6SM2BbN6NSyGWt3HUWlUnFb/66MT3S9MUqVSsW8foP4R/JWcsvK6RgayjP9Bzo7rOsWEurL/Q8NZsXnO9HrjcT1acPEyfHODsvhiisqeGffDsqMRkZGt2dAVCu71WUwmXhv3U7yy/QktmvJiO4Na1a1aLokCTdSr/53A7+nngfg+LkLfL/7KEazBYDUrHzaRobQqVWYM0O8qiFtohncug0GsxkPbcP9eI6dEMfN43tiNlvQaptet6TJYuHBH75hR8ZZANacOs6i4WPtloifWbqWjYdOAbB231EqTCbG9epkl7pEw1L5FCXXnZjlOm1yUW+MJjMZuYW2+y4mYIDCMj37T59zdFi1plKpGnQCvkSlUjXJBAyQWpjPrszL98jm6cv5KeWkXeoqqzCyL/Xy57nMYGTbkVS71CVEfZMk3Ai5aTW0bBZUbd8lAd6e9Ihu7uiwmhx9uQF9ucHZYThFgIcnPu6298h6u7nbpS4PNw0+nrZ16WSpSHGRRcaEhTO8OHU4b32zlfyScjq2CCM6Ioi1u44CcNuAri7ZFd1YKIrC4gVr2fLj76hQMXx8T2Y8MaLO19u9+zRf/HcnBoOJhIQY7p7Wvx6jtY8QnTcP9+jD+/t3U2IwkBAZxaO9EuxSl0at5sFhCbzzfTJ5JWV0bhHG7BGJdqlLNDyKosLiQrOhq5Ik3Ei1CgvkzYcm2OybNLC7k6JpWrb8eJDVX+3CYlYA+ObzHXTr3Zo+A65/slBxsZ433/yB81mVj8c7fvw8YWH+jBjepV5jtodZcQlM7tiVYkMFUb7+aNT2+yK8tU9nhndtS0GpnvAAX7Qa1/3SFeJK8kkVop5lZxZYEzCA0WAmMz2/TtdKScm2JmAAk8nCqZNZNxyjowR56WjlH2jXBHyJj6cHUcH+koCFDQUwo6q2uQr5tApRz/omdaBZhL+1HNkiiMSk2Dpdq3XrUJo187OWtVo1bdqE3nCMouEyKxYURbn2gaJBkO5oIepZy+hmPPv3SaxdsQdUMOGuvoRekZSvh5+fF08+MYIv/ruTigozCQkxjBrVrZ4jbtpKS/X842+rSDmZTWCQjpE3d2fDjwcpKCijfYcI5j57M+7uzv+qNFjMPLvna/blnsXXzZPZsYMYHeX6wxLOVnmLkuu2N53/yRKiEerUvSWduresl2slJLQlIUHWnbaXf7+7kR3bjgNwPrOA0yezqagwAXA2LZeAQB0PPzbcmSECsPjIFr4/dxiALH0xrx/cQFJ4e3Ra+8w6F47huj8PhBDCAfJyS2zKBoPJppx9xZi8M12oKLUpZ+uLyTeUOSmahsNylfFgGRMWoo7MipmTJac4WyoPS69JcYmezduOcvzEeWeH0iC06xBuU/b2tr3nuHW0a4zB9wpugYf6cudlx4Bwwjz9/uAMAZdvUaq6uQrpjhYNhsFi5M1j/+Rg0WG0aEhqNoD72tzt7LBcytlz+bzwt5WkpF7Aw0PLXZMSuG+q699X7EzT7huIxaxw/GgGAYHeTLi9F18u20lRYTmxHSOYdp9rrF9+S6selJuNJGefxlvrzmOdhqB1wKxzYV+ShEWD8eP59RwsqhwTM2Fmc/bP9A9JpJ2vjJde8uX/fiEl9QIAFRUmvlvzK5Mn9sHLU8YNa6JWq7jvgSSbfS+8MtFJ0fyxu6L7cFd0H2eH0aAoqFzq0YVVuW5kQlRhsNguAWnGTImptIajmyaz2fbWFYPRjMloqeFoIYSzSRIWDUa/4ATCPJpZy+18YujsL0/KudLwmzoRHORtLfeLj8HXt+E9k1mI+qJQOTmr6uYqpDtaNBgRXuE82f5Rtl1IRqvSMiZiBO5qWaj/SnHdW/G3F25jxy8nCQz0ZsLNPZ0dkhBOZVGu3h3tKs83kyQsGpTmukjubHm7s8NwaR07RNCxQ4SzwxBC1IIkYSGEEI1W5YpZ1bufpSUshBB1ZLEorF6zn/y8Evr1a0e7duHXPkkIFyRJWAjR4Cx4bQ3r1h0EYPWaAzz/3Hh69KifZUJF46KgwnyVOciuMptEZkcLIRqUkhI9ycknreXc3BLWrfvdiREJV6YolZOzqm6uQpKwEKJB0WjUaLW2X10arXyViYZJuqOFcDBFUfh48UZ2bj+Bu7uGSVMSGTzcPo+kMxhMvPfeRs6l5xEW5s+s2UPx8mrYq2d5ebkzbmwPvvxqFwaDidatQ7hzcoKzw6qTQ6mZfPzjbiqMJgZ3i+H2Qd2dHVKjo6DC4sLtTUnCQjjYujUHWPH5DuvqVovfWkeXHq0ICfWt97reeecn1qzeby2X6w08//yEeq/H0aZPH0TfvjFkZBQQHx+Nr6+Xs0O6biXlFTz/6Y+kZecDsO9kBkF+Om7q0c7JkQlHct2fB0I0UmdTL9gsL5mfV0rKyWy71JWakmNTTru4rnRj0LFjc4YO7dwgEzDA0bPZ1gQMUG4wsv9UhhMjapwsqDAr1TdXIUlY1MhoMjs7hEapY5coPDwvz80MjwygQ6dIu9QVHOzzh2XhPK3DAgny1VnLKiAqxN95AQmnkO5oUc1vKZn846uNZOWX0CY8iFfuGUV4UP13lTZV/QfHMjPnJpK3HkPrrmXy1ET8/O3Tmpv9yDDKygycO5dHWHgAjzw6zC71iOsX4u/DUxMHseSnPVQYTfTt2IpJMiZc7y7NjnZVkoRFNW+t3MbhM5Xdo7nFZbz1zTb+fv8YJ0fVuEyYFM+ESfF2ryc01I9/LLjT7vWIuhkT35Ex8R2dHUajVrlilut2+rpuZMJp8kvKbMvFZTUcKYQQ4kZIEhZA5W0zBfpy9CYjHaJCbV7r1CrMSVGJhk5RFIyW+plboCgKBeV6LIpy7YOFuKhyxazqm6uQ7miB0Wzm8R/WkHw2DU83d+7p2oO7fHVkF5TQNjKEB0f3dXaIogFak3aEt3/fSpGhgriQ5rzZfwIemrp95ZzMyeXZb38kNS+fSH8//jpmKHEt7DOZTQhHkiQs+PDXPfxw6gQARQYD7+/bxZop9xLl5+fkyERDpTcZWbB/I2dLCwH4If0YrX/fxrweQ+p0vdc3/szvmVkAHMu+wJubtrP0nkn1Fq9ovFx9YpZ0Rwvyy8ttysUGA+dLip0UjWgMCgx6cvSlNvsuVClfj6LyCptyYbm+ztcSTculiVlVN1fhOpEIpxnSOoZAT09ruWuzMLo0a+bEiBqH8+fyOHM6B4vFckPXMRpMfPPxFv77zk9cyCyop+jsK9TTm86Bl+cSuKs19AltUefr9WgRYVPu1lweXSgaB+mOFiS2aMFrw0fzw8njeGg0zO6TgKeGgUTqAAAgAElEQVTWVR701TC9++p3/LhyL2aThYTBsTy/8E402ut/jLjZbOGVBz5i96YjAGxauYeXlzxEWFRQfYdcrzRqNe8MuI3Xf9tCibGC/mGtmRRT93tgn7ppAH6eHpzIziUq0I9HByXWY7SiMatcO/ra3dGHDh1i5syZ9O/fn4ULF1r3l5WVMXbsWOuP6ePHj/Pyyy8zc+ZMduzYwVNPPYVGo6F169YsXbr0uuOTJCwAuKlNNDe1iXZ2GI3C/l9O8f2K3ZhMlf9okzccZvVXu5gw5foTx+E9KdYEDHDmRBbff76D++bdXG/x2kuYzpcFfcfWy7XUKhUP9bf/fdWiaTKZTLz11luMHz+erKwsm9d0Oh0bN24EwGKx0LFjR0aNGkVFRQV33nknq1evpmvXrnWuW7qjhahn+bnF1gR8SWlx3cYwte4aVGrbX/EajetOMhHC1VgUrrl2tFar5YMPPiAiIqKGq1T65ptv6NGjB1FRUSQnJ9O1a9cbSsAgSViIepeQFEvbjpdvnwlvEcigkXX7hxrboxVDbumF6uJ3RoeeLbl15uB6iFLYy55DaSz4eD1vL91EYXH5tU8QdlWfE7PeeOMNHn/8cQAyMjLIzMzk5ptvZsCAASxYsKBO15TuaCHqmc7bk5f/dQ8rPt2Gxaww6rbeRLUOqdO1VCoVT785lYE396C0sIz+Y7rjqfOo54hFfdl35Cwv/HMt+UWVq8wdPpXFO3+ehFsd5gMI17J79270ej39+vUDwMPDgw4dOrBs2TIMBgOdOnXirrvuokWL65uAKElYCDsICvHlwafrZ71tlUpF3+Fd6uVawr627j5pTcAAvx0/x8kzOXSMltnczqIoqnq5T/jNN9+0toIBunTpwuuvvw5U/htVqVR4eV3/g1gkCQshRD3x1rnblHVe7gT762o4WriS8ePHk5KSQnFxMampqUydOhWj0ciUKVPIysoiOTmZTz/91Hp8bGwsAwcOpE+fPnh6evKnP/2JkJDr7/GSJCyEEPXknvEJHDl1nt2HzqLzdOOuMb1pFiwrzzmTArW6Rem7776r8bWwsDBSU1Or7a/rOPCVJAkLIUQ98XDXsvCZ28jKLUbn6Ya/r32eEy0aD0nCQghRj9RqFRGh0vp1FZWzo133tj5JwkIIIRqtyolZrns3riRhIWqwa38qm3cex8vDjXsn9cXP53LX4uFdp/jkla8pyiuhfVwbHn9jGlq3y/+c9qw7wA+fbkatVjPhkZF0TmzvjLcgGpj3f9jB5oOn8XDTcu9NvRnSNcbZIQk7s2sSTktL4+GHH6akpASTycTnn3+Oh4cHd999N0ajEX9/f5YuXUpAQADLly9n/vz56HQ6xo8fz7x58zAajTz88MMcPnwYjUbDu+++S/fudV9/Voja2n0gjZfeXkNhUeViC0dOnWfRX+9Aq9WgKArvzFvG6YPpAKQdyySomT/T/3IrACf3p/LazPfJP1/5sIWju07y2vq/ENaybvcKi6Zh9a7DfLBuF0Zz5WprGbkb6d4mgiAfmV19IxSa8KMMfX19+eKLL9i2bRtDhgxh0aJFPP300zzzzDNs3bqVgQMH8uabb1JYWMiLL77I9u3b2bJlC8uWLeP06dN88sknhIaGsmPHDt544w3mzJljz3CFsPp590lrAgY4ePQcZzLyASgrLic7Pc/m+PNnLlj/f/ePB6wJGOB8ag67vv/VzhGLhu54xgVrAgbIKizhWHq2EyMSjmDXJBwUFERAQABQuTanTqdj69atDB8+HIChQ4eyY8cOkpOTSUhIwMvLC41Gw6BBg9i5cycbNmxg9OjRAMTHx3Po0CF7hitukKIozg6h3lS939Pb24MAv8ruaJ2vF1ExYTavt+7YHKj8M4ho0wyN9vI/LQ+dO606Rdk5YtHQdWzRDPcrPjcRgb50bBH2B2eI2rj0FKWqm6twyJhwXl4eS5cuZcOGDSxevBiNpnIJt8DAQAoLC8nPz8fHx8d6fE37dTod5eXldVqVRNjPzgu/sCpzNXpzBbF+scxocx9qletOhKiNeyf25dipLH49dBadlztTJvQhKMAbuLiU5LvT+fiVlZQUlNKhVxui+rVjxhNLKCnV0zk2kvGPjuLnFTtRa9SMvG8w3QZ2dPI7Eq5udK9YMvKK2HLoNO5aLffd1JsAb/muu1GK4trd0XZPwhUVFUycOJFXX32V6Oho1Go1JpMJrVaLyWQiODgYHx8fCgsLrefUtN/d3V0SsIspNZXx37PLyTdWdtX+fGE7EZ7hjI2snyUbncXDw43Xnp9ITl4xXp7u+Pl42rwe1TacF5bMAsBoNHPvYx+TnlnZBZ2ZXcTdkxL4bP4UVCpQqxv2DxLhODOGxzNjuDyysSmx67eDyWTijjvuYOLEiUyZMgWA3r17s379egC+//57hg4dSu/evdm+fTt6vR6z2cymTZtISkoiPj6edevWAbB//346dpTWhKspMORbE/AleYa8Go5uWNRqFWEhftUScFVFJeVcyCux2ZdzoQSNRi0JWAgnu3SfcNXNVdi1JfzVV1+xYcMGCgsLWbFiBSEhISxcuJAHHniA1157jcjISP7973/j5eXFI488wqBBg/D29ubBBx8kLCyMWbNmMW3aNAYPHoyHhweLFy+2Z7iiDkI9m9HSqwVnys8C4K5yp71v07odJ9Dfm+jWoRw+lgmAVqOmU4c/fi4pgL6sAjcPNzQaxybqQmMJaSWZtPFpjq+bzLwFKDGVkK3PIdIrAk/NH//oEqI+2TUJT5kyxdoCvtKOHTuq7Zs7dy5z58612efv78+qVavsFp+4ce5qN2bFPMg3GauoMFfQPaArfYPrrzvtox17WHvoGGqVmqm9u3NL9071du36olareOGpsXz4n58pLaugZ9eW3DK6Z43H68sqmD/rU47uT8PH34u7nxrD4Am9HBJrcs5v/OvUCnINhYR7BPNE+zvpHti0fjRVtSt3D5+f+YJ8YwERnuE8ED2DGJ9oZ4cl6kmTHxMWjV+kLpLZbR+q9+tuOn6aRVt2oDeaADi7voDOEWG0axZc73XdqMiwAF54amytjl26cC27Nh4GoCivlCUL1tBvVHfcPez/z3HFuY3kGirnWZyvyOWrs+ubfBJek7mWfOPF8Xz9eVZlrOGJ9o85OSpRX1x92UoZsBIu68j5bGsCBsgv07M/PcOJEdWPovwym3JhbgllxeU1HF2/KswGm7LBYqrhyKajwlL1z8RQw5FC1D9JwsJldW8egY+7m7Uc4qOjT6vmToyofnTr1w53j8vvq223FvgH+/zBGfWnT1AnNBf/2burtMQHd3ZIva6si//lIQ43lRvdA7o5MRrnURSFYmMxpkb2w+zSowyb9H3CQtRF/5hWPDNsEKsOHkWjUjO1T3daBwc5O6wbNnxSPGaTif0/n8Dbz4t7nhmNSuWYL4V7W99MhGcwaWXnaefbgiHNejukXlc2teVdhHmGkV2eTbRvNInBCc4OyeHyDQW8e/J9zpSl4+/mx10tJxEX2MPZYTUJkoSFS7uzVzfu7NX4Wiaj7urHqLv6ObxelUrFyIhEh9frylQqFcPDhjo7DKdafnYlx4pPAFBuLmf52ZX0DOjusB+H9qS42C1JVUl3tBBCNHFlZtt5CiWmUixYajha1CdJwkII0cR19OuAFo21HOPTBo1K8wdnNByXnqLUJBfrEEII4fpGhg9Dq9Jyovgk/m5+TGxxi7NDqjeu3h0tSVgIIQRDwwYzNGyws8NociQJCyGEaLQudUe7KknCQgiHOHg8g/+u2ovJbGb4gI4M7dfB2SEJ4XSShIUQdpebX8rLi9ZyLqtyycwDR84R5K+jZ+cWTo5MNHqKCsWFW8IyO1oIYXd7D56xJmCAohI9ew+ecWJEoqlw9RWzJAkLIewuplUIPjp3a1mtVhEZ5u/EiIRwDdIdLYSoF+dzi1j4n03k5JfQKiKQP903HK+La2THtAxl5uT+/O+HXzGZLAyMb8vopMaxbvXP24+zfv1htFo1k++I5/eD6ezZk4JO58H0+wbQvHnDX2q1IXP1pyhJEhZC1Iv/9+l6dv6eCsCR1Cw83LU8P32E9fVJY+K4/eJzlhvDcogAv/12ltdf/57CwsqnYP36ayolJXqMxsrVptLT83jnn3ej1TaOhS9E/ZPuaCFEvcjIKbQtZxdVO0alUjWaBAywd1+qNQED5OeXWRMwwOnT2Vy4UOyM0MRFilK5YEfVzVVIEhZC1IvmzWzHeCOb+TkpEscJa+bLlb8p3LS2X6khIb74++scHJW40qXuaFm2UgjRqD137zBeX7aZ7PxiWkcG8dSUIc4Oye5Gj+7OiRPZJO84gVarZtzYnqSk5PDb72fx9vbg7rv74+Xlfu0L1UGZwcjfvt1I2oV8IgP9eH7cEAK8vexSl7AfScJCiHoRFuzHgjnjnR2GQ6lUKh5/fASPPjoMlUqFWl3ZwlIUxe7d7n/7dgPf7jsCwP4zmRhNZt6aNs6udTZIF7ujXZUkYSGEuEEajW03tCPGvc/kFtiU06qURcMgY8KiUakwmjiReYGiMr2zQxHCriL8bcfcIwJ8nRSJa5MxYSEc5ExOPs98uoaj6TmE+ul4YvxAxvbp5OywhLCL58cPwWA2cya3gMhAX/48fqizQxJ1IElYNBrv/bCTo+k5AOQUlfHhT7u5uXfHRnVLjBCXBHh78baMAV9T5S1Kzo6iZpKERaNRXmGsUjZgURQ0koSFaLIurR3tqmRMWDicyWgm53whRqPphq9lNJo5dfw8eRdKGNQ5Gi83N+trvdu2QKOWj7gQwnVJS1g41NGD6bz90rekp+USERXInD+Po0tc6zpdqzC/jBef+ZLDv53Fx9eTu6YP4OWpw9l76hyhfj7cN7R3/QYvhGiAXGuFrKqkmSAc6j+LN5JyIgujwcSZ0zksfW9T3a/18RYO/3YWgJJiPSv+s4P+7Vrzp9tvYuaIeLQa+XgLIVybfEsJhyotqbAtF9X9VqIKvW13dlmZAb3eUOfrCSEaH0XBpW9RkiQsHKpLXCubtXa7xLWq87UGDInFz//yMn09erUmKNjnRsITjURRuZ4V+w6y+dhpFBeaGqsoCutW72fZx1s5dfy8s8NpEhQuz5C+cnMVMiYsHOr+OcMJDPIh9WQ2Ua2Duf3e/nW+Vny/djz3t4ns/Pk4Pr6e3HnvALkdSZBVVMJD//mGo1k5aFQqbu3Zmb9NGO7ssABY9I+1rP1mL4oFVv9vL8++fCvde7V2dljCiSQJC4dSqVTcdne/erteXEI0cQnR9XY90fAt2bGPo1mV94ubFYVvDxxmZv/etA4JdGpc5eUGtm86gnLxSYe5OcX8uGq/JGF7c7FHF1Yl3dFCiEbFgm1fo8WiYDRbajjacTRqNeoqkwXVcgtdkyefgCbColhYfW4nS1PWk1aa5dC6czILWfbvzaxY8jOGKgtqiOoURWFb9q+sOLuB9DLH/l3VJ6PBxIp3f+LT+d9y8uDZaq+bjCa+XvQ9n73yNWlH0uut3sm9utE6+HKrd1jHtrRtFlRv168rdw8toyf0xM29sgMyqmUQt92V4OSoGr/KMWFVtc1VSHd0E6AoCvMP/5efzu8DYG3GLl7qeg/t/aLsXnd2ZgF/fmQpZ05Xdg/u3XGKV/45Da2bxu51N1TvnfqaNRk/Y8bCqoxtPBt7Lx392zg7rOtisVh4deYH/LLudwA2LN/Fnz98gA4X7wm3WCz8beo/2bGq8jO5ftk2Xvjycdp2b33DdbcJCeTDu29l7cFj+Hl6MqlXF5eZK3DvQ0PolRBN+pk84vu3k4mEDqC42GzoqqQl3ATkGorYmv27tZypz2P1uV8cUvcPK/daEzDArztPsX/XaYfU3RDpzRVsydmHmcru0+yKfNZkJjs5qut3/kwuezYetpYvZOSz8X+7rOVzJ7PY9f1+azkr7QI//efneqs/KtCfBwfGc2efbi63alqXHq0YNb6nJGABSBJuEjQqNRqVuso+x/wy1GqrPGdVrcLD062GowWoqq1y6yKNuOvi6eWOu6dtR9uVvR8eXm5o3au8rpWvI2EfrnyLknzqm4BAd1/GRMajVVV+CbbxDmdSqySH1H3L1ESCWwZh9tRi9tTSZ1CHG7o3uLHz1LgzNCweLZV/VxGewYyPHOSw+k9n5vLyf37ir0vX8evJc9c8fv3KvSyc9yXvvfod5aWXF2IJCvPn5nsH4eZRmWjb92jFpEdHWF9v1iKE0dMHo3WrfL1dXBtuf3JsPb+bhiflVDb/XLCWdxd+z7kzec4ORziAjAk3EY+2n0BiSCey9fn0D+mMn7u3Q+r9ce3vXCgoA21lUsktKcdiUdBoGmDzzkFmRk+gm39bMvW59A3uQpinYyYV5RaV8vT7q0jNygfglyNpLJp9C+2iQq96/E9f7+Gdv36DQV852S79dA5/+2iG9fUZf7mVIbf1IT+rkC6J7fDwcrc5f9bCu0m6vS+5mQX0GdEVT29PO72zhiEzI5+X533FubOVyXf/7lQWLL6bwCDptr4RlS3fq33fuEZzWJJwE9IrqJ3D6zx14jwWy+UPe1rKBfLzSwkJ8XV4LA1JfHBnh9f58+8p1gQMkJVfwpbfTteYhH9NPmlNwADHfjtLSVE5Pn6XVzGL7hwFnWueANipr+M/k65q24Yj1gQMkJaSw7YNRxg/qY8To2oMapoN7RpJWLqjhV01C/evVvb31zkpGvFHokL9cddeHrdVAaEBNbfCvH1tW66+fl546dxrOFpcS0Cgzmb8X61WERwqP1YbO2kJC7uaes8AcrKKOHQwHR8fT+6bOQg3B92edOZ0Np+8vY7igjJiu7Xk/idHyOIIf6BX+xZMvakn3yQfwmyxMKRHW8Yndqrx+HufHEn6qWyOH0zHN0DHtDnD0Whr/3d74tfTLF+4CrPRxE1TBtL/lnib11d+tYs9O0/h5eXG1OmDaNO2WZ3fW0Mw/ObuHNibxtb1h1CpVQwb3Y1+SR2cHVaDp+Aqbd6rkyQs7EqjVfPUnxw/4UZRFN74y9ccvfiow4P70tD5eDDloSEOj6UheezWgcwYnYDZYsFX98djtD5+Xvy/JQ9QmF+KzscTd/faf50UXCji/931FunHMwHYv/kQfsE+dB1YmfTX//AbH727HoPBDED6mVwWfTgDd4/G+5WlUql4+oXx3PfwENQaFcEyZNMkSLNANEplJRWkp+bY7Es72XBXn3Iknaf7NRPwJSqVioAgn+tKwAD7Nx20JmCAogvF7P3pN2v5yO/nrAkY4PTJbDIz8mnsVCoVoWF+koDrk+LaK2ZJEnYxRksp5aZsFMX5a902ZDofD5pFBNjsC49y7gL+4rJWHaPQXTGmrFKraNYyxFoOqzqXIMyPkAY6PpqfX8qevSkUFpY5vO6c4lJO5+Rhtsj3iatqvH07DdDJwuUcyf8Qg7mYEK9u9At7HTeNY24lamxUKhVzXpjAx2+to6igjPado7h79jBnhyUuatOlJdNevINV//oBo8FM4rjejJ4x1Pr6xCl9OZeex4F9qXh6uTPlvgF4+zS8W5iSd5zg7bfXkZNTTLNQX56cO4qE+BiH1L1oYzKf7fwVvdFIQpuWvHvXeDzdmuhXvgsPCjfRvxHXY7SUciT/Q/TmCwBkl+/mcP6/6R7ypJMja7hiu7VkwccznR2GqMGkueOY+MTNKBal2oQujUbNk06YS1DfvvpqFzk5xQBk5xTz5Ze/OCQJp+cX8umOfZQZKm8h234qjQ9/3s2jQxLtXrercbXu56okCbsIo6UEg7moyr5SJ0UjhGOo1epGPShmMJhsysYrxrntKa+0zJqALyk1GBxSt7g+jfjj37B4aUIJ8exhLWtVPkTo+jsxIiHEjUrs29a6Jrabm4a+iY7pio4Nb0Zci0hrOdRHx4iOzlsYxaIYKTWmY7I4flwcXHvtaGkJuwiVSk1i+Gsczv83JkspEboBNPeR22mEaMjuvrs/YeH+nD6dTduYMIYNc8xKaO5aDe9NncD723ZTYTIxqnN7eraMvPaJdlBiOMOu7D9RUHEML20zugY/SZTPcKfE4ookCbsQd40PPULmOjsMIUQ9GjG8i1Pq9fPy5JkRA51S95UO5y+moKLysZblpkyO5n/g0CRc89rRrkGSsBBCCLup2gVtspSgKAoqhz2jUwUunIRlTFiIOrBYFJZt2sdr/9vMhv0nHFq32WLhm/2HWfrLrxSV6x1a941QFIUDPx9j09e7KS9pOHGLGxPmlYiay2uKB3vGOTABuz5pCQtRB//vyw2s2P47ACu3H6SwTM9t/bravV6LReHJ5WtYd/hkZd2/HuaDabcQ7OP695P/6/++4PulP2M2WWjXvSUvffEogaF+zg5L2FlMwJ1o1Tpy9b/hpQ2lQ+CMa59Un1xsIlZV0hIW4jopisKOo2nWcpnByObfTjuk7n1nM6wJGOBwZjbLdh1wSN03IiMlmx8/T8Zsqly56cSBM3y9eL2ToxKO0spvPHHN/kzHoIdQq6TtdyX50xDiOqlUKjyqrDzk4aAnQ13tJ73iyj/zLzJWmDCbbO+RNRsdc8+sEK68Ypa0hIVLO30ii7fmr+bt+atJOZXt7HCs7hvem2C/yucix0QE88CoBIfUG9eyOcNjL99r2jE8lKkJPf7gDNfQon04vYdeniUcGhXEyKkDnBiR+DkljTkrV/PEN2vYc/acs8OxG8XFH+AgLWHhss5nFPDy/31FRnrl03P270nl7+/eXW1xf2cYn9CZAZ3acDankHbNg9F5OOZh9mq1ircmj2Xlr4cpMxgY160jgd5eDqn7RpjNFkrd3FD8vMFiwRLoC47qPRDVHMnK5v/W/Eh2SeWqfL9mZLL0rttpGRhwjTNFfZOWsHBZP286Yk3AAOfO5pG8+agTI7IV5Kuje3SEwxLwJRq1mtt7deGexLgGkYABDu1L4+C+M+DtBb7eXLhQyqY1rj+W3VhtO51mTcAAGUXFbD2d4sSI7Ey5yuYiJAkLlxXSzBe1+nK3kVqtarCPs2vqfHw9rcs3XuLuIR1xztIyMAA39eW/D0+tlpjgkD84Q9iLJGHhspKGdWbMLXF4+3jg4+vJzbf1YsBNHet8PUOFkaXvb+L9N37g0P4z9RipuJa2nZozamIftBe7oLv3iea2e2RtdGcZFduO6X3iiPD1IdLPlwf79iGxdQtnh2UfVxkPljFhIWpBpVIx59mbue/hyjW0/QJ0db6WxWLhpaf/y57kytt7Nv3wO8/Nn0S3Xq3rI1RRC48+P45RE3tRVlJBp+4trQlZOMczQwYyN6k/KpUKtSye4TSShIXLu5Hke0lWRgH7d10e88rPLWHrTwclCTtY21jnPERAXJ1G3UQ6Q11oDLiqJvI3IJo6L507nl5uNvs8PN1qOPrqFEWhrNzQIO7LFUJcSXWVzTX8YUt4/fr1rFy5kuTkZE6dOkVZWRk+Pj60bduW+Ph4xo8fz8iRI2Ud0EbGYlF4Z8VWDp4+T4CvF49NGkiLZoHODuuGBAT5cOuUviz/LBl9uYFO3VsweXrtnzCTdiaXf7z9PekZ+YSH+TN39nBi20fYMWIhRFNQYxLu1asXnTt3ZtiwYcyaNYvWrVuj0+koLS0lNTWVffv2sXTpUubNm8dvv/3myJiFnX20aidLf9hrLecVlvHhc3c6MaL6Me3BIQwb24OC3FJiYsNxc6v9aMz7n27h0NEMAAqLyvn3p1t54/9NtleoQoj64mK3JFVV47dQcnIyHh4e1fb7+vrStWtXunbtyr333oteL09DaWxOnbtgU07LykdvMOLpfn3dt64oPDKQ8Mjrb9UXFZfblAuKymo4Ugghaq/GJFw1AZ88eZL09HQslsoF2DUaDUlJSXh6eto3QuFwEcG2T7aJCPKttlZyU9OhfQS/H768tF+HduHVjlEUhZKSCnx8PGSI5qKM8wVkZhXSqUMEXp6OXdTEVZ3OzSO/rJyukeG4a2SGuEM0xJbwlWbPns3KlStp37496ouz6bRaLUlJSdc899ChQ8ycOZP+/fuzcOFCTp48yUMPPYSiKOj1epYuXUpMTAzLly9n/vz56HQ6xo8fz7x58zAajTz88MMcPnwYjUbDu+++S/fu3W/sHYtrmnVbf/KLyzmSlkWAjxdzJg1s8kll9v2D8dG5k3Y2l4gwf2bcbTuefPR4Jgvf+ZHMrEKaRwQyb85I2kaHOSdYF/Hlt7v59L87KC2roG2bUF55dgLNwxv23IIbtXDzzyzZ8yt6k4neUc359+3j8ZWGjH0pqsrNRdUqCa9Zs4aUlJTrbvWaTCbeeustxo8fT1ZWFgCLFi1i6tSp3H///bz22mssWbKEp556ihdffJG9e/fi7u5OXFwct99+O+vXryc0NJQdO3awa9cu5syZw5YtW67/XYrr4u6m5a8zRzk7DJei0aiZ/gcPHPjgs62cuPiAiWMnzvPvJVtZ8NIkR4XnckwmM8u/20tpWQUAJ1NyWLr8F559rOl+rjIKi1i6dz96kwmAPenneH/nbp4eXPsJgqLxqdUtSiNGjCA5Ofm6L67Vavnggw+IiLg8i7RPnz6sWLGCU6dOsW7dOkaNGkVycjIJCQl4eXmh0WgYNGgQO3fuZMOGDYwePRqA+Ph4Dh06dN0xCOEIRcX6Pyw3NUaTmYoKk+0+o6mGo5uGYoMBvdFos6/CJI9ztDeFS09Sst1cRa2S8JgxYxg9ejQtW7akZcuWtGjRgujo6DpVOGXKFDQaDf379yc8PJz4+Hjy8/Px8fGxHhMYGEhhYWG1/TqdjvLy8qtdVtQDg9FUp3tgTUYzPyzbzqqPt6AvrajT+YV5JXW+/9ZgMGE2W+p0bn3p2MF2EYrOHZr2ohRenu7Ex7W2lv18PRnSv4Nd6zSbzJQWue73Q9vgIPq1bmkth/v6MLaTff9MhOurVXf0448/zueff05cXJ/5HwgAACAASURBVJx1TFhTxwkFjz32GGPHjuW7777jL3/5C3PmzGHUqFEUFhZajzGZTAQHB+Pj42Oz393dHS+vhvHUmIYkO6+YFz/8gdMZuYQG+vDMlCF0b9e8VueajGZeuvc99mys7KXY9PUuXv3yMby8azd0sW3Nfj75x2oKc0uI6dyc5xbfR0Bw7R7SYDZbeOsv/2Pvzyfw8HTjtnsHMG5qYq3OrW+PPzSUQH8v0s/l0yIqiLsnOycOV/Lc42NoHxNGXn4Z8XGt6dWtld3q2vz1bj77xyqKC0pp37M1f/7ogVp/Bh1Fo1azeOIEPvxlN2VGE2Ni29MlomnPG3AIF79FqVYt4bZt2zJixAjatGlDq1ataNWqFVFRUXWqMDk5mUGDBqFSqRg9ejQnT56kd+/ebN++Hb1ej9lsZtOmTSQlJREfH8+6desA2L9/Px071n3xflGzRcu3sufoWfKKyjiWls2ir7bW+tzta361JmCAI3tSWP1p7c5XFIXPFq4lM/UCZcV6ft95iiWvra113d/9J5n13/xK/oUSzqfns/Sf68k5X1Dr8+vTpTHjv8wbx31T+qPRyGJ0Go2ayRP6MOu+JLsmYKPBxJK/f0tmag4lBWXs23SYJfO/s1t9N8LTTcujAxKZN2SgJGBHujQ568rNRdSqJdy5c2eSkpIYO3asdZ9Go+HFF1+85rnjx48nJSWF4uJiUlNTeeWVV7j77rvx9fXFaDTy5ptvEhkZySOPPMKgQYPw9vbmwQcfJCwsjFmzZjFt2jQGDx6Mh4cHixcvrvs7FTUqKLEdv8wrrv09sFfvQa7dB9xYYaKk0LausuvoTszNKbYpFxeVk3kmj9BweTB5U1Jeoqcor9RmX9XPlRCuqlZJODw8nAkTJtjsq2139HffVf9FOm7cuGr75s6dy9y5c232+fv7s2rVqlrVI+quc5swfjmUZi13aFX7X+gDxvZk/Zcd2bv5SOW5ca0Zc0/Ns4iv5O7pRvvurdi1obIl7eaupWtiu1rXHZcYww/Ld1N6cRJU6/bhtOtcu2500Xj4BnrTvkcr9m87BoC7hxvd7Tz+LBoOlVK5uapaJeHnn3/e3nEIJ3rolv64abUcP5tNaIAPj91e+1smVCrQuF/+QRbaPBAv7+orrdXkT/+6l89eW0tRXgldEmIYdVftx1Lj+rdnzku3sH3dIdzctUx+cPB11S0aB5VKxXMfPsBnf19FaVEZPQZ0YPidMiYvGoYak/CoUaO4/fbbGTlyJC3+P3v3HR1FuT5w/Lt9s+mk0kMLBAiht0AoUgVCb4INFL3ixXb1iuXay0+9XtSLDa9iV1QEpEjvhN4JLQQSIL33ZMv8/oguLpAQwiYz2byfc3KOMzuz8xg2+8y85XmbXrvY86lTp/j1119Zu3YtW7ZsqckYhRqmVqu4L7p3tc5d+/0u9q47bt/e8dshto86yIBx3at0vtFNz5x/javWtQGiRnQiakSnap8vuAZPH3fmvln365sLNaQuPgl/+umnfPbZZ4wePZrk5GSaNGmCm5sbhYWFJCYmEhQUxJgxY/jss89qM15BYYqvMx9WydNEBEEQlKTCJNysWTNefvllXn75ZQoLC4mPjycvLw+TyUTLli3x9vauzTgFhRo0uRcbftrDhZPlKwy17tSMAeO6yRxVObPZyv++2k5Scg4hzfy4e4YYtewqbDYby7+JISM1l869WtEjyrEPWJIkdv52kKy0XPqP7Y5vgFcF7yS4PIWNhr5alfqE3d3dCQ8Pr+lYhDqoQYAXL387l1VfbkOlUjHu/kG4e5nkDguAdz9Yy5r15U3l23ZCUXEZDz9wm8xRCc7w3r9+Zd3S8uU2f/9pH3Ofj2bwmC721xc88hXrvtuJZJNY+flWXvruYRqGBMgVriA3BTdHi8cC4ZYFNPLlnvljufvpaLyrWGijNpyOS3XcPpsiUySCM1ktVg7sPGvfLiooJWZjrH079WIGm3/eg2Qr/+ZNPJXEsk821nqcglAV9Xt9OsGl+Xg7Vlfz8hLV1lyBWqPGYHD86tLpr2yX1wa+6tFHScWChdrlChWz+vfvz++//862beWVkNauXUu3bt0YPXo08fHxNRqgIFTXg7MH0S40GB9vE+EdGvPgrIFyh1QrjqSl8NKOTby+awvZJa43SE6lUjFxVn88/7jJat46kKlzBtpfD2rqR/9x3flz9c0mrYMZPXvgtW8kCApQpSfhnTt38vXXX1NYWMjq1atZtWoV7777LpcvX+b+++9n40bR1CMoT2jrID5+7y5KyywY9PWj0Sc2I42/rV1OUkF5NbG9yZf4YexUjFqdzJE51+1TetFrYDuSErJo3b6Rw/xwlUrFPxbeS9cBYeRk5DNgfHcCGvvJGK0gOwU/CVfpm6lhw4Z8++23SJJEx44dsdlsDB06FIA33nijRgMUhFtVXxIwwKpzp+0JGOBwWgoxly8yqHn1Vj1TMr9Ab/wCrz9LQ61WM2Ra31qOSFAkibo/Orpx48asXbuWsrIygoKCOH36NKmpqRgMBiyW+r1GqCAoiade77Ct12gINHlUcHT9JkkSKpVyv5yF+qFKSXjx4sW88soraDQagoODadGiBf369SM7O5tnnnmmpmMUBKGKZnXqzr7ky2xNPI9eo2Vmhwg6BATKHZbTpVzM4sN//UzqpSwahQTwyJtT8fGr2s2GzWbjvUe+4tC2kxhNBibPG8HQ6dV7apYkic8+38aBQxdwM+q4c0ZfunYJqdZ7CTWnzteObt++Pd9//73DPpvNRklJCSaTMuaECoJQ/uS7aOR44nOy8NDpCfZQzpQxZ/rohV/Yt7l80ZDEs6kYjDqe/uCuKp277OMNrP12p337fy/9Qo+h4fj43/zvavlvB/l+yW774Ov0zHUs+vAe3Nz0lZ8oCH+oUhI+ceIEb7/9NpcvX8ZqtQLlqyitX7++RoMTBOHmqVUqWvu69kCk9OTsSrcrPfdSlsN2TloeqYkZ1UrCFxIyHWY/Xb6cTWpaHiHN/W/6vYQa4gpTlCZPnkxoaChPPvkkzz77rP1HEARBDo1bOFa/atKi6stvdowMxfiX0dQt2jemebtG1YqjZYsA1Oor/cpNm/oSFChKZNZFJ06coE+fPvzjH/+45rUXX3yRAQMG0KFDB95///0qnVNVVXoSdnNzE32/giAoxrzXp6A36Em7nEXjFgE8+OKEKp8bOaorBW8UEbPmMEaTnumPj8Joqt4SmGNGdSYjI58DBxMwuum4a0akaIqugywWCwsWLCA6OprU1NRrXp85cyYvvvgiaWlptGnThnnz5t3wnKqqUhIeP348s2fPZsaMGfZ9Go2GAQMGVPvCgiAI1eXp486T/5lx4wMrMHxmP4bP7HfLcahUKmbdE8Wse275rYQaouLGA7O0Wi2LFi1i8eLF102orVu3BuDo0aO0bdu2SudUVZWS8LZt2zCbzbzyyisOQYskLAiCILi6PXv28Mgjj3Dy5EkWLlzo1PeuUhJet26dUy8qCELdYLHaWHn0JBabxJhO7TDo6k/hk7qgxGIhITebYA9PvA1GucNRJicsZdirVy92797N5cuX6dmzJ126dKFDhw5OCa/Sv6jPPvuM++67j1deeQWz2ezwmkaj4YUXXnBKEIIgKI/FamPud8vZdvYCAMsPx7LorgkYRSJWhLNZGczbsIqTmekEu3vwbN+BjGndTu6wXFJaWhqBgYEEBwfj7e1NSUmJ09670r8mm80GlCfcP//7TxqNxmlBCIKgPL+fOGNPwAD7Ey7z474j3N23m3xBCXbv7Y/hZGY6ACmFBSw8sFsk4YpUYYpSdHQ058+fJz8/nwsXLjBjxgzMZjPTpk3jnnvuIT8/n9LSUqZMmUK3bt2ue86yZctuOrRKk/CcOXMAxMhoQaiHrFfdeJfvU/CEy3qmyOLYOlloNotSnNdTxXnCK1asqPC11atX3/Q5VSWKdQiCAhw/msjemHM08Pcgenx3h7mnchnZMZRfD51gz/lLAIQ3DmZyt44yR+Ucpw4ncvJwImGdm9GuczO5w6mWQc1asPNSAqV/fCf3btxUJOA6qEpJePLkycycOZM77rjD3gwtmqMFwTn27DrL26+vJDenCIAzp5J56tlomaMCvVbLJzPHs+TAMaw2G5O6dsTDWL35tEqy9pf9fPbWKgrySvD0duO+p25n2ITucod10+7s2AUvg4EDyUk08vTk/ogecoekWHW+drTBYBBN0oJQQzatP2FPwAB7Y85RUmLGaJR/DWCDTsudvbvIHYZTrf1lHwV55QNr8nOLWfvz/jqZhAHGtmnP2Dbt5Q5DuAUVJuHt27fbR0QPGDBAFOsQhBqi0zm2Kun0GrTaKlWUFapDunpTwY9Jwq1TeO3oCpPwq6++SllZmcM+UaxDcCVmq5XFWw+QU1hMVFgLerV2Xt9gWamF1csPYLHYGBHdBQ+PiudwTp3Zh9Mnk7hwPgM3Nz1jJ3RHq624uyclKZufvtiB1WLjttERhHcLcVrcriLxbApLFm7AYrYwILorfYZ3sr82dHw3EuPTKcwrxt3LjWETxGhvl1cXk/DatWtrMw5BqFWSJPHkt6vZcDwOgBUHY3ltygiiwlrc8nuXlVl47onvOXzgAgBbN8TyxvszKkzETZv5s+Djezh6KIFGTRrQPKTiFXgKC0p46dHvOH+mvEze3h1neGnBHbTp0PiW43YVhXnFvPbAFySeTQHg4PYzuHub6NS7vPTgyCk9adYqkFNHE2kX0YwOXUNkjFao76rU5hUVFXXNvv79+zs9GEGoLTlFJew+m2jfzi4sYd2xs055720bY+0JGOD0ySRWLz9U6Tnu7gb69AutNAEDHNp9zp6AAbLS89m15dQtxetqju05Z0/AAPnZhRzYctLhmA7dQph4b5RIwPWASrr+j1JUOjArNjaWmJgYzpw5w+eff27fn5aWRkJCQo0HJwg1xaDTYtRrKSgt+8s+54z411ynP9dZM46CGvuiN2gpK7XY9/k0cHfOm7uIxi0CMHkYKCoote/zb+gjY0SCULFKk/CFCxfYsWMHxcXFbN++3b7fZDLx448/1nhwSlBQVsprhzaRUpRHW+9A/tF5AFq1GDRT15n0Ou6O6sanG/eQX1JGxyZBPDC4l1PeO2pQezb2PcbeXeVN3eGdmzHaSaNv24Q1YtLdkaz6aT9ms4We/UIZPaVntd+vuLiMr7/ZSVFhGZGRbejRo6VT4pRT09ZBzHhsJMs/34q5zEKPwe0ZdWek3GFVy/JNR1m74yQqtYrxt3ViSB9REatabrF2dE2qNAnffvvt3H777XTq1InHHnustmJSlCd3r+L3i6cB2JIUj0Wy8ly3ITJHJTjDvQO6M6JTKGl5hbRrFOC0xQk0WjUv/d9UNq8/jsViZfCwcPQG59Vbvuuh25h4VyQWsxUvH1O1CzTYbBLP/2spBw9eAGDL1pPMnx9Nr551PxFPmDOIsbOisFolp/7ua9OBE4l88O1WCovLW2vOX8qkZRN/WjatvMtCuIrCR0dX6ZHuv//9b03HoVhnctIdtk9mp8kUiVATGvp6EdG8odNXB9Jo1QwZ2YkRY7rUSBJw9zDi7et+SxWS0tLzOHr0on07L6+EXTvPOCM8RdBoNXU2AQMcP5tsT8AA2XlFHDp1ScaIhJpQpSQcERHBqlWrajoWRQp083DY9jeK/jfBNbibDHh4OFbAMrrJXyBEKNeuZSBuhiv/Ht4eRsLbNJIxospZbTYs16k3rgR1dmDWnxITExkzZgxNmjQByqd36HQ64uPjazQ4JXiu2228uG89yUV5tPb259lut8kdkuJcKNjH4ayl2LAQ6jmIjr63yx2SUAWenkamT+vNt9/FUFhYQnh4U+6qo32nrqhXpxbMmRLJup0nUatVjLstgtCQQLnDuq4P9sfww8lj2Gw2Rrdux7ORA+UOqc6oUhJetmyZfeGGP9WX2tEdGgTz0/A7xeokFcg1p7Ix5V0KLH8sqVZ8Gk9dIM096mYZwPpm8uSeDBsWTkFBMcHBPmg0YtChkkwb2Y1pI5VdTGRP0kU+PLiHYkv5iP3FRw8SERTMaKUsq+gKfcJNmjQhLS2NpUuXsnz5cgoLC+1PxfWFSMDXl1R01J6AAcxSEUnFJ2SMSLhZ3t5uNG7cQCRgoVrisrLsCRjALNlIzM2VMaK6pUp/dW+99RZ33HEHly5dIiEhgXHjxvHVV1/VdGxCHRBoCMWo8bJvq9HhbwiRLyBBUCCzxXrjg+qogc1CaOLpbd8OcDMxoNmtV55zFhUu0Ce8aNEi9u/fj7d3+S/6+eefp1+/ftx11101GpygfH7G5vQLuJ8j2cuxSVZae/WnjZeoKS7c2NKdx/htdyyoYHyfjkT36SB3SE6XkpnHv75Yy4XkLAJ9PfjnHYMJb9lQ7rCcqrGXN+8NuZ3Fxw5hkySmhHWkQ4CC+q4V3hxdpSSs0WjsCRjAx8eHoqKiSs4Q6pMOPiPp4DNS7jCEOuRg3GUW/LqN/D/nwCZn0aqRHx2aB8scmXO99/N2Dpwun1aUmVfEgp+28b9/TpU5Kufr1rAx3RqK+uXVUaUkfPvttzNu3DjuueceVCoV3377Lf369avp2ARBcFFHzyfbEzBAblEJR+KTXS4J5xQUO2xn5xdXcKRQoxT8JFzlPuGhQ4fy9ddf88UXX9CtWzc+/vjjmo5NEAQXFR4SjKeb3r7tbTLSqYXzErAkSaSl5JKZke+096yOsJAgh+22zQJkikRQqio9CWu1WubOncvcuXNrOh5BEOqBbm2a8MjY/qzYE4tapWJc3450DHFOX6nVauP/XlrGzq2n0WrVDB/dmYceG+6U975ZD4/vh0Gn5eylDIJ8Pfj7RLH6XK1T2ECsq1WahG+77TbMZvN1X/Pw8GD16tU1EpQgCK5vYv9OTOzfyenvu3r5QbasL58mZy6DlUv306d/KF261/6IXbVaxQPRfWr9ukLdUWkSfuWVV65Jwjt37uTdd9+97hrDgiAIcsvLdRw0arHYSE/LkykaQahcpUm4b9++9v8+cOAAzz77LHl5efz0008MGjSoxoNTAkmSiD2TRGmZlU7tm6AVBQ2cLjk1h8NHLxLWtiEhzeRbIeb8ySQ2Ld2HzqBl4gODcfd0ky0WofqiBrdn7cojpCTlANCydRB9+4fKHJUgq7raHA1w/PhxnnvuOS5cuMDLL79MdHR0bcSlCJIk8dqC1WzYGotNgl5dW/D6M+PROWnxdwFi9p3jnQ/WkpFZgJeXkb/NGsTtQ8NrPY4Lp5N5afYiUi9mAXB0Vxyvfz+3Tq/CU181be7Pv16fxO+/HUatUTHpjj54iBsqQaEq/YZ54IEH2LlzJ08//TTTpk1DrVZj+2OVDJVK5fKlHGMOxLN+ayzSH3dRew6e59c1h5kSrexarnXJLysOkJFZAJQvpffryoOyJOEtyw7YEzDAiX3xnNgbR5f+Cql/K9yU1m0b8nBb1yqKIVSP0ipkXa3SttXvv/+e2NhY7rrrLgwGAzqdDp1Oh1arxd3d9Zf0KygotSfgP5WWXn+gmlA9Fovj0mdWqzxLoRn+Ml0GQKvT4OnrUcHRgiDUKdJ1fhSi0ifhvLz6PZghqk8b2q9qSOyZZABCmvoxfFD7Wrt+dukpLhVsQKf2INRnBmqV6631OmRAGKfjUigqKkOv1zCwX9tKj09Pz2Pp0v0AjB7dhcaNfZ0Sx8Q5gzi2O45D20+j02u5fWYkrTvWzUVKzl3KYOnmo6hUMG1oF5oEOed35Gy52QXsWH0UHz8P+g4Pd/mWtdqWll3AV6v2YrZYGRnZns5tql7RSpIklizfz+WUHNq0DGTMsAinxnbgaAIx++Lx8jQyfXzPet3FJzq8KmE06HjnxUksWXEAq9XGmGGdCPT3uvGJTpBZfJyY1KcotqYCkFFymMjgd1GpXGtg2OgREQQFenPi1GVaNg8gKrLiATS5ucU8M38J8fHlqzbFxMTxztvT8A+49X8TvVHHK189wLnjlzB5udGkpYJq396EpPQcnvxgBZdSywcl7Y9N5MOnJ9PAS1ktV6mXsnjx/i+4cDoZlVrF0IndefTNKSIRO0lxqZl/LFjGqYQ0AHYcjuftR8bSvooFUd7/bBNLVx1EkkCrUZOdU8RdU5wz1SpmfzyvL1hFbn4JAKfjUnh1/ria+7ev6MlXIR811/pGrwEe7kZmTY/k/pn9CQ70vvEJTpJYsNqegAFSimIoMF+utevXph5dQ7jnjshKEzDAli2x9gQMcDExk02bYp0Wh0arIbRz8zqbgAE27DtrT8AA8UlZbDlwTsaIrm/54h1cOF3ewiTZJDYtO0hyYqbMUbmOg6cu2hMwQHpOIVsPxlX5/EPHEu1dcRarjf1HEpwW25Zdp+0JGODAsURy8+pvOU+RhBXq6qZntUqPVm2UKRpl8PYycfXNsodn/f6dXM33qlHAarUKP2+TTNFUTLpqsIUkSdhkGg/gioIaeGLQOTZ0erlX/W/FdNUYCTej87rC9HrHuNwMOgw1PAtByUsZiiSsUO187qGBoSMAavS09p6Mm7Z+150dMLAdQ4d2RKtVo9GoGDgwjBEjqt5XlZSYweuPfcf8Wf/j6w/WX5MI/qrMbOGVxeu4+7XvePyDZVxKz6nwWCUZFdmBkX3DMOg0GA1axg0IJ6pLK7nDusaoGX1p2upKi0P/kZ1o3EIZn+/s0iJeOrSKJ/ctZXnCEbnDqZbWTQOYMbIbXu5GDHoNg7u3YerQrlU+/64pfWkU7INKBS2a+XH31Jtrij59LoU33lvDG++vJvZMksNrd0/pQ1ib8mZxD3cDU8f2wM2ov97b1AsqqbJvIgV78cUXefHFF+UOo0ZZbSWkFx/GqPXDx9BG7nAUQZIkLiZmYpMkmjf3r3I/kiRJPDHjE04eTgRApYJ7Hx/B5NnXr/z29neb+HHTYft2j3ZN+egfk2/9f6AWSJJEWlYBarWKAAWP8E5PymHrykN4+pgYOqkHarX8zwQ2SeKeHV+yL6O8+dVNo+PlLmMY3bT2p805Q25BMWVmK/4+7jfd51pUVEZaZj4NA70wGKr+JJySlsujzy8hKaX8xjU40It/vzSFpo2uDBAsLbMQn5COfwMPAvw8byquylwvL/hHDSdgwIhrjp2qyVVEDhEDsxRMozYS7N5b7jAURaVS0az5zVfVKsgr4WL8lT4ySYJzp5IqPP7qJ99LaTlIklQnBg6pVCqCnPjFVlMCGvkwaY6yKu9llBRwNPvK2Itiq5mY9Pg6m4S9PapfpMRk0hNi8rvp83bujbMnYICUtDx27DnL9PE97fsMei1hbWpnHrcKZTU/X03+W09BqAXungb8gxwH1gUG+1R4fOMAx9caB/rUiQQs3BpPnZEGesc+dB+d8vrUlaxhkDc67ZXUotGoCHLCDAZXJZ6EbyA/r5hF728gJ7uQth0acce9/cWXcR2kVqt56LkxfPGfdeTlFNE2vAl3/n1Ihcc/MjmKopIy4i5n4O/lzmNTB9RitIJc3LQ65oUNYuGpbeSZi+nSoCkPhw2UO6w6pW+P1kwY1ZW1W2KRbBJDosIYFFn5/P8apbDiHFcTSfgG3vzXr+zbVT60f8+OMwDMmCVWkKqLwnu05N3vHqzSsQadlhdnXduPJLi+cc07M7ppJ0qsZjx0BrnDqZPmzhrErDsikaRrR1oLjkQSroTVaiP+7JW5upIEZ08myxiRIAi1QatW46EWCfhWKGrEs4KfhEWfcCU0GjW+fo6jS318lVV5SBAEQajEdeYIK2mglngSvoEHHhnKp++tJzurkJatA7nv77fJHZIgOM22nWf45ocYiorL6BrRnEfnDkWtrpkxD/EnL/O/V5eTl1VI6/AmzH1tCtp6XDNYEEAk4RuK6BbCwq/ux2aTauzLSRDkUFBYyn8/3URaWvlCLZcuZ9O4kQ9TJ/a8wZnV88E/f+TUwQsAxB27iIe3idnPjq2RawmCAwU9+V5NNEdXkUjAgqtJT8+3J+A/Jafk1si1ykrNpCRmOF7rQkYFRwtC/eGSSViSJOKyM4nLzqy0NGFtsVisiohDqJzZYmX17pOs3n0Ss8UqdzhOk5lVwJr1xzh11nFQYaOG3rQMuVL4RKfTENb21gooFBWUkJyQgcXs+PvTG3Q0DHEsS9m4pTLKVArykySJvLKSmvmevN5awgr6Ona55mibJPHYltWsjD8FwOiW7fjPwNtRyzC3t6CwhFfe/I2z59Lw8TbxwKwB9OrRstbjEG7MbLHy6PvL2H2yvKzl6t0n+c/DY9Fp63af5akzKbz8fytISsrBzU3HXdP7Mn1yLwAMBh3P/GMUX34fQ0lxGT26tWD4kI7VvtaW5Qf4/LUVZKXl0qpjU+Z/dA/BTa9UXHr0nen875Xl5GYV0Ca8KXc9OeqW//+Euu9YZjLP7l3DxYIcmnv68mavUbTzdd5KZkqvmOVySXh5XCzLzl1Z3m7ZuVgGNglhfJsOtR7Loi+2sWf/eQCysgv5+PMt9OzeQhT7UKAVO0/YEzBAzIkEVsbEMr5/3SxX+KclS/eSlFReQrC42MzyVYeZOrGnvXuldasgXnlu3C1fR5IkvvvP76QnZQNw5nAC3/x7Nf9YcKf9mGatg3npywdu+VqCa3nr8GaOZ6cAcDQrmbePbOZ/A6fKHFXtcbnm6OzSkirtqw05uUWO2zlFmM2u08zpSszXWUbP4gJL69lsjo8AVqsNWw00+dmsNooKHP/OigtKnX4dwfXklBVfte3k72uFN0e7XBIe07ItoT5X+rlCffwZ01Kekmkdwho5DOhq3TLwmrU0XUFZqYXlX+3gl8+2UpBXdM3rFrOVuBOXSPvjKUmJxkZ2oFOrK/2hnVo1ZHSf9jJG5BzDbuuAzx/rCavVKqIi26DVOO/PPvViJpuX7Sc9KYcOPa90teiNOroNDHPadQTXFeHXyGG75dwcQAAAIABJREFU81Xbrq7GM8KJEye47777iIyM5J133gHg9ddf5+eff8ZkMvH3v/+dqVOn8tNPP/HGG29gMpmIjo7mqaeewmw28+CDDxIbG4tGo2HhwoVERFS+fmyAyYPPh0/gq9iDANzVvisBJnmWc5s8vgeSBMdPJuHt6cace12v3KW5zMIL93/O4Zjy0p5bVh3mtS/uw8unvKhJYX4JLz24mGN743EzGRg/uz93zhsmZ8jX5WbQsfDRCfyy9SgAkwZG4HYTy7cpVd9erXn9hQnsO3ie4CBvht3mvG6ZXb8f4b/P/kR2ah6+AZ7MeXECTVoFkZ2WR8derRg8oYfTriW4rhe7D8fXYCIhP4tWXv483DHS6deot33CFouFBQsWEB0dTWpqefnH5cuXs27dOvbs2YNOV/4ll5ubywsvvMCBAwfQ6/V07dqVSZMmsWHDBgICAoiJiWHv3r3MmzePrVu33vC6TT29ebaX/EukqVQqpk7siSv3buxad9yegAHijl9m9Q97mPbgYAB+/HgTx/bGA1BcVMqvX2xn1LTeNAhU3qoqJqOeO4d3lzsMp2sf1oj2Yc5/ulj22RayU8unOGWn57Ny8Xbe+fVRp19HcG1atZonIurvAik12hyt1WpZtGgRDRteaeZbunQpDz/8sD0BA+zatYtevXrh5uaGRqMhKiqK3bt3s3HjRkaOHAlAz549OXHiRE2GK1SD5joVjzR/aYIvKzE7vFZaXEZ+bvHVpwh1kMXs2GducaFpXYILEX3CjpKSkvj2228ZMmQIQ4YMYd++fWRnZ+PhcaXJ2NfXl9zc3Gv2m0wmiovFF7iS9BnSgW5RofbtRiH+jJl5pTmp/8hONAi4ssB8p56taNLCn9pgs0n8+Os+3v1wHavXHauVa9Yn/UZ3xvDHCjl6o47I2yvvKhIE2Sg4Cdf6KCGDwcCDDz7I8OHDWb9+PfPnz2fevHnk5l6p1GOxWPDz88PDw8Nhv16vx83NrbZDFiqhVqscCjNkpuaya/1xBo/tCkCH7i149oM72bH2GEaTnikPDEJTS3NvP/h0I7+uPIgkwWrtUXJyC7ljcu9auXZ9MOH+QTRq7s+ZI4m0Dm9C3xEiCQvCzar1J+GOHTsSH1/eR2iz2fD09KR79+7s3LmTkpISrFYrmzdvZsCAAfTs2ZN169YBcPjwYcLCxGhLpcnNKuTM0Yv27dJiM4d2nnU4pn23EOY8M4a7Hh2OsRbXFj14JIE/Z+OYLTb2/VG3WHCe3sPCuevJUSIBC4qlquBHKWr8STg6Oprz58+Tn5/PhQsX+Oijj5gyZQo//vgjOp2OhQsX0qhRI+bOnUtUVBTu7u7MmTOHoKAg/va3vzFz5kwGDhyIwWDgo48+qulwhZtk8jDi6eNOcWGZfZ+7p1HGiK64ejFxN7e6P9pZEATXUuNJeMWKFdfs2759+zX7Hn/8cR5//HGHfd7e3vz22281Ftv1FBSU8Ml/N5CZUUDLVoHce/9ANFr5p1NLksTCX3dy6OxlPE1GHhrXl9CmNVN712aT+O+yHRw9l4yXu5G/T4ikRbDfdY/VG7TMnDeUb99fT35eEaHhzbjz0eE1EtfNunNqH97/ZCPJqbmENPPjzml95Q5JEITaprA+4Ku5XuWIW/T2a7+xa/sZAPbtPofNJjFnrvxrCH+97gCL1+yzf5bScwr46tnpaNTOv0H4/Pe9LF67376dmVfI4n9Oq7Dc5tAJ3RkwqjMlRaV4+pgUU5azb6/WRIQ3JT09n+Agb4xG8SQsCPVRvZ0nXBedj09z2I6PS5UpEkfnLmc43MwlpGaTnV+Mv7d7jVzrrxJSsykoKcPTzVDhOXqDFr1BeR8nd5MB9+YVxy0IgiAn+dtZFcbPz9Nh27eB85NcdTQJ8HHYbuzvhY9HzYwUb+TnWEijYQMvPIy1N6BKEATBaRQ+T1h5jy4ye+iRoSxcsI6M9Hyah/jzwMND5A4JgFmjepKeW8Cxc8l4mIw8NK6PU2sA/9WDY/uSlV/MiYRUvE0GHh7fTzFNzIIgCK5EJOGrtGnbkAUf3Y0kSYpKPBq1mmdm1s4NgU6j4YW7lVffWRAEoVoU9OR7NZGEK6CkBCwIgiBUj0pS9sAs0ScsCIIgCDIRT8KCUMdtSN3M1rQdgES/gL4MD1bGOAZBUAwFPwmLJCwIddiZ/DiWXFxKsbV8YZOUi2k0cWtEB+/2MkcmCEJViOZoQajD4gvO2xMwQImthAuFiTJGJAgKI13pF/7rj1KIJCzUCkmSkCQFffJdRDuvUDw1V+a2e2jdCfVsI2NEt8Zqs934IEFwIaI5Wqhx/1u4ka0bYtFoVIye0J2Jd4jlBJ0lxL05d7eYzqa0bUiSxIDAfrTxbCV3WDetpMzC/G/WcDwhGW+TkYdHRTKwY937/xAUSsH3/yIJCzVqy/oT/PJtDBZL+RPO14u20LFzM9q2byRzZK6jl19Pevn1lDuMW/Lhml1sOhoHQFpuIf9ZsY1+YS1qrCCNUH8orfn5auITLtSoy4mZ9gQMUFRYxrkzKTJGJChRZn6Rw3ZGXhEFJaUyRSMItUckYQWx2Wys+Xo7372zivOxl+QOxym69GyBt4/Jvh3U0IfuvVvW2PX2pVzis2P7OZyWVGPXqM9sNonDBy9w5HACNpvzHi86NQ9G+5cVwUIbBeBtUsa61Dcro6SAXy4c4mCmGCDnDJIkcSQ5mcNJSdUfVyJqRwtV8Z95X7L+hxgA1ny9jfmL7qd9z9YyR3Vr2oc35ZH5o9iw6ghqjZpx03oRGOxz4xOr4dtTh3ljzzbyzaX4GIw812sQk0M71si16iOrxcZLz/9CzB9LffYf0I5nXx6PxglNxlP7d8ZstbH/3CU8jXrmja6b9cpP56byyJ4lJBRmYVRrmRUayd/DBsodVp1lkyQeWbmKNWfKP3MjQkN5f/Qo1HXws1ERkYQVIi+rgF2rD9u30y9ns/a7XXU+CQP0GxRGv0FhNX6dJaePk28ub8LMKS3hx9NHRRJ2orVrjtgTMMD2rafYuO44w0Z2csr7zxzYlZkDuzrlveSy+GwMCYVZAJTYLPx0/gAPhPZDrxFftdXx64lYVp+58plbc+YMv55owcSOHar+Jgp78r2a+GQohEarvuaJQq1xnbs9OdTFJyklKy4qu84+0W/7V7arvu2tkg2rmJpXbfml136+rrfvRsTArHpsy6oj/Pvpn/jo1RUU5BVXeJy7l4kRM/uhM5TfFzVr25BJD4mVjG7GlNCOeOkMAPjojUwJDZc5ItcydEQnWrUJsm+3aRvMbcNr7ncccymRF3Zs5O3d2ykyX3sDoETjm0UQaPQAyr9chzYKw02rkzeoOmxs+zDCAvzt2+38/RnbvuZb1WqTeBKuQdvWHOX9fy21P0FcOJvKG1/MRq2+/r3PrBcm0mNYOMkXMug1PBzvBp7XPU64vhlhnQn19edwWjJdgxrTLUhMg3ImL283Xn9nGr/9egBUKsZN6o6HR80Mntp5KYF561eRWVI+avpwWgpfjZ6IpoK/HaXoHdiSj3rfwbbUswS5eTGuWYTcIdVpvm5ufD5hAl8dOgSouKtLZ3zd3G7uTURzdP11YOdZhya800cSyUrPxz/Iu8JzwvuEEt4ntDbCc0k9gpvQI7iJ3GG4rAZ+Htx934Aav87v8WftCRhgT1IiF/NyCfHxrfFr36r2vg1p79tQ7jBcRqCHB//o31/uMGqMSMI1yMPL8SnB09eEh+dN3sUJQj1k0jk24brr9HgaDDJFI9RlKkCl4H55kYRr0J0PD+XC2VROHUrAw8uNOx66DaNJb3/97IlL/PLZNsxmKwNHR9B/hHNGmQqCElgtNn7+NoaszAK69WpJz8iq17R+uGtvjqSlsDfpIh56PXM698TPzXTjEwXhaqI5uv4ymvS8+um95OUU4WbSozdcubvPyy7k/x7/gcsXMgA4ujceL193InqJermCa3jrpWVsXnscgLW/HebR+aMYOKxqU8Y8DQa+GTOZ8zlZ+BjdCDC512SogiAbZY9ycAEqlQpvX3eHBAxwdO95ewIGKMgt5vCuuNoOTxBqRFmphcP7ztu3iwpL2bn19E29h1atpk0Df5GAhVsmljIUrtG8dSDunlf6jFVqFYFNlD/oRBCqQqNVY3DTO+zT6zUyRSMIyiWao2XStFUgdz4ylBVf78JcZqXX4DBGTOohd1jXlZKcw/v//p2U5BwaNfbl8adG0cDfQ+6wBAXTaNRMntmHLz/ZTF5OMa3bBjNjdpTcYckuOT+PN3ZsI7ukmI4BQTwZ2d+lSjACWCxW3vxyIyfiU/DxMPLQ5H6Et5ZxuqDoExYqMvbOSKJn9kWSpArnDivBwv+sY9/ucwBcTMhk4Xtref6ViTJHJSjdmInd6TewHWmpeYS0DMBgFEUrnlj/O7svXQRg58VE9BoNj/WJlDkq5/pk6S5WbDtu337rq0189dIM2SrYlY+OluXSVSKSsMxUKpXiyyumpeU6bKen5csUiVDX+Pp54OsnWk0AzFYr8VlZDvtOZ2ZUcHTdlZTu+H2RnJFHYUkZHm5iitn1KPfxq57Kzy3ickIGVotV7lDsmjbzq3RbEIQb02k0BHk43pAEe7jeDUqLxv4O280b+uJu1FdwdC243jKGCnoyFk/CCrLih9189+kW8nOKCevUlOcXTMfbR/6RoY/983YMRi2pyXk0auzDQ48Od3h9yaoDrN5c3vw0alA4k0fV7ZVwBKGmvDRwMG/s2EZmURFhAYH8M9L1+slnRfeiqKSU4+dS8PYwMndy5ctS7t8bz5efbaWwoIQO4U159Mnb0Wjrz/OhSMIKUVpiZsnn28nJLATg+KEEvvloM3Pnj5Y5MnB3N/LkM9HXfe3g8UQW/bCDomIzAJ/+uIM2LQLp3F6UjhSEq3UObsSPk6bJHUaNUqtVzJtWtdKmZWUWFi5Yx6XETAAuJmYRGOTNnbOcW6ZSyX3C9ed2Q+GKi8ooKihx2FdSrPyVY07Hp9oTMEBRURln4lNkjEgQhLoiN6eItFTHPuSUlByZopGHSMKAJEnklhVjsdlki8Hb10SHriH2bXdPI72i2lZ6TklxGefPplJ4VfKuTZ3bN8HH60o9bB8vNyLaN5UtHiWRpBIkSfk3UoIglwYNPGgecqUPWa1WEdrWyYtfiD5hZcspLeaxfT8Rm5OMr97E4x1uY1jj9rUeh0ql4rl3pvH9oi3k55XQI7INvQe2q/D4U8cv8e7zS0mMTye4sQ9/e3r0DZN2TQhr3ZDHZt/GbxuPggTRQzrRtmXQjU90YZIkIeW9CqVrABWScQJqryfkDksQFEejVfPP58eyeNEWigpLCe/cnOgJ3Zx6DTFFSeEWxG5kd3p5eb08cwkLYjcxpFGYLBPoDUYd9/x9aJWO/f6TLSTGpwOQcjmH7xdtlSUJAwyJbMeQyIpvGOobqWQ1FH8H/DHCvegLbIaeqA2uuxybIFRX8xB/XnhtktxhyKbeN0fnmh2bcnPLiim1WmSKpuqu7i8uKRLNnophvYw9AQNQBpZEuaIRhPpNkq7/oxD1Pgn38g9Br7pS07ajbyPctMqv7NO1T2u0uvK4VSqI6NFC5ogEO+MgUP+lX0vdDAyD5ItHEK7DJtk4kH2QPZl7sdiU/+Dhqup9c/S0lt1Rq2BfRiK+ehOPdBgsd0hVMnV2FD4N3Ik7lUyjpg0YN6OP3CEJf1Br22Dz+TcU/QBowHQXaq2MtXMF4So2ycbCuI/Yn30QgDDPdjze9lH0auU/gFSH6BNWuCktujOlRXe5w7hpw8d3Y/iNDxNkoNZ3B33d+0wJ9cOBrIP2BAxwMv8Um9I2MyJ4mIxR1RCFjYa+Wr1vjhYEQahvyq4zdU40SctDJGEBgPSiAr6OPcTGxDi5QyH2bDKrNh0jPatA7lAEwSX1aNCDtp5t7NvNTc0YEOCao/dVEqhs1/4ohWiOFjiXk8mcDcuJy8lEq1IxM6wLL/W9TZZYvlm2l8U/x1BcYqZxsA8vPjKK9m2cPHlfEOo5vVrHE6GPsSFtE1ablYGBA/DUecodVr0knoQFPj9+kLic8tqtFknil7jjZJUU1XocNpvEsrVHKC4pL4N5OSWHH1ceqPU4BKE+MGgMjGo4kujGo/Fy9QQsKmYJSiZd9Ym02iSsttr/lNokCavNcQlHq4ylRAVBcAGSGB0tKNyMdp3ZfvkCifm5qIDRLdsSYKr9JRS1GjUDeoXyy++HsNkkGnibGDmwQ7Xfz2qzsfD3XZxMSsff08QTo6PwcXe78Yku5uejx9l+PgF3vZ7H+vclwMN5/7bL9p7g5OVUmgf4Mj2yc6VL1gn1R25uMb8u219eyja6Cw0auN66yc4ikrBAB/9Avh4xmTUXzhBocmdC6+onvlv1yL2DaNcyiOT0PHp2bk6HNtWfX7vw910s2rTPvp1dWMx/Z41zRph1xrITJ3lx/WZKLOUjX89lZvHdHZPRqG+9J+p/m/axcO0uzFYbahVczsrjyeiqLWEnuK6CghL++fSPnDlTvprazl1neeft6Xh7y3gDXEmFLEmSeOaZZ9i4cSNarZaXX36ZIUOG2F9fuXIlb775JgCBgYH88MMP6PV6YmJieOKJJ9BoNISEhPD111/z5ptv8sMPP+Dj44NarWbFihV4eFR+AyKSsABAiLcvf4voJXcYqFQqRtzC0+9fnUxKd9g+m5yBJEn16mltd8JFewIGOJGaRnphEcGet/5ksuPUeczW8u4CmwS7z4rSnAJs2hRrT8AA586lsX79cSZN6iFjVBVbv34958+fZ+/evVy+fJlBgwZx+vRp+/fEM888w7Jly2jZsiWjRo1ix44dREZGMm3aNFauXEl4eLj9vUpKSliwYAEDBw6s8vXFwCzBZfl7Oja7Bnp51KsEDNDA5Pj04Wdyw9tocMp7u+kdqysZdeKeXgA307WfL6NRvkpcKun6P3/auHEjI0eOBKBx48YYDAaSk5Ptr/fo0YMvvviCQ4cOkZ6eTpcuXdi1axfh4eEOCbi6XP6vxmyx8tyP6ziakIyXycDDw/vSP0zUWa4Pnhjdn6yCIuJSMgjwcuepsfWvqfSRfn2Iz8rm4OUkPPUG5vXrjZvOOV+Ic4b04lJmLufTs2no68mc23o65X2V5rf9sfyw8wgWm43hEaHMGqzMJzqluG1we3buPMO2bacB6Nu3NSNHdqrw+OTkHL5YvJ3iolK6dWvBuHHOXcrwRrKzsx2ajH19fcnNzaVRo/KusDfeeIPIyEg++ugjXn31VXx9fUlKSiI5OZlRo0aRm5tLdHQ0Tz31FIGBgcyfPx+Anj178vbbb6PX6yu9vssn4Y/WxbD60KnyjSx4a8UWeoc2Q6fRVH6iUOf5uLuxcPa4etcE/VcGrZaPJ0RTbDZj0GqdukRn55BG/PDoHSRk5NC4gRdebkanvbdSxKdm8vaKbWQXFgNwLjWLkIAGDA5vJXNkyqVWq/jX8+M4diwRq1UiIqIZGs31G12tVhuvvLqcU6fKnzz37T+Pm5ue4cNv/QnTQSWjoz08PMjNzbVvWywW/Pz8ALDZbIwePZolS5bQpk0bJk6cSIMGDTAYDLRt25Zvv/2WsrIy2rdvz/Tp03nooYd46KGHsFgszJo1i8WLFzNnzpxKQ3P55uiUXMeqS2l5heQWlVRwtOCK6msC/is3na5G1sg2GfSENQ50yQQMcCQh2Z6AAUrNFk4npckYUd2gVquIiGhO164hFSZggPSMfOLiUu3bZWVWjh+/5NxgbtAc3bNnT9atWwdAcnIyFouFwMBAAHJzc0lKSiI8PBwPDw+ioqKIi4ujY8eOnD9fvg69SqVCpVLh5uZGWVl5OVCtVovJZMJkMt0wPJd/Eu7YLJjVh07Z5722DvKngfuNfzGCIAjdWzYhwMud9LxCAEwGHZ2aiwpuzuLt5YafnwepqXn2fX5+tTudaeLEiaxevZr+/fuj0+n48MMP+emnnzCbzdxxxx088MAD9O7dGw8PD0wmE9988w0NGjSgf//+9OjRA6PRyPz58/H392fWrFmcO3cOSZKIiIhg+vTpN7y+yyfh6X0jKC41sz/+El5GA38f2Re1WjwZCYJwY039fXhuwm18t+MQFquNEV3aEtkuRO6wXIabm54HHxzMl1/uoLCwlPCOTZg5s6+TryJVOkVJp9Px1VdfOezr3v3KCmjPP/88zz///DXnvfXWW9fs+/zzz286OpdPwiqVitmDezBbDKYQBKEaBoe3En3ANWhAVDsGRLXDZpPq5QOSyydhQRAEQflqKgFf3QesNC4/MKs25OcVU1RUKncY9VJBWSk/nD7KqvhT2CppchIEpSi1mrlUlE2p1Sx3KPWHWMDBNdlsEu+88Rsx28+i0aoZObozsx8cJHdY9UZOSTF3rf2JIxnl1XnGXGjHB4PGiNHQgmIdzkrkpaPLSSzMoKm7H8+HR9PNL0TusAQZiSfhW7Dmt0OsX3OMgoIScnOKWLpkL8ePXpQ7rHrj8xMH7AkYYOX5U+xPvSxjRIJQuY/PbCa+IA2LZON8QTofn9ksd0iu7wZTlOQmkvAtyM4qdNguK7OQkpQjUzT1j+2qNiUJMF+1FKIgKEmBxbFGQaFFdGPVdyIJ34I+/UPx978yp61Zcz969BGjKGvLzHadaevrb98e3KQlvYKbyhiRIKeM/EIS07MVvQZ1D78WqCnvLlEB3fyayxtQfSBRvsLI1T8KIfqEb0GTpg3w8jGRkVFelcvNpMfNWHmdUMF5gt09+XbEFH6JO4FJq2N6uwinLNEn1D0fro3hm+0HKS610KtNUxbcG63IBSXmtRuKv8GTuPw0Wnr6M7OFs+fECtelnJx7DeV9SuuQ31ceIT7uSgm70yeTWbvmCGNquQB5fRZg8uDBTvIvwSjI53JmLl9uOUBRWflo452nE/h80z4eGt5H5siupVKpmNFSeXEJ8hFJWBCEOi2roMiegP9UWFomUzSC0ihtINbVRNvdLRg5ujOdOjezb0d0ac7w2yNkjEgQ6p92jQPp0qKRfdvf08SQ8NY39R4pSdmsXLKXYwcuODk6QaiceBK+BXqDljf+PZ3fVx0GlYoRt0eg14tfqSDUJp1Ww39njWXRxr2UmC2MiGhLlxaNq3z+8YMJvPXsL6Ql52B00zF1dhTT76t/a0+7rsprR8tNZIxbpDdoiZ7Q/cYHCoJQY7xMRp4YE1Wtc1d8v5u05PKphSXFZn7/9SDTZkeJoi9CrRBJWBCEes121XQVyWZDkiSRhF2FwvuERRIWhGranr6TswVxNND5MqbxKDQqjdwhOd3a34+ycWMsGo2KceO706uX682DHz6+KycOJ5KdWYBWp2HAsHDUYqrbTUmIS+XH/23FXGZlwIhw+g3tKHdIjupzEj5x4gT33XcfkZGRvPPOO/b999xzD2vWrCE1NRWAn376iTfeeAOTyUR0dDRPPfUUZrOZBx98kNjYWDQaDQsXLiQiQgx8EuS3NmU9SxJ/xoIFgNTSNB5odZ/MUTnXwYMXWLhwA4WF5VWd4uPTee/9mQQH+8gcmXP16BfKK/+dyYHd52jctAH9hnSQO6Q6pbCghNf+8QOJ58qnax7aHYeHpxude7veDVtNqNEkbLFYWLBgAdHR0fZkC3Dw4EFKS0spLS3/487NzeWFF17gwIED6PV6unbtyqRJk9iwYQMBAQHExMSwd+9e5s2bx9atW2syZEGokhO5sfYEDHAm/6zLNWEeP3bJnoABMjLyOXQogZEjXSsJA7QOa0TrsEY3PlC4xomDCfYEDFCQV8LB3XGKScLlU5SU+yhco20uWq2WRYsW0bBhQ4f98+fP56WXXrJv79q1i169euHm5oZGoyEqKordu3ezceNGRo4cCUDPnj05ceJETYYrCFVm1Bgdtt00RpdKwAAhLfzR6680sXt6GmnXrmElZwj1UZMQP9w9DQ77/IO8ZIqm7qn1jo/Vq1fTvHlzQkND7fuys7Px8LhSg9nX15fc3Nxr9ptMJoqLi2s1XkG4nvGNxxLiHoIaNQEGf8Y2jpY7pCopKCrlnx+sYNI/P+fBN5dwOjGtwmOjotoR2b0lHmYbnhYbgyNDadEisBajFeqCRs38uWvuEIIb++IX6MnwCd0wuul57YnveOfZn0m+nCVvgBJgu86PQtTqwCybzcYLL7zAsmXLHPZ7eHiQm5tr37ZYLPj5+V2zX6/X4+bmVmvxCkJFGroF83zYfLLKsvHWeWHQGG58kgK89+M2Nh+IAyAxNYd3vtnEomemXffY0ycucWxHHKV55Sv/7FpzjFFju9KqrXgaFhyNndGXMdN7Y7NKxGw+yb+f/5mSovIqZonxafz7qznoZKvlLdXf5uirZWVlUVJSwuzZsxkxYgSFhYXce++9dO/enZ07d1JSUoLVamXz5s0MGDCAnj17sm7dOgAOHz5MWFhYbYYrCJXSqrUEGgPqTAIGSMvKr3T7r47tv0Bu9pXlOnOyCjmyN77GYhPqNrVajVan4ei+eHsCBjh74jLJFzNljEzZavzWJDo6mvPnz5Ofn8+FCxc4duyY/bWAgAC++OILAObOnUtUVBTu7u7MmTOHoKAg/va3vzFz5kwGDhyIwWDgo48+ckpMZrOVjMx8fH3cMRp1TnlPQagLWjb2I+bYBft2i0Z+FR7bqm1DjCad/QvVzaSntegTFm7Ax8/DYdvX3xNfPxn7iCXq9xSlFStWVPhaenq6/b8ff/xxHn/8cYfXvb29+e2335waz4XEDF5/exXx59MJDPRk7pzbiOx9c3VmBaGumju5PzZJIu5iBn7eJh6dXnF5xi69W3Hng4NZ/9thAIaN7UKnHi1rK1TZnI5N4tMPN5KTXUTr0CCeeHo0eoMoqVBVU2cPIDE+naP74nFzNzD9voF4eotuxIrUu0/WF9/s4Exc+XSppORcvvx2p0jCQr2h1ah5bPrAKh8/8e47D3JpAAAdY0lEQVR+TLy7X80FpEAf/Od3Tp9MBiAxIQNvHxMPPTJM5qjqDq1Ow/y3pmI2W9BqNcqYNaDgPuF6l4SLihyXOCsqFkueCYJQzmKxkpaa57AvLTW3gqOFysg3EMuRWMpQYbp1bo5Wc+V/O7xDExmjEQRBSbRaDU2bOfaTt2gppmVVh6Tgp08lUcatSi2aNqkXHh5GTp1OJjDAixlTe8sdkiAIwDc/7+b3TccBGDG4IzMnyfO3+dRz0XyycAN5OcWEtmvIzHv7yxJHXWW12Xht8Xr2n7qIyaDjnlE9GdFbzpktYilDxRk9IoLRI0QNakFQin2HzvPljzGUlpWXAv3yxxjatg6mR+eQWo8lKNibf70ysdav6yq+W3uA33ZcqW74/pJt9A1vgZe7sZKz6q961xwtCILynEvIsCdggNIyC/EJGTJGJFRXylVzz9NzCknJrHg+eo2TQGW79kcp6mUSTkzL5sv1+1mz75TotxCqLSO/kNhLqZSaLTc+WEZ7jiewcMl2ft545Iafd6vFyi8frmfxq79y6kDtFeboHtGcBr7u9u0Gvu50i2hWa9evj2w2icuZueQWltz0uQXFpfxvzR4+XbWbzLxCh9ciWjfCoL/SyNqmqT/NXGzlLWeqd83Rxy8k88/PVpOUlYdGpeLwuSTmTxssd1hCHfP9zsMsXBtDblEJHZsG8e7do2noo7yi9ev3nOLNxRvJLypFBZxNTGP+vUOve6wkSfzfA/9j2/L9AKz9fidPf3wfEf3b1XicrVsE8uRDw1i57igAo4d1onWIGBBVU4pLzTy2aAX74y7hYdRz75Ae3D2ke9XOLTPz8Ae/cjS+fBrXxkNn+XDeBPy8ym+ihvVqR35RKbuOnceg13LfmN4Y9TIWRZIQfcJK8sv2YyRllU9BsEoSvx84zcNjI/F0qzulBwV5ma1Wvtiyn9yi8ieI4xdT+XT9Hl6YfP3kJqcNe86QX1S+HKEE7Dp6HqvNhuY6i9bnZhaw/4+BUQDZqXls/mVvrSRhgMierYnsKebs14ZP1+5m9+lEAHIKS/h8/T7G9GpPA0/TDc/dcOCMPQEDnL2UwW8xsdwzvId938RBEUwcpKBxN8rNwfWwOfqqeePqa3cJQqXMFivFZWaHfcUKbZLW6TQO23qdFnUFxRP0Bh16g+MTi/aq8wXXUFji+PktKC4lr4rN0kbDtU+1f21+Fm5OvUvC0wd0pmlAef+EVqNmdK/2eNSRp+Ds0svsy/ie49mrkSQFjSyoZ0wGPT1aNbVvexoN3NZBmU9wd4/uSUjDBgB4mgxMH9a1wgpGJk8jo2cPxOimB6Blx6ZMfXRErcWqZHHHLrL5l71kp+fd+OA6YFB4S3z+Mlq5S6vGNAmoWr/tbZ3bMKRrG/t2n7DmTOwX7vQYnUUlSdf9UYp6d/sS2jSQTx6dyNYj52js503/8LpRCze9JJ6Vl14gz1LeDHSp6AjDGz2tjJJw9dD/zRhJu0YBZBcW07dtc/q3ayF3SNfVpmkAn/9rOkfOXqZ5wwY0Caz8i3bmk2OIHNWFtEtZdOobipuHmFby88L1fPef1RQXlNKoRQBPfzKbNp3q9qCxPmEhvH73SDYdicNk0HH/8F4ORYwqo1arePO+Uew+mYDZYqVvxxB0GtFiUl31LgkDNPT1YtrALnKHcVOO56yyJ2CAc/nbKbDch6cuQMao6i+dRsOcIb3kDqNKPEwGIiOqfrPZon0TWrQXleSgfA30VV9uo7igvF896Xw6Sz/ewD8/nCVzZLeub1gIfcNCqnWuWq2ib4fqnSsLBT35Xq1eJuG66OonXpVKg1pV73oTBKFWSTYJi9nqsM9qEV1BdYoEKPifTCThm3Q8K4Wl8cfQqTU82KEPvobaWaKrS4NJXCo8TGbZBVRoCPMehru24rVgBUFpSkvN/PrTPsrKLAy/PYKgYO9au/aKtYfZuO0UGq2aSWO60bd7qyqdp9Fq6DMygpWLtyHZJLz9PBg8sWe147BarPzyySayUnPp3L8tvYcqty9VqB0iCd+EU9lp/G3rUi4Vla+qsjctke+HzMCorfk5cN66YCY2+zfnCnbioQ0gxKPHjU8SBIUwm608+48fOXIoAYAtG2N5/Z1pBDes+SIOew6e58MvtlBUXD4i+MLFTD588w6CA6t2E/C316bQOrwZqZcy6T6wPWG3sKbyu499w6al5fOw1y3Zw8OvTWHwRPG3XKMUNhDraqI98yasSjhpT8AAhzOTiUlNqLXru2m96ehzu0jAQp1zaH+8PQEDXEzIZO2qI7Vy7ROnk+wJGCAjs4AjJy5V+XyVSsWw6X2488nRt5SAzWUWDu04Y98uzi9h9/pj1X4/wTWIJHwTPHR6h22dWoO/0b2CowVB+JPRTc/VQxh0tTS3tGVzf/R/me/s7elGuzbBtXLtv9Jo1RhNjt8hBqOMlaTqE0m69kchRHP0Tbi3XU/2pF1ka9I59Botd7TpQrhfQ7nDEgTFC49oxpDh4WxYewzJBp06N2PClOr3rd6MgX3bknApiy27TqPTaJgwqgvNm9TOeApJkvj6ndWcPngeTx93hk/rw6+LNpObWUCrDk2YXsk8bEmS+Pq9dRzYfgaDUceUBwbRPaptrcTtUipMusqY3imS8E3QazR8NnAycbkZmLQ6mniIouSCUBUqlYonnxnDsJGdKCospXuvVuhrscrS3VP6cPeUPrV2vT8t+e96fnjvd3sOyEzN5eNNz5B2OYvmoQ0xuOkrPHftz/tY8slm+2js9KRsPlj+KB5etTMYVKgdIgnfJLVKRaiPmJsrCDdLpVLRuWuI3GHUqvjYyw4PYYlnUjC46QmNaH7DcxPOpDhMh0q5lM2l82m0q8K5wl8ofIqS6BMWap0kSUi2ArGMpODy/IIcR2D7N/S5pl+4Ii3aBqPRXvmKDm7agKYtxcpSrkY8CQu1ylZ2EPJeAmsyaFsgeb+FWivu7F3RyeOX+P7rXZSWWegT2YZxk+rfqP57nh5DVloeZ48m4N3Ag3ufGetQeGfnuuMs/3onVouNqNv/v717D4+qOvc4/p1MJpcJNxOCgAQhgLlAkMQYGgQCokeBAnJRCgTOwQIFVKTqA1J7FPGKpdbSioJoK5yISq03EAtyvwSBgFzlYoBIACMCiSHXmcw+f0RHEgIkwGRPkt/neeZ5WHv2zH5dhrystfZeb0cGjLzN/d5/DUng+xM5bF9/AL9AG0PH305QfU1FV5UFLvKIktaEpS7KnQXOr0v/7NgBuX+C6/5ubkxyzeXnF/HS85+S+e0ZAPbsOsZ1wUEk3R5tcmTVyy/AxuNz/qfC9747doY5Mz7mzKlcANK/PkHTFsF07hnlPid50p0kT/K+Eply7Wg6WqqX62y5drY5cYhHHU7/3p2AAYqLnOzbc9zEiLzPnrQj7gQMUFTg4ODuyj+/LJVU0eNJXrQUpiQs1cvWoVzbiwp/yzUTFhZCcON67rbFB25ocZ2JEXmf6NgbaRTySx/5BdhoG928St9hGAZOZ8nlT6zrvDgJazpaqpWl4XMYPqFQcgJ8b8JSb7zZIXmcyzCY/p9VbDqaQYDNxu8Sb6VfdKTZYXlUw0Z2Jjx0J++npFJY5ODWhHD6DbzF7LC8SvMbGzP+iX588n+bcDpK6Hp3DIl3tK/057enHWXunJWcPZtHeHgT/vDkABro8aUaR0lYqpXF4oelwRSzw6hWb2/dwaKvdrnbL65ax22tbiTYXrt/YfboFU2PXnVrDbiqkvreTFLfK5sNmvfaStK/yQLgzOlzvPH6Kh6d0vdahlc7eNnItzxNR4t42PEfc8q0vz+Xx7FsrYXLlSspcZGdnVfmWE52vknRyNVQEhbxsLgbmhNo+2XSqV3jENo1bmxiRFJVa748yDufbiXj+GmzQwHAavUhvM317rbFB6KquJ5cZ/y8WUf5l5fQdLSIh/WJiuBsQQGrvzlCoM3GhC4J2P20cX9N8bcFa1i8bAfOEheLl+3gqYd60ykqzOyw+MP/DmD+3NVkn80jMro5Q4dV/7acNYU3lzJUEhapBiPiOjEirpPZYUgVOZwlLN+4H2dJ6dAp64dcPlqx2yuScL16AUx+tLfZYchVUhIWEbkES7mdlSwW7x1VSQV0Y5aIeFL2uQJSVm3no017KXGZu9i1JuMwf9uWyvpjRz16HcMw2LzjCMvW7CEvv8hj17H5WumTFI3Nt7QecbPQBgy+K9Zj15O6RyNhkRrsVPY5Hnz1Qw4e/wGA9XsO86cxv8bHp/r3xX1z5zZe3rKRPIeDIJuNaYlJJHfwzBT8S68vZ8nqPbhcBhHhO5j1h0EENwryyLXGD+9GTOQNZJ48S9f4NtxwvUqY1igG4NJIWEQ84L11O90JGGDVV9+w6+hJU2L5+NB+8hwOAPIcDj46uM8j1zmaeZrP1u7F9dMv1gOHs3hvyTaPXOtnt8WFM7TvLUrANZK2rRQRDyk/3rVUcKy6+FrKXtlq8cyvF6ezBFdJ2Wn3Ei8e6YhcipKwSA32mx6diGgR6m7fEXcTHVs3MyWWke07ERxQugtY40A7o2I8s3Ya3jKUrvFt3e2wZtfR7/YYj1xLagEDrx4Ja01YpAYLaRDE3EmDWbbtAPUD/eh9a1SZerXVaWBke6JDm5D23QkSmrWgbXCIR67j42Phmcf689F/dpJXUMQdXSM1TSw1lpKwSA3XsF4gv+nhHc8gR4SEEhESevkTr5Kv1YchfXSXslSGd418y9N0tIiIiEk0EhapZVZsP8jaXenYA/z4XZ9fEdLAM4/uiNQIXv6IkpKwSC2yemc6M1JWcK6gGIBDx3/gjcn34mvVpJfUUYYBhhdVbChHfzNFapFN+464EzDA7iMnOXE65xKfEBEzaSQsUos0tAeUaTcKCqRRvUCTohHxBt59Y5aSsEgtMrbPrzh04jTbv8mkXoA/9991Kw3KJWYR8R5KwiK1iL/Nl1fG9ye3oIgAmy9+Nv0VlzpON2aJSHWyWCwa/Yr8zMt2yCpPN2aJiIiYRCNhERGp3TQSFhERkfI0EhYRkdrLy9eElYRFRKT2MgxwaccsERERKUcjYRERqd28eDpaI2ERERGTaCQsYoKiEifrvjtEA78AEhq3wmKxmB2SXERuYRHzt2yjyFnCgPZRRF0fanZIUhW6MUtEzneuuJBxqe+w48wxrFi458ZOPBPbT4nYCxWXlPC7Dz5m67HjAHx+4CCvDx5AZBMlYrk2NB0tUs3+mb6ZHWeOAVCCwccZX7Ev+6TJUUlFdmSecCdggOM5uSzZd8DEiKTKft47uvzLS2gkLFLNHK6SMm0nBvklDpOikUtpZA/EZvXBUfLLIy4BKopRwxgYhh5REpGfDL4xllZBIe72bU3CiQ0OMzEiuZiI0MaMiL0Zf18rAF1uDGNMQrzJUV2cq2AprtxXcBWuNjsUqST9k06kmrWsF8y8LsP55Ngu/K02RrXpjK+P/j3srZ64owfDYjuSW1RE+6bXe+3/K1fu3yHvdaAYsOOq/xg+Qclmh2U+w7umn8tTEhYxQVi9YB6I6mF2GFJJ4SHBZodweUUrKE3AAPlQ+DkoCXs9JWERkVrBVrZpsVV8Wl1joEeURGq7Qz9+w0cnPqawpJCoBlEMbjFQjxxJ9QoaBT++AMZp8GkG9v8xOyLv4OV7RysJi1ylYlcxbx39JycKSx8zSs87QkNbQ+5s2svkyKQu8Qnsj8vWCRx7wBaHj29Ts0OSSlASFrlKZ4vP8l1hlrttYJBZcPwSnxDxDB/fluDb0uwwvIx375jlnbf5idQgjWyNaBLQpMyxZgEahYjI5Xk8Ce/du5fExEQee+wxABYsWECPHj2IjY3loYcecp+3ePFi4uLi6Nq1Ky+99BIADoeD3/72tyQmJtK1a1d27tzp6XDlCrhyX8F1qjeuH/rjyn/P7HCqnb/Vn9GtRhJVP5JW9lb0bnoXdzW90+yw5BpZuOsrfv3uQvouWsD8HdvMDkeqyDAMDJfrgpe38Oh0tNPp5JVXXqF///5kZZVO13Xt2pVRo0bhdDpp06YNY8aMoVWrVjz11FOkpaXh5+dHXFwcQ4YM4YsvviA0NJTU1FS2bNnCpEmTWLt2rSdDlipyFSyFvDeAn3Z8yp2Fy3YLPra2psZV3SIbRPJ4g0izw5BrLO3kcf6Uup7c4tJHfzK+3ERESGO6tWxlbmBSeV5+d7RHR8K+vr688cYbNGvWzH0sPDzc/R6A3W5n06ZNdO7cmcDAQKxWK927d2fz5s2sXLmS3r17A5CQkMDevXs9Ga5cCWc67gQMYOSAc49p4YhcS7uystwJGCDP4WDfqe9NjEhqG9PWhFNSUoiMjKRdu3acPXuWevXqud+77rrryMnJueC43W6noKDAjHDlYvzigfq/tH2aga2zaeGIXEu/ahFG40C7u90oIIBbm99gYkRSZUYFxRu8aActU+6O3rZtGy+++CLLly8HoF69euTk5LjfdzqdhISEXHDcz8+PwMDAao9XLs7HvwuuBk9A4RLAF4L+Gx/fZpf9nEhNENU4lBk9evHOnp0YBgyJbk9cMyVhuXaqPQnv27eP0aNH8+9//9s9TR0fH8/vf/97CgsLsdlsrF69mocffpjDhw+zfPlybr/9dr766iuioqKqO1ypBB/7ILAPMjsMEY/o3fYmere9yeww5IoZ4MVVlDyehPv378+RI0fIzc3l6NGj5OfnU1BQwNixYwFITk5mzJgxPPDAA3Tv3p2goCDGjRvH9ddfz4QJE0hOTqZHjx74+/vz2muveTpcERGRauPxJPzJJ59U6rxHHnmERx55pMyxhg0b8umnn3oiLBERqQsMMLxoDbg87ZglIiK1l+Hd09HaMUtERMQkSsIiIlJrle6YdeHr/PenTZtGQkICXbp04YsvvrjgO4qLi4mIiGDo0KEAZGRk0Lt3b7p160ZiYiJHjhy54viUhEVEpM5asWIFR44cYcuWLSxevJiJEydilNth69VXX6VTp04UFRUBUL9+fRYtWsT69evp2bMns2fPvuLrKwmLiEgt9tOacPnXT87fmfGGG27A39+fkydPut/Pzs5m0aJFjBkzxn0sODiYRo0aAaW7P9rtv2zoUlU19saszMxMpk+fbnYYIiLiJTIzMy84Nu6FURQWFl5wPCAgAOCiOzY2b94cgOeff57Jkydjs9ku+I4zZ86wcOFCVq5cecUx19gkPH/+fLNDEBERL/f4449f8v2L7dgIcOzYMdauXcvMmTMvKB5UVFTE4MGDee6559w1Ea6EpqOvUkUL9CdOnKBXr150796dfv36kZ2dDVRcrrGuWrlyJTabjc8++0z9VQmpqal06dKFbt26MXLkSPXZZRiGwQMPPEBSUhJxcXG8+eab6rOLKF9utir9VBvKzSYkJLi3UD558iROp5MmTUrrg+/atQvDMOjduzdTpkwhNTWVv/zlLzidTu677z4GDx7M8OHDry4AQ67K6dOnjbNnzxqGYRjTpk0zJk+ebAwbNsxYtmyZYRiGMXPmTOPJJ580srOzjaioKCM/P99wOp1Gx44djfT0dDNDN43L5TK6detm9OnTx/jwww/VX5dRWFhotGzZ0ti1a5f7mPrs0tLS0oxu3boZhmEY33//vREWFqY+q4DD4TDGjBljPP/888bDDz9sGEbVfrbmzp1rTJ061TAMw/jyyy+N7t27m/bfcqWKi4uNkSNHGl27djV69uxpbN261Xj//feNlJSUMuetXbvWGDRokGEYhpGSkmIEBQUZSUlJRlJSkjF48OArvr5GwlepogX6devWceedpUXde/XqRWpq6kXLNdZFCxcupGfPnoSGhgKovy5j06ZNxMTEEBMT4z6mPru0Vq1akZWVxebNm3n77bcZMGCA+qwCFZWbrUo/1YZyszabjQULFrB+/XpWrVpFfHw899577wUj3O7du/PBBx8AMHz4cM6dO8eaNWtYs2YN//rXv674+jV2TdjbnL9A/9prr2G1WoGLl2X8+XhdU1hYyF//+ldWr17NpEmTAMjPz1d/XcKJEyc4efIkffv2JScnh/79+6vPLiM4OJgZM2bQv39/AgMDWbVqFSkpKeqzSqjKz9bFys2q2l3lKQlfA+UX6H18fHA6nfj6+l60LOP5i/91yezZsxkxYgQNGjRwH1N/XZq/vz8RERGkpKRQXFxMdHQ0OTk56rNL2LRpE//4xz84evQoX3/9NXfddZd+ziqpKv2kcrNXT9PRV6miBfr4+Hj3rivLli2jV69exMfHs3HjRgoLCykpKWH16tUkJSWZGbop9u/fz9KlS7n77rtZsWIFM2bMoFOnTuqvS+jQoYN7Rx6LxYLFYqFnz57qs0vYsmULN998M3a7ndjYWHx9fYmNjVWfVUJVfn+df1OTys1eGYthGN5bXqIGeOeddxg3bhzx8fEANG7cmOnTpzN27FjsdjvNmzdn3rx5BAYG8vLLL/Puu+8SFBTEqFGjGD16tMnRm+v+++9n4MCBtG7dWv11GVOmTGHVqlUEBAQwevRoOnfurD67hFOnTpGcnEx+fj4Oh4MBAwbQr18/9VkFzi8326lTJ5599tlK91NOTg7Jycnk5ua6y81ezeM6dZGSsIiIiEk0HS0iImISJWERERGTKAmLiIiYRElYRETEJErCIiIiJlESFqljNmzYUKY26uVkZma6tzEUkWtLSVjEQ4KCgmjbti1t27alXbt2fPvtt1X6/MaNG9mwYYO7/eSTT7JkyZKrjmvu3Lncd999lT6/RYsWFBYWcvDgwau+toiUpeeERTzEYrHgcDjw9b2y3WGnT5+O0+nk2WefvWYxlZSUEBYWxuHDh91FzSvj6aefJiAggKlTp16zWEREI2GRarV7924SExNp164dHTt2JC0tDYC1a9fSvn17WrduzZw5c3jvvfeYPXs2c+bMISIiguLiYsaOHcvbb78NQEREBDNnziQmJob4+HgyMzPd39OxY0dCQ0OxWq0MGzaszPX3799P06ZN3Ql47NixTJs2jc6dOxMeHs6KFSsYNGgQYWFhZZJ/XFwcqamp1dFFInWKkrCIB3Xo0IHIyEjGjRsHQNOmTfn88885dOgQU6ZMcRdHnzVrFrNmzeLw4cPcc889DB06lEmTJjFx4kQOHDiAn58fDocDh8MBlBYNycvLY/fu3cTGxjJv3jwApk6dyrvvvst3331HYmIif/zjH8vEk56eTsuWLd1th8PB9u3bWbduHRMnTmTIkCFMnz6drVu38sILL5CXlweUlgZMT0/3eH+J1DWqoiTiQXv27CkzHR0UFMSjjz7Khg0byMnJoU2bNgD07duX6dOnA7jrs17OqFGjgNIN97ds2QKA1WrF4XDgdDopLCzEbreX+cyPP/54QZWbe++9F39/f+Li4oiJiaFjx44AhISEkJWVRXh4OHa7nXPnzlW9A0TkkjQSFqlGf/7znykoKGDHjh289dZb/HxLxvjx45k/fz7z5s3jwQcfrNJ3Wq1WSkpKgNJScg8++CC33HILI0aMoHXr1mXOtdvt7tFteT4+Phe0XS4XAHl5edSvX79KcYnI5SkJi1SjM2fO0LRpU6xWK0uXLnUf37VrFzExMTzxxBNs3LgRKE2Yp06dqtL3WywWMjIyKCoqIiUlxb3m/LPw8HAyMjKqHPfRo0dp27ZtlT8nIpemJCziIXa7ncjISPdjSgcOHGDChAksWbKE6OhoAgICsNlsALz66qu0bt2a5ORknnnmGQAGDhzIypUruemmmygqKsJqtWK1WgEq/PPp06fJzMxk3759HDp0iHHjxvHcc8+Vial9+/YcO3bMXYj9Yt8J4Ovr625v376dhIQED/aWSN2kR5REaol9+/YxcOBAdu7cSUBAAE8//TRZWVnMmTOnzHlDhgxh0KBBDB8+vNLf3aVLF+bPn090dPS1DlukTtNIWKSWiI6Opm/fvkRGRtKuXTvS0tLcN3udb+LEiWWmwi8nIyMDm82mBCziARoJi4iImEQjYREREZMoCYuIiJhESVhERMQkSsIiIiImURIWERExiZKwiIiISf4fGz0toDuSWeYAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {
"tags": []
}
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "Jkur6P4hgaJ2",
"colab_type": "text"
},
"source": [
"### Normal score transform"
]
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:33.515489Z",
"start_time": "2020-06-03T02:32:33.400546Z"
},
"id": "tOVXuHEugaJ3",
"colab_type": "code",
"colab": {}
},
"source": [
"unscore = gs.Program(exedir+'unscore')"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:33.700496Z",
"start_time": "2020-06-03T02:32:33.525504Z"
},
"id": "QoOK00UCgaKK",
"colab_type": "code",
"colab": {},
"outputId": "be4d9b77-a12c-46f6-846b-772809d60996"
},
"source": [
"parstr = \"\"\" Parameters for NSCORE\n",
" *********************\n",
"\n",
"START OF PARAMETERS:\n",
"{datafl} -file with data\n",
"{nvar} {varcols} - number of variables and columns\n",
"{wtcol} - column for weight, 0 if none\n",
"0 - column for category, 0 if none\n",
"0 - number of records if known, 0 if unknown\n",
"-1.0e21 1.0e21 - trimming limits\n",
"0 -transform using a reference distribution, 1=yes\n",
"../histsmth/histsmth.out -file with reference distribution.\n",
"1 2 0 -columns for variable, weight, and category\n",
"1000 -maximum number of quantiles, 0 for all\n",
"{outfl} -file for output\n",
"{trnfl} -file for output transformation table\n",
"\"\"\"\n",
"pars = dict(datafl=dat.flname, nvar=dat.nvar,\n",
" varcols=dat.gscol(dat.variables),\n",
" wtcol=dat.gscol(dat.weights),\n",
" outfl = os.path.join(outdir, 'nscore.out'), trnfl = os.path.join(outdir,'nscore.trn'))\n",
"unscore.run(parstr=parstr.format(**pars))"
],
"execution_count": null,
"outputs": [
{
"output_type": "stream",
"text": [
"Calling: ['../pygeostat/executable/unscore', 'temp']\n",
" UNSCORE Version: 1.1.1\n",
" data file = Output/declus.out\n",
" columns = 4 5 7 0\n",
" number of records = 0\n",
" trimming limits = -1.000000000000000E+021 1.000000000000000E+021\n",
" consider a different reference dist = 0\n",
" file with reference distribution = ../histsmth/histsmth.out \n",
" columns = 1 2 0\n",
" maximum number of quantiles = 1000\n",
" file for output = Output/nscore.out\n",
" file for transformation table = Output/nscore.trn\n",
" Determining size of Output/declus.out\n",
" Reading Output/declus.out\n",
" Building transform table 1 category 1\n",
" Building transform table 2 category 1\n",
" Computing normal scores 1 category 1\n",
" Computing normal scores 2 category 1\n",
" Generating output file\n",
" Total execution time 0.100 seconds.\n",
"\n",
" UNSCORE Version: 1.1.1 Finished\n",
"\n"
],
"name": "stdout"
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "CQ2QVQ1vgaLU",
"colab_type": "text"
},
"source": [
"### Load and inspect the normal score transformation result\n",
"\n",
"Use the notvariables kwarg leads to isolation of the normal score variables as the dat_ns.variables attribute. Since DataFiles with iniitialized wts are passed to histplt, wt may be True."
]
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:33.759525Z",
"start_time": "2020-06-03T02:32:33.704493Z"
},
"scrolled": false,
"id": "JP-Vf6E1gaLW",
"colab_type": "code",
"colab": {},
"outputId": "be68f8ed-0347-4bb0-a028-974b2961b9de"
},
"source": [
"dat_ns = gs.DataFile(os.path.join(outdir, 'nscore.out'), \n",
" notvariables=['Base Elevation']+dat.variables)\n",
"print('variables = ', dat_ns.variables)"
],
"execution_count": null,
"outputs": [
{
"output_type": "stream",
"text": [
"variables = ['NS_Top Elevation', 'NS_Thickness']\n"
],
"name": "stdout"
}
]
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:33.805490Z",
"start_time": "2020-06-03T02:32:33.766491Z"
},
"id": "5_OzhxMngaLc",
"colab_type": "code",
"colab": {},
"outputId": "feebbf1a-7484-4c61-f273-c5110c7a3f43"
},
"source": [
"dat_ns.data"
],
"execution_count": null,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" HoleID | \n",
" X | \n",
" Y | \n",
" Top Elevation | \n",
" Thickness | \n",
" Base Elevation | \n",
" Declustering Weight | \n",
" NS_Top Elevation | \n",
" NS_Thickness | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" 3.0 | \n",
" 405.63 | \n",
" 2135.75 | \n",
" 376.69 | \n",
" 47.98 | \n",
" 328.71 | \n",
" 1.314446 | \n",
" -0.97569 | \n",
" -0.42111 | \n",
"
\n",
" \n",
" | 1 | \n",
" 5.0 | \n",
" 235.89 | \n",
" 1865.70 | \n",
" 379.69 | \n",
" 51.00 | \n",
" 328.69 | \n",
" 1.973776 | \n",
" 0.22327 | \n",
" 0.46519 | \n",
"
\n",
" \n",
" | 2 | \n",
" 6.0 | \n",
" 325.03 | \n",
" 2055.81 | \n",
" 376.86 | \n",
" 49.34 | \n",
" 327.52 | \n",
" 1.025138 | \n",
" -0.85014 | \n",
" 0.04552 | \n",
"
\n",
" \n",
" | 3 | \n",
" 7.0 | \n",
" 675.54 | \n",
" 2195.25 | \n",
" 381.49 | \n",
" 48.75 | \n",
" 332.74 | \n",
" 0.912454 | \n",
" 0.85265 | \n",
" -0.21400 | \n",
"
\n",
" \n",
" | 4 | \n",
" 8.0 | \n",
" 355.73 | \n",
" 1995.74 | \n",
" 376.97 | \n",
" 48.94 | \n",
" 328.03 | \n",
" 1.758901 | \n",
" -0.79922 | \n",
" -0.14376 | \n",
"
\n",
" \n",
" | ... | \n",
" ... | \n",
" ... | \n",
" ... | \n",
" ... | \n",
" ... | \n",
" ... | \n",
" ... | \n",
" ... | \n",
" ... | \n",
"
\n",
" \n",
" | 225 | \n",
" 296.0 | \n",
" 245.18 | \n",
" 1655.43 | \n",
" 378.33 | \n",
" 49.45 | \n",
" 328.88 | \n",
" 1.973776 | \n",
" -0.22713 | \n",
" 0.08391 | \n",
"
\n",
" \n",
" | 226 | \n",
" 297.0 | \n",
" 715.53 | \n",
" 1845.27 | \n",
" 380.60 | \n",
" 47.48 | \n",
" 333.12 | \n",
" 0.516043 | \n",
" 0.58374 | \n",
" -0.58766 | \n",
"
\n",
" \n",
" | 227 | \n",
" 298.0 | \n",
" 685.11 | \n",
" 1605.33 | \n",
" 374.90 | \n",
" 50.95 | \n",
" 323.95 | \n",
" 0.548732 | \n",
" -1.79667 | \n",
" 0.43325 | \n",
"
\n",
" \n",
" | 228 | \n",
" 300.0 | \n",
" 835.21 | \n",
" 1455.71 | \n",
" 376.89 | \n",
" 49.25 | \n",
" 327.64 | \n",
" 1.093967 | \n",
" -0.83368 | \n",
" -0.03194 | \n",
"
\n",
" \n",
" | 229 | \n",
" 301.0 | \n",
" 125.06 | \n",
" 1455.92 | \n",
" 379.69 | \n",
" 52.50 | \n",
" 327.19 | \n",
" 1.098478 | \n",
" 0.20613 | \n",
" 0.84717 | \n",
"
\n",
" \n",
"
\n",
"
230 rows × 9 columns
\n",
"
"
]
},
"metadata": {
"tags": []
}
}
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "b2DTREdIgaLt",
"colab_type": "text"
},
"source": [
"## 4. Variogram Calculation and Modeling\n",
"Normal score variograms are calculated and modeled, before being input to sequential Gaussian simulation in the next section.\n"
]
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:36.333771Z",
"start_time": "2020-06-03T02:32:36.318769Z"
},
"id": "CQlmzg4SgaLw",
"colab_type": "code",
"colab": {}
},
"source": [
"var_calc = gs.Program(program=exedir+'varcalc')"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:36.448769Z",
"start_time": "2020-06-03T02:32:36.336767Z"
},
"id": "5f_Co0UTgaME",
"colab_type": "code",
"colab": {},
"outputId": "0990a4e9-e3d0-41ea-f09c-2ab0a909c833"
},
"source": [
"parstr = \"\"\" Parameters for VARCALC\n",
" **********************\n",
" \n",
"START OF PARAMETERS:\n",
"{file} -file with data\n",
"2 3 0 - columns for X, Y, Z coordinates\n",
"2 8 9 - number of variables,column numbers (position used for tail,head variables below)\n",
"{t_min} 1.0e21 - trimming limits\n",
"{n_directions} -number of directions\n",
"0.0 90 10000 0.0 22.5 1000 0.0 -Dir 01: azm,azmtol,bandhorz,dip,diptol,bandvert,tilt\n",
" {n_lag} {lag_length} {lag_tol} - number of lags,lag distance,lag tolerance\n",
"{output} -file for experimental variogram points output.\n",
"0 -legacy output (0=no, 1=write out gamv2004 format)\n",
"1 -run checks for common errors\n",
"1 -standardize sills? (0=no, 1=yes)\n",
"2 -number of variogram types\n",
"1 1 1 1 -tail variable, head variable, variogram type (and cutoff/category), sill\n",
"2 2 1 1 -tail variable, head variable, variogram type (and cutoff/category), sill\n",
"\"\"\"\n",
"\n",
"n_directions = 1\n",
"varcalc_outfl = os.path.join(outdir, 'varcalc.out')\n",
"\n",
"var_calc.run(parstr=parstr.format(file=dat_ns.flname,\n",
" n_directions = 1,\n",
" t_min = gs.Parameters['data.tmin'],\n",
" n_lag=10,\n",
" lag_length = 30.0,\n",
" lag_tol = 20.0,\n",
" output=varcalc_outfl),\n",
" liveoutput=True)"
],
"execution_count": null,
"outputs": [
{
"output_type": "stream",
"text": [
"Calling: ['../pygeostat/executable/varcalc', 'temp']\n",
"\n",
"varcalc version: 1.400\n",
"\n",
" data file: Output/nscore.out\n",
" x,y,z columns: 2 3 0\n",
" number of variables: 2\n",
" Variable columns: 8 9\n",
" tmin,tmax: -998.000000000000 1.000000000000000E+021\n",
" number of directions: 1\n",
" direction parameters:\n",
" azm,azmtol,bandhorz 0.000000000000000E+000 90.0000000000000 \n",
" 10000.0000000000 \n",
" dip,diptol,bandvert 0.000000000000000E+000 22.5000000000000 \n",
" 1000.00000000000 \n",
" tilt 0.000000000000000E+000\n",
" nlags,lagdist,lagtol 10 30.0000000000000 \n",
" 20.0000000000000 \n",
" output file: Output/varcalc.out\n",
" legacy output? 0\n",
" run checks? 1\n",
" attempt to standardize sills? 1\n",
" number of variogram types: 2\n",
" Variogram tail,head,type 1 1 1\n",
" standardizing with sill = 1.00000000000000 \n",
" Variogram tail,head,type 2 2 1\n",
" standardizing with sill = 1.00000000000000 \n",
" Reading data file\n",
" Setting up final parameters for variogram calculation\n",
" Calculating variograms...\n",
" working on direction 1\n"
],
"name": "stdout"
}
]
},
{
"cell_type": "code",
"metadata": {
"ExecuteTime": {
"end_time": "2020-06-03T02:32:36.666769Z",
"start_time": "2020-06-03T02:32:36.452767Z"
},
"id": "v8CsteF7gaMg",
"colab_type": "code",
"colab": {},
"outputId": "cbf914a2-66ba-4e2d-e947-7cacb6ab5a4e"
},
"source": [
"varfl_exp = gs.DataFile(varcalc_outfl)\n",
"varfl_exp.head()"
],
"execution_count": null,
"outputs": [
{
"output_type": "execute_result",
"data": {
"text/html": [
"\n",
"\n",
"
\n",
" \n",
" \n",
" | \n",
" Variogram Index | \n",
" Lag Distance | \n",
" Number of Pairs | \n",
" Variogram Value | \n",
" Variogram Number | \n",
" Calculation Azimuth | \n",
" Calculation Dip | \n",
" Variogram Type | \n",
" Variogram Tail Index | \n",
" Variogram Head Index | \n",
"
\n",
" \n",
" \n",
" \n",
" | 0 | \n",
" 1.0 | \n",
" 13.055639 | \n",
" 44.0 | \n",
" 0.013837 | \n",
" 1.0 | \n",
" 0.0 | \n",
" 0.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
"
\n",
" \n",
" | 1 | \n",
" 1.0 | \n",
" 34.531689 | \n",
" 407.0 | \n",
" 0.100780 | \n",
" 1.0 | \n",
" 0.0 | \n",
" 0.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
"
\n",
" \n",
" | 2 | \n",
" 1.0 | \n",
" 61.017428 | \n",
" 790.0 | \n",
" 0.262916 | \n",
" 1.0 | \n",
" 0.0 | \n",
" 0.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
"
\n",
" \n",
" | 3 | \n",
" 1.0 | \n",
" 90.311682 | \n",
" 1250.0 | \n",
" 0.447996 | \n",
" 1.0 | \n",
" 0.0 | \n",
" 0.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
"
\n",
" \n",
" | 4 | \n",
" 1.0 | \n",
" 120.364562 | \n",
" 1541.0 | \n",
" 0.610774 | \n",
" 1.0 | \n",
" 0.0 | \n",
" 0.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
" 1.0 | \n",
"
\n",
" \n",
"
\n",
"