{ "cells": [ { "cell_type": "markdown", "id": "83ff2437", "metadata": {}, "source": [ "

\n", " \n", "\n", "

\n", "\n", "## Spatial Data Analytics \n", "\n", "### Interactive Demonstration of Variogram h-Scatterplots \n", "\n", "#### Michael Pyrcz, Professor, The University of Texas at Austin \n", "\n", "##### Contacts: [Twitter/@GeostatsGuy](https://twitter.com/geostatsguy) | [GitHub/GeostatsGuy](https://github.com/GeostatsGuy) | [www.michaelpyrcz.com](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446)\n", "\n", "This a simple demonstration of the variogram h-scatterplot for a 1D datasets with variable spatial continuity and visualization.\n", "\n", "* we will see the correlogram (equal to the covariance function when the sill, variance is 1.0) is the correlation coefficient of the h-scatterplot. \n", "\n", "* there is some deviation due to the lag effect, the edge effect with variogram calculation that excludes some of the data (e.g., at large lags only the samples at the edges of the area of interest are included in the pairs)\n", "\n", "* we will perform the calculations in 1D for fast run times and ease of visualization.\n", "\n", "Please cite this code as:\n", "\n", "Pyrcz, Michael J. (2021). PythonNumericalDemos: Educational Data Science Demonstrations Repository (1.0.0). Zenodo. https://doi.org/10.5281/zenodo.5564991\n", "\n", "#### Load the required libraries\n", "\n", "The following code loads the required libraries." ] }, { "cell_type": "code", "execution_count": 1, "id": "812f4f78", "metadata": {}, "outputs": [], "source": [ "supress_warnings = True\n", "import numpy as np # arrays\n", "import matplotlib.pyplot as plt # plotting\n", "from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) # control of axes ticks\n", "import pandas as pd # dataframes\n", "from geostatspy import GSLIB # affine correction\n", "from ipywidgets import interactive # plotting widgets and interactivity\n", "from ipywidgets import widgets \n", "from ipywidgets import Layout\n", "from ipywidgets import Label\n", "from ipywidgets import VBox, HBox\n", "cmap = plt.cm.inferno # default color bar, no bias and friendly for color vision defeciency\n", "plt.rc('axes', axisbelow=True) # grid behind plotting elements\n", "if supress_warnings == True:\n", " import warnings # supress any warnings for this demonstration\n", " warnings.filterwarnings('ignore') " ] }, { "cell_type": "markdown", "id": "f6c226de", "metadata": {}, "source": [ "#### Declar Functions" ] }, { "cell_type": "code", "execution_count": 2, "id": "f1a3abb0", "metadata": {}, "outputs": [], "source": [ "def add_grid():\n", " plt.gca().grid(True, which='major',linewidth = 1.0); plt.gca().grid(True, which='minor',linewidth = 0.2) # add y grids\n", " plt.gca().tick_params(which='major',length=7); plt.gca().tick_params(which='minor', length=4)\n", " plt.gca().xaxis.set_minor_locator(AutoMinorLocator()); plt.gca().yaxis.set_minor_locator(AutoMinorLocator()) # turn on minor ticks " ] }, { "cell_type": "markdown", "id": "c46a4b7c", "metadata": {}, "source": [ "#### Set the working directory\n", "\n", "I always like to do this so I don't lose files and to simplify subsequent read and writes (avoid including the full address each time). Also, in this case make sure to place the required (see below) data file in this working directory. " ] }, { "cell_type": "code", "execution_count": 3, "id": "90b55ab1", "metadata": {}, "outputs": [], "source": [ "#os.chdir(\"C:\\PGE337\") # set the working directory" ] }, { "cell_type": "markdown", "id": "63fed3ff", "metadata": {}, "source": [ "#### Load the Dataset\n", "\n", "It is a small 1D dataset available on my GitHup [GeoDataSets](https://github.com/GeostatsGuy/GeoDataSets) repository.\n", "\n", "Datasets are cited as: \n", "\n", "* Pyrcz, Michael J. (2021). GeoDataSets: Synthetic Subsurface Data Repository (1.0.0). Zenodo. https://doi.org/10.5281/zenodo.5564874" ] }, { "cell_type": "code", "execution_count": 4, "id": "e4439445", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
DepthNporosity
00.25-1.37
10.50-2.08
20.75-1.67
31.00-1.16
41.25-0.24
\n", "
" ], "text/plain": [ " Depth Nporosity\n", "0 0.25 -1.37\n", "1 0.50 -2.08\n", "2 0.75 -1.67\n", "3 1.00 -1.16\n", "4 1.25 -0.24" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "df = pd.read_csv('https://raw.githubusercontent.com/GeostatsGuy/GeoDataSets/master/1D_Porosity.csv')\n", "npor = df['Nporosity']\n", "npor = GSLIB.affine(npor,0.0,1.0) # ensure variance is 1.0 for results to work below\n", "df.head()" ] }, { "cell_type": "code", "execution_count": 5, "id": "38c10646", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAHKCAYAAADsGyoAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACXcklEQVR4nO3deXwTdfoH8M8kve+D3i1tAeWQQw4VKuVyLSKwxW5VPEFR17UqCKyKugqK4gEIrvX4rQseux5rjbiorHhQrAKKUIRyCdjSuwV6pG16Jt/fH9NJmzZpMpNrkjzv14tXSzKT+U4nnT75Hs/DMcYYCCGEEEI8hMLZDSCEEEIIcSQKfgghhBDiUSj4IYQQQohHoeCHEEIIIR6Fgh9CCCGEeBQKfgghhBDiUSj4IYQQQohHoeCHEEIIIR6Fgh9CCCGEeBQKfojbefvtt8FxnP6fl5cXEhMTcccdd6CiosLZzTNKaHNJSYn+sffffx+bNm2y+bHq6uqwcOFCREdHg+M4LFiwwObHsKXjx4/jtttuw5AhQ+Dn54dBgwZhwoQJuP/++6FWq/XbLV68GCkpKXZvT0pKChYvXmz163R2diImJgaTJ082uY1Op8PgwYMxduxYq48HAHv27MHq1avR0NDQ77kZM2ZgxowZNjmOMZZcn/vvvx8cx6G6utrg8bq6OigUCnh7e6O5udngufLycnAch+XLl4tqT9/rWFJSAo7j8Pbbb4t6HeKaKPghbmvr1q3Yu3cvvv76a9x999344IMPkJ6ejpaWFmc3rZ+5c+di7969iIuL0z9mr+DnmWeewaeffoqXX34Ze/fuxYsvvmjzY9hKYWEhJk6ciGPHjuHJJ5/E//73P7zxxhuYO3cuvvrqK9TV1em3/dvf/oZPP/3Uia0Vx9vbG7fddht++uknHDt2zOg233zzDcrKyrBkyRKbHHPPnj1Ys2aN0eDntddew2uvvWaT40g1c+ZMAEB+fr7B47t374aXlxc4jsMPP/xg8NyuXbsM9iXEEl7ObgAh9jJ69GhMmjQJAH9j1Gq1eOaZZ7Bt2zbccsstVr12a2sr/P39bdFMAEBUVBSioqJs9noDKSoqwtChQ63+GQgYY2hra7Ppz0OwadMmKBQK5OfnIzg4WP94dnY2nnnmGfQuTTh06FCbH9/elixZgg0bNmDLli1Yv359v+e3bNkCHx8f3HrrrVYdp7W1FX5+fgNuM2rUKKuOYQszZswAx3HIz8/HwoUL9Y/n5+fjsssuA2MMu3btwjXXXGPwnEKhwLRp05zRZOKiqOeHeAxheOHs2bMAgLa2NqxatQqpqanw8fFBQkICcnJy+n0qTklJwbx586BSqTB+/Hj4+flhzZo1APhAIjMzE+Hh4fDz88Oll16Kd955x2B/nU6HtWvXYvjw4fD390dYWBjGjh2LzZs367fpO+w1Y8YMfPHFFzh79qzBEB5jDBdddBFmz57d7/yam5sRGhqKnJwco+cvdOt/8803OH78uP41hU/ZdXV1uO+++5CQkAAfHx8MGTIEjz/+ONrb2w1eh+M43H///XjjjTcwcuRI+Pr69jvn3j766CNkZGQgLi4O/v7+GDlyJB599FGLeuAuXLiAkJAQBAUFGX2e4zj998aGVYS2vvfeexg5ciQCAgIwbtw4fP755/1e67PPPsPYsWPh6+uLIUOGYPPmzVi9erXBMUxRq9VYuXKlwXtp2bJlZs9x5MiRmDJlCt577z10dXUZPNfQ0IDPPvsMmZmZiIyMBAD88ssv+OMf/4iIiAj4+flh/Pjx+M9//mOwn/Be2rlzJ+68805ERUUhICAAq1atwl//+lcAQGpqar/rb2zYq729HU8//TRGjhwJPz8/REZGYubMmdizZ49+m9zcXEybNg3R0dEIDAzEmDFj8OKLL6Kzs9Psz62vyMhIjBkzpl/PT35+PmbMmIHp06fre3p6PzdhwgSEhoYCkH4tiGehnh/iMU6fPg2A72VhjGHBggX49ttvsWrVKqSnp+Pw4cN46qmnsHfvXuzduxe+vr76fQ8ePIjjx4/jiSeeQGpqKgIDA3Hy5EmkpaUhOjoar7zyCiIjI/Gvf/0LixcvRk1NDR5++GEAwIsvvojVq1fjiSeewLRp09DZ2YkTJ04YHXoQvPbaa7jnnntw5swZg6EcjuPwwAMPYNmyZTh16hQuuugi/XPvvvsu1Gq1yeAnLi4Oe/fuxX333YfGxkb8+9//BsB/4m9ra8PMmTNx5swZrFmzBmPHjkVBQQHWrVuHQ4cO4YsvvjB4rW3btqGgoABPPvkkYmNjER0dbfJcTp06hWuvvRbLli1DYGAgTpw4gRdeeAE///wzvvvuO5P7AcCUKVPwxRdf4JZbbsGf//xnXH755aJ7mL744gvs378fTz/9NIKCgvDiiy/iuuuuw8mTJzFkyBAAwP/+9z9kZWVh2rRp+Oijj9DV1YX169ejpqbG7OtrNBpMnz4d5eXleOyxxzB27FgcPXoUTz75JI4cOYJvvvlmwABqyZIluOuuu/DFF18gMzNT//j777+PtrY2/ZCX0ONxxRVX4I033kBoaCg+/PBD3HjjjdBoNP3mId15552YO3cu3nvvPbS0tGDSpEnQaDT4+9//DpVKpR9iNdXj09XVhTlz5qCgoADLli3DrFmz0NXVhX379qG0tBRpaWkAgDNnzuDmm2/WBxu//vornn32WZw4cQJbtmwx+/Pra+bMmdi8eTOqqqoQFxeHCxcu4MiRI3jppZeg0+nw0ksvQa1WIyQkBGVlZfj999/xpz/9ySbXgngQRoib2bp1KwPA9u3bxzo7O1lTUxP7/PPPWVRUFAsODmbV1dXsf//7HwPAXnzxRYN9P/roIwaA/d///Z/+seTkZKZUKtnJkycNtl24cCHz9fVlpaWlBo/PmTOHBQQEsIaGBsYYY/PmzWOXXnqpRW0uLi7WPzZ37lyWnJzcb1u1Ws2Cg4PZ0qVLDR4fNWoUmzlz5oDHYYyx6dOns0suucTgsTfeeIMBYP/5z38MHn/hhRcYALZz5079YwBYaGgoq6urM3usvnQ6Hevs7GS7d+9mANivv/464PZtbW1swYIFDAADwJRKJRs/fjx7/PHHWW1trcG2ixYt6vfzAsBiYmKYWq3WP1ZdXc0UCgVbt26d/rHLLruMJSUlsfb2dv1jTU1NLDIykvW9TSYnJ7NFixbp/79u3TqmUCjY/v37DbbLy8tjANiXX3454Dk2NTWxoKAg9sc//tHg8YkTJ7KkpCSm1WoZY4yNGDGCjR8/nnV2dhpsN2/ePBYXF6ffTngv3X777f2O9dJLL/V7nwmmT5/Opk+frv//u+++ywCwf/zjHwO2vzetVss6OzvZu+++y5RKpcF7xNj1MWbbtm0MAHv//fcZY4x98sknzMvLizU1NTG1Ws2USiX7/PPPGWOMvfPOOwY/YzHXou91LC4uZgDY1q1bLT5f4rpo2Iu4rcmTJ8Pb2xvBwcGYN28eYmNjsWPHDsTExOh7HPp+Wr7++usRGBiIb7/91uDxsWPH4uKLLzZ47LvvvsNVV12FpKQkg8cXL14MjUaDvXv3AgAuv/xy/Prrr7jvvvvw1VdfGaxQkiI4OBh33HEH3n77bX1X/nfffYdjx47h/vvvl/Sa3333HQIDA5Gdnd3vXAD0+3nMmjUL4eHhFr3277//jptvvhmxsbFQKpXw9vbG9OnTAfAruQbi6+uLTz/9FMeOHcPLL7+MhQsX4ty5c3j22WcxcuRInDx50uzxZ86caTBfKCYmBtHR0frhz5aWFvzyyy9YsGABfHx89NsFBQVh/vz5Zl//888/x+jRo3HppZeiq6tL/2/27NkGw0qmBAUF4YYbbsCXX36p72kqKirCgQMHsHjxYigUCpw+fRonTpzQz9PqfZxrr70WVVVV/X4WQm+IVDt27ICfnx/uvPPOAbcrLCzEH//4R0RGRuqv7+233w6tVovffvtN9HGnT5+un+cF8MNakyZNQlBQEIKDgzFhwgT90Fd+fj68vLwwdepUANZfC+I5KPghbuvdd9/F/v37UVhYiMrKShw+fBhXXnklAH4uiZeXV79JxhzHITY2FhcuXDB4vPcqLMGFCxeMPh4fH69/HgBWrVqF9evXY9++fZgzZw4iIyNx1VVX4ZdffpF8bg888ACampr0Q1evvvoqEhMTDYZNxLhw4QJiY2P7DQlER0fDy8vLop+HMc3NzUhPT8dPP/2EtWvXIj8/H/v374dKpQLAT8S1xMiRI7Fs2TL861//QmlpKTZu3IgLFy7gb3/7m9l9hfkyvfn6+uqPXV9fD8YYYmJi+m1n7LG+ampqcPjwYXh7exv8Cw4OBmMM58+fN/saS5YsQVdXF9577z0A/ERnjuNwxx136I8BACtXrux3nPvuuw8A+h3H0mtkyrlz5xAfHw+FwvSfidLSUqSnp6OiogKbN29GQUEB9u/fj9zcXACWX9/ewsLCcOmll+oDnF27dumDZYAPjoQgZteuXZg0aZI+uLXFtSCegeb8ELc1cuRI/WqvviIjI9HV1YVz584ZBECMMVRXV+Oyyy4z2N7YPIHIyEhUVVX1e7yyshIAMGjQIACAl5cXli9fjuXLl6OhoQHffPMNHnvsMcyePRtlZWUICAgQfW7Dhg3DnDlzkJubizlz5uC///0v1qxZA6VSKfq1hHP56aefwBgzONfa2lp0dXXpz0Vg6byJ7777DpWVlcjPzzf4AzbQfCdzOI7DQw89hKeffhpFRUWSX0cQHh4OjuOMzu/pm2/GmEGDBsHf39/k/Ja+Pztj0tLSMHLkSGzduhVLly7Fv/71L8yaNQupqakGr7Fq1SpkZWUZfY3hw4cb/N/auS1RUVH44YcfoNPpTAZA27ZtQ0tLC1QqFZKTk/WPHzp0yKpjz5w5Exs2bMDhw4dx9OhRg3QM06dPx8aNG3H48GGUlJTgpptu0j9ni2tBPAP1/BCPdNVVVwEA/vWvfxk8/sknn6ClpUX/vLnXEP649/buu+8iICDAaPK6sLAwZGdnIycnB3V1dQZJDfvq3TthzNKlS3H48GEsWrQISqUSd999t9k2D3Quzc3N2LZtm8Hj7777rv55KYQ/wL0njwPAm2++adH+xoJLgA8w1Wq1vpfNGoGBgZg0aRK2bduGjo4O/ePNzc1GV4X1NW/ePJw5cwaRkZGYNGlSv3+WJl688847cezYMTzxxBM4d+6cwXDT8OHDcdFFF+HXX381eozevR8DEa6DJT0yc+bMQVtb24BJ/4xdX8YY/vGPf5h9/YEIOXvWrFkDhUKhH9YCoP9eWHHZO7+Pra4FcX/U80M80tVXX43Zs2fjkUcegVqtxpVXXqlf7TV+/HjcdtttZl/jqaeewueff46ZM2fiySefREREBP7973/jiy++wIsvvqhfejt//nx9zqGoqCicPXsWmzZtQnJyssFqrb7GjBkDlUqF119/HRMnToRCoTDoybr66qsxatQo7Nq1C7feeuuAK67Muf3225Gbm4tFixahpKQEY8aMwQ8//IDnnnsO1157Lf7whz9Iet20tDSEh4fj3nvvxVNPPQVvb2/8+9//xq+//mrR/vfccw8aGhrwpz/9CaNHj4ZSqcSJEyfw8ssvQ6FQ4JFHHpHUrr6efvppzJ07F7Nnz8bSpUuh1Wrx0ksvISgoyCCRojHLli3DJ598gmnTpuGhhx7C2LFjodPpUFpaip07d2LFihW44oorzLbh9ttvx2OPPYaXXnoJYWFh/Xp43nzzTcyZMwezZ8/G4sWLkZCQgLq6Ohw/fhwHDx7Exx9/bPYYY8aMAQBs3rwZixYtgre3N4YPH240cLrpppuwdetW3HvvvTh58iRmzpwJnU6Hn376CSNHjsTChQtx9dVXw8fHBzfddBMefvhhtLW14fXXX0d9fb3Ztgxk2rRpUCqV+PTTT/sFdmFhYRg3bhw+/fRTeHt764eyAdtdC+IBnDnbmhB7EFa79F3x0Vdrayt75JFHWHJyMvP29mZxcXHsL3/5C6uvrzfYLjk5mc2dO9foaxw5coTNnz+fhYaGMh8fHzZu3Lh+q0U2bNjA0tLS2KBBg5iPjw8bPHgwW7JkCSspKenX5t6rcOrq6lh2djYLCwtjHMf1W3XEGGOrV6/Wr2yzlLHVXowxduHCBXbvvfeyuLg45uXlxZKTk9mqVatYW1ubwXYAWE5OjsXH27NnD5syZQoLCAhgUVFR7K677mIHDx60aGXNV199xe688042atQoFhoayry8vFhcXBzLyspie/fuNdjW1GovY23tu9KHMcY+/fRTNmbMGP01ev7559mDDz7IwsPDze7b3NzMnnjiCTZ8+HDm4+PDQkND2ZgxY9hDDz3EqqurBzzH3q677joGgN13331Gn//111/ZDTfcwKKjo5m3tzeLjY1ls2bNYm+88YZ+G3Pv/1WrVrH4+HimUCgYALZr1y7GWP/VXozxvyNPPvkku+iii5iPjw+LjIxks2bNYnv27NFvs337djZu3Djm5+fHEhIS2F//+le2Y8cOg9dmzPLVXoLLL7+cAWArV67s99yyZcsYAHbllVf2e87Sa0GrvTwbx1ivFKmEEJcyadIkcByH/fv3O7spbqezsxOXXnopEhISsHPnTmc3hxBiQzTsRYiLUavVKCoqwueff44DBw64VD0rOVuyZAmuvvpqxMXFobq6Gm+88QaOHz9ukImbEOIeKPghxMUcPHgQM2fORGRkJJ566inZV2V3FU1NTVi5ciXOnTsHb29vTJgwAV9++aXk+U6EEPmiYS9CCCGEeBRa6k4IIYQQj0LBDyGEEEI8CgU/hBBCCPEoNOG5D51Oh8rKSgQHB1udHp4QQgghjsEYQ1NTk9madAAFP/1UVlb2q9JNCCGEENdQVlaGxMTEAbeh4KcPIY16UVGRxUFQc3MzACAoKMji44jdxxHHAPiCkz/++COuvPJKhIWFyaZdjthH7Lk7ql1yvfZSjkPXXp7vF0++9oD73Pfkeu2l7CPlGFVVVRgxYoRFde4o+OlDGOoKDg5GSEiIRfsI3WtiLpLYfRxxDIAf9gsICEBISIhF5++odjliH7Hn7qh2yfXaSzkOXXt5vl88+doD7nPfk+u1l7KPlGMIAZMlU1ZowjMhhBBCPAr1/BBCCCEuTqvVoqCgANXV1UhNTUV6ejqUSqWzmyVbFPwQQgghLkylUmHF0qUoKS/XP5aSmIgNmzcjKyvLiS2TLxr2IoQQQlyUSqVCdnY2xpSXYy+AJgB7AYypqEB2djZUKpWTWyhP1PNjQltbm37ylDkajUb064vdxxHHAIDW1lb9Vy8v828PR7XLEfuIPXdHtUuu117Kcejay/P94snXHnDd+55Wq8XyBx/EPMawDT29GZMBbGMMCzgOKx58EFdddZXJITC5vvetuY6W8OjgJzc3F7m5uQaPabVaJ7WGEEIIsdyePXtwtqICH6L/MI4CwCrGkFZRgT179iA9Pd0JLZQvjw5+cnJykJOTY/CYWq1GaGgo/Pz8RC2xA8QtyZO6j72P0dXVBQDw9/cXtZ8jzt3e+0g9d3u3y1HHcOT507WX1/vFk6894Lr3vYaGBgDAaBP7je61nanXkPt7X8z2arXa4m1pzg8hhBDiguLi4gAARSaeL+qzHelBwQ8hhBDigtLT05GSmIjnOA66Ps/pAKzjOKQmJdGQlxEU/BBCCCEuSKlUYsPmzfgcwAKOM1jttYDj8DmA9Zs2Ub4fIyj4IYQQQlxUVlYW8vLycCQhAWkAQgCkAShKTEReXh7l+THBoyc8E0IIIa4uKysLmRkZ+HrLFlRfuICU+Hik33EHlD4+zm6abFHwQwghhLg4pVaLaePHAwACAgIAXd9ZQKQ3GvYihBBCXF17u+H/Ozud0w4XQcEPIYQQ4urCw6EbOhTgOP7/FPwMiIIfQgghxNV5ewMhIWBhYfz/KfgZEAU/hBBCiLvw9ua/UvAzIJrwTAghhLi6mhqgrQ0sIgJISgJopdeAKPgxgaq6u1Z1Y1vs407Vjamyt3tceyn70LX3wPueTgfFb7+ho60NXZdcwq/0amuzaFe5vvepqrsdUVV3QgghLk8Y4lIoAAsDGE/n0T8lquren6tWN7bFPu5U3VjKPp5c2Vvu117KPnTtPei+xxgQEAD4+UHn54eg5magqwtITDT72nJ/71NVd0IIIYT0153jhwmTnSsq+DlANJJhEgU/hBBCiCvr6OC/+vjwQ19CIVNa8WUSBT+EEEKIK+sd/PT+KjxO+qHghxBCCHFl3UEOE4IeyvVjlkdPeCaEEEJcXnIyP+9HmONDwY9Z1PNDCCGEuDJfXyAkpCfooeDHLAp+CCGEEHdCwY9ZNOxFCCGEuKr2dqCuDvD370lwGB4OBAVRiYsBUM8PIYQQ4qpaWoDKSj6vj8Dbm096SNmeTXLp4GfdunW47LLLEBwcjOjoaCxYsAAnT550drMIIYQQx+i7zJ1YxKWDn927dyMnJwf79u3D119/ja6uLmRkZKClpcXZTSOEEELsTwh+fH0NH6+pAcrKKMuzCS7dJ/a///3P4P9bt25FdHQ0Dhw4gGnTpjmpVYQQQoiDmOr5qa7m63sNGsTPByIGXDr46auxsREAEBERYfT59vZ2tHfXQDFFKIzW3NyMhoYGi46r0WgA9BSIs8c+jjgGADQ1NRl8lUu7HLGP2HN3VLvkeu2lHIeuvTzfL5587QHXvu8pzp8H19YGrUYDTXcg1NXVBUVrK//4uXP8MngT5Pret+Y6WoJjjDGLt5YxxhgyMzNRX1+PgoICo9usXr0aa9assej13n//fQQEBNiyiYQQQohNhfz+OzidDk2DB0PXq/cnsKoKXi0taI2ORscAwY870Wg0uPnmm9HY2IgQM+fsNsFPTk4OvvjiC/zwww9ITEw0uo2lPT9JSUk4evQo4uPjLTq2EKGKCZbE7uOIYwB85Hzw4EFMmDABwcHBsmmXI/YRe+6Oapdcr72U49C1l+f7xZOvPeDC9z2tFsojR/hvx46Fpq1Nvz1XWgpFXR10sbFgsbEmjyHX976UY9TU1GDEiBEWBT9uMez1wAMP4L///S++//57k4EPAPj6+sK376QwE4KCghAWFmbRtl7dywmDgoIs2l7KPo44Rm/BwcEWnb+j2uXI87f03B3VLrleeynHoWsvz/eLwJOvPeCC9z3GgMmT+Xk/ISHwam7u2V6j4ZMcBgQAFpyT3N77Uo4hBEwWvb7FW8oQYwwPPPAAPv30U+Tn5yM1NdXZTSKEEEIcg+MAPz/+X1+U5XlALh385OTk4P3338dnn32G4OBgVFdXAwBCQ0PhT7PbCSGEeCoKfgbk0sHP66+/DgCYMWOGweNbt27F4sWLHd8gQgghxFHq64G2NiA0lB/e6i04GBg5kpIfmuDSwY+bzNUmhBBCxKuv5/8plf2DH2OPET2XzvBMCCGEeCwqbSEZBT+EEEKIKxJSt5gKfs6d40tcmEnx4oko+CGEEEJcjU7Hl68ATAc/588DtbX8vCBigIIfQgghxNUIQ14KBeBlYvourfgyiYIfQgghxNWYqubem9AjJGxL9Cj4IYQQQlyNJZOdqefHJJde6k4IIYR4pIgIICiIL3FhCgU/JlHwY0JbWxuau+ukmCOmnojUfRxxDABobW3Vf/UyNY7shHY5Yh+x5+6odsn12ks5Dl17eb5fPPnaA25w3+v+W9Vv+/Z2KDQaMMbATPw9k+t735rraAmPDn5yc3ORm5tr8JhWq3VSawghhBAb6u754bq6QCmBDYkOfsrLy/H6669jz549qK6uBsdxiImJQVpaGu69914kJSXZo512kZOTg5ycHIPH1Go1QkND4efnJ7oqsJQqwnI7Rlf30kl/f39R+zni3O29j9Rzt3e7HHUMR54/XXt5vV88+doDLnrfq6jgV3pFRfVb7aXfPiCAL3Ph7W1yRZjc3/titler1RZvKyr4+eGHHzBnzhwkJSUhIyMDGRkZYIyhtrYW27Ztw9///nfs2LEDV155pZiXJYQQQoilGANqavivgwaZ3k6hAKjIt1Gigp+HHnoId911F15++WWTzy9btgz79++3SeMIIYQQ0kdnJx/4cJzpHD9kQKKWuhcVFeHee+81+fyf//xnFBUVWd0oQgghhJjQe5k7xw287YULfImLlhb7t8uFiAp+4uLisGfPHpPP7927F3FxcVY3ihBCCCEmiClo2tDAl7iQsHrKnYnqL1u5ciXuvfdeHDhwAFdffTViYmLAcRyqq6vx9ddf46233sKmTZvs1FRCCCGEiAp+KNePUaKCn/vuuw+RkZF4+eWX8eabb+qXhSuVSkycOBHvvvsubrjhBrs0lBBCCCEQF/xQiQujRM+UuvHGG3HjjTeis7MT58+fBwAMGjQI3kJ0SQghhBD7aW/nv1LPj2SSp4l7e3vT/B5CCCHE0YYO5XtyLFnpRcGPUTYtbPqHP/wBQ4YMseVLEkIIIaQ3hQLw86Pgxwo2TRBw3XXX6YfCCCGEEOJkQvDT1dWTG4jYNvjpWyqCEEIIITbU3g5UV/OZm6OjzW/v5QWMGsUHQRT46Nl02IsQQgghdtTaCpw/D9TVWb6Pvz9lgu7Dpj+NsrIyPPXUU9iyZYstX9Yp2tra0NzcbNG2GgnJo8Tu44hjAEBra6v+q5cFvyyOapcj9hF77o5ql1yvvZTj0LWX5/vFk6894Fr3Pa6+HpxGA+bjA9bnb5Rcr72Ufaw5F0vYNPipq6vDO++84zLBT25uLnJzcw0eE3IXEUIIIbIjZpm7oKEBXHMzWHAwEBpqn3a5GFHBz3//+98Bn//999+taoyj5eTk9JunpFarERoaCj8/PwQFBYl6PbHbS9nH3sfo6uoCAPj7+4vazxHnbu99pJ67Pdul1WpRUFCA6upqpKamIj09HUql0i7tcuT507WX173Ck6894GL3PW9vICAACA8HTLxWv2PU1/O1vYKD++0j9/e+mO3VarXF24oKfhYsWACO48AYM7kNRxOqCLEJlUqFFUuXoqS8XP9YSmIiNmzejKysLCe2jBDiNEKmZl9fy/eh5e79iC5s+sknn0Cn0xn9d/DgQXu1kxCPolKpkJ2djTHl5dgLoAnAXgBjKiqQnZ0NlUrl5BYSQpxCTGkLAZW46EdU8DNx4sQBAxxzvUKEEPO0Wi1WLF2KeYxhG4DJAIK6v25jDPMArFy2jOanEeJpdDo+Xw8gLvihnp9+RA17/fWvf0VLS4vJ54cNG4Zdu3ZZ3ShCPFlBQQFKysvxAfp/OlEAWMUY0srKUFBQgBkzZji+gYQQ51AogPHj+R4cEXP/KPjpT1Twk56ePuDzgYGBmD59ulUNIsTTVVVVAQBGm3h+dJ/tCCEeRChtIYYQ/Gi1fO+RglL80U+AEJkRCgYXmXi+qM92hBAyIKWyJ+Ch3h8AFPwQIjvp6elISUzEcxwHXZ/ndADWcRxSk5KQfvnlzmgeIcRZLlwAzp4FRCzp1hs5Ehg3TtwqMTdGwQ8hMqNUKrFh82Z8DmABxxms9lrAcfgcwPpnn4Xy+HGAhr4I8RxqNV/aQkQmYz1Lq8B7CMnBz48//oj27kyTvb8nhFgvKysLeXl5OBITgzQAIQDSABQlJiIvLw9ZV13Fb1hZCfz+Oz+OTwhxb1KWuROjJIeBc+bMwaFDhzBkyBCD7wkhtpGVlYXMzEzs/OILVJ8/j9QhQwwzPPv4AKWlfPbW9nZg2DDA29vqrNCEEJmyJvhpauLvFYGBQGSkbdvlgiQHP73z+VBuH0LsQ6lUIn3WLABG0rwPGsR3ZZ85A2g0wPHjUBUVYcUjj1BWaELcDWPWBT8aDXDuHL/ii4IfmvNDiGxZ8qEiKAgYMQLw94fqq6+QfeutlBWaEHckrNLiuJ6l62JQrh8DNPvJhLa2NjQ3N1u0rUajEf36YvdxxDEAoLV7Il1rayu8LJgc56h2OWIfsedu13bpdFAcPQoWEABNdLTZhGba2Fgs37BBnxVa+FQjZIVewHFY8eCDuOqqq0wOgTni/Onay/Ne4cnXHnCR+15TExQaDeDrC52Jv00DHqOjg99fqzXYX67vfWuuoyU8OvjJzc1Fbm6uwWNUMoDIQnMz0NUFrrXVokyue376CWdravAhBsgKXVGBPXv2mE1WSgiRoe4eGyal1wegnp8+PDr4ycnJQU5OjsFjarUaoaGh8PPz6z/Hwgyx20vZx97H6OquG+Pv7y9qP0ecu733kXrudmlXYyMQEABERkIXEGB2+4aGBgDms0I3NDSYfB1Hnj9de3ndKzz52gMuct8LCgISEviVnWYCIKPHCAjg/wGAv7/+Q5XYc7d2QYU93y9qEfmPaM4PIXIk/BKHhFi0OWWFJsQDKJXS5vsAfIZnIUiR2PujUqkwLCUF1157Le68807MnDkTw1JSXHI+IQU/hMhNVxe/MgMAgoMt2sXirNA05EWI57Ji6EulUiE7O9ttFlRIDn4ee+wxRERE9PueEGKlpib+q7+/xZ/yLMoKvWkT5fshxFWdPcv/E5a7S3HRRXyJCws/VAm0Wi1WLF2qX1AxGUAQehZUzAOwctkyl5ozKzn4WbVqFcLCwvp9TwixksghL4E+K3RCgmFW6Lg4Pis05fkhxHXV1fGlLazJq+fjI6nERUFBAUrKy/EYTC+oKC4rQ0FBgfS2OZhHT3gmRJYCA/mMzSKDH6BXVuidO1H9229I9fPjJySOGmWHhhJCHKKrq6eEjRNKW1R11xA0t6CiyoVqDVLwQ4jcDBrE/5NIqVTyc3uuuAJBxcV8EcSuLipqSIirEmpnenvzSQ6lamnhK8P7+gIxMRbv1ntBxWQjz7viggqa8EyIu/Lx6Zk3ZIfCw8KS148//hj5+fkuNd5PiEuxVUHTjg6+xEV3agxLueOCCpsEP2+++aYtXoYQ0txs2yRkF10EjB3LD6XZkDsteSVE9mwV/Ehc7WWwoAJwiwUVNgl+9u7diwceeAC67jHJkydP4rbbbrPFSxPiORgDTp8GDh/uWepuLak5QQbgbkteCZE9Ifjx9bXudYTgScKKsaysLOT98584Eh1tuKAiMdElF1TYZBLA22+/jY0bN+Laa69FWFgYTp06hb/+9a+2eGlCPEd33R0olfxwlS0xxk+YtPKTWd8lr8ZqiK1ctgyZmZku9SmQEFnrzsJss54fxiTNA8y67TZkzpuHnfn5qK6oQGpQENKvvRbK+Hjr2uUENgl+Dh48iB9//BE1NTX47bffsGvXLiQnJ9vipQnxHEJ+n+Bg6yY19nX+PFBRAUREAElJVr2UsOT1AwxQQ6x7yeuMGTOsOhYhpFtqKjB4sPX3BY7jA56uLn7oS+wiCC8vKKOikD5nDrjaWgTW19uul9rBbBL83HfffXjyySdx7bXX4ueff0ZmZiZyc3Nx5ZVX2uLlnYKqusu4urGd9nF2dWOuqgqcRgNdRAQ/98dWx2hrg0KtBtrboQsPN7mPJedfXFwMwPyS1+LiYkyaNGngdlnIE669Lfehqu6ee9+zZHuusxNcayt09fWAViv92iuV4DQaoK0NupgYs4GZW1Z137dvn/77yy+/HF988QWuv/567NmzxxYvbzdU1Z3Ihk4HrqWF/15k9lWzQkL4G1N7O7/s3YohtdjYWADml7wK2xFCZEYIcMQurLhwAVxXF5iQ0Njfv6cXqaWFL7zqQiQFP4sWLcKMGTNwxx13AADOnj2LY8eOIS0tDaGhoUhISMB3331n04baA1V1788lqhvbaR+nVvZWq/mbiY+PyRw/Vh0jNpavFK/VmrxJWXL+GRkZ/JLXigpsY8xg6Eu/5DUxERkZGQPO+aFrL697BVV1l/F9r70dQVVV/P0hMdH6Y1xySc/wF0Sce3k5H+SEhuqrwwfExAD19fwcIgvPy6Wrun/11VcYMWIEAKC+vh4TJkxAVlYWRo0ahZMnTwIA/Pz8pLw0IZ5JYkkLi4WG8l8bG616GXdc8kqIrLW18fcHEX/YB+TtLX6uj1bbM7en9z1K6KUW5iu6EEnBT2NjIxK7I9D//Oc/iI+PR2NjI26++WasWrXKpg0kxCPExPCTGq3I7DwgIfhpbu5ZOSJR1h/+gLwXXsCRmBi3WPJKiKzZapm7NZqa+N4dPz/DFWdC8NPc3FN+w0VIGvZKSkpCcXExkpKSoFKpsGjRIvj4+ODuu+92qQyPhMiGtze/GstefHz4rmqNhu/9iYyU/lq1tciaNQuZ11+Pnfn5qCkuRsqYMUi/7jrq8SHExjhhbo6tanq1tQG1tXzai4QEy/Yx1TPt5weEhfFDctYUXHUCScHP4sWLcf/992Pu3Ln47rvv9JOGtVqtxSukCCEONmgQ/ynSmmzP7e36oTNlXBymTZsGbuhQBCQkWJ1DiBBihK2yOwu0Wr7EhY+P9cEPAAwdapt2OZik4GfVqlVgjGHnzp14/vnnMWzYMADA/v37MXjwYJs2kBC3V1PDfw0Pt2/F5qgo61/j3Dn+a0gI4OsL5ucHDuBXkRFCbE+oy2ere4PYEhft7fw/jnO5FV0DkRT8cByHxx9/HI8//rjB4zU1Nbj55ptt0jBCPEZtLf/pTljtJVc6HZ8wEQCio/mv/v5gISF84EYIsTmus5MPWGwd/AhZns1pawMUCn7Y3FTvrlbLzwsKDnaZHmCb5PkRUEkLQkRqa+MDH0d9qtLp+JuUVit+jlFLC7+/r29P97e3N9jQoW71iZAQ2WCsZy6NrYKf3lmeLanxFRoKXHrpwD1FJ07w97KhQ/k5QC7ApsEPIUQkYYloUBD/6cre1GrgzBn+Rio2+AkO5ivEC13ghBD74jjoxozhe4Vt2aPi49NT4sLCdgwYfAUH88FPU5PLBD8OuNsSQkyyd36fvkJC+CCro0PaPB0vL+MTpjs7e+YmEEJsy9ZDSZbO+7F0BZcL5vtx6eDn+++/x/z58xEfHw+O47Bt2zZnN4kQyzHWc7NwVPCjUPTcqBoaLN9voJvkhQvA4cNAaalVTTNGq9WioKAAH3/8MfLz86n8DCG2YGnwU10NHD3K/44PRLintLaKL5vhJC4d/LS0tGDcuHF49dVXnd0UQsTTaPi5N0qlVfW2RBO6pS3N9tzZCRw5Avz2m/FEZkLyNRuv+FKpVBiWkoJrr70Wd955J2bOnIlhKSlQqVQ2PQ4hsnXuHLjTp80HH2IlJPDzeOLiBt5OreaHs8z1AHl59dzDXKT3x6WDnzlz5mDt2rWUUZa4pvZ2vidGKDzqKEK255YWyz6lnT/fM/HS2Lwk4abX2ckHczagUqmQnZ2NMeXlBiU0xlRUIDs7mwIg4hE4jQZcU5NlE5PF8PIyP5Sm1fL3CMCyYstC77WLBD82n/CsUCgwY8YMvPTSS5g4caKtX94q7e3taDczL0EojNbc3IwGC4cFNN01T7pElA0Qu48jjgEATd1v3CYL38COapcj9hF77la3KyAASE7mbzIDvNfsce6Kri5wra3QlZaCdWd7Nnr+jEFx5gy4ri7oIiLA+rRTOI6urQ1cZye0VVUDrvyy5Fy0Wi0eeuABzGMM29DzCW0ygG2MYQHHYfmDD2L69OkmM0rL/trL7F7hiPOX688LkOd9T6vVYve33+JcZSUSxo7FlNmzzWZQt+W5c42NUKjVYD4+0LW2GvTsGj2OVgtlUxNYezt0wgcsK9pmzblYgmPMtjmp3377bZw9exY7d+7Ejz/+aMuXHhDHcfj000+xYMECk9usXr0aa9assej13n//fQR0V64lxN341tXBr64OHaGhaB0g+aF3czMCqqvBvLygTk422UMVWFUFr5YWtEZFocPIjU+MI0eO4G9/+xv2gg94+toLvpbYM888gzFjxlh1LELkaO/evXj7rbdQ02u4KyYyEovvugtTpkyx+vW5ri741dcDOh1aY2KMbuN/7hx8GhvN3iP0tFr4tLSgy98fOmFOkYNpNBrcfPPNaGxsRIiZeZQ2D36cxZLgx9Ken6SkJBw9ehTx8fEWHVuIUMUES2L3ccQxAD5yPnjwICZMmIBgC7o6HdUuR+wj9tytahdjCLCwzIRdzl0YovLz0z9k7PwVp06Ba2mBLiYGzMj8AOE4gQ0NUNTWQhcZCZaUZNW55OXl4e6770YTAGN9SE3gi6n+4x//QHZ2tuTjGLymI6+9nfdxxO+9o9rlife97du3Y9GiRZjHGB4DMBpAEYDnOA6fA3jnnXcwf/5869rV2Qnl0aMAgIYhQ3CwsLDfuSuOHwfX3g5daipYnw80cn0f19TUYMSIERYFPx6V58fX1xe+FlbGDQoKQpiF+Qq8vLz0+1hK7D6OOEZvwcHBFp2/o9rlyPO39NytaRdXXo7AxkZ+wqGZ7MhOu/atrT1zkoYN61khYuw4fn789j4+A+b5sKRdQrmcIhjv+SnqtZ2p6yTna2/vfRzxe++odnnafU+r1eLJxx4bcMj3qccfx80332x0CMzidjHG/14zBm33nD2Dc+/o4H+XfX2BxMR+84Pk+j4WAiZLiJ7wXF5ejscffxwzZ87EyJEjMWrUKMycOROPP/44ysrKxL4cIZ6puVn+9bCEUhZhYUYDHwOBgXzJCxvUD0tPT0dKYiKe4zj0XVumA7CO45CalIT09HSrj0WInBQUFKCkvByPof8fZwWAVYyhuKwMBQUF1h1IyPIMGF/0wBhfCDksTFyOIZ2OL9dTUmJd+xxAVM/PDz/8gDlz5iApKQkZGRnIyMgAYwy1tbXYtm0b/v73v2PHjh248sor7dVeA83NzTh9+rT+/8XFxTh06BAiIiKowCqRJa1Wi4Jdu1Dzyy9IiY9H+ujRcFolnM5OoKyMD8IuuaT/84mJfFDTa2jMJF9fYIDhLjGUSiU2bN6M7OxsLOA4rGJM3/W/rrvrP2/TJrOTPwlxNVVVVQD4oS5jRvfZzio+PkBnJ187rC9fX34xhhTl5XzwFBtr2b3DSUQFPw899BDuuusuvPzyyyafX7ZsGfbv32+Txpnzyy+/YObMmfr/L1++HACwaNEivP322w5pAyGWUqlUWLF0KUrKy/WPpaxdiw2bNzsnXYNSyef60en4nEN9cZz4Ehg2kpWVhbyPPsKK++9HWm2t/vHU2Fjkvfoqpbcg7qezE3Hd0zLMDfnGmcvPYwmx1d0toVDwqz2bmvh/7hL8FBUV4V//+pfJ5//85z/jjTfesLpRlpoxYwbcZL42cXNC3pp5jOED9JrE2J23Ji8vz/F/0IX5PA0NfBBkbaJFrZbvRbJR0sasjAxkbt+Orw8eRJVOh1SlEumXXQbl+PFWvzYhssEYcO4cUFmJ9MREpMTH47mqKmxjzGDoSz/km5homyFfU8FPRwf/WECAtPxjISE9wY8NhsHtRdScn7i4OOzZs8fk83v37rVNREqIG9FqtVixdKl+EuNk8KuYhEmM8wCsXLbMOaUbhFUcvbM9NzYCRUX8DVmM2lrg5EmgpsY2bWtuhlKpRPq0abj+ttswY/JkKHuXBCHERZgs09LUBBw/zg8/a7VQhoRgw/r1+BzAAo4zSPC5oHvId72thnxNBT8XLvBV2qXO23GROl+ien5WrlyJe++9FwcOHMDVV1+NmJgYcByH6upqfP3113jrrbewadMmOzWVENckTGL8AKYnMaZ1T2KcMWOGYxtnJNuz4sIF/tOo2EKlQm+PrSZyd2eXZcIn0PBwPiBranJcLTRCrGR0uDshARsefxxZl1/OP+DlxZecGDQIWSNGIM/XFyuWLkVar31SExORt2mT7XqIo6OB6GiwpiaguLjncaHYsoTVcgD4HiOlkq8ar9Hw/5chUcHPfffdh8jISLz88st488039dGrUqnExIkT8e677+KGG26wS0MJcVUOncQolrc3P6m5pQWcWg1FZyc4tZr/9Ca2y1oIfoRaQNaU7GCsJ7W+cBOOieHb5Mg6aIRYweRwd2Ulsu+7D3kvvoishQuB+HiDVVVZWVnIzMzEzp07UV1djdTUVKSnp9t2kr+x1+pd0kLqBwyO439nGxv5DyruEPwAwI033ogbb7wRnZ2dON+9FHbQoEHwdlJGR0LkThgKdsgkRilCQ6FVq/HDd9/hx/37oRw1CrPnzIHSwpxYej4+/DwinY7vNbJmsqPQ66RU9ryO2PYQ4kR9h7uN5exZ+coryFy+3GhQo1Qq9XN7pOQskqS5mf/g4etr3e9bcDDfgyTjCu+Skxx6e3vT/B5CLKDPW1NRYf9JjBKo8vOxYvlylFRXAwBeAJDy0kvY8Pe/i+ti5zg+UNFo+KEva4IfPz++6rQw+bIva3uWCLEzi4a7y8udM9wN8B9Syst7hrmBniEva4eVBw3ie2mNFUKWCZu27A9/+AOGDBliy5ckxOUJeWscMolRJJVKhexbbsGY6mrD6ulVVdKqp/ce+rIWx/X/9KnV8vMTiop6btiEyJCsh7sBPjA5fx5cYyM4oXiorYIfpVLWgQ9g4+Dnuuuuw6JFi2z5koS4haysLOTl5eFITAzSwNemSgNQlJjonGXusNMqNFtPeu5LqeTnEXR09NyoCZGh3sPdxjh9uBvQr/hSaLX875TwocXC+m4WkemHFJvW9srJybHlyzlVW1sbmpubLdpWTD0Rqfs44hgA0Nr9R6u1tVVfW0UO7XLEPmLPXewxMjIy8Ou2bcgvKECVlxcGDx+OtLQ0KJXKAd9r9jp3S1eh7dy50+SQXL/jKBTgIiL4FVpGzsmic+nshOL0abDAQLDBg/vtw/n4gGtsBCsvBzPRWya3a+/IfRzxe++odrnyfW/8+PFITkjAc5WVJoe7U+LjMX78eJO///b+GXMdHWhrbwfX1YXWri54JSTwH1zMfHix6DgtLeDKygClEuyiixz6PraERxU27Ss3Nxe5ubkGjzkl1wrxGF6MYeq4cei66CIEDBrk1LZUd8/xMdctL2xnEX9/MGtXY7W08J9ATczpYWFh4M6d4wMgnU723evEMymVSjz7wgu47dZbsQDAKqBfmZb3XnjBuWVafHwAAIqurp5VWraaXO3lBa61FeA4/vdUZkQHP+Xl5Xj99dexZ88eVFdXg+M4xMTEIC0tDffeey+SbFTfxxFycnL69Vap1WqEhobCz89P9Ax7KTPy5XaMru6xX39/f1H7OeLc7b2P1HMX1S4fH8DPD7rQUKdf+9TUVADmV6GlpqaaPbZNz6WhgV8eGxVlcCPW7xMUxOf76ejgJ20OMD9BVtfewfs44vfe3u1y1D72uu/d8sc/wv/FF7Fi/XrDMi0ic/bY7WccFgbdhQtQaLW2v/ZBQXxh1I4O/v/dS94tOYaQFFLsMn+1iKFwly5sSohL6erqGf+2cGjBnuy2Cq2tje+98feXluNDyDMSGGh6m/BwPpN0fT1/gyVEjmpqkDVrFjKvvx47jx+3X84eqbrn/HhpNPwQFcf1JD61heBgPmN0U5PFr2s0KWRios1rILp0YVNCXIqwokKplMUybbtVT6+t5XtmYmLEBz/GkhsaIwQ/DQ187w8NfRG56erSz3tTxsUhvXuY22E5eyzRHfwo29v5Je8BAU4NfhxZA1HUHaOoqAj33nuvyef//Oc/o6jI1Nx2QjyckK9GRglB9avQEhJstwrNmuXuGg0fAHl5DZxkLTCQv5nGxsp2NQnxcF5ewJgxwNChsk3QqQ0Jwe7GRuw8ehQFhw5BO1BvqxTCqrGWFj5NxUBtcXANRCpsSoijKJVARASYzOpSZWVl4XRJCbZv347ly5dj+/btOFVcLP0TljXL3YVVL5bchIcNA+LijKfpJ0QOlErZDsuqVCoMGzIE8xcswLqtWzH3oYcwbNIk8bm9BtI9xxGA2UKnwurTx2B69Wlxdw1EW6DCpoQ4SkAAkJoKZmEKBUdSKpWYOnUqtFotpk6dat18BCH46ejgP+2JeS2Fgr9ZymlogBCxOjr0K6nkaMCaYzYeXkJ4OLQaDQp+/hnVjY0m5zw5OikkFTYlhNiWUsnf+Ds6+N4fMYFMVJS4gqo6HT/vx8eHAiYiD4wBJ07ww9tDhshuyMuimmPLliEzM9Mmk7JV+/aZnsB83XX83Khe5bIcVQNR9CzBG2+8Efv27YNGo0FFRQUqKiqg0Wiwb98+CnwIGYhW6znzU+yd6VlQVcWXu+i1jJgQp6qr4+f3dXbKsvfHkcNLQg/TmPJyw/I53ROYVS+/DJw5A6DX6lOOQ9+sQPrVp0lJNquBKHmJhFDYNC4ujiq6E2KJ4mLg4EF+9YO7kxL8SAkOw8P5r42NfC8QIc5WU8N/jY6WxarOvhw1vGR2AjNjWPnii9C2tABdXQ6vgUjrQwlxlN5L3d1dRARw0UX8hGRLVVYChw6J68UJCOCHFYThL0KcqamJD/gVCr6yuQw5quaY2R4mAMU1NShQq/V5z0yuPk1IsHkNROdnWiPEU8hwqbvd+Pv39P5YqrmZD2LEJoAMDweqq/mEhxER4vYlxJaEXp/ISFkkMjXGbslN+7C4h6lP+ZysrCxkZmZi586ddk0KKc+rQ4g7Enp+ZHpTdCqdrmeITOzE5YgIPvhpbBS/uowQW2lr49+DAD/kJVN2S27ahzUTmJVKpT74sldSSMnDXj/++CPa29v7fU8IMUKn65mT4gk9PwCgVvNDWZZUZxaSG3p7i58k6u/PL49nrOePDyGOVlfHfw0N7cltI1N2SW7ah80mMGu1/FD42bNWt6k3yR9B58yZg0OHDmHIkCEG37uLtrY2NFuYj0Vjyc3dyn0ccQwAaO3+9N3a2govC3ooHNUuR+wj9txFHaO9HQqNBlAooBGZ+Viu197ccbizZ8E1NIAlJIB1fxI2tT1XUwNOowELC+uXB8mSc+F8fMDV1YGdOwfm4yOva+/gfeRw7W3VLrm+940eIySkJ4A38rdDbtc+IyMDh4qK8M0332DPnj1IS0vDH/7wByiVSrN/+yw9ztrnn8dtt91msofpvXXr9D97k8fo6oLi5EkAgC4oaMDUAaZeyxjJwQ/rtSqDuejy3dzcXOTm5ho8ZqvU2YQYEOb7eNKQl78/PwnZkhuSmMzORrBBg8AiI2W5tJh4EFvWxXIApVKJtLQ0eHl54fLLL7f5vJrMzEy89957ePyRR5BWUaF/PCU+Hu+98AIyMzPNv4iXF1hQELjmZnCNjfoPUtbyoDtxfzk5OcjJyTF4TK1WIzQ0FH5+fqLHGqWMTcrtGF3d81L8/f1F7eeIc7f3PlLP3aJjKJVAUhKgVELXXezTXa69yeMMGsQPQykU/ebxGN0+IIAvhmoiAHLZa++kfZx67W24vSP2seq+xxj/z8Liup527W+55RYsXLhQ8gTmoKAgIDERKC3lE6cOcEy1Wm1x2z06+CHEYfz9gZQU/nsZlrewi94FThkznfOEMX6CaEuL+CrwxlC+H+JI9fVAWRmf1kHGE52dyeoJzGFhfPDT0mKz0iGU54cQYh8+PvynYZ2Ov2GZwnH8H45hw6xLCqfTAadPA7/+araCNCE2U1PDr+SkoNt+vL17KsTX19vkJannhxBH0Gr5QECGGV/thuP4VS8aDT/vx941jhQKoKMD2s5OFOzYgerWVrvlCCEEAJ/UsHshg1yTGrqN8HD+511fzw+PW4mCH0Ic4exZ/pd28GDxyf9cmb9/T/ATFmZ8G7Wa384GKQBUe/dixWOPoaRXlmh9EUUbZoclnkur1aKgoICfv+Ljg/SLLoIyKsqzFjM4Q1gYP7zo5TXwMLqF6GoR4gieuNoL4Iez4uJMj9FrtcCpU/z3Y8daFQCpVCpk33035jGGDwD9strnuoso2jo9PvE8KpWqf4Xy6Gg+uB482Ikt8wDe3sCll1o8sdwcya/y2GOPIaI7lXzv7wkhRnhqdmdfX/6fqU9pLS38Vx8fqwIfs0UUAaxctoxSWRDJTFYoP3cO2TffDJVK5eQWegAbBT6AFcHPqlWrENbdjd37e+I5hO7fjz/+GPn5+fSHZSCeVNdLDCH4sTKFvdkiioyhuKwMBQUFVh2HeCYKrmWmo8PqRQ0WBz80YZD0plKpMCwlBddeey3uvPNOzJw5E8NSUujTjzGM9fyielrPD8Cnpi8u5pe892VlckOBxUUUu7cjRAwKrmWktBQ4cgS4cMGql7E4+OmdxbmpqcmqgxLXZrL7t3tuBQVAfQhDXhznmcFPfT1f90jo5enNRj0/vYsoGjNQEUVCzKHgWkaEumlWLnm3OPjheo3Zp6eno7pPGXriGaj7VwJPnews6J3ssLe2tp4UAFaugLOoiGJCgvkiioQYQcG1jAhTbJqbe+6tEkia8zNp0iRcccUVOHHihMHjhYWFuPbaayU3hsgfdf9KoFAAkZGml3q7OyGw6Vvjq/eQl5XLVpVKJTZs3ozPASzgOIMeyQXdRRTXL10KZa8l8IRYymYVyon1fHx6hsmt6P2R9FH0rbfewpo1azB16lRs27YN0dHReOKJJ/DJJ5/gj3/8o+TGyAlVdTde3bi4uBiA+e7f4hMnMGnSJJu1y+WrugsJ0JqbXfbaSz6OVstXtO/qgiYkpOdxpRKIiuK/DvC7Zum5ZGRkmCyi+K9HH8U148ZBc/o0WEMDWGKiQcBFVd2pqrvw1dT5W1OhXGq73OXaS9lnoO05Hx9w586BlZeD9SqJ45Cq7k899RR8fHxw9dVXQ6vVYvbs2di/fz8mTJgg9SUdjqq6ixcbGwuA/6WfbOR5ffdvRwe44mKw1FSD57VaLX788UfU1NQgOTkZaWlpNJne3Qk9P8IKDeF6e3sDNk6RkZmZiXnz5uG7777r9x5j586BKy8Hd/480NkJlpJi06WzxL3pK5T/9a9I6zW3R1SFcmITLCwMXEUFuJYWsM5OSatoJQU/VVVVWLduHd566y2MGjUKJ06cwMKFC10q8AGoqrsx5ir8ZmRk8N2/FRXYxpjB0Je++zc2FldPngxlRITBRFajCcJEZt91ycreJkpbuNq1t+o4YWFARwd0HAcEBNj9XK6++ur++wQF8Snyf/+dnytQWQkMGwYtxyE/Px/ff/89lEol5syZIyogl2OVcrHbU1V3y87/lltuwcJp0/D1p5+iqr0dqZddJrp8iidfeyn7mNw+KopfMNHVxf9eQ1xVd0kfe4YMGaLP73LgwAGoVCrcd999eOGFF6S8HHEhFs2tePVVKMeNM6i/ovr3vz13hVhpKXDwIL/k21P5+fHBn1DgVKMBqqv5r44UFgZcfDHf+9TSAtXrr2NYSgrmz5+PjRs3Yv78+ZSygQxI2dWFaePH44Y//QkzZsygnmtniYsDUlMl9x5LCn62bt2KwsJCzJ07FwAwe/Zs7Nq1C5s3b8Z9990nqSHEdWRlZSEvLw9HYmKQBiAEQBqAosTEnhICvr761U1arRYrVq703BViwlJ3T75JpqYC48f33KgaGoCKCucEhEFBwIgRUBUUIHvpUs8MyIl03asWmamSLcQxQkP5+4nE+6rFw1698/wsXLiw3/MTJkzAnj17aLWXh8jKykJmcjK+/vFHVPn4IHXECJPdvwUFBSiprsYHML1CLK17hdiMGTMc0HoHo+zO/Zf5CxOcrczvI5XW2xsrNmzQB+TC+1IIyBdwHFYuW4bMzEz6ZE8MCSkbhHwzxCVZHPzodH0X+PWXkpKCH3/80aoGERfR1gYlgGkTJ0I3diyCgoNNburxCcI8ta6XKYz1JDe0MrOzVELKBo8NyIk0nZ2A8LfQ19e5bSH8vfX8eX44XWRhWZsvdQjvnnhE3Fz3XA3m7282R4vHJwgTgh9P7vkBgLNnwZ04AajV/B8QpdJpn549PiAn0nRXFtcNH251bipiA1otP3x+7lzPfdZCtM6TSCN8cu+VY8EUj04Q1tXF93QA1PPT3AyutRWcUJPHBskNpfL4gJxIp1RadN8jDuDr23MtGhpE7UrBD5GmO/hhFtwELFohtmmTe86tEOb7KJX0SbE73w/X2Mj/30lDXoCHB+SEuBNhtElktmcKfoh4jPUsUbbwD5h+hVhCgukVYu5IKG1Bw8HQ+vjg+8JC/Oebb5D/yy/QWlnPyxoDBuQAH5C//LJ7BuREuqoqoKysf6kW4jzCvbWpSdTQl4f3wxNJOA4YPZoPgEQM5WRlZSEzMxM7d+5EdXU1UlNT+RViwtCQO/aM+PoCKSnOboXTqVQqrHjgAZRUVuofS1m7VlSCS1sTAvIVS5cirVfizdSYGOStWIGsadOc0i4iY3V1/Gqv+Hhnt4QIhKEvjQYQepUtYPOeH4VCgVmzZuHAgQO2fmkiJz4+kgp1KpVKpKen4/rrr+cThJWVAUVFot60xLWoVCo+wWVlpezy6WRlZeF0SQm2b9+O5cuXY/v27Th16BCyZs1yfAJGIm+MAe3t/Pe00ktewsOh1WqxNz/f4l1s3vOzZcsWnD17Fg8++CAteyfmCSugamrcs+q5idIWnkKr1WLF0qWyzqejVCoxdepUaLVaTJ06FcrQUD7/kJNyEBGZ6ujgAyCFgv/wR2RDtXs3Vjz4IEpEJE21ec/P4sWL8dRTT1Hg485KS/mxb5FLC42KjuYDg+Zm9/ykXV7Ol7aornZ2S5xCyKfzGEzn0ynuzqcjGxxHgQ/pT0huSL0+sqJSqZB9000YU1uLb0TsR3N+TGhra0OzkIXWDI2EP9pi93HEMQCgtXsiX2trK7yMzefRaqE4exYAoPP3h0ao1WRFuzgfH3D19WC//85X2rZgHynHMcfsuUs4BtfQAE6jga6tTZ/V2GWvvYTjFBcXAzCfT6e4uBiTJk2SdAxb7GPy3Ds7wVVXgyUk9KsA76jr4oj3iz2uvS3aJaf3PldXB06jAfPxke25yPXaS9nHku21Wi2WP/igvmfZsr/YPEk9P4sXL8b3338vZVdZyc3NxahRowz+XXbZZc5ulrwJb0gfH5vlrWHR0QD4QAESgilZ8/DszrGxsQDM59MRtpMb7swZcOfPg6Nkh4TKWsjOnj17cLaiwmjPsjmS7shNTU3IyMhAUlIS7rjjDixatAgJCQlSXsqpcnJykJOTY/CYWq1GaGgo/Pz8ECSy61vs9lL2sfcxurr/WPv7+xvfr7mZn1kfHm4wNGBVu4KC+ARVTU38EtIBqvTa8/zNnruUYwhzA8LD+6UFcLlrL+E4GRkZfD6digpsY8zgBqXPp5OYiIyMDLNzfpxy7S+6CDh9ms9rpVTqcxU5ql3W7OPsa2+r7R2xj0Xn7+3N3/siIsC6h77keC5it3fktZeyz0DbN3QnNjTVszwQST0/n3zyCSoqKnD//ffj448/RkpKCubMmYO8vDx0CkndiHuyV02m7t4fsYmqZM/De35cPsFlaCgfuDIGnD3bk62beJ5hw4Bx49xzYYaLMpepfSCSJzxHRkZi6dKlKCwsxM8//4xhw4bhtttuQ3x8PB566CGcOnVK6ksTObNX8BMayufDGTXKtq/rTDpdTxFED67r5fIJLpOS+F6flha+iCLxXF5e/HuByMJAmdrNsXq1V1VVFXbu3ImdO3dCqVTi2muvxdGjRzFq1Ci8/PLL1r48kZPOzp5yDbaubcNxfCZkhRslHRd+VgqFe52XBEI+nS+//BJbtmzBrl27cKq4WP6BD8AHrkJSu4qKnutKCHGqvj3LP4vYV1JffGdnJ/773/9i69at2LlzJ8aOHYuHHnoIt9xyC4KDgwEAH374If7yl7/goYceknIIIkdtbfwfcV9f+/8x7+py/aEiIaAjAHoSXALS5go4VVQUn923pYVP8zB4sLNbRBxJrQZqa4GQkJ4heiILvTO1X90rU7s5kv66xMXFQafT4aabbsLPP/+MSy+9tN82s2fPRhiNjbqX4GDg0kvt+8m3pYWfW+HlBVx8sf2O4wg+PlTawl1wHB/w1NUB8fHQarUoKCgwLNNCwyHuq6WFz0Lv6h/I3JRQOmnbtm3Izs62aB9JV3Lp0qVYsWIFAvoMfTDGUFZWhsGDByM8PFyf44O4EY6zb3ZTb2++h0konmrr4TVCpAoIAAIC+DplS5eipNenzJTERKfWKSN2JpS1oGXusqVUKjFlyhSLt5c0drF69WqjCQDr6uqQmpoq5SUJ4fn49FTpralxbluspdX2THgmbkFfp6y8XHZ1yogdUXZntyMp+GEmlns2NzfDjyJj99TWBhw9ype2sLeYGP5rfb1rTy6tqAAKC/k5IsTl9a1TNhlAEHrqlM0DsHLZMmi1Wmc2k9gDJTh0O6KGvZYvXw4A4DgOTz75pMGwl1arxU8//WR0/g9xAy0t/A3AEWPeAQF84sPmZn6SoQsm0ATQk+OH5oK4BaFO2QcwXacsrbtO2YwZMxzfQGIfXV18Ly5APT9uRFTPT2FhIQoLC8EYw5EjR/T/LywsxIkTJzBu3Di8/fbbdmqqca+99hpSU1Ph5+eHiRMnyqtAojsRylrYOr+PKULvz7lzrjt0JPRaeXCOH3dS1d2DZ65OWRX19LkXodfHx8fjU1a4E1Ef43ft2gUAuOOOO7B582aEhITYpVGW+uijj7Bs2TK89tpruPLKK/Hmm29izpw5OHbsGAbTUlTbEpIbOmoCcmgo/ymrvZ0vfWHPSdb24uHZnd1N72yyk408X9RnO+Imurr43lsa8nIrksLYrVu3Oj3wAYCNGzdiyZIluOuuuzBy5Ehs2rQJSUlJeP31153dNPcirLwCHNfzw3FAYiJw0UXQhoaioKAAH3/8MfLz811nTgX1/LiVgbLJ6uuUJSXpcxkRNxEWxqf4GDrU2S0hNmTxR9Lly5fjmWeeQWBgoH7ujykbN260umHmdHR04MCBA3j00UcNHs/IyMCePXuM7tPe3o52YcmiCWq1GgA/eVsommaOpjswEArE2WMfRxwD4IvW9v7a/UJQqtVgSiV0ra188VEHtWv7Z5/hiUcfRWllpf6xwfHxWPv885g/f77NjgOYOHczTB6DMSi73z/a5uaernMJ7XLqtbfxcRx1Lva69k8/9xwWLVqEBRyHVYxhNPgen3UAPgfwzrPPDvgacr1XePK1B8Sfv1zPRa7XXso+1pyLJSwOfgoLC/VFSwsLC01ux3GcxQe3xvnz56HVahEjzA3pFhMTg+rqaqP7rFu3DmvWrLHo9X/99VePrk928OBB/fc+jY3wP3cOXQEBaKmrc1gb9u7dixdfeAFzAXwE6P/QPFtZiUW3346HH3lEVF4HS/U+d6m4zk6EnD0LcBwaGxtt0CrHscX5uypz5x4SEoKHH34Yb7/1FtIuXNA/Hh8WhkfuugshISHYvXu3vZtpN5587QHPPn93OHchYLIEx0ytW5e5yspKJCQkYM+ePQZ/AJ999lm89957OHHiRL99LO35SUpKwtGjRxEv1PMxQ/iB9036aMt9HHEMgI+cDx48iAkTJuhLlXAXLoCrrQULCwPrM5/BXu3SarWYMHYsxlVWYhsMx2d14Ou4HI6Px4FffzWZWVds24yduzkmj9HZCUVVFcAYdMnJVrXLmdfe1sdx1LnY+9prtVrs2rULtbW1GBwbiynTp1uU4Vmu9wpPvvbAAOfPGBQnTwK+vtANHqxfuSnXc5HrtZeyj5Rj1NTUYMSIEWhsbDQ7NUfSTMzW1lYwxvSNOnv2LD799FOMGjUKGRkZUl5StEGDBkGpVPbr5amtre3XGyTw9fWFr4VLFYOCgiwuz+HVPaFVTL0isfs44hi9BQcH95x/WJjJ8W57tSs/Px+llZX4CAMsK66owJEjR0wuK5Z6/gbnbsaAx4iKskm7nHrtbXwcR52LI679Nddc45B2OeJnLPDkaw8YOf+ODn7OHmNARAQ/F9GB7XKXay9lHynHENPzI2nCc2ZmJt59910AQENDAy6//HJs2LABmZmZDpts7OPjg4kTJ+Lrr782ePzrr79GWlqaQ9pA7IeWFROXo9PxealcNTUD6a93ZmcHTekgjiEp+Dl48KB+RUNeXh5iY2Nx9uxZvPvuu3jllVds2sCBLF++HG+99Ra2bNmC48eP46GHHkJpaSnuvfdeh7XB7el0/KceB+u9rNgY2S8rptIWnufkSaCsDDh/3tktIbZCmZ3dlqRhL41Gox8b3LlzJ7KysqBQKDB58mScPXvWpg0cyI033ogLFy7g6aefRlVVFUaPHo0vv/wSyX3mWBAr1Nby5RliYgAL50DZgn5ZcUUFtjHWb87POo5DamKifJcVV1XxtcliY103QzURZ9AgvvxLdTX/PSXEc31U0NRtSfrtHDZsGLZt24aysjJ89dVX+nk+tbW1Ds//c99996GkpATt7e04cOAApk2b5tDju72WFr4Hw8ElGpRKJTZs3ozPwU9u7l1EcgHH4XMA6zdtsmiSqVMIOX4owaHnGDSIT8bZ2clnJieujwqaui1Jwc+TTz6JlStXIiUlBVdccYV+tdXOnTsxfvx4mzaQOJmQ2dlRyQ17ycrKQl5eHo4kJCANQAiANABFCQnIy8tDVlaWw9tkMcru7Hk4rqd3tLraKcOeWq3WNROCyhUNe7ktSXfm7OxsTJ06FVVVVRg3bpz+8auuugrXXXedzRpHnKyzs6cHw1FlLfrIyspCZmYmdu7cierDh5E6ZAjS58+HUu43I8ru7JkiIvghz/Z2fsg4NtZhh1apVFixdClKysv1j6UkJmLD5s3y/qAgV4zxH166uij4cUOSP5bGxsYits8v9uWXX251g4iMCL0+/v5Onb+gVCr5uT3p6ZKWcDoF9fx4JqH3p7iY7/2JinLIkLFKpUJ2djbmMYYP0JMQ9LmKCmRnZ8u/p1SOOA4YOdLZrSB2IvnO3NDQgH/+8584fvw4OI7DyJEjsWTJEoSGhtqyfcSZnDjk5fKE4Id6fjxPeDjf++Pv75D5clqtFiuWLsU8xgwSgk4GsI0xLOA4rFy2DJmZmfKdI0eIg0n6OP/LL79g6NChePnll1FXV4fz58/j5ZdfxtChQ90iRTbpJsfgR6cD1Go++ZhcdXX1pAegnh/PI/QYDBnikOC3oKAAJeXleAymE4IWl5WhoKDA7m0hxFVIujM/9NBD+OMf/4h//OMf+iyMXV1duOuuu7Bs2TJ8//33Nm0kcZLQUH64S07BT3Ex0NDAV3w3kclbFiIj+UCNEqN5JgcOE1NCUDspLQWam4G4OL43j7gVyT0/jzzyiD7wAfhU1A8//DB++eUXmzWOOFlMDDBsGN99LxdC7Rk5Fwv18gJSUvhP/sSztbcDZ8/2DIPagcsnBJWr1lb+H3FLknp+QkJCUFpaihEjRhg8XlZWZnFhNLlra2tDc3OzRduKqScidR9HHAPg67YJX70sGLJxVLv0+ygUUGg0QGsrdLGxA37CFnscsecu5RhS9pHrtZdyHIe/Xyxkr3PnTp4Ep9GAtbWBJSTY5VzGjx+P5IQEPFdZaTIhaEp8PMaPH2/ynubJ1x4wfv6Kujqgqwu6zk6+B8gJ7XLEz9id7nutIoJVST0/N954I5YsWYKPPvoIZWVlKC8vx4cffoi77roLN910k5SXdIrc3FyMGjXK4N9ll13m7GbZhFarxY8//giVSoWCggLx+T5aW3uWa8uJnx+fSI4xoKnJ2a0xjkpbkG6se0Usd/683X6flEolnn3+eXzOGBYARhOCrn3hBZrsLIZW29NbR8vc3ZKknp/169eD4zjcfvvt6Op+g3h7e+Mvf/kLnn/+eZs20J5ycnKQk5Nj8JharUZoaCj8/PxEL6uWsgzbHsewJt+HcD0Dz59HUH09P3RjwXi3I85dv09sLF8/SacDLHgNS48jnLu/v79116Wigl/mHBPDz02ysl1Stxe7j83O3w7b23sfu517UBAfpLe0QNvYiILjx1FdXY3U1FSkp6eLCkgGOs4tCxbA/6WXsGL9eqTV1OgfT01MRN6mTRb/3nvitQeMnH9LC5/bzNsbGKBqgRzPRez2jrz2UvYRs71arbZ4W9HBT2dnJ2bPno0333wT69atw5kzZ8AYw7BhwxDgpER4pIdN8n0wBq61lb9xy/GahoTwwY+IN7pDCZ8Y6ZM2AYD4eKjefBMr1q9HSW2t/mGbJiAMCkLWgw8i8557sPPHH1FdVobUhASkz54NJa04FI8yO7s90cNe3t7eKCoqAsdxCAgIwJgxYzB27FgKfGSgb76PyQCC0JPvYx6AlcuWmR0CU7a392Q3lWNNG+GTWFubPJe8U3Zn0ovqm2+Q/cgjGFNbazAkNab7A4lKpbLNgXx8oAwJQXp6Om4cPRozYmOhtONEa7dGBU3dnqQ5P7fffjv++c9/2rotxEq2yvehFH7x5bTEvTelEkhNBUaN4uf/yA1ldybdbPWBxKS6OqM9oEz4vRB+l4k4SiW/ylVOK12JTUm6O3d0dOCtt97C119/jUmTJiGwzx/JjRs32qRxRBxb5ftQCl2+cu7Ni4hwdgtMo54f0k34QPIBTH8gSev+QDJjxgxxL97Zyeei0WqBiy4ynJsiVJen4EeamBh55xEjVpMU/BQVFWHChAkAgN9++83gOY6SujlN73wfk408b2m+Dy+59/zIHfX8kG52TUBYVsYHPoGBPfmvBL6+fPAjx2FhQmRA0t15165dtm4HsYH09HSkJCbiuYoK4/k+AKTGxCD9oov4OT3GAlWtFgrhhin34Ke+ns/2HBMjn16q3svcqefH49nqA0k/jY38+5/jgMGD+/8u07AXIQNyXqluYnNKpRIbNm/G5+Dze/TL98FxWL9iBZQtLabLLnAcNDEx0MXEyL/nor6en/Mgt2zPgwYBYWEOLXFA5En/gYTj0Dfzk5CAMDUxEenp6Za/qE7HD3cBQHS00cBfP+eHen7Ea2oCCguBM2ec3RJiR5Lvzg0NDdiwYQPuuusu3H333di4cSMa5fZHyANlZWUh7+OPcSQqCmkAQgCkAShKTOSXud95Jz9ZWCAU4AQ/OfOHPXvwbWEhCs6ckT4J01GEOQ5yWvKuVALJycDQoc5uCZEBsx9IGOM/kIj5Xauq4oMaHx8gPt74NsIqTer5Ea+tjQ8we90bifuhqu5uKGvuXJzevh07XnkFW/75T+zatQuniov5fCKRkYbDMWVlQEkJVHl5GJaSgvnz52Pjxo2YP38+hqWk2G4Zrj0IwU9LCz/cRIgMZWVlIS8vD0cSEgw/kCQkIG/jRmSlpwMnTgCWpPNvbweEJIaDB5vuXRQCo6QkG52FB6Fl7h6Bqrq7o9ZWKJVKpKelgQ0fbjpDZns7cP48VN9+i+xHHrEuMaIz+PjwN6i2Nr73Rw6Vl4X5PjTkRXrJyspCZmYmdu7caZjhWacDTp3iy8mcPMn3GA6QURi+vnzPbXMzEBpqejuFgq9GTsQTVrvKMccZsRlJwc8vv/xiEPgAPVXdJ02aZLPGEYmE4m7mclT4+kI7dChWXHutPg+J8CdbyEOygOOwctkyZGZmyrM2UEiIvIKfmhqgshKIiuI/mRPSTalU6uf26D+QKJXA8OH8/JKmJr4ndtQo03PyAP59Lof3urui7M4eQdLHU6Gqe1/uVNXdpQnd5xYk6CooLERJTY3ViRGdRm7zfoQcP3KfLE7kQ6nk8/RERQHDhhkEPlqtFgUFBfj4ww+R/+234ubhdXbyvxciKl17PMZ6JolT8OPWJN2hharu69evR1paGjiOww8//IC//vWvLlXVfSBtbW1obm62aFuNJWP1Vu4jZnvFhQtARwc0ERFgZvYrLi4GYD4PSXFxsdFePUec+4D7cBwUbW1gjIE1NhrU0xJ7nNbuPxKtra0GvZpi2sU1NoLTaMDa28FMvH/see2t2ccW52+PdjliH0ecu9l9IiL4gKU7gP7so4/w+JNP4mxlpX6T5Lg4PPvSS8jMzDR7DK6qClx1NdigQWBm5v548rUHep1/QwN8WloAhQK69naTE8blei5y/b2Xso8152IJj67qnpubi9zcXIPHZL/CyZyuLv0nF2bBJ5fY2FgA5vOQCNvJjkIB3ejR8iki2v2Hi1HPD7HCZx98gNvuuQfzAHyIXvPwqqtx22234b333hswAAIA5u0NDqAVXyKx0NCBhx2JW5B0h/bx8cHmzZtdvqp7Tk4OcnJyDB5Tq9UIDQ2Fn5+f6YnCJojdXso+Zrfv6uInTXZ2Qtc9BDnQPhkZGQMnRuzOQ5KRkTHgnB9HnLu99xECeX9/f+nXxdeX7zoPCwPMvIbc3l82OX87bW/vfRx57ub20Wq1eGL1aswDTM7D+9uqVVi4cOHAv5P+/sCFC/zqTjNt9ORrD/Q6//BwBPZOBeLkdkndR66/91L2EbO9WsT0B1FzfjQaDXJycpCQkIDo6GjcddddiIuLo6rucuLlxS9xTU62aHOzeUgArN+0SZ6TnfvSap2fm4PqehErFRQUoKSy0vp5eJTlmRCTRAU/Tz31FN5++23MnTsXCxcuxNdff42//OUv9mobcRCTeUiExIhyXObe1++/A7/+yi8BdhbGevIN0bAXkchm9cB8fPjhG8Z6gnIyMFef9kAsJuoOrVKp8M9//hMLFy4EANx666248sorodVqXaNnwBM0N/NDLyJ7HoQ8JDt27MCuXbswc+ZMzJkzx3Wuq0LB3+TV6v5FHh1Fp+NLW3R1UfBDJLNZPTCO4+8DHR187w/1RpqlOH6cHyIcPtyi1bLEdYnq+SkrKzOoQXP55ZfDy8sLlb1WIxAnYgz47Tfg8GFJXd1KpRJTp07FtGnTMHXqVNcJfICeJe/OLLFCpS2IDVhUDywpybJ6YFTmwnJaLbiuLr73RxgyJG5LVPCj1Wrh0+dN4eXlpZ8wRZysvZ0PgBQKz8tOKgQ/ra3UxU9cmk3n4UVHAykpzusNdSFK4e+Yl5d8Vo8SuxHVN88Yw+LFi+Hb6w9rW1sb7r33XgQGBuofk3U9KHcmIrmh2/Hy4qtbazT80FdkpOPbQKUtiI0I8/BWLF2KtPJy/eOpiYnI27TJ8nl4YWH2aaAbUlByQ48iKvhZtGhRv8duvfVWmzWGWElI8OSpK+9CQpwb/NTWAhUV/LwfC1fbEWKKyXpg1CthFwqhx5iCH48gKvjZunWrvdpBbMGTe34APviprnZeqQuh25z+OBEbMVoPTAydjl8E0dnpnA8ELkQpBD+eNmXAQ9GSFHfi6T0/QUF8wcegIOfk+6EcP0RutFq+ajzAl8+gzMUm0bCXZ6Hgx110dfX88fXUnh+OA4YMcd7xe0+YJEQOvL35OWg6Hb/knXo1TOoKDAQLCfHc+6eHobu0u+A4flVHRwdNuHUW6vkhcuTjA7S18atBKfgxqS0iArohQ+hn5CEo+DHBJau6+/ry/7rb7YgKv7KsbtzRATQ2QuPvLyoQtLa6sUKt5muqtbUNeFy5VsP25Mrecq1sLWWfvttzXV3gNBroLlww+b705GsPuMl9T+Ix5Prel2VVd3fhllXdCRSnT0Or0WBPRQWqm5qQnJyMtLQ0+6+SEYa9qOeHyEl3TwbX2QknV76Tr64ucHTv9ygWBz/Lly+3+EU3btwoqTGO5lZV3evq+D+6gYH9Pt05osKvnKobqwoLseLJJ1FSW6t/LCUxERs2bx4wP4pV1Y0DAoCkJH7oKyzMoomlcnt/eXJlb7lXtpayj3778HCgpYWfi2biNTz52gOA9uxZhBQXI2D4cASZKxviwHZJ3Ueuv/dS9rFXVXeLg5/CwkKD/x84cABarRbDhw8HAPz2229QKpWYOHGixQcnNsIYUFLCfx0zxqNTs6tUKmT/5S+Yxxg+AF8EsgjAcxUVyM7Otl+hVoWCcvsQeaISF2Zxwkovmu/jMSwOfnbt2qX/fuPGjQgODsY777yD8PBwAEB9fT3uuOMOy+rNENtqa+MDH6XSowMfrVaLFUuXYh5j2Iae2i2TAWxjDAs4DiuXLUNmZiYliiOeIzCQXwxBS7iN0mq1KNizB/sKC9GcmIg5qal0f/AAkpYFbdiwAevWrdMHPgAQHh6OtWvXYsOGDTZrHLGQpyc37FZQUICS8nI8hv5vbAWAVYyhuKwMBQUFtj+4TsfnVCFEbry9+QSHvUoQEZ5KpcKwlBTMu/9+rP33vzF/4UIMS0mhEk0eQFLwo1arUVNT0+/x2tpaNDU1Wd0oIpKnJzfsVlVVBYAf6jJmdJ/tbOr8eeDQIaC42PavTQixOZVKhezsbIwpLzcoHjume4icAiD3Jin4ue6663DHHXcgLy8P5eXlKC8vR15eHpYsWWKf+RRkYELw4+E9P3HdExWLTDxf1Gc7mxJy/FCCQyJHLS18gN7W5uyWyELfIfLJAILQM0Q+D8DKZcto9a8bkxT8vPHGG5g7dy5uvfVWJCcnIzk5GbfccgvmzJmD1157zdZtJObQsBcAID09HSmJiXiO46Dr85wOwDqOQ2pSkn3mpVF2ZyJnNTXA2bPOq3snM04dIieyICn4CQgIwGuvvYYLFy6gsLAQBw8eRF1dHV577TUE0riyY3V29vzh9fDgR6lUYsPmzfgcwAKOM+jKXsBx+BzA+k2b7DOZkbI7EzmjFV8GnDpETmTBqjoIgYGBGDt2LMaNG0dBj7N4eQHDhwOpqVTWAkBWVhby8vJwJCEBaQBCAKQBKEpMtN8yd4B6foi8CatAbRT8aLVaFBQU4OOPP0Z+fr7LDQ85dYicyILkO3VBQQHefPNNnDlzBnl5eUhISMB7772H1NRUTJ061ZZtJAPhOJOJyzxVVlYWMjMzsfOrr1B98iRSIyORftNNUNqzV4Z6foicCT0/Qj4bK6hUKqxYuhQl5eX6xyxJIion6enpSElIwHOVldjGmEEvgH6IPDGRUre4MUldBZ988glmz54Nf39/FBYWor3700RTUxOee+45mzaQECmUSiXSp03DjWlpmHHJJVCKqPkiCfX8EDmzUc+Pu6yQUra1YcPSpc4ZIieyICn4Wbt2Ld544w384x//gHevT7ppaWk4ePCgzRpHLFBVxa/icLFuZ0dhwcH8N/ac6MkYEBHBl7Wgnh8iR0LPj07XE6iLZM0KKVkNkzEGnD2LrFmzkJeb6/ghciILkj6mnjx5EtOmTev3eEhICBoaGqxtE7GUTscHP4wBISF8hmdigAUH86kA7Bn8cByVtiDyxnF8709HB9/7I6GHUlgh9QFMr5BK614hNWPGDP1zshsmq67m7wleXsi6+25k3nMPduzYgV27dmHmzJmYM2cO9fh4AEnBT1xcHE6fPo2UlBSDx3/44QcMGTLEFu1yura2NjQ3N1u0rUZYai6C2H2Mbq/RQNFdsFDX0dFvPF9Ku1q7h4daW1vhZcEN0hHnbtU+SiUUGg2g0UBXXz9gz4zYc7eqXXbcXuo+jjh/ub5f5HrtpexjcvuICH5RRFcX0OfeZsn5F3cn8DS3Qqq4uBiTJk0CNBp89sUXuG3JEpO19t577z1kZmaKP5cBDLhPaysUZ84AjEGXnKzPezRx4kRotVpMnDhR/7OQfAwn7iPX33sp+1hzLpaQNOz15z//GUuXLsVPP/0EjuNQWVmJf//731i5ciXuu+8+KS/pFLm5uRg1apTBv8suu8zZzbJc94VmVLPHNC8vMCEFgL2yj1NpC+IKQkOB4GDJPcSxsbEAzK+QErbTlZbi8UceGXCY7IlHHnHcEBhj4EpLAcbAQkL4YJB4LEk9Pw8//DAaGxsxc+ZMtLW1Ydq0afD19cXKlStx//3327qNdpOTk4OcnByDx9RqNUJDQ+Hn54cgkauoxG4vZR+D7evr+ZIWgwYNuOJLzDG6uucD+Pv7i9rPEecudZ/A2Fg+yZtON+DPSeq549w5BFVUAOHhgIieT7m9vySfv8jjSNne3vs48twdtY+tr31GRgafRLSiYsAVUhkZGVAqlcg/eRJnz53DhxhgmKyiAoWFhQbDZNaei8l9amv5r8HBwMiRBkWg3e2+J9ffeyn7iNleLWJ6g+TEMM8++yzOnz+Pn3/+Gfv27cO5c+fwzDPPSH05IgXV9LKMMB/KTnmQOOGTK80TIHLW2ckvjjh/XtLuYpOIVnX/vlmUSFDXNye7HSZJ+/jww94JCQaBD/FMkv4azJo1C2vWrEFAQAAmTZqEyy+/HEFBQaivr8esWbNs3UZiCpW1sExwMDBuHDB4sH1en3L8EFfQ3s6XuKiulvwSYpKIWpxIMDIS+PVXvm3d9zSh2vq1116LO++8EzNnzrS+2npYGHDJJUBUlPTXIG5DUvCTn5+PV199FQsWLEBLS4v+8Y6ODuzevdtmjSMD6Ozk55lwHAU/5nAc/89eKMcPcQVCb0dHB79CVKKsrCycLinBl19+iS1btmDXrl04VVzcb+WWxbX2LrmE7/k5fx44fhyqV16xbS6h3udKvbOkm+RxgG+++QbV1dWYPHkySkpKbNgkYhFvb2DsWODii+37h93dCL009nhN6vkhcubjw98rGLPu96ClBUqOQ3p6Oq6//nrMmDHD6NJwi4fJEhL4Ej0REdDqdFjx7LO2q7be2QkcPQrU1Uk/X+KWJAc/cXFx2L17N8aOHYvLLrsM+fn5NmwWsYi3N5W2sFRHB3DkCH8jtOJTrzH6OT/U80PkztoCp11dwIkT/DCVBQGIxcNkQUFAaioKGhtRUltru2rrpaX8udbU2Pz3nrg2SXdrrrunwdfXF//+97+xdu1aXHPNNXjkkUds2jhCbMbbm79Za7X8vAJbFuLt7OT/qFDPD5E7Hx8+t017Oz8XTiwhXYSvr8VDSPpaezt3orq6GqmpqUhPTzfaW1R17hwACyZJ90qY2JswSbq6uhqpgwYhPSYGSi8vICWFesiJAUnBD+sTQT/xxBMYOXIkFi1aZJNGEQucPcv3NERH0x9dS3Acv+qrvp7P9mzD4IeFh/esJCFEzqwtcCokRxTZ46xUKvVFQgdautx7kvRkI8/rJ0lrNEBJCT95uft32Wgm6ehobHj2WWRNmCCqvcT9SRr2Ki4uxqBBgwwe+9Of/oR9+/Zhy5YtNmkYGYAwOdCKVRseyU51vlhSEjB0KE2mJPJnbYFToedHSq+RBSyaJB0bi/RLLwUuXOCH4I4fh+qdd4xPkj53Dtn33OMyBVeJ40gKfpKTk6EwkjNl9OjR1PvjCEJ+Hy8v6m0QIySE/9rSYjSvCCFuLzycXySRmCh+366unnuPneYaWjRJOjcXyksuASIjAY6DtqkJKx5+2HaTpIlHsHjYa/ny5XjmmWcQGBiI5cuXD7jtxo0brW6YOc8++yy++OILHDp0CD4+Pp5VUJWSG0rj68v/a2/nP8GGhlr/mjodTaQkrkP4HZBCGPLy8+M/dEntPTJDmCS9YulSpPUawkpNTETepk09k6QDA4HERBT8978oqa0VXXCVeDaLg5/CwkJ0di+PLCwsNLkd56BJZR0dHbj++usxZcoU/POf/3TIMWWDkhtKFxzM37TVatsEP/X1UJSWAnFxwLBh1r8eIXIlcb6PFBZPkvbyQlV3EGZRJmlCulkc/Ozatcvo986yZs0aAMDbb7/t3IY4g9DzQ8GPeOHhfJmLsDCbvBxHCQ6Jq6mr4+8hUVHiyjxERvI9Pg7qcbb5JOnu7QgBJK72clXt7e1oN9NVKxRGa25utngoTdPdEyMUiLPHPr23V9bWAlottB0dwABtlNKupu4JjU0WVkB3xLnbfJ/gYH7Je5+fndhzB4C2+nooW1rQ2dICZqf3i6N+XlLOX67nInYfR5y7o/Yxt73it9/AtbVB29Wl7/20+Px9ffW/O3K59mPGjMHg+Hg8V1VlsuBqcnw8xowZY/Ke7i73Pbn+3kvZx5pzsYSoOT+WcsScHynWrVun7zEy59dff8WpU6fs3CLxuK4uBJeWgmMMjQ0NdstdcfDgQbu8riuw9Ny1Wi1+//FHNNbUIDgxEcMuv9xo7hJXQ9fevQVUVcG7pQWtZ8+io08PqKue/0233YYXX3gBmQAeAz/UVQTgOQBfMIaHb70VP/zwg9nXcdXztwV3OHchYLKEqDk/vR04cABarRbDhw8HAPz2229QKpWYOHGixQfva/Xq1WaDk/3792PSpEmSXn/VqlVmgzi1Wo2kpCSMGzcO8fHxFr2u8AMPENEdLHafftt3dppd6SWlXU1NTTh48CAmTJiAYAuWszri3G2+D2NAczM4rRas181fzLlv374dTzz6KEorK/WPDY6Px9rnn8f8+fNtei6O+nmJvfZSjuP0a2+CI87dUfuY256rqIDi3DnooqLAEhIAmD9/rrER0GrBgoP19x05Xfvp06fjkksuwROPPoq0Xr+TyQkJeGfdOrO/k+5y35Pr772UfaQco6amxuJtJc352bhxI4KDg/HOO+8gPDwcAFBfX4877rhDP0Yrxf3334+FCxcOuE1KSork1/f19YWvhSsdgoKCEGbhvBCv7vkeA41LW7uPI47RW3BwsEXn76h22XSfpiagtpa/iRt5P5k7d5VKhUWLFmEeY/gIvT5lVlVh0aJF/apbW3sucr32Uo7j9Gtvhj3P3VH7mN2+o4PP8uzn12/um8nzr6nhJzyHhur3kdu1v+2223DzzTdblEnaFFe/78n1917KPlKOYZeen942bNiAnTt36gMfAAgPD8fatWuRkZGBFStWSHlZDBo0qF/yREJsLjCQn/Tc2clP/BQxcVyr1WLF0qX6nCLC/AIhp8gCjsPKZcuQmZnpFkNgxA2JzfKs0/G5sQDZ1xK0dJI0IZKSHKrVaqPdS7W1taImHFmjtLQUhw4dQmlpKbRaLQ4dOoRDhw6hWViO6a5OnQJ+/116enrCBz7CjVHk+7WgoAAl5eW2K7xIiKOJLW7a0sIPFXt7S88RRIjMSOr5ue6663DHHXdgw4YNmDyZX1y4b98+/PWvfx2wu9+WnnzySbzzzjv6/48fPx4APzzntomstNqe0gyDBzu3La4uJIT/WarVfH00Cwm5QiinCHFZwvJ2rZbP2mwuTYOdS1oQ4gySgp833ngDK1euxK233qpPfOjl5YUlS5bgpZdesmkDTXn77bc9L8ePkN/H25vyylhLuJE3NfGfai1cNUc5RYjLUyiA4cMtv484MLkhIY4i6S9oQEAAXnvtNbz00ks4c+YMGGMYNmwYAm1YKZsYQWUtbCcggL/xd3Xx3foW3tj1hRcrKkzmFElNTLRq4j8hdmdpIMNYz3wf6vkhbkT0nJ/Ozk7MnDkTv/32GwIDAzF27FiMGzeOAh8H4Nra+G8os7NtCIVORcz7USqV2PDkk/i8e3Kz0cKLmzbRZGfiHjQafsKzlxe/OowQNyE6+PH29kZRUZHDaniRXqiml23FxgIjR/JfLdXYiKwJE5D34os4EheHNAAhANIAFCUmml3mTogsaDRAZSVw4cLA2wUGAmPHUt064nYkrfa6/fbbPa+YqBNptVoUFBTg4+3bkf/LL9DSigvb8Pfnh78sDeQ7O4GSEgBA1sKFOF1aii+//BJbtmzBrl27cKq4mAIf4ho0GqCqiq/zZY63Nx8EEeJGJM356ejowFtvvYWvv/4akyZN6jfkJdfyFq5IpVJhxdKlKCkv1z+WsnYtNmzeTH9oHYkxoLiYnyPk7w8kJECpUFBOEeKaxOb6IcTNSAp+ioqKMGHCBAB8WYveaDjMdlQqFbKzszGPMXyAXpmEKyqQnZ1NQyy20NICnDvH/zEYaDixpoafG6RQAEOG8F8JcVXCcveBcv1oNEBFBZ/VWUQ6CEJcgaTgp3epC2IflEnYQTo7+XkPfn5Ad50jo7y9+YBn8GCa+Elcn48PP9zLGP87YExTE58HS6Gg4Ie4HauSxRw7dgylpaXo6NV1ynGc2SJyrqCtrc3ibNFi6olYuo+QSfgDmM4knFZWhp07d5pcVi2lXa3dy+lbW1v1tVUGYo9zd+g+HAeFRgNoNGjtHroyeu6+vkByMv9Ho9f7whHn4qifl9hrL+U4srr2vTji3B21j6XbK7q6gPZ26C5cQCtjAAzPn6uuBqfRgIWHgxm5F7rLtQfc574n1997KftYcy6WkBT8/P7777juuutw5MgRcBwH1v2LIwx5abVaKS/rcLm5ucjNzTV4TC5tr66uBmA+k7CwHZFIqQQLCACn0YAzFuxqtYDQsyYMFRDiBpi3N7j2dn7oq+97mzFw3fl9GM1nI25IUvCzdOlSpKam4ptvvsGQIUPw888/48KFC1ixYgXWr19v6zbaTU5ODnJycgweU6vVCA0NhZ+fn+hJrFImvZraJzU1FYD5TMKpqalmjyumXV1dXQAAf39/UfvZ8twdvk9cHFBV1f/cz50DqquB1FSzSeEccS72PobUay/2OFK2t/c+jjx3R+1jdvuICD6Hj48Purrnu+nPX6PhezuVSiAqasAVke7w83K3+55cf++l7CNme7VQ/skCkmZt7t27F08//TSioqKgUCigUCgwdepUrFu3Dg8++KCUlyR96DMJcxx0fZ7TZxJOSqJMwrbQu9SFoLUVKCvjV8MIGW4JcSdxccDo0cbzXPUuaUGLWIgbkhT8aLVafTQ2aNAgVFZWAgCSk5Nx8uRJ27XOgymVSmzYvBmfA5RJ2N6CgqBlDD/88gsKvvsOP3z/PbSnTvGTQUNDgZgYZ7eQENvz8eF7d4wFN8IHARryIm5KUvAzevRoHD58GABwxRVX4MUXX8SPP/6Ip59+GkOGDLFpAz1ZVlYW8t59F0eioiiTsB2pPv0Uw+bPx9yHHsKGV17B/MxMDLv6aqh27wZSUpzdPEIcT6Hg/1E9L+KmJM35eeKJJ9DSPRSwdu1azJs3D+np6YiMjMRHH31k0wZ6uqz0dGRu346dx46hmjGkpqYiPT2denxsxGQupdpaZK9cibyUFAoyiXtijM/y3N4OhIUZPpeayj9PiJuSFPzMnj1b//2QIUNw7Ngx1NXVITw8nJIc2hJjQEMDlEol0jMygJAQyiRsQwPmUgKwAKBcSsR9cRyfvFOn48u8GHueEDdlszS1ERERFPjYQ1ISvyqDup9tTsil9BhM51IqLitDQUGBE1pHiAMYy/Qsk3QfhNiTxT0/y5cvt/hFqbaXjXAcH/hERBgk1iO2UVVVBcB8LiVhO0Lcjq8v0NYGrneNr6NH+fk+F13UUwOMEDdjcfBTWFho8P8DBw5Aq9Vi+PDhAPgaX0qlEhMnTrRtCwmxk7i4OADmcykJ2xHidvoWOG1v58tdcBxf0oUQN2Vx8NO7ntfGjRsRHByMd955B+Hh4QCA+vp63HHHHZR3xlY0Gr6uTng4ffqyE30upYoKbGPMYOhLn0spMZHe08R9dQ97CT0/+izngYFUvJe4NUnv7g0bNmDdunX6wAcAwsPDsXbtWmzYsMFmjXOmffv2ObfUxfnzfEVlGnKxG8qlRDye8MGqe86PUNKC5hgSdycp+FGr1aipqen3eG1tLZp6Z8l1YTfeeCOGpaRApVI5/uDdq7wA8D0/xG6ysrKQl5eHIwkJlEuJeJ6+w16U3JB4CEnBz3XXXYc77rgDeXl5KC8vR3l5OfLy8rBkyRK3+WPxDYAxFRXIzs52fADU0sKPuyuVQEiIY4/tgbKysnC6pATbt2/H8uXLsX37dpwqLnab9zIhJvn5AaNHQzd6NLjOTnDCfB8Kfoibk5Tn54033sDKlStx6623orOzk38hLy8sWbIEL730kk0b6CyXAdjGGBZwHFY8+CCuuuoqk8MfGo1G9OsPtA9XXg5OowGLiADr7oa29TFMaW1t1X/18jL/9nBUuxyxz8SJE6HVajFx4kT9z0EO7ZLrtZdyHLlee0ecu6P2kXTt29rg1dqKtvZ2KBgDs+A13OXaA+5z35Pr772Ufaw5F0tICn4CAgLw2muv4aWXXsKZM2fAGMOwYcMQGBgo5eWcJjc3F7m5uQaP9Z7nI+R6SauowJ49exw28ZXrHvJifbOuEkKIneh8fKATUmsQ4uYkBT+CwMBAjB071lZtcbicnBzk5OQYPKZWqxEaGqr/v5DrpaGhwWx2ZSnZl/vt09zMLzH18wPi4/tlWbXJMQbQ1dUFAPD39xe1n73b5Yh9pJ67vdvlqGM48vzp2svn/aI9fx67v/wS3x88CAQGYs6ECaIm+bvDz8vd7nty/b2Xso+Y7dVqtcXbSg5+vv32W3z77beora2FTqczeG7Lli1SX1Z2HJ7rpb2dX2IaFkbp5QkhdqVSqbDi/vtR0r2q9KXXXkNKYiI2bN5Mc96IW5M04XnNmjXIyMjAt99+i/Pnz6O+vt7gn7vQ53pJSnJcrpfISGDcOCAhwTHHI4R4JKGo75iqKoM0D05b6EGIA0me8Pz222/jtttus3V7ZONnAK9053rJc3SuF4WCEowRQuxmwKK+3Qs9qKgvcWeS/sJ2dHQgLS3N1m2RlasBFCUkODbXS/fYKyGE2BMV9SWeTlLwc9ddd+H999+3dVtk5eN163Dq558dO+59/Dhw7BjQ1ua4YxJCPA4V9SWeTtKwV1tbG/7v//4P33zzDcaOHQvvPgXw3KGq+5QxY6AUkTPAai0tfJZVhUJfb4cQQuyBivoSTycp+Dl8+DAuvfRSAEBRUZHBc5ybrFBiwcGOza4sTBQPC6P5PoQQu6KivsTTSQp+eld4d1e6lBRg0CDHHVAIfqiWFyHEzoSivtnZ2VjAcVjFGEaD7/FZ56yFHoQ4kFVJDo8dO4bS0lJ0CEXxwPf8zJ8/3+qGeZTeQ15Uy4sQ4gBCUd8VS5cirbxc/3hqYiLyNm2iPD/ErUkKfn7//Xdcd911OHLkCDiOA2MMQM+QV+8SES6to4MvMGrvsh1CBffQUBryIoQ4TFZWFjIzM7Fjxw7s2rULM2fOxJw5c6jHh7g9SX9ply5ditTUVNTU1CAgIABHjx7F999/j0mTJiE/P9/GTXSSpibgyBGguNj+x6IhL0KIkyiVSkydOhXTpk3D1KlTKfAhHkFSz8/evXvx3XffISoqCgqFAgqFAlOnTsW6devw4IMPorCw0NbtdLg2hYKvKqvRQFdfz9fbMsHqCrfh4eAaGsCUSr62lz2OYSF3qW4sZR93qm4s1+rOdO3l+X7x5GsPuM99T67XXso+sqzqrtVq9cXGBg0ahMrKSgwfPhzJyck4efKklJd0igGruiuVYP7+4Fpb+YDEnr0yoaFgvYqpEkIIIcR+JAU/o0ePxuHDhzFkyBBcccUVePHFF+Hj44P/+7//w5AhQ2zdRrsZqKq7n58fAv38gHPn+CcsqCwrtwq3UvZxt+rGYvZxp+rGUvahqu7yvfZS9qFr77n3Pbleeyn7yKqq+xNPPIGWlhYAwNq1azFv3jykp6cjMjISH330kZSXlKegID74MTEUZbXWVqCxEYiIAPz97XMMQgghhBiQFPzMnj1b//2QIUNw7Ngx1NXVITw83G2SHALo6e3RaACtFrDxRECuvp6fWN3WBgwdatPXJoQQQohxklZ7lZaW6pe3CyIiIsBxHEpLS23SMFnw8ekpNdHd02VLnLDEnVZ5EUIIIQ4jqecnNTUVVVVViI6ONnj8woULSE1NdZ88PwAQH8/n3gkIsNlLarVaFHzzDWoKC5GSmIj0MWNAi0sJIYQQx5AU/DDGjA5vNTc3w8/Pz+pGyUpkpE1fTqVSYcXSpSjplVE15emnsWHzZsqoSgghhDiAqOBn+fLlAPhMzn/7298Q0Ks3RKvV4qefftIXPCX9qVQqZGdnYx5j+ADQ19J5rqIC2dnZyMvLowCIEEIIsTNRwY+QvJAxhiNHjsBHmA8DwMfHB+PGjcPKlStt20I5aGnhJyaHhkpelaXVarFi6VLMYwzb0DPZajKAbYxhAcdh5bJlyMzMpAyrhBBCiB2JCn6Eau533HEHNm/ejBBPKcJZXd1Tf0ti8FNQUICS8nJ8gP6zzBUAVjGGtLIyFBQUYMaMGdLbSgghhJABiVrt9dNPP2HHjh3YunWrPvB59913kZqaiujoaNxzzz1ob2+3S0OdSljybkW+n6qqKgD8UJcxo/tsRwghhBD7EBX8rF69GocPH9b//8iRI1iyZAn+8Ic/4NFHH8X27duxbt06mzfS6WwQ/MTFxQHg5/gYU9RnO0IIIYTYh6jg59ChQ7jqqqv0///www9xxRVX4B//+AeWL1+OV155Bf/5z39s3kinCwjgl7trtXxWZgnS09ORkpiI5zgOuj7P6QCs4zikJiUhPT3d6uYSQgghxDRRwU99fT1iYmL0/9+9ezeuueYa/f8vu+wylJWV2a51csFxQGAg/73E3h+lUokNL72Ez7snN+8F0ARgL4AFHIfPAazftIkmOxNCCCF2JmrCc0xMDIqLi5GUlISOjg4cPHgQa9as0T/f1NQEb29vmzfSGdra2tDcK9DhFApwGg1YTQ1Yn0nPGo3Gotecfeml+Pczz2DVq68iraZG/3hKfDzee+EFZGRkGBxTyjGs3ae1u2ertbUVXl7m3x6Oapcj9hF77o5ql1yvvZTj0LWX5/vFk6894D73Pbleeyn7WHMulhAV/FxzzTV49NFH8cILL2Dbtm0ICAgwGKY5fPgwhrpQjarc3Fzk5uYaPGYqOzULDAQHgGttBTO6hRk6Hbi6OmROn45rFy3Cd7/8gpqaGiQnJyMtLY16fAghhBAHERX8rF27FllZWZg+fTqCgoLwzjvvGOT62bJlCzIyMmzeSHvJyclBTk6OwWNqtRqhoaHw8/NDkDDRGeDn/YSE8F9NFG812L6v8+cBX1/+X1ISru6u5zXgPmKPYYN9urq6AAD+/v6i9rN3uxyxj9Rzt3e7HHUMR54/XXt5vV88+doD7nffk+u1l7KPmO3VarXF24oKfqKiolBQUIDGxkYEBQX16634+OOPJf0wXIJC0TPvR4pz5/ivUVG2aQ8hhBBCJJFU2ys0NNTo4xEREVY1xm21tAAaDd9jZONaYYQQQggRR1Lw47E6OoCqKqCzExg2zPL9hOzQERGAhRPKCCGEEGIf9JdYDIWCn7sDAF1dlgcyCQl8XTAKfAghhBCnE5Xnx+N5eQF+fvz3YvP9BAX17EsIIYQQp6HgRyyxpS50ffM5E0IIIcSZKPgRKziY/2pJ8KNWA4cPAxUV9m0TIYQQQixGwY9YQs+PRmO+V+fcOb4eGPX+EEIIIbLhksFPSUkJlixZgtTUVPj7+2Po0KF46qmn0NHRYf+D+/jw/xjjl7Cb0tkJNDby3w8aZP92EUIIIcQiLrn86MSJE9DpdHjzzTcxbNgwFBUV4e6770ZLSwvWr19v/wYEBfE9PyZKYQDge30Y47ftUwuMEEIIIc7jksHPNddcY1BNfsiQITh58iRef/11xwQ/KSkmS1wA4IMeYUk8ZXQmhBBCZMUlh72MaWxsdFyG6YECH74x/LCXlxfQXcOLEEIIIfLgkj0/fZ05cwZ///vfsWHDhgG3a29vR3t7+4DbCIXRmpub0SBkZjaFMYAxaNraAPQUiFOcOQOuqQm66GgwYd5PHxqNxmAfc8RuL3WfpqYmg69yaZcj9hF77o5ql1yvvZTj0LWX5/vFk6894D73Pbleeyn7WHMuluAYY8zire1s9erVWLNmzYDb7N+/H5MmTdL/v7KyEtOnT8f06dPx1ltvWf36gvfffx8BAQEmn/erq4NvQwPaIiLQHhZm8JyisxM+ajXaQ0LAvL0tOh4hhBBCpNNoNLj55pvR2NiIkJCQAbeVVfBz/vx5nBfmypiQkpICv+5MyZWVlZg5cyauuOIKvP3221AoBh7Fs7TnJykpCUePHkV8fLzJ7biaGiiqqsBCQ9EcEwMAAwZLfQlRraX7iN1e6j5NTU04ePAgJkyYgGAhp5EM2uWIfcSeu6PaJddrL+U4dO3l+X7x5GsPuM99T67XXso+Uo5RU1ODESNGWBT8yGrYa9CgQRhk4bLwiooKzJw5ExMnTsTWrVvNBj4A4OvrC19fX4tePygoCGF9enQMeHvziQ4VCii6f8hBQg4gC3h11/mydB+x20vdRxAcHDzw+Tu4XY48f0vP3VHtkuu1l3IcuvbyfL8IPPnaA65/35PrtZeyj5RjCAGTRa9v8ZYyUllZiRkzZmDw4MFYv349zp07p38uNjbWMY0ICOALnXZ1AW1tfN2uCxf4Cu4xMT3JEAkhhBAiKy4Z/OzcuROnT5/G6dOnkZiYaPCcw0bxOA4IDASamvgeID8/PrdPSwv/OAU/hBBCiCy55FL3xYsXgzFm9J9DdQc4XEsLn/SwpYUPiiijMyGEECJbLtnzIxtBQdBqtSj4/ntUNTcjNSQE6VddBaUX/VgJIYQQuaK/0lZQ7dyJFQ88gJLqav1jKQkJ2PDKK8jKynJiywghhBBiiksOe8mBSqVC9g03YEx1NfYCaAKwF8CYykpkZ2dDpVI5uYWEEEIIMYaCHwm0Wi1WLF2KeYxhG4DJAIK6v25jDPMArFy2DNqBCp8SQgghxCko+JGgoKAAJeXleAz9f4AKAKsYQ3FZGQoKCpzQOkIIIYQMhIIfCaqqqgAAo008P7rPdoQQQgiRDwp+JIiLiwMAFJl4vqjPdoQQQgiRDwp+JEhPT0dKYiKe4zjo+jynA7CO45CalIT09HRnNI8QQgghA6Cl7ia0tbWhubnZ5PNrn38et912GxZwHFYxhtHge3zWcRw+B/DeunVobW01ub+YGiRStpe6j9Dm1tZWfW0VObTLEfuIPXdHtUuu117Kcejay/P94snXHnCf+55cr72Ufaw5F0t4dPCTm5uL3Nxcg8csXaGVmZmJ9957D48/8gjSKir0j6fEx+O9F15AZmamTdtKCCGEENvw6OAnJycHOTk5Bo+p1WqEhobCz8/PbDXZW265BQsXLsTOnTtRXV2N1NRUpKenQ6lUWtwGsdV3pVTrFbNPV1cXAMDf31/UfvZulyP2kXru9m6Xo47hyPOnay+v94snX3vA/e57cr32UvYRs71arbZ4W48OfmxBqVTq5/ZIeSMQQgghxLFowjMhhBBCPAoFP4QQQgjxKBT8EEIIIcSjUPBDCCGEEI9CwQ8hhBBCPAoFP4QQQgjxKBT8EEIIIcSjUPBDCCGEEI9CwQ8hhBBCPAoFP4QQQgjxKFTewgRzVd17k2uFW0+ubixlH3eqbizX6s507eX5fvHkaw+4z31Prtdeyj5U1d2OrKnqTgghhBDX5NHBj7VV3fuSW4VbKfu4W3VjMfu4U3VjKft4cmVvuV97KfvQtffc+55cr72UfexV1Z3m/BBCCCHEo1DwQwghhBCPQsEPIYQQQjwKBT+EEEII8SgU/BBCCCHEo3j0ai9jGGMAgKamJotnjgv5gHQ6ncXHEbuPI44B8LPlNRoN1Go1FArzsbGj2uWIfcSeu6PaJddrL+U4dO3l+X7x5GsPuM99T67XXso+Uo7R1NQEoOfv+EAo+OnjwoULAIDRo0c7uSWEEEIIEaupqQmhoaEDbkPBTx8REREAgNLSUrM/vN4uu+wy7N+/X9SxxO7jiGOo1WokJSWhrKwMISEhsmmXI/aRcu6OaJejjuGo86drL7/3iydfe8C97ntyvfZS9hG7PWMMEydORHx8vNltKfjpQ+j2Cw0NFfVGUCqVoraXso8jjiEICQmxeD9HtctR+4g5d0e1S67XXspx6NrL8/0CePa1B9zjvifXay9lHynH8PHxsWj4jiY8E0IIIcSjUPBDCCGEEI9CwQ8hhBBCPAoFPzbSt0CqPfZxxDGkcFS7PPn85XruUo5D116e7xcp6NrLbx+5Xnsp+9jzGByzZEG8BxGqujc2NkqaNObqPPn8PfncAc8+f08+d4DO35PP31PPnXp+CCGEEOJRKPghhBBCiEeh4KcPX19fPPXUU/D19XV2U5zCk8/fk88d8Ozz9+RzB+j8Pfn8PfXcac4PIYQQQjwK9fwQQgghxKNQ8EMIIYQQj0LBDyGEEEI8CgU/hBBCCPEoFPz08dprryE1NRV+fn6YOHEiCgoKnN0ku1u9ejU4jjP4Fxsb6+xm2c3333+P+fPnIz4+HhzHYdu2bQbPM8awevVqxMfHw9/fHzNmzMDRo0ed01g7MHf+ixcv7vd+mDx5snMaa2Pr1q3DZZddhuDgYERHR2PBggU4efKkwTbuev0tOXd3vvavv/46xo4dq69ePmXKFOzYsUP/vLted8D8ubvzdTeFgp9ePvroIyxbtgyPP/44CgsLkZ6ejjlz5qC0tNTZTbO7Sy65BFVVVfp/R44ccXaT7KalpQXjxo3Dq6++avT5F198ERs3bsSrr76K/fv3IzY2FldffTWampoc3FL7MHf+AHDNNdcYvB++/PJLB7bQfnbv3o2cnBzs27cPX3/9Nbq6upCRkYGWlhb9Nu56/S05d8B9r31iYiKef/55/PLLL/jll18wa9YsZGZm6gMcd73ugPlzB9z3upvEiN7ll1/O7r33XoPHRowYwR599FEntcgxnnrqKTZu3DhnN8MpALBPP/1U/3+dTsdiY2PZ888/r3+sra2NhYaGsjfeeMMJLbSvvufPGGOLFi1imZmZTmmPo9XW1jIAbPfu3Ywxz7r+fc+dMc+69owxFh4ezt566y2Puu4C4dwZ87zrzhhj1PPTraOjAwcOHEBGRobB4xkZGdizZ4+TWuU4p06dQnx8PFJTU7Fw4UL8/vvvzm6SUxQXF6O6utrgfeDr64vp06d7xPtAkJ+fj+joaFx88cW4++67UVtb6+wm2UVjYyMAICIiAoBnXf++5y7whGuv1Wrx4YcfoqWlBVOmTPGo69733AWecN1783J2A+Ti/Pnz0Gq1iImJMXg8JiYG1dXVTmqVY1xxxRV49913cfHFF6OmpgZr165FWloajh49isjISGc3z6GEa23sfXD27FlnNMnh5syZg+uvvx7JyckoLi7G3/72N8yaNQsHDhxwqyywjDEsX74cU6dOxejRowF4zvU3du6A+1/7I0eOYMqUKWhra0NQUBA+/fRTjBo1Sh/guPN1N3XugPtfd2Mo+OmD4ziD/zPG+j3mbubMmaP/fsyYMZgyZQqGDh2Kd955B8uXL3diy5zHE98HghtvvFH//ejRozFp0iQkJyfjiy++QFZWlhNbZlv3338/Dh8+jB9++KHfc+5+/U2du7tf++HDh+PQoUNoaGjAJ598gkWLFmH37t365935ups691GjRrn9dTeGhr26DRo0CEqlsl8vT21tbb9PA+4uMDAQY8aMwalTp5zdFIcTVrnR+6BHXFwckpOT3er98MADD+C///0vdu3ahcTERP3jnnD9TZ27Me527X18fDBs2DBMmjQJ69atw7hx47B582aPuO6mzt0Yd7vuxlDw083HxwcTJ07E119/bfD4119/jbS0NCe1yjna29tx/PhxxMXFObspDpeamorY2FiD90FHRwd2797tce8DwYULF1BWVuYW7wfGGO6//36oVCp89913SE1NNXjena+/uXM3xp2uvTGMMbS3t7v1dTdFOHdj3P26A6DVXr19+OGHzNvbm/3zn/9kx44dY8uWLWOBgYGspKTE2U2zqxUrVrD8/Hz2+++/s3379rF58+ax4OBgtz3vpqYmVlhYyAoLCxkAtnHjRlZYWMjOnj3LGGPs+eefZ6GhoUylUrEjR46wm266icXFxTG1Wu3kltvGQOff1NTEVqxYwfbs2cOKi4vZrl272JQpU1hCQoJbnP9f/vIXFhoayvLz81lVVZX+n0aj0W/jrtff3Lm7+7VftWoV+/7771lxcTE7fPgwe+yxx5hCoWA7d+5kjLnvdWds4HN39+tuCgU/feTm5rLk5GTm4+PDJkyYYLAM1F3deOONLC4ujnl7e7P4+HiWlZXFjh496uxm2c2uXbsYgH7/Fi1axBjjlzs/9dRTLDY2lvn6+rJp06axI0eOOLfRNjTQ+Ws0GpaRkcGioqKYt7c3Gzx4MFu0aBErLS11drNtwth5A2Bbt27Vb+Ou19/cubv7tb/zzjv19/aoqCh21VVX6QMfxtz3ujM28Lm7+3U3hWOMMcf1MxFCCCGEOBfN+SGEEEKIR6HghxBCCCEehYIfQgghhHgUCn4IIYQQ4lEo+CGEEEKIR6HghxBCCCEehYIfQgghhHgUCn4IIR6jpKQEHMfh0KFDovf97rvvMGLECOh0OsnHb29vx+DBg3HgwAHJr0EIsR4FP4QQh1i8eDE4jgPHcfD29kZMTAyuvvpqbNmyxaqAYqDjLViwwGav9/DDD+Pxxx+HQiH9tunr64uVK1fikUcesVm7CCHiUfBDCHGYa665BlVVVSgpKcGOHTswc+ZMLF26FPPmzUNXV5ezm2fSnj17cOrUKVx//fVWv9Ytt9yCgoICHD9+3AYtI4RIQcEPIcRhfH19ERsbi4SEBEyYMAGPPfYYPvvsM+zYsQNvv/22frvGxkbcc889iI6ORkhICGbNmoVff/1V//zq1atx6aWX4s0330RSUhICAgJw/fXXo6GhQf/8O++8g88++0zf25Sfn6/f//fff8fMmTMREBCAcePGYe/evQO2+8MPP0RGRgb8/Pz6tWHLli0YPHgwgoKC8Je//AVarRYvvvgiYmNjER0djWeffdbgtSIjI5GWloYPPvhA+g+SEGIVCn4IIU41a9YsjBs3DiqVCgDAGMPcuXNRXV2NL7/8EgcOHMCECRNw1VVXoa6uTr/f6dOn8Z///Afbt2/H//73Pxw6dAg5OTkAgJUrV+KGG27Q9zRVVVUhLS1Nv+/jjz+OlStX4tChQ7j44otx0003Ddjz9P3332PSpEn9Hj9z5gx27NiB//3vf/jggw+wZcsWzJ07F+Xl5di9ezdeeOEFPPHEE9i3b5/BfpdffjkKCgqs+rkRQqTzcnYDCCFkxIgROHz4MABg165dOHLkCGpra+Hr6wsAWL9+PbZt24a8vDzcc889AIC2tja88847SExMBAD8/e9/x9y5c7FhwwbExsbC398f7e3tiI2N7Xe8lStXYu7cuQCANWvW4JJLLsHp06cxYsQIo+0rKSlBfHx8v8d1Oh22bNmC4OBgjBo1CjNnzsTJkyfx5ZdfQqFQYPjw4XjhhReQn5+PyZMn6/dLSEhASUmJ9B8YIcQqFPwQQpyOMQaO4wAABw4cQHNzMyIjIw22aW1txZkzZ/T/Hzx4sD7wAYApU6ZAp9Ph5MmTRgOe3saOHav/Pi4uDgBQW1trMvhpbW01GPISpKSkIDg4WP//mJgYKJVKg0nRMTExqK2tNdjP398fGo1mwDYSQuyHgh9CiNMdP34cqampAPjelLi4OIM5OoKwsDCTryEET8LXgXh7e/fbb6AVZ4MGDUJ9ff2AryO8lrHH+r52XV0doqKizLaTEGIfFPwQQpzqu+++w5EjR/DQQw8BACZMmIDq6mp4eXkhJSXF5H6lpaWorKzUD0ft3bsXCoUCF198MQDAx8cHWq3WJm0cP348jh07ZpPXAoCioiKMHz/eZq9HCBGHJjwTQhymvb0d1dXVqKiowMGDB/Hcc88hMzMT8+bNw+233w4A+MMf/oApU6ZgwYIF+Oqrr1BSUoI9e/bgiSeewC+//KJ/LT8/PyxatAi//vorCgoK8OCDD+KGG27QD3mlpKTg8OHDOHnyJM6fP4/Ozk7J7Z49ezZ++OEH606+l4KCAmRkZNjs9Qgh4lDwQwhxmP/973+Ii4tDSkoKrrnmGuzatQuvvPIKPvvsMyiVSgD8MNGXX36JadOm4c4778TFF1+MhQsXoqSkBDExMfrXGjZsGLKysnDttdciIyMDo0ePxmuvvaZ//u6778bw4cMxadIkREVF4ccff5Tc7ltvvRXHjh3DyZMnpZ98t71796KxsRHZ2dlWvxYhRBqOMcac3QhCCBFj9erV2LZtm6QyFVI9/PDDaGxsxJtvvmnV61x//fUYP348HnvsMRu1jBAiFvX8EEKIBR5//HEkJydbNY+ovb0d48aN089vIoQ4B/X8EEJcjjN6fggh7oOCH0IIIYR4FBr2IoQQQohHoeCHEEIIIR6Fgh9CCCGEeBQKfgghhBDiUSj4IYQQQohHoeCHEEIIIR6Fgh9CCCGEeBQKfgghhBDiUSj4IYQQQohH+X/rI4nruHe+WwAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.plot(npor,color='red',ls='--',alpha=0.2,zorder=1)\n", "plt.scatter(np.arange(0,len(npor),1),npor,color='red',edgecolor='black',zorder=1)\n", "plt.xlabel('Depth (m)'); plt.ylabel(r'Standardized Porosity, $\\overline{x} = 0.0$, $s_x = 1.0$'); plt.title('Porosity for a Single Vertical Well')\n", "plt.xlim([0,39]); plt.ylim([-2.5,2.5]); add_grid()" ] }, { "cell_type": "markdown", "id": "d61e9599", "metadata": {}, "source": [ "Notice that we ensured that the dataset variance is 1.0 as we assume this to calculate the correlogram below.\n", "\n", "#### Interactive Interface\n", "\n", "Here's the interactive interface. I calculate the variogram, plot the h-scatterplot and calculate and annotate the correlogram / h-scatterplot correlation coefficient. \n", "\n", "* the user specifies lag to investigate" ] }, { "cell_type": "code", "execution_count": 6, "id": "ba7230a6", "metadata": {}, "outputs": [], "source": [ "l = widgets.Text(value=' Variogram h-Scatterplot Demonstration, Prof. Michael Pyrcz, The University of Texas at Austin',\n", " layout=Layout(width='930px', height='30px'))\n", "\n", "lag = widgets.IntSlider(min=1,max = 40,value=5,step = 1,description = 'Lag',orientation='horizontal',style = {'description_width': 'initial'},layout=Layout(width='940px',height='30px'),continuous_update=False)\n", "\n", "ui2 = widgets.VBox([l,lag],)\n", "\n", "def run_plot(lag): \n", " size = 0.25\n", " gamma = np.average(np.square((npor - npor.shift(lag)).dropna()))*0.5\n", " gamma_all = []; num_pairs_all = []\n", " for ilag in range(0,40):\n", " num_pairs_all.append(len((npor - npor.shift(ilag)).dropna()))\n", " gamma_all.append(np.average(np.square((npor - npor.shift(ilag)).dropna()))*0.5) \n", " \n", " npor_shift = npor.shift(lag)\n", " correl = np.round(np.corrcoef(npor[~np.isnan(npor_shift)],npor_shift[~np.isnan(npor_shift)]),2)[0][1]\n", " \n", " plt.subplot(131)\n", " scatter = plt.scatter(np.arange(0,40)*size,gamma_all,s=num_pairs_all,color='red',edgecolor='black')\n", " plt.scatter(lag*size,gamma,color='darkorange',edgecolor='black',s=40,zorder=10)\n", " plt.plot([lag*size,lag*size],[0,gamma],color='black',ls='--',zorder=1)\n", " plt.plot([lag*size,lag*size],[1.0,gamma],color='red',ls='--',zorder=1)\n", " if gamma < 1.0:\n", " plt.plot([lag*size,lag*size+0.1],[gamma,gamma+0.1],color='red',zorder=1)\n", " plt.plot([lag*size,lag*size-0.1],[gamma,gamma+0.1],color='red',zorder=1)\n", " else:\n", " plt.plot([lag*size,lag*size+0.1],[gamma,gamma-0.1],color='red',zorder=1)\n", " plt.plot([lag*size,lag*size-0.1],[gamma,gamma-0.1],color='red',zorder=1)\n", " plt.plot([0,lag*size],[gamma,gamma],color='black',ls='--',zorder=1)\n", " plt.plot([0,10],[1.0,1.0],color='black')\n", " plt.xlim([0,10]); plt.ylim([0,3]); plt.title('Experimental Variogram')\n", " plt.annotate(r'$\\gamma(\\bf{h}) =$ ' + str(np.round(gamma,2)),[lag*size+0.2,(gamma)/2])\n", " plt.annotate(r'$\\bf{h} =$ ' + str(lag*size),[lag*size-0.3,0.03],rotation=90)\n", " plt.annotate(r'$\\sigma^2 - \\gamma(\\bf{h}) =$ ' + str(np.round(1.0-gamma,2)),\n", " [lag*size+0.2,(gamma+1.0)/2],color='red')\n", " plt.xlabel(r'$\\bf{h}$ Lag Distance'); plt.ylabel(r'$\\gamma(\\bf{h})$ Variogram')\n", " add_grid()\n", " legend = plt.legend(*scatter.legend_elements(\"sizes\", num=4),loc='upper left')\n", " legend.set_title('Number of Pairs')\n", " \n", " plt.subplot(132)\n", " #plt.plot(npor,color='red',ls='--',alpha=0.2,zorder=1)\n", " n = 0\n", " for itail in range(0,len(npor)):\n", " if itail + lag < len(npor)-1:\n", " plt.plot([itail,itail+lag],[npor[itail],npor[itail+lag]],color='gray',alpha=0.4,zorder=1) \n", " n = n + 1\n", " itail = 0.0 \n", " while itail < len(npor)-1:\n", " if itail + lag < len(npor)-1:\n", " plt.plot([itail,itail+lag],[npor[itail],npor[itail+lag]],color='black',lw=2,zorder=2)\n", " itail = itail + lag + 1\n", " plt.scatter(np.arange(0,len(npor),1),npor,color='red',edgecolor='black',zorder=1)\n", " plt.xlabel('Depth (m)'); plt.ylabel(r'Standardized Porosity, $\\overline{x} = 0.0$, $s_x = 1.0$'); plt.title('Standardized Porosity for a Single Vertical Well and Pairs for Lag = ' + str(lag))\n", " plt.annotate('Number of Pairs = ' + str(n),[25,-2.3])\n", " plt.xlim([0,39]); plt.ylim([-2.5,2.5]); add_grid()\n", " \n", " plt.subplot(133)\n", " plt.scatter(npor,npor.shift(lag),color='darkorange',edgecolor='black',s=20)\n", " plt.plot([-3,3],[-3,3],color='black')\n", " plt.xlim([-3,3]); plt.ylim([-3,3]); \n", " plt.title(r'h-Scatter Plot, lag = ' + str(lag) + r', $\\bf{h} =$ ' + str(lag*size))\n", " plt.xlabel(r'$Z(\\bf{u})$ Tail'); plt.ylabel(r'$Z(\\bf{u}+ \\bf{h})$ Head')\n", " plt.annotate(r'$\\rho_{Z(\\bf{u}),Z(\\bf{u} + \\bf{h})}$ = ' + str(correl),[1.0,-2.5],fontsize=12)\n", " add_grid()\n", " plt.subplots_adjust(left=0.0, bottom=0.0, right=2.6, top=0.9, wspace=0.2, hspace=0.3); plt.show()\n", " \n", "# connect the function to make the samples and plot to the widgets \n", "interactive_plot = widgets.interactive_output(run_plot, {'lag':lag})\n", "interactive_plot.clear_output(wait = True) # reduce flickering by delaying plot updating \n", " " ] }, { "cell_type": "markdown", "id": "7df92474", "metadata": {}, "source": [ "### Interactive Variogram h-scatterplot Demonstration \n", "\n", "#### Michael Pyrcz, Professor, The University of Texas at Austin \n", "\n", "Change the number of sample data, train/test split and the data noise and observe overfit! Change the model order to observe a specific model example.\n", "\n", "### The Inputs\n", "\n", "* **lag** - the lag number to calculate, h = lag $\\times$ data spacing" ] }, { "cell_type": "code", "execution_count": 7, "id": "6aa493e8", "metadata": { "scrolled": false }, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "82ee4ad9b4054e17b7e9a3eef1134b08", "version_major": 2, "version_minor": 0 }, "text/plain": [ "VBox(children=(Text(value=' Variogram h-Scatterplot Demonstration, Prof. Michae…" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "40d66cab08a94d3f93a6beb91ff2e120", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Output()" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "display(ui2, interactive_plot) # display the interactive plot" ] }, { "cell_type": "markdown", "id": "5267d3ae", "metadata": {}, "source": [ "#### Comments\n", "\n", "This was an interactive demonstration of the variogram h-scatter plot. \n", "\n", "I have many other demonstrations on simulation to build spatial models with spatial continuity and many other workflows available [here](https://github.com/GeostatsGuy/PythonNumericalDemos), along with a package for geostatistics in Python called [GeostatsPy](https://github.com/GeostatsGuy/GeostatsPy). \n", " \n", "We hope this was helpful,\n", "\n", "*Michael*\n", "\n", "***\n", "\n", "#### More on Michael Pyrcz and the Texas Center for Geostatistics:\n", "\n", "### Michael Pyrcz, Professor, The University of Texas at Austin \n", "*Novel Data Analytics, Geostatistics and Machine Learning Subsurface Solutions*\n", "\n", "With over 17 years of experience in subsurface consulting, research and development, Michael has returned to academia driven by his passion for teaching and enthusiasm for enhancing engineers' and geoscientists' impact in subsurface resource development. \n", "\n", "For more about Michael check out these links:\n", "\n", "#### [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig) | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)\n", "\n", "#### Want to Work Together?\n", "\n", "I hope this content is helpful to those that want to learn more about subsurface modeling, data analytics and machine learning. Students and working professionals are welcome to participate.\n", "\n", "* Want to invite me to visit your company for training, mentoring, project review, workflow design and / or consulting? I'd be happy to drop by and work with you! \n", "\n", "* Interested in partnering, supporting my graduate student research or my Subsurface Data Analytics and Machine Learning consortium (co-PIs including Profs. Foster, Torres-Verdin and van Oort)? My research combines data analytics, stochastic modeling and machine learning theory with practice to develop novel methods and workflows to add value. We are solving challenging subsurface problems!\n", "\n", "* I can be reached at mpyrcz@austin.utexas.edu.\n", "\n", "I'm always happy to discuss,\n", "\n", "*Michael*\n", "\n", "Michael Pyrcz, Ph.D., P.Eng. Associate Professor The Hildebrand Department of Petroleum and Geosystems Engineering, Bureau of Economic Geology, The Jackson School of Geosciences, The University of Texas at Austin\n", "\n", "#### More Resources Available at: [Twitter](https://twitter.com/geostatsguy) | [GitHub](https://github.com/GeostatsGuy) | [Website](http://michaelpyrcz.com) | [GoogleScholar](https://scholar.google.com/citations?user=QVZ20eQAAAAJ&hl=en&oi=ao) | [Book](https://www.amazon.com/Geostatistical-Reservoir-Modeling-Michael-Pyrcz/dp/0199731446) | [YouTube](https://www.youtube.com/channel/UCLqEr-xV-ceHdXXXrTId5ig) | [LinkedIn](https://www.linkedin.com/in/michael-pyrcz-61a648a1)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "032a7648", "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.11.4" } }, "nbformat": 4, "nbformat_minor": 5 }