Skip to content

Instantly share code, notes, and snippets.

@nickodell
Last active July 31, 2023 02:16
Show Gist options
  • Save nickodell/490d22cb3077fde89dceb55cd32db38a to your computer and use it in GitHub Desktop.
Save nickodell/490d22cb3077fde89dceb55cd32db38a to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"id": "1ab0e3df",
"metadata": {},
"source": [
"# Quantizing normally distributed floats in Python and NumPy\n",
"https://stackoverflow.com/questions/76798643/quantizing-normally-distributed-floats-in-python-and-numpy"
]
},
{
"cell_type": "markdown",
"id": "eaf966fd",
"metadata": {},
"source": [
"Let the values in the array `A` be sampled from a Gaussian\n",
"distribution. I want to replace every value in `A` with one of `n_R`\n",
"\"representatives\" in `R` so that the total quantization error is\n",
"minimized.\n",
"\n",
"Here is NumPy code that does linear quantization:\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "b2c36039",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from scipy.stats import norm\n",
"import math\n",
"import sklearn.preprocessing\n",
"import jenkspy\n",
"import scipy.optimize"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "07717fe5",
"metadata": {},
"outputs": [],
"source": [
"n_A, n_R = 1_000_000, 256\n",
"mu, sig = 500, 250\n",
"A = np.random.normal(mu, sig, size = n_A)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "5284bf49",
"metadata": {},
"outputs": [],
"source": [
"# Linear method\n",
"lo, hi = np.min(A), np.max(A)\n",
"linear_R = np.linspace(lo, hi, n_R)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "a35a82b4",
"metadata": {},
"outputs": [],
"source": [
"def measure_loss(R):\n",
" lhits = np.clip(np.searchsorted(R, A, 'left'), 0, n_R - 1)\n",
" rhits = np.clip(np.searchsorted(R, A, 'right') - 1, 0, n_R - 1)\n",
" ldiff = R[lhits] - A\n",
" rdiff = A - R[rhits]\n",
" I = lhits\n",
" idx = np.where(rdiff < ldiff)[0]\n",
" I[idx] = rhits[idx]\n",
"\n",
"# plt.plot(A - R[I])\n",
" L = np.mean(np.abs(A - R[I]))\n",
" return L"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "00efe98a",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Linear method: 2.446699788910965\n"
]
}
],
"source": [
"print('Linear method:', measure_loss(linear_R))"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "02da61f1",
"metadata": {},
"outputs": [],
"source": [
"from scipy.stats import norm\n",
"\n",
"dist = norm(loc = mu, scale = sig)\n",
"bounds = dist.cdf([mu - 3*sig, mu + 3*sig])\n",
"pp = np.linspace(*bounds, n_R)\n",
"gaussian_R = dist.ppf(pp)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "cfe549ed",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Gaussian method: 1.6443746592274302\n"
]
}
],
"source": [
"print('Gaussian method:', measure_loss(gaussian_R))"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "94336abe",
"metadata": {},
"outputs": [],
"source": [
"A_samp = np.random.choice(A, size = 10000)\n",
"breaks = jenkspy.jenks_breaks(A_samp, n_classes=n_R-1)\n",
"jenks_R = np.array(breaks)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "62c120c0",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Jenks method: 1.2877975273144864\n"
]
}
],
"source": [
"print('Jenks method:', measure_loss(jenks_R))"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "213f6453",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABLgElEQVR4nO3dd3gU1dfA8e9NDwkhEEIPvRM6oVeliQqKDVQURUEFRQFfQQERG8WOhSZFsWDDX0CQ3msIndBDIAk1AQLpm+x9/5hFIyUsYUuSPZ/n2Se7s7MzZ7LJnp07956rtNYIIYRwXW7ODkAIIYRzSSIQQggXJ4lACCFcnCQCIYRwcZIIhBDCxXk4O4DbVbJkSV25cmVnhyGEEAVKZGRkgtY6+EbPFbhEULlyZbZv3+7sMIQQokBRSp242XPSNCSEEC5OEoEQQrg4SQRCCOHiCtw1ghsxmUzExcWRnp7u7FAKJB8fHypUqICnp6ezQxFCOEGhSARxcXEULVqUypUro5RydjgFitaaxMRE4uLiqFKlirPDEUI4gd2ahpRSs5RS55RS+27yvFJKfaGUOqqU2qOUapLXfaWnpxMUFCRJIA+UUgQFBcnZlBAuzJ7XCOYA3XN5/h6ghuU2EPjmTnYmSSDv5HcnhGuzW9OQ1nqdUqpyLqv0Ar7TRh3sLUqpQKVUWa31aXvFJITIv8xmTUaWmdTMLNJM2aSbsknLNJNmyiYzy0yW2Uy2WWPK1mSb9T+Ps7I1WWZNttmMKVujMZo8AbQGY8nV+9cvw7L+P/dzrJNz2R27Scl/rc2YdBoZ5mQydAqZ5mQyzMlk6QyydaZxw4Rf2nE613+cBxt0tEU0/+HMawTlgdgcj+Msy65LBEqpgRhnDVSsWNEhwd0uf39/kpOTOXXqFK+88gq//fabs0MSwqm01lxOy+JUUhpnktI5lZTG6UvpJKZkkpSWyaVUExdTTSSlZnIpzURqZrazQ7YTM8rzEm5eCbh7JaK8EnDzSsTNMxHlkQJuaShlXaqpujsNClkisJrWejowHaBZs2b5eiadcuXK2T0JZGVl4eFRIN464SIupGSyLz6JI+eSOXruCkfOJnPkXDJJaab/rOemoISfF4FFvAj09aR8oA/1ygVQzNeToj4e+Hi642u5+Xj9e9/Lww0Pd4WHm8LdTeHp7oa72/WP3ZXCzc1o6lQKFEbT59XGT2OZImdraM5l/66nLK+1vulUa01cchx7zu/h0IVDxFyO4cTlE8ReicVk/vf34OvhS0jRECoWbUiQbxABXgEU8y5GgFcAAd4BFPMqRoB3AEXcffE5thKvFe/inXYJr7bDUO1H5O0NugVnfprEAyE5HlewLCvQYmJiuO+++9i3bx9z5swhPDyc1NRUjh07xoMPPsikSZMAWLZsGW+//TYZGRlUq1aN2bNn4+/vz/jx41m4cCFpaWm0bt2aadOmoZSiY8eONGrUiA0bNtC3b1+GDx/u5CMVrux0UhrrDp8nIuYiO05cJDoh5Z/nihfxpEbpotzXoCxVSvpRtpgvZQN9KFvMh2B/bzzcC/7wJa018cnx7Dm/hwMXDnAg8QBRF6K4knkFAE83TyoFVKJKsSp0COlA5YDKhBQNoVJAJYJ9g2+dXJLiYdFwOLwEyjWBp8KhdD27HY8zE0E4MEQp9TPQAkiyxfWBdxbuJ+rU5TsOLqe65QJ4+/68vQm7du1i586deHt7U6tWLV5++WV8fX157733WLFiBX5+fkycOJFPPvmEsWPHMmTIEMaOHQtAv379WLRoEffffz8AmZmZUmdJOIXWmt1xSfy97wxrDp3j4BnjA6+EnxdNKgbycLMKNKoQSK0yRQny93ZytLZn1maiL0UTeTaSyHOR7Di7g7OpZwHjQ79m8Zp0q9yNOiXqEFoylJrFa+LhloePV7MZdsyF5WMh2wRd34eWL4Kbu42P6L/slgiUUj8BHYGSSqk44G3AE0BrPRVYDPQAjgKpwDP2isWZ7r77booVKwZA3bp1OXHiBJcuXSIqKoo2bdoAxgd8q1atAFi9ejWTJk0iNTWVCxcuUK9evX8SwWOPPeacgxAuK+5iKn/ujOePHfFEJ6Tg4aYIq1yCN3vUpmOtUtQo5V9oe51lZmeyMX4jS08sZWP8Ri5lXAIg2DeYJqWb0KRUExqXakz1wOp4uttgMGbiMVg4FGLWQ5X2cP/nUKLqnW/XCvbsNdT3Fs9rYLCt95vXb+724u3977cjd3d3srKy0FrTpUsXfvrpp/+sm56ezksvvcT27dsJCQlh3Lhx/+nf7+fn57C4hWvbcfIi09YeY1nUWbSG5lVKMKhDVbqHlqWYb+EdgW7KNrH59GaWxixl1clVJJuSCfAKoGNIR8LKhNG0VFMqFK1g2+SXnQVbvoLVH4C7N/ScAo37gQMTrFxxdIKWLVsyePBgjh49SvXq1UlJSSE+Pp5SpUoBULJkSZKTk/ntt994+OGHnRytcBVaa9YcPs/Xq48SEXORYr6evNSxGn3CKhJSooizw7ObbHM2W09vZfHxxayKXcWVzCsU9SpK50qd6Va5Gy3KtsDTzU7J78xe+N8QOL0Lat0L934MAWXts69cSCJwguDgYObMmUPfvn3JyMgA4L333qNmzZo8//zzhIaGUqZMGcLCwpwcqXAVe+Iu8eHig2yOTqR8oC9j76vLY2Eh+HkX3o+Iy5mX+Xbvt4QfCychLYGinkXpVLET3Sp3o1XZVrZp7rmZrAxYNxk2fAq+xeGROVD3AYeeBeSk9E0GOeRXzZo109deMD1w4AB16tRxUkSFg/wOXVNSqomJSw/y49aTBPl5MbRzDfo2r4hnIejZczOpplT+OPIHM/fO5GLGRTpU6MD91e6nfYX2eLs74EL3ya0QPgQSDkPDx6Hb+1CkhN13q5SK1Fo3u9FzhTfdCyFy9fe+M4z+cx8XUjIY0LYKr3auQVGfwtv+vz9hP78e/pUlx5eQmpVKs9LN+Drsa+oG1XVMABnJsHI8bJsOxSrAk79D9c6O2fctSCIQwsWkZGQxfmEU87fHUq9cAHOeCSO0fDFnh2UXWeYsVp5cybyoeew6vwtfD1+6Ve7GwzUfpmFwQ8cFcnQFLHwNkmKh+UC4ewx4F3Xc/m9BEoEQLuTouWQGfr+d4wkpvNixGq91romXR+FrBkrKSOKPI3/w48EfOZNyhgr+FXgj7A16Ve9FUS8HfgCnXoClb8HuH6FkTXh2KVRs4bj9W0kSgRAuYkXUWV6dvwtvDzd+eK4FrauVdHZINncm5Qyz981mwdEFpGWl0bxMc95s/ibtK7TH3c6Dsv5Da4j6Exa/DmkXof3r0G4EePo4LobbIIlACBcwc3007/11gPrlizGtX1PKBfo6OySbOpNyhpl7Z/LHkT/QWnNv1XvpV7cftUrUcnwwl0/D4hFwcBGUbQT9FkCZ+o6P4zZIIhCiEDObNRP+Psj0ddH0qF+GTx5thI+nA78Z29m1CeCBGg/wXP3nKO9f3vHBaA07v4eloyE7A7qMh5aDwT3/f8zm/wjFTY0dO5b27dvTuXP+6Hkg8pdss2bk73v4NTKOp1pV4u376+HuVjjKQZxOPm0kgKN/APBg9Qd5rv5zlPMv55yALkQb5SGOr4NKbaHnFxBUzTmx5IEkggJs/Pjxzg5B5FPmHElg6N01eLVzjUJREyjVlMq3+75lzr45mDHTu3pvnqv/HGX9HT8aFwBzNmz5Bla9B+6ecN9n0ORpcCtYF+ALVrT53LvvvkutWrVo27Ytffv25aOPPmLGjBmEhYXRsGFDHnroIVJTUwHo37//f+Yt8Pf3B+D06dO0b9+eRo0aERoayvr168nOzqZ///6EhoZSv359Pv300+u2MX78eMLCwggNDWXgwIH/zNDUsWNH3njjDZo3b07NmjVZv369I38lwgm01ry5YO8/SeC1LjULfBLQWvP38b/p+WdPpu+ZTudKnVn84GLGtBrjvCRwNgq+7QLL3oKqHeClLdDsmQKXBKAwnhEsGWnU77ClMvXhngm5rhIREcHvv//O7t27MZlMNGnShKZNm9K7d2+ef/55AEaPHs23337Lyy+/fNPt/Pjjj3Tr1o233nqL7OxsUlNT2bVrF/Hx8ezbtw+AS5cuXfe63MpXZ2VlsW3bNhYvXsw777zDihUr8vJbEAXExL8P8XNELC/fVZ1XO9dwdjh37EDiASZGTCTybCR1StRhUvtJNCndxHkBZWXA+k9g/cfgEwAPfQuhDzmtPIQtFL5E4CQbN26kV69e+Pj44OPj88+H8L59+xg9ejSXLl0iOTmZbt265bqdsLAwnn32WUwmEw888ACNGjWiatWqREdH8/LLL3PvvffStWvX616XW/nq3r17A9C0aVNiYmJse+AiX5mz8ThT1x7jiRYVGVbAzwQS0hKYsnMKC44soLhPcca2Gkvv6r0d2w30WrERRnmI8wehwWPQ7UPwC3JePDZS+BLBLb65O1r//v35888/adiwIXPmzGHNmjUAeHh4YDabATCbzWRmZgLQvn171q1bx19//UX//v0ZNmwYTz31FLt372bp0qVMnTqVX375hVmzZv2zj1uVr75aCvtqGWxROP297wzvLIqia93SjO8VWmCTQGZ2Jj8c+IFpe6aRkZXBU3WfYlDDQY4dCHZdUCnGdYAt30BAeXj8V6h5/ReygqrgNWblU23atGHhwoWkp6eTnJzMokWLALhy5Qply5bFZDLxww8//LN+5cqViYyMBCA8PByTyZjT9MSJE5QuXZrnn3+e5557jh07dpCQkIDZbOahhx7ivffeY8eOHf/Z99UP/Zzlq4VrOXjmMsN+2UWDCoF80bdxgewdpLVm5cmVPPC/B/gk8hOalW7Ggl4LGBE2wrlJ4Nhq+LolbPkawgbAS5sLVRKAwnhG4CRhYWH07NmTBg0aULp0aerXr0+xYsV49913adGiBcHBwbRo0YIrV4wp/p5//nl69epFw4YN6d69+z+TzqxZs4bJkyfj6emJv78/3333HfHx8TzzzDP/nEF8+OGH/9l3YGCglK92YZdSMxn4XST+3h5M79e0QI4TOHLxCBMjJrL19FaqFavGtM7TaF2+tXODSrsIy0bDznkQVB2eWQKVnByTnUgZahtKTk7G39+f1NRU2rdvz/Tp02nSxIkXtW5DfvkdittjNmuenr2NrdEX+HlQS5pULO7skG5LRnYG03ZPY/a+2RTxLMLgRoN5tNajeZvv15aiwo3RwSkJ0GYodHgj35aHsJaUoXaQgQMHEhUVRXp6Ok8//XSBSQKi4Ppm7THWH0nggwfrF7gkcPjiYV5f+zrRSdH0rNaT15u9TqBPoHODunLWSAAHwqFMA3jiVyjrwCqlTiKJwIZ+/PFHZ4cgXMj2mAt8svww9zcsR9/mIc4O57asPLmSUetH4e/pz9TOU2lTvo1zA9Iadv0IS0eBKR06j4NWQ4xBYi5AEoEQBVBSqolXftpJ+UBfPniw4PQQ0lozfc90vtz1JfVL1uezTp9Rqkgp5wZ1MQYWvgrRq6Fia6M8RMmCP/7idkgiEKIAejt8H+euZPD7i60LzKxiVzKvMG7TOJadWMa9Ve9lXKtx+Hg4sd3dnG3MFrZyPCg3Y+L4ps8WyJHBd0oSgRAFzN/7zvDnrlMMvbsGDUMCnR2OVQ4kHuC1Na9xJuUMrzV9jWfqPePcs5hzB42BYXERUKMr3PsJBBas5jVbkkQgRAGSmJzBWwv2Uq9cAEPuqu7scKyy8NhC3tn8DoHegczpPodGpRo5L5isTNjwKaybbEwV2Xsm1H+4QJeHsAXXOweyo9at89bHuGPHjlzbJVaIG3k7fD9X0rP45NFGeLrn739fU7aJD7Z+wJsb3qRBcAPm3zffuUkgPhKmd4A1H0DdXjAkAho84vJJAOSMwKY2bdrk7BBEIbb60DkW7TnNa51rUqtM/pn4/EYupl/ktTWvEXk2kqfqPsVrTV9z3tiAzFRY/b4xMti/DPT9GWrd45xY8qn8/ZWigLlaSnry5MmEhYXRoEED3n77bQBiYmKoU6cOzz//PPXq1aNr166kpaX95/Vms5n+/fszevTom5aeFq4pLTObMX/uo1qwHy90rOrscHJ19OJR+v7Vl73n9zKh3QReD3vdeUkgei180wo2f2nMEzB4iySBGyh0ZwQTt03k4IWDNt1m7RK1eaP5G1atu2zZMo4cOcK2bdvQWtOzZ0/WrVtHxYoVOXLkCD/99BMzZszg0Ucf5ffff+fJJ58EjFLRTzzxBKGhobz11ltERkbesvS0cB2frzxC3MU05g9sibdH/i0hsT5uPa+vex1fD1/mdJ9D/WAnzdWbdgmWj4Ed30GJqtD/L6jc1jmxFAByRmBjy5YtY9myZTRu3JgmTZpw8OBBjhw5AkCVKlVo1KgRcH1J6EGDBv2TBID/lJ7++++/CQgIcPShiHzi6LlkZq6P5pGmFWhRNX+WPNZa88OBHxiyagghRUP46d6fnJcEDv4FX7UwagS1GQovbpIkcAuF7ozA2m/u9qK1ZtSoUQwaNOg/y2NiYv4pBw1GSeicTUOtW7dm9erVDB8+HB8fH4oXL55r6WnhOt7/KwpfT3feuKe2s0O5oSxzFhO2TWD+ofl0CunEhHYTKOJZxPGBJJ+DJf8H+xdA6VDo+xOUlzIv1pAzAhvr1q0bs2bNIjk5GYD4+HjOnTt3y9cNGDCAHj168Oijj5KVlXXL0tPCNaw5dI7Vh87z8t3VKenvfesXONjlzMsMXjmY+Yfm80zoM3zW6TPHJwGtYddP8FVz42zgrtEwcI0kgdtg1zMCpVR34HPAHZiptZ5wzfMVgblAoGWdkVrrxfaMyZ6UUnTt2pUDBw7QqlUrwLiAPG/ePNzdb92uO2zYMJKSkujXrx8jR47MtfS0KPyyss2899cBKgUV4enWlZ0dznXOppzlhRUvEJMUwzut36F3jd6OD+LSSaM8xLGVENICen4JwTUdH0dBp7W2yw3jg/0YUBXwAnYDda9ZZzrwouV+XSDmVttt2rSpvlZUVNR1yxwtISFBV6xY0dlh5Fl++B2K/5q76biu9MYi/fe+084O5ToxSTG6669ddYsfWugtp7Y4PoDsbK23TNP6vbLGbcs0Y5m4KWC7vsnnqj3PCJoDR7XW0QBKqZ+BXkBUzjwEXL0KWgw4Zcd47ObUqVN07NiRESNGODsUUUgkpZr4ZPlhWlcLomvd0s4O5z8OJB7ghRUvADCr2yzqBtV1bADnD0P4yxC7BardDfd/BoEVHRtDIWPPRFAeiM3xOA5occ0644BlSqmXAT+g8402pJQaCAwEqFgx/73h5cqV4/Dhw84OQxQiX64+wuU0E2Puq5uvKotGno1kyMohFPUqyvQu06lcrLLjdp5tgo2fw9qJ4OUHD04zJpDPR7+fgsrZF4v7AnO01hWAHsD3SqnrYtJaT9daN9NaNwsODr7hhnQBm2ktP5HfXf5yJimduZtP0LtJBeqUzT/dhtfHreeF5S9Q0rck393znWOTwKmdML0TrHoXat8Lg7dBwz6SBGzEnmcE8UDOcn4VLMtyGgB0B9Bab1ZK+QAlgVt3s8nBx8eHxMREgoKC8tW3p4JAa01iYiI+PgV7Gr7CZMqqI2itGXp3/qmJ/3fM34xaN4oaxWvwTedvCPJ10HgGUxqs+RA2TQG/UvDYD1DnPsfs24XYMxFEADWUUlUwEkAf4PFr1jkJ3A3MUUrVAXyA87e7owoVKhAXF8f587f9UoGRSCtUqODsMAQQeyGV+RGx9GkeQkgJJ/TFv4HwY+GM2TiGRsGN+PLuLynq5aA6RzEbjGsBF6KN8hBdxoNvoGP27WLslgi01llKqSHAUoweRLO01vuVUuMxrl6HA8OBGUqp1zAuHPfXeWin8PT0pEqVKrYMXwin+GzFEdzdFC/flT/OBv448gfjNo2jeZnmfHHXF44ZI5CeBMvfhsjZULwyPBUOVTvYf78uzK7jCLQxJmDxNcvG5rgfBTh5slIh8oej55JZsDOOZ9tUoXSA85vq/jjyB29veps25dvwWcfPHDOb2KElsGgYJJ8x5gzu9BZ45Y8zo8Ks0JWYEKKg+nTFYXw83XmxYzVnh8LCYwsZt2kcbcq14fNOn+PtbudRzSkJsOQN2PcblKoHfeZB+ab23af4hyQCIfKBg2cu89ee0wzpVJ0gJ5eSWBS9iNEbR9O8THM+6/SZfZOA1rD3VyMJZFwxzgDavAoeXvbbp7iOJAIh8oGvVx/Dz8ud59o571qX1poZe2cwZecUwsqE8cVdX9i3OSgpDha9BkeWQYUw6DkFStWx3/7ETVmdCJRSRbTWqfYMRghXFJOQwqI9p3i+XVUCizjnm7DJbOLdze+y4OgC7q16L+Nbj8fL3U6xmM0QOQuWjwOdDd0nQPOB4JZ/51ko7G6ZCJRSrYGZgD9QUSnVEBiktX7J3sEJ4QqmrTuGh7sbA9o652zgcuZlRqwZwebTmxnYYCBDGg2x33ichCMQ/gqc3ARVOxnlIYpXts++hNWsOSP4FOgGhANorXcrpdrbNSohXMSZpHR+i4yjT1hFSjmhp9DxpOO8suoV4q7EMb71eB6s8aB9dpRtMgaFrZkAnj7Q62to9LiMDM4nrGoa0lrHXvMNIds+4QjhWmasj8asYWB7x89DvCl+E8PXDsfL3YuZ3WbStLSdeumc3g3/GwJn9kCdntDjIyiavwrpuTprEkGspXlIK6U8gaHAAfuGJUThdyElkx+3nqRXo3IOH0W8KHoRYzaMoWpgVb6860vK+pe1/U5M6UaBuI2fg19JePR7qNvT9vsRd8yaRPACxuQy5TFKRSwDBtszKCFcwZyNx0nPyuYlB48bmBc1j4kREwkrE8bnnT63T8mIE5uM8hCJR6Hxk9D1PfAtbvv9CJu4ZSLQWicATzggFiFcRmpmFt9tOUGXOqWpXspBtXuA+QfnMzFiIp0rdmZC+wm2HyOQfhlWvgMRM405Avr9CdU62XYfwuZumgiUUlMw6v/ckNb6FbtEJIQL+H1HPJdSTQ69NrA4ejHvb32fjhU6MqnDJDzdPG27g8PLjHEBl+Oh5UvG3MFefrbdh7CL3M4ItjssCiFciNmsmbXhOA1DAmlayTHNJRviN/DWhrdoWropkztMtm0SSEmEpaNgz3wIrg0DlkNImO22L+zupolAaz3XkYEI4SpWHTzH8YQUpvRt7JD5M6ISoxi2ZhjVi1dnyl1TbDdaWGvY97tRHiI9CTqMhHbDwMO5JTLE7cutaegzrfWrSqmF3KCJSGstl/+FyIOZG6IpH+jLPaFl7L6v+OR4XlrxEoHegXx999f4e/nbZsNJ8fDXcDi8BMo1gV5fQul6ttm2cLjcmoa+t/z8yBGBCOEK9sUnsSX6Am/2qI2Hu31nik3KSOKF5S+Qac5kVrdZBBe58TSvt8Vshh1zYflYY5BY1/eh5YtSHqKAy61pKNJyt5HW+vOczymlhgJr7RmYEIXRtxuO4+flzmNhFe26n/SsdF5e9TLxyfHM6DqDqoE2uCideAwWDoWY9VClPdz/OZRw/EA4YXvWfCV5+gbL+ts4DiEKvTNJ6SzcfYpHw0Io5mvjHjs5ZJmzeHPDm+w8t5MP23145yOGs7OMQWHftIbTe4wqoU+FSxIoRHK7RtAXY47hKkqp8BxPFQUu2DswIQqbH7aeIFtrnmltv+JyJrOJUetHsfzEcl5v9jrdKne7sw2e2WuUhzi9C2rfZ5SHCLDDKGThVLldI9gEnAZKAh/nWH4F2GPPoIQobDKzzPy0LZa7apWiYpB9ykmkmlIZvnY4G+I3MKzpMJ6q91TeN5aVAesmw4ZPjRHBj8yFur2kSFwhlds1ghPACaCV48IRonBasu80CckZ9GtVyS7bT8pIYvDKwexN2MvYVmN5pOYjed/Yya0QPgQSDkPDx6Hb+1CkhO2CFfmONfMR9AYmAqUAZblprXWAnWMTotD4fvMJKgcVoX0NG/TcucbZlLO8sOIFTlw+wUcdPqJLpS5521BGMqwcD9umQ7EK8OTvUL2zbYMV+ZI1RecmAfdrraXiqBB5EHXqMttPXGT0vXVwc7Nt08rxpOMMWj6IpIwkvun8DS3Ktsjbho6ugIWvQVIstBgEd40BbxuNORD5njWJ4KwkASHy7vstMfh4uvFI0xCbbnd/wn5eXPEiSilmdZ9FvaA8DOhKvQBL34TdP0HJmvDsUqiYx2QiCixrEsF2pdR84E8g4+pCrfUf9gpKiMIiKc3EnztP0atheYoVsV2X0U2nNvHa6tco7lOcaV2mUSngNq89aA1Rf8Li1yHtIrR/HdqNMGYPEy7HmkQQAKQCXXMs04AkAiFu4bfIONJM2Ta9SHx1UpkqgVWY2nkqpYqUur0NXD4Ni0fAwUVQthH0WwBl6tssPlHwWDMfwTOOCESIwsZs1szbcoImFQMJLV/MJtucu38uH23/iLAyYXzW6TMCvG6jz4bWsOM7WDYGsjOgy7tGuWh3q2asFYWYNb2GfIABQD3gn/NGrfWzdoxLiAJvw9EEjiekMPSxRne8La01n0Z+yuz9s+lWuRsftP0AL3cv6zdwIdooD3F8HVRqCz2/gCDHzowm8i9rSkx8D5QBumHUF6qAMahMCJGL7zafIMjPi3vq31mV0ZxJoE+tPkxqP8n6JGDOhk1fwtet4dQuuO8zeHqhJAHxH9acE1bXWj+ilOqltZ6rlPoRWG/vwIQoyE5dSmPVwbO82LEa3h55r8yptWbKzinM3j+bx2o9xpst3rR+DoOzUcbAsPhIqHkP3PsxFCuf51hE4WVNIjBZfl5SSoUCZzAGlwkhbuKX7bFooM8dVhmdunsqM/bO4KEaD1mfBLIyYP3HsP4T8CkGD8+Cer2lPIS4KWuahqYrpYoDY4BwIApjkNktKaW6K6UOKaWOKqVG3mSdR5VSUUqp/ZazDSEKtGyz5tftcbStXpKQEnmrK6S15pvd3/D17q/pVa0XY1uNxU1Z8e8aGwHT2sPaiRDaGwZvg9CHJAmIXFnTa2im5e5awOq6s0opd+AroAsQB0QopcK11lE51qkBjALaaK0vKqXkTEMUeBuOJhB/KY1RPWrn6fVXm4Nm7J1Bz2o9eaf1O7dOApkpsOo92PINBJSHx3+Fml1zf40QFtb0Ghp7o+Va6/G3eGlz4KjWOtqynZ+BXhhnFFc9D3yltb5o2eY5a4IWIj+bH3GSEn5edKlb+rZfq7Xmk8hPmLN/Dg/XfJgxLcfcOgkcWw0LX4FLJyHseej8NngXzWP0whVZc40gJcd9H+A+wJqSE+WB2ByP44Brx67XBFBKbQTcgXFa67+v3ZBSaiAwEKBiRfvO7CTEnUhIzmB51FmeblU5TxeJp+6eypz9c+hTq8+trwmkXYSlo2HXPAiqDs8sgUqt7yB64aqsaRrKORcBSqmPgKU23H8NoCNGt9R1Sqn6WutL18QwHZgO0KxZM22jfQthc3/siMOUrXks7PbrCs3dP5evd3/NA9UfYFSLUbkngahwY3RwSgK0HQYd3pDyECLP8jKksAjGh/atxAM5/xsqWJblFAds1VqbgONKqcMYiSEiD3EJ4VRaa36OiKVppeLUKH17TTO/Hf6Nj7Z/RNdKXRnXatzNm4OunDUSwIFwKNMAnvgVyja0QfTClVlzjWAvRm0hMJpvgoFbXR8A48O8hlKqCkYC6IMx9WVOfwJ9gdlKqZIYTUXRVkUuRD6z/cRFos+nMOnh2xusteT4EsZvHk/b8m2Z0G4C7m43aFLSGnb9CEtHgSkdOo+DVkPA3X5zHwvXYc0ZwX057mdhlKXOutWLtNZZSqkhGM1I7sAsrfV+pdR4YLvWOtzyXFelVBSQDbyutU687aMQIh/4eVss/t4e3NfA+jl918Su4c31b9K0dFM+7fgpnjf6YL8YAwtfhejVULG1UR6iZA1bhS2EVYng2nISATnbLrXWN53IXmu9GFh8zbKxOe5rYJjlJkSBdTndxF97T9G7SQWKeFnX4ropfhPD1gyjdonaTLlrCj4e17Txm7ON2cJWjgflbowMbvosuFkz/EcI61nzF7sDo63/IsY0lYHASctzmtsYWyBEYRW+6xTpJjN9rLxIHHEmgldWv0K1wGpM7TIVf69rZgM7d9AoDxEXATW6wn2fGtNHCmEH1iSC5cACy7d7lFL3AA9orQfZNTIhCpCfI05Sp2wA9a0oN33wwkGGrBxCSNEQpnWZRjHvHK/JyoQNn8K6ycZYgN4zof7DMjJY2JU155gtryYBAK31EkA6KwthsS8+iX3xl+kTFnLLWkBnUs4weMVginoVZVqXaZTwKfHvk/GRML0DrPkA6vaCIRHQ4BFJAsLurDkjOKWUGg3Mszx+Ajhlv5CEKFjmR8Ti5eHGA41yr+yZnJnM4JWDSclK4bt7vvt3ZrHMVFj9Pmz5GvzLQN/5UKu7AyIXwmBNIugLvA0swLgmsM6yTAiXl5aZzZ+74ukRWibXOYlTTakMXjmY6EvRfHn3l9QsXtN4InqtUR7iYgw0e9boFupjm9nMhLCWNSOLLwBDHRCLEAXOkn2nuZKexWO5lJtOz0rn5VUvs+v8Lia1n0Sb8m0g7RIsH2NMHVmiKvT/Cyq3dVzgQuQgk5UKcQd+3hZL5aAitKxa4obPm8wmhq0ZRsSZCD5o9wHdKneDg3/BomGQch7aDIWOo8DT18GRC/EvSQRC5NGx88lsi7nAG91r3/Aisdaa8ZvHsz5+PWNajuG+Us3h1/6wfwGUrg+P/wzlGjs+cCGuIYlAiDz6JSIWdzfFQ01vfJH4691f8+fRP3mhwSAeTcuGL8PAlAp3jTHOBKQ8hMgnrKk1FIwxb0DlnOtrrZ+1X1hC5G+ZWWZ+3xHH3bVLUaro9VU/5+ybw9TdU3mwYhde2rsSoldBSAvo+SUE13RCxELcnDVnBP/DmKx+BUY9ICFc3qqDZ0lIzqRP8+tHEv9y6Bc+jvyY7kVrMHbTz0az0T2TIew5KQ8h8iVrEkERrfUbdo9EiALk54hYygT40L5G8H+Wb4rfxAdb36ed2ZsP96zEo3pnozxEoEyoJPIva76eLFJK9bB7JEIUEPGX0lh7+DyPNquAh/u//0JnLsfyxqpXqJqRwUfnE/B4cBo88ZskAZHvWXNGMBR4UymVCZgsy7TWOsB+YQmRf/263ZiB9ZFm/zYLmeIiGL7seUzKxCcBjSjS9wvwL+WsEIW4LdYMKJNZsIWwyDZrft0eR9vqJQkpUQRMabDmQz45+B17AorycfV+VG4z0tlhCnFbrOo+qpTqCbS3PFyjtV5kv5CEyL82HE0g/lIao3rUhpgNEP4yv5nOMq9kEE/WeISurSUJiILnltcIlFITMJqHoiy3oUqpD+0dmBD50fyIk4T4ZnJPzESYcy+r3LN5r2Qwbcq3YXjLN50dnhB5Ys0ZQQ+gkdbaDKCUmgvsBEbZMzAh8puE5AyyDyzmL9+5uO9MZHOTxxiRtJ26JWrzcYeP8XCT8ZmiYLK2U3NgjvtSGlG4npQELn3/FNM8PsK7aBC7Hv6GoVd2USmgEt90/gY/Tz9nRyhEnlnzFeZDYKdSajXGVJXtAWkIFa5Ba9j7K3rJG1RKu8xPfv1o3GcQL60YRLBvMDO6zvjvDGNCFEDW9Br6SSm1BgizLHpDa33GrlEJkR8kxcGi1+DIMlKCG/PApb483r0Vg9e+ShHPIszoOoOSviWdHaUQd+ymiUApVVtrfVAp1cSyKM7ys5xSqpzWeof9wxPCCcxm2P4trBgH2gzdJzLuRHPOnDvNysTJXM64zPc9vqecfzlnRyqETeR2RjAMGAh8fIPnNHCXXSISwpkSjkD4K3ByE1TtBPd/xmXf8iz6awXV6ixnT8IuJrefTO0StZ0dqRA2c9NEoLUeaLl7j9Y6PedzSqnryy0KUZBlm2DTFFgzATx9oNfX0OhxUIr/bY4h228zJ00reTb0WbpXkfmEReFizcXiTUATK5YJUTCd3g3/GwJn9kDdXkal0KKlAWNymTmRa/ApE07b8m15pfErTg5WCNvL7RpBGaA84KuUaozRYwggACjigNiEsC9TGqydCBu/AL+S8Oj3ULfnf1ZZc+wIZ3ynEeRVmontJ+Lu5u6kYIWwn9zOCLoB/YEKwCc5ll8BZAilKNhObILwlyHxKDR+Erq+B77F/7NKRnYGYza/jnLLZMrdXxDgJXUWReGU2zWCucBcpdRDWuvfHRiTEPaTfhlWvgMRMyGwEvT7E6p1um41rTVjN4wjyRxNE5/XaFi6luNjFcJBrLlGEKqUqnftQq31eDvEI4T9HF5mjAu4HA8tB8Ndb4HXjUcEzzswj8Uxi8g435mhjzzg2DiFcDBrEkFyjvs+wH3AAfuEI4QdpCTC3yNh7y8QXBsGLIeQsJuuvuX0Fj7e/jFFTA0p7daTJhWL33RdIQoDa0YW/2ccgVLqI2Cp3SISwla0hn2/w5L/M5qEOoyEdsPAw/umLzmQeIDha4ZTtkgIBw72ZlCPysacw0IUYnmZSbsIxgXkW1JKdVdKHVJKHVVK3bQ+kVLqIaWUVko1y0M8QlwvKR5+6gu/D4DilWHQOug0KtcksD9hPwOWDcDP04+6bq/h5eZL7yblHRezEE5yyzMCpdRejJHEAO5AMHDL6wNKKXfgK6ALRnmKCKVUuNY66pr1imLMd7D19kIX4gbMZtgxF5aPNQaJdfsAWrwAt+j2uff8XgYtH0SAdwBfdZrOA18coEdoGQKLeDkocCGcx5prBPfluJ8FnNVaZ1nxuubAUa11NIBS6megF8bkNjm9C0wEXrdim0LcXOIxozzEiQ1QpT3c/wWUqHLLl+06t4sXV7xIoHcgs7rNYuOhbK6kZ9GnuUw6L1zDLZuGtNYngCCMD/HeQH0rt10eiM3xOM6y7B+WgnYhWuu/ctuQUmqgUmq7Umr7+fPnrdy9cBnZWbDxc/imNZzZCz2nwFPhViWBHWd3MGj5IEr4lGB299mU9S/Lz9tOUrWkHy2qlHBA8EI4nzVTVY4F5mIkg5LAHKXU6DvdsVLKDWOg2vBbrau1nq61bqa1bhYcHHynuxaFyZm9MPNuoymoemcYvBWaPAVWXOCNOBPBCyteoFSRUszuPpsyfmU4fPYK209cpE/zELlILFyGNU1DTwANrxaes8xhvAt47xaviwdCcjyuYFl2VVEgFFhj+YcrA4QrpXpqrbdbFb1wXVkZsG4ybPjUGBH8yFyjTpCVH95bT29lyMohlPMvx8yuMwkuYnzB+HlbLJ7uioeaWNUfQohCwZpEcApj/MDVCqTe/PcD/WYigBpKqSqW9fsAj199UmudhHGGAYBl8psRkgTELZ3cCuFDIOEwNHwcur0PRaxvxlkXt45ha4YRUjTkP5PLpJuy+WNnHF3rliHI/+a9i4QobHIrOjcFo7dQErBfKbXc8rgLsO1WG9ZaZymlhmCMOXAHZmmt9yulxgPbtdbhtjgA4UIykmHleNg2HYqFwJO/G81Bt2FR9CJGbxhNzeI1mdplKiV8/k0gi/ee5lKqib5ykVi4mNzOCK5+M48EFuRYvsbajWutFwOLr1k29ibrdrR2u8IFHV0BC181po9sMQjuGgPe/re1iR8O/MCEbRMIKxPGF52+wN/rv6//bvMJqpb0o3W1IBsGLkT+d6uic0I4V+oFWPom7P4JStaEZ5dCxRa3tQmtNVP3TOXrXV9zV8hdTOowCW/3/zb97I1LYlfsJcbeVxc3N7lILFxLbk1Dv2itH71mQNk/tNYN7BqZcG1aQ9SfsPh1SLsI7f8P2o/IdWTwjTej+Xr310zdPZWe1XryTut38HC7/s/++y0x+Hq681BTuUgsXE9uTUNDLT/vy2UdIWzv8mlYPAIOLoKyjYxS0WVCb3szOZPAA9Uf4J3W7+Cmru8xfSk1k//tOkXvJhUo5ut55/ELUcDk1jR02lImYo7W+vqC7ULYmtaw4ztYNgayM6DLu9DyJXC3pnPbtZuyLgkA/Lo9jowsM0+1qnSnRyBEgZTrf5jWOlspZVZKFbN09xTCPi5Ew8KhcHwdVG4H938OQdXytKnbSQJms2be1hOEVS5OnbIyA5lwTdbOR7DX0n005epCrbXM4i3unDkbtnwDq94Dd08jATR+CtzyUhjXSAJTdk5hxt4Zt0wCAOuOnOdEYirDu8oMZMJ1WZMI/rDccrru4rEQt+3sfmPe4PhIqHkP3PcJBJTL8+a01ny0/SO+i/qOh2o8xNhWY3NNAgDfbz5BSX9vutcrk+f9ClHQWZMIArXWn+dcoJQaerOVhbilrAxY/7Fx8wmEh2dBvd5Wl4e4EbM288HWD5h/aD59a/dlZPORt0wCsRdSWXXoHEM6VcfLI29nIEIUBtb89T99g2X9bRyHcBWxETCtPaydCKEPweBtxs87SALZ5mzGbhzL/EPzeabeM4xqPuqWSQBg3tYTuCnF4y1kJLFwbbmNI+iLURuoilIqZzmIAOCCvQMThUxminEdYMs3EFAenvgNanS5482azCbeWv8WS2KW8GLDF3mx4YtWVQ1NN2XzS0QsXeqUpmwx3zuOQ4iCLLemoU3AaYzCcDnnLb4C7LFnUKKQObYaFr4Cl05C2PPQ+W3wLnrHm83MzuT/1v0fK0+u5NUmrzKg/gCrX7tgZzwXU00806byHcchREGX2ziCE8AJpVRnIE1rbVZK1QRqA3sdFaAowNIuwtLRsGseBFWHZ5ZApdY22bTJbGLE2hGsjl3NyOYjeaLOE1a/1mzWfLvhOKHlA2guk88IYdU1gnWAj1KqPLAM6AfMsWdQohCICoevWhg1gtoOgxc22iwJmLWZ0RtG5ykJAKw9cp6j55J5rm1VmXxGCKzrNaS01qlKqQHA11rrSUqpXXaOSxRUV84a5SEOhEOZBsa1gLK2K0ultebdLe+y+PhihjYZettJAODb9ccpE+BDj/plbRaXEAWZVYlAKdUKY6ayq42w7vYLSRRIWsOuH2HpKDClQ+dx0OrlPJWHuPkuNB9s/YDfDv/Gc/Wf47n6z932Ng6eucyGown8X/da0mVUCAtr/ktfBUYBCywTy1QFVts1KlGwXIwx5gqIXg0VWxuTx5esbtNdaK2ZFDGJnw/9zNN1n+aVxnkb2P7t+uP4errzuEw+I8Q/bpkItNZrgbU5HkcDUl5CGOUhtk03Zg1T7nDvx9D02TyXh7gZrTWfRH7CvAPzeLLOkwxvNjxPbfvnrqTzv12neCwshMAiXjaNUYiCLLdxBJ9prV9VSi3kxvMR9LRrZCJ/O3fQmDc4LgJqdDPKQxSzfS3/q7WD5uyfQ59affi/sP/L8wXe7zefwGQ2S5dRIa6R2xnB95afHzkiEFFAZGXChk9h3WRjLEDvmVD/4TsaGZybaXumMWPvDB6u+TBvtngzz0ngSrqJuZti6Fq3NFWDb2+KSyEKu9zGEURafq5VSgVb7p93VGAiH4qLNM4CzkVB/Ueg+wTwK2m33c3aN4uvdn1Fr2q9GNNyzB119fxx60kup2fxUkfbXrsQojDItTFXKTVOKZUAHAIOK6XOK6VuOPm8KMQyU2HpW/BtZ0i7BH3nw0Mz7ZoEvo/6nk8jP+WeKvfcspT0raSbspmx/jjtapSkYUig7YIUopC46X+XUmoY0AYI01qX0FoXB1oAbZRSrzkqQOFk0Wvhm1aw+Uto2h8Gb4Fa3e26y/kH5zMpYhJdKnXhg7Yf4O52Z72Vf42MIyE5Q84GhLiJ3K4R9AO6aK0Tri7QWkcrpZ7EGGH8qb2DE06UdgmWjzGmjixRDfr/BZXb2nWXWmtm7p3JFzu/oGOFjkxsN/GGE83fDlO2mWlrj9G4YiAtq0o5CSFuJLf/Ms+cSeAqrfV5pZTM8F2YHVgEfw2HlPPQ5lXoOBI87VuhMz0rnbGbxrLk+BJ6VOnB+Dbj8XS/8z+zhbtPEXcxjXH315NyEkLcRG6JIDOPz4mCKvkcLH4dov6E0vXh8Z+hXGO77/ZsylleWf0KBxIPMLTJUAaEDrDJh7bZrPlmzTFqlynKXbVL2SBSIQqn3BJBQ6XU5RssV4CPneIRzqA17P4Z/h4JpjS4awy0GWrMIWxne87vYejqoaSaUvniri/oGNLRZttevO80R84l83mfRri5ydmAEDeTW/dRqSfkCi6dNMpDHFsJIS2N8hDBNR2y64XHFjJu0ziCiwQzvct0ahSvYbNtZ2Wb+WT5YWqW9ue+BnmfB1kIV2C7imCiYDGbIWImrBhnDAbr8RE0G2Dz8hA3km3O5vOdnzN732zCyoTxcYePKe5T3Kb7+HPXKaLPpzD1ySa4y9mAELmSROCKzh+C8JchditU7wz3fQqBjinCZso2MXL9SJadWMajNR9lZIuReLrZtgkqM8vM5ysPE1o+gG71yth020IURpIIXEm2CTZ+BmsngZcfPDgNGjxmt/IQ10o1pTJszTA2ntrI8KbD6R/a3y77+TUyltgLaYx/JlR6CglhBbsmAqVUd+BzjPkLZmqtJ1zz/DDgOSALOA88a5kiU9jaqZ3wvyFwdh/U6w33TAL/YIftPikjicErB7M3YS/jW4/nwRoP2mU/6aZspqw8StNKxelY03HHJ0RBZrdEoJRyB74CugBxQIRSKlxrHZVjtZ1AM8sMaC8Ck4DH7BWTSzKlwZoPYdMU8C8NfX6E2vc6NITzqecZtGIQMUkxfNThI7pU6mK3ff2w9SRnLqfzyWMN5WxACCvZ84ygOXDUMn8BSqmfgV7AP4lAa51zgpstwJN2jMf1HF8PC1+BC9HQ5GnoMh58Ax0aQuyVWAYuG0hieiJf3f0Vrcq1stu+klJNTFl1hDbVg2hdzX51kIQobOyZCMoDsTkex2HUKrqZAcCSGz2hlBoIDASoWFFmlrql9CRY/jZEzobiVeDphVClvcPDOHLxCIOWDyIjO4OZXWfSINh2cxffyOcrj5CUZuKtHnXtuh8hCpt8cbHYUr+oGdDhRs9rracD0wGaNWt23SQ5IodDS2DRMEg+A61fho5vglcRh4ex89xOhqwcgo+7D3O7z6V6cfsWfIs+n8x3m2PoExZC3XIBdt2XEIWNPRNBPBCS43EFy7L/UEp1Bt4COmitM+wYT+GWkgBL3oB9v0GpetBnHpRv6vAwtNb8cugXJkRMoJxfOaZ1mUaForafuexaHyw+gI+nO8O61LL7voQobOyZCCKAGkqpKhgJoA/weM4VlFKNgWlAd631OTvGUnhpDXt/NZJAxhXo9JZRKM7D8XPyZmRn8O7md/nfsf/Rrnw7Pmz3IcW8i9l9vxuOJLDiwDne6F6b4KLedt+fEIWN3RKB1jpLKTUEWIrRfXSW1nq/Umo8sF1rHQ5MBvyBXy09PE7KXMi3ISkOFr0GR5ZBhTDo+SWUqu2UUE4nn+bVNa8SlRjFCw1f4MWGL97RZDLWyso28+6iKEJK+MpcxELkkV2vEWitFwOLr1k2Nsf9zvbcf6FlNsP2b43yENoM3SdC8+fhDidwyastp7fwf2v/D5PZxJS7pti0cNyt/LD1JIfOXuHrJ5rg4ynlsYTIi3xxsVjchoQjEP4KnNwEVTvB/Z9D8UpOCcWszXy791u+3PUlVQKq8Fmnz6hcrLLD9n/2cjqTlx6iXY2S3BMqpSSEyCtJBAVFtskYFLZmgjFJzAPfQMO+DisPca3LmZd5a8NbrIldwz1V7mFcq3EU8XRs76R3Fu7HlG3mvQeklIQQd0ISQUFwerdRHuLMHqjbC+6ZDEVLOy2cQxcO8dqa1zidfJqRzUfyeO3HHf5BvOrgWRbvPcPr3WpRKcjPofsWorCRRJCfmdJg7UTY+AX4lYRHv4e6zruWrrVmwdEFfLj1QwK8ApjdfTaNSjVyeBxJaSbeWrCPGqX8eb5dVYfvX4jCRhJBfnVik1EqOvEoNO4HXd8FX9vW7L8dKaYUxm8ez+Lji2lRtgUT2k2gpK9zyji8tyiKc1cymPpkU7w87N8zSYjCThJBfpN+GVa+Y0waE1gJ+v0J1To5NaSDFw4yYu0IYq/E8nLjlxkQOgB3J/VQWnngLL9GxjG4UzUahgQ6JQYhChtJBPnJ4WXGuIDL8dByMNz1ljFvgJNorZl/aD6TIyYT6B3It12/pVmZZk6L50JKJqP+2EvtMkV55W7bTWsphKuTRJAfpCQaE8fv/QWC68CA5RAS5tSQrmRe4e1Nb7P8xHLalm/L+23fp4RPCafFo7VmxK+7uZRqYvYzYXh7yJgBIWxFEoEzaQ37focl/2c0CXUcBW2HOaU8RE77EvYxYu0IzqScYVjTYTxd72mHjBLOzayNMaw6eI53etajXjn7l60QwpVIInCWpHj4azgcXmIUh+v5JZR2bvlkrTXzDszjk8hPCPYNZk73OU7pFXStPXGXmLDkAF3qluapVs4ZPCdEYSaJwNHMZtgxx5gvINsE3T6AFi84rTzEVbGXY3l/6/tsPLWRjiEdea/New4pGHcrCckZvPB9JKWK+jD54QYycEwIO5BE4EiJx4zyECc2GBPF3P8FlKji7KhYFL2I8ZvH46bceLPFm/Sp1SdffOCass28NG8HF1Iz+e2F1gQWcW6TmRCFlSQCR8jOgi1fweoPwN3baAZq/KTTykNclWpKZfL2yfx2+DealGrCxPYTKeOXP2r2aK15Z+F+tsVc4PM+jQgt7/yzEyEKK0kE9nZmr1Ee4vQuqH0f9PgIAso6NSStNYuiF/Fp5KecTzvPgNABDGk8BA+3/PPnMG1dNPO2nGRQh6r0alTe2eEIUajln//8wiYrA9ZNhg2fGiOCH5lr1Aly8llA9KVo3tn8DjvO7SA0KJRPO31Kw+CGTo3pWn/ujGfCkoPc37Acb3RzzvwKQrgSSQT2cHKLUR4i4TA0egK6vgdFnNcHHyAzO5OZe2cyY+8M/Dz9GN96PL2q93J6t9BrrT54jtd/202rqkF89EgD3Nycf61CiMJOEoEtZSTDyvGwbToUC4En/4Dqdzs7Kraf2c47m98h5nIM91a9l9ebvU6Qb5Czw7rO2sPnGfR9JLXKFGVqv6YyaEwIB5FEYCtHV8DCV43pI1sMgrvGgLe/U0M6cfkEX+78kr9j/qa8f3mmdp5Km/JtnBrTzaw7fJ7nv9tO9VL+zBvQgmK+ns4OSQiXIYngTqVegKVvwu6foGRNeHYpVGzh1JDOppxl6p6pLDiyAC93L56v/zzP1X/O4RPHWCt89ymG/7KL6qWK8sNzLaSbqBAOJokgr7SGqD9h8euQdhHa/x+0HwEe3k4L6dilY8zdP5dF0YvQaB6t9SgDGwx0Wrloa8zeeJx3FkbRvEoJZjzVTM4EhHACSQR5cfm0UR7i0F9QrrFRKrpMqFNC0Vqz/ex25uyfw7q4dfi4+9C7Rm/61+tPhaIVnBKTNTKzzIxbuJ8ft56kW73SfN6nsUw+L4STSCK4HVrDju9g2RjIzoAu70LLl8Dd8b/GzOxM/o75m3lR8zhw4QAlfErwUqOX6FOrD8V9nDeBjTXOXU5n8I87iIi5yEsdqzG8ay3cpXeQEE4jicBaF6Jh4VA4vg4qt4P7P4egag4P43zqeX45/Au/HPqFC+kXqFqsKmNajqFntZ74ePg4PJ7btTzqLG/8vofUzCym9G3M/Q3LOTskIVyeJIJbMWfDlq9h1fvg7mkkgMZPgZtj+9/vS9jHvAPzWBqzlGxzNu0rtOfxOo/TqmyrfFEX6FaS0kxMWHKQn7adpF65AD7v05jqpZzbq0oIYZBEkJuz+43yEKd2QK0ecO/HEOC4b7Ams4kVJ1Yw78A89pzfg5+nH31q9aFv7b5UDKjosDjuhNaaxXvPMG7hfhKTMxjUvirDutaUMQJC5COSCG4kKwPWf2zcfALh4dlQ70GHlYcwmU0sOLKAaXumcS71HBWLVmRk85H0qtYLf6+C8y1658mLfLj4INtiLhBaPoBZT4dRv4IUjxMiv5FEcK3YCAgfAucPQoM+0P1Dh5SH0FqzP3E/C48t5O+Yv7mQfoHGpRrzdqu3aVu+bb4rBZGbXbGX+HLVUVYcOEtJf2/eeyCUPmEheLgXnGMQwpVIIrgqIxlWvQdbp0JAeXjiN6jRxa67TMpIIuJMBFtOb2Hzqc2cvHISLzcvOoR04IHqD9CufLsC0f4PRnfQZVFn+GHLSTZHJ1LM15NXO9fguXZV8feWPzMh8jP5DwU4tsroEXTpJIQ9D53fBu+iNt9NqimVned2svX0Vrac3sLBCwfRaHw9fGlauinPhj5Ll8pdCPAKsPm+7UFrzYHTV/jfrnh+i4wjMSWT8oG+jLqnNk+0rCQJQIgCwrX/U9MuwtLRsGseBNWAZ/6GSq1stvmrzT3r49az5fQW9iTsIcuchYebBw2DG/JioxdpWbYloSVD8XQrGCNqM7Ky2XHiEisOnGVZ1BliL6Th4aboXKc0fZqH0K5GsIwJEKKAcd1EEBUOi0dASgK0G26UiPC8s374WmtOpZxi7/m97Dy3k9WxqzmdchqFom5QXfrV7UfLMi1pXLoxvh6+NjoQ+zp3JZ2oU5fZcfISW6MT2Rl7icwsM17ubrStUZLBHavTuW5pSvo7r7SGEOLO2DURKKW6A58D7sBMrfWEa573Br4DmgKJwGNa6xh7xsSVs0YCOBAOZRoY1wLKNritTWRmZxKfHE/slViOJx3n2KVjHEs6RvSlaJJNyQB4u3vTsmxLBjcaTIcKHQj0CbTDwdhGuimb2AupnEhM5cSFVE4mpnA8MZUDpy9z/koGAG4K6pUrRr+WlWhRpQStq5eUph8hCgm7/ScrpdyBr4AuQBwQoZQK11pH5VhtAHBRa11dKdUHmAg8ZpeAtIZdPxiVQrMyoPM4aPUyuHtg1maSTckkZSRxOfMylzMucznzMpfSL5GYnkhiWiIJaQkkpidyJuUM51LPodH/bDrIJ4hqgdW4r+p9VA+sTv3g+tQoXsMuzT1aa7LNmiyzxpRtxpStyco2k5ltJivbWJaSmU1KRpZxy8wiOePfxxdSMklIziAx+d+fVzKy/rMPf28PKgUVoX2NYOqVC6BuuQDqlQugqE/BaL4SQtwee36law4c1VpHAyilfgZ6ATkTQS9gnOX+b8CXSimltdbY2JR5T7IsPZK0koFcdAsg6/ACOPwHWmWiVRqom+xSK5T2Q2UXRZmL4pZdCa/sJqjsINyySuBmKk2m9iPKcmBag+YUcIqrR3F1y/8elb7Bc/qG62ptpByzWWPK1pjMZvL621EKihfxoqS/F0F+3tSvEEiQnxfBRb0pH+hLxaAiVCpRhBJ+XgWmt5IQ4s7ZMxGUB2JzPI4Dri3U/886WusspVQSEAQk5FxJKTUQGAhQsWLeRtSq8nfje/wCqZ7lCcQdhRsKN9yUJ57446n88MAPT2W54Yen8sdLBeCm3FEAls9GZbmj/nnMNY/Vv/f/+Ty9/jXqltv798PYw03h4e6Gl7vx09PdDU93hae7Gx6Wn1cfF/Fyx8/LAz/vqzd3/L098PV0lw94IcR1CkQjr9Z6OjAdoFmzZnn6PjzkrmcZwrM2jUsIIQoDew71jAdCcjyuYFl2w3WUUh5AMYyLxkIIIRzEnokgAqihlKqilPIC+gDh16wTDjxtuf8wsMoe1weEEELcnN2ahixt/kOApRjdR2dprfcrpcYD27XW4cC3wPdKqaPABYxkIYQQwoHseo1Aa70YWHzNsrE57qcDj9gzBiGEELmTcpBCCOHiJBEIIYSLk0QghBAuThKBEEK4OFXQemsqpc4DJ/L48pJcM2q5kJPjLbxc6VhBjtcWKmmtg2/0RIFLBHdCKbVda93M2XE4ihxv4eVKxwpyvPYmTUNCCOHiJBEIIYSLc7VEMN3ZATiYHG/h5UrHCnK8duVS1wiEEEJcz9XOCIQQQlxDEoEQQrg4l0kESqnuSqlDSqmjSqmRzo7H1pRSMUqpvUqpXUqp7ZZlJZRSy5VSRyw/izs7zrxSSs1SSp1TSu3LseyGx6cMX1je6z1KqSbOizxvbnK845RS8Zb3eJdSqkeO50ZZjveQUqqbc6LOG6VUiFJqtVIqSim1Xyk11LK8UL6/uRyv895frXWhv2GUwT4GVAW8gN1AXWfHZeNjjAFKXrNsEjDScn8kMNHZcd7B8bUHmgD7bnV8QA9gCcasny2Brc6O30bHOw4YcYN161r+pr2BKpa/dXdnH8NtHGtZoInlflHgsOWYCuX7m8vxOu39dZUzgubAUa11tNY6E/gZ6OXkmByhFzDXcn8u8IDzQrkzWut1GHNW5HSz4+sFfKcNW4BApVRZhwRqIzc53pvpBfystc7QWh8HjmL8zRcIWuvTWusdlvtXgAMY85kXyvc3l+O9Gbu/v66SCMoDsTkex5H7L74g0sAypVSkUmqgZVlprfVpy/0zQGnnhGY3Nzu+wvx+D7E0h8zK0dRXaI5XKVUZaAxsxQXe32uOF5z0/rpKInAFbbXWTYB7gMFKqfY5n9TGOWah7Stc2I/P4hugGtAIOA187NRobEwp5Q/8Dryqtb6c87nC+P7e4Hid9v66SiKIB0JyPK5gWVZoaK3jLT/PAQswTh3PXj1ltvw857wI7eJmx1co32+t9VmtdbbW2gzM4N/mgQJ/vEopT4wPxR+01n9YFhfa9/dGx+vM99dVEkEEUEMpVUUp5YUxN3K4k2OyGaWUn1Kq6NX7QFdgH8YxPm1Z7Wngf86J0G5udnzhwFOW3iUtgaQcTQwF1jXt4A9ivMdgHG8fpZS3UqoKUAPY5uj48koppTDmLz+gtf4kx1OF8v292fE69f119hV0R90wehocxrji/paz47HxsVXF6FWwG9h/9fiAIGAlcARYAZRwdqx3cIw/YZwumzDaSAfc7PgwepN8ZXmv9wLNnB2/jY73e8vx7LF8OJTNsf5bluM9BNzj7Phv81jbYjT77AF2WW49Cuv7m8vxOu39lRITQgjh4lylaUgIIcRNSCIQQggXJ4lACCFcnCQCIYRwcZIIhBDCxUkiEIWaUqqMUupnpdQxS/mNxUqpmkqpyjkre9p4ny8opZ66zdesUUq5zOTsIn/xcHYAQtiLZeDOAmCu1rqPZVlDjJo1sbm99k5orafaa9tC2IOcEYjCrBNgyvnBrLXerbVen3Mly9nBeqXUDsuttWV5WaXUOktt+H1KqXZKKXel1BzL471Kqdeu3amlrvwIy/01SqmJSqltSqnDSql2luW+ljOVA0qpBYBvjtd3VUpttsTyq1LKXylVyVKXv6RSys0Sb1f7/NqEq5EzAlGYhQKRVqx3DuiitU5XStXAGNXbDHgcWKq1fl8p5Q4UwSgIVl5rHQqglAq0YvseWuvmlolG3gY6Ay8CqVrrOkqpBsAOy/ZKAqOBzlrrFKXUG8AwrfV4pdREjMJk24AorfUy634NQuROEoEQ4Al8qZRqBGQDNS3LI4BZlgJhf2qtdymlooGqSqkpwF+ANR/GV4uoRQKVLffbA18AaK33KKX2WJa3xJiIZKPRsoUXsNmy3kyl1CPACxgJSQibkKYhUZjtB5pasd5rwFmgIcaZgBf8MzlMe4xKj3OUUk9prS9a1luD8YE804rtZ1h+ZnPrL18KWK61bmS51dVaDwBQShXBqDwJ4G/FfoWwiiQCUZitArxzTNSDUqrB1Xb6HIoBp7VR/rcfxtSmKKUqAWe11jMwPvCbWJpu3LTWv2M04eR1vtx1GE1PKKVCgQaW5VuANkqp6pbn/JRSV89QJgI/AGMxyhQLYRPSNCQKLa21Vko9CHxmaWtPx5jb+dVrVv0a+N3S5fNvIMWyvCPwulLKBCQDT2HMDDVbKXX1S9SoPIb3jWU7BzCmKoy0xHxeKdUf+Ekp5W1Zd7SlRHEY0EZrna2Uekgp9YzWenYe9y/EP6T6qBBCuDhpGhJCCBcniUAIIVycJAIhhHBxkgiEEMLFSSIQQggXJ4lACCFcnCQCIYRwcf8PmoeQi6Zo0M8AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Plot quantiles of R within distribution against index of R\n",
"\n",
"def quantile_plot(R, name):\n",
" x = np.arange(len(R))\n",
" dist = norm(loc = mu, scale = sig)\n",
" y = dist.cdf(R)\n",
" plt.plot(x, y, label=name)\n",
" plt.xlabel('Class index')\n",
" plt.ylabel('Distribution quantile')\n",
" plt.legend()\n",
"\n",
"lo, hi = np.min(A), np.max(A)\n",
"linear_R = np.linspace(lo, hi, n_R)\n",
"quantile_plot(linear_R, 'linear')\n",
"\n",
"dist = norm(loc = mu, scale = sig)\n",
"bounds = dist.cdf([mu - 3*sig, mu + 3*sig])\n",
"pp = np.linspace(*bounds, n_R)\n",
"gaussian_R = dist.ppf(pp)\n",
"quantile_plot(gaussian_R, 'gaussian')\n",
"\n",
"quantile_plot(np.array(breaks), 'jenks')\n",
"\n",
"# plot synthetic sigmoid curve\n",
"# synthetic_R = sigmoid_classes(7.15067561, n_R)\n",
"# quantile_plot(synthetic_R, 'synthetic')"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "4bc1984b",
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Estimated curve strength of linear: 16.7025\n",
"Average error 0.235003\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAr4klEQVR4nO3deXxU9b3/8ddnsgfCEnYJkIAoOwphE7e6XAEX2qoV3KvV1tZbe229pbWXenttb7XL7+p1K7gU1KsiVKUWxYq4AxIQZBMIGEjYEhIIgeyZ7++PMwwJsgRIcjKT9/PxyGPOfM+ZyefrhLdnzvme7zHnHCIiEvkCfhcgIiINQ4EuIhIlFOgiIlFCgS4iEiUU6CIiUUKBLiISJWL9LkCkoZlZDvA9YCTQ2zn3PX8rEmkaCnSJWs653/ldg0hT0iEXkQZmZjF+1yAtkwJdopaZPWBmL4SW083MmdktZrbVzHab2f21tg2Y2RQz22RmhWY2y8xSa61/1cx2mlmxmX1oZgNrrfurmT1pZvPM7ADwjSbtqEiIAl1amnOBM4GLgalm1j/U/q/AN4ELgNOAPcDjtV73FtAX6AwsB1487H2vB34LpAAfN1LtIsekQJeW5j+dc2XOuZXASmBoqP0HwP3OuTznXAXwAHCNmcUCOOeedc6V1Fo31Mza1nrfN5xznzjngs658ibrjUgtOikqLc3OWsulQOvQci/gNTML1lpfA3Qxs514e9/XAp2Ag9t0BIpDy7mNVrFIPWkPXcSTC4x3zrWr9ZPonNuGdzhlInAJ0BZID73Gar1e05aK7xToIp6ngN+aWS8AM+tkZhND61KACqAQSAY0HFKaJQW6iOcRYC7wjpmVAIuBUaF1M4EtwDZgbWidSLNjusGFiEh00B66iEiUUKCLiEQJBbqISJRQoIuIRAnfLizq2LGjS09P9+vXi4hEpGXLlu12znU60jrfAj09PZ2srCy/fr2ISEQysy1HW6dDLiIiUUKBLiISJRToIiJRolnNtlhVVUVeXh7l5dE5+2hiYiJpaWnExcX5XYqIRKFmFeh5eXmkpKSQnp6OmR3/BRHEOUdhYSF5eXlkZGT4XY6IRKHjHnIxs2fNLN/MVh9lvZnZo2aWbWZfmNmwky2mvLycDh06RF2YA5gZHTp0iNpvHyLiv/ocQ/8rMO4Y68fj3ZqrL3An8OSpFBSNYX5QNPdNRPx33EMuzrkPzSz9GJtMBGY6b9rGxWbWzsy6Oed2NFSRIhJBgkGoLoeq0tBP2aHH6goI1kCwGoJVocfQ85qDz0M/zgHusEe+1hZ0UBMMhn+CwSA1Qec9r/GWXTBI0DkcDufwlp13KDSI93iwHSDoarVxaNuDJdSepdZxqDSHC9/qxB0s9bB2gNRhEzlj2AUN/p++IY6hd6fu7bfyQm1fC3QzuxNvL56ePXs2wK9ueI8++ihPPvkkO3fu5Oc//zlTpkzh9ddf54wzzmDAgAF+lyfS9JyD8r1QvA32bYd9ed5y6W4o2+P9lO45tFx1oEnLC4R+mnKoQdCd2rftpW26QTMN9Hpzzk0DpgFkZmY2y4nYn3jiCd59913S0tLCba+//jpXXHGFAl2i34FC2LECCtZDwZeHHsv31t3OApDcAZJSIak9tE2DbkNwCW0psWSKKgMUVcRSVBlDYWUM+WUx5Jcb+aWO3WWOymCAGgJUEUsNAaqJocZ5zwMxsSTEx5EYH09ifIDkuFiS4mNJSoglMT6W5NgYEuNjiY0JEB8bIC42lrgYIy42hvjYGOJiAsTFxBAfGyA2NvQYiCEmJkDAjJiAERMgvHyozVsOGEds9x7BMDAw74GAWWjZewTqPA+YYQfbQhuMonE0RKBvA3rUep4Waos4P/jBD9i8eTPjx4/ntttuY9OmTVx//fXMnTuXDz74gAcffJA5c+bQp08fv0sVaRjF22DTAti6GHKXQGH2oXVJqdC5Pwz6NqT2gbbdoU0atDmNiqSObCwoZ+32fWzYVUJOYSlbNh9gS1EpldXBOr+ifXIcXdok0rldIl16JNAnJYH2yfG0TY6jfXI87ZLjaJ8cR9ukeNomxREfq8tjTlZDBPpc4G4zexnvfzzFDXH8/D//voa12/edcnG1DTitDb++cuBR1z/11FO8/fbbLFy4kDfffBOAc845h6uuuoorrriCa665pkHrEWlyzsG25bBuLmz8J+Sv8dqTO0DaSDjrBug+HLoMhFYdQy9xbCrYz5Kvili2eg9rd2wiO38l1UHvS3ZiXIBeqa3o3akVF/XrTK8OrUjvkEyP1GQ6pSSQGBfjV29bnOMGupm9BFwIdDSzPODXhA5XOeeeAuYBE4BsoBT4bmMVKyInae9W+OIVWPmytxceiIWeY+DS/4K+l0KnflBrFFbO7gMsXPEVizcXsjRnD0UHKgHo2DqBwd3bcHH/zgzo1pYBp7WhV2oygYBGcDUH9RnlMvk46x3wowarKORYe9IiUk+5S+GT/4Ev/wE46DUWxt4D/a+CpHbhzZxzLN+yh/lrdvLuul1sLvBObPZMTeaifp0ZmZ7KyIxUenVI1vDbZqxZXSnaXKWkpFBSUuJ3GSL14xxkvwsf/Qm2LoLEdnDevTDsFmjfq86mX+0+wGufb+ONFdvYUlhKXIwxuncHbh7di4v7d6FHarI/fZCTokCvh0mTJnHHHXfw6KOPMnv2bJ0UleZr23L451TI+Qja9oRxv4ezb4KE1uFNgkHHe1/m89dPc/g4ezdmcE6fDtz9jdMZN6grKYmaayhSWe0B8k0pMzPTHX6Di3Xr1tG/f39f6mkqLaGP4oOyPfDuf8Ky5yC5I1w4BYbfCjGHwrmiuoZXlubyzMdfsaWwlK5tErlxdE+uGd6Drm0T/atdToiZLXPOZR5pnfbQRSLdur/Dm/d6F/qM/pEX5oltwquraoLMXpbH/y7YyPbicob1bMd9l53JZQO7EhejIYLRRIEuEqkq9sPbU+Dz56HrELjhVTjtrPBq5xzvrN3F7+atY0thKWf1aMfD1wxl7OnROQGeKNBFIlPBBnj5em8I4rn/Bhf+EmLjw6tzi0r59dw1vPdlPmd2SeGZWzK5qF9nBXmUU6CLRJr1b8GcOyA2AW6ZCxnnh1fVBB3TPtzMIws2EDDj/gn9uXVsug6ttBAKdJFIsuhxmP9L6HYWTHrRm0MlJLeolH97ZQVZW/Zw2cAu/PrKgZzWLsm/WqXJKdBFIkEwCO9OhU//FwZMhG/9BeIOhfVbq3Zw3+wvMOB/rjuLb57d3b9axTf6HlbL3r17eeKJJ/wuQ6SuYA3M/VcvzEfcAdc8Fw7z6pog//3WOu56cTmnd27NvHvOU5i3YAr0Wo4W6NXV1T5UI4K3Zz73x7DiBbhgCkz4AwS8ya5Kyqu4bUYWf/lgMzeM6skr3x+tKztbOB1yqWXKlCls2rSJs846i7i4OBITE2nfvj1ffvkl77zzDldccQWrV3u3Vv3jH//I/v37eeCBB9i0aRM/+tGPKCgoIDk5menTp9OvXz+feyMRzzl4855DYf6NX4RX7dpXznefW8r6XSX897cHM3lk87xhjDSt5hvob02Bnasa9j27Dobxvz/q6t///vesXr2aFStW8P7773P55ZezevVqMjIyyMnJOerr7rzzTp566in69u3LkiVL+OEPf8h7773XsLVLy/PuA7B8Jpx/n3exUEjO7gPc8PQS9pRW8swtmVx4Zmf/apRmpfkGejMwcuRIMjIyjrnN/v37+fTTT7n22mvDbRUVFY1dmkS7JX/xZknMvA2+cX94attNBfu5fvpiqmocr9w5hsFpbf2tU5qV5hvox9iTbiqtWrUKL8fGxhIMHroTS3l5OQDBYJB27dqxYsWKpi5PotW6v8NbP4d+V8CEP4bDPDu/hEnTlgCOl+4YzZldU/ytU5odnRSt5VjT5Hbp0oX8/HwKCwupqKgI39GoTZs2ZGRk8OqrrwLe5dYrV65sspolyuxaA3/7PnQfBlc/HT4BmrenlBuf/gyAl+9UmMuRNd89dB906NCBsWPHMmjQIJKSkujSpUt4XVxcHFOnTmXkyJF07969zknPF198kbvuuosHH3yQqqoqJk2axNChQ/3ogkSy0iLvcv6EFLjuxfDQxML9Fdz8zGccqKxm1vfHcHpnhbkcmabPbWItoY9yEoJBePFqyPkYbp0HPUYAUFZZw3XTFrFhVwnP3z6KEempPhcqftP0uSLN3cd/hk3vwRX/Ew7zYNDxs1dXsmpbMdNuylSYy3HpGLqI37YuhoW/g0FXezelCHlkwUb+sWoHvxjfj0sHdDn660VCml2g+3UIqClEc9/kJJXtgdm3Q7se3t55aETL26t38siCjVw7PI07zuvtb40SMZpVoCcmJlJYWBiVweeco7CwkMRE3epLapn377B/J1zzbPguQ1sKD3Df7JUMTWvLg98apDnMpd6a1TH0tLQ08vLyKCgo8LuURpGYmEhaWtrxN5SWYd3fYdUs77L+7sMBKK+q4Uf/t5yAGY9dP4yE2Bifi5RI0qwCPS4u7rhXZopEhQO74e8/8W4dd/7Pws3/PW8dq7ft4+mbMzXRlpywZhXoIi3GvPugYh986+8QEwfAhxsKmLFoC7eNzeASnQSVk9CsjqGLtAgb/wlr/gbn/Qy6DABgb2kl981eSd/Orfn3cWf6XKBEKu2hizSlylL4x73Q8Qw49yfh5qlvrKFwfyXP3DKCxDgdN5eTo0AXaUofPAR7t3pXg8YmAPDOmp3MXbmdey89g0HdNXuinDwdchFpKgUbYNFjcNaNkD4W8O46NPWNNfTrmsJdF/bxuUCJdNpDF2kq79wPcclwyQPhpj+9s4FdJeU8eeMw4mK0fyWnRn9BIk1h47uw8R3v7kOtOwGwIncvMxblcMuYdM7u2d7nAiUa1CvQzWycma03s2wzm3KE9T3NbKGZfW5mX5jZhIYvVSRC1VTD/F9C+wwY9X3Am3hr6hur6ZySwM8u06gWaRjHDXQziwEeB8YDA4DJZjbgsM1+Bcxyzp0NTAKeaOhCRSLWsudg93r4lwfDJ0JnL8vji7xifjG+P60TdORTGkZ99tBHAtnOuc3OuUrgZWDiYds4oE1ouS2wveFKFIlgZXtg4W8h43zodzkA+8qreHj+lwzv1Z6JZ53mc4ESTeoT6N2B3FrP80JttT0A3GhmecA84F+P9EZmdqeZZZlZVrTO1yJSx4d/hPJiuOy/wzMpPvruRgoPVPLAlQM18ZY0qIY6KToZ+KtzLg2YADxvZl97b+fcNOdcpnMus1OnTg30q0WaqX3b4bPpMHQydB0EQG5RKTMW5XDt8DQGp2nMuTSs+gT6NqBHredpobbabgdmATjnFgGJQMeGKFAkYn34B3BBuODn4aY/vbOemIBx76U6ESoNrz6BvhToa2YZZhaPd9Jz7mHbbAUuBjCz/niBrmMq0nLtyYHlM2HYzdC+FwCrtxXz+ort3DY2g65tNS++NLzjBrpzrhq4G5gPrMMbzbLGzH5jZleFNvspcIeZrQReAm510XiXCpH6ev8hCMR6485DHnr7S9olx/H9C3RFqDSOeo2Xcs7NwzvZWbttaq3ltcDYhi1NJEIVbIAvXobRP4Q23QBYtKmQjzbu5leX96dtUpzPBUq00pWiIg3t/d9BbBKc+2/hpkcWbKBzSgI3ju7lY2ES7RToIg1p1xpY8xqMvgtaeeMClmwuZPHmIn5wQR9NjSuNSoEu0pA++jPEt4YxPwo3PbJgIx1bJ3D9qJ4+FiYtgQJdpKEUbvLuRJR5GySnArA0p4hPNxXygwt6a+9cGp0CXaShfPIIBOLq7p2/u5GOreO5YZSOnUvjU6CLNIR922HF/8GwmyClKwDLthTxcfZu7jy/N0nx2juXxqdAF2kInz7mXRV6zo/DTY8uyKZDq3iNbJEmo0AXOVUHCr0pcod8J3xV6PqdJXywoYDvjk0nOV7T40rTUKCLnKolT0FVWZ1x59M/2kxSXIz2zqVJKdBFTkXlAVg63ZvrvJM34daufeW8sWIb38lMo11yvM8FSkuiQBc5FStf8m5icc6hWwDM+DSHmqDjtnMzfCxMWiIFusjJCgZh0RPQfTj0GAXAgYpqXli8hcsGdqVXh1Y+FygtjQJd5GRtnA9Fm7xx56E7D83KymVfeTV3nN/b5+KkJVKgi5ysRY9D2x7Q37vFbnVNkGc/+YrMXu0Z1rO9z8VJS6RAFzkZ21dAzkcw8k6I8YYl/nPtLnKLyvjeedo7F38o0EVOxuInvEm4ht0cbpqxKIe09klcOqCLj4VJS6ZAFzlR+7bD6jlw9k2Q1A7wLiRavLmIG0f3IiZg/tYnLZYCXeREZT0LwRoY9f1w0/OLc4iPDfCdzB7HeKFI41Kgi5yI6kpYNgPOuAxSvXHmJeVVvLZ8G1cOOY3UVrqQSPyjQBc5EevmwoF8GHFHuOlvy7dxoLKGm8foMn/xlwJd5ER8Nh1Se0OfiwBwzvH84i0M7dGOoT3a+VubtHgKdJH62rkKchdD5u0Q8P7pLNpUSHb+fm7WJFzSDCjQRerrs+kQmwRn3xBumrEoh9RW8Vw+pJuPhYl4FOgi9VG2F1a9CoOvgSTvKtBd+8p5d10+12am6X6h0iwo0EXqY8X/QVUpjDx0MvTVrFxqgo7JI3r6WJjIIQp0keMJBmHp05A2EroNDTU5XsnKZUzvDqR31KyK0jwo0EWOZ/NCb1bFWnvnn2zaTW5RGZNHae9cmg8FusjxLH0akjvCgInhppc/y6V9chyXDdS8LdJ8KNBFjqU4Dza8DcNvgdgEAHbvr+CdtTv59rA0EmJ1MlSaDwW6yLF8/gI4B8NuCTfNWZZHVY1j8kjN2yLNS70C3czGmdl6M8s2sylH2eY7ZrbWzNaY2f81bJkiPgjWeIHe+0Jo71045JzjlaW5ZPZqz+mdU/ytT+Qwxw10M4sBHgfGAwOAyWY24LBt+gK/AMY65wYCP2n4UkWa2OaFUJxbZ87zJV8VsXn3ASaN1MlQaX7qs4c+Esh2zm12zlUCLwMTD9vmDuBx59weAOdcfsOWKeKD5TMhuQP0uzzcNCsrl5SEWC4frCtDpfmpT6B3B3JrPc8LtdV2BnCGmX1iZovNbNyR3sjM7jSzLDPLKigoOLmKRZrC/gL4ch4MnRw+Gbq/opq3Vu3kiqHdSIrXyVBpfhrqpGgs0Be4EJgMTDezdodv5Jyb5pzLdM5ldurUqYF+tUgjWPkSBKu8uxKFvLVqB2VVNVwzPM3HwkSOrj6Bvg2ofTo/LdRWWx4w1zlX5Zz7CtiAF/Aikcc573BLj1HQuV+4efayPDI6tmJYz/Y+FidydPUJ9KVAXzPLMLN4YBIw97BtXsfbO8fMOuIdgtnccGWKNKGti6FwY52ToblFpSz5qoirh3XHTPcMlebpuIHunKsG7gbmA+uAWc65NWb2GzO7KrTZfKDQzNYCC4H7nHOFjVW0SKNaPhPiU2Dgt8JNc5bnYQbfGqbDLdJ8xdZnI+fcPGDeYW1Tay074N7Qj0jkKi+GNa/B0EkQ7026FQw65izP45w+HejeLsnnAkWOTleKitS2ajZUl9U53LI0p4jcojKu1t65NHMKdJHals+ALoPhtLPDTbOX5dEqPoZxg7r6WJjI8SnQRQ7avgJ2rPT2zkMnPksrq5m3agcTBncjOb5eRyhFfKNAFzno8+chJgGGXBtuenv1Tg5Uauy5RAYFughAZSl88ao353nSoXHmc5bn0SM1iRHpqT4WJ1I/CnQRgHVzoaK4zsnQbXvL+HRTIVcPSyMQ0Nhzaf4U6CIAy2ZAam9IPzfc9NryPJxDo1skYijQRXZvhK2f1jkZ6pxjzvJtjMpIpUdqss8FitSPAl1k+UywGBh6/aGmrXv4avcBrtbJUIkgCnRp2aorvZkVzxwPKYdu+Dx7WR5JcTFM0LznEkEU6NKybXgbDhTUORlaXlXDmyt3MH5wV1onaOy5RA4FurRsy2dCymnQ5+Jw0/w1OympqOYanQyVCKNAl5Zrby5kvwtn3wgxh/bE5yzfRvd2SYzu3cHH4kROnAJdWq4VL3qPZ98YbtpZXM7HGwv49rDuGnsuEUeBLi1TsAY+fwF6Xwjte4WbX/t8G0GNPZcIpUCXlmnzQijOrXMy1DnH7GW5ZPZqT3rHVj4WJ3JyFOjSMi2fCUmp0O/ycNPKvGI2FRzQRFwSsRTo0vLsL4Av58HQyRCbEG6evSyXhNgAE4Zo7LlEJgW6tDwrX4Jg1dfGnv995Q7GDepKm8Q4H4sTOXkKdGlZnPMOt/QYBZ37hZsXrMunuKxKJ0MloinQpWXZuhgKN9bZOwfvcEu3tomMPb2jT4WJnDoFurQsy2dCfAoM+Ga4ade+cj7Y4I09j9HYc4lgCnRpOcqLYc1rMPhqSGgdbtbYc4kWCnRpOb6YBdVlRxh7nsfwXu3p3an1MV4s0vwp0KVlcA6Wz4Cug+G0YeHmlXnFZOfv51qNPZcooECXlmH757BzFQy7JXxXIvBOhibGaey5RAcFurQMy/4KsUkw5DvhpvKqGuau2M64gRp7LtFBgS7Rr2I/rJ4Dg74NiW3Dzf9cu4t95dVcM7yHj8WJNBwFukS/1XOgcr93uKWW2cvyOK1tImP6aN5ziQ4KdIl+y/4KnfpBj5Hhpp3F5Xy0sYBvD0vT2HOJGgp0iW47V8H25TD81jonQ8NjzzW6RaJIvQLdzMaZ2XozyzazKcfY7mozc2aW2XAlipyCZTMgJgGGXBducs7x6rJcRqS3J0PznksUOW6gm1kM8DgwHhgATDazAUfYLgW4B1jS0EWKnJTKUu9iogETITk13Px57l42a95ziUL12UMfCWQ75zY75yqBl4GJR9juv4CHgPIGrE/k5K19HSqKYfjXT4YmxgWYMFhjzyW61CfQuwO5tZ7nhdrCzGwY0MM5949jvZGZ3WlmWWaWVVBQcMLFipyQZTOgw+nQa2y4yZv3fDvjB3UjRWPPJcqc8klRMwsAfwZ+erxtnXPTnHOZzrnMTp06neqvFjm6/HWQu/hrV4bOW7WDkvJqrs3U4RaJPvUJ9G1A7Ssv0kJtB6UAg4D3zSwHGA3M1YlR8dXymRCIg7Our9P88me5pHdIZkxvjT2X6FOfQF8K9DWzDDOLByYBcw+udM4VO+c6OufSnXPpwGLgKudcVqNULHI8VeXebeb6XQ6tDt2wIjt/P5/lFHHdiJ6Yaey5RJ/jBrpzrhq4G5gPrANmOefWmNlvzOyqxi5Q5IStfR3K9nhjz2t5ZelWYgOm0S0StWLrs5Fzbh4w77C2qUfZ9sJTL0vkFHw23TsZmnFBuKmiuoY5y7dxSf8udEpJ8LE4kcajK0Ulumz/HLZlwYjvQeDQn/c/1+6i6EAlk0ZqIi6JXgp0iS5Ln4a4ZBg6uU7zy5/l0r1dEuf11egqiV4KdIkepUWwarY353lSu3Dz1sJSPs7ezXcye2giLolqCnSJHitehOpyGHFHneaXl24lYPCdEToZKtFNgS7RIRiEpc9AzzHQdVC4uaK6hleW5nJRv850a5vkY4EijU+BLtFh03uw5yvvZGgtb63aSeGBSm4ek+5PXSJNSIEu0WHpdGjVGfrXvTRi5qIcMjq24tzTOx7lhSLRQ4EukW9PDmyY782qGBsfbl69rZjlW/dy4+heBHQyVFoABbpEvqxnwQIw/Lt1mmcuyiEpLkZXhkqLoUCXyFZVBsufh34ToO2hWZ33llbyxortfPPs7rRN0jS50jIo0CWyffEKlBXBqLvqNL+alUdFdZCbx/TyqTCRpqdAl8gVDMKiJ6DbUOh1Tq1mxwtLtjAivT39u7XxsUCRpqVAl8i1aQHsXg9j7q5zE4uF6/PZUljKTRqqKC2MAl0i16LHIOU0GPDNOs3TPtzMaW0TGT+oqz91ifhEgS6Radca2Pw+jLyjzlDFL/L2suSrIm47N4O4GP15S8uiv3iJTIue8GZVPOwmFtM/+oqUhFiuG6FpcqXlUaBL5CnZBatmefcLTU4NN+cWlTJv1Q4mj+pJSqKGKkrLo0CXyLN0OtRUfW2o4nOf5GDAreek+1KWiN8U6BJZyvfBZ9O8G0B3PD3cXFxWxStLt3LFkG6c1k6zKkrLpECXyLLsOSgvhvPurdP80mdbOVBZw/fO6+1TYSL+U6BL5Kgqh08fg97fgO7Dw83lVTU88/FXjD29A4O6t/WxQBF/KdAlcqx4AQ7kw3k/rdP80mdbKSip4O5v9PWpMJHmQYEukaGmCj55BNJGQPq54ebyqhqe+mATIzNSGdOng48FivhPgS6RYfUc2LvV2zuvdZn/rKxcdu2r4CcXa+9cRIEuzV8wCB//P+g8EPpeFm6uqK7hyfc3MSK9vfbORVCgSyRY+zoUfOmNbAkc+pOdlZXHjuJyfnxxX8x0RyIRBbo0bzXVsPB30Kk/DPxWuLmiuoYnF2YzrGc73S9UJESBLs3bqllQuBG+8UsIxISbZy3NZXtxOfdccob2zkVCFOjSfFVXwvu/925g0f/KcPOBimoeWZDNyIxUzu+rvXORgxTo0nyteAH2boGL/qPOyJanP/qK3fsrmDK+n/bORWqpV6Cb2TgzW29m2WY25Qjr7zWztWb2hZktMDPdyFFOTVU5fPAH6DEKTr8k3Lx7fwXTPtzEuIFdGdazvY8FijQ/xw10M4sBHgfGAwOAyWY24LDNPgcynXNDgNnAww1dqLQwWc9CyXa46Fd19s7/d8FGyquD3DfuTB+LE2me6rOHPhLIds5tds5VAi8DE2tv4Jxb6JwrDT1dDKQ1bJnSopTtgQ8fhowLIOP8cPP6nSW8sGQrk0f2oE+n1j4WKNI81SfQuwO5tZ7nhdqO5nbgrSOtMLM7zSzLzLIKCgrqX6W0LB88DGV74bLfhpucc/zmzTW0io/h3ku1dy5yJA16UtTMbgQygT8cab1zbppzLtM5l9mpU6eG/NUSLXZne/OdD7sZug4ON89fs4tPsgv56b+cSWqr+GO8gUjLFVuPbbYBtW/QmBZqq8PMLgHuBy5wzlU0THnS4rzzK4hN8o6dh5RX1fDgP9ZyZpcUbhjV08fiRJq3+uyhLwX6mlmGmcUDk4C5tTcws7OBvwBXOefyG75MaRE2LYQNb8H5P4XWncPNjy7YSN6eMn591QBiYzTSVuRojvuvwzlXDdwNzAfWAbOcc2vM7DdmdlVosz8ArYFXzWyFmc09ytuJHFlNNcz/JbTrVedeoet27GPah5u5Znga5/TRRUQix1KfQy445+YB8w5rm1pr+ZKvvUjkRGQ9C/lr4doZEJcIQE3Q8Yu/raJNUhz3T+jvc4EizZ++v4r/9u2ABb+BPhfBgEMjYmd8msOK3L1MvWIA7XUiVOS4FOjiv7d/DsEquPxP4YuIsvP389DbX/KNMzsx8azTfC5QJDIo0MVfG+bD2jfg/PsgtTcAVTVB7p21guT4GB66eojmaxGpp3odQxdpFGV74c1/g0794Jwfh5sfey+bL/KKeeKGYXRuk+hffSIRRoEu/pl/P5TshOueh1jvGPnizYX873sb+dbZ3ZkwuJvPBYpEFh1yEX+sf9ubHvfcn0D34QAUlFTw45c+J71DK/7rm4P8rU8kAmkPXZregUL4+z3eTZ8v+DngDVH8ySufU1xWxYzbRtI6QX+aIidK/2qkaTkHr98FZUVww6sQmwDAw29/ySfZhTx09WD6d2vjc5EikUmBLk1r8ZOwcT6M/wN0GwLAq1m5/OXDzdw0uhfXjdBcLSInS8fQpelsWw7/nApnXg4j7wBgaU4Rv3xtFWNP78DUKw+/b4qInAgFujSN/QXwyk2Q0hUmPgZmrN9ZwvdmZJHWPpknrh9OnCbeEjklOuQija+mCmbdDKWFcPt8SE4lt6iUm55ZQkJsgJm3jaRtcpzfVYpEPAW6NC7n4K2fw9ZP4epnoNtQdhaXc9MzS6ioDjLr+2PokZrsd5UiUUHfcaVxffIIZD0DY++BwdewbW8Z101bxO79lTz33RGc2TXF7wpFoob20KXxfDEL3v01DLoaLn6A3KJSJk9fTHFpFTNvH8mwnu39rlAkqijQpXFseAde/yGknwfffJJV20v47l+XUlUT5MU7RjEkrZ3fFYpEHR1ykYaX/S68cgN0GQDXvcDC7GKum7aIhNgAc+4aozAXaSTaQ5eGlb0AXroeOp2Ju+l1nli8mz++s54B3drw3K0jNHuiSCNSoEvDWTUbXvsBdOrHvu/M4b45m5m/ZhdXDj2Nh64eTHK8/txEGpP+hUnDWPyUd+ehXmNZOuZx7vnLanaVVPCry/tz+7kZukmFSBNQoMupqa6Et/4dlj1H9RmX83Dr+5g+Yy29UpOZc9c5nNWjnd8VirQYCnQ5eSU74dVbYesiNp15JzdtvpTtJdu5YVRPfjmhP600Ba5Ik9K/ODk5X86DuXcTrDzAUx3v5+GVA+nfLYnHbhqh8eUiPlGgy4kp2+tdLLTsr+Qm9OXW0inkV/XiP644g1vG9CJWE2yJ+EaBLvXjHKx9nao37yOmrJDp1VfwRHASN17QlzvO60275Hi/KxRp8RToclyVW5ZQMvcXdChcxpfBdH7DvQwdfSELLuxDx9YJfpcnIiEKdDmiYNCxLmsh9vGfGLDvY4KuLX9O+AGp532PZ0ekk5Ko6W5FmhsFuoSVV9WweOMOdi2Zw+lbZzHcrWava8VbnW6l9YX38JMBGQQCGk8u0lwp0Fuw8qoavsgrZulXhWzfkEXv7f9gon3AhbaP3bFdWHXGT+k97l8Z30ajVkQigQK9haisDpKdv5+1O/axdvs+1uXlE7M9iwtdFlcGsugZKKAmEENR2sVUnvs9Op5xCR0DMX6XLSInQIEeRYJBx66ScnJ2l7Kl8AA5hd7jV7sPUFyQxxkuh2GBjYyLWccUyyY+popgII7q9Ath4JXEnDmBTq07+d0NETlJ9Qp0MxsHPALEAE87535/2PoEYCYwHCgErnPO5TRsqS1TVU2Q4rIqisuq2FdWxd6yKgpKKsjfV05+SQW7Qo/5+yrYV1JCl+BOetkuelk+GYF8LorP5wy20C6uCABnAeg6BOt1J6SPJZBxPvEJumuQSDQ4bqCbWQzwOHApkAcsNbO5zrm1tTa7HdjjnDvdzCYBDwHXNUbBTc05R9BBdTBIdY2jusZRFVquqglSHXRU1wSpqnHeNkEX2i5IVXhdkPKqIGVVNZRV1niPFdWUV1VSWVlJVUUFVVUVVFVWUFlZTnlZBWUV5QTL9xNbc4BWlNOKclpbmbdsZaRSQv/YErrElNDRimkbLCYp7kDd2uNTsNQM6DIOug2BroOxrkMgsY1P/zVFpDHVZw99JJDtnNsMYGYvAxOB2oE+EXggtDwbeMzMzDnnGrBWAJb+7RE6r54G4Xd2WPiJ92ju0DKA4XChx/Cy+/q6us/rthFqO/g8BojB1Wk72rYHHw0IECSOauKoIWDH+c8TE/o5jMMgORVr1RladYTWA6BVJ++nXU9onwGpGVhyB9AshyItRn0CvTuQW+t5HjDqaNs456rNrBjoAOyuvZGZ3QncCdCzZ8+TKjgupROFyX1Cbxg4+MahALVaAWaYeeHnNYW2MS9aDy6H6vra+3ivDXhviWEBIwBYIEDAwCz0GAgQMCNghpkRCBgBC2AGgdC2gUCA2IARGxMgJjaOmLh4LDYeYuIhJg4CcaHlWO8xEOe1x7eC+NaQ0Np7DC1bXLKCWkS+pklPijrnpgHTADIzM09q7/2sS6+HS69v0LpERKJBfWZS2gb0qPU8LdR2xG3MLBZoi3dyVEREmkh9An0p0NfMMswsHpgEzD1sm7nALaHla4D3GuP4uYiIHN1xD7mEjonfDczHO0X3rHNujZn9Bshyzs0FngGeN7NsoAgv9EVEpAnV6xi6c24eMO+wtqm1lsuBaxu2NBERORG6G4GISJRQoIuIRAkFuohIlFCgi4hECfNrdKGZFQBbTvLlHTnsKtQop/5Gr5bUV1B/G0Iv59wRp0X1LdBPhZllOecy/a6jqai/0asl9RXU38amQy4iIlFCgS4iEiUiNdCn+V1AE1N/o1dL6iuov40qIo+hi4jI10XqHrqIiBxGgS4iEiUiLtDNbJyZrTezbDOb4nc9Dc3McsxslZmtMLOsUFuqmf3TzDaGHtv7XefJMrNnzSzfzFbXajti/8zzaOiz/sLMhvlX+ck5Sn8fMLNtoc94hZlNqLXuF6H+rjezy/yp+uSYWQ8zW2hma81sjZndE2qPys/3GP317/N1zkXMD970vZuA3kA8sBIY4HddDdzHHKDjYW0PA1NCy1OAh/yu8xT6dz4wDFh9vP4BE4C38G7ROhpY4nf9DdTfB4CfHWHbAaG/6QQgI/S3HuN3H06gr92AYaHlFGBDqE9R+fkeo7++fb6RtocevmG1c64SOHjD6mg3EZgRWp4BfNO/Uk6Nc+5DvDnzazta/yYCM51nMdDOzLo1SaEN5Cj9PZqJwMvOuQrn3FdANt7ffERwzu1wzi0PLZcA6/DuNxyVn+8x+ns0jf75RlqgH+mG1cf6DxiJHPCOmS0L3VQboItzbkdoeSfQxZ/SGs3R+hfNn/fdocMMz9Y6hBY1/TWzdOBsYAkt4PM9rL/g0+cbaYHeEpzrnBsGjAd+ZGbn117pvO9uUTvWNNr7F/Ik0Ac4C9gB/MnXahqYmbUG5gA/cc7tq70uGj/fI/TXt8830gK9PjesjmjOuW2hx3zgNbyvZLsOfhUNPeb7V2GjOFr/ovLzds7tcs7VOOeCwHQOfe2O+P6aWRxeuL3onPtbqDlqP98j9dfPzzfSAr0+N6yOWGbWysxSDi4D/wKspu5NuG8B3vCnwkZztP7NBW4OjYYYDRTX+uoesQ47TvwtvM8YvP5OMrMEM8sA+gKfNXV9J8vMDO/+wuucc3+utSoqP9+j9dfXz9fvM8UncWZ5At7Z5E3A/X7X08B96413FnwlsOZg/4AOwAJgI/AukOp3rafQx5fwvoZW4R1DvP1o/cMb/fB46LNeBWT6XX8D9ff5UH++CP0j71Zr+/tD/V0PjPe7/hPs67l4h1O+AFaEfiZE6+d7jP769vnq0n8RkSgRaYdcRETkKBToIiJRQoEuIhIlFOgiIlFCgS4iEiUU6CIiUUKBLiISJf4/g/nPUTE5KlwAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Estimated curve strength of gaussian: 0.0005\n",
"Average error 0.000154\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAoL0lEQVR4nO3dd3hUZf7+8fcnIRBAOkgLgQChNzEColIUkI64KN3d1d9iX+ziWhbbWrGwioiuunRFFFFQUVRARSBIEUILASTUEHoLKc/vj4x+IwskwCQnM3O/rouLmTlnZu6HE29PnplzjjnnEBGRwBfmdQAREfEPFbqISJBQoYuIBAkVuohIkFChi4gECRW6iEiQUKGLnAMzG2xmc7zOIZKT6XvoIiLBQXvoIiJBQoUuAcXMWprZMjM7ZGbTzOx9M3vKzMqZ2WdmlmJm+3y3o3I8b7OZdcpxf6SZTfTdjjSziWaWamb7zWyJmVX2LfuLmSX53m+TmQ3O8fj3OV7vVTPbamYHzWypmV1x0nt9YGbjfa+z2sziCuLfS0KLCl0ChpkVBT4G3gPKA1OAvr7FYcC7QE0gGjgGvJbHl/4zUAaoAVQAbgGOmVlJYDTQzTlXCmgLLD/NaywBWvhyTQammVlkjuW9galAWWDmWWQTyTMVugSSNkARYLRzLt059xGwGMA5l+qcm+6cO+qcOwQ8DbTP4+umk13kdZ1zmc65pc65g75lWUATMyvunNvhnFt9qhdwzk30Zchwzo0CigH1c6zyvXNutnMuE5gAND/LsYvkSoUugaQasM398ZP8rQBmVsLM3jSzLWZ2EJgPlDWz8Dy87gTgS2CqmW03s+fNLMI5dwToT/Ye+w4zm2VmDU71AmZ2n5mtMbMDZraf7D3+ijlW2Znj9lEg0syK5G3YInmjQpdAsgOobmaW47Eavr/vJXuPuLVzrjTQzvf4b+seAUrkeF6V32749vYfd841IntapSdwg2/Zl865zkBVYC3w1smhfPPlDwDXA+Wcc2WBAzneW6RAqNAlkCwEMoE7zKyImfUBWvmWlSJ73ny/mZUH/nnSc5cDA8wswveBZL/fFphZRzNr6tubP0j2FEyWmVU2sz6+ufQ04DDZUzAnKwVkAClAETN7DCjtnyGL5J0KXQKGc+4EcC1wE7AfGAJ8RnbZvgIUB/YAPwFfnPT0R4E6wD7gcbI/uPxNFeBDsst8DTCP7GmYMOAeYDuwl+w5+VtPEe1L3/utB7YAx/FNBYkUJB1YJAHNzBYBY51z73qdRcRr2kOXgGJm7c2sim/K5c9AM/53b1wkJOlTdgk09YEPgJJAEtDPObfD20gihYOmXEREgoSmXEREgoRnUy4VK1Z0tWrV8urtRUQC0tKlS/c45yqdaplnhV6rVi3i4+O9ensRkYBkZltOt0xTLiIiQUKFLiISJFToIiJBolB9Dz09PZ3k5GSOHz/udZR8ERkZSVRUFBEREV5HEZEgVKgKPTk5mVKlSlGrVi3+eEK9wOecIzU1leTkZGJiYryOIyJBKNcpFzN7x8x2m9mq0yw3MxttZolmttLMWp5rmOPHj1OhQoWgK3MAM6NChQpB+9uHiHgvL3Po7wFdz7C8GxDr+zMMeON8AgVjmf8mmMcmIt7LtdCdc/PJPnXo6fQBxrtsP5F9lZiq/gooIhIsjp3IZMrU99i2e0++vL4/vuVSnT+e+znZ99j/MLNhZhZvZvEpKSl+eGv/Gz16NA0bNqRcuXI8++yzAMyYMYOEhASPk4lIIFu8ej3zn+vLwLXD2fX16Hx5jwL9UNQ5Nw4YBxAXF1cozwo2ZswYvv76a6Kion5/bMaMGfTs2ZNGjRp5mExEAtHBYyf4bMrrdNnyEmXtCMnN7qRl74fz5b38Uejb+L/rOgJE+R4LOLfccgtJSUl069aNG2+8kY0bNzJo0CBmzpzJvHnzeOqpp5g+fTp16tTxOqqIBID58Stg1j0McvFsv6ARmYPGEVW9ab69nz8KfSbZ13icCrQGDvjj/NSPf7qahO0HzztcTo2qleafvRqfdvnYsWP54osv+Pbbb/nss88AaNu2Lb1796Znz57069fvtM8VEfnNnkPH+GrC8/TY9QZFLZPtrR+h2tX3QFh4vr5vroVuZlOADkBFM0sm++K7EQDOubHAbKA7kAgcBf6aX2FFRAoz5xxff7+QsnPvYyCr2Vr2EioPfpNqFxbMb/W5FrpzbmAuyx1wu98S+ZxpT1pEpLDZvvcQ3094nN573yMzrCi72r9AjfZ/gwL8unKhOlK0sCpVqhSHDh3yOoaIFEJZWY7ZX39FzA8Pcr0lseXCDkQNeYOSZaoVeBadnCsPBgwYwAsvvMBFF13Exo0bvY4jIoXEpl17+filW7n6hwFEhe9lT7dx1LxtBuEelDl4eE3RuLg4d/IFLtasWUPDhg09yVNQQmGMIsEuIzOLT2fNoNnSh6lj29kU1YdaA1/GSlbI9/c2s6XOubhTLdOUi4jIWVizZQeJU+6nz7HP2BdRiX29pxLTrJvXsQAVuohInqRlZPLp9Am0SXiSHpbKlrqDqXXdM1hkaa+j/U6FLiKSi+XrN7Fr2j30S/+GXcWiOdzvv8TUu9zrWP9DhS4ichpH09KZ9f5YOm58niZ2hM2Nb6PWNf+EiEivo52SCl1E5BQWr1jN8U/u5rqsRWwrUZ8TA8dRK7qF17HOSIUuIpLDgSMnmDPpBa7e9hrFLINfLx5BdPf7Ibzw16W+h57D/v37GTNmjNcxRMQj835azPoXr+S67c+zv3QDuPUHons9FBBlDir0PzhdoWdkZHiQRkQKSsqBo0x/bQStPu9BI7eRbZc/Q/TdcylWuZ7X0c5KYPxvp4CMGDGCjRs30qJFCyIiIoiMjKRcuXKsXbuWOXPm0LNnT1atyr606osvvsjhw4cZOXIkGzdu5PbbbyclJYUSJUrw1ltv0aBBA49HIyK5cc7x9XffUWXeffyJRDZVuIKooW9QvVyN3J9cCBXeQv98BOz8xb+vWaUpdHv2tIufffZZVq1axfLly/nuu+/o0aMHq1atIiYmhs2bN5/2ecOGDWPs2LHExsayaNEibrvtNr755hv/ZhcRv9q2Zz+Lxz9CjwOTORZWkp2dxhDTdlCBnkzL3wpvoRcCrVq1IiYm5ozrHD58mB9//JHrrrvu98fS0tLyO5qInKOsLMcXX35KvZ8eoq8ls7FaD2IGj6bMBRW9jnbeCm+hn2FPuqCULFny99tFihQhKyvr9/vHjx8HICsri7Jly7J8+fKCjiciZylp+24SJt5P9yOfsLdIRVJ6TKROy15ex/IbfSiaw5lOk1u5cmV2795NamoqaWlpv1/RqHTp0sTExDBt2jQge05uxYoVBZZZRHKXnpnFzI8mEfFmW3oenUFSTH8q3L+USkFU5lCY99A9UKFCBS677DKaNGlC8eLFqVy58u/LIiIieOyxx2jVqhXVq1f/w4eekyZN4tZbb+Wpp54iPT2dAQMG0Lx5cy+GICInSUjaQvLUe+l94it2FY1i37WfULdhB69j5QudPreAhcIYRQqD4+mZfDFtHG3XPUt5O8iW+v+POv2eLLSH7eeVTp8rIiFlecI6Dnw0nGsyFpJcPJZj/adTJ+aUHRhUVOgiEjQOH0/nq8kvceWWV4i0dJJa3E/tXg9CeITX0QpEoSt05xwWwN8DPROvprdEQsGin3/GPruLvlkr2FKqOZUGvUntaqE1vVmoCj0yMpLU1FQqVKgQdKXunCM1NZXIyMCevxMpbPYfPsZ3E56my85xYMaWNk9Ss8sdEBZ6X+IrVIUeFRVFcnIyKSkpXkfJF5GRkURFRXkdQyRoLPhhAWW+uodrWE9SucuoNuQNalas6XUszxSqQo+IiMj1yEwRkd37DvLT+EfpuncCx8JKktxxNLWvuCGgD9v3h0JV6CIiZ+KcY+7cz6n5/QP0ZisbKnclZsi/KVP6Qq+jFQoqdBEJCMm79rBy4oNcfXA6+8MrsLPru8S2utbrWIWKCl1ECrXMLMdXsz6gUfyjdLddrK/Rj7qDRhFWoqzX0QodFbqIFFpJW5PZOPleuh77gp0R1Ui55iPqNbnK61iFlgpdRAqd9Mws5kz/D3Grn+ZK28+6ujdS7/qnsaIlvI5WqKnQRaRQWbM+kd3T7qJH+gKSi9Xm4HXvUz+2tdexAoIKXUQKheMnMvh66mgu2ziKOnacDU3uIrbvIyFz2L4/5KnQzawr8CoQDrztnHv2pOXRwH+Bsr51RjjnZvs3qogEq2UrV5L+yd/pmbmMzSWbUGTQW8RGNfI6VsDJtdDNLBx4HegMJANLzGymcy4hx2qPAB84594ws0bAbKBWPuQVkSBy6Fga8yY+Q8fkMZgZGy/5J3W63RWSh+37Q1720FsBic65JAAzmwr0AXIWugNK+26XAbb7M6SIBJ9FSxYSOfsuerq1JJZpTfWhb1Knko4UPx95KfTqwNYc95OBkz+hGAnMMbM7gZJAp1O9kJkNA4YBREdHn21WEQkC+w4e4ccJj9Fp93ukhRVnc7uXqNvxxpA/bN8f/PV7zUDgPedcFNAdmGBm//Pazrlxzrk451xcpUqV/PTWIhIInHMsmPcVu19qS4+Ut9lcsQPFhi+h1pU3qcz9JC976NuAGjnuR/key+kmoCuAc26hmUUCFYHd/ggpIoFtV+o+lo1/kM77P2B/eDm2dnqL+m2v9zpW0MnLHvoSINbMYsysKDAAmHnSOr8CVwGYWUMgEgjOc+CKSJ455/jmi+kc/3cbuh54n/XVr6HMvT9TQ2WeL3LdQ3fOZZjZHcCXZH8l8R3n3GozewKId87NBO4F3jKzu8n+gPQvTpfnEQlpW7fvZN3Eu+l0dDY7w6uys9cHNGxxtdexgpp51btxcXEuPj7ek/cWkfyTmeWY+8l/abb8cSqxj3W1h9JgwLOEFSvpdbSgYGZLnXOnvOK1jhQVEb/ZuHkz26b8nS5p89haNIZ9/SbRqH5br2OFDBW6iJy3E+mZfDPtdVqte45oO8aaBnfSoN+jWJFiXkcLKSp0ETkvCWsSODz9TrpmxLOpeCPCBoylYa3mXscKSSp0ETknx9LSmTf5OS7b/BpFLIt1LR6mfu97ISzc62ghS4UuImdt2c9LCPvs73TNSmBDqUuoMmQs9avU9TpWyFOhi0ieHTxylB8njKTjjnc4YUXZ0PY5YjvfrCM9CwkVuojkyaIfv6XsV3fT1W1iTfmO1Bo6htjy1byOJTmo0EXkjFL3H2Dp+Ie4MnUKB8LKsOmqsTS8YqDXseQUVOgickrOOb6fO5Ma34+gC9tZXaU3sUNepUKp8l5Hk9NQoYvI/9i5ezcJE+7lykMz2RVWma3dp9A4rrvXsSQXKnQR+V1WlmP+rEnUj3+MDuzll+jBNBr8POGRF3gdTfJAhS4iAPy69Ve2TB5Oh2PfkBwRza6+79G0cTuvY8lZUKGLhLiMjEzmffwmLVY9QxuOsCr2Vhr3H4lFRHodTc6SCl0khCUmrmPvB3dw1YnFbCpWH9d/LE3qtPQ6lpwjFbpICEpLT2f+1FG0SXyFKMtkdZMHadT3ASxclRDItPVEQszqVctI//hOOmf+wvqSLblw0FgaR9X3Opb4gQpdJEQcPX6cHyc+weVbx5FuEay55Gkadr9dh+0HERW6SAhYvmQBkZ8Pp1PWRhLKXkH0kDE0rBTtdSzxMxW6SBA7cOgw8eP/QbvdEzkUdgHr279Gow5DtFcepFToIkFq0bzZXPjtfVzFNlZW6k69oaOpV6aS17EkH6nQRYLMnr2prBp/L+32zSAlrCJJV4+nWZs+XseSAqBCFwkSzjl+nPM+MQsfoZ3bw8rq/Wk89AUqFy/tdTQpICp0kSCwY8c2kiYO57IjX5EcXoNtvT+mRfOOXseSAqZCFwlgWZlZLPj0PzRZ9iStOMzy2n+j6cAnCS9a3Oto4gEVukiA2rI5kV1T7qR92o8kFY0lvd90WtS/xOtY4iEVukiAycjI5PtpL9Ny7SgqWzorGt5Ls34PYeERXkcTj6nQRQLIhrUrOfrh7XTIWMm64s2pOPBNmtds6HUsKSRU6CIBIO3ECRZOforWm94g04rwy0WP06TXnVhYuNfRpBBRoYsUcquXLyRs5t/pkLWeVaXaEjXkDZpWqeV1LCmEVOgihdSRI0dYMuER2u74L0esJKvbvkKTzn/RYftyWmF5WcnMuprZOjNLNLMRp1nnejNLMLPVZjbZvzFFQsvyhXPY/WJrOux8hzXlOxExPJ7GXf6qMpczynUP3czCgdeBzkAysMTMZjrnEnKsEws8BFzmnNtnZhfmV2CRYHZg/35WjL+Py1M/JCWsAus7vUPzy//kdSwJEHmZcmkFJDrnkgDMbCrQB0jIsc7fgNedc/sAnHO7/R1UJNgtmfsR1RY8SDt283OVfjQaOorKF5T1OpYEkLwUenVga477yUDrk9apB2BmPwDhwEjn3Bcnv5CZDQOGAURH61zMIgApKTvZMGE4bQ9+QXJYdZK6T6NlXBevY0kA8teHokWAWKADEAXMN7Omzrn9OVdyzo0DxgHExcU5P723SEByzrFw1nvExo+klTvI0ui/0mzIv4goVsLraBKg8lLo24AaOe5H+R7LKRlY5JxLBzaZ2XqyC36JX1KKBJntyZtJnnQHbY8tIKlIHY71fZ+LG7fxOpYEuLx8y2UJEGtmMWZWFBgAzDxpnRlk751jZhXJnoJJ8l9MkeCQlZnF99NepeRbbWl+9CeWxg6n1oM/Ea0yFz/IdQ/dOZdhZncAX5I9P/6Oc261mT0BxDvnZvqWdTGzBCATuN85l5qfwUUCzZaNCex//zYuP7GMtZFNKdt/DBfXbuZ1LAki5pw3U9lxcXEuPj7ek/cWKUjp6en8NPUZLk58jSwLY33T+7io7906bF/OiZktdc7FnWqZjhQVyUcbVi0h4+PbuSJzHatKtqbK4DdoWb2O17EkSKnQRfLB8ePHWDLxUVpvfYejVoIVrV6gebe/6UhPyVcqdBE/S4j/lsjZw7kiawvLynaiztDXaF6xqtexJASo0EX85PDhg6wY/wBtdk0l1cqzuv04LurY3+tYEkJU6CJ+sGL+TCp+ex+XuV0sqdSXRkNf4sIy5b2OJSFGhS5yHg7sTWHN+OG02T+LZKvKum5TuKR1d69jSYhSoYucA+ccS+dMoubCR4hzB1hcfSjNhjxLVIkLvI4mIUyFLnKWUnb8ypaJdxJ35DuSwmM42HsCrZpf4XUsERW6SF65rCwWzxxDg+X/oplLY1Ht27l44D8pUrSY19FEABW6SJ5s37yOPVNuo3VaPGsjGlHyujdoXa+F17FE/kCFLnIGmZmZLPngeZqufZkyGIsbPURcv/sJC9dh+1L4qNBFTmPLumUc/fA22qQnsLL4JVw4aAytout5HUvktFToIidJP5FG/KR/cvHmtzhqkcS3fIaLe96CheXpmuoinlGhi+SwYfkCwmfewaVZm1laqgO1hr5GXOUauT9RpBBQoYsAx48eZtmEB2m1fRKpVo5lbV/n4i5DvI4lclZU6BLyEn6cTemv7uFSt4NF5XvR4IZXuKhcRa9jiZw1FbqErEP7U1kz4W5apX7CNqvMqk4TaH15b69jiZwzFbqEpJVzp1BlwT+42O1jYZWBNL/heaqXLO11LJHzokKXkLIvZTtJE+7g4oNzSQqryd4e73LpxR28jiXiFyp0CQkuK4tls9+idvyTNHVH+aHmzcQNfpxixYp7HU3Eb1ToEvRStm1kx8RbaXlsEWuLNKDota9zWaNTXmNXJKCp0CVouaxMlk4fRcPVo6jrHD/Wu59W1z9IkYgIr6OJ5AsVugSl7YkrOfDBrcSdWMWKoi2pMGAMbWs39DqWSL5SoUtQyUw/wdKpT9A8cSwlKcrCZk/S+po7CAvXYfsS/FToEjQ2r1pI5ozbaZWxkSUl2xE95DUurVbT61giBUaFLgHvxLEjrJj0EBdtncA+K83iVq9ySbc/Y2ZeRxMpUCp0CWiJS+ZQ7PO7uCRrGz+V7U69oa/QqmJlr2OJeEKFLgHp2KF9rJ5wD3G7P2Iblfm5/bu06Xit17FEPKVCl4CTMH8aFb4dQcusVL6v1J9mNzxP9dJlvY4l4jkVugSMQ3t3kjj+Ti7aP4dNVoNV3aZxeZvOXscSKTRU6FL4OccvX75D1E8jaeyOsKD6TcQNeYqYEiW8TiZSqOTpy7lm1tXM1plZopmNOMN6fzIzZ2Y6rlr8Yt+OTfzyYnea/nQPu8KrsOnaz7li2EsUV5mL/I9c99DNLBx4HegMJANLzGymcy7hpPVKAcOBRfkRVEKLy8pkxczR1Fn+HHVdJvNr302bgQ9TtKgO2xc5nbxMubQCEp1zSQBmNhXoAySctN6TwHPA/X5NKCFn9+bV7J1yCy3SVrIyojmlrh9Du9gmXscSKfTyMuVSHdia436y77HfmVlLoIZzbtaZXsjMhplZvJnFp6SknHVYCW5ZGen8PGUkpd9tT7Xjicxv+BiNR3xHjMpcJE/O+0NRMwsDXgL+ktu6zrlxwDiAuLg4d77vLcFj29rFHJ9+Gy3TN7CkeFuqDXyddjVrex1LJKDkpdC3ATVy3I/yPfabUkAT4DvfodZVgJlm1ts5F++voBKcMtKOsXLKIzTd9C4H7QJ+aDmKtj1vxMJ0Mi2Rs5WXQl8CxJpZDNlFPgAY9NtC59wB4PdLpJvZd8B9KnPJzeZl3xD+2Z20zEzmx1JdqDt0NJdVrup1LJGAlWuhO+cyzOwO4EsgHHjHObfazJ4A4p1zM/M7pASXtKMHSJhwP823f8BOq8iiy97i0k7X6WRaIucpT3PozrnZwOyTHnvsNOt2OP9YEqw2/DiD0l/dR/OsPXxf/lqa3vAircuV9zqWSFDQkaJSII7u382G8X+n+d7P2WzVWdZ5Ku0u7+p1LJGgokKX/OUca7+ZwIULHqGRO8y3Vf5C3NCnqXXBBV4nEwk6KnTJNwd3/8qvE26jyaEFrAurQ3KPyXS8+HKvY4kELRW6+J9zrP7sNaKX/ou6Lp1vat5J28GPElmsmNfJRIKaCl38KnXrWvZMvoXGx5axokhTIvu+zpWNm3sdSyQkqNDFL1xmBr989Cyxq1+lmgvnm3oPc0X/e4gooh8xkYKi/9rkvO1K/JnDH9xKsxNrWVK0NRUHvMaVtet5HUsk5KjQ5ZxlnTjOL+8/RqPEtylCSb5r9hxXXDOM8HAdti/iBRW6nJNtv8wj65M7aJ7xK9+XvIqYwa/SoXqN3J8oIvlGhS5nJePYIVZPeoCmW6ewy8qzoNUbXN5toA7bFykEVOiSZ1sWz6LYF3fTPGsX35bpQ5MbXuKKihVzf6KIFAgVuuTq+MFUNkwYTtOUT9lCNX5qP5GOHXt5HUtETqJClzPaOG8yZb/7Bw2zDvB1xSHE/fkZapYu7XUsETkFFbqc0tHUbWyecBuN9n/HOqvNpu7/pVPr9l7HEpEzUKHLHznHui/fpOpPT1DHneCr6rdy6dCRXFA80utkIpILFbr87uD2RHZMupn6R+JZGd4Iev+bzs3jvI4lInmkQhfIymTNJy9Sa8Uoqrkw5sQ8QLtBDxBZNMLrZCJyFlToIS518woOTLmFhmkJLI64hDL9/k2X+g29jiUi50CFHqJcRhqrP3ic+uvHEuZKMKfhU3TsdxsRRcK9jiYi50iFHoJ2JfzAiY9vo0n6ZhZEdqD6wNF0qVnT61gicp5U6CEkK+0ICZMepOGWiaRQjrktR9Ox1w2EhemwfZFgoEIPEdt+/oLwWcNpkrmTuaV60mDIKK6qUsXrWCLiRyr0IJd+ZB/rJ9xF450z2EJV5l36Hld2uUYn0xIJQir0ILblhw+44OsHqZ+1ny/LD6Dl0OdoX76s17FEJJ+o0IPQ8X072Dzhdhrsnct6arGu0ziuvqKz17FEJJ+p0IOJcyTNfZuK348kxqXxeeW/0XboE9QrVcLrZCJSAFToQeLI7k1sn3AzsYcWsTKsAendR9MtrrXXsUSkAKnQA11WFutnvUzU0uep6mB29L20HzyCkpFFvU4mIgVMhR7ADvy6ir1TbqbesVUsDm9J8WtH071xU69jiYhHVOgByGWcYP1HTxOT8BplXSSzY0dyVf87KRahzSkSyvLUAGbWFXgVCAfeds49e9Lye4D/B2QAKcCNzrktfs4qwN4Nizgy7Rbqn0hiftF2VO7/Kt3r1PY6logUArkWupmFA68DnYFkYImZzXTOJeRYbRkQ55w7ama3As8D/fMjcKhyJ46yburDxCa9S7oryxdNX6JT379SJDzM62giUkjkZQ+9FZDonEsCMLOpQB/g90J3zn2bY/2fgCH+DBnqdq/8mqyZd9IgYztfF+9G3cGj6BpV3etYIlLI5KXQqwNbc9xPBs70fbibgM9PtcDMhgHDAKKjo/MYMXRlHt3Phkn30mDbh/zqKvNVq7e5qls/nUxLRE7Jr5+imdkQIA445dWEnXPjgHEAcXFxzp/vHWy2L/qIol/eR2zmXj4vcx0thj5H50oVvI4lIoVYXgp9G1Ajx/0o32N/YGadgIeB9s65NP/ECz0nDuxm08Q7qJ/yJeuJZlX7MXTteLVOpiUiucpLoS8BYs0shuwiHwAMyrmCmV0EvAl0dc7t9nvKUOAcv857jzLzHiUm6yifVbyRNjc8Qb0ypbxOJiIBItdCd85lmNkdwJdkf23xHefcajN7Aoh3zs0EXgAuAKb59iR/dc71zsfcQeV46haSx99C3QM/8ovV41DXV+h56WVexxKRAJOnOXTn3Gxg9kmPPZbjdic/5woNWVls+mI0lRc/QzXn+LT6cNoN/gdlSkZ6nUxEApAOLfTI4W0JpEy+hZgjK1gS1pyw3qPp1aKF17FEJICp0AtaZjobP3mGqJWjKeeKMjPmEToNvIsSxSK8TiYiAU6FXoAOJMVz8P1bqJO2gflF2lKu3yv0blDf61giEiRU6AXAnThK4oePEbP+P6S50nza8Dm69PsbxYqEex1NRIKICj2fpa7+jvQZtxObnsxXxbpQc9AoetXUUbIi4n8q9Hzijh9gw+T7qffr+2x1FzLrojfo2nsg4TpsX0TyiQo9H+xaOpPw2fdQN2MPsy7oS5MhL9CjaiWvY4lIkFOh+1Hm4T1smngndXfOZoOrwdJLJ9CtS0+dTEtECoQK3R+cY/sPkyg59yGis47wSbkbaHXDU1xdvozXyUQkhKjQz9OJvVvZOvFW6uxdwC/UZU+nUfS+vL1OpiUiBU6Ffq6ystg69w3K//AU1VwGH1e+nfZDH6VpqeJeJxOREKVCPwfHd65n18Rh1Dy8jCXWlBM9XqbvJZd4HUtEQpwK/WxkZrB51vNU/fllyrkiTK8xgs6D76V08aJeJxMRUaHn1eEty9g/9WZqHVvHgvDWFL/mFf7UtJHXsUREfqdCz01GGknTRxK95k0iXUk+iv0X3a67meLF9E8nIoWLWukMDqybz7Hpt1P7xK/MKXoV1a5/iWvr1vI6lojIKanQT8GlHSJp6oPEbJrMIVeBGU1fo/s1gylaJMzraCIip6VCP0nq8tm4T4cTk5HCrBK9aDD4Ba6JquJ1LBGRXKnQfbIOp7Jp8nDqbP+Uja4aCy95l+7dr9HJtEQkYKjQnWPXovcpNudBojMP8VHpQcQN/Re9LizndTIRkbMS0oWesX8bWyfeRsye71hFHZLbvU3fKzvpsH0RCUihWejOsePbcZRa8DhVs04wrcLNtLvhMZqUvcDrZCIi5yzkCj1tdyK7Jt1M9IF44q0xh65+iX6XttZeuYgEvNAp9KxMtn4+ikpLXqCsC+eDqvfRecj9lLsg0utkIiJ+ERKFfix5JXsnD6PG0TUsCLuEsF4vcf1FzbyOJSLiV8Fd6Blp/PrJk1T7ZQzFXEner/UE3QfcSimdTEtEglDQFvrhxIUcnnYL0WmbmVOkAxWuHUX/RnW9jiUikm+Cr9DTDrPlw4eosWECB115pjV4mV79/kxkRLjXyURE8lVQFfr+X74k45O/UzNjJzOL9aBO/xe4rnZ1r2OJiBSIoCh0d3QfW6bcTa2tH5PkqjH/orfp1etPRITrZFoiEjoCvtBT4z8kfPb9RGXu58OS/bloyNNcW62S17FERApcngrdzLoCrwLhwNvOuWdPWl4MGA9cDKQC/Z1zm/0b9Y+yDu5k66TbqbnraxJcLRIvHUPfLl11Mi0RCVm5zkmYWTjwOtANaAQMNLOTr712E7DPOVcXeBl4zt9Bf+ccu+f/h6MvX0zlnfOYWuYmSt25gN5du6nMRSSk5WUPvRWQ6JxLAjCzqUAfICHHOn2Akb7bHwKvmZk555wfswKwavI/aLJhDPE0JKXjC/Rvf7kO2xcRIW+FXh3YmuN+MtD6dOs45zLM7ABQAdiTcyUzGwYMA4iOjj6nwBkthjBpfzE6D32AuNIlzuk1RESCUYF+KOqcGweMA4iLizunvfcWjRvTonFjv+YSEQkGefle3zagRo77Ub7HTrmOmRUBypD94aiIiBSQvBT6EiDWzGLMrCgwAJh50jozgT/7bvcDvsmP+XMRETm9XKdcfHPidwBfkv21xXecc6vN7Akg3jk3E/gPMMHMEoG9ZJe+iIgUoDzNoTvnZgOzT3rssRy3jwPX+TeaiIicDR0bLyISJFToIiJBQoUuIhIkVOgiIkHCvPp2oZmlAFvO8ekVOeko1CCn8QavUBoraLz+UNM5d8pTynpW6OfDzOKdc3Fe5ygoGm/wCqWxgsab3zTlIiISJFToIiJBIlALfZzXAQqYxhu8QmmsoPHmq4CcQxcRkf8VqHvoIiJyEhW6iEiQCLhCN7OuZrbOzBLNbITXefzNzDab2S9mttzM4n2PlTezr8xsg+/vcl7nPFdm9o6Z7TazVTkeO+X4LNto37ZeaWYtvUt+bk4z3pFmts23jZebWfccyx7yjXedmV3tTepzY2Y1zOxbM0sws9VmNtz3eFBu3zOM17vt65wLmD9kn753I1AbKAqsABp5ncvPY9wMVDzpseeBEb7bI4DnvM55HuNrB7QEVuU2PqA78DlgQBtgkdf5/TTekcB9p1i3ke9nuhgQ4/tZD/d6DGcx1qpAS9/tUsB635iCcvueYbyebd9A20P//YLVzrkTwG8XrA52fYD/+m7/F7jGuyjnxzk3n+xz5ud0uvH1Aca7bD8BZc2saoEE9ZPTjPd0+gBTnXNpzrlNQCLZP/MBwTm3wzn3s+/2IWAN2dcbDsrte4bxnk6+b99AK/RTXbD6TP+AgcgBc8xsqe+i2gCVnXM7fLd3ApW9iZZvTje+YN7ed/imGd7JMYUWNOM1s1rARcAiQmD7njRe8Gj7Blqhh4LLnXMtgW7A7WbWLudCl/27W9B+1zTYx+fzBlAHaAHsAEZ5msbPzOwCYDpwl3PuYM5lwbh9TzFez7ZvoBV6Xi5YHdCcc9t8f+8GPib7V7Jdv/0q6vt7t3cJ88XpxheU29s5t8s5l+mcywLe4v9+7Q748ZpZBNnlNsk595Hv4aDdvqcar5fbN9AKPS8XrA5YZlbSzEr9dhvoAqzijxfh/jPwiTcJ883pxjcTuMH3bYg2wIEcv7oHrJPmifuSvY0he7wDzKyYmcUAscDigs53rszMyL6+8Brn3Es5FgXl9j3deD3dvl5/UnwOnyx3J/vT5I3Aw17n8fPYapP9KfgKYPVv4wMqAHOBDcDXQHmvs57HGKeQ/WtoOtlziDedbnxkf/vhdd+2/gWI8zq/n8Y7wTeelb7/yKvmWP9h33jXAd28zn+WY72c7OmUlcBy35/uwbp9zzBez7avDv0XEQkSgTblIiIip6FCFxEJEip0EZEgoUIXEQkSKnQRkSChQhcRCRIqdBGRIPH/Ac/BXKMj7ipCAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Estimated curve strength of jenks: 6.7242\n",
"Average error 0.020398\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAq50lEQVR4nO3dd3hUZd7G8e8vjUBIKEmoARIQpFgQAyIo2AVkwYYCYgERG253RX1l0dVdXV0VKwsoKmBBQUFAwIJiAaQjVQEpoUYIgQCp87x/TNyNLCXCJGdmcn+uKxcz5xwm9+OE25NnTjHnHCIiEvoivA4gIiKBoUIXEQkTKnQRkTChQhcRCRMqdBGRMKFCFxEJEyp0CVtmttLMLjjJ1xhmZuMCk0ikbEV5HUCkrDjnWnmdQaQ8aQ9dRCRMqNAlbJnZRjO7xMwizGyIma03s91mNsHMahZvk2pmzsxuNrPNZvaTmT14lNeLNrO3zGyimcWYWTszW2hm+8xsp5k9Xb4jFPklFbpUBPcAVwKdgXpAFvDiYducB5wKXAwMNbMWJVeaWWXgAyAPuM45lw8MB4Y75xKAJsCEshuCyPGp0KUiuAN40DmX4ZzLA4YB15pZyc+QHnbOHXLOLQOWAWeWWJcAzADWA/2dc0XFywuAU8wsyTmX45ybV+YjETkGFbpUBI2A981sr5ntBVYDRUDtEtvsKPH4IFC1xPP2wBnA4+6XV7O7FWgGrDGzBWbWvSzCi5SWjnKRimALMMA59/XhK8wstRR/fxawHPjUzC5wzu0EcM79APQxswjgauA9M0t0zh0IXHSR0tMeulQEI4DHzKwRgJklm1nPX/MCzrl/Am/iL/Wk4tfpZ2bJzjkfsLd4U1/gYov8OtpDl4pgOGDALDOrB+wC3gEm/5oXcc79zcwqAZ+Y2UVAF+BpM6sCbAJ6O+cOBTa6SOmZbnAh4crMNgP9nHNzvM4iUh405SJhycySgWRgo8dRRMqNCl3Cjpm1BX4AnnfObfY6j0h50ZSLiEiY0B66iEiY8Owol6SkJJeamurVtxcRCUmLFi36yTmXfKR1nhV6amoqCxcu9Orbi4iEJDPbdLR1mnIREQkTKnQRkTChQhcRCRNBdep/QUEBGRkZ5Obmeh2lTMTGxpKSkkJ0dLTXUUQkDAVVoWdkZBAfH09qaipm5nWcgHLOsXv3bjIyMkhLS/M6joiEoeNOuZjZq2a2y8xWHGW9mdlzZrbOzJabWZsTDZObm0tiYmLYlTmAmZGYmBi2v32IiPdKM4f+Gv6ryh1NV6Bp8dcg4OWTCRSOZf6zcB6biHjvuFMuzrk5x7kJQE/gjeI7ucwzs+pmVtc5tz1QIUVEgp7PB/n74VAWHNrr/zN3L+QfpCDvEHuys9mzbz/79uWQdHYPmrTuFPAIgZhDr4//jjA/yyhe9j+FbmaD8O/F07BhwwB868B77rnnePnll9mxYwf33XcfQ4YM4YMPPqBZs2a0bNnS63gi4hVfEWRvgd3rYc8G/9fu9ZD1IxzIhNxscEe+v0k0/vsd/nzPw3nxtYO20EvNOTcSGAmQnp4elFcFe+mll/jkk09ISUn5z7IPPviA7t27q9BFKgLnIGsjZCyEHcv9pb17nb+4i/L/u1l0FQ5Vbciu6AZsq9yKDIthQ04Mu4sqk+3iyKYqVaslUbdWEg1q1aRhrZqk1vZ/ta9UNtUbiFfdCjQo8TyleFnIueOOO9iwYQNdu3ZlwIABrF+/nr59+zJlyhS++OILHn30USZOnEiTJk28jioigeAc7N3kL+/tS2H7Mv9XbrZ/fWQM1GyCL7Epe+pfxAZfHZYeqMmXexKYuyuawv3+zWrGxdC8TjzNWyTQtk48zevG07RWPJVjIst1OIEo9CnAYDN7GzgHyA7E/PnDH65k1bZ9Jx2upJb1Evjrb1oddf2IESOYMWMGs2fPZurUqQB06NCBHj160L17d6699tqA5hGRcubzwU9rYdPXsGkubJ4L+4r3PyNjoHYraHUV+clnsNqaMDu7Ft9u2seSVXs5VFAEQEJsFGekVGdQ82qc2aA6Z6ZUp3ZCpaA46OG4hW5mbwEXAElmlgH8Ff+UEM65EcB0oBuwDjgI9C+rsCIiv1phHqz7FFa+D+s+gUN7/Mur1oFG50LDDhyofTbzc2oxf3MO327cw3dzsyn0HcDsR1rUSeD6tg04q2F1zkipTmpilaAo7yMpzVEufY6z3gF3ByxRsWPtSYuIHFNhPmyY7S/xNdMgbx/EVodTu0Lq+RQ1OJcVB2sw54ef+HLJTyzevJNC3w5iIiM4I6Uat3VqTLu0mpzdqAYJsaFzZndQnSkarOLj49m/f7/XMUTkWHxFsOFzWDER1kz1z4PHVoMWPaDVVexKPofP1+1lzupMvv5wHVkHCwA4rX4Cgzo15rymSbRpWIPY6PKd9w4kFXop9O7dm9tuu43nnnuO9957Tx+KigSTQ3vhq2dg2VuQsxMqVYPmV0Crq/ixWltmrtnDrI93sGTLHJyD5PhKXNi8Fp2bJdPxlCSSqlbyegQB49k9RdPT093hN7hYvXo1LVq08CRPeakIYxQpF3k5sGQsfPk0HPwJmnXFnXEdq+PPZdrqLGat3MkPu3IA/174ZS3rcHGLWrSsmxC0c+ClYWaLnHPpR1qnPXQRCS1bF8OiMbBiEuTnQKPz2HbOG7y3PYnJM7ayPnMhkRHGOWk1ueGchlzaqg71q1f2OnW5UKGLSPArKoQ1H8K8l2HLfIiuwqFmPZlVuQuvbExi+Rt7MNtDu9SaDDgvjW6n1aVGXIzXqcudCl1EgtehLFj8BswfCfsycDVSWXfWg7yQ1Y5pSw5S6HOcVt/xYLcWdD+zLnWrVYw98aNRoYtI8MneCl8P98+RFxwkN6UjM+r9nn+ub8S2uQUkxhUw4Lw0rktP4ZRa8V6nDRoqdBEJHtkZ/iNWFr+Bcz52pvbkpUOX8sa6BCIMOjWrztC2DbioeW1ionQHzcOp0EXEe9kZ/qNVlozFOR8/1LuSYVmX8c2qqiRVrcTvLm7I9W0bUK+CfLh5olToJezdu5c333yTu+66y+soIhXD3i3w1dOweCwOWJLYnfszL2HtuhqcXr8aT1+XyhVn1KVSVOie7FOeVOgl7N27l5deeul/Cr2wsJCoKP2nEgmY/AP+qZWvn8M5H/Ord+MvOy9h69YkupxWh793TKVNwxohfby4F9RSJQwZMoT169fTunVroqOjiY2NpUaNGqxZs4ZZs2bRvXt3Vqzw31r1qaeeIicnh2HDhrF+/XruvvtuMjMzqVKlCqNGjaJ58+Yej0YkCDkHKyfBrIdg31YWJlzC7zN7sCu/Fte3bcCgTo1pULOK1ylDVvAW+kdDYMd3gX3NOqdD18ePuvrxxx9nxYoVLF26lM8//5wrrriCFStWkJaWxsaNG4/69wYNGsSIESNo2rQp8+fP56677uKzzz4LbHaRULd9Gcy4HzZ9zeaYU/hj3lBW7W1Fv/MbMfC8NGolxHqdMOQFb6EHgXbt2pGWlnbMbXJycvjmm2/o1avXf5bl5eWVdTSR0JGzCz59BLdkHDmR1fh7wa18xKXcdFETRnVIrZAnAJWV4C30Y+xJl5e4uLj/PI6KisLn++/9AnNzcwHw+XxUr16dpUuXlnc8keBWmAfzR+D74p+4/EOMKerGiKJruO78VnzRqQnVqoTOZWlDhQ7kLOFYl8mtXbs2u3btYvfu3eTl5f3njkYJCQmkpaXx7rvvAuCcY9myZeWWWSToOAdrpuF74Rz4eChf5DXj8vwn2Xj2/Uy/9wr+0qW5yryMBO8eugcSExPp2LEjp512GpUrV6Z27dr/WRcdHc3QoUNp164d9evX/8WHnuPHj+fOO+/k0UcfpaCggN69e3PmmWd6MQQRb+1chZsxBPvxC34khWH5Q6h5RhdeubQZjRLjjv/35aTo8rnlrCKMUSqgwjyY8yS+L5/hALE8mX8Nq+v34qEep3NGSnWv04UVXT5XRMrOzpUUvnMLUXu+Z1JRJ0ZVvpW7erbl4TPr6TjycqZCF5ETVrTqQ3wTbyOrMJYhvvs5rdPVvN+5MVViVC1eCLr/6s65sP2/ulfTWyIB5xw7pz5K7UVP8Z2vCa+kPMrQqzuTmqR5ci8FVaHHxsaye/duEhMTw67UnXPs3r2b2FidPCGh7dC+LH4cM4CWWZ8x3c6nqOdwnmvTOOz+zYaioCr0lJQUMjIyyMzM9DpKmYiNjSUlJcXrGCInbOWiL6kx9Vaa+TL5qO6ddLzxEarpxKCgEVSFHh0dfdwzM0Wk/OUWFDHzzWe5fMM/yLYEVnV5h67nXuZ1LDlMUBW6iASfJT/uZOObf+CqgmlsiD+L2gPeonbNul7HkiNQoYvIERUW+Xhl1iJazx3MVRFr2NriVhpf+0+IVG0EK70zIvI/tuw5yL/GT+aPmQ9RNzKbg91HUP/sPl7HkuNQoYvIL0xeupUZ74/lSZ4lqnIc0TfOIDrlbK9jSSmo0EUEgEP5Rfx18ndUWfoKL0SPoyipJTE3vgPVdGRWqFChiwgbMnMYPG4BvXe/yE3RH+Nr1o2Ya0ZBpapeR5NfQYUuUsFNW76dRyfO5V/2LB2ilkHH3xFx8TCI0NW1Q02pCt3MugDDgUhgtHPu8cPWNwReB6oXbzPEOTc9sFFFJJDyC338ffpqPvpmMe/GPUkDtxW6Pw9tbvI6mpyg4xa6mUUCLwKXAhnAAjOb4pxbVWKz/wMmOOdeNrOWwHQgtQzyikgA7M7J487xi9m1cSUzE56kGgew3pOgcWevo8lJKM0eejtgnXNuA4CZvQ30BEoWugMSih9XA7YFMqSIBM7KbdkMemMRyTlrmBn/JJWiIqDfVKjX2utocpJKM0lWH9hS4nlG8bKShgH9zCwD/975PUd6ITMbZGYLzWxhuF6vRSSYTV2+jWte/obTi1YysfLfqVQ5DgbMVJmHiUB96tEHeM05lwJ0A8aa2f+8tnNupHMu3TmXnpycHKBvLSLH4/M5npq5lsFvLqFf4ve87B4jMqEODJgBSad4HU8CpDRTLluBBiWepxQvK+lWoAuAc26umcUCScCuQIQUkRN3IK+Q3729lE9W7+Sxpj/Qd+vfsFotod8kqKodq3BSmj30BUBTM0szsxigNzDlsG02AxcDmFkLIBbQnIqIx3bty+X6kXOZvXYXY9M30DfjYSylLdwyVWUeho67h+6cKzSzwcBM/IckvuqcW2lmjwALnXNTgD8Bo8zsD/g/IL3F6fY8Ip76fud++o9ZQNbBfKae9yMtvn0I0s6HPm9DjO4sFI5KdRx68THl0w9bNrTE41VAx8BGE5ET9c36n7h97CJioyOZdcFmUuY8AKdcAtePg+jKXseTMqJTwUTCzPtLMrj51W+pkxDLzIu2kTLnL9DkYrh+vMo8zOnUf5Ew4ZzjxdnreGrW95zbOJFX2mygyrTf+adZeo+HaN3PNtyp0EXCgM/neGTqKl77ZiNXta7Hk3U+JmrqY5BaPGeuPfMKQYUuEuIKinz85b3lvL9kK4M6NOB+37+xz8fB6ddBzxcgqpLXEaWcqNBFQlhuQRGD31zMJ6t38cBFdbltx1Bsw2zodC9c+CCYeR1RypEKXSRE7c8tYODrC/l24x6GXxJHz9W3QdZG6PECtLnR63jiARW6SAjanZPHLWMWsHr7PsZfeIAOC+6AyBi4+UNo1MHreOIRFbpIiNm5L5e+o+aRkXWIyZ220mreEEhu7v/ws3qD47+AhC0VukgI2Z59iD4j55G5P4+Z564kdd7f/Eey9B4PsdW8jiceU6GLhIite/1lnnUgn+kdv6fR3L9Bi9/A1aN1jLkAKnSRkLBlz0H6jJpH9qECPrxwB40+HwrNusK1YyAy2ut4EiR06r9IkNu0+wC9R85jf24hUy7PJXXOH6FRR+ilMpdf0h66SBD78acD9B01j9yCIiZeWYW0qbdArRbQ5y2d/Sn/Q4UuEqQ2/nSA3iPnUlDkeLd3fU6ZfCVUqQl934XYhOP+fal4VOgiQSgj6yA3jJ5PfqGPCTc155QPr4aiPP9x5gl1vY4nQUqFLhJkdmTn0nfUfPbnFvBW/9Y0/eQm2LsJbvwAajX3Op4EMRW6SBD5KSePG0bPY8+BfMb2b0OreX+ELfOg12uQqnvIyLHpKBeRILH3YD79Rs9n695DjLmpNWctuBdWT4HL/w6trvI6noQA7aGLBIF9uQXc+Mq3bPjpAK/d0Iq2c++GdR/DpY/AuXd7HU9ChApdxGMH8grpP2YBa3bs45XrmtDh61th60Lo/iyk9/c6noQQFbqIh3ILirj19QUs3bKXUVfWo9NXN8Oe9f4585Y9vY4nIUaFLuKRwiIfg99czPwf9zCqazUu+rofHMqCG96Dxp29jichSIUu4gGfz3HfxO/4ZPUuXrwALpl3E2Bwy1Sod5bX8SRE6SgXkXLmnOMfH61m4uIMnjl7D1csHgQxcXDrLJW5nBQVukg5G/HFBkZ9+SNPnbqGK1f/Hqo3ggGzILGJ19EkxGnKRaQcvf3tZp6YsYZnGnzJVZte9t+c4vpxULm619EkDKjQRcrJjBXbeeD95byYPJkrMif4Txa66t8QVcnraBImNOUiUg6+Wf8Tv31rCf+qMYkr9k+AtrfBNa+qzCWgtIcuUsZWbM1m0BsLGVb1fa46OBHaDoRuT4KZ19EkzGgPXaQMbdt7iAGvLeD3UZPomzcB2twMXVXmUjZKVehm1sXM1prZOjMbcpRtrjOzVWa20szeDGxMkdCzL7eA/q9+y8357zCw6B1ofYP/dP4I7UdJ2TjulIuZRQIvApcCGcACM5vinFtVYpumwP1AR+dclpnVKqvAIqGgoMjH3eMW0TPrVe6K/ADO7AM9nleZS5kqzU9XO2Cdc26Dcy4feBs4/CITtwEvOueyAJxzuwIbUyR0OOf4v0nf0XHj8/4yb3Mz9HwJIiK9jiZhrjSFXh/YUuJ5RvGykpoBzczsazObZ2ZdjvRCZjbIzBaa2cLMzMwTSywS5F76fD11lg7njqip/g9AfzNce+ZSLgL1UxYFNAUuAPoAo8ys+uEbOedGOufSnXPpycnJAfrWIsFj8tKt7Pnkaf4QPRHXuq8+AJVyVZpC3wo0KPE8pXhZSRnAFOdcgXPuR+B7/AUvUmF8++MeFrz3DA9Fj6eoRQ/sN5ozl/JVmp+2BUBTM0szsxigNzDlsG0+wL93jpkl4Z+C2RC4mCLBbUNmDhPfGM4jUaMpSLuIyGtegUid5iHl67iF7pwrBAYDM4HVwATn3Eoze8TMehRvNhPYbWargNnAvc653WUVWiSY7M7JY9ToF3nMPU9+vXOI7jMeomK8jiUVkDnnPPnG6enpbuHChZ58b5FAyS0o4okXXuL+vcMoSG5J3MBpEJvgdSwJY2a2yDmXfqR1muATOUE+n2PEa69z396Hya3elLgBk1Xm4ikVusgJGjNpCgMz7udg1YYkDJoKVWp6HUkqOBW6yAmY9Pl8rvjudxTFJFDj9qkQl+R1JBFdbVHk15qzYj0tPhtItcg8ovt/iCXU8zqSCKBCF/lVVm/aRty7fWkasZXC694mqt7pXkcS+Q9NuYiU0vbdezjwWi9a2/fkXPEysS0u8zqSyC+o0EVKYf+Bg2x5uRdtfCvZceEzVG97vdeRRP6HCl3kOAoKi1jy0i20K1zI+naPUL/zLV5HEjkiFbrIMTjnmDPyT3Q6MJMVp9xB0yt+63UkkaNSoYscw7xxw7h41xiWJ/+G02543Os4IsekQhc5iuWTn+Xc9c+yKP5CTrt9jC6DK0FPhS5yBN9//QEtFz/Mopi2tBr8DhFR0V5HEjkuFbrIYbb8+D3JHw9mU2QD0u6cQGylSl5HEikVFbpICXv25bBvbD+iKaRSn3HUrKHrs0joUKGLFMstKOKbEXfRyreWnRc+RUrTM7yOJPKrqNBF8F8K9/3Rf6f7wclsaHIjTTr38zqSyK+mQhcBPnh7JL12PM3mmh1o3PdZr+OInBAVulR4n01/hyvWPsD2KqfSYNA7uheohCwVulRoS7+eSfv597ArJoW6d0/DdMchCWEqdKmwNnw3l8azbiErsiY175hOVNVEryOJnBQVulRIu7ZtpurEvhyyykT3n0JcYn2vI4mcNBW6VDg5Bw+S+cr1xLsccq4eR60GzbyOJBIQKnSpUAqLfCx4aSCtilaxoeM/aXJGB68jiQSMCl0qDOccU8c8xoU501jZeACtLuvvdSSRgFKhS4Uxecokum15hvXVzqVVv6e8jiMScCp0qRA+nreYjov/QFZMXdJufwsiIr2OJBJwOoNCwt6CddupM/1WqkbmE9F/AhFVangdSaRMaA9dwtq6nfvZPu52To/YQNGV/6ZSvVZeRxIpMyp0CVuZ+/OYNnooPfiCvefcS9Uze3odSaRMqdAlLB3ML2T4qFHcnf8a2Y0up/rlD3gdSaTMqdAl7BQW+XhyzAT+nP13cqs1plrfVyBCP+oS/kr1U25mXcxsrZmtM7Mhx9juGjNzZpYeuIgipeecY+Q7E/n9tj8RVTmBqre8B5XivY4lUi6OW+hmFgm8CHQFWgJ9zKzlEbaLB34HzA90SJHSen/qFPqt/S2uUjWq3j4TaqZ5HUmk3JRmD70dsM45t8E5lw+8DRzp06W/AU8AuQHMJ1JqX302jUsW3k5eTHUS7pgJNRp5HUmkXJWm0OsDW0o8zyhe9h9m1gZo4JybdqwXMrNBZrbQzBZmZmb+6rAiR7Nq3kxafzGAA1E1iL9jBhE1VeZS8Zz0J0VmFgE8DfzpeNs650Y659Kdc+nJyckn+61FANi8+GNSZ9xIVkRNqgyaQWyiylwqptIU+lagQYnnKcXLfhYPnAZ8bmYbgfbAFH0wKuUhc/kskqfcwA6SiRgwnWq1VeZScZWm0BcATc0szcxigN7AlJ9XOueynXNJzrlU51wqMA/o4ZxbWCaJRYrt/24a8ZNuIINaFN44hfoN9AGoVGzHLXTnXCEwGJgJrAYmOOdWmtkjZtajrAOKHEne4reoPPFG1rn67LtuEs2aNPE6kojnSnVxLufcdGD6YcuGHmXbC04+lsjRFc0dQaWZ9zHX15IDV4/lkpaneB1JJCjo9DkJKb7ZjxM58z5mFqWzuesbXNJaZS7yMxW6hI7Z/yDii3/wXlEnfuj8Atef29TrRCJBRYUuoWH2P+CLx5lQ2Jnvzn6Muy9u7nUikaCjG1xI8CtR5nOaD2V4j9MxM69TiQQdFboEt5/LvKgz09MeYGTvNkRGqMxFjkSFLsHJOZj9GMx5kneLOjOx/n28dmNbYqI0SyhyNCp0CT7Owaz/g7kvMMF3EW/W+gNjb2lH5Rjd2FnkWFToElx8PvjoXlgwmvHucsZWv4u3BrQnPjba62QiQU+FLsHDVwRTfgtLx/EqPXgjrj8TBranRlyM18lEQoIKXYJDUQG8fzusmMi/rRevR/dmwm3tqZUQ63UykZChQhfvFebBewNgzVRejOzHGLuSCbe1J6VGFa+TiYQUFbp4q6gA3u0Pa6fxdOStvOG68NbAc2icXNXrZCIhR8eAiXd8PvjgTn+ZR/nL/M2B7WlRN8HrZCIhSYUu3nAOpv0BvnuXf0fdwBu+LowfeA4t66nMRU6UCl3Kn3Mw/V5Y9BpjI6/mpaIrGXfrObSqV83rZCIhTYUu5cs5mPkALBjFm5E9eaqoN+MHnsNp9VXmIidLhS7lxzn4+CGY9xITIrvzRNENjL+tvcpcJEBU6FI+nIPP/gbfPM+7EV34u+8mxg1UmYsEkg5blPIx50n48l9MtEt5wgbyzqD2nFon3utUImFFhS5l76tnYfZjfMAFPB19O+8O6kBaUpzXqUTCjgpdytbcl+CTvzLVdeS5uHt457YOOgNUpIyo0KXsLBgNM+9nhq8dz1f7M2/f1kHXZhEpQyp0CTzn4Mt/wWd/41NfG0Yk3c9bt55HTV01UaRMqdAlsAoO4SYPxla8x+SiDryXMoQ3bu5Agq5nLlLmVOgSOPu24d7qA9uX8c+C69nU8nZGX9+aSlG605BIeVChS2BkLMS93Ze8A/sYnP9HGpx7Dc9f0ZII3dBZpNyo0OXkLXsHN+UeMqlBv9xhXNPlUgZ1aoyZylykPKnQ5cT5iuDTh+Hr4SyPOp2BBwdzf6/zuLpNitfJRCokFbqcmMJ8mHQbrPqA9yIu47GCWxh+Szs6NUv2OplIhaVCl18v/wC8cyOs/5Qnim5gapVrmXBzW5rW1qn8Il4q1cW5zKyLma01s3VmNuQI6/9oZqvMbLmZfWpmjQIfVYLCoSzc2KvwrZ/NvQWDWFi/Hx/c1VFlLhIEjlvoZhYJvAh0BVoCfcys5WGbLQHSnXNnAO8B/wx0UAkC+3fiG9ONwozF3JX/W3xn9mPcwHNIrFrJ62QiQun20NsB65xzG5xz+cDbQM+SGzjnZjvnDhY/nQfoU7Fwk7WRwtGXkbdrA/3z/kzry2/iqV5n6BhzkSBSmjn0+sCWEs8zgHOOsf2twEdHWmFmg4BBAA0bNixlRPHczlXkv9aTQ4cOcod7iP43XMtlrep4nUpEDhPQD0XNrB+QDnQ+0nrn3EhgJEB6eroL5PeWsuE2zyPvjV7sLYjkwaqP87ebr+KUWlW9jiUiR1CaQt8KNCjxPKV42S+Y2SXAg0Bn51xeYOKJZ5wjd95oomYOYbsvkVcbP83wPl2oWkkHRokEq9L861wANDWzNPxF3hvoW3IDMzsL+DfQxTm3K+AppXwV5JL17j3U+H4CnxedycYLhvPIRa115qdIkDtuoTvnCs1sMDATiARedc6tNLNHgIXOuSnAk0BV4N3if/SbnXM9yjC3lBG3dzN7Xr2exH2rGB1xLa36Pc4tp+hkIZFQUKrfn51z04Hphy0bWuLxJQHOJR7IWf0p7t3+RBfl82yth7nx5jt1SKJICNGEqIDPx9apj1Fn8dNs8NVl8bnP89vLL9KVEkVCjAq9gsvPyWLTKzfSNOtLPok8j+Qb/831TXQagUgoUqFXYBtXzCdm0s2kFu1icr3fctFNDxFfWbeJEwlVKvQKqKjIx1fvPkO71Y+z36qy5KKx9Ox8hdexROQkqdArmB+2bGfbuDvonPc5q6q0ofYtr9Outs7aFQkHKvQKIr/Qx7tTp9FhyZ85z3aysvk9tOw1DIvUj4BIuNC/5gpgyaY9fPXWEww69AqHoqux/5r3adXiQq9jiUiAqdDDWNaBfJ77aBFtlw3lnshv+aleJ5L6jYG4JK+jiUgZUKGHIZ/P8c7CLUz9aCr/8D1DSuRuci/8K0nn/x4iSnVPExEJQSr0MLNsy16GTl7B2dvf5vXot3DxtYm4fgaxDdp5HU1EypgKPUxs2XOQp2atZcmyJTwR+zrnRi/FndoV6/kSVKnpdTwRKQcq9BCXfbCAFz9fx2tfb6RHxFd8UvkVoqMi4ZKnsLYDQVdIFKkwVOghKq+wiLFzN/HC7HUUHNrP63Xe59ysDyGlA1wzGqrV9zqiiJQzFXqIyS/0MWlxBi/MXkdG1kHuq7eMgbGvEZ21C877A1z4f6Bjy0UqJP3LDxG/LPJDXFF3H9MajKZa5kKo1wb6vAkN2nodU0Q8pEIPcocX+dkpcYxtMpfUVSOwmDjo8Ty07qfDEUVEhR6s9ucW8Pa3Wxjz9Y9sy87lzAbVea7DIc5adi+24ns4vRdc/g+oqrsJiYifCj3IbNt7iDFf/8jb325hf14h7RvX5OlL4zln4wjs04lQvSH0mwin6CZRIvJLKvQg8V1GNq98tYGpy7fjgG6n1+XONpVp+f3LMG0sRFWC8/8M5/8RYuK8jisiQUiF7qGD+YV8uGwb4+dvZnlGNnExkdzcIZXbmudRZ+VomDABnIO2t/rLPL6215FFJIip0D3w/c79vDl/MxMXZ7A/t5Bmtavy8G9acm3SRuIWDoNxMyGqMrS5CTrcAzVSvY4sIiFAhV5O9h7MZ+ry7by/ZCuLNmURExlBt9Pr0K9tHc7ePxub/xBsXwZVkuCCB6DtQIhL9Dq2iIQQFXoZyi/0MXvtLt5fvJXP1uwiv8hHs9pVeaBbc3qdGk2NVeNg0qtwIBOSToXuz8CZfSC6stfRRSQEqdADrLDIx7wNe5i+YjsffbedrIMFJFWNoV/7Rlzdpj6t3Dps/uPw+fvgK4Rml8M5t0PjC3XdFRE5KSr0AMgv9PH1+p+Y8d0OZq3aQdbBAqrERHJxi9pcfVZ9zm9cjai1H8JHv4WMBRAT759SaXcbJDbxOr6IhAkV+gnKOpDPnB8y+WzNLj5bs4v9uYVUrRTFJS1q0fX0unRulkxshA+WjIXnn4T926BmY+jyBLTuC7EJXg9BRMKMCr2UnHOs3r6f2Wv9Bb5kcxY+BzXjYrisZR26nV6H85omUSkyArYtho+fhhUT4eBP0KA9/OZZOOVSnaIvImVGhX4Uzjk27znI3PW7+ab466ecPABOr1+NwRc15cJTkzkjpTqRuVmw8SuY+QWs/wz2bIDISnBqF/91VppeqvlxESlzKvRizjm2Zecyf4O/vOeu383WvYcAqBVfifNOSaTDKUlc0CyZWrFFsHkurBkLH83xH26Ig+g4aNQBOv4eWvaEytW9HJKIVDAVttBzC4pYuS2bxZv2snhzFos3Z7Fzn38PvEaVaNo3TuSOzo05t0kSTZKqYNuXwA/jYeIc2PIt+AogIhoatIML7ofGnaH+2RAZ7fHIRKSiqhCFXljkY33mAVZuy2bF1n0s2ZLFyq37yC/yAdCgZmXaN06kTcMapKfWoEXteCL2Z/iPSFk4H9ZOh+wtgEG91nDuXZDWGRqeCzFVPB2biMjPSlXoZtYFGA5EAqOdc48ftr4S8AZwNrAbuN45tzGwUUvnQF4h63blsHLbPn+Bb9vHmu37yCv0l3dsdASn169G/46ptKkfR3r1/STmb4PMryFzDaxaC5lrIS/b/4JRsdD4ArjwAWjWRTdcFpGgddxCN7NI4EXgUiADWGBmU5xzq0psdiuQ5Zw7xcx6A08A15dF4J/lFBf39zv3s25XDj/s3M+6nfvI3ruHBDtANQ5Qr9IhLqjh4/bGBaRVPkDdyP0kFO0h4kAmrNkK324D3H9fNC4ZkpvDGb2gVguonw61W2kaRURCQmn20NsB65xzGwDM7G2gJ1Cy0HsCw4ofvwe8YGbmnHME2PyJw6m94t/4ioqogY8O5uN8fMRFFFDVHSAi1vfLv5BV/AVQJRHiavlvCpHWyX/Rq+qNoEYj/6n3unaKiISw0hR6fWBLiecZwDlH28Y5V2hm2UAi8FPJjcxsEDAIoGHDhicUOLZaMllVmxIXW4nYyjHEV44lLjaGiJjKULkGxFb3H13yi8c1IS5Je9oiEtbK9UNR59xIYCRAenr6Ce29n3lJX7ikb0BziYiEg9KctrgVaFDieUrxsiNuY2ZRQDX8H46KiEg5KU2hLwCamlmamcUAvYEph20zBbi5+PG1wGdlMX8uIiJHd9wpl+I58cHATPyHLb7qnFtpZo8AC51zU4BXgLFmtg7Yg7/0RUSkHJVqDt05Nx2YftiyoSUe5wK9AhtNRER+DV36T0QkTKjQRUTChApdRCRMqNBFRMKEeXV0oZllAptO8K8ncdhZqGFO4w1fFWmsoPEGQiPnXPKRVnhW6CfDzBY659K9zlFeNN7wVZHGChpvWdOUi4hImFChi4iEiVAt9JFeByhnGm/4qkhjBY23TIXkHLqIiPyvUN1DFxGRw6jQRUTCRMgVupl1MbO1ZrbOzIZ4nSfQzGyjmX1nZkvNbGHxsppm9rGZ/VD8Zw2vc54oM3vVzHaZ2YoSy444PvN7rvi9Xm5mbbxLfmKOMt5hZra1+D1eambdSqy7v3i8a83scm9Snxgza2Bms81slZmtNLPfFS8Py/f3GOP17v11zoXMF/7L964HGgMxwDKgpde5AjzGjUDSYcv+CQwpfjwEeMLrnCcxvk5AG2DF8cYHdAM+AgxoD8z3On+AxjsM+PMRtm1Z/DNdCUgr/lmP9HoMv2KsdYE2xY/jge+LxxSW7+8xxuvZ+xtqe+j/uWG1cy4f+PmG1eGuJ/B68ePXgSu9i3JynHNz8F8zv6Sjja8n8IbzmwdUN7O65RI0QI4y3qPpCbztnMtzzv0IrMP/Mx8SnHPbnXOLix/vB1bjv99wWL6/xxjv0ZT5+xtqhX6kG1Yf6z9gKHLALDNbVHxTbYDazrntxY93ALW9iVZmjja+cH6/BxdPM7xaYgotbMZrZqnAWcB8KsD7e9h4waP3N9QKvSI4zznXBugK3G1mnUqudP7f3cL2WNNwH1+xl4EmQGtgO/AvT9MEmJlVBSYCv3fO7Su5Lhzf3yOM17P3N9QKvTQ3rA5pzrmtxX/uAt7H/yvZzp9/FS3+c5d3CcvE0cYXlu+3c26nc67IOecDRvHfX7tDfrxmFo2/3MY75yYVLw7b9/dI4/Xy/Q21Qi/NDatDlpnFmVn8z4+By4AV/PIm3DcDk71JWGaONr4pwE3FR0O0B7JL/Ooesg6bJ74K/3sM/vH2NrNKZpYGNAW+Le98J8rMDP/9hVc7554usSos39+jjdfT99frT4pP4JPlbvg/TV4PPOh1ngCPrTH+T8GXASt/Hh+QCHwK/AB8AtT0OutJjPEt/L+GFuCfQ7z1aOPDf/TDi8Xv9XdAutf5AzTescXjWV78j7xuie0fLB7vWqCr1/l/5VjPwz+dshxYWvzVLVzf32OM17P3V6f+i4iEiVCbchERkaNQoYuIhAkVuohImFChi4iECRW6iEiYUKGLiIQJFbqISJj4f9YNo0IuYeuoAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# Define sigmoid function, determine curve strength of methods so far\n",
"\n",
"def sigmoid(x, strength):\n",
" scaler_input = sklearn.preprocessing.MinMaxScaler(feature_range=(-0.5, 0.5))\n",
" # Avoid having points at inf\n",
" epsilon = 1e-5\n",
" x = scaler_input.fit_transform(x.reshape(-1, 1)).flatten()\n",
" if strength != 0:\n",
" output = 1 / (1 + np.exp(-(x * strength)))\n",
" else:\n",
" # 0-strength line. Make sure final output has some variance\n",
" output = np.linspace(0, 1, len(x))\n",
" scaler_output = sklearn.preprocessing.MinMaxScaler(feature_range=(0 + epsilon, 1 - epsilon))\n",
" return scaler_output.fit_transform(output.reshape(-1, 1)).reshape(-1)\n",
"\n",
"def sse(y_true, y_pred):\n",
" return np.sum((y_true - y_pred) ** 2)\n",
"\n",
"def sigmoid_classes(strength, n_classes):\n",
" x = np.arange(n_classes)\n",
" return dist.ppf(sigmoid(x, strength))\n",
"\n",
"def est_curve_strength(R, name):\n",
" y = dist.cdf(R)\n",
" x = np.arange(len(R))\n",
" res = scipy.optimize.minimize(lambda x0: sse(y, sigmoid(x, x0[0])), x0=[10], options={'gtol': 1e-10})\n",
" print(f\"Estimated curve strength of {name}: {res.x[0]:.4f}\")\n",
" print(f\"Average error {res.fun:.6f}\")\n",
" plt.title(name)\n",
" plt.plot(x, sigmoid(x, *res.x), label='fit')\n",
" plt.plot(x, y, label='true')\n",
" plt.legend()\n",
" plt.show()\n",
"\n",
"est_curve_strength(linear_R, 'linear')\n",
"est_curve_strength(gaussian_R, 'gaussian')\n",
"est_curve_strength(jenks_R, 'jenks')"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "1809fa9d",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
" fun: 1.3198519807543092\n",
" hess_inv: array([[0.04675019]])\n",
" jac: array([6.63101673e-06])\n",
" message: 'Optimization terminated successfully.'\n",
" nfev: 44\n",
" nit: 7\n",
" njev: 22\n",
" status: 0\n",
" success: True\n",
" x: array([6.48649694])"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"initial_strength = 16 # linear strength\n",
"res = scipy.optimize.minimize(lambda x0: measure_loss(sigmoid_classes(x0[0], n_R)), x0=[initial_strength])\n",
"sigmoid_R = sigmoid_classes(res.x[0], n_R)\n",
"res"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "48be5c0c",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7fdc8859fb80>]"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD4CAYAAAAEhuazAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcpklEQVR4nO3de3xUZ73v8c8vIQmEOyRACKQBGqC0oRRSaLXSVloLRIttbTd1n23VHtmeY/fZe3v21vrS81Lrrkd9WT11W90Hj9XWrVTtRWlJb/Ri1dpCoBDuNA0Ucg+Ea26TmXnOH7MCIw3XZLJmZn3fr9e8MvPMYub3dE2/efKsZ9Yy5xwiIhIsGX4XICIiA0/hLyISQAp/EZEAUviLiASQwl9EJIAG+V3AucjLy3PFxcV+lyEiklI2bNhwwDmX39tzKRH+xcXFVFZW+l2GiEhKMbN3T/ecpn1ERAJI4S8iEkAKfxGRAFL4i4gEkMJfRCSAFP4iIgGk8BcRCSCFv4hIknpiQy2r1u1LyGsr/EVEktSTb9Xy+IbahLy2wl9EJEm1hyLkZmcm5LUV/iIiSaojFGFIls/hb2YPm1mzmW2Na/uamdWZ2SbvtjTuuS+ZWbWZ7TKzm+LaF3tt1WZ2b/91RUQkvbSHIgzNScwp2M5n5P9zYHEv7d93zs3xbhUAZjYLWA5c6v2bH5lZppllAg8BS4BZwJ3etiIicor2UIQhCZr2OedfKc6518ys+Bw3XwY85pzrAvaYWTUw33uu2jlXA2Bmj3nbbj/3kkVEgqEjFCbX72mfM7jHzKq8aaHRXlshsD9um1qv7XTt72FmK8ys0swqW1pa+qFMEZHU4ZyjvTt5D/j+GJgGzAEagAf6WlAP59xK51yZc64sP7/XaxGIiKStzu4ozsGQ7MTM+ffpVZ1zTT33zewnwDPewzpgctymk7w2ztAuIiKe9lAYIDlH/mZWEPfwFqBnJdBqYLmZ5ZjZFKAEWAesB0rMbIqZZRM7KLy6LzWIiKSj9lAESFz4n/PI38xWAdcBeWZWC3wVuM7M5gAO2Av8PYBzbpuZ/YbYgdww8DnnXMR7nXuA54FM4GHn3Lb+6oyISLro6O4Jf5+nfZxzd/bS/NMzbH8/cH8v7RVAxbm+r4hIECV65K9v+IqIJKH2rticf6LW+Sv8RUSSkEb+IiIB1J7gOX+Fv4hIEupI5qWeIiKSGJr2EREJoJ7w1wFfEZEAaQ+FycwwsjMTE9MKfxGRJNQeipCblYmZJeT1Ff4iIkmoIxQhNycxUz6g8BcRSUqx6/cmZpknKPxFRJJSewKv3wsKfxGRpNQeCidsmSco/EVEklIir98LCn8RkaTUEYowVHP+IiLB0t6taR8RkUBxznGkvZuhORr5i4gExjstbRztDHNZ4YiEvYfCX0Qkyazb0wrA/CljE/YeCn8RkSSzfm8recNyKB6bm7D3UPiLiCSZdXtamT9ldMLO6wMKfxGRpPJOy3HqDncwv3hMQt9H4S8ikkQe31BLZoaxtLQgoe+j8BcRSRLhSJQnNtRy/Yx8xo0YnND3UviLiCSJF7c30Xysi4/Nm5zw91L4i4gkgUjU8f21u5maP5QbLhmX8PdT+IuIJIHVm+vY3XScf75hOoMSdOnGeAp/ERGfHW4Pcf+aHcyeNJLyBB/o7ZG4E0eIiMhZOef4+tPbOdTezaOfXkBGRuLW9sfTyF9ExEe/razlqbfq+IcPXsysiYk7l8+pFP4iIj55vfoAX/ndVt5/8Vj+4YMlA/reCn8RER9sqT3CZx6tZEreUH708XlkDtB0Tw+Fv4jIANvVeIxP/mwdo3KzeeTT8xmZmzXgNSj8RUQG0Js1B7n9P14nM8P4xd3zmTAysd/kPR2Fv4jIAFlT1cDfPbyOvOE5PPnf38fU/GG+1aKlniIiCdYdifKd53bykz/uYW7RKH5615WMHprta00KfxGRBNrf2s7nf7OJ9XsPcdfVF/Hl8llkD/J/0kXhLyKSAM45fltZy33PbAfgweVzWDan0OeqTlL4i4j0s5qW43x19Tb++PYBFkwZw3dvv5zJYxJ3ScYLcc5/e5jZw2bWbGZb49rGmNmLZva293O0125m9gMzqzazKjObG/dv7vK2f9vM7urf7oiI+KezO8L3XtjF4v/zRzbtP8x9yy5l1WeuSrrgh/Nb7fNzYPEpbfcCLznnSoCXvMcAS4AS77YC+DHEflkAXwUWAPOBr/b8whARSVXRqON3b9Vxw/f+wA9erqZ8dgEv/c9r+cTVxQN2rp7zdc7TPs6518ys+JTmZcB13v1HgFeBL3rtjzrnHPCGmY0yswJv2xedc60AZvYisV8oqy68CyIi/nDO8equFr793E52Nh5jVsEIfvWZ2bxvWp7fpZ1VX+f8xzvnGrz7jcB4734hsD9uu1qv7XTt72FmK4j91UBRUVEfyxQR6V8b3m3lO8/t4s09rRSNyeXB5XP4yOyJSTvSP1W/HfB1zjkzc/34eiuBlQBlZWX99roiIhfKOcefqw/yw1fe5o2aVvKGZXPfsktZfmVRUizfPB99Df8mMytwzjV40zrNXnsdEH8RykleWx0np4l62l/tYw0iIgkVjTrW7mjioVffYfP+w4wfkcNXyi/hzvlFDM1JzUWTfa16NXAX8C3v5+/j2u8xs8eIHdw94v2CeB74ZtxB3g8BX+pjDSIiCdEeCvPExjp+9uc91LS0UTQml2/eUspt8wrJGZTpd3l9cs7hb2ariI3a88ysltiqnW8BvzGzu4F3gTu8zSuApUA10A58CsA512pm3wDWe9vd13PwV0QkWdQd7uDR1/eyat0+jnaGmT1pJA8un0N5acGAXF93IFhsQU5yKysrc5WVlX6XISJpzDnHm3taefQve3l+WxMAiy+dwKevKWZu0WjMUuNAbjwz2+CcK+vtudScrBIR6SetbSGe2FDLqnX7qDnQxsghWfzXD0zhE1cXUzhqiN/lJYzCX0QCp2eUv2rdPp7d0kgoEmXeRaN54PqLWVpawJDs1J7PPxcKfxEJjOajnfxuUx2Prd9PTUsbwwcP4uMLirhzfhEzJgz3u7wBpfAXkbTW2R1h7Y4mHt9Qy2u7W4g6mFs0iu/efjnlARnl90bhLyJpxznHxn2HeWJjLc9srudoZ5iCkYP5b9dN49a5k5jm4xW0koXCX0TSRv3hDp56q44nNtRSc6CNwVkZLLmsgNvmTuLqaWPJTJFTLwwEhb+IpLRDbSEqtjawelM96/a24hzMnzKGz147jSWlExg+OMvvEpOSwl9EUk5bV5i1O5pYvameP+xuIRx1TM0fyj8tms4tVxRSNDb5zp+fbBT+IpISQuEor+1u4feb61m7vYmO7ggFIwdz9zVT+MjlE7l04oiU/CKWXxT+IpK0IlHHuj2trN5cR8WWRo50dDM6N4tb5xZy8+UTubJ4TMqcQjnZKPxFJKk459hce4RnNtfzdFU9TUe7yM3O5EOzxrNsTiHXlOSRlSbn1/GTwl9EfNcT+Guq6qnY0kjd4Q6yMo1rp4/jK+UTueGS8YFdj58oCn8R8cXpAv8DJfn8843TuXHWeEYO0UqdRFH4i8iAcc6xaf9hKrY0KPB9pvAXkYRS4Ccnhb+I9DsFfvJT+ItIvzhT4H/+xuncoMBPKgp/EblgCvzUpfAXkfPSE/hrqhp4duvJwF+owE8pCn8ROSsFfvpR+ItIr3oL/OzMDD5QkqfATwMKfxE5wTnHW/sPU6HAT3sKf5GAi0a9wN/SwLNbGqg/0qnADwCFv0gAxQL/EGuqGnl2awMNXuAvnJ7Hv9w0gxtmjWeELoKS1hT+IgERjTo27jvEmi0NPLe1MS7w8/nC4hksukSBHyQKf5E01hP4z1TFAr/xaCfZgzK4dno+X1w8k0WXjNNlDgNK4S+SZqJRx4Z9h1jTS+DfW6rAlxiFv0gaiEYdle8eih203dpA09EusgdlcN30fL40eyYfnKnAl7+m8BdJUZGoo3Jvqxf4jTQf6yJnUAbXzchnaWkBiy4Zz7Ac/S8uvdMnQySFRKKO9XGB3+IF/vUzxrF0dgEfnDlOgS/nRJ8SkSTXE/hrqhp4btt7A3/RzHEMVeDLedInRiQJRaKOdXtOjvAPHO9icJYX+KWxEb4CX/pCnx6RJBGJOt7cc5CKLQ08t7XpROB/cGYs8K+focCX/qNPkoiPwpEo6/a0smZLA89va+TA8dCJwC8vncj1M/PJzdb/ptL/9KkSGWDhSJQ3ewJ/ayMH20IMyco8OcJX4MsA0CdMZACEI1HeqIkF/gvb4gL/knGUlxZw3QwFvgwsfdpEEuRk4Nfz/LYmWttC5GZnelM6BVw3YxxDsjP9LlMCql/C38z2AseACBB2zpWZ2Rjg10AxsBe4wzl3yMwMeBBYCrQDn3TObeyPOkT8Fo5E+UtN7KDtqYH/4dkFXDtdgS/JoT9H/tc75w7EPb4XeMk59y0zu9d7/EVgCVDi3RYAP/Z+iqSk7kiUv7zTE/iNHGrvJjc7k0WXjKe8dALXzRjH4CwFviSXRE77LAOu8+4/ArxKLPyXAY865xzwhpmNMrMC51xDAmsR6VfdkSivv3OQiqoGnt/eyOH2boZ6gb/Um8NX4Esy66/wd8ALZuaA/+ucWwmMjwv0RmC8d78Q2B/3b2u9tr8KfzNbAawAKCoq6qcyRS5cdyTKn6sPULGlgRe2N50I/BtmxQL/2ukKfEkd/RX+1zjn6sxsHPCime2Mf9I557xfDOfM+wWyEqCsrOy8/q1If+kJ/DVVscA/0tHNsJxB3HBJbFnmQgW+pKh+CX/nXJ33s9nMngLmA0090zlmVgA0e5vXAZPj/vkkr00kKYTCUf78zgEq4gJ/eM6gEyP8D5TkKfAl5fU5/M1sKJDhnDvm3f8QcB+wGrgL+Jb38/feP1kN3GNmjxE70HtE8/3it1DYG+F76/CPdoZPBH55aQEfmJ5HziAFvqSP/hj5jweeiq3gZBDwK+fcc2a2HviNmd0NvAvc4W1fQWyZZzWxpZ6f6ocaRM5bKBzlT9UtrKlq5MXtJwP/xp4RvgJf0lifw985VwNc3kv7QWBRL+0O+Fxf31fkQnSFI/zp7dgI/8XtTRzrDDN8cCzwy0sLuKZEgS/BoG/4Sto7EfhVDby442Tgf2jWBMpnT+D9FyvwJXgU/pKWusIR/rg7tizzxe1NHOsKM2LwIG66dALlpQW8/+I8sgdl+F2miG8U/pI2Orsj/PHtWOCv9QJ/5JAsFl82gaWzC3j/NAW+SA+Fv6S0zu4Ir+1uiQX+jmaOe4G/pHQCS0sLeJ8CX6RXCn9JOZ3dEf7gBf5LXuCPys1iaekEymdP5H3TxpKVqcAXOROFv6SEzu4Ir+7qCfwm2kIRRuVmUV5awNLZBQp8kfOk8Jek1RP4a7Y08LIX+KNzs/jI5RNZWlrA1Qp8kQum8JekEgv8ZtZsafyrwL95Tizwr5qqwBfpDwp/8V1HqCfwG3h5ZzPtoQhjhmZz85xCyksLuGrqGAYp8EX6lcJffNERivCKF/iveIE/dmg2H70iFvgLpijwRRJJ4S8Dpj0U5pWdsYO2L+9spqM7Qt6wbG7xAn++Al9kwCj8JaF6An/Nlnpe2dlyIvBvnVtI+ewCFkwZS2aG+V2mSOAo/KXftYfCvLyz+cQIv7M7St6wHG6bV8jSUgW+SDJQ+Eu/aOs6Gfiv7DoZ+LfPm8xSb0pHgS+SPBT+csHausK8tLOZiqoGXt0dC/z84TncURYL/CuLFfgiyUrhL+elZ4S/pio2wu8KRxk3PIe/8QK/TIEvkhIU/nJWvU3pjBuew/IrFfgiqUrhL73q7aBtvkb4ImlD4S8nxC/LjA/8O8omU67AF0krCv+AO/FN26r4L17poK1IulP4B1D8qRVe3nHym7YfmzdJyzJFAkLhHxCd3RFe2RkL/JfiAl9fvBIJJoV/GuuORPlT9QGe3lTPC9ubON4VZuxQnVpBRBT+aScadVS+e4jVm+uo2NJIa1uIEYMHUV5awM1zJupsmSICKPzTgnOObfVHWb25nmc211N/pJPBWRncOGsCN18+kYXT88gZlOl3mSKSRBT+KeydluOs3lTP05vrqTnQxqAM49rp+XxxyUxuuGQ8Q3O0e0Wkd0qHFNNyrIvfb6rjd5vq2Fp3FDO4aspYPrNwKksum8Co3Gy/SxSRFKDwTwGd3RHW7mjiyY11/GF3C5GoY/akkfyvD8/iw7MLGD9isN8likiKUfgnKeccG/cd4vENdTxTVc+xzjATRgxmxcKp3HpFISXjh/tdooikMIV/ktnf2s5Tb9Xx5MZa9h5sZ0hWJosvm8Btcydx9TQtzRSR/qHwTwKhcJS1O5pYtW4ff6o+gHNw9dSxfO76i1lSWsAwHbgVkX6mVPFRTctxfr1+P49vqOVgW4iJIwfzj4tK+Ni8SUwanet3eSKSxhT+AywUjvLs1gZWrdvHGzWtZGYYi2aO484FRSwsyde0jogMCIX/AGk+2skv39zHr9bto+VYF0VjcvnXm2Zw+7xJjNNqHREZYAr/BHtr3yF+/vpeKrY00B1xXD8jn7veV8zCknwyNMoXEZ8o/BMgEnW8sK2R/3iths37DzMsZxD/5aqL+MTVxUzJG+p3eSIiCv/+1BWO8NTGOla+VkPNgTYuGpvL12++lNvmTdKKHRFJKr4lkpktBh4EMoH/55z7ll+19FVnd4T/fONdVr5WQ/OxLi4rHMEPP34FSy4r0AFcEUlKvoS/mWUCDwE3ArXAejNb7Zzb7kc9FyoUjvLr9fv495eraT7WxdVTx/LAHZdzzcV5mCn0RSR5+TXynw9UO+dqAMzsMWAZkBLhH45EeeqtOh586W1qD3VwZfFo/v3OK1gwdazfpYmInBO/wr8Q2B/3uBZYEL+Bma0AVgAUFRUNXGVn8Xr1Ab7+9HZ2NR2jtHAk//bRy7h2er5G+iKSUpL2KKRzbiWwEqCsrMz5XA61h9r5ZsUOKrY0Mmn0EH70t3NZctkEhb6IpCS/wr8OmBz3eJLXlnQiUcfP/ryH776wC4DP3zidFQunMjhLV8YSkdTlV/ivB0rMbAqx0F8OfNynWk5rz4E2/vW3m6l89xCLZo7jGx+9jImjhvhdlohIn/kS/s65sJndAzxPbKnnw865bX7U0hvnHP/5xrvcX7GD7MwMHrj9cm6dW6gpHhFJG77N+TvnKoAKv97/dNq6wtz75Bae3lzPtdPz+fZts5kwUufeEZH0krQHfP2wv7WdT/98Pe+0HOcLi2fw2YXTdP4dEUlLCn/P1rojfOrn6wmFo/zi7gW8/+I8v0sSEUkYhT/wZs1B7n6kkhGDB/Grz16t6+OKSNoLfPhvrTvC3Y9UMmHkYH5x93wKRmo1j4ikvwy/C/DT3gNtfPJn6xg5JEvBLyKBEtjwP9LRzV0/W0fUwaMKfhEJmECGv3OOLzy+mbpDHfzkE/OYlj/M75JERAZUIMP/yY11PL+tiS8unsm8i8b4XY6IyIALXPgfPN7FN9ZsZ95Fo7n7mil+lyMi4ovAhf8PX6nmWGeY/31rqb7AJSKBFajwrz/cwS/f2MdtcwuZrrX8IhJggQr/h/+0B4fjfywq8bsUERFfBSb8Q+EoT75Vx42zxjNpdK7f5YiI+Cow4f/SjiZa20LcXjb57BuLiKS5wIT/U2/VMX5EDgtL8v0uRUTEd4EI/2jU8eaeVq6dnk+mVviIiAQj/N9uPs6Rjm6uLNYXukREICDhv27PQQAWTBnrcyUiIskhGOG/9xDjR+QweYxO3iYiAgEJ/037D1F20RhdgF1ExBOI8D/S3k3esGy/yxARSRqBCP+O7ghDsgN/0TIRkRPSPvy7I1G6I46h2Zl+lyIikjTSPvzbQxEAhij8RUROSPvw7/DCP1fTPiIiJ6R9+LeHwgDkauQvInJCAMJf0z4iIqcKTPgP1bSPiMgJAQj/2LSPRv4iIielffifPOCr8BcR6ZH24d+u8BcReY8AhH/Pah/N+YuI9AhA+GvkLyJyqsCE/5Ashb+ISI+0D/+O7giDszLI0OUbRUROSPvwbw+FNd8vInKKAIR/RPP9IiKnSP/w71L4i4icqk/hb2ZfM7M6M9vk3ZbGPfclM6s2s11mdlNc+2KvrdrM7u3L+5+Ldl3IRUTkPfojFb/vnPtufIOZzQKWA5cCE4G1Zjbde/oh4EagFlhvZqudc9v7oY5edYTC5Gqlj4jIX0nUtM8y4DHnXJdzbg9QDcz3btXOuRrnXAh4zNs2YTTnLyLyXv0R/veYWZWZPWxmo722QmB/3Da1Xtvp2t/DzFaYWaWZVba0tFxwcR2hCLk5mvYREYl31vA3s7VmtrWX2zLgx8A0YA7QADzQX4U551Y658qcc2X5+fkX/DptmvYREXmPsw6JnXM3nMsLmdlPgGe8h3XA5LinJ3ltnKE9IdpDEZ3OWUTkFH1d7VMQ9/AWYKt3fzWw3MxyzGwKUAKsA9YDJWY2xcyyiR0UXt2XGs6mQ3P+IiLv0dfJ8O+Y2RzAAXuBvwdwzm0zs98A24Ew8DnnXATAzO4BngcygYedc9v6WMNphcJRwlGn8BcROUWfwt8593dneO5+4P5e2iuAir6877k6eSEXHfAVEYmX9t/wLZ9dwLRxw/wuQ0QkqaT1kHhkbhYPfXyu32WIiCSdtB/5i4jIeyn8RUQCSOEvIhJACn8RkQBS+IuIBJDCX0QkgBT+IiIBpPAXEQkgc875XcNZmVkL8G4fXiIPONBP5SS7IPUV1N90FqS+QmL6e5Fzrtdz4qdE+PeVmVU658r8rmMgBKmvoP6msyD1FQa+v5r2EREJIIW/iEgABSX8V/pdwAAKUl9B/U1nQeorDHB/AzHnLyIify0oI38REYmj8BcRCaC0Dn8zW2xmu8ys2szu9bueRDCzvWa2xcw2mVml1zbGzF40s7e9n6P9rvNCmdnDZtZsZlvj2nrtn8X8wNvfVWaWUlfyOU1fv2Zmdd7+3WRmS+Oe+5LX111mdpM/VV84M5tsZq+Y2XYz22Zm/+i1p93+PUNf/du/zrm0vBG7QPw7wFQgG9gMzPK7rgT0cy+Qd0rbd4B7vfv3At/2u84+9G8hMBfYerb+AUuBZwEDrgLe9Lv+fujr14B/6WXbWd5nOgeY4n3WM/3uw3n2twCY690fDuz2+pV2+/cMffVt/6bzyH8+UO2cq3HOhYDHgGU+1zRQlgGPePcfAT7qXyl945x7DWg9pfl0/VsGPOpi3gBGmVnBgBTaD07T19NZBjzmnOtyzu0Bqol95lOGc67BObfRu38M2AEUkob79wx9PZ2E7990Dv9CYH/c41rO/B87VTngBTPbYGYrvLbxzrkG734jMN6f0hLmdP1L131+jzfN8XDcFF5a9dXMioErgDdJ8/17Sl/Bp/2bzuEfFNc45+YCS4DPmdnC+Cdd7G/ItF3Pm+79A34MTAPmAA3AA75WkwBmNgx4Avgn59zR+OfSbf/20lff9m86h38dMDnu8SSvLa045+q8n83AU8T+NGzq+XPY+9nsX4UJcbr+pd0+d841Oecizrko8BNO/umfFn01syxiYfhL59yTXnNa7t/e+urn/k3n8F8PlJjZFDPLBpYDq32uqV+Z2VAzG95zH/gQsJVYP+/yNrsL+L0/FSbM6fq3GviEtyrkKuBI3PRBSjplTvsWYvsXYn1dbmY5ZjYFKAHWDXR9fWFmBvwU2OGc+17cU2m3f0/XV1/3r99HwRN5I7Y6YDexI+Vf9rueBPRvKrEVAZuBbT19BMYCLwFvA2uBMX7X2oc+riL253A3sXnPu0/XP2KrQB7y9vcWoMzv+vuhr7/w+lLlBUJB3PZf9vq6C1jid/0X0N9riE3pVAGbvNvSdNy/Z+irb/tXp3cQEQmgdJ72ERGR01D4i4gEkMJfRCSAFP4iIgGk8BcRCSCFv4hIACn8RUQC6P8DtWpNAFMkQagAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(sigmoid_R)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "111ba97d",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAz80lEQVR4nO3dd3hU1fb/8fdKrxBI6BASegchgIAFGyIqiAUVULF7r71dGyo/b7FhwQuigFwQESygRqSpF8WrdJQWqtRQQyAJBELKrN8fM/CNCGGATE4ys17PM09mzjkz89kZmJWzzzl7i6pijDEmcAU5HcAYY4yzrBAYY0yAs0JgjDEBzgqBMcYEOCsExhgT4EKcDnC6EhISNCkpyekYxhhToSxZsmSvqlY70boKVwiSkpJYvHix0zGMMaZCEZEtJ1tnXUPGGBPgrBAYY0yAs0JgjDEBrsIdIziRgoIC0tPTycvLczpKuRYREUHdunUJDQ11Oooxphzxi0KQnp5ObGwsSUlJiIjTccolVSUzM5P09HSSk5OdjmOMKUd81jUkImNFZI+IrDzJehGRd0Rkg4gsF5H2Z/peeXl5xMfHWxEogYgQHx9ve03GmD/x5TGCcUDPEtZfATT23O4BRp7Nm1kRODX7HRljTsRnXUOqOldEkkrYpA/wobrHwZ4vInEiUktVd/oqkzHGlDsuF+QfgMP74XCW+2deFuQfwlWQx/6cHDKzctiffYCEDr1p2O6CUo/g5DGCOsC2Yo/TPcv+VAhE5B7cew0kJiaWSbjScNddd/HYY4/RokULn71Hr169+Pjjj4mLi/vD8iFDhhATE8MTTzzhs/c2xnjBVQTZ2yDzd9i30X3L/B32b4LcDMjLBnWd8KlBQLznBjA/tobfFQKvqeooYBRASkpKhZlJZ8yYMT5/j+nTp/v8PYwxXlCF/ZshfTHsWu7+ss/c4P7CL8r/v+1CoyiqkkxWZDK7wtuy7XA4Gw+GsSk3hGyNJktjkKg46lRPoF71eJJqVqVRrXga1orn3HDfnPHnZCHYDtQr9riuZ1mFlJubS79+/UhPT6eoqIjnn3+ekSNHMnToUFJSUvjggw949dVXiYuLo23btoSHhzN8+HAGDRpEZGQkv/76K3v27GHs2LF8+OGHzJs3j86dOzNu3DgAJk2axL/+9S9UlSuvvJJXX30V+L8hNxISEvjnP//J+PHjqV69OvXq1aNDhw4O/kaM8WOqkLXF/aW/8zfYucx9y8t2rw8Og6oNIaEx2qQnGeH1SMtLYF52HD9uD2LdtoO4PH/S1omLpGXdSrSqU5metd0/q8eGl+kxPScLQSrwgIhMBjoD2aVxfOD/fb2KtB05Zx2uuBa1K/Hi1S1L3GbmzJnUrl2bb775BoDs7GxGjnQf/96xYwd///vfWbp0KbGxsVx88cW0bdv22HP379/PvHnzSE1NpXfv3vz888+MGTOGjh078ttvv1G9enWeeuoplixZQpUqVejRowdffvkl11xzzbHXWLJkCZMnT+a3336jsLCQ9u3bWyEwprS4XLB3LWz5GbbMg63zIMfzd2twGNRoCS37Qq22FNU6h7SiROZtzmL+xn0snb+frEMFAMSG59MuMY7LW9Wiff0qtKlTmSrRYQ42zM1nhUBEJgHdgQQRSQdeBEIBVPU9YDrQC9gAHAJu91WWstC6dWsef/xxnnrqKa666irOP//8Y+sWLlzIhRdeSNWqVQG44YYbWLdu3bH1V199NSJC69atqVGjBq1btwagZcuWbN68mS1bttC9e3eqVXMPHDhgwADmzp37h0Lw008/0bdvX6KiogDo3bu3r5tsjH8rPAIbvodVX8CG7+DwPvfymJpQvwskdoXEzrgSmrEu8wjzfs/kl7RMFkzbQ07eDgAaJETTo0UN2idWoX39KjSqFkNQUPk7e8+XZw3dfIr1Ctxf2u97qr/cfaVJkyYsXbqU6dOnM3jwYC655BKvnxseHg5AUFDQsftHHxcWFtqVwMaUlcJ82DjH/eW/5hs4kgMRcdD0Ckg6310AqiSTdbiAH9dlMOfHPfy0fi6Zue5jAIlVo7iiVS26NIynS8N4alSKcLY9XqoQB4srgh07dlC1alUGDhxIXFzcHw4Ud+zYkUceeYT9+/cTGxvLlClTjv3V741OnTrx0EMPsXfvXqpUqcKkSZN48MEH/7DNBRdcwKBBg3jmmWcoLCzk66+/5t577y219hnjt1xFsPEHWDkF1kxz9/NHVIbmvd3dPQ0uRINCWL3zAHOW7WHOmnks3bofl0LV6DAubFKNrp4v/rpVopxuzRmxQlBKVqxYwZNPPklQUBChoaGMHDny2KmbderU4dlnn6VTp05UrVqVZs2aUblyZa9fu1atWrzyyitcdNFFxw4W9+nT5w/btG/fnhtvvJG2bdtSvXp1OnbsWKrtM8bvHM6C/70FyybBwd0QXhmaXen58u9OUVAoCzftY8a0tXybtpud2e6r8lvVqcQDFzXiombVaVM3juBy2NVzusTdQ1NxpKSk6PET06xevZrmzZs7lMg7Bw8eJCYmhsLCQvr27csdd9xB3759yzxHRfhdGeNTRw7CrxPgpzfh0F5ocgW0vREaX05hUBjzN+5j+sqdzF61i70H84kIDeLCJtW4pFkNujetRvUK0t1zPBFZoqopJ1pnewRlZMiQIXz33Xfk5eXRo0ePPxzoNcaUge1LYcl/YOVUyD8I9c+Dyz9Da7Vl6db9TJ22nukrdrL/UAFRYcFc3Kw6vVrXonvTakSF+fdXpX+3rhwZOnSo0xGMCTxFhbDma5g/ErYtgNAoaHktdBjE1qiWTP01nS8m/sCWzENEhAZxWYuaXOn58o8IDXY6fZmxQmCM8T+H98PSD2HBKMhJhypJ0PMVDre4ia/XHeTTadtYvGUOItClQTwPXNSIK1rXIiY8ML8SA7PVxhj/lL0dfh7mPgZQcMh9ymev11lXuSsfL9rOlJkLOZBXSMNq0fytZ1OuaVeH2nGRTqd2nBUCY0zFl53uPgNo6YfuAdza3Eh+x3uZkZHAxB+2snDzz4QFB3FF65oM6FyfjklVbFj2YqwQGGMqrux099k/v05wF4BzBpLd4UEmrFHGj9tCxoF0kuKjeLZXM67vUI+q5WA4h/LICkEp6tq1K7/88stpP6979+7HBqczxnghaxv8701YOsH9+JyBbGv5F0YtL+CzkevJK3BxQZNqvHFDMuc1SiiXwzqUJ1YIStGZFAFjzGnIz3V3Af38jnsPoP0trGxwJ8MW5/Hd6A2EBgVxzTm1ufO8BjStGet02grDl1NVBpyYmBgAXn/9dTp27EibNm148cUXAdi8eTPNmzfn7rvvpmXLlvTo0YPDhw//4fkul4tBgwYxePBgioqKGDRoEK1ataJ169a89dZbZd4eY8oNVfcQEMM7wtzXoUVvll83h1t238RVH25h8eZ9PHhRI/739EW8dn1bKwKnyf/2CGY8DbtWlO5r1mwNV7zi1aazZ89m/fr1LFy4EFWld+/ezJ07l8TERNavX8+kSZMYPXo0/fr1Y8qUKQwcOBCAwsJCBgwYQKtWrXjuuedYsmQJ27dvZ+XKlQBkZWWVbpuMqSh2LoOZz7iHgK7ZhlVd3uRfq+L4ecJW4qPDeOaKZgw8tz7RAXrqZ2mw31wpmz17NrNnz+acc84B3ENLrF+/nsTERJKTk2nXrh0AHTp0YPPmzceed++999KvXz+ee+45ABo0aMDGjRt58MEHufLKK+nRo0dZN8UYZx3cA9+/BL9+BFHxbO7yT57e1I75X2WTEHOQwVc2p3/nRL+/6rcs+N9v0Mu/3H1FVXnmmWf+NPLn5s2b/zDEdHBw8B+6hrp27cqcOXN4/PHHiYiIoEqVKixbtoxZs2bx3nvv8emnnzJ27Ngya4cxjik8Agvegx9fh8LD7G97Dy9m9yJ1Ti7VYvN44aoW9O+cGFBX/vqa/xUCh11++eU8//zzDBgwgJiYGLZv3+7VfAJ33nknc+fOpV+/fkydOpWsrCzCwsK47rrraNq06bEuJGP8liqsnQ6znoP9mzjS4DKGhw5ixEIhKuwIT/Rowh3nJdsegA/Yb7QUiQg9evRg9erVdOnSBXAfQP7oo48IDj71Xy+PPfYY2dnZ3HLLLTz99NPcfvvtuFwuAF5++WWfZjfGUbvTYObTsOlHihKaMrX5Ozy/sjpFLuXWLvV58OJGxMeEn/p1zBmxYahLSWZmJu3bt2fLli2O5jiV8vC7MuaYwiPus4D+9xYaFsOKJvdzb1pbdh4spHfb2jzRoymJ8RVzspfyxoah9rEdO3bQvXv3YxPRGGO8sHsVfHY77F1LVpPreSyrH/9dWEibutGMuLUl7ROrOJ0wYFghKAW1a9f+w2T0xphTWD0Npt6DKyyGjxq+yZAVNakcKbx8bWtuTKlnVwKXMb8pBKpqg0idQkXrBjR+SBXmDoU5/2BfXGtuznmQ9WkxDOhcn8d7NCEuysYCcoJfFIKIiAgyMzOJj4+3YnASqkpmZiYRERVzmj3jB/KyIfUhSPuSX6Iv4fZdt9C8XnVSr2lFqzrez+FtSp9fFIK6deuSnp5ORkaG01HKtYiICOrWret0DBOIdi5DP7kFzU7nLVd/xuT05qmrmnFb1yS/mPy9ovOLQhAaGkpycrLTMYwxJ7JsMq7Uh9mnMdyT9zzRjboyu29r6lW1s4HKC78oBMaYcqgwH531LLJoNAtdLXgm6DHuv74L17WvY1245YwVAmNM6cvN5MjH/QnfPp/Rhb1Y3PgRPr22HdVi7aKw8sgKgTGmdO1ZTe646wk5tJsnXQ+S0uce3kupZ3sB5ZgVAmNMqcldOYPgqXeQWxTG0Kqvcf+AG0lKiHY6ljkFKwTGmLOnytYZb1Jn4T9Y46rHgnNH8K/LuxESbHNfVQRWCIwxZ0WLCkgb+xdabv+Mn4I7UfnWcdzRoI7TscxpsEJgjDlj2fv2sm30jbQ6vJhZcTfR5d53qBRpB4QrGp/ut4lITxFZKyIbROTpE6xPFJE5IvKriCwXkV6+zGOMKT2r1qwh89/daXroV/7X4kV6PPyeFYEKymd7BCISDIwALgPSgUUikqqqacU2Gwx8qqojRaQFMB1I8lUmY8zZU1WmfvsjnX++kyqSy5ZeEziv85VOxzJnwZddQ52ADaq6EUBEJgN9gOKFQIFKnvuVgR0+zGOMOUt5BUWM+HgKt218nLDgIFwDv6ZRg45OxzJnyZeFoA6wrdjjdKDzcdsMAWaLyINANHDpiV5IRO4B7gFITEws9aDGmFPbmX2YYWPH82zWEIioTMxd0wiq1tjpWKYUOH1u183AOFWtC/QCJojInzKp6ihVTVHVlGrVqpV5SGMC3ZIt+3h52L95Met5givVpNJfv7ci4Ed8uUewHahX7HFdz7Li7gR6AqjqPBGJABKAPT7MZYw5DZ8s2srPX43mzZARFFZrTvSgLyHG/iDzJ77cI1gENBaRZBEJA24CUo/bZitwCYCINAciABtL2phyoLDIxYtfrWT+F+/ydshwqJtCxF3TrQj4IZ/tEahqoYg8AMwCgoGxqrpKRF4CFqtqKvA4MFpEHsV94HiQ2jRaxjju4JFC7p+4lOq/f8abYaPRpPMJ6T8Zwmy4CH/k0wvKVHU67lNCiy97odj9NKCbLzMYY07PzuzD3DFuMW0yUnk1dBQ0uhS58SMIjXQ6mvERpw8WG2PKkbQdOfQd8Qtt983klZDR0PASuHGiFQE/Z4XAGAPAD2v3cMN7v3C560deDnoXST4fbpoIoTbPtb+zsYaMMUxauJXBX67ghUozuC1vAiSdDzdPtj2BAOF1IRCRKFU95MswxpiypaoM+349w79bzX/iP+aC3JnQuh/0GQ4hNm5QoDhl15CIdBWRNGCN53FbEXnX58mMMT7lcilDUlcx9rvf+Cb+HXcRuOBJuHaUFYEA480ewVvA5XiuAVDVZSJygU9TGWN8Kr/QxROfLWPl8sXMiXuHqod3QO/h0P4Wp6MZB3jVNaSq246bb7TIN3GMMb52KL+Qv3y0FNeG75kRNZzwoAi47Wuo39XpaMYh3hSCbSLSFVARCQUeBlb7NpYxxheyDuVzx7hF1N8+jTfCRxGU0Mx9UDiu3qmfbPyWN6eP3gfcj3s00e1AO89jY0wFsicnj37vz6P9zsm8FfouQfW7wO3TrQiYU+8RqOpeYEAZZDHG+MjO7MP0H72A7jmpDA7+EJpfDdeOsWsEDFBCIRCRf+Me/+eEVPUhnyQyxpSqbfsO0X/MfLrkzuGFoLHQ5Aq4/j8QHOp0NFNOlLRHsLjMUhhjfGLT3lwGjJ5P2yOLeTXoXSSxG9xgRcD80UkLgaqOL8sgxpjStWHPAfqPXkCjot8ZEfI2ktAcbp5kVwubPympa+htVX1ERL7mBF1Eqtrbp8mMMWds9c4cBo5ZQB3J4MPw1wkKi4f+n0FEpVM/2QSckrqGJnh+Di2LIMaY0pG2I4f+Y+ZTPfgQn0cPJSQvHwZMg0q1nI5myqmSuoaWeO62U9VhxdeJyMPAj74MZow5fWt3HWDgBwuICykkNX44YXu2wS1fQvVmTkcz5Zg31xHcdoJlg0o5hzHmLK3ffYD+o+cTEeTim3ofEbFzkXvcoCSb+8mUrKRjBDcD/YFkESk+13AssM/XwYwx3vs94yA3j15AqBQxu/5HRG/4Bi7/F7Ts63Q0UwGUdIzgF2AnkAC8UWz5AWC5L0MZY7y3eW+ue09ADzOrzgdEb5gDl70EXWwAAOOdko4RbAG2AF3KLo4x5nRszTzEzaPnE1mYw4xqw4nc9itc9Tak3O50NFOBeDMfwbUisl5EskUkR0QOiEhOWYQzxpzctn3uIhCbv4eZlV8mcu8KuGGcFQFz2rwZffQ14GpVtRFHjSkndmYfpv+Y+cTnbeHzmNcJy82GAZ9DgwudjmYqIG8KwW4rAsaUH5kHjzBwzAJq565hYuTrhLiCYNA0qH2O09FMBeVNIVgsIp8AXwJHji5U1am+CmWMObGcvAJu+89CErMWMCZ8GMER8e7rBOIbOh3NVGDeFIJKwCGgR7FlClghMKYMHc4v4q5xi2myewZDQ98nqGpTGDjFrhg2Z82b+QjsyJMxDssvdPGXiUtomz6B50ImQv3z4caPIDLO6WjGD5yyEIhIBHAn0BI4NouFqt7hw1zGGI8il/LoJ79y7u/DuC9kmvsisb7vQ0i409GMn/BmiIkJQE3gctzjC9XFfVGZMcbHVJXBXyyn9eo33UWg491w3VgrAqZUeVMIGqnq80CuZ46CK4HOvo1ljFFVXp6+mtq/Hi0Cd0Gv1yHIm/+2xnjPm39RBZ6fWSLSCqgMVPddJGMMwLs//E7kL0N5MORLtP1tcMXrIOJ0LOOHvCkEo0SkCvA8kAqk4b7I7JREpKeIrBWRDSLy9Em26SciaSKySkQ+9jq5MX5swrzNHPnuZR4NnYK27Y9c9bbtCRif8easoTGeuz8CDbx9YREJBkYAlwHpwCIRSVXVtGLbNAaeAbqp6n4RsT0NE/BmrthJzjcv8FjoV7ja3ERQn+FWBIxPeXPW0AsnWq6qL53iqZ2ADaq60fM6k4E+uPcojrobGKGq+z2vuceb0Mb4qyWbM0n/7EnuD/mawna3EtJ7mBUB43Pe/AvLLXYrAq4Akrx4Xh1gW7HH6Z5lxTUBmojIzyIyX0R6nuiFROQeEVksIoszMjK8eGtjKp7fMw6ycPzT3BX0NXntbiekzztWBEyZ8KZrqPhcBIjIUGBWKb5/Y6A77tNS54pIa1XNOi7DKGAUQEpKipbSextTbuzJyWP6+8/xoH7Kweb9iOn9ph0YNmXmTP7ciML9pX0q24F6xR7X9SwrLh1IVdUCVd0ErMNdGIwJGAePFDL5vb/zYOE4spJ7EXP9SNsTMGXKm/kIVojIcs9tFbAWeNuL114ENBaRZBEJA27CfdZRcV/i3htARBJwdxVt9Dq9MRVcQZGL8aOG8kDuCDJrXUDcgPEQ7M0QYMaUHm/+xV1V7H4h7mGpC0/1JFUtFJEHcHcjBQNjVXWViLwELFbVVM+6HiKShvv4w5OqmnnarTCmAlJVxo8byT17X2NvfAeq3/4JhIQ5HcsEIFEtuctdRKqWtF5Vy3Qi+5SUFF28eHFZvqUxPvHZJ+PpnfYYWbFNqPHALIio5HQk48dEZImqppxonTd7BEtx9/XvBwSIA7Z61imncW2BMcbt228+5+q0x8mMTKbWX7+xImAc5c0RqW9xT1WZoKrxuLuKZqtqsqpaETDmNM37eQ5dFt7P3rA6VP/rdCSqxJ1uY3zOm0JwrqpOP/pAVWcAXX0XyRj/tTwtjQazbycvOIb4+6YRUskupjfO86ZraIeIDAY+8jweAOzwXSRj/NOm7TuJ+PQmYiSP/FtmEBlf79RPMqYMeLNHcDNQDfgC9/SU1TzLjDFeysjMJOeDvjQgnQNXf0CVZJto3pQf3lxZvA94uAyyGOOXcg8eYMfIa2hVtIZtFw8nqcOVTkcy5g/s8kVjfKgg/wjrh19L64IVrOnyGkkXDnQ6kjF/YoXAGB9Rl4tf3x1Eu7yFLGk9mJY973E6kjEnZIXAGB+ZP/ZJOmVNZ369u+h4/RNOxzHmpLyZj6Aa7nkDkopvr6p3+C6WMRXb0kkv0SV9DIvirqTz7a87HceYEnlz+uhXwE/Ad7jHAzLGlGD118Nov/YNFkR1p/394xAbSdSUc94UgihVfcrnSYzxAxvmfUXjxUNYFJZCqwc/ITTUBpEz5Z83f6pME5FePk9iTAW3ddM64mfdz5ageiTd9ynRkRFORzLGK94UgodxF4M8ETngueX4OpgxFUlG1gFyJgwklELC+n9Etfh4pyMZ4zVvLiiLLYsgxlRUuUcKmffe/fR2rWXTRe+S3LiN05GMOS1eTYUkIr2BCzwPf1DVab6LZEzFUVjk4tP3/8HteV+xtfFtJF84wOlIxpw2b6aqfAV391Ca5/awiLzs62DGlHeqysTxI7klcxg7ErqReNObTkcy5ox4s0fQC2inqi4AERkP/Ao848tgxpR3Uz//iJu2vEBGbHNq3/2JzTVsKixvT3COK3a/sg9yGFOhfDsrlStWPs6+iHrU/Os0CLdDaabi8uZPmJeBX0VkDu6pKi8AnvZpKmPKsUXzfqDzL/eSExpPwl9n2AxjpsLz5qyhSSLyA9DRs+gpVd3l01TGlFOr1q2n3szbORIcRaV7phNauZbTkYw5ayftGhKRZp6f7YFaQLrnVtuzzJiAsnn3fgo+HkicHCR4wCdEVU92OpIxpaKkPYLHgHuAN06wToGLfZLImHJo78EjLB99L71Zw+4e71KjYYrTkYwpNSctBKp6dPD0K1Q1r/g6EbFr503AyD1SyKcjX+KvhbPY1fo+ana1awWMf/HmrKFfvFxmjN8pKHLx1tgPuevgSPbWvICaff/ldCRjSt1J9whEpCZQB4gUkXNwnzEEUAmIKoNsxjhKVXnlk++5b9cQ8qLrknDbBAgKdjqWMaWupGMElwODgLpA8UsmDwDP+jCTMeXCWzOX03vN36gcUkDYoE8hMs7pSMb4REnHCMYD40XkOlWdUoaZjHHchF82Uf+X52gbvBG9YSJUb+Z0JGN8xpsLylqJSMvjF6rqSz7IY4zjZq7cyabpb/BCyE8UXfgMwc2vcjqSMT7lTSE4WOx+BHAVsNo3cYxx1sJN+5j0yQTGhkykqMmVBF/4N6cjGeNz3lxZ/IfrCERkKDDLZ4mMcci63QcYOv5TxgQPQ+MbE3Ld+2DzDZsAcCb/yqNwH0A+JRHpKSJrRWSDiJx0fCIRuU5EVETsKh3jiJ3Zh3l5zMeM4SWiYisTMuATG0jOBIxT7hGIyArcVxIDBAPVgFMeHxCRYGAEcBnuoSkWiUiqqqYdt10s7vkOFpxedGNKR/bhAl4eNZFh+UOIqBRPyB3fQJX6Tscypsx4c4yg+JGyQmC3qhZ68bxOwAZV3QggIpOBPrgntynu78CrwJNevKYxpepwfhGvjf6Qfx4cTEhsNcLu/AbiEp2OZUyZOmXXkKpuAeJxf4lfC7T28rXrANuKPU73LDvGM3hdPVX9pqQXEpF7RGSxiCzOyMjw8u2NKVlBkYu3x47nmcxnkZhqRN49w4qACUjeTFX5AjAedzFIAMaJyOCzfWMRCcJ9odrjp9pWVUepaoqqplSrVu1s39oYXC5l5PgPeWjn07hiahBz7yyo7NWhL2P8jjddQwOAtkcHnvPMYfwb8I9TPG87UK/Y47qeZUfFAq2AH0QEoCaQKiK9VXWxV+mNOQOqyoRJH3LXlr+RF12bqvfNhNiaTscyxjHenDW0A/f1A0eF88cv9JNZBDQWkWQRCQNuAlKPrlTVbFVNUNUkVU0C5gNWBIzPffXZOPqte4KDkXWo8pdZVgRMwCtp0Ll/4z5bKBtYJSLfeh5fBiw81QuraqGIPID7moNgYKyqrhKRl4DFqppa8isYU/rmfj6cK1e9wK7IhtS5fzoSa12NxpTUNXT0L/MlwBfFlv/g7Yur6nRg+nHLXjjJtt29fV1jzsTyKa9ywcp/sSayLY0eSiUoKs7pSMaUC6cadM4Yv7Dp8+dps/IdFkZ0pc3DnxMSGe10JGPKjZK6hj5V1X7HXVB2jKq28WkyY0rJ9i9eIHnlO3wffgmdHppIRGSk05GMKVdK6hp62PPThl40Fdbur16kzrJhzAi5hE4PTCQ2yoqAMccrqWtop2eYiHGqelEZZjKmVGR8PYQav77NN8EXc879E4iPtSJgzImUePqoqhYBLhGpXEZ5jCkV+6YNodqSt5gWdDFt/jqB2lXsmIAxJ+PtfAQrPKeP5h5dqKoP+SyVMWdKlf3fDKHq4rdJlYtpdd946sXHOJ3KmHLNm0Iw1XMr7k8Hj41xnCoHUp+iyq/vM1UuodW9/6FB9UpOpzKm3POmEMSp6rDiC0Tk4ZNtbIwjXC5yv3qU2GXjmERP2t71Pk1qWo+mMd7wZoiJ206wbFAp5zDmzLmKODzlL0QvG8dY7U2LO96jRZ04p1MZU2GUdB3BzUB/IFlEig8HUQnY5+tgxnilqIAjn91N5JovGKHX03nQa7RNrOJ0KmMqlJK6hn4BduIeerr4vMUHgOW+DGWMVwqPcGTyIMI3TGeo62a63fZPUpLjnU5lTIVT0nUEW4AtInIpcFhVXSLSBGgGrCirgMacUFEBRybdSvjvM/mnaxCXDHqBcxtYETDmTHhzjGAuECEidYDZwC3AOF+GMqZELheHP7vnWBG47PYXrQgYcxa8KQSiqodwT1P5rqreALT0bSxjTkKVQ188ROSaqbztuonL73iRTslVnU5lTIXmVSEQkS64Zyo7OrdwsO8iGXMSquR++ShRKyYwRq/h/DtfISXJioAxZ8ubQvAI8AzwhWdimQbAHJ+mMuZ4qhz46kmil/2HcXoVHe54iw717ewgY0rDKS8oU9UfgR+LPd4I2PASpuyosv/Lp6iybDQT6cU5dw63U0SNKUUlXUfwtqo+IiJfc+L5CHr7NJkxAKrs+XIw1Ze9z6dBPUm5+32a1rJhI4wpTSXtEUzw/BxaFkGMOZFtX/0/6i0bzlfBPej8lw+on2ADyBlT2kq6jmCJ5+ePIlLNcz+jrIIZ8/sX/6DhsreYGXIxnR4YR604G0raGF8o8WCxiAwRkb3AWmCdiGSIyAknnzemNK2a+jINl73OD2EX0vGhiVYEjPGhkxYCEXkM6AZ0VNWqqloF6Ax0E5FHyyqgCTzzP3mNlstfYX54N9o/8gnxlaKcjmSMXytpj+AW4GZV3XR0geeMoYHArb4OZgKPq8jFj2Oe5tzV/+S3qC60fWQKlWyOYWN8rqSDxaGquvf4haqaISKhPsxkAlD+4VxWjLyFC3O+Z1ncpbT6ywRCwq0IGFMWSioE+We4zpjTcjBjK3tGXcs5+RuY1+B+zr3lH0iQN9c6GmNKQ0mFoK2I5JxguQARPspjAkzm2p+RyQOo7jrMvE7v0O1K63U0pqyVdPqojSdkfGr7j+NImPMEe7QKG6+cSLdO3ZyOZExA8mbOYmNKl6uITZ/8jeS1Y1girYi65SNSGiY7ncqYgGWFwJQpLTzC7+8PoFHGt8yI6EX7e0dRo0qs07GMCWhWCEyZOXIoh00jrqNZ7kK+SLiPK+79FxGh1gNpjNN8emqGiPQUkbUiskFEnj7B+sdEJE1ElovI9yJS35d5jHMy9+5m01uX0/jgIr5t9DzX3P+KFQFjygmfFQIRCQZGAFcALYCbRaTFcZv9CqSoahvgc+A1X+Uxzklbt479Iy6jQf46fjv3bS4b+AQi4nQsY4yHL/cIOgEbVHWjquYDk4E+xTdQ1TmeaTAB5gN1fZjHOODrH34heuJV1NbdpPcaT4crBjkdyRhzHF8WgjrAtmKP0z3LTuZOYMaJVojIPSKyWEQWZ2TYAKgVQV5BEW989AWd5txMfPAhCgZ8SYPOVzkdyxhzAuXiYLGIDARSgAtPtF5VRwGjAFJSUv40SY4pX7btO8SwcR8xOHsIweFRRN4xg+Cax/cKGmPKC18Wgu1AvWKP63qW/YGIXAo8B1yoqkd8mMeUgTlrdvPz5Nd5mbEUxNYl6s5UqJLkdCxjTAl8WQgWAY1FJBl3AbgJ6F98AxE5B3gf6Kmqe3yYxfhYfqGLN2cso8GCFxkc8iOH619M1E1jIdLmFjamvPNZIVDVQhF5AJgFBANjVXWViLwELFbVVOB1IAb4zHMWyVabC7ni2ZhxkL9PnMWj+/5Om5BNFJ73BJEXPwc2cJwxFYJPjxGo6nRg+nHLXih2/1Jfvr/xLVXl8yXpzEidzBtBw4gNU7h+EiHNejkdzRhzGsrFwWJT8eTkFTB46nLqrnqPMaGfU1S1EaH9P4aExk5HM8acJisE5rT9vGEvL302jycOv8VloUtwtbyO0N7vQHiM09GMMWfACoHx2sEjhbw8fTVLF/7E2Ihh1A7eC5e/SlDne8GuFDamwrJCYLzyy4a9PPnZMs47OIPUyA8JiaqK9PsGEs91Opox5ixZITAlyj1SyMszVvPF/LUMixnPpaFzof6FcO1oiK3hdDxjTCmwQmBO6ru03byYuorKOauZW3kkVfO3w0WD4fzHIMhGDjXGX1ghMH+yPeswQ1JX8W3aLh6L+4kHIsYSFBYPN38NSec5Hc8YU8qsEJhjCopcfPC/TQz7bj0xHOL7ehNpmPEdNLoM+r4H0QlORzTG+IAVAgPAwk37eP7LlazdfYC7Guzn6YOvErJ3O1z2EnR50K4SNsaPWSEIcFszD/HKzNVMX7GLOpUjmNl5Bc1WDIXYmnDHTKjXyemIxhgfs0IQoLIPFzBizgbG/byZ4CBhSLcobtn/DsHL/gtNe0GfERBV1emYxpgyYIUgwBQUuZi0cCtvfbuOrMMFXN++LoMTV1D5u7+BBEGvodDxLrtAzJgAYoUgQBS5lGnLdzDsu/Vs3JtLlwbxPN8jkRYrXoUZ4yCxK1w3BiqXNImcMcYfWSHwcy6XMmvVLt76bh3rdh+kaY1YRt/SgUsLf0A+vw0O7oLzHnVfHxBs/xyMCUT2P99PqSr/XbOHN2avI21nDg2qRfPvm8/hypo5BH0zCLbOg9rt4caPoF5Hp+MaYxxkhcDPFHn2AN778XeWp2eTWDWKN25oS59W8YT88jZ89QaERUPvf0O7gXZaqDHGCoG/OFJYxNSl2xk1dyOb9uaSFB/FK9e25roOdQlNnw+j+8LeddD6Brj8ZYip5nRkY0w5YYWggsvJK+DjBVv54H+byDhwhNZ1KjOif3t6tqpJ8P6N8MVdsHIKxCXCwCnQyCaFM8b8kRWCCmr97gN8OG8LU5emk5tfxPmNE3j7xnZ0bRiPHNgJ3zwCSydASDic/4R7oLiwaKdjG2PKISsEFUhhkYvvVu9m/C9bmLcxk7CQIK5uU5vbuyXRqk5l2LMGUl+E5Z+CKnS8010EbLhoY0wJrBBUANuzDjNlSTqTFm5lZ3YedeIi+VvPptyYUo/46DDY/D+Y+G9YPwtCIqH9rdD1QaiS5HR0Y0wFYIWgnMorKGJ22m4+W7yN/23Yiyp0axTPkN4tuaRZdUK0AFZOhQUjYecyiEqA7s+6rwqOjnc6vjGmArFCUI6oKsvSs5myJJ2vfttOTl4hdeIiefDixtzQoS71qkbBgV0w9xVYPBZyMyChKVz1FrS9GUIjnW6CMaYCskLgMFUlbWcO05bvZNryHWzbd5jwkCB6tqrJDR3q0bVhPEFBAtuXwJT3YNUX4CqEJpdD53uhwUU2LpAx5qxYIXCAqrJu90Gmr9jJ18t3sDEjl+AgoVujBB66uDE9WtakcmQoFBXAqimw4D1IXwRhse6un053Q3xDp5thjPETVgjKSGGRi8Vb9vNt2m6+TdvN1n2HEIFzk+O587xkrmhVi6rRYe6NiwrcXT8/vg4HdkDVBtDzVWjXHyIqOdsQY4zfsULgQ/tz8/n59738d/Ue/rt2D1mHCggLCaJbw3juu7AhlzavTvVKEe6NVd3dP8s+cV8Admgv1DsXrn7bPVWkDQVhjPERKwSlKL/QxdKt+/lpfQY/rd/Liu3ZqEJcVCgXN63OZS1qcEGTakSHe37th/ZB2mzY9CP8/l/YtxGCw6FpT/c4QI0vs/5/Y4zPWSE4C4VFLtJ25rBw0z5++T2T+RszOZRfRHCQ0D4xjkcvbcL5jRNoXacyIcFBkJ8LW3+AjT/Cprnu0z5RCI2G+l2h2yPQog9ExjnbMGNMQLFCcBryCor4bVsWizbtY+HmfSzdsp/c/CIAkuKjuK59Xc5vnECXhvHERoS6u3t2LIWfvnV/8W9bCK4CCAp1zwXc/RlocCHU6QDBoQ63zhgTqKwQnITLpWzKzGV5ehbLtmWzPD2LFduzKShSAJrVjOXa9nXpmFyVTklVqVk5wv3Fn70NNqTC1gWwdrr7MQK120GXv0LyhZDYBcKiHG2fMcYc5dNCICI9gWFAMDBGVV85bn048CHQAcgEblTVzb7MdCIul7J13yHW7MphWbr7S395ejYH8goBiAwNplWdStzRLZmOSVVJqRtFXP5u2L8JMn6AH9ZAxlr37Ui2+0VDIqBBd7joWWjS0yaCN8aUWz4rBCISDIwALgPSgUUikqqqacU2uxPYr6qNROQm4FXgRl9lAsg6lM+aXQdYu+sAa3blsHrnAdbvziY4/yCVJJf4oEO0quriyvoumsbmUz88l6pkEZSbATv2wJrtkLMD0P970ehqUK0ZtLkBqjeHOilQo6V19xhjKgRf7hF0Ajao6kYAEZkM9AGKF4I+wBDP/c+B4SIiqqqUsgVThlFj5fu4ioqohYu64uJyUaIln+igXIIiXP+38QHP7aioeIiu7p7MJfkC92BucfWhSn33EA82to8xpgLzZSGoA2wr9jgd6HyybVS1UESygXhgb/GNROQe4B6AxMTEMwoTGVeNrJjGREeGEx0ZTqXoCCLCQpHQSIisAhFx7rN1/nC/KkQn2F/2xhi/ViEOFqvqKGAUQEpKyhntLbS5pD9c0r9UcxljjD/w5eWq24F6xR7X9Sw74TYiEgJUxn3Q2BhjTBnxZSFYBDQWkWQRCQNuAlKP2yYVuM1z/3rgv744PmCMMebkfNY15OnzfwCYhfv00bGqukpEXgIWq2oq8AEwQUQ2APtwFwtjjDFlyKfHCFR1OjD9uGUvFLufB9zgywzGGGNKZkNaGmNMgLNCYIwxAc4KgTHGBDgrBMYYE+Ckop2tKSIZwJYzfHoCx1217Oesvf4rkNoK1t7SUF9Vq51oRYUrBGdDRBaraorTOcqKtdd/BVJbwdrra9Y1ZIwxAc4KgTHGBLhAKwSjnA5Qxqy9/iuQ2grWXp8KqGMExhhj/izQ9giMMcYcxwqBMcYEuIApBCLSU0TWisgGEXna6TylTUQ2i8gKEflNRBZ7llUVkW9FZL3nZxWnc54pERkrIntEZGWxZSdsn7i94/msl4tIe+eSn5mTtHeIiGz3fMa/iUivYuue8bR3rYhc7kzqMyMi9URkjoikicgqEXnYs9wvP98S2uvc56uqfn/DPQz270ADIAxYBrRwOlcpt3EzkHDcsteApz33nwZedTrnWbTvAqA9sPJU7QN6ATMAAc4FFjidv5TaOwR44gTbtvD8mw4Hkj3/1oOdbsNptLUW0N5zPxZY52mTX36+JbTXsc83UPYIOgEbVHWjquYDk4E+DmcqC32A8Z7744FrnItydlR1Lu45K4o7Wfv6AB+q23wgTkRqlUnQUnKS9p5MH2Cyqh5R1U3ABtz/5isEVd2pqks99w8Aq3HPZ+6Xn28J7T0Zn3++gVII6gDbij1Op+RffEWkwGwRWSIi93iW1VDVnZ77u4AazkTzmZO1z58/7wc83SFji3X1+U17RSQJOAdYQAB8vse1Fxz6fAOlEASC81S1PXAFcL+IXFB8pbr3Mf32XGF/b5/HSKAh0A7YCbzhaJpSJiIxwBTgEVXNKb7OHz/fE7TXsc83UArBdqBescd1Pcv8hqpu9/zcA3yBe9dx99FdZs/PPc4l9ImTtc8vP29V3a2qRarqAkbzf90DFb69IhKK+0txoqpO9Sz228/3RO118vMNlEKwCGgsIskiEoZ7buRUhzOVGhGJFpHYo/eBHsBK3G28zbPZbcBXziT0mZO1LxW41XN2yblAdrEuhgrruH7wvrg/Y3C39yYRCReRZKAxsLCs850pERHc85evVtU3i63yy8/3ZO119PN1+gh6Wd1wn2mwDvcR9+eczlPKbWuA+6yCZcCqo+0D4oHvgfXAd0BVp7OeRRsn4d5dLsDdR3rnydqH+2ySEZ7PegWQ4nT+UmrvBE97lnu+HGoV2/45T3vXAlc4nf8023oe7m6f5cBvnlsvf/18S2ivY5+vDTFhjDEBLlC6howxxpyEFQJjjAlwVgiMMSbAWSEwxpgAZ4XAGGMCnBUC49dEpKaITBaR3z3Db0wXkSYiklR8ZM9Sfs/7ROTW03zODyISMJOzm/IlxOkAxviK58KdL4DxqnqTZ1lb3GPWbCvpuWdDVd/z1Wsb4wu2R2D82UVAQfEvZlVdpqo/Fd/Is3fwk4gs9dy6epbXEpG5nrHhV4rI+SISLCLjPI9XiMijx7+pZ1z5Jzz3fxCRV0VkoYisE5HzPcsjPXsqq0XkCyCy2PN7iMg8T5bPRCRGROp7xuVPEJEgT94evvm1mUBjewTGn7UClnix3R7gMlXNE5HGuK/qTQH6A7NU9Z8iEgxE4R4QrI6qtgIQkTgvXj9EVTt5Jhp5EbgU+AtwSFWbi0gbYKnn9RKAwcClqporIk8Bj6nqSyLyKu6ByRYCaao627tfgzEls0JgDIQCw0WkHVAENPEsXwSM9QwQ9qWq/iYiG4EGIvJv4BvAmy/jo4OoLQGSPPcvAN4BUNXlIrLcs/xc3BOR/Ozu2SIMmOfZboyI3ADch7sgGVMqrGvI+LNVQAcvtnsU2A20xb0nEAbHJoe5APdIj+NE5FZV3e/Z7gfcX8hjvHj9I56fRZz6jy8BvlXVdp5bC1W9E0BEonCPPAkQ48X7GuMVKwTGn/0XCC82UQ8i0uZoP30xlYGd6h7+9xbcU5siIvWB3ao6GvcXfntP102Qqk7B3YVzpvPlzsXd9YSItALaeJbPB7qJSCPPumgRObqH8iowEXgB9zDFxpQK6xoyfktVVUT6Am97+trzcM/t/Mhxm74LTPGc8jkTyPUs7w48KSIFwEHgVtwzQ/1HRI7+EfXMGcYb6Xmd1binKlziyZwhIoOASSIS7tl2sGeI4o5AN1UtEpHrROR2Vf3PGb6/McfY6KPGGBPgrGvIGGMCnBUCY4wJcFYIjDEmwFkhMMaYAGeFwBhjApwVAmOMCXBWCIwxJsD9f3GOxE2tjn30AAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"# plot Jenks vs. optimized\n",
"quantile_plot(sigmoid_R, 'sigmoid')\n",
"quantile_plot(jenks_R, 'jenks')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9fc381f5",
"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.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