pyerrors/examples/01_basic_example.ipynb

353 lines
86 KiB
Text

{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Basic pyerrors example"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Import numpy, matplotlib and pyerrors."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"import pyerrors as pe"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"plt.style.use('./base_style.mplstyle')\n",
"import shutil\n",
"usetex = shutil.which('latex') not in ('', None)\n",
"plt.rc('text', usetex=usetex)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We use numpy to generate some fake data"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"test_sample1 = np.random.normal(2.0, 0.5, 1000)\n",
"test_sample2 = np.random.normal(1.0, 0.1, 1000)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From this we can construct `Obs`, which are the basic object of `pyerrors`. For each sample we give to the obs, we also have to specify an ensemble/replica name. In this example we assume that both datasets originate from the same gauge field ensemble labeled 'ensemble1'. **For correct error propagation it is crucial that observations from the same Monte Carlo chain are initialized with the same name.**"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"obs1 = pe.Obs([test_sample1], ['ensemble1'])\n",
"obs2 = pe.Obs([test_sample2], ['ensemble1'])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now combine these two observables into a third one:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"obs3 = np.log(obs1 ** 2 / obs2 ** 4)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"`pyerrors` overloads all basic math operations, the user can work with these `Obs` as if they were real numbers. The proper error propagation is performed in the background via automatic differentiation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we are now interested in the error of `obs3`, we can use the `gamma_method` to compute it and then print the object to the notebook"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"1.367(20)\n"
]
}
],
"source": [
"obs3.gamma_method()\n",
"print(obs3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the method `details` we can take a look at the integrated autocorrelation time estimated by the automatic windowing procedure as well as the detailed content of the `Obs` object."
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Result\t 1.36706932e+00 +/- 2.04253682e-02 +/- 1.02126841e-03 (1.494%)\n",
" t_int\t 5.01998002e-01 +/- 4.47213595e-02 S = 2.00\n",
"1000 samples in 1 ensemble:\n",
" · Ensemble 'ensemble1' : 1000 configurations (from 1 to 1000)\n"
]
}
],
"source": [
"obs3.details()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As expected the random data from numpy exhibits no autocorrelation ($\\tau_\\text{int}\\,\\approx0.5$). It can still be interesting to have a look at the window size dependence of the integrated autocorrelation time"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAGJCAYAAAC5Lib1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAss0lEQVR4nO3dbWxc133n8d/hsx7IGVG2bMWyaw2TIHHVbDKW4gZJsUg83A2QpovsSlYXzr4oUJObAt0mQEEu/cbImxDUiz5ssciSKbAvGqOQxK2xTQt4V7SLtClSxzabTZWkSKxR4iiRS1vUzFAPHJLDsy94Lz0kZ6ghOXfOffh+gAF5n4Z/jcdzf3POuecaa60AAADQXG2uCwAAAIgjQhYAAEAACFkAAAABIGQBAAAEgJAFAAAQAEIWAABAAAhZAAAAAehwXYArxhgj6T2SFlzXAgAAIqVX0i/sPSYbTWzI0lrAuua6CAAAEEnHJP18ux2SHLIWJOlnP/uZ+vr6XNeCJnu7tKgLr13TUyeP6f6+HtflAABiolQq6eGHH5Ya6AlLcsiSJPX19RGyYmhRXeo5cFC9fX3qI2QBABwIzcB3Y0zaGDNkjLm0y+N3dRziqbuzXSceSqm7s911KQCAhApFS5YxJivppKS0pP5dHH9aUq7JZSHCUvs6NfjYA67LAAAkWChasqy1s9baKUn5nR5rjElrF8EM8bZcWdU7t8parqy6LgUAkFChCFl79JSkC/fayRjTbYzp8x9au/wSMXXz9pL+7Ns/1c3bS65LAQAkVKRDljEmJ2mmwd3HJBWrHkzfAAAAAhPpkCUpba1ttItxXFKq6nEssKoAAEDihWLg+24YY4a8cVwNsdaWJZWrjg+kLgAAACmiLVne1Yivua4D4dbeRpAGALgT1ZasfklZb0yWJA1IkjFmRFLeWjvtrDKEwpG+Hv2XJ9/nugwAQIKFLWTVnIrBGJORdNpae06SrLUzqhrw7rVsDfnbAQDRMlda1NxCecv6I73dOsJdGxBRoQhZfoiSdFZrLVQTkl6tapHKSRqWtCVEeRORnvV+n5B0yQthSLAbt8p68ftv6dO//KAOH+x2XQ6Ae3j+lTf1xy/9eMv633vyffrS4PsdVATsnbHWuq7BCW+urGKxWOTehTE0V1rU86+8qaefeIRvwUAE+C1Zb8zd0hfPf1d/dPbDeu+Rg7RkIXRKpZJSqZQkpay1pe32DUVLFgAg2Y709WwIU+89clAnHko5rCjcZmZmNDo6qv7+fl26VPvWvbOzsxodHVU+n9eVK1daXKF07ty7nU83btzQ8PCwpqenNTIy0vJaXCFkAQBCobJq9b1rBUnS964V9MGjfVwlXEcul9PY2JjGx8fr7pPNZjU6Oqrh4eEdPXehUNCFCxd08eLFugHuXoaHhzU8PKxsNru+7syZM7t6rs388OYHx8nJyYaOSafTktb+fa0KepGcwgEAEC8vXr6uT0y8rGdfuCxJevaFy/rExMt68fJ1x5WFlx8attPfv7Nb+87OzurChQsqFAqan5/fZWXShQsXNgQsSfra17626+fzjY6OamRkRCMjI+vhanBwcNtj/FA2NDSkoaEhZbPZHQfP3SJkIZb69nXqMx86qr59na5LAXAPL16+ri98fVbXi4sb1r9VXNQXvj5L0GqhbDaroaEhZTKZPT1PoVBQPr/xhizpdFqnTp3a03POzs6qUCisrxseHtbMzMyWv1VtfHxcQ0ND68u5XE5TUw3PZb4ndBcilno62/X+B7gHOBB2lVWrL3/jB6p1CZaVZCR9+Rs/0OBjD7as6/DcuXPKZDLK5/PKZDI6ffr0+hgoaa1FJp/PK5/P68aNG5qYmFg/dmpqSplMZj1kpNPp9RN8o887Pz+v119/XZOTk5qamlJ/f7/Onz+vsbGxLa1DkjQ9vXYh/vz8fMNdYbVqaaZsNqvBwUFNTk4ql8utr99rN91rr72mfD6//jr4YbA6eFXL5/MqFAo1W/1mZmY21BYEQhZi6XZ5Rf/81oI+8GCvDnTzNgfC6jtX57e0YFWzkq4XF/Wdq/P62MDhwOs5c+aMzp49ux46BgcHlclklMvlNDExoeHhYc3Pz69vHxgY0NmzZ5XNZtfDjn/izufzmpmZ2fXzjo6ObghwzzzzjF5//fUN9c7OziqXy62HiKmpKQ0PD287TqleLbUC3G5dvHhRg4OD6115uVxOo6Ojewo16XRaN2/e3LDOf33rtbzVa+FKp9N1g1kzcfZBLN0ur+hvf/S2Hj60j5AFhNjcQv2AtZv99iKfz2t6eloXL15cX3fmzBlNTk5qcnJS/f39yufzG4KC3xrkB5SLFy/qqaeeUjqdViaT0cmTJ3f9vNWy2WzNwJDNZje00gwNDckYo9HR0ZrB4161NEsmk9GVK1c0MzOjS5cuaWZmRoODg7p48eJ6uGt0XNTjjz++obuv2vj4uCYnJxsan1atv79/T2POGsXZBwDgzJHexubAanS/vZiZmVE6nV5vHZHWrmCrDjebg0s6nV4/WZ8+fVqTk5M6dOiQstmszp49q5GREU1NTe3qeQcGBnb178hkMpqdna0Zshr5NzZTLpdbD4+jo6N65pln1kPWXkPd6Oiozp49WzeAbacVAUsiZAEAHPro8X4dTfXoreJizXFZRtKDqR599PjOrpLbjUKhsN6F59tp99alS5c0OzurmZmZDSFir8+7E9sFiGb8G++lUChoZmZmyziviYkJnTt3ru4YqZ2Ynp7WwMDAPQNWvW5E/3UIGiELAOBMe5vRc599TF/4+qyMtCFo+cPcn/vsYy0Z9J7NZmvOO9VoKJiamlqfIsC/Su/JJ5/UxMTEnp53pwqFQt3xVXv9Nzbq1VdfrTmYPpPJrP+d3XYX+q1w/jp/uolaocn/e/4A/2pBD3qXmMIBMdXV0abM/QfU1cFbHAi7T584qq9+PqsHUxu7BB9M9eirn8/q0yeOtqSOXC6nkydPrg9g9124cKHuMdWDpwuFwpapAfxWo70873b8q+d8ftCr10rTaC31WsPy+XxD0x9MTU1t6JKUtKV1yx8Hdq9HdcCanZ3V7Ozs+hg1vx5/PrB8Pr9hpnlJGhsb21DL9PT0rroYd4N7F3LvQgAIhcqq1flX39SzL1zWVz53QmdPPeJkxvfR0VENDAysn7hPnz6t2dlZjY+Pa3p6WhMTExoZGdG5c+c0Pj6uTCajsbGx9WBSfcIfGhpab7nZ6fOePHly/dY5/j4jIyMaGxtTOp3W7Ozs+rQN0tYpHKqfe2RkZMOVirVq8Wuenp7W+fPnNTs7q5GREZ06dWp9+9TUlEZHR3X16tW6LV/+jPH+2LBqe5nCoVAo6Pjx4zUDqJ9lpqenNTo6uuU2Qv6UFdJaK1v1a7FTO7l3ISGLkBVLlVWr8kpF3R3t3JYDiJDLPy/q1//kW/qr3/0E9y4MKX+QfCvGNIXRTkIWfSmIpRu3ypr8Zl43bpVdlwIAsVLvykVsxcB3AIBzc6VFzS2U9cbcLUla/3mkt1tH+oKfvgGNa9X0B3FAyAIAOPf8K2/qj1/68fryF89/V5L0e0++T18afL+jqrBZPp/XyZMnXZcRGYQsAIBzTz/xiAYfe2DL+iO93Q6qQT10E+4MIQsA4NyRvh66BRE7hCzE0n0Hu/U7nxxQZxvXdgAA3CBkIZba2oy629pdlwEASDC+5iOWbt5e0l/MXtPN20uuSwEAJBQhC7G0XFnVT2/c0XJl1XUpAICEImQBAAAEgJAFAAAQAEIWAABAAAhZiKWDPR365AeO6GAPF9ACANzgDIRY2t/VoQ8/nHZdBgAgwWjJQiwtLlf0w+slLS5XXJcCAEgoQhZiqXR3WS9efkulu8uuSwEAJBQhCwAAIACELAAAgAAQsgAAAAJAyEIsdbS36WiqRx3tvMUBAG6EZgoHY0xa0lOSzlhrBxs8ZsT7dUCSrLXDwVSHqOk/0KXf/OgjrssAACRYKEKWMSYr6aSktKT+Bo+ZsNaOVi1PGmMuNRrQAAAAghSKvhRr7ay1dkpSvpH9vVavrPfTNykpZ4zJNL9CRM1caVF/eOlHmistui4FAJBQoQhZu3RSUnWg8gNauvWlAAAAbBSK7sKdstYWJB3atDrn/azZGmaM6ZbUXbWqt/mVAQAArIlyS9ZmY5KGvQBWb3ux6nGtRXUBAIAEikXIMsZMSDrvjeuqZ1xSqupxrBW1AQCAZIpkd2E1Y8xpSVfuEbBkrS1LKlcdF3RpcKj/QJd+6+OP6mB35N/iAICIinRLljEmJ0l+wDLGpLm6ENLaZKTp/V1MRgoAcCZsZ6Cac2QZYzJVE4/667KSspJmve0ZSUOS5oMvE2FXvLOsFy9fV/HOsutSAAAJFYqQVRWihrU2/9WE1w3oy3nb/P3Tkl6SNCHpStVjYpuB70iQ8kpFP7y+oPJKxXUpAICECsWAFWttXtI571Fr+5SkqarlgrZO4QAAABAaoWjJAgAAiBtCFgAAQAAIWYil/d0d+tXMYe1nCgcAgCOcgRBLB7s79LGBw67LAAAkGC1ZiKXySkU/eec2VxcCAJwhZCGWineW9cI//px5sgAAzhCyAAAAAkDIAgAACAAhCwAAIACELMRSW5tRen+n2tqM61IAAAnFFA6IpfsOduu3Pn7cdRkAgASjJQsAACAAhCzE0tsLZf2Pb17R2wtl16UAABKKkIVYstbq7lJF1lrXpQAAEoqQBQAAEABCFgAAQAAIWQAAAAEgZCGW0vu7dPbUw0rv73JdCgAgoZgnC7HU1dGm96T3uS4DAJBgtGQhlhYWl/XNH72thcVl16UAABKKkIVYurtU0exPb+ruUsV1KQCAhCJkAQAABICQBQAAEABCFgAAQAAIWYilnq52/auHU+rpanddCgAgoZjCAbHU19OpT33gAddlAAASjJYsxNJyZVVzpUUtV1ZdlwIASChCFmLp5u0lPf/Km7p5e8l1KQCAhCJkAQAABICQBQAAEABCFgAAQAAIWYgns3aTaBnXhQAAkspYa13X4IQxpk9SsVgsqq+vz3U5AAAgAkqlklKplCSlrLWl7falJQsAACAAoZmM1BiTlvSUpDPW2sEGjxmRVPAW09bac8FUh6i5causv/6n6/rMrxzV4YPdrssBACRQKEKWMSYr6aSktKT+Bo8ZkSRr7ZS3nDPGTFprh4OqE9FRWbW6cWtJldVkdocDANwLRXehtXbWC0v5HRw2Jmmq6jlmJA01uzYAAIDdCEXI2iljTEZr3YOFGttyra8IAABgo1B0F+5Cps76gta6HLcwxnRLqh6c09vckgAAAN4VyZasbcyr/piuMUnFqse1VhWF1uvb16nf+PB71Lev03UpAICEilvI2m7Q/LikVNXjWEsqghM9ne0auP+gejrbXZcCAEioqIasegPk0/W2WWvL1tqS/5C0EFRxcO92eUXfuTqv2+UV16UAABIqkiHLWpuXVPAGwG/eNuOgJITM7fKK/v6NdwhZAABnwhayanb3GWMy/rxYVcYl5ar2Oa2qKR0AAABcCkXIqgpRw5KyxpgJLzT5ct62dd7s7mljzGlv31NMRAoAAMKCG0Rzg+hYmist6vlX3tTTTzyiI309rssBAMQEN4hG4nV3tOt9DxxUdwdXFwIA3IjqZKTAtlL7O/XrH3qP6zIAAAlGSxZiqbJqtbC4zA2iAQDOELIQSzdulfWnf3dVN26VXZcCAEgoQhYAAEAACFkAAAABIGQBAAAEgJAFAAAQAKZwQCzd39ut3/3Ue9XeZlyXAgBIKEIWYskYo452AhYAwB26CxFLN28v6eJrP9PN20uuSwEAJBQhC7G0XFnVtZt3tVxZdV0KACChCFkAAAABIGQBAAAEgJAFAAAQAEIWYqm3p1ODjz2g3p5O16UAABKKKRwQS/u62nXioZTrMgAACUZLFmLp7lJFl39e1N2liutSAAAJRchCLC0sLuvSD/5FC4vLrksBACQUIQsAACAAhCwAAIAAELIAAAACQMhCLHW2t+nYoX3qbOctDgBwgykcEEuHDnTpzMmHXZcBAEgwvuYjlqy1WqmsylrruhQAQEIRshBLby+U9Scvv6G3F8quSwEAJBQhCwAAIACELAAAgAAQsgAAAAJAyAIAAAgAUzgglg4f7NZv/9px7e/iLQ4AcIMzEGKpvc2ot6fTdRkAgASjuxCxVLyzrL/63i9UvLPsuhQAQEIRshBL5ZWKfvwvt1ReqbguBQCQUKHqLjTGjEgqeItpa+25Bo4ZkpT2jhuQNG6tLWxzCAAAQOBCE7K8gCVr7ZS3nDPGTFprh+9xzJQfqowxaUlfk3Qm8IIBAAC2EabuwjFJU/6CtXZG0tA9jhmsbrXyfk8HUBsAAMCOhCJkGWMyWuseLNTYltvm0IIx5pLXguU/T77O3+g2xvT5D0m9e68cYXWgu0Mff+99OtAdmsZaAEDChCJkScrUWV/Q9i1Tz3jH3jTGTEjKbdO9OCapWPW4tqtKEQkHujv00eP9hCwAgDNhCVn1zEvqr7fRa/makDQtaUTSGb9Vq4ZxSamqx7FmFopwWVyu6Mrbt7S4zNWFAAA3wh6y6gYsSfJar/LW2jNau7KwX9Lrtfa11pattSX/IWmh6dUiNEp3l/WX3/2FSneZJwsA4EZYQlbNcVRa6yqsN8bKH8c1I0nW2ry19nGtjdM6HUiVAAAADQpFyLLW5rUWjraMzfJDVA0ZvTunVrXJJpYGAACwK6EIWZ5xSetXEnqtUVNVyxl/Li1pPXxla4zBetxaOx1wrQAAANsKzaVX1tpzxpiRqq6+U5uuFMxJGpZUPQv8GUljxpgbevdKxNEWlIuQa28zOnywS+1txnUpAICEMtZa1zU44c2VVSwWi+rr63NdDgAAiIBSqaRUKiVJKe9CurrC1F0IAAAQG00PWcaYR6t+/4gx5j9UrwNaYW5hUf/9b97Q3MKi61IAAAkVREvW+uB1a+0/Wmv/V/U6oCWstLSyKiWzNxwAEAJNGfhujElJekprp7RBYzYMNk5LOiXpT5vxtwAAAKKgKSHLWls0xsxo7cq+AUk3qzYXJP3XZvwdAACAqGjaFA7W2quS/rMx5klr7UvV2xiTBQAAkqbp82RZa18yxnxKa92EvrPeA2iJQwe69PQTj+jQgS7XpQAAEqrpIcsYc0FrAatQtfojzf47wHY629t0pK/HdRkAgAQLYsb3S9bar1WvMMY8GcDfAeoqLS7rtZ/M6+Sj/err6XRdDgAggYKYwuFGg+uAwCwuVfT/flbU4lLFdSkAgIQKoiVrwBjzfyTNVq3LaW0aBwAAgERoasjy5ss6K+n85k3N/DsAAABh19SQ5c2X9Yy19h+r13tzaAEAACRG08dkbQ5Ynps11gGB2dfVruwvHdK+rnbXpQAAEqpZt9X595JmrLUlY8zv19hlUNK/bcbfAhrR29Opf/3++12XAQBIsF23ZBljfrtq8VlJJ73ff1NrY7CqH4d3+3eA3VhaWdUvCnfXbhINAIADxlq7uwONqVhrt/TFGGM+UmNM1pZ1rhlj+iQVi8Wi+vr6XJeDJpsrLer5V97U0088wqSkAICmKZVKSqVSkpSy1pa223cvY7JqXjFYK0yFLWABAAAEbS8ha3dNYAAAAAmwp5YsY8zvG2M+7HW9AQAAwLOXqwutpGlJj0t61hhzXNK81mZ6f1Xe1YZ7LxHYOWOM9nW1yxjmwQUAuLGXge+rktLVQcoY8xGtzfie09qAsPc1pcoAMPAdAADs1E4Gvu+lJWtGa61Yf+Ov8Aa4M8gdAAAk3l7GZJ2R9AVjzKNNqgVomndulfU///6q3rlVdl0KACChdh2yrLVFa+1TkgaaWA/QFKurVoU7y1pd5SJYAIAbe753obX2pWYUAgAAECdNv0E0AAAACFkAAACBIGQhllL7O/W5jzyk1P5O16UAABJqL1M4AKHV3dGuR+874LoMAECCEbIQS7fKK/qna0X9yrGUDnbzNgeAuJsrLWpuYeu0PUd6u3Wkr8dBRYQsxNSd8or+IX9DA/cfIGQBQAI8/8qb+uOXfrxl/e89+T59afD9DioiZAEAAIWzJWgnnn7iEQ0+9oDemLulL57/rv7o7If13iMHdaS321lNoQpZxpgRSQVvMW2tPdfgcROSrniL89ba6QDKAwAgtsLYErQTR/p6NoTB9x45qBMPpRxWFKKQ5QUsWWunvOWcMWbSWju8zTFpSS9JetJaWzDGZCW9Lsm0oGQAANbREoTNQhOyJI1JOu4vWGtnjDGXJNUNWZImJJ231ha8Y2aNMYOBVolI6O5o1weP9qq7o911KQASgpYgbBaKkGWMyWite7BQY1vOWjtT59AhSQPe8Rlr7cw2+yJBUvs79ekTR12XASBBaAnCZqEIWZIyddYXJKVrbfCClSRlJeUl5Y0xk5Iu1gpaxphuSdXv9N7dFovwW6ms6lZ5RQe7O9TRzpy7AIJHSxA2C0vIqmdeUn+dbX7IKlhrZyXJGDMq6aqkQzX2H5P0XNMrRCjN317S86+8qaefeCQSYyEARH9ME7BZ2ENWvYBV7TX/F2/we7pOF+O4pD+oWu6VdK0JNQIAmiDqY5qAzcISsvJ11qe32VZvfUE1uh+ttWVJ61+RjOECRAAIE8Y0IW5CEbKstXljTMEYk7HW5jdtqzmQ3Tsmr7VANVu1Ka2q1i0AQDQwpglxE4qQ5RmXlJPkz5N12v/dW85IOr1pgtJRSWflhSzvmBl/jBYAJAljmoBwCU3IstaeM8aMeEFJkk5tmog0p7U5s85VHTNtjOn3JzKVdNhayzxZ0JG+HsZwIHEY0wSES2hClrQWtKoWpzdtm1JVy9am9QCwZ1FvCWJMExAuoQpZcC/qJxnf/O0l/d/vv6V/88sPqv9Al+tyEiPq75+otwQxpgkIF0IWNoj6Sca3UlnV9eKiViqrrktJlKi/f2gJAtBMhCxswEnGrai3BEX9/UNLEIBmImRhgyN9PTp8sFvfu1aQJN1ZWtEHj/apvY15xVoh6i1BhBQAeBchCxu8ePm6vvyNH+h6cVGS9OwLl/UnL7+h5z77WCRuuOy3BM3fXtLcwqL++a0FzS2UaQkCALQcIQvrXrx8XV/4+qzspvVvFRf1ha/P6qufz4Y+aG1uCfrz7/xMEi1BAJAElVW73hPzvWsF5z0xhCxIWntjfvkbP9gSsCTJSjKSvvyNH2jwsQdD3XVISxAAl8J2kt+pKNcfxp6YNid/FaHznavz62/MWqyk68VFfefqfOuK2oUjfT068VBKxw7tkyQdO7RPJx5KRaKrEEC0vXj5uj4x8bKefeGypLWT/CcmXtaLl687rqwxUa7f74nZfB7ze2Jc/RsIWZAkzS3UD1i72c+1O0uVDT8BRMPmlpTKaq329fAJ60m+UVGu/149MdJaT4yL9xIhC5KkI72NtfQ0uh8A7FRUW1LCfJJvRNTrD3NPDCELkqSPHu/X0VSP6vW8G0lHUz366PH+VpYFICGi3JIS5pN8I6Jef5h7YghZAaisWn37yg397+/+XN++ciO06b9ae5vRc599TJK2BC1/+bnPPhaJAZCVVasfXi9Jkn54vRSJ179aVLtLfNTvVhTrj3pLSphP8o2Iev1h7okx1obzTRs0Y0yfpGKxWFRfX1/Tnnfz1Q3SWgtQVOaZon63qN8t6nfj21du6D9+7R/uud+fP/Or+tjA4RZUtDPU71Zl1eoTEy/rreJizaBuJD2Y6tG3Rj/VlIaCUqmkVColSSlrbWm7fWnJaqIoN3f7Pn3iqL41+il95XMnJElf+dwJfWv0U6H+gPZF/fWnfreo352ot6REfbhF1OsPc08MIatJot7cXa29zehDx9KSpA8dS0emizDKrz/1u0X9boW5u6cRYT7JNyLq9UtrDQRf/XxWD6Y2vkceTPU4nUibkNUkUR846JsrLeryz4t6Y+6WJOmNuVu6/POi5krh/Abpi/rrT/1uUb9bUW9JkcJ7km9U1OuXwtkTw4zvTRL15m7f5tvSfPH8dyWF/7Y0UX/9qd8t6nfLb0n5wtdnZaQNLXJRaUmR1k7yg489qPOvvqlnX7isr3zuhM6eeiT0dfuiXr8Uvp4YQlaTRL252+fflmazsN+WJuqvP/W7Rf3u+S0pmwfuPxiBgfvVwnaS36mo1x82hKwm8Zu773V1Q5ibu6WtNyiOiqi//tTvFvWHQxxaUoBqjMlqkjgMHIyyqL/+1O8W9YcHLSmIE0JWE8Vh4GCURf31p363qB9As9Fd2GR+c/d3rs5rbmFRR3rXmuj5NtYaUe9uoH63qB9AMxGyAtDeZkI5K25StLcZHb/vgCTp+H0HIneCiXp3CfW7FfX6gTihuxCxtFyxG34CANBqhCwAAIAA0F2IWJkrLWpuoayr79yWJF1957b6D3TpSG93JKemAABEFyELsbJ5xvrn/vL7ksI/Yz0AIH4IWYgVf8b6u0sV/eTGbT16+ID2dbWHfsZ6n98SV33vSEm0xAFABBGyECvVM9afCvns1rVE9d6RAICtCFmIpcXlit6cv6NH+verp7PddTkNi+q9I320xAHAuwhZiKXS3WX99feu6+knHolUyIrqvSN9tMQBwLsIWQCahpY4JBnvH2xGyALQNLTEuRX1k3zU6+f9g80IWQDgiXpLXNRP8lGvn/ePW2EMicba8Nx2xBgzIqngLaatted2ePwla+1gg/v2SSoWi0X19fXtrFCE3o1bZb34/bf06V9+UIcPRuMDDtgr/ySzWVRaIqJef9RF/fX/w0s/2hASfc0OiaVSSalUSpJS1trSdvuGJmR5AUt+sDLG5CSdsdYON3j8aUkXrbUN3Q2VkAUAQHy0KiTuJGSFqbtwTNJxf8FaO2OMuSTpniHLGJOWFL1JkQAAQFOEcUxoKG4QbYzJaK17sFBjW66Bp3hK0oV7/I1uY0yf/5DUu6tiEQlzpUX9t5d+rLnSoutSAAAJFYqQJSlTZ31BUnq7A70QNtPA3xiTVKx6XGu8PERRZTUcXeEAgGQKS8iqZ1737gZMW2vzDTzXuKRU1ePYHmsDAACoK0xjsmrZNmAZY4astVONPJG1tixpfUScMQ2NjwcAANiVsLRk1WuJStfbZozJSnotqIIAAAD2IhQtWdbavDGmYIzJbO76s9bWG2/VLylbNTB+QFqfCiJvrZ0OrmKE3aEDXfpPH/slpfZ1ui4FAJBQYZsnq+B3/3nzXg3682R5VyCerjdBqdey9TrzZAEAgKDsZJ6ssHQX+pOQpo0xp72AdWrTRKQ51Zkzy9t/zPt9osFpHxBjxbvLuvSDf1Hx7rLrUgAACRWalqxWoyUr3uZKi3r+lTf19BOPhG5yOgBAdEWyJQsAACBOCFkAAAABIGQBAAAEgJCFWNrX1a5Tj/ZrX1e761IAAAkVinmygGbr7enUJ953n+syAAAJRksWYqm8UtHP5u+ovFJxXQoAIKEIWYil4p1lTb9+TcU7zJMFAHCDkAUAABAAQhYAAEAACFkAAAABIGQhlowx6u3pkDEN3S8cAICmYwoHxNL9vd367V/LuC4DAJBgtGQBAAAEgJCFWHp7oaw//bu83l4ouy4FAJBQhCzEkrVWC4srsta6LgUAkFCELAAAgAAQsgAAAAJAyAIAAAgAIQuxlNrfqdOPH1Nqf6frUgAACcU8WYil7o52Pdy/33UZAIAEoyULsbSwuKxv/fgdLSwuuy4FAJBQhCzE0t2lil79ybzuLlVclwIASChCFgAAQAAIWQAAAAEgZAEAAASAkIVY6u5s14mHUurubHddCgAgoZjCAbGU2tepwccecF0GACDBaMlCLC1XVvXOrbKWK6uuSwEAJBQhC7F08/aS/uzbP9XN20uuSwEAJBQhCwAAIACELAAAgAAQsgAAAAJAyEJstbcZ1yUAABLMWGtd17DOGDMiqeAtpq215xo8RpIGJMlaO9zg3+qTVCwWi+rr69tFtQAAIGlKpZJSqZQkpay1pe32Dc08WX5YstZOecs5Y8zkdqHJGDNhrR2tWp40xlyy1g4GXzEAAEB9YeouHJM05S9Ya2ckDdXb2RiTlpT1fvomJeWMMZmAakRE3LhV1vOv/FQ3bpVdlwIASKhQhCwvFKWttYUa23LbHHpSUnWgyns/000rDpFUWbWaK5VVWQ1PdzgAIFnC0l1Yr+WpoDqByQtkhzat9gNZftN6GWO6JXVXrerdSYEAAAA7EYqWrG3MS+rfwf5jkoZrtYh524pVj2t7rg4AAKCOsIeshgOWMWZC0nl/4HwN45JSVY9jey8PAACgtrB0F27p3vOkt9m2zhhzWtKVbQKWrLVlSeWqY3ZYIqKkb1+nPvOho+rb1+m6FABAQoWiJctam5dUqHVVoHeVYV3+wPiqqR/SXF2Ins52vf+BXvV0trsuBQCQUKEIWZ5xvTtw3W+dmqpazlRNPOqvy0rKSpr1tme0Nu3DfGtKRljdLq/o9Z/e1O3yiutSAAAJFZqQ5c3unjbGnPYC1qlNE5HmJK0ve/NjvSRpQtKVqsdEnYHvSJDb5RX97Y/eJmQBAJwJy5gsSetByze9aduUNk5WWtDWKRwAAABCITQtWQAAAHFCyAIAAAgAIQux1NXRpsz9B9TVwVscAOBGqMZkAc2S3t+lf/fhh1yXAQBIML7mI5Yqq1Z3lla4QTQAwBlCFmLpxq2yJr+Z141b5XvvDABAAAhZAAAAASBkAQAABICQBQAAEABCFgAAQACYwgGxdN/Bbv3OJwfU2cb3CACAG4QsxFJbm1F3W7vrMgAACcbXfMTSzdtL+ovZa7p5e8l1KQCAhCJkIZaWK6v66Y07Wq6sui4FAJBQhCwAAIAAELIAAAACQMgCAAAIACELsXSwp0Of/MARHezhAloAgBucgRBL+7s69OGH067LAAAkGC1ZiKXF5Yp+eL2kxeWK61IAAAlFyEIsle4u68XLb6l0d9l1KQCAhCJkAQAABICQBQAAEABCFgAAQAAIWYiljvY2HU31qKOdtzgAwA2mcEAs9R/o0m9+9BHXZQAAEoyv+QAAAAEgZCGW5kqL+sNLP9JcadF1KQCAhCJkAQAABICQBQAAEABCFgAAQAAIWQAAAAFgCgfEUv+BLv3Wxx/VwW7e4gAANzgDIZY62tuU3t/lugwAQIKFKmQZY0YkFbzFtLX2XBDHIP6Kd5b17fw7+ljmPqX2d7ouBwCQQKEZk+WFJVlrp6y1U5JmjTGTzT4GyVBeqeiH1xdUXqm4LgUAkFChCVmSxiRN+QvW2hlJQwEcAwAAELhQhCxjTEZrXX2FGttyzToGAACgVcIyJitTZ31BUroZxxhjuiV1V63qlaRSqdRAeYiahdKi/u77b+o3PphWj5ZclwMAiImd5IawhKx65iX1N+mYMUnPbV758MMP76IsRMULX3JdAQAgpnolbZu4wh6ydhqwtjtmXNIf1Nh3fhd/oxG9kq5JOiZpIaC/gfp4/d3i9XeL198tXn+3WvH690r6xb12CkvIytdZn95m246OsdaWJZU3rQ6sr9AY4/+6YK2lT7LFeP3d4vV3i9ffLV5/t1r0+jf0vKEY+G6tzUsqeIPZN2+badYxAAAArRKKkOUZl7R+VaAx5rSqpmcwxmT8ebEaPQYAAMCV0IQsb6b2tDHmtBeWTllrh6t2yUka3uExLpUlfVlbuyjRGrz+bvH6u8Xr7xavv1uhef2NtdZ1DQAAALETmpYsAACAOCFkAQAABICQBQAAEABCFgAAQADCMhlprHhTTRS8xbR3FSRapGqqjwFJCtEVp4ljjLlkrR10XUfSGGMmJF3xFuettdMu60kSY8yQ1ibFLmjtM2jcWltwWFJsGWPSkp6SdKbW50wYzsWErCbzT/DW2ilvOWeMmeRE3xrGmAlr7WjV8iQneje8aVVy99wRTeOddF6S9KS1tmCMyUp6XZLZ9kA0hff5P+WHKu+/x9cknXFYVix57+2TWgu0W26nF5ZzMVM4NJkx5qak49XfXIwx1lrLh1zAvA+0i1r7VlPw1vknmQHvLgFogapvmJO891vHGDMp6Ur1N3ZjTI67YLRGrS90fMkLlvdlbsxa+/im9aE4FzMmq4m8W/ykazUNG2P4Rt8aJyVV32rJD1bp1peSaE9JuuC6iAQakjTt3SEjJ3GbsRYrGGMueV8y/HMCX+5aLEznYkJWc225j6KnIE7ygbPWFqy1h6y1s1Wr/f+h+KBrEe9DjBN7i1XdxzWrtc+bvNddzhe81nlGa+eBm964uBxDRZwIzbmYkNUa86rRZ4yWGJM0zMDTlkrTNeuEf2IpWGtnvf8Go1rrQkcLeJ8zE5KmJY1IOuO3aiEUWn4uJmS1BgHLAe+b5Hl/4COCZ4wZ4ko2517zf/FO+mlas1rD+8zJW2vPaO3Kwn6tjQlFOLT8XEzIaq56397T22xDALzBkBsGACNY3kUGr91zRwSl3mdMQfW7T9AkVeOAZiTJWpv3BmMXvM8jtE5ozsVM4dBE1tq8MaZgjMls7i5h8GnrVA349S/dTUvqpwsrcP2SslWtJgPS+qXUeVq4guV9/uS1FqiqxyWmRfhthYzenZOp2mSL60i8MJ2LaclqvnFVzQ3kfYOhu6pFvNaUrKRZ7wqrjNauuJp3W1n8WWtnrLXn/Ie8k4u3TMBqjVFJZ/0F7/NnZtPFIAiAd/LO1hiD9Tjv/0DV6wIMxbmYebIC4H9z9xZPVU+OieB4H25XVePqEeZqai3vA+2spNOSzkm6RGtua1TNOC5Jh/n8aR3vM2hM0g29eyXb+uSkaB7vC7T/OZPV2ufMq9WBNgznYkIWAABAAOguBAAACAAhCwAAIACELAAAgAAQsgAAAAJAyAIAAAgAIQsAACAAhCwAAIAAELIAAAACQMgCkDjGmBFjzBVjjDXGXKy636KMMUPGmNe9bZe8GdSrj53wtl3ZvA0AqjHjO4BEMsZMShqqdcsl73YcE5IO1bolijHmorX2TPBVAogyWrIAJFWhgX223HzWa/XifoAA7omQBSCprkjrN/Vd5y2f8hY3bPNkrLX5GusBYANCFoCkmvd+bm6tekrSeK1txpgha+1U0IUBiAdCFoCkKng/0/4KY0xW0mt1tmUk0YIFoGGELABJ5bdkZarWnbTWzqp2K9dpa+1MSyoDEAuELABJVfB+9kuSMea0pAuSVHVF4YC3LSdpurXlAYg6QhaARKoavJ72B79vmq6hoHe7CxnsDmDHCFkAku6wpKestZtbquYl9TPYHcBuEbIAJFlBUk5rg91rbcuKwe4AdqnDdQEA4NC8pNe8we61ts0z2B3AbtGSBSDJZlV/9vZZScMtrAVAzHDvQgAAgADQkgUAABAAQhYAAEAACFkAAAABIGQBAAAEgJAFAAAQAEIWAABAAAhZAAAAASBkAQAABICQBQAAEABCFgAAQAAIWQAAAAEgZAEAAATg/wMgaj+3lh+p5AAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 640x395.55 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"obs3.plot_tauint()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This figure shows the windowsize dependence of the integrated autocorrelation time. The vertical line signalizes the window chosen by the automatic windowing procedure with $S=2.0$.\n",
"Choosing a larger windowsize would not significantly alter $\\tau_\\text{int}$, so everything seems to be correct here."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Correlated data"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now generate fake data with a given covariance matrix and integrated autocorrelation times:"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"cov = np.array([[0.5, -0.2], [-0.2, 0.3]]) # Covariance matrix\n",
"tau = [4, 8] # Autocorrelation times\n",
"c_obs1, c_obs2 = pe.misc.gen_correlated_data([2.8, 2.1], cov, 'ens1', tau)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and once again combine the two `Obs` to a new one with arbitrary math operations"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Result\t 3.27194697e-01 +/- 1.53249111e+00 +/- 2.49471479e-01 (468.373%)\n",
" t_int\t 4.75187177e+00 +/- 1.33949719e+00 S = 2.00\n",
"1000 samples in 1 ensemble:\n",
" · Ensemble 'ens1' : 1000 configurations (from 1 to 1000)\n"
]
}
],
"source": [
"c_obs3 = np.sin(c_obs1 / c_obs2 - 1)\n",
"c_obs3.gm() # gm is a short alias for gamma_method\n",
"c_obs3.details()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This time we see a significant autocorrelation so it is worth to have a closer look at the normalized autocorrelation function (rho) and the integrated autocorrelation time"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAmYAAAGfCAYAAAD1WR7GAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCXklEQVR4nO3dfXAc52Hn+d8DvgAg8TIEKZAwKUoERVmEKVkCrbeYsXwW4Cgbs3K+o8ykpFzqyisy8l3O0lWWtLR1x2ivSgqYeNc++1Yhld3UlqVyKDHHxPLGzImS7axtWZIJMRJMxSYJyhQpiLAEDgC+gC/Ac390NzQYzMvTwDSme/D9VE21Zrp75sE0NfOb59VYawUAAIDyqyp3AQAAAOAhmAEAAMQEwQwAACAmCGYAAAAxQTADAACICYIZAABATBDMAAAAYoJgBgAAEBMEMwAAgJiYW+4CAEgeY0xKUpekDkmtkrol/SzjkCZJe6y1e7POa/fP+4Sk7dba3TNS4DIwxnRI2mqtvbcEz9UqabukdMbDu+S9/westb3TfQ3/ddolbZaUkn9drbXbHc9Nybu2B/3zF7ueC+BDhiWZAEyVMWaTpOckdVprD2Tte0FSk7V2fY7zzqgCg1lGOJG80NRrre2c5nO2ynuP77bWpjMe3yXpC5LWlyKYBaEsM0wVuoY5zj/mlyXt39/kP9+0gykwm9CUCWA60gX23Sup3Q8Q2QaiKU55WWvT1tqt1tqt8moRS6FL0hOZocx/ra0lev7A1hw1XFvlXcMthU40xmyTV7uWzijfXkkdfuAD4IhgBiASGV/SHeUsRwXoUP4AXMoaxy3ZITqjJq5Yrd9mSa/leLxXXH8gFIIZgEj4zXpS6WqOZqteef3Lctmj0tU+HpB0bIrntit3eOxV8VAHIAOd/wFEpUveF/MD+Q7w+yE1yessfqukB7Kb7Pw+Vlv1YWhYHaZTeUaH9mOSVks6FvRtyxqMEJQzb3n8Dv2t8sJQcJystTtdyzMFT0h6zu/DtV1eZ/+0/7olC725+sL5772UuzYs+5hcAXFAXmgD4IhgBqAU7s34gl4tL7y8VqQfVKe8AQC9kmSMeU5eSBo/xw9OT2V2PjfGtGZ3NM/HD1K7rLWrMx47aIyRtXa3H2w6/cEInZK68pUnCIiZndn954+0Rshau9cYs90vy3P+6/b6f1eUgVDy/vb0NF+nqVSFAWYDghmAUnguGJXpB5hd8mqcCsoaTfiaMkJZ8Lz6cJTj+DnGmG5lhbg8dmWfL68GqksT+2cNyBt9WKg8rf4tsywHMgJpZKy1O40xu+WNwuyU12+ryxizOSu07souo8Nz5wyWfijeJqnoiMwCCGVASAQzACXlh5tOY8wZY8yuArVm2c1j6cw7fjBoldf3KdsLKhLMMs7/WdaubuUOLwXLE4Qwv3btWUkvWGv3ztSUH37t4G7/FoyE7DLGbAnKUOKRmk/JmwalWHNp0ISZL4RV5AhcICp0/gcQlQPyanjySRc5v1DfJUlKZQwwKHR+hzFmS3CTV9uUK8AUK48krZIXjDrk9fs6E/V0EH4/vEn85sUDiqAp1a95eyJ7bro85Uj7/5nKsTslr58hAEfUmAGIUsoYkyrWFyyP4Au9SZNDU0qaEAoKnX+gFJ3kjTGtfm3gdknbMyaTfUrTa+4rplPS3jz7XlBGk/EUmjLT2RPA+jVxL2Su2mCMaS/yHnZLWpzj8SZ5I0cBOCKYAYhKUNMVLNmkMCHNWtttjEnLq53Kbi68VfnDSvb5n1DWlB3GmA6X2qAs7X5A2es/f1rSVr9pM0pfMMZsz/O+rZYXzuSXaVpNmX7tXDp7KS1516BQMNuj3DV37SowKhfAZDRlApiO9qxtpiAwdEgTpsYIpLKOz74veasHbM1ssvSbDl2/8O/Vh7Vbwfkp5Q4RLuV5JPOO/1z5Al5KefpdGWNS/ujQgjPqZ3gqu8nUv/+JUvVx858vGIG6JePWlXVcUPbMJtbdklozB0L4f9vuUtRWArMJa2UCCC2jGS+d8XBK3vQX6YzjtsgLR8/5D/1M3pf/FnlNjXuttdv9L/9N8mrX9srr3xTUsgWBIZjHbLFyLFFUoKyZ56clKWseM6fy+P8dzF0WvHZKXvjI/Ju7/MeD0LXbP368zP77d1xeM2vBtSSDART+exlM+5GSV7NVskXC/Zq/VJ7d9wa1aBlln7DWacZC6yxiDkwDwQwAyiRzRCUASDRlAkA5pcpdAADxQjADgDLwm1HpfwVgAoIZAJTHVEaGAqhw9DEDAACICWrMAAAAYoIJZkMwxhhJH5E0XO6yAACARKmX9K4t0lRJMAvnI5JOlrsQAAAgkVZIOlXoAIJZOMOS9M4776ihoaHcZQEQd8PD0qFD0s03S/X15S4NgDIZGhrS1VdfLTm0uBHMpqChoYFgBqA4Y6SFC6WGBoIZACd0/geAqMybJ7W0eFsAcECNGQBEpaZG+uhHy10KAAlCjRkARGVsTDp3ztsCgAOCGQBE5dw56bXXvC0AOCCYAQAAxATBDAAAICYIZgAAADHBqMxpGB2zevX4gPqHR9RcX6PbVjVpTpUpd7EAxIUxUlWVtwUAB4kOZsaYlKQvSLrXWtvpeM42SWn/bspau3Mqr72/p0+PPX9YfYMj44+1NNZox8Y23bOuZSpPCaDS1NVJn/pUuUsBIEES25RpjGmXF8pSkpocz9kmSdba3dba3ZK6jTG7wr72C4ff04NPd08IZZL03uCIHny6W/t7+sI+JQAAQHKDmbW22w9XvSFOe0TS7oznOCBpS9jX/rPv/YtyLQ0fPPbY84c1OlZw8XgAs8H589LPfuZtAcBBYoNZWMaYVnlNl+kc+zrynFNtjGkIbpLqJen00MW8r2Ml9Q2O6NXjAyUpN4AEGx2Vzp71tgDgYNYEM0mteR5Py2sOzeURSYMZt5OuL9Y/PFL8IAAAgAyzKZjlM6D8fdSekNSYcVvh+qTN9TXTLxkAAJhVEj0qs0TyDhyw1l6UNN5uafwh70sbqvX+ReXsZ2YkLWv0ps4AAAAIYzbVmOUbJJAqsC+nr/z2DZK8EJYpuL9jYxvzmQGQamulj33M2wKAg1kTzKy1vZLS/iCA7H0HwjxXZ9syPXl/u5Y1TmyuXNZYoyfvb2ceMwCeuXOlq67ytgDgoBI+LXK2GfoBbFPWBLJPSOqQP2WGMWaTMqbPCOOedS3qbFvGzP8A8rt0STp9Wlq6VJo/v9ylAZAAiQ1mQfCStFlSuzGmS9Jr1tq9/iEdkrZKGg9m1tqdxphtfiCTpFuttVunWoY5VUZ3rl481dMBVLqLF6Vjx6RUimAGwImxlolQXflzmQ0ODg6qoaGh3MUBEHfDw9LBg9L69VJ9fblLA6BMhoaG1NjYKEmN1tqhQsfOmj5mAAAAcUcwAwAAiAmCGQBEZe5cafFiRmUCcManBQBEpbZWuvHGcpcCQIJQYwYAUbFWunzZ2wKAA4IZAETl7Fnpxz/2tgDggGAGAAAQEwQzAACAmCCYAQAAxATBDAAAICaYLgMAolJXJ23YIM2ZU+6SAEgIghkARMUYJpcFEApNmQAQlQsXpDfe8LYA4IBgBgBRuXJFGhjwtgDggGAGAAAQEwQzAACAmCCYAQAAxATBDACiUlMjrVnjbQHAAeO4ASAq8+ZJy5eXuxQAEoQaMwCIyuXL0unT3hYAHBDMACAqIyPSW295WwBwQDADAACICYIZAABATBDMAAAAYoJgBgBRmTNHamjwtgDggOkyACAqCxZI7e3lLgWABKHGDAAAICYIZgAQleFh6Qc/8LYA4IBgBgAAEBMEMwAAgJggmAEAAMQEwQwAACAmmC4DAKKycKF0++1SdXW5SwIgIQhmABCVqiqptrbcpQCQIDRlAkBURkakt97ytgDggBqziPUPjah/+OKkx5vrq9XcUFOGEgGYMZcvS6dPSytWSDX8/w6guMQHM2PMNklp/27KWrvT4ZwtklL+easlPWGtTRc4ZcqeeeWEvv7ikUmPf/nuNXq48/ooXhIAACRUooOZH8pkrd3t3+8wxuyy1m4tcs7uIIgZY1KSnpJ0bxRlvO/2lepsW6qj/Wf10J5D+trmm3Vdc52a6+kMDAAAJkp6H7NHJO0O7lhrD0jaUuSczszaMf+/UxGUTZLU3FCjtS0NOn/piiTp/KUrWtvSQDMmAACYJLHBzBjTKq/pMp1jX0eBU9PGmBf8mrLgeXojKaSk/T192tD1kh7d1yNJenRfjzZ0vaT9PX1RvSSAuKiulq69lukyADhLbDCT1Jrn8bQK14A94J97xhjTJakjX9OnMabaGNMQ3CTVhyng/p4+Pfh0t/oGJ47Iem9wRA8+3U04Ayrd/PleMJs/v9wlAZAQSQ5m+QxIasq3069h65K0V9I2SfcGtWc5PCJpMON20rUQo2NWjz1/WDZXGfztY88f1uhYriMAVIQrV6SBAW8LAA4qMZjlDWWS5NeS9Vpr75U3IrNJ0sE8hz8hqTHjtsK1EK8eH5hUU5bJSuobHNGrxwdcnxJA0ly4IL3xhrcFAAdJDmb5+oWl8u3L6Jd2QJKstb3W2vXy+p1tyj7eWnvRWjsU3CQNuxauf9htQknX4wAAQOVLbDCz1vbKC1ST+poFwSuHVn0451mmXSUsmiSpud5t1KXrcQAAoPIlNpj5npA0PgLTr/XanXG/NZjrTBoPbO05+pStt9buLWXBblvVpJbGGpk8+42klsYa3baqYMsrAACYRRIdzPxZ/lPGmE1+KLs1a4Rlh6TsEZf3SnrEGLPNGLPFD27bS122OVVGOza2SdKkcBbc37GxTXOq8kU3AIkXLGJeleiPWgAzyFjLqEBX/pQZg4ODg2poaHA6Z39Pnx57/vCEgQAtjTXasbFN96xriaikAAAgLoaGhtTY2ChJjX6f9bwIZiFMJZhJ3tQZe147oUf39ejxz6/T5ltXUlMGAMAsESaYUb8esf6hEb3VN6QF871lSRfMn6u3+obUP8RoTKDinT0r/fjH3hYAHCR6EfMkeOaVE/r6i0fG7z+055Ak6ct3r9HDndeXqVQAZoS10uXL3hYAHBDMInbf7SvV2bZ00uPN9aydBwAAJiKYRay5oUbNDcxVBgAAiqOPGQAAQEwQzAAgKgsWSO3t3hYAHNCUCQBRmTNHCjG1DgBQYwYAUbl4UTp61NsCgAOCGQBE5dIl6eRJbwsADghmAAAAMUEwAwAAiAmCGQAAQEwQzAAgKvPmScuXe1sAcMB0GQAQlZoaac2acpcCQIJQYwYAURkdlYaHvS0AOCCYAUBUzp+XDh70tgDggGAGAAAQEwQzAACAmCCYAQAAxATBDACiYoy3kLkx5S4JgIRgugwAiEpdnfSbv1nuUgBIEGrMAAAAYoIasxjpHxpR//DFSY8311eruaGmDCUCMC3nzkmHD0ttbdLCheUuDYAEIJjFyDOvnNDXXzwy6fEv371GD3deX4YSAZiWsTEvnI2NlbskABKCYBYj992+Up1tS3W0/6we2nNIX9t8s65rrlNzfXW5iwYAAGYAwSxGmhtqJjRZXtdcp3XLG8tYIgAAMJPo/A8AABATBDMAiEptrbRunbcFAAc0ZQJAVObOlZYsKXcpACQINWYAEJVLl6QTJ7wtADggmAFAVC5elHp7vS0AOCCYxczomNUbJ9OSpDdOpjU6ZstbIAAAMGMIZjGyv6dPG7pe0qP7eiRJj+7r0Yaul7S/p6/MJQMAADOBYBYT+3v69ODT3eobHJnw+HuDI3rw6W7CGQAAswDBLAZGx6wee/6wcjVaBo899vxhmjWBpJk7V7rqKm8LAA4IZjHw6vGBSTVlmaykvsERvXp8YOYKBWD6amulj32MecwAOEv8zzhjzDZJaf9uylq70/G8LknH/LsD1tq9ERTPSf9w/lA2leMAxMTYmHT5sjRvnlTF72AAxSU6mPmhTNba3f79DmPMLmvt1gLnpCS9KOlua23aGNMu6aAkMwNFzqm5vqb4QSGOAxAT585JBw9K69dL9fXlLg2ABEj6T7hHJO0O7lhrD0jaUuScLkl7rLVp/5xuSZ1RFdDFbaua1NJYkzcZGkktjTW6bVXTTBYLAADMsMQGM2NMq7ymy3SOfR0FTt0iaa8xpjU4zg90ZTOnymjHxjZJk6vtgvs7NrZpTlXZKvUAAMAMSGwwk9Sa5/G0pFSuHX6Yk6R2/5heY8yufEHOGFNtjGkIbpIia4u4Z12Lnry/XcsaJzZXLmus0ZP3t+uedS1RvTQAAIiJRPcxy2NAUr42vyCYpf0mTBljtks6LmlRjuMfkbSj5CXM4551LepsW6Y9r53Qo/t69Pjn12nzrSupKQMAYJZIco1ZPi4dsX4W/IffFJrKU2v2hKTGjNuKUhSwkDlVRjetSEmSblqRIpQBSVZXJ33qU94WABwkucasN8/jqQL78j2eVo6mUWvtRUnjqw8bQ0gCEIIx3g0AHCW2xsxa2yspndFvLHNfzs78/jm9mhzCUsqoRQOAkjh/Xjp0yNsCgIPEBjPfE5LGmyCNMZuUMX2GP/JyW9Y52yVtzjrnQNDnDABKZnRUSqe9LQA4SHJTpqy1O40x2/xwJUm3Zk0u2yFpq6SdGefsNcY0ZQS2xdbass5jBgAAICU8mEleOMu4uzdr325l1KBlPQ4AABAriQ9mlaR/aET9wxd1tP+sJI1vm+ur1dzAckwAAFQ6glmMPPPKCX39xSPj9x/ac0iS9OW71+jhzuvHHw8CXDYCHBAzNTXSRz/qbQHAAcEsRu67faU625ZOery5vnrC/ewAF8gOcADKbN48qYVVOwC4M9bacpchMfxlmQYHBwfV0NBQtnJkNnk+tOeQvrb5Zl3XXEeNGRA3ly9L778vLVnihTQAs9LQ0JAaGxslqdFaO1ToWGrMEqi5oWZCALuuuU7rljeWsUQAchoZkX7xC2/mf4IZAAdJn8cMAACgYhDMAAAAYoJgBgAAEBMEMwCIypw5UirlbQHAAZ3/ASAqCxZIN99c7lIASBBqzAAgKtZKY2PeFgAcEMwAICpnz0r/9E/eFgAcEMwAAABigmAGAAAQEwSzhBods3rjZFqS9MbJtEbH6MMCAEDSEcwSaH9PnzZ0vaRH9/VIkh7d16MNXS9pf09fmUsGAACmg2CWMPt7+vTg093qGxyZ8Ph7gyN68OluwhkQJwsXSnfe6W0BwAHBLEFGx6wee/6wcjVaBo899vxhmjWBuKiqkqqrvS0AOGCC2QR59fjApJqyTFZS3+CIXj0+oDtXL1b/0Ij6hy9OOq65vlrNDTURlhSAJOnCBam3V2ptlWpry10aAAlAMEuQ/uH8oSzXcc+8ckJff/HIpP1fvnuNHu68vqRlA5DDlSvSr38trVxZ7pIASAiCWYI017vVcgXH3Xf7SnW2LdXR/rN6aM8hfW3zzbquuU7N9dVRFhMAAEwRwSxBblvVpJbGGr03OJKzn5mRtKyxRretapIkNTfUTGiyvK65TuuWN85MYQEAQGj0SE2QOVVGOza2SfJCWKbg/o6NbZpTlb0XAAAkAcEsYe5Z16In72/XssaJzZrLGmv05P3tumddS5lKBmCS6mqv43813QcAuKEpM4HuWdeizrZl2vPaCT26r0ePf36dNt+6kpoyIG7mz6fjP4BQqDFLqDlVRjetSEmSblqRIpQBcXTlivT++94WABwQzAAgKhcuSD093hYAHBDMAAAAYoJgBgAAEBMEMwAAgJggmAFAVKqqpIULWcQcgLNpT5dhjPkfJG2W1CqpV9L/Z639T9N9XgBIvIULpVtvLXcpACTItH7GGWOelfRX8iaeP+5v/9wY848lKBsAAMCsMuUaM2PMn0jaY639Qo59Dxhj/sRa+xfTKh0AJNnZs9Lrr0u33CLV1ZW7NAASYDo1ZoPW2r/NtcNa+5QmL+eIMhgds3rjZFqS9MbJtEbHci1/DiAS1kqjo94WABxMp49ZsU+aM9N4bhTQPzSi/uGLOtp/VpLGt8311Wpu+HANzf09fXrs+cPqGxyRJD26r0ffeOmodmxsY01NAABiaDrB7Lpp7scUPfPKCX39xSPj9x/ac0iS9OW71+jhzusleaHswae7J6Xn9wZH9ODT3RMWPA+CXrbsoAcAAKI1nWC2x+/k/2eSDlprh4wxDZI6JD0i6YFSFLAYY8w2SWn/bspauzPk+S9YaztLXrAI3Xf7SnW2LZ30eHN9tSSv+fKx5w/nrNK08tqYH3v+sDrblmlOlZkU9AKZQQ8AAERvysHMWvu6MebPJT0laZUx413K0pK2WGsPTbt0RfihTNba3f79DmPMLmvtVsfzN8kLkonS3FBTsCbr1eMD482XuVhJfYMjevX4gO5cvXg86B3tP6uH9hzS1zbfrOua68aDHoApWrBAWr/e2wKAg2nNY2atPSDpOmNMu6T1knqttS+WpGRuHpG0KrM8xpgXJBUNZsaYlKSm6IpWPv3D+UNZruOyg951zXVat7wxkrIBs8qcOVJ9fblLASBBSjIdtbW221r71EyGMmNMq7ymy3SOfS61YF+Q9GyR16g2xjQEN0mJ+IRtrnfrF+Z6HIApGhmRjhzxtgDgILJ1QowxT0b13L7WPI+nJaUKnegHtwMOr/GIpMGM20n34pXPbaua1NJYk3e+EiOppbFGt62qyApDID4uX5ZOnfK2AOBgWk2Z/nJMudYbScmrkXpwOs8/RQMq3kSZstb2+s2ZhTwh6d9n3K9XAsLZnCqjHRvb9ODT3TKaOK9JENZ2bGzTnCqmmgMAIE6mXGNmjPkzecsxrZe0OsetXAqGMmPMFmvtXpcnstZetNYOBTdJwyUp4Qy4Z12Lnry/XcsaJzZXLmusmTBVBgAAiI/pdv7PG4L84Bal3jyPp/Lt8wcp/CyqAsXNPeta1Nm2THteO6FH9/Xo8c+v0+ZbV1JTBgBATE0nmB0rsv+JaTx3UX5TZNoY02qt7c3al6//WJOk9ozBAaul8Wk3el1r0pJkTpXRTStSkqSbVqQIZcBMmj9fWrHC2wKAg+n2MWvwm/hyuVdeU2eUnpA3D1kwj9mm4L/9+62SNgWTzvqB7UDG/nZ5c66FmpQWAJxUV0vXsQgKAHfOwcwY85msh45J6jLGpCW9luOUrYo4mFlrdxpjtvmBTJJuzZpctsMvx6Tg5Z+z2f/vLkkvFKhpA4DwRkelc+ekhQu9Oc0AoIgwNWZ75fXfSufYl2tC1xmZoTSrtmtv1r7dyqhBy9q3N/t4ACip8+el7m5v9n8mmgXgIEww+5m19rOuBxtj/nIK5UGMsdg5AADRChPMtod87l0hj0fMsdg5AADRcg5m1trXM+8bYx6Q1G6tzTmJbPbxSD4WOwcAIFrTGZXZqfxziaECsdg5EJIx0rx53hYAHExnrczXrLVfybfTGBPpPGYAEHt1ddInP+ltAcDBdGrMnjPG/In/393y1qjM1CFvEXAAAAA4mE4wK9aMaYvsR4yMjlm9cTItSXrjZFprWxpYJQCYrnPnpJ4ead06by4zAChiOk2Z3ZIWWWurct0k/W2JyoiI7e/p04aul/Tovh5J0qP7erSh6yXt7+krc8mAhBsbky5c8LYA4GA6wWy7tXawwH6my0iA/T19evDpbvUNjkx4/L3BET34dDfhDACAGTTlYGatfXE6+xG9/qER9Zwa1NH+s5Kko/1n1XNqUP1DXggbHbN67PnDOducg8cee/6wRsdolQYAYCZMaxFzxFv2hLAP7Tkk6cMJYV89PjCppiyTldQ3OKJXjw/oztWLIy4tAAAgmFWwYELYbMGEsP3D+UNZJtfjxo9n6SbAU1sr3XSTtwUABwSzCpY9Ieyk/fVuIcn1uABLNwG+uXOlpqZylwJAghDMZrHbVjWppbFG7w2O5OxnZiQta6zRbavCfbGwdBPgu3RJevdd6SMfkebPL3dpACTAdEZlIuHmVBnt2NgmyQthmYL7Oza2hZ7PrLmhRuuWN+q6Zm+282DpJpoxMetcvCi9/ba3BQAHBLNZ7p51LXry/nYta5wYmpY11ujJ+9t1z7qWMpUMAIDZh6ZM6J51LepsW6Y9r53Qo/t69Pjn12nzrSuZ+R8AgBlGjRkkec2aN61ISZJuWpEilAEAUAYEMwCIyrx50tKl3hYAHNCUiVBKudg5852h4tXUSGvXlrsUABKEYAZn+3v69Njzh8dXC3h0X4++8dJR7djYNqVBAsx3hoo3NuaNyKyulqpooABQHMEMToLFzrPnOwsWO5/KCE7mO0PFO3dOOnhQWr9eqq8vd2kAJADBDEUVW+zcyFvsvLNtWahmzeyVCYL5zgAAmK2oW0dRYRY7BwAAU0cwQ1FRLXYOAAAmIpihqKgWOwcAABPRxwxFRbXYeRhMrYFEqq+XPv3pcpcCQIIQzFBUsNj5g093y0gTwtl0FjsPg6k1AACzAcEMToLFzjPnMZO8mrKpzmMWBlNrIJHOn5f+5V+kG26QFiwod2kAJADBDM7Kudg5U2sgkUZHpaEhbwsADuj8j1BcFzvPXrppdCxX7zQAAJCJYIaS29/Tpw1dL+nRfT2SvKWbNnS9pP09fWUuGQAA8UYwg/qHRtRzalBH+89Kko72n1XPqUH1D4WflyxYuil7Qtpg6SbCGQAA+dHHDJNGPD6055Ck8CMeo1q6KQym1UCs1NRIa9d6WwBwQDDD+IjHbGFHPIZZuunO1YvDFtMJ02ogVubNk5ZO/n8LAPIhmGHSiMepisPSTUyrgVi5fFnq75eam72QBgBFJD6YGWO2SUr7d1PW2p2O50jSakmy1m6NpnSzSxyWbnKdVoMmT8yIkRHpyBGpoYFgBsBJooNZELCstbv9+x3GmF2FgpYxpstauz3j/i5jzAvW2s7oS1zZ4rB0kyuaPAEAcZT0UZmPSNod3LHWHpC0Jd/BxpiUpHZ/G9glqcMY0xpRGWeNYOkm6cOlmgL5lm4q13xn992+Ut/94w362uabJUlf23yzvvvHG3Tf7Stn5PUBAMglsTVmfpBKWWvTOfZ1+CEtl09IapXU7d/v9bepUpdxNgqzdNP+nr4Jxz26r0ffeOnojCzxRJMnACCOEhvM5IWrXNLKE7L8ELco6+EOf9ub9biMMdWSMnuN14cp4GzlsnRTMN9Zdv1YMN/Zk/e3Rx7OXNDkiWmZO1dqavK2AOCgEj8tBiSF6cT0iKStuWre/H07SlGo2abQ0k1xmO/MlesoT2rWkFNtrXTTTeUuBYAEqcRg5hzKjDFdkvYEgwdyeELSv8+4Xy/p5DTKBsVjvjNXrk2e1KwhJ2u9BcznzJFMeX9kAEiGJAezSU2PvlSBfeOMMZskHSsQymStvSjpYsY5IYuIXOIw31mpUbOGnM6elQ4elNavl+rpCQGguMQGM2ttrzEmbYxptdb2Zu3L1/Ffkjc4wD8umGYjJakp+3kQjanOd5Y9gnNtS0PZmzoD1KwBAEohscHM94S8zvtBwNqkjOkz/JGbmzInnTXGtEtql7Q3Y4qMCechWlOZ76ycIzhLiZUJAACFJHoeMz9wpYwxm/xQdmvW5LIdksbv+zVjL0rqknQs49aVp/M/IhB2vrNgBGd2v7RgBOf+nr6IS1w6zQ01Wre8Udc110n6sGaNZkwAgJTwYCZ54cxau9e/bc/at9tauzrjftpau8haa7JvM1/y5OkfGlHPqUEd7T8rSTraf1Y9pwbVPxS+L1gw39myxomBZFljzYSpMoqN4JS8EZwzNTEtAABRSnpTJmZQdv+oh/YckjT1/lEu850laQRnqTFQoALU1Umf/CTzmAFwxqcFnAX9o7JNp39UofnOpMocwemKgQIVwBgWLwcQCsEMzrJHHs7Ia05hBGecR2+GwUCBCnDhgnT0qHTddd5kswBQBMEMsRZ2BGeljN6U3KfgQIxduSJ98IF07bXlLgmAhEh8539UtjAjOCtp9CYAYHYimCH2XEZwTmX0ZnaTJyM7AQRGx6xePvaB/v7QKb187AM+HzBjaMpEIhQbwRl29GYlNXkCKK3szwdJamms4fMBM4JghsQoNIIzzOjNoMkz+/dv0OSZOY9aEjCtRoxVV0urV3tbJELYz4fRMatXjw+of3hEzfVef9ckDjZCfBDMUBFcR28uWVitP9n7z3mbPI28Js/OtmWJ+XBlWo0Ymz9fuvrqcpcCGQoFqWJdIrI/H6hZQxQIZqgIrqM3ZRR6wtq4T7/BtBoxduWKdOaMtGgRk8zGQLEgFaZLxOCFSxVV8474oPM/Sq6USze5ch29+f7ZyU1+uQRNo/t7+rSh6yU9uq9HktcXbUPXS7Ea4cn6mzF24YL08597W5SVy6ht1y4R7w1eYKk4RIZghpJ75pUT+tw3fjS+ZNNDew7pc9/4kZ555USkr+syejPMhLVMvwFUBtdR20vq3GqZB85dcq5ZA8Kibh0lF8XSTa6Kjd50bfJcf80i3fXn3w/VFy3uTZ7AbOXaRCkrp8+HJscAV4lLxSF61Jih5IKmtezbTDWtFRq96drkefBXZ0L9Ik5CkycwW7kGpPfPXXT6fFjm+FnmWkMPZCKYYdZxafKcyvQbrk2eTGw7i8yZI9XVeVuUTZguDC6fD0HNe776cCOv5i1YKk5iwlq4oykTs1KxJs+opt9gYttZZsEC6ROfKHcpZr2wa+4Gnw/5ptUIat4ffLpbRprwnNlLxUlMWItwCGaYtQo1eUYx/Ua5h9czEe308R4mU9ggFZwTTJmTS1Czlh24lmUFrkqb0BrRI5ihbIIvucxpNaR4fMm5fpC7Tr/x3uAF7fzHX5R1Ylsmos3PNXCFfg/PnpW6u6X2dq9JE2XjGqTCPmehmrWwE9YCEsEMZZT9JRdMrxGXoODyQf7ysQ+cnivM8PrMX+mlHOkZ94loo6iNKnXgcn0Pg9etOjus2tNDunAqrbG60UmvSw3czCoWpKaiUM1a2DV8AYlghjIq57Qarko1/cZUhteXuj9ac0PNhC/7YCLauIiiRq/Ugcv1PQxed+HF87rx9FG9+ZMLOle9YNLrUos584o1UZZSmEFEQIBghrLJ/pKLK5fpN4o1eTbWznd6rWDQQRL6pbjW9rgeF6ZGr9TPWerQGrzu28ff09P/8ai6/sebdO2qZZNeN+61mJieMKNBgQDBDJgmlybP0THrPCosKf1SXGt7XI8LE46ieM5SCl636uywJKn1qoVqy/G6ca/FxPSEHQ0KSAQzJEScBwpIxZs8w4wKe/nYB2Xtl1Lq2qgoaoWSUtM0tmCh/nnZ9RpbsLDcRal4o2O2pH3HSmEqo0EBghkSIe4DBaTCTZ6S+6iwcvdLKXVtVBS1Qompaaqq0oX5NVIVc3lHKc7zhEUxGhSVjWCGREjCQAEXxWrWpPL3S0lKbVQSmJERtX5wUmZkRFIMg2MFSEJ/zDCjQeNY84eZRTBDIiRloICLYjVrU+mXUsppNRJTG5UA5splNZ8bkLlyudxFqUhJ6Y8puY0GjXPNH2YO9etAzLgutJ653AsLqGM2CjNPWNyFXXMXlYtghorSPzSinlODEwYJ9JwaVP9QsuYJcllIWeLDHLNbuftjlkqxmj/Jq/lj4fP4K8Vi9TRloqIkYZCAq2L90ZLUjANEodz9MUuFFQIqQ6Gm6N9Y6T4ym2CGilIpgwQChfqj8WEef3befJ1qaJad5zbBMMKplHnCKqXmbzYrNgjlL/77Nc7PRTBDRXEdJBD3edFc8GEef7a6Wu+klslWJ/OHQdxVyjxhU635YwRnPLi0XvzZ9/7F+fkIZpiVKqHJs1KacSra6KgaRs5Ko6PlLknFqoR5wqZS88cIzvhwab04PTR50u58CGaYlSqhybNSmnEqWdWF82rr71XVhfOSuA5RCTNPWByFrflLwtxts0mpWyUIZpiVKmFetEppxgFKwWWesDhzrfmbyqAfmjyjVepWCYIZkGCV0IwD5DPbAoVLzV/YQT9hmjxn2/tdKi6tF0sbqvWO4/MRzICEc1nmCUia2dqHqljNX5hBP2GaPGfr+10KLq0XX/ntG7Tp37k9HxPMAgUkZcLaYss8BbKXbmLCymjZqipdmjNPlkXMQ2Hi5Pxcm82WLKx2nrSW93v6ik0K3tm2zPm5El9jZozZJint301Za3dGcQ5mp0oYvRnI/kX86L4efeOlo/wijpBdsFDdy9fKLnCfXHK2Y+LkwlwH/cjIqcnzp8c+oM9aiRRqih4aGnJ+nkQHMz9gyVq727/fYYzZZa3dWspzMHu5jt6M+7xojOJCUjBxcmGug37eP+s2PcPLve9H1mdtNirFIJSk168/Iml3cMdae0DSlgjOwSzV3FCjdcsbJ92yw9Yzr5zQ577xo/EatYf2HNLnvvEjPfPKiTKUeiLW4XPj2swbpjnYnD+n9lNvyZw/V+riJlqh9QSZOLk4l7V03UcKutV0ZfZZm41NnqVYA9NVYmvMjDGt8poh0zn2dfiBa1rnGGOqJWVWjdRL0qFDh1RXVzf+4KJFi7Rq1SqNjIzo8OHDk8ra3t4uSfrFL36hc+cmfkBfe+21ampq0q9//Wu9887EMRv19fVas2aNRkdH9c///M+TnvfGG2/UvHnzdOzYMQ0ODk7Yt3z5ci1dulRnzpzR8ePHJ+yrra3V2rVrJUmvv/66rJ34D2zt2rWqra3Vr371K33wwQcT9i1dulTLly/X8PCwjhw5MmHfvHnzdOONN0qS3nzzTV2+fHnC/jVr1qi+vl6nTp3S6dOnJ+xbvHixrrnmGl24cEFvvfXWhH3GGN1yyy2SpLfeeksXLlyYsH/VqlVatGiRTp8+rVOnTk3Y19jYqNWrV+vy5ct68803le3jH/+45syZoyNHjmh4eHjCvquvvlpXXXWVBgYG9Pbbb0/Yt3DhQn30ox+VJHV3d+vG6ova+ekP/02sXnODqmtqdOGDPnV3d48/PnD2oubUNWnQ1Gns4nn94w9f1uGmWjUtmK+mumpVV1frYx/7mCTpjTfe0JUrV3S0f1gX3zuqw2/W6drG9aqrq9PJkyfV398/oUwfXPGW/blw4by6u4+NP97Tf9H5F/GqhZf17rvvStL4677zqxVat/zjunTpknp6eiadP/eqVZKkt3uP6tLpiR/yK1eu1JIlS/T+++/r8Js/H/87Lp2uV11dna6//nqNjY3p0KFDkrwPvwNH0rqcfk+vv/2+1rY06FdvH1c6nZ7wvB/5yEe0bNkyDQ2mJzynJNXU1KitrU2S9//q5Suj48/5/D/9TNdsvEP1dQt14sQJvf/++5Kkn568oP/0+pA+uDAmyWvm/er+w/riLQ26Y0Wt93fOnat3q66a1BwcHHf/p29UQ0OD+vr61NfnfUH1vn1aZ3/9tk6deFsfW/ORvJ8R85euliQdP3ZEl05P/L1caZ8Rp8ySSTUui2ur9MVbGvQH/91NzoFi8L0T6u72PkeS8BmRra2tTTU1NTp+/LjOnDkzYV9LS4taWlo0NDSko0ePTtgXfEbcs65FS6/06833zuvMyJgW1VRp7ZL5Wnut9//BR+ad1+LaqvF/07ksrq3SktEP8u7PVGMv60+fP1a4yfM7P9fikXcnNWvefPPNqqqq0i9/+UudPXt2wr7Mz4gTJyb+kM31GZFp3bp1mj9/vnp7e/N+RqTTafX29k78W7I+I8bGJr5HN9xwgxYsWDD+GZH9+SBJS+vn6w9vXDD++SB5nxE33XSTJOnnP/+5Ll78sNYy++8uyFqbyJukDq/4kx4/I2lTKc6R9Kfy/s0VvN13333WWmuPHDmSc3/gjjvumLTvW9/6lrXW2m9+85uT9n32s5+11lo7ODiY83n7+/uttdZu3Lhx0r6vfvWr1lprn3322Un7brnllvEyzZ8/f9L+np4ea621X/ziFyft+8pXvmKttfb73//+pH3Lly8ff97ly5dP2v/973/fWmvtV77ylUn7vvjFL1prre3p6Zm0b/78+ePPe8stt0za/+yzz1prrf3qV786ad/GjRuttdb29/fnfA8HBwettdZ+9rOfnbTvm9/8prXW2m9961uT9t1xxx3jZcr1vEeOHLHWWnvfffdN2tf4yd+312z/rm2+97FJ+1avXj3+vEuWLJm0/yc/+Ym11tqHH3540r7Nf/iv7TXbv2v3fO8HEx5fsPZT9prt3y16+7vXT9rHH3980vN2/s7vWmutfeedd3L+rQePnbbXbP+u/cQdn5y076mnnrLWWvvUU09N2nfXXXdZa60dGRmxkmzt9Xfa5Q/+9YQy3fH4AfuZP/zfJ537+OOP2yujY/bffHOPXbD2U7b66hutTJWVZNva2sbfw6aPf2bSc7b/6T/Y7735rv3Sl740/rortz1vV257fsJxK7d9x67c9rytvf5OK8kuXf9b9toc71tw3OPf+gdrrbU7duzwymmqbPXVN9oFaz9lP735QXtldCz3Z4Spss/89G17zfbv2ra7N43/HYn8jDBV9n/e/oT9u9dP2v+498Ckv2XFHZ8r+B7++bdfsFdGx+wN2//Wrtz2nbz/Vpc/+NcTnrvSPiN27NhhrbV2//790/qM+PDf9uT3cvzftqmyyx/867zv98pt37HLH/xru/O//L3T50j11TdOKtPIyIi11tq77rprWp8R2bd33nnHWmvtpk2bcn5GWGvt3//930/al/kZUV9fP/kz7eBBa621X/rSl/J+Pkx4D/3zlixZMv68q1evzllmSQ22SL4xNuuXUFIYYzokvWCtNVmPH5PUZf0+ZNM5J0+N2ckf/vCH1JhRYzalX8NBjdlVS5fp7PCQTrzt/ZIrVmP25b85pK//3s36VxsK15g98P/+Ss/963bNGfqwSaGn/6L+zx8MTCpjtm8/cMekGrMv/80h/ectn9Zv/0buGrPRMatfXlmsf/t3P9f/dsci/eaK6gm/loNfw6f7f63//OIb+uYP39b/ete16liTUmND/fiv4b/87sva+ZP0pDIFfWi2/UZqwi/TX5yv1f/9397NWevy6esWqa2tTft7+vRHT0++NkHp/t1vrdQN9Zf0R/+1v2itwv/zr67S//K99/XB+fxLKy1rqNaPv3K3+k+/p32v9U76hd3SWKNHfmuNVuj98cdy/RIP/o7g703KZ0Sxv2V0zOqP/uHXTu/h3/y3w/q333t70v7g2v2brH8PlfQZIbnVmEkffkZkuv766yd8RuS6LssaqvU/rasdfw9/evKCdv4kPanPWmDbb6S0aEmzHvnOL3PsnejhO1L6zZW1Ex5Lao3Z8bd/pU3/5XDRz4e//J1mzakyRWvM7rrrLklqtNYWHAlQicHsjKTtIYNZ3nOyjmuQNDg4OKiGhoZp/w2Ai55Tg/rcN36k7/7xBq1b3hj6uNExqw1dLxUdxfWj7Z+ZEKqKva5rJ+BixwXly9fcml2+fAMZgpIHQ9NdnvMv7v247vurV3Iek+n/+J21+r/+61tFj/v2A3do8MKlouW7Z12L09+RlM7ULn9LY+18/f5TPy36XN9+4A46mZeYyyjKYu/3y8c+CHX9XF83zqbyN+czNDSkxsZGySGYJbaPmaTePI+nCuybyjlAWZRqpOdUlm7K7uC+tqVhwn7XUZ4uxzXWznfuA3fbqianof31NfOcnvPlY279a341cN7puPcGL2jnP/6iaPk+c8PSipkSwnV6i2333OD0fEGn/qSvfxknLiMFi73fYdfmrYRgXa6BKIkdlWmt7ZWU9jv0Z++b1PF/qucA5VLKkZ4uo7gC+3v6tKHrJT26z2u2fHRfjzZ0vTQ+4sp1lOelK2NOx73nOFlv//CI81QKroErd8PNZNc0LXA6buDcJafyfevlt53DaKaZHBnmyvWaDDhO35DZ+T8IFL9783LduXoxoSxihd7v4AeeNHkcZ/YPvEoZvek6EIW1Mid6Ql6H/mBOsk3KmArDD2Cb7MQJZAueA8SF6xxqrlyWbiplDZdr+Ajzhe3+y9QtsNzZukR/232qaC3AH9x5rf7qR8eLHtdU53ZtXGvgMv/euNZAuF6TpoXzQ9W4IH5c1uatpAmCw9YSlkpia8wkyQ9cKWPMJj9g3WonThTbIWlryHOAWHCdQy2MQks3udaEudZwuYaP4As730e0kRdAblvV5PzL9M7WJU7PecfqxU61APPnVjkdt8zx2rjWwAV/b5xrIFyvybLGWucaF8TXPeta9KPtn9G3H7hDX/+9m/XtB+7Qj7Z/ZvzHQZgJguMuTC1hKSU6mEle0LLW7vVv27P27bbWrg5zDjBblbpJyjV8hPnCDn7Blipwzakyzs28Lse5lu8P7rzWOYzGfYJg17/5tlVNoZrUEV+Fmjyn0i/LtYk+iqb8Ys9Zjn+zSW/KBFAipW6Scm3+CzoYF2sikcINZHBpdgm4NPO6HOdavqAGzuXvePnYB2VfoqjQ6Lqwg0vo1F/ZwvbLKtXo7qlwfc6Z/jeb+BozAKVR6iYp1+a/zC/sH23/jB7//DpJ0uOfXzehiSQQ5hes63NKhZt5wxxXyho4qfxLFAWDQX7/qZ/qy39zSL//1E8nDAaRwtcq0Km/coWpQXVtoo+iKT/sc87kv1lqzABICtfR1bWGK0ytleQejlxruMI8ZykF5Xv2x0f1H/7mZT38e3fqC5+8LnQNnFS+kWGS+7QoEjVh8LjWoEpyGiQw1allCtXyxn2AAsEMSLhyzXdWqua/qSpH4ApjTpXRja3N6q9frBtbm6dcAzeVkWGlmNhzKl9eLvNlofK5/CBzbaIPM7VM8G+vWBNlmAEK5fj3TDADEu6ZV07o6y9+uDxWMO/Zl+9eo4c7rw/1XFHVcMU9REVmbEy1l0aksfxLuhQTNjCH6YtTKMDF/csL8VasBtW16T3s1DIutbwXr7j9/xhV94BiCGZAwpVjvjO4qTp/Th9/75eqOn+7pEVTfh7XwBym6bFYgCt33zYkX6EaVNem9zBTy7jW8v7FvR93fs5yIJgBCdfcUDOtuc1ymbU1XDFWrAYiTNPjC4ffKxrgytm3DZXPtYk+zOhu11peWcV6smNGZQJAQhQaGeb6pfTTYx84zYu2/ppFzqPrgLBcJ28NM7rbtfb2/XMXYz3ZMcEMACqA65fSy73vOwW4g786E+svLyRfqaeWCVPLG+fJjmnKBICoGKMxUyWZ6MOLe5OiW1n6h0f0uzcvDzUYBAjLdZoVl+PCjmCO6xQvBDMAiMjYwjq9evU6jS2si/y1XL+U7ly9WN/8/tGizxcEvbh+eaFyuE6zUuy4sCOYw7z2TKIpEwAqgGufnTtaF4fuO8ZM/UiKODdRuqLGDAAiYi6c1419R2QutEtqjPz1XKfVCFurACRJ0mt5CWYAEBEzOqqFly/IjI7O2Gu6fCmFnUgYSJo4NlG6IpgBQIVx+VJKeq0CUKkIZsAsUao1NVE5klyrAFQqghkwS5RyTU0AQDQIZsAsUeo1NVHcWE2tfrnkGo3V1Ja7KAASgmAGzBJRrKmJIubO1cCCRmkuH7UA3PBpAWDWi6r/nbl0SS1Dv5a5dKkk5QRQ+QhmACqWa+By7X8XNsCZSxd1TbpP5tLF0v5hACoWwQxALIQJPaUOXK797xhAASBqBDMAsRAm9JQ6cLn2v2MABYCoEcwARMq1ditM6Cl14HJV6udjbjkA2QhmACZwDQulbk4ME3riPsI0eG96B0Z0prZBRwdGNHZqcMrvDYDZw1hrix8FSZIxpkHS4ODgoBoaGspdHCAS/+GFX04IC4HssOB6XBBSslVyrRDvDYBMQ0NDamxslKRGa+1QoWMJZiEQzDAbuIYFQkV+4++NtdKVK948Zsbw3gCzFMEsIgQzAKEMD0sHD0rr10v19eUuDYAyCRPMqmamSAAAACiGYAYAABATBDMAAICYIJgBAADEBPOYAUBU6uqkDRukOXPKXRIACUEwA4CoGONNlQEAjmjKBICoXLggvfGGtwUABwQzAIjKlSvSwIC3BQAHBDMAAICYSHTnB2PMNklp/27KWrvT8RxJWi1J1tqt0ZQOAAAgnMQGsyBgWWt3+/c7jDG7CgUtY0yXtXZ7xv1dxpgXrLWd0ZcYAACgsCQ3ZT4iaXdwx1p7QNKWfAcbY1KS2v1tYJekDmNMa0RlBDCb1dRIa9Z4WwBwkMhg5geplLU2nWNfR4FTPyEpM4T1+ttUntepNsY0BDdJrEIMwN28edLy5d4WABwkMphpYrjKlFaekGWtTVtrF1lruzMeDkJcb65z5NXKDWbcToYuKYDZ6/Jl6fRpbwsADpIazPIZkNQU4vhHJG3NVfPme0JSY8ZtxbRKB2B2GRmR3nrL2wKAg1h0/jfGbJK02eHQJ7JqvLI5hzJjTJekPcHggVystRclXcw4x/XpAQAAQotFMLPW7pW0N8Qp+ZoeUwX2jfOD4LFCoQwAAGCmJbIp01rbKymdazSlPzozr2BwQMY0GylGZQIAgDhIZDDzPaEPO+8HtWC7M+63ZkwmGzzWLqldUre/v1XeFBsDM1NkALPKnDlSQ4O3BQAHxlpb7jJMmR+8gqbLW7Mmj90iabu1drV/PyXpuHKM2rTWOnUe86fMGBwcHFRDQ8P0Cg8AAGaFoaEhNTY2SlKjtXao0LGJDmYzjWAGAADCChPMktyUCQDxNjws/eAH3hYAHBDMAAAAYoJgBgAAEBMEMwAAgJggmAEAAMRELGb+B4CKtHChdPvtUnV1uUsCICEIZgAQlaoqqba23KUAkCA0ZQJAVEZGpLfe8rYA4IBgBgBRuXxZOn3a2wKAA4IZAABATBDMAAAAYoLO/1MwNFRwmSsA8AwPSz/9qbRmjcS6xMCsFSY3sIh5CMaY5ZJOlrscAAAgkVZYa08VOoBgFoIxxkj6iKTMFYnr5YW1FVmPo7y4LvHDNYknrkv8cE3iabrXpV7Su7ZI8KIpMwT/zZyQdL2sJkkattbSxhkTXJf44ZrEE9clfrgm8VSC6+J0Dp3/AQAAYoJgBgAAEBMEs+m7KOkxf4v44LrED9cknrgu8cM1iacZuS50/gcAAIgJaswAAABigmAGAAAQEwQzAACAmCCYAQAAxAQTzE6DMWabpLR/N2Wt3VnG4sxKxpiUpC9Iutda25ljP9eoDPz3XZJWS5K1dmuO/Wn/LtdlBmT8vyJ516VV0gPW2nTGMVyXMjLGvJD9OcY1mXnGmA5JWyW9IKlXUqek16y1ezOOiey6UGM2RcEXj7V2t7V2t6RuY8yuMhdrVjHGtMv7oklJasqxn2tUBsaYLmvtTv+21X/shYz9XJfy6JJ0wH/ft0sakPRcsJPrUl7GmE2SOrIe45qUR0retdjl347lCGWRXRemy5giY8wZSauyfm1aa63Jfxai4H+gPWKtXZ/1ONdohvm1Ms/Jq8FM+4+1SzooabW1tpfrUh5+OH4h+GXvf7k8Yq1d5N/nupRJRm3mrsz3m2tSHv53yoHM9z1rf6TXhRqzKTDGtMqrukzn2Ncx+QzMNK5RWX1CXjNZoNffprgu5WOt7cxqbrlV0gGJ/19i4AuSns18gGsSTzNxXehjNjWteR5Py6sCRflxjcrA/7BalPVw8GHVKy+05ZIW12XG+DUCKUn3+g/x/0uZ+F/mB3Ls4pqU1xeMMQPyusms9pv/pRm4LgSz0gouIuKLazTzHpG01VqbNiZvTT/XZQZkNJmlJD2Xr6kmA9cleim/iT/leDzXJHrdkmSt7ZUkY8wWY8xz1tp7C5xTsutCMCst/meJP67RDDLGdEna43eQLYTrMgP8ILZbGv+yOSNpVYFTuC4RMsZscfh/IxvXJGJBIMvwrKRdRcJzya4LfcymJvuiBVIF9mFmcY3KzG8uO5bVr4nrUgbGmJQxpivri+WAPhx9xnWZYf6gmJ8VOIRrUib+Z9e4jJrlVs3AdSGYTYGfptN+J8Dsfbn6CmCGcY3KK+gEG9QG+MGgletSNq2Stmnir/qUv01zXcqiSVKHMWabP0K2S/JGyxpjNnFNyiMYWZ75vmf8oOmdietCMJu6J5Qx54yfsMNWSaM08lUhc43KwK8JaJc3t0+r/wG2RV4fDInrMuOstd2SdmY10WyW1J3xZcJ1mUHW2gMZ8/3tlDdflvz7wZxZXJMZ5teOZf+/skXS3oyas0ivC/OYTYP/Kye4eLdmjNrADPC/8DfJ+4Jpl7RTuWdn5hrNEP+X5XHlGJ2UNT8T12WG+ddmS8ZDqyVtzzHzP9dlhvlf7JvlfZ7tlDffXDCVCddkhuX4f2Vx9vse5XUhmAEAAMQETZkAAAAxQTADAACICYIZAABATBDMAAAAYoJgBgAAEBMEMwAAgJggmAEAAMQEwQwAACAmCGYAUIQxZosx5gVjjDXGHDPG7MrYt8kY85y/70zmvoxzj/n7npv50gNIEmb+BwBHxhgrb828e3PsOyZvkePOHPs6JLX7ayICQF7UmAGAu73KWLw4S7ekDn+dvWythDIALghmAOBul6SUXwOWzxdyPJaKpjgAKg1NmQAQgjHmjKQDmc2Zxph2/z+fU1ZzpjGmVV6N2YGZLSmAJKLGDADCeVbSpqzHOqy13fKbOrOaMzsIZQBcEcwAIJznpPEO/YG0vw1GZGY2Z6aiLxKASkFTJgCElNmcGTRj+jVmwejMtLV2Pc2YAMKixgwAwstszgyaMQN7JbX7zZk0YwIIhWAGAOEFzZmb9GEzZiBXcyYAOKEpEwCmwG/OHJB0b1aNWdCcKUlbqTEDEAY1ZgAwNc9KH/Yty7JLUhOhDEBYc8tdAABIqOc0uRkzsFfS6pkrCoBKQVMmAABATNCUCQAAEBMEMwAAgJggmAEAAMQEwQwAACAmCGYAAAAxQTADAACICYIZAABATBDMAAAAYoJgBgAAEBMEMwAAgJggmAEAAMQEwQwAACAm/n+SxcaP8hXDGAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 640x395.55 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlQAAAGJCAYAAABM0K1FAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAyeElEQVR4nO3dXXBb533n8d/Dd0kkAVE2LcWyGlG2m6hK69CyXU+zs5uY7LrTajvpUFJn7N2ZThuq6U6aZKYVK99ofBOW2pkmmU7HpdKrrT07krj17jqdais63bTZdS3baJoySppYkK3IpkNbFADqhW/gsxc4oAHwkDggDnBwgO9nBgPx4CH5Fw8J/PC8HWOtFQAAADavKegCAAAAwo5ABQAAUCYCFQAAQJkIVAAAAGUiUAEAAJSJQAUAAFAmAhUAAECZWoIuoFqMMUbSRyTNBV0LAAAIlS5J79oNNu9smEClTJi6FnQRAAAglHZLeme9BxspUM1J0k9+8hN1d3cHXQuAGvd+al5nX7+mIwd36+7ujqDLARCQVCql++67TyoywtVIgUqS1N3dTaACUNS82tSxrVNd3d3qJlABKIJJ6QDgor21WQfujai9tTnoUgCEQMP1UAGAF5EtrRrcf0/QZQAICXqoAMDFUnpFH9xc0FJ6JehSAIQAgQoAXNy4tai/fOVt3bi1GHQpAEKgZgKVMSZqjBk2xlzw0LZoGwAAgGqpiTlUxph+SQclRSX1FGk7JGmgCmUBAAB4UhOBylobkxRzwtK6jDFRFQlcAAAA1VYzQ34eHZF0NugiADSG5iYTdAkAQqImeqi8MMYMSJosoX27pPacQ12+FwWgbvV2d+j3n3gg6DIAhERoApWkqLU27gz7eXFC0skK1gMAqGMzqXnNzC2sOd7b1a5eds9HgVAEKmPMsLX2dImfNirpT3I+7hIXRwbg0fWbCzr//ff05M/t1I7O9uKfgLrzwqtX9fWXf7zm+BefeEBfHnwwgIpQy2o+UDkrAF8v9fOstQuSVt9aGMNcCADepVesZlILSq/YoEtBQJ56bI8G99+jN2du6ktnvquvHX1I9/d2qreLgI21aj5QKbOqr9+ZQyVJ+yTJGHNcUtxaOxFYZQCAutXb3ZE3tHd/b6cO3BsJsKLaderUqdV/X79+XceOHdPExISOHz8eYFXVVWuBas2WCNbaSeVMRnd6rIattacK2wIA4Kf0itX3riUkSd+7ltDHd3XX1erPRCKhs2fP6ty5c7pwYXN7Zh87dkzHjh1Tf3//6rHDhw/7Ul82qF2+fFmSND4+7ulzotGopMz/r1qhria2TTDG9Dk9TseU6Y0ac9uTyjl2wvn3WE6vFQAAvjo/Na1PjX1Lz7w4JUl65sUpfWrsWzo/NR1wZf6IxWI6e/asEomEZmdnN/11zp49mxemJOkb3/hGueVpZGREx48f1/Hjx1eD1ODg4Iafkw1gw8PDGh4eVn9/v44dO1Z2LV4YaxtjfoAxpltSMplMqru7O+hyANS4+aW0rs7e1p6erepobQ66HFTZ+alpff75mApfIbN9U8893a8nD+yqdlkVMTExodHRUb3xxhub+nxjjC5fvqy+vr6846dOndp071AikdDhw4d17ty51d6mWCymhx9+2PV7ZW3fvl1XrlxZ/ZxsfeVknVQqpUgkIkkRa21qvXa1NuQHADWho7VZD97D9nWNKL1i9exLl9aEKUmyyoSqZ1+6pMH9O6s2/Hfq1Cn19fUpHo+rr69PQ0NDmpyc1MjIiKRMj1A8Hlc8Htf169c1Nja2+rmnT59WX1+fEomE4vG4otGohoeHfautv79fg4ODGh8f18DAhwNH5Q61vf7664rH46u9X9kQlUgkXNvH43ElEom8MJU1OTmZV1sl1MSQHwDUmlsLy3rj7Ru6tbAcdCmosotXZjWdnF/3cStpOjmvi1c2P0xWisOHD6+GqOzwVywW08DAgMbGxlaH7LKPT0xMKBaLScr0PknSwMCAhoaGNDS04RXeNuXcuXOSMsNxxhgNDg5qctLzPtyuotGobty4kTeUmP2a6/VOxePxdb/WeiHMT/RQAYCLWwvL+vsfva/7tm/RtnaeKhvJzNz6YWoz7coRj8c1MTGxGlqkTMAaHx/X+Pi4enp6FI/H83pfsj1Z2TBy7tw5HTlyRNFoVH19fTp48KCvNfb19eny5cuanJzUhQsXNDk5qcHBQZ07d241wHmdx/Twww+v23s2Ojqq8fFx1x6ojfT09JQ1R8wrniUAAMjR2+VtF3Sv7coxOTmpaDSa1+Nz+fLlvN6Ywh6baDS6GiCGhoY0Pj6u7du3q7+/X0ePHq3YqreBgYHVYDcyMqLPfe5zq4HKy+q8jYyMjOjo0aObGqqsRpiSCFQAAOR5dG+PdkU69F5y3nUelZG0M9KhR/eu2enHd4lEQn19fXk9UKXOBbpw4YJisZgmJydXg41foSqRSGhycnLNUOLY2JhOnTq17pymUkxMTGjfvn1Fw9R6Q4HZn2GlEagAAMjR3GR08tB+ff75mIyUF6qyU9BPHtpflQnp/f39Gh0dXXPca1A5ffr06vYB/f39Gh4e1hNPPOFrL9Vrr73mOjerr69vtcbNDvlle+ayx7LzxdwCUvb7ZSfu56r0hHSJQAUArtpamtR39za1tbB2pxE9eWCXnnu6X8++dClvgvrOSIdOHtpftS0TBgYGdPDgQU1MTOSFlrNnz67bY5M7ATuRSKyGqqzCsLHekFg8Htfk5GTRnqHTp09rcHAwL7QU9lptZsgvFospFotpaGhodYhzYmJitZ7s/LLccHjixIm8mnPbVxr7UAEAsI70itWZ167qmRen9JXPHtDRR/YEslP6yMiI9u3bp56ezDDj0NCQYrGYRkdHNTExobGxMR0/flynTp3S6Oio+vr6dOLEidWwlP28eDyu4eHh1Z6ciYkJnTlzRrFYTMePH9cjjzyyGoROnz6tkZGRNfs65crutN7X17e6sjCrnF6wRCKhvXv3uq7Oy+aWiYkJjYyMrO6inpXdYkLK9J7lbiGxGV73oSJQAYCL9IrVwnJa7S3NdXWpEZRu6p2kfu1Pv6NvfuFTDXctv2zPUDXmINUqr4GKvmwAcHH95oLGvx3X9ZsLQZcCBCYWizV0mCoFc6gAAHAxk5rXzNyC3py5KUmr971d7ertrvyWCbWgWlsO1AMCFQAALl549aq+/vKPVz/+0pnvSpK++MQD+vLggwFVVT3xeNz3TUDrGYEKAAAXTz22R4P771lzvLerPYBqqo+hvtIQqAAAcNHb3dEwQ3soH4EKAFzc1dmu3/v0PrU2sXYHQHEEKgBw0dRk1N7UHHQZAEKCt14A4OLGrUX9VeyabtxaDLoUACFAoAIAF0vpFb19/baW0itBlwIgBAhUAAAAZSJQAQAAlIlABQAAUCYCFQC46Oxo0ac/1qvODhZDAyiOZwoAcLG1rUUP3RcNugwAIUEPFQC4mF9K6wfTKc0vpYMuBUAIEKgAwEXqzpLOT72n1J2loEsBEAIEKgAAgDIRqAAAAMpEoAIAACgTgQoAXLQ0N2lXpEMtzTxNAiiObRMAwEXPtjb95qN7gi4DQEjw1gsAAKBMBCoAcDGTmtdXL/xIM6n5oEsBEAIEKgAAgDLVzBwqY0xU0hFJh621gy6PH3f+uU+SrLXHqlcdAADA+moiUBlj+iUdlBSV1OPy+Ji1diTn43FjzAW34AUAAFBtNTHkZ62NWWtPS4oXPub0XPU791njkgaMMX3VqRAAAGB9NRGoPDgoKTc8ZYNXtPqlAGgEPdva9Fu/9FH1bGsLuhQAIVATQ34bsdYmJG0vODzg3K/p0coyxrRLas851OVvZQDqWUtzk6JbCVMAvAlLD1WhE5KOOWFrozbJnNu1KtQFoE4kby/p/NS0kreXgi4FQAiELlAZY8YknXHmXG1kVFIk57a70rUBqB8Ly2n9YHpOC8vpoEsBEAI1P+SXyxgzJOmyhzAla+2CpIWcz61kaQAAoIGFpofKGDMgSdkwZYyJssoPAADUgloLVGv2oJJW96nqlxQzxvQ5QWpY0mw1iwMAAHBTE0N+TkAaknRUmT2nxiS9Zq2dcPafelmZLRLGcj/PWnuqyqUCaBBb21v0i307tLW9Jp4mAdQ4Y60NuoaqMMZ0S0omk0l1d3cHXQ4AAAiBVCqlSCQiSRFrbWq9drU25AcANWFhOa23PrjFKj8AnhCoAMBF8vaSXvynd9iHCoAnBCoAAIAyEagAAADKRKACAAAoE4EKAFw0NRlFt7aqqYmrLAAojg1WAMDFXZ3t+q1f2ht0GQBCgh4qAACAMhGoAMDF+3ML+vNvX9b7cwvFGwNoeAQqAHBhrdWdxbQa5WoSAMpDoAIAACgTgQoAAKBMBCoAAIAyEagAwEV0a5uOPnKfolvbgi4FQAiwDxUAuGhradJHoluCLgNASNBDBQAu5uaX9O0fva+5+aWgSwEQAgQqAHBxZzGt2Ns3dGcxHXQpAEKAQAUAAFAmAhUAAECZCFQAAABlIlABgIuOtmb9wn0RdbQ1B10KgBBg2wQAcNHd0arPfOyeoMsAEBL0UAGAi6X0imZS81pKrwRdCoAQIFABgIsbtxb1wqtXdePWYtClAAgBAhUAAECZCFQAAABlIlABAACUiUAFAG5M5gLJMkEXAiAMjLU26BqqwhjTLSmZTCbV3d0ddDkAACAEUqmUIpGIJEWstan12tFDBQAAUCYCFQC4uH5zQf/1lbd0/eZC0KUACAECFQC4SK9YXb+5qPRKY0yLAFAeAhUAAECZauZafsaYqKQjkg5bawddHj8uKeF8GLXWnqpedQAAAOuriUBljOmXdFBSVFKPy+PHJclae9r5eMAYM26tPVbNOgEAANzU1LYJxpghSSestQ8XHL8haa+1NpFzzFprPe8Qw7YJAEoxv5TWO4k7uje6RR2tzUGXAyAgdbNtgjGmT5khvoTLYwPVrwhAI+hobda+uzsJUwA8qflAJalvneMJZYYIXRlj2o0x3dmbpK4K1AagTt1aWNbFK7O6tbAcdCkAQiAMgWo9s3KZb5XjhKRkzu1aNYoCUB9uLSzr/775AYEKgCdhDlQbhSlJGpUUybntrnhFAACgIdXEKr8i4uscj27wmKy1C5JWtzg2hiucAgCAyqj5HiprbVxSwpmcXvjYZAAlAQAA5Km1QLXeMN6opNUVfc72CqerUhGAhtTe0qwH7ulUewur/AAUVxP7UDm9T0OSjkrql3RK0mvW2omcNsf14RDfI9bakRK/B/tQAQCAknjdh6omAlU1EKgAlCK9YnV7cVlb21rU3MQcTKBR1c3GngAQhOs3F/QX/3BF128uFG8MoOERqAAAAMpEoAIAACgTgQoAAKBMBCoAAIAyhWGndACouru72vWFz9zPCj8AnhCoAMCFMUYtzYQpAN4w5AcALm7cWtS513+iG7cWgy4FQAgQqADAxVJ6Rddu3NFSeiXoUgCEAIEKAACgTAQqAACAMhGoAAAAykSgAgAXXR2tGtx/j7o6WoMuBUAIsG0CALjY0tasA/dGgi4DQEjQQwUALu4spjX1TlJ3FtNBlwIgBAhUAOBibn5JFy79VHPzS0GXAiAECFQAAABlIlABAACUiUAFAABQJgIVALhobW7S7u1b1NrM0ySA4tg2AQBcbN/WpsMH7wu6DAAhwVsvAHBhrdVyekXW2qBLARACBCoAcPH+3IL+9Ftv6v25haBLARACBCoAAIAyEagAAADKRKACAAAoE4EKAACgTGybAAAudnS263f+zV5tbeNpEkBxPFMAgIvmJqOujtagywAQEgz5AYCL5O0lffN77yp5eynoUgCEAIEKAFwsLKf145/e1MJyOuhSAIQAgQoAAKBMBCoAAIAyhWpSujFmWFJUUkLSPkmj1tpEgCUBAACEJ1AZY45LOp0NUMaYqKRvSDocYFkA6tS29hb90v13aVt7aJ4mAQQoTEN+g7m9Uc6/o0EVA6C+bWtv0aN7ewhUADwJU6BKGGMuOD1TMsb0SYoHWxKAejW/lNbl929qfolVfgCKC1Og+pykPkk3jDFjkgastcfWa2yMaTfGdGdvkrqqVSiA8EvdWdL/+u67St1hHyoAxYUmUDlDfGOSJiQdl3Q421u1jhOSkjm3axUuEQAANKjQBCqnVypurT2szAq/HklvbPApo5IiObfdFS8SAAA0pFAEKme+VNRaOylJ1tq4tfZhZeZVDbl9jrV2wVqbyt4kzVWxZAAA0EBCEaiUmTuVcDk+XuU6ADSI5iajHZ1tam4yQZcCIARCEaicnql+lzlTD1trJwIoCUCd29HZrv/0+Ee1o7M96FIAhECYNlg5LOmEMea6Mr1VUUkjQRYEAAAgVSBQGWM+aq19y/n3J5UZrnsje2yznFV+BCgAVTEzN69zr1/T4YO71dvVEXQ5AGpcJYb8BrL/sNb+k7X2v+ceA4BQsNLi8opkgy4EQBj40kNljIlIOqLMU8+gMXmTOKOSHpH0F358LwAAgFrjS6Cy1iaNMZPKDMntk3Qj5+GEpD/y4/sAAADUIt/mUFlrr0j6XWPME9bal3MfM8Z81K/vAwAAUGt8n5RurX3ZGPMZZYb6so46NwAIhe3b2vTUY3u0fVtb0KUACIFKrPI7q0yYSuQc/qTf3wcAKqm1uUm93azuA+BNJfahumCt/UbuAWPMExX4PgBQMan5Jb3+1qwOfrRH3R2tQZcDoMZVYtuE6x6PAUDNml9M659/ktT8YjroUgCEQCV6qPYZY/63pFjOsQFltk4AAACoO74GKmc/qqOSzhQ+5Of3AQAAqCW+BipnP6rPWWv/Kfe4s0cVAABAXfJ9DlVhmHLccDkGADVrS1uz+n9mu7a0NQddCoAQ8OvSM78hadJamzLG/IFLk0FJ/96P7wUA1dDV0ap/++DdQZcBICQ23UNljPmdnA+fkXTQ+fdvKjNnKve2Y7PfBwCCsLi8oncTdzIXSAaAIoy1m7uUujEmba1d0xdujPmkyxyqNceqzRjTLSmZTCbV3d0dZCkAQmAmNa8XXr2qpx7bwwafQANLpVKKRCKSFLHWptZrV84cKteVe27BKegwBQAAUEnlBKrNdW0BAADUmbJ6qIwxf2CMecgZTgMAAGhI5azys5ImJD0s6RljzF5Js8rskP6anFV/5ZcIANVnjNGWtmYZw77EqA8zqXnNzC2sOd7b1Z43T9BrO+QrZ1L6iqRobmgyxnxSmZ3SB5SZvPWAL1X6gEnpAICwqET4+eqFH+nrL/94TdsvPvGAvjz4YMntGiV4eZ2UXk6g+ltJo9bav9tcidVFoAIAhIXf4Uf6MAC9OXNTXzrzXX3t6EO6v7dz3ZBWrF2jBC+vgaqcIb/Dkr5hjLlirX2rjK8DADXng5sLeumf39WhX/iI7upsD7oc1Amv4eKpx/ZocP89rqEml9d2ktTb3ZH3Pe7v7dSBeyObbuf1e7/w6lXPoS8oG50Xr5Fv04HKWpuUdMQY84Sktzb7dQCgFq2sWCVuL2llhQXNKM5rUPIaLvwOP5Xgd/AKsidro/Py24/t9PQ1yr70jLX25XK/BgAAYeY1KJXSo1QvvAYvrz/DUoKXPz2Ci57+n75cyw8AgHrk9xBdkD1Kta4SQ4h+9AimUgQqAADK4vcQHdbn9xBiqW3LRaACABeRra367CfvVWRra9ClIECNOERX60oJr9UMugQqAHDR3tKsj961LegyUCFeh/LoeYJXBCoAcHFzYVn/ci2pT+yOqLOdp8p6E4al/AgXniUAwMXthWX9Y/y69t29jUBVhxjKg994lgAANByG8uC3pqALAAAACDt6qAAAdSPs141DeIUuUBljxiRddj6ctdZOBFkPgPrU3tKsj+/qUntLc9CloARMNkdQQhOojDFRSS9LesJamzDG9Et6Q5IJtDAAdSmytVVPHtgVdBkoEZPNEZTQBCpJY5LOWGsTkmStjRljBoMtCUC9Wk6v6ObCsjrbW9TSzHTTsGCyOYISpkA1LGmfMaZPUp+1dtJaOxl0UQDq0+ytRb3w6lU99dge5t4EjHlRCINQBConRElSv6S4pLgxZlzSufVClTGmXVJuH29XZasEAFQC86IQBqEIVJKygSphrY1JkjFmRNIVSdvX+ZwTkk5WoTYAQAUxLwphELaJAa9n/+HMpYoaYwbWaTsqKZJz213x6gAAvuvt7tCBeyO6v7dT0ofzohjuQy0JSw9VfJ3jCX3Ye5XHWrsgaXXQ3RgWAwKNzOs8nGy72VuLmpmb1w/fm9PM3MK67bx+vWLtGhE/G9STUAQqa23cGBNXJjzFch6KKqfXCkB9KOWF1mtbr/NwCtv9t4s/8dTO69dbr10jhgvmRqGehCJQOUYkHZUTqIwxQ5Ims3OqAHjnd++K3+1KeaH12tbrPJyg2jViuGBuFOpJaAKVtXbCGNNjjDnuHNphrWUfKtS9WuqtqVa7Ul5ovbb1uj9Rtl3yzpIk6e6u9g3bef16xdp5/X8E2ZPld3BmzyjUk9AEKkmy1p4OugbAL/TWlB9+Sm1bivSKzbuvNK//j1J+H4LqOWzE3jYgVIEKqCdh6q2pdjusr5Tfh6B6DhnKQyMiUAE+8/puP0y9NagdpZzjoHoO+T1EIyJQAT7z+m6fFx1UGgEIqB4CFeAzhjvqw7b25rx7ANgIgQrwiJVLjaW9pTnvHgA2QqACPGLlUmOZX0rn3QPARghUgEcM5TWW24vpvHsA2AiBCg2PoTwAQLkIVGh4DOUBAMpFoELDYygPAFAuAhUaHkN5cNPSZPLuAWAjBCrUrSAvIovw697SmncPABshUKFuMTcK5VixNu8eADZCoELdYm4UypG4vZR3D/glvWL1vWsJSdL3riX08V3danYZWva7XaltUZqmoAsAKqW3u0MH7o3o/t5OSR/OjWK4D0AlFIaV9Mra3s3zU9P61Ni39MyLU5KkZ16c0qfGvqXzU9MVbVdqWy//F+QjUCF0ZlLzmnonueY2k5oPujQAdcivoHR+alqffz6m6WT+c9V7yXl9/vnYalu/222mLcGrdAQqhM4Lr17Vr/3pd9bcXnj1atClAagzfgWl9IrVsy9dklvcyB579qVLWlxe8bVdesV6/t7pFUvwKgOBCqHz1GN79M0vfEpfO/qQJOlrRx/SN7/wKT312J5gCwMQKsVe5P0MSv8Yv77m6xS2nU7O6y9fecvXdhevzOrilVlPbf/x8vWGDl7lfm8CFUKnnuZGef0DLuUPvRHeCVZDxNkuIcK2CXWp2Iu830HplcvXPdX19uxtX9vNzM1rZs7bdIhX4h/UZfDyc27bRghUqAn1Ni8qiMmppbatRJirJ9mVT6yAqj9eXuS99up4DUpyjR9r/UzPVl/b9XZ1qLfL65tNb7/rYQpefgzZXrj03gY/jQ8RqFATwjIvqlYnp26mrd9hzu+AFnSQm5tfyrtHeGz0u+O15+k9z2/mvP1ePt53l3ZFOtaNLEbSrkiH/uPjH/W13aN7e/To3h5PbR/ft8PT/yVMwcuPIds//psfFv2/SgQq1IgwzIuq1cmplZp0GmRAK7X7vTJhLum0SzZMr1w9KPa747Xnafbm2qssuPEalH5x3w6dPLR/9VhhG0k6eWi/2lqafG3X3GTU3GQ8tf3Fvh11E7z8HLL9acrb7wKBCjUh6HlRYZ6cWolJp6WEuSCXgmfbVyLMjTrvSkf/5ocMm4aEl98dr/OJera1+RqUmpuMnjywS8893a+dkfzntZ2RDj33dL+ePLBLknxv57VtvQSvUp4TvQ/ZFkegQsML++TUSkw69Rrm/A5opQQ5Kfh9fSrR2wZ3xUKp17/Tuzq9XSlhZ2SL70FJygSb74x8Rl/57AFJ0lc+e0DfGflMXptKtPPath6CVynPiV6HbL0gUKGian2yeT1MTq3EpFOvYc7vgFZKr1yQ+/pUqreNnix3XkKp179TWXmee1SJoCRlFjr8/O6oJOnnd0fXXfjgdzuvbcMevEp5TvQyZHtPt7cQTqBCRQU92bwRJqdWYtKp1zDnd0ArpVfO6wuo32GuEsOmUnj26qk2r6HUa4/EB7cWPPc8SZUJSmEQ5uBVynOilyHbP/qVj63zVfIRqFBRQU42b5TJqZWYdOo1zPkd0ErplfP6Aup3mPO7V+7ildmKLRkPAz/e9KRXrOceid6ujpJ6nqT6Ckp+q8XgVcpzopch28H9Oz39LAhUqKhKTTb3YxJ5PU1O9drW65OM1zDnd0ArpVfO6wuo32HO716595J3KrJkPCuo7Sz82GLE65uei1dmPfdIPLq3R1JpPU8oXzWD12bblvv7QKBC6Pg1ibzeJqd6betnmPM7oJXSK+f1BdTvMOd3r9zsrUXfl4xnBbWdhV9bjHh90zMzN19Sj0QWPU+1x6/gtZm25f4+EKiwKUFNNvdzEnk9Tk712tbPMBfUUvCgwpzfvXI9HoN9KUvGpeBWQPq5xYjXNz3Z3spSh/IQXn4/J/qBQIVNqcRkc7+WRXudRN7ok1P9DHNBLQUPIsz5HeR2eh7+9r5kPKgVkF7bed1ipJQ3PVkM5SEoBCpsit+Tzf1cFu11EjmTU/0V1FLwSoW5E87KnhO/8rGK9sr5vUKzt6sjsBWQXtt53WKk1Dc9WfydIggEKmyKn5PN/V4W7XUSOZNT60clwtyBeyOSpAP3RiraK1eJJeNBrYD02s7rFiObedMDBCW0gcoYcyHoGuBNtZdFlzKJPIt3tCgU2dKad+/GryDn95LxoFZAem3ndYsR3vQgTEIZqIwxQ5IGgq4DxQW1LJp3tQgbPxcABLUC0mu7UrYYyeJND2pd6AKVMSYqqadYO2yOn6v3gl4WzbtalGP21mLefTX4tQAgqBWQpWx7wZse1JvQBSpJRySdDbqIeuXX6r1aWRbNu1rUI6/BK4jtLKq9mSJQK1qCLqAUxpgBSZMe27ZLyn217qpIUXXmqcf2aHD/PXpz5qa+dOa7+trRh3R/b6d6u9YGn8K5UR/f1b36xF7qXlDvJeddw5dR5om4cFn04P6dOvPaVT3z4pS+8tkDOvrIHsISUMDr30pQ7STe9KB+hK2HKmqtjXtse0JSMud2rWJV1RGvq/eKzY2q1IVKs3gSBrwJajsL/kbRaEITqIwxw9baiRI+ZVRSJOe2uyKFNSAvc6MqeaFSAABqTSiG/Iwx/ZJeL+VzrLULklZ3eDSGd0d+KDY3yigzN+rbf/jpkobyGMZDrenuaMm7B4CNhKWHqkfSgDHmuDHmuKQxSXI+Hgq2tPqz0b5RXudGvfH2DZZFI9Rampvy7gFgI6F462WtnVTOZHSnx2rYWnsquKrCZSY1r5m5tZdk6e1qz5sfdX5qWs++dGk1ND3z4pT+9Ftv6uSh/XrywK6Stjn49Yfu1XNP9+d9PSnTM5X9ekCturmwnHcPABsJRaDK5fRIHXX+PSbpghO4sIEXXr2qr7/84zXHv/jEA/ry4IOSPpwbVThEl50b9dzT/SXNjZIYykN4LS6v5N0DwEZCF6iciemlTE6Him+HUKm5UVLjDeVlewPfnLkpSav3hb2BXtsBAGpf6AIVNqe3uyPvRTq7HUJWqXOjPv98TEb5lzjdaJuDeuA1ABX2Bn7pzHcl5fcGltKO4AUAtY9A1UA22oizHudG+d1T5DUAZXsDCxVujuq1ndfvCwAIDoGqQRSbbB6muVFB9RR5DUCFvYHr8drO6/eV6M3yU0drU949AGyEQNUAvEw2H9y/MzRzo4LqKfIagPxWyvdlGNE/W9ta8u4BYCM8U9Q5r5PNB/fvDHRuVCkv8EH1FIWB38OIjRy8ltIrefcAsBECVcgV21/K62Tzi1dmVy8BE8TcqFLmCdVTAPKb38OIjRy85uaX8+4BYCMEqpArtr9UKZPNJf/nRnl9oS1lnhDKF1TwAoB6RaAKuace26PPfKxXf3vpPf3Z313Wf/70Pv2yMx9KUsmTzSV/50Z5faGl16k2+R286rEnCwAkAlXoxa7eyBui+7O/u6y/ir2zOkT36N6ekiebe0HPE3J5DV6VGEIkpAGoBQSqEPOyeu/JA7sqMtmcnidsRiWGECs1z8uY/HsA2AiBKqRKWb1XymRzep5QSZXYe6tSE+yv31yUJF2/uaipd5L0eAHYEIEqpEpZvff4vh2eJ5vT84RaUMrvFysbAdQCAlVIlbp6T/I22ZyeJ9SrUoPXjduL+j//+r7+3c/ere1b21jZCGBDBKoaVGxvqcy/S1+956k9PU9ocNm/gZnUvF5/64Z+9p4u178JVjYCyEWgqkHF9paSVPLqPZ7UAX8FubIR1eP3RdY5z/WLQFWDiu0tJWWG70pZvcfwBBAMNketTUFdZN3vgE1Aqx0EqhpUbG+prFJW7zE3CggGm6NWl99Bye+LrPsdsEsJ4vzuVBaBqsZ43Vsqy+vqPeZGAaWJbG3V0MO7FdnaWpXv5/cQotSYL6B+ByW/L7Lud8Au5c0yw8+VRaCqIaXsLZUNTNlf/K1tmVO5ta1FP5hO8YsPlKm9pVn39WwNuow1GvEFtJT6/A5KQfE7oElsHVJpBKoaUureUhLzLoBKmZtf0j//JKlfuC+iro7q9FJ5EYYXUL/blfI8V+tBKUjs2VZZBKoaspm9pZgbBVTGncW0XntrVg/e01lTgaoUQb2A+t2O57nqInhtDoGqhmxmbynejQEoV1DzeupleK5REbzyEahqSKl7SwFANQU1QRvhFlTwkvwZfvb6G2qsdXvprj/GmG5JyWQyqe7u7qDLWVd2lZ/kvrdU4So/AJUxk5rXC69e1VOP7eFFH6ghXq4mkvXVCz8qulF2sXa//dhORSIRSYpYa1Pr1UWgqqKZ1Lymk/P6/rtJ3bi9pO1bW/VzH4loVyQ/vZ+5eFX/5W//VR84V7uXpLs62/SHv/yzOvroniBKBxpO8s6SLl6Z1aN7exTZEs45VECj8xq+NmrXoUVPgYohvyp69qVL+ut/mV5z/Fc/sUt/9lT/6sfvJufzwpQkfXBzUe9usAIQgL8iW1pdhx8AhIcfw8+p1KLr8UIEqio5PzXtGqYk6a//ZVqHpqZXh/JY0QIEbym9ouSdJUW2tKq1uSnocgDUOAJVFWQ37FxP4YadTNQEgnfj1iJzqAB4xtuuKihlw04AABA+BKoq2MyGnQAAIDwIVFWwmQ07AQBAeBCoqiC7YedGdrFhJ1BzshchB4BiCFRV0NxkdPLQ/g3bnDy0nydvoIb0dnfo9594gAnpADwJVaAyxhx3buPGmPGg6ynFkwd2aew3PqG7Otvyjt/V2aax3/gEu58DABBiodkp3RgzZq0dyfl4XFKftXbQ4+cHvlO6lNlC4eKVWc3Mzau3KzPMR88UUHuu31zQ+e+/pyd/bqd2dLIHHNCoUqlU/eyUboyJSuo3xkSttQnn8LikN4wxfdbaeGDFObwGpeYmo8f37QigQgClSK9YzaQWlF4Jx5tOAMEKRaByHJTUJynmfJwNUdFAqslxfmpaz750KW+vqV2RDp08tJ+hPAAAGkAo5lBZaxPW2u3W2ljO4QHn3rV3yhjTbozpzt4kdVWitvNT0/rd52NrNu6cTs7rd5+P6fyU++VmAABA/QhFoFrHCUnHcoYA3R5P5tyu+V1AsUvKSJlLyjBkAABAfQtloDLGjEk6Y609vUGzUUmRnNtuv+sodkkZiUvKAGHVvaVVv/rzu9S9pTXoUgCEQJjmUEmSjDFDki4XCVOy1i5IWsj5PN9r4ZIyQP3qaG3Wg/dUZKYAgDoUqh4qY8yAJGXDlDEmaozpC6oeLikD1K9bC8t64+0burWwHHQpAEIgNIHKGNMvqV9SzBjT5wSpYUmBjadxSRmgft1aWNbf/+h9AhUAT0IRqJx9qF6WNCbpcs5tbINJ6RXHJWUAAIAUkkCVs22CKbwFXduTB3bpz5/uX9NTtSvSoT9/up99qAAAaAChm5Rei548sEuD+3dySRkAABoUgconXFIGqC9tLU3qu3ub2lpC0ZEPIGAEKgBwEd3apl9/6N6gywAQErz1AgAX6RWr24vLXOkAgCcEqiLSK1avXL6u//ndd/TK5es8uQIN4vrNBY1/O67rNxeKNwbQ8Bjy28D5qWk9+9KlvMvL7Ip06OSh/azeAwAAq+ihWsf5qWl9/vnYmmv1vZec1+efj+n81HRAlQEAgFpDoHKRXrF69qVLchvcs87t2ZcuMfwHAAAkEahcXbwyu6ZnqtB0cl4XrwR21RsAAFBDmEPlYmZu4zBVajsA4XNXZ7t+79P71NrE+04AxRGoXPR2bXzB41LbAQifpiaj9qbmoMsAEBK89XLx6N6eNdfmK7Qrkrm8DID6dOPWov4qdk03bi0GXQqAECBQuWhuMjp5aL+MpMKr8WWPnTy0n2v1AXVsKb2it6/f1lJ6JehSAIQAgWodTx7Ypeee7tfOgp6qnZEOPfd0P/tQAQCAVcyh2sCTB3ZpcP9OXbwyq5m5efV2ZYb56JkCAAC5CFRFNDcZPb5vR9BlAACAGsaQHwC46Oxo0ac/1qvODt53AiiOZwoAcLG1rUUP3RcNugwAIUEPFQC4mF9K6wfTKc0vpYMuBUAIEKgAwEXqzpLOT72n1J2loEsBEAIEKgAAgDIRqAAAAMrUsJPS0yuW/aUAAIAvGjJQnZ+a1rMvXdJ0cn712K5Ih04e2s8O6AAkSS3NTdoV6VBLMx35AIoz1tqga6gKY0y3pOTEK/+qP/gfP1bh/zrbN8VlZQAAQFYqlVIkEpGkiLU2tV67hnvr9cd/88M1YUrS6rFnX7qk9EpjhEwAAOCPhgtUP00trPuYlTSdnNfFK7PVKwhATZpJzeurF36kmdR88cYAGl7DBSovZuZ4AgUAAN4RqFz0dnUEXQIAAAiRhgtU93S3a73NEYwyq/0e3dtTzZIAAEDINVyg+qNf+ZgkrQlV2Y9PHtrPflQAAKAkDbdtQjKZ1P+7eot9qABsaDm9opsLy+psb2EvKqCBed02oSEDVXd3NzulAwCAorwGqlDtlG6MOS4p4XwYtdae2uzXam4yenzfDl/qAlB/kreX9Er8Az3ed5ciW1uDLgdAjQtNP7YTpmStPW2tPS0pZowZD7gsAHVqYTmtH0zPaWE5HXQpAEIgNIFK0glJp7MfWGsnJQ0HVw4AAEBGKAKVMaZPmSG+hMtjA9WvCAAA4ENhmUPVt87xhKSo2wPGmHZJ7TmHuqTM5DIAKGYuNa9/+P5V/YePR9WhxaDLARAQr7khLIFqPbOS1tuF84Skk4UH77vvvooWBKC+vPjloCsAUCO6JNXHKj8XG21pPirpT1za5175uEvSNUm7Jc35Wxo2iXNSmzgvtYdzUps4L7Wp3PPSJendjRqEJVDF1zkeXe8xa+2CpIWCw3nJ0pjVfafmNtpbAtXDOalNnJfawzmpTZyX2uTDeSn6OaGYlG6tjUtKOJPTCx+bDKAkAACAVaEIVI5RSasr+owxQ8rZRgEAACAooQlUzq7oUWPMkBOmHrHWHivzyy5IelZrhwYRHM5JbeK81B7OSW3ivNSmip+XhrmWHwAAQKWEpocKAACgVhGoAAAAykSgAgAAKBOBCgAAoExh2djTd8aY48pcC1DKXHj5VIDlNBxjTFTSEUmHrbWDLo9zfgLi/OwlaZ8kFa6m5dxUV87fipQ5J32SPpd7sXjOSfCMMRcKn8s4L9VljBmQdEzSBWU2/R6U9Jq1diKnTcXOSUP2UGVfMKy1p621pyXFjDHjAZfVMIwx/cq8QETlcvkgzk9wjDFj1tpTzu2Yc+xCzuOcm+obkzTp/MxHlLl81rnsg5yT4Dlb+QwUHOO8VF9UmfMw7twuu4Spip2Thtw2wRhzQ9Legnd41lpr1v8s+M15EjphrX244DjnJwBOT8g5ZXoNE86xfklvSNpnrY1zbqrPCbQXsu+knReFE9ba7c7HnJMA5fQgjuf+zDkv1ee8pkzm/swLHq/oOWm4Hirn8jVRtx+4012IAHF+AndQmSGlrOy1MqOcm2BYawcLhiUekTQp8fdSI45IOpt7gPNSe6pxThpxDtWa6wE6Esp0FyJYnJ+AOE802wsOZ59o4sqELTcJcW6qwnkHHpV02DnE30uAnBdit+vJcl6Cc8QYM6vMdJJ9zjC5VIVz0oiBaj3ZE4DaxPkJxglJx6y1iZyrtRfi3FRYzrBSVNK59YY0cnBOqiPqDIVHPbbnvFRWTJKstXFJMsYMG2POWWsPb/A5vp0TAtWH+CWvbZyfKjPGjEk640ze3AjnpsKcAHVaWn2RuCFp7wafwjmpMGPMsIe/jUKclwrKBqkcZyWNFwm8vp2ThptDpQ/nhBSKbvAYqofzUwOcoaXLBXN3ODdVZoyJGmPGCl4QJvXhaibOSQCcxRqvb9CE8xIA53lrVU5Pbp+qcE4aLlA5CTbhTFArfMxtLBxVxPkJXnaCZvbdt/Oi3se5CUSfpOPKfxcdde4TnJPA9EgaMMYcd1ZdjkmZFZjGmCHOS/VlVynn/sxz3ojEq3FOGi5QOUaVs2eIk2pL7bpF+dbrauX8BMR5592vzP4sfc6Tz7Ay8wwkzk1VWWtjkk4VDGUclRTLeRHgnFSZtXYyZ7+2U8rseSTn4+y+R5yXKnJ6owr/VoYlTeT0VFX0nDTkPlTS6l4u2R/8IzkrAVBhzov0kDIvDP2STsl9N1vOTxU57+auyGXFS8H+OpybKnLOy3DOoX2SRlx2SuecBMB5UT6qzHPaKWX2DMtua8F5qSKXv5UdhT/zSp6Thg1UAAAAfmnUIT8AAADfEKgAAADKRKACAAAoE4EKAACgTAQqAACAMhGoAAAAykSgAgAAKBOBCgAAoEwEKgB1z7nG2mVjjDXGnMter9B5bNgY84bz2AVjzHDB5445j10ufAwAstgpHUBDMMaMSxrOvYxOzmPZC9xuz72kS87j56y1hytfJYCwoocKQKNIeGiz5oLdTm8W12ADsCECFYBGcVlavYDqKufjR5wP8x5z9BVcwR4A1iBQAWgUs859YS/UEUmjbo8ZY4attacrXRiA8CNQAWgUCec+mj1gjOmX9Po6j/VJomcKgCcEKgCNIttD1Zdz7KC1Nib33qsha+1kVSoDEHoEKgCNIuHc90iSMWZI0llJylnZt895bEDSRHXLAxBmBCoADSFnYnk0OzG9YIuEhD4c8mMiOoCSEKgANJodko5Yawt7oGYl9TARHcBmEKgANJKEpAFlJqK7PdYvJqID2ISWoAsAgCqalfS6MxHd7bFZJqID2Ax6qAA0kpjW3/U8JulYFWsBUEe4lh8AAECZ6KECAAAoE4EKAACgTAQqAACAMhGoAAAAykSgAgAAKBOBCgAAoEwEKgAAgDIRqAAAAMpEoAIAACgTgQoAAKBMBCoAAIAyEagAAADK9P8Bt/EXMB94LdsAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 640x395.55 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"c_obs3.plot_rho()\n",
"c_obs3.plot_tauint()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now redo the error analysis and alter the value of S or attach a tail to the autocorrelation function via the parameter `tau_exp` to take into account long range autocorrelations"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Result\t 3.27194697e-01 +/- 1.67779862e+00 +/- 2.08884244e-01 (512.783%)\n",
" t_int\t 5.69571763e+00 +/- 2.09295390e+00 tau_exp = 20.00, N_sigma = 1\n",
"1000 samples in 1 ensemble:\n",
" · Ensemble 'ens1' : 1000 configurations (from 1 to 1000)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlQAAAGJCAYAAABM0K1FAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAA9hAAAPYQGoP6dpAABEZUlEQVR4nO3deXxcV33///eZ0b7NWLZleY0tOZswWWQnISwtSaw2FBwKeGm/CfRBCzYppUDbr01o+zUpFD/klq2Upnb4lv5o0jaxgW8IoS5WCFsIsR0lpIqz2XIWO3JkWxqt1jZzfn/MkpE00lx5rmbRvJ6Phx6y7r2aOb6S5r7nnM8511hrBQAAgAvnyXQDAAAAch2BCgAAIEUEKgAAgBQRqAAAAFJEoAIAAEgRgQoAACBFBCoAAIAUFWS6AelijDGSlkjqy3RbAABATqmU9JqdZvHOvAlUCoepk5luBAAAyEnLJJ2aamc+Bao+SXr11VdVVVWV6bYgh5zpHdL9R05q87plWlhVkunmAADSqLe3V8uXL5eSjHDlU6CSJFVVVRGoMCNDKlJJeYUqq6pURaACACRAUTqQRHGhV2uW+lRc6M10UwAAWSrveqiAmfKVFqqpYVGmmwEAyGL0UAFJjAZDOts/rNFgKNNNAQBkKQIVkET3wIj+7bGX1T0wkummAACyVNYEKmOM3xiz1Rhz0MGxSY8BAABIl6yooTLGNEpaJ8kvqTrJsRslrU9DswAAABzJikBlrW2V1BoJS1MyxviVJHABAACkW9YM+Tm0WdL9mW4E8o/XYzLdBABAFsuKHionjDHrJbXM4PhiScVxmypdbxTyQk1Vif70posz3QwAQBbLmUAlyW+tbY8M+zlxh6Sds9geAECO6OwdUmff8KTtNZXFquEOCHBBTgQqY8xWa+3eGX7bLklfjvu6UtwcGRfgXP+wDjxzWje/qVbzK4qTfwOArHPv46/oaw+/OGn7J2+6WJ9uuiQDLcJck/WBKjID8MhMv89aOywp9nbEGGpgcGGCIavO3mEFQzbTTQFwgW69boWaGhbpWGe/PnXfU/rqlqu0uqZCNZW8SYI7sj5QKTyrrzFSQyVJ9ZJkjNkuqd1auz9jLQMA5ISaqpJxQ3urayq0Zqkvgy1yXyAQ0P333699+/bp4MG5uVzj7t27JUnHjx+XJO3Zs2fSfr/fLyl8PrZv3562tmVboJq0JIK1tkVxxeiRHqut1trd6WwYACC3BUNWT58MSJKePhnQ5Yur5swM3tbWVh05ckSBQEBdXV2Zbs6s2LFjh5qbm2Nfb9u2TU1NTbHwGA1bW7dulSS1tLRo27Ztk0LXbDHWZn4YwxhTJ2mjpC2SGiXtlnR4Yu9TZJ2qLZFjd0s6GAlcTp6jSlJPT0+Pqqqq3Gw+5rjO3iHd+/gruvW6FRSvAjnqQFuH7nzwqDp6hmLbFvtKtHNDg25esziDLXPX/v37tWvXLj3xxBOZboqrAoGANm3apH379sV6oFpbW7V27VodP35cdXV1mjdvnk6cOBHbL4XLfVLNOb29vfL5fJLks9b2TnVcVvRQWWvbFQ5I0/Y6RQIWQ3xIq6rSQr37isWqKi3MdFMAXIADbR26/Z5WTbysnu4Z0u33tOqu2xrTGqp2796turo6tbe3q66uThs3blRLS4t27NghSbr77rvV3t6u9vZ2nTt3blyvzN69e1VXV6dAIKD29nb5/f5Yj4wbor06XV1dCgQC48LJxLCSbkeOHFF7e7saGxslSXV1dZIUOxcT2xvV0tKi9etn/wYrWRGogGxWUujVJYtYxgzIRcGQ1Z0PHp0UpiTJSjKS7nzwqJoaatMy/Ldp0yZt2bJFGzeGbwzS1NSkuro6rV+/Xs3NzbEwE91fX1+vLVu2qLGxUfv3h/sTouGgvb1dLS2Ol2dMKhAI6ODBg7H6pN27d6e1Bmk6fr9f3d3d47ZF/+91dXU6ciTx3DW/369AIDDbzZNEoAKSGhge03On+3RZbaXKi/mTAXLJoRNd44b5JrKSOnqGdOhEl66vnz+rbWlvb9f+/fu1b9++2LZNmzZpz5492rNnj6qrq9Xe3j6uNyXakxXtldm3b582b94sv9+vuro6rVu3ztX2RXvDZjuEbNu2zdFxa9eunbIHbteuXdqzZ8+0vWbV1dVpqynj6gAkMTA8pp+9cEbL55USqIAc09k3dZi6kONS0dLSIr/fP65X6fjx42pvb499HR3GivL7/bFAsHHjRu3Zs0fz5s1TY2OjtmzZ4moPUjS0Rdsa/7XbUi0U37Fjh7Zs2ZJ0uDOdBfpcHQAAc1ZNpbOJJE6PS0UgEIgN70XNtLbn4MGDam1tVUtLSyyUzMaw3H333ae7777b9cd1w/79+1VfXz8uTE0MolHRc54OBCoAwJx17apqLfaV6HTPUMI6KiOp1leia1dNWrXHdY2Njdq1a9ek7VMVU0+0d+9ebd26VY2NjWpsbNTWrVt10003zUqgam1tndUC9Asd8ov27kW3RZeJqKurk9/vjxX6x0tHQbpEoAIAzGFej9HODQ26/Z5WGWlcqIqWoO/c0JCWgvT169dr3bp12r9/f6zoXJLuv//+KYeu4muZAoFALFRFTQwPUw1xRQvYncwIjM4enGjHjh265pprYnVeXV1d2rZtm5qbm7V+/XqtXbtWzc3N8vv92rRpk5qbm1VdXa2DBw9qx44d49p6IUN+ra2tam1t1caNG2PDpPv374/9n+64445x/8f4femQFetQpQPrUOFCBQZH9NMXzug3L1kof1lRppsD4AJk0zpUO3bsUH19vaqrw71iGzduVGtrq3bt2qX9+/erublZ27dv1+7du7Vr1y7V1dXpjjvuiIWl6Pe1t7dr69atsZ6Z/fv367777lNra6u2b9+ua665Jhbc9u7dqx07djha+mD//v1qb28f1/O1d+/ecSuPRxfUbG9v17Zt27Rt2zatX78+9thr167Vww8/HGtbU1NTbPbghQgEAlq1alXCYvn4HBNdkkKSDh8+PG7JiQvldB0qAhUAIC8EQ1b3HX5Fn/1em774vjXacs2KObNSuhPRXp0LqSnatm2b/H6/mpqaJIVDXbRofffu3Tp8+PC42Ytr164dt7hookU3c4XTQOVJX5OA3BQMWQ2OjHFzZCDHeT1GVyzzS5KuWObPqzAlhYfMLrRAe+3atZLCw5bxNUnRmqW6urrYOllR0d6kQCCg6urqnAxTM0ENFZDEuf5hbj0DIOelsoTA1q1btXv3bu3du1fV1dWqq6vT3r171dzcrCeeeEKNjY1au3aturq6YnVL999/v6qrq3X48OE5e7PmeAQqAMCc19k7pM6+YR3r7Jek2OeayuK8eKPU3t6e8iKgE2cTRmcaSolXMo/uiy/An8sIVACAOe/ex1/R1x5+Mfb1p+57SpL0yZsu1qebLslQq9InXWsxxXO6HMRcQaACAMx5t163Qk0NiyZtr6kszkBr5rbW1la1t7dr7969WXMvwHQgUAEA5ryaqpK8GNrLBo2NjZOG//IBgQpIYkFFsf74hnoVepgUCwBIjEAFJOHxGBV7vJluBgAgi/GWG0iie2BE3209qe6BkUw3BQCQpQhUQBKjwZBePjeo0WAo000BAGQpAhUAAECKCFQAAAApIlABAACkiEAFJFFRUqAbLqtRRQmTYgEAiXGFAJIoKyrQVcv9mW4GACCL0UMFJDE0GtSzHb0aGg1muikAgCxFoAKS6D0/qgNtp9V7fjTTTQEAZCkCFQAAQIoIVAAAACkiUAEAAKSIQAUkUeD1aLGvRAVe/lwAAImxbAKQRHV5kX7v2hWZbgYAIIvxlhsAACBFBCogic7eIX3l4Avq7B3KdFMAAFmKQAUAAJCirKmhMsb4JW2WtMla25Rg//bIP+slyVq7LX2tAwAAmFpWBCpjTKOkdZL8kqoT7G+21u6I+3qPMeZgouAFAACQblkx5GetbbXW7pXUPnFfpOeqMfI5ao+k9caYuvS0EAAAYGpZEagcWCcpPjxFg5c//U1BvqkuL9KH37ZS1eVFmW4KACBLZcWQ33SstQFJ8yZsXh/5PKlHK8oYUyypOG5TpbstQ74o8HrkLyNMAQCmlis9VBPdIWlbJGxNd0xP3MfJNLQLc1DP4KgOtHWoZ3A0000BAGSpnAtUxphmSfdFaq6ms0uSL+5j2Wy3DXPT8FhQz3b0aXgsmOmmAACyVNYP+cUzxmyUdNxBmJK1dljScNz3zmbTAABAHsuZHipjzHpJioYpY4yfWX4AACAbZFugmrQGlRRbp6pRUqsxpi4SpLZK6kpn4wAAABLJiiG/SEDaKGmLwmtONUs6bK3dH1l/6mGFl0hojv8+a+3uNDcVeaisuEBvqZuvsuKs+HMBAGQhY63NdBvSwhhTJamnp6dHVVVVmW4OAADIAb29vfL5fJLks9b2TnVctg35AVlneCyol84OMMsPADAlAhWQRM/gqL735CnWoQIATIlABQAAkCICFQAAQIoIVAAAACkiUAFJeDxG/rJCeTystg8ASIyFdYAkFlQU68NvW5XpZgAAshg9VAAAACkiUAFJnOkb1j//9LjO9A0nPxgAkJcIVEAS1lqdHwkqX+4qAACYOQIVAABAighUAAAAKSJQAQAApIhABSThLyvSlmuWy19WlOmmAACyFOtQAUkUFXi0xF+a6WYAALIYPVRAEn1Do/rpC2fUNzSa6aYAALIUgQpI4vxIUK0vd+v8SDDTTQEAZCkCFQAAQIoIVAAAACkiUAEAAKSIQAUkUVLk1ZXLfSop8ma6KQCALMWyCUASVSWFuvGyRZluBgAgi9FDBSQxGgyps3dIo8FQppsCAMhSBCogie6BEd37+CvqHhjJdFMAAFmKQAUAAJAiAhUAAECKCFQAAAApIlAByZjwDZJlMt0QAEC2MtbaTLchLYwxVZJ6enp6VFVVlenmAACAHNDb2yufzydJPmtt71TH0UMFAACQIgIVkMS5/mF9+7GXdK5/ONNNAQBkKQIVkEQwZHWuf0TBUH4MjwMAZo5ABQAAkKKsuZefMcYvabOkTdbapgT7t0sKRL70W2t3p691AAAAU8uKQGWMaZS0TpJfUnWC/dslyVq7N/L1emPMHmvttnS2EwAAIJGsWjbBGLNR0h3W2rUTtndLWmWtDcRts9ZaxysDsWwCLtTQaFCnAue11F+qkkJvppsDAEijObNsgjGmTuEhvkCCfevT3yLkm5JCr+oXVhCmAABTyvpAJaluiu0BhYcIEzLGFBtjqqIfkipnoW3IAwPDYzp0oksDw2OZbgoAIEvlQqCaSpcS1FvFuUNST9zHyXQ0CnPPwPCYHj12lkAFAJhSLgeq6cKUJO2S5Iv7WDbrLQIAAHkpK2b5JdE+xXb/NPtkrR2WFFva2hjubAsAAGZH1vdQWWvbJQUixekT97VkoEkAAADjZFugmmoYb5ek2Iy+yPIKe9PSIuS94gKvLl5UoeICZvkBABLLinWoIr1PGyVtkdQoabekw9ba/XHHbNcbQ3zXWGt3zPA5WIcKAADMiNN1qLIiUKUDgQoXKhiyGhwZU1lRgbweavEAIJ/MmYU9gUw71z+sb/78hM71Dyc/GACQlwhUAAAAKSJQAQAApIhABQAAkCICFQAAQIpyYaV0IKMWVhbrEzeuZoYfAGBKBCogCWOMCryEKQDA1BjyA5LoHhjRviOvqntgJNNNAQBkKQIVkMRoMKST3ec1GgxluikAgCxFoAIAAEgRgQoAACBFBCoAAIAUEaiAJCpLCtXUsEiVJYWZbgoAIEuxbAKQRGmRV2uW+jLdDABAFqOHCkji/EhQbad6dH4kmOmmAACyFIEKSKJvaFQHj76uvqHRTDcFAJClCFQAAAApIlABAACkiEAFAACQIgIVkESh16Nl80pV6OXPBQCQGMsmAEnMKy/SpnXLM90MAEAW4y03kIS1VmPBkKy1mW4KACBLEaiAJM70DevrPz6mM33DmW4KACBLEagAAABSRKACAABIEYEKAAAgRQQqAACAFLFsApDE/IpifeQdq1RWxJ8LACAxrhBAEl6PUWVJYaabAQDIYgz5AUn0DI7qB0+/pp7B0Uw3BQCQpQhUQBLDY0G9+Hq/hseCmW4KACBLEagAAABSRKACAABIUU4VpRtjtkrySwpIqpe0y1obyGCTAAAAcidQGWO2S9obDVDGGL+kuyVtymCzkAfKiwv0ttULVF6cM38uAIA0y6Uhv6b43qjIv/2ZagzyR3lxga5dVU2gAgBMKZcCVcAYczDSMyVjTJ2k9sw2CflgaDSo42f6NTTKLD8AQGK5FKg+KqlOUrcxplnSemvttqkONsYUG2Oqoh+SKtPVUMwtvedH9f2nXlPvedahAgAkljOBKjLE1yxpv6TtkjZFe6umcIeknriPk7PcRAAAkKdyJlBFeqXarbWbFJ7hVy3piWm+ZZckX9zHsllvJAAAyEs5Eagi9VJ+a22LJFlr2621axWuq9qY6HustcPW2t7oh6S+NDYZAADkkZwIVArXTgUSbN+T5nYgD3k9RvMriuT1mEw3BQCQpXJiHri1tsUYs8MY45+wkOfa6QrTATfMryjWh65fmelmAAAuQGfvkDr7hidtr6ksVk1ViWvPkxOBKmKTpDuMMecU7q3yS9qRyQYBAIDsdu/jr+hrD784afsnb7pYn266xLXnMdZa1x5MkowxK621L0X+fbXCw3VPRLdlSmTphJ6enh5VVVVlsinIMZ19Q9p35KQ2rVummkr33s0AAGZftIfqWGe/PnXfU/rqlqu0uqbCcQ9Vb2+vfD6fJPkiNdkJzUYN1froP6y1T1prvxO/Dcg5VhoZC0nuvvcAAKRBTVWJ1iz1aXVNhSRpdU2F1iz1uTrcJ7k05GeM8UnarPAlp8mYccW7fknXSPqmG88FAACQbVwJVNbaHmNMi8I1TfWSuuN2ByR9xo3nAQAAyEauFaVba09I+pgx5iZr7cPx+4wxK916HgAAgGzj+iw/a+3DxpgbFR7qi9oS+QByzrzyIt163QrNKy/KdFMAAFnK9UBljLlf4TAViNt8tdvPA6RLodfjevEiAGBumY11qA5aa++O32CMuWkWngdIi96hUR15qUvrVlarqqQw080BAGSh2Vg24ZzDbUBOGBoJ6tev9mhoJJjppgAAstRs9FDVG2P+W1Jr3Lb1Ci+dAAAAMOe4Gqgi61FtkXTfxF1uPg8AAEA2cTVQRdaj+qi19sn47ZE1qgAAALLOdDdQdjolaTaWTXgywebuBNuAnFBa5FXjRfNUWuTNdFMAADMUfy8/SbHP8ffym+4Gyn90Xa2j53Hl5sjGmPdLarHW9hpj/iLBIU3W2t9O+YlSwM2RAQDIP185+MKUYenTTZdImv4GyiUacXRz5AvuoTLGfMRaG70/32cVXnfqx5J+T5NrqOZf6PMAmTYyFtLZ/mEtqChWUcFsTIwFAFyI6Ybqor1Pt163Qk0Ni/TIc5360sEX9OdNl+iGy2pUU1n8xvFVJePWG4zeQFmSentHHLUllSG/PYrc8Nhauy5uOzVUmFMCgyO67/CruvW6FSzwCQBZZLqhumjvUzQsRYf6lleXxcKSm1IJVAln7iWqoZqirgoAAGASJz1P0hu9T4mG6tItlUCVevEVAADABE56nqTph+rSLaUeqkgBeouk9ukKtQAAAJzKpp4np1Ltodovaa2kzxpjVknqUniF9MOKzPpLvYlAZhljVFrklTGsTwtgak6GqZwOZc1lTs5BNvU8OZVqDVWXtfY7kr4jScaYqxVeKf2zkpolXZxyC4EMW1hZrI/9Zn2mmwFgFjgNOE6OczJM5eSYuR66nA7n5ZpUAlWLwr1Tj0Q3RIrPKUAHAOQEpxd3J8c5GaZycozTNuVq8MrF4TwnUglUmyTdbYw5Ya19yaX2AFnnbP+wHvz1a9pw5RItqMjtP3ggn8xkjaJkF3cnxzkZpnJyjNM2ZWNPz1wdznPiggOVtbZH0mZjzE2SXnKtRUCWCYWsAoOjCoWY2ApkC7eG4Jxe3NMZApw+l5Pgle5erGwMeemS8r38rLUPu9EQAACccmsILpc5CV5u1my52eM3F7l+c2QAAFLh1oV7rg4tzYSbNVtu9vjNRQQqAEBW4cLtHjdrtvK598kJAhWQhK+sUO+7eql8ZYWZbgqQF7hwp1c21pHlIgIVkERxgVcrF5RnuhlAznNaq8OFG7mIQAUk0T88pv852aM3L/Opopg/GeBC5fMMMMx9XB2AJAaHx/Sr9nOqX1hOoAJSwFAe5jKuDgCAtGAoD3OZJ9MNAAAAyHX0UAEAUpar95UD3JJzgcoY0yzpeOTLLmvt/ky2B3NfcYFXly+uVHGBN9NNAbIWBefIdzkTqIwxfkkPS7rJWhswxjRKekKSyWjDMOf5ygp185rFmW4GkNUoOEe+y5lAJalZ0n3W2oAkWWtbjTFNmW0S8sFYMKT+4TFVFBeowEvZIZAIBefId7kUqLZKqjfG1Emqs9a2WGtbMt0ozH1dAyO69/FXdOt1K6gFQV6iPgpILicCVSRESVKjpHZJ7caYPZL2TRWqjDHFkuL7mitnt5UAMDdRHwUklxOBSlI0UAWsta2SZIzZIemEpHlTfM8dknamoW0AMKdRHwUkl2sFIUei/4jUUvmNMeunOHaXJF/cx7JZbx0AzEE1VSVas9Sn1TUVkt6oj2K4D3hDrvRQtU+xPaA3eq/GsdYOS4oN+hvDZEAg2zmp1XHrmJkcN5dxDpAPgiGrFzv7JEkvdvYpGLLyetzNBTkRqKy17caYdoXDU2vcLr/ieq2A2VBTVUKdyDTcDC9OanXcOsbpcXM9cFAfhbnuQFuHvvDQszrZfV6S9I1HjuuBp17TX737cleXxMmJQBWxQ9IWRQKVMWajpJZoTRWQD9wML24d42Z4cVKr49YxTo+b64GD+ijMZQfaOnT7va266bIa/cPvX61LF1Xq+df79E+PHNPt97bqrlsbXQtVOROorLX7jTHVxpjtkU3zrbWsQ4VZ1zUwoh89c1q/9aZaVZcXXdBjZGN4cesYN8OLk7WM3DrG6XFO2p3uXiw3hz1ZPwpzVTBk9YWHntVNl9Vo7wfXyRMZ4mtcMU97P7hOW//tiP72h8+qqaHWleG/nAlUkmSt3ZvpNiD/jAVD6ugZ0lgwlHB/Ooey3Awvbh3jZnjJRk7a7TTopjNYz/WeNSCZQye6dLL7vP7h96+Ohakoj8fo9neu1gfu+qUOnejS9fXzU36+nApUQDZK51CWm+HFzZ6efOc06KYzWDOUh3zX2TckSbp0UeJlKC+trRx3XKoIVMAUoj0FXQMj6uwb0nOn+9TZNzxpyCSdQ1nITk5/dukM1vw+Id/VVIZ//59/vU+NKyYvWfn86b5xx6WKQAVMYWJPwX8celXS5CETLlxwiiAEt4RCVh6PkbVWgcERSdLrvUOqKC7QaDCkVQvKVeD16Fhnv871D+vFzn5J0qET51RRXKCVC8p1sntQj7d3aSwU0ljI6tWuwXHPcddPjmtkLKRgKKSgtRoLWf3R21fF9n/jkWPylxUqGLIKhqR3XrpQG65cohdf79Ou/3pWkvSFHxxVRUmBqkoL9eXNV0mSPvmfT+pM37D6hkYlSdv3/1q7N16pNUt9uudXL2v/EydlrZWVNDA8Fnu+V84Nats9T8haK0myVvJ6jH74yXfEjvnje59QSaE3st9qYWWx/umRY+NqqKLn766fHNPy6lJdu6rajR8JgQqYSrSn4Ohrvdr+nae1+wNXqGFJFUMmQB7rj1zgO3qGVFzg0fBYSMury+QrLVT7mX698HqfjkXCy0P/06Gz/cN656U1CgyOaM/P2jUyFtLIWEine8JT+EOhcDj4/A+O6uhrvRoJhvePBkP6QGN4PeqHn31dH/qXQxodC2kkGN53ff183fuRt+j8aFC3/d9DkqQ/+v/eWEXo0F/epJrKEu0+8Jx+dPT12Pa/+cGzGg1abfvNerWd6tGf7/u1JMljNKnO6J5fvayxUEheY+T1GhV4PNq8bnls//On+1RRUiCPMfJ6jAKD4eAftFaDw8HYv8OP/8Zj+0oLFbJSoTe8bZGvRMUF4XXGF1QU6bLaShkTXj+ye2BEx88MSJJKijy6dmW4pym6tuTEYvKrls/Tgoo3XqNrq4q168Bz+ui/HdEfv3O1Lq2t1POn+3TXT47p4ec6ddetja6tR0WgQl5yUvg7saegYUkVPQVAlolOFjk3MKxTgfNa6i+VJP2q/Zz6h8Y0OBrUsdfDQztdA+G/+YNHX9ejx85qaDSo86NBnR8Jqj6yCvwr5wb0Z/c/paHRkIZGgxoeC6nQ69GRvwrflOMv7n9KkvTRb78RXu7+0Do1NSzSf7Wd1t/99/Ox7Xt/elwdgfN656U1Gh4L6aGnO1RU4FGR16OgDbc7GjjKiwu0sLI4vD9yjL+8UJJ00fxyfeQdq1To8ajQa1RY4NESX/j/WVzg1V/+zuX62x8+q8+/901aXVOpQq+RrzT8vZ+75U36zLsuU/uZAX3k20f07T+8RtesDBdgNzXU6oUvvEsFHiOPx6jtVI/e8/VfxNr/6GduTHjO2071SJL+4fevTviaeFltlT7/u2v0nq//Qjs3vGnSMX/z3jWxx/npC7/QnzddqosjdU43r1k8bhmDtlM9+q+205LCQ3N3Rr53Klt/o27S862YX6YvPPSsPnDXL2PblleXurpkgkSgQp6ayQyoodHguM8AZm5kLKTBkTH1D4/pdM8bRcAHj76uwZExDQwHY58/eP1Fqi4v0g/bOiRJf/1Amwo8RoMjQW1cu0wfftsq/fL4Wf3BvxzSaDAcSP7gXw5rqb80FgI+8R/hYaV40edtP9OvX7WfU0mhVyWFHpVGhogkqay4QG+tXxDbV1LoVXnRG/s/9s7V+qv/16a//d01umxxpYoLvFoxv0yS9OG3rdSt163Qsc5+bfznx/TAn7w9dnFfVFWin22/IfY40fBS6A33zPxZgpmX0eCyuqZCv3v10oTn1esxsRlqV6+YNylMLIkEzMGR8OtXdXmxSiP/H6/HuL5aeDa6ec1iNTXU6ssHn9c3Hjmuj99Qrz9rujQ/V0oH3DaTGVDRF6LoZyAfWGs1GrQqKvAoGLJ65rUe/c/J8AX+x8+9ridf6dbGtctVWuTVvY+/rNaXAxoYHtPrveHQ8vCzr2vNUp9+9Mxp/cm/P6mRuGVHVkYCiCR94j9aNTQakjFSeVGByoq8evcV4TXfhiJ/cyUFXi32lai0yKtl88LfW7egQn/9ngad6x/W1x4+pp0bGvTmuDDxwMffpkKvR2VFXh3r7Nd7v/GoGpaE92/7zXpt+836cf/ftlM9uusnx7Wgolifu+VNU56Xq5b7JUlXLvdPCi9lRQUqK1KshgfZw+sxurgm3At2cU3lrARJAhXmnAsZzqPwF3PJWDCk/uEx9Q2Nxep5JOknz3fqpbMDsX29Q2O65colur5+vn7Vfk6S9OFvHdZI5PvXXjRP92+7XsGQ1S3/+Gjscb588EUVeIxuvHyRlhaV6tWu83rp3IDKiwtUVhwOE9Ehp0trK/WX775c5cUFKi/yqry4QF0DI/rUfU9Jkn75mZtUGukNmnjP1fc3LtO/PPqS/vLdl0/6+6z1lehD169U26kefe3hY7pmZfW4Y6I9M9LkOhtgNhCoMOewoCFy3VgwpN7z4RlQL77ep8DgqFbXVKjWV6KnXg3o4NHT6j0/ple6wsW633r0JX1p85XqHhjR25t/rIEpelP/9Zcv6ZfHz6myuECVJQWqLCnUOy5eIEnyl4UD0A2XLVTdwgpVFBdo6bxwKCkq8OihP327Xguc10e//YS+c/v1alwxLxaAPvOuy2LPER3KWrcyPHPqovnl+oO3lo9rR3QoS9IF330AyDYEKsw5LGiIbBAKWfUNjSlwfkQ950cVGBxVVWmhrlruj00X/1rLC/J4jHrOj6r3/Jh+8Im3y+Mx2rznMbW+EpAkffr+8Cysv9t4hTatW67jnf36f0++pqrSQkUmRqk80itUUVKgTzddoqqSwlhgOtM3FHuMb35onQoiNTsTXVZbJUn60PUrE/bWvmmJT5H6aRUXeCf1JgH5jkCFOcft4byCyHBBAcMGeSkajF4LhKe5t77cHft9+vZjL+l4Z78C50fVPTiqjsgxkrTnZ+1qPvDcuMdaf/kiffMP1ikyU14vdQ1qsa9UCytLtLqmQCPBkEo8Xn1q/SV6rqNXX/yv5/TVLVdp7UXztDDyhuADa5fpA2vD0+mjvUHRqeyFXo8+8o66cc8Z3xs0VZgCkDoCFXJKum9CK0lVkVqQ6GfkrmDIquf8qLoGRrSoqliVJYV66tWAvtd6SpL01ZYXJBldvcKvj9+wWi+dHdCNX/pJLABJ0v/5/jO69S0Xyesx+vmLZ/Vq16B8pYWaV1akVQvKYwsoNjUs0kXzy+QvLZSvrFD+siLNiwyrReuLvrL5qoRh/zcuWRgbCltdU6Hl1WWTjgGQXQhUyCmZqI8KRcY5op+RHWzcz2NgeEyHXupSV/+IugdHdG5gROcjU+wl6f880Kb2swPqOT8aG7ba88G1+u031erQiXP6TutJSdKr3YNa6i+LzdKqqSrW37x3jarLi9Q9OKK//F6b/vXD1yjaWXn3h9aNa1PbqZ7YIoqrayq0OrK2EYDMCoasXuwMr0f2YmefgiE7abJCMGT19MmAJOnpkwFdvrhqRhMaCFTIKZmojwoMjo77jNlzfiSos/3Dej6yEOOPjp5WdXmRlvhL9YOnX9P+J07qXP+IugZGdKb/jbWMXu8d0oe/dViSVF7kVXVFkarLi/X+yNo916ys1k2XL1J1eZHmlRWpurxIlywKh52tv1Gvt9Yv0Hu+/gt9adP4HqOyogLd9paLJL0xdLagopj6oRzi5CLp9ELq5mPBPcnO+YG2Dt354FF1RNYh+8Yjx/Xd1lPauaEhtrDnxGM++702ff3H4eU43rpi/KSKqRCokFNY7iC3hEJWgfOjClmrBRXF6hoY0QNPhYfX/uHhFxWyVudHg7r3I2+RJL3n6z+P3WYifMwxrVni0xJ/qUJWKvJ6tGZplarLizQ8FtI3f35CkrS8uky//MyNqi4vGrcGUDQEbbhyCb8nc9BML6TxF0knF9L4VbTdfCxClzNOzlOyc36grUO339OqieMLp3uGdPs9rbrrtkZJmvaYv7tltaP2EqiQFTJRG4UL1z0wos6+YZ3pG9bZ/vDHRfPD7+J+9Mxp/e/9T+tc/7C6BkY0FrJ6f+NSfXnzVeo9P6p//eVLkqQTZwe0vLpMS3ylsRu97twQXlAxMDiqP/3PJ/W9P36rro7cJf6WK5foliuXxNrQdqonFqgKvZ5x6w4h96Ualty4kN51W+OsPBahKzmnAXa6c/6N/3W1Pv/Qs5P2S5KVZCR97vvPSDLTHvP5h446ajOBClmBtaMy7/xIUIHBEfnLitQd6Uk60x8OTcfP9I87dvOex2LF11J4mO0vfvtSSeH1jNZe5NeCimLNryjWwooi1S0MD69dNL9M3739rdrwj4/qK1smF2T/xiULJb3Rs1TIrLS8lGpYcutCeueDR3XjZYt054NHXXmsUEj6+L+7F7qk3A1e07XbSYBtaqhN+nP5qwfa1DUwdamGlXS6d/Ib+YnHdA+OOfo/EaiQFVg7anZEC7dfOjeg7sERdfYO6y3187XUX6oHnjqle3/1il7tHpQkbdrzmN539VJ9ZctV6hsa0xd/+JwWVhZrYWVx7E7woch0t79935vl9RgtrCjWgsoilRUVxELQtavm6w/fXpegNeE7xFN/lN+c9DylGpbcupB29Azp3x57KRZsUn2sv3qgzbXQJWXvEGMqvYtOgtKdDx5VZUlh0p/LdL8Ds4FAhayQzbVR0SnuvmmWTXCrWHWmL3xDo0H9/MUzer13WK/3Dqmzd0iB86P62u9dLUn65H8+JUn6k39/MvY9/3Rro5b6S1Vc4NUSf4mWzSvVd588pb/4rUv0zktrJIXvxP78F26OhZ/oekeeSFuuXVU9ZZuAqSQLAMGQdaXXwc0L6ctdg649VtfAyJT7ZhK6mhpqdfDo6YwMMc72UOyn1l+cNCh19AzpsePnpjwmUwhUmHW5XB8VvSmsJD3zWo/WLPXNuCgy1WNuvW6FugdHI4FpWB0952M9d79qP6e//9ELksJDbYsqS1RTVazRYEgeY/TW+vlqPzugj/3GKv3etRep1lcSK9q+eU2tbl5Tq7ZTPfruk6f0zktrYiGWXiS4zckwjq+0KOt6HS5K8xpgTkLXr46fc9SL4/YQYzqGYr/16EvTnZ4J35FcdXm4hCHR0UbSoqpiSUav9w5Ne8yrDp6LQIVZl631UW7NELrQYtWOuGP+49Cr+ukLZya18XTPkP7+Ry9oUWVxuA7JhF9wX+0Or8j99z96QQsri/XX775ct1y1dMq2//PPTuiBX3dc0LtSporDqal+D5z0PN354FFtv/myBEdcGDcupLW+En3w+pX65i9O6HTP9BfcZI81r7zQtTD4WPtZR704bg4xpmsoNnDe2Tm6vm6BvtN6atqfS62vRH/97gZ9/N9bZTQ+gkVfnT53S3gizO33TH3MZ951mTb+TfI2UfGJWXfrdSv0g0+8XV/dcpUk6atbrtIPPvF23Xrdioy16UBbh97e/GN99nttksJh6e3NP9aBto7Y/tvvaZ30ohV98TjQ1pH0ImElbd//tD73/cTHRH3u+8/oiZe7E+6LvhB5PEYffMtFOnyia9KNb8/2DeuT//nUjNru5By4eUzUxAtuMDT5zLh1DNJrut+DQye6HAWArv7pa5GiqsuLNFVcN5IW+0r0hfeuiX09cb8UvpB+7paGaY/ZuaFBRQUe7dww/XFOHusL712jxb6SadtdXe70bgzO3qykOsQohUPXyFho2tc6KRyW3Opd9JcWJv35vqV+ftKfy84NDfqdKxbrrtsaVesbPxpS6yuJhcWb10x/TFNDraN2E6gw62qqSrRmqS+2anS0Pmo2h/umu+AmCxw/fPq1pC8edz54VL9qPzftC4gk9Q6N6XRv8oLW/uGpZ5G4/cJ354NH9cOnk4cuJ8HMaXiT0h/gnIYuwlnqkv0etBw97ehxqsuLkoYOJ2HJrQtptDfXyXHJjvmdK5YkDQBOQtdiX4mur58/xREz5yR0OSnMd3Mo9sNvWyVp+p+v12Nm9PP7xY4b9fEb6iVJH7+hXr/YceO43vroMV98X/h364vvWzPpmGQY8sOck+oMEifd0pkoinTrhc/JUICTaeBOp4o7LaCVkq/j4+QYN2tC8IZUhvO+F1nMNZlaX6l2bmiYdvgl+rO5y9M47mcX/v6ScT+7m9csVlNDre47/Io++702ffF9a7TlmhXjhqOdHOPWY0UDwHTt9nhM0v//W+rma7GvZNrhLjeHGN0szE82FFvrK9Gf3Lhal9ZWJP35Ss5/fl6P0cU1lZKki2sqE5YkeD1GVyzzS5KuWOafcdkCgQopybaCczdmkDh/Ecq+ngynL3zJwpmTaeBOp4o7KaB1M8C5VRMS/w4332vEpgueTgvJnVxIr11VLa/HJA0d0swupMkukk4vpG48lhuhS1LS4PmF967R5x961pXQ5bQw38nPOFlNU3zvk5Ofr5R6EHILQ35Iyb2Pv6L3fP0Xkz7uffyVWXvOqYZokr1TlmYygyS56+sWJO2er60qVm1V6sc4ra1I94wkJ5wU0J7uHXY0NJrsGLeHRoMh6/oQY65xazjvd68Kr3KfbBhHcj78ki0X0plyErqS/f/TOcT4wetXpnUoNirXfr4EKqQk3QXnqRa+Op1B4qTo1UlRpJNiVbcKWp2+8DkvfHVLel8E3Rwa/ccfH3OtwD+bpfImxelwXlNDreMLqZR7F1O3Ofn/JwteboQup4X5MwlLbtQrZSMCFVLiZsF5snf4br1TdjKDxMk7LadFkekqaHX6wucknLnVs+Z2Aa1bnA6NfuvRE64U+MdL50xHJ8ek+iYlOpyX7Pfg2lXVc/ZCmkmp9na5WZjv5PmctjsXEaiQFZK9w3fznbKTGSQz6ZZ22j2f6jFuvfA5CWdu9azFF9Dm4tDodD2aMxlijAaZdM50dHpMuofz5uKFNNu5McQ4k+Py9WdMoEJCnb1DajvVM+mjc5oalgvlZOq9m++U/+TG1TN6p3XHu8KLDd7xrsumfDftVuFrul740tWzdvOaxfJ6TE4OjfqnudVQPCdDjIdOdLm2DIVbx2RyOA/Zx83C/HzFLD8k5Pbq5ulaQfl3r1qibz36kmszSLweo8sXV0lS1szscuuFz60p5U6PcTJzyY1jnEw7jw6NTnfch9+2Ul9pmfw3cCFO95zX7v9+Pi0zHZ3OhnR6c9mZzM5zOisLmIsIVEjo1utWqKlhkY519utT9z2lr265SqtrKmL3kJuJVKdcz2QF5aaGWl27qtrR+iW800pfz5qUvgDnNLwlO66poVb/efhVV6addw2MOJrpOB23l7Nwuo6a0zcpEn9TyG8EKiRUU1UyrrA8Wmw+U8nW+vnDt6109DjRFZST3bdpNt4px18skJp0Do26sVijG2v91PpKVF0x8zcis8/Zsg4zeZMC5LOcraEyxhzMdBsQlo4p19EVlKX0F776IrU0Poc1NcgObgyNujX7sjYDi9wm42QdNWbnAc7lZKAyxmyUtD7T7UD6p1xT+Ip0c6PA/9pV1Wmb6eh0OQunN5dlOA9wJucClTHGL6k60+3IZW7N4MvmFZTdFJ3FNd1sLsxtqc6+TOdMR6fLWczk5rIAksvFGqrNku6XtCfTDclVbszgc/OGqDOt0eCdMrKRk9CVrpmOTo+JtovZeUDqcipQGWPWS2pxeGyxpPhK0MpZaVQOcjqDb7obws5kOI8p10BYupeqyLWbywK5LNeG/PzW2naHx94hqSfu4+SstSrHOLldTLJVljv7nA0PsoIyMF46l6rgbwpIn5wJVMaYrdba/TP4ll2SfHEfy2alYXOQk1WWayqdzVpiBWUAQD7IiSE/Y0yjpCMz+R5r7bCk2Ap3xvDOzAmnK5f/9H/fkLF1odKtqqRg3GcAACbKlR6qaknrjTHbjTHbJTVLUuTrjZltWm6aau0oJ7VRHT1DeuLl7ryZcl3g9Yz7DADARDnxltta26K4YvRIj9VWa+3uzLUqO3X2Dqmzb/KtJ2oqi2M1UtPdCmZ4LOTsefqG9N6rljqeSZTL+ofHxn0GAGCinAhU8SI9Ulsi/26WdDASuKDkSyIkuxXMp9Zf7Oh5ojVU+TDleiQSMkcchk0AQP7JuUAVKUyfSXF6XpluSQQn9VH/cegV1VaV6PXe5LVRUbk8nDedaG/fibMDkqQTZwdUXV40rrcv/rhjnf2SFPs88TgAwNyVc4EK05vupsaPHT/n6I73n15/ib7a8oKju8vnKichaGJv387vPyNp8gKoE4/71H1PTTqO0AUAcxuBag6aakFOp2tHrVxQlrW1UU6CyYWEpUQhKNrb1zUwov9q69C71iyO9VDFix43UfxxTp4PAOC+6DXh1a5BSdKrXYNqO9Xj+htaAtUcM13BudO1o2oqS3R9/fy01kY57cFxEkxmEpYmig9B0d6+vqFReYzRlct9qiwpnPw9E3oFE3HyfPRiAYD7Jl4TvnTwBX3p4Auuv6ElUM0hyQrOv/G/rna8dpSU3toopz04ToLJTMKSE5UlhXr7xQscHTsVJ8/H0CEAuM/p7dZSRaCaI5wUnH/+oWf11+9u0Mf/vTWt9VFOQoCTECQ5CyYzCUtODI8F1dk7rJqqYhUXeF173IncHDokeAFA2HS1xW4iUOWQ6daYOn5mwNGCnPPKi9JeH+UkBLgdgtzUMziq/U+c1K3XrVBN1ewFKreGDiV6uwAg3QhUOWS6NabqFpY7eozogpxu1Ue52fuE5JwGTwrlASC9CFQ55NbrVujGy2r0o6On9Y1HjuvjN9TrtxpqtdhXouNnBhw9RrQw3a36qFzvfZqrKJQHgPQiUOWQ1le6xw3VfeOR4/pu6ynt3NCgpkiwclpwnozTCym9T7kr3YXyhDMAcxmBKkckm8F3122N2rmhQbff407BudPhoHzofTLGqLKkQMbk9mKmF8LNoUPqugDMZQSqHOBkBt+dDx7VL3bc6KjgnLqnmVlYWayPvKMu083ICDcL5d0KZ4QuANmIQJUDDp3ocjSD79CJLkc3K6buCW5y+rviVjgjdAHIRgSqHOD0ljHR45IVnNP7NDNn+ob1wFOn9N6rlmoh52hWpTN0AYCbCFRZYLr1pWqqSmZ0yxhHx9H7NCPWWvUNjcnaRIOuSDc3hyHpyQLgFgJVFphufalPN12ia1dVO57BxwUCcP6mgeHD7OTWTdD52SGdCFRZYLr1paTwEJ7TGXwMdQDOMXyYXum+CTrLfiCdCFRZYLr1paIz825es9jRDD7qowDnWADVPU7OQbpvgs6yH0gnAlWGOVlfKj5UJZvBR32U+3xlhdq4dpl8ZYWZbgoywK0FUKW5fcF1cg7SfRN0lv1AOhGoMsjp+lJNDbXyekzsD7GsKPxjKysq0LMdvfwhzrLiAq+WV5dluhnIYnP9ptVurV2XjW/4WPYDbiFQZdBM1pe6vn4+tRwZ0jc0ql+/2qMrl/tUWUIvFSZL902r3SzIdmuoLhvDUrqle9kPgld2IVBl0EzXl6I+KjPOjwR1+KUuXbKogkCFlKTzgutmbRCvPe5xcxiS3q7sQqDKoJmuL8U7QGDuc+uC62ZtEK896ZWNPZ7p5mYPa6rHOD0DBKoMmsn6UgAQ5VZB9kyOQ/bJxiFGt45xs4c11WP+6LraSecvEZMvqz8bY6ok9fT09KiqqirTzYmJzvKTEq8vFT/LD5nR2Tukex9/Rbdet4ILD4A5J9ndOqK+cvCFaRehdvMYp21yclyqx5RoRD6fT5J81treSQdFEKhmWWfvkDp6hvTMaz3qHhzVvLJCvWmJT4t9b7yzuO/QK/q7Hz2vs/0jse9bUFGk//1bl2rLtSvS1lYk1nN+VIdOdOnaVdXylVJDBSA/pSO8ZOOb1t7eXkeBiiG/WXbng0f10P90TNr+7jcv1jdubZQkvdYzNC5MSdLZ/hG9Ns0MQKSPr7QwYZc5AOSTdK4PlosIVLPoQFtHwjAlSQ/9T4c2tHXo5jWLmUGT5UaDIfWcH5WvtFCFXk+mmwMAyEIEqlkSXbRzKvGLds7VtD5XdA+MUEMFAJgWb7dnyUwW7QQAALmNQDVLZrpoJwAAyF0Eqlky00U7AQBA7iJQzZLoop3TWcyinTnD6zHJDwIA5C0C1Szxeox2bmiY9pidGxq4UOeAmqoS/elNF1OQDgCYUk4FKmPM9sjHHmPMnky3J5mb1yxW8/vfrAUVReO2L6goUvP738wK6AAAzBE5s1K6MabZWrsj7us9kuqstU0Ovz9jt54JhqwOnehSZ9+QairDw3z0TOWOc/3DOvDMad38plrNr2BtMADIJ3NqpXRjjF9SozHGb60NRDbvkfSEMabOWtueqbY5CUtej9H19fMz1EKkKhiy6uwdVjCUG28+AADplxOBKmKdpDpJrZGvoyHKn5HWKLwS+p0PHh233tRiX4l2bmhgOA8AgDySEzVU1tqAtXaetbY1bvP6yOeEvVPGmGJjTFX0Q1Klm2060Nahj93TOmnxzo6eIX3snlYdaEt8yxkAADD35ESgmsIdkrbFDQEm2t8T93HSrSdOdlsZKXxbGYaIAADIDzkZqIwxzZLus9buneawXZJ8cR/L3Hr+ZLeVkbitzFxSVVqod1+xWFWlhZluCgAgS+VSDZUkyRizUdLxJGFK1tphScNx3+daG7itTH4pKfTqkkWujhgDAOaYnOqhMsasl6RomDLG+I0xdeluB7eVyS8Dw2N64uVuDQyPZbopAIAslTOByhjTKKlRUqsxpi4SpLZKSvu4GreVyS8Dw2P62QtnCFQAgCnlRKCKrEP1sKRmScfjPpqnKUqfNdxWBgAAxMuJQBW3bIKZ+JGpNt28ZrH++bbGST1Vi30l+ufbGlmHCgCAPJJzRenZ5OY1i9XUUMttZQAAyHMEqhRxW5m5r6jAo7qF5SoqyIkOXQBABhCogCT8ZUV671VLM90MAEAW4y03kEQwZDU4MsbK9wCAKRGophEMWT12/JweeOqUHjt+jgtqnjrXP6w9P23Xuf7h5AcDAPISQ35TONDWoTsfPDruFjOLfSXauaGBGXwAAGAceqgSONDWodvvaZ10v77TPUO6/Z5WHWjryFDLAABANiJQTRAMWd354FElGtyzkY87HzzK8B8AAIghUE1w6ETXpJ6piTp6hnToRNrveAMAALIUNVQTdPZNH6Zmehxy34KKYv3xDfUq9PD+AwCQGIFqgprK6W96PNPjkPs8HqNijzfTzQAAZDHeck9w7arqSffnm2ixL3yLGeSH7oERfbf1pLoHRjLdFABAliJQTeD1GO3c0CAjaeId+aLbdm5o4H59eWQ0GNLL5wY1GgxluikAgCxFoErg5jWLdddtjaqd0FNV6yvRXbc1sg4VAAAYhxqqKdy8ZrGaGmp16ESXOvuGVFMZHuajZwoAAExEoJqG12N0ff38TDcDAABkOYb8gCQqSgp0w2U1qijh/QcAIDGuEEASZUUFumq5P9PNAABkMXqogCSGRoN6tqNXQ6PBTDcFAJClCFRAEr3nR3Wg7bR6z49muikAgCxFoAIAAEgRgQoAACBFeVmUHgxZ1pcCAACuybtAdfDoaf39I4fV0TMU27bYV6KdGxpYAR0JFXg9WuwrUYGXDl0AQGLGWpvpNqSFMaZKUs+KT90vU1w2fl/kM7eVAQAA8Xp7e+Xz+STJZ63tneq4vHvLnSg+Rrfd+eBRBUP5ETABAIB78i5QTcVK6ugZ0qETXZluCrJMZ++QvnLwBXX2DiU/GACQlwhUE3T2cdEEAAAzQ6CaoKayJNNNAAAAOSbvZvlNtTiCkVTrCy+hAAAAMBN52UM1MVRFv965oYH1qAAAwIzlXaD68pYrVesbP6xX6ythyQRMqbq8SB9+20pVlxdluikAgCyVd+tQ9fT0qLyikpXSAQBAUk7XocqpGipjzHZJgciXfmvt7gt5HK/H6Pr6+a61C3Nbz+CoHms/q+vrFshXVpjp5gAAslDODPlFwpSstXuttXsltRpj9mS4WcgDw2NBPdvRp+GxYKabAgDIUjkTqCTdIWlv9AtrbYukrZlrDgAAQFhOBCpjTJ3CQ3yBBPvWp79FAAAAb8iVGqq6KbYHJPkT7TDGFEsqjttUKYWLy4CZ6Osd0s+feUW3XO5XiUYy3RwAQBo5zQ25Eqim0iVpqpU475C0c+LG5cuXz2qDMHd979OZbgEAIIMqJc2NWX4JTLes+S5JX05wfPTux5WSTkpaJqnP/aYhAc55+nHO04vznX6c8/TLx3NeKem16Q7IlUDVPsV2/1T7rLXDkoYnbI4lS2Ni6071TbeuBNzDOU8/znl6cb7Tj3Oefnl6zpP+P3OiKN1a2y4pEClOn7ivJQNNAgAAiMmJQBWxS1JsRp8xZqPillEAAADIlJwJVJFV0f3GmI2RMHWNtXZbCg85LOlOTR4WxOzhnKcf5zy9ON/pxzlPP855AnlzLz8AAIDZkjM9VAAAANmKQAUAAJAiAhUAAECKCFQAAAApypWFPV1ljNmu8H0ApfBNl3dnsDlzjjHGL2mzpE3W2qYE+zn/syByXiWpXpImzoLlvLsn7ndcCp/vOkkfjb+BO+d7dhljDk58feGcu8cYs17SNkkHFV5Au0nSYWvt/rhjON9x8q6HKnrRsdbutdbuldRqjNmT4WbNGcaYRoUvNH4luDUQ5392GGOarbW7Ix/bItsOxu3nvLurWVJL5HzuUPiWVvuiOznfsyuydM76Cds45+7yK3yO90Q+jicIU5zvOHm3bIIxplvSqgnvJK211kz9XZipyAveHdbatRO2c/5dFukt2adwj2Agsq1R0hOS6q217Zx3d0XC6sHoO/LIxeUOa+28yNec71kS1zu4J/58cs7dFXkNb4k/nxP2c74nyKseqsita/yJfkEi3ZuYRZz/WbVO4WGnqOg9Lv2cd/dZa5smDG9cI6lF4vc8DTZLuj9+A+c8vTjfieVbDdWkewFGBBTu3sTs4vzPgsiL2rwJm6Mvau0Kh61EAuK8pyzyTt4vaVNkE7/nsyRysU50/1bO+ezYbIzpUrh8oz4yvC1xvhPKt0A1legvDDKD8+++OyRts9YG4u4MPxHnPQVxQ09+SfumGhqJw/lOnT8yhO13eDzn/MK1SpK1tl2SjDFbjTH7rLWbpvmevD7fBKqwvP0FyBKcfxcZY5ol3RcpFJ0O5z0FkQC1V4pdbLolrZrmWzjfKTDGbHXwOz0R5/wCRYNUnPsl7UkSZvP6fOdVDZXeqCuZyD/NPriH8z/LIsNPxyfU93DeXWSM8RtjmidcWFr0xqwozrfLIpMsjkxzCOfcZZHXkpi4Htg6cb4TyqtAFUncgUhB3cR9icbl4SLO/+yKFoNG38VHLvx1nHfX1UnarvHvxv2RzwHO96yolrTeGLM9MqOyWQrPrjTGbOScuys6czj+fMa9gWjnfCeWV4EqYpfi1i+JpPCZdiMjuam6fjn/syDyDr5R4bVg6iIvdFsVrmmQOO+usda2Sto9YUhki6TWuIsJ59tF1tqWuHXWdiu8LpIiX0fXRuKcuyTSGzXxd3yrpP1xPVWc7wnybh0qKbZmTPQX5Zq4mQtIUeRCvlHhC0yjpN1KvLou598lkXeOJ5Rgds2EdXo47y6JnPOtcZvqJe1IsFI659tlkQv3FoVfZ3YrvB5YdMkKzrlLEvyOz594Pjnf4+VloAIAAHBTPg75AQAAuIpABQAAkCICFQAAQIoIVAAAACkiUAEAAKSIQAUAAJAiAhUAAECKCFQAAAApIlABmPMi93w7boyxxph90fseRvZtNcY8Edl30BizdcL3Nkf2HZ+4DwCiWCkdQF4wxuyRtDX+djxx+6I33J0Xf/uYuP37rLWbZr+VAHIVPVQA8kXAwTGTbuod6c3K63uUAUiOQAUgXxyXYjd9jYl8fU3ky3H7Iuqste0JtgNADIEKQL7oinye2Au1WdKuRPuMMVuttXtnu2EAch+BCkC+CEQ++6MbjDGNko5Msa9OEj1TABwhUAHIF9Eeqrq4beusta1K3Hu10VrbkpaWAch5BCoA+SIQ+VwtScaYjZLul6S4mX31kX3rJe1Pb/MA5DICFYC8EFdY7o8Wpk9YIiGgN4b8KEQHMCMEKgD5Zr6kzdbaiT1QXZKqKUQHcCEIVADySUDSeoUL0RPtaxSF6AAuQEGmGwAAadQl6UikED3Rvi4K0QFcCHqoAOSTVk296nmrpG1pbAuAOYR7+QEAAKSIHioAAIAUEagAAABSRKACAABIEYEKAAAgRQQqAACAFBGoAAAAUkSgAgAASBGBCgAAIEUEKgAAgBQRqAAAAFJEoAIAAEgRgQoAACBF/z/vAcn4GlayxQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 640x395.55 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"c_obs3.gamma_method(tau_exp=20)\n",
"c_obs3.details()\n",
"c_obs3.plot_tauint()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.10"
}
},
"nbformat": 4,
"nbformat_minor": 4
}