{
"cells": [
{
"cell_type": "markdown",
"id": "8cdf9df1",
"metadata": {},
"source": [
"# 2-D prior toy example (RBF covariance)"
]
},
{
"cell_type": "markdown",
"id": "b92e9179-4463-4fda-95fb-191e97796811",
"metadata": {},
"source": [
"> Demonstration of prior error bounds for a 2-D toy example"
]
},
{
"cell_type": "markdown",
"id": "6b817ae4",
"metadata": {},
"source": [
"In this notebook we will demonstrate the error bounds for the statFEM prior for the toy example introduced in [twoDim](01_twoDim.ipynb#2-dimensional-case-(PDE)). We first import some required packages."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "354f84b8",
"metadata": {},
"outputs": [],
"source": [
"from dolfin import *\n",
"set_log_level(LogLevel.ERROR)\n",
"import numpy as np\n",
"import numba\n",
"import matplotlib.pyplot as plt\n",
"plt.rcParams['figure.figsize'] = (10,6)\n",
"# import required functions from twoDim\n",
"from statFEM_analysis.twoDim import mean_assembler, kernMat, cov_assembler\n",
"from scipy.linalg import sqrtm\n",
"import sympy; sympy.init_printing()\n",
"from tqdm import tqdm\n",
"# code for displaying matrices nicely\n",
"def display_matrix(m):\n",
" display(sympy.Matrix(m))"
]
},
{
"cell_type": "markdown",
"id": "e2531b90",
"metadata": {},
"source": [
"Since we do not have ready access to a solution operator for this problem we will utilise a different approach to demonstrate the appropriate rate for the error bounds. This approach involves comparing the statFEM distributions on finer and finer meshes as follows. We first assume we have the following error bound:\n",
"\n",
"$$W(\\nu_{\\star},\\nu_{h})\\leq\\mathcal{O}(h^{p}) \\text{ as } h\\rightarrow 0$$\n",
"\n",
"From this we have:\n",
"\n",
"$$W(\\nu_{\\star},\\nu_{h})=Ch^{p}+\\mathcal{O}(h^{p+1})$$ \n",
"\n",
"The triangle inequality then yields:\n",
"\n",
"$$W(\\nu_{h},\\nu_{h/2})\\leq W(\\nu_{h},\\nu_{\\star}) + W(\\nu_{\\star},\\nu_{h/2})$$\n",
"\n",
"And so with the above we obtain:\n",
"\n",
"$$W(\\nu_{h},\\nu_{h/2})=Ch^{p}+C(h/2)^{p}+\\mathcal{O}(h^{p+1})$$\n",
"\n",
"Similarly we have:\n",
"\n",
"$$W(\\nu_{h/2},\\nu_{h/4})=C(h/2)^{p}+C(h/4)^{p}+\\mathcal{O}(h^{p+1})$$\n",
"\n",
"Dividing the two above equations yields:\n",
"\n",
"$$\\frac{W(\\nu_{h},\\nu_{h/2})}{W(\\nu_{h/2},\\nu_{h/4})}=\\frac{Ch^{p}+C(h/2)^{p}+\\mathcal{O}(h^{p+1})}{C(h/2)^{p}+C(h/4)^{p}+\\mathcal{O}(h^{p+1})}=\\frac{1-2^{-p}+\\mathcal{O}(h)}{2^{-p}-2^{-2p}+\\mathcal{O}(h)}=2^{p}+\\mathcal{O}(h)$$\n",
"\n",
"From which it follows:\n",
"\n",
"$$\\operatorname{LR}(h):=\\log_{2}\\frac{W(\\nu_{h},\\nu_{h/2})}{W(\\nu_{h/2},\\nu_{h/4})}=p+\\mathcal{O}(h)$$\n",
"\n",
"Note that this ratio above does not rely on having access to the true prior $\\nu_{\\star}$.\n",
"\n",
"We will thus compare the errors on a succesively refined mesh and take the base-2 log of the above ratio. We hope to see these logarithms converge to $p=2$."
]
},
{
"cell_type": "markdown",
"id": "242e557f",
"metadata": {},
"source": [
"We now set up the mean and kernel functions for the random forcing term $f$."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "6ddcbd15",
"metadata": {},
"outputs": [],
"source": [
"# set up mean and kernel functions\n",
"f_bar = Constant(1.0)\n",
"\n",
"l_f = 0.4\n",
"σ_f = 0.1\n",
"\n",
"@numba.jit\n",
"def k_f(x):\n",
" return (σ_f**2)*np.exp(-(x**2)/(2*(l_f**2)))"
]
},
{
"cell_type": "markdown",
"id": "66640463",
"metadata": {},
"source": [
"We now set up a function to get the statFEM prior for a FE mesh size $h$, using functions from twoDim."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "d450240b",
"metadata": {},
"outputs": [],
"source": [
"def fem_prior(h,f_bar,k_f,grid):\n",
" J = int(np.round(1/h))\n",
" μ = mean_assembler(h,f_bar)\n",
" Σ = cov_assembler(J,k_f,grid,False,True)\n",
" return μ,Σ"
]
},
{
"cell_type": "markdown",
"id": "a38cac40",
"metadata": {},
"source": [
"We now set up a reference grid on which we will compare the covariance matrices."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "6cfe93ea",
"metadata": {},
"outputs": [],
"source": [
"N = 41\n",
"x_range = np.linspace(0,1,N)\n",
"grid = np.array([[x,y] for x in x_range for y in x_range])"
]
},
{
"cell_type": "markdown",
"id": "1c83e55a",
"metadata": {},
"source": [
"We now set up a function to compare the covariance functions on the reference grid."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "44442a43",
"metadata": {},
"outputs": [],
"source": [
"def compute_cov_diff(C1,C2,tol=1e-10):\n",
" N = np.sqrt(C1.shape[0])\n",
" C1_sqrt = np.real(sqrtm(C1))\n",
" rel_error_1 = np.linalg.norm(C1_sqrt @ C1_sqrt - C1)/np.linalg.norm(C1)\n",
" assert rel_error_1 < tol\n",
" \n",
" C12 = C1_sqrt @ C2 @ C1_sqrt\n",
" C12_sqrt = np.real(sqrtm(C12))\n",
" rel_error_12 = np.linalg.norm(C12_sqrt @ C12_sqrt - C12)/np.linalg.norm(C12)\n",
" assert rel_error_12 < tol\n",
" \n",
" hSq = (1/(N-1))**2\n",
" return hSq*(np.trace(C1) + np.trace(C2) - 2*np.trace(C12_sqrt))"
]
},
{
"cell_type": "markdown",
"id": "cc481005",
"metadata": {},
"source": [
"We now set up a function to compute the Wasserstein distance between the statFEM priors."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "3a7317e1",
"metadata": {},
"outputs": [],
"source": [
"def W(μ_1,μ_2,Σ_1,Σ_2):\n",
" mean_error = errornorm(μ_1,μ_2,'L2')\n",
" cov_error = compute_cov_diff(Σ_1,Σ_2)\n",
" cov_error = np.sqrt(np.abs(cov_error))\n",
" error = mean_error + cov_error\n",
" return error"
]
},
{
"cell_type": "markdown",
"id": "29670c5d",
"metadata": {},
"source": [
"In the interests of memory efficiency we will now create a function which will compute the ratios of errors mentioned above by succesively refining the FE mesh."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "d4f43c4f",
"metadata": {},
"outputs": [],
"source": [
"def refine(h,n,f_bar,k_f,grid):\n",
" # set up empty lists to hold h-values and errors (this being the ratios)\n",
" h_range = []\n",
" errors = []\n",
" # get the statFEM prior for h and h/2\n",
" μ_1, Σ_1 = fem_prior(h,f_bar,k_f,grid)\n",
" μ_2, Σ_2 = fem_prior(h/2,f_bar,k_f,grid)\n",
" # compute the distance between these and store in numerator variable\n",
" numerator = W(μ_1,μ_2,Σ_1,Σ_2)\n",
" # succesively refine the mesh by halving and do this n times\n",
" for i in tqdm(range(n), desc=\"inner loop\", position=1, leave=False):\n",
" # store mean and cov for h/2 in storage for h\n",
" μ_1, Σ_1 = μ_2, Σ_2 \n",
" # in storage for h/2 store mean and cov for h/4\n",
" μ_2, Σ_2 = fem_prior(h/4,f_bar,k_f,grid)\n",
" # compute the distance between the priors for h/2 and h/4\n",
" # and store in denominator variable\n",
" denominator = W(μ_1,μ_2,Σ_1,Σ_2)\n",
" # compute the ratio and store in error\n",
" error = numerator/denominator\n",
" # append the current value of h and the ratio\n",
" h_range.append(h)\n",
" errors.append(error)\n",
" # store denominator in numerator and halve h\n",
" numerator = denominator\n",
" h = h/2\n",
" # return the list of h-values together with the ratios for these values\n",
" return h_range,errors"
]
},
{
"cell_type": "markdown",
"id": "a27d149b",
"metadata": {},
"source": [
"
\n",
"\n",
"Note:\n",
"\n",
"The above function will compute the ratio $\\frac{W(\\nu_{h},\\nu_{h/2})}{W(\\nu_{h/2},\\nu_{h/4})}$ for the values $h,h/2,\\cdots,h/2^{n-1}$.\n",
" \n",
"
"
]
},
{
"cell_type": "markdown",
"id": "61be823b",
"metadata": {},
"source": [
"We now set up a list of starting $h$ values and number of refinements $n$ to get a decent number of ratios to approximate $p$."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "0d6b6a9c",
"metadata": {},
"outputs": [],
"source": [
"my_list = [(0.25,4),(0.2,3),(0.175,3),(0.22,4),(0.176,3),(0.177,3),(0.178,3),(0.179,3),(0.18,3),(0.21,3),(0.215,3),(0.1,2),(0.3,4),(0.31,4),(0.315,4),(0.24,4),(0.245,4),(0.25,4),(0.14,2),(0.16,2),(0.15,2),(0.18,2)]"
]
},
{
"cell_type": "markdown",
"id": "4929a34c",
"metadata": {},
"source": [
"
\n",
"\n",
"Warning: \n",
"\n",
"As we are now in 2-D the matrices needed to compute the statFEM prior get very large in size as $h$ decreases. Due to memory constraints we must thus make sure $h$ does not get too small.\n",
"\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "49f4d257",
"metadata": {},
"source": [
"We now compute the results:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "576e550e-fe49-46a5-abb8-cfa291c1a39d",
"metadata": {},
"outputs": [],
"source": [
"set_log_level(LogLevel.ERROR)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "bcd8874d-7277-45be-b39d-c95f083146b6",
"metadata": {},
"outputs": [],
"source": [
"n_outer=len(my_list)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "e4311ad1",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"outer: 0it [00:00, ?it/s]\n",
"inner loop: 0%| | 0/4 [00:00"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(hs,np.log2(es),c='r')\n",
"plt.xlabel('$h$')\n",
"plt.ylabel('$\\operatorname{LR}(h)$')\n",
"plt.title('Base-2 log plot of the ratios against $h$')\n",
"plt.hlines(2.0,hs[0],hs[-1],colors='black',linestyles='--')\n",
"plt.xlim(hs[-1],hs[0])\n",
"plt.grid()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "e9eca82a",
"metadata": {},
"source": [
"We can see from the above plot that the logarithms seems to be approaching $p=2$ as $h$ gets small just as predicted. However, the results aren't that smooth and they haven't seemed to settle on $p=2$ yet. This can be due to memory constraints meaning we cannot use small enough $h$. We thus smooth the above results by taking a cumulative average of the ratios and then applying $\\log_2$. We take the rolling average starting with large $h$. We choose a cut-off point of $h=0.15$, i.e. we discard any results for $h>0.15$.\n",
"\n",
"The results are shown below:"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "92e582a7",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmoAAAGHCAYAAAAA4H6+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdd5xU9fX/8deBpQqCiK4NxN4QG1Y0Lmps0VhjjVGjIbEk9liSqImar/5MjBorUURjwRZjRWNbAcECCgKxBHCXKkVAWKQInN8fn7sy7M7uzuzOnTu7834+HvcxM/dz595z7swOh8+993PN3RERERGRwtMq6QBEREREJD0VaiIiIiIFSoWaiIiISIFSoSYiIiJSoFSoiYiIiBQoFWoiIiIiBUqFmoiIiEiBUqEmIiIiUqBUqEnBMbMKMzskge1ONLOyetoTias+ZradmY01s8Vm9psM31NweTRFU/Jp6DOXpmkJ3884viON2S8ZrHOame2ei3VJYVGh1sJFP3pLzazKzBaY2ctm1iOP229nZg+aWWX0ozTWzI7I1/az4e47uXs5FN4/FvX4LfC2u3d29ztrNjajPPKi5v5I/cwlFs3q+5kunpi+I/Xul2yZ2XrAJsCnTY5s7fV2M7PnzGxJ9Bt+WgPLX2hmo81suZkNrtFWbmbLon+Lqszs81zG2pKpUCsOR7t7J2BjYDbw9zxuuwSYBhwIdAF+DzxlZr3yGEO9zKwk6RiaYHNgYtJBSP41k+9twXw/C2x/NWq/1JPDzkCFuy9tUlS13Q2sAEqB04F7zWynepafCdwIDKqj/UJ37xRN2+U21BbM3TW14AmoAA5JeX0k8EXK66uAycBi4L/AcSltVwIzorbPgYOj+ZsAzwJzgS+B32QZ0yfACZnEDOwAlAMLCT9sP05Zbnfg4yi+p4EngRszzK0iyu8TYDmhoKwADgH+CawGlgJVwG9T3nN59J5vou21r7HOK6L2JcCDhB+4oVEMbwDr1bdv0+yL+vJ/C1gFLIvi3LbGe7POI5vPFugB/Cta9mvgrpQ2B7ZOeT24xmeTzb7KZF3V35f6PvNa+6PGe68EnqmR4x3AnY3YN2njaOo2SP+9rS/nhv5G6tweTfyOEsP3sxGfQ8b7q4F4Gvw9ynSfpdsv9a03XQ5p1nkh8BpwZ7RfZgI/zOC3uCvwE0Jh9XSNtnUIRdq2KfP+CdycwXpvBAbXmFcOnNvQezWl2Z9JB6Ap5g947R+ZjsDDwCMp7T+JfvhaAScT/tHcGNiO0BO2SbRcL2CraLkxwLVAW2BLYApwWIbxlEY/UNs3FDPQBpgEXBNt66DoB3C76HUlcFG03PHRj8qNDeWWso2xhGKjQ5p99f3zGnF9EK2zG+Eww69qtL8X5bgpMAf4CNgNaE/4gb6urn2bZj/UmX/KMvX++GWTRzafLdAaGAf8jfCD3h7YP6U9k+KqwX2VxbqqP7c6P/N0+6PGezcHvgU6p+Q4C9gnm33TwN9Vk7ZB+u9tXduq92+kvu2Ro+8oOfx+NhRzPevPaH81EE+9v0fRchnts5r7JYP11sohzfruJ/xn6bAor+uAN+pYdhfgamAYoRB+CTgf2LzGcrsB39aYdznwYga/83UVanOBecC7QFlD69EUJh36LA7/NrOFhD/KHwK3Vje4+9PuPtPdV7v7k8D/gL0I/+NrB+xoZm3cvcLdJwN7Ahu4+5/cfYW7TwH+AZzSUBBm1gZ4DHjY3T/LIO59gE6E/8GtcPe3CD8qp0ZtJYSeiO/c/V+EH/fv1ZNbtTvdfZpnd7jgzmid84EXgV1rtP/d3We7+wxgOPC+u3/s7suA5wg/fnXt22zyb6p0eWTz2e5F+MfuCndf4u7L3H1EljFksq+yksFnXt97KwnF4nHRrIMI/1C9R5bf+7riyNE21vre1pNzQ38j9W0v6e9oXX9njfn9yXR/NaShXDPdZ9mut1YOafQBbnL319x9NaGnsBYze4FQ9B0L3AyUuvtR7n5P9N1M1QlYVGPeN0DnDHJK50pCYb0pMBB40cy2auS6iooKteJwrLt3JfRUXAi8Y2YbAZjZz6IT/BdGxVxvoLu7TwIuBq4H5pjZEDPbhNAjsEn18tF7riH0jGBmp6ecLDq0OgAza0XoNl8RxZCJTYBp0Q9PtUrCH/omwAz38F+1yLTUN9eVW13LZ+irlOffEn7MUs1Oeb40zetO9ezbmurLv6nS5VHvZ1tDD6DS3Vc2IYYG91W2K8zgM2/I46z5B/K06DVkt28aiqOp28j0e97Q30id2yuA72hdf2dZfQ6RbH8X6lJvrlnss6zWmy6HGvlYlMOLKbN7k75Y+yNwE6EX72FgkJmdYWYbplm2Cli3xrx1Cb19WXP39919sbsvd/eHCb1qRzZmXcVGhVoRcfdV0f+qVwH7m9nmhP+NXgisHxVzEwCLln/c3fcn/Dg6cAvhB+NLd++aMnV29yOj9zzma04WPQK+/yGpPgfpBHf/LsOQZwI9oiKvWk/COSCzgE2jdVf7/mrWhnKr3iX17a4MY2yUOvZtTfXln/Gmsli23s82zbI96zm5+VvCofZqG2URR6PWlYPPHMJ5XGVmthmh16u6iMp432QQR1O38X0ODWyr3r+RhraXp+9otn9n2XxHa20jg8+mvngazDXDfZb1ehuIa4vocVLKvN0IPWdrcfcx7v57d9+dcAHC68DRwKdm9kiNxb8ASsxsm5R5u5C7i0Octf82pQ4q1IqIBccA6xHO+1iH8McyN2o/m/A/sepxfg4ys3aEc8qWEk60/QBYbGZXmlkHM2ttZr3NbM96Nn0v4WTZo+vpuk/nfcI/0r81szYWxjI6GhgCjCIUnBeaWUmUV+rhizpzy9BsQjd9ztWzb2uqL/9MZZNHNp/tB4RC4GYzW8fM2ptZv5T2scBp0ToOJ1z121iZriuTz7ze/eHucwnn0jxEKAiqhzvIZt/UG0eOtpHJthr6G6lze3n8jmb7d9aYfZSqoe9IffHUm2sW+yyr9WagDzC+Rs/pboRzSOvk7l+5+0PufhKwIaGnLbV9CeFioT9Ff+P9gGMIR0bSir5n7QnnXraOfhdKzKyrmR2W8vp04AfAqxnmWNRUqBWHF82sinC+wU3Ame4+0d3/C/yV8IM+m/A/rHej97QjnMMwj3AYYkPgandfBRxFOGfky6j9AcLQG7VE/4P9ZbT8V7bmsOjpDQXt7isIP1hHRNu5B/iZu38WtR0PnEO4UuqnhPM6lkfvrS+3TPwf8Pvo8MjlWbwvE2n3bc2F6ss/i21lnEc2n2207NHA1sBUYDrhxOxqF0XtCwmX9f87i5hrymhdGX7mmeyPxwknj1f3dGW7bzKJo0nbyGRbGfyN1Le9fH1Hs/o7a8w+qvH+hj6bOuPJINeM9lmamJq6D3cmpSgzs+6EXucJNRe0MJaZ15yAlYSjHjWdD3QgXOzzBHCeu09MWd9QM7smZfnfEwrUqwjft6XRvDaECwyqLyb4NeGUnC8yzLGo2dpFuEjzZWbvA/e5+0NJxyJSiPQ3ItL8qEdNmi0zO9DMNoq60s8kHAJQV7pIRH8jIs1fIY3ULJKt7YCnCOedTAFOdPdZyYYkUlD0NyLSzOnQp4iIiEiB0qFPERERkQKlQk1ERESkQLXYc9S6d+/uvXr1inUbS5YsYZ111ol1G4WqmHOH4s6/mHOH4s5fuRdn7qD882HMmDHz3H2DmvPzUqiZWQ/gEcLI9A4MdPc7aixjwB2EW0p8C5zl7h9FbauA8dGiU939xw1ts1evXowePTp3SaRRXl5OWVlZrNsoVMWcOxR3/sWcOxR3/sq9LOkwElPs+eeDmdW83yqQvx61lcBl7v6RmXUGxpjZ69Hgg9WOALaJpr0Jo9nvHbUtdfeaN78WERERadHyco6au8+q7h1z98WE2xfVvGnvMcAjHrwHdDWzjfMRn4iIiEghyvvwHGbWCxgG9Hb3RSnzXwJudvcR0es3gSvdfbSZrSTc729ltEzaW8iY2QBgAEBpaekeQ4Zkc0vE7FVVVdGpU6dYt1Goijl3KO78izl3KO78lXtx5g7KPx/69+8/xt371pyf14sJzKwT8CxwcWqRloHN3X2GmW0JvGVm4919cs2F3H0gMBCgb9++Hvfx9GI+Zl/MuUNx51/MuUNx56/cy5IOIzHFnn+S8jY8h5m1IRRpj7n7v9IsMgPokfJ6s2ge7l79OAUoB3aLNVgRERGRApCXQi26ovNB4FN3v62OxV4AfmbBPsA37j7LzNYzs3bReroD/YD/1rEOERERkRYjX4c++wFnAOPNbGw07xqgJ4C73we8QhiaYxJheI6zo+V2AO43s9WEwvLmGleLioiIiLRIeSnUogsErIFlHLggzfyRwM4xhSYiIiJSsHQLKREREZECpUJNREREpECpUBMREREpUC32puwiIiICrF4NCxbA3LlhWrkSNtwwTN26QevWSUco9VChJiIi0hx9+y2MHg1ffbWmCEs3ff01rFqVfh2tWkH37lBauqZ423DDWq9bf/ttfnOT76lQExERaQ7c4fPPYehQePVVeOcdWL587WW6dYMNNgjTtttCv35rXldPJSWhgJs9G+bMWXv64IPwuHjxWqvdv1Ur2GMPOPDAMO2/P3Ttmsfki5cKNRERkUK1eDG89VYozF59FSoqwvwddoDzz4dDDoGePUMBtv76oQjLhaVLQzE3Zw7MmkXl00/Tq6IC7rwT/vIXMINdd11TuB1wQNi+5JwKNRERkULhDhMmhKJs6FAYMQK++w46dYKDD4arroLDDoNeveKNo0OHUAD27AlARefO9CorCwXc+++H3rx33oH77oPbbw/v2XnnNYXbD34QDptKk6lQExERScrChTBuXJg++gjeeANmzAhtO+8Ml1wChx8eDmG2bZtsrBAKuLKyMEE49Prhh2sKt0GD4K67Qlu/fnDiiWHabLOkIm72VKiJiIjEzT0cthw7NhRl1Y/VhzIhHL78wQ/giCNCcbbppklFm7l27cL5avvvD7/7Xej9GzMGXn8dnnkmFJqXXAL77rumaIt66SQzKtRERKQ4rF4No0bBv/8derKysO2sWfD006Ewads2TA09nzdv7aJs0aKwslatwon+e+8Nv/wl7LJLON9ro43CuV/NWZs2sM8+YfrDH+CLL0LB9vTTcNllYdp77zVFW9yHcFsAFWoiItJyuYdC6Ykn4MknYerUUEhleeL7+suXw3vvwYoVYap5tWVdOnWCPn3gpz9dU5D17g0dOzYimWZo223hmmvCNGlSKNqeeQauuCJMe+65pmjbcsukoy1IKtRERKTl+fzzUJwNGRKel5TAoYfCTTfBMcdA585ZrW5UeTll1edlQSgAV65cU7SlFnDVzzt3DsVHK90ECICttw4XQ1x1FUyZsqZou/LKMO20Uyjcdt89TLvsEgrdIqdCTUREWoapU0Ov2RNPwMcfh8OIBx4Il14Kxx8fBnbNFbNwmK9NG1hnndytt1hsuSX89rdhqqgIBdtbb8Err8DgwWEZM9h++zWF2+67w267QZcuSUaedyrURESk+Zo7F556KhRn774b5u21F/ztb/CTnzSPE/KLXa9ecPnlYXKHWbPCFbBjxoTHd96Bxx5bs/xWW4XBd3ffPVx9uueeLbrXUoWaiIg0L6tXh96Xf/wDnnsuXGnYu3c4rHnyyeEfcmmezGCTTcJ01FFr5s+ZE4q26unDD0OBDmHZY48NU1lZ6OVsQVSoiYhI8zBrFjz0EDz4YDjHqVs3uOAC+PnPw5hj0nJtuGEYsuTww9fMmzcvDAr83HPhe3HPPeG2VkcdBccdFwYGbgGHpVWoiYhI4Vq1Cl57LfSevfhieF1WBjfeGP4xbt8+6QglKd27wxlnhOnbb8PYbc89F74njz4aBuc99NDwPTn66FDYN0Mq1EREpPBMmxZGuX/wwfB8ww3DGFznngvbbJN0dFJoOnYMV/Mec0y4GnfYsFC0/fvf8Pzz0Lp1uLDkuOPgpJOa1e2tVKiJiEjy3MOFASNHht6zV18N8374Q7jtNvjxjwvjFkpS+EpK4KCDwnTnnTB6dCjYnnsOfv3rMLbboYcmHWXGVKiJiEh+rFgRhtCYPDmcYzZ58prnU6ZAVVVYbuON4eqr4ZxzYIstko1ZmjezcFXonnuGi00+/7zZfadUqImISO7NnQuPPw4TJ64pxqZODVdsVmvfPoynteWW0L9/eNxxx3AOWon+eZIYbLdd0hFkTX8JIiK5Nncu/Otf4f6GkyeHwyzHHBMOxbT0k99Hj4a//z3cEWDFinCj8a22gv32Cyd9b7lleL3VVuHeli14/CuRXFChJiKSC3PnhnNgnnoK3n479Bxtu224z+Pjj8PAgWGogMMOC0Xbj36U9f0mC9by5WFk+b//Hd5/P9z25xe/CENn7LBD0tGJNGsq1EREGmvevLWLs1WrwhWJV18dRsXv0yecI7N8eWh//nl44YXQ29aqFey//5or1ZrjIK0zZsD994cidPbsUJjeeSeceSasu27S0Ym0CCrURESy0Oabb8JViU8/HUbHX7Uq3Gz6yivDZf/VxVmqdu3WDNZ5993h1jgvvBAKt8suC9OOO64p2nbaac06UtdV17xWrfI3Grs7jBgBd90VCs5Vq8IAoxdeCIccokOZIjmmQk1EisvKlfDyy/DAA2F8rlat1kytW9f/eskS9hs5MhzWrC7OfvIT2GWX2sVZXVq1WnMV2g03hJPsX3ghTP/v/8H//V/j8tp443CYcfvtw1T9fNNNM4+tLu4wfz4bv/wyXHwxjBsH660Xnp93XjjvTERikZdCzcx6AI8ApYADA939jhrLGHAHcCTwLXCWu38UtZ0J/D5a9EZ3fzgfcYtICzJzZijO/vEPmD493B9wjz1CEbJ6de1p1apQ1KXOa9WKqaecwuZXXJFdcVafLbcMBc/FF8P8+WH8sBkzQpv72o91Pf/uu1DwffZZGJF90aI1y3TqtHbhVv18k03C9ubMWTPNnbv269R5K1eyHYQew3/8A047LQwyKiKxyleP2krgMnf/yMw6A2PM7HV3/2/KMkcA20TT3sC9wN5m1g24DuhLKPLGmNkL7r4gT7GLSHO1ejW8+Sbcd184zLhqVbgC8847wy1lGjEExJfl5Wy+664xBEu4xc1ppzVtHe7w1VehaPv00zWPb78N//xnw+9fZ50wavuGG0LPntC3b3i+wQZ8XFLCbhdemJsCVUQykpdCzd1nAbOi54vN7FNgUyC1UDsGeMTdHXjPzLqa2cZAGfC6u88HMLPXgcOBJ/IRe30uvvhiunbtuta8k046ifPPP59vv/2WI488stZ7zjrrLM466yzmzZvHiSeeWKv9vPPO4+STT2batGmcccYZtdovu+wyjj76aD7//HN++ctf1mr//e9/zyGHHMLYsWO5+OKLa7X/+c9/Zr/99mPkyJFcc801tdpvv/12dt11V9544w1uvPHGWu33338/2223HSNHjuT666+v1f7Pf/6THj168OSTT3LvvffWan/mmWfo3r07gwcPZvDgwbXaX3nlFTp27Mg999zDU089Vau9vLwcgL/85S+89NJLa7V16NCBoUOHAnDDDTfw5ptvrtW+/vrr8+yzzwJw9dVXM2rUqLXaN9tsMx599FEgfLZjx45dq33bbbdl4MCB32+/Zv677rort99+OwA//elPmT59+lrt++67L/8XHdY64YQT+Prrr9dqP/jgg/nDH/4AwBFHHMHSpUvXaj/qqKO4/PLLASgrK6OmvHz3DjyQuWPHUpbms8/Xd+/FF1/kr3/9a6321O/eY3feyeFffcWPZ81is6VL+aakhDbnnUfHiy5i8IgRDL7jDrhjrU79jL97Tz75ZK3PPp/fvQEDBvDFF1+s1b7Wd++MM9J/915/HRYv5rfHHEOXWbPotmIF37Rpw8I2bdi6Xz/Ovfpq2GADjjjhhDXfvUWLYNEijtp+ey6/5BK+KS+nrH//WvumGH73gIy+ey35d6/B715L/d07+uha8/Mp7+eomVkvYDfg/RpNmwLTUl5Pj+bVNT/dugcAAwBKS0u//3LHZdWqVSxcuHCteV988QXl5eUsW7asVhvAZ599Rnl5Od98803a9okTJ1JeXs6cOXPSto8fP57OnTszderUtO3jxo2jpKSESZMmpW3/6KOPWLFiBRMmTEjbPnr0aBYuXMi4cePStr///vvMmjWrzvxGjRrF5MmTmThxYtr2d999ly5duvDZZ5+lbR82bBjt27fniy++SNte/ZlOnjy5VvvSpUu/b//yyy9rta9evfr79nT7r02bNt+3T58+vVb7zJkzv29fuXJlrfbp06d/3z579uxa7VOnTv2+fe7cuSxKPTwVxVzdPn/+fJYvX75W++TJk79vT7dv4vzutV+1ik533cV3p53GiUuW0Gq99bhro41YkHICe76+e+PHj0//3Rs5krlz57LdwIE8PWEC7dz5uGNH7unZkze6dOF3hxxCl+nTm/zdW758eaLfvZkzZzbpu/feypUsatcuXOAQ2X3FCrb+8kv48st6v3tVVVV5/+5BYfzudevWre7vXhH87lVVVTX5u9fcfvdgzXcvUe6etwnoBIwBjk/T9hKwf8rrNwmHOy8Hfp8y/w/A5Q1ta4899vC4vf3227Fvo1AVc+7uRZT/8uXud93lXlrqDu5HHeXTjjvOvaTEfd113W+9NSyTdIwPPujep0+IsXNn9/PPd//kk1g2VzSffRrKvXgVe/75AIz2NPVM3q6jNrM2wLPAY+7+rzSLzAB6pLzeLJpX13wRicuqVfDII+F2KxdeGB5HjIAXX2TSb34D48dDv35wxRVhKIkXX1z7JPd8WLIEbr89jD92zjlh3v33h4sG7r4bdt45v/GIiMQgL4VadEXng8Cn7n5bHYu9APzMgn2Abzyc2/YacKiZrWdm6wGHRvNEJNfcwwCuffqEQUu7dQtXIZaXh8Ks2vbbwyuvhKmkBH784zDi/sSJ8ce4YEEY1mLzzeGSS8JVk0OHwtixMGBAuMpRRKSFyFePWj/gDOAgMxsbTUea2a/M7FfRMq8AU4BJwD+A8wE8XERwA/BhNP0pmiciufTGG7D33nD88aFH7emnw30bDzus7qv8jjgCPvkk9Gx9+GEYsuLCC6HGCcM5MWtW6MHr2ROuvRb22Sf08r3zThhIVlciikgLlK+rPkcA9f6KRsdnL6ijbRAwKIbQROS99+B3vwuj7PfsCYMGhZtnZzp0RZs2cNFFcPrpcN11cO+98NhjcP31cP75TR8xf8qUMBDsQw+Fcc1OPhmuuir0+omItHC614dIsZoyBY49FvbdN5xzdvvt8MUXcPbZjRpfjO7dw7lh48aFsbcuvjgUU9HQAVkbPz4Uf9tsE4q0s86Czz8PNzhXkSYiRUK3kBIpNitXhqLs2mtDQXbDDaGoytW5Xb17w3/+Ay+9BJdeCkceGQZRrR4Sol07aNt27dc158+fHw7FrrNOOA/t0kvDSPoiIkVGhZpIMRk7Nlwh+dFHYWT+e+6BzTbL/XbMwvoPOywcSv3f/2DFCli+PP20YgUsXrzmtRn88Y/hfLdu3XIfn4hIM6FCTaQYLF0aCp+//AXWXx+eegpOPDH+E/DbtoVf/arh5UREJC0VaiIt3dtvh2ErJk2Cn/8cbr1VvVQiIs2ELiYQaakWLIBf/AIOOiiMj/bGG/DggyrSRESaERVqIi2NOzzzDOywQ7ha8re/DWOdHXxw0pGJiEiWdOhTpCWZMQMuuACefx522y3cOWD33ZOOSkREGkk9aiItgTs88ADsuCO89loYIPaDD1SkiYg0c+pRE2nuli+H884Lhzn794eBA2HrrZOOSkREckA9aiLN2cyZcOCBoUi79tpwwYCKNBGRFkM9aiLN1fvvw3HHwaJF8Oyz4WbqIiLSoqhHTaQ5eugh+MEPoH17GDVKRZqISAulQk2kOfnuO7joojBw7QEHwIcfws47Jx2ViIjERIWaSHMxb164d+add4Yblb/6argdlIiItFg6R02kOfjkEzjmGJg1Cx5+GH72s6QjEhGRPFCPmuTH88/Do48mHUXz9PTTsO++sGIFDBumIk1EpIioR03it2gRnHkmfPMNdO8Ohx+edETNw+rVYciNm24Khdqzz8LGGycdlYiI5JF61CR+998firTNN4ef/hSmTUs6osL3zTfhUOdNN8E558Dbb6tIExEpQirUJF7Ll8Pf/gYHHQT/+U94fcop4epFSW/OnDCI7auvwt13wz/+Ae3aJR2ViIgkQIWaxOvRR8MJ8FdeCdtuG4qOkSPhd79LOrLCNGNGKNK++AJefhnOPx/Mko5KREQSokJN4rN6Ndx6K+y2G/zwh2HeKaeE+1Leeiu89FKy8RWaioowiO2MGeHG6ocemnREIiKSMBVqEp/nn4fPP4ff/nbtXqHbbgvF289+BpWVycVXSL74Igxgu2BBuF/nAQckHZGIiBQAFWoSD3e45RbYYgs48cS129q3D0NOrFoFJ50Uhp0oZuPHh560FSugvBz22ivpiEREpECoUJN4DBsWbhp++eVQkmYUmK22gkGD4IMPwvlrxWr0aCgrC/vonXegT5+kIxIRkQKiQk3iccstsMEGcPbZdS9zwgnw61/D7bfDc8/lL7ZCMWJEuBq2SxcYPhy23z7piEREpMCoUJPc++QTGDoUfvMb6NCh/mVvvRX23DMUdFOm5Ce+QvDGG+G+nZtsEnoft9gi6YhERKQA5aVQM7NBZjbHzCbU0b6emT1nZp+Y2Qdm1julrcLMxpvZWDMbnY94pYn+3/+DddYJQ0s0pF07eOqpcLHBSSeFcdZauhdfhKOOgq23Doc7N9ss6YhERKRA5atHbTBQ332DrgHGunsf4GfAHTXa+7v7ru7eN6b4JFcqKmDIEBgwALp1y+w9vXrB4MEwZgxcdlmMwRWAp56C448P56K9/TaUliYdkYiIFLC8FGruPgyYX88iOwJvRct+BvQyM/0L1hzddlvoHbvkkuzed8wxcOmlYST+p56KJ7akDR4Mp54a7tv5xhuZF7IiIlK0CuWm7OOA44HhZrYXsDmwGTAbcOA/Zk3zTJQAACAASURBVObA/e4+MLkwpV7z5sEDD8Dpp0OPHtm//+abYdQoOPfcMM7aNtvkPsZqCxeGc+K+/DKMXda6dbjysqQk/fMa8zp/+il07rx2e33TkCHhwolDDw0XTnTsGF9uIiLSYpi752dDZr2Al9y9d5q2dQmHO3cDxgPbA79w97Fmtqm7zzCzDYHXgV9HPXTptjEAGABQWlq6x5AhQ2LJpVpVVRWdOnWKdRuFKl3uvQYPptfDD/PBQw/xba9ejVpvuzlz6PuLX7B8gw346O67Wd3Ie1zaypW0mzOHDrNm0X7mzLUfZ82izeLFjVpvU8zr14+J116Lt22b923nUjF/76G481fuxZk7KP986N+//5h0p3gVRKFWYzkDvgT6uPuiGm3XA1Xu/peGtte3b18fPTreaw/Ky8spKyuLdRuFqlbuS5ZAz57Qrx+88ELTVv7KK/CjH4Xz3O6/Pwyeu2QJzJ8fer8WLEj/fP58mDs39JJNmxYG1K3Wpk04F27LLWtP668fll21ClauDFO65ynzPvn4Y/rstNOa9zU0de4Mp50W4mjmivl7D8Wdv3IvSzqMxBR7/vlgZmkLtYI49GlmXYFv3X0FcC4wzN0Xmdk6QCt3Xxw9PxT4U5KxSh0efDAUSrkYvPbII8N6brklHCZcsCAUSHVp3RrWWy9M3bvDfvutXYhtsQVsumlYLkfmt28fBqoVERGJUV4KNTN7AigDupvZdOA6oA2Au98H7AA8HJ2HNhE4J3prKfBc6GSjBHjc3V/NR8yShe++g7/+NfSm9euXm3XeeCO0bQtz5oQCrFu3NcVY6vP11gu9Van3EhUREWkh8lKoufupDbSPArZNM38KsEtccUmOPPkkTJ0Kd92Vu3WWlMCf1HkqIiLFTXcmkKZxDwPc7rhjOK9MREREcqYgzlGTZmzoUBg/PowR1kp1v4iISC7pX1ZpmltuCWOmnVrv0W0RERFpBPWoSeO99164ofjf/hZO/BcREZGcUo+aNN4tt4SrLs89N+lIREREWiQVatIoHadOheefhwsvBI1WLSIiEgsVatIoPYYMgfbtw/0rRUREJBYq1CQ77vDhh5S+/jr8/OewwQZJRyQiItJi6WICadi8efDGG/Daa/Cf/8DMmXj79nDZZUlHJiIi0qKpUJPavvsuXNH52mthGjMm9KSttx4ccggcdhjvd+nCfltskXSkIiIiLZoKNQm+/HJNYfbWW7BoURjAdp994Prr4bDDoG/f729svqK8PNFwRUREioEKtWK3bBkcfDCMHBle9+wJJ58cCrODD4auXZONT0REpIipUCt2o0aFIu2KK+Ccc2DbbcEs6ahEREQEFWoyYkQozK65Rr1nIiIiBUbDcxS74cNh551VpImIiBQgFWrFbOXKcOjzgAOSjkRERETSUKFWzMaNg6oq2H//pCMRERGRNFSoFbPhw8OjCjUREZGCpEKtmI0YAb16wWabJR2JiIiIpKFCrVi5hx41nZ8mIiJSsFSoFatJk2DOHBVqIiIiBUyFWrHS+WkiIiIFT4VasRo+HNZfH7bfPulIREREpA4q1IrViBGhN023ixIRESlYKtSK0VdfhXPUdH6aiIhIQVOhVoxGjAiPOj9NRESkoKlQK0bDh0PHjrD77klHIiIiIvVQoVaMhg+HffaBNm2SjkRERETqkZdCzcwGmdkcM5tQR/t6ZvacmX1iZh+YWe+UtsPN7HMzm2RmV+Uj3hZt0aJwj08d9hQRESl4+epRGwwcXk/7NcBYd+8D/Ay4A8DMWgN3A0cAOwKnmtmO8Ybawo0aBatX60ICERGRZiAvhZq7DwPm17PIjsBb0bKfAb3MrBTYC5jk7lPcfQUwBDgm7nhbtBEjoHXrcOhTRERECpq5e342ZNYLeMnde6dp+zPQwd0vMbO9gJHA3sAWwOHufm603BnA3u5+YR3bGAAMACgtLd1jyJAhcaTyvaqqKjp16hTrNnJt14svptWyZXx0331NWk9zzD2Xijn/Ys4dijt/5V6cuYPyz4f+/fuPcfe+NeeXJBFMGjcDd5jZWGA88DGwKtuVuPtAYCBA3759vaysLJcx1lJeXk7c28ipFSvg88/hV79qctzNLvccK+b8izl3KO78lXtZ0mEkptjzT1JBFGruvgg4G8DMDPgSmAJ0AHqkLLoZMCPvAbYUY8bAsmU6P01ERKSZKIjhOcysq5m1jV6eCwyLircPgW3MbIuo/RTghaTibPZ0I3YREZFmJS89amb2BFAGdDez6cB1QBsAd78P2AF42MwcmAicE7WtNLMLgdeA1sAgd5+Yj5hbpBEjYNttYcMNk45EREREMpCXQs3dT22gfRSwbR1trwCvxBFXUVm9OhRqxx+fdCQiIiKSoYI49Cl58OmnsGCBDnuKiIg0IyrUikX1+Wm6kEBERKTZUKFWLIYPh402gi23TDoSERERyZAKtWIxYkToTTNLOhIRERHJkAq1YjB1aph02FNERKRZUaFWDEaMCI+6kEBERKRZybpQM7N1zKx1HMFITIYPh86doU+fpCMRERGRLDRYqJlZKzM7zcxeNrM5wGfALDP7r5ndamZbxx+mNMnw4bDfftBa9bWIiEhzkkmP2tvAVsDVwEbu3sPdNwT2B94DbjGzn8YYozTF/PkwcaLOTxMREWmGMrkzwSHu/l3Nme4+H3gWeNbM2uQ8MsmNd98NjyrUREREmp0GC7XqIs3MDgJOBxYCE4BPgAnuvjxdIScFYsQIaNMG9twz6UhEREQkS9nc63MQcDHhZup9gGOBnQCdo1bIhg8PRVqHDklHIiIiIlnKplCrdPd/R8+fjiMYybGlS2H0aLjkkqQjERERkUbI5KrPR8zsYuA9M7s0DzFJrnzwAXz3nc5PExERaaYyuepzMGBAKXCGmVWa2QtmdoOZ/STW6KRpqm/Evt9+ycYhIiIijZLJxQRvAW9VvzazEmAHYBdgL3QYtHCNGAG9e0O3bklHIiIiIo3QYKFmZubuXv3a3VcC46Pp0XTLSAFYtQpGjoSfaog7ERGR5iqjAW/N7Ndm1jN1ppm1NbODzOxh4Mx4wpNGGzcOFi/W/T1FRESasUyu+jwc+DnwhJltQRhHrT3QGvgPcLu7fxxfiNIo1Tdi14UEIiIizVYm56gtA+4B7onuQNAdWOruC+MOTppg+HDo2RN69Eg6EhEREWmkTA59fs/dv3P3We6+MLpZ++lxBSZN4B561NSbJiIi0qxlMo7aumZ2tZndZWaHWvBrYApwUvwhStYmT4avvtL5aSIiIs1cJueo/RNYAIwCzgWuIYyrdqy7j40xNmksnZ8mIiLSImRSqG3p7jsDmNkDwCygZ3TumhSi4cPD2Gk77JB0JCIiItIEmZyj9l31E3dfBUxXkVbghg+Hfv2gVVanIIqIiEiByaRHbRczWxQ9N6BD9NoAd/d1Y4tOsjd7Nvzvf/CLXyQdiYiIiDRRJsNztK6rzcwuyW040mTV56fpQgIREZFmr6nHxjIq1MxskJnNMbMJdbR3MbMXzWycmU00s7NT2laZ2dhoeqGJ8bZ8I0ZA+/awxx5JRyIiIiJN1NRCzTJcbjDhDgd1uQD4r7vvApQBfzWztlHbUnffNZp+3OhIi8HMmfDEE+H8tLZtG15eREREClpTC7WMbsTu7sOA+Q2sp7OZGdApWnZlE2MrLsuXwwknQFUV3H570tGIiIhIDjR4jpqZLSZ9QWZAhxzFcRfwAjAT6Ayc7O6ro7b2ZjaaULjd7O7/ztE2Ww53uOACeO89eOYZ6N076YhEREQkB8w9o06xpm/IrBfwkrvXqiLM7ESgH3ApsBXwOrCLuy8ys03dfYaZbQm8BRzs7pPr2MYAYABAaWnpHkOGDIkll2pVVVV06tQp1m1kYpPnn2fb22+n8vTT+fLcc/OyzULJPSnFnH8x5w7Fnb9yL87cQfnnQ//+/ce4e99aDe6elwnoBUyoo+1l4ICU128Be6VZbjBwYibb22OPPTxub7/9duzbaNDw4e4lJe5HHOG+cmXeNlsQuSeomPMv5tzdizt/5V68ij3/fABGe5p6plBGRJ0KHAxgZqXAdsAUM1vPzNpF87sTet3+m1iUhWb6dDjxROjVCx5/HFrXOZKKiIiINEOZDHjbZGb2BOFqzu5mNh24DmgD4O73ATcAg81sPOHctyvdfZ6Z7Qfcb2arCRc+3OzuKtQAli0LFw8sWQJvvglduyYdkYiIiORYXgo1dz+1gfaZwKFp5o8Edo4rrmar+uKBDz6Af/0Ldtop6YhEREQkBplc9Xlpfe3uflvuwpGM3HsvDBoEv/89HHdc0tGIiIhITDLpUescPW4H7EkYRgPgaOCDOIKSegwfDhddBD/6Efzxj0lHIyIiIjHK5F6ffwQws2HA7u6+OHp9PeFqTcmX6osHttgCHn0UWhXKtSAiIiISh2zOUSsFVqS8XhHNk3xYtgyOPx6WLoXycl08ICIiUgSyKdQeAT4ws+ei18cCD+c+JKnFHc4/Hz78EP79b9hhh6QjEhERkTzIuFBz95vMbChwQDTrbHf/OJ6wZC333AMPPQTXXgvHHJN0NCIiIpInGZ/kFN0wfUegi7vfAXxtZnvFFpkEw4bBxRfD0UfDddclHY2IiIjkUTZno98D7AtUj4m2GLg75xHJGtUXD2y1Ffzzn7p4QEREpMhkc47a3u6+u5l9DODuC8ysbUxxCcDf/w4LF4ZetS5dko5GRERE8iybLprvzKw14ABmtgGwOpaoJHjlFTjgANh++6QjERERkQRkU6jdCTwHlJrZTcAI4M+xRCUwbRpMmABHHJF0JCIiIpKQbK76fMzMxgAHR7OOdfdP4wlLePXV8HjkkcnGISIiIonJuFAzs3bA7kCX6H0/MTPc/U9xBVfUXnkFevbUmGkiIiJFLJuLCZ4HvgHGAMvjCUcAWLEC3ngDTj8dzJKORkRERBKSTaG2mbsfHlskssa770JVlc5PExERKXLZXEww0sx2ji0SWWPoUGjTBg46KOlIREREJEEN9qiZ2XjCkBwlwNlmNoVw6NMAd/c+8YZYhIYODcNydO6cdCQiIiKSoEwOfR4VexSyRvWwHLfemnQkIiIikrAGD326e6W7VwLnVz9PnRd/iEWmelgOnZ8mIiJS9LI5R+2Haeapmsi16mE5dtwx6UhEREQkYZmco3YeoedsSzP7JKWpM/BuXIEVJQ3LISIiIikyOUftcWAo8H/AVSnzF7v7/FiiKlYalkNERERSNFioufs3hIFuTzWzXYADoqbhgAq1XNKwHCIiIpIi43PUzOw3wGPAhtH0qJn9Oq7AipKG5RAREZEU2dyZ4Fxgb3dfAmBmtwCjgL/HEVjR0bAcIiIiUkM2V30asCrl9aponuTC0KHhUeeniYiISCSbHrWHgPfN7Lno9bHAg7kPqUgNHQo9emhYDhEREflexj1q7n4bcDbhAoL5wNnufnum7zezQWY2x8wm1NHexcxeNLNxZjbRzM5OaTvTzP4XTWdmus1mo3pYjiOP1LAcIiIi8r1setRw94+Ajxq5rcHAXcAjdbRfAPzX3Y82sw2Az83sMaATcB3Ql3DP0TFm9oK7L2hkHIVHw3KIiIhIGlkVajWH53D3cZm+192HmVmv+hYBOpuZEYqz+cBK4DDg9eox28zsdeBw4IlsYi9oGpZDRERE0shmeI6LiHd4jruAHYCZwHjgIndfDWwKTEtZbno0r+XQsBwiIiKShrl7ZguG20ftmzI8xzrAKHfvk/HGQo/aS+7eO03biUA/4FJgK+B1YBdgANDe3W+MlvsDsNTd/5JmHQOi5SktLd1jyJAhmYbWKFVVVXTq1KlJ62g3Zw77nnwyk3/1K6adfHKOIotfLnJvzoo5/2LOHYo7f+VenLmD8s+H/v37j3H3vjXnZ3PoM+7hOc4GbvZQOU4ysy+B7YEZQFnKcpsB5elW4O4DgYEAffv29bKysnSL5Ux5eTlN3sbAgQBsdeGFbLXTTk0PKk9yknszVsz5F3PuUNz5K/eypMNITLHnn6TGDs9hhOE5BuUwlqnAwcBwMysFtgOmAJOAP5vZetFyhwJX53C7ydKwHCIiIlKHjAs1d7/NzMoJhycBznT3sZm+38yeIPSMdTez6YQrOdtE674PuAEYbGbjCYXgle4+L3rvDcCH0ar+1GJuBl89LMfpp2tYDhEREakl40LNzPoCvwN6Re8bYGae6Tlq7n5qA+0zCb1l6doGkdveu8KgYTlERESkHtkc+nwMuIJwRebqeMIpMhqWQ0REROqRTaE2191fiC2SYqRhOURERKQe2RRq15nZA8CbwPLqme7+r5xHVQymTYMJE+DWW5OORERERApUNoXa2YThMtqw5tCnAyrUGmPo0PCo89NERESkDtkUanu6+3axRVJsNCyHiIiINCDjW0gBI81MVUUuVA/LccQRGpZDRERE6pRNj9o+wNjojgHLCWOdZTw8h6SoHpbjyCOTjkREREQKWDaF2uGxRVFsNCyHiIiIZCCbQ597AfPdvRI4A/gb0C2WqFo6DcshIiIiGcimUPuDuy82s/2BQ4AHgfviCasFqx6WQ1d7ioiISAOyKdRWRY8/Aga6+8tA29yH1MJpWA4RERHJUDaF2gwzux84GXjFzNpl+X4BDcshIiIiGcum0DoJeA04zN0XEs5PuyKWqFoqDcshIiIiWcj4qk93/5aUuxC4+yxgVhxBtVjVw3LosKeIiIhkQIcu86l6WI6DD046EhEREWkGVKjlk4blEBERkSyoUMuXr74Kw3IcrnGDRUREJDMq1PLlf/8Lj310xy0RERHJjAq1fKmoCI+9eiUZhYiIiDQjKtTypbIyPPbsmWwcIiIi0myoUMuXigooLYUOHZKORERERJoJFWr5UlkJm2+edBQiIiLSjKhQyxcVaiIiIpIlFWr5sHp1KNR0IYGIiIhkQYVaPsyeHe7zqR41ERERyYIKtXzQ0BwiIiLSCCrU8qF6aA71qImIiEgWVKjlQ3WPmgo1ERERyUJJPjZiZoOAo4A57t47TfsVwOkpMe0AbODu882sAlgMrAJWunvffMScU5WV0K2bbsYuIiIiWclXj9pgoM67kbv7re6+q7vvClwNvOPu81MW6R+1N78iDTQ0h4iIiDRKXgo1dx8GzG9wweBU4IkYw8m/igpdSCAiIiJZM3fPz4bMegEvpTv0mbJMR2A6sHV1j5qZfQksABy4390H1vP+AcAAgNLS0j2GDBmSs/jTqaqqolOnTvUv5M4BRx7JzKOOYvIFF8QaTz5llHsLVsz5F3PuUNz5K/fizB2Ufz70799/TLojh3k5Ry0LRwPv1jjsub+7zzCzDYHXzeyzqIeulqiIGwjQt29fLysrizXY8vJyGtzG3LmwbBk99t+fHjHHk08Z5d6CFXP+xZw7FHf+yr0s6TASU+z5J6nQrvo8hRqHPd19RvQ4B3gO2CuBuBpPQ3OIiIhIIxVMoWZmXYADgedT5q1jZp2rnwOHAhOSibCRNNitiIiINFK+hud4AigDupvZdOA6oA2Au98XLXYc8B93X5Ly1lLgOTOrjvVxd381HzHnjHrUREREpJHyUqi5+6kZLDOYMIxH6rwpwC7xRJUnlZWw7rrQtWvSkYiIiEgzUzCHPlusiorQmxZ6BUVEREQypkItbhrsVkRERBpJhVrcNNitiIiINJIKtTgtXAiLFqlHTURERBpFhVqcNDSHiIiINIEKtThpaA4RERFpAhVqcaou1NSjJiIiIo2gQi1OFRXQoQN07550JCIiItIMqVCLU2Vl6E3TGGoiIiLSCCrU4lQ92K2IiIhII6hQi5MGuxUREZEmUKEWl6oq+PprXUggIiIijaZCLS4amkNERESaSIVaXDTYrYiIiDSRCrW4qEdNREREmkiFWlwqK6FtW9hoo6QjERERkWZKhVpcKiqgZ09opV0sIiIijaMqIi4amkNERESaSIVaXCoqdCGBiIiINIkKtTgsXQqzZ6tHTURERJpEhVocpk4Nj+pRExERkSZQoRYHDc0hIiIiOaBCLQ7VhZp61ERERKQJVKjFoaICWreGTTZJOhIRERFpxlSoxaGyEnr0gJKSpCMRERGRZkyFWhwqKnR+moiIiDSZCrU4aLBbERERyYG8FGpmNsjM5pjZhDrarzCzsdE0wcxWmVm3qO1wM/vczCaZ2VX5iLdJVqyAGTN0IYGIiIg0Wb561AYDh9fV6O63uvuu7r4rcDXwjrvPN7PWwN3AEcCOwKlmtmM+Am606dPBXT1qIiIi0mR5KdTcfRgwP8PFTwWeiJ7vBUxy9ynuvgIYAhwTQ4i5o6E5REREJEcK6hw1M+tI6Hl7Npq1KTAtZZHp0bzCVVERHtWjJiIiIk1UaONHHA286+6Z9r6txcwGAAMASktLKS8vz2FotVVVVdXaRq933mFzM4ZNmYJPm5b+jS1AutyLSTHnX8y5Q3Hnr9zLkw4jMcWef5IKrVA7hTWHPQFmAD1SXm8WzUvL3QcCAwH69u3rZWVlMYS4Rnl5ObW2MXgwbLIJB/7wh7FuO2lpcy8ixZx/MecOxZ2/ci9LOozEFHv+SSqYQ59m1gU4EHg+ZfaHwDZmtoWZtSUUci8kEV/GKit1fpqIiIjkRF561MzsCaAM6G5m04HrgDYA7n5ftNhxwH/cfUn1+9x9pZldCLwGtAYGufvEfMTcaBUVsN9+SUchIiIiLUBeCjV3PzWDZQYThvGoOf8V4JXcRxWDVavC8By6kEBERERyoGAOfbYIM2fCypU69CkiIiI5oUItlzQ0h4iIiOSQCrVc0mC3IiIikkMq1HKpuketZ89EwxAREZGWQYVaLlVWQmkpdOiQdCQiIiLSAqhQy6WKCp2fJiIiIjmjQi2XKitVqImIiEjOqFDLldWrYepUXUggIiIiOaNCLVdmz4bly9WjJiIiIjmjQi1XNDSHiIiI5JgKtVzRYLciIiKSYyrUcqW6R02FmoiIiOSICrVcqaiAbt2gc+ekIxEREZEWQoVarlRW6vw0ERERySkVarmiMdREREQkx1So5YK77kogIiIiOadCLRe+/hq+/VaHPkVERCSnVKjlgobmEBERkRioUMsFDXYrIiIiMVChlgvqURMREZEYqFDLhcpKWHdd6No16UhERESkBVGhlgvVV3yaJR2JiIiItCAq1HJBg92KiIhIDFSo5YIGuxUREZEYqFBrqoUL4ZtvVKiJiIhIzqlQayoNzSEiIiIxUaHWVBqaQ0RERGKiQq2p1KMmIiIiMclLoWZmg8xsjplNqGeZMjMba2YTzeydlPkVZjY+ahudj3izUlEBHTpA9+5JRyIiIiItTEmetjMYuAt4JF2jmXUF7gEOd/epZrZhjUX6u/u8eENspOqhOTSGmoiIiORYXnrU3H0YML+eRU4D/uXuU6Pl5+QjrpzQ0BwiIiISE3P3/GzIrBfwkrv3TtN2O9AG2AnoDNzh7o9EbV8CCwAH7nf3gfVsYwAwAKC0tHSPIUOG5DiLtVVVVXHY6acz98AD+eLSS2PdVqGpqqqiU6dOSYeRmGLOv5hzh+LOX7kXZ+6g/POhf//+Y9y9b835+Tr02ZASYA/gYKADMMrM3nP3L4D93X1GdDj0dTP7LOqhqyUq4gYC9O3b18vKymINevjQobRZtIhN9tuPTWLeVqEpLy8n7v1byIo5/2LOHYo7f+VelnQYiSn2/JNUKFd9Tgdec/cl0blow4BdANx9RvQ4B3gO2CuxKGto99VX4YkOfYqIiEgMCqVQex7Y38xKzKwjsDfwqZmtY2adAcxsHeBQoM4rR/Ot/ezZ4YmG5hAREZEY5OXQp5k9AZQB3c1sOnAd4Zw03P0+d//UzF4FPgFWAw+4+wQz2xJ4zsIVlSXA4+7+aj5izkR79aiJiIhIjPJSqLn7qRkscytwa415U4gOgRai9rNnQ9u2sNFGSYciIiIiLVChHPpsltrPng09e0Ir7UYRERHJPVUYTdD+q690fpqIiIjERoVaE7SbPVvnp4mIiEhsVKg11rJltJs/X4WaiIiIxEaFWmNNnRoedehTREREYqJCrbEqKsKjetREREQkJirUGquyMjyqR01ERERiokKtsSor8VatYJNNko5EREREWigVao1VUcGyDTeEkkK5r72IiIi0NCrUGquykuWlpUlHISIiIi2YuoMaq3dvFixbRtek4xAREZEWS4VaY917L5Xl5WyRdBwiIiLSYunQp4iIiEiBUqEmIiIiUqBUqImIiIgUKBVqIiIiIgVKhZqIiIhIgVKhJiIiIlKgVKiJiIiIFCgVaiIiIiIFSoWaiIiISIFSoSYiIiJSoFSoiYiIiBQoFWoiIiIiBUqFmoiIiEiBMndPOoZYmNlcoDLmzXQH5sW8jUJVzLlDcedfzLlDceev3ItXseefD5u7+wY1Z7bYQi0fzGy0u/dNOo4kFHPuUNz5F3PuUNz5K/fizB2Uf5J06FNERESkQKlQExERESlQKtSaZmDSASSomHOH4s6/mHOH4s5fuRevYs8/MTpHTURERKRAqUdNREREpECpUEvDzA43s8/NbJKZXZWm/Qdm9pGZrTSzE2u0rTKzsdH0Qv6izp0m5v+qmS00s5fyF3H2MsixnZk9GbW/b2a9ovltzOxhMxtvZp+a2dUp77nIzCaY2UQzuzhl/vVmNiPle3FkPnJsjDj2S6FpbI5RWx8zGxV9xuPNrL2ZdTSzl83ss2j+zSnLn2Vmc1M++3Pzk2V6Tfh8e5nZ0pQ87ovmF3zuTfm8o/aeZlZlZpfXmN/azD5O/a0zs8eibU0ws0Fm1iaaX2Zm36Tsi2vjybZuTfjs90qJe5yZHRfN3y5l/lgzW1T9u9ecfvOaBXfXlDIBrYHJwJZAW2AcsGONZXoBfYBHgBNrtFUlnUPC+R8MHA28lHQuTczxfOC+6PkpwJPR89OAIdHzjkBFtD96pV13DAAAByJJREFUAxOieSXAG8DW0XLXA5cnnXcS+yXpnHKcYwnwCbBL9Hr9aH0dgf7RvLbAcOCI6PVZwF1J552D3HsBE9Kss6Bzb0rOKe3PAE/X/BsGLgUeJ+W3DjgSsGh6Ajgvml9Ggr+JTfzsOwIl0fONgTnVr2us/yvCOGDQTH7zmsukHrXa9gImufsUd18BDAGOSV3A3Svc/RNgdRIBxqxJ+bv7m8DivETaeA3mGL1+OHr+DHCwmRngwDpmVgJ0AFYAi4AdgPfd/Vt3Xwm8Axwffyo5Fcd+KTRNyfFQ4BN3Hwfg7l+7+6roM387mrcC+AjYLA+5ZKspuafVDHJvUs5mdizwJTAx9Q1mthnwI+CB1Pnu/opHgA8onH3R6P2Q8psG0J7wt17TwcBkd497kPmipEKttk2BaSmvp0fzMtXezEab2XvRH3lz09T8m4NMcvx+mehH6htCD8ozwBJgFjAV+Iu7zyf0ph1gZuubWUfC/6x7pKzvQjP7JDocsl4MOeVCHPvl/7d3b6FWVHEcx78/0jKULNLIMK8ZhuWpTC1TgxQkKEnswVv0UA89CAmVoUFZRFEQSVfIgqDrQw9RYGbYRTGK0LwdFS0LSUIrirBCyvPvYa1t4/EcL3vvc85s9+8Dmz2zZs2c9Z/Ze+a/Z82cKZtaYrwUCEkfKXX9L26/cEnnks4orykUz87b/l1JF7efpxvVEjvA8NzV97mkKe0XXtLYq45ZUj/gAeCRDpa7HFhMJz/Wc5fn7cCqQvF1uevwQ0ljqoilFjVte0kTJbUCW4G7C4lbxRzSGcSiRtjnNQQnavU3NNJ/b54HLJc0sqcbZHU1ATgMXAQMB+6VNCIidgBPAqtJO+dNuR7AS8BI4EpSIvN0dze6G3S4Xnq2SXXXC5gMzM/vsyRNq0zMZxPfBp6NiD25+ANSF/BY4GP+P2PRaH4ChkTEVeQuP0nnVCaeprEvA56JiIPFQkk3AwciYsNx5n0RWBsR6/L4RtKxoQV4DnivC9rbZSLiq4gYA4wHlkjqU5km6UxgJql7uKIZ9nndxonasfZx9JmQwbnspETEvvy+B/gMuKqejesGNcXfIE4mxiN18kGoP/ArKQFfFRH/RMQBYD1wDUBEvBoR4yJiKvAbsCuX789dZG3AClJSU0Zdsl5KppYYfyQdfH+JiL+AlcDVhfleBnZHxPJKQe4ePZRHXwHG1TGWU1V17BFxKCJ+BcgJynekM4wVZY29lu09EXhK0g/AImCppIXA9cDMXP4OcKOkNyoLk/QwMJCU0AIQEX9UEr6IWAn0ljSgfmGeUC3r4Yj8g/Qg6ZrcipuAjRGxv1CvUfZ5DcGJ2rG+BkZJGp5/KcwBTuruTUnnSTorDw8gfaG3d1lLu0bV8TeQk4nxfeCOPHwb8Em+7mQvcCOApL7AtcDOPH5Bfh9Cuj7trTw+qLDcWaRu0jLqkvVSMrXE+BFwhdKdjr2AG8jfb0mPkQ5si4oLarftZwI76hzPqag6dkkDJZ0BkM+UjgL25PEyx151zBExJSKGRcQwUlfn4xHxfEQsiYjBuXxOrr8AQOnO1hnA3JykkMsvLFz3NoF07D0qCepitWz74fnzjqShwGjSzUIVc2nX7dlA+7zGcKK7DZrxRbq+aBfpV+ODuexRYGYeHk/6df0n6cvWmssnkfrwN+f3O3s6lu6MP09bB/wM/J3rzOjpeKqMsQ/pVP63pIuCR+Tyfrm8lXSQvr9d7Nvz9p9WKH89fx62kHaGg3o6/u5cL2V7VRtjnrYgx7gNeCqXDSZdYL2D1OW9CbgrT3si198MfAqMbsTYgdk5jk2kbrxbGiX2WrZ3YRnL6OAuRtrdzQn8m/9OZV08lMsXFtbFl8CkRlkPpGvtitv+1sIy+5KOAf3b/a2G2ec1wstPJjAzMzMrKXd9mpmZmZWUEzUzMzOzknKiZmZmZlZSTtTMzMzMSsqJmpmZmVlJOVEzMzMzKyknamZmZmYl5UTNzOw4JH1RxTzTJL3eFe0xs+biRM3M7DgiYlIVs7UA39S7LWbWfJyomVlTkjRM0k5Jr0naJelNSdMlrZe0Oz+TEUkHC/V3SFohqVXSaklnd7L4FuBCSWsl7ZU0vbviMrPTixM1M2tmlwBPkx40PRqYB0wG7gOWdlB/FPBCRIwBfic9A7MjLcDPETEVuAeYX+d2m1mTcKJmZs3s+4jYGhFtpAdPr4n0AOStwLBO6m/Kwxs6qiOpN3A+KQEE6E1K6szMTpkTNTNrZocKw22F8Tag1wnqH+6kzmXA5pz8AYwFttXYTjNrUk7UzMzqqwXYXBgfC2zpobaYWYNzomZmVl8tHJ2YXY7PqJlZlZQuxzAzMzOzsvEZNTMzM7OScqJmZmZmVlJO1MzMzMxKyomamZmZWUk5UTMzMzMrKSdqZmZmZiXlRM3MzMyspJyomZmZmZXUf83J/An8RTsBAAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# get the cumulative average of es:\n",
"i = -19\n",
"ES = es[:i]\n",
"HS = hs[:i]\n",
"pos = np.arange(0,40,step=5)\n",
"es_av = np.log2(np.cumsum(ES[::-1])/np.arange(1,len(ES)+1))\n",
"plt.plot(es_av,label='log_of_cum_av',c='r')\n",
"plt.xticks(pos,HS[::-1][pos])\n",
"plt.hlines(2.0,0,len(es_av),colors='black',linestyles='--')\n",
"plt.title('Base-2 logarithms of the cumulative averages of the ratios for $h\\leq 0.15$')\n",
"plt.xlabel('min $h$')\n",
"plt.ylabel('smoothed $\\operatorname{LR}(h)$')\n",
"plt.grid()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "b1849b79",
"metadata": {},
"source": [
"From the smoothed results above we can see more clearly that the ratios are converging to around $p=2$. In fact, discarding the values corresponding to $h>0.15$ seems to result in the the rolling average converging to a value slightly greater than $2$."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.9"
}
},
"nbformat": 4,
"nbformat_minor": 5
}