Skip to content

Instantly share code, notes, and snippets.

@kwinkunks
Created November 8, 2021 18:47
Show Gist options
  • Save kwinkunks/9fb0093496699dbeebd0411e16169701 to your computer and use it in GitHub Desktop.
Save kwinkunks/9fb0093496699dbeebd0411e16169701 to your computer and use it in GitHub Desktop.
Publishable graphics with Python
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"id": "4a2559a5",
"metadata": {},
"source": [
"# \"Publishable graphics with MATLAB\"... in Python :P\n",
"\n",
"I'm trying to replicate [this nice tutorial by Martin Trauth](https://mres.uni-potsdam.de/index.php/2017/02/09/create-publishable-graphics-with-matlab/). \n",
"\n",
"It generates data and fits a model to produce this plot:\n",
"\n",
"<img src=\"http://141.89.112.21/wp-content/uploads/2017/02/publishable_graphcs_vs3.png\" width=600 />\n",
"\n",
"I'm going to stick somewhat closely to Martin's code, but with some exceptions:\n",
"\n",
"- I'll use [recommended best practice](https://towardsdatascience.com/stop-using-numpy-random-seed-581a9972805f) for NumPy's random number generation pattern (not just `np.random.seed(0)`, which is not safe).\n",
"- I'll use `matplotlib`'s object-oriented API, because it's more flexible than the procedural MATLAB style of the PyPlot interface.\n",
"- I'm not going to put all the variables into one array, because it's not my usual practice; I tend to see arrays as homogenous 'lumps' of data, not tables of data (I'd use `pandas` for that). This saves a lot of indexing to 'reach' into the `data` array."
]
},
{
"cell_type": "markdown",
"id": "08558328",
"metadata": {},
"source": [
"## Make some data"
]
},
{
"cell_type": "code",
"execution_count": 99,
"id": "9f427f45",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy.optimize as so\n",
"\n",
"rng = np.random.default_rng(999)\n",
"\n",
"# How many data points do we want?\n",
"n = 26\n",
"\n",
"# Generate synthetic variables x, y and error e.\n",
"x = np.linspace(0.5, 3.0, n) + 0.2 * rng.standard_normal(n)\n",
"y = 3 + 0.2 * np.exp(x) + 0.5 * rng.standard_normal(n)\n",
"e = np.abs(rng.standard_normal(n))"
]
},
{
"cell_type": "markdown",
"id": "21187b6f",
"metadata": {},
"source": [
"Notice that `y` is generated here by the function \n",
"\n",
"$$ \\hat{y} = f(x) = c_0 + c_1 \\mathrm{e}^{x}$$\n",
"\n",
"where $c_0$ and $c_1$ are the coefficients of the function, equal to 3.0 and 0.2 respectively.\n",
"\n",
"This is the function we're going to try to fit. But we don't **know** the function, so I have not defined a Python function for it yet. We're pretending it's hidden &mdash; buried in the data, if you like."
]
},
{
"cell_type": "markdown",
"id": "bfe32977",
"metadata": {},
"source": [
"## Fit a model\n",
"\n",
"We'll use [`scipy.optimize.curve_fit()`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html) because it's convenient and I think it does exactly what we want.\n",
"\n",
"First we define the model as a Python function expressing the function $f$ above."
]
},
{
"cell_type": "code",
"execution_count": 100,
"id": "c4a1c9fc",
"metadata": {},
"outputs": [],
"source": [
"def model(t, c0, c1):\n",
" \"\"\"The non-linear model we're fitting.\"\"\"\n",
" return c0 + c1 * np.exp(t)\n",
"\n",
"# Define a domain.\n",
"t = np.linspace(x.min(), x.max(), n)\n",
"\n",
"# Unweighted solution.\n",
"(c0, c1), pcov = so.curve_fit(model, x, y, p0=(0, 0), method='lm')\n",
"y_hat_fit = model(t, c0, c1)\n",
"\n",
"# Weighted solution.\n",
"(c0, c1), pcov = so.curve_fit(model, x, y, p0=(0, 0), sigma=e, method='lm')\n",
"y_hat_fit_weighted = model(t, c0, c1)"
]
},
{
"cell_type": "markdown",
"id": "546b9e68",
"metadata": {},
"source": [
"We can take a look at these coefficients; we know they should be 3.0 and 0.2:"
]
},
{
"cell_type": "code",
"execution_count": 101,
"id": "6093168b",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(2.954089846079857, 0.21970627478563665)"
]
},
"execution_count": 101,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"c0, c1"
]
},
{
"cell_type": "markdown",
"id": "88350bd6",
"metadata": {},
"source": [
"## Make a plot"
]
},
{
"cell_type": "code",
"execution_count": 102,
"id": "c2e65142",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfMAAAFJCAYAAACPXsRYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABa20lEQVR4nO3dd3hUVfrA8e8hhdBraAFpUkMJVTqDKIigoqKAsUQFRM2yYkH0pxBkXVxkV0RURJSgorBiBURxEaQqAsYSQDpCpISSQIBAMjm/P86dYZLMJJOQZDLJ+3meeSa3zL3n3juZ955zT1Faa4QQQgjhv8r4OgFCCCGEuDISzIUQQgg/J8FcCCGE8HMSzIUQQgg/J8FcCCGE8HMSzIUQQgg/J8FcFDtKqVillFZKxfg6LZ4opdoqpbYopS5ZaW3j6zTlRikVZaV1TR4+U2TXQil1wNrX0MLel8s+Heckroj2t8baX1QePqOtV6PCSxkopWzWfpJyWKeyUuoLpdRZa91of/h/LQ0kmJcQSqmeSqmlSqmTSqlUpdRepdRrSqlgX6ctH1YCrwI/+DohOfgX0AnYgknrCd8mxyvbMWldUpAbdQnCtoLcbgm1BHMNthfUBl2C8IGC2mYOxgI3AyeBWUAcWf5f83PTKK5coK8TIK6cUmoE8AEQAPwC/AQ0wvzjPQ9c8lni8kgpFai1/hD40NdpyUVz6/05rfV3Pk2Jl7TWm4HNvk5Haaa1nu3rNFwhx/f+Pa31JJf5xf3/teTTWsvLj19AecxdsgbeB8q4LGsKBFt/twO+xuQgE4GlQAuXdQ9Y25gCxAMpwEygNbAVOAt8BJS11o+y1l+HuUM/A+wDIl22+SSwGzgHXMTcaAxzWR5rbeMt4FvMTYfNZX6MtV5Haz9nrHT9Djzssp1bMTcwZ4GDwOtAVWtZI2tbGngA+BM4DbySy3nNaZsHXLapzb+R22041rNlOWdrrGmbNX0AeBY4br2espaHWsuTMaVo1YEM4Ki1vIG1/DigrHkPWOc5xTr3zwKB7vZvzXsEOIT5XkxwSfPQLNfoNeBT4DzwKxDh6VwAUdaymzE3D2esc/hvoLybfScCT2Xdt5vzeT3ws3U+0qxtTnFZ7ji+9cArQBKQQObvZD1MTvIc5js1xfpMnId9PmUtn21NP25N/8uanmRNv2xN18B8nw9gvjsbgN4u21uT5RxVARZb5+hXl+0nuXzGcV5HATus7X4ABHP5O5Tt+4j5bXgJ2GMd7zbXc5tl37+423eWcxHrZl82l/kxLtfA9XXA17+TpeHl8wTI6wovoPmBc/zTtPCwTl1MANPAMuvHTANHgGrWOo4f0rPAAkxg1cApYCHmx14Do6z1Hf+0GcAmTPGhY7qdtc7rmADwOvBfIB1IBRpZy11/HNYA72ICt/PHwVpvvTW9BHgb8yM8z1o2yFp20Ur379b019byRi77+NM6Frs13d/D+cptm5OsH0BHmmZ62I7jnNqynLM11rTNJW3brXPlOIfNrHXirXntgCEu618NjHCkwVr3IZfjjAV2WtOTc9l/BiZn9ZvLuRma5RplAJ+5pGedp3MBdAUGWvMSMTeZW63p+Vn2bbeuyS9Z9+3mfEZhbvrewnxXTlnrj8hyfBpzE+H4nqcAla111ljzdgDvYb6POQXzLtbybda043u+3pr+2poegrnhcnxX11rHdRa4gPW/SfZgvsCaPgzMdzmXSS5pcBxTonU9LljTD2K+B440nbHO/0zrcx9Z87dax5poXUdbln0fss5nctZ9ZzkXd2G+pxpTpD7T2n8sl4N5V5fzfthaZ5KvfydLw8vnCZDXFV5AiHT5Zw/xsM4Ea/lql3k/W/PGWNMHrOnnrGnHj85/rel/W9OvW9NR1vQxIMia9xmZcykVgHuBqZic0hFr+V3WcsePwPdZ0uv8cbCmf+RyzroNEAQEWMu+InPAqonJtWlMkWAjl/PTxVrne2v6SQ/nK8dtZjlfthyuTaZ18BxM04E61ryD1rxh1vQb1vRYYBqwH3OTEYXJLWvgb9a6jkDrCKofWNNHPex/HpkDbKjLcQ7Nci2WW9P9rOkUT8dpzVtuzVtppcVxHBmYHKNj3+9Y61fn8g3kUA/nswxwI/Ac5vv0k7X+3CzHdxIIwXxP0q15nYH6Lt+FBlm+13Ee9hmACZLpQEXMd/gPzE1AOUzu3w5U5XLgdwZVTG5YAy9l+b+KsrZ90Zruay0fj+dgfoc17QjCjtICG1lywFwu1bFjviczXfa9KMu+e1uf+VvWfbs5H47vQ0wO/6+O67DG17+Ppeklz8z933GXvxtifmiyamS973CZtxOIsD7jyrFOkvXu2N5Z671ClvX3aq3TXLYJUN+qePcDJvhmFZpleqObdVw9jgkG8wCFyWlNwvygN3JNt9b6hFLqBFAHc2y7Xbbzs/WeZL1X9LC/3La5K5f0ehLgYf5RrfVRl7Rd5ZK2NcDDQA8rXauBVkBPTCmGYx3XdN+eZfu1lVLujjXMenccZ6LLcWaV9dxl/R5k5UjL9dbLQQFNXPb9h7XvU0qpkx727fAmMMbN/Kzfpx1a61QApdQ5oDLmfDr2eUFrfcj6O8drqbW2K6XWY0prRljpm4wpHYjCFFVv01onudQ2rwT8Pcumrnaz+ZqYonK4/H+XU8U4b7+/cPn8lwGi3aTFdd+O//H8fq9FMSC12f3fRkwROsBzSinnNVVKNVRKBWFyTgAtXT7Xwno/mGV79lyms2pq7cN1+4cxz9rbWJ9vhvmuOX6oVJZtXMxlH1u01u2BaphcSBDwklIqkCzHppSqgfmhgizHprVOd/yZy/683mYuzlnvla13T83X0l3+zpq27633Ppic3wbrdR3QHpML/T1Lum/WWivHC2iitU5xs98E670ZgFKqJpeP01Ma3Z07x3fE9ffEkZZxWdLSVGv9u8u+W1j7ro553pyT4dZ7FObG6E1rOuv3ydP5dOyznFKqgfV3c3LnuAZPYEouPsA8dnoiy/ID1vtfmFIyxzGXJ3tAxdqGo3JqM+u9pZv1HDxdg5zO/yUg1CUtwZj6IK77dvwWeHMuvOEuPaKQycn2c1rrc5jisQzgbmCbUmquUmoZ5k67AubHJxnop5T6Uin1NdABU0R+pc2UagLfK6WWAEMxPzSOZ+wZmB/d/2CedTbzsI3cLFVKrQJexlSaKospKbBjnscDPKuUisXkUgOBb7XW+c1pFNQ2HTmpfyilXrXSnida62OYEo+GmKJjRzBvhDm3a7VVtgk4akp/YLX9fU8ptR3zLNad9633+5VSC4HvyN9vgiOX+4JSaqYVKB1pma6U+kgp9a5SaguwyprvqP0cZe17Nbm3rjlmvY+z0h6Vl0RqrQ9jnmUDrFRKvYf7IJvVGuu9JSYXfh5zE900y/KtmPoj9YCflFJzlFKfY4L7DW7SY8c81wb4SCn1LvBCHg7JwXH+6yul5imlntZaJ2LqqQQDP1pp+dha90Fr34utz31o7fvFfOw7p/R0Ukq9oZQaXUDbFTmQYF4CaK0XYp5lfoUpor0PUxT7NnBea/2XtXwlpni2M+aZZj+t9akr3P0GzA/b9Zhc631a6zjrh/NvmB/gvpgfutyK0z1Zg/mBjAQGY56VDtfGcuBOzPPiYZhiz7e4nIvLswLc5nOYH/cmmCLx/DZLcuT8TmKKRDe4WQYwB1PjeT8m3Tdibqrmuduo1vp74FHMc+AbMDdhjoCZW2mJqxhMjenumOLl2lrrFZgc4C9WOm7D3Ny9au37O0xQPoIpwv4EU3EvJ6MwNzatMUXZb+UhjQ6RwP8wN0fNMTeaudmKebQDl8+9411jKmSitc4AbsFch8qYm40OmP9LT30m/B34GFPq1BnTfwHk4fxrrQ8AMzA37A8C91iLHsTUZs+w0tIT83382lo+DnMzXxVTcW2at/vMxVrMzZod84jolgLarsiBunxTL4T3rB6s5mMqr9l8mxqRX0qpKlrrZOvv+pgbsjLA1VrrvT5NXCmglKqEqUyorelngH9iasv39mnihF+RCnBClG4/K6W+wuT6R2AC+VcSyItMf0xdlxWYOgP3W/Nn+S5Jwh8VWjG79YzsuFLqd5d51ZVS3yqldlvv1Qpr/0IIr2zDBPGnMRXJZmDaE4ui8Sem7sMTmPP+C3Cn1vpjn6ZK+J1CK2ZXSvXBPGd6T2vdxpo3HTiltX5JKTUR02HJ04WSACGEEKKUKNRn5la7y2UuwfwPTMcSR5RSdTGdCrTIaRtCCCGEyFlR12avrbU+AmC91yri/QshhBAlTrGsAFexYkXtWmJQs2ZNQkOzdvIkhBBClBxbt249obXOV7Ar6mB+TClV16WY/bi7lVq2bMmWLVuKOGlCCCGE7yil8tLDZCZFXcz+JaZDE6z3L4p4/0IIIUSJU5hN0z7C9DbUQil1WCnl6I3oeqXUbkyPYS8V1v6FEEKI0qLQitm11iM9LOpfWPsUQgghSqNiWQHOnbS0NA4fPkxqaqqvkyJEsRYSEkL9+vUJCgrKfWUhRIngN8H88OHDVKpUiUaNGqFU1hEPhRAAWmtOnjzJ4cOHady4sa+TI4QoIn4zalpqaio1atSQQC5EDpRS1KhRQ0qwhChl/CaYAxLIhfCC/J8IUfr4VTD3NaUUTzzxhHN6xowZxMTE5GtbBw4coE2bNgBs2bKFcePGFUQSvRYTE8OMGTMKbHsHDhzgww8/LLDtCSGE8J4E8zwoW7Ysn376KSdOnCjQ7Xbu3JlZswp3xEO73V6o2y8uwVxrTUZGhq+TIYQQRUqCeR4EBgYyZswYXnnllWzLDh48SP/+/WnXrh39+/fnzz//BCAqKopx48bRo0cPmjRpwpIlS7J9ds2aNQwZMgQwOeYHHngAm81GkyZNMgX5Dz74gK5duxIREcFDDz3kDNAPP/wwnTt3Jjw8nMmTJzvXb9SoES+88AK9evXi4489j6i4d+9ebrjhBjp16kTv3r3ZuXMnAEuXLuWaa66hQ4cOXHfddRw7dgyA77//noiICCIiIujQoQNnz55l4sSJrFu3joiIiGzn58iRI/Tp04eIiAjatGnDunXrAJg/fz7Nmzenb9++jB49mujoaOc5cz1PFStWBCAlJYX+/fvTsWNH2rZtyxdfmD6HDhw4QKtWrXjkkUfo2LEjhw4d4uWXX6ZLly60a9fOeU7OnTvH4MGDad++PW3atGHx4sUez4kQQvgTv6nNnsljj0FcXMFuMyICZs7MdbVHH32Udu3aMWHChEzzo6Ojuffee7nvvvt49913GTduHJ9//jlggtn69evZuXMnN998M8OGDctxHzt37mT16tWcPXuWFi1a8PDDD7Nnzx4WL17Mhg0bCAoK4pFHHmHhwoXce++9vPjii1SvXh273U7//v359ddfadeuHWCaKa1fvz7H/Y0ZM4Y5c+bQrFkzfvzxRx555BG+++47evXqxQ8//IBSinnz5jF9+nT+/e9/M2PGDF5//XV69uxJSkoKISEhvPTSS8yYMYNly5Zl2/6HH37IwIED+b//+z/sdjvnz5/nyJEjTJ48ma1bt1KlShX69etHhw4dckxnSEgIn332GZUrV+bEiRN069aNm2++GYA//viD+fPn88Ybb7By5Up2797N5s2b0Vpz8803s3btWhITE6lXrx7Lly8HIDk5Ocf9CSGEv/DPYO5DlStX5t5772XWrFmUK1fOOX/Tpk18+umnANxzzz2Zgv3QoUMpU6YMrVu3duZuczJ48GDKli1L2bJlqVWrFseOHWPVqlVs3bqVLl26AHDhwgVq1TKDzv33v/9l7ty5pKenc+TIEbZv3+4M5sOHD89xXykpKWzcuJE77rjDOe/ixYuAaQ44fPhwjhw5wqVLl5xNnXr27Mnjjz9OZGQkt912G/Xr189xH126dOGBBx4gLS2NoUOHEhERwapVq7DZbM4BdIYPH86uXbty3I7WmmeffZa1a9dSpkwZEhISnOezYcOGdOvWDYCVK1eycuVK581BSkoKu3fvpnfv3jz55JM8/fTTDBkyhN69e+e4PyGE8Bf+Gcy9yEEXpscee4yOHTty//33e1zHtUZx2bJlnX97M3686/oBAQGkp6ejtea+++5j2rRpmdbdv38/M2bM4KeffqJatWpERUVlapZUoUKFHPeVkZFB1apViXNT0vG3v/2Nxx9/nJtvvpk1a9Y4K/tNnDiRwYMH89VXX9GtWzf+97//5biPPn36sHbtWpYvX84999zDU089ReXKlT3Wug4MDHQ+99Zac+nSJQAWLlxIYmIiW7duJSgoiEaNGjmP1fU4tdY888wzPPTQQ9m2vXXrVr766iueeeYZBgwYwKRJk3JMuxBC+AN5Zp4P1atX58477+Sdd95xzuvRoweLFi0CTNDp1atXge6zf//+LFmyhOPHzUBzp06d4uDBg5w5c4YKFSpQpUoVjh07xooVK/K03cqVK9O4cWPnM3WtNb/88gtgiqHDwsIAWLBggfMze/fupW3btjz99NN07tyZnTt3UqlSJc6ePet2HwcPHqRWrVqMHj2aBx98kG3btnHNNdewZs0aTp48SVpaWqZn+o0aNWLr1q0AfPHFF6SlpTnTU6tWLYKCgli9ejUHD7ofYGjgwIG8++67pKSkAJCQkMDx48f566+/KF++PHfffTdPPvkk27Zty9O5EkKI4so/c+bFwBNPPMHs2bOd07NmzeKBBx7g5ZdfJjQ0lPnz5xfo/lq3bs0//vEPBgwYQEZGBkFBQbz++ut069aNDh06EB4eTpMmTejZs2eet71w4UIefvhh/vGPf5CWlsaIESNo3749MTEx3HHHHYSFhdGtWzf2798PwMyZM1m9ejUBAQG0bt2aQYMGUaZMGQIDA2nfvj1RUVGMHz/euf01a9bw8ssvExQURMWKFXnvvfeoW7cuMTExdO/enbp169KxY0dnhb7Ro0dzyy230LVrV/r37+/MdUdGRnLTTTfRuXNnIiIiaNmypdvjGTBgADt27KB79+6AqUD3wQcfsGfPHp566inKlClDUFAQb775Zp7PlRBCFEfKm2Lfota5c2eddTzzHTt20KpVKx+lSBS22NhYtmzZkukGSeSf/L8I4X+UUlu11p3z81kpZhdCCCH8nBSzi2IhKiqKqKgoXydDCCH8kuTMhRBCCD8nwVwIIYTwcxLMhRBCCD8nwVwIIYTwcxLMvTR+/HhmuvQ8N3DgQEaNGuWcfuKJJ/jPf/7j8fOTJk3Ktac0T8OSJiUl8cYbb+Q5zQU9zKk7NpsNRzPCG2+8kaSkpELdnyvXYWQLysyZMzl//nyBblMIIQqbBHMv9ejRg40bNwKmC9QTJ04QHx/vXL5x48YcO2x54YUXuO666/K17/wG86L21VdfUbVq1ULbfnp6eqFt26G4BPOiOFYhRMkhwdxLPXv2dAbz+Ph42rRpQ6VKlTh9+jQXL15kx44ddOjQga1bt9K3b186derEwIEDOXLkCJB5WM+vvvqKli1b0qtXL8aNG+cc/hRg+/bt2YY/nThxInv37iUiIoKnnnoKwO0QnwAvvvgiLVq04LrrruOPP/5weyyehmXVWvPUU0/Rpk0b2rZt6xwidM2aNdhsNoYNG0bLli2JjIx028d8o0aNOHHihHNI0tGjRxMeHs6AAQO4cOECkPfhVmNiYhgzZgwDBgzg3nvv9Xh97HY7Tz31lPOcvPXWW4DnYVPdDYc6a9Ys/vrrL/r160e/fv2y7WPixIm0bt2adu3a8eSTTwKmb/zu3bvTpUsXnn/+eedwra7D2oIZVS82NhYwN3ZdunShTZs2jBkzxnkubTYbzz77LH379uXVV1/1+F2aNWuWMx0jRozweE6EEKWI1rrYvTp16qSz2r59++WJo3/X+kDfgn0d/Xu2fWbVsGFDffDgQT1nzhz95ptv6ueee04vX75cr1+/Xvfu3VtfunRJd+/eXR8/flxrrfWiRYv0/fffr7XW+r777tMff/yxvnDhgq5fv77et2+f1lrrESNG6MGDB2uttZ48ebLu3r27Tk1N1YmJibp69er60qVLev/+/To8PNyZjm+++UaPHj1aZ2RkaLvdrgcPHqy///57vWXLFt2mTRt97tw5nZycrJs2bapffvnlbMdx33336WHDhmm73a7j4+N106ZNtdZaL1myRF933XU6PT1dHz16VDdo0ED/9ddfevXq1bpy5cr60KFD2m63627duul169ZprbXu27ev/umnn5znJzExUe/fv18HBATon3/+WWut9R133KHff/99rbXW1157rd61a5fWWusffvhB9+vXT2ut9alTp3RGRobWWuu3335bP/74485z0rFjR33+/Plsx+F6Xt566y09depUrbXWqampulOnTnrfvn06LS1NJycna621TkxM1E2bNtUZGRl6yZIletSoUc5tJSUlZTqGrE6ePKmbN2/uTOPp06e11lrfdNNNesGCBVprrWfPnq0rVKigtdZ69erVzuuqtdaPPvqonj9/vnNbDnfffbf+8ssvnefy4Ycf1lrrHL9LdevW1ampqZnSkVWm/xchhF8Atuh8xk3pNCYPHLnzjRs38vjjj5OQkMDGjRupUqUKPXr04I8//uD333/n+uuvB0xusW7dupm2sXPnTpo0aeIcTnTkyJHMnTvXudzd8KdZeRri8+zZs9x6662UL18ewDnWtzvuhmVdv349I0eOJCAggNq1a9O3b19++uknKleuTNeuXZ1DnUZERHDgwIEcB5Np3LgxERERAHTq1IkDBw7ka7hVx3G4DjfrzsqVK/n111+dpQzJycns3r2b+vXrux02tW3btnkaDrVy5cqEhIQwatQoBg8e7Mx1b9iwgU8++QQwQ98+/fTTOW4HYPXq1UyfPp3z589z6tQpwsPDuemmm4DLQ9bm9F1q164dkZGRDB06lKFDh+a6PyFEyeefwbz2TJ/s1vHc/LfffqNNmzY0aNCAf//731SuXJkHHngArTXh4eFs2rTJ4zZ0Ln3huxv+1N023A3xOXPmTI/Diua0H0eackqbN+nKaf0LFy7ka7hVyH0YV0faX3vtNQYOHJhpfmxsrNthU5s3b56n4VADAwPZvHkzq1atYtGiRcyePZvvvvsOwO05dx3GFXAO1ZqamsojjzzCli1baNCgATExMW6HrM3pu7R8+XLWrl3Ll19+ydSpU4mPjycw0D//lYUQBUOemedBz549WbZsGdWrVycgIIDq1auTlJTEpk2b6N69Oy1atCAxMdH5A5yWlpapkhxAy5Yt2bdvHwcOHABwPpfOSdbhRT0N8dmnTx8+++wzLly4wNmzZ1m6dGmejq9Pnz4sXrwYu91OYmIia9eupWvXrnnaRk7yM9yqtwYOHMibb77pHC51165dnDt3zuOwqZ6GQ/U0lGtKSgrJycnceOONzJw503lD0rNnz0xD3zo0bNiQ7du3c/HiRZKTk1m1ahVwOajXrFmTlJQUZ0lCVp6+SxkZGRw6dIh+/foxffp0kpKSnN8DIUTpJbfzedC2bVtOnDjBXXfdlWleSkoKNWvWBGDJkiWMGzeO5ORk0tPTeeyxxwgPD3euX65cOd544w1uuOEGatas6VWwrFGjBj179qRNmzYMGjSIl19+2e0Qnx07dmT48OFERETQsGHDXIuOs7r11lvZtGkT7du3RynF9OnTqVOnjrOSWkHI63Cr3ho1ahQHDhygY8eOaK0JDQ3l888/9zhs6m+//eZ2ONQxY8YwaNAg6taty+rVq53bP3v2LLfccgupqalorXnllVcAePXVV7nrrrt49dVXuf32253rN2jQgDvvvJN27drRrFkz5yORqlWrMnr0aNq2bUujRo3o0qWL2+MJDg52+11q3rw5d999N8nJyWitGT9+fKG2IBBC+AcZAtUHUlJSqFixIlprHn30UZo1a5Zp/G/hvypWrFgscsol6f9F5E9sbKyzBcXRo0epU6cOIIMaFWdXMgSq5Mx94O2332bBggVcunSJDh06ZHv2LYQQV8o1aNtsNtasWePT9IjCJcHcB8aPHy858RKqOOTKhRClj1SAE0IIIfycBHMhhBDCz0kwF0IIIfxciXxmLrU4hRBClCYlMmceFRXFwoUL6d67H6fOp9G9dz8WLlx4xYE8ICCAiIgIwsPDad++Pf/5z38y9fLlzoEDB/jwww+vaL/uOIYbzTqiWtYBPoQQQpR8JTKYL1y0mK79hzB/expBgyYyf3saXfsPYeGi3Htby0m5cuWIi4sjPj6eb7/9lq+++oopU6bk+JnCCuaO4UYLc3jUrF22ejsspwzfKYQQRavEBfOEhAQmxLxE4JBJhDTrTlC1eoQ0607gkElMiPkXCQkJBbKfWrVqMXfuXGbPno3WmgMHDtC7d286duxIx44dncOlTpw4kXXr1hEREcErr7zicT1X06dPdw5/On78eK699loAVq1axd133w1cHm7U3fCoKSkpuQ5X6mko0qioKB5//HH69evH008/nW06Li6Obt260a5dO2699VZOnz4NZB++8+OPP6ZNmza0b9+ePn36FMg5F0II4V6Je2Y+e8487OGDCQoMzjRfBQZjDx/E7DnzmDZ1sodP502TJk3IyMjg+PHj1KpVi2+//ZaQkBB2797NyJEj2bJlCy+99BIzZsxg2bJlAJw/f97teq769OnDv//9b8aNG8eWLVu4ePEiaWlprF+/PlsXrS+99BK///67s6/wNWvW8PPPPxMfH0+9evXo2bMnGzZsyDbC2ZgxY5gzZw7NmjXjxx9/5JFHHnEOHLJr1y7+97//ERAQQFRUVKbpdu3a8dprr9G3b18mTZrElClTmDlzJgBJSUl8//33gOnm9ptvviEsLIykpKQCOd9CCCHcK3HBPH7XbgJq9nW7LKBGQ3bsWleg+3PketPS0oiOjiYuLo6AgAB27drldn1v1uvUqRNbt27l7NmzlC1blo4dO7JlyxbWrVvnzLHnJLfhSnMaihTgjjvuICAgINt0cnIySUlJ9O1rzu99992XaRuO4TvBDEASFRXFnXfeyW233ZZrmoUQQuRfiQvm4c2bsXn7QYKq1cu2zH7yIK3Cry6wfe3bt4+AgABq1arFlClTqF27Nr/88gsZGRmEhIS4/cwrr7yS63qOoTrnz59Pjx49aNeuHatXr2bv3r1e9bed23ClOQ1FCtmHHPVmCNKs682ZM4cff/yR5cuXExERQVxcHDVq1PBqO0IIIfKmxD0zjx47ioD45ej0S5nm6/RLBMSvIPqhUQWyn8TERMaOHUt0dDRKKZKTk6lbty5lypTh/fffx263A9mH1PS0XlZ9+vRhxowZ9OnTh969ezNnzhwiIiKyjZ3tacjOnOQ0FGlOqlSpQrVq1Vi3zpRuvP/++85celZ79+7lmmuu4YUXXqBmzZocOnQoT2kUQgjhvRIXzMPCwpgeM5H0ZVNJ3b2RtFMJpO7eSPqyqUyPedo5ZnZ+XLhwwdk07brrrmPAgAFMnmyevz/yyCMsWLCAbt26sWvXLmcutV27dgQGBtK+fXteeeUVj+tl1bt3b44cOUL37t2pXbs2ISEhboc0dR0e1VEBzhsLFy7knXfeoX379oSHh/PFF1949bkFCxbw1FNP0a5dO+Li4pg0aZLb9Z566inatm1LmzZt6NOnD+3bt/c6bUIIIfKmxA6BmpCQwOw583hnwfs8eN89RI8ddUWBXAh/IkOgClcyapp/uJIhUEtczhxMD3CRkZFsWrea6uUD2bRuNZGRkc5e4YQQQoiSpMRVgAPptlUIIUTpUiJz5kIIIURp4pNgrpQar5SKV0r9rpT6SCnlvh1XFsXx+b4QxY38nwhR+hR5MFdKhQHjgM5a6zZAADAit8+FhIRw8uRJ+aESIgdaa06ePOmxnwMhRMnkq2fmgUA5pVQaUB74K7cP1K9fn8OHD5OYmFjoiRPCn4WEhDh7ABRClA5FHsy11glKqRnAn8AFYKXWeqXrOomJiXTufLl2/pgxYxgzZgyNGzcu2sQKIYQQfqDIg7lSqhpwC9AYSAI+VkrdrbX+wLFOaGhotsFHhBBCCOGeLyrAXQfs11onaq3TgE+BHj5IhxBCCFEi+CKY/wl0U0qVV6aj8f7ADh+kQwghhCgRijyYa61/BJYA24DfrDTMLep0CCGEECWFT2qza60nA5N9sW8hhBCipJEe4IQQQgg/J8FcCCGE8HMSzIUQQgg/J8FcCCGE8HMSzIUQQgg/J8FcCCGE8HMSzIUQQgg/J8FcCCGE8HMSzIUQQgg/J8FcCCGE8HMSzIUQQgg/J8FcCCGE8HMSzIUQQgg/J8FcCCGE8HMSzIUQQgg/J8FcCCGE8HMSzIUQQgg/J8FcCCGE8HMSzIUQQgg/J8FcCCGE8HOBvk6AEEKI7GJjY4mNjQXg6NGj1KlTB4CoqCiioqJ8lzBRLEkwF0KIYsg1aNtsNtasWePT9OSJ1qCUr1NRqkgxuxBCiIKTfgT+7AMpX/s6JaWK5MyFEKIYi42NJS4uDpvNVvyL2y/8AAm3gT0Z9AVfp6ZUkWAuhBDFWFRUFLGxsaxZs6Z4F7cnvQPHHoHAMGj4DYS09XWKShWPwVwp9asXn0/UWvcvwPQIIYTwJ/oSHBsPSW9A+eshbBEEVPd1qkqdnHLmAcCNOSxXwJcFmxwhhBB+I/04JAyDC+ug+pMQOg2UFPj6Qk5n/SGt9cGcPqyUeqSA0yOEEMIfpG6Fw0PBfhLqfQiVR/o6RaWax9rsWuv1uX3Ym3WEEEL4h9jYWGw2GzabjZYtWzr/drR3d0p+Dw72BMpAww2ZA/nZs3DXXbB5c1EmvdTLtTxEKdUMmAa0BkIc87XWTQoxXUIIIYpYrm3bdTocfwpOz4Ty/aDeYggMvbx850647TbYtQuuvRa6di2qpJd63rQznw+8CaQD/YD3gPcLM1FCCCGuXEJCAs88P4Xt+w7xzPNTSEhIyP/G0k/AoYEmkFf7OzT4JnMg//xzE7xPnIBvv4VRo640+SIPvAnm5bTWqwCltT6otY4Bri3cZAkhhLgSCxctpmv/IczfnkbQoInM355G1/5DWLhocd43lhoHBzrDhQ1QNxZqzwQVZJbZ7fB//we33gotW8LWrdCvXwEeifCGN9UOU5VSZYDdSqloIAGoVbjJEkIIkV8JCQlMiHmJwCGTCAoMBiCoWj10405MiJmKrXcvwsLCPH529px5ztx8dGQoYepJ09zsqnVQrsvllU+dMs/Hv/nG5MRfew1CQtxuVxQub3LmjwHlgXFAJ+Ae4N5CTJMQQogrMHvOPOzhg1FWIHdQgcHYwwcxe848t5/Lnpu/SNeb/snClU2g0dbMgTwuDjp3htWrYe5cePttYhct8q4CnShw3uTMM7TWKUAKcD+AUuqmQk2VEEKIfIvftZuAmn3dLguo0ZAdu9Zlm+85N9+FCa+/gO22dJyZ+Q8+gNGjoUYNWLsWrrkG8PPBYfycNznzt5VSzn75lFIjgOcKL0lCCCGuRHjzZthPuO8mxH7yIK1aXJ1tfs65+Rt5KPrv9O/Th0/q14d77iGubFluveoqYnfsKJRjEHnjTTAfBixQSrVSSo0GHgUGFG6yhBBC5Ff02FEExC9Hp1/KNF+nXyIgfgXRD2WvaR6/czsBNRu63V5AjYYE6jK8l5rKluMn6VK9Dosf/TuzP/64+A32UkrlGsy11vuAEcAnmMA+QGudXNgJE0IIkT9hYWFMj5lI+rKppO7eSNqpBFJ3byR92VSmxzydvfLb+bWE11mNPXGf2+3ZTxwkY/X/6Lr3CPOHPM5fw15g/k57/mvHiwLnMZgrpX5TSv1qDbiyBKgONAJ+9HIQFiGEED4SOWI4m1ct5f7WwaR9PZ37WwezedVSIkcMv7ySTofEGPizH9EjKhAQv8xtbp7V77I5uCKB988ipHkPgqrVI6RZdwKHTGJCzL+urP26KBA5VYAbUmSpEEIIUeDCwsKYNnUym9atZtrUyZkXph2EvyJN2/HK9xLWbDbTp3zFhJip2MMHEVCjIfbE/QR8/z5dK5Vjc+eROdaOz7Z9UaQ8BvPcBlkRQgjhp84sgaOjATvUXQhV7gJMbt7Wuxez/28ScR+9SUTaRaL/Hs3D+w8TUK2R2015qh0vilZO45lv01p3zOnD3qwjhBDFVWxsrLMN9NGjR6lTpw6QuYlViZJxzow9nvw2hHSFeh9BsMswG2lphL3+OtPem8+hkBAarPseunQh/PkpbN5+kKBq9bJt0n7yIK3Cs9eOF0Urp2L2Vrk8G1dAlQJOjxBCFJlS1S469Wf46y649AdUnwihL1zukhVg927Tm9uWLTBqFGN27mRFF9NJTPTYUbzXfwi6cadMRe3O2vGzlhb10YgscgrmLb34vL2gEiKEEKIQZFzkwZv3w4GuZmCUBt9Chf6Xl2sN774L48aZrlg/+QRuu40LNptzFUft+EzP008eJCB+hfva8aLI5TSe+UEvXoeLMrFCCFEaXbx4MX+jn134AQ504J4b/4TKkdD498yB/ORJGDbM9KvevTv8+qsZwtQNr2rHC5/xptOYAqeUqqqUWqKU2qmU2qGU6u6LdAghRHG3cNFiftl/NG+jnzmejR/sARkpTJjVFurFmsFSHFatgnbtYOlSePllWLmS2G+/dfanfvTo0Ux9q8fGxhIZGcmmdaupXj6QTetWExkZKf2uFxPe9M1eGF4FvtZaD1NKBWMGchFCCOHC0V966N3/cT6rznX0s3PfmZrqafug6iMQ+hKbt7sMp3HxIjz/PMyYAS1awLJl0KEDkHvFvxJZKbCEyDVnrpT6lzfzvKWUqgz0Ad4B0Fpf0lon5Xd7QghRUuVp9DN7MhwZA4f6A2XgqjVQ53UIqHR5nR07THH6yy/D2LFm7HErkAv/5k0x+/Vu5g26gn02ARKB+Uqpn5VS85RSFVxXSExMpHPnzs7X3Llzr2B3Qgjhn8zoZ577S9+xa4+ZSFkG+8Mh+R2o/iQ0/gXKXx41TWkNs2dDp05w6BB8+SW88QaUl0LRkiKnduYPA48ATbI0UasEbLjCfXYE/qa1/lEp9SowEXjesUJoaChbtmy5gl0IIYT/C2/eLOf23a3qwV93w5mFULYNhH2WecxxgPh4XouLM0OV3nCDqblet27RHIAoMjk9M/8QWAFMwwRbh7Na61NXsM/DwGGt9Y/W9JIs2xdCCAHcfstgXn97GCm7fiCoWj0qtO1PYKWapn33758R/dgxOHMGakyGms+CcimOv3gRpk2Df/6TBlrD++9DZCQo5bsDEoUmp+5ck4FkYKRSKgCoba1fUSlVUWv9Z352qLU+qpQ6pJRqobX+A+gPbM/PtoQQoqRauGgxE2JeolyPSAJDG3Hp+H4SP3+Jcg3CqXBqG9MfOkRYgwio8y6EtM384U2bTHOz7dvhrru4d/9+vrj7bp8chygaudZmV0pFAzHAMSDDmq2Bdlew378BC62a7PuA+69gW0IIkW8JCQnMnjPP2YY7euwon3eC4qjFHjhkEkEutdjLN+1C0qJxfPFqIp1tL0L1x0G5/IyfPQvPPguvvw7168Py5XDjjSS7dAAjSiZvKsA9BrTQWodrrdtarysJ5Git47TWnbXW7bTWQ7XWp69ke0IIkR8LFy2ma/8heWvDXQRyqsVe7pq7+WRbNNSYkDmQL18O4eEmkP/tbxAfDzfeWMQpNzci+ergRlwRb4L5IUxxuxBClAixsbF0796d0Y9NJHDIJEKadS9WY3TnXIu9MTv2HLs84/hxGDkShgyBSpVgwwZ49VXzdxErrjdHpYE3ncbsA9YopZYDFx0ztdb/KbRUCSFEIYqKiuKPvQfZvz2tWI7RHd6sAZt37M95lDJHpbbx4yElBaZMgYkTITjYzRYLn6dHAzl2cCMKjDc58z+Bb4FgTLM0x0sIIfyW1224i5LOgKR3iR74NgG/vodOv5R5sWOUskEDYOBAuO8+aNkSfv4ZJk3yWSCHPHZwIwpcrjlzrfUUAKVUBa31ucJPkhBCFL5c23AX9RjdFzbDsWhI/Ymwq3ow/fmbmfDPzKOUnd3wIXN7dSHs2mshMNA8Hx87Fsr4ZJiNTMzNUV+3y8zN0boiTlHp4k13rt2VUtuBHdZ0e6XUG4WeMiGEKETRY0cREL/cc+73oVFFk5D0Y3DkATh4DaQfhrrvw1Xribz3aZdRyv7F/WoPG0//ReQnH8PQoabZ2SOPFItADubmyH7ioNtl9pMHadWiiG+OShlvvgUzgYHASQCt9S+YvtWFEMJvOcboTvzgcVJ3byTtVAKpuzeSvmxq0YzRnXEOTrwI+5pD8gdQfQI0/gOq3O3s2CUsLIxptwxmVfIJpn3yX0KCg2HdOvjoI2jQoHDTl0fF5uaolPJq1DSt9SGVudcge+EkRwghioZjWM/A8ydJ3xBLUuolateszmOPji3cMbr1JUiaCyf+AfZjUPFmCJ0OZVtkXu/oUdNmPDaWsMBAeOcdHl6wgO969Sq8tF0Bx83RhJjMjwYC4lcUzc1RKedV0zSlVA9AK6WClVJPYhW5CyGEv4qKimLNmjV07NiRNi2b0S0inIN7/mD8+PGFs0NtNznwfS3h2N8guCU03Aj1v8gcyC9ehOnToXlz+OADePJJ7u7aFR54gIxi3hVr5IjhLo8GpnN/62A2r1pauDdHAvAumI8FHgXCMP2qR1jTQgghcqM1nF0KByLgyD1QpirU/xquWg3lumde78svTccvTz8NNpvp+GX6dM4HelWIWiyEhYUxbepkWjepz7SpkyVHXkS8qc1+AogsgrQIIUTJcn4tJD4DFzZCUDOotwgq3QEqSz4qPt60F//2W2jVCr7+2jQ9E8JL3vTN3hjTl3oj1/W11jcXXrKEEMKPpcZB4rNwbgUE1oM6c6FKFKigzOsdPw5Tp8Kbb5oe2159FR5+GIKC3G1VCI+8Kbv5HHgHWMrlgVaEEEJkdWkPJD4PZxdBmWqmYlu1aChTLvN6SUkwYwbMnAkXLsBDD8ELL0DNmr5ItSgBvAnmqVrrWYWeEiGE8FeXdsHJf0Hye2ZM8Rr/B9WfhICqmdc7dw5mzTIV3JKS4M47TTesLVv6ItWiBPEmmL+qlJoMrCRz3+zbCi1VQghRyFyHPi0XHEi90Op530jqL3Dyn3D2Y1BlodrDUONZCKyTeb2LF+Gtt+DFF03R+pAhpng9IqJAjkUIb4J5W+Ae4Foyj2d+bWElSgghCtPCRYuZEPMS9vDBBA2ayIVj+/jlh0UsXLTYu2ZU5zeaIH5uOZSpBNWfhuqPQWDtzOulp0NsrClCP3QI+vWDzz+H7t3dbFSI/PMmmN8KNNFaX8p1TSFEsefoLAXg6NGj1KljcpFRUVFERUX5LmFFxNPoXuWu7prz6F5aw/n/mSB+fg0E1ICaU80z8azF6RkZsHixGfxkzx645hqYPx/69y/04xOlkzfB/BegKnC8cJMihCgKrkHbZrOxZs0an6anqDlG9wryduhTnQEpX5ognvqTqZ1e6xWoOhrKVMi8ca1h6VJ47jn47Tdo2xa++AJuusnZRasQhcGbYF4b2KmU+onMz8ylaZoQwu94PbqXToczi+DkNLi0HYKamCZmle+FMmUzf1Br0zY8JgY2b4ZmzUz/6XfeWWwGQhElmzfBfHLuqwghhH/IaejTtOP7+X3HD2xacS/dm6+HtP1Qtg3U+9Dq7CXLT2Z6Onz8Mbz0Evz6K1x1FcybZ8YZt3ptu9LHGrGxscTFxWGz2Th69Cg2my1Pnxelgzc9wH1fFAkRQoiiED12FO/1H4Ju3AnlUtSu0y8RtP0Tvn/nGGGh70NAV6g1EyoOyd5j24ULpmLbyy/D/v2m17bYWBg5EoIzF99f6WONqKgoYmNjS93jEJE3HoO5Umq91rqXUuospva6cxGgtdaVCz11QgifKakV5dyO7pW4j4DfP2D62BOENbkVqo+DkG7Zn3MnJZne2mbONE3MrrkGXnnFPBOX4nThQx6Duda6l/VeqeiSI4QoLkpyRbnIO2/C1u4As+e8zo69Z2hSL50nPnmAsFZPQ5CbmuxHjpgA/uabcPYs3HADTJwIffpIxTZRLOSUM8+xBwWt9amCT44QQhSitINw+g1IepuwMqeZNj6Cl+aW57ufahHWzk1Hl3v2mKL02FjzfPzOO82IZtLZiyhmcioX2gpssd4TgV3AbuvvrYWfNCGEKABaw/nv4fDtsLcJnPo3VLgOrlpHQtBS3vxYE7crgWeen0JCQoL5zJYtMHw4tGgBCxbAAw/Arl2mhroEclEMeQzmWuvGWusmwDfATVrrmlrrGsAQ4NOiSqAQouAlJJjgtX3focxBrCSxn4bTs8044n/aTEcv1SdA030Q9l8WfplA1+tu4kLb2wkaNJH529Po2ut6FrZoBV26mKZmEybAgQOmeL1pU98ejxA58KbGRhet9VeOCa31CsB9I00hRLG3cNFiuvYfwvztaZeDWP8hLFy02NdJu3I6A859B3/dBXvqwrG/AYFQZx5cfRhqTYOgqzL1Ale+RU+CqtUjpFl3Am9/kQmnz5MwZQr8+SdMmwZ16uS629yUipsn4VPeBPMTSqnnlFKNlFINlVL/B5ws7IQJIQqeaxALadb9chAbMokJMf8q8iDjCHI3D7/7yoJc2mE48Q/YdzUc6g8pK6DKaGj0MzTeClUfzDQM6ew5b2NvfWOmpmlg9QLX+25mp2moUuVKDs2pRN88iWLDm2A+EggFPrNeodY8UQrExsZis9mw2Wy0bNnS+bejyZLwL46uTN0GMasr06LiGuS2Vu2b9yCnL8HZT+HQYNjbEE48D0GNoe5CuPovqPMahERk/sy5c/D228TPe5eA0EZuN2t6gdtzRcfmUNxunkTJ5U2nMaeAvyulKmqtU4ogTaIYKcnNk0ojr7syLWQJCQmMmxhDxWH/zDTYiW7ciXETn/U82Anw2aJ/cinxTQZcc4xqldNIPB3M15vqU67uowwbOcH9DvfuhTfegHffhaQkwus3ZHPiAbe9wNlPHqRV+NUFcpx57gfez2Xtm0B6qys6uQZzpVQPYB5QEbhKKdUeeEhr/UhhJ04IUbBy6sq0IINYbmbPmUdQp9vdlhAEdbote5BLPwFnl8CZ97k1YiMQCBVv5umXdvGv1+K4p3tA9p1cugTLlsE778CKFRAQALffDtHRRDdqxHvX3YRu0jlbL3AB8SuInrW0QI6zuNw8FRUJ2r7jTTH7K8BArOfkWutfgD6FmSghROGIHjuKgPjl6PTMIxo7g9hDo4okHSbINXS7zFnMbT8LyR+YYvQ9deHYw6aGeujLpjJb/U/48fcaoLIE8l9/hfHjISzMBO+ff4bnnze10hctgl69CKtfn+kxE0lfNpXzO9eTdiqB1N0bSV82lekxT3ssFcir8ObNsJ846HaZ/eRBWrUompsnUfJ5M9AKWutDKnMvR/bCSY4Q/q24d4HqtivTkwcJiF9RoEEsN+HNm/FjvIdi7hP7adXwd9hTG/QFCLwKqj8Ole+Csu3c97h2+jR8+KEZM3zrVggKgltugfvvhwEDnIOeuIocMRxb71706HstFy7ZefC+e4ietbRAz0FO/cAXZAmAEN4E80NWUbtWSgUD44AdhZssIfyTP9QxcASx2XPm8c6C6YUSxHLjDHLuirl/e5/o8eehyv1QeSSU65F9oBMAu53Op07BiBHw+edw8SK0bw+vvgp33QU1a+aajrCwMBrXrwtQKM+ui8vNkyj5vClmHws8CoQBh4EIa1oI4afCwsKYNnUyrZvUZ9rUyUUeVMLq1WP6s3eQ/uUEUnetM8Xcu9aR/uVTTP+/+wjrdQzqvA7le2UP5Hv3mmLzRo2Y8dtvsHIljB4N27ZBXByMG5drIHdtpXH06FHnEKOF0UojcsRwNq9ayv2tg0n7ejr3tw5m86qlRI4Ynqc0SksSkRNvarOfACKLIC1CiJJMp8P59ZDyOZz9nMiuB7HNC2L2ZyvY8WcorVp3I/q1793fWJw6BZ9+Ch98AN9/b4raBw5kcuXKTNm2DcqWzVNSsj72KOxSFMfN06Z1q70uASguj2aEf/CmNvt04B/ABeBroD3wmNb6g0JOmxDC32VcgHPfmgCe8iXYT4IqC+Wvh5qTCWt2K9P6VHX/2TNn4IsvTKW1lSvNQCfNmsGLL8K990L9+nxvs+U5kAtREnnzzHyA1nqCUupWTDH7HcBqQIK5ECI7+2lIWW4F8K9Bn4MyVaDiEKg4FCreAGUquv/s+fOwfLkJ4MuXm+fgV11laqePGAEdOsiQo0K44U0wD7LebwQ+0lqfUvLPJESJ5a5GflxcHLGxsZ6LfdMSIOULOPuZGdCEdAisC1XuhUpDobwNVLD7z168CN98YwL4l1+aXtrq1IGHHjIB/JproIw31XuEKL28CeZLlVI7McXsjyilQoHUwk2WKI2Ke7Ou0sJdjXybzZb5Gmg7pP5k+kA/t8L8DRDcHKo/YQJ4SFf3tdAB0tLgu+9g8WLzLDw5GWrUgMhIE8D79DGdvAghvOJNBbiJSql/AWe01nal1HnglsJPmiht/KFZV6mWfgzOfWMF8JWQcQooY4J2zReh0q0Q3NJzMfjZsyYH/sUXpgj99GmoXBluvdUE8P79TftwIUSeedtpzGmXv88B5wotRUKI4kGn07ZpMiQ+x1vPbIU91lCgAbWh4k3m2XeF6yGghudtHDliis6/+AJWrTJdrFavDjfdZIL4DTdASEjRHI8QJZhXwVwIUXLkOBhG5ACT+z63As59y2tPJcHJ37h4qYLJfVccBGXbey4+1xp27DDB+4sv4McfzfwmTeDRR02vbD17uu2RraSSwUdEUfD4H6WU6qm13qCUKqu1vliUiSqN5HmxKCqO71RsbCyfffw27Zsn0qLBX/Rs+CjsPW9WCqwHlW5n8submPLyBsaNHcqaNc+636DdDhs3Xg7ge6zhQzt3hqlTYehQCA8vtbXQffU/LDcRpUtOt8ezgE7AJqBj0SSn9JLnxaJI2JPg/Fo4v5qovt8R1e1XAC6klqFcjYFQvh9UGAhl24JSfL/NBgFVs2/nyBHT9vvrr+Hbb+HkSfO8+9pr4fHH4eabzUAnwmckaJcuOQXzNKXUfCBMKTUr60Kt9bjCS5YQokBkpJhe186vhvPfQeo2IANUCJTraYrOK/RjyOCnWfXdVx43E5SRYWqff/ONCeC/mpsAateGG2+EwYNh0CBToU0IUeRyCuZDgOuAa4GtBb1jpVQAsAVI0FoPKejtC1Eq2U/ChU1wfgNcWAsXNgPpQBCU6w41nocK/SCkG5S53HOaPcPNM/A9e7g1IQFuuokvN2y4XNu8Z0+YNs1UXmvXTtqAC1EMeAzmVp/si5RSO6wxzAva3zGjr8mtvBD5oTWk7bYCt/W6tNNaGAQhnaDGU1D+WjPyWJnyOW8vKcn0e/7NN+a1bx9/t/bzdZ063Prmm9CvH1SqVLjH5YWEhARmz5nH9n2HeOb5KUSPHSUjkIlSzZtb6pNKqc+UUseVUseUUp8opepfyU6tzw8G5l3JdoQojhISEnjm+SnOQJOQkFAwG864COc3wsmX4fBQM973vhZw9AE4+ykENYXQf8JV30PzZGi0yUxXuM59IE9OhmXL4MkneWvrVtNpy9Ch8N57psLa7NlEdu0Ke/bwarNm5jl4MQjkCxctpmv/IczfnkbQoInM355G1/5DWLhosa+TJoTPeNM+ZD7wIaZPdoC7rXnXX8F+ZwITAN//MghRgBYuWsyEmJewhw+2As1B3us/hOkxE3Md8jITrSHtoOlZLfUnuLARUreAo2FJ0NVQ4UYo39M8+w5u6bm5mENyMqxfD2vWmNe2bZCRAcHBXChXDiZNApsNuneHYNP1asLHH+fnNBSahIQEJsS8ROCQSQRZ46AHVauHbtyJCTFTsfXuJTl0USp5E8xraa3nu0zHKqUey+8OlVJDgONa661KKZu7dRITE+ncubNzesyYMYwZMya/uxSiSFxRoEk/ZgXtny4HcPsJs0wFQ9mOUC3aBO5yPSCwdu4JOnMmc/DeutUZvOnWzYwJbrPBNdfw2KBBrJns3dCcvjR7zjxzoxSYuZ93FRiMPXwQs+fM83qIUWm6JUoSb4J5olLqbuAja3okcPIK9tkTuFkpdSMQAlRWSn2gtb7bsUJoaChbtmy5gl0IUfS8DjT2ZEjdagXvzeY9/ZC1dhko29r0sBbSxXq18zxIiYPW8Oefpr234/XLL6YNuCN4P/ecCd7dukG5cjluLusz6YsXi0dXE/G7dhNQs6/bZQE1GrJj1zqvtyVBW5Qk3gTzB4DZwCuABjZa8/JFa/0M8AyAlTN/0jWQC+Gvcg00v74K+z6CS39cXhDUxOS0Q7pAua4Q0sHz8KCuLl6En382QXvTJvP+119mWYUKZqSxZ5+9HLzL51L5zYW7RwVn9h9l4aLFeXtUUAjCmzdj8/aDBFWrl22Z/eRBWoVf7YNUCeF73gy08idwcxGkRRQwR+4qftduwps38+sav8W+hzxtJ7xpNTb/ccB9oDmxj1ZXnYTgblD5bit4d865X3NXx45dDtobN8KWLSagAzRubIJ2jx7m1bZtvrtL9fSooGwxeSYdPXYU7/Ufgm7cCeVSAqLTLxEQv4LoWUt9ljYhfMmnHSRrrdcAa3yZhpLqoUce5f1Pv6Ji95EE1erLD7/t49VOvbjntht5643X87St4tAMqFj1kGc/CRfj4eLvkBoHF3+Bi78RPfAC730chm7cOXug2f410a9t9a5XtORkUzntp59M0N6yBfbvN8uCg6FTJ4iONoG7e3eoW7fADi2nRwVprQYybcYrzH5lRoHtL6/CwsKYHjORCTFTsYcPIqBGQ+wnDxIQv4LpMU/77c2qEFeq9Ix2UIokJCSw7LuNhN79H2dQCapWj3JXd2XZsqkkJCR4/aNXYLWz88FTbvzo0aN52k6+SyjsSSZoX4q3grf1srvsv0w1UzRedSxhdSKY/vwxJvzzBezhN2YJNM+432dKiikudwTtLVtg167Lyxs1Mn2cP/yw6aylY8dCHWUsp0cFQaGNmf/hi3S/potPi9sjRwzH1rsXs+fM450F03nwvnuInrVUArko1ZTW2tdpyKZz5866NFeAu9Kc5zPPT2H+9jRCmnXPtix190bubx3sVY3fhIQEuvYfQuCQSdlymunLprJ5VeH9gGY9B67TeTk/rjcjATUbYj9xkID45ZlvRuxn4NL2y8HaEbzTXdqHqwqmYlrZcAgON+9l20Bg/WwDiDhuHt5Z8L4JNI6bh3Pn4LffMgfuHTtMDXOA+vVN4Ha8OnWCmjXzeQa9czmt7/HgffdyNuUsS/4McfvdOf/HRi4lHqDCiR2Feu3zwuelNEIUIKXUVq1159zXzC7XnLlSqjbwT6Ce1nqQUqo10F1r/U5+digKX0HV+M1vM6Di8nw7x6Zik57C1vgVwqofAPuxyx9S5SC4lek1zRGwg8Mh6Krc23FbwurVY9ro++Gzj5lWNgAee8zULN+zx9Q6B6hVC7p0gWHDLgfuAiwu94a7UhfiVpCWoSnr5pn0mS2fU/PmCaQfbZSnJmBCiMLnTTF7LKaTmP+zpncBiwEJ5sVUQdX4ze9Ngc+eb2sN9kS4tAsu/cHs/7yPvfUg9zcjbYYx+8NPmDbxRghuYTpdKdsGghqBCvB+n6mpEB9vgrXj9euvcPo008C05W7SBNq3h8hI8965s3l27sMhQXO60bm0eAInFj5J+W7DCQ5tRFriQc5s+ZxKnW4msFJNdNrFPDUBE0IUPm+CeU2t9X+VUs8AaK3TlVL2Qk6XuAIFVeP3Sm4KCq3SnE6nTo0LcO5/cGkvpO29/J6214wSZonfWY2A+pFuNxNQswk7jnaDuu96t9+0NNi9G7ZvN6/4ePj9d/jjD9OWG0yTsLZt4c47oX17Hp07l9e//75YjiSWU6lLcNc7afzXWv7Y8yOp+7YSWKU2NW+eQGAlU+QvTcCEKH68CebnlFI1MG3MUUp1A5ILNVXiihRUjd/83hRccaU5ezJNw1Lg7BfOYD39b7+SsKExsz88zvnUYJ6ZcBPRd6YSVicYghqbfsnL9zHvwc0huAXhHRfk/Wbk4sXLQTs+/nLw3rUL0tPNOkqZ5mDh4XDbbSa33b49NG2aaQSx+MWLi2Ugh9xLXWqm1+Ho6T1u60tIEzAhih9vgvnjwJdAU6XUBiAUGFaoqRJXLC31AldVK8tfvywh6ew5qlaqQL3Q6qSlXvB6G/m5Kci1S9NePQmrW870O552wP17RhLvPA8kDDUbLVOFtVvtvLKkDPY2jxFQvxGbjx7gvYeWMX3yRCJHjnSb/uixYzzfjPy+nOj7XjCDiuzaZSqibd9uArkjp62UCdCtW8Mtt5j38HBo0SJPnbAUR7mVunRs34Z7Rw5j9GOPU7nnXdIETIhizqva7EqpQKAFoIA/tNZphZkoqc1evGroeqyd7cYzzz3H/B2akGY9si1L3b2e++u+wbToM5kXlKlonlUHNjTvQQ2Z/GIsU16cD8FNSThygWade2dqagde1Kq/cIGFs19nwhvzsUfcQkBoI+zH9xOw/gOmJx0l0m7ltAMCzHPt8PDLAbt1axO0c+n2NCfF7Tq68ralQvfu3bFdd4NX194XivM5FiKvCrs2+21ZZjVXSiUDv2mtj+dnp8K/hIWFMW3qZOI2/49pz4+EtD8geRWkHTZ9ijve0w8T/7OdgPrT3W4noEZjdhxuCrXucQZtghqZttpWZbCEhARmvz6Pd5akcqnCCqLHjmL2W+9QsfvITEEHXGrV/2Ma0wZeB/v2wd69Jne9axf8+SeRWmMDZifsY0elKrSqVoXoG/oS1qEDNG9uXo0bO0cJKy28LXUpW7Ys06ZOZtO61VJ7XYhizJti9geB7sBqa9oG/IAJ6i9ord8vpLSJopKRAulHIP0v6+X+7xWvppjxs10F1DBtrQMbQLnuhLfaz+b9Hro0PXmQVu1ugerj3SbD07P2OlUqEhRxr9vPBNRoyI73JsAcq1e7KlVMgO7Vyxmsw5o3Z1qzZsViLO7iRDpfEaLk8CaYZwCttNbHwNnu/E3gGmAtIMG8OMq4YIbVtB9zeT/qZt4RyDib/fOqHATWhcB6UDYCKgxizjtfMvbRKSZ4B9U372UyF0NHj0+wnlO76dLUXcUpux2OHCFh61YmPDOFwNtfzPasfffchyhfd5/7G4TEA7QaPAgmPGGeb1erdsWnrjRxlLpIzlsI/+ZNMG/kCOSW40BzrfUppVShPjsXLjIumvGt7YnW6wSkJ16eTk/MHKTdBWiAMlXNWNgBtaFse6gwyARsR+B2/F2mSrZ20ItWxjH22XtyTKbb4tsTBwn45QumD+hD2Ny5cPDg5dfhw5CWxuygEOw3Pua2qVTZHsM5s2EhZa/umv0GYcfXRK9a6l2f50IIUUJ5E8zXKaWWAR9b08OseRWApMJKWImWcdEM1uF4ZZzKND3hnj/g0E0ugTvRc3BGmaLugFATpEM6QWCdywE703stKFO24I7jzBkTjLO8Ig8fxqYuMPt/b7IjPYNWaReJTkslbE68uUGoVw8aNjRDczZsCFddRfx/PyOgdhO3uwm+qi2Xtn1G+jIZXEMIIdzxJpg/CtwG9MLUZl+gtV5iLetXWAkr9nQGZJwB+2nIOG3enX8nucw7CfbMwRp93vN2VVm6tNaQftgE6OCmEFDT/B0QCoGhmacDquWtxzJvpKTAkSPZXs/s3AnXXw8JCSZwn3Vzg1GrFtSvb55TX3stc1esYExMDDRoYAJ3/fpuK5uF/3U8x6ZSVStXYOOqpfJ8Vwgh3PBmPHMNfGK9UEr1Ukq9rrV+tLATV2yc/x5OTM0SrJOw+tHxINAE2oDqlyuJlW1v5aJrXJ6fdbpMee4ojOY2aWmQmGjGxT5+3LwfO+Y2aLsN0sHBtFfKBPrWrWHAABOYw8LMe/36JsddtuzlvtmTkzkaHMyH75pe1nLqmz23DmrqhVaX57tCCOGBV0OgKqUigJHAcGA/8Gkhpqn40RkmNx1YBwJamaZUAdarTDUIqOryt/VSFQq3722tzbjXJ06YIO0uULv+feqU++2UL28G+KhXDyIiYNAgM531Vb06I/r18+omIz8DquTWVOrtOW/maXtCCFGaeAzmSqnmwAhMED+JGVxFaa1LX9F6hX5QYWPhbV9rMzzmqVNw8iSdTp+GRYtMgHYE6xMnsv/t6F40qypVTHF37domF22zmb9r17483/F3pUo+HfDDVU5NpfwlmGcdMc5mswFFP2KcEKJ0ySlnvhNYB9yktd4DoJRy30BYGBkZplJYUhKcPm3erQCd6d3dvEuXnJv5N4Cji1KloHp1M651aChcfbWpOBYaenlezZrm5QjQISE+OPiC4e9F6SUhaLu7IYmLiyM2Ntbvj02IkiqnYH47Jme+Win1NbAIUwGu9Nm7F7755nKA9vSenGwCuichIVCjhgnO1aub7kKrV788z3ofN2UKsz76yATo6tVNd6NCFBF3NyQ2m00CuRDFmMdgrrX+DPjMaoI2FBgP1FZKvQl8prVeWTRJ9B1HDqVvYiJTtm8H4GKZMtgrV6Z8vXpQtap5ntyqlemspFo1M8/13RG4a9Twup/vX1991WxTCJGNPMoQIjtvarOfAxYCC5VS1YE7gIlAiQ/mzh+H8+e57brr+PS77yjrx0XY/sz1R1t+wEs3ueZCZOdVbXYHrfUp4C3rVXqUL8+p4GC/fhbt7+rUqSOjYwkhhAdlfJ0AIXKSkJDAM89PYfu+Qzzz/BQSEhJ8naQSIzY2FpvNhs1mc5Z22Gw2ZxG2EMJ/5ClnLkRR8jSK2vSYiUSOGO7r5Pk9Ka4WouSQYC6KpYsXLzIh5iUCh0zKNorahJip2Hr3kq5chRDCIsXsolj6K/EU9vDBmbp2BTOKmj18ELPnzPNRyoQQoviRYC6KpQuX0gmo2dDtsoAaDdmxa08Rp0gIIYovCeaiWCoXHIj9xEG3y+wnD9KqxdVFnCIhhCi+JJiLYqleaHUC4pej0y9lmu8YRS36oVE+SpkQQhQ/EsxzIU2jfKNs2bJMj5lI+rKppO7eSNqpBFJ3byR92VSmxzwtld+EEMKFBPMcLFy0mK79hzB/e5rVNCqNrv2HsHDRYl8nrVSIHDGczauWcn/rYNK+ns79rYPZvGqpNEsTQogspGmaBwkJCdI0qhjw91HUhBCiKEjO3IPZc+ZJ0yiklzAhhPAHkjP3IH7XbgJq9nW7zDSNWlfEKfIN6SVMCCGKP8mZexDevJk0jRJCCOEXJJh7ED12lDSNEkII4RckmHsQFhZW5E2jpBmcEEKI/JBgnoOibBolzeCEEELkl1SAy0VRNI2SZnBCCCGuhOTMiwFpBieKI2mWKIT/kJx5MSDN4ERxJM0ShfAfEsyLgfDmzdi8/SBB1eplW2Y/eZBW4aWjGVxsbKwz1+fICQISUIQQIhcSzIuB6LGjeK//EHTjTpmK2p3N4GYt9WHqik5OOUEp2hVCCM/kmXkx4ItmcEIIIUoOCebFhIwQJoQQIr+kmL0YkRHChBBC5EeR58yVUg2UUquVUjuUUvFKqb8XdRq8JU1zhBBC+ANf5MzTgSe01tuUUpWArUqpb7XW232QlhxJ0xwhhBD+oMhz5lrrI1rrbdbfZ4EdgNTwEkIIIfLJp8/MlVKNgA7Aj67zExMT6dy5s3N6zJgxjBkzpmgTJ4QQQvgJnwVzpVRF4BPgMa31GddloaGhbNmyxTcJE0IIIfyMT5qmKaWCMIF8odb6U1+kQQghhCgpfFGbXQHvADu01v8p6v0L/yGtCYQQwju+KGbvCdwD/KaUirPmPau1/soHaRHFmLQmEEII7xR5MNdarwdUUe9XCCGEKKmkO1chhBDCz0kwF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz0kwLyZiY2Ox2WzYbDaOHj3q/Ds2NtbXSRNCCFHMKa21r9OQTefOnfWWLVt8nQwhhBCiyCiltmqtO+fns5IzF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz0kwF0IIIfycBHMhhBDCz/kkmCulblBK/aGU2qOUmph1eWJioi+SVWzMnTvX10nwKTl+Of7SqjQfO8jxAzXz+8EiD+ZKqQDgdWAQ0BoYqZRq7brOiRMnijpZxUpp/0LL8cvxl1al+dhBjh8Ize8HfZEz7wrs0Vrv01pfAhYBt/ggHUIIIUSJEOiDfYYBh1ymDwPXuK5w/vz5i0opu8usRKA0ZddrKqVK0/FmJccvx19aj780HzvI8bfI7wd9EcyVm3mZ+pTVWocUUVqEEEIIv+eLYvbDQAOX6frAXz5IhxBCCFEi+CKY/wQ0U0o1VkoFAyOAL32QDiGEEKJEKPJgrrVOB6KBb4ADQD3gC3dN1JQxy2rC9qtSqmPRprZw5dZETyllU0olK6XirNckX6SzMCil3lVKHVdK/e5heUm/9rkdf0m+9g2UUquVUjuUUvFKqb+7WafEXn8vj78kX/8QpdRmpdQv1vFPcbNOSb7+3hx/3q+/1tonLyAA2As0AYKBX4DWWda5EViBec7eDfjRV+n10fHbgGW+TmshHX8foCPwu4flJfbae3n8Jfna1wU6Wn9XAnaVsv99b46/JF9/BVS0/g4CfgS6laLr783x5/n6+7IHOG+aqN0CvKeNH4CqSqm6RZ3QQlKqm+hprdcCp3JYpSRfe2+Ov8TSWh/RWm+z/j4L7MC0cnFVYq+/l8dfYlnXNMWaDLJeOstqJfn6e3P8eebLYO6uiVrWL7Q36/grb4+tu1Ucs0IpFV40SSsWSvK191aJv/ZKqUZAB0zuxFWpuP45HD+U4OuvlApQSsUBx4Fvtdal6vp7cfyQx+vvy2CeaxM1L9fxV94c2zagoda6PfAa8HlhJ6oYKcnX3hsl/torpSoCnwCPaa3PZF3s5iMl6vrncvwl+vprre1a6whMa6auSqk2WVYp0dffi+PP8/X3ZTD3polaSW7Gluuxaa3POIpjtNZfAUFKqXz33etnSvK1z1VJv/ZKqSBMIFuotf7UzSol+vrndvwl/fo7aK2TgDXADVkWlejr7+Dp+PNz/X0ZzL1povYlcK9Vs7EbkKy1PlLUCS0kuR6/UqqOUkpZf3fFXK+TRZ5S3yjJ1z5XJfnaW8f1DrBDa/0fD6uV2OvvzfGX8OsfqpSqav1dDrgO2JlltZJ8/XM9/vxcf1/0AAeYJmpKKUcTtQDgXa11vFJqrLV8DvAVplbjHuA8cL+v0lvQvDz+YcDDSql04AIwQltVHf2dUuojTI3Nmkqpw8BkTEWQEn/twavjL7HXHugJ3AP8Zj03BHgWuApKxfX35vhL8vWvCyxQZtCtMsB/tdbLSstvP94df56vvyo53w8hhBCidPJlMbsQQgghCoAEcyGEEMLPSTAXQggh/JwEcyGEEMLPSTAXQggh/JwEcyHySCllt0Yyire6W3xcKZXv/yWl1LMufzdSHkZSy/KZsUqpe/O7z7ywRnBaZv19s3Izwl8h7K+Hh2UxSqkEpdQLedzmQqXUKaXUsIJJpRDFi8/amQvhxy5YXTGilKoFfAhUwbQVz49ngX/m5QNWW9Qip7X+kuydOxU0G5ACbPSw/BWt9Yy8bFBrHamUir3CdAlRbEnOXIgroLU+DowBoq3eqgKUUi8rpX5SZhzmh8CZ21yrlPpMKbVdKTVHKVVGKfUSUM7K6S+0NhuglHrbyvmvtHqJysTKoT5p/b1GKfUvZcZI3qWU6u1m/brW/uOUUr871lFKDVBKbVJKbVNKfaxMf+EopW5QSu1USq0HbnPZTpRSarb1d6xS6k1lxubep5Tqq8w47TtcA2cO+ziglJpizf9NKdVSmYFHxgLjrbRmOxY352GBdZ4OKKVuU0pNt7b3tTLdpgpR4kkwF+IKaa33Yf6XagEPYrqe7AJ0AUYrpRpbq3YFngDaAk2B27TWE7Fy+lrrSGu9ZsDrWutwIAm43YtkBGqtuwKP4b6E4C7gG6tEoT0Qp0xfz88B12mtOwJbgMeVUiHA28BNQG+gTg77rQZcC4wHlgKvAOFAW6VUhKd9uHz+hDX/TeBJrfUBYA4m9x2htV7nxbE3BQZjhs38AFittW6L6TlrsBefF8LvSTG7EAXDMcrTAKCdy7PZKpjgfAnYbAV+R3euvYAlbra1X2sdZ/29FWjkxf4dg3V4Wv8n4F0rp/q51jpOKdUXaA1sUKYb6GBgE9DSSsNuK60fYEof3FmqtdZKqd+AY1rr36zPxFvpqO9hH+7SfRv5s0JrnWalIQD42pr/G96dOyH8ngRzIa6QUqoJYMeMTayAv2mtv8myjo3sQzh66kv5osvfdiBbMXsOn7Hj5v9aa71WKdUHk1N9Xyn1MnAaM5byyCxpjcghbZ72m5El3RlWOuzu9uFtuvOSBq11hlIqzaUPa0cahCjxpJhdiCuglArFFAvPtoLIN5gBEoKs5c2VUhWs1bsqM0peGWA4sN6an1bYz3aVUg2B41rrtzEjdnUEfgB6KqWuttYpr5RqjhnBqbFSqqn1cU+B2Bue9pGTs0ClK9inEKWOBHMh8s5RYS0e+B+wEphiLZsHbAe2KdPE7C0u5w43AS8BvwP7gc+s+XOBX10qwBUGG+Y5+c+YZ/Cvaq0TgSjgI6XUr5jA21JrnYopVl9uVYA7mN+detpHLh9bCtzqTQU4IYQho6YJUQSsYvYntdZDfJwUv6aUigFS8to0zfpsLLBMa+2unoIQfk1y5kIIf5ICjFH56DQG6AukFkqqhPAxyZkLIYQQfk5y5kIIIYSfk2AuhBBC+DkJ5kIIIYSfk2AuhBBC+DkJ5kIIIYSfk2AuhBBC+Ln/B6WSdyKHdbPvAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 576x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 5))\n",
"\n",
"params = {'elinewidth': 0.75,\n",
" 'capsize': 3,\n",
" 'markersize': 15,\n",
" 'markerfacecolor': 'C0',\n",
" 'markeredgewidth': 0.75,\n",
" }\n",
"\n",
"ax.errorbar(x, y, yerr=e, fmt='k.', label='Data with errors', **params)\n",
"ax.plot(t, y_hat_fit, 'r', label='Nonlinear least squares')\n",
"ax.plot(t, y_hat_fit_weighted, 'gold', label='Weighted nonlinear least squares')\n",
"ax.tick_params(direction=\"in\")\n",
"ax.set_title('Comparison of unweighted and weighted fit', fontweight='bold')\n",
"ax.set_xlim(0, 3.5)\n",
"ax.set_xlabel('Depth in sediment [m]')\n",
"ax.set_ylim(0, 10)\n",
"ax.set_ylabel('Age of sediment [ka]')\n",
"ax.legend()\n",
"\n",
"plt.savefig('Figure.svg')\n",
"plt.savefig('Figure.png', dpi=300)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "c8c7d1fd",
"metadata": {},
"source": [
"## Another way, using `least_squares`\n",
"\n",
"I'm not smart enough to know why this yields a different result from the unweighted least squares solution we get from `curve_fit()`, above."
]
},
{
"cell_type": "code",
"execution_count": 103,
"id": "fcf2e633",
"metadata": {},
"outputs": [],
"source": [
"def residual(phi, y, t):\n",
" return y - model(t, phi[0], phi[1])\n",
"\n",
"solution = so.least_squares(residual, x0=(0, 0), args=(y, t), method='lm')\n",
"\n",
"c0, c1 = solution.x\n",
"\n",
"y_hat_lsq = model(t, c0, c1)"
]
},
{
"cell_type": "code",
"execution_count": 106,
"id": "5d6a98d7",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfMAAAFJCAYAAACPXsRYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABmgUlEQVR4nO3dd1xV9f/A8ddHpnuBCwdqThykuBemubVlpVmGZmZmllZmyzQr+6mlmZaVKVqWlt+W2rBMcqaBooZ7oeICFRQFZHx+f5zDFfACFwQuF97Px+M+Lmfc8/mccy73fc7nfIbSWiOEEEIIx1XC3hkQQgghxO2RYC6EEEI4OAnmQgghhIOTYC6EEEI4OAnmQgghhIOTYC6EEEI4OAnmotBRSgUqpbRSaqq985IZpVRzpVSwUuqGmddm9s5TdpRSAWZeg3LwmQI7F0qpE2Za9+Z3WmnSTD0moQWUXpCZXkAOPqPNl3f+5QyUUv5mOtFZrFNOKfWTUuqque44R/h/LQ4kmBcRSqlOSqnVSqmLSql4pdRRpdRHSilXe+ctF9YBHwL/2DsjWfg/oDUQjJHXKPtmxyb7MPK6Ki83miYI++fldouoVRjnYF9ebTBNED6RV9vMwhhgEHARmAeEkuH/NTcXjeL2Ods7A+L2KaWGAF8BTsBu4F/AG+Mf7w3ght0yl0NKKWet9dfA1/bOSzYamu+va63/smtObKS13gHssHc+ijOt9Xx75+E2pX7vl2mtp6SZX9j/X4s+rbW8HPgFlMK4StbAl0CJNMvqA67m3y2A3zDuICOB1UCjNOueMLcxDQgDYoG5QFMgBLgKfAO4mesHmOtvwrhCvwIcA4al2eaLwGHgGpCAcaExOM3yQHMbnwJ/YFx0+KeZP9Vcr5WZzhUzX/8BT6fZzn0YFzBXgXBgAVDBXOZtbksDI4GTwGVgTjbHNattnkizTW38G1ndRup6/hmOWZA57W9OnwBeBS6Yr5fM5Z7m8hiMUrRKQApwzlxey1x+AVDmvJHmcY41j/2rgLO19M15Y4FTGN+LSWnyfG+Gc/QR8D1wHdgD+GZ2LIAAc9kgjIuHK+YxfB8oZSXtSOCljGlbOZ53A7vM45FobnNamuWp+7cZmANEAxGk/07WwLiTvIbxnZpmfiY0kzRfMpfPN6cnmtP/Z05PMadnmdOVMb7PJzC+O1uALmm2F5ThGJUHVprHaE+a7Uen+UzqcR0F7De3+xXgys3v0C3fR4zfhveAI+b+7kx7bDOkvdta2hmORaCVtPzTzJ+a5hykfZ2w9+9kcXjZPQPyus0TaPzApf7TNMpkneoYAUwDa8wfMw2cBSqa66T+kF4FlmIEVg1cApZj/NhrYJS5fuo/bQqwDaP4MHW6hbnOAowAsAD4FkgC4gFvc3naH4cgYDFG4Lb8OJjrbTanVwGfY/wILzKX9TWXJZj5/s+c/s1c7p0mjZPmviSb0z0yOV7ZbXOK+QOYmqe5mWwn9Zj6ZzhmQea0f5q87TOPVeoxbGCuE2bOawEMSLP+HcCQ1DyY6z6VZj8DgQPm9JvZpJ+CcWe1N82xuTfDOUoBfkiTn02ZHQugLdDbnBeJcZEZYk4vyZB2snlOdmdM28rxDMC46PsU47tyyVx/SIb90xgXEanf81ignLlOkDlvP7AM4/uYVTBvYy7faU6nfs83m9O/mdMDMC64Ur+rG839ugrEYf5vcmswX2pOnwaWpDmW0WnykLpPkeb5iDOnn8D4HqTm6Yp5/Oean/vGnB9i7mukeR79M6R9yjyeMRnTznAsHsH4nmqMIvW5ZvqB3AzmbdMc99PmOlPs/TtZHF52z4C8bvMEwrA0/+zumawzyVy+Ic28Xea80eb0CXP6dXM69UfnW3P6fXN6gTkdYE6fB1zMeT+Q/i6lNDAcmI5xp3TWXP6IuTz1R+DvDPm1/DiY09u5eWfdDHABnMxlv5A+YHlg3LVpjCJB7zTHp425zt/m9IuZHK8st5nhePlncW7SrUPmwTQJqGbOCzfnDTanPzanxwAzgOMYFxkBGHfLGnjWXDc10KYG1a/M6XOZpL+I9AHWM81+3pvhXKw1p7ub07GZ7ac5b605b52Zl9T9SMG4Y0xN+wtz/UrcvIC8N5PjWQLoB7yO8X3611z/swz7dxFwx/ieJJnz/ICaab4LtTJ8r0MzSdMJI0gmAWUwvsMHMS4CSmLc/ScDFbgZ+C1BFeNuWAPvZfi/CjC3nWBOdzOXTyDzYP6gOZ0ahFNLC/zJcAfMzVKdZIzvydw0aa/IkHYX8zPPZkzbyvFI/T5MzeL/NfU8BNn797E4veSZueO7kObvOhg/NBl5m+/708w7APian0krdZ1o8z11e1fN99IZ1j+qtU5Ms02AmmbFu38wgm9Gnhmmt1pZJ62JGMFgEaAw7rSmYPyge6fNt9Y6SikVBVTD2LfDabazy3yPNt/LZJJedts8lE1+M+OUyfxzWutzafJWO03egoCngY5mvjYATYBOGKUYqeukzfcDGbZfVSllbV+9zPfU/YxMs58ZZTx2Gb8HGaXm5W7zlUoB9dKkfdBM+5JS6mImaaf6BBhtZX7G79N+rXU8gFLqGlAO43imphmntT5l/p3ludRaJyulNmOU1gwx8/cmRulAAEZR9U6tdXSa2uZlgecybOoOK5v3wCgqh5v/d1lVjLP1+ws3j38JYJyVvKRNO/V/PLffa1EISG12x7cVowgd4HWllOWcKqXqKKVcMO6cABqn+Vwj8z08w/aSs5nOqL6ZRtrtn8Z41t7M/HwDjO9a6g+VyrCNhGzSCNZatwQqYtyFuADvKaWcybBvSqnKGD9UkGHftNZJqX9mk57N28zGNfO9nPmeWfO1pDR/Z8zb3+Z7V4w7vy3mqyfQEuMu9L8M+R6ktVapL6Ce1jrWSroR5nsDAKWUBzf3M7M8Wjt2qd+RtL8nqXkZnyEv9bXW/6VJu5GZdiWM581Zedh8D8C4MPrEnM74fcrseKamWVIpVcv8uyHZSz0HL2CUXHyF8djphQzLT5jvZzBKyVL3uRS3BlTMbaRWTm1gvje2sl6qzM5BVsf/BuCZJi+uGPVB0qad+ltgy7GwhbX8iHwmB9vBaa2vYRSPpQCPAjuVUp8ppdZgXGmXxvjxiQG6K6V+Vkr9BtyJUUR+u82UPIC/lVKrgHsxfmhSn7GnYPzofoDxrLNBJtvIzmql1HpgFkalKTeMkoJkjOfxAK8qpQIx7lKdgT+01rm908irbabeSb2tlPrQzHuOaK3PY5R41MEoOk4N5t4Yx3ajNss2gdSa0l+ZbX+XKaX2YTyLteZL832EUmo58Be5+01Ivct9Syk11wyUqXmZqZT6Rim1WCkVDKw356fWfg4w095A9q1rzpvv4828B+Qkk1rr0xjPsgHWKaWWYT3IZhRkvjfGuAu/jnERXT/D8hCM+iM1gH+VUguVUj9iBPc+VvKTjPFcG+AbpdRi4K0c7FKq1ONfUym1SCn1stY6EqOeiiuw3czLd+a6T5hprzQ/97WZ9ju5SDur/LRWSn2slHoyj7YrsiDBvAjQWi/HeJb5C0YR7eMYRbGfA9e11mfM5eswimf9MJ5pdtdaX7rN5Ldg/LDdjXHX+rjWOtT84XwW4we4G8YPXXbF6ZkJwviBHAb0x3hW+rA2rAUewnhePBij2PNTbt7F5VgebvN1jB/3ehhF4rltlpR653cRo0h0i5VlAAsxajwfx8h3P4yLqkXWNqq1/ht4BuM5cB+Mi7DUgJldaUlaUzFqTHfAKF6uqrX+FeMOcLeZj/sxLu4+NNP+CyMon8Uowv4fRsW9rIzCuLBpilGU/WkO8phqGPAnxsVRQ4wLzeyEYDzagZvHPvVdY1TIRGudAtyDcR7KYVxs3Inxf5lZnwnPAd9hlDr5YfRfADk4/lrrE8BsjAv2J4DHzEVPYNRmTzHz0gnj+/ibuXw8xsV8BYyKazNsTTMbGzEu1pIxHhHdk0fbFVlQNy/qhbCd2YPVEozKa/72zY3ILaVUea11jPl3TYwLshLAHVrro3bNXDGglCqLUZlQm9OvAO9i1JbvYtfMCYciFeCEKN52KaV+wbjrH4IRyH+RQF5gemDUdfkVo87ACHP+PPtlSTiifCtmN5+RXVBK/ZdmXiWl1B9KqcPme8X8Sl8IYZOdGEH8ZYyKZLMx2hOLgnESo+7DCxjHfTfwkNb6O7vmSjicfCtmV0p1xXjOtExr3cycNxO4pLV+Tyk1GaPDkpfzJQNCCCFEMZGvz8zNdpdr0gTzgxgdS5xVSlXH6FSgUVbbEEIIIUTWCro2e1Wt9VkA871KAacvhBBCFDmFsgJcmTJldNoSAw8PDzw9M3byJIQQQhQdISEhUVrrXAW7gg7m55VS1dMUs1+wtlLjxo0JDg4u4KwJIYQQ9qOUykkPk+kUdDH7zxgdmmC+/1TA6QshhBBFTn42TfsGo7ehRkqp00qp1N6I7lZKHcboMey9/EpfCCGEKC7yrZhdaz00k0U98itNIYQQojgqlBXgrElMTOT06dPEx8fbOytCFGru7u7UrFkTFxeX7FcWQhQJDhPMT58+TdmyZfH29kapjCMeCiEAtNZcvHiR06dPU7duXXtnRwhRQBxm1LT4+HgqV64sgVyILCilqFy5spRgCVHMOEwwBySQC2ED+T8RovhxqGBub0opXnjhBcv07NmzmTp1aq62deLECZo1awZAcHAw48ePz4ss2mzq1KnMnj07z7Z34sQJvv766zzbnhBCCNtJMM8BNzc3vv/+e6KiovJ0u35+fsybl78jHiYnJ+fr9gtLMNdak5KSYu9sCCFEgZJgngPOzs6MHj2aOXPm3LIsPDycHj160KJFC3r06MHJkycBCAgIYPz48XTs2JF69eqxatWqWz4bFBTEgAEDAOOOeeTIkfj7+1OvXr10Qf6rr76ibdu2+Pr68tRTT1kC9NNPP42fnx8+Pj68+eablvW9vb1566236Ny5M999l/mIikePHqVPnz60bt2aLl26cODAAQBWr15Nu3btuPPOO+nZsyfnz58H4O+//8bX1xdfX1/uvPNOrl69yuTJk9m0aRO+vr63HJ+zZ8/StWtXfH19adasGZs2bQJgyZIlNGzYkG7duvHkk08ybtw4yzFLe5zKlCkDQGxsLD169KBVq1Y0b96cn34y+hw6ceIETZo0YezYsbRq1YpTp04xa9Ys2rRpQ4sWLSzH5Nq1a/Tv35+WLVvSrFkzVq5cmekxEUIIR+IwtdnTef55CA3N2236+sLcudmu9swzz9CiRQsmTZqUbv64ceMYPnw4jz/+OIsXL2b8+PH8+OOPgBHMNm/ezIEDBxg0aBCDBw/OMo0DBw6wYcMGrl69SqNGjXj66ac5cuQIK1euZMuWLbi4uDB27FiWL1/O8OHDeeedd6hUqRLJycn06NGDPXv20KJFC8BoprR58+Ys0xs9ejQLFy6kQYMGbN++nbFjx/LXX3/RuXNn/vnnH5RSLFq0iJkzZ/L+++8ze/ZsFixYQKdOnYiNjcXd3Z333nuP2bNns2bNmlu2//XXX9O7d29ee+01kpOTuX79OmfPnuXNN98kJCSE8uXL0717d+68884s8+nu7s4PP/xAuXLliIqKon379gwaNAiAgwcPsmTJEj7++GPWrVvH4cOH2bFjB1prBg0axMaNG4mMjKRGjRqsXbsWgJiYmCzTE0IIR+GYwdyOypUrx/Dhw5k3bx4lS5a0zN+2bRvff/89AI899li6YH/vvfdSokQJmjZtarm7zUr//v1xc3PDzc2NKlWqcP78edavX09ISAht2rQBIC4ujipVjEHnvv32Wz777DOSkpI4e/Ys+/btswTzhx9+OMu0YmNj2bp1Kw8++KBlXkJCAmA0B3z44Yc5e/YsN27csDR16tSpExMnTmTYsGHcf//91KxZM8s02rRpw8iRI0lMTOTee+/F19eX9evX4+/vbxlA5+GHH+bQoUNZbkdrzauvvsrGjRspUaIEERERluNZp04d2rdvD8C6detYt26d5eIgNjaWw4cP06VLF1588UVefvllBgwYQJcuXbJMTwghHIVjBnMb7qDz0/PPP0+rVq0YMWJEpuukrVHs5uZm+duW8ePTru/k5ERSUhJaax5//HFmzJiRbt3jx48ze/Zs/v33XypWrEhAQEC6ZkmlS5fOMq2UlBQqVKhAqJWSjmeffZaJEycyaNAggoKCLJX9Jk+eTP/+/fnll19o3749f/75Z5ZpdO3alY0bN7J27Voee+wxXnrpJcqVK5dprWtnZ2fLc2+tNTdu3ABg+fLlREZGEhISgouLC97e3pZ9TbufWmteeeUVnnrqqVu2HRISwi+//MIrr7xCr169mDJlSpZ5F0IIRyDPzHOhUqVKPPTQQ3zxxReWeR07dmTFihWAEXQ6d+6cp2n26NGDVatWceGCMdDcpUuXCA8P58qVK5QuXZry5ctz/vx5fv311xxtt1y5ctStW9fyTF1rze7duwGjGNrLywuApUuXWj5z9OhRmjdvzssvv4yfnx8HDhygbNmyXL161Woa4eHhVKlShSeffJInnniCnTt30q5dO4KCgrh48SKJiYnpnul7e3sTEhICwE8//URiYqIlP1WqVMHFxYUNGzYQHm59gKHevXuzePFiYmNjAYiIiODChQucOXOGUqVK8eijj/Liiy+yc+fOHB0rIYQorBzzzrwQeOGFF5g/f75let68eYwcOZJZs2bh6enJkiVL8jS9pk2b8vbbb9OrVy9SUlJwcXFhwYIFtG/fnjvvvBMfHx/q1atHp06dcrzt5cuX8/TTT/P222+TmJjIkCFDaNmyJVOnTuXBBx/Ey8uL9u3bc/z4cQDmzp3Lhg0bcHJyomnTpvTt25cSJUrg7OxMy5YtCQgIYMKECZbtBwUFMWvWLFxcXChTpgzLli2jevXqTJ06lQ4dOlC9enVatWplqdD35JNPcs8999C2bVt69OhhueseNmwYAwcOxM/PD19fXxo3bmx1f3r16sX+/fvp0KEDYFSg++qrrzhy5AgvvfQSJUqUwMXFhU8++STHx0oIIQojZUuxb0Hz8/PTGccz379/P02aNLFTjkR+CwwMJDg4ON0Fksg9+X8RwvEopUK01n65+awUswshhBAOTorZRaEQEBBAQECAvbMhhBAOSe7MhRBCCAcnwVwIIYRwcBLMhRBCCAcnwVwIIYRwcBLMbTRhwgTmpul5rnfv3owaNcoy/cILL/DBBx9k+vkpU6Zk21NaZsOSRkdH8/HHH+c4z3k9zKk1/v7+pDYj7NevH9HR0fmaXlpph5HNK3PnzuX69et5uk0hhMhvEsxt1LFjR7Zu3QoYXaBGRUURFhZmWb5169YsO2x566236NmzZ67Szm0wL2i//PILFSpUyLftJyUl5du2UxWWYF4Q+yqEKDokmNuoU6dOlmAeFhZGs2bNKFu2LJcvXyYhIYH9+/dz5513EhISQrdu3WjdujW9e/fm7NmzQPphPX/55RcaN25M586dGT9+vGX4U4B9+/bdMvzp5MmTOXr0KL6+vrz00ksAVof4BHjnnXdo1KgRPXv25ODBg1b3JbNhWbXWvPTSSzRr1ozmzZtbhggNCgrC39+fwYMH07hxY4YNG2a1j3lvb2+ioqIsQ5I++eST+Pj40KtXL+Li4oCcD7c6depURo8eTa9evRg+fHim5yc5OZmXXnrJckw+/fRTIPNhU60Nhzpv3jzOnDlD9+7d6d69+y1pTJ48maZNm9KiRQtefPFFwOgbv0OHDrRp04Y33njDMlxr2mFtwRhVLzAwEDAu7Nq0aUOzZs0YPXq05Vj6+/vz6quv0q1bNz788MNMv0vz5s2z5GPIkCGZHhMhRDGitS50r9atW+uM9u3bd3Pi3HNan+iWt69zz92SZkZ16tTR4eHheuHChfqTTz7Rr7/+ul67dq3evHmz7tKli75x44bu0KGDvnDhgtZa6xUrVugRI0ZorbV+/PHH9Xfffafj4uJ0zZo19bFjx7TWWg8ZMkT3799fa631m2++qTt06KDj4+N1ZGSkrlSpkr5x44Y+fvy49vHxseTj999/108++aROSUnRycnJun///vrvv//WwcHBulmzZvratWs6JiZG169fX8+aNeuW/Xj88cf14MGDdXJysg4LC9P169fXWmu9atUq3bNnT52UlKTPnTuna9Wqpc+cOaM3bNigy5Urp0+dOqWTk5N1+/bt9aZNm7TWWnfr1k3/+++/luMTGRmpjx8/rp2cnPSuXbu01lo/+OCD+ssvv9Raa33XXXfpQ4cOaa21/ueff3T37t211lpfunRJp6SkaK21/vzzz/XEiRMtx6RVq1b6+vXrt+xH2uPy6aef6unTp2uttY6Pj9etW7fWx44d04mJiTomJkZrrXVkZKSuX7++TklJ0atWrdKjRo2ybCs6OjrdPmR08eJF3bBhQ0seL1++rLXWeuDAgXrp0qVaa63nz5+vS5curbXWesOGDZbzqrXWzzzzjF6yZIllW6keffRR/fPPP1uO5dNPP6211ll+l6pXr67j4+PT5SOjdP8vQgiHAATrXMZN6TQmB1Lvzrdu3crEiROJiIhg69atlC9fno4dO3Lw4EH+++8/7r77bsC4W6xevXq6bRw4cIB69epZhhMdOnQon332mWW5teFPM8psiM+rV69y3333UapUKQDLWN/WWBuWdfPmzQwdOhQnJyeqVq1Kt27d+PfffylXrhxt27a1DHXq6+vLiRMnshxMpm7duvj6+gLQunVrTpw4kavhVlP3I+1ws9asW7eOPXv2WEoZYmJiOHz4MDVr1rQ6bGrz5s1zNBxquXLlcHd3Z9SoUfTv399y171lyxb+97//AcbQty+//HKW2wHYsGEDM2fO5Pr161y6dAkfHx8GDhwI3ByyNqvvUosWLRg2bBj33nsv9957b7bpCSGKPscM5lXn2iXZ1Ofme/fupVmzZtSqVYv333+fcuXKMXLkSLTW+Pj4sG3btky3obPpC9/a8KfWtmFtiM+5c+dmOqxoVumk5imrvNmSr6zWj4uLy9Vwq5D9MK6pef/oo4/o3bt3uvmBgYFWh01t2LBhjoZDdXZ2ZseOHaxfv54VK1Ywf/58/vrrLwCrxzztMK6AZajW+Ph4xo4dS3BwMLVq1WLq1KlWh6zN6ru0du1aNm7cyM8//8z06dMJCwvD2dkx/5WFEHlDnpnnQKdOnVizZg2VKlXCycmJSpUqER0dzbZt2+jQoQONGjUiMjLS8gOcmJiYrpIcQOPGjTl27BgnTpwAsDyXzkrG4UUzG+Kza9eu/PDDD8TFxXH16lVWr16do/3r2rUrK1euJDk5mcjISDZu3Ejbtm1ztI2s5Ga4VVv17t2bTz75xDJc6qFDh7h27Vqmw6ZmNhxqZkO5xsbGEhMTQ79+/Zg7d67lgqRTp07phr5NVadOHfbt20dCQgIxMTGsX78euBnUPTw8iI2NtZQkZJTZdyklJYVTp07RvXt3Zs6cSXR0tOV7IIQovuRyPgeaN29OVFQUjzzySLp5sbGxeHh4ALBq1SrGjx9PTEwMSUlJPP/88/j4+FjWL1myJB9//DF9+vTBw8PDpmBZuXJlOnXqRLNmzejbty+zZs2yOsRnq1atePjhh/H19aVOnTrZFh1ndN9997Ft2zZatmyJUoqZM2dSrVo1SyW1vJDT4VZtNWrUKE6cOEGrVq3QWuPp6cmPP/6Y6bCpe/futToc6ujRo+nbty/Vq1dnw4YNlu1fvXqVe+65h/j4eLTWzJkzB4APP/yQRx55hA8//JAHHnjAsn6tWrV46KGHaNGiBQ0aNLA8EqlQoQJPPvkkzZs3x9vbmzZt2ljdH1dXV6vfpYYNG/Loo48SExOD1poJEybkawsCIYRjkCFQ7SA2NpYyZcqgteaZZ56hQYMG6cb/Fo6rTJkyheJOuSj9v4jcCQwMtLSgOHfuHNWqVQNkUKPC7HaGQJU7czv4/PPPWbp0KTdu3ODOO++85dm3EELcrrRB29/fn6CgILvmR+QvCeZ2MGHCBLkTL6IKw125EKL4kQpwQgghhIOTYC6EEEI4OAnmQgghhIOTYC6EEEI4OAnmQgghhIOTYJ4DqSNi5ZXAwEDOnDmTp9ssaD/++CNvvfUWAPPnz2fJkiX5kk7GoUlv91zk9bkUQgh7kmBuR4UtmOdmDO2ZM2cyduxYAEaOHGkZtjWvFZZxxoUQojBy2Hbm/oH+t8x7yOchxrYZy/XE6/Rb3u+W5QG+AQT4BhB1PYrB3w5OtywoIChH6c+aNYtvv/2WhIQE7rvvPqZNmwYYo5GdOnWK+Ph4nnvuOUaPHk1ycjJPPPEEwcHBKKUYOXIktWrVIjg4mGHDhlGyZEm2bdtmdWSwyZMn8/PPP+Ps7EyvXr2YPXs2x48f55FHHiEpKYk+ffowZ84cYmNjCQoKYvbs2axZswYwxtD28/MjICCAt956i9WrVxMXF0fHjh359NNPUUrh7+9Px44d2bJlC4MGDcLf35+JEydauqgNDAykevXqzJs3j4ULF+Ls7EzTpk1ZsWIFhw4dws3NzdKVbalSpfD29mbHjh23dFMbEBBAyZIlOXDgAOHh4SxZsoSlS5eybds22rVrZ+mpat26dbz55pskJCRQv359lixZwuLFiy3jjHt4eFi6WX3ttddYs2YNJUuW5KeffqJq1aqEh4czcuRIIiMj8fT0ZMmSJdSuXfuWYyaEEEWJ3Jnnwrp16zh8+DA7duwgNDSUkJAQNm7cCMDixYsJCQkhODiYefPmcfHiRUJDQ4mIiOC///5j7969jBgxgsGDB+Pn58fy5csJDQ21GsgvXbrEDz/8QFhYGHv27OH1118H4LnnnuPpp5/m33//tXTRmJ1x48bx77//8t9//xEXF2cJ+ADR0dH8/fffjB8/nmeffZZVq1YREhLCyJEjee211wB477332LVrF3v27GHhwoWAMfxnq1at0qXj5+fHpk2brObh8uXL/PXXX8yZM4eBAwcyYcIEwsLC2Lt3L6GhoURFRfH222/z559/snPnTvz8/Pjggw8YP348NWrUYMOGDZZAfu3aNdq3b8/u3bvp2rUrn3/+uWU/hw8fzp49exg2bBjjx4/P9TETQghH4bB35lndSZdyKZXlco9SHjm+E08rs/HEu3btyrx58/jhhx8AOHXqFIcPH6ZRo0YcO3aMZ599lv79+9OrVy+b0insY2ifPXsWT0/PdOlUqVIl04FZBg4ciFKK5s2bU7VqVZo3bw6Aj48PJ06c4PTp0+zbt49OnToBcOPGDctgMhm5urpajkfr1q35448/ANi2bRvff/+95dhMmjQp18dMCCEchcMGc3vKbDzxoKAg/vzzT7Zt20apUqXw9/cnPj6eihUrsnv3bn7//XcWLFjAt99+y+LFi7NNp7CPoV2yZEliYmLSrRcfH2+1lAFujnFeokSJdOOdlyhRgqSkJJycnLj77rv55ptvsj02Li4ulmOQ1fjqaY+TrWO9CyGEo5Fi9lzIbDzxmJgYKlasSKlSpThw4AD//PMPAFFRUaSkpPDAAw8wffr0bMfOTlXYx9Bu0qQJR44cSbeNQ4cO0axZsxwdz1Tt27dny5Ytlm1ev36dQ4cOAdkfq1QdO3ZMd2w6d+4MZH7MhBCiKJA781zo1auX1fHE+/Tpw8KFC2nRogWNGjWiffv2gBHsR4wYYblrnjFjBmBUChszZkymFeAK+xjaXbt25YUXXkBrbbnr3bJlC2+++WaujqunpyeBgYEMHTqUhIQEAN5++20aNmyY6TjjGc2bN4+RI0cya9YsSwW4rI6ZEEIUBTKeeRFgzzG0n3vuOQYOHEjPnj3ZtWsXH3zwAV9++aVd8iJukv8XkZYMgeoYbmc8cylmF7fl1VdftbT/joqKYvr06XbOkRBCFD9SzF5I3HfffRw/fjzdvP/7v/+jd+/e2X7WnmNoV61alUGDBgFYasELIYQoWHYJ5kqpCcAoQAN7gRFa6/isP1W0pTZnE0IIIXKqwIvZlVJewHjAT2vdDHAChhR0PoQQQoiiwl7PzJ2BkkopZ6AUUHg6KBdCCCEcTIEXs2utI5RSs4GTQBywTmu9Lu06kZGR+PndrNA3evRounTpYnMagYGBlr6+z507Z+m+MyAggICAgNvcAyGEEKJwsUcxe0XgHqAuUAMorZR6NO06np6eBAcHW16jR4/OURoBAQEsX76cDl26c+l6Ih26dGf58uW3HcidnJzw9fXFx8eHli1b8sEHH6Trcc2aEydO8PXXX99Wutb069eP6OhooqOj+fjjjy3zg4KCLN2cCiGEKB7sUczeEziutY7UWicC3wMd8zKB5StW0rbHAJbsS8Sl72SW7EukbY8BLF+x8ra2W7JkSUJDQwkLC+OPP/7gl19+sYyWlpn8Cua//PILFSpUuCWY56WMXaTaOkRqboZSFUIIkXv2COYngfZKqVLK6DasB7A/rzYeERHBpKnv4TxgCu4NOuBSsQbuDTrgPGAKk6b+HxEREXmSTpUqVfjss8+YP38+WmtOnDhBly5daNWqFa1atWLr1q2AMYTppk2b8PX1Zc6cOZmul9bMmTMt44JPmDCBu+66C4D169fz6KNGIYa3tzdRUVFMnjyZo0eP4uvry0svvQQYTdUGDx5M48aNGTZsGNY6Bjp69Ch9+vShdevWdOnSxTI4SkBAABMnTqR79+68/PLLt0yHhobSvn17WrRowX333cfly5cBo1OKV199lW7duvHhhx/y3Xff0axZM1q2bEnXrl3z5JgLIYSwzh7PzLcrpVYBO4EkYBfwWV5tf/7CRST79MfF2TXdfOXsSrJPX+YvXMSM6bnrbjSjevXqkZKSwoULF6hSpQp//PEH7u7uHD58mKFDhxIcHMx7772Xbozx69evW10vra5du/L+++8zfvx4goODSUhIIDExkc2bN99Sd+C9997jv//+s/TbHhQUxK5duwgLC6NGjRp06tSJLVu2WPooTzV69GgWLlxIgwYN2L59O2PHjrUM4nLo0CH+/PNPnJycCAgISDfdokULPvroI7p168aUKVOYNm0ac+fOBW4OpQrQvHlzfv/9d7y8vIiOjs6T4y2EEMI6u7Qz11q/CeRNRM0g7NBhnDy6WV3mVLkO+w9ZH2s7t1LvehMTExk3bhyhoaE4OTlZBgjJyJb1WrduTUhICFevXsXNzY1WrVoRHBzMpk2bLHfsWWnbti01a9YEwNfXlxMnTqQL5rGxsWzdupUHH3zQMi+1L3SABx98ECcnp1umY2JiiI6Opls34/g+/vjj6baROpQqGAObBAQE8NBDD3H//fdnm2chhBC5V+R6gPNp2IAd+8JxqVjjlmXJF8Np4nNHnqV17NgxnJycqFKlCtOmTaNq1ars3r2blJQU3N3drX5mzpw52a7n4uKCt7c3S5YsoWPHjrRo0YINGzZw9OhRm/rbTju8qLXhQVNSUqhQoYLlbj6j1CFRM5vOTNr1Fi5cyPbt21m7di2+vr6EhoZSuXJlm7YjhBAiZ4pc3+zjxozCKWwtOulGuvk66QZOYb8y7qlReZJOZGQkY8aMYdy4cSiliImJoXr16pQoUYIvv/yS5ORk4NahOzNbL6OuXbsye/ZsunbtSpcuXVi4cCG+vr63jMlt69CgaZUrV466devy3XffAUbpwu7du7P9XPny5alYsSKbNhmlG19++aXlLj2jo0eP0q5dO9566y08PDw4depUjvIohBDCdkUumHt5eTFz6mSS1kwn/vBWEi9FEH94K0lrpjNz6st4eXnlettxcXGWpmk9e/akV69eluE+x44dy9KlS2nfvj2HDh2y3KW2aNECZ2dnWrZsyZw5czJdL6MuXbpw9uxZOnToQNWqVXF3d7fa1r5y5cp06tSJZs2aWSrA2WL58uV88cUXtGzZEh8fH3766SebPrd06VJeeuklWrRoQWhoKFOmTLG63ksvvUTz5s1p1qwZXbt2pWXLljbnTQghRM4U2SFQIyIimL9wEV8s/ZInHn+McWNG3VYgF8KRyBCoIi0ZAtUxyBCoGQQGBjJs2DC2bdpApVLObNu0gWHDhll6hRNCCCGKkiJXAQ6k21YhhBDFi0PdmRfGRwJCFDbyfyJE8eMwwdzd3Z2LFy/KD5UQWdBac/HixUybRgohiiaHKWavWbMmp0+fJjIy0t5ZEaJQc3d3t3QaJIQoHhwmmLu4uFC3bl17Z0MIIYQodBymmF0IIYQQ1kkwF0IIIRycBHMhhBDCwUkwF0IIIRycBHMhhBDCwUkwF0IIIRycBHMhhBDCwUkwF0IIIRycBHMhhBDCwUkwF0IIIRycw3TnKoQQxUlgYCCBgYEAnDt3jmrVqgEyxLOwToK5EEIUQmmDtr+/P0FBQXbNT45oDUrZOxfFihSzCyGEyDtJZ+FkV4j9zd45KVbkzlwIIQqxwMBAQkND8ff3L/zF7XH/QMT9kBwDOs7euSlWJJgLIUQhFhAQQGBgIEFBQYW7uD36Czg/lgtJ1Tjk9jGdy95n7xwVK5kGc6XUHhs+H6m17pGH+RFCCOFI9A04PwGiPyYsoQO9158gRb/Csecext3Z3d65KzayujN3AvplsVwBP+dtdoQQQjiMpAsQMRjiNkGlF6lV7hVa/Pco7/Z4VwJ5AcsqmD+ltQ7P6sNKqbF5nB8hhBCOID4ETt/LtYQoZpy8l9fqv0U5l5L8MuwXe+esWMq0NrvWenN2H7ZlHSGEEI4hMDAQf39//P39ady4seXv1PbuFjHLILwTB68k0y7Ii3d3/MT64+uNZVevwiOPwI4dBZ7/4izbCnBKqQbADKApYCk30VrXy8d8CSGEKGDZtm3XSXDhJbg8l/9d8GHElnDcnBNZ99g6etbrCQcOwP33w6FDcNdd0LZtge9DcWVLO/MlwCdAEtAdWAZ8mZ+ZEkIIcfsiIiJ45Y1p7Dt2ilfemEZERETuN5YUBad6w+W5zDvdjcEbwmji2ZSdo3cagfzHH43gHRUFf/wBo0bl2X6I7NkSzEtqrdcDSmsdrrWeCtyVv9kSQghxO5avWEnbHgNYsi8Rl76TWbIvkbY9BrB8xcqcbyw+FE74QdwWqB5Iv1aLmNh+IhsDNlKrTA147TW47z5o3BhCQqB79zzfH5E1W9qZxyulSgCHlVLjgAigSv5mSwghRG5FREQwaep7OA+YgouzKwAuFWug67Zm0tTp+HfpjJeXV6afnb9wkeVuftwwT7zUi2y6WJpvLw5iXsPh3KEU7/d+Hy5dMp6P//67cSf+0UfgLrXY7cGWO/PngVLAeKA18BgwPB/zJIQQ4jbMX7iIZJ/+KDOQp1LOriT79GX+wkVWP3fr3XwCbQe+y6NLKtL9z8v8dmIXUdejjJVDQ8HPDzZsgM8+g88/J3DFCtsq0Ik8Z8udeYrWOhaIBUYAKKUG5muuhBBC5FrYocM4eXSzusypch32H9p0y/zM7+bbsHL5s9z96t2sGLGC8u7l4auv4MknoXJl2LgR2rUDHHxwGAdny53550qp5qkTSqkhwOv5lyUhhBC3w6dhA5KjrHcTknwxnCaN7rhlflZ385XaDMfp91Lc32sg/6tZEx57jFA3N+6rXZvA/fvzZR9EztgSzAcDS5VSTZRSTwLPAL3yN1tCCCFya9yYUTiFrUUn3Ug3XyfdwCnsV8Y9dWtN87AD+3DyqGN1e84e3jjjxLL4eIIvXKRNpWqsfOY55n/3XeEb7KWYyjaYa62PAUOA/2EE9l5a65j8zpgQQojc8fLyYubUySStmU784a0kXoog/vBWktZMZ+bUl2+t/HZ9Iz7VNpAceczq9pKjwknZ8Cdtj55lyYCJnBn8FksOJOe+drzIc5kGc6XUXqXUHnPAlVVAJcAb2G7jICxCCCHsZNiQh9mxfjUjmrqS+NtMRjR1Zcf61Qwb8vDNlXQSRE4lJdwfd78kLgV/afVung2L2eFaBucR83Bv2BGXijVwb9AB5wFTmDT1/26v/brIE1lVgBtQYLkQQgiR57y8vJgx/U22bdrAjOlvpl+YGA5nhhFxcQsjdlbnj4iz+N57B+dWT0M3649T5TokRx7H6e8vaVu2JDv8hmZZO/6W7YsClWkwz26QFSGEEA7qyio49ySQzOlS09ge9T4L+y9kdOvRnDlzhvmvTSH0m0/wTUxg3HPjePr4aZwqelvdVGa140XBymo8851a61ZZfdiWdYQQorAKDAy0tIE+d+4c1apVA9I3sSpSUq7B+QlER37Omsh6PNrlD9q51uPk888ZTc4SE/FasIAZy5Zwyt2dWpv+hjZt8HljGjv2heNSscYtm0y+GE4Tn1trx4uClVUxe5Nsno0roHwe50cIIQpMsWoXHb8LzjxC0OkDDP+3HGevn6RTK0VdV4xAfviw0ZtbcDCMGsXoAwf4tU0bwKgdv6zHAHTd1umK2i214+etttdeCVNWwbyxDZ9PzquMCCGEyAcpCTwx6DgJR9vwepg77+9X3FGpKpsf/pK6FeuC1rB4MYwfb3TF+r//wf33E+fvb9lEau34SVOnk+zT13iefjEcp7BfrdeOFwUuq/HMw214nS7IzAohRHGUkJCQu9HP4v6BE3fyaN+T+G+sxOz913iq9VPsemoX7Wq2g4sXYfBgo1/1Dh1gzx5jCFMrbKodL+zGlu5c85xSqgKwCGgGaGCk1nqbPfIihBCF2fIVK9l9/BzHLf2lh7OsxwBmTp2ceSBNuQaRr3PlwlzKutfk5Y+a88yYSUwpWZm+Dfoa66xfD8OHQ2QkzJoFEycSuGxZujoE/ubdeeqjiNRllUo5s23TBrZt2lB06xc4GKW1LvhElVoKbNJaL1JKuQKltNbRqcv9/Px0cHBwgedLCFF8FcZn5hEREbTtMQDnAVNueVadtGY6O9avvrWI+9pf6LOj+O7YcZ7bXYr3erzPkokrbu5bQgK88QbMng2NGsHXX8Odd+Z53gvj8SzslFIhWmu/3Hw22x7glFL/Z8s8WymlygFdgS8AtNY30gZyIYQQhhyNfpYcA2dHE36gBwM2nOfhbeBVvgktqre7uc7+/UZx+qxZMGaMMfZ4PgRyUfBs6Zv9bivz+t5GmvWASGCJUmqXUmqRUqp02hUiIyPx8/OzvD777LPbSE4IIRyTMfqZ9f7SjfbdR4yJ2DVw3IfFuxbR9Ddn/o5UzOk9h39G/cOd1e9EaQ3z50Pr1nDqFPz8M3z8MZQqVYB7I/JTVu3MnwbGAvUyNFErC2y5zTRbAc9qrbcrpT4EJgNvpK7g6emJFLMLIYo7n4YNsm7f3aQGnHkUriwHt2Z41JzAXTFBLOi3gNrlaxsrhoXxUWioMVRpnz5GzfXq1Qt2R0S+y6oC3NfAr8AMjGCb6qrW+tJtpHkaOK213m5Or8qwfSGEEMAD9/RnweeDiT30Dy4Va1C6eQ+cy3oY7bv/+4GAcWeZ8HcMHhW681qf3xikXBnYfCJKKePZ+IwZ8O671NIavvwShg0Dpey9WyIfZNWdawwQAwxVSjkBVc31yyilymitT+YmQa31OaXUKaVUI631QaAHsC832xJCiKJq+YqVTJr6HiU7DsPZ05sbF44T+eN7lKzlQ+lLOxny8CnuDinB6WvwXLuWoIzn6kop2LbNaG62bx888gjDjx/np0cftfMeifxkSwW4ccB54A9grflac5vpPgssN4vvfYF3b3N7QgiRKxEREblrw53PeZo09T2cB0yhpDlKWelGnag29F1STm2n4cMn+MA5ifKlG7Bl5Bbm9JljfPDqVXj2WejUyfh77VpYvpwYV9esExQOz5YKcM8DjbTWPlrr5uarxe0kqrUO1Vr7aa1baK3v1Vpfvp3tCSFEbixfsZK2PQawxNKGO7FQjNGdVS32Uu0eZds/ZXnnrncIGb2TDrU6GAvXrgUfH1iwwAjoYWHQr1+B570wXhwVB7YE81MYxe1CCFEkBAYG0qFDB558fjLOA6bg3qBDoRqjO8ta7B516VGhH692eRVXJ1e4cAGGDoUBA6BsWdiyBT780Pi7gBXWi6PiwJYe4I4BQUqptUBC6kyt9Qf5lishhMhHAQEBHDwazvF9iYVyjG6fBrXYsf94prXYW/g0NvpU//JLmDABYmNh2jSYPBnsVKSe9tGAi3lMXSrWQNdtzaSp0/Hv0ln6cM9HttyZn8R4Xu6K0Swt9SWEEA7L5jbcBUmnQPRi7uu6kNiQZeikG+kXp45S1rcX9O4Njz8OjRvDrl0wZYrdAjnksIMbkeeyvTPXWk8DUEqV1lpfy/8sCSFE/su2DXdBj9EdtwPOj4P4f6lVyw/Xrse48sNk3HwHW0Ypu7rlaz7r3Aavu+4CZ2fj+fiYMVDClvuy/GVcHHWzusy4ONpUwDkqXmypzd5BKbUP2G9Ot1RKfZzvORNCiHw0bswonMLWZn73+9SogslI0nkSIwL46I923LN+D7raMqo33sG5JWf5b+Nf5ihl/8cIdYStl88w7H/fwb33Gs3Oxo4tFIEcjIuj5Khwq8uSL4bTpFEBXxwVM7Z8C+YCvYGLAFrr3Rh9qwshhMNKHaM78quJxB/eSuKlCOIPbyVpzfSCGaM75Ro68m3WbKpL85VLGb8Lrjl3IMZ9ICiFq5MrXl5ezLinP+tjopjxv29xd3WFTZvgm2+gVq38zV8OFZqLo2LKpiFQtdanVPpeg5LzJztCCFEwAgMDCQwMxPn6RZK2BBIdf4OqHpV4/pkx+TtGt74B0Z9xMnwaj2yJYksUNKzkzc9D5jGg4QAsv7XnzsGrr0JgIF7OzvDFFzy9dCl/de6cf3m7DakXR5OmTifZp6/l0YBT2K8Fc3FUzNnUNE0p1RHQSilXpdSLmEXuQgjhqAICAggKCqJVq1Y0a9yA9r4+hB85yIQJE/InQZ0MMV9xaX8DOP8sVco2IdGlKQv6LWDv2IMMbDTwZjesM2dCw4bw1Vfw4os82rYtjBxJSiHvinXYkIfZsX61+WhgJiOaurJj/er8vTgSgG135mOADwEvjH7V1wHP5GemhBCiyNAaYtew48DzTNl5jANXXTn05Grcy/Xnn3rcvBPXGlavhokT4ehRGDgQ3n8fGjTgur+/XXchJ7y8vJgx/U22bdpgt6Z9xZEttdmjgGEFkBchhCharm9k14FneTNkD6vPQGX3srzc+XVSSvcEpbDcZ4eFGe3F//gDmjSB334zmp4JYaNsg7lSqi5GX+readfXWg/Kv2wJIYQDiw+FyFf5J/xXOqyHCm6leLv7y4xvN4Gybmm66bhwAaZPh08+MXps+/BDePppcHGxW9aFY7KlmP1H4AtgNZCSr7kRQghHduMIBw8/R9jZX7i/dkXaNf4/PnJx5tGWI6ngXuHmetHRMHs2zJ0LcXHw1FPw1lvg4WGnjAtHZ0swj9daz8v3nAghhKO6cYhDR1/lnR3f81W4pmrJMgzsehAXV0/GpY3P167BvHlGBbfoaHjoIaMb1saN7ZVzUUTYEsw/VEq9iVHxLW3f7DvzLVdCCJHPIiIimL9wEfuOnaKkqzM1PCvlfCPxu9l/ZDKvb/+NH06Dm5MTE9o+yaQu03Bx9by5XkICfPopvPOOUbQ+YIBRvO7rm2f7I4o3W4J5c+Ax4C5uFrNrc1oIIRzO8hUrmTT1PZJ9+uPSdzJx54+x+58VLF+x0qZmVCnXNnPt/NuUvfE7166WYkOkO692GsOz7SdTtUzVmysmJUFgoFGEfuoUdO8OP/4IHTrk276J4smWYH4fUE9rfSPbNYUQhV5qZykA586do1q1aoDR7jogIMB+GSsgmY3uVfKOtlmP7qU1N67+ytf/TmT2noO08XBjSd/p+DUYR4SfGyVdSt5cNyUFVq40Bj85cgTatYMlS6BHjwLaS1Hc2BLMdwMVgAv5mxUhREFIG7T9/f0JCgqya34KWuroXi62Dn2qU7hycQWf/fMyc/edJiIOmleuTu+W74JHAAAlnVLXNduKv/467N0LzZvDTz8ZbcYLeYcvwrHZEsyrAgeUUv+S/pm5NE0TQjgcm0f30klwZQVcnMG07fv44BDcVbMxX3SdRa87+pOui2utjbbhU6fCjh3QoIHRf/pDDxWagVBE0WZLMJcufIQQRUZWQ58mXjjOf/v/4evvB7H++gaG146lW61mTOj6EUO7+uHn1T79B5KS4Lvv4L33YM8eqF0bFi0yxhl3Nn5eb/exRmBgIKGhofj7+3Pu3Dn8zd7gistjEWEbW3qA+7sgMiKEEAVh3JhRLOsxAF23NSpNUbtOukHK3pVUGRXOsL1HKOlUgjbeY+jmvYCaqgQ1024kLs6o2DZrFhw/bvTaFhgIQ4eCa/ri+9t9rBEQEEBgYGCxexwicibTYK6U2qy17qyUuopRe92yCNBa63L5njshhN0U1YpyVkf3ijzGlZ1fEdfuPGdc3ZnRfjhPtJmGZ5kq6T8cHW301jZ3rtHErF07mDPHeCYuxenCjjIN5lrrzuZ72czWEUIUXUW5otwjDw7A1XM97378MbUiNHfUSKLeh52o0/Bt+jV+DKcSTuk/cPasEcA/+QSuXoU+fWDyZOjaVSq2iUIhqzvzLHtQ0FpfyvvsCCFE/om5upcvd0zk4z1/sf9KCpX8nPC/Wov/dtTig3t/v/UDR44YRemBgcbz8Ycegpdfls5eRKGT1TPzEIzidQXUBi6bf1cATgJ18ztzQghx27SGuI3sP/Y2bX76k2tJ0MazIoF9n6aT5wh63N2fuBsRvPLGNMaNGWW0MQ8ONoL4qlXGoCcjR8KLL0L9+vbeGyGsyvQhj9a6rta6HvA7MFBr7aG1rgwMAL4vqAwKIfJeRIQRvPYdO8Urb0wjIiLC3lnKc1evn2Tx5mF8uq4WnPSnkUsIT/u0ZUfAz+wYewnnYy3o1vdB4po/gEvfySzZl0jbznezvFETaNPGaGo2aRKcOGEUr0sgF4WYLTU22mitf0md0Fr/ClhvpCmEKPSWr1hJ2x4DWLIv8WYQ6zGA5StW2jtrty0lJYk/983iseV1qfpBHZ5Y/zXLjkajq35OiQYRzLp3O23qDEzXC1ypRp1wqVgD9wYdcH7gHSZdvk7EtGlw8iTMmAFmxb/bURwunoR92RLMo5RSryulvJVSdZRSrwEX8ztjQoi8lzaIuTfocDOIDZjCpKn/V+BBJjXIDXr40dsLcomnIeptJn5Xmbu/m8Tq8BMMb+jDtkeXsvnpq6iKo6DEze5W5y/8nOSm/dI1TQOzF7gujzI/UUP58rezaxZF+eJJFB62BPOhgCfwg/nyNOeJYiAwMBB/f3/8/f1p3Lix5e/UJkvCsaR2ZWo1iJldmRaUtEEupEK3HAe5mLhIFm0dS+dPKrHz39oQ9QYBDRuxou84zr5wkYUP/Uf7+sPT99R27Rp8/jlhixbj5OltdbtGL3BH8mAPC9/Fkyi6bOk05hLwnFKqjNY6tgDyJAqRotw8qTiyuSvTfBYREcH4yVMpM/jddIOd6LqtGT/51UwHO0lOSeatJU+y/sIKQpLiiE+BhqXh281VOFbrBQYPnYSvtQSPHoWPP4bFiyE6Gp+addgRecJqL3DJF8Np4nNHnuxnjvuBd3AZ+yaQ3uoKTrbBXCnVEVgElAFqK6VaAk9prcfmd+aEEHkrq65M8zKIZWf+wkW4tH7AagmBS+v70wU5rTXnY/ZRTW8i/mIgH5zbjrOCkY3qcXUrLH3jIKqElZ+yGzdgzRr44gv49VdwcoIHHoBx4xjn7c2yngPR9fxu6QXOKexXxs1bnSf7WVgungqKBG37saWYfQ7QG/M5udZ6N9A1PzMlhMgf48aMwilsLTop/YjGliD21KgCyYcR5OpYXeZUuQ77Dh0m5PRGXl4zgHoflKbnkmZw/mlKl7jCX/eN5+zEEyx46CgnQ2vdGsj37IEJE8DLywjeu3bBG28YtdJXrIDOnfGqWZOZUyeTtGY61w9sJvFSBPGHt5K0Zjozp75sfQjUXPBp2IDkqHCry5IvhtOkUcFcPImiz5aBVtBan1LpezlKzp/sCOHYCnsXqFa7Mr0YjlPYr3kaxLLj07AB28OsF3MnRh5nU/xq/L5YjrOCu6u781CjvqTUeYcS7r60sdbj2uXL8PXXxpjhISFG2/B77oERI6BXL8ugJ2kNG/Iw/l0607HbXcTdSOaJxx9j3LzVeXoMsuoHPi9LAISwJZifMovatVLKFRgP7M/fbAnhmByhjkFqEJu/cBFfLJ2ZL0EsO5YgZ6WYO2X3lzQffZ3h7f25t/mLVK7UF5SVQsTkZPwuXYIhQ+DHHyEhAVq2hA8/hEceAQ+PbPPh5eVF3ZrVAfLl2XVhuXgSRZ8txexjgGcAL+A04GtOCyEclJeXFzOmv0nTejWZMf3NAg0qWmsuOV2k3QMeXPhuPNcObiLxUgRxhzaR9PNLfDjlcf5++TpP+G+gcuX+twbyo0eNYnNvb2bv3Qvr1sGTT8LOnRAaCuPHZxvI07bSOHfunGWI0fxopTFsyMPsWL+aEU1dSfxtJiOaurJj/WqGDXk4R3mUliQiK7bUZo8ChhVAXoQQRdSN5BtciD1DTecTXIr6Bt+ln6FdocMYKLdvBVyog2+zLoz76G/rFxaXLsH338NXX8HffxuDm/TuzZvlyjFt505wc8tRfjI+9sjvUpTUi6dtmzbYXAJQWB7NCMdgS232mcDbQBzwG9ASeF5r/VU+500I4cAuxV3i10M/8fO+L/jt+A7aVII//ROprNz4/i4/2tUdRrVqAeBUwfoGrlyBn34yKq2tW2cMdNKgAbzzDgwfDjVr8re/f44DuRBFkS3PzHtprScppe7DKGZ/ENgASDAXQtwq+TLPrRnCgt1/kKw1Vd3hodou3HdHJ6jxDJTpwz2Nylj/7PXrsHatEcDXrjWeg9eubdROHzIE7rxThhwVwgpbgrmL+d4P+EZrfUnJP5MQRZa1GvmhoaEEBgbeUuybmJzI9ojtrN7/Nb8d+ZmNvRtQPmkzrV2TeLlpaQY16E2b+k9SovRdoFxvTQyMgP3770YA//lno5e2atXgqaeMAN6uHZSwpXqPEMWXLcF8tVLqAEYx+1illCcQn7/ZEsVRYW/WVVxYq5Hv7++f7hzsObeLN/96nr/Ct3PlRgLOCvyrQGSsC+VrvMDwOveCe1vrtdABEhPhr79g5UrjWXhMDFSuDMOGGQG8a1ejkxchhE1sqQA3WSn1f8AVrXWyUuo6cE/+Z00UN47QrKs4iomPIcojiqdXP87AWh708zyD05lf2BVxhSE1oVftBvS8YwjlPYaCa+PMi8GvXjXuwH/6yShCv3wZypWD++4zAniPHkb7cCFEjtnaaczlNH9fA67lW46EEHaXlJLEOxunc7BjMJVnViS5uSZ8TxhNNFC6Kk1r3Mvx0X1QZXqBU+XMN3T2rFF0/tNPsH690cVqpUowcKARxPv0AXf3AtsvIYoqm4K5EKLoyPg4o5t/N66Vvkbj3o3xa92A531q43ztV5bv/B+1q6cwqpqiflIpht01GZfyA8CtJSqz4nOtYf9+I3j/9BNs327Mr1cPnnnG6JWtUyerPbIVVTL4iCgImf5HKaU6aa23KKXctNYJBZmp4kieF4uCkvqdenrB0xw8vJeT5Y8QVyKZkMshHAyB5ysDzjXY+9Bw3v1gB9PGbsG/x70EDH3D+gaTk2Hr1psB/Ig5fKifH0yfDvfeCz4+xbYWur3+h+UionjJ6vJ4HtAa2Aa0KpjsFF/yvFjkl6SUJHae3cnfJ/5m26mNrOwzApeETbh7rqJ0/CX6eUKHCoq7G3SjTpV+ULo3uDXHTSn+3ulvvR342bNG2+/ffoM//oCLF43n3XfdBRMnwqBBxkAnwm4kaBcvWQXzRKXUEsBLKTUv40Kt9fj8y5YQ4nZtPrmZdza+xeaTm4lNjAOgYVmIOLoG7zLuzGzbgTl39YTS3enR/2VG/bUh0225pKQYtc9//90I4Hv2GAuqVoV+/aB/f+jb16jQJoQocFkF8wFAT+AuICSvE1ZKOQHBQITWekBeb1+I4kBrTXhMONtObWPb6W1sO7WJdzreR68qcSSdWc3JyDAeqw3dqjjRtZYf1T36QOnu4N4elxI3e05LTrHyDPzIEe6LiICBA/l5y5abtc07dYIZM4zKay1aSBtwIQqBTIO52Sf7CqXUfnMM87z2HMboa3IpL4SN4hLjuJZ4DY9SHpyMDqfdojacuxYJQCknRdvKGhUVCs4u+FdrTdijr0Cpu6BkRyhRKuuNR0cb/Z7//rvxOnaM5wC05rdq1bjvk0+ge3coWzaf9zJ7ERERzF+4iH3HTvHKG9MYN2aUjEAmijVbqpReVEr9AHQCNLAZeE5rfTq3iSqlagL9gXeAibndjhCFUV4GmoNRBwk+E8z2iO1sO7WV0POhjPZpy4K2VfC6toW+VaLwqwgdPMvSvHpnnMt0gZKdwL0NlCiZ9cZjYmDTJggK4tOQEKPTlpQUKF3a8ux72LJlLN++nQ/9/blv0KBc7UNeW75iJZOmvkeyT39c+k5myb5wlvUYwMypk7MdiUyIosqWYL4E+BqjT3aAR815d99GunOBSYD9L/GFyEO5DTQpOoXDFw8TcjaEhKQERvgGQGI4vZZ15uTVKEo5l6BNRc2LjTR9PbZBwh04le3P4kGdjODt2jjz3tZSxcTA5s0QFGS8du40grerK3ElS8KUKeDvDx06gKvR9WrEd9/l1aHJExEREUya+h7OA6bgYo6D7lKxBrpuayZNnY5/l85yhy6KJVuCeRWt9ZI004FKqedzm6BSagBwQWsdopTyt7ZOZGQkfn5+lunRo0czevTo3CYpRIGwNdBorUkd3+Cj7R+xav8qdp3dydUbsQA0LFeaEWUmQXIUi1uDp7szTT1b41ymsxG4S3YE56rZZ+jKlfTBOyTEErxp394YE9zfH9q14/m+fQl607ahOe1p/sJFxoWSc/p+3pWzK8k+fZm/cJHNQ4xK0y1RlNgSzCOVUo8C35jTQ4GLt5FmJ2CQUqof4A6UU0p9pbV+NHUFT09PgoODbyMJIQpeVoEmqWkf7p30MKXuduJg1AFOP/UVzjd2En76CxJiw3ms9g38KkLrSoomnt5Qqi24t6FHnTbg3iLzQUpSaQ0nTxrtvVNfu3cbbcBTg/frrxvBu317KJl1EXzGRwUJCYWjq4mwQ4dx8uhmdZlT5TrsP7TJ5m1J0BZFiS3BfCQwH5iD8cx8qzkvV7TWrwCvAJh35i+mDeRCOCKtNcFhu3Hy7Gt1ubOHN7uDF9K6YxwPVI/n+olelHOB2b71oP19xjPukm3B/U4okcnwoGklJMCuXUbQ3rbNeD9zxlhWurQx0tirr94M3qWyqfyWhrVHBVeOn2P5ipV2fybt07ABO/aF41Kxxi3Lki+G08TnDjvkSgj7s2WglZNA4aj5InIk9e4q7NBhfBo2cOgav4Wph7yo61G4OblR1q0sW09tZdIfk9hzfg+x0YlUppH1QBN1jBc6uTLjnq7g7mcGb7+s+zVP6/z5m0F761YIDjYCOkDdukbQ7tjReDVvnuvuUjN7VOBWSJ5JjxszimU9BqDrtkalKQHRSTdwCvuVcfNW2y1vQtiTXTtI1loHAUH2zENR9dTYZ/jy+18o02EoLlW68c/eY3zYujOP3d+PTz9ekKNtFYZmQPbqIe9S3CWW71nOvsh97I/az77IfURej2RZv1d57A4vSkWvp0TCHobXuU6tx5N5/60v0fXa3hpo9v3GuI9CbOsVLSbGqJz2779G0A4OhuPHjWWurtC6NYwbZwTuDh2gevU829+sHhUkNunNjNlzmD9ndp6ll1NeXl7MnDqZSVOnk+zTF6fKdUi+GI5T2K/MnPqyw16sCnG7is9oB8VIREQEa/7aiuejH1iCikvFGpS8oy1r1kwnIiLC5h89ezYDyuxu/Ny5cznaTlYlFFprjkcf50DUAfZF7rME7UeaPcKzfo8Rf3UT438bTwVXN5pWcOee6jdoUhba6nfhPPi6VWRj/zbg1hLcfampzjPp3bdI9umXIdC8Yv2Yx8YaxeWpQTs4GA4durnc29vo4/zpp43OWlq1ytdRxrJ6Ju3iWZclX79Dh3Zt7FrcPmzIw/h36cz8hYv4YulMnnj8McbNWy2BXBRrSmtt7zzcws/PTxfnCnC3e+f5yhvTWLIvEfcGHW5ZFn94KyOautpU4zciIoK2PQbgPGDKLXeaSWums2N9/v2AZjwGaadzcnwsFyNN++Hk6U1S5AniQ7/jgVH+LH71M5JSkij5TkmSUpIAqFqyFE0ruBHgDcNrX0ZrOBcP1UqWQrn7gJsPuJrvbs3AueYtA4ikXjx8sfRLI9CkXjxcuwZ796YP3Pv3GzXMAWrWNAJ36qt1a/DwuM0jmbWbeV3GE48P52rsVVaddLf63bl+cCs3Ik9QOmp/vp77nJBxDERRopQK0Vr7Zb/mrbK9M1dKVQXeBWporfsqpZoCHbTWX+QmQZH/8qrGb26bAdnj+XaKTuFc7DmOXT7G4YuHKaFK0NOzp9Xnv+71/Fj1+fNM79Eer0onWN4uiRoloWk5qOSuwbWuGax9UG7NqO7qAy61s2/HbfKqUYMZT46AH75jhpsTPP+8UbP8yBGj1jlAlSrQpg0MHnwzcOdhcbktrJW6EPoriSkaNyvPpK8E/4jHoEkknfPOURMwIUT+s6WYPRCjk5jXzOlDwEpAgnkhlVc1fnN7UZBfz7evJlzlePRxojyiWBq6lMd9HwfgsR8eY9W+VcQnxVvWbVqpBgdOLCG5aV+rFyOl7hzK/K//x4zJ/XiodSOj0xW3ZuDiDcrJ9kzFx0NYmBGsU1979sDly8wAoy13vXrQsiUMG2a8+/kZz87tOCRoVm3ib6ycRNTyFynV/mFcPb1JjAznSvCPlG09COeyHujEhBw1ARNC5D9bgrmH1vpbpdQrAFrrJKVUcj7nS9yGvKrxezsXBbmpNHftxjVOxpzkZMxJzlY/a+lc5Y2/3mBrp62Ue8/sxr85PPPLGIbXjkMlHaNtqb1Ua1iWeqU0dUslcEcZ8C59hvtfisOp5jCraTl51GP/ufZQfbFNx4LERDh8GPbtM15hYfDff3DwoNGWG4wmYc2bw0MPQcuWPPPZZyz4++9COZJYVqUurm0fou6ZjRw8sp34YyE4l6+Kx6BJOJc1ivylCZgQhY8twfyaUqoyRhtzlFLtgZh8zZW4LXlV4ze3FwWZVZp7dfIY2tx9JydjThIeHc5Tfk9RxrUM87bP462/3+JiXJq+iBpDTNTXVHA6T32nrfSvqmhcpgIHNiRw6ZwbrerHE7HnaWpWd+XZBnXBpQ241geX+uDaEFwb4dNqac4vRhISbgbtsLCbwfvQIUgynqujlNEczMcH7r/fuNtu2RLq1083gljYypWFMpBD9qUuHknVOHf5iNX6EtIETIjCx5ZgPhH4GaivlNoCeAKD8zVX4rYlxsdRu6IbZ3avIvrqNSqULU0Nz0okxsfZvA1bLwqSUpI4F3uO01dOs+fwHl5+7QPK3TfjluLb5994lqS9EZZx8npVK0HzCnCHCuLBOmWpU9KF2u5XqFPyOrVLQ7moR0FBQPXyRP3uypxVlUhuNgSnet4cjDzBt2PWMPPNyQwbOtRq/seNGZ35xch/axn3+FuwbJkRqPfvN4L24cM377SVMgJ006Zwzz3Gu48PNGqUo05YCqPsSl1atWzG8KGDefL5iZTr9Ig0AROikLOpNrtSyhloBCjgoNY6MT8zJbXZC1cN3QPHDjB1zgx+/W09rTo3x7t7NaLdonmp40t0rNWR1QdXM2iF0a+Q+tudyhWep3TDzrdsJ+7QJrqWXcDrY2KpXRoquJiPjUuUMZ5VO9cx3l3q8OY7gUx7Zwm41ifibBwN/Lqka2oHNtSqj4tj+fwFTPp4Ccm+9+Dk6U3yheM4bf6KmdHnGJZs3mk7ORnPtX18bgbspk2NoJ1Nt6dZKWznMS1bWyp06NAB/559bq2ZX0gU5mMsRE7ld232+zPMaqiUigH2aq0v5CZRYV/XblzjwrULlHIpRdUyVbkcd5k5/8zh7NWznLt2jrNXz3I29izT/KcxqtUoEkolsNJjGTwKQURQ/kxpapapQPT5T8HlS1rpwyxsX5OabjHM2uLEYc96VtN19qhH4oUGtGz0mCVo4+INJSpaKoNFREQwf8EivlgVz43SvzJuzCjmf/oFZToMTRd0IE2t+rdnMKN3Tzh2DI4eNe6uDx2CkycZpjX+wPyIY+wvW54mFcszrk83vO68Exo2NF5161pGCSsubC11cXNzY8b0N9m2aYPUXheiELOlmP0JoAOwwZz2B/7BCOpvaa2/zKe8iRw4H3ueC9cuWF7nr52nsUdj+tzRh7jEOLov7c75a8Y61xOvA/B6l9eZftd0kpNjeHvj21QpVYHqpctRrWRJWnhVpG7y13ByOQ3jT7O5Z0mqu8VR1R1KO18DrgHL4GplvNxq8lRTX3CpyeYWxzlw/ETmz6lb3AOVJljdh8yetVcrXwYX3+FWP+NUuQ77l02ChWavduXLGwG6c2dLsPZq2JAZDRpAWRlxNy3pfEWIosOWYJ4CNNFanwdLu/NPgHbARkCCeR5K0SkkOidy+OJhNJqGlRsC8PG/HxMeHU7U9Sii4qKIuh5Fx5odmdVrFgB3fHQHseYQmqkCmg+mj1dF3BPPUcn5Gg2qlKeKe3mquCVTxTUBv/JfwqEPqZx8lRsPgnOJy8Bl48OqJDjHga5BydKt6FS+Pwu/+Jkxz0wzOkpxqWm8l0hfDD1uQoT5nNrPtopTyclw9iwRISFMemUazg+8c8uz9sOfPUWp6sesXyBEnqBJ/74w6QXj+XbFirdz+IsdLy8vufMWogiwJZh7pwZy0wWgodb6klIqX5+dO7KklCSu3bhGeffyAGw7tY1jl49xKe4SF+MucinuEpVKVmKq/1QABn0ziK2ntnI5/jIpXVJoOL8h/t7+bHjcKBCZt/1DjkefwLNkBTxKlsHDvSQV9X6IfAOSI5nfvhGlSsRSxSWWKi6Xqep6nQquqyB8FQr4pb2ZsRIVjLGwnaob7841UM41cHauDs41zFd1KFH+lnbQK9aFMubVx7Lcb6vFt1HhOO3+iZm9uuL12WcQHn7zdfo0JCYy38Wd5H7PW20q5dbxYa5sWY7bHVb6PN//G+PWr7atz3MhhCiibAnmm5RSa4DvzOnB5rzSQHR+Zayw2Re5j0MXDxEdH010fDSX4y6TmJLIuz3eBeD1v17n54M/G8viLxN7I5ba5WsT/nw4AG9tfIvfjvxm2V4F9/K0rtoYrvlDyiXaVFLUcmtARRc4fzSMri3rU7fUBThaH5Ij2eV/FXcnUOoCxvUUQBhc/AWcKvN4HU8jCDtVBedqZsCumuG9CpRwy7uDcuWKEYwzvIadPo2/imP+n5+wPymFJokJjEuMx2thmHGBUKMG1KljDM1Zpw7Urk3Ytz/gVNX6s3bX2s25sfMHktbI4BpCCGGNLcH8GeB+oDNGbfalWutV5rLu+ZWxwmb+jvl8EvxJunkepSrzTudRKB1NOXWB+mXdqVC5FhVdvangClXdFUQ8BMmX+KjZGVKaVKOScwwVnONwLhEDbIdTxiF8o465UeVGpIfGsxrg5AXOnuDkQUknT3DytEyTOu1UMWc9ltkiNhbOnr3l9cqBA3D33RARYQTuq1dv/WyVKlCzpvGc+q67+OzXXxk9dSrUqmUE7po1rVY28zlzIcumUhXKlWbr+tXyfFcIIaywZTxzDfzPfKGU6qyUWqC1fia/M1doXP+bF71382TVxlRwiaWi01XKlojBqcRFOF4fgEnVgXRdazsbgTbhPDhV5o5K9Y2xqy2vStanS5TiwfxobpOYCJGRxrjYFy4Y7+fPWw3aVoO0qystlTICfdOm0KuXEZi9vIz3mjWNO243t5t9s8fEcM7Vla8XG72sZdU3e3Yd1NTwrCTPd4UQIhM2DYGqlPIFhgIPA8eB7/MxT4WPTqFeGQ3l6hkBukRF493yd4Vb56vS+dv3ttbGuNdRUUaQthao0/596ZL17ZQqZQzwUaMG+PpC377GdMZXpUoM6d7dpouM3Ayokl1Tqc8XfpL9RoQQopjKNJgrpRoCQzCC+EWMwVWU1rrYFK1blO4Opbfm3/a1NobHvHQJLl6k9eXLsGKFEaBTg3VU1K1/p3YvmlH58kZxd9Wqxl20v7/xd9WqN+en/l22rF0H/Egrq6ZSjhLMM44Y5+/vD+TviHFCCJHVnfkBYBMwUGt9BEApZb2BsDCkpBiVwqKj4fJl490M0Onerc27ccOymfcBUrsoVQoqVTLGtfb0hDvuMCqOeXrenOfhYbxSA7S7ux12Pm84elF6UQja1i5IQkNDCQwMdPh9E6KoyiqYP4BxZ75BKfUbsAKjAlzxc/Qo/P77zQCd2XtMjBHQM+PuDpUrG8G5UiWju9BKlW7OM9/HT5vGvG++MQJ0pUpGd6NCFBBrFyT+/v4SyIUoxDIN5lrrH4AfzCZo9wITgKpKqU+AH7TW6womi/aTeofSLTKSafv2AZBQogTJ5cpRqkYNqFDBeJ7cpInRWUnFisa8tO+pgbtyZZv7+d7z4YfGNoUQt5BHGULcypba7NeA5cBypVQl4EFgMlDkg7nlx+H6de7v2ZPv//oLNwcuwnZkaX+05Qe8eJNzLsStbKrNnkprfQn41HwVH6VKccnV1aGfRTu6atWqyehYQgiRiRL2zoAQWYmIiOCVN6ax79gpXnljGhEREfbOUpERGBiIv78//v7+ltIOf39/SxG2EMJx5OjOXIiClNkoajOnTmbYkIftnT2HJ8XVQhQdEsxFoZSQkMCkqe/hPGDKLaOoTZo6Hf8unaUrVyGEMEkxuyiUzkReItmnf7quXcEYRS3Zpy/zFy6yU86EEKLwkWAuCqW4G0k4edSxusypch32HzpSwDkSQojCS4K5KJRKujqTHBVudVnyxXCaNLqjgHMkhBCFlwRzUSjV8KyEU9hadNKNdPNTR1Eb99QoO+VMCCEKHwnm2ZCmUfbh5ubGzKmTSVoznfjDW0m8FEH84a0krZnOzKkvS+U3IYRIQ4J5FpavWEnbHgNYsi/RbBqVSNseA1i+YqW9s1YsDBvyMDvWr2ZEU1cSf5vJiKau7Fi/WpqlCSFEBtI0LRMRERHSNKoQcPRR1IQQoiDInXkm5i9cJE2jkF7ChBDCEcideSbCDh3GyaOb1WVG06hNBZwj+5BewoQQovCTO/NM+DRsIE2jhBBCOAQJ5pkYN2aUNI0SQgjhECSYZ8LLy6vAm0ZJMzghhBC5IcE8CwXZNEqawQkhhMgtqQCXjYJoGiXN4IQQQtwOuTMvBKQZnCiMpFmiEI5D7swLAWkGJwojaZYohOOQYF4I+DRswI594bhUrHHLsuSL4TTxKR7N4AIDAy13fal3goAEFCGEyIYE80Jg3JhRLOsxAF23dbqidkszuHmr7Zi7gpPVnaAU7QohRObkmXkhYI9mcEIIIYoOCeaFhIwQJoQQIrekmL0QkRHChBBC5EaB35krpWoppTYopfYrpcKUUs8VdB5sJU1zhBBCOAJ73JknAS9orXcqpcoCIUqpP7TW++yQlyxJ0xwhhBCOoMDvzLXWZ7XWO82/rwL7AanhJYQQQuSSXZ+ZK6W8gTuB7WnnR0ZG4ufnZ5kePXo0o0ePLtjMCSGEEA7CbsFcKVUG+B/wvNb6Stplnp6eBAcH2ydjQgghhIOxS9M0pZQLRiBfrrX+3h55EEIIIYoKe9RmV8AXwH6t9QcFnb5wHNKaQAghbGOPYvZOwGPAXqVUqDnvVa31L3bIiyjEpDWBEELYpsCDudZ6M6AKOl0hhBCiqJLuXIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8G8kAgMDMTf3x9/f3/OnTtn+TswMNDeWRNCCFHIKa21vfNwCz8/Px0cHGzvbAghhBAFRikVorX2y81n5c5cCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHZ5dgrpTqo5Q6qJQ6opSanHF5ZGSkPbJVaHz22Wf2zoJdyf7L/hdXxXnfQfYf8MjtBws8mCulnIAFQF+gKTBUKdU07TpRUVEFna1Cpbh/oWX/Zf+Lq+K87yD7D3jm9oP2uDNvCxzRWh/TWt8AVgD32CEfQgghRJHgbIc0vYBTaaZPA+3SrnD9+vUEpVRymlmRQHG6XfdQShWn/c1I9l/2v7juf3Hed5D9b5TbD9ojmCsr89L1Kau1di+gvAghhBAOzx7F7KeBWmmmawJn7JAPIYQQokiwRzD/F2iglKqrlHIFhgA/2yEfQgghRJFQ4MFca50EjAN+B04ANYCfrDVRU4Z5ZhO2PUqpVgWb2/yVXRM9pZS/UipGKRVqvqbYI5/5QSm1WCl1QSn1XybLi/q5z27/i/K5r6WU2qCU2q+UClNKPWdlnSJ7/m3c/6J8/t2VUjuUUrvN/Z9mZZ2ifP5t2f+cn3+ttV1egBNwFKgHuAK7gaYZ1ukH/IrxnL09sN1e+bXT/vsDa+yd13za/65AK+C/TJYX2XNv4/4X5XNfHWhl/l0WOFTM/vdt2f+ifP4VUMb82wXYDrQvRufflv3P8fm3Zw9wtjRRuwdYpg3/ABWUUtULOqP5pFg30dNabwQuZbFKUT73tux/kaW1Pqu13mn+fRXYj9HKJa0ie/5t3P8iyzynseaki/nSGVYryufflv3PMXsGc2tN1DJ+oW1Zx1HZum8dzOKYX5VSPgWTtUKhKJ97WxX5c6+U8gbuxLg7SatYnP8s9h+K8PlXSjkppUKBC8AfWutidf5t2H/I4fm3ZzDPtomajes4Klv2bSdQR2vdEvgI+DG/M1WIFOVzb4sif+6VUmWA/wHPa62vZFxs5SNF6vxns/9F+vxrrZO11r4YrZnaKqWaZVilSJ9/G/Y/x+ffnsHcliZqRbkZW7b7prW+kloco7X+BXBRSuW6714HU5TPfbaK+rlXSrlgBLLlWuvvraxSpM9/dvtf1M9/Kq11NBAE9MmwqEif/1SZ7X9uzr89g7ktTdR+BoabNRvbAzFa67MFndF8ku3+K6WqKaWU+XdbjPN1scBzah9F+dxnqyife3O/vgD2a60/yGS1Inv+bdn/In7+PZVSFcy/SwI9gQMZVivK5z/b/c/N+bdHD3CA0URNKZXaRM0JWKy1DlNKjTGXLwR+wajVeAS4DoywV37zmo37Pxh4WimVBMQBQ7RZ1dHRKaW+waix6aGUOg28iVERpMife7Bp/4vsuQc6AY8Be83nhgCvArWhWJx/W/a/KJ//6sBSZQy6VQL4Vmu9prj89mPb/uf4/Kui8/0QQgghiid7FrMLIYQQIg9IMBdCCCEcnARzIYQQwsFJMBdCCCEcnARzIYQQwsFJMBcih5RSyeZIRmFmd4sTlVK5/l9SSr2a5m9vlclIahk+M0YpNTy3aeaEOYLTGvPvQcrKCH/5kF7HTJZNVUpFKKXeyuE2lyulLimlBudNLoUoXOzWzlwIBxZndsWIUqoK8DVQHqOteG68Crybkw+YbVELnNb6Z27t3Cmv+QOxwNZMls/RWs/OyQa11sOUUoG3mS8hCi25MxfiNmitLwCjgXFmb1VOSqlZSql/lTEO81NgudvcqJT6QSm1Tym1UClVQin1HlDSvNNfbm7WSSn1uXnnv87sJSod8w71RfPvIKXU/yljjORDSqkuVtavbqYfqpT6L3UdpVQvpdQ2pdROpdR3yugvHKVUH6XUAaXUZuD+NNsJUErNN/8OVEp9ooyxuY8ppbopY5z2/WkDZxZpnFBKTTPn71VKNVbGwCNjgAlmXm/ZFyvHYal5nE4ope5XSs00t/ebMrpNFaLIk2AuxG3SWh/D+F+qAjyB0fVkG6AN8KRSqq65alvgBaA5UB+4X2s9GfNOX2s9zFyvAbBAa+0DRAMP2JANZ611W+B5rJcQPAL8bpYotARCldHX8+tAT611KyAYmKiUcgc+BwYCXYBqWaRbEbgLmACsBuYAPkBzpZRvZmmk+XyUOf8T4EWt9QlgIcbdt6/WepMN+14f6I8xbOZXwAatdXOMnrP62/B5IRyeFLMLkTdSR3nqBbRI82y2PEZwvgHsMAN/aneunYFVVrZ1XGsdav4dAnjbkH7qYB2Zrf8vsNi8U/1Rax2qlOoGNAW2KKMbaFdgG9DYzMNhM69fYZQ+WLNaa62VUnuB81rrveZnwsx81MwkDWv5vp/c+VVrnWjmwQn4zZy/F9uOnRAOT4K5ELdJKVUPSMYYm1gBz2qtf8+wjj+3DuGYWV/KCWn+TgZuKWbP4jPJWPm/1lpvVEp1xbhT/VIpNQu4jDGW8tAMefXNIm+ZpZuSId8pZj6SraVha75zkgetdYpSKjFNH9apeRCiyJNidiFug1LKE6NYeL4ZRH7HGCDBxVzeUClV2ly9rTJGySsBPAxsNucn5vezXaVUHeCC1vpzjBG7WgH/AJ2UUneY65RSSjXEGMGprlKqvvnxzAKxLTJLIytXgbK3kaYQxY4EcyFyLrXCWhjwJ7AOmGYuWwTsA3Yqo4nZp9y8O9wGvAf8BxwHfjDnfwbsSVMBLj/4Yzwn34XxDP5DrXUkEAB8o5TagxF4G2ut4zGK1deaFeDCc5toZmlk87HVwH22VIATQhhk1DQhCoBZzP6i1nqAnbPi0JRSU4HYnDZNMz8bCKzRWlurpyCEQ5M7cyGEI4kFRqtcdBoDdAPi8yVXQtiZ3JkLIYQQDk7uzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHASzIUQQggHJ8FcCCGEcHD/D4e2qu0zC/wmAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 576x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import matplotlib.pyplot as plt\n",
"\n",
"fig, ax = plt.subplots(figsize=(8, 5))\n",
"\n",
"params = {'elinewidth': 0.75,\n",
" 'capsize': 3,\n",
" 'markersize': 15,\n",
" 'markerfacecolor': 'C0',\n",
" 'markeredgewidth': 0.75,\n",
" }\n",
"\n",
"ax.errorbar(x, y, yerr=e, fmt='k.', label='Data with errors', **params)\n",
"ax.plot(t, y_hat_fit, 'r', label='Nonlinear least squares')\n",
"ax.plot(t, y_hat_fit_weighted, 'gold', label='Weighted nonlinear least squares')\n",
"ax.plot(t, y_hat_lsq, 'g--', label='least_squares() method')\n",
"ax.tick_params(direction=\"in\")\n",
"ax.set_title('Comparison of unweighted and weighted fit', fontweight='bold')\n",
"ax.set_xlim(0, 3.5)\n",
"ax.set_xlabel('Depth in sediment [m]')\n",
"ax.set_ylim(0, 10)\n",
"ax.set_ylabel('Age of sediment [ka]')\n",
"ax.legend()\n",
"\n",
"plt.savefig('Figure.svg')\n",
"plt.savefig('Figure.png', dpi=300)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "5cc82d8e",
"metadata": {},
"source": [
"---\n",
"[Thanks to Martin Trauth](http://mres.uni-potsdam.de/index.php/2017/02/09/create-publishable-graphics-with-matlab/)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.9.7"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment