{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<p align=\"center\">\n",
    "    <img src=\"https://github.com/GeostatsGuy/GeostatsPy/blob/master/TCG_color_logo.png?raw=true\" width=\"220\" height=\"240\" />\n",
    "\n",
    "</p>\n",
    "\n",
    "### Boostrap vs. Spatial Bootstrap in Python \n",
    "\n",
    "#### Michael Pyrcz, Professor, The University of Texas at Austin \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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Bootstrap\n",
    "\n",
    "Uncertainty in the sample statistics\n",
    "* one source of uncertainty is the paucity of data.\n",
    "* do 200 or even less wells provide a precise (and accurate estimate) of the mean? standard deviation? skew? P13?\n",
    "\n",
    "Would it be useful to know the uncertainty in these statistics due to limited sampling?\n",
    "* what is the impact of uncertainty in the mean porosity e.g. 20%+/-2%?\n",
    "\n",
    "**Bootstrap** is a method to assess the uncertainty in a sample statistic by repeated random sampling with replacement.\n",
    "\n",
    "Assumptions\n",
    "* sufficient, representative sampling, identical, idependent samples\n",
    "\n",
    "Limitations\n",
    "1. assumes the samples are representative \n",
    "2. assumes stationarity\n",
    "3. only accounts for uncertainty due to too few samples, e.g. no uncertainty due to changes away from data\n",
    "4. does not account for boundary of area of interest \n",
    "5. assumes the samples are independent\n",
    "6. does not account for other local information sources\n",
    "\n",
    "The Bootstrap Approach (Efron, 1982)\n",
    "\n",
    "Statistical resampling procedure to calculate uncertainty in a calculated statistic from the data itself.\n",
    "* Does this work?  Prove it to yourself, for uncertainty in the mean solution is standard error: \n",
    "\n",
    "\\begin{equation}\n",
    "\\sigma^2_\\overline{x} = \\frac{\\sigma^2_s}{n}\n",
    "\\end{equation}\n",
    "\n",
    "Extremely powerful - could calculate uncertainty in any statistic!  e.g. P13, skew etc.\n",
    "* Would not be possible access general uncertainty in any statistic without bootstrap.\n",
    "* Advanced forms account for spatial information and sampling strategy (game theory and Journel’s spatial bootstrap (1993).\n",
    "\n",
    "Steps: \n",
    "\n",
    "1. assemble a sample set, must be representative, reasonable to assume independence between samples\n",
    "\n",
    "2. optional: build a cumulative distribution function (CDF)\n",
    "    * may account for declustering weights, tail extrapolation\n",
    "    * could use analogous data to support\n",
    "\n",
    "3. For $\\ell = 1, \\ldots, L$ realizations, do the following:\n",
    "\n",
    "    * For $i = \\alpha, \\ldots, n$ data, do the following:\n",
    "\n",
    "        * Draw a random sample with replacement from the sample set or Monte Carlo simulate from the CDF (if available). \n",
    "\n",
    "6. Calculate a realization of the sammary statistic of interest from the $n$ samples, e.g. $m^\\ell$, $\\sigma^2_{\\ell}$. Return to 3 for another realization.\n",
    "\n",
    "7. Compile and summarize the $L$ realizations of the statistic of interest.\n",
    "\n",
    "This is a very powerful method.  Let's try it out.\n",
    "\n",
    "#### Load the Required Libraries\n",
    "\n",
    "We will also need some standard Python packages. These should have been installed with Anaconda 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "supress_warnings = True                                   # supress warnings?\n",
    "import numpy as np                                        # ndarrys for gridded data\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt                           # for plotting\n",
    "import math\n",
    "from matplotlib.ticker import (MultipleLocator, AutoMinorLocator) # control of axes ticks\n",
    "import scipy as sp                                        # lower upper matrix decomposition\n",
    "from scipy.interpolate import make_interp_spline          # smooth curves\n",
    "import statsmodels.formula.api as smf\n",
    "from scipy.stats import norm                              # Gaussian distribution PDF and sampling\n",
    "import geostatspy.GSLIB as GSLIB\n",
    "import geostatspy.geostats as geostats\n",
    "from ipywidgets import interactive                        # 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",
   "metadata": {},
   "source": [
    "If you get a package import error, you may have to first install some of these packages. This can usually be accomplished by opening up a command window on Windows and then typing 'python -m pip install [package-name]'. More assistance is available with the respective package docs.  \n",
    "\n",
    "#### Declare Functions\n",
    "\n",
    "I just added a convenience functions to:\n",
    "\n",
    "1. add major and minor gridlines\n",
    "2. calculate a isotropic variogram between 2 points\n",
    "3. calculate the azimuth between 2 points to rotate a label"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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    \n",
    "\n",
    "\n",
    "    \n",
    "    \n",
    "def n_effective(df,xcol,ycol,seed,nreal,vario):\n",
    "    \"\"\"Calculate the number of effective data from spatial locations and spatial continuity model\n",
    "    Used in bootstrap to account for spatial continuity, use n effective instead of number of data\n",
    "    :param df: source DataFrame\n",
    "    :param xcol: column with the X locations\n",
    "    :param ycol: column with the Y locations\n",
    "    :param seed: random number seed for the random sampling\n",
    "    :param nreal: number of realizations to sample the variance of the average \n",
    "    :param vario: variogram model as a dictionary, see the GeostatsPy Package's GSLIB.make_variogram() function \n",
    "    :return: n_eff as effective number of data\n",
    "    \"\"\" \n",
    "\n",
    "# Set constants\n",
    "    np.random.seed(seed)\n",
    "    PMX = 9999.9\n",
    "    \n",
    "# load the variogram\n",
    "    nst = vario['nst']\n",
    "    cc = np.zeros(nst); aa = np.zeros(nst); it = np.zeros(nst)\n",
    "    ang = np.zeros(nst); anis = np.zeros(nst)\n",
    "    c0 = vario['nug']; \n",
    "    cc[0] = vario['cc1']; it[0] = vario['it1']; ang[0] = vario['azi1']; \n",
    "    aa[0] = vario['hmaj1']; anis[0] = vario['hmin1']/vario['hmaj1'];\n",
    "    if nst == 2:                                   # include 2nd structure if present (optional)\n",
    "        cc[1] = vario['cc2']; it[1] = vario['it2']; ang[1] = vario['azi2']; \n",
    "        aa[1] = vario['hmaj2']; anis[1] = vario['hmin2']/vario['hmaj2'];\n",
    "    \n",
    "# Set up the rotation matrix\n",
    "    rotmat, maxcov = geostats.setup_rotmat(c0,nst,it,cc,ang,PMX)\n",
    "    \n",
    "# Load the data\n",
    "    nd = len(df)\n",
    "    x = df[xcol].values\n",
    "    y = df[ycol].values\n",
    "    \n",
    "# Calculate Symmetric Covariance Array - assuming variogram with spherical structure with range specified\n",
    "    cov = np.zeros((nd,nd))\n",
    "    var_range = 100.0\n",
    "    for i in range(0, nd):\n",
    "        x1 = x[i]; y1 = y[i]\n",
    "        for j in range(0, nd):\n",
    "            x2 = x[j]; y2 = y[j]\n",
    "            cova = geostats.cova2(x1, y1, x2, y2, nst, c0, PMX, cc, aa, it, ang, anis, rotmat, maxcov)\n",
    "            cov[i,j] = cova\n",
    "            \n",
    "# Lower and upper deconvolution            \n",
    "    P, L, U = sp.linalg.lu(cov) \n",
    "    \n",
    "# Build realization and calculate the average    \n",
    "    realizations = np.zeros((nreal,nd))\n",
    "    average_array = np.zeros(nreal)\n",
    "    rand = np.zeros((nd)) \n",
    "    for l in range(0, nreal):\n",
    "        rand = np.random.normal(loc = 0.0, scale = 1.0, size = nd)\n",
    "        realizations[l,:] = np.matmul(L,rand)\n",
    "        average_array[l] = np.average(realizations[l])\n",
    "        \n",
    "# Back out the number of effecitve data useing the standard error in the average\n",
    "    var_average = np.var(average_array)\n",
    "    n_eff = max(min(1.0/var_average, nd),1.0)    # filter n effective less than 1.0 or greater than number of data\n",
    "    return n_eff, realizations\n",
    "    \n",
    "def percentile(n):\n",
    "    def percentile_(x):\n",
    "        return np.percentile(x, n)\n",
    "    percentile_.__name__ = 'percentile_%s' % n\n",
    "    return percentile_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### n Effective Demonstration\n",
    "\n",
    "Now let's set up our first dashboard. Here we visualize the impact of spatial continuity range and number of samples on the effective number of data. \n",
    "\n",
    "* How many independent peices of information do you have?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "max_L = 100\n",
    "\n",
    "l = widgets.Text(value=r'                                                      Number of Effective Data, Michael Pyrcz, Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "\n",
    "n = widgets.IntSlider(min=1, max = max_L, value = 1,step = 1,description = '$n$:',orientation='horizontal',layout=Layout(width='300px', height='50px'),continuous_update=False)\n",
    "n.style.handle_color = 'gray'\n",
    "\n",
    "vrange = widgets.FloatSlider(min=0.1, max = 200.0, value = 10.0,step=5,description = r'$\\gamma_{range}$',orientation='horizontal',layout=Layout(width='300px', height='50px'),continuous_update=False)\n",
    "vrange.style.handle_color = 'gray'\n",
    "\n",
    "seed = widgets.IntSlider(min=100, max = 999, value = 1, description = '$s$:',orientation='horizontal',layout=Layout(width='300px', height='50px'),continuous_update=False)\n",
    "seed.style.handle_color = 'gray'\n",
    "\n",
    "window = widgets.Checkbox(value=False,description='Conv. Fit',disabled=False,layout=Layout(width='200px', height='30px'))\n",
    "poly = widgets.Checkbox(value=True,description='Poly. Fit',disabled=False,layout=Layout(width='200px', height='30px'))\n",
    "\n",
    "ui1 = widgets.HBox([n,vrange,seed,],kwargs = {'justify_content':'center'})\n",
    "ui2 = widgets.HBox([window,poly,],kwargs = {'justify_content':'center'})\n",
    "ui_all = widgets.VBox([l,ui1,ui2],)\n",
    "\n",
    "def f_make1(n,vrange,seed,window,poly): \n",
    "    np.random.seed(seed=seed)\n",
    "    values = np.random.rand(max_L*2)*100\n",
    "    X,Y = np.split(values,2)\n",
    "    df = pd.DataFrame(np.vstack((X[:n],Y[:n])).T,columns = ['X','Y'])\n",
    "    \n",
    "    ax1 = plt.subplot(121)\n",
    "    plt.scatter(df['X'],df['Y'],color='red',edgecolor='black',zorder=10,label = 'Spatial Data')\n",
    "    plt.gca().add_patch(plt.Rectangle((0, 0), 100, 100, fill=False,edgecolor='black',lw=2))\n",
    "    \n",
    "    for ipts in range(0,n):\n",
    "        if ipts == 0:\n",
    "            circle1 = plt.Circle((df['X'][ipts],df['Y'][ipts]),vrange,fill=True,color='red',edgecolor=None,\n",
    "                lw=0.0,alpha=0.1,zorder=1,label='Spatial Continuity')\n",
    "        else:\n",
    "            circle1 = plt.Circle((df['X'][ipts],df['Y'][ipts]),vrange,fill=True,color='red',edgecolor=None,\n",
    "                lw=0.0,alpha=0.1,zorder=1)\n",
    "        plt.gca().add_patch(circle1)\n",
    "    \n",
    "    plt.xlim([-10,110]); plt.ylim([-10,110]); add_grid(); plt.xlabel(\"X (m)\"); plt.ylabel(\"Y (m)\") \n",
    "    plt.title('Random Spatial Dataset with Spatial Correlation Range Indicated')\n",
    "    \n",
    "    vario = GSLIB.make_variogram(nug=0.0,nst=1,it1=1,cc1=1.0,azi1=0.0,hmaj1=vrange,hmin1=vrange)\n",
    "    neffective = n_effective(df,'X','Y',seed=seed,nreal=1000,vario=vario)[0]\n",
    "    \n",
    "    neffect_array = []; l_array = []\n",
    "    for l in np.arange(1,max_L,2):\n",
    "        l_array.append(l)\n",
    "        df = pd.DataFrame(np.vstack((X[:l],Y[:l])).T,columns = ['X','Y'])\n",
    "        neffect_array.append(n_effective(df,'X','Y',seed=seed,nreal=1000,vario=vario)[0])\n",
    "    \n",
    "    df_effective = pd.DataFrame(np.vstack((l_array,neffect_array)).T,columns = ['l','e'])    \n",
    "    olsres2 = smf.ols(formula = 'e ~ l + I(l**2)+ I(l**3)-1', data = df_effective).fit()    \n",
    "      \n",
    "    plt.annotate(r'$n$ = ' + str(np.round(n,2)),(0,-5))\n",
    "    plt.annotate(r'$n_{effective}$ = ' + str(np.round(neffective,2)),(15,-5))\n",
    "    neffective_model = olsres2.params[2]*np.power(n,3) + olsres2.params[1]*np.power(n,2) + olsres2.params[0]*n\n",
    "    plt.annotate(r'$\\hat{n}_{effective}$ = ' + str(np.round(neffective_model,2)),(42,-5))\n",
    "    plt.legend(loc = 'upper left')\n",
    "        \n",
    "    ax2 = plt.subplot(122)\n",
    "    plt.plot([1,max_L],[1,max_L],color='black')\n",
    "    plt.scatter(n,neffective,c='red',edgecolor='black',s=20,zorder=10)\n",
    "    plt.scatter(n,neffective_model,c='darkorange',edgecolor='black',s=20,zorder=10)\n",
    "    plt.scatter(l_array,neffect_array,c='grey',edgecolor='grey',lw=0,s=30,alpha=0.4,zorder=5)\n",
    "    kernel = np.array([0.076923,0.230769,0.384615,0.230769,0.076923])\n",
    "    #plt.plot(np.arange(1,max_L+1,1),poly(np.arange(1,max_L+1,1)),c='red',zorder=1) \n",
    "    if window:\n",
    "        plt.plot(l_array[2:-2],np.convolve(neffect_array,kernel,mode = 'valid'),c='red',zorder=1) \n",
    "    # plt.plot(np.arange(1,max_L,3),max_pool,c='red',zorder=1) \n",
    "    if poly:\n",
    "        xs = np.arange(1,max_L,1)\n",
    "        neff_poly_model = olsres2.params[2]*np.power(xs,3) + olsres2.params[1]*np.power(xs,2) + olsres2.params[0]*xs\n",
    "        plt.plot(xs,neff_poly_model,color='darkorange',zorder=10)\n",
    "        plt.fill_between(xs,neff_poly_model,np.zeros(len(xs)),color='darkorange',alpha=0.2,zorder=1,label='Effective')\n",
    "        plt.fill_between(xs,xs,neff_poly_model,color='grey',alpha=0.2,zorder=1,label='Ineffective')\n",
    "        plt.legend(loc='upper left')\n",
    "        \n",
    "    plt.plot([0,n],[n,n],c='grey',zorder=3)\n",
    "    plt.plot([n,n],[0,n],c='grey',zorder=3)\n",
    "\n",
    "    if n > 10: \n",
    "        plt.annotate('No Spatial Correlation = ' + str(np.round(n,2)),[1.2,n+0.4],c='grey',zorder=3)\n",
    "        if abs(n - neffective) > 1.5:\n",
    "            plt.annotate('With Spatial Correlation = ' + str(np.round(neffective_model,2)),[1.2,neffective_model+0.4],c='darkorange',zorder=5)\n",
    "            plt.plot([n,n],[0,neffective_model],c='darkorange',zorder=5); \n",
    "            plt.plot([0,n],[neffective_model,neffective_model],c='darkorange',zorder=5) \n",
    "    \n",
    "    plt.ylim([1,max_L]); plt.xlim([1,max_L]); plt.xscale(\"linear\"); plt.yscale(\"linear\")\n",
    "    plt.xlabel(\"Number of Samples\"); plt.ylabel(\"Number of Independent Samples\"); plt.title(r'Number of Independent Samples ($n_{effective}$) vs. Number of Samples')\n",
    "    add_grid()\n",
    "    \n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.2); \n",
    "    plt.savefig('SpatialBootstrap.jpg',dpi=600,bbox_inches='tight')  \n",
    "    plt.show()\n",
    "    \n",
    "interactive_plot1 = widgets.interactive_output(f_make1, {'n':n,'vrange':vrange,'seed':seed,'window':window,'poly':poly})\n",
    "interactive_plot1.clear_output(wait = True)                # reduce flickering by delaying plot updating  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1cef6f901e0545a69d4601501f9b4c93",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                                                      Number of Effective Data, Mic…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f3487015831b434490c9c0e9c9e1f152",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(ui_all, interactive_plot1)                            # display the interactive plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (586471281.py, line 1)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;36m  Cell \u001b[1;32mIn[5], line 1\u001b[1;36m\u001b[0m\n\u001b[1;33m    UNDER CONSTRUCTION\u001b[0m\n\u001b[1;37m          ^\u001b[0m\n\u001b[1;31mSyntaxError\u001b[0m\u001b[1;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "UNDER CONSTRUCTION"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "nreal = 100\n",
    "\n",
    "l = widgets.Text(value=r'                                                      Number of Effective Data, Michael Pyrcz, Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "\n",
    "n = widgets.IntSlider(min=1, max = max_L, value = 1,step = 1,description = r'$n$:',orientation='horizontal',layout=Layout(width='300px', height='50px'),continuous_update=False)\n",
    "n.style.handle_color = 'gray'\n",
    "\n",
    "vrange = widgets.FloatSlider(min=0.1, max = 200.0, value = 10.0,step=5,description = r'$\\gamma_{range}$',orientation='horizontal',layout=Layout(width='300px', height='50px'),continuous_update=False)\n",
    "vrange.style.handle_color = 'gray'\n",
    "\n",
    "ireal = widgets.IntSlider(min=1, max = nreal, value = 1,step = 1,description = r'$\\ell$:',orientation='horizontal',layout=Layout(width='300px', height='50px'),continuous_update=False)\n",
    "ireal.style.handle_color = 'gray'\n",
    "\n",
    "seed = widgets.IntSlider(min=100, max = 999, value = 1, description = r'$s$:',orientation='horizontal',layout=Layout(width='300px', height='50px'),continuous_update=False)\n",
    "seed.style.handle_color = 'gray'\n",
    "\n",
    "window = widgets.Checkbox(value=False,description='Conv. Fit',disabled=False,layout=Layout(width='200px', height='30px'))\n",
    "poly = widgets.Checkbox(value=True,description='Poly. Fit',disabled=False,layout=Layout(width='200px', height='30px'))\n",
    "\n",
    "ui1 = widgets.HBox([n,vrange,ireal,seed,],kwargs = {'justify_content':'center'})\n",
    "ui2 = widgets.HBox([window,poly,],kwargs = {'justify_content':'center'})\n",
    "ui_all2 = widgets.VBox([l,ui1,ui2],)\n",
    "\n",
    "def f_make2(n,vrange,ireal,seed,window,poly):\n",
    "    \n",
    "    nreal = 100;\n",
    "    maxf = 50; maxfboot = 15\n",
    "    dx = -5; dy = 3; pdfy = 13.0\n",
    "    \n",
    "    np.random.seed(seed=seed)\n",
    "    values = np.random.rand(max_L*2)*100\n",
    "    X,Y = np.split(values,2)\n",
    "    df = pd.DataFrame(np.vstack((X[:n],Y[:n])).T,columns = ['X','Y'])\n",
    "    \n",
    "    ax1 = plt.subplot(121)\n",
    "    plt.scatter(df['X'],df['Y'],marker='x',s=30,color='red',edgecolor='black',lw=2,zorder=10,label = 'Spatial Data')\n",
    "    plt.scatter(df['X'],df['Y'],marker='x',s=40,color='white',edgecolor='black',lw=4,zorder=9,label = 'Spatial Data')\n",
    "    plt.gca().add_patch(plt.Rectangle((0, 0), 100, 100, fill=False,edgecolor='black',lw=2))\n",
    "    \n",
    "    xx = np.arange(-3,3,0.1)\n",
    "    pdf = norm.pdf(xx, loc=0, scale=1) * pdfy\n",
    "    px = np.linspace(0,10,len(xx))\n",
    "    \n",
    "    for ipts in range(0,n):\n",
    "        if ipts == 0:\n",
    "            circle1 = plt.Circle((df['X'][ipts],df['Y'][ipts]),vrange,fill=True,color='red',edgecolor=None,\n",
    "                lw=0.0,alpha=0.05,zorder=1,label='Spatial Continuity')\n",
    "        else:\n",
    "            circle1 = plt.Circle((df['X'][ipts],df['Y'][ipts]),vrange,fill=True,color='red',edgecolor=None,\n",
    "                lw=0.0,alpha=0.05,zorder=1)\n",
    "        plt.gca().add_patch(circle1)\n",
    "        \n",
    "        plt.plot([df['X'][ipts]+dx,df['X'][ipts]+10+dx],[df['Y'][ipts]+dy,df['Y'][ipts]+dy],color='black')\n",
    "        plt.plot([df['X'][ipts]+dx,df['X'][ipts]+dx],[df['Y'][ipts]+dy,df['Y'][ipts]+5+dy],color='black')\n",
    "        plt.plot(px + df['X'][ipts] + dx,pdf + df['Y'][ipts] + dy,color = 'black',zorder=1)\n",
    "    \n",
    "    plt.xlim([-10,110]); plt.ylim([-10,110]); add_grid(); plt.xlabel(\"X (m)\"); plt.ylabel(\"Y (m)\") \n",
    "    plt.title('Random Spatial Dataset with Spatial Correlation Range Indicated')\n",
    "    \n",
    "    vario = GSLIB.make_variogram(nug=0.0,nst=1,it1=1,cc1=1.0,azi1=0.0,hmaj1=vrange,hmin1=vrange)\n",
    "    neffective,realizations = n_effective(df,'X','Y',seed=seed,nreal=nreal,vario=vario)\n",
    "    \n",
    "    for ipts in range(0,n):\n",
    "        y = df['Y'][ipts]+dy; x = df['X'][ipts]+dx+(realizations[ireal-1,ipts]-(-3.0))/(3.0-(-3.0))*10\n",
    "        plt.plot([x,x],[y,y+5],color='black',lw=3,zorder=1); plt.plot([x,x],[y,y+5],color='red',lw=1,zorder=3)\n",
    "    \n",
    "    ax2 = plt.subplot(122)\n",
    "    plt.axis('off'); plt.xlim([0, 10]); plt.ylim([0, 10])\n",
    "    \n",
    "    plt.plot([2,8],[6,6],color='black'); plt.plot([2,2],[6,9],color='black')\n",
    "    \n",
    "    xpmin = 2; xpmax = 8; ypmin = 6; ypmax = 9\n",
    "    xprange = xpmax - xpmin; yprange = ypmax - ypmin\n",
    "    xmin = -3.0; xmax = 3.0; nbin = 16\n",
    "    xhalf = (xmax-xmin)/(nbin-1)/2.0; xsize = xhalf*2.0\n",
    "    xrange = xmax - xmin\n",
    "    xbins = np.linspace(xmin,xmax,nbin)\n",
    "    ybins = np.linspace(0.0,maxfboot,11)\n",
    "    xcents = np.linspace(xmin + xhalf,xmax-xhalf,nbin-1)\n",
    "    \n",
    "    for ibin, xbin in enumerate(xbins):\n",
    "        xx = ((xbin-xmin)/xrange*xprange)+xpmin\n",
    "        plt.plot([xx,xx],[ypmin,ypmax],color='grey',lw=0.2,zorder=1)\n",
    "        if ibin % 2 == 0:\n",
    "            plt.plot([xx,xx],[ypmin-yprange*0.1,ypmin],color='black',zorder=3)\n",
    "            plt.annotate(np.round(xbin,1),(xx,ypmin-yprange*0.18),ha='center')\n",
    "        else: \n",
    "            plt.plot([xx,xx],[ypmin-yprange*0.05,ypmin],color='black',zorder=3)\n",
    "            plt.annotate(np.round(xbin,1),(xx,ypmin-yprange*0.13),ha='center')\n",
    "    \n",
    "#     for xbin in xbins:\n",
    "#         xx = ((xbin-xmin)/xrange*xprange)+xpmin\n",
    "#         plt.plot([xx,xx],[ypmin-yprange*0.05,ypmin],color='black',zorder=3)\n",
    "#         plt.annotate(np.round(xbin,1),(xx,ypmin-yprange*0.15),ha='center')\n",
    "#         plt.plot([xx,xx],[ypmin,ypmax],color='grey',lw=0.2,zorder=1)    \n",
    "    \n",
    "    for ybin in ybins:\n",
    "        yy = (ybin)*3/maxfboot+6\n",
    "        plt.plot([xpmin-xprange*0.04,xpmin],[yy,yy],color='black',zorder=3)\n",
    "        plt.plot([xpmin,xpmax],[yy,yy],color='grey',lw=0.2,zorder=1)\n",
    "        plt.annotate(np.round(ybin,1),(xpmin-xprange*0.05,yy),ha='right')\n",
    "    \n",
    "    histboot = np.histogram(realizations[ireal-1,:], bins=xbins, weights=None)[0]\n",
    "    average = np.average(realizations[ireal-1,:])\n",
    "    \n",
    "    for ibin, prop in enumerate(histboot):\n",
    "        xx = ((xcents[ibin]-xmin)/xrange*xprange)+xpmin\n",
    "        yy = histboot[ibin]*yprange/maxfboot\n",
    "        ax2.add_patch(plt.Rectangle((xx-xhalf,ypmin), xsize, yy, lw=1, fc = 'darkorange',color='black', ))\n",
    "        xavg = ((average-xmin)/xrange*xprange)+xpmin\n",
    "        #plt.plot([xavg,xavg],[6,9],color='red')\n",
    "        \n",
    "    for ipts in range(0,n):\n",
    "        xx = ((realizations[ireal-1,ipts]-xmin)/xrange*xprange)+xpmin \n",
    "        plt.scatter(xx,ypmin,color='red',edgecolor='black',s=20,alpha=1.0,zorder=100)\n",
    "        \n",
    "    plt.annotate('Spatial Bootstrap Distribution Realization',(xpmin+xprange*0.5,ypmax+yprange*0.08),ha='center') \n",
    "    plt.annotate('Feature (NSCORE)',(xpmin+xprange*0.5,ypmin-yprange*0.25),ha='center') \n",
    "    plt.annotate('Proportion',(xpmin-xprange*0.18,ypmin+yprange*0.5),va='center',rotation = 90.0) \n",
    "    \n",
    "    plt.plot([2,8],[1,1],color='black'); plt.plot([2,2],[1,4],color='black')\n",
    "    \n",
    "    xpmin = 2; xpmax = 8; ypmin = 1; ypmax = 4\n",
    "    xprange = xpmax - xpmin; yprange = ypmax - ypmin\n",
    "    xmin = -3.0; xmax = 3.0; nbin = 16\n",
    "    xhalf = (xmax-xmin)/(nbin-1)/2.0; xsize = xhalf*2.0\n",
    "    xrange = xmax - xmin\n",
    "    xbins = np.linspace(xmin,xmax,nbin)\n",
    "    ybins = np.linspace(0.0,maxf,11)\n",
    "    xcents = np.linspace(xmin + xhalf,xmax-xhalf,nbin-1)\n",
    "    \n",
    "    for ibin, xbin in enumerate(xbins):\n",
    "        xx = ((xbin-xmin)/xrange*xprange)+xpmin\n",
    "        plt.plot([xx,xx],[ypmin,ypmax],color='grey',lw=0.2,zorder=1)\n",
    "        if ibin % 2 == 0:\n",
    "            plt.plot([xx,xx],[ypmin-yprange*0.1,ypmin],color='black',zorder=3)\n",
    "            plt.annotate(np.round(xbin,1),(xx,ypmin-yprange*0.18),ha='center')\n",
    "        else: \n",
    "            plt.plot([xx,xx],[ypmin-yprange*0.05,ypmin],color='black',zorder=3)\n",
    "            plt.annotate(np.round(xbin,1),(xx,ypmin-yprange*0.13),ha='center')\n",
    "                 \n",
    "    for ybin in ybins:\n",
    "        yy = (ybin)*yprange/maxf+1\n",
    "        plt.plot([xpmin-xprange*0.04,xpmin],[yy,yy],color='black',zorder=3)\n",
    "        plt.plot([xpmin,xpmax],[yy,yy],color='grey',lw=0.2,zorder=1)\n",
    "        plt.annotate(np.round(ybin,1),(xpmin-xprange*0.05,yy),ha='right')\n",
    "    \n",
    "    average = np.average(realizations[ireal-1,:])\n",
    "    \n",
    "    averages = np.average(realizations,axis=1)\n",
    "    #hist_avg = np.histogram(averages[:ireal+1], bins=xbins, weights=None)[0]\n",
    "    hist_avg = np.histogram(averages[:ireal], bins=xbins, weights=None)[0]\n",
    "    \n",
    "    for ibin, prop in enumerate(hist_avg):\n",
    "        xx = ((xcents[ibin]-(-3.0))/(3-(-3))*6)+2\n",
    "        yy = hist_avg[ibin]*3/maxf\n",
    "        ax2.add_patch(plt.Rectangle((xx-xhalf,ypmin), xsize, yy, lw=1, fc = 'gold',color='black',zorder=10))\n",
    "        xavg = ((average-(-3.0))/(3.0-(-3.0))*6)+2\n",
    "        plt.plot([xavg,xavg],[6,9],color='red')\n",
    "    \n",
    "    #hist_avg_prior = np.histogram(averages[:ireal], bins=xbins, weights=None)[0]\n",
    "    if ireal > 1:\n",
    "        hist_avg_prior = np.histogram(averages[:ireal-1], bins=xbins, weights=None)[0]\n",
    "    \n",
    "        for ibin, prop in enumerate(hist_avg):\n",
    "            xx = ((xcents[ibin]-(-3.0))/(3-(-3))*6)+2\n",
    "            yy = hist_avg_prior[ibin]*3/maxf\n",
    "            ax2.add_patch(plt.Rectangle((xx-xhalf,ypmin), xsize, yy, lw=1, fc = 'red',color='black',zorder=15))\n",
    "            xavg = ((average-(-3.0))/(3.0-(-3.0))*6)+2\n",
    "            #plt.plot([xavg,xavg],[6,9],color='red')\n",
    "        \n",
    "    plt.annotate('Spatial Bootstrap Average Realization',(xpmin+xprange*0.5,ypmax+yprange*0.08),ha='center') \n",
    "    plt.annotate('Bootstrap Average (NSCORE)',(xpmin+xprange*0.5,ypmin-yprange*0.30),ha='center') \n",
    "    plt.annotate('Frequency',(0.8,2.5),va='center',rotation = 90.0) \n",
    "    \n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.2)\n",
    "    plt.savefig('SpatialBootstrap.jpg',dpi=600,bbox_inches='tight')  \n",
    "    plt.show()\n",
    "    \n",
    "interactive_plot2 = widgets.interactive_output(f_make2, {'n':n,'vrange':vrange,'ireal':ireal,'seed':seed,'window':window,'poly':poly})\n",
    "interactive_plot2.clear_output(wait = True)                # reduce flickering by delaying plot updating     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d09a1a2f52c94ed9acc03d8332c16575",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                                                      Number of Effective Data, Mic…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3648f48c68714bf5a0737d1f1da96b95",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(ui_all2, interactive_plot2)                            # display the interactive plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.57744207, -0.12475772,  0.97571096,  1.59494036, -1.20194452,\n",
       "       -1.28620972,  0.82640545,  0.78747084,  0.65703186,  1.50348689])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nreal = 100\n",
    "\n",
    "l = widgets.Text(value=r'                                                      Number of Effective Data, Michael Pyrcz, Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "\n",
    "n = widgets.IntSlider(min=1, max = max_L, value = 1,step = 1,description = r'$n$:',orientation='horizontal',layout=Layout(width='300px', height='50px'),continuous_update=False)\n",
    "n.style.handle_color = 'gray'\n",
    "\n",
    "vrange = widgets.FloatSlider(min=0.1, max = 200.0, value = 10.0,step=5,description = r'$\\gamma_{range}$',orientation='horizontal',layout=Layout(width='300px', height='50px'),continuous_update=False)\n",
    "vrange.style.handle_color = 'gray'\n",
    "\n",
    "ireal = widgets.IntSlider(min=1, max = nreal, value = 1,step = 1,description = r'$\\ell$:',orientation='horizontal',layout=Layout(width='300px', height='50px'),continuous_update=False)\n",
    "ireal.style.handle_color = 'gray'\n",
    "\n",
    "seed = widgets.IntSlider(min=100, max = 999, value = 1, description = r'$s$:',orientation='horizontal',layout=Layout(width='300px', height='50px'),continuous_update=False)\n",
    "seed.style.handle_color = 'gray'\n",
    "\n",
    "window = widgets.Checkbox(value=False,description='Conv. Fit',disabled=False,layout=Layout(width='200px', height='30px'))\n",
    "poly = widgets.Checkbox(value=True,description='Poly. Fit',disabled=False,layout=Layout(width='200px', height='30px'))\n",
    "\n",
    "ui1 = widgets.HBox([n,vrange,ireal,seed,],kwargs = {'justify_content':'center'})\n",
    "ui2 = widgets.HBox([window,poly,],kwargs = {'justify_content':'center'})\n",
    "ui_all2 = widgets.VBox([l,ui1,ui2],)\n",
    "\n",
    "def f_make2(n,vrange,ireal,seed,window,poly):\n",
    "    \n",
    "    nreal = 100;\n",
    "    maxf = 50; maxfboot = 15\n",
    "    dx = -5; dy = 3; pdfy = 13.0\n",
    "    \n",
    "    np.random.seed(seed=seed)\n",
    "    values = np.random.rand(max_L*2)*100\n",
    "    X,Y = np.split(values,2)\n",
    "    df = pd.DataFrame(np.vstack((X[:n],Y[:n])).T,columns = ['X','Y'])\n",
    "    \n",
    "    ax1 = plt.subplot(121)\n",
    "    plt.scatter(df['X'],df['Y'],marker='x',s=30,color='red',edgecolor='black',lw=2,zorder=10,label = 'Spatial Data')\n",
    "    plt.scatter(df['X'],df['Y'],marker='x',s=40,color='white',edgecolor='black',lw=4,zorder=9,label = 'Spatial Data')\n",
    "    plt.gca().add_patch(plt.Rectangle((0, 0), 100, 100, fill=False,edgecolor='black',lw=2))\n",
    "    \n",
    "    xx = np.arange(-3,3,0.1)\n",
    "    pdf = norm.pdf(xx, loc=0, scale=1) * pdfy\n",
    "    px = np.linspace(0,10,len(xx))\n",
    "    \n",
    "    for ipts in range(0,n):\n",
    "        if ipts == 0:\n",
    "            circle1 = plt.Circle((df['X'][ipts],df['Y'][ipts]),vrange,fill=True,color='red',edgecolor=None,\n",
    "                lw=0.0,alpha=0.05,zorder=1,label='Spatial Continuity')\n",
    "        else:\n",
    "            circle1 = plt.Circle((df['X'][ipts],df['Y'][ipts]),vrange,fill=True,color='red',edgecolor=None,\n",
    "                lw=0.0,alpha=0.05,zorder=1)\n",
    "        plt.gca().add_patch(circle1)\n",
    "        \n",
    "        plt.plot([df['X'][ipts]+dx,df['X'][ipts]+10+dx],[df['Y'][ipts]+dy,df['Y'][ipts]+dy],color='black')\n",
    "        plt.plot([df['X'][ipts]+dx,df['X'][ipts]+dx],[df['Y'][ipts]+dy,df['Y'][ipts]+5+dy],color='black')\n",
    "        plt.plot(px + df['X'][ipts] + dx,pdf + df['Y'][ipts] + dy,color = 'black',zorder=1)\n",
    "    \n",
    "    plt.xlim([-10,110]); plt.ylim([-10,110]); add_grid(); plt.xlabel(\"X (m)\"); plt.ylabel(\"Y (m)\") \n",
    "    plt.title('Random Spatial Dataset with Spatial Correlation Range Indicated')\n",
    "    \n",
    "    vario = GSLIB.make_variogram(nug=0.0,nst=1,it1=1,cc1=1.0,azi1=0.0,hmaj1=vrange,hmin1=vrange)\n",
    "    neffective,realizations = n_effective(df,'X','Y',seed=seed,nreal=nreal,vario=vario)\n",
    "    \n",
    "    for ipts in range(0,n):\n",
    "        y = df['Y'][ipts]+dy; x = df['X'][ipts]+dx+(realizations[ireal,ipts]-(-3.0))/(3.0-(-3.0))*10\n",
    "        plt.plot([x,x],[y,y+5],color='black',lw=3,zorder=1); plt.plot([x,x],[y,y+5],color='red',lw=1,zorder=3)\n",
    "    \n",
    "    ax2 = plt.subplot(122)\n",
    "    plt.axis('off'); plt.xlim([0, 10]); plt.ylim([0, 10])\n",
    "    \n",
    "    plt.plot([2,8],[6,6],color='black'); plt.plot([2,2],[6,9],color='black')\n",
    "    \n",
    "    xpmin = 2; xpmax = 8; ypmin = 6; ypmax = 9\n",
    "    xprange = xpmax - xpmin; yprange = ypmax - ypmin\n",
    "    xmin = -3.0; xmax = 3.0; nbin = 6\n",
    "    xhalf = (xmax-xmin)/(nbin-1)/2.0; xsize = xhalf*2.0\n",
    "    xrange = xmax - xmin\n",
    "    xbins = np.linspace(xmin,xmax,nbin)\n",
    "    ybins = np.linspace(0.0,maxfboot,11)\n",
    "    xcents = np.linspace(xmin + xhalf,xmax-xhalf,nbin-1)\n",
    "    \n",
    "    for xbin in xbins:\n",
    "        xx = ((xbin-xmin)/xrange*xprange)+xpmin\n",
    "        plt.plot([xx,xx],[ypmin-yprange*0.05,ypmin],color='black',zorder=3)\n",
    "        plt.annotate(np.round(xbin,1),(xx,ypmin-yprange*0.15),ha='center')\n",
    "        plt.plot([xx,xx],[ypmin,ypmax],color='grey',lw=0.2,zorder=1)    \n",
    "    \n",
    "    for ybin in ybins:\n",
    "        yy = (ybin)*3/maxfboot+6\n",
    "        plt.plot([xpmin-xprange*0.04,xpmin],[yy,yy],color='black',zorder=3)\n",
    "        plt.plot([xpmin,xpmax],[yy,yy],color='grey',lw=0.2,zorder=1)\n",
    "        plt.annotate(np.round(ybin,1),(xpmin-xprange*0.05,yy),ha='right')\n",
    "    \n",
    "    histboot = np.histogram(realizations[ireal,:], bins=xbins, weights=None)[0]\n",
    "    average = np.average(realizations[ireal,:])\n",
    "    \n",
    "    for ibin, prop in enumerate(histboot):\n",
    "        xx = ((xcents[ibin]-xmin)/xrange*xprange)+xpmin\n",
    "        yy = histboot[ibin]*yprange/maxfboot\n",
    "        ax2.add_patch(plt.Rectangle((xx-xhalf,ypmin), xsize, yy, lw=1, fc = 'darkorange',color='black', ))\n",
    "        xavg = ((average-xmin)/xrange*xprange)+xpmin\n",
    "        #plt.plot([xavg,xavg],[6,9],color='red')\n",
    "        \n",
    "    for ipts in range(0,n):\n",
    "        xx = ((realizations[ireal,ipts]-xmin)/xrange*xprange)+xpmin \n",
    "        plt.scatter(xx,ypmin,color='red',edgecolor='black',s=20,alpha=1.0,zorder=100)\n",
    "        \n",
    "    plt.annotate('Spatial Bootstrap Distribution Realization',(xpmin+xprange*0.5,ypmax+yprange*0.08),ha='center') \n",
    "    plt.annotate('Feature (NSCORE)',(xpmin+xprange*0.5,ypmin-yprange*0.25),ha='center') \n",
    "    plt.annotate('Proportion',(xpmin-xprange*0.18,ypmin+yprange*0.5),va='center',rotation = 90.0) \n",
    "    \n",
    "    plt.plot([2,8],[1,1],color='black'); plt.plot([2,2],[1,4],color='black')\n",
    "    \n",
    "    xpmin = 2; xpmax = 8; ypmin = 1; ypmax = 4\n",
    "    xprange = xpmax - xpmin; yprange = ypmax - ypmin\n",
    "    xmin = -3.0; xmax = 3.0; nbin = 16\n",
    "    xhalf = (xmax-xmin)/(nbin-1)/2.0; xsize = xhalf*2.0\n",
    "    xrange = xmax - xmin\n",
    "    xbins = np.linspace(xmin,xmax,nbin)\n",
    "    ybins = np.linspace(0.0,maxf,11)\n",
    "    xcents = np.linspace(xmin + xhalf,xmax-xhalf,nbin-1)\n",
    "    \n",
    "    for ibin, xbin in enumerate(xbins):\n",
    "        xx = ((xbin-xmin)/xrange*xprange)+xpmin\n",
    "        plt.plot([xx,xx],[ypmin,ypmax],color='grey',lw=0.2,zorder=1)\n",
    "        if ibin % 2 == 0:\n",
    "            plt.plot([xx,xx],[ypmin-yprange*0.1,ypmin],color='black',zorder=3)\n",
    "            plt.annotate(np.round(xbin,1),(xx,ypmin-yprange*0.18),ha='center')\n",
    "        else: \n",
    "            plt.plot([xx,xx],[ypmin-yprange*0.05,ypmin],color='black',zorder=3)\n",
    "            plt.annotate(np.round(xbin,1),(xx,ypmin-yprange*0.13),ha='center')\n",
    "                 \n",
    "    for ybin in ybins:\n",
    "        yy = (ybin)*yprange/maxf+1\n",
    "        plt.plot([xpmin-xprange*0.04,xpmin],[yy,yy],color='black',zorder=3)\n",
    "        plt.plot([xpmin,xpmax],[yy,yy],color='grey',lw=0.2,zorder=1)\n",
    "        plt.annotate(np.round(ybin,1),(xpmin-xprange*0.05,yy),ha='right')\n",
    "    \n",
    "    average = np.average(realizations[ireal,:])\n",
    "    \n",
    "    averages = np.average(realizations,axis=1)\n",
    "    #hist_avg = np.histogram(averages[:ireal+1], bins=xbins, weights=None)[0]\n",
    "    \n",
    "    for ibin, prop in enumerate(hist_avg):\n",
    "        xx = ((xcents[ibin]-(-3.0))/(3-(-3))*6)+2\n",
    "        yy = hist_avg[ibin]*3/maxf\n",
    "        ax2.add_patch(plt.Rectangle((xx-xhalf,ypmin), xsize, yy, lw=1, fc = 'gold',color='black',zorder=10))\n",
    "        xavg = ((average-(-3.0))/(3.0-(-3.0))*6)+2\n",
    "        plt.plot([xavg,xavg],[6,9],color='red')\n",
    "    \n",
    "    hist_avg_prior = np.histogram(averages[:ireal], bins=xbins, weights=None)[0]\n",
    "    \n",
    "    for ibin, prop in enumerate(hist_avg):\n",
    "        xx = ((xcents[ibin]-(-3.0))/(3-(-3))*6)+2\n",
    "        yy = hist_avg_prior[ibin]*3/maxf\n",
    "        ax2.add_patch(plt.Rectangle((xx-xhalf,ypmin), xsize, yy, lw=1, fc = 'red',color='black',zorder=15))\n",
    "        xavg = ((average-(-3.0))/(3.0-(-3.0))*6)+2\n",
    "        #plt.plot([xavg,xavg],[6,9],color='red')\n",
    "        \n",
    "    plt.annotate('Spatial Bootstrap Average Realization',(xpmin+xprange*0.5,ypmax+yprange*0.08),ha='center') \n",
    "    plt.annotate('Bootstrap Average (NSCORE)',(xpmin+xprange*0.5,ypmin-yprange*0.30),ha='center') \n",
    "    plt.annotate('Frequency',(0.8,2.5),va='center',rotation = 90.0) \n",
    "    \n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.2, hspace=0.2)\n",
    "    plt.savefig('SpatialBootstrap.jpg',dpi=600,bbox_inches='tight')  \n",
    "    plt.show()\n",
    "    \n",
    "interactive_plot2 = widgets.interactive_output(f_make2, {'n':n,'vrange':vrange,'ireal':ireal,'seed':seed,'window':window,'poly':poly})\n",
    "interactive_plot2.clear_output(wait = True)                # reduce flickering by delaying plot updating     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-2.8, -2.4, -2. , -1.6, -1.2, -0.8, -0.4,  0. ,  0.4,  0.8,  1.2,\n",
       "        1.6,  2. ,  2.4,  2.8])"
      ]
     },
     "execution_count": 255,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xmin = -3.0; xmax = 3.0; nbin = 16\n",
    "xhalf = (xmax-xmin)/(nbin-1)/2.0\n",
    "xbins = np.linspace(xmin,xmax,nbin)\n",
    "ybins = np.linspace(0.0,maxf,11)\n",
    "xcents = np.linspace(xmin + xhalf,xmax-xhalf,nbin-1)\n",
    "\n",
    "xcents"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 0, 3, 1, 0, 0], dtype=int64)"
      ]
     },
     "execution_count": 206,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hist_avg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-3., -2., -1.,  0.,  1.,  2.,  3.])"
      ]
     },
     "execution_count": 207,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xbins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 210,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0, 0, 3, 1, 0, 0], dtype=int64),\n",
       " array([-3., -2., -1.,  0.,  1.,  2.,  3.]))"
      ]
     },
     "execution_count": 210,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "averages[:ireal]\n",
    "\n",
    "np.histogram(averages[:ireal], bins=xbins, weights=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0, 0, 2, 1, 0, 0], dtype=int64),\n",
       " array([-3., -2., -1.,  0.,  1.,  2.,  3.]))"
      ]
     },
     "execution_count": 212,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.histogram(averages[:ireal-1], bins=xbins, weights=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 211,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.10180635,  0.3473193 , -0.08361122, -0.4094051 ])"
      ]
     },
     "execution_count": 211,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "averages[:ireal]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.10180635,  0.3473193 , -0.08361122])"
      ]
     },
     "execution_count": 213,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "averages[:ireal-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 214,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.24737254766389993"
      ]
     },
     "execution_count": 214,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.average(realizations[ireal,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 217,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.24737254766389993"
      ]
     },
     "execution_count": 217,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "averages[ireal]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 215,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 215,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ireal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxwAAAMqCAYAAAAIJEe4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAADIwElEQVR4nOzdeVhU1f8H8PfMADPACO4sRmq44r4LSmpR5JotVtrXrbINK6WyLNPU0rQsWzDLMsuyRTPzl6aZaZqSG2rumvsGrgjDzsz9/THMxMAA9w5zmXuH9+t5eNA75577udxZzmfOPedoBEEQQEREREREJAOtpwMgIiIiIiLvxYSDiIiIiIhkw4SDiIiIiIhkw4SDiIiIiIhkw4SDiIiIiIhkw4SDiIiIiIhkw4SDiIiIiIhkw4SDiIiIiIhkw4SDiIiIiIhk4xUJx6hRo9CoUSNPh6F4vXv3Ru/evV3aV6PR4PXXX3drPFS+U6dOQaPRYNGiRaLLvvPOO/IHVoUq89pu1KgRRo0a5dZ4lEKO12Nl3h/IM5y9R7z++uvQaDSeC6oS+DlD5L0kJRyLFi2CRqOx//j4+KBBgwYYNWoUzp8/L1eMqrRv3z7cf//9aNiwIQwGAxo0aIA77rgDH374oazHPXjwIF5//XWcOnVK1uM4Y/vws/34+vqibt26iImJwSuvvIIzZ864XPeFCxfw+uuvY8+ePe4LuBJWr17tkQ9GOY976tQpjB49GpGRkTAYDAgNDcWtt96KKVOmyHI8GyVc29zcXLz33nvo1q0bgoODYTAY0KxZM4wdOxZHjx71WFzu5sn3h7Js3LjR4X1Dp9Ohfv36uP/++3Ho0CFPh+c2ts/PnTt3ejqUSlmyZAnmzp3r6TCISGV8XNlp2rRpaNy4MXJzc/H3339j0aJF+Ouvv7B//34YDAZ3x6g6W7duRZ8+fXDzzTdjzJgxCA0NxdmzZ/H333/j/fffxzPPPCPbsQ8ePIipU6eid+/epb4Z/u2332Q7bnFDhw5Fv379YLFYcP36dezYsQNz587F+++/j88//xwPPfSQ5DovXLiAqVOnolGjRmjfvr37g5Zo9erVSEpKkjXpaNiwIXJycuDr6yv7cf/991906dIF/v7+eOSRR9CoUSNcvHgRKSkpmDVrFqZOnerW4xVX3rVdsGABLBaLbMcGgCtXruCuu+7Crl27MGDAAAwbNgxGoxFHjhzBd999h08//RT5+fmyxlBVlPD+UJZnn30WXbp0QUFBAf755x/Mnz8fGzduxP79+xEaGurR2NRk0qRJePnll2Wrf8mSJdi/fz/GjRsn2zGIyPu4lHD07dsXnTt3BgA89thjqFu3LmbNmoWVK1figQcecGuAavTmm28iODgYO3bsQM2aNR0eu3TpkmeCAuDn51clx+nYsSP+97//OWw7ffo07rzzTowcORItW7ZEu3btqiQWNdNoNFWWwL/33nswmUzYs2cPGjZs6PCYJ5+zxZMtuYwaNQq7d+/GsmXLcN999zk8Nn36dLz66qtuOU5WVhYCAwNLbRcEAbm5ufD393fLcVxVVe8PZYmNjcX9999v/3/z5s3x1FNP4auvvsKECRM8GJm6+Pj4wMfHpY92IiLZuGUMR2xsLADg+PHj9m35+fmYPHkyOnXqhODgYAQGBiI2NhYbNmxw2Lf4veeffvopIiMjodfr0aVLF+zYsaPUsVasWIHWrVvDYDCgdevW+Omnn5zGlJWVheeffx4RERHQ6/Vo3rw53nnnHQiC4FBOo9Fg7NixWLp0KaKiouDv74/o6Gjs27cPAPDJJ5+gSZMmMBgM6N27t6hbEY4fP45WrVqVSjYAoH79+k6P/80336B58+YwGAzo1KkTNm3a5FDu9OnTePrpp9G8eXP4+/ujTp06GDJkiEM8ixYtwpAhQwAAffr0sd+isHHjRgCl79EWe43coWHDhli0aBHy8/Mxe/Zs+/Zr167hhRdeQJs2bWA0GhEUFIS+ffti79699jIbN25Ely5dAACjR4+2n5ftvuXNmzdjyJAhuPnmm6HX6xEREYHx48cjJyfHIYbU1FSMHj0aN910E/R6PcLCwnD33XeXuqa//vorYmNjERgYiBo1aqB///44cOCA/fFRo0YhKSkJABxuBSlLYmIi6tSp4/Dce+aZZ6DRaPDBBx/Yt6WlpUGj0eDjjz8GUPr+bLHHFfM6Kun48eO46aabSiUbQOnnbKNGjTBgwAD89ttvaN++PQwGA6KiorB8+XKHcu64ts7GcLzzzjuIiYlBnTp14O/vj06dOmHZsmUVnqMz27Ztw6pVq/Doo4+WSjYAQK/XlxoX88cff9ifHzVr1sTdd99d6tYf2330Bw8exLBhw1CrVi307NnT4e+3du1adO7cGf7+/vjkk08AAOnp6Rg3bpz9fatJkyaYNWtWhb08crw/ANZk89FHH0VISAgMBgPatWuHL7/80qGM1PdwsZx9rgDir7/tvdX2maHX69GqVSusWbOmVNmNGzeic+fOMBgMiIyMxCeffFLmWIivv/4anTp1gr+/P2rXro2HHnoIZ8+edekcR40aBaPRiPPnz2Pw4MEwGo2oV68eXnjhBZjNZoey6enpGDVqFIKDg1GzZk2MHDkS6enppeosL+6uXbsiICAAtWrVwq233urQq/Xzzz+jf//+CA8Ph16vR2RkJKZPn+4QR+/evbFq1SqcPn3a/vwp/vrMy8vDlClT0KRJE/t78YQJE5CXl+cQS15eHsaPH4969eqhRo0aGDRoEM6dO+fS35CI1MEtX4PYPtRq1apl35aRkYHPPvsMQ4cOxZgxY5CZmYnPP/8c8fHx2L59e6lbJ5YsWYLMzEw88cQT0Gg0mD17Nu69916cOHHC/i3nb7/9hvvuuw9RUVGYOXMmrl69am9AFicIAgYNGoQNGzbg0UcfRfv27bF27Vq8+OKLOH/+PN577z2H8ps3b8bKlSuRkJAAAJg5cyYGDBiACRMmYN68eXj66adx/fp1zJ49G4888gj++OOPcv8eDRs2RHJyMvbv34/WrVtX+Pf7888/8f333+PZZ5+FXq/HvHnzcNddd2H79u32/Xfs2IGtW7fioYcewk033YRTp07h448/Ru/evXHw4EEEBATg1ltvxbPPPosPPvgAr7zyClq2bAkA9t8lSb1GlRUdHY3IyEisW7fOvu3EiRNYsWIFhgwZgsaNGyMtLQ2ffPIJevXqhYMHDyI8PBwtW7bEtGnTMHnyZDz++OP2hkhMTAwAYOnSpcjOzsZTTz2FOnXqYPv27fjwww9x7tw5LF261H6s++67DwcOHMAzzzyDRo0a4dKlS1i3bh3OnDlj/9BcvHgxRo4cifj4eMyaNQvZ2dn4+OOP0bNnT+zevRuNGjXCE088gQsXLmDdunVYvHhxhecdGxuL9957DwcOHLBfz82bN0Or1WLz5s149tln7dsA4NZbb3Vaj5jjinkdOdOwYUP8/vvv+OOPP3DbbbdVeE7Hjh3Dgw8+iCeffBIjR47EF198gSFDhmDNmjW44447ALjn2jrz/vvvY9CgQXj44YeRn5+P7777DkOGDMEvv/yC/v37Vxh7cStXrgQADB8+XFT533//HX379sUtt9yC119/HTk5Ofjwww/Ro0cPpKSklEqOhgwZgqZNm2LGjBkOCeeRI0cwdOhQPPHEExgzZgyaN2+O7Oxs9OrVC+fPn8cTTzyBm2++GVu3bsXEiRNx8eLFcu+bl+P9IScnB71798a///6LsWPHonHjxli6dClGjRqF9PR0PPfccw7lXX3ulcXZ5wog7fr/9ddfWL58OZ5++mnUqFEDH3zwAe677z6cOXMGderUAQDs3r0bd911F8LCwjB16lSYzWZMmzYN9erVKxXTm2++iddeew0PPPAAHnvsMVy+fBkffvghbr31Vuzevdvpl0wVMZvNiI+PR7du3fDOO+/g999/x5w5cxAZGYmnnnoKgPUz7e6778Zff/2FJ598Ei1btsRPP/2EkSNHijrG1KlT8frrryMmJgbTpk2Dn58ftm3bhj/++AN33nknAGtCajQakZiYCKPRiD/++AOTJ09GRkYG3n77bQDAq6++ihs3buDcuXP2z1Gj0QgAsFgsGDRoEP766y88/vjjaNmyJfbt24f33nsPR48exYoVK+zxPPbYY/j6668xbNgwxMTE4I8//pD82iUilREk+OKLLwQAwu+//y5cvnxZOHv2rLBs2TKhXr16gl6vF86ePWsvW1hYKOTl5Tnsf/36dSEkJER45JFH7NtOnjwpABDq1KkjXLt2zb79559/FgAI//d//2ff1r59eyEsLExIT0+3b/vtt98EAELDhg3t21asWCEAEN544w2H499///2CRqMR/v33X/s2AIJerxdOnjxp3/bJJ58IAITQ0FAhIyPDvn3ixIkCAIeyzvz222+CTqcTdDqdEB0dLUyYMEFYu3atkJ+fX6osAAGAsHPnTvu206dPCwaDQbjnnnvs27Kzs0vtm5ycLAAQvvrqK/u2pUuXCgCEDRs2lCrfq1cvoVevXvb/i71GtjinTJlS7nnbruXbb79dZpm7775bACDcuHFDEARByM3NFcxmc6l69Hq9MG3aNPu2HTt2CACEL774olSdzv42M2fOFDQajXD69Gn7eVUUW2ZmplCzZk1hzJgxDttTU1OF4OBgh+0JCQmC2JfPpUuXBADCvHnzBEEQhPT0dEGr1QpDhgwRQkJC7OWeffZZoXbt2oLFYhEE4b+/Z/FzLuu4Ul5Hzuzfv1/w9/cXAAjt27cXnnvuOWHFihVCVlZWqbINGzYUAAg//vijfduNGzeEsLAwoUOHDvZt7ri2I0eOdHhtC0Lp652fny+0bt1auO2220rFOXLkyHLP+5577hEACNevXy+3nE379u2F+vXrC1evXrVv27t3r6DVaoURI0bYt02ZMkUAIAwdOrRUHba/35o1axy2T58+XQgMDBSOHj3qsP3ll18WdDqdcObMGfu2kq9HOd4f5s6dKwAQvv76a/u2/Px8ITo6WjAajfb3xso+9zZs2CAAEBYuXChcvnxZuHDhgrBmzRqhSZMmgkajEbZv3+5QXuz1ByD4+fk5vN/v3btXACB8+OGH9m0DBw4UAgIChPPnz9u3HTt2TPDx8XF4rZ06dUrQ6XTCm2++6XCcffv2CT4+PqW2l2T7/NyxY4d928iRIwUADq8HQRCEDh06CJ06dbL/3/aZNnv2bPu2wsJCITY2ttRrx/bcK34uWq1WuOeee0q9Hm3vNYLg/Dn0xBNPCAEBAUJubq59W//+/Uu9JgVBEBYvXixotVph8+bNDtvnz58vABC2bNkiCIIg7NmzRwAgPP300w7lhg0bJupzhojUyaVbquLi4lCvXj1ERETg/vvvR2BgIFauXOnQ06DT6ez3BFssFly7dg2FhYXo3LkzUlJSStX54IMPOnyTZfum88SJEwCAixcvYs+ePRg5ciSCg4Pt5e644w5ERUU51LV69WrodDr7N8c2zz//PARBwK+//uqw/fbbb3f4ZrJbt24ArN+I16hRo9R2W0xlueOOO5CcnIxBgwZh7969mD17NuLj49GgQQP7N6rFRUdHo1OnTvb/33zzzbj77ruxdu1ae3d28fu7CwoKcPXqVTRp0gQ1a9Z0+vcUQ+o1cgfbt2GZmZkArLesaLXWp6HZbMbVq1dhNBrRvHlz0TEU/9tkZWXhypUriImJgSAI2L17t72Mn58fNm7ciOvXrzutZ926dUhPT8fQoUNx5coV+49Op0O3bt1cvtWsXr16aNGihf02uS1btkCn0+HFF19EWloajh07BsDaw9GzZ89KTWlZ0euoLK1atcKePXvwv//9D6dOncL777+PwYMHIyQkBAsWLChVPjw8HPfcc4/9/0FBQRgxYgR2796N1NRUAO65ts4Uv97Xr1/HjRs3EBsb61KdGRkZAODwOi+L7T1o1KhRqF27tn1727Ztcccdd2D16tWl9nnyySed1tW4cWPEx8c7bFu6dCliY2NRq1Yth+dfXFwczGZzqdssi5Pj/WH16tUIDQ3F0KFD7dt8fX3x7LPPwmQy4c8//3Qo7+pzz+aRRx5BvXr1EB4ejrvuugs3btzA4sWL7bfc2Ui5/nFxcYiMjLT/v23btggKCrLHZDab8fvvv2Pw4MEIDw+3l2vSpAn69u3rUNfy5cthsVjwwAMPOFyf0NBQNG3atFK3opZ8nsTGxjr83VavXg0fHx97jwdgff8WMwHJihUrYLFYMHnyZPvr0ab4e03xv2tmZiauXLmC2NhYZGdn4/DhwxUeZ+nSpWjZsiVatGjh8Pex9Zja/j6210nJz2cOQifybi7dUpWUlIRmzZrhxo0bWLhwITZt2gS9Xl+q3Jdffok5c+bg8OHDKCgosG9v3LhxqbI333yzw/9tH1y2xuHp06cBAE2bNi21b8kGzOnTpxEeHl6qEWG7dcBWV1nHtiU0ERERTreX1WAtrkuXLli+fDny8/Oxd+9e/PTTT3jvvfdw//33Y8+ePQ5JkrNzatasGbKzs3H58mWEhoYiJycHM2fOxBdffIHz58873J5x48aNCuMpi5Rr5A4mkwnAfw08i8WC999/H/PmzcPJkycd7he23fJQkTNnzmDy5MlYuXJlqWtj+9vo9XrMmjULzz//PEJCQtC9e3cMGDAAI0aMsM+AY2v4l3VLUVBQkIQzdRQbG2v/oN28eTM6d+6Mzp07o3bt2ti8eTNCQkKwd+9eDBs2zOVjABW/jsrTrFkzLF68GGazGQcPHsQvv/yC2bNn4/HHH0fjxo0RFxdnL9ukSZNSiVGzZs0AWG+FCQ0Ndcu1deaXX37BG2+8gT179jjcG+5Koma7ppmZmRXeDmN732jevHmpx1q2bIm1a9eWGhhe1uvI2fZjx47hn3/+cXorD1D+4H053h9Onz6Npk2blmqkin0flfLcA4DJkycjNjYWJpMJP/30E7777rtSxwakXf+SMdnissV06dIl5OTkoEmTJqXKldx27NgxCILg9P0acH2CA4PBUOqaF48RsP6tw8LC7F/Y2Dh7LpZ0/PhxaLXaUl/MlXTgwAFMmjQJf/zxhz0RtxHzHDp27BgOHTpU4fP39OnT0Gq1DokgIO5ciEi9XEo4unbtap+lavDgwejZsyeGDRuGI0eO2N8Qv/76a4waNQqDBw/Giy++iPr160On02HmzJmlBgEC1m9rnBFKDPKWQ1nHdkdMfn5+6NKlC7p06YJmzZph9OjRWLp0qeS1DZ555hl88cUXGDduHKKjoxEcHAyNRoOHHnrI5WlDpV4jd9i/fz/q169vb+jNmDEDr732Gh555BFMnz4dtWvXhlarxbhx40Sdl9lsxh133IFr167hpZdeQosWLRAYGIjz589j1KhRDnWMGzcOAwcOxIoVK7B27Vq89tprmDlzJv744w906NDBXnbx4sVOp+GszMwvPXv2xIIFC3DixAls3rwZsbGx0Gg06NmzJzZv3ozw8HBYLBb7t8KucsdzVqfToU2bNmjTpg2io6PRp08ffPPNNw4JhxiVvbbObN68GYMGDcKtt96KefPmISwsDL6+vvjiiy+wZMkSyfW1aNECgHXdnMr+7Z0pa+YpZ9stFgvuuOOOMmdksiV0zsjx/iBVZZ97bdq0sT/HBg8ejOzsbIwZMwY9e/a0f/kj9fq783PFYrFAo9Hg119/dVpvyWRArLJirErp6eno1asXgoKCMG3aNPtaPCkpKXjppZdEPYcsFgvatGmDd9991+njJb/AI6LqpdKDxm0N1D59+uCjjz6yz/+9bNky3HLLLVi+fLnDN0+uLiJmmz3H9i10cUeOHClV9vfff0dmZqZDL4etW9jZTDxVwZakXbx40WG7s3M6evQoAgIC7N8WLVu2DCNHjsScOXPsZXJzc0vNUiLlW153X6OKJCcn4/jx4w5T5i5btgx9+vTB559/7lA2PT0ddevWtf+/rPPat28fjh49ii+//BIjRoywby8+ML24yMhIPP/883j++edx7NgxtG/fHnPmzMHXX39t/8atfv36FTaupX6bbmvMrlu3Djt27LC/Tm699VZ8/PHHCA8PR2BgoMOtde44bmWV9Zz9999/IQiCQzy2BfJstydW9to68+OPP8JgMGDt2rUOvapffPGF6DqKGzhwIGbOnImvv/66woTD9r5R8v0GsL631K1b1+m0t2JFRkbCZDJJTuwAed4fGjZsiH/++QcWi8Whp6Gq3kffeust/PTTT3jzzTcxf/58AO6//vXr14fBYMC///5b6rGS2yIjIyEIAho3blxu8ieHhg0bYv369TCZTA6JjbPnYkmRkZGwWCw4ePBgmROBbNy4EVevXsXy5csdJq04efJkqbJlPYciIyOxd+9e3H777eU+zxo2bAiLxYLjx4879GqIORciUi+3TIvbu3dvdO3aFXPnzkVubi6A/761Kf5N0rZt25CcnOzSMcLCwtC+fXt8+eWXDt2769atw8GDBx3K9uvXD2azGR999JHD9vfeew8ajabUvbnutmHDBqffoNluqSnZdZycnOxwS9jZs2fx888/484777T/HXU6Xak6P/zww1JTJ9oaPM6mSyzJ3deoPKdPn8aoUaPg5+eHF1980SGGkue1dOnSUivXl3Vezs5BEAS8//77DuWys7Ptz02byMhI1KhRw35bRnx8PIKCgjBjxgyH28tsLl++XGE8ZWncuDEaNGiA9957DwUFBejRowcAayJy/PhxLFu2DN27d6+wF0XqccXavHmz03Mu6zl74cIFhympMzIy8NVXX6F9+/b23qHKXltndDodNBqNw/P+1KlTDjPgSBEdHY277roLn332mdM68vPz8cILLwBwfA8qHuv+/fvx22+/oV+/fi7FYPPAAw8gOTkZa9euLfVYeno6CgsLy9xXjveHfv36ITU1Fd9//719W2FhIT788EMYjUb06tWrwjoqIzIyEvfddx8WLVpkHxfk7uuv0+kQFxeHFStW4MKFC/bt//77b6mxfvfeey90Oh2mTp1a6m8tCAKuXr3qUgxi9OvXD4WFhfYpswFr7+6HH35Y4b6DBw+GVqvFtGnTSvVU2M7D2ftofn4+5s2bV6q+wMBAp7dYPfDAAzh//rzTMV85OTnIysoCAPvnb/EpwQFw9XIiL+e21YFefPFFDBkyBIsWLcKTTz6JAQMGYPny5bjnnnvQv39/nDx5EvPnz0dUVJT9Pn6pZs6cif79+6Nnz5545JFHcO3aNXz44Ydo1aqVQ50DBw5Enz598Oqrr+LUqVNo164dfvvtN/z8888YN25cqXtH3e2ZZ55BdnY27rnnHrRo0QL5+fnYunUrvv/+ezRq1AijR492KN+6dWvEx8c7TIsLwGF15wEDBmDx4sUIDg5GVFQUkpOT8fvvv5e6F759+/bQ6XSYNWsWbty4Ab1ej9tuu63UWgq2Ot19jQAgJSUFX3/9NSwWC9LT07Fjxw78+OOP0Gg0WLx4Mdq2besQw7Rp0zB69GjExMRg3759+Oabb3DLLbc41BkZGYmaNWti/vz5qFGjBgIDA9GtWze0aNECkZGReOGFF3D+/HkEBQXhxx9/LHXf+NGjR3H77bfjgQceQFRUFHx8fPDTTz8hLS3NvvJ5UFAQPv74YwwfPhwdO3bEQw89hHr16uHMmTNYtWoVevToYU9ibT0Rzz77LOLj46HT6SpcQT02Nhbfffcd2rRpY7+/vWPHjggMDMTRo0dFjd9w5bhizJo1C7t27cK9995rvz4pKSn46quvULt27VIDOps1a4ZHH30UO3bsQEhICBYuXIi0tDSHb5ore22djXPo378/3n33Xdx1110YNmwYLl26hKSkJDRp0gT//POPS+f+1Vdf4c4778S9996LgQMH4vbbb0dgYCCOHTuG7777DhcvXrSvxfH222+jb9++iI6OxqOPPmqfFjc4OLjSq7+/+OKLWLlyJQYMGIBRo0ahU6dOyMrKwr59+7Bs2TKcOnXKoWeoODneHx5//HF88sknGDVqFHbt2oVGjRph2bJl2LJlC+bOnStqoH1lvfjii/jhhx8wd+5cvPXWW7Jc/9dffx2//fYbevTogaeeesr+ZVXr1q2xZ88ee7nIyEi88cYbmDhxIk6dOoXBgwejRo0aOHnyJH766Sc8/vjj9uTU3QYOHIgePXrg5ZdfxqlTp+zr3ogZW9GkSRO8+uqrmD59OmJjY3HvvfdCr9djx44dCA8Px8yZMxETE4NatWph5MiRePbZZ+3v1c6+OOvUqRO+//57JCYmokuXLjAajRg4cCCGDx+OH374AU8++SQ2bNiAHj16wGw24/Dhw/jhhx/s6860b98eQ4cOxbx583Djxg3ExMRg/fr1TnuZiMiLSJnSytm0fjZms1mIjIwUIiMjhcLCQsFisQgzZswQGjZsKOj1eqFDhw7CL7/8Umqay/KmUoWTKfJ+/PFHoWXLloJerxeioqKE5cuXO506MzMzUxg/frwQHh4u+Pr6Ck2bNhXefvtth2kAbcdISEhw2FZWTLbpG5cuXVru3+nXX38VHnnkEaFFixaC0WgU/Pz8hCZNmgjPPPOMkJaW5vT4X3/9tdC0aVP736rktJXXr18XRo8eLdStW1cwGo1CfHy8cPjwYadTfy5YsEC45ZZbBJ1O5zAFZslpL8VeI1ucYqfFtf34+PgItWvXFrp16yZMnDjRPkVtcbm5ucLzzz8vhIWFCf7+/kKPHj2E5OTkUrEKgnWazaioKPt0lbapIA8ePCjExcUJRqNRqFu3rjBmzBj79Je2MleuXBESEhKEFi1aCIGBgUJwcLDQrVs34YcffigV04YNG4T4+HghODhYMBgMQmRkpDBq1CiHqYsLCwuFZ555RqhXr56g0WhETZGblJQkABCeeuoph+1xcXECAGH9+vVO/57Fp7ws67hSX0clbdmyRUhISBBat24tBAcHC76+vsLNN98sjBo1Sjh+/LhD2YYNGwr9+/cX1q5dK7Rt21bQ6/VCixYtSr0u3HFtnT0XP//8c/trpUWLFsIXX3xRaipQW5wVTYtrk52dLbzzzjtCly5d7K/Zpk2bCs8884zDtKqCIAi///670KNHD8Hf318ICgoSBg4cKBw8eNChjC2ey5cvlzqW7e/nTGZmpjBx4kShSZMmgp+fn1C3bl0hJiZGeOeddxym1S55TeV4fxAEQUhLS7PX6+fnJ7Rp06bU9MWVfe5V9L7au3dvISgoyD4dutjr7+y9XRCcPy/Wr18vdOjQQfDz8xMiIyOFzz77THj++ecFg8FQav8ff/xR6NmzpxAYGCgEBgYKLVq0EBISEoQjR46Ue55lTYsbGBhYqqyz87l69aowfPhwISgoSAgODhaGDx8u7N69u8JpcW0WLlwodOjQQdDr9UKtWrWEXr16CevWrbM/vmXLFqF79+6Cv7+/EB4ebp/OvfhzRBAEwWQyCcOGDRNq1qxZakr6/Px8YdasWUKrVq3sx+nUqZMwdepU+1TogiAIOTk5wrPPPivUqVNHCAwMFAYOHCicPXuW0+ISeTGNIFTBqGwqk0ajQUJCQqnbv4iUqlGjRmjdujV++eUXT4dCJJvBgwfjwIEDTsfYERGRNG4Zw0FERKRWOTk5Dv8/duwYVq9ejd69e3smICIiL+O2MRxERERqdMstt2DUqFG45ZZbcPr0aXz88cfw8/Mrc4piIiKShgkHERFVa3fddRe+/fZbpKamQq/XIzo6GjNmzChzkT8iIpKGYzhUYNOmTXj77bexa9cuXLx4ET/99BMGDx5c7j4bN25EYmIiDhw4gIiICEyaNAmjRo2qkniJiIiIiGw4hkMFsrKy0K5dOyQlJYkqf/LkSfTv3x99+vTBnj17MG7cODz22GNO5/cnIiIiIpITezhURqPRVNjD8dJLL2HVqlXYv3+/fdtDDz2E9PR0rFmzpgqiJCIiIiKy4hgOL5ScnIy4uDiHbfHx8aUWbysuLy/PvuI2YF1R+NChQ4iIiIBWy44wIvIMQRBgMpkQFhbG9yIiIpkJgoDMzEyEh4e79T2XCYcXSk1NRUhIiMO2kJAQZGRkICcnB/7+/g6PJSUl4fXXX8eVK1eqMkwiIiIiUqCzZ8/ipptuclt9TDgICQkJeOyxxxx6OM6ePYvWrVvj7NmzCAoKqrAOk8kEADAajaKOmZ6eji1btqBHjx6oWbOm2+t3ZR+lnYPa43dlH6WdQ3WL35V95D6HixcvokWLFjhz5gyCg4NF7UNERK7JyMhAREQEatSo4dZ6mXB4odDQUKSlpTlsS0tLQ1BQUKneDRu9Xg+9Xm//v+2DPSgoSFTCYet2E9uIsFgsCAgIkK1+V/ZR2jmoPX5X9lHaOVS3+F3ZR+5zsCUowcHBouonIqLK02g0bq2PN8R6oejoaKxfv95h27p16xAdHe2hiIiIiIioumLCoQImkwl79uzBnj17AFinvd2zZw/OnDkDAJg4cSJGjBhhL//kk0/ixIkTmDBhAg4fPox58+bhhx9+wPjx4z0RPhERERFVY0w4VGDnzp3o0KEDOnToAABITExEhw4dMHnyZADWe5xtyQcANG7cGKtWrcK6devQrl07zJkzB5999hni4+M9Ej8RERERVV8cw6ECvXv3RnnLpSxatMjpPrt375YxKiIiIiKiirGHg4iIiIiIZMOEg4iIiIiIZMOEg4iIiIiIZMOEg4iIiIiIZMOEg4iIiIiIZMOEg4iIiIiIZMOEg4iIiIiIZMOEg4iIiIiIZMOEg4iIiIiIZMOEg4iIiIiIZMOEg4iIiIiIZMOEg4iIiIiIZMOEg4iIiIiIZMOEg4iIiIiIZMOEg4iIiIiIZMOEg4iIiIiIZMOEg4iIiIiIZMOEg4iIiIiIZMOEg4iIiIiIZOPj6QCIiKiIxQKYzYAgWP+v1QI6HaDReDYuIiKiSmDCQUTkSYIA5OUBJtN/iUZJOh2g11dtXERERG7CW6qIiDwlPx+63FwgP7/sZAOw9npkZ1t/LJaqi4+IiMgN2MNB5TKZTNBqK85Ls7OzJdWbk5Nj/+3jU/HTUGr9ruyjtHNQe/yu7KO0c5CtvCAAubnIvXEDGgC5ubkwmUwV7paTkwNkZVmTDoVcM7mvga1+IiJSLyYchKSkJCQlJTlsKygo8FA0RF5OECrXUyEIQE4OYDAAvr7ujY2IiEgGTDgICQkJSEhIcNh27tw5REREwGg0wmg0iq5LbNnCwkIAgL+/vyz1V2YfpZ2D2uN3ZR+lnYPbytuSjYAAAP/FbzAYXItfrxeVdMh5zeS+BhkZGaLrJCIiZeIYDiKiqpKT494xGLm5QFGDn4iISKmYcBARVYXcXOvgbznq5UByIiJSMCYcRERyKyiw/sihaAA6ERGRUjHhICKSk22dDTmZzdapdYmIiBSIg8aJiOSUl1fmGhsWiwU3MjJw5epVFBSNxdAAeG/ePAT4+yNx7Fj4+/uLO05+vnUAOVclJyIihWHCQUQkF7MZQn6+fe2JgIAAaIolBNeuX8fDTzzhsEvXTp2wfdcuAMCRY8fw1aefijuWrSfFYHBP7ERERG7CW6qIiOSSl4fs7GwYw8JgDAsrc9E7fwAdin7bkg0AWPzdd/hn/37xxysokGdgOhERUSUw4SAikoOExn8LAClFvwHgpfHjMeSeewAAn335pbTjyj1ehIiISCImHEREcqjEIO7R//sfRg0bBgBY9vPPsEiZ9tZs5tocRESkKEw4iIjcrbDQ5bUxmkZGonnTpri9d28YjUZcTE3F3n37pFXCGauIiEhBmHAQEblbJRr8t/XqBQDQ6/W47dZbAQC//fGHtErMZi4GSEREisGEg4jInSyWSg3cvjUmxv5vW8KxYdMm6RWxl4OIiBSCCQcRkTtVckXxbp072//dOzYWALBl2zYUSh2XUVBQ5vofREREVYkJBxGRO1Uy4QgKCrL/u3VUFIKCgmAymbDvwAHplXHwOBERKQATDiIidyksdGuvgk6nQ/eiHo+t27ZJr6CSyQ8REZE7MOEgInIXGRr4Md26AQCSt2+XvjMHjxMRkQIw4SAicgdBkOUWpu5dugAA/t6xw7UK2MtBREQexoSDiMgdZBov0bVTJwDA8ZMncfnqVekVVGLGLCIiIndgwkFE5A4yNexr1aqFFs2aAQB27d4tvQKzmbNVERGRRzHhICJyBxlnhLLdVrV91y7XKmAvBxEReRATDiKiyrJYZO1FiO7aFUAlEg5Oj0tERB7EhIOIqLJcbNCnXbokqpwt4diZkiJ9AUCAPRxERORRTDiIiCrLxQb9rj17RJVr1bIlgoODkZWdjf0HD0o/kMXC6XGJiMhjmHAQEVWWiwnHTpGDwLVaLWKKejm2urIeB8BeDiIi8hgmHB62adMmDBw4EOHh4dBoNFixYoXD44IgYPLkyQgLC4O/vz/i4uJw7NgxhzLXrl3Dww8/jKCgINSsWROPPvooTCZTFZ4FUTVWiVmgpIzJuLVHDwDAX8nJLh2L4ziIiMhTmHB4WFZWFtq1a4ekpCSnj8+ePRsffPAB5s+fj23btiEwMBDx8fHIzc21l3n44Ydx4MABrFu3Dr/88gs2bdqExx9/vKpOgah6c7Ehf/nKFRwp8eVBeXr37AnAmnBYXLk9ij0cRETkIT6eDqC669u3L/r27ev0MUEQMHfuXEyaNAl33303AOCrr75CSEgIVqxYgYceegiHDh3CmjVrsGPHDnTu3BkA8OGHH6Jfv3545513EB4eXmXnQlQtudiQX79xo6TynTt2RA2jEdfT07Fr9250KVoQUDRBsI7j0PJ7JiIiqlpMOBTs5MmTSE1NRVxcnH1bcHAwunXrhuTkZDz00ENITk5GzZo17ckGAMTFxUGr1WLbtm245557StWbl5eHvLy8co+dmZkJAEhPTxf1bWp2djYAiJ5B59Zbb8WFCxfg6+sLrYgGkC0GMWVd3ceV8gUFBbKdg9rjdzUmqeU9fQ0CLBZoyilf/PXTrl07+2Op6elO6+vWvTt8yjieJT8fAHDXwIGoYzRKjj9Po0GhRlPeLi5fg/DwcGzatKnC8lLfK2zvRUREpF5MOBQsNTUVABASEuKwPSQkxP5Yamoq6tev7/C4j48PateubS9T0syZMzF16lRRMWzZsgUBAQFSQ6/QhQsXcPXqVbfXS1TVKmr2Fx/dkXrpEjRF27LKKH+pqIwztiZ6usmEfBfGaeUX/cjlzz//dHudtgSFiIjUiwlHNTRx4kQkJiaWW+b8+fOIiopCjx49EBQUVGGdtkaB2OTE19cXgPVb1NDQ0ArL89t1ecqzh6Ny5bWCAP9iA8bL6uE4XrTeRmj9+tBqtcjKy0PW9ev25KO4+vXrl9nDkV9YiNNXrsACoE7duvDzKf8tvGQ8Zo0GuW7u4UhNTYXFYoGvry969epVYXmp7xVpaWmiyhERkXIx4VAwW0M8LS0NYWFh9u1paWlo3769vcylEouHFRYW4tq1a2U25PV6PfR6fbnHzsjIAADUrFlTVMLhU9TwMYq4zQP4rzETGhqK8+fPV1jeNuuW2Ppd2Udq+fT0dPz555/o1asXatas6fF4lBa/K/so7RwqLF9QABSbwMFZ+aysLBiLXr979+5FYGAghj3yCL5dtgwPDRmCb5cudahy299/o17dumXGM3joUKz/80+MeOQRTJs0SVr8Wi0QGChtnwo0aNAAFy5cgFarFXUNpL5XsIeDiEj9OHpQwRo3bozQ0FCsX7/evi0jIwPbtm1DdHQ0ACA6Ohrp6enYVWx6zT/++AMWiwXdunWr8piJqhUXBoxfuHgRy37+GQAwpGgyCCmGP/QQAODTRYscZqsThYv/ERGRBzDh8DCTyYQ9e/ZgT9GKwydPnsSePXtw5swZaDQajBs3Dm+88QZWrlyJffv2YcSIEQgPD8fgwYMBAC1btsRdd92FMWPGYPv27diyZQvGjh2Lhx56iDNUEcnNhQb8m2+/jYKCAvTo3h1tWrWSvP+gfv0QcdNNSLt0CR9/9pnk/Tk9LhERVTUmHB62c+dOdOjQAR06dAAAJCYmokOHDpg8eTIAYMKECXjmmWfw+OOPo0uXLjCZTFizZg0MBoO9jm+++QYtWrTA7bffjn79+qFnz5749NNPPXI+RNWKxIRj09atmL9wIQBgegW3Q5XF19cXk196CQAwecYMSWt5AGAvBxERVTmO4fCw3r17QyhnlWKNRoNp06Zh2rRpZZapXbs2lixZIkd4RFQWQZC8wviDo0bBYrFg5LBh6HPrrfj3+HGXDv3I8OFY/N132LRlC+IGDcJXn3yC3rGx0FQwIBwAEw4iIqpy7OEgInKFyIb7Z19+af93ZmYmekZH4+3p05GVlYXsnJxS5bOys5GVleX8JzsbWdnZyMnJwaL589GsSROcO38etw0YgNbduuHGjRtui5uIiMhd2MNBROQKkQ33kj2YfyUno35kZJnlG7du7VI4Bw8fts8AVS4mHEREVMXYw0FE5AqRDff4uDiZA5GICQcREVUx9nAQEblC5PiN5k2bIq1orEaAv7/DOIvjJ0+iXUwMDgPoCOAwgJP796NenTpO6zJlWdcnN5axlobYxfQgCICY8R5ERERuwISDiNTBbHac0lWrlTxo261EHlur1aJ+vXpOH7spPBzffPIJunbtiuDgYABAndq1y1zl23Z7VmAFi/dViAkHERFVISYcRKRcFguQnw8UFjpv4JtM0ObleWZtCTckO1qtFsFBQahbp46oVbrdxpOJGhERVTtMOIhImfLyrMlGBbQWC5CdDQQEAAZD1X1zr+ZGu5pjJyIi1eGgcSJSFkGwJhAikg0HhYVAVlbV9XaoudGu5tiJiEh1mHAQkXJYLJVLGmzJSmGhe+Mq61hqpebYiYhIdZhwEJEy2G6NckdjOCdH3qRD7Q12tcdPRESqwoSDiDxPEKxJgjsbwjk58t1epfYGu9rjJyIiVWHCQUSel5srz4J0ubnyNK7VvngeEw4iIqpCTDiIyLNs097KwWKxJh3upvYGu9rjJyIiVWHCQUSeY7FYp7+VU2EhUFDg3jrV3mBXe/xERKQqXIeDiDxHZLJhsVhw9do1XLh4EfMXLkRI/fp45fnnpR3Hh293dkw4iIioCvETmIg8o7AQQkEBsrOzAQABAQFlFr167Rrq33KLw7ZLly9jxuTJ4o4lCNLX9fAQs9mMzVu34mJqKsJCQxEbEwOdTufpsIiIiFzGW6qIyDPy8pCdnQ1jWBiMYWH2xKM8/gA6FP3+9IsvcP7iRfHHy89X/GDv5StXokmrVujTvz+GPfoo+vTvjyatWmH5ypWeDo2IiMhlTDiIqOq52PhvASCl6LfZbMa3y5ZJq0Du8SKVsHzlStw/fDjaXLiAZACZAJIBtLl4EfcPH86kg4iIVIsJBxFVLRdubzpx8qTT7T+vWiXt2IWF7lmbw81jIMxmM56fMAEDBAErAHQHYCz6vUIQMADACy+9BLNc64oQERHJiAkHEVWtwkLJDfb1Gzc63b7/4EFcuXZN2vEVOJYjeft2nLpwAa+g9JuyFsBEQcDJ8+exdds29xyQg8aJiKgKMeEgoqrlQoN/45YtpbZ1aNsWALBn3z5plbmQ8Mgt9dIlAEDrMh63bU9NS3PPATUa99RDREQkAhMOIqo6ZrPksRtmsxnJ27eX2n5rjx4AgH0HD0qPo7K9HG5usIfWrw8A2F/G47btoSEhbj0uERFRVeC0uCqRlJSEt99+G6mpqWjXrh0+/PBDdO3atczyc+fOxccff4wzZ86gbt26uP/++zFz5kwYDAZJxzWZTNBqK85LxcwwVJxQ9A2zIAgwmUxur9+VfaSWz8nJsf/2EbHGg9zxKC1+p/vk5DisKp5V7HFTVha0Thrye/75x+lzpEO7dgCA/YcPIzc3V9TzyHbOuH4dMBqlx2+Tn+90ALq9fpFyi1ZB79C2LRqGhWFGaipWCILDN0EWADM1GjQKC0OHtm0lHaPMsuUkTEp7LUv9mxIRkfKwh0MFvv/+eyQmJmLKlClISUlBu3btEB8fj0tFt2GUtGTJErz88suYMmUKDh06hM8//xzff/89XnnlFaflk5KSEBUV5fDTp08fOU+JqiNBcEg2xPp7xw6n2zu2awetVou0S5ek32rkYixy0el0eHP6dPwCYLBG4zBL1WCNBr8AeGPaNK7HQUREqsQeDhV49913MWbMGIwePRoAMH/+fKxatQoLFy7Eyy+/XKr81q1b0aNHDwwbNgwA0KhRIwwdOhTbyhhwmpCQgISEBIdt586dQ0REBIxGI4wivgm2EVtWU/QNq0ajkaX+yuwjtnxhUYPV399f1nNQe/z2fZzcxqQp9k27MTDQ/m158WPs2rPHaZ2hISGIatEC+w8exD8HDqB9UY+H6Hh8fAB/f/Hli8vPB3x9xZcvg+0aGAwGPPzAA/A3GPD8hAmIuXDBXqZxeDiWzZqFewcNsvcgVPo5odFU2MOjlNdyRkaG6DqJiEiZ2MOhcPn5+di1axfi4uLs27RaLeLi4pCcnOx0n5iYGOzatQvbi+57P3HiBFavXo1+/fqVeZy8vDxkZGTYfzIzM917IkQFBS7tllxGDwcAdOnYEQCwfdcu6RUrcPD4vYMG4d8DB7Bh1Sos+fxzbFi1Csf278e9gwa590AcNE5ERFWIPRwKd+XKFZjNZoSUGCwaEhKCw4cPO91n2LBhuHLlCnr27AlBEFBYWIgnn3yyzFuqAGDmzJmYOnWqW2MnsrNYXFro72JqKk6fOVPm4107dsQXX3+Nnbt3uxZXQQHg5yd9Pxkb7DqdDr1jY2Wrn4iIqKqxh8MLbdy4ETNmzMC8efOQkpKC5cuXY9WqVZg+fXqZ+0ycOBE3btyw/xx0ZeYforK4OF5i286dAIAWzZo5fbxzUQ/H3n37kO/KzFOuLqSn9h4CERNBEBERuQt7OBSubt260Ol0SCsxKDYtLQ2hoaFO93nttdcwfPhwPPbYYwCANm3aICsrC48//jheffVVp7NO6fV66PV6+/953zS5lYsNe9uA8c4dO+Lw0aOlHo9s3BjGwECYsrLwz/799gRE7rhU32BXe8JERESqovJPTe/n5+eHTp06Yf369fZtFosF69evR3R0tNN9srOzSyUVttltBIXds07VRCUTjk5lDAjXaDRo3qQJADhdq6NCguBabGpvsKs9fiIiUhUmHCqQmJiIBQsW4Msvv8ShQ4fw1FNPISsryz5r1YgRIzBx4kR7+YEDB+Ljjz/Gd999h5MnT2LdunV47bXXMHDgQE6rSVXPbHZpcHZhYSF2pKQAQLk9Fy2aNgVQ9vS5Ig4kfR+1N9jVHj8REakKb6lSgQcffBCXL1/G5MmTkZqaivbt22PNmjX2geRnzpxx6NGYNGkSNBoNJk2ahPPnz6NevXoYOHAg3nzzTU+dAlVnhYUuDczeu28fsrOzUbNmTTSNjCyzXMui8R1bXenhACo3jkOtPYZMOIiIqAox4VCJsWPHYuzYsU4f27hxo8P/fXx8MGXKFEyZMqUKIiOqgIsN+i1//w0AiO7SpdzV7ptFRkKr1eLU6dO4cPEiwsPCpMcnCNIb4Uw4iIiIROEtVUQlCYJ1utTcXCAr67+f3FzrdrU2Mj3B1TESADYXrTPTo3v3cssFBASgVcuW1n22bnXpWNVuHIeaYyciItVhwkFkIwhAXp5jcmFbP8JicUxC8vKYeIjhYrIhCAI2bdkCAOjVs2eF5Xt06wYA+POvv1w6XrUbx6Hm2ImISHWYcBAB1oZxVhaQn19xIiEI1nJZWdC4ev9/deHi+hsHDh3CpcuX4e/vb19NvDw9i2Zs+2PTJpeOxx4OIiIi+TDhIMrPB7KzpfdYCAJ0+fnW3g5yzsWEbN0ffwAAYqOjHdaHKUvP7t2h1Wpx5NgxnDl7VvoBLRbp11/NjXY1x05ERKrDhIOqt7y8yicMTDrKZrG4tNuv69YBAOLj4kSVDw4ORrfOnQEAa37/3aVjSk6O1Lr4n1rjJiIi1eInD1Vf+fnWH3fVxaTDkYvJxo2MDGwsGovRPz5e9H62sj+vWuXScSXHq9aGu1rjJiIi1eInD1VPhYXuTxDy860Dy8nKxYTjlzVrUFBQgJbNm6N50aJ+YtwzcCAAYN2GDbh27Zr0AzPhICIikgU/eaj6sViss03JIS/P5Ya213Fx/MaSpUsBAEPvv1/SflEtWqBt69YoKCjAt8uWST+w1Oum0ahzLAQTDiIiqmL85KHqJzdXviltBQHIyZGnbrVxIfE6cOgQNm3ZAq1Wi5HDhkne/5H//Q8AkLRgASxSj+9KoqjGxrsaYyYiIlXjJw9VL/n59m/eLRYLLl+58t/P1avSG6nOWCzuGxuiZi78Ld985x0AwL2DBuHmiAjJ+4/+3/8QHByMQ0eO4NuinhLRBEF6IqrGxrsaYyYiIlXjJw9VG4LFgqzr15GVlQVBEHD12jXUv+UW+88tbdrg6vXrDvvMTUpCrZtvRs8778Sly5fFH0zMeh7eTmLC8ePKlfi/X3+FTqfD6xMnunTIoKAgTHjuOQBA4iuv4MLFi9IqkHobmE4nrbynabXqvA2MiIhUjQkHVRvZ16/DGBoKY1gYsrOz7dv9AXQo+l3cth07MH7iRKSnp2PL339j7AsviD+YIMg3TkQNXOjdeGzsWADAuKefRquWLV0+9PPPPIM2rVrh0uXL6N2vH7bt2CF+Z28fOK62eImIyCvw04eqB7O5zBmkWgBIKfpd3Oz33wcAtGzeHBqNBkt/+gnHT5wQf8zCQpcHTqueyIb7rPfes/+7oKAAA+66C+MTEpCVleX4UyxBtMnKzkZWdjZyc3Ot/y4qW1hYiG8//xw3NWiAY8ePo/vtt6Nr794oEDODGBMOIiIit/PxdABEVULiFLjXr1/HytWrAQA/fPklXnj1Vaxdvx5ff/89pki53Sc3FwgMlHRsryAy0TKXKPfLmjW4qUXJ1M+5xq1biw5nR0oK8sWMq3F1piq13D7HhIOIiDyAnz7k/cxmyT0Nq9auRWFhIVpHRaF1VJR9itaffvlF2rEtFmtPR3UjsgF+/+DB8sYhlbfPVKWmWImIyGuwh4O8nwszRv26bh0AYGDfvgCAfvHx0Gg02LtvHy5cvIjwsDDxlRUUAD7V7KUmMuFo16YN0o4fBwAE+Pvbb50ylugVunz1Khq3bo3DADoCOAzg5P798PXxwd/JyegeHY2awcHlHisgIABZWVluiduBTqeOW+c0GvUNciciIq9QzVpBVO0IguQeBkEQ8MemTQCAO2+7DQBQr25ddGzfHrt278aGTZvw8IMPiq+wsLD6LQYosuGu1WpRv169UtsDSyQc/v7+uFRi/Eyd2rWRkZEBg8GAwICAUvu4TBCkzeTk46OOaZCZbBARkYewf528mwsNwaP//ovUtDQYDAZ079LFvr1PbCwAYMPmzdLjEDNg2Zu4eUyDVqtFvbp1HX60ct0eJDV2nU4dU81Wt142IiJSDCYc5N1caOhv3roVANC1UycYDAb79lt79AAA/JWcXCVxqJqae3Rcva1K6dQQIxEReSV+5UXlMplMor5JznYybWl5hKJGnSAIMJlMbq8fALIzMqDJzbU3IItPrWrKykJ2Tk6pfXJzcrCxqAeja6dODrG1bdUKAHDk2DGcOn0aBr3euk9urrhzEARJ3zJLPWep5XOKzj8nJwc+IuISXb8gAEV/jxwnf2MxMYmVW7TWidhrIKp+s9l+nUSfc34+kJenjPidlZew4J/SXstSz5mIiJSHCQchKSkJSUlJDttErVmgdC6ew/ZduwAA3Tp3dthep3ZtNGvSBEf//Rc7UlIQGx0tPZ7qcFuLWqaILYsr8fv4SJ56uUqxd4OIiDyoGrR+qCIJCQlISEhw2Hbu3DlERETAaDTCaDSKrktsWU3Rt60ajUaW+iEIQGamwz6aYt/wGgMDnfbc5Obl4VjRrEl9YmNLHS+mWzcc/fdf7Nm3D3f06QMAMBgM4uIymQB/f0nnWzx+d5cvLBpM7y8xpgrLms2lvk2X+xxEXwMx9ev1gJ+f9HiKPZ88Gr+z8jVrSk52lfJazsjIEF0nEREpE8dwkHdyce2LPfv2AQCaNWmCOnXqlHo8umtXAMC2nTtdi8sbeo4qUh17OABl916xh4OIiDyICQd5JxcTjr3//AOg9O1UNrZZq7bt3FlqlWxR1LBeQ2VV14RDqY16tcyiRUREXosJB3knFxv2u4sSjuLT4RbXqmVLBAYGwmQy4fCxY1UWl6pU14RDqT0cSk2EiIio2mDCQd7HYnG50VhRwqHT6dC1UycAwM6UFOkHEATvTzqqa8Kh0Sgz6fD19XQERERUzTHhIO/j4u1UAJBpMiEgIABtW7cus4xtHIdtNivJKhEfVYHKJExKa9zrdA6D2YmIiDyBn0TkfSrZoO/aqVO561JEF/V+uJxweHsPh4zMZjM2bt6Mb5cuxcbNm10bRyMnHx9ljZdQWgJERETVkgL7/4kqqZKN0B7dupX7eEzR4/+eOIEbrkzZaTZbv0VXUsNUBZavXInnJ0zAqQsX7NtuDgvD/x56CD2KVoFXBCXdVqWkWIiIqNpiDwd5FzfcrtSzggX9ateujVYtWwIADhw54tpBlPbNvDvJMIZj+cqVuH/4cLS5cAHJADIBJANol5qKme+9h//79Ve3H9NlJdbw8BhfXya1RESkCEw4yLtUsiGv1WrRo3v3CsvdGhMDANh38KBrB+I4DtHMZjOenzABAwQBKwB0B2As+r1CENAfwGuvv66c26u0WmXMDKWUxIeIiKo9JhzkXSrZkG/bujVq1KhRYbnbevUCAOzdv9+1AymlcawCm7duxakLF/AKSr9haQG8AuD0xYvYvHVr1QdXFk839jlYnIiIFISfSORdLJZK7V7R+A2b2269FVqtFmfOncP5ixelH6gSU/dWNxdTUwEAZc0b1rpEOUXw8fFsg9/TCQ8REVExTDjIe1Qy2QCAXj17iipXu3ZtdGrfHgDw2/r1rh3MDfFWB2GhoQCAsvqS9pcopxieavRrtRwsTkREisKEg7yHGxrwHdu1E132rrg4AMCqtWtdO5i3JhxuHqgcGxODRuHhmKHRoORfzAJgBoCGYWGILRpXoxi+vp4Zy2EwVP0xiYiIysGEg7yHG8ZFaCU0EAf16wcA2PjXX7h0+bL0g3lrwuFmOp0Oc2bPxi8ABms0DrNUDdZosArA9Ndfh04JA7VL0uur9ng+PsoYsE5ERFQMEw7yHlXcgG9yyy1oGhkJs9mMr779VnoFTDhEu3fQICxbvBj7wsIQAyAIQAyAf0JDMXH8eAzs29fDEZZBp6u6xfc0GvZuEBGRIjHhIO/hgQZ839tvBwB8MH8+8vLypO3MmaokuXfQIPx74AA2rFqFJZ9/jg2rVmHX1q2I6drVvQdy99oVen3VrIfh58d1N4iISJGYcJD3kJhwXLl6tdKH7N2jB8JCQnD23Dm8//HH0nYWBO+cqUrGRq9Op0Pv2FgMHTIEvWNj5bmNyt3xazTy31ql03FmKiIiUiwmHOQdJCYbubm5eHLcuEof1s/PD5MmTAAAvPbGG9i0ZYu0Crzxtiq1f8suR/y+vvLdWsVbqYiISOGYcJB3EHF7klCsN+HWvn2dJgfZ2dnIysoS95OdjdzcXAzq1w8D+/ZFfn4+4gYNwguTJmH/gQMOxysTEw7lkSt+g0GetTnkqpeIiMhNOFk7eQcRDffsnBz7vw8ePuy0TGuRC/+VpaCgAHM++ABzPvgAposXERgYWP4OTDiUR874AwKA7Gz3XXeDgWtuEBGR4vFrMfIOam24qzXu8jDhKL/ugAD39EgYDFU3AxYREVEl8Ksx8g4ibl+qW6cOLhw9Cp1OhwB/f1y5dg2NW7fGYQAdARwGsH/bNjS6+WZRh0y/cQN/Jyeje3Q0agYHlwhHQEBAQMWVMOFQHrnjtyUdubmu7+/vz54NIiJSDX5ikXcQkXBotVqEhYba/x8QEIBLJ07Y/2/KykKdWrUqvg2qSEFBAQwGAwIDAkTvUwpnqVKeqojfljRInb1Kp+NtVEREpDr81CLv4ELDXavVol7duvb/+3tiph9vTDgA6y1Dau29qcqEyc8PhQbDf2tolPV8sE17q/ZkjoiIqiUmHOQd1NxwFwTva0iq+XyqOnbbOh1Go3W2teIzrmm11mRDzX9PIiKq9phwkPqpOdkAmHAojSdj1+msP0RERF6Es1SR+nlDwuFt1JpwqDVuIiIiBWMPB5XLZDJBK2IKz+zsbEn12hbFEwQBJpOpcvWbzda1DUrIKbbuhhhSy+cWzTKUm5sr6hzKrN9sdjoIWOrfVGp5Wzw5OTnwETEIWVL9eXlAfr56roFNid4FVV8DF8q7so9bX8tOSH1OEBGR8jDhICQlJSEpKclhW0FBgYeicYHaewjUHr8zar0tiCt2ExERuR0TDkJCQgISEhIctp07dw4REREwGo0wGo2i6xJbVlN064pGo6l8/fn55TZwpdQvpXxhYSEAwGAwVO4c9HrrDESVjEdqeVv8/v7+7r/GFovDNVH8NbAp41qo8hpUoryUfdz6WnYiIyNDdJ1ERKRM/DqP1E/tPQRqj98ZtfYUqDVuIiIiBeOnK6mf2hvsao+/LGpsvKv1VjAiIiIFU2GLgMjLMOFQBo2Gs1QRERHJgGM4SP1kbLCbzWZs3roVF1NTERYaitiYGOj4Lbg4aks41BYvERGRSjDhICrD8pUr8fyECTh14YJ9W6PwcMyZPRv3DhrkwchUQm2JGRMOIiIiWfATlsiJ5StX4v7hw9HmwgUkA8gEkAygzcWLuH/4cCxfudLDEaqA2hrwaouXiIhIJfgJS1SC2WzG8xMmYIAgYAWA7gCMRb9XCAIGAHjhpZdgNpvdc0CO4VAGtcVLRESkEryliqiEzVu34tSFC/gWpTNyLYCJgoCY8+eRvH07OMS4AlqtdU0ONVDbLWBEREQqwa/0iEq4mJoKAGhdxuO27amXLlVJPKqmlka8VssZqoiIiGTChIOohLDQUADA/jIet20PrV+/SuJRNR+VdKKqJTEiIiJSISYcRCXExsSgUXg4Zmg0KHkzkAXATI0GjRs0QHTXrp4IT13U0pBXS2JERESkQkw4iErQ6XSYM3s2fgEwWKNxmKVqsEaDXwC8M2uW+9bj8OZbeTQadSQdaoiRiIhIpZhwEDlx76BBWLZ4MfaFhSEGQBCAGAD7w8OxbPFirsMhhdIb8zqddyd9REREHsb7CEj9ZGos3jtoEO7u358rjVeW0m9X4vUkIiKSlcJbAkSepdPp0Ds2Vt6DePu367YGvVLXG1F6QkRERKRyvKWK1E/tDXa1xy+GUnsR1DLGhIiISMWYcJD6qb3Brvb4xVBqo16pcREREXkRJhykfmpvsKs9fjGUetsSEw4iIiLZMeEg9VN7g13t8Yuh1Vp/lMbX19MREBEReT0FtgCIJFJ7g13t8YultMa9j0/1+dsTERF5EBMOUj+1NxrVHr9YSks4lBYPERGRl2LCQeqn9ga72uMXS6NRzlgOJcVCRETk5ZhwkPqpvcGu9vil8PPzdARW7N0gIiKqMkw4yDuotdGu1rhdpdMpY/C4UhIfIiKiakABn/xUHrPZjNdeew2NGzeGv78/IiMjMX36dAjFVm0WBAGTJ09GWFgY/P39ERcXh2PHjnkwag9QQiPWFdUt4QA839jnYHEiIqIqpdJWWvUxa9YsfPzxx/joo49w6NAhzJo1C7Nnz8aHH35oLzN79mx88MEHmD9/PrZt24bAwEDEx8cjNzfXg5FXMbU2INUad2X4+nr2vD2d8BAREVUzHDWpcFu3bsXdd9+N/v37AwAaNWqEb7/9Ftu3bwdg7d2YO3cuJk2ahLvvvhsA8NVXXyEkJAQrVqzAQw895LHYq5RaezjUGndl+fkBeXlVf1wfHy72R0REVMWYcChcTEwMPv30Uxw9ehTNmjXD3r178ddff+Hdd98FAJw8eRKpqamIi4uz7xMcHIxu3bohOTnZacKRl5eHvAoae5mZmQCA9PR0WCyWCuPMzs4GABQWFoo6L1udFosF6enpla+/oAAo0aOTk5MjKSap5U0mk8Nvl+rX60vFbSP1byq1vO0a2367u/4K98nKAko8t+S+BtkaDVBQINvfVGnXwO3XzAm3v5ZLEPu3ISIi5WLCoXAvv/wyMjIy0KJFC+h0OpjNZrz55pt4+OGHAQCpqakAgJCQEIf9QkJC7I+VNHPmTEydOlXU8bds2YKAgIBKnIFzBQUF9t9//vln5Ss0m+GTn1/5elywd+9el/ct9PPz+DfuKSkpHjmuxmyGzk3XTMw1sPj4wKLQ2ak8dQ3cwe2v5RJsCQoREakXEw6F++GHH/DNN99gyZIlaNWqFfbs2YNx48YhPDwcI0eOdKnOiRMnIjExsdwy58+fR1RUFHr06IGgoKAK67Q1CsQmJ75FDT9fX1/06tWr8vULAlDiW27bt+X+/v6iYpJa3mQyYe/evWjXrh2MRqNr9RuNZY5nkPo3lVo+MzMTKSkp6NixI2rUqOH2+kXtk5MDFPumW7ZroNEAgYHILqpfrr+p0q6BLNesBLe/lktIS0sTVY6IiJSLCYfCvfjii3j55Zdx9epVDBw4EKmpqahbty6mTJmCkSNHIjQ0FID1QzksLMy+3/nz55GZmYmwsDBcu3YNDRs2xNy5c9GvXz/o9Xro9fpyj5uRkQEAqFmzpqiEw6doETUxDW8A0BaNXdBqtahZs6Z76vf1tSYeLsYktbyN0Wh07Rw0GmvC4aZ4XI2/Ro0a7rsGUvcJCrLeWlWJY9jKl3sOBgPg6wufosaxXH9TpV0DWa5ZCbK8lothDwcRkfpV0xGr6pGdnY2UlBQkJiZiypQpSElJQXh4OM6ePYtLly6hcePGCA0Nxfr16+37XLlyBcnJydDpdFi2bBmOHDmCBQsWoEGDBh48kyqgtsHA1XXAeHFarXUci5x8fLjQHxERkQexh0PhBg4ciMWLF+POO+9Enz59sHv3bpw8eRIBAQFYuHAhXn75ZYwbNw5vvPEGmjZtisaNG2P48OHQ6XTYtGmT/TaNRo0aefZEqoLaGvBqi1cufn6A2exwa5XbaLXW3g0iIiLyGCYcCjdnzhx89dVX2L59O1q2bInw8HA88cQTOHXqFJKTkwEAEyZMQFZWFh5//HGkp6fDaDSib9++eP755/Hzzz+jXr16GDZsGF566SXoyugFKDlzlSpnhlFbA15t8crJYADkuHXGYKiea50QEREpCFs8CpeXlwdBEPDzzz8jJycHx48fxxtvvIHw8HD7LFQajQbTpk1DamoqcnNzUbduXfz2228wm81YvXo1XnvtNcyZMwdvvPFGmceZOXMmgoOD7T9RUVFVdYruo7YGvNpuAZOTRgP4+7s3OfD359+YiIhIAVTWQiMxLBYL6tevj08//RSdOnXCgw8+iFdffRXz588vc5+JEyfixo0b9p+DBw9WYcRuoraEQ23xyk2rBQIC3JN0+Ptbx24QERGRx/ETWeHq1q0LnU5XamrItLQ0+wxVJYWFhcHX19fh9qmWLVsiNTUV+fn58PPzK7VPyZmrbLNUqYpGY220ilio0OM0Gt7q44xWCwQGWnsmzGbp+9t6StizQUREpBj8ilXh/Pz80KlTJ4dZqCwWC9avX4/o6Gin+/To0QP//vuvwwrhR48eRVhYmNNkw6uopaHJb9/LptFYezqkPld9fP5LVoiIiEgxmHCoQGJiIhYsWIAvv/wShw4dwlNPPYWsrCyMHj0aADBixAhMnDjRXv6pp57CtWvX8Nxzz+Ho0aNYtWoVZsyYgYSEBE+dQtVRS0OejeKK6fXWBMLXt9zeIIvtVix3jwEhIiIit1BJ66x6e/DBB3H58mVMnjwZqampaN++PdasWYOQkBAAwJkzZ+yLbwFAREQE1q5di/Hjx6Nt27Zo0KABnnvuObz00kueOoWqo5aGvFoSI08rPq2t2ex4m5VWCxQWwqLXq+e6ExERVUNs9ajE2LFjMXbsWKePbdy4sdS26Oho/P333zJHpUAajev3/1cVrZbfxLtCpyudWPDvSEREpHi8pYq8j9K/7WbvBhEREVUjTDjI+yi9Qa/0hIiIiIjIjZhwkPdReoNe6fERERERuRETDvJOSu3l0Ok47oCIiIiqFSYc5J2U2oug1ESIiIiISCZMOMg7KbVhr9S4iIiIiGTChIO8k1arvF4Ordb6Q0RERFSNsPVD3svX19MROFJaPERERERVgAkHeS+l3b7EhIOIiIiqISYc5L00GuU08n18ODsVERERVUtMOMi7KSXh8PPzdAREREREHsGEg7ybTuf5weNKHMBOREREVEWYcJD383QvB3s3iIiIqBpjwkHez9fXc+MnlDSOhIiIiMgDmHBQ9aDXV6/jEhERESkEEw6qHnx9q34chU7H3g0iIiKq9hS2UAEpjclkglbE6tjZ2dmS6hUEwf7bZDK5vX6n+1gsQFZWmeVzcnIk1Z+bm2v/7fQcAgKAYtulnoPc5W3nm5OTAx8Ra5a45Rq4ubzc51Dd4ndlH7lfy1Jfl0REpDxMOAhJSUlISkpy2FZQUOChaGSk1Vp7HKri3DzRo0JERESkQEw4CAkJCUhISHDYdu7cOURERMBoNMJoNIquS2xZTdEgbo1GI0v9Ze4TGGjt5Sj6VrYyxygsLAQAGAwGx300GutxyhioLvUc5Cpvi9/f379qr4Eby1fVOVS3+KXsI/drOSMjQ3SdRESkTBzDQdWLRgMYDPIew2DgquJERERERZhwUPXj4yPf2hh+ftb6iYiIiAgAEw6qrvR69ycGOh2nwSUiIiIqgQkHVV/+/u4b2K3TWesjIiIiIgdMOKh68/evfE+HLdnguA0iIiKiUphwUPWm0ViTBRfHdFh8fKzrbTDZICIiInKKCQcRYB17ERAg/harol4NC1cSJyIiIioXp9MhstHprEmHwQAUFgJms3V1chuNxnr7le0nP99zsRIRERGpBBMOopJsCQVgXSBQEKzJBm+bIiIiIpKMCQdReZhoEBEREVUKx3AQEREREZFsmHAQEREREZFsmHAQEREREZFsmHAQEREREZFsmHAQEREREZFsmHAQEREREZFsmHAQEREREZFsmHAQEREREZFsmHAQEREREZFsmHAQEREREZFsmHAQEREREZFsmHAQEREREZFsmHAQEREREZFsmHAQEREREZFsmHAQEREREZFsmHAQEREREZFsmHAQEREREZFsfDwdACmbyWSCVltxXpqdnS2pXkEQ7L9NJpPb63dlH6nlc3Jy7L99fCp+Kckdj9Lid2UfpZ1DdYvflX3kfi3b/kZERKReTDgISUlJSEpKcthWUFDgoWiIiIiIyJsw4SAkJCQgISHBYdu5c+cQEREBo9EIo9Eoui6xZTUajf23HPVXZh+x5QsLCwEA/v7+sp6D2uN3ZR+lnUN1i1/KPnK/ljMyMkTXSUREysQxHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsmHEREREREJBsfTwdAymYymaDVVpyXZmdnS6pXEAT7b5PJ5Pb6XdlHavmcnBz7bx+fil9KcsejtPhd2Udp51Dd4ndlH7lfy7a/ERERqRcTDkJSUhKSkpIcthUUFHgoGiIiIiLyJkw4CAkJCUhISHDYdu7cOURERMBoNMJoNIquS2xZjUZj/y1H/ZXZR2z5wsJCAIC/v7+s56D2+F3ZR2nnUN3il7KP3K/ljIwM0XUSEZEycQwHERERERHJhgkHERERERHJhgkHERERERHJhgmHSiQlJaFRo0YwGAzo1q0btm/fLmq/7777DhqNBoMHD5Y3QCIiIiIiJ5hwqMD333+PxMRETJkyBSkpKWjXrh3i4+Nx6dKlcvc7deoUXnjhBcTGxlZRpEREREREjphwqMC7776LMWPGYPTo0YiKisL8+fMREBCAhQsXlrmP2WzGww8/jKlTp+KWW26pwmiJiIiIiP7DhEPh8vPzsWvXLsTFxdm3abVaxMXFITk5ucz9pk2bhvr16+PRRx8VdZy8vDxkZGTYfzIzMysdOxEREREREw6Fu3LlCsxmM0JCQhy2h4SEIDU11ek+f/31Fz7//HMsWLBA9HFmzpyJ4OBg+09UVFSl4iYiIiIiAphweJ3MzEwMHz4cCxYsQN26dUXvN3HiRNy4ccP+c/DgQRmjJCIiIqLqgiuNK1zdunWh0+mQlpbmsD0tLQ2hoaGlyh8/fhynTp3CwIED7dssFgsAwMfHB0eOHEFkZGSp/fR6PfR6vf3/XN2XyiQI1h+NxvpDREREVA72cCicn58fOnXqhPXr19u3WSwWrF+/HtHR0aXKt2jRAvv27cOePXvsP4MGDUKfPn2wZ88eREREVGX45C3MZiA3F8jKAkym/36bTEBODlBQ4OkIiYiISKHYw6ECiYmJGDlyJDp37oyuXbti7ty5yMrKwujRowEAI0aMQIMGDTBz5kwYDAa0bt3aYf+aNWsCQKntRBUym4G8POtvZ0mFIACFhdaf/HzAzw/w9a36OImIiEixmHCowIMPPojLly9j8uTJSE1NRfv27bFmzRr7QPIzZ85Aq2VnFblZfr412RDLYrH2ghQWAgYDb7ciIiIiAEw4VGPs2LEYO3as08c2btxY7r6LFi1yf0Dk3XJzXb9NqrAQyM4GAgLcGxMRERGpEr8WJyJHlUk2bCwWa9IhCO6JiYiIiFSLCQcR/ccdyYYNkw4iIiICEw4isikocP9sU7ZxHURERFRtMeEgov+mvZWDbQYrIiIiqpaYcBCR/L0Qtql1iYiIqNrhLFVE1V1+PqDTOX3IYrHg6rVrAABTVhYAwN/fH3v37UPtWrXQqGFD8cfJy+PMVURERNUQEw6i6kwQIOTlIctiAQAEBARAU2z9jKvXrqH+Lbc47DJq2DAsWrIEPj4+WPrVVxg8YIC4Y9kWD+TCgERERNUKb6kiqsa0+fnIzs6GMSwMxrAwZGdnOy3nD6BD0e9FS5YAAAoLC/HU+PHIyckRf8C8PM5aRUREVM0w4SCqrsxmaIt6NirSAkBK0W8AGDF0KG6OiEBqWhqW/vST+GMKAgeQExERVTNMOIiqq7w8l3cd9/TTeGzECADAdz/+KG3n/Hz2chAREVUjTDhU4Pz58/jf//6HOnXqwN/fH23atMHOnTvtjwuCgMmTJyMsLAz+/v6Ii4vDsWPHPBgxKZ7Z7PKsUY1uvhnt27bFkHvuAQD8vnEjTCaTtErYy0FERFRtMOFQuOvXr6NHjx7w9fXFr7/+ioMHD2LOnDmoVauWvczs2bPxwQcfYP78+di2bRsCAwMRHx+PXC64RmWpRIM/rndvaDQaNG/aFI0aNkRBQQE2bt4srRJ3LzBIREREisWEQ+FmzZqFiIgIfPHFF+jatSsaN26MO++8E5GRkQCsvRtz587FpEmTcPfdd6Nt27b46quvcOHCBaxYscKzwZMyCYJ1MT4X9YqNBQBoNBrc0acPAGCD1IRDEJh0EBERVROcFlfhVq5cifj4eAwZMgR//vknGjRogKeffhpjxowBAJw8eRKpqamIi4uz7xMcHIxu3bohOTkZDz30UKk68/LykFfB/fuZmZkAgPT0dFhEDCy2zW5UKLIha6vTYrEgPT3d7fW7so/U8ra/ke23p+MRXT4vD8jPt98GlVXsdqj0GzdQUCwRuHHjRqndo5o1s1+zrp06YcGiRVi/caPT62ibwcppTFotEBjo2jkUUe01KKK0+F3ZR+7Xsti/DRERKRcTDoU7ceIEPv74YyQmJuKVV17Bjh078Oyzz8LPzw8jR45EamoqACAkJMRhv5CQEPtjJc2cORNTp04VdfwtW7YgQIbF2myN2oKCAvz5559ur78qpaSkeDoESXQ5OdAU+/++ffvs//47ORkGg8H+/xsZGaX2P3z4MC5euAAA8C1as2PfgQNY9/vvCPD3lxRLoV5vTTwqSW3XoCQ1xy/3a7msqZqJiEg9mHAonMViQefOnTFjxgwAQIcOHbB//37Mnz8fI0eOdKnOiRMnIjExsdwy58+fR1RUFHr06IGgoKAK67Q1CsQmJ75Fi7/5+vqiV69ebq/flX2kls/MzERKSgo6duyIGjVqeDweUeULC4GiXgeTyYS9e/eiTZs29oe7R0cjsNj+V65eLVVF165dUbdOHfv/I956C2fPnYOfwYAePXo4lLX1cPiXlYj4+QF6vbRzKEaV16AYpcXvyj5yv5bT0tJElSMiIuViwqFwYWFhiIqKctjWsmVL/Fg0FWloaCgA64dyWFiYvUxaWhrat2/vtE69Xg99sUaeMxlF32zXrFlTVMLh42N9KhmNxgrLAoC26FttrVaLmjVrur1+V/Zx5RgAUKNGDVnOQZbyOTkODXwACCxWvmZwMAKL3eZU4OS2l+DgYIfz7dGtG747dw77Dh7E3SVWHa8wJo0GKPZYtbgGTiglflf2kfu1zB4OIiL146BxhevRoweOHDnisO3o0aNo2LAhAKBx48YIDQ3F+vXr7Y9nZGRg27ZtiI6OrtJYSeEqOVi8LN27dAEA/L1jh/SdZYqJiIiIlIM9HAo3fvx4xMTEYMaMGXjggQewfft2fPrpp/j0008BWGcKGjduHN544w00bdoUjRs3xmuvvYbw8HAMHjzYs8GTsri47kZFbAnHtp07IQgCNBpNBXuUUFgI+PCtiIiIyFvxU17hunTpgp9++gkTJ07EtGnT0LhxY8ydOxcPP/ywvcyECROQlZWFxx9/HOnp6ejZsyfWrFnjMPiXSK6ehPZt20Kv1+PK1as4fuIEmhRN2ezpuIiIiEgZmHCowIABAzCgxL3xxWk0GkybNg3Tpk2rwqhIdWTq4dDr9ejYrh2St29H8vbt0hMOQQAsFrfMVkVERETKw094orLYGsIi1iFRPJnPo1LjOAD2chAREXkx9nAQFWc2A7m51t8lG+g6nXWsga+vdXYlNXGxd+NiGWu5lGRLOJJdTThk6n0hIiIiz2PCQQRYk4ucHOs37WXd2mM2W3/y863rR6iJiz0I23fuFFUupls3AMDeffuQmZkpak0JB0w4iIiIvBZvqSIqKACyssQ3ygUByMsDsrOt/1YDFxv023btElXupgYN0PDmm2GxWFyfHpdJBxERkVdiwkHVW26u9ccVZjN0ubnKH+NhNrucGCVv3y667K0xMQCATVu2uHQsjuMgIiLyTkw4qPrKzbX2blSCBrD2dCg56XCx5+Dq1as4cOiQ6PK9evYEAGzYvNml47GHg4iIyDsx4aDqKS+v0smGnSAo+/YqF3sOpCYOt916KwDrAoCZmZnSD1iJnhgiIiJSLiYcVP0UFFgHfruTIFgHnSuRi70vv65bJ6l840aN0OSWW1BYWIj1f/7p0jEV3VNERERELmHCQdWLxWLt3ZCD2ezeui0Wa3Jk+3HlliNBcKnXwGKxYNXatZL363vHHQCAX9askbxv0YFd24+IiIgUiwkHVS85OfLetpOfX7nBz4JgrSMry/pjG9Sem2u9bSsz03oOYpMPFxvwW7dtQ9qlSzAajZL2u7t/fwDAz6tWocCVW9aYcBAREXkdJhxUfRQUSGrQCoKA7Tt3YvPWrbh85Yr9x1JRHa72ctim583LKz/OwkJr8iEmeXJxIPY3P/wAAOhX1GMhVq+ePVGvbl1cuXoVv61fL/3ATDiIiIi8DhMOqh4EAUJuLrKyspCVlQVBRC/Hq9Omodttt+HWu+5C/Vtusf9cvXat/B0tFmljRAThv14MKb0vhYXWBKW8pMKFBnxGRoY94Xjwvvsk7evj44OHH3gAADB/4ULJx+ZMVURERN6HCQdVD3l5yM7KgjEsDMawMGRnZ5dbfP/Bg3jr3Xft//cH0KHotyj5+eKSB9tgc1dnzLLNkFXW/i4kHEkLFiAzMxMtmzdHj+7dJe//1KOPQqPR4Jc1a7DvwAFpO7s45oSIiIiUiwkHeT/b4GsJkhYsgCAI6HfnnQCAFgBSin6LYluNvKIyUsZjlCc31/nYEYkJx5mzZzFjzhwAwCvPPw+tVvpbRLOmTfHAPfcAAJ5/9VWYpZ4fb6siIiLyKkw4yPtJHFNRWFiIpT/9BAAY/b//uX7cisaMuCvZKKs+F3oLhowcCZPJhJhu3TCs6NYoV7w1dSoCAwORvH07nnvpJWkDyJlwEBEReRUmHOTdLBbJs0Zt3bYNV69dQ53atV26pchBWWM58vLkGa9QfCC5yPr/3r7d/u8Dhw6hfv36+OzDD5GTk4MsJ7eeZWVn28fClPVTr25dfPzeewCAL5csQcvOnfHm228jX8zYFiYcREREXsXH0wEQySo/H9DrJe2ytmh2pfjbb4ePTyVfIgUF1uNrNP9tKyx0/8KDNrbbtAICRDfcV65e7fD/S5cuIapr1zLLN27dWnJYx0+exKTp0/HcU0/Bz8+v/MIcOE5ERORV2MNB3ksQXFoTY8OmTQCAuD593BNH8duJbDNSyclstiY0IhOO+NtvlzeeYjTFE6+ycNA4ERGRV2EPB3mvggLJjdes7GzsSEkBAPTu2dNpGTFT6jrIzwds3+oX/7cTFovFPu2uqej2pFVr1yIwIACPjx5dce9A8WOKHPB91x13IO34cQBAgL+/Q1Jw+epVNG7dGocBdARwGMDJ/ftRr04dUXWbsrIAAMbAQGv9AQEV78RbqoiIiLwKEw7yXi5MNbtr924UFhbipgYN0KhhQ1y5erVUmXPnz6N+vXriK7X1tFgsEPLykFUUV0BAQKlv/K9eu4b6t9zitJo9+/bhs48+En/M3FzA17fColqtFvVr13b6mL+/Py6dOAHgv+Th5ptuEj17lS05CyxKOEQTBMfb0IiIiEi1eEsVeSez2aVvyrft3AkA6NGtW5m3/9h6QCQpKLCuBZKTI2otEGfrfnz+1Vc4cOiQ+GPm51d6PIRWq0W9unWtP3XqoF6dOi5NlSsZb6siIiLyGkw4yDu5uJDe30UJR0y3bmWW2elKwpGbKymm4ut+PDZiBO4ZOBAA8OkXX4g/ppi1QJSKCQcREZHX4C1VVC6TySTqG+2KVu4uyXarjSAIMJlMbq8fJhNycnIcNhWf4tWUlVVqLEZ2djZ27NoFAGjbujVMJpP9NqLitqekwGQyIbdo8Hdubm7F55CdjZzCQuQUSzqcxeDseIMHDEB2Tg5++r//w9IVKzB90iSn16Tk+cIWU34+4GS2LUnxO6tfBKn72MtnZjqNuazyOTk5omYUk/o8kru80uJ3ZR+5X8uuPO+IiEhZmHAQkpKSkJSU5LBN0kJtSmM2u/QN+clTp3A9PR16vR5tW7Uqs9yBQ4ekNYJst3e5MGMWAES1aIGgGjUQGBCAi6mp2HfgANq1aVP+TsXPv4yEQ9HYw0FEROQ1VNYKITkkJCQgISHBYdu5c+cQEREBo9EIo9Eoui6xZW3jIzQajfvrz8tzmKHJtk/xMRnGwMBSA5n/OXgQANCxXTvULhpEneNkCluz2Yyjx4+jVYsWAACDwVB+XDk51ng0GlgqiMHZ8YxGI+rUqYPbevXC//36K/76+2/0iI4u83BGo7H0+JWAgFKzVhUWJUAVxu+sfomk7mMMDCx3Ni8b2zn4+/vLeg5ylVdq/FL2kfW1DCAjI0N0nUREpEwcw0Hex8WB0rbbqbp36VJh2b937BBXafGeDUGo1CDu23v1AvDfOiHlKtlDINdCg3JhDwcREZHXYMJB3qUSjfrtciQcJW9Nc/G2KgDoVbQuyJZt2+zfjJepZIO9sFBdjXg1xUpERETlYsJB3sXFZCMrKwv7i6acje7atcLyydu3i1sAsGTCUYlF7dq0aoWgoCCYTCb8s39/+YVLxubiqusew4SDiIjIazDhIO/iYqN6286dMJvNuCk8HBE33VRuWZ1Oh/MXLuDs+fMVx1Ky4VyJW6p0Oh2ii3pfkrdvL7+wswa7miYCYMJBRETkNZhwkHdxsUH/V3IyAHG3U9lmsKqw0S9DA9/W++JSwuHiYogewYSDiIjIazDhIO9hsbjcoN60dSsAoEf37hWW7W5r9G/bVnYhF29hqmhshi3hED2GpCQ3JEFmsxkbN2/Gt0uXYuPmzTBXcjVzIiIi8m5MOMh7uNjwzcvLw5a//wYgLuGwrUL+V9E+7ozl0JEj5T7etVMnAMDxkydx+coV6QeoZHLw8+rVaNKqFfr0749hjz6KPv37o0mrVli+cmWl6iUiIiLvxYSDvIeL4ze2btuG3Nxc1K9XDy2aNauwfEzXrtDpdDh+8iQuXb7s1lh27dlT7uM1a9ZEy+bNAQDbXOnlcHFRRMCabAwfMwZtLlxAMoBMAMkA2ly8iPuHD2fSQURERE4x4SDv4eLtVGvXrwdgXeei+OKAZalRo4Z9rMeuvXudF3KxJ2HX7t0VlrEdO7m8hKO8pMKF2MxmM1597TUMEASsANAdgLHo9wpBwAAAL7z0Em+vIiIiolKYcJD3cDHh+L9ffwUA3NGnj+h9+t5xBwBgu7MEoRJjSXZISDhcHsfhQu/L5q1bcfriRbyC0m8aWgATBQEnz5/H5qKxMJXGQeNERERegwkHeQcXG/iHjhzBwcOH4evriztuu030fnf37w8A2P3PP7iRkeH4oIvf8l+5ehUnTp6ssJxt4Pi2nTsrXgDQGRfiu5iaCgBoXcbjrUuUqzQRPU1ERESkDkw4yDu42Mj/+vvvAQDxt9+OmsHBovdr1bIlWjRtisLCQqxcvdrxwUqMJRF77KCgIGRlZZW9AGB5DXYXemDCQkMBAGUtN7i/RDkiIiIiGyYcEh06dAhTpkzBbbfdhsjISISFhaFt27YYOXIklixZgry8PE+HWD250MORl5eHhYsXAwBGDB0qaV+NRoMH77sPAPDF1187rjruYvIj9nYkrVaLmKJeDpdvYZIYY2xMDBqGhWGGRoOSf2kLgJkaDRo3aIDYmBjX4iEiIiKvxYRDpJSUFMTFxaFDhw7466+/0K1bN4wbNw7Tp0/H//73PwiCgFdffRXh4eGYNWsWE4+q5kLCseibb5CaloYG4eH2W6Sk+N+DD8LP1xe7//kH6/74w7pRyixQJWL+86+/RB/71h49rPts2SJ6HwcSe2F0Oh3enD4dvwAYrNE4zFI1WKPBLwDemTULOp3OtXiIiIjIa/l4OgC1uO+++/Diiy9i2bJlqFmzZpnlkpOT8f7772POnDl45ZVXqi7A6k5iwpF2+TImTZ8OAJjw3HPw8/NDfn6+pDrq1qmDu26/HSvXrMG4l1/Grk2b4C+lwV0s5hs3blQ4JW5xvXv2BABs3LwZFosFWq3E7w5c6IW5u18/LF6wAJMmT0bMhQv27Y3Dw7Fs1izcO2iQ5DqJiIjI+zHhEOno0aPw9fWtsFx0dDSio6NR4IYVnUkCiQnHgyNH4srVq2jTqhWeeuwxlw879L778PeuXTh05AgeGj0aSz76CIF6vbidi8W8oShxaNyoEU6eOlXhrl06dUKNGjVwPT0dKXv2oHPHjtICFwSXeoXu7tcPD913HzZv3YqLqakICw1FbEwMezaIiIioTLylSiQxyUZlylMliPy2vvhCeTt370bN4GAsmj8f+fn5yMrKQlZ2tvUnK+u/n+zsUvXYy2Rnw9fHB0nvvANfX1+sXL0abbt3x2tvvolMk6nigIo1+H9dtw4AcFtsrKhz8fHxwe29ejns60DMLE8ujjXR6XToHRuLoUOGoHdsrDzJBmepIiIi8hrs4XDRjh07sGHDBly6dAmWEt8Uv/vuux6KqpoS+U39+j//dPh/+o0b6CSygV9c49ZlTQ4LnLh0CW/Mm4cXnnmm4oqK4rZYLPa1QOL69MHnRQPZK9I/Ph4rfvkFP69ejddeesnxQTENdhenEq4STDiIiIi8BhMOF8yYMQOTJk1C8+bNERIS4rA6tZiVqsnNRDac4+PiMGPOHJmDsfLxEfHSKop7y99/42JqKoKCgtAzOlr0MQb27QutVotdu3fjxMmTqF+v3n8PinkeKnlxPb6OiIiIvAYTDhe8//77WLhwIUaNGuXpUAgQ3XDuGR2NtOPHAQAB/v6lkkNTVhYAwBgYaN92+epVNG7dGocBdARwGMDJ/ftRr04dpN+4gb+Tk9E9Otq6hkdBAVA0O1mAv7/o2L/69lsAwD0DBkAvdvwHgJD69XFbr174fcMGfPXtt3jh2Wf/e5A9HERERKQQTDhcoNVq0aNoWlJSAJEJh1ardewFKFWNtZ7AYgmHv78/Lp044VCuTu3a0Gq1KCgogMFgQGBAgHWfvDxA4niGq1evYsnSpQCAR4YPl7QvAIx++GH8vmEDPl20CM888cR/CYvYhEOpvRxMOIiIiLwGB427YPz48UhKSvJ0GGQjY6NZq9WiXt26Dj9lTkHrQo/B2++/j+zsbHRo186lRfPuHzwY4WFhuJiaii++/vq/B8Q22JlwEBERkczYw+GCF154Af3790dkZCSioqJKzUi1fPlyD0VWTSml0Sxx1qfd+/fj3XnzAABTX3nFpfE/fn5+mPTii3g6MRFvvP027rrjDrSOihLdYNco9bYqqeuKEBERkWLxU90Fzz77LDZs2IBmzZqhTp06CA4OdvihKqaEhEMQJMdxz6OPoqCgAPcMHIgBd93l8qEfHz0a3bt0wY2MDAx68EHsTEmxPiAi6dAo4W/nDHs4iIiIvAZ7OFzw5Zdf4scff0T//v09HQoBykg4RPYULP3lF/u/r1y7hnatW+Ojd95BdtF6H2Wt+xFQNKC9LF998gnuuPtunDx9Gl1690a3zp3x6+LFqBUUVH5ASvjbOcOEg4iIyGsw4XBB7dq1ERkZ6ekwCFBOg1lkwpFZInHYu38/GjRvXu4+5a37UZZtO3fC18+vwnKKvaWKCQcREZHX4C1VLnj99dcxZcoU+7fS5EEqSzjievaUOZD/iBoTopS/X3FMNoiIiLwKezhc8MEHH+D48eMICQlBo0aNSg0aT7HdQ+8FTCZT2bMyFSM1+bJNQSsIAkwmk+v1m81AGY/l5ORIiklq+dzcXPtvU2EhUFhYfvmcHDQMD8fx5GQAQIDBAI2vL2Aw2MtcuXoVrbt1c1j3Y/+2bahbp47o+P2L1gCx4L+1RZzJs8WfkwOTiIUKpf59XNknJycHgkYDiHhOFK8/JydH1GKLUp+ncpdXWvyu7OO213IZXHneERGRsjDhcMHgwYM9HYJbJSUllZrmt6CgwEPRSKSUb+ilrAVSPHnQ6YBiiwT6Gww4sW8fAGsiAAANwsNFJX3aop4BW8KBgoIKkyApsVcZzlBFRETkVZhwuGDKlCmeDsGtEhISkJCQ4LDt3LlziIiIgNFohNFoFF2X2LK22300Gk3l6i8oqHCxPSn1SylfWNSYNxgMMPr5iZ4Wt/hK5tBqgeL/BxBUNNDb9m2xy/GX0/sDFItfr5flGru8j5+f5Gvg7+8v6znI/RxSWvxS9nHba7kMGRkZouskIiJl4leJIglK+xaYrJRyXVyNQ874xfYUKOVvaCNxtXYiIiJSNiYcIrVq1Qrfffcd8vPzyy137NgxPPXUU3jrrbeqKDKSi9lsxsbNm/Ht0qXYuHkzzOX1YCit0Q5YB1+rceA4b6kiIiLyKrylSqQPP/wQL730Ep5++mnccccd6Ny5M8LDw2EwGHD9+nUcPHgQf/31Fw4cOICxY8fiqaee8nTIVAnLV67E8xMm4NSFC/ZtjcLDMWf2bNw7aFDpHZTWaLfR6Soex6G02JlwEBEReRUmHCLdfvvt2LlzJ/766y98//33+Oabb3D69Gnk5OSgbt266NChA0aMGIGHH34YtWrV8nS4VAnLV67E/cOHY4Ag4FsArQHsBzDj4kXcP3w4li1e7Jh0KK3BXpyYxruS4meyQURE5HWYcEjUs2dP9KzCtRSoapnNZjw/YQIGCAJW4L97DrsDWCEIGKzR4IWXXsLdxVeZV1KDvSQmHERERORh/HQndXNzY3nrtm04deECXkHpF4cWwERBwMnz57F561b3xCB3Y19tDXi1xUtEREQV4qc7UTGpaWkArLdROWPbfjE19b+NSuohKElMA17kKulVgjNUEREReR0mHCqRlJSERo0awWAwoFu3bti+fXuZZRcsWIDY2FjUqlULtWrVQlxcXLnlVU3MLEwShIaEALCO2XDGtj0sNPS/jZVJONwcv9P6KzqGkhIm9nAQERF5HX66S3Ch2IxFVen7779HYmIipkyZgpSUFLRr1w7x8fG4dOmS0/IbN27E0KFDsWHDBiQnJyMiIgJ33nknzp8/X8WRq09Mt25oFB6OGRoNSn7vbwEwU6NB4wYNEBsT898DSmqwO6OWXgONhgkHERGRF+KnuwStWrXCkiVLqvy47777LsaMGYPRo0cjKioK8+fPR0BAABYuXOi0/DfffIOnn34a7du3R4sWLfDZZ5/BYrFg/fr1VRy5+uh0OsyZPRu/ABis0SAZQCaA5KL//wLgnVmzoKuiRryktUDK4qOSuSHUkhgRERGRJEw4JHjzzTfxxBNPYMiQIbh27VqVHDM/Px+7du1CXFycfZtWq0VcXBySk5NF1ZGdnY2CggLUrl27zDJ5eXnIyMiw/2RmZlY6drW6d9AgLFu8GPvCwhADIAhADID94eGlp8SV0c+rV6Nd167o078/hj36KPr0748mrVph+cqV0ipSS0NeLXESERGRJEw4JHj66afxzz//4OrVq4iKisL//d//yX7MK1euwGw2I6RobIFNSEgIUosPXC7HSy+9hPDwcIekpaSZM2ciODjY/hMVFVWpuNXu3kGD8O+BA9iwahWWfP45NqxahWP791dZsrF85UoMHzMGbS9edOhlaVO0FoikpEOrVcetSmrpiSEiIiJJ+AkvUePGjfHHH3/go48+wr333ouWLVvCp0RDKSUlxUPRlfbWW2/hu+++w8aNG2EwGMosN3HiRCQmJtr/f/78eXUkHTIOutbpdOgdGytb/QCcxi9lLRDRt3bpdGXPRqWEMSgcv0FEROS1mHC44PTp01i+fDlq1aqFu+++u1TC4U5169aFTqdDWtF0rTZpaWkILT5TkhPvvPMO3nrrLfz+++9o27ZtuWX1ej30er39/xkZGa4HXZXknuVJbk7i37x1K05duIBvUfZaIDFFa4GIToh8fICCgspGKx/2bhAREXktfspLtGDBAjz//POIi4vDgQMHUK9ePVmP5+fnh06dOmH9+vUYPHgwANgHgI8dO7bM/WbPno0333wTa9euRefOnWWN0aO8MOGwrfEhaS2Qiih9fITS4yMiIiKXMeGQ4K677sL27dvx0UcfYcSIEVV23MTERIwcORKdO3dG165dMXfuXGRlZWH06NEAgBEjRqBBgwaYOXMmAGDWrFmYPHkylixZgkaNGtnHehiNRhiNxiqLu0p4YcJhW+NjP6y3UZXkdC0QMcfR6QBns1wp4W/IHg4iIiKvxU95CcxmM/755x/cdNNNVXrcBx98EJcvX8bkyZORmpqK9u3bY82aNfaB5GfOnIG22P3vH3/8MfLz83H//fc71DNlyhS8/vrrVRm6/JTQWK4MJ/HHxsRY1wK5eBErBMHhtir7WiDh4Y5rgYhRVsLhaVqt+q8jERERlYkJhwTr1q3z2LHHjh1b5i1UGzdudPj/qVOn5A9IKdTeUHUSv20tkPuHD8dgjQYTBQGtYe3ZmFm0FsgyV9YC8fEB8vPdEbV7sXeDiIjIq3FaGFI/NScdZczMdO+gQVi8YAH+CQ1131ogOp0y/1ZMOIiIiLwaP+lJ/TQaz07tWplGfDn73t2vHwbEx2P3P//gYmoqwkJDERsTU7lVzn19ldXLodVywDgREZGXY8JB6ufpb+1lSjgAGdYCcZZwePLv5+vruWMTERFRleAtVaR+nl4wTsaEw+2c9Sgw4SAiIiIZMeEg9fPiHg5ZlGzke+rv5+Pj+WtHREREsmPCQern6Uarq8ev7g19Pz9PR0BERERVgAkHqZ9ab6nyVNwajWMvhyfi4GBxIiKiaoMJB6mfpxMOwLWkw5NxF084PNHbwd4NIiKiakMBLTWiSmLC4dqxbUlHVSccJXtYiIiIyKtxWlwql8lkglZEwzg7O1tSvULRuhmCIMBkMlW+/pwcwGwusSlHUkxSy+fm5tp/m/LzAYul/PIl67dYyl0TQ/b4LRYIAHLz8kRdA6n1l7mPwQCUcTypzyNb/Tk5OfARsYCg1PrlLq+0+F3Zx+2v5RJced4REZGyMOEgJCUlISkpyWFbQUGBh6JxkVZbKuGoUmrr4QAAjQaWqh5ArtOxd4OIiKiaYcJBSEhIQEJCgsO2c+fOISIiAkajEUajUXRdYstqihq5Go3GPfX7+QF5eZWKSWr5wsJCAIDBYIDRYABEJmnGwEBrI1/kcWSN39cXhsBAWa6x030CA0UlWlKvgb+/v6znIPdzSGnxS9nH7a/lEjIyMkTXSUREyqSAm9+J3EABvQWSeDre4gyGqjmOr6+yzpuIiIiqBD/9yTt4uiEr9fiejreIAFjX5RAxfqBSNBpAr5f3GERERKRIymj1EFWWVuvZxeykJhAKWYNCsMVtMMibBPn7K2OxQSIiIqpyTDjIe3h6mlk5y8vFFodGI9+tVXq9YhIsIiIiqnoKafUQuYEnG/EajbRv8BWScAjFY9bp3J90+PpykT8iIqJqThmtHiJ3kHscQkXEJhE6nWJuLxJKxuHr676kw511ERERkWox4SDv4enbdqQkHErhLGZ3JApMNoiIiKgIEw7yHhqNZxvzYo/t6Z4Ym/JuA/P1Fb1mRqk6DQYmG0RERGSnkJYPkZvodJ5bcVxM49zTSVFxFcWh1VqTjvx8609FfH2tA8QVcrsYERERKQMTDvIuPj7iGsdyEJNwKCXZAMT3Xvj5WX8KC63JnE4HCIL1MZ3O+uPjw0SDiIiInGLCQd7FNiDb1iCuSrZblMo7tkJmpwIgPRbbAoEBAfLEQ0RERF5JQa0fIjfxZC9CRY14pYzfAJSV/BAREZHXYouDvI9SB45rNMpp5CtpLAkRERF5NYW0fojcyJO9COUdW0m9G0w2iIiIqIow4SDvo9V6riehvEX9lNTIV1IsRERE5NWYcJB38vX13LGdNeY1GmX1cCgpFiIiIvJqTDjIOykt4fBkPCUpaSwJEREReT22Osg7ebJHwdlxlZRwsHeDiIiIqhATDvJenmrka7WO4zh0OmX1KHD8BhEREVUhBbWCiNzMx8dzDf3ivQh+fp6JoSzs4SAiIqIqxISDvJunejlsvQhKGyxesveFiIiISGZMOMi7+fp6poFtSzLYu0FERETVHBMO8m4ajWd6OWzHVdJgcYAJBxEREVU5Jhzk/fz8PNPLERCgrNuXtFoOGCciIqIqx4RDZd566y1oNBqMGzfOvi03NxcJCQmoU6cOjEYj7rvvPqSlpXkuSKXRaAC9vmqPqdUCgYHKSjiU1ttCRERE1QITDhXZsWMHPvnkE7Rt29Zh+/jx4/F///d/WLp0Kf78809cuHAB9957r4eiVChf36r9dt9gsP5W0i1MTDiIiIjIA5hwqITJZMLDDz+MBQsWoFatWvbtN27cwOeff453330Xt912Gzp16oQvvvgCW7duxd9//+3BiBWoqno5fHz+S26UMmjcx0dZvS1ERERUbSjo61cqT0JCAvr374+4uDi88cYb9u27du1CQUEB4uLi7NtatGiBm2++GcnJyejevXupuvLy8pCXl1fu8TIzMwEA6enpsFgsFcaXnZ0NACgsLBR1PrY6LRYL0tPT3V5/mfvk5gIFBU7L5+TkSDqGyWRy+A3A2qgPDASK/32zswGzWXL9bi0fEOAYE/67xrbfFXHbNXBjebnPobrF78o+cr+Wxf5tiIhIuZhwqMB3332HlJQU7Nixo9Rjqamp8PPzQ82aNR22h4SEIDU11Wl9M2fOxNSpU0Ude8uWLQgICJAcc0UKihr9BQUF+PPPP91ef3l0ubnQCILb6tu7d6/932Y/Pwglbt3SmM3Q5ee77XhSCVotzOX07qSkpFRhNPJQ+zmoOX65X8u2BIWIiNSLCYfCnT17Fs899xzWrVsHg21cQCVNnDgRiYmJ5ZY5f/48oqKi0KNHDwQFBVVYp61RIDY58S0aT+Dr64tevXq5vf5y97FYrL0OJZIOWw+Bv7+/qPpNJhP27t2Ldu3awWg0Wm+fKqthn5WFnKwsSfVLjafM8gaD0/EbmZmZSElJQceOHVGjRo0K63frNXBTebnPobrF78o+cr+WOQEGEZH6MeFQuF27duHSpUvo2LGjfZvZbMamTZvw0UcfYe3atcjPz0d6erpDL0daWhpCQ0Od1qnX66GvYDxDRkYGAKBmzZqiEg6fosHRRqOxwrIAoNVq7b9L9s64o/4K9wkOtiYdlTyGrXzNunWB8hKDwED4XLkiqX6p8Tgtb5stqxw1atTwzDVw0zEA+c6husXvyj5yv5bZw0FEpH5MOBTu9ttvx759+xy2jR49Gi1atMBLL72EiIgI+Pr6Yv369bjvvvsAAEeOHMGZM2cQHR3tiZDVQaezJghFvQKV4uNTfrIB/DdLltlc+eNJUdXTARMRERGVwIRD4WrUqIHWrVs7bAsMDESdOnXs2x999FEkJiaidu3aCAoKwjPPPIPo6GinA8apGFuikJtb6vYqsSy2xEUMvb5Ur4qsfHyUNS0vERERVUtsjXiB9957D1qtFvfddx/y8vIQHx+PefPmeTosdfDxsc7glJsrfV+DARYp097qdFW7FgZ7N4iIiEgBmHCo0MaNGx3+bzAYkJSUhKSkJM8EpHZa7X/TxoqZTcrX1zpA3JXbo/R669S5bpwlyyk/P+t5EREREXkYEw4iG19f64/BYE0mzGZrYqDRWBvvOl3lF9DTaKxJhys9KmJptcpZcJCIiIiqPSYcRCXJPfbB19eazJSxAGGlaDTWMSVcVZyIiIgUgvdcEHmCXi/PLU9y1UtERETkIrZMiDxBo7GOG3FncqDXV+2gdCIiIiIRmHAQeYrt9id3JB16PcdtEBERkSIx4SDyJNsMWTqd63UYDEw2iIiISLGYcBB5mu32KqlJg1YLBAbyNioiIiJSNM5SRaQUtjEY+fnlz2Blm/aWiQYRERGpABMOIiXRaq23SOn1gMVinT7X1vNhWwuEs1ARERGRijDhIFIijcaaXOh01uSDiIiISKX4VSkREREREcmGCQcREREREcmGCQcREREREcmGCQcREREREcmGCQcREREREcmGCQcREREREcmGCQcREREREcmGCQcREREREcmGCQcREREREcmGCQcREREREcmGCQcREREREcmGCQcREREREcnGx9MBkLKZTCZotRXnpdnZ2ZLqFQTB/ttkMrm9flf2kVo+JyfH/tvHp+KXktzxKC1+V/ZR2jlUt/hd2Ufu17Ltb0REROrFhIOQlJSEpKQkh20FBQUeioaIiIiIvAkTDkJCQgISEhIctp07dw4REREwGo0wGo2i6xJbVqPR2H/LUX9l9hFbvrCwEADg7+8v6zmoPX5X9lHaOVS3+KXsI/drOSMjQ3SdRESkTBzDQUREREREsmHCQUREREREsmHCQUREREREsmHCQUREREREsmHCQUREREREsmHCQUREREREsmHCQUREREREsmHCQUREREREsmHCQUREREREsmHCQUREREREsmHCoRJJSUlo1KgRDAYDunXrhu3bt5dbfunSpWjRogUMBgPatGmD1atXV1GkROQSQQAKC4G8vP9+Cgut24mIiFSMCYcKfP/990hMTMSUKVOQkpKCdu3aIT4+HpcuXXJafuvWrRg6dCgeffRR7N69G4MHD8bgwYOxf//+Ko6ciCpkNgPZ2YDJBOTkAPn5//3k5Py33Wz2dKREREQuYcKhAu+++y7GjBmD0aNHIyoqCvPnz0dAQAAWLlzotPz777+Pu+66Cy+++CJatmyJ6dOno2PHjvjoo4+qOHIiKpMgANnZ8MnPrziZKCy0JiU5OezxICIi1fHxdABUvvz8fOzatQsTJ060b9NqtYiLi0NycrLTfZKTk5GYmOiwLT4+HitWrCjzOHl5ecjLy7P/PyMjo3KBE1HZzGbXei1siUdAAKDTyRObQvE9iYhIfrb3WsHNX24x4VC4K1euwGw2IyQkxGF7SEgIDh8+7HSf1NRUp+VTU1PLPM7MmTMxderUygdMROUrLLQmG66yWKz7+/tXq6QjIiLC0yEQEVUbV69eRXBwsNvqY8JBAICJEyc69IqcPXsWrVu39mBERF6ossmGjSBUu6Tj7NmzCAoK8nQY1UJGRgYiIiL4N6+meP2rtxs3buDmm29G7dq13VovEw6Fq1u3LnQ6HdLS0hy2p6WlITQ01Ok+oaGhksoDgF6vh16vt//fnVktEeG/26jcxZZ0BAQAWu8fjhcUFMTGTxXj37x64/Wv3rRu/lzx/k8plfPz80OnTp2wfv16+zaLxYL169cjOjra6T7R0dEO5QFg3bp1ZZYnIpkJApCbq556iYiI3Ig9HCqQmJiIkSNHonPnzujatSvmzp2LrKwsjB49GgAwYsQINGjQADNnzgQAPPfcc+jVqxfmzJmD/v3747vvvsPOnTvx6aefevI0iKqv3Fzr2As5mM3W+g0GeeonIiKqJCYcKvDggw/i8uXLmDx5MlJTU9G+fXusWbPGPjD8zJkzDl1fMTExWLJkCSZNmoRXXnkFTZs2xYoVKySNySh+exURVUJhofVHBIvFgqvXrjlsq1O7dsVd2wUFgK9vtRnPQfLS6/WYMmUKPweqKV7/6k2u68+EQyXGjh2LsWPHOn1s48aNpbYNGTIEQ4YMcfl4fKMhqjzBYkH21auAICAgIAAajabc8levXUP9W25x2HbpxAnUq1u34oPl5gKBgZUJlwiA9f3/9ddf93QY5CG8/tWbXNefYziIiGSSnZ4OY2gojGFhyM7OrrD8iZMnAQD+ADoU/RbNYrH2dBARESkMEw4iIjlYLEB+vqRd3v3oIwBACwApRb+LL8hZobw8rkRORESKw4SDiEgOUhIFWOe+/3n16lLb15aYca5cgiD5uERERHJjwkFE5G4Wi+iB4jarf/vNaW/GL2vWSDt2QQF7OYiISFGYcFRjSUlJaNSoEQwGA7p164bt27d7OiQi7yDxVirAmnA4s2HTJhRKTF5cOT5VL1Lf/9PT05GQkICwsDDo9Xo0a9YMq530yJE6SL3+c+fORfPmzeHv74+IiAiMHz8euVwDSJU2bdqEgQMHIjw8HBqNBitWrKhwn40bN6Jjx47Q6/Vo0qQJFi1aJPm4TDiqqe+//x6JiYmYMmUKUlJS0K5dO8THx+PSpUueDo1I3QRB8uBtQRCwbsMGp4/dyMjArt27pcXAweNUDqnv//n5+bjjjjtw6tQpLFu2DEeOHMGCBQvQoEGDKo6c3EHq9V+yZAlefvllTJkyBYcOHcLnn3+O77//Hq+88koVR07ukJWVhXbt2iEpKUlU+ZMnT6J///7o06cP9uzZg3HjxuGxxx7D2rVrJR2X0+JWU++++y7GjBljXzxw/vz5WLVqFRYuXIiXX37ZXs5kMola3l7MDDzFCUW3fAiCAJPJ5Pb6XdlHavmcnBz7bx+fil9KcsejtPhd2Udp5+BS+fx8+3oYWcX2N2Vl2Z/3NrZvCA8cOoTUtDT4+fkh30nvxNr169GqZUv7+YpSUIBsFxIPpb2WJZ0ziSL2/d9m4cKFuHbtGrZu3QpfX18AQKNGjaoyZHIjqdd/69at6NGjB4YNGwbAeu2HDh2Kbdu2VWnc5B59+/ZF3759RZefP38+GjdujDlz5gAAWrZsib/++gvvvfce4uPjRdfDHo5qKD8/H7t27UJcXBwAa9dq69atkZWVhTfffBNRUVHo0qWLh6MkUikXGvnbdu4EALRv08bp41tdud2RvRzkRMn3fwDQarWIi4tDcnKy031WrlyJ6OhoJCQkICQkBK1bt8aMGTNgNpurKmxyE1euf0xMDHbt2mW/7erEiRNYvXo1+vXrVyUxk2clJyc7PF8AID4+vsznS1nYw1ENXblyBWaz2b5SeUJCAhISEjBhwgT8+eef2LZtGzIyMhAcHAyj0Qij0Si6brFlbQugaTQaWeqvzD5iy9vuq/f395f1HNQevyv7KO0cRJc3mwGz2V6++EJ/xsBABJZYmM8W/649ewAAMd26YfuuXaWq3bZzJwICAqTHIwiAVivrNZP7tZyRkSG6TqpYyfd/m5CQEBw+fNjpPidOnMAff/yBhx9+GKtXr8a///6Lp59+GgUFBZgyZUpVhE1u4sr1HzZsGK5cuYKePXtCEAQUFhbiySef5C1V1URqaqrT50tGRgZycnLg7y9uxSj2cBARuYuLvQq2JKNrx46lHjMYDEhPT8eRY8eqLB6i4iwWC+rXr49PP/0UnTp1woMPPohXX30V8+fP93RoVAU2btyIGTNmYN68eUhJScHy5cuxatUqTJ8+3dOhkYow4aiG6tatC51Oh7S0NIftaWlpCA0N9VBURConCJKnwgWAjMxMHC9aYbxThw6lHm/fti0AIJm3VZEbuPL+HxYWhmbNmkFXNDYJsN7HnZqa6nTMESmXK9f/tddew/Dhw/HYY4+hTZs2uOeeezBjxgzMnDkTFoulKsImDwoNDXX6fAkKChLduwEw4aiW/Pz80KlTJ6wvtqCYxWLB+vXrER0d7cHIiFSssNCl9S+O/PsvAKB506aoVatWqcc7t28P4L9xHpIIgvU2L6Iirrz/9+jRA//++69D4/Lo0aMICwuDn5+f7DGT+7hy/bOzs0tNHmNLPktOhEHeJzo62uH5AgDr1q2T3F5kwlFNJSYmYsGCBfjyyy9x6NAhPPXUU8jKyrLPWvH44497OEIilXGxYW9LOLp17uz0cVuvh0s9HIBLvS7k3Sp6/x8xYgQmTpxoL//UU0/h2rVreO6553D06FGsWrUKM2bMQEJCgqdOgSpB6vUfOHAgPv74Y3z33Xc4efIk1q1bh9deew0DBw506PUidTCZTNizZw/2FI0dPHnyJPbs2YMzZ84AACZOnIgRI0bYyz/55JM4ceIEJkyYgMOHD2PevHn44YcfMH78eEnH5aDxaurBBx/E5cuXMXnyZKSmpqJ9+/ZYs2aNfWDQuXPnPBwhkcq42LC3jc2I7trV6eOdixKOA4cOIdNkQg2pA8CZcFAJFb3/nzlzxuEb7YiICKxduxbjx49H27Zt0aBBAzz33HN46aWXPHUKVAlSr/+kSZOg0WgwadIknD9/HvXq1cPAgQPx5ptveuoUqBJ27tyJPn362P+fmJgIABg5ciQWLVqEixcv2pMPAGjcuDFWrVqF8ePH4/3338dNN92Ezz77TNKUuAATjmpt7NixGDt2rNPHVq9ejeDg4CqOiEilLBaXbqeyWCz2Ho7uZUxFHVK/PhrefDNOnzmDXf/f3p3HR1Xd/x9/ZyGTRCCAMQsxLC4IFAybwUBVsKkR+VKpSykiIEVa7MQK6U8kVUmtlaCV5atGUhBEWxFcABUQpaGBIkFkSQsCUQEJIAmmSBISSCA5vz+A+TJkIdudmSSv5+NxH3XOnHPv53QW5p277dihQbfcUrfaLrpiFlDd9396enqFtpiYGG3evNniquAqtXn9fX19lZSUxBXJmohBgwZVeyhcZXcRHzRokHbU9ga0l+CQKgCorzruRdj79dcqPnVKVwQGqkf37lX2u/n84VZfbN9ep+2wlwMA4E4EDgCorzqev7Hl/IngfXr1qvZO6RcOt9pclxPHJU4cBwC4FYEDAOqrjnsQLlx5qqoTxi8YePPNks4FlDpdhpI9HAAANyJwAEB91GPvwYUrT10ucPS68UZdccUVOpGfrz1ZWbXfkDHnzuUAAMANCBwAUB913Htw6PBhHTx0SN5eXpcNHL6+vhpw/rCqjRkZddoeezkAAO5C4ACA+qjjHo5169dLkq675poaXep28K23SpLWf/ZZnbbHeRwAAHchcABAXdXjTt6frlsnSYrq0aNG/WMHDZIkbfjsM52ty94KAgcAwE0IHABQV3U8L6KsrEyfpKVJkvpGRdVoTJ9evdS2bVvlFxRo0+ef136jxtTpXiEAANQXgQMA6qqOgWPT55/rv8ePq01QkLpef32Nxvj4+Cju9tslSR+sWlWn7XLiOADAHQgcAFBXdfwBv/T99yVJd8bGVnv/jUsNu+suSdK7K1bU7fK4HFYFAHADAgcA1FUdfvSfPn1ab7/3niTp3rvvrtXYOwYPVutWrXTo8GHHSee1wh4OAIAbEDgAoK7q8AP+rXfe0fEfflDk1Vdr8C231Gqsv7+/RtxzjyTppdTUWm+bwAEAcAcCBwDURR1upnfq1Cn9+S9/kSQ99sgj8vHxqfVmH3n4YXl5eemjjz/W5vM3DqwxAgcAwA0IHABQF7X88W6M0eNPPaVvDx5URPv2emT8+Dpt9vprr9VDo0ZJkn5ltys/P782RXClKgCAyxE4AKAuahA4zEU/7n8zaZJS5s+XJP3vjBkyxqiouFinT59WUXGxioqKVFRcXGEdjufOP19UXKykxESFhYZqT1aWbh0yRBs3bVJRUZHT9qrEieMAABer+eVRAJz7kVlWJp09KxUXy6ekRDp9WjpzRvLxkbzJ8M1GDQJH8alTjv9e/M47jv++b8yYGm+m82VuDPifXbt0y513SpJOHj2qK664ovoVclgVAMDF+HUE1ERZmXTqlFRUdC5gnD0rlZXJq7z8XNg4ffrcc6dO8YOuuWisr3NjrduFNmzYoGHDhql9+/by8vLSihUrLjsmPT1dffr0kc1m03XXXadFixZZXicANBbs4QAup6REKi2tWd+zZ88tfn6SzWZtXXCvGvxwD77ySmWkpcnm56frr71WXl5eTs+fyM/X5owM3RwTozZBQfr+v/9V5x49tFdSH0l7JR3YtUtXXXmlJOlkUZEkqeUlezH+e/y4vL29FRgYePm6OYfjsoqKihQVFaVf/epXuuf8VcGqc+DAAQ0dOlQTJ07UW2+9pbS0ND388MMKDw9XXFycCyoGAM9G4ACqYsy5PRZ1Oea9tPTcuIAA6ZIfmWgiavDD3dvbWzffdFOVz585c0b+/v66IjBQV1xxhQICAnRs/36nPle2ayfv84fqXThH49LDpi57GNXF2MNxWUOGDNGQIUNq3D81NVWdO3fWzJkzJUndunXTxo0bNXv2bAIHAIjAgSpcuItxTa+Ac/LkSadxNV1/eXm5CgoKGnz9dRnj1N8Yqbi42h9nBYWFKi4uVkFhoeMHYQX5+VJgoOTlVb96LOhfUFBwrv6Cgqrrr8f66zLG0+ZQbf/Cwor9z++BqHH9lbyHbH5+ldZQl/VX2t/Lq9oQ7Wmf5QvrrNOd1V0kIyNDsbGxTm1xcXGaNGlSlWNKSkpUUlLieFxeXq7jx4/ryiuvrLAnDABcxRijwsJCtW/fvkb/rtYUgQOVOnr0qCSpQ4cOlm4nJydHQUFBlm4DgPWs/iwfPXpUbdq0sWz99ZGTk6PQ0FCnttDQUBUUFOjUqVMKCAioMCY5OVnPPPOMq0oEgFo5dOiQrr766gZbH4EDlWrZsqWkc2+41q1bu7kaAM1VQUGBIiMjHd9JTUViYqISEhIcj/Pz89WhQwe+cwG41YXv3FatWjXoegkcUEpKilJSUpzaSs+fJN26dWv+8QPgdp58mFFYWJhyc3Od2nJzc9W6detK925Iks1mk62SC0vwnQvAEzT0dy6BA7Lb7bLb7U5thw8fVmRkpJsqAoDGIyYmRqtXr3ZqW7t2rWJiYtxUEQB4Fu7DAQDARU6ePKnMzExlZmZKOnfZ28zMTGVnZ0s6dzjUmItu3jhx4kTt379fU6ZM0d69e/Xqq6/qnXfe0eTJk91RPgB4HAIHAAAX2bp1q3r37q3evXtLkhISEtS7d29NmzZN0rkT2C+ED0nq3LmzVq1apbVr1yoqKkozZ87Ua6+9xiVxAeA8AoeHS05O1k033aRWrVopJCREw4cPV1ZW1mXHvfvuu+ratav8/f3Vs2fPCrv7AQCVGzRokIwxFZYLdw9ftGiR0tPTK4zZsWOHSkpKtG/fPj300EMurxsAPBWBw8OtX79edrtdmzdv1tq1a3XmzBndcccdKjp/ff3KbNq0SSNHjtT48eO1Y8cODR8+XMOHD9euXbtcWDkAAAAgeRlTg9vlwmN8//33CgkJ0fr163XrrbdW2mfEiBEqKirSypUrHW0333yzevXqpdTU1Bpt58JJ4/n5+VwxBYDbFBQUKCgoqMGvCe9pLsyT71wA7mTVdxFXqWpkLtz5u127dlX2ycjIcLq+u3TurrcrVqyocsyld70trOQuynC97Oxs5eXl1XpccHCw5TdtBAAAqAkCRyNSXl6uSZMmaeDAgerRo0eV/aq6621OTk6VY7jrrefJzs5Wt643qPjU6VqPDQzw1569WYQOAADgdgSORsRut2vXrl3auHFjg6/70rveHjlyRN27d2/w7aDm8vLyVHzqtP7+gNQtpObj9hyTHlx8Wnl5eQQOAADgdgSORiI+Pl4rV67Uhg0bLnscc1V3vQ0LC6tyzKV3vS0oKKhfwWgw3UKkPk330HUAANDEcZUqD2eMUXx8vJYvX65169apc+fOlx0TExOjtLQ0pzbuegsAAAB3YA+Hh7Pb7Vq8eLE++OADtWrVynEeRlBQkAICAiRJY8aMUUREhJKTkyVJjz32mG677TbNnDlTQ4cO1ZIlS7R161bNmzfPbfMAAABA88QeDg83d+5c5efna9CgQQoPD3csS5cudfTJzs7W0aNHHY8HDBigxYsXa968eYqKitJ7772nFStWVHuiOQAAAGAF9nB4uJrcJuXSO95K0v3336/777/fgooAAACAmmMPBwAAAADLEDgAAAAAWIbAAQAAAMAyBA4AAAAAliFwAAAAALAMgQMAAACAZQgcAAAAACxD4AAAAABgGQIHAAAAAMsQOAAAAABYhsABAAAAwDIEDg+XnJysm266Sa1atVJISIiGDx+urKwspz6nT5+W3W7XlVdeqZYtW+ree+9Vbm6umyoGAAAA/g+Bw8OtX79edrtdmzdv1tq1a3XmzBndcccdKioqcvSZPHmyPvroI7377rtav369vvvuO91zzz1urBoAAAA4x9fdBaB6a9ascXq8aNEihYSEaNu2bbr11luVn5+vBQsWaPHixbr99tslSa+//rq6deumzZs36+abb3ZH2QAAAIAkAkejk5+fL0lq166dJGnbtm06c+aMYmNjHX26du2qDh06KCMjo9LAUVJSopKSkmq3U1hY2IBVAwAAoLkicDQi5eXlmjRpkgYOHKgePXpIknJycuTn56c2bdo49Q0NDVVOTk6l60lOTtYzzzxjdbkAAAAAgaMxsdvt2rVrlzZu3Fiv9SQmJiohIaHaPkeOHFH37t3rtR0AAACAwNFIxMfHa+XKldqwYYOuvvpqR3tYWJhKS0t14sQJp70cubm5CgsLq3RdNptNNput2u0VFBQ0SN0AAABo3rhKlYczxig+Pl7Lly/XunXr1LlzZ6fn+/btqxYtWigtLc3RlpWVpezsbMXExLi6XAAAAMAJezg8nN1u1+LFi/XBBx+oVatWjvMygoKCFBAQoKCgII0fP14JCQlq166dWrdurUcffVQxMTFcoQoAAABuR+DwcHPnzpUkDRo0yKn99ddf10MPPSRJmj17try9vXXvvfeqpKREcXFxevXVV11cKQAAAFARh1R5OGNMpcuFsCFJ/v7+SklJ0fHjx1VUVKRly5ZVef4GAODyUlJS1KlTJ/n7+6t///7asmVLtf3nzJmjG264QQEBAYqMjNTkyZN1+vRpF1ULAJ6NwAEAwEWWLl2qhIQEJSUlafv27YqKilJcXJyOHTtWaf/Fixdr6tSpSkpK0p49e7RgwQItXbpUf/jDH1xcOQB4JgIHAAAXmTVrliZMmKBx48ape/fuSk1NVWBgoBYuXFhp/02bNmngwIF64IEH1KlTJ91xxx0aOXLkZfeKAEBzQeAAAOC80tJSbdu2TbGxsY42b29vxcbGKiMjo9IxAwYM0LZt2xwBY//+/Vq9erXuuuuuKrdTUlKigoICpwUAmipOGgcA4Ly8vDyVlZUpNDTUqT00NFR79+6tdMwDDzygvLw8/fjHP5YxRmfPntXEiROrPaQqOTlZzzzzTIPWDgCeij0cAADUQ3p6uqZPn65XX31V27dv17Jly7Rq1So9++yzVY5JTExUfn6+Yzl06JALKwYA12IPBwAA5wUHB8vHx0e5ublO7bm5uVVe/e/pp5/W6NGj9fDDD0uSevbsqaKiIv3617/Wk08+KW/vin/bs9lsstlsDT8BAPBA7OEAAOA8Pz8/9e3bV2lpaY628vJypaWlKSYmptIxxcXFFUKFj4+PpHOXNgeA5o49HAAAXCQhIUFjx45Vv379FB0drTlz5qioqEjjxo2TJI0ZM0YRERFKTk6WJA0bNkyzZs1S79691b9/f33zzTd6+umnNWzYMEfwAIDmjMABAMBFRowYoe+//17Tpk1TTk6OevXqpTVr1jhOJM/Oznbao/HUU0/Jy8tLTz31lI4cOaKrrrpKw4YN03PPPeeuKQCARyFwAABwifj4eMXHx1f6XHp6utNjX19fJSUlKSkpyQWVAUDjwzkcAAAAACxD4AAAAABgGQIHAAAAAMsQOAAAAABYhsABAAAAwDIEDgAAAACWIXAAAAAAsAyBAwAAAIBlCByNwIYNGzRs2DC1b99eXl5eWrFiRbX909PT5eXlVWHJyclxTcEAAADAeQSORqCoqEhRUVFKSUmp1bisrCwdPXrUsYSEhFhUIQAAAFA5X3cXgMsbMmSIhgwZUutxISEhatOmTcMXBAAAANQQeziasF69eik8PFw//elP9dlnn1Xbt6SkRAUFBY6lsLDQRVUCAACgKSNwNEHh4eFKTU3V+++/r/fff1+RkZEaNGiQtm/fXuWY5ORkBQUFOZbu3bu7sGIAAAA0VRxS1QTdcMMNuuGGGxyPBwwYoH379mn27Nn629/+VqF/SkqKlixZoi5dujjazpw5owMHDrikXgAAADRdBI5mIjo6Whs3bqz0ObvdLrvd7tR2+PBhRUZGuqI0AAAANGEcUtVMZGZmKjw83N1lAAAAoJlhD0cjcPLkSX3zzTeOxwcOHFBmZqbatWunDh06KDExUUeOHNGbb74pSZozZ446d+6sH/3oRzp9+rRee+01rVu3Tp9++qm7pgAAAIBmisDRCGzdulWDBw92PE5ISJAkjR07VosWLdLRo0eVnZ3teL60tFS///3vdeTIEQUGBurGG2/UP/7xD6d1AAAAAK5A4GgEBg0aJGNMlc8vWrTI6fGUKVM0ZcoUi6sCAAAALo9zOAAAAABYhsABAAAAwDIEDgAAAACWIXAAAAAAsAyBAwAAAIBlCBwAAAAALEPgAAAAAGAZAgcAAAAAyxA4AAAAAFiGwAEAAADAMgQOAAAAAJYhcAAAAACwDIEDAAAAgGUIHAAAXCIlJUWdOnWSv7+/+vfvry1btlTb/8SJE7Lb7QoPD5fNZlOXLl20evVqF1ULAJ7N190FAADgSZYuXaqEhASlpqaqf//+mjNnjuLi4pSVlaWQkJAK/UtLS/XTn/5UISEheu+99xQREaGDBw+qTZs2ri8eADwQgQMAgIvMmjVLEyZM0Lhx4yRJqampWrVqlRYuXKipU6dW6L9w4UIdP35cmzZtUosWLSRJnTp1cmXJAODROKQKAIDzSktLtW3bNsXGxjravL29FRsbq4yMjErHfPjhh4qJiZHdbldoaKh69Oih6dOnq6ysrMrtlJSUqKCgwGkBgKaKwAEAwHl5eXkqKytTaGioU3toaKhycnIqHbN//3699957Kisr0+rVq/X0009r5syZ+vOf/1zldpKTkxUUFORYIiMjG3QeAOBJCBwAANRDeXm5QkJCNG/ePPXt21cjRozQk08+qdTU1CrHJCYmKj8/37EcOnTIhRUDgGtxDgcAAOcFBwfLx8dHubm5Tu25ubkKCwurdEx4eLhatGghHx8fR1u3bt2Uk5Oj0tJS+fn5VRhjs9lks9katngA8FDs4QAA4Dw/Pz/17dtXaWlpjrby8nKlpaUpJiam0jEDBw7UN998o/LyckfbV199pfDw8ErDBgA0NwQOAAAukpCQoPnz5+uNN97Qnj179Mgjj6ioqMhx1aoxY8YoMTHR0f+RRx7R8ePH9dhjj+mrr77SqlWrNH36dNntdndNAQA8CodUAQBwkREjRuj777/XtGnTlJOTo169emnNmjWOE8mzs7Pl7f1/f6+LjIzUJ598osmTJ+vGG29URESEHnvsMT3xxBPumgIAeBQCBwAAl4iPj1d8fHylz6Wnp1doi4mJ0ebNmy2uCgAaJwIHlJKSopSUFKe20tJSN1UDAACApoTAAdnt9grHGh8+fJjrwgMAAKDeOGkcAAAAgGUIHAAAAAAsQ+AAAAAAYBkCBwAAAADLEDgAAAAAWIbAAQAAAMAyBA4AAAAAliFwAAAAALAMgQMAAACAZQgcAAAAACxD4AAAAABgGQIHAAAAAMsQOAAAAABYhsDRCGzYsEHDhg1T+/bt5eXlpRUrVlx2THp6uvr06SObzabrrrtOixYtsrxOAAAA4FIEjkagqKhIUVFRSklJqVH/AwcOaOjQoRo8eLAyMzM1adIkPfzww/rkk08srhQAAABw5uvuAnB5Q4YM0ZAhQ2rcPzU1VZ07d9bMmTMlSd26ddPGjRs1e/ZsxcXFWVUmAAAAUAGBownKyMhQbGysU1tcXJwmTZpU5ZiSkhKVlJQ4HhcWFlpVnkfIzs5WXl5erccFBwerQ4cOFlTU8Pbs2VPrMSUlJbLZbC4Z58ptSXV/7er6XnF1nQAAeCoCRxOUk5Oj0NBQp7bQ0FAVFBTo1KlTCggIqDAmOTlZzzzzjKtKdKvs7Gx163qDik+drvXYwAB/7dmb5dE/CI8WSN5e0oMPPljrsT5eUpmp/TbrMs6V25Lq9trV573iyjoBAPBkBA5IkhITE5WQkOB4fOTIEXXv3t2NFVknLy9PxadO6+8PSN1Caj5uzzHpwcWnlZeX59E/Bk+clsqNaj2/1Xulp9e4ZpwrtyXV/bWr63vF1XUCAODJCBxNUFhYmHJzc53acnNz1bp160r3bkiSzWZzOvyjoKDA0ho9QbcQqc/V7q7COrWd355jrhvnym01hMZSJwAAnoirVDVBMTExSktLc2pbu3atYmJi3FQRAAAAmisCRyNw8uRJZWZmKjMzU9K5y95mZmYqOztb0rnDocaMGePoP3HiRO3fv19TpkzR3r179eqrr+qdd97R5MmT3VE+AAAAmjECRyOwdetW9e7dW71795YkJSQkqHfv3po2bZok6ejRo47wIUmdO3fWqlWrtHbtWkVFRWnmzJl67bXXuCQuAAAAXI5zOBqBQYMGyZiqL3dT2V3EBw0apB07dlhYFQAAAHB57OEAAAAAYBkCBwAAAADLEDgAAAAAWIbAAQAAAMAynDQOpaSkKCUlxamttLTUTdUAAACgKWEPB2S327V7926nJT093d1lAYDbpKSkqFOnTvL391f//v21ZcuWGo1bsmSJvLy8NHz4cGsLBIBGhMABAMBFli5dqoSEBCUlJWn79u2KiopSXFycjh07Vu24b7/9Vv/v//0/3XLLLS6qFAAaBwIHAAAXmTVrliZMmKBx48ape/fuSk1NVWBgoBYuXFjlmLKyMo0aNUrPPPOMrrnmGhdWCwCej8ABAMB5paWl2rZtm2JjYx1t3t7eio2NVUZGRpXj/vSnPykkJETjx4+v0XZKSkpUUFDgtABAU0XgAADgvLy8PJWVlSk0NNSpPTQ0VDk5OZWO2bhxoxYsWKD58+fXeDvJyckKCgpyLJGRkfWqGwA8GYEDAIA6Kiws1OjRozV//nwFBwfXeFxiYqLy8/Mdy6FDhyysEgDci8viAgBwXnBwsHx8fJSbm+vUnpubq7CwsAr99+3bp2+//VbDhg1ztJWXl0uSfH19lZWVpWuvvbbCOJvNJpvN1sDVA4BnYg8HAADn+fn5qW/fvkpLS3O0lZeXKy0tTTExMRX6d+3aVTt37lRmZqZj+dnPfqbBgwcrMzOTQ6UAQOzhAADASUJCgsaOHat+/fopOjpac+bMUVFRkcaNGydJGjNmjCIiIpScnCx/f3/16NHDaXybNm0kqUI7ADRXBA4AAC4yYsQIff/995o2bZpycnLUq1cvrVmzxnEieXZ2try9OUAAAGqKwAEAwCXi4+MVHx9f6XPp6enVjl20aFHDFwQAjRh/ogEAAABgGQIHAAAAAMsQOAAAAABYhsABAAAAwDIEDgAAAACWIXAAAAAAsAyBAwAAAIBlCBwAAAAALEPgAAAAAGAZAgcAAAAAyxA4AAAAAFiGwAEAAADAMgQOAAAAAJYhcAAAAACwDIGjCUlJSVGnTp3k7++v/v37a8uWLe4uCQAAAM0cgaOJWLp0qRISEpSUlKTt27crKipKcXFxOnbsmLtLAwAAQDNG4GgiZs2apQkTJmjcuHHq3r27UlNTFRgYqIULF7q7NAAAADRjvu4uAPVXWlqqbdu2KTEx0dHm7e2t2NhYZWRkVOhfUlKikpKSatdZWFjY4HUCAACg+SFwNAF5eXkqKytTaGioU3toaKj27t1boX9ycrKeeeYZV5UHAACAZoxDqpqhxMRE5efnV7vs3r3b3WUCAACgCWAPRxMQHBwsHx8f5ebmOrXn5uYqLCysQn+bzSabzVbtOgsKChq0RgAAADRP7OFoAvz8/NS3b1+lpaU52srLy5WWlqaYmBg3VgYAAIDmjj0cTURCQoLGjh2rfv36KTo6WnPmzFFRUZHGjRvn7tIAAADQjBE4GomUlBT95S9/UU5OjqKiovTyyy8rOjra8fyIESP0/fffa9q0aTpy5IjKysokyXFIlc1m0+nTp91SOwAAAJovDqlqBGp6U7/4+HgdPHhQr732mlq3bq2jR486loMHD7qpegAAADRnBI5GoC439fPy8lJYWJhjufSSuQAAAIArEDg83IWb+sXGxjraqrup3wUnT55Ux44dFRkZqbvvvltffvlltdspKSlRQUGBY+HGfwAAAGgIBA4PV91N/XJyciodc8MNN2jhwoX64IMP9Pe//13l5eUaMGCADh8+XOV2kpOTFRQU5Fi6d+/eoPMAAABA80TgaIJiYmI0ZswY9erVS7fddpuWLVumq666Sn/961+rHHPpzQC58R8AAAAaAlep8nC1valfZVq0aKHevXvrm2++qbLPpTcD5MZ/AAAAaAjs4fBwDXFTv7KyMu3cuVPh4eFWlQkAAABUisDRCCQkJGj+/Pl64403tGfPHj3yyCNON/UbM2aMEhMTHf3/9Kc/6dNPP9X+/fu1fft2Pfjggzp48KAefvhhd00BABqVlJQUderUSf7+/urfv7+2bNlSZd/58+frlltuUdu2bdW2bVvFxsZW2x8AmhsOqWoELr6pX05Ojnr16qU1a9Y4TiTPzs6Wt/f/ZccffvhBEyZMUE5Ojtq2bau+fftq06ZNnAgOADVw4d5Hqamp6t+/v+bMmaO4uDhlZWUpJCSkQv/09HSNHDlSAwYMkL+/v55//nndcccd+vLLLxUREeGGGQCAZyFwNBLx8fGKj4+v9Ln09HSnx7Nnz9bs2bNdUBUAND0X3/tIklJTU7Vq1SotXLhQU6dOrdD/rbfecnr82muv6f3331daWprGjBnjkpoBwJNxSBUAAOfV9d5HFysuLtaZM2fUrl07q8oEgEaFPRwAAJxX3b2P9u7dW6N1PPHEE2rfvr1TaLlUSUmJSkpKHI+5MiCApow9HAAANJAZM2ZoyZIlWr58ufz9/avsd+nNViMjI11YJQC4FoEDAIDz6nPvoxdffFEzZszQp59+qhtvvLHavpfebPXQoUP1rh0APBWBAwCA8+p676MXXnhBzz77rNasWaN+/fpddjs2m02tW7d2WgCgqeIcDgAALpKQkKCxY8eqX79+io6O1pw5cyrc+ygiIkLJycmSpOeff17Tpk3T4sWL1alTJ+Xk5EiSWrZsqZYtW7ptHgDgKQgcAABcpLb3Ppo7d65KS0t13333Oa0nKSlJf/zjH11ZOgB4JAIHAACXqM29j7799lvrCwKARoxzOAAAAABYhsABAAAAwDIEDgAAAACWIXAAAAAAsAyBAwAAAIBlCBwAAAAALMNlcaGUlBSlpKQ4tZWWlrqpGgAAADQlBA7IbrfLbrc7tR0+fFiRkZFuqggAAABNBYdUAQAAALAMgQMAAACAZQgcAAAAACxD4AAAAABgGQIHAAAAAMsQOAAAAABYhsABAAAAwDIEDgAAAACWIXAAAAAAsAyBAwAAAIBlCBwAAAAALEPgAAAAAGAZAgcAAAAAyxA4AAAAAFiGwAEAAADAMgQOAAAAAJYhcAAAAACwDIEDAAAAgGUIHAAAAAAsQ+AAAAAAYBkCBwAAAADLEDgAAAAAWIbAAQAAAMAyBA4AAAAAliFwNBIpKSnq1KmT/P391b9/f23ZsqXa/u+++666du0qf39/9ezZU6tXr3ZRpQDQ+PGdCwANh8DRCCxdulQJCQlKSkrS9u3bFRUVpbi4OB07dqzS/ps2bdLIkSM1fvx47dixQ8OHD9fw4cO1a9cuF1cOAI0P37kA0LAIHI3ArFmzNGHCBI0bN07du3dXamqqAgMDtXDhwkr7/+///q/uvPNOPf744+rWrZueffZZ9enTR6+88oqLKweAxofvXABoWL7uLgDVKy0t1bZt25SYmOho8/b2VmxsrDIyMiodk5GRoYSEBKe2uLg4rVixosrtlJSUqKSkxPE4Pz9fklRQUFCrenNycpSTk1OrMdK5OZWXl7tkXFZWliRp2xHpZMllOl88Lu/c/27btk0nT56s1TZdWeee83+E9eRxrq6xrq+dq18DV7/H3DEuLCxMYWFhNe5/4TuoLtuqi8b2nQsADenCd5AxpmFXbODRjhw5YiSZTZs2ObU//vjjJjo6utIxLVq0MIsXL3ZqS0lJMSEhIZX2f+WVV0xwcLCRxMLCwuKRy5YtWxrmS/UyXPGda4wxSUlJbv//lIWFhaWqZd++ffX/Qr0Iezggu92uhx9+2OmvbWfPnlXv3r31r3/9S97eNTvybtCgQUpPT69R38LCQnXv3l27d+9Wq1atGnz9NRlTWQ2eNgcr+7uifqu30dxfg5qMb+jPzaWsfg3Ky8t16623qnfv3jXq31gkJiY67RU5ceKEOnbsqOzsbAUFBbmxMmsVFBQoMjJShw4dUuvWrd1djqWay1yZZ9OSn5+vDh06qF27dg26XgKHhwsODpaPj49yc3Od2nNzc6s8LCEsLKxW/SXJZrPJZrM5tV1xxRXq0KFDjWv18/PT1VdfXaO+F3bZRURE1PiDW5v112RMZTV42hys7O+K+q3eRnN/DWoyvqE/N5dyxWsQGBgoX1/X/HPlzu9cSQoKCmrSP2YuaN26dbOYp9R85so8m5aa/rG5xutr0LWhwfn5+alv375KS0tztJWXlystLU0xMTGVjomJiXHqL0lr166tsj8A4By+cwGg4bGHoxFISEjQ2LFj1a9fP0VHR2vOnDkqKirSuHHjJEljxoxRRESEkpOTJUmPPfaYbrvtNs2cOVNDhw7VkiVLtHXrVs2bN8+d0wCARoHvXABoWASORmDEiBH6/vvvNW3aNOXk5KhXr15as2aNQkNDJUnZ2dlOu74GDBigxYsX66mnntIf/vAHXX/99VqxYoV69Ojhril4JJvNpqSkpEoPa0DDsNlsCg4O5v9ji/AetoY7vnOby2vZXOYpNZ+5Ms+mxbJ5Nugp6GhSXnnlFcv65+fnG0kmPz/fsnrqMsbT5tDY66/LGE+bQ3Orv7ZjPHUOAADP4WVMQ19oF7i8goICBQUFKT8/v9GefNXY59DY65ca/xwae/1S05gDAMBanDQOAAAAwDIEDgAAAACWIXDALZrCyVeNfQ6NvX6p8c+hsdcvNY05AACsReCAW9hsNv3xj3/0mB8p3377rcaPH6/OnTsrICBA1157rZKSklRaWlrlGE+bQ21ZXX9KSoo6deokf39/9e/fX1u2bGnwbTTH1yA5OVk33XSTWrVqpZCQEA0fPlxZWVkWVlm9xv4auFJtPxPvvvuuunbtKn9/f/Xs2VOrV692UaX1U5t5zp8/X7fccovatm2rtm3bKjY21pLvCivU9TtuyZIl8vLy0vDhw60tsAHVdq4nTpyQ3W5XeHi4bDabunTp0ijev7Wd55w5c3TDDTcoICBAkZGRmjx5sk6fPu2iautmw4YNGjZsmNq3by8vLy+tWLHismPS09PVp08f2Ww2XXfddVq0aFHtN+zus9YBT/Dxxx+bhx56yHzyySdm37595oMPPjAhISHm97//vbtLa5SWLFli/Pz8zMKFC82XX35pJkyYYNq0aWNyc3PdXVqjFxcXZ15//XWza9cuk5mZae666y7ToUMHc/LkSXeXhmrU9jPx2WefGR8fH/PCCy+Y3bt3m6eeesq0aNHC7Ny508WV105t5/nAAw+YlJQUs2PHDrNnzx7z0EMPmaCgIHP48GEXV147df2OO3DggImIiDC33HKLufvuu11TbD3Vdq4lJSWmX79+5q677jIbN240Bw4cMOnp6SYzM9PFlddObef51ltvGZvNZt566y1z4MAB88knn5jw8HAzefJkF1deO6tXrzZPPvmkWbZsmZFkli9fXm3//fv3m8DAQJOQkGB2795tXn75ZePj42PWrFlTq+0SOIAqvPDCC6Zz587uLqNRio6ONna73fG4rKzMtG/f3iQnJ7uxqqbp2LFjRpJZv369u0tBNWr7mfjFL35hhg4d6tTWv39/85vf/MbSOuurvp/9s2fPmlatWpk33njDqhIbRF3mefbsWTNgwADz2muvmbFjxzaawFHbuc6dO9dcc801prS01FUlNojaztNut5vbb7/dqS0hIcEMHDjQ0jobUk0Cx5QpU8yPfvQjp7YRI0aYuLi4Wm2LQ6qAKuTn56tdu3buLqPRKS0t1bZt2xQbG+to8/b2VmxsrDIyMtxYWdOUn58vSbxXPVhdPhMZGRlO/SUpLi7Ooz9DDfHZLy4u1pkzZzz6/VzXef7pT39SSEiIxo8f74oyG0Rd5vrhhx8qJiZGdrtdoaGh6tGjh6ZPn66ysjJXlV1rdZnngAEDtG3bNsdhV/v379fq1at11113uaRmV2mo7yLuNA5U4ptvvtHLL7+sF1980d2lNDp5eXkqKytz3JX5gtDQUO3du9dNVTVN5eXlmjRpkgYOHFiru1rDterymcjJyam0f05OjmV11ldDfPafeOIJtW/fvsIPHE9Sl3lu3LhRCxYsUGZmpgsqbDh1mev+/fu1bt06jRo1SqtXr9Y333yj3/72tzpz5oySkpJcUXat1WWeDzzwgPLy8vTjH/9YxhidPXtWEydO1B/+8AdXlOwyVX0XFRQU6NSpUwoICKjRetjDAbdwxQnFkjR16lR5eXlVu1z6ZXLkyBHdeeeduv/++5WXl3fZE3RPnz4tu92uK6+8Ui1bttS9996r3NxcS+ZTXzNmzJCXl5cmTZrkaGsM9R85ckQPPvigrrzySgUEBKhnz57aunWr43ljjKZNm6bw8HAFBAQoNjZWX3/9tRsr/j9lZWV6+umnnS5I8Oyzz8pcdM/VutZvt9u1a9cuLVmypEFrvtxJhTWp9/jx4xo1apRat26tNm3aaPz48Tp58mSD1ommZcaMGVqyZImWL18uf39/d5fTYAoLCzV69GjNnz9fwcHB7i7HcuXl5QoJCdG8efPUt29fjRgxQk8++aRSU1PdXVqDSk9P1/Tp0/Xqq69q+/btWrZsmVatWqVnn33W3aV5JAIHXG7p0qVKSEhQUlKStm/frqioKMXFxenYsWMNvq3f//732rNnT7XLNddc4+j/3XffafDgwRowYIDmzZun9evXy263a/PmzVq7dq3OnDmjO+64Q0VFRY4xkydP1kcffaR3331X69ev13fffad77rmnwedSX1988YX++te/6sYbb3Rqb+j6g4OD5ePjUyG05ObmKiwsrNbr++GHHzRw4EC1aNFCH3/8sXbv3q2ZM2eqbdu2jj4vvPCCXnrpJaWmpurzzz/XFVdcobi4OI+4Wsjzzz+vuXPn6pVXXtGePXv0/PPP64UXXtDLL7/s6FOX+uPj47Vy5Ur985//1NVXX92gNRcVFSkqKkopKSmVPl+TekeNGqUvv/xSa9eu1cqVK7Vhwwb9+te/btA6G4u6fCbCwsIa7DPkKvX57L/44ouaMWOGPv300wrfUZ6mtvPct2+fvv32Ww0bNky+vr7y9fXVm2++qQ8//FC+vr7at2+fq0qvtbq8puHh4erSpYt8fHwcbd26dVNOTk61V350p7rM8+mnn9bo0aP18MMPq2fPnvr5z3+u6dOnKzk5WeXl5a4o2yWq+i5q3bp1jfduSOIqVXA9Tz2h+PDhw+b66683v/zlL83Zs2cr7XPpCbonTpwwLVq0MO+++66jz549e4wkk5GR4ZK6a6KwsNBcf/31Zu3atea2224zjz32mDHGuvqjo6NNfHy843FZWZmJiIio02v8xBNPmB//+MdVPl9eXm7CwsLMX/7yF0fbiRMnjM1mM2+//Xatt9fQhg4dan71q185td1zzz1m1KhRxpja119eXm7sdrtp3769+eqrr6wt3lQ8qbAm9e7evdtIMl988YWjz8cff2y8vLzMkSNHLK/ZE9X2M/GLX/zC/M///I9TW0xMTKM4aby2n/3nn3/etG7d2qO+My+nNvM8deqU2blzp9Ny9913m9tvv93s3LnTlJSUuLL0Wqvta5qYmGg6duxoysrKHG1z5swx4eHhltdaH7WdZ58+fcyUKVOc2hYvXmwCAgKq/A3haS79fq/MlClTTI8ePZzaRo4cWeuTxgkccKmSkhLj4+NT4Q0+ZswY87Of/cw9RZlzYeO6664zP/nJT8zhw4fN0aNHHcvFvv76ayPJcWnKtLQ0I8n88MMPTv06dOhgZs2a5aryL2vMmDFm0qRJxhjjFDisqn/JkiXGZrOZRYsWmd27d5tf//rXpk2bNiYnJ6fW6+rWrZuZNGmSue+++8xVV11levXqZebNm+d4ft++fUaS2bFjh9O4W2+91fzud7+r8xwaynPPPWc6duxosrKyjDHGZGZmmpCQEPP3v//dGFP7+h955BETFBRk0tPTnd6nxcXFltR/6T9INal3wYIFpk2bNk7Pnzlzxvj4+Jhly5ZZUqenu9xnYvTo0Wbq1KmO/p999pnx9fU1L774otmzZ49JSkpqNJfFrc08Z8yYYfz8/Mx7773n9H4uLCx01xRqpLbzvFRjukpVbeeanZ1tWrVqZeLj401WVpZZuXKlCQkJMX/+85/dNYUaqe08k5KSTKtWrczbb79t9u/fbz799FNz7bXXml/84hfumkKNFBYWmh07dpgdO3YYSWbWrFlmx44d5uDBg8YYY6ZOnWpGjx7t6H/hsriPP/642bNnj0lJSeGyuPB8R44cMZLMpk2bnNoff/xxEx0d7aaqjHn99deNpEqXC8rKyszQoUOdLnn31ltvGT8/vwrru+mmmyr85cNd3n77bdOjRw9z6tQpY4xz4LCy/pdfftl06NDB+Pn5mejoaLN58+Y6rcdmsxmbzWYSExPN9u3bzV//+lfj7+9vFi1aZIw598NMkvnuu++cxt1///0e8cVfVlZmnnjiCePl5WV8fX2Nl5eXmT59uuP52tZf1fv09ddft6T+SwNHTep97rnnTJcuXSqs66qrrjKvvvqqJXU2BtV9Jm677TYzduxYp/7vvPOO6dKli/Hz8zM/+tGPzKpVq1xccd3UZp4dO3as9P2clJTk+sJrqbav58UaU+AwpvZz3bRpk+nfv7+x2WzmmmuuMc8991yj+Kt/beZ55swZ88c//tFce+21xt/f30RGRprf/va3Ff6A52n++c9/VvqZuzC3sWPHmttuu63CmF69ehk/Pz9zzTXX1OnfGwIHXMpTA0dNTJw40XTs2NEcOnTI0ebpgSM7O9uEhISYf//73442VwWOhtKiRQsTExPj1Pboo4+am2++2Rjj+YHj7bffNldffbV5++23zX/+8x/z5ptvmnbt2jWawETgAADUFyeNw6Ua+oRiV6nqBN2wsDCVlpbqxIkTTv09ZT7btm3TsWPH1KdPH8fJiuvXr9dLL70kX19fhYaGenT90rkTELt37+7U1q1bN2VnZ0uSo05PfU89/vjjmjp1qn75y1+qZ8+eGj16tCZPnqzk5GRJnl//pWpSb1hYWIWLQJw9e1bHjx/3yDkBAKxF4IBL+fn5qW/fvkpLS3O0lZeXKy0tTTExMW6srHLGGMXHx2v58uVat26dOnfu7PR837591aJFC6f5ZGVlKTs72yPm85Of/EQ7d+5UZmamY+nXr59GjRrl+G9Prl+SBg4cWOFSxF999ZU6duwoSercubPCwsKc5lBQUKDPP//cI+ZQXFwsb2/nr1ofHx/HVUw8vf5L1aTemJgYnThxQtu2bXP0WbduncrLy9W/f3+X1wwAcDN372JB89OQJxRbrSYn6E6cONF06NDBrFu3zmzdutXExMRUOATIk1x8SJUxnl//li1bjK+vr3nuuefM119/bd566y0TGBjoOOnamHMnnrZp08Z88MEH5j//+Y+5++67TefOnR3nrbjT2LFjTUREhFm5cqU5cOCAWbZsmQkODnY6ZM3T6r/cSYU1qffOO+80vXv3Np9//rnZuHGjuf76683IkSPdMh8AgHsROOAWDXVCsdVUgxN0T506ZX7729+atm3bmsDAQPPzn/+8wtWtPMmlgaMx1P/RRx+ZHj16GJvNZrp27ep0lSpjzl2q9emnnzahoaHGZrOZn/zkJ46rQrlbQUGBeeyxx0yHDh2Mv7+/ueaaa8yTTz7pdClMT6v/cicV1qTe//73v2bkyJGmZcuWpnXr1mbcuHEef+UhAIA1vIy56Ha3AAAAANCAOIcDAAAAgGUIHAAAAAAsQ+AAAAAAYBkCBwAAAADLEDgAAAAAWIbAAQAAAMAyBA4AAAAAliFwAAAAALAMgQNAg8rKylJYWJgKCwvrtZ6bb75Z77//fgNVBQAA3IXAAcBJWVmZBgwYoHvuucepPT8/X5GRkXryySerHZ+YmKhHH31UrVq1qlcdTz31lKZOnary8vJ6rQcAALiXlzHGuLsIAJ7lq6++Uq9evTR//nyNGjVKkjRmzBj9+9//1hdffCE/P79Kx2VnZ+u6667TgQMHFBERUa8aysrKFBERoQULFmjo0KH1WhcAAHAf9nAAqKBLly6aMWOGHn30UR09elQffPCBlixZojfffLPKsCFJ77zzjqKiopzCxqJFi9SmTRutXLlSN9xwgwIDA3XfffepuLhYb7zxhjp16qS2bdvqd7/7ncrKyhzjfHx8dNddd2nJkiWWzhUAAFjL190FAPBMjz76qJYvX67Ro0dr586dmjZtmqKioqod869//Uv9+vWr0F5cXKyXXnpJS5YsUWFhoe655x79/Oc/V5s2bbR69Wrt379f9957rwYOHKgRI0Y4xkVHR2vGjBkNPjcAAOA6BA4AlfLy8tLcuXPVrVs39ezZU1OnTr3smIMHD1YaOM6cOaO5c+fq2muvlSTdd999+tvf/qbc3Fy1bNlS3bt31+DBg/XPf/7TKXC0b99ehw4dUnl5uby92SELAEBjxL/gAKq0cOFCBQYG6sCBAzp8+PBl+586dUr+/v4V2gMDAx1hQ5JCQ0PVqVMntWzZ0qnt2LFjTuMCAgJUXl6ukpKSeswCAAC4E4EDQKU2bdqk2bNna+XKlYqOjtb48eN1uWtMBAcH64cffqjQ3qJFC6fHXl5elbZdekWq48eP64orrlBAQEAdZwEAANyNwAGgguLiYj300EN65JFHNHjwYC1YsEBbtmxRampqteN69+6t3bt3N1gdu3btUu/evRtsfQAAwPUIHAAqSExMlDHGccJ2p06d9OKLL2rKlCn69ttvqxwXFxenjIwMp6tN1ce//vUv3XHHHQ2yLgAA4B4EDgBO1q9fr5SUFL3++usKDAx0tP/mN7/RgAEDqj20asiQIfL19dU//vGPetdx5MgRbdq0SePGjav3ugAAgPtw4z8ADSolJUUffvihPvnkk3qt54knntAPP/ygefPmNVBlAADAHbgsLoAG9Zvf/EYnTpxQYWGhWrVqVef1hISEKCEhoQErAwAA7sAeDgAAAACW4RwOAAAAAJYhcAAAAACwDIEDAAAAgGUIHAAAAAAsQ+AAAAAAYBkCBwAAAADLEDgAAAAAWIbAAQAAAMAyBA4AAAAAlvn/unT6fiTS1TcAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.subplot2grid((3,2), (0,0), colspan=2, rowspan = 2)\n",
    "\n",
    "seed = 103; n = 10; max_L = 100; ireal = 2; nreal = 100\n",
    "dx = -5; dy = 3; pdfy = 13.0\n",
    "\n",
    "np.random.seed(seed=seed)\n",
    "values = np.random.rand(max_L*2)*100\n",
    "X,Y = np.split(values,2)\n",
    "df = pd.DataFrame(np.vstack((X[:n],Y[:n])).T,columns = ['X','Y'])\n",
    "\n",
    "ax1 = plt.subplot(121)\n",
    "plt.scatter(df['X'],df['Y'],color='red',edgecolor='black',zorder=10,label = 'Spatial Data')\n",
    "plt.gca().add_patch(plt.Rectangle((0, 0), 100, 100, fill=False,edgecolor='black',lw=2))\n",
    "\n",
    "xx = np.arange(-3,3,0.1)\n",
    "pdf = norm.pdf(xx, loc=0, scale=1) * pdfy\n",
    "px = np.linspace(0,10,len(xx))\n",
    "\n",
    "for ipts in range(0,n):\n",
    "    if ipts == 0:\n",
    "        circle1 = plt.Circle((df['X'][ipts],df['Y'][ipts]),vrange,fill=True,color='red',edgecolor=None,\n",
    "            lw=0.0,alpha=0.05,zorder=1,label='Spatial Continuity')\n",
    "    else:\n",
    "        circle1 = plt.Circle((df['X'][ipts],df['Y'][ipts]),vrange,fill=True,color='red',edgecolor=None,\n",
    "            lw=0.0,alpha=0.05,zorder=1)\n",
    "    plt.gca().add_patch(circle1)\n",
    "    \n",
    "    plt.plot([df['X'][ipts]+dx,df['X'][ipts]+10+dx],[df['Y'][ipts]+dy,df['Y'][ipts]+dy],color='black')\n",
    "    plt.plot([df['X'][ipts]+dx,df['X'][ipts]+dx],[df['Y'][ipts]+dy,df['Y'][ipts]+5+dy],color='black')\n",
    "    plt.plot(px + df['X'][ipts] + dx,pdf + df['Y'][ipts] + dy,color = 'black',zorder=1)\n",
    "\n",
    "plt.xlim([-10,110]); plt.ylim([-10,110]); add_grid(); plt.xlabel(\"X (m)\"); plt.ylabel(\"Y (m)\") \n",
    "plt.title('Random Spatial Dataset with Spatial Correlation Range Indicated')\n",
    "\n",
    "vario = GSLIB.make_variogram(nug=0.0,nst=1,it1=1,cc1=1.0,azi1=0.0,hmaj1=vrange,hmin1=vrange)\n",
    "neffective,realizations = n_effective(df,'X','Y',seed=seed,nreal=nreal,vario=vario)\n",
    "\n",
    "for ipts in range(0,n):\n",
    "    y = df['Y'][ipts]+dy; x = df['X'][ipts]+dx+(realizations[ireal,ipts]-(-3.0))/(3.0-(-3.0))*10\n",
    "    plt.plot([x,x],[y,y+5],color='black',lw=3,zorder=1); plt.plot([x,x],[y,y+5],color='red',lw=1,zorder=3)\n",
    "\n",
    "plt.subplot2grid((3,2), (2,0))\n",
    "\n",
    "plt.hist(realizations[ireal,:],color='darkorange',edgecolor='black',bins=np.linspace(-3,3,20))\n",
    "\n",
    "plt.subplot2grid((3,2), (2,1))\n",
    "\n",
    "plt.subplots_adjust(left=0.0, bottom=0.0, right=1.0, top=1.5, wspace=0.2, hspace=0.2); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5e86c6eb67f24fbfa39644106f441258",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(Text(value='                                                      Number of Effective Data, Mic…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4ef19a72599f40bd9400780ab338cca6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output(outputs=({'output_type': 'display_data', 'data': {'text/plain': '<Figure size 640x480 with 2 Axes>', 'i…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(ui1, interactive_plot1)                            # display the interactive plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "can only concatenate list (not \"int\") to list",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[1;32mIn[23], line 1\u001b[0m\n\u001b[1;32m----> 1\u001b[0m a \u001b[38;5;241m=\u001b[39m [\u001b[38;5;241m12\u001b[39m,\u001b[38;5;241m12\u001b[39m]\u001b[38;5;241m+\u001b[39m \u001b[38;5;241m2\u001b[39m\n\u001b[0;32m      2\u001b[0m a\n",
      "\u001b[1;31mTypeError\u001b[0m: can only concatenate list (not \"int\") to list"
     ]
    }
   ],
   "source": [
    "a = [12,12]+ 2\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "olsres2 = smf.ols(formula = 'Y ~ X + I(X**2)-1', data = df).fit()\n",
    "print(olsres2.summary())\n",
    "xs = np.arange(1,max_L,1)\n",
    "plt.plot(xs,olsres2.params[1]*np.power(xs,2) + olsres2.params[0]*xs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "olsres2.params[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l = widgets.Text(value='                                 Monte Carlo Method to Estimate Mean Demonstration, Michael Pyrcz, Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "\n",
    "L = widgets.IntSlider(min=1, max = 40, value = 1, description = '$L$:',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "L.style.handle_color = 'gray'\n",
    "\n",
    "dist1 = widgets.Dropdown(\n",
    "    options=['Uniform','Triangular','Gaussian'],\n",
    "    value='Gaussian',\n",
    "    description='$X_1$:',\n",
    "    disabled=False,\n",
    "    layout=Layout(width='200px', height='30px')\n",
    ")\n",
    "min1 = widgets.FloatSlider(min=0.0, max = 100.0, value = 10.0, description = 'Min',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "min1.style.handle_color = 'gray'\n",
    "max1 = widgets.FloatSlider(min=0.0, max = 100.0, value = 90.0, description = 'Max',orientation='horizontal',layout=Layout(width='230px', height='50px'),continuous_update=False)\n",
    "max1.style.handle_color = 'gray'\n",
    "\n",
    "ui = widgets.HBox([L,dist1,min1,max1],kwargs = {'justify_content':'center'})\n",
    "ui1 = widgets.VBox([l,ui],)\n",
    "\n",
    "def make_sample_cdist(dist,zmin,zmax,L):\n",
    "    cdf = np.linspace(0.001,0.999,1000)\n",
    "    pvals = np.random.random(size = L)\n",
    "    if dist == 'Triangular':  \n",
    "        distribution = stats.triang(loc = zmin,c=0.5,scale = zmax-zmin)\n",
    "    if dist == 'Uniform':        \n",
    "        distribution = stats.uniform(loc = zmin,scale = zmax-zmin)\n",
    "    if dist == 'Gaussian':\n",
    "        mean = (zmax + zmin)*0.5; stdev = (zmax - zmin)/6.0\n",
    "        distribution = stats.norm(loc = mean,scale = stdev)\n",
    "    \n",
    "    cdfx = distribution.ppf(cdf)\n",
    "    sample = distribution.ppf(pvals)\n",
    "    \n",
    "    return sample, pvals, cdfx, cdf \n",
    "        \n",
    "def f_make1(L,dist1,min1,max1): \n",
    "    np.random.seed(seed = 73073)\n",
    "    sample, pvals, cdfx, cdf = make_sample_cdist(dist1,min1,max1,L)\n",
    "    truth_mean = np.average([min1,max1])\n",
    "    counts = np.arange(1,len(sample)+1,1)\n",
    "\n",
    "    L_all = 40\n",
    "    np.random.seed(seed = 73073)\n",
    "    sample_all, pvals_all, cdfx_all, cdf_all = make_sample_cdist(dist1,min1,max1,L_all)\n",
    "    prop_all = np.cumsum(sample_all, axis=0, dtype=None, out=None)/np.arange(1,len(sample_all)+1,1)\n",
    "    counts_all = np.arange(1,len(sample_all)+1,1)\n",
    "    \n",
    "    ax1 = plt.subplot(121)\n",
    "    plt.plot(cdfx,cdf,'--',color='black',linewidth = 3,zorder=10); add_grid(ax1)\n",
    "    plt.xlim(0,100); plt.ylim(0,1.0); plt.xlabel(\"$X_1$\"); plt.title(\"Cumulative Distribution Function, $X_1$\"); plt.ylabel('Cumulative Probability')\n",
    "     \n",
    "    for l in range(0,L):\n",
    "        alpha = max(0.02,(6 - (L-l) )/5)\n",
    "        dhead = 1;lhead = 0.02        \n",
    "        if l == L-1:\n",
    "            plt.plot([0.0,sample[l]],[pvals[l],pvals[l]],color='red',alpha = alpha,lw=1,zorder=1)\n",
    "            plt.plot([sample[l],sample[l]],[pvals[l],lhead+0.02],color='red',alpha = alpha,lw=1,zorder=1)\n",
    "            head = plt.Polygon([[sample[l],0.02],[sample[l]-dhead,lhead+0.02],[sample[l]+dhead,lhead+0.02],[sample[l],0.02]], color='red',alpha = alpha,zorder=1)\n",
    "            plt.annotate(int(sample[l]),[sample[l]+1,0.07],rotation=90.0,color='red',alpha = alpha)\n",
    "            plt.scatter(sample[l],0.01,s=30,color='red',edgecolor='black',alpha=1.0)\n",
    "            plt.gca().add_patch(head)\n",
    "        else:\n",
    "            plt.plot([0.0,sample[l]],[pvals[l],pvals[l]],color='darkorange',alpha = alpha,lw=1,zorder=1)\n",
    "            plt.plot([sample[l],sample[l]],[pvals[l],lhead+0.02],color='darkorange',alpha = alpha,lw=1,zorder=1)\n",
    "            head = plt.Polygon([[sample[l],0.02],[sample[l]-dhead,lhead+0.02],[sample[l]+dhead,lhead+0.02],[sample[l],0.02]], color='darkorange',alpha = alpha,zorder=1)\n",
    "            plt.annotate(int(sample[l]),[sample[l]+1,0.07],rotation=90.0,color='darkorange',alpha = alpha)\n",
    "            plt.scatter(sample[l],0.01,s=10,color='yellow',edgecolor='black',alpha=1.0)\n",
    "            plt.gca().add_patch(head)            \n",
    "\n",
    "    ax2 = plt.subplot(122)\n",
    "    plt.hist(sample,density= False,bins=np.linspace(0,100,20),weights=None,color='red',alpha=1.0,edgecolor='black',zorder=1)\n",
    "    if l > 0:\n",
    "        plt.hist(sample[:-1],density= False,bins=np.linspace(0,100,20),weights=None,color='darkorange',alpha=0.8,edgecolor='black',zorder=5)\n",
    "    \n",
    "    for l in range(0,L):\n",
    "        if l == L-1:\n",
    "            plt.scatter(sample[l],0.1,s=30,color='red',edgecolor='black',alpha=1.0,zorder=l+10)\n",
    "        else:\n",
    "            plt.scatter(sample[l],0.1,s=10,color='yellow',edgecolor='black',alpha=1.0,zorder=l+10)\n",
    "    \n",
    "    plt.vlines(truth_mean,0,10,colors='black',zorder=1000)\n",
    "    plt.annotate('Population Mean',[truth_mean+1,7],rotation=270.0)\n",
    "    \n",
    "    plt.plot(prop_all,counts_all/L_all*10,color='gray',lw=2,ls='--',zorder=999)\n",
    "    plt.plot(prop_all[:l+1],counts_all[:l+1]/L_all*10,color='blue',lw=2,zorder=1000)\n",
    "    plt.scatter(np.average(sample),L/L_all*10,c='blue',edgecolor='black',zorder=1000,label='Sample Mean')\n",
    "    plt.annotate('Sample Mean',[np.average(sample)+1,L/L_all*10+0.2],color='blue',zorder=1000)\n",
    "    \n",
    "    plt.xlabel(r\"$x_{\\alpha}, \\alpha = 1,\\ldots,L$\"); \n",
    "    plt.title(r\"Monte Carlo Arithmetic Average\"); plt.ylabel('Frequency')\n",
    "    plt.xlim([0,100]); plt.ylim([0,10]); add_grid(ax2)\n",
    "    plt.legend(loc='upper left')\n",
    "    \n",
    "    ax2 = plt.gca().twinx()\n",
    "    ax2.set_ylabel('Number of Realizations',color='blue'); ax2.set_ylim([1,L_all])\n",
    "    ax2.spines['right'].set_color('blue'); ax2.tick_params(axis='y', colors='blue')\n",
    "    \n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.3, hspace=0.2)\n",
    "    plt.show()    \n",
    "\n",
    "interactive_plot1 = widgets.interactive_output(f_make1, {'L':L,'dist1':dist1,'min1':min1,'max1':max1})\n",
    "interactive_plot1.clear_output(wait = True)                # reduce flickering by delaying plot updating  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(ui1, interactive_plot1)                            # display the interactive plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Monte Carlo Method to Calculate a Geometric Ratio of the Area of Circle Inscribed in a Square Demonstration\n",
    "\n",
    "Now let's set up our second dashboard. This is an example where we solve a geometric problem, the ratio of circle inscribed in a square by random sampling. \n",
    "\n",
    "* All we need to calculate is whether a point is in the circle, easy to do given the circle centroid and radius. \n",
    "\n",
    "#### Monte Carlo Method to Calculate Geometric Ratio of the Area of Circle Inscribed in a Square\n",
    "\n",
    "1. randomly sample a point within the square $v_1$\n",
    "    * sample independent uniform random values over x and y.\n",
    "2. check if the point is within the circle with a distance to center of the circle\n",
    "    * inside if less than the circle's radius\n",
    "3. calculate the proportion of points within the circule\n",
    "4. go to 1., continue until enough samples for the required accuracy\n",
    "\n",
    "From geometry we can calculate the exact solution:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\pi \\cdot r^2}{\\left( 2 \\cdot r \\right)^2} = \\frac{\\pi \\cdot r^2}{4 \\cdot r^2} = \\frac{\\pi}{4} \n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "title = widgets.Text(value='                     Monte Carlo Method to Estimate Ratio Circle Incribed in a Square Demonstration, Michael Pyrcz, Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "\n",
    "l = widgets.FloatSlider(min=0, max = 5, step = 0.1, value = 0.3, description = '$log(l)$:',orientation='horizontal',layout=Layout(width='750px', height='50px'),continuous_update=False)\n",
    "l.style.handle_color = 'gray'\n",
    "\n",
    "ui = widgets.HBox([l],kwargs = {'justify_content':'center'})\n",
    "ui2 = widgets.VBox([title,ui],)\n",
    "\n",
    "def f_make2(l): \n",
    "    L = 100000; radius = 50; cx = 50; cy = 50; seed = 73071\n",
    "    \n",
    "    i = int(10**l)\n",
    "    np.random.seed(seed=seed)\n",
    "    pts = np.random.rand(L,2)*100\n",
    "    pts_l = pts[:i]\n",
    "    dist = (pts[:,0]-50)**2 + (pts[:,1]-50)**2\n",
    "    dist_l = (pts_l[:,0]-50)**2 + (pts_l[:,1]-50)**2\n",
    "    inside = dist < (radius * radius)\n",
    "    inside_l = dist_l < (radius * radius)\n",
    "    \n",
    "    counts = np.arange(1,len(inside)+1,1)\n",
    "    prop = np.cumsum(inside, axis=0, dtype=None, out=None)/np.arange(1,len(inside)+1,1)\n",
    "    truth = math.pi/4\n",
    "    \n",
    "    ax1 = plt.subplot(121)\n",
    "    \n",
    "    plt.scatter(pts_l[inside_l,0],pts_l[inside_l,1],color='red',edgecolor='black',zorder=10)\n",
    "    plt.scatter(pts_l[inside_l==False,0],pts_l[inside_l==False,1],color='gray',edgecolor='black',zorder=10)\n",
    "    plt.scatter(50,50,marker='x',c='black',s=80,zorder=1)\n",
    "    \n",
    "    circle1 = plt.Circle((50, 50),radius,fill=False,edgecolor='black',zorder=1)\n",
    "    plt.plot([52,50+radius-1],[50,50],c='black',lw=1)\n",
    "    plt.gca().add_patch(circle1)\n",
    "    \n",
    "    plt.gca().add_patch(plt.Rectangle((0, 0), 100, 100, fill=False,edgecolor='black',lw=2))\n",
    "    \n",
    "    head = plt.Polygon([[cx+radius-1,cy-1],[cx+radius-1,cy+1],[cx+radius,cy],[cx+radius-1,cy-1]], color='black',alpha = 1.0,zorder=1)\n",
    "    plt.gca().add_patch(head)\n",
    "    plt.xlim([-10,110]); plt.ylim([-10,110]); add_grid(ax1); plt.xlabel(\"X (m)\"); plt.ylabel(\"Y (m)\") \n",
    "    plt.title('Ratio of Circle Inscribed Within Square')\n",
    "    \n",
    "    ax2 = plt.subplot(122)\n",
    "    plt.plot(counts,prop,c='gray',zorder=1)\n",
    "    plt.plot([0,L],[truth,truth],c='red',zorder=3); plt.annotate('Exact Solution',[2,truth+0.01],c='red')\n",
    "    plt.scatter(counts[:i-1],prop[:i-1],c='gray',edgecolor='black',s=30,zorder=9);\n",
    "    plt.scatter(counts[i-1],prop[i-1],c='red',edgecolor='black',s=30,zorder=10); add_grid(ax2)\n",
    "    plt.ylim([0.0,1.3]); plt.xlim([1,100000]); plt.xscale(\"log\"); plt.xlabel(\"Number of Samples\"); plt.ylabel(\"Cumulative Proportion\"); plt.title('Monte Carlo Geometric Ratio')\n",
    "    \n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.3, wspace=0.3, hspace=0.2); plt.show()\n",
    "    \n",
    "interactive_plot2 = widgets.interactive_output(f_make2, {'l':l})\n",
    "interactive_plot2.clear_output(wait = True)                # reduce flickering by delaying plot updating      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(ui2, interactive_plot2)                            # display the interactive plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Monte Carlo Method to Calculate Volume Integrated Variograms (Gamma Bar) Demonstration\n",
    "\n",
    "Now let's set up our third dashboard. This is an example of an expensive integration that we can approximate with a limited number of samples. \n",
    "\n",
    "\\begin{equation}\n",
    "\\overline{\\gamma}(v_1,v_2) = \\frac{1}{v_1 \\cdot v_2} \\int_{v_1} \\int_{v_2} \\gamma\\left( \\bf{u}_{v_1},\\bf{u}_{v_2} \\right) dv_1 dv_2\n",
    "\\end{equation}\n",
    "\n",
    "For a 2D problem this is a fourfold (quadruple) integral and for a 3D problem this is a sixfold (sextuple) integral, see below: \n",
    "\n",
    "In 2D:\n",
    "\\begin{equation}\n",
    "\\overline{\\gamma}(v_1,v_2) = \\frac{1}{v_1 \\cdot v_2} \\int_{{v_1}_X} \\int_{{v_1}_Y} \\int_{{v_2}_X} \\int_{{v_2}_Y} \\gamma\\left( \\bf{u}_{v_1},\\bf{u}_{v_2} \\right) dv_{1_X} dv_{1_Y} dv_{2_X} dv_{2_Y}\n",
    "\\end{equation}\n",
    "\n",
    "and in 3D:\n",
    "\n",
    "\\begin{equation}\n",
    "\\overline{\\gamma}(v_1,v_2) = \\frac{1}{v_1 \\cdot v_2} \\int_{{v_1}_X} \\int_{{v_1}_Y} \\int_{{v_1}_Z} \\int_{{v_2}_X} \\int_{{v_2}_Y} \\int_{{v_2}_Z} \\gamma\\left( \\bf{u}_{v_1},\\bf{u}_{v_2} \\right) dv_{1_X} dv_{1_Y} dv_{1_Z} dv_{2_X} dv_{2_Y} dv_{2_Z} \n",
    "\\end{equation}\n",
    "\n",
    "given the semivariogram calculation computational cost, combinatorial of variogram values to average to calculate 1 gamma bar value, and combined with the large number of gamm bar values often required in geostatistical spatial estimation and simulation, this is a too high computational complexity and we need a speed up!\n",
    "\n",
    "#### Monte Carlo Method to Calculate Gamma Bar Values\n",
    "\n",
    "1. randomly sample a lag vector between volumes $v_1$ and $v_2$\n",
    "    * uniform random value over x and y for each the lag vector tail (in $v_1$) and the lag vector head (in $v_2$).\n",
    "2. calculate the variogram value over the lag vector\n",
    "    * in this example, I assume a isotropic spherical variogram model\n",
    "3. calculate the average variogram over all random lag vectors\n",
    "4. go to 1., continue until enough samples for the required accuracy\n",
    "\n",
    "Note, the 'Exact Solution' for the gamma bar was calculated with 100k samples and hard coded in the dashbaord below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "title = widgets.Text(value='                Monte Carlo Method to Estimate Volume Integrated Variograms, Gamma Bars, Demonstration, Michael Pyrcz, Professor, The University of Texas at Austin',layout=Layout(width='950px', height='30px'))\n",
    "\n",
    "L_dia = widgets.FloatSlider(min=0.25, max = 2.0, step = 0.25, value = 0.3, description = '$log(l)$:',orientation='horizontal',layout=Layout(width='750px', height='50px'),continuous_update=False)\n",
    "L_dia.style.handle_color = 'gray'\n",
    "\n",
    "ui = widgets.HBox([L_dia],kwargs = {'justify_content':'center'})\n",
    "ui3 = widgets.VBox([title,ui],)\n",
    "\n",
    "def f_make3(L_dia): \n",
    "    L = int(10**L_dia)\n",
    "    seed = 73073; r = 200.0; L_all = 100\n",
    "    \n",
    "    gamma_truth = 0.81038       # calculated with 100k samples\n",
    "    np.random.seed(seed=seed)\n",
    "    pts_v1_all = np.random.rand(L_all,2)\n",
    "    pts_v1_all[:,0] = pts_v1_all[:,0]*50.0+25.0; pts_v1_all[:,1] = pts_v1_all[:,1]*50.0+35.0; \n",
    "    pts_v2_all = np.random.rand(L_all,2);\n",
    "    pts_v2_all[:,0] = pts_v2_all[:,0]*100.0+125; pts_v2_all[:,1] = pts_v2_all[:,1]*100.0+10.0; \n",
    "    gamma_all = isotropic_spherical_variogram(pts_v1_all[:,0],pts_v1_all[:,1],pts_v2_all[:,0],pts_v2_all[:,1],r)\n",
    "    \n",
    "    pts_v1 = pts_v1_all[:L,:]; pts_v2 = pts_v2_all[:L,:]; gamma = gamma_all[:L]\n",
    "    \n",
    "    counts_all = np.arange(1,L_all+1,1)\n",
    "    prop_all = np.cumsum(gamma_all, axis=0, dtype=None, out=None)/np.arange(1,L_all+1,1)\n",
    "    \n",
    "    f, (ax1, ax2) = plt.subplots(1, 2, gridspec_kw={'width_ratios': [3, 1]})\n",
    "    \n",
    "    for l in range(0,L):\n",
    "        alpha = max(0.02,(11 - (L-l) )/10)        \n",
    "        if l == L-1:\n",
    "            rotate = rotation(pts_v2[l,1]-pts_v1[l,1],pts_v2[l,0]-pts_v1[l,0])\n",
    "            ax1.scatter(pts_v1[l,0],pts_v1[l,1],color='red',edgecolor='black',alpha = alpha,s=20,zorder=150)\n",
    "            ax1.scatter(pts_v2[l,0],pts_v2[l,1],color='red',edgecolor='black',alpha = alpha,s=20,zorder=150)\n",
    "            ax1.plot([pts_v1[l,0],pts_v2[l,0]],[pts_v1[l,1],pts_v2[l,1]],color='red',alpha = alpha,zorder=130)\n",
    "            ax1.annotate(r'$\\gamma(v_1^{\\ell},v_2^{\\ell})$ = ' + str(np.round(gamma[l],2)),\n",
    "                         [(pts_v1[l,0]+pts_v2[l,0])*0.5,(pts_v1[l,1]+pts_v2[l,1])*0.5+2],\n",
    "                         rotation = rotate,color='red',ha='center',zorder=100)\n",
    "            ax1.annotate(r'$\\overline{\\gamma}(v_1,v_2)$ = ' + str(np.round(prop_all[l],2)),\n",
    "                         [5,5],color='red',zorder=100)\n",
    "        else:\n",
    "            ax1.scatter(pts_v1[l,0],pts_v1[l,1],color='gray',edgecolor='black',alpha = alpha,s=20,zorder=50)\n",
    "            ax1.scatter(pts_v2[l,0],pts_v2[l,1],color='gray',edgecolor='black',alpha = alpha,s=20,zorder=50)\n",
    "            ax1.plot([pts_v1[l,0],pts_v2[l,0]],[pts_v1[l,1],pts_v2[l,1]],color='black',alpha = alpha,zorder=30)\n",
    "    \n",
    "    ax1.add_patch(plt.Rectangle((25,35), 50, 50, fill=False,edgecolor='black',lw=2))\n",
    "    ax1.annotate(r'$v_1$',[26,86])\n",
    "    ax1.add_patch(plt.Rectangle((125,10), 100, 100, fill=False,edgecolor='black',lw=2))\n",
    "    ax1.annotate(r'$v_2$',[126,111])\n",
    "    add_grid(ax1)\n",
    "    ax1.set_xlim([0,240]); ax1.set_ylim([0,120]); ax1.set_xlabel('X (m)'); ax1.set_ylabel('Y (m)'); ax1.set_title('Volume Integrated Variograms (Gammabar)')\n",
    "    \n",
    "    counts = np.arange(1,L+1,1)\n",
    "    prop = np.cumsum(gamma, axis=0, dtype=None, out=None)/np.arange(1,L+1,1)\n",
    "    \n",
    "    ax2.plot(prop_all,counts_all,color='gray',zorder=10); ax2.scatter(prop,counts,color='gray',edgecolor='black',s=20,zorder=100)\n",
    "    ax2.scatter(prop[l],counts[l],color='red',edgecolor='black',s=50,zorder=120)\n",
    "    ax2.plot([gamma_truth,gamma_truth],[1,L_all],color='red',zorder=20); \n",
    "    plt.annotate('Exact Solution',[gamma_truth,3],c='red',rotation=270)\n",
    "    ax2.set_xlim([0.5,1.0]); ax2.set_ylim([1,L_all]); ax2.set_xlabel(r'$\\overline{\\gamma}(v_1,v_2)$'); ax2.set_ylabel('Number of Samples'); ax2.set_title('Monte Carlo Gammabar')\n",
    "    add_grid(ax2); ax2.set_yscale('log')\n",
    "    \n",
    "    plt.subplots_adjust(left=0.0, bottom=0.0, right=2.0, top=1.2, wspace=0.3, hspace=0.2); plt.show() \n",
    "    \n",
    "interactive_plot3 = widgets.interactive_output(f_make3, {'L_dia':L_dia})\n",
    "interactive_plot3.clear_output(wait = True)                # reduce flickering by delaying plot updating      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "display(ui3, interactive_plot3)                            # display the interactive plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Comments\n",
    "\n",
    "This was an interactive demonstration of Monte Carlo methods to solve statistical and mathematical problems. I have other lectures and inteactive dashboards on Monte Carlo Simulation and I wanted to expand our coverage with general Monte Carlo methods here. I also have lectures on Markov chain Monte Carlo (McMC) methods!\n",
    "  \n",
    "I hope this was helpful,\n",
    "\n",
    "*Michael*\n",
    "\n",
    "#### The Author:\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. Professor, Cockrell School of Engineering and 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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}