Commit d50aeb5f authored by Maximilian Schanner's avatar Maximilian Schanner
Browse files

Added notebook for testdata generation

parent c800cb7e
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Isoline for 1 sigma in polar wander plot"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Imports\n",
"import sys\n",
"import os\n",
"# relative import\n",
"sys.path.append(os.path.abspath('') + '/../')\n",
"\n",
"import numpy as np\n",
"import pandas as pd\n",
"\n",
"from matplotlib import pyplot as plt\n",
"\n",
"from scipy.stats import uniform, gamma\n",
"\n",
"import pyfield\n",
"from corbass.utils import load, nez2dif"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# seed for reproducability\n",
"np.random.seed(161)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Generate synthetic data for tests\n",
"\n",
"This notebook serves the purpose of generating synthetic data for testing the `CORBASS` algorithm. We therefore generate data from a given set of Gauss coefficients and add synthetic errors from the Fisher-von Mises and the gamma distribution."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# some basic parameters\n",
"# the name of the output file\n",
"out = '../dat/synth_data_clean_complete.csv'\n",
"# switch for using locations and incompleteness structure from the example data file\n",
"real_locs = False\n",
"# the number of records to be generated, is only used if real_locs is False\n",
"n_points = 412\n",
"# the fraction of incomplete records, works as a switch if real_locs is False\n",
"r_inc = 0.\n",
"# switch for corrupting the data by noise\n",
"noise = False\n",
"# the error levels to be stored\n",
"ddec = 4.5\n",
"dinc = 4.5\n",
"dint = 8250\n",
"# the average concentration parameter from GEOMAGIA for the interval [750, today] is 650\n",
"kappa = 650\n",
"\n",
"header = f\"# This file was produced using the notebook Gen_Data.ipynb with the following parameters:\\n\" \\\n",
" + f\"# real_locs={real_locs}, n_points={n_points:d}, r_inc={r_inc:.2f}, noise={noise}, ddec={ddec:.2f}, \" \\\n",
" + f\"dinc={dinc:.2f}, dint={dint:d}, kappa={kappa:.1f}\\n\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Sampling from the Fisher-von Mises distribution"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# The sampling process involves some rotations, thus we first define some convenience functions\n",
"def angles(vec):\n",
" return np.arctan2(vec[1], vec[0]), \\\n",
" np.pi/2 - np.arctan2(vec[2], np.sqrt(vec[0]**2 + vec[1]**2))\n",
"\n",
"\n",
"def rot_z(ang):\n",
" return np.array([[np.cos(ang), np.sin(ang), 0],\n",
" [-np.sin(ang), np.cos(ang), 0],\n",
" [0, 0, 1]])\n",
"\n",
"\n",
"def rot_y(ang):\n",
" return np.array([[np.cos(ang), 0, np.sin(ang)],\n",
" [0, 1, 0],\n",
" [-np.sin(ang), 0, np.cos(ang)]])\n",
"\n",
"\n",
"def rotator(vec):\n",
" vec = np.asarray(vec)\n",
" p, t = angles(vec)\n",
" return np.dot(rot_y(t).T, rot_z(p))"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# This cell contains the sampling procedure\n",
"def sample_Fisher(n, mu=(0, 0, 1), kappa=20):\n",
" \"\"\" Generate samples from the Fisher distribution\n",
" \n",
" Parameters:\n",
" -----------\n",
" n : int\n",
" The number of samples to be generated\n",
" mu : array-like of length 3, optional\n",
" A vector pointing towards the center of the distribution. Its length is ignored.\n",
" kappa : float, optional\n",
" The concentration parameter.\n",
" \n",
" Returns:\n",
" --------\n",
" numpy array of shape (3, n) containing the sampled vectors\n",
" \n",
" Reference:\n",
" ----------\n",
" [1]: W. Jakob, \"Numerically stable sampling of the von Mises \n",
" Fisher distribution on S^2 (and other tricks)\",\n",
" http://www.mitsuba-renderer.org/~wenzel/files/vmf.pdf,\n",
" 2015\n",
" \"\"\"\n",
" if kappa <= 0:\n",
" raise ValueError(f\"The concentration parameter has to be positive, but kappa={kappa} was given.\\n\"\n",
" f\"For kappa=0 use a uniform sampler on the sphere.\")\n",
" trafo_mat = rotator(mu)\n",
" \n",
" # sample from the uniform circle, V in [1]\n",
" angles = uniform.rvs(scale=2*np.pi, size=n)\n",
" vs = np.array([np.cos(angles),\n",
" np.sin(angles)])\n",
" # sample W in [1] via inverse cdf sampling\n",
" def inv_cdf(x):\n",
" return 1 + np.log(x + (1-x)*np.exp(-2*kappa))/kappa\n",
" \n",
" unis = uniform.rvs(size=n)\n",
" ws = inv_cdf(unis)\n",
" ret = np.sqrt(1-ws**2)*vs\n",
" res = np.einsum(\"i...,ij->j...\", np.array([ret[0], ret[1], ws]), trafo_mat)\n",
" if n == 1:\n",
" return res.flatten()\n",
" else:\n",
" return res"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We draw some samples and plot the result, to eye-check whether the procedure works:"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlsAAAE9CAYAAAA8gKerAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzdWZBk53XY+f9391xr765q9N5YG9wFQRQJwLJFagvZVIwcQ/FBDzMP8DjGD5rQhGMeZmzPzMs8UDGWQw6LDI9CITpCYngoCVKEPBRpLiBAgmRjIQk00Hv1VvuW692/bx5uVqGqutburu5G9flFNKq7MvPmrcxE3FPnO985yhiDEEIIIYTYG9b9PgEhhBBCiP1Mgi0hhBBCiD0kwZYQQgghxB6SYEsIIYQQYg9JsCWEEEIIsYck2BJCCCGE2EPO/T6BvTA8PGyOHz9+v09DCCGEEA+J119/fc4YM7LRbfsy2Dp+/Dhnzpy536chhBBCiIeEUurqZrfd12VEpdSfKKVmlFJvb3L7LymlGkqpt3p//tW9PkchhBBCiDtxvzNbfwr8EfBnW9zne8aY37w3pyOEEEIIcXfd18yWMeZlYOF+noMQQgghxF76IOxG/EWl1E+UUv9FKfX0/T4ZIYQQQojduN/LiNt5AzhmjGkrpX4D+GvgsY3uqJR6EXgR4OjRo/fuDIUQQgghtvBAZ7aMMU1jTLv3978DXKXU8Cb3/bIx5hljzDMjIxvuvBRCCCGEuOce6GBLKTWqlFK9vz9Lcb7z9/eshBBCCCF27r4uIyql/hz4JWBYKXUD+NeAC2CM+WPgnwL/XCmVASHwO8YYc59OVwghhBBi1+5rsGWM+cI2t/8RRWsIIYQQQogPpAd6GVEIIYQQ4oNOgi0hhBBCiD0kwZYQQgghxB6SYEsIIYQQYg9JsCWEEEIIsYck2BJCCCGE2EMSbAkhhBBC7CEJtoQQQggh9pAEW0IIIYQQe0iCLSGEEEKIPSTBlhBCCCHEHpJgSwghhBBiD0mwJYQQQgixhyTYEkIIIYTYQxJsCSGEEELsIQm2hBBCCCH2kARbQgghhBB7SIItIYQQQog9JMGWEEIIIcQekmBLCCGEEGIPSbAlhBBCCLGHJNgSQgghhNhDEmwJIYQQQuwhCbaEEEIIIfaQBFtCCCGEEHtIgi0hhBBCiD0kwZYQQgghxB6SYEsIIYQQYg9JsCWEEEIIsYck2BJCCCGE2EMSbAkhhBBC7CHnfp+AEEIIIR4MV2bbvHJpnulGxMG+gOdODXFipHq/T+sDTzJbQgghhODKbJuvnrlBJ8oY7QvoRBlfPXODK7Pt+31qH3gSbAkhhBCCVy7N019yqZdcLKWol1z6Sy6vXJq/36f2gfdQLSNKelQIIYTY2HQjYrQvWPO9auAw1Yju0xntH/c1s6WU+hOl1IxS6u1NbldKqX+nlLqolPqpUuoTt/tckh4VQgghNnewL6AdZWu+144yDq4LwMTu3e9lxD8Ffm2L238deKz350XgP9zuE0l6VAghhNjcc6eGWApTmmGKNoZmmLIUpjx3auh+n9oH3n0NtowxLwMLW9zlc8CfmcJrQL9Saux2nmu6EVEN1q6aVgOHaUmPCiGEEJwYqfL5Zw5T6S0dVgKHzz9zWMpt7oIHvWbrEeD6qn/f6H1vcv0dlVIvUmS/OHr06C0HWk6P1kvuyvckPSqEEEK878RIVYKrPXC/lxG3ozb4ntnojsaYLxtjnjHGPDMyMnLL7ZIeFUIIIcT98KAHWzeAI6v+fRiYuJ0DSXpUCCGEEPfDg76M+DfAv1BK/QXwC0DDGHPLEuJOSXpUCCGEEPfafQ22lFJ/DvwSMKyUugH8a8AFMMb8MfB3wG8AF4Eu8N/dnzMVQgghhLg99zXYMsZ8YZvbDfA/3qPTEUIIIcQ+8KA1MX/Qa7aEEEIIIXbsQWxiLsGWEEIIIfaNB7GJuQRbQgghhNg3HsQm5hJsCSGEEGLfeBBnPEqwJYQQQoh940FsYv6g99kSQgghhNjQZrsOP//MYV65NM9U7/u/+vTB+7obUYItIYQQQnzgLO867C+5jPaWDr965sbKdJi7GVzdaSsJCbbugXvd7+NB6y8ihBBC3G2rdx0CK19fuTR/1wOtrYK6nZCarT12r/t9PIj9RYQQQoidujLb5iuvXeWLXz/HV167uun1617tOrwbrSQks7UDd5IpuleR9/16PiGEEGK9271u7iaLtLzrcPk6B3uz63C6ETG67pjVwGFqF0HdvsxsxZneNCreacS8+v53kim61/0+HsT+IkIIIR4ed3Ld3E0W6V7tOrwbrST2ZbC12Ek2fJNv5wNwp+nDe93v40HsLyKEEOLhcSfXzd0kDJZ3HVZ6WaZK4Oyqjmqn7kZQty+XEW1LbbiMBux6ie1O04fPnRriq2durDyuHWUshSm/+vTBXf5UO3Ovn08IIYRY7U6um7tdGly963B56fJrb9y8q5vD7kYriX2Z2bIttebfy1Hx7Syx3Wmm6F5F3vfr+YQQQojV7uS6ebtZpL3cHHY3dvjvy8xWrs2af69+k3dbTLfTTNFWb8bd7vexnXv9fEIIIcSyO1lhud0s0p1sDtvq+r1Rwf6XXr7MaN0n0+w4+Nq3wVYzTDd8k3f7AdjJG7+T3RPS+0oIIcTD4E6X3W4nYXC7S5fbXb/XB3FJljM+12Ghm/DCYyNr7r+VfRlsDVS8lWW09W/y7XwAtnvjt4uo70ZDNCGEEOKDYq9XWNYnMGxr9ytXsP31e30Qd3G2Q73skmR6pfh/+f5b2ZfBlu9Y/O4nj2142158ALaLqKX3lRBCCHF3bJTAmGnFGBNzfKiyq6XL7a7f6wv2m1GKb1vUSy5zrYiLsx2aYYrBoBzP3+x59mWB/L22vhhwrh3x8oVZfnZzia+8dpX3JprS+0oIIYS4CzZqLXFssMJo3d/15rDtivnXF+x7tsVSmDJUdjhzdZE40/iOhe/aWOW+4c2eR4Ktu2D1mzHbCnn14jzNbsrHj/TTiTKuLXa5NtdZ8xjpfSWEEELs3madBTJdXI8P9gVMNyJeuTS/7W7E7XY/rt/h/9RYjZMjFW4uxQSejQLCVPOhQ3XQebbZ8+zLZcR7bXUx4A/OztNXcvnQoTrDtSKYemq0xtnJFoNVX3pfCSGEEHdgs15cjsWu66NXX7/fnWjSiNI1DViXS4/Wb4r7N397FjD0lT1OH6oxXA0wRuebnbMEW3fJ8puxvP5rqfd7fR0dqhAm+aZF+0IIIYTYmc1aS/iO2rY+eqvOADcXQ44MlleOuVmgdmKkymdOH6SzLuBTyrI3O2cJtu6yzSLux8fqmxbtCyGEEGJnNmst8bU3bm64vLhc7L5VZ4DdbmR77tQQX3r5MgudhCTTeI6Fst1NC+Ql2LpNm0XHD8q4HOnrJYQQYr/aqLPAdqN+tgqobqdPl1KAAoMBtendgH0abE03I77y2tU9CzC265t1pzOUbveclj8wjgVTzWILrPT1EkII8TDYLtmxVUC125mMr1ya59hghQ8/0r/yvT/J03izc9uXwZZrWytzkXYTYOw0G7RduvFej8tZH/y9fGGWZjflUF+ApVzp6yWEEGLfW53sODfZZDFM6QveL3bfKqDa7arURoHbQ1kgv9sAY3XA4trwvQuz/NXrN3j+8WF+62OPrDnGnUw036ndLAPeOk5AUy+7XJztrOyIvNvnJ4QQQjxoVhe7PzKwttj90ycHefXyAnBrQLXbVamNAreHtkA+zjJ+cHZ3AUuS57x+tUHZsxmu+7w72SLO1mbIdptu3K3lwE/nmslmxBtXF/n621P8s+dP8PwTB265//rgr15yiZOcVpTuyfkJIYQQe2E3iYbN7rvZ6tP4YrhlQLWbVamNMmFY9qYx1b5tajrXjvj+pQV8x2K0L1hZVtyswdlyk7SLsx3Knk3g2pQcmzTXa3puwPZN0O7UK5fmWWhFfOvcDO9MNGnHKUma8+VXrmx4/us74D46UmEpTHFta0/OTwghhLjblhMNnSjb9rq91X03a3o63Yg4MVLldz95jP/5V5/gdz957LZKa5aDvFaY8vZEg3OTTSqBg+425jZ7zL7MbLXjjL964yZGwemx6i3DIpeHQ68uKH97YonXr2nmWglHB0sErk2UaWqBe8sS3F4Xwb830eT160tYSlHzbDJtWAwTKtrecFl0fYTt2TYnRyocqPnS10sIIcQHwmYZqZfeuslQLViTwdqqdno3q0+73bm/uuToqUP1laXI504NYbLk4SqQz7XBsS0O1jwuzYYMVHyGq8FK0LS+Puv7lxaIkgzXsQDDpZk2RwbLKKV4emxgwzdpL4vgG1FKmmkGKh4KhWsr0lzTjXO+eXbqlg/FRsHfi8+flOBKCCHEB8ZG9dBxlvHyhTl+5fTomt31rTDlqUP1Nfddvsb/9ice2VGx+3adBTYKxLYK8rayL4OtvpLLU2N14kyjoCgUrwYrQdPqF+ub787TilKiNMdOFUMVj+sLIbPthF97+iCeY9/zPllFEGgRJjklzybLDVGaE2U5x0eqG34o7vUOSCGEEOJOrA9mHItbMlLvTLYYrvq3BDc3FrubZq92uvq0XeC0USC2VZC3lX0ZbEFRt3Tm6iKBZ9PoJit1S8tdZkf7AubaERen21QDh6rv0I5zPMfmM08dYHyhS6ahEjj3fAnuibE6aZ7z+rUGrTil6rtYCsqew4cO1TdcFhVCCCE+KDbKKk01Y5SKOTZYWclIzTVjXnh8eOVxc+2ICzNtxuc6LHVTTo/WODpcuSV7tZMExOpM2lwr4uJsh2aYYjDMt6KVQGyuXdw224qYb8WcnVqiHWkAjgyWeHS4ypHhypbPtW+DreFawDPHBnh7ogmoNUHT8nruxdkOrqNY7CYkmcF3LYw2TDZjPnN69L6N13nu1BA3F0P+0RMek82IuXbMYgd++Ylhhqvvp1ilnYMQQogPoo2ySseHKoRptmaO8POPDxO4Ragy1444M76EpeDEcIWxvoB3p1p0kpwnD9U3TIxsVZO1HAskWc6Zq4uUPAffsTAKXr4wxz94fJi5ds6Z8SXKnk3JsZhrJ9xc0vQFNrmBH4+HnJ1o8vufeWzLn/e+BltKqV8D/hCwgf9ojPm/1t3+S8BLwJXet/7SGPN/bHfcXBu0MXiOzcmRKp8+Ocj4YsjX3rjJwb6A4wMlXr28wPX5LmmmiRKNZUHJsbi+2KXsO5vu3LsXY3BWp0A9x+bZE0PMtSJK7tq369p8h4lGxBe/fk5G8gghhPjA2KxfZTvO1iQ6ljNgABdm2lgKtDE8dqDKcC1gsOJTCZwNkyPb1WQtby67PNsm8GwU0E01zxzv5+2bTd6ZbBG49kqHggszTTJtUArmuikVz6HiOfiO4qWfTqEcb9PZiPet9YNSygb+PfDrwGngC0qp0xvc9XvGmI/1/mwbaAGEac5fvXmDtycanBoqAqvV20NfvbzAp08OkuYalOJgn89Q1ce2LVzHouoXu/6++PVzfOW1qyvbTnezLfVOrd+e+lsfe2RNu4nxuTY/Gl/kUD3Y83MRQggh7qb1LYtg4x2Dy8mHSuAwuRRRCxyeOTawpmH39CYrPKuzZ8vlN6tbOS0fO840cZrjuRbPHO9nuBpweqzGXDNmthXhOYr5dsRMM8YCLAWOZeHaFv0lF9+x0cag/Ep9wxPh/vbZeha4aIy5bIxJgL8APnc3DmwM1HyXJNP8Ra856PKLnWQ5l2fb/MdXxym5NiNVn2NDVZ4arXNiuEo9cGlH+YYB1XZv3F5a/YGbakScn25RdW3Oz7b50ZUFkiy/Z+cihBBC3Ind9KtcTj587uOP8PShvpVAC7Zu2L1Vvy14f6UKwHdsHh2urJTqBK7D848PM1TxmWsnLEUpB2o+geeQ5gaLYuPa9aWQZpTh2wpl2S6buJ/B1iPA9VX/vtH73nq/qJT6iVLqvyilnt7pwftKLkopJpshF2ZaQFEAd+bqYm86t2Gg6pFpTZrlNKMU37EoOTaDFZezU02++e40Z6eaaK1Xlg63euP22vIH7rc/8QjNMKNccugLXOJMc+bqIlGa3bNzEUIIIW7X+gRCJXC2nWW8WYB2fKDEV167estq1HL2bK4V8drlef7+7BQvn5/FttauVH3iaB+NMOX7l+aZaYYrx/2tjz3C53/uMLalmG8n2JbCGE2Wa9LcoI0mzzS2pVjoJgBms3O/nzVbaoPvrT/RN4Bjxpi2Uuo3gL8GNqxCU0q9CLwIUBseQylFnmtyDW9cX6JW8mhHKSXPQQF9ZY9HRyo0uimObfHC4yO0o4y/mblJybVJUo2t4OJ0izevLjDWV+JTp4b2dEzPTr1yaZ7huo+FQilFyS3GMZ2dbPH84yP39FyEEEI83G63lnmrHYObHXN9S4enx2q8enlhw7qs504N8aWXLzM+16FedvFti6UwZaYV89dv3Vy12zCn6tu8N9nk4mybE8NlTo/W+dNXr/CziSYl16LqFY3OfcdhoALNKEUbRb3scHK4QpIbtoi17mtm6wZwZNW/DwMTq+9gjGkaY9q9v/8d4CqlhtmAMebLxphnjDHPlGr9tKOUqwtdAsdCAc0w5cJMm06U0E1yHh0p0oWfOjVInOmVyLrqOzSjlHNTTX56o0GSaQLXptnblnp1obNnY3p2aroR8fRYjW6SE6U5BgPGMNeOZSSPEEKIe2Yvapm3Oub6eubxxXDT8p4TI1VG6z71clFW5Hs2n350iGODFd66vkQ1cFZ2OMaZpuQ5lD0breGtGw2+e3GOsmtRDTz6Si4Hqj4nRyrUSx7PHB3gUH+J40NVBqo+nzo1CKiNkkjA/Q22fgw8ppQ6oZTygN8B/mb1HZRSo0oVJ6+UepbifLctStIGZloxCrCU4smxOvWyi2Mp5jrpSgHcXDvinckWYFZ2KbajlFQbwjQHYLGbEmeaeuBwfKgYgbObtOdeONgX4DsOzxzvx3MtGmGKUfDCY8OyG1EIIcQ9sxe1zLs55lblPVdm27x5vUGc5tRLLkNll4uzHX5weY6bSxHX5jsr85CbUYZSiv6yR5jmpLnGAppxRsm1GawG1MsutcAhyTRhorEUTDVDrsx2mFgKMTpPN/uZ7tsyojEmU0r9C+DrFK0f/sQY845S6n/o3f7HwD8F/rlSKgNC4HeMMZvn6VZZ7MZUA5e+ksenTg4xXAs4NVzmu+fn8Gyb2VbI9y8toIBPnhykE2V86XtXGKh4DCrF+ekWSilcS+Eqi9H+0obbUu+H5e2q/SWXZ48PrjRz+9zHNip5E0IIIfbGZi0c7qQH5G6OudkcRMcqOsD7jgUKbi50+fa7bWwbLBQVz+Z7F+bwHZujQyXaUYal4EDV59pCFwNUfIdOVCReAsciyXKePtSHpQw/uLxI1bPpK7l0opRvn5sDrfPNfqb72mertzT4d+u+98er/v5HwB/t9riBa+G7DlGaM9YXrFSH+Y7DC48NUwkcfnB2nr6Sy4cO1Vd2NmhtCFyLTCuGq712GaYYbP3oSOW+1GdtZK8HYQshhBA7sVGwc6c9IHczSHo5+QDFHMV3JlvMNWPqJYcDNZ8s17w72WQxTNC5wVM2vq3QpujH2U1S5to2Fd+mP3CpBi6OXQQN/SWXRjfl/FSLVpSQaMN0IyTVhr7AoZPmLC2GBJ7NwZqH5Zf7NvuZ9mUH+W6SU/Nt5js5F2daXJhucWywzJGhysqA5uXI2Vq1xDpU9ZjrxDx5sMZUM2RqKcJ3LU6NVPDsO5uReLebod7OLMR70ZBVCCHEw2N1sFMNHK7Nd/jR+CLPHhvYsJHoThwfKPGl711Ba8NQ1WOsHmDZFk+P1fi33zjHm9cbgOFjR/qLHYPPHOalt27y8oU5hqs+Lzw+zHfOzXJxpsNo3SfTmiQ1oMDKDCMVlyjLmW3FVAOHY0Nlap7Du9MtlroJJdcm00Vddr3k0I4SFrspQ1WXsu9webZDrg0jNZ+Rqk+aG7ppjrLdTbMx97Nma88oYLDskWtDO9YYbZjrpKxegNyoodpYPaAVpnz73CxpnjNS8/AdRTNMCdPstuuz7mUz1Af5HIQQQuwv61s4TDQinj02wPGR6m3VcF2ZbfPq5QVOj9UYqhUJkHenW5waKvF370zxg8sLBI6F79r88PICX3r5MgBDtYBfOT3KLz1xgAP1ErkpdgdeXexiUNgKMBBnmhuLIUvdhCjN0VpjDJR9m6MDZaqBw+mxOv/wiREGqh6ubZPk0F/26Cv7lBybJNOkmWaqGTPfTopj6K0rnPZlZsu2FLOdhL6SR38JbBsGKz7HhyorOxTWR+PtKGOiETLZCAkTjemt6R6s+zx9qM6l2c7KuJ/dZoS2mix+rzJLuzkHyYAJIYTYqdUrLV/8+rk7quFafa06PlwcsxmmfPfiPEmuGSh7BK5NK0pphBk/vDxPO87oC1yeOvR+A/eyazOZ5mhjcCwLA+jebVobMhRKaSzLou47DNUCfu+zT6w5l//1r35GnndpeBZV3yXLDVcXOqRaowCjNVGaMd+JURhg85ryfRlsRZnmxkKXkmdT9R2i1DBShXcmG1yeafPDy8WbOVz1CNOMdlwU012ebROlBttSaAPdNGeqEdGMMw4PlPjkyaHbSomem2zSiFJavTXoR0cqDFb8PR8ivTpo+tnNBp842ge8vwa+0f8A282SEkIIITazm3qrjWxWHD/ZCKn6Dv1lj1aUcm2hizaadpzx5rVFAtcmzXI+dmywdx4lri+GAERpjuMoTGbQBlBgMNiqqN36/uV5PprcWtu+FKbYlqLqu+R5UeO10E3JdRHApblmKTS4tmKo4nPhQdyNuJe0NnhO8ULMtGICx6IRZsw0IiaWQiaWIgyGU8NVHh+r8eLzJ3nl0jzNKC9akilwrWLnwWKY0opT+ksel2ZazHdTZlsR7001+Njh/qKNf5jSF7g8eah+Sxboymyb8fkujqXoK7vEqebM+BJPjlY5MlTZs9dgfdB0frrF9y8t8OlHh1bGEWz0P8CDkIUTQgjxwXR8oMSXX7lCpjXDVb+ot7IsfvXpgztaNVkfrM21It6eaNIIMzpxjtaGRpShjWGpm2FbRaATuBbfv7wAQJhpri906SQZxwbLRKkm1y43lyKyXIMqyo0MRbww34m5tthdKatZPsfpRogB+gKHq/NdFrtFn05tINNFnixwwLMVmTaYNO5s9rrsy2DLcyzqJYeZRrzyos60QqYaMRaGWslDKcX4YpdUG1566yaZhjjL8XsvWqrz4sUzkGvwHcU33p3hxFCZkmNxbqrNYifFsRS1wKURplQ8m68uhmuyQK9cmuf0aI33pltEqSZwLKI0592pFl949uievQbrg6YPHarz/UvzvD3R5IXH/JV2EesL/vdiG68QQoj9b7ne6qmDNSabEfPthMVOyj97/gTAjlZNVpf4RGnGa5cXMMBzpwZ58/oSF2c7OAoSrdHGUHXtov+l79COU16/vsSBus+RoTJjdZ+fTRS9NGuBy+H+Eu0kY7GbYCuF61jEmSHXmrG6z1+/dZOZVsy1hQ43FiJmmhFJpnHson/n8h8A14JMQ5QZMJqBqo3RecYm9mWw5ViKR/rLHB4o04lS3p5oUQ0cbKWwLYso1VT8Ykp3mue8eb3BZ04fxHfsIrByoRNlGGOwFFR8myQ3lF2bRpSR5pp6ySPNNWkORwYrhGnOVCvm9Fh9TRbovYkmzSglSnMWuwll1+ZAPaAeuHuaKVofNA3XAj55sviwbtUu4k5TwEIIIR5Oa+qtRt6vtxpfDNd0eofNV01Wtzb6weU56uVei6ZqwGDF50dXFnhvukWWGQ7WfY4MlqkGLmGaY1kWB+suv/mRQyvHOzJY4fxMi2aYMVhxSRZyPMciyw1xmpMow3DVYbYdc2WuS5zlLHQS2lFGrjU5gC6qsZZ7Fzhq7bxBbSBONo2zisfc6Yv7IKr4Dp89PQrAd87PcGQwJ8k0nSjDshRKFb2zAK7Md5hqRhwfLFEPbMI0I8kMcVakCF1bFUMm2xEDFZ9OlJJow2MjVa4vhisveJblXJlt0wxTDGZlbM61xS62Uoz1lYgyTZhkjNUDjgzv3RIibBw0Ba7DZ06PbtmUdaONA3fS8kIIIcTDYbuVkZ2umiwX3K9v0TRcC/i1D49xbKiycm2t+A5hmhMmRe21byteuzxPM0qxKGqzWnHG848Po4BGNyNKNV2dEfguFtBNNG9cW8KzLWwLWlFOrg1KKTwb8l46S1MEWdm6MYjaGMJUoxzX3+y12ZetH3JtVmYXzjVjPnliAEspPNsi14Y013TjnDTLUSgGSi7vTbep+Q7GUGS0rGId1nMsHukLiHLDbCui7Ds8NlLFdWxcW+HYinaUcnmug21Z+L0tqV89c4OX3rrJU6O1IurNNIFrYSnF2anWjmcYXpltbzjNfDubTUff7nlvZxK7EEIIsVFLpeWVka1u2+3xnjxU58XnTpBpw2w7xrMVTx6sYSuLhW5CnBWjdq7Md7k406YeOJRchygz/PLpAzx+sEp/2SfNNFGmaccp7Thjrp2w2C3qwYwxZL0gy7YoOtFz66hpBeQGbKVQluNt9rPsy8yWZSnenSh2ANo2TDRiPvJIlTPXcpJ2XMw8UmBbFn0lh5MjNRzbYny+w+mxOo04ZbGTEGeGimdjWTYnhypMLEX84okhyr7Na5cX0MbQiTNeb0QoBYMVlzDVxcxC2+YHl+f57OmDVAOHt64vcWm2jQEGy+62PwPc2c7AO+kyfzsNU4UQQjzctlsZ2e2qyVbHOzFS5fBgeU3BfZzlvDfdRgFznRjXUmCslX5fUNSBdRONQpNpTZZrcgOusogo+mcZwFJFYJVrQ+DaADgWpHrtOSqg7NkMVj1Mnsab/Sz7Mtiq+g4aQ5JrFIrzUy36yg7/8PERLs10ePPGErHOSXLNYjfjpxNLnByq0IlTRusBQxWfk0NVOknGTCum0U0ZrpYZqLicm2kBitG6R6Y1JdcmTEN8BxpRzkcOlxmuBuheQ7V2lIGBLDecHKmCMRjFjoKmne4M3GyHhwRNQggh7pXtfsnfbQJgu+Otv8Z98evn+NSpQS7NdVnqpigMYab58f8xoKQAACAASURBVPgiM+2YnzvSTzXweP7xYf7Ta1cJXIfM1ni2RTMslh2XlwrzVSksW0GS5zi2umVJ0XcUjwyUaMfZgzmIei/dXOzyg0vzHBkoMdZfwncsLs11+Pb5OUbrHoGrSHOF25t/tNhOeL2TUnYtmlFCqg3X5jvFIOvAYaTq0Y5zxvpKvPD4CO0o45WLc3zi6ADHh6u8dnmeOCuanM13Ux6jSHV+7Eg/S2HK5dk2gWejgO6qzNdy0LRZsLSTnYHSF0sIIcSDYqtf8m8nAbCbxxzsTUf55Ikhphsh7062cG2Lqu/Q7Kb85ZsTHBsq88tPHsBd7iQQZaRZRtKLriwFng1JBkoVhfGWpfCUw2DZpRVlhGlGnBX3z4xhvhVT8m1M3Gludm77MtjKtEFrzXtTLSp+SOBYZFlObDQLHUWYaLLcoLV5f/1VG1LbcHGmTdW16WaaRpgybSmqvs1imOHZislGyDNH+9HaMNmMOD5c5dGRCmeuLhJ4No1uslIf9flnDgPwL7/2ExY7KUopjg6WgCJoOjfZ5P/+xjm+d36O4brP02O1lTE6n3/m8I52BkpfLCGEEB8Eez2dZPWy42KnGKMTakOmNWlusNBEaU43zlFKMdmIyLQh6iVLDEU3A1B4TrE8+OiBKv/oyYP8zU9uYqkiYZJFveyWAZ1DlGs+PNLHq1uc274MttJcM99J8R2LpBc0NcMUWxWjfJaL3vJ1lW55rqn7Hlop+kounTijEabMtBL6Sw59FY8o1Xzj3VkO9flcX+igUDSjlCTPuT7VJTcGBRwbLPO1N27iWKA1HB4oY1sw2Yi4MH2D4aqPZRVjhIbrPhaK1682eOZ4/8ocqa2awy2TvlhCCCEeVMsB1rnJJuPzXU6P1jg6XNmTVZjVy47znQTXtrBtiyjJepvkHBY7Ce9NtzgyWOb6QpeFTrzSVcAC6BXGe26x03G0v0S95HJ4oMzZyQaWokh/5UXneM+2GK56RLnGrg4d2uzc9mWwpZQiznLiXNNNMrJeVOVYRUPRZFWU5VgKSxXNS9Mcyr6Da1scGShzdaFLmGlsneM6RXZrsOJR8mwmGiFg0V9KidKMc1NtEq05VPe4ONPm/EybJw9WSbKiAG++HdGOc0quje9ajC90ONxfohEmjNZLqN7W1ouzHZ49Psi5ySY3F8MNm8Ot/mBKXywhhBAPotVlLo2oaAL+3nTR93K4Vlyj7vYqzPKy41d/dI1unFEveUw0Qjzboh2ltLOcdyaaWBbEaU7JczBxRtbro6UsC2U0NtCJU356fZF3J5qUHJhrxRwZLK/UbcWpphLYZHlx3VWuv+kPsi+DLQw4SqGBNC+WCpcL32y19q5aG3LAtYudBq0oY6y/xEw7xrUt8lzj2u93yGhFGfXAphnl/PqHhplrJ5ybbuM6FlXH5mYjwbcthqouU82YdpxxdKDEdDum5DorRfWeXfTemmiGRL2AbLYdsdRNiZKcJMv55MnhDZvDPb/q/KUvlhBCiHttqyXB5du+eXYa37H40KE6rSijr+wSpZqLsx2Ga8FdX4VZfU5KQSfJ0BTTYRrdmDQvklKuVcQDSW5wLYPv2ljaUPVs4lzTjorYoWQ0lsqAnGtxitGGMMlQgGVZeA4kqcagSfKclazJBvZlsGWpIkOV5Jo019iWIkk1jl1M/nZXbd/UgA1FgOZYtOMU1ypxfSEkB3JjKNmKwYpXFMYlOWXPYqDs8uHDA3zrvWns3uvbCIuOs37gEKYGyzJUfYf5TkKcak4f7UOhCNNiBqNRUHJs5lsR890Ux4KKZ9MIU1phkTFbnbGqBg7vTjT5ymtX13zAb7fFgxBCiP1pL+ujttqYBe+P5YEiXXTm6iKOrYh7I+uaUbFp726twlyZbfPSWzd5+cIcw1Wf02M1xvpKtMOURpj26rMVjl2sYkWpJjMGZYpxOzaaasnDdyxSXZQCOTY4tqLZm8NoemP75jspdd8hN5o0zwmTvNgAp2x42Mb1KKWwbcVw4JHkGs9W3FiK8G1FO8nIe4GWU9TBYalixlFgKYyGqWZRNGdbUA9cwiSnmxQZrbJnUfFdPnWyzrW5Dhdm2jh20ac2zTWm1+Yh0zlDVY+Ka3N2somy4LVLc7iOReDafORQjYlGQn/FI81yGlFKnGmOHqjy0SP9vD3R5OxkiwP1oqB+rhXxwyvzXFsIme3EtxTTb9UVXgghxMNjr3epb7UxC1i5ra/skaSakueQ5jndJCdKc2qBs7KRbPUqzO0EiMs/69s3FgmTjIszCZfn2nzkUA2jFENVj9Njffzwyjxz7RijoRln1H0Hz1UshhkpYMKU3LOxlWKk5tGJM+JU94rmLTKjiXNDyVXYdlGM3UlyMg06yfEcA5a9aUy1L4OtwLV54mCNqUZEX6nY2jla91FKsRSmWBYEtqKbmiKjZUFfyaG/5NFJMgLHZrjic7MRUfZsyp5NJ8pIMs1TYzX++0+d4PBgmf/tpXcIHAtlMuY6GWlucCxoRgnVwKPiWtxY6hK4FgroJDmuNhyoetxsJAxXPR47UOV7F+d4cqzOYweqDFeLKP/pMcN3z8/R7GW4Xru8wFQz4thg6ZZi+pfeuslQLdizHR5CCCE+OPZ6l/pOx/I8OlLhzPgSJddCa3hytMq7Uy3qgUslcNaswiwHTVprJpsRr19b4O/fmeLF506saV5qW0Vt1Vw7YSlMmW6ERXun2Q4V3ybONN0k5+WLCwyWXTQwPtem1et5qXqb1jpJhoHekiC4toVlWViq6NUZphpL65USJG2KxIxSCs9RdOMi2wVFsiZLNEpZ9mav2b4MtjKtKXsOv/mRMQLXYXy+Q5SkfP3sTO+NUiuTu4v7Q5TmBHVFMyra/6e9wdNhptE6x7Utfv8zjxFp+OH4IuOLIUNllzTLODsZETgWga2KN9DAkX6fqWZMrg0lz8FSiv5ykWlrx5rjNZdHD1T5vc8+wVCt6A2yesnQdxxeeGyYSuCsDOM0GIaqwZpi+lPDZV6+MMevnB6VPltCCCH2fJf6dhuzlm8brgY8c7xYqTEYjgxV+MKzR9cEWMtlMePzHcquxVQzoezZHKgFNLopf/itCzx6oEY9cLgw0+LcdNFYfLBcHH+qGeNaEKY57ThDKYXWhtgYkizn1EiFVpRTC2wyrVf6Yy2XEjkUK1hKKQ7WfCYbEXOdhCjNihoviuVEVNF3KzeG6UZEnBvWNZOnCOU2ti+DrYO1AM+xeOXiPKN9AacPVvjW+RZZplFWUfhm9aZ2L/fW0AYuznYpeTZxltOKMowxeK6FNopOkvMH3zzPyeEaQ1WPJM2Z76ak2nCg7tOKcgAOBA6ebdFf8fGcHINhoZNQ9W1yDbY29JddXnh8ZOWDv7rIPc4y3plsMdeMef7xYY4PlABFnOZ005yFTsxQNcB3LRphyjuTLYarvvTZEkIIAez9LvXdjOXxbJuTw9VbEgDrlzpfv7bAu82Ysb6AwC1GDPaVXK4udAgWuliWxY3FLloXq0dL3ZR64FINHGabEYoi4Cq5DpYFGEWaGyYbEZ5jYykL11LE66cbWqAwxKnm8lybNDd4DpRcG0VOkkOcGyyKdlE61dgWaxI2O7Evg61GmJJkmqrvsNRJ+H9encVzbALfQWuD1jm5eX+gpALSrIhS0zAr0oaOJtWQJxrXViRZMd6n5FkkueHN64u045RLs136Sg6DFR9joJvmPP/YMNXAY3y+QzsqOtPmucG1LZJcs9RNefn8LHGm+bffOIcBWmHKuakG082YwwNlPvxIndfHF/nzH12n6tuM9Qf0BS5X5roAlD0Hz7GYa8a88Pjwmp9f+mwJIcTDa693qd/pWJ4rs23+4BvnWWgnDNd8Hh2pMFILuLnYpRmlK+U0UVZscJtuxRys+yx2i/6ZVi97dWG2w6PDZaYaRWNSW0GaF3VUvqMYrnpkeZGEuDrXIUz1SkcC3YsBMg3z3Yyab2OhUMqQZMUORotiQ53Tu3ZnGjCQ57t/zfZlsBVnGtVrTHp+pkU7ySkZg2fbLMXJmmamhrVTvFXvP504x3UsMmNIMkOcg60MV2Y7nBypMtdO6KY5Jc9GKcVsK+ZA3efEUJmLM236Kz7jsx3CLGeg7LIU56S5Jkoz4kzTCD2eGq3yg8sLKOCTJweZ7cTUA5dH+nwuzXWZ7yb0BQ7aGK7OhRwfLnN8qMRsO6Hs5bzw2DBmtEbgrn0bpc+WEEI8vLYLhu7Wc9zOWJ7ljNZ8J2ak5hNnmjNXFzk1Usa1bZa6KcYUwVOYZBys+0wsRcy0InJtaEUZuTbYlsIYw3Qrpq/s0ZjvYgyUPUU1sDFGkWvwHMXp0Ro3FopEhW2plZZQq7Xj9yMotzcDMTdFUJZqfcv9N2Y2vduWwZZSatM5P8t3ASaNMY/v6DzuEdtSlHpTunNtsFC045z+so1tKfL1reNXce3ixc00kBuMMVi9Gilj4OpCl0ac0R84xGnOgZpPUQVWFNgZY7gw0+FzH+1DAeNzHcbnu1hK4VjFGzhY8fn0o0NcnO0wUPYwwOW5LkmmqZddzlxb4lB/mVxHlPzieQ71+SyFKUMVn4N1i3/zj59emasofbaEEEKsdjtzCO+F5eL9kVpQ7FTsXavnOxnPHOvnx+NLzLZjhioexw7WaCUZi+2E64shnmOR98bxpTmkWU4rTHn8YI1KrwSoGKFjsGxFK0oZrPhcmu1wpL/EuZk2WS/QKroRvD9RZvVK1+oSo+XvLxfK367tMluXjDEf3+oOSqk37+D574lUF+m/pW5Cprd+0dK8iIxNWrTsL6rdiq/LL/xCO6ETpRiKVOZIzaebZix1U9Lc8OiBCsdHqqR5zo+vxL3u9DmmN1k8cBX/9ewME82IsmtRL7mUPZuDfSXiJKcZppwcLlpERL0eHoNVH8ex+IUTQ1QCZ83Uc+mzJYQQ4oNguXh/eaciFEt+c+2Yjxzu53//x6cZXwzX7K7/01evMNmcxBjwbYXWFrp3XQfDxbk2JceiGvhEqSZKcnxVFLaHac7NpZCjAyWODZa4Mtft9dFSpJnBsRSuVcxUTnVxjY83SMjcSaAF2wdbv72DY+zkPvdUrk1vLE9OoxdgAStfN3vRloMwz3HwbEMjzLDtYn12ufmpAjKKRmglBxbDoj/W8aEyNd/h2kJI3Xd47fI8040QbQxRmhGlGt+xyHPDQjctjq0U7d4sxTDRnB6rc7YRUfIcwiynL3BY6CQcrPuEWbEjcqOs1YP6G4wQQgix2nLx/vJOxYuzHWZbEUMVf6WI/vl1j6mVPD730TFeubTAUjembDkkWU4nyRksu1hWsXrl2kXwlNmqt94ExhTLgRNLEf/tzx/h1UtzvDvZIsk0tgMjFY92nJPGm/Yj3YXNO8hvuk2xd5KXtzv0Tu5zr+XGMD7f4d2JBu04u2VEz2pq1VdDEVD1l1wOD5apBjauVbxEjgLLVit9OXxbYdkWg2WXTpzz3nS7aEo6VGKg4hFnmncmGyR5UT8WuBbaGLSCdpSjUNiWQueGZpwzWve42Yg4OVLh2eMDzDVjSr7DZ54coeQ7zDVjnhqrSUsHIYQQH1jPnRpiKUxphsUS3+nROh89PMDvf/bxTa9ttgUTjZihilu0burVUJVdG9e28Wwb31FFv6vcUHKX2zzo3q5Bw3w34UeXF/jFk0OcGC4zUPE40lfCtRVKbVVttRu3WbO1FaXUl40xL97u4/dSxSvWgGfbCUZD2bdXit/WvxIGKHsWaarp9ThluhWz0E2o+EWXW9uCkle0bgjznP6Sg2NbhGmObVmM9vl0kpzPPDVKkuecGV+i7Nloo0hyjdZgWYrcFE9gKBqq+Y7Nob6AKMtxbJs4zXnx+ZMrtVjLTdyeP1i/pVHpXo5iEEIIIfbCbktfrsy2mWnFNMKU/pLHcDXjxmIXTFHC00k0eS/Asi1FK8ppxSkKRcV38B2bJNfYVs47Uw3CPGeg7HFyuMLF2Q7jc10GKg59JZe5dkwnub0Fwy1yOsCd7Ub80h08dk8Vy3qVYqunLrrJbhW0RoleWVoMXIuSa9GJUrQxeLaFRTHd23MtHBtcxyI38MhAiadG++imGT+9vsg7Ew1acYZjQ5LnOMszGR2rKNrrFeLZqmj/byh6dA1UfE6P1W+pxdrqw/fl713uzVzMOT9jc3aisRKoCSGEeLjtxS/kd+uYuyl9eeXSPMcGK4z1BVyc7VALXFzbItOGMNUkWU6WQ9Abo3N4wOfSbI6tDEmu6cQ5mdZobchtqAcOcWqYbad8/Eg/v3hiiKlWzHuTzd5K1u0FW2bNl1vddrBljHn9dh97L4zPtQl7RVpbbD4sWvX3GpzS2+bZjXNSDUmUE3iK/pKP1ppuqsnzokA+cC0GSxXCNGe6EeLaVjHVvOSubFl9dKTCdCumE2d0krxYzTXFtlJtNBaKiUbEU6P1Xe0gfOmtm1zu7WTsL3tFM7bZDi+9dZPf++wTd/rSCSGE+ADbi9mI2x3zTgOx1Y93rCJqyTX87GaDTxztY6RWWum/debKHF8/O0MnySi7NmUXGlHO4arPZ0+PMvn9qyyFKXaWFTsOe93es0zzwysLPD5Sob8S8O1zs/z2Jx7hkyNDDFUc/vOZGyslRbdp0wTXdq0f/nar5zXG/JPbP6e9E6U51xfDHW3VXNll2Fvei7K1j0gSQ8PEaIqJ5b5T1FpZlsWV+S6PWha2ZfGhQ3UuzHS5ttClGtgEtg1K8fEjRQHgjYUOyirmQxVLjEXH2qlGyLfOz/Cpk0ObnuP6D/Grl+YYLPsEvS2zaa5phil/+WYxI1GWFIUQ4uG1F7MRtxs+vToQuzbf4V+9M8XRgTJPHrq1DGa91YGca8P3L73ff9J3LL5/aYFPPzq0Emy1U83HD/fjeRY3FkIWuwnDVZfDAyVGaiWGqj5xVjQ3Nawdq5PlhrNTLX7uqItjKc5OtjhQLzHfyRgo+4RJyE56lq4PymwFxphNQ47tMltf7H39b4BR4D/1/v0FYHwH53Mf3dq0bCM5xYuw2YubA2FqcG3oLzv4blGM59hFgNZf8RjEsBhmjPV6YbWjjI7K+dBYjd/77BO89NZNvvbGTQLXZqjiMtuKmWrGWBYM+C6OZfFf353hnYkmv/fLj/H8EwdWnn+j3yamlmJ826LkObSilGsLXaBY8uzIbEQhhHio7cVsxK2OuToQm2tHvDfVxlaKZpTu6Jq0+vGvXWmu6T/5SJ/P+ekWf/XGTT50uI+xerAyOeVAvQTA35+doh44TDdjXrsyTytOgaJ9U74u/FluVvrG9QWGyh4/Gp9nthXRjnMMujfq5/0O88ss9f6IHlsVddxpplfaRPSVHK6bPN3s9dsy2DLGfBdAKfV/GmNeWHXT3yqlXt7qsffT8lruju+/RVS2fNNy47NWlBO4RWBzuL/E8aEK4/MdbKXor/gM9SLvpW5CjuLESDFs+ueODvCH37rAtcUQrYs2EGmuaUQZnudQL7k0o5Qvv3KFw4PllQ/lRr9NHB0ucWMxohJ4zLaK/3my3HBipCKzEYUQ4iG3F7MRtzrm6kDs4myHsmfjOxbNKN3RNWn145thSl/JBQOTjZClrs0jAwHXFyLevtngratLjPX5dJMiRTLXjpjvxJyfapHkGs9W1HyXdpShuLVee3m/YJzROz+fZpQRpTmtKMe1FXF2a7LG7u1YtOltcovXpmjCJAdl2Zu9flu2flhlRCl1cvkfSqkTwMgOH3vPWQos647WXW+RG1gKM5IsJ0pyFrsJE40Q2ypaRSz39jIUX3NtGFj1oXz+iQN8+tFhPnSoD98pdk2UfJvAc4hSjedYKBStMOMPvnGeL379HF957SrnJptUg7Ux8bPHBym7Fu2o6Ko7347JjeHYYBkoftuYltmIQgjxUFrdXkEbQzNMWQpTnju1ebnKnRxzORCDIljyXYso09SC4hq43TXpYG/p8bUr89xcDLkw3WKhHReDpT2HKCnqoBUKzylqnV96a4KvnbnOd8/NUrJtuklG4FhMNiKiNEcpxYGaj7Wuimp1GqaTGBY6CUmWc3SgTDVwcFQRFq3vmKVNEXDp3m2OtfY+YWZQdyHY+p+A7yilvqOU+g7wbeD3dvjYTSmlfk0pdU4pdVEp9b9scLtSSv273u0/VUp9YodHpuI5ODv96baxHLSt7ChUoI1hsZvw7XdnWApTaoHNzaWQH12Z56c3lphvxyyGKVdm20CxHPjW9aViJE/JpRY4aA2upchyTZobtDa04pSFTsJoX0Anyhif73JtvrPmfHzH4eePD2JZFv1ll4N9AccGy1ya7TLXjmQ2ohBCPMSW2ytUest8lcC549KSrY753Kkhri50+M75GW4sdDk70WBisUs7Tvn7s1O8fGF2y+vx8YESPxpfpNlNOTIQ0IpSLvWunZ0o4dJsh4rv4NmKxW5KK854pN/nRiNkqhFRLTkcGShR9ovZirPtmMGy20tibM2yFPOdlPH5Dkf7S4S5xlJFS1SLojbLBmqBw2jdx7KKfl6pfn9Z8X1q059yR7sRjTH/n1LqMeDJ3rfeM8bEO3nsZpRSNvDvgc8CN4AfK6X+xhhzdtXdfh14rPfnF4D/0Pu69bGBOM3Jbm8H55YyU6QLA9diuOoTZpqya/Gd8wscrPmUPQejDfPtmMCz+P3//BM+fKhOnBt810YZ6A9cJpdCWlFKWyk8x6YZJliWYrQeMFz1sZSiXnI5PVrj3akWgxV/zezDgYrHp/vLJFnOmauLlDwHjOHtiSYnh6sb7myU3lxCCPFw2IvJIlsd0/Tm2Q1WPG4uhSx2U/pKDuXApdlNmWrGXJltrzx+9fVofL7D4yMVotzQilIeO1gjznKmmzFpbqgFy32wEhzbwlOKzIDvWDw5VidwbKjDlbkuB+o+cy1wHZuFbsKh/oBmmNGKsg03zOVaY1sWBmjGKbk2+I6Fr4o+mWleZLC0gZLv8nTV42c3W2uOtV1AB7tr/fAY8AQQAB9VSmGM+bNdPH69Z4GLyx3olVJ/AXwOWB1sfQ74M2OMAV5TSvUrpcaMMZNbHdiwUcR5dyjAsy0Cz8HptaaPcsPBqsfVhS6WKgZSWwrCVHN0oMwPxxfpL7k8NVbl0mwICny3aLQWpxowWKpoxuq7No+OVFae7+hwhU6Sr/w2sdwA7mtv3KQaOFjK5ZljA1yc7dAMUwxmw99g9mIrsBBCCPHKpXmOD1X4yOF+AL55dorriyHXF0PqpQwFvHV9iX/5//6Ef/LxwxwfKPHq5YWV69Hr1xZoWBY/f2yA4VqxKqON4YeX5njrRoNOnNFNUtphRq7AsyyiNKPs2fzk2iJhqgncYkh1J1Z0k4yFbkqmDQ2LXiPyouhqdc9SAyTZ/8/em8ZYmp33fb9z3vWutW+97z3dMySHwxE5JKdHkiVKtiyFBORElhPGMYSQMSQgDmQgRhLb+RLEHyTE2QCJTgQIhB0ztiNKtiVRlCx4FqlNDskROdMz3dP7VnvVrbu9+zn5cO69XVVd1V3V3TPT7Dk/YNC13Hpv3WVwnnqe//P/a3ynIE4Fl7opoSNQGuKesr5fT5R9yf7REs1uNhDLrw+r7oUDPfA2IgBCiH8I/BhwGvh9TMfpVeBhiq29wI11n9/k7q7VVrfZC9y72NKarLch4Ep23eFyhdFoba7XRO97Sa7QUcpaN6MWuLx+ZRm3Z356YLTMQislKwrTAfPNLPnAWNmkmh8a5pvnFvBdwXS9xDN7ary71GGlnbLUTgHNv/6+mW3vHy1xbLzKU3vqfPGFgxt+l/VixfFayHgtpBllG4xR1/NerAJbLBaLxbJ5U1EB+0dC3l3sMFzymF2LEQLmmjk3ljp84805Ts/UBufQRC2k2c24uNgZFFvXlzo04pwDoyWuLXeYW0soNAS9qqWdmG6VAMqBSzcpiPICXwqyQqP7xZCGdqIIHEiKOwXSeusGhRknaq2RjqC7yUVeA44jubXa5fpyNGjm7Kans1NV018DfgKY01r/LeBjQLCL+9mKrTpvW9U397uNuaEQXxJCvC6EeL3org2eTAH4u9RubbWd2J/dqt5/SQESzVjVY7mdMt+MSZXm9lpMM0rp9kT0K+2EoZKP0NCKM8arIWNVn+NTdfaNlljq5OwdKnNyqorSmusrMUme40p4d77FKxeXODRSuuv32a0Acn4tvktob4X0FovFYnlY1gvkAeqhx43VmGrgshabaDrPcYxVRCtBKc1s887Zc2yiQqE0S+1kcJ6dm2txarrGickqca5xHWGaJ8Udb0wJuI4gV5puWiCATGlKvrmvki9xHEE1cFEIhBBUfEnoShxpmjGBJ/AdiScFZc9hvZFB/9zXwGor5vJiF43aNm5a3yMbcadlSNQz68qFEHVgAThyn5+5HzeB/es+3wfcfoDbAKC1/orW+nmt9fPl+iihK3AEOAN7+N2hgcDpResIcB3wXTMe7Adb18s+jjCz3izX5EXBaiclzhRFrwV5ZbnLickKjSjDc0wYte9Kmt2MTpKz1Eq4vNjm3GwLgImqT65MZ60WeoxVPK6uRnf9frsVQG7+nwEefhXYYrFYLJbNf/xP1wLWooyxskeUmi39rFDsqRvx+1jVZ6l9R/Y9Xg05NVNjtOIPzrNDY2XKvsOlpS4l18F1JJ4rcR1JuWfonWuIc02aFyiM47wJo1ZoNI4QlFxnYLVU9hym6yWm6wGjVZ966KK0oOS7nJiqMVwJKJQeNGqEMLWABFJlPLqUFgNri/W4joA8u/uw7n9/h8/l60KIYeCfAN8B2sC3dvxKbM23geM9G4lbwF8H/sam2/we8Cs9PdengLX76bXAzHoDzyHOc/Jc4wjT3cpUr1DSkN/vIpjuVR9HQaUkETh00wKJEQKudlITv9MzT+t3vzINFSFAK/70wiK+IzkwWuL8bJPT0zUuzLd54+YqQyUfRwrzhtQavyLxEb1xZML1lYg/KQfOewAAIABJREFUPje/pZh9NwLIF4+O8bXXbwJsENrvNCLIYrFYLJat2BwuvX+8wl95eorZVgKtFICDo2UcRxK4ktARLDQT/vX3bzFeDZiph0gp+YVP7OPqasT8Wkwjyri62EEIKNDGekEKKr5DoczHRaEHo0K4M/aKc41GI4Q2XTVXMlIOOD5ZRUppiiVPstY1puBDJZfRasB8K6HsO3SSgkKbpofEFF3mPs0ZH2emuOuXW1LCWMXnqlYPrtkSQgjgf9ZaN4DfEEL8IVDXWn//wV6W3pOidS6E+BXgG5jNyt/SWr8lhPivet//DYw+7GeAi0AX+Fs7vX7Zd+gmOVmvSyQETFQ9skIjBDTjnFzdbbm//S/c02ppE4/jStHbeDQVe6GNv1foO6DNi5HmBUqZ9uexiQqZglxrPnFghL+42TBhmHnCaMVnpGI2LYyg0OXCfEHoOr3UcvnQYvbdJq1bLBaLxXIvttpwB6MFLno5w8cmylyY73BhvoXnSk5MVLi+kvKJA8NEmWK5nbLayfj8R6d57fIKqlC8u9ji4nyH2WZM6ApGKh6ehDhTOJh/XSlIC2M+Ks3CHgrweo2VrFBIAa5QLLYSnpmpsW+0TD10mW3GLLRiXCn52y8d5ne/P8dalBH6ktB1kcKI7LOi56vV03+50hie9r27+pWVxLjlC+l6dz1JPYTefsR450ZCfEdr/YmHeVHeT4YOPKV/5L/+DW6sdECA5ziEnmSyVmKhGVHyXEAx10rRWu9KQC8ESG0qWUcK0ty8wP2WI0DgmtmwQLBnuETJdzgyXjUOtUmKFJJcKUqu5MpSl7RQaKVppTlJWlAJPbTWVHyHPcMlzhyfwHcdKqF7l1DeYrFYLJb3m/Ub7v1pydVl04k6OFqhGrr84OYqr11coRY6OI7ElSaF5cRkBc9zmWtERHlBnisa3ZyhskM3VRRKUw89ri13SAtF2XcYrfgIoBFlRGnBZC3kdiPa0CzxXRM6nSuzPZgrjRSCw2NlPvf0NF94du+GsOuVTsq1lS4rHbPwVg4kjW5OK87opmqD9htMQ2UrTbcjzBjxyld+OUkXrmypzdnpGPGsEOJHtNbf3vlL8cGhlGalkyGEma1KIWjFBWU/o1AMDM/2DAWsdnPa2/hvbIXuPftam4ic9c97fy6cFRA6mlxr0JqJqtklCDzJxcWEqu9RaM3NbkSSK1pRhpAwXgmQFUEzypFSUAk9zhyfYLwWorR+qFwri8VisVgeFVttuK90UhDwkb3GAiLKNccmq9RLHi8cMV2vf/H6DS7Md5gZDlntpmS5YqmdkBaKRmSCCQPPoRZ6hJ6D0sbIuxZ67Bku001z4jTn4wdG+J3v3aIVZ/g9XVaSKwqVM1Hz+cxRE3ITZQW+IzZkJDa6KedmmxRKMTVUYs9QCSkEhVJUPEleOCgFaaFA3ymytov2KzSoXIOU29ZUOy22fhz4shDiGtChv1Gp9Ud3+PPvLwKm6yGNSPQyDQWuLLixYqpgY7mvcYXk1HSN1W7KrUZEViiE7mmuei+M07OBWI8jwHUd8kKZqqv3IqS5aVsqpUmloB64TA2Fg8iCJFPEaYHWMFL2aMUZrjRTX1cKFPATJydY7mastGIWWinfvdGgHnqEjqCbK37tG+etIanFYrFYPlC2CqZOcyNM79OMMuoll1Z8J585KxSZUjTjHN9xaMU5Re+czgtFoSDNc9K8gxQCKaEoCppRxuFxycGRGqXAIck1P35ygn9/YWkg5wGj0dq3boM/dCWL7YQjE5VBJ24tzmjFOXmhmKhqqqHHzFCJq8sdjk7VSQszrvz+rTXirOiFVN8b8wi221PcebH1V3Z4u8cCVwj2jZVwGtCKMla6KZ1UGT8OTyKFoJ0UKA2jWvPFzxwmlPAvv3uLt+daxFkBGOHdVuaoCshyo56Xoqf6Mk0shBT4rjEtfXpPnSTXdLMcoU37U0rBWMUnzpVxoM8U3TRHA4fHyix3zQbHt6+sUA1c6qFptV5d6fLjJ8ffE0NS6yxvsVgslt2wVTC178oN2//1knGPX38bV0o8R9KOc6qBKbbWH7P9j6NUUQvNZEohODJe4fRMnUaU8YVn9wKmu5ZkBZeWugPzrP3DIZ3MZBQHnmQtygbxO/1OXN+fK/QcFnpbkQuthFurXQql2T9awnUcjk9W+f6ttcHvJjb9jlvwcMWW1vraTm73uFANXdBGJB+lCs+RuLLAEQJHSgJP0k0KumlOIzImay9fXERrzXjVZ6WTopJ8Sy1XX3znSuMUXyh1Z8sRcIQgcCVHJip8/MAIGuOcC4JPHxnl8mKb22sxt1YjYyMhJUMlF40R+J27tUahNI7QjFY8mnFOVBQcGisT5XoQ47PSSfj1b17g0FjloQok6yxvsVgslt2y1Yb7aMU3C2hRRjV0ma4FXF/ucmLS+Ei245xayeXEZIW351t0ktzooKXAlxJHCNbizDQuhCnMlDa2SpeWunxsf37X2bTcimnEpkFyaLRMUmhqSnFursVCM0Zp+OieOt+7scbhsRLnZnNurUa04hzPgU6z4MZKF88RplgEFtYiljsZaWFsJfrcq9i6n8PUPYstIcR3tdb3DH/eyW3eb4bLPkcmqgyXPJI857deu4rWRsBWCRyzyeAYm34BzLUSuqlpUzpSEnqmmm5EGaq3+ikFJg5Ag1IKRwqGQpckN+3Q/pujEri40uiuXrm4xBc+vo//8eeeHrw5/uHXf8Dr1xrEmWlLljwIXcF41efKcpc4UyitGCl5KA3PHRjmu9cb1EOXZmRasUvtmLdnWxRa88KRsYcqkL7+xi0uL7ZJC0U99Dg2UWG45FlneYvFYrFsy1Yb7l9+ydhvDiwgxiqcOTbO1dVocJsvvXiY1y6v8Omyz9vzLXKlmG3EhKGD4wjKnuzJf6CTFuwbCTk4WibONfE60dSVxTZfeeUyl3su9VrA+fk2gSvopAWukD2vSp9MQydK+cZck2OTVcYrHrONLu1UITHC+gSB70m01ix3cjpJjiPFYLol4Z7a7l7DZdum1/06W6eEEPeyeBDA0H2u8YEQuIJ/9848K70co0roIrUxJYuy3IjbtbH8f3u2SSvOSAtFxZF0kpy4Z/PgO1APfaKsICtMIqXjwEjomSieTBnzVMc41SJgtZvQiFKW2g7/9Ow1/uStWT52YISlVsrL75oO2nQ9YLWboZRGC8FyJ0VpwemZGmtxRpwqlloJf9HTbDWjjHrZtGIvLnZwpGCk5A86XbD76J0ri21eubDEeD1gKPSIc8Xr11Z57sAw7WQnTmQWi8Vi+bCyndfj5q+d2fT9faNlXr20bIoiyTqZj9ObQAlKnsR3HU5M1YmygmooNzQCXr20zHInZaTsE/ZMTkVZcG25w4GxMqHnkGaK0HOIsmKQlvL2XIuiUGSFHhRQSa4JzCWI0tx005RxGuj7bg3if7Z5LqQArVWxzbfvW2w9dZ/vA2x78Q+K1U5KlBZGKDfsglbcXI2Mx1ZidFgKKHlmFTXOcppRTqEhzlIExtZBChMNEOeFEe5pKDmAcOhkBVKC5wi056AxovzVTkpWgJSaqXpAnOZ8+1qXZpSzd7Tce7GMzUOhTcGnlCbKNE/vqTFWDfFdh+srXVxHcG25y186Ocn1lS4nepE+iz1/kPWB1dWek/xuePXSMuP1ANmLMSj13rDnZlucOTHxqF4Oi8VisVgGbC7SXjm/wFdevUKuFKErWemkLLZSDoyWibKCKM15emZkwzk3vxaTZAXDZX9wndCVdNOcJCtodFLivCBOTUj1ajch9ByacU7oO2hR4Dn08hYdBKbAmmsmVHwXraHsu/iOJCsS4nU5PlsVXerOt7bknsXWD5tWq0+SK/70/CJeTw9VKBiv+qx2s8ET5kqo9ErZRjcbtAr7onhVGOd5LYwxW985NsrBdxSB69KICwql8B2J6wjWooyVboYjoOK7lDyPdpIQupKbjYh6xSfwJK04J/AcDo2WyJSiEWXUSlDubS3WQo8Do2VuNSK6ac7+8Qpnjt9pxY5VAvbUw0FgJzxY9M78WszTMzW+c80IAANPgjb5VNtlLFosFovFshvut4R15uTkoNs1vxbjSPiLmw3iXFMNJU/PjDBeC2lG2eCcmxoKubDgkGSKvFAstBNWOwmNKCe+vYZCUPUdhko+3aygm6meG7zA6Qnm854TfV4oMmXOfKWhk5oRotCavFDkPTP09QHU/TQa1fscDTqN1tiGnW4j/hCiQWguLrWp+S5l32ig8iKh0EaM3oxyVjvG42qs4rLSzTfYPBTaFGVaCJOzJI0gPs6gUDmVwCUrjFu8LAR+KAcBmVkvVDNKczxXkOXGpO2G0mitaSUmQwo0riMZLvm8dXONODdjyWrgUQsdfvrjewdGpv1WbF/U3hchPmj0ztRQSCfOef7QMBcXO6xFGb4reen4uNVrWSwWi2XAg26tv3J+gd985QpKacaqPmlW8LXV6C6N8eZu12bT1GaU0Ygyxsouv/LPvsvbtxusdMw56kiJK424XgBZbgKnG5GCnrt8URQspgWBK1HKFFyJ1sie87zQprlS8iVJZsKmW4nZlDTmqBsf12ZLKAQIr1Tb7nl4Youtku9SFBqloJsVBFqz0EpMJSu0ifDRJm+JAtppgRRm7trzMUNjPg5dQQGDUSJAWmhIzfqo03Opr4Uuy53UXBtNoTRJrkhzzZ7hMscmKnz7yjJKm1lxocwQeMiTNLpGMyaFyVhc6SaAz/MHRu56bI8qeqe/TTJc8vjkodFB0fb53lqtxWKxWCy73VrvF2bv3G5y9soyYxWfmaESca54Z77FU1O1uzTGWxVzm885VeT8r//uErkqyJWxmmhEGaDIco0rIfAdOnHBamEW4VRbUQk8XMfBdYzdQyvOe0J8RW6k2IS+JHAdqqHLYjMe5CuCyRMEUx9sZQcF4EmBLrJk6+8+ocWW0pqh0GV2LUEKTTspiDJBrnrFVO/JWm/tEGfmi31rftn7pNCQK1NkrX+SNUZUJzHjNyk0raSgGjhEmUJj3Gh9KUiU5tl9Q2g0aa5IckU1dBmv+ewfLnNtpUNSFDw1XWMtzmlGZtyZFpqvfecm+0bLDxVCvR02L9FisVgs92Mrt/j+1zefF+sLs2ackRUFK92MeuhR7UllZpsxvuds+TObi7n+ZOfKYptf+u3X8R2BRpLmBXGhUco0TXxHICR3YnY0+I5DkiuOTYW4rYQoU0zWQoZLOe3E6K5NuoxD6DmmeOtmvbg9fUeb1UuNcSWk2ynkBQgh5Tbfva/1Q4vtLSW01rp+r5//oAhcSSPOjFmoBq0VWc6Gmet29DtaRnflmGKrl13oCsjW/Xw/DzFwJJ7rUPEdpusBC60EKQSeI5msu0zVAvaMlPj3F5YIPUklDAldh6LQJHlBnCkEmtFqQOA5RGlB2XPIlWalk97X1uFhTEkfRdFmsVgslieXrdzit1vKWl+YtZKcoZJPkhlNVTX0Bo7unzwytuXPwNbF3KuXlklyxVjFY60Rk+YKIQSa/tKbJs9N0aWEoOgVYUIKbjViaoHLZDUg05AXkloo+fGTe3h3oUVaKC4vdFjtpkRJQa7vNF4ExjbKc0wc33ZFRJxpcNxthdP3E8hvO398nIkzxcxQiT31Eo0o4/pym8VWuu3OZuCI3qqnuYnvmNVThMAVoJUgQ6MFuL0xY7/idXoVV8V36KQ5zdjMkJ87MMLH9g/jOyZAGuCnTk9zbrbJSiehGed04pxGnHFgrMRyOzNvyFaC75ji2JGC8WpwT9+rBzEltY7xFovFYtkpW7nFb7eUtb4wq4ceRaGYTRPaSY5GDxzd1y9h7aSYm1+LGSl7RJkiVwohQSs9kP2oXh2ktEYIgesYR4HAkcRZwc88PcWlpS4l3wVtzvNWkrN/pMzr11bJlSZKi0H+oQZ81zgSqELjetBJ9Lb2DxoQQmzb2dr2Gz/MOFKQ5orAdzg9UyPJNY4jcOT2D7gauFQDB98xL1KWK5K8YLjsMVEP8RwTzVPyXDzHJCBpjPPtsYkKviNoxTmZ0hweK+NKyWsXl7m20uHFo2PMr8VUQ5djExWkEMwMlXh6b51q4HJgtMLRiQqr3ZRWnKLQdLOCim/sHaqhO/AI2cz6vwj6nlv94mwr+sVZJ86Z7gnkv/b6Ta4sth/Nk2+xWCyWJ4oXj47RiDKakWkm9MXqW22t9wszgGMTFYQQjFY8yp7DQium0Jovnzm84Q/89T/TZ3MxNzUUcnq6RpQWKKXJckWUGeepknunAMoKKAojmv/o3iGOT9UYqwRMDZV57sAwWmsW2wkztQCtYc9wmWpgxo2bU2O0NnWDApLMFFoO26O13tb39InUbNVCl586PQ3A2SvLxpQs7b2Qm8pSV0ChNaXAxettG841YyNcL3kcGDE+H4vtmKIXdKn1HdGcBhbaCZ2kYGYo5EePj7PczWnFGUMlj8lawOGJ6uDNNF4Lef7gCBcXOyy1E8YqAV86Y1x3f/eNW/ybH8yRZAWHxyo8u98kp798YZEkV3z17LXBm7vfmfrBrQYf3z+84S+Oe3lu7Wb2brFYLBbLbvS962N8RqsBT03VODfX4uRUjRMz9S0nKVtF/2zesH/x6Bi3ViM+c2SUVy4us9yJkb1zGgEyy8lzTa7BdQX1wEhx0kLzt186TKyMifmpmRqCGt+70aAVpXz3+irXljoU+u5+VdoTeEtgqOwR5wVZpigydVd3q2cr8dBB1D9UFOtmqs0oI3DNWmia38kjF0DZF4xWA5ZaCY1uxmjFoxa6dFOPSugyUw8ZqQb8yESFhXZMM8qRQhhzVKWRUiMxK6RSwGjFZ6QScHzamOorrQdFz+Y34GnXoRFlG8Z9f+dzJ/n8s3sHY8E4yzl7eQUNfOboKJ045zdfvowQcHC0wnTPZ+Ts5RU+c3Rs4Lt1L8+t3czeLRaLxWKBnet7Nxdm+8cr/OKnDtzzZ3dSzK2/jec6vPzuIgdHS4zVQt682aChoVZzSTJFLfTopAWNbsoLR8a4tmqu+alDI7x2eYXhksdKO+byUhelTVF1Lzm3EDBdD5hrmsZK33l+fT0hhBGyb3eNJ7bY6ntQ5XnBQishTtWGTUQjhBeEjsNwyacVG80U2vhnffrwKM8eGB1cc/9omVurMSXfAQ2L7YROkhN6DqFvthiqgcfFxc6WRc9O/zJYf7s/v7xEvezxzJ4641VznZVOCgI+std0vZ7ZU+e1i8u8ebvJSyeC+3pu7Wb2brFYLBbLbnmQxaud/Mz62/zjb57n3FzLaMCk5NkDw6x1M243Y+LcyHAKZWyX+t2yr7x6hVNTNeolj7lWSqHAcyRSqLt9s9ZR8R2uLHfxHEnVc2jpnPUDQ0fS8+H6kI0RRyo+ldDl/GzT2DBojecIMnXH3sE1Zuksd1KGSx4TtYCZ4TJRmjNdC5htJhxZZxp6YLTCUOjzzlwTV0rQMFENODFdo+K7zDUjip6RaT/dfHPRs5u/DA5PVAddKLmuM5nmivU1+Hg15DNHR/nu9bUd2TfspF1rsVgsFssHwU4XuD7/7F7i3hTorXCNq4sdri53GK34DJU8bjVMRN98s0u9NES95JnQ62bMoYkqSaZQWhFl+p6FlgRyDTpXhK5DrIx/Z8kTaKV7wnhJXijgQyaQD1zJF184iBAmNTzOFXFhxG1+TyjvOCbdu5vmJsYnVxSFouS7RLnqabVyvnlujj+/vMxULeAXf2Q/p/cMUQ4kvivRGhZaCcudhKl6iVMzNUYrPnNrMZXQvedG4E7YSjTou5LA2yjRC1yXnzw9xd/96ZN88YWDO2rXrn9sobvtmNlisVgslveF3Sxw9c+ySugyHHrMNmOGQpfhsk+uTKB0Xij++es3+X++dY1359YYrwYsd1KWWsY6Ii3MJGzzCeg7gpIr8Xvi+CgtiHpuAZ3U+HO5QpIW9MxVBcNmWvThGiMmueIffP0H/Mvv3qISOLjS2DUU2ri8Sox9f5praqHDaMXDlZJrK10OjJSYb5pu19uzTcbrAU/P1Ahcl9cur/Cjx8b4rT9r4zmw1k3pZjkLjuSlYyFSSn71cycAIzj/V9+99cDWClcW2yy3Yl5+d4nxasDpmRqh5zJa8RGCh47qSXLNp4+MD65xP7sIi8VisVgeJZu7WEuteFcLXOunRd+73kAIE58XZwWFNuM9lCbOFN98e5GP7KkRZwV/dmnZuJ5iJl2OY2qC/jaiIwWBK+km5gv9Cqr/b66g1UudUdrYTWWyH9a3NU9ksbXaSfnW1VUqvoPo2ZJ5UpMp86S4vflqyZMcn6zSjAsE4DmCdxdbKCWo+i5lX3Jzpcs7t5scm6pybKLK7781h+dIKoGHEMa/QynNYifhb//4MYBtfa/gzhahK+/EAW0uyNZ7Z/3oiXHemm3x8oUlzpwY58svHRlc50Fd3+1GosVisVg+SLbyiHzlwhIvnRjf8Xb9eg5NVGjHOcNln29dXabiu8R5gec6lH2HolCcX+hw5tgYs60ER0gCR/VMUU0aTMWXhK7DR/cP8Rc311CAL0BLQb6FiL6XuNcLsO45qG/DE1lsRWnBfK+l2M0UniMphKZQCgXUQjO7LXsOUgim6wGNKKMdFzTjnJ8+PcXV5S7L7ZTAdfBdwbnbTd6db7EW5Xx0b53xWnXd/eXEueLwRJWvnr22ZSHz9TdukeSa4ZKH58CfXVpBAC8cMVuGX3nlMpO1gELB1eUOe4bC3s96/FitRDPKqITuoBh6kE7Zg9pFWCwWi8WyHQ9ilN3/oz8tCr51tUkzyuhmBS+/u8i+0QrNKKNe8ii5gk6q+LVvnL/ntX/hE/v49T9+FzA1gCvAlZLxik+UFcZcXEOt5HN8us5sI6YVZWRKk+YFQgimqj5aCJ7eM4zScGG+RT10ubLUxZF3cpPXowHXcZiqB1xRKr/rF+vxRGq2XMe0ANfinGpvUzDOFUioBA6n99SZrJU4MFahVvIYqQSMVwOOT9U4PlnlI/tG6GYFUggKrVlup7TinHaSE2U5l5e7tJNscH/abCEADMxL11MNXd640RgUYZeWuoyUfYbKPpeXusw3u/z5pWV+743bXFvucGu1y9uzLZba8YZrbGdsej82z8EDz9hFLLXuXM9uJFosFotltzyoUfb8WkyS57x+tUGaKYZKHmVP8s5cm7lGRD10mWtE/On5JcqeZHoo5MZyh3/we2/x3/9/P+CrZ69tuI99o2U+trfOzUZEkhkn+Gf21Pno/hFOzwwxVgk4OF7BlfDyu4sINI0ow3ME9dBFa81CO+W5A8N86tAItdAjzTUrnQwhjKn5+kLLEXf+zZWi7DsgxLaep09ksaU0HB6rUPSyBbtpARoEAt8RzDVjPrq3anIIs4JPHh7lU4fHODJR5fR0nZffXWSxGbPaTZlrRMS5QghwpSBwJK0o48pSG601UVbQ7GYDA9LtnHBBDIqwZpQReJLQlcytRbxycRnPEXiOIC00rTgnyQsuLnY2XMOV8NWz1/i1b5y/6412Lza7zD+zp44G3rzdvK8bsMVisVgs27HbFJM+U0Mhb822BiHQAkGUKyaqAVFW0OrpkUNP8L0ba/y7d+b53o0GjhA042xDUdcv+GaGy/wXnznEz35kBs+R5EqjtKLRTWlEGT96bIy5ZkKzmzE9XOLIeJlGJ+XWWkw1cPnLT08Rug6//sfvUvEl0/WAuOdSv1mOVWhTaJliS3NrLQYpP2zFlmao7HN4tEQ7zUlzhSOh7AlqocdMLSTXklM9G/+3bzd58/Yas6td3rzdZL4RMVIJKHmSdpqjlKbiuwyVfPaOlil5DrONmO/dWOXacoexqs8Xnt0LbB9r8PH9Q4MirF7ySDJF3IsbUErjOZLQdyl5DtNDIcudlMVWPLjGtZWOMVR7gJidzd22vl1EkqtHtjlpsVgslg8f201z7jeJefHoGEvNBNVLZYmygnaSc3SizHg14OP7h4nSgopvrj23FrPUSnCkcYJfX9RtLviePTjKjx0fZ7Wb8saNBjcbXT66t87NtZhDYxU+c3SM0HUYKvtMD5f4xIFh/suXjnJyZoi5VsJwySPKNT/9zDRP7xmi5LogxaBgcsWdzlaqTDxQyRXoPI22e7xPpGarEjj4nmSpm+FKyd7hEM+VzDZiFlsJrThjuOxTPjzGL3xiH69dXqESOPzZpWWacQoIRsoeHUzSd8lzqIUua1GO7wrivMB3JUfGq/iupBbe0T5tZ156c6XLV169Qq4UJVey1M4o+w4lT5IXkigtODpRAmCsEtBOcsYqweAak7WAkuc+kKh9KyPTvl3EF184+IiffYvFYrF8WHhQo+zDE1XOnBjn7dkWzTijFnocn6yilDnfLi52qAYehVKUfIc4Kyj5DleWOlQClz86N0ctcKmHHiMV/65kFPMzin2jJcarAaHnDAT4k/XSwHz8G2/N9hbpDM04o15yaUYZ49WQ/+jZvfzsxzTfvrLMaxeXaUQpUgiU0khHoHONIyEvPoQCeVdKTk/X+d61VfYMhYSeS6EUjiMotKCbFEzWHYSA16+vopTinbk23bRguBwQpwWdRPETpyZ57eISN1YjVrsZI2WXdlIgENRDn+cODjNeDWlG2YaiZ7N56ZXFNq9dXuHUVI3ZZsxyJ8V3BSemqlxb6eI7kjjXuI4cpKJXA49f/dyJwXV+7Rvnt/zrYSeidmtkarFYLJb3goc5X77w7F6S3GwkVkOX68sdvnV1lROTVS4stBmruFxdjpiqB2itWYtSVjsZe4dKDIUea73J0XjN31DwLbVi/v27S0ih6SQF88013p5tUQ9czs22mKyXWGrFXFzscH2lS+A6LLVixmsh9dCjGWVIabKVm1FGnhufrX0jJcq+ZKWbkReaaugSpTm+IymHLkK63naP9YkstrJCUQldpusBhVLcaMR0kgJHQuAIPMflZz4yje84JhKn5FFw1oaKAAAgAElEQVT2HRPvU2hKnkM7yVnuZvzYiUn+4NwcoStJc0UzznCFwJOCN240+MlT0/cteta3OA/1iqf+duHf+uxhvvb6TVTP2XahFeNKuW0q+oPE7OwmRNRisVgslp3yMOfL4YkqR8dK/PbZ6yy1E8arAT/z9CS1Ssj5hRbl0ONzp6ssd4wc6NZqwWjFY7QW9Bzg4dR0DQE0IrO0Vg1d3rzdpJMUlHxJUWhqgUeUFiy1TZTP1cU2b8+3cISg7EkWWgn/96tXGC57jJV9GrGZYgki4kyxFmdMVH329qJ/4rxFoTHmplLiOQ7VwLunQP6JLLaGyz5gzM0uL3VxBGS5IkOTCsGpmRrj1dD4YiBYaidM1kImqwHXVrpkBZQDh8VWzJ6hEs/uHaKd5lxa6FAPPcqeQ+BL3l1o8+y+GN917ln09GN3+pV0M84G7c8vvnBw8Eb1XYdPHh574FT0e/EgWVUWi8VisdyPBz1fXjm/wD/99k3GqwFHJio0o5x//YN5fuTQCPuGy1xf7eI5Dp88XKcd57SSnFrgcO52E4D9oyXKvkOu2FDwJbmi7AuitCDOFL4jqQQOCMFUPeDPLy9xaalDnquedYNEYxJhcqVxBMQZRkJUcskKRZJpVjoJ7aSgUCbyL0kV1dDFlea+0LrY7rE+kcXWaielE+cINGnPSRborW/e8cloxzkf3z/Ed643WOtmDJW8Xi5ijCMFY5Vg8AK+cmGRU3uGKArFhYUWS62EAvj9N+d4dv/wwGx0K6aGQq4vdXhnvkXJdze0P68stnccwGm7UxaLxWJ5Uvjad8wIsd8gcR1BO8n41tVV/rMXDlLxHV6/tsobNxqUfJdq4BDniiMTVcKepdPZyyt86sjohnP0f/nmed6eXaPQGk8KCqWZbybUQoduknNxsUvZc0gQdNKCTCmGSy6uI/EdyWonZWa4xMnpeu83bbPWSbm20uXpPUM4UnC7l714cjjEcx3m1mK0Vh8uzZYjBWlRcKMRG5FcriiKO9b8b882+bffv8WBsQpfOnOETxwY4TdfucJiO2Gs4jMzHDLbTBgKzabDoZESv9NOmKgGKDRKGVv/WuiQZAXiPtGCLx4d4++/OYcrBaEnN7Q/d+PabrtTFovFYnlSmF2L2DtcGny+0EqoBi7dNEcKY5fkO5J62eOl4xP84Q9mubEaUw89Qi9AYJon64/gK4ttvn+zQZIrnJ6QXaGRAiaqAddXYso9/812miOEcXVYi3KGSh6OhKQoiHPTpGrHGUleMN+KkULgu5LD41XqJY/lVsKtnnvBWMUHVWRswxNbbF1c7CClJM3Ni6YkoDW51jhCcH6uxaX5NhcXWnz26DhfPnOYq6sRF2abXFnucnq6xoFxY///2uUVPrKnzmwr4cZKRK3k8dRMCVdKfE9ycLRyz6Lp8ESVQ2Nl1uKMtZ4r7uk9NUZ724YWi8VisXzYqIUub882EUIQeg7NvsloT5t8cbFDveyR5gopBK4rOTpephHlOI6kXvL4zNFR+lZYfb+tOFecmq5xYb5DUiimaiGTNZ+VbkYzTvEdQSsuKJRG9UZdGmNOutROcaVEa1huxcw2Y4QQyF5X5dtXV5kZCvns0VG0rvEHb84zVgmol9x7+mw9kcVWoYw3Vc13aHTSgR+G7Nntu44mU5qR3gbDn19eYb6V8KUzZhS4d6R8l8VC4ApKgUc7yhmv+aS5ppsWnN5T2yCQ3y624ORMnU6ck/bMSr97vYHvSE7N1N7/J8hisVgslg+QK4ttKr7RH1dDlywvWIsyXCn42Y8Yg+1mnBH0iiqAeugR55K0MGdzM8p4a7bF6WlzjvaX0SZqIWmmeO6gx61GZGwZMBqsauiRF4qSI+gkmyRWGrppzlDocXqmzuWlDoXWCK2RQjBe9Sn5Lgi4tNilUIp9wyWWOymXl9oI4WxbUz2RpqaF0viuRErjyq6ECaM21SloJCXPpey7KA0jZZ/lTjookrayWCh6AjzfE/zg5hqXFtt4vRr2+lKHq8sd/off+QF//3ff4sZS5y7j0RePjnF1ucOfXVomzgsCR7IWZSy0kh07wVssFovF8iTw6qVlPrJ3mL/8zBSh79BJC0bKHhNVn6mhMkprfEfSiDKOTVQAODZRYX4tYqEVk6TG77LZzZhrmnO0f34fm6jQTQtcKTk+WWWiFiCk5BMHRvjskRHS3Eh51hdAvmvi+aQQKDQnpqrUQocDo2UqoctT0zU811hGdZOc2bWYN26ssdhJ8B0YLvmP3zaiEGIU+BpwCLgK/Cda69UtbncVaAEFkGutn9/J9UcqPqena3z/RoPxqk8nVSS5SZAMPEgKzVDokhWawHMIPEmjmw66UdtZLNxc6dKMcpJcEXpGr/XH5+bpJAXTQwHvzDYRApLMVOp907RXLy3zxRcOEjiCRpQxuxZTL3k8f3CYqXppV7oti8VisVgeV3YaSt3f0q+Xhjg+NQSY9Jfzs00qvWnRqZkaC60E33FM8eU6OFKyb7hEWihqocczR+v4rsOrl5YH5/d4NeT5Q8NcXOyw2IoZqwQMhR4HxisD+6WzVxogwAEmqj6p0sRpQeg7VAOf0HMJXJeTUzUuLLQZCj2G05wri23WopyxikBpxVI7YXYtJnAECLFtA+uDGiP+PeBPtNb/SAjx93qf/7fb3PbHtdZLu7l44Er+zudOAnBurkWjk9GIUlqxMSLzHY0CVKHYM1QmyRSB5wzeGFtZLDw9U+Mrr16hHpqIgNm1hOsrEVprCqAW+iy2UlxHsNROeePmHQ+u87NN/vE3z/ONc/OMlH2OTlbwHYdLi12Gyx7tNRPj8yDJ6RaLxWKxPA70NVPDJY/pXuHztddvbhkHt1Vj4/pyh9Uo23AGAhu28E/P1Hlqpj7QUIEp0ubWYn7+ub2D83u0EnDacWgMlQab/P37e/7wBM8fnuDffv8WV5c7+I4DaYEbCrLcxAelRcGeesCfvLOAUqYxM171EUKyr7dEVyij80JDigax/brcBzVG/Dzw272Pfxv4wntyJ8/u5ch4lb/01CT/+acP8XMf3WOcaIVmoRVTKEWc5ax2U8Yq/qC4+YXn9w0q635u4OvXV7m52uV2I2Kpk7JnOOTUniHiXBG4RtwX+g4gKHkON1dMRNL15Q5Xlrucm2sxUvbIlebGSkRWmJTwt2ZbTA2FD5ycbrFYLBbL48BuQqk35whfXWrzraur7KmHG85AgC++cJC/+9Mn+eILBzk5Ux/kDPfpT5+2O78PT1QH9/fGtRX+2X+4xv/5pxe5sRLhYNJb0sKI8Gsll4NjZd640WC2lVALPA6Nl2lGGTcbEY7Q3G7GKK2MXAmzDWmE9o9fsTWltZ4F6P07uc3tNPBHQojvCCG+tNs72fzEl32HZ/YO8/PP7ee5A8NoDddWIp6aqvKlM0c2xO2sf3EBXrmwhASaccqt1Yjv32qw1k1IC0XomQbhZDWgk2TMNWNurnT5Jy9f4g/emmPPUECaK/asW3FdbMUoNEvNhBePjj1wcrrFYrFYLI8Duwml3nw+316L+eTBEQ5NVO95Bm4u0po9z8oXj47dczrUd6v/s8srNOOU0YrpvnUyzUjJZbwaMFLxOTFZY7QamEJMCA6Ml/nc6Rn++if3c2KqxnI3QyuYqocEroMjxcBJ3phIbM17VmwJIf5YCPHmFv99fheX+azW+jngrwC/LIR46R739yUhxOtCiNcXFxe3vM3FxQ5aK+ZaCRrBR/YP85mjoyx3M/7Vd2/x1bPXtuwkvXppmbJv8guzgp6NP8ak1JP4jiDKCjSaOCvopDmeKwk8hyQtuLTYwRHgOQ4HR8uEnhHHo+HMiXEOT1QfODndYrFYLJbHgf5ocD33ipVb39g4NFbhwHhlw/e3OgO3614B95wOXVls89tnryMFDJcDpmohe4bLjFZ8mknBz350holqgONI4rQgV4pC6YE4f7wa8tLxCSZrPqNVH6WMZKnkOWhAKQV6+2LrPdNsaa1/crvvCSHmhRAzWutZIcQMsLDNNW73/l0QQvwO8Eng5W1u+xXgKwDPP/+8hrvnx69dXGShlXBkvMJoNWClnfCdpQ6TtYAXjoxtO1+eX4vxXFNUOcJhLc5IsgIQPLd/iMl6mZVOyvWVCKVhvBrw1HSdWuhxfs4kmg+VPbppQdl3mBkuMTNU4shElS88uxfYPvvQkfDVs9esjstisVgsjzUPEyu3m/zfrQy+v3r22mA6BMYaYqWT8OvfvMBQ6HF9tctiK2GyHpAXmmsrXQ6Oltk/UuLcbAvfdXjuwDDnZlsstRMOjZU5Pmmi/db/PkcmapQ9yXevr+FI4/015vkURgi+ranpBzVG/D3gb/Y+/pvA726+gRCiIoSo9T8Gfgp4czd3snk0l2ttLBfiHIFgLc4JHEneW/fcrm05NRQS54rpekiUFfiOZLjss2c4pNEt+KvPTHPmxARDJY+hssfJ6Rq10LzgM0MBcaaIsoJPHBwajA5PzdQ2FHVbtUavLndYaCVWx2WxWCyWx557aabux3Zn4HIr5te+cX7byVOfzdOhpXbM27MtVjopzTjDEcbaoRVl+I7EcyQL7QSlBZO1gH/z/dv882/f4HYj4hef38c/+NmnkVJu+H2urXQYLXtcXOgwXfc5Nllh33CZaujx2aOjFHF7W93PB7WN+I+A/1cI8UvAdeA/BhBC7AH+L631zwBTwO/0xP0u8M+01n+4mzvpr5b2KbkOkSzoxBlaa9pxhpSCknvHGmO9QWmfF4+O8UdvzTHXihivBkghSAvFVD2gHLhcXY3uaLveXUSuCw/wXIcDoyVGKgFZAWeOTwy6U1cW2xu6Vp89MsrV1WiwdTFdDwg99y6DVWsVYbFYLJbHkQeNlduc/+tIIzcPPXfQJdtusxHu7oxdXOzgSMFIyTfTpZLHwdESlxa7+O4dy6d2nIPWHBirDMKw/+m3bzK9botxbi3GlUaR1R899jtgP3ZynM8/u5fDE1X+t/90dWsNEx9QsaW1XgZ+Youv3wZ+pvfxZeBjD3M/m5/86eESaIjygmacUQlcSq5jvt5jq7bl4YkqX3rxMP/d19/Ek5pqyWWqEiCF5OmZ2mCm/OLRMc7dXuPyYofhkocW0OxmHBqv8OWXjmx4g2y1Ivva5ZUNb6Rf+8b5LXVcNuLHYrFYLOt50qyDLi122DMU7rjZsHmEudiKcaXk2ESFi4sd4lyxb7RCmhs305V2ykjZp+xLyv6dMOz+v1/7zk3+j7/x3OC+vnr22qD5US95TNZLNKOMSuju6Hl+Ih3k+2xuS07XAlKl+fThMX7i1BSfPjxGqszXN281bObMyUl+/rm9HJs2mYaj1YDnDw0TuO6gODs8YbYan5qqcmmpwzuzTUJf8lefmb7rxeiPONOi4FtXVzh7ZZnLS21+941bg9vsVmxosVgslg8fj5t1UH9qs5Px3/qfWf8YVtopb8+2WGrfaS7ca2ls8whzrBJwaqrGeC3k2ESFKM1Z62YcmqjwV5/Zw899bC//+y9+nEJjcg3XUS9tjOD76tlrfP17N3lrdm3Hv89mnshsxPUEruDPLy8Bgo/vH+JXf/L4YFS3f7zCmePjG0Z3P/301LZV6uef3Uvc60bdS/xXCjz+2nP7Brd57fIK+0bLd4vuHfjOtTXKvsNQySNOC15+d2nQknwYsaHFYrFYPhys1yfDBys52Y2x6Xo2P4bxWkAzyri42BmI1O/XbFg/wnzl/AK/+coV3rzdZKzqM10LuN1MGA49KqE7OOtnhko0o3zQ0QJoRqbgW/9Yhkou78w2+YvrDY5NVXl2/zC+4+y4+fHEFlvrn6TPnZ4eFCr7RsucObnR1uvMDq+5eaa8VXG20zf91FDIK+8uUvYdwn7IohCMV4PBbXdyfxaLxWL5cLNZnwwfnOTkQQu/zY/h2ESFb19bNZ6UWu+q2XBlsc1rl1c4PVNjthmz1E5wpeTLZw7fdf7/wif28et//G7vd3VpRuZ+fumzhzZMoJrdglxB4ElmGxGdpODIRIUvnTmyo+fliS223qtK/37iv52+6V88OsbvfOcm4/UArTVxrojSnOcODG9oSz6o2NBisVgsHw52Y5vwXvOghd/mxzBeCzk1VeN2M951s2H9+X9o3Ny+GWVcXY02NFeuLLa5uhpxcKTEpSWTo3h4osovffYQZ05O8h++cZ7poZBvXW0yXgsYKnsstmLWooz9ox6TtWDH5/MTW2zN97YHzs02acYZ9dDjyHiZdpLf/4cfgp2+6Q9PVDlzYpy3Z40PVy30eHpmBN91qIRP7MtisVgslkfM4yQ5edDCb6vHIB3Jr37uxK4bDjsp+NZPv37kyBin9gzRiLIN487+Y2lGZpsx9BxcR3JoXPLJw6O76hw+sQJ5R8LZyyskuWIo9EhyxdnLKzi7fMS7FfrdK0pgM194di9HJqp86vAYnzw8iu86297WYrFYLJateBh/q0fNbs7A9TzKx9AvkpZaMWcvL/NH5+Z4+cLihvN/JxF5/cfiu8ZVPsoKojTn2ERl153DJ7aFIjDBiqL3yYbPe9xvVfZBhX6bRfnb3d5qsiwWi8XyKHivJSc7tZZ4mHPtUT2GF4+O8ZsvX+bqUod62SNwJI0oY6GVcGWxPYjIu1/3q/9YfveNW7z87hLj1YDnDgwPGiO76Rw+scVWruAzR0e5tNRlLcqolzw+c3SUrDDf30khtVvd13ai/HthNVkWi8VieZzZbePhgz7XDk9Uma4HrHRT0lxRL3l8dm8d33EG5/duJD9/53Mn+fyzewfF5vptxp3yxBZbUz2vkRcO32ldNqOM0ap5yIMtg7zgW1eMrst3JF9/4xb/zedOArsX+j1O67cWi8VisTwKfhjPtlzBS8cnkOLOPEtpPTi/d6tze9gC8onVbN1vbjy/FhNnOa9fWx3ouhDwyoWlgS5rt6aim7OZYHemZxaLxWKxPG78MJ5t9zu/32+d2xPb2brf3NiR8IdvzZPmimroMlEN8BzJeP2Oz9VuK9/Haf3WYrFYLJZHwQ/L2bZeV+ZKmGsmHBqr3HV+b9af/fxze9/zDt0TW2zB9m2/K4ttFlqJCacMXfJCcXmpw1jZ58yJsUG1vluh3+O0fmuxWCwWy6PgvT7bHkWu41a6MiESoiynneSD8xt4oMW3h+WJLra249VLyxwcrXBqOmZ2LSZXmtCT1MsugesOdF2wuzmt3S60WCwWy5PGe3m2PejW/2a20pUdHK1QCV2++MLBwe2+evbaB6I/+1AWW33h+8f2D5MVDcq+g+8KltrpQ1frH/QWhsVisVgsj5r36mx7VOL7nS60fVDRRh/KYqs/fx6vhjx/aJiLi8amf6wSfCBGcI+ihfpeXMtisVgslveSR1X87FRX9iD6s0dxrj6x24j3Yv2m4mgl4PR0nY/tG3mgWICHpd9C7cQmZbzTa6Hez6n+vb6WxWKxWCzvNbvd+t+OnTrX79bh/lGdqx/KYutxijbYSWTAB3Eti8VisVjeax403mczOz3Xd3v+P6pz9UM5RoTHR1v1KOfHH9Qs2mKxWCyWB+FRiu93eq7v5vx/VOfqh7bYelx4lP4lPyxeKBaLxWKx9Hlcmh9b8ajO1Q/lGPFx4lG1UB/1tSwWi8Vi+bDzqM5VobV+j37FD47Dpz6qf/kf/4sfmm08u41osVgsFsvjyU7PVSHEd7TWz291jSey2Dp48iP6f/rtfzNwuf2gxO8Wi8VisVg+HNyr2HoiNVuOFIOtAXi8k8ktFovFYrG8/7yfk6AnstjKC8XZy8s044xa4FIPvfv/kMVisVgslg8FjyomaKc8kcVWJy1IcsVQ6LHWE7NdWWzb7pbFYrFYLI8570fH6VHFBO2UJ3IbUQChJ0lyhdJwarpmjT0tFovFYnnMeb+SUObXYqrhxn5TNXSZf498KZ/IYktKwVqU4XuS5w8Nc2Cs8p49gRaLxWKxWB4N71cSyqOKCdopT+QYseQ5/NTp6cHnzSh7ZE+gtVawWCwWi+W94f1KQnnx6P/f3t2GaFbWcRz//lrXB7asyE1dXXsAK3Mpo2FRzDBaSyQ0BcGIqHyx+EIqQigSygqhCKTIICWDih7BLEtLkzQNMne0NVfXDbHCaaVdMzUzsNV/L+ZMDDWjO3vf15y9z3w/MOx9zn3mnP9cDPf89rrOua6X8f3pmf+ef272gncef/hYrzNnkD1bzzxbTSb2dKFnSZLaWY4ep7lOk3/8699s2/k4Ox5+ovkayYMMWy9dc2CTRaZd6FmSpHZar4Qyv9PkuHWHsmHdi1lz8Ormo1SDHEY86IAX8L4TXzH287rQsyRJ7YxzYeqFLPdTiHMGGbZacaFnSZLaarkwdV+dJoMcRmzFhZ4lSZpcy/0U4hzD1hLMdW+2uB9MkiS11VenicOIS9Sye1OSJLXT+p6wxfQStpKcC1wCHAdsrKrpRY47HfgSsAr4WlV9btmKlCRJg9NHp0lfPVvbgHOAKxY7IMkq4CvAacAMsCXJtVV13/KUOH5OiCpJ0srTyz1bVbW9qnY8z2EbgQeq6sGqehr4HnBW++racEJUSZJWpv35BvmjgIfmbc90+xaUZHOS6STTu3fvbl7cUjkhqiRJK1OzsJXkpiTbFvja296pLLCvFju4qq6sqqmqmlq7du2+Fd3Qcq8wLkmS9g/N7tmqqk0jnmIGWD9v+2hg54jn7I0TokqStDLtz8OIW4Bjk7wqyYHAecC1Pde0z5wQVZKklamXsJXk7CQzwEnAdUlu6PavS3I9QFXtAS4EbgC2Az+oqnv7qHccnBBVkqSVKVWL3gY1saampmp6esGpuyRJksYuyZ1VNbXQe/vzMKIkSdLEM2xJkiQ1ZNiSJElqyLAlSZLUkGFLkiSpIcOWJElSQ4YtSZKkhgxbkiRJDRm2JEmSGjJsSZIkNWTYkiRJasiwJUmS1JBhS5IkqSHDliRJUkOGLUmSpIYMW5IkSQ0ZtiRJkhoybEmSJDVk2JIkSWrIsCVJktSQYUuSJKkhw5YkSVJDhi1JkqSGDFuSJEkNGbYkSZIaMmxJkiQ1ZNiSJElqyLAlSZLUkGFLkiSpIcOWJElSQ4YtSZKkhgxbkiRJDRm2JEmSGjJsSZIkNdRL2EpybpJ7kzybZOo5jvtTknuSbE0yvZw1SpIkjcMBPV13G3AOcMVeHPu2qnqkcT2SJElN9BK2qmo7QJI+Li9JkrRs9vd7tgq4McmdSTb3XYwkSdJSNevZSnITcMQCb11cVT/ey9OcXFU7k7wc+EWS+6vq1kWutxnYDHDMMcfsU82SJEnj1ixsVdWmMZxjZ/fvriTXABuBBcNWVV0JXAkwNTVVo15bkiRpHPbbYcQka5K8aO418A5mb6yXJEmaGH1N/XB2khngJOC6JDd0+9club477HDg10nuBu4Arquqn/dRryRJ0r7q62nEa4BrFti/Ezije/0g8MZlLk2SJGms9tthREmSpCFI1fDuJU+yG/jzEr7lMMCJU8fPdm3Htm3Ddm3Htm3Ddm1nqW37iqpau9AbgwxbS5VkuqoWXTZI+8Z2bce2bcN2bce2bcN2bWecbeswoiRJUkOGLUmSpIYMW7Ou7LuAgbJd27Ft27Bd27Ft27Bd2xlb23rPliRJUkP2bEmSJDVk2Ook+WyS3yfZmuTGJOv6rmkIknwhyf1d216T5CV91zQUSc5Ncm+SZ5P4NNKIkpyeZEeSB5J8vO96hiLJ15PsSuJya2OUZH2Sm5Ns7z4HPtx3TUOQ5OAkdyS5u2vXT4/lvA4jzkpyaFU90b3+EPD6qrqg57ImXpJ3AL+sqj1JPg9QVR/ruaxBSHIc8CxwBXBRVU33XNLESrIK+ANwGjADbAHeU1X39VrYACR5K/Ak8M2q2tB3PUOR5EjgyKq6q1tH+E7g3f7OjiZJgDVV9WSS1cCvgQ9X1e2jnNeerc5c0OqsAUyhY1BVN1bVnm7zduDoPusZkqraXlU7+q5jIDYCD1TVg1X1NPA94KyeaxqEqroVeLTvOoamqh6uqru61/8AtgNH9VvV5KtZT3abq7uvkfOAYWueJJcmeQh4L/DJvusZoPOBn/VdhLSAo4CH5m3P4B8uTYgkrwTeBPy230qGIcmqJFuBXcAvqmrkdl1RYSvJTUm2LfB1FkBVXVxV64FvAxf2W+3keL527Y65GNjDbNtqL+1N22osssA+e7e130vyQuBq4CP/M0KjfVRVz1TVCcyOxGxMMvLw9wGjlzU5qmrTXh76HeA64FMNyxmM52vXJO8H3gW8vbxJcEmW8Dur0cwA6+dtHw3s7KkWaa909xRdDXy7qn7Ydz1DU1WPJbkFOB0Y6QGPFdWz9VySHDtv80zg/r5qGZIkpwMfA86sqqf6rkdaxBbg2CSvSnIgcB5wbc81SYvqbuS+CtheVZf1Xc9QJFk799R8kkOATYwhD/g0YifJ1cBrmX2668/ABVX1l36rmnxJHgAOAv7W7brdpzzHI8nZwJeBtcBjwNaqeme/VU2uJGcAXwRWAV+vqkt7LmkQknwXOBU4DPgr8KmquqrXogYgyVuA24B7mP27BfCJqrq+v6omX5I3AN9g9nPgBcAPquozI5/XsCVJktSOw4iSJEkNGbYkSZIaMmxJkiQ1ZNiSJElqyLAlSZLUkGFLkiSpIcOWpImX5JkkW5Pcm+TuJB9N8n+fb0lOTfJ4kpHmIkrygSSXL7D/lCT3JRlptmlJw2LYkjQE/6qqE6rqeOA04AwWX27rtqo64393Jlk1ahFVdVt3bUn6L8OWpEGpql3AZuDCbkmTRXU9XTcn+Q6zM3GT5EdJ7ux6yTbPO/aDSf6Q5FfAyS1/BknDsqIWopa0MlTVg90w4suZXSLmuWwENlTVH7vt86vq0W5dtC3dUl4HAp8G3gw8DtwM/K5N9ZKGxrAlaaies1drnjvmBS2AD3XrTgKsB44FjgBuqardAEm+D9M0p6AAAADvSURBVLxmbJVKGjTDlqTBSfJq4Blg114c/s9533cqsAk4qaqeSnILcHD3tgvJSton3rMlaVCSrAW+ClxeVUsNSC8G/t4FrdcBJ3b7fwucmuRlSVYD546vYklDZ8+WpCE4JMlWYDWwB/gWcNk+nOfnwAVJfg/sAG4HqKqHk1wC/AZ4GLgLWAWQ5Exgqqo+OeoPIWmYsvT/+EnSZOqGCS+qqnc1vMYrgZ9W1YZW15A0WRxGlLSSPA1sGHVS08UkOQX4CfBIi/NLmkz2bEmSJDVkz5YkSVJDhi1JkqSGDFuSJEkNGbYkSZIaMmxJkiQ19B99/H+Rlvn8hgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 720x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"SF = sample_Fisher(1000, mu=(-1, 0, 0), kappa=10);\n",
"fig, ax = plt.subplots(1, 1, figsize=(10, 5))\n",
"\n",
"D = np.arctan2(SF[1], SF[0])\n",
"I = np.arctan2(SF[2], np.sqrt(SF[0]**2 + SF[1]**2))\n",
"\n",
"ax.set_aspect(1)\n",
"ax.set_xlim((-np.pi, np.pi));\n",
"ax.set_xlabel('D [rad.]')\n",
"ax.set_ylim((-np.pi/2., np.pi/2.));\n",
"ax.set_ylabel('I [rad.]')\n",
"ax.scatter(D, I, alpha=0.4);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Generate synthetic data\n",
"We first need a set of coefficients for the field. At the time of writing this, IGRF-13 \\[1\\] had just been released and we take the reported coefficients as a reference model. You can get it [here](https://www.ngdc.noaa.gov/IAGA/vmod/coeffs/igrf13coeffs.txt). Place the file in the `/dat/` folder."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"IGRF = pd.read_csv('../dat/igrf13coeffs.txt', header=0, delim_whitespace=True, skiprows=3)\n",
"coeffs = IGRF[['2020.0']].to_numpy().flatten()\n",
"# retrieve the maximal degree using pyfield and the index of the last entry in coeffs\n",
"l_max = pyfield.i2lm_l(len(coeffs)-1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then we generate random points as points of observation:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"def ran_sph(npoints, r=1., ndim=3, nocluster=True):\n",
" \"\"\" Sample random points on the ndim-sphere.\n",
"\n",
" Parameters:\n",
" -----------\n",
" npoints : int\n",
" The number of points to sample\n",
" r : float, optional\n",
" The radius of the sphere to be sampled on\n",
" ndim : int, optional\n",
" The dimension of the sphere\n",
" nocluster : bool, optional\n",
" Whether toexclude points sampled outside of the ball, i.e. sample \n",
" uniformly distributed\n",
"\n",
" Returns:\n",
" --------\n",
" Array of shape (ndim, npoints) including the sampled points.\n",
"\n",
" References:\n",
" -----------\n",
" [1] http://mathworld.wolfram.com/SpherePointPicking.html\n",
" [2] http://mathworld.wolfram.com/HyperspherePointPicking.html\n",
" \"\"\"\n",
"\n",
" vec = np.random.randn(ndim, npoints)\n",
" if nocluster:\n",
" n = np.linalg.norm(vec, axis=0)\n",
" vec = vec[:, n <= 1.]\n",
" while vec.shape[1] < npoints:\n",
" ad = np.random.randn(ndim, npoints)\n",
" n = np.linalg.norm(ad, axis=0)\n",
" ad = ad[:, n <= 1.]\n",
" vec = np.concatenate((vec, ad), axis=1)\n",
" if vec.shape[1] > npoints:\n",
" vec = vec[:, 0:npoints]\n",
" vec /= np.linalg.norm(vec, axis=0)\n",
" if npoints == 1:\n",
" vec = vec.flatten()\n",
" return vec*r"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"if real_locs:\n",
" # load reference data\n",
" pars = load('../examples/Example_Parfile.py')\n",
" # ungroup the data\n",
" ref_data = pars.data.filter(lambda x: True)\n",
" ref_data.reset_index(inplace=True)\n",
" n_points = len(ref_data)\n",
"\n",
" x_obs = np.zeros((3, n_points), order='F')\n",
" x_obs[0] = ref_data['co-lat']\n",
" x_obs[1] = ref_data['lon']\n",
" x_obs[2] = ref_data['rad']\n",
"else:\n",
" x_obs = ran_sph(n_points, r=pyfield.REARTH)\n",
" # transform x_obs from cartesian coordinates to co-lat, lon, rad\n",
" pyfield.mapLoc(fromSys=pyfield.SYS_GEO,\n",
" fromForm=pyfield.COOR_CAR,\n",
" toSys=pyfield.SYS_GEO,\n",
" toForm=pyfield.COOR_CLR,\n",
" t=0,\n",
" x=np.asfortranarray(x_obs))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we use the `pyfield` library to get the basis functions and generate a field from the coefficients."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"dspharm = np.empty((len(coeffs), 3*n_points), order='F')\n",
"pyfield.dspharm(src=pyfield.SOURCE_INTERNAL,\n",
" gSys=pyfield.SYS_GEO,\n",
" atSys=pyfield.SYS_GEO,\n",
" atForm=pyfield.COOR_CLR,\n",
" bSys=pyfield.SYS_GEO,\n",
" bForm=pyfield.FIELD_NED,\n",
" lmax=l_max,\n",
" R=pyfield.REARTH,\n",
" t=0.,\n",
" at=x_obs,\n",
" B=dspharm)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The fields $N, E, Z$ components are then easily calculated using a dot product with the coefficients."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"field_obs = np.dot(coeffs, dspharm).reshape(-1, 3).T"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using a utility function from `CORBASS`, we can transform these to $D,I,F$."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"decs, incs, ints = nez2dif(*field_obs)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we set up error levels for the components and add some synthetic noise."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"if noise:\n",
" mus = np.array([np.cos(np.deg2rad(incs))*np.cos(np.deg2rad(decs)),\n",
" np.cos(np.deg2rad(incs))*np.sin(np.deg2rad(decs)),\n",
" np.sin(np.deg2rad(incs))])\n",
"\n",
" # the angular components retrieve an error from the Fisher-vonMises distribution\n",
" # the intensities retrieve an error from the gamma distribution\n",
" for it, mu, d, i, f in zip(np.arange(n_points), mus.T, decs, incs, ints):\n",
" samp = sample_Fisher(1, mu=mu, kappa=kappa)\n",
" decs[it] = np.rad2deg(np.arctan2(samp[1], samp[0]))\n",
" incs[it] = np.rad2deg(np.arctan2(samp[2], np.sqrt(samp[0]**2 + samp[1]**2)))\n",
" b = ints[it]/dint**2\n",
" a = ints[it] * b\n",
" ints[it] = gamma.rvs(a=a, scale=1./b)\n",
"\n",
"# mimic some incompleteness in the data\n",
"if r_inc != 0.:\n",
" if real_locs:\n",
" # incompleteness of reference data\n",
" decs[ref_data.query('not (D==D)').index] = np.nan\n",
" incs[ref_data.query('not (I==I)').index] = np.nan\n",
" ints[ref_data.query('not (F==F)').index] = np.nan\n",
" else:\n",
" if r_inc != 0.:\n",
" # generate indices of missing points\n",
" ind = np.arange(n_points)\n",
" np.random.shuffle(ind)\n",
" mnum = np.int(r_inc*n_points)\n",
" mind = ind[0:mnum]\n",
" # generate a boolean array of triples of True and False for missing\n",
" # values at each point\n",
" mdist = np.zeros((3, n_points), dtype=bool)\n",
" bools = [True, True, False, False]\n",
" for j in mind:\n",
" np.random.shuffle(bools)\n",
" mdist[:, j] = bools[0:3]\n",
" # set the missing values to nan\n",
" decs[mdist[0]] = np.nan\n",
" incs[mdist[1]] = np.nan\n",
" ints[mdist[2]] = np.nan"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally we build a pandas `DataFrame` from the synthetic data and store it."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"synth_data = pd.DataFrame({'co-lat': x_obs[0], \n",
" 'lat': 90-x_obs[0],\n",
" 'lon': x_obs[1], \n",
" 'rad': x_obs[2],\n",
" 't': 2020,\n",
" 'dt': 0,\n",
" 'D': decs, \n",
" 'I': incs, \n",
" 'F': ints,\n",
" 'dD': ddec,\n",
" 'dI': dinc,\n",
" 'dF': dint})\n",
"\n",
"out_frame = synth_data.to_csv()\n",
"outfile = open(out, \"w\")\n",
"outfile.write(header + out_frame)\n",
"outfile.close()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References\n",
"\\[1\\] P. Alken, E. Thebault, C. Beggan, H. Amit, J. Aubert, J. Baerenzung, T.N. Bondar, \n",
" W. Brown, S. Cali, A. Chambodut, A. Chulliat, G. Cox, C. C. Finlay, A. Fournier, \n",
" N. Gillet, A. Grayver, M. Hammer, M. Holschneider, L. Huder, G. Hulot, T. Jager, \n",
" C. Kloss, M. Korte, W. Kuang, A. Kuvshinov, B. Langlais, J.-M. Leger, V. Lesur, \n",
" P. W. Livermore, F. J. Lowes, S. Macmillan, W. Magnes, M. Mandea, S. Marsal, \n",
" J. Matzka, M. C. Metman, T. Minami, A. Morschhauser, J. E. Mound, M. Nair, \n",
" S. Nakano, N. Olsen, F. J. Pavon-Carrasco, V. G. Petrov, G. Ropp, M. Rother, \n",
" T. J. Sabaka, S. Sanchez, D. Saturnino, N. Schnepf, X. Shen, C. Stolle, \n",
" A. Tangborn, L. Tner-Clausen, H. Toh, J. M. Torta, J. Varner, P. Vigneron, \n",
" F. Vervelidou, I. Wardinski, J. Wicht, A. Woods, Y. Yang, Z. Zeren and B. Zhou, \n",
" \"International Geomagnetic Reference Field: the thirteenth generation\", \n",
" submitted to Earth, Planets and Space, see also: \n",
" https://www.ngdc.noaa.gov/IAGA/vmod/igrf.html"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment