Skip to content

Instantly share code, notes, and snippets.

@santiago-salas-v

santiago-salas-v/synthetic_ramp_method.ipynb

Last active May 4, 2025 08:48
Show Gist options
  • Save santiago-salas-v/36a044cf0e7fe16dcab2d9af4247cd89 to your computer and use it in GitHub Desktop.
Save santiago-salas-v/36a044cf0e7fe16dcab2d9af4247cd89 to your computer and use it in GitHub Desktop.
Download ZIP
calculate synthetic ramps given a series of time points and coefficients of unit steps at each time. \n convolution ramps
Display the source blob
Display the rendered blob
Raw
synthetic_ramp_method.ipynb
{
"cells": [
{
"cell_type": "markdown",
"id": "d7ca8601",
"metadata": {},
"source": [
"synthetic ramp method\n",
"\n",
"calculate synthetic ramps using convolution for a given series of time points:\n",
"* define matrix $A$ with the coefficients $a_{i,j}$ of the unit steps in $n$ different series $i$ for each time $t_j$.\n",
"* define matrix $\\sigma(t-t_j)*\\sigma(t)=(t-t_j)\\cdot \\sigma(t-t_j)$ which is lower triangular with permutations of $(t-t_j)\\cdot \\sigma(t-t_j)$\n",
" * example with $t_0<t_1$: $(t_0-t_1)\\cdot \\sigma(t_0-t_1)=0$, and $(t_1-t_0)\\cdot \\sigma(t_1-t_0)=(t_1-t_0)$ due to step definition\n",
"* the matrix $\\sigma(t-t_j)\\cdot A$ has then the slopes $\\frac{dy_i}{dt}$ at times $t_j$ and $(t-t_j)\\cdot \\sigma(t-t_j)\\cdot A$ has the synthetic ramps by convolution.\n",
"$$\n",
"\\underbrace{\\left[\\begin{array}{cccc}\n",
"1 & 0 & 0 & ...\\\\\n",
"1 & 1 & 0 & ...\\\\\n",
"1 & 1 & 1 & ...\\\\\n",
"... & ... & ... & ...\n",
"\\end{array}\\right]}_{\\sigma(t-t_j)}\\cdot \n",
"\\underbrace{\\left[\\begin{array}{ccc}\n",
"a_{1,0} & a_{2,0} & ...\\\\\n",
"a_{1,1} & a_{2,1} & ...\\\\\n",
"a_{1,2} & a_{2,2} & ...\\\\\n",
"... & ... & ...\\\\\n",
"\\end{array}\\right]}_{A} = \n",
"\\underbrace{\\left[\\begin{array}{ccc}\n",
"a_{1,0} & a_{2,0} & ...\\\\\n",
"a_{1,0}+a_{1,1} & a_{2,0}+a_{2,1} & ...\\\\\n",
"a_{1,0}+a_{1,1}+a_{1,2} & a_{2,0}+a_{2,1}a_{2,2} & ...\\\\\n",
"... & ... & ...\\\\\n",
"\\end{array}\\right]}_{\\mathrm{slopes}}\n",
"$$\n",
"\n",
"$$\n",
"\\underbrace{\\left[\\begin{array}{cccc}\n",
"t_0-t_0 & 0 & 0 & ...\\\\\n",
"t_1-t_0 & t_1-t_1 & 0 & ...\\\\\n",
"t_2-t_0 & t_2-t_1 & t_2-t_2 & ...\\\\\n",
"... & ... & ... & ...\n",
"\\end{array}\\right]}_{(t-t_j)\\cdot \\sigma(t-t_j)}\\cdot \n",
"\\underbrace{\\left[\\begin{array}{ccc}\n",
"a_{1,0} & a_{2,0} & ...\\\\\n",
"a_{1,1} & a_{2,1} & ...\\\\\n",
"a_{1,2} & a_{2,2} & ...\\\\\n",
"... & ... & ...\\\\\n",
"\\end{array}\\right]}_{A} = \n",
"\\underbrace{\\left[\\begin{array}{ccc}\n",
"a_{1,0}\\cdot(t_0-t_0) & a_{2,0}\\cdot(t_0-t_0) & ...\\\\\n",
"a_{1,0}\\cdot(t_1-t_0)+a_{1,1}\\cdot(t_1-t_1) & a_{2,0}\\cdot(t_1-t_0)+a_{2,1}\\cdot(t_1-t_1) & ...\\\\\n",
"a_{1,0}\\cdot(t_2-t_0)+a_{1,1}\\cdot(t_2-t_1)+a_{1,2}\\cdot(t_2-t_2) & a_{2,0}\\cdot(t_2-t_0)+a_{2,1}\\cdot(t_2-t_1)+a_{2,2}\\cdot(t_2-t_2) & ...\\\\\n",
"... & ... & ...\\\\\n",
"\\end{array}\\right]}_{\\mathrm{synth\\ curves}}\n",
"$$\n",
"\n",
"example below with two series $i=\\{1,2\\}$"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "ffd86ec0",
"metadata": {
"lines_to_next_cell": 2
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"lower triangular (t-tj)sigma(t-tj)) \n",
" [[ 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]\n",
" [ 3. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]\n",
" [ 4. 1. 0. 0. 0. 0. 0. 0. 0. 0. 0. ]\n",
" [ 6. 3. 2. 0. 0. 0. 0. 0. 0. 0. 0. ]\n",
" [ 8. 5. 4. 2. 0. 0. 0. 0. 0. 0. 0. ]\n",
" [ 8.5 5.5 4.5 2.5 0.5 0. 0. 0. 0. 0. 0. ]\n",
" [11. 8. 7. 5. 3. 2.5 0. 0. 0. 0. 0. ]\n",
" [14. 11. 10. 8. 6. 5.5 3. 0. 0. 0. 0. ]\n",
" [16. 13. 12. 10. 8. 7.5 5. 2. 0. 0. 0. ]\n",
" [17. 14. 13. 11. 9. 8.5 6. 3. 1. 0. 0. ]\n",
" [20. 17. 16. 14. 12. 11.5 9. 6. 4. 3. 0. ]]\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAosAAAHrCAYAAACn9tfQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACUWElEQVR4nOzdeVxUVf8H8M8wDAwgjCsCsgiWooGmuONWGogplbllmZbZZhlRv7SnJ8XMNH1My7LMVCp92hR7tEUl93LJBRNF0RSDFMQVUGSZmfP7Y2QSZWCAmblzZz7v12tewMy5Z75nzsydL+fee45CCCFARERERFQFF6kDICIiIiL7xWSRiIiIiExiskhEREREJjFZJCIiIiKTmCwSERERkUlMFomIiIjIJCaLRERERGSSq9QBSEGv1+Ps2bPw9vaGQqGQOhwiIiKbEEKgqKgIAQEBcHHheBGZxymTxbNnzyIoKEjqMIiIiCSRk5ODwMBAqcMgmXDKZNHb2xuA4cPi4+MjcTRERES2UVhYiKCgIOP3IJE5nDJZrDj07OPjw2SRiIicDk/BotpwymSRiByHTi/we9Yl5BeVwNdbja6hjaF0cewvQmdsMxFJh8kiEcnW+sO5mL4uA7kFJcb7/DVqTBvSDgMj/CWMzHqcsc1EJC1ZXgqVlJQEhUJR6ebn5yd1WERkQ+sP5+K5FQcqJU0AkFdQgudWHMD6w7kSRWY9zthmIpKeLJNFALjrrruQm5trvKWnp0sdEhHZiE4vMH1dBkQVj1XcN31dBnT6qkrIkzO2mYjsg2wPQ7u6upo9mlhaWorS0lLj34WFhdYKi4hs4PesS7eNrt1MAMgtKMH45N/h66O2XWBWlF9YYlabf8+6hB6tmtguMCJyeLJNFk+cOIGAgAC4u7ujW7dueOeddxAWFlZl2VmzZmH69Ok2jpCIrCW/yHTSdLOtxy9YORL7Y+5rQ0RkLoUQQnbHLH7++WcUFxejdevWOHfuHN5++20cO3YMR44cQZMmt/9HXdXIYlBQEAoKCjh1DpEM7Tp5EY8s2V1juVFdghDcxNMGEVlf9sVifL03p8ZyX03ozpFFMqmwsBAajYbff1QrshxZjIuLM/4eGRmJHj16oFWrVvj888+RmJh4W3l3d3e4u7vbMkQisqKmDdzgogBMnZ6nAOCnUWPmQ5EOM6WMTi+w7fh55BWUVHneYkWbu4Y2tnVoROTgZHuBy828vLwQGRmJEydOSB0KEVlZZl4RHlmyp9pEEQCmDWnnMIkiAChdFJg2pB2Af9p4MwHgjfvbOlSbicg+OESyWFpaiqNHj8Lfn3OMETmyQ39fwchPd+HC1VKE+3ljzrD28NdUvoDFT6PGx491csg5BwdG+OPjxzrB75Y2V6SHu05ehAzPLCIiOyfLw9CvvvoqhgwZguDgYOTn5+Ptt99GYWEhxo4dK3VoRGQle09fwhPL9+JqqRYdghri8ye6oKGnGx7uFOhUq5kMjPDHfe38KrW58Ho5nl25Hyv3ZOMO3wZ4IjpU6jCJyIHIMln8+++/8cgjj+DChQto1qwZunfvjt27dyMkJETq0IjICnacOI8JX+xDSbke3UIbY+m4Lmjgbth9KV0UTndBR1VtnjIwHLN+PoYZP2SgZVMv3NPGV6LoiMjRyPJq6Pri1WBE8pGacQ4TVx5AmU6Pvq2b4ZPHouDhppQ6LLsjhMBrqw7hu/1/o4G7K1Ke74nWzb2lDovsDL//qC4c4pxFInJMa/84i2dX7EeZTo/Yu5rj08eZKJqiUCgw86FIdA1tjKulWoz/fC8uXi2teUMiohowWSQiu/Tt3hy89HUadHqBhzq2wEejO8HdlYliddxcXfDJY1EIbuyJnEvX8cyX+1Gq1UkdFhHJHJNFIrI7y3/LwmurD0EI4JGuwZg3vANcldxdmaOxlxuWjesMb7Ur9v11Ga+npPMKaSKqF+59iciufLTlT0xflwEAeKpXKN55KAIuDnx1szXc4euNj0Z3gtJFgZQDZ/DxtpNSh0REMsZkkYjsghACczccw9wNmQCASf3vxBv3t4VCwUSxLvq0bmacxHvO+kysP5wncUREJFdMFolIcnq9wPR1Gfhoi2EE7PW4cCTe15qJYj093qMlHu9hmFLs5W8O4vCZAokjIiI5YrJIRJLS6QVeT0lH8s7TAIAZD9yFZ/q2kjYoBzJ1cDv0vrMprpfr8NTn+3CusETqkIhIZpgsEpFkynV6JHxzEN/sy4GLAvjP8A4Y06Ol1GE5FFelCz4c3Qmtmnkhr7AEE77Yh+tlvEKaiMzHZJGIJFFSrsNzKw5g3R9n4eqiwMJHOmFYVKDUYTkkjYcKy8Z1QSNPFQ79XYBXV/0BvZ5XSBOReZgsEpHNFZdpMeGLffjl6Dm4ubrg08ejcH97f6nDcmghTbzwyWNRUCkV+PFQLhZsOiF1SEQkE0wWicimCkvKMXbZ79hx4gI83ZRIHtcF94Y3lzosp9AtrAlmPhgJAPhg0wn87+AZiSMiIjlgskhENnP5WhkeXbIHe09fhrfaFV+O74aedzSVOiynMqJLEJ7uEwYA+L9Vh3Ag+7LEERGRvWOySEQ2kV9UglGf7kb6mQI09nLDVxO6IyqkkdRhOaXJA8MxoG1zlGn1ePqL/Thz5brUIRGRHWOySERWd+bKdYxcvBuZ54rg6+2Ob57ujogWGqnDclpKFwXeH3U3wv28ceFqKcYn78XVUq3UYRGRnWKySERWdfrCNYz4ZBeyLlxDi4Ye+O7ZHrizubfUYTk9L3dXLB3XBU0buONYXhESvk6DjldIE1EVmCwSkdWcOFeEEYt34cyV6whr6oXvnu2BkCZeUodFN7Ro6IFPH4+Cm6sLfjmaj3fXH5M6JCKyQ0wWicgqDp8pwIjFu5BfVIo2zb3xzTM9ENDQQ+qw6Badghth7rD2AIBPt5/CN3uzJY6IiOwNk0Uisrj9f13CI5/uxuXicrQP1ODrp7ujmbe71GGRCQ/c3QKT+t8JAHhjzWHsPnVR4oiIyJ4wWSQii9r55wWMWfo7ikq16NqyMVY+1Q2NvNykDotqkND/Ttzf3h9avcCzK/bj9IVrUodERHaCySIRWcymo+cwLnkvist06H1nU3z+ZFd4q1VSh0VmcHFRYN7wDugQqMGV4nKM/3wvCq6XSx0WEdkBJotEZBE/HsrFM1/uR5lWj/vaNcdnYzvDw00pdVhUC2qVEkse7wx/jRonz1/DC/89AK1OL3VYRCQx2SeLs2bNgkKhQEJCgtShEDmt7/bl4MWvDkCrF3jg7gAserQT3F2ZKMqRr48aSx7vDA+VEjtOXMD0dRlSh0REEpN1srh37158+umnaN++vdShEDmtL3adxv+tOgS9AEZ1CcJ7I+6GSinrXYvTi2ihwYJRd0OhAL7c/Re+2HVa6pCISEKy3aNfvXoVjz76KJYsWYJGjbhkGJEUPtl2ElP/dwQA8ER0S8waGgmli0LiqMgSYu/yw2ux4QCA6esysP34eYkjIiKpyDZZnDhxIu6//34MGDCgxrKlpaUoLCysdCOiuhNC4L2NmZj9s2ES5xfvvQNTB7eDQsFE0ZE82zcMD3cKhE4vMHHlAfyZXyR1SEQkAVkmi19//TUOHDiAWbNmmVV+1qxZ0Gg0xltQUJCVIyRyXEIIvP3jUXyw+U8AwGsD2+CVmDZMFB2QQqHAO0Mj0KVlIxSVavFk8j5culYmdVhEZGOySxZzcnLw0ksvYcWKFVCr1WZt8/rrr6OgoMB4y8nJsXKURI5Jpxf415p0LP01CwAwPf4uPN/vDomjImtyd1Xik8eiENTYA9mXivHsCsMV70TkPBRCCFmtHP/999/joYceglL5z5WWOp0OCoUCLi4uKC0trfRYVQoLC6HRaFBQUAAfHx9rh0zkELQ6PV797g98f/AsXBTA7IfbY0RnjtI7i+PnivDwop0oKtVieFQg5gxrz9FkGeL3H9WF7EYW+/fvj/T0dBw8eNB469y5Mx599FEcPHiwxkSRiGqvVKvD8ysP4PuDZ+HqosD7ozoyUXQyrZt7Y+HojnBRAN/t/xufbj8ldUhEZCOuUgdQW97e3oiIiKh0n5eXF5o0aXLb/URUf9fLdHhmxX5sP34ebq4uWDS6Ewa0ay51WCSBfm188ebgdpi+LgOz1x9DWLMGuI/vBSKHJ7uRRSKynaKScoxd/ju2Hz8PD5USy8d1YaLo5Mb1bIlHuwVDCOClr9Nw5GyB1CERkZXJ7pxFS+A5G0Q1u1JchrHL9+KPnCvwdnfF8ie6oHPLxlKHRXagXKfHuOW/47c/LyJAo8b3L0TD19u8Cw5JWvz+o7rgyCIR3eZ8USlGfbobf+RcQSNPFf47oTsTRTJSKV2waHQUwpp64WxBCSZ8sR8l5TqpwyIiK2GySESV5BZcx8jFu3AsrwjNvN3xzTM9EBmokTossjMaTxWWjusCjYcKf+Rcwf+tOgQnPFBF5BSYLBKRUfbFYgz/ZBdOXbiGFg098N0zPdC6ubfUYZGdCm3qhU8ei4KriwLr/jiL9zedkDokIrICJotEBAD4M78IwxfvxN+Xr6NlE098+2wPtGzqJXVYZOd6tGqCtx80zESx4JcTWPfHWYkjIiJLk93UOURUfzq9wO9Zl5BfVAJfbzW83JUYt3wvLl0rQ+vmDbBifDf4+vCCBTLPqK7B+DP/Kj77NQuvfvcHghp7IrKFptJ7rGtoYyhdOIk3kRwxWSRyMusP52L6ugzkFpQY71MAEAAiW2jwxZNd0cjLTbL4SJ5eH9QWWReuYdOxfDy+dA/UKiXyi0qNj/tr1Jg2pB0GRvhLGCUR1QUPQxM5kfWHc/HcigOVEkXAkCgCwBPRLZkoUp0oXRR4/5GOCNCoUViirZQoAkBeQQmeW3EA6w/nShQhEdUVk0UiJ6HTC0xflwFT16sqAMzdkAmdnle0Ut14qJQoN/H+qbh3+roMvseIZIbJIpGT+D3r0m0jijcTAHILSvB71iXbBUUO5fesSzh/y4jizfgeI5InJotETiK/yHSiWJdyRLfie4zIMTFZJHIS5i7HxmXbqK74HiNyTEwWiZzEob+vVPu4AoYrVruGclk/qpuuoY3hr1HD1AQ5fI8RyROTRSIHJ4TAgl+OY9bPx4z33fplXvH3tCHtOBce1ZnSRYFpQ9oBuP09VoHvMSL5YbJI5MCEEHjnp6NY8IthGbb/i22DTx7rBD9N5cOAfho1Pn6sE+fAo3obGOGPj6t4j3m6KfkeI5IpTspN5EBuXpmlWQN3rDt0Fl/9ngMAmDq4HZ7sFQoAuK+dH1fXIKsZGOFvfI/tOnkBH2z+E3q9QI+wplKHRkR1wGSRyEFUtTJLhXcfjsTILsHGv5UuCvRo1cSW4ZGTqXiPdQ9rjI0Z53Asrwjf7c/BU73DpA6NiGqJh6GJHICplVkqaDxUNo6IyEChUGBMjxAAwIrdf0HPCbmJZIcji0QyZ87KLNPXZeC+dn481EySePDuFpj98zGcvliM7SfOo18bX6lDojrQ6XQoLy+XOgyyEJVKBaVSaVZZJotEMleblVl46Jmk4OXuimFRgVj+22l8uesvJosyI4RAXl4erly5InUoZGENGzaEn58fFIrqBxKYLBLJHFfNIDkY0z0Ey387jc2Z+ci5VIygxp5Sh0RmqkgUfX194enpWWNiQfZPCIHi4mLk5+cDAPz9q5+lgMkikcy5u5p36jFXzSAphTVrgN53NsWOExewYvdfeH1QW6lDIjPodDpjotikCY9MOBIPDw8AQH5+Pnx9fas9JC3LC1w+/vhjtG/fHj4+PvDx8UGPHj3w888/Sx0Wkc3lFZRgzvpj1ZbhqhlkLx7v0RIA8M2+HJSU66QNhsxScY6ipydHgh1RRb/WdC6qLJPFwMBAzJ49G/v27cO+fftw77334oEHHsCRI0ekDo3IZnIuFWP44p04daEYjTwNVztzZRayZ/eG+6JFQw9cKS7H2j/OSh0O1QIPPTsmc/tVlsnikCFDMGjQILRu3RqtW7fGzJkz0aBBA+zevVvq0Ihs4uT5qxj+yS7kXLqOkCaeWPdiL67MQnZP6aLAY90N0+h8ses0hOA0OkRyIPtzFnU6Hb777jtcu3YNPXr0qLJMaWkpSktLjX8XFhbaKjwiizuaW4gxS/fgwtUy3OnbACue6obmPmoENvLkyixk90Z2CcL8X47j8JlCpOVcQafgRlKHREQ1kOXIIgCkp6ejQYMGcHd3x7PPPos1a9agXbt2VZadNWsWNBqN8RYUFGTjaIks42DOFYz6dDcuXC3DXQE++OaZHmju889oYsWqGQ/c3QI9WjVhokh2p7GXG4a0DwAAfLnrL4mjIVvR6QV2nbyI/x08g10nL0Inw8nZk5OT0bBhQ6vVn5SUhLvvvttq9deHQsj0OEBZWRmys7Nx5coVrF69Gp999hm2bdtWZcJY1chiUFAQCgoK4OPjY8uwiepsz6mLeDJ5L66V6dApuCGWP9GVK7OQLB36+wriP/wNbkoX7Hz9XjRt4C51SE6jsLAQGo3G7O+/kpISZGVlITQ0FGp13WZUqGopUn+NGtOGtJPVKTLXr19HUVERfH2tM0/o1atXUVpaWuerzlNSUrB48WLs378fFy9eRFpaWo3Jp7n9K9uRRTc3N9xxxx3o3LkzZs2ahQ4dOuD999+vsqy7u7vxyumKG5GcbDt+HmOX/45rZTr0bNUEX47vxkSRZKt9YEN0CGqIMp0e3+zNkTocsiJTS5HmFZTguRUHsP5wrkSR1U55eTk8PDysligCQIMGDeo1PdG1a9cQHR2N2bNnWzAqA9kmi7cSQlQaPSRyFOsP5+Gpz/eipFyPe8N9sWxcF3i5y/50Y3JyY29aL1qr00scDZlLCIHiMq1Zt6KSckxbe6TKpUgr7ktam4GiknKz6qvNgdBVq1YhMjISHh4eaNKkCQYMGIBr164ZH1++fDnatm0LtVqN8PBwLFq0yPjY6dOnoVAo8O2336Jfv35Qq9VYsWJFlYeh161bh6ioKKjVaoSFhWH69OnQarXGx5OSkhAcHAx3d3cEBARg0qRJJmOu72HoMWPGYOrUqRgwYECd6zBFlt84//rXvxAXF4egoCAUFRXh66+/xtatW7F+/XqpQyOyqO/TzuCV7/6ATi9wf6Q/5o+8G25mTsJNZM8GRfrj7R+PIregBL8czcfACD+pQyIzXC/Xod3UDRapSwDIKyxBZNJGs8pnvBULT7ea05bc3Fw88sgjmDNnDh566CEUFRVhx44dxmRzyZIlmDZtGj788EN07NgRaWlpmDBhAry8vDB27FhjPZMnT8a8efOwfPlyuLu7Y+PGynFu2LABjz32GD744AP07t0bJ0+exNNPPw0AmDZtGlatWoX58+fj66+/xl133YW8vDz88ccfZr46wMqVK/HMM89UW2bx4sV49NFHza6zrmSZLJ47dw5jxoxBbm4uNBoN2rdvj/Xr1+O+++6TOjQii/nvnmy88X06hAAe7hSIdx+OhKuSiSI5BrVKiVFdgrBo60l8ses0k0WymNzcXGi1WgwdOhQhIYYR7MjISOPjM2bMwLx58zB06FAAQGhoKDIyMrB48eJKyWJCQoKxTFVmzpyJKVOmGLcJCwvDjBkz8Nprr2HatGnIzs6Gn58fBgwYAJVKheDgYHTt2tXsdsTHx6Nbt27VlmnevLnZ9dWHLJPFpUuXSh0CkVV9tuMU3v7xKADg8R4hSBpyF1x4ZTM5mEe7h+CTbSex8+RF/JlfhDt8vaUOiWrgoVIi461Ys8r+nnUJ45bvrbFc8hNdzFphykNlejm6m3Xo0AH9+/dHZGQkYmNjERMTg2HDhqFRo0Y4f/48cnJyMH78eEyYMMG4jVarhUajqVRP586dq32e/fv3Y+/evZg5c6bxPp1Oh5KSEhQXF2P48OFYsGABwsLCMHDgQAwaNAhDhgyBq6t5qZe3tze8ve3jM8FhCiI7IoTAwk0njIniM33DMD2eiSI5phYNPTCgrWFk5AtOoyMLCoUCnm6uZt1639kM/hr1bStLGeuC4aro3nc2M6s+c1cbUSqVSE1Nxc8//4x27dph4cKFaNOmDbKysqDXG86PXbJkCQ4ePGi8HT58+LaFPby8vKp9Hr1ej+nTp1eqJz09HSdOnIBarUZQUBAyMzPx0UcfwcPDA88//zz69OlT49J6FVauXIkGDRpUe1u5cqVZddWXLEcWiRyREAKz1x/D4m2nAACv3NcaL9x7B5fZIof2eI+W2JhxDqv3/43/i20DbzWv8ncUShcFpg1ph+dWHIACqHShi7WXIlUoFIiOjkZ0dDSmTp2KkJAQrFmzBomJiWjRogVOnTpV73P9OnXqhMzMTNxxxx0my3h4eCA+Ph7x8fGYOHEiwsPDkZ6ejk6dOtVYPw9DE1Eler1A0rojxtGVf9/fFk/1DpM4KiLri76jCcKaeeHU+WtYk3YGj/doKXVIZEEDI/zx8WOdbptn0c+K8yzu2bMHmzZtQkxMDHx9fbFnzx6cP38ebdu2BWC46njSpEnw8fFBXFwcSktLsW/fPly+fBmJiYlmP8/UqVMxePBgBAUFYfjw4XBxccGhQ4eQnp6Ot99+G8nJydDpdOjWrRs8PT3x5ZdfwsPDw3geZU1qexj60qVLyM7OxtmzhnXXMzMzAQB+fn7w86vfOcFMFokkptXpMSUlHav2/w2FApj5YCRGdwuWOiwim1AoFHi8ewiS1mXgi11/YUz3EI6mO5iBEf42XYrUx8cH27dvx4IFC1BYWIiQkBDMmzcPcXFxAICnnnoKnp6emDt3Ll577TV4eXkhMjISCQkJtXqe2NhY/PDDD3jrrbcwZ84cqFQqhIeH46mnngIANGzYELNnz0ZiYiJ0Oh0iIyOxbt26es2lWJ21a9fiiSeeMP49atQoAIYrs5OSkupVt2xXcKmP2s5gT2QtZVo9Xv7mIH5Mz4XSRYF5wzvgwY4tpA6LyKYKS8rR/Z1NKC7T4b9PdUPPO5pKHZLDkmIFF7JfDr+CC5HclZTr8OyK/fgxPRduShd8NLoTE0VySj5qFYZ2Mrz3eaELkf2xWbJYsRRfZmam8WokImd1rVSLJ5P3YvOxfLi7umDJ2M6cZ46cWsW5ihsz8nD2ynVpgyGiSmx2zmJERAQA4OWXX8aff/6JBg0a4K677kJERAQiIiJw//332yoUIkkVXC/HE8t/x4HsK/ByU2LpuC7oHmadc1iI5KJ1c290D2uM3acu4b97svFqbBupQyKiG6w+spiTY1gkvn///gCAn376CcePH8fWrVvx3HPPoVGjRkhNTbV2GER24dK1MoxeshsHsq9A46HCygndmSgS3VAxuvjV79ko1eqkDYaIjKyeLIaHh+PNN9+stIA3YLhaqWfPnnj66aexYMECa4dBJLlzhSUYuXgXjpwtRNMGbvj66e64O6ih1GER2Y372jWHn48aF6+V4ef0PKnDIaIbrJ4spqamYuPGjbjzzjuxfPlyaz8dkV36+3IxRizehRP5V+Hno8Y3z/RAW39eiU90M5XSxTht1Oe7TksbDBEZWT1Z7NmzJ/bs2YPZs2dj6tSp6NixI7Zu3WrtpyWyG6fOX8WIT3bhr4vFCG7sie+e7YFWzRpIHRaRXRrVNQgqpQJp2Vdw+EyB1OEQEWx4NfTjjz+O48ePY8iQIbj//vvx0EMP4c8//7TV0xNJ4lheIUYs3o2zBSVo1cwL3z7TA0GNPaUOi8hu+XqrEXdjVY8vOLpIZBdsOs+iEAIxMTF4+umnsXbtWkREROCVV15BUVGRLcMgsolDf1/BqE9348LVUrTz98E3z/SAn4aT2hLVZGxPw3Jo/zt4FpevlUkcDZFBcnIyGjZsaLX6k5KScPfdd1ut/vqwerL4ySefYPz48Wjfvj00Gg0GDBiA3377DRMnTsSiRYtw8OBBtGvXDvv27bN2KEQ2s/f0JYxesgdXisvRMbghvprQHU0buEsdFpEsdApuhHb+PijV6vHd/hypwyFL0OuArB1A+irDT738rnYfOXIkjh8/brX6X331VWzatKlO25aXl2Py5MmIjIyEl5cXAgIC8PjjjxvXia4vqy/3FxQUhO7duxtvnTt3hrt75S/Nd955B//9739x+PBha4ZixOX+yJp2nDiPCV/sQ0m5Ht3DGuOzsV3QwJ3LsBPVxjd7szF5dTqCGntg66v3WG0dYWcjyXJ/GWuB9ZOBwpsSF58AYOC7QLv4utVpY+Xl5VCpVFKHYVJBQQGGDRuGCRMmoEOHDrh8+TISEhKg1WqrHYyzm+X+cnJy8N133+GVV15BdHT0bYkiAIwfPx5Hjx61dihEVrfxSB7GJxsSxX5tmiH5ia5MFInqIL5DC2g8VMi5dB3bjudLHQ7VVcZa4NvHKyeKAFCYa7g/Y61VnnbVqlWIjIyEh4cHmjRpggEDBlSawm/58uVo27Yt1Go1wsPDsWjRIuNjp0+fhkKhwLfffot+/fpBrVZjxYoVVR6GXrduHaKioqBWqxEWFobp06dDq9UaH09KSkJwcDDc3d0REBCASZMmmYy5PoehNRoNUlNTMWLECLRp0wbdu3fHwoULsX//fmRnZ9epzpvZxdrQvr6+2Lx5s9RhENXL/w6ewXMrD6BMp0dchB8+HdMZapVS6rCIZMnDTYnhUYEAgM93cr1ouyEEUHbNvFtJIfDzawCqOoB54771kw3lzKnPzAOhubm5eOSRR/Dkk0/i6NGj2Lp1K4YOHYqKA6lLlizBG2+8gZkzZ+Lo0aN455138Oabb+Lzzz+vVM/kyZMxadIkHD16FLGxsbc9z4YNG/DYY49h0qRJyMjIwOLFi5GcnIyZM2cCMCSs8+fPx+LFi3HixAl8//33iIyMNPulXrlyJRo0aFDtbeXKlSa3LygogEKhsMh5lnYx5KFQKNC3b1+pwyCqs2/2ZmNKSjqEAIZ2bIE5w9rDVWkX/4sRydZj3UOw9LcsbDt+HqcvXEPLpl5Sh0TlxcA7ARaqTBhGHGcHmVf8X2cBt5rfA7m5udBqtRg6dChCQgwXS92cpM2YMQPz5s3D0KFDAQChoaHGZG/s2LHGcgkJCcYyVZk5cyamTJli3CYsLAwzZszAa6+9hmnTpiE7Oxt+fn4YMGAAVCoVgoOD0bVrV/PaCiA+Ph7dunWrtkzz5s2rvL+kpARTpkzB6NGjLXK6nV0ki0RytuzXLLz1QwYA4NFuwZjxQARceH4VUb21bOqFvq2bYWvmeXy5+y+8Obid1CGRDHTo0AH9+/dHZGQkYmNjERMTg2HDhqFRo0Y4f/48cnJyMH78eEyYMMG4jVarhUajqVRP586dq32e/fv3Y+/evcaRRADQ6XQoKSlBcXExhg8fjgULFiAsLAwDBw7EoEGDMGTIELi6mpd6eXt7w9vbuxYtNygvL8eoUaOg1+srHV6vDyaLRPXw0ZY/MXdDJgDg6T5heD0uHAoFE0UiSxnboyW2Zp7Hd/ty8EpMa3i68WtLUipPwwifOf7aCawcVnO5R1cBIT3Ne24zKJVKpKamYufOndi4cSMWLlyIN954A3v27IGnp6GOJUuW3DZqp1RWPm3Iy6v6UUy9Xo/p06dXOfqoVqsRFBSEzMxMpKam4pdffsHzzz+PuXPnYtu2bWZdLLNy5Uo888wz1ZZZvHgxHn30UePf5eXlGDFiBLKysrB582aLXcTLT1096fQCv2ddQn5RCXy91ega2tgprtqzdLvl8DpWjtEd246fxyfbTgEAEgbciZf638lEkcjC+rZuhuDGnsi+VIz5qccR0UJjt/sIS7PL/aJCYdahYABAq3sNVz0X5qLq8xYVhsdb3Qu4WPb8boVCgejoaERHR2Pq1KkICQnBmjVrkJiYiBYtWuDUqVOVkqy66NSpEzIzM3HHHXeYLOPh4YH4+HjEx8dj4sSJCA8PR3p6Ojp16lRj/bU9DF2RKJ44cQJbtmxBkyZNzG9MDWSXLM6aNQspKSk4duwYPDw80LNnT7z77rto06aNzWNZfzgX09dlILegxHifv0aNaUPaYeCNFQgckaXbLYfXsaoYK7wxqC0m9AmTICoix+fiokCXlo2QfakYS3ZkGe+3t32Epclhv1gjF6VhepxvHwegQOWE8UbSO3C2xRPFPXv2YNOmTYiJiYGvry/27NmD8+fPo23btgAMVx1PmjQJPj4+iIuLQ2lpKfbt24fLly8jMTHR7OeZOnUqBg8ejKCgIAwfPhwuLi44dOgQ0tPT8fbbbyM5ORk6nQ7dunWDp6cnvvzyS3h4eBjPo6xJbQ5Da7VaDBs2DAcOHMAPP/wAnU6HvLw8AEDjxo3h5uZmdruqIrsz8Ldt24aJEydi9+7dSE1NhVarRUxMTKVL4m1h/eFcPLfiwG3JQ15BCZ5bcQDrD+faNB5bsXS75fA6moqxQlBjDxtHROQ81h/ORcqBM7fdb0/7CEuTw37RbO3igRFfAD63JLg+AYb7rTDPoo+PD7Zv345BgwahdevW+Pe//4158+YhLi4OAPDUU0/hs88+Q3JyMiIjI9G3b18kJycjNDS0Vs8TGxuLH374AampqejSpQu6d++O9957z5gMNmzYEEuWLEF0dDTat2+PTZs2Yd26dRYd8avw999/Y+3atfj7779x9913w9/f33jbuXNnveu3+qTc1nb+/Hn4+vpi27Zt6NOnj1nb1HdSbp1eoNe7m00mDwDQwN0V46JbwsWBDkvqhUDyb6dxtVRrskxt2m3p+qyhphgVAPw0avw6+V7pDw8RORhn3Ndae58jyaTcgGHFlr92AlfPAQ2aG85RtPCIItWeuf0ru8PQtyooKABgGGY1pbS0FKWlpca/CwsL6/Wcv2ddqnbnBQBXS7X4cPOf9XoeObJ0u+39dRQAcgtK8HvWJfRoZfn/FomcGfe1t5PtPsdFCYT2ljoKqiNZJ4tCCCQmJqJXr16IiIgwWW7WrFmYPn26xZ43v6j6nVeF3nc2RagDzQuWdeEadpy4UGM5c9tt6fqswdwYzX1PEJH5nHFfy30O2SNZJ4svvPACDh06hF9//bXacq+//nqlk1YLCwsRFGTmJKBV8PU2byj++X53yOs/vxrsOnnRrJ2Yue22dH3WYG6M5r4niMh8zriv5T6H7JHsLnCp8OKLL2Lt2rXYsmULAgMDqy3r7u4OHx+fSrf66BraGP4aNUydLaKA4aq1rqGmD43LkaXbLYfXUQ4xEjkqZ/z8OWObyf7JLlkUQuCFF15ASkoKNm/eXOurlyxB6aLAtCGGlQRu/UBX/D1tSDuHu+DB0u2Ww+sohxiJHJUzfv6csc1k/2SXLE6cOBErVqzAf//7X3h7eyMvLw95eXm4fv26TeMYGOGPjx/rBD9N5UMBfho1Pn6sk3zmwaolS7dbDq+jHGIkclTO+Pmzxzbr9XqbPydZn7n9Krupc0ytkLF8+XKMGzfOrDrqO3XOzexyhn0b4Aou9hkjkaNyxs+fNdpc2+8/vV6PEydOQKlUolmzZnBzc+NKVQ5ACIGysjKcP38eOp0Od955J1xcTI8fyi5ZtARLJotERERyUZfvv7KyMuTm5qK4uNjK0ZGteXp6wt/fv8YVXmR9NTQRERFZl5ubG4KDg6HVaqHT6aQOhyxEqVTC1dXVrJFiJotERERULYVCAZVKBZVKJXUoJAHZXeBCRERERLbDZJGIiIiITOJhaCIisn96HfDXTuDqOaBBcyCkp2G9YUfmjG0mu8RkkYiI7FvGWmD9ZKDw7D/3+QQAA98F2sVLF5c1OWObyW7xMDQREdmvjLXAt49XTpoAoDDXcH/GWmnisiZnbDPZNY4sEhGRfdLrDKNrqGo64Bv3rXsJ0OuBaiYUlhW9HvjxZZhuswJYPwUIv5+HpMlmmCwSEZF9+mvn7aNrt7p+CVg11jbx2AUBFJ4xvDahvaUOhpwEk0UiIrJPV8+ZV67JnYBXU+vGYivXLgAXT9RcztzXhsgCmCwSEZF90pu5Wsjg+Y4zypa1A/h8cM3lGjS3fixENzjISR5ERORQTqQCP71aQyEF4NPCMKWMowjpabjqGaaWYHPANpPdY7JIRET2Q68Hts0FVg4HSguBJnfAkDjdmjzd+HvgbMe60MNFaZgeB4DTtJnsHpNFIiKyDyWFwLdjgC1vAxBA5yeB53YBI74AfPwrl/UJMNzviHMOtot3vjaTXVMIIaq6Pt+hFRYWQqPRoKCgAD4+PlKHQ0RE548D3zwKXDgOKN2A+98DOo3553FnXM3ECm3m9x/VBS9wISIiaR39AVjzLFBWZDgfb8SXQGBU5TIuSse5iMVczthmsktMFomISBp6HbDlHWDHfwx/h/QChicDDZpJGhYRVcZkkYiIbO/6ZWD1BODPVMPf3Z8H7nsLUKqkjYuIbsNkkYiIbOvcEeDrR4HLWYCrBxD/AdB+hNRREZEJTBaJiMh2Dq8G/vcCUF4MNAwGRq4E/NtLHRURVYPJYn054xV6gPO221E4Uv85UlvMJYc23xpjYFdgywxg50LD42H3AMOWAZ6NpY2TiGoky2Rx+/btmDt3Lvbv34/c3FysWbMGDz74oO0DyVgLrJ9ceaF7nwDDhKqOPA+Ws7bbUThS/zlSW8wlhzZXFaPSDdCVGX7v9TJw75v2l+ASUZVkOSn3tWvX0KFDB3z44YfSBZGxFvj28co7QwAozDXcn7FWmriszVnb7Sgcqf8cqS3mkkObTcVYkSj2eBEYkMREkUhGZDmyGBcXh7i4OOkC0OsM/zWjqvnMBQCFYU1T37aOtUPU626s1VpNu9dPAcLvd6x2OwpHet+a816US1vMJYc2VxvjDUdSgPumO06/EDkBWSaLtVVaWorS0lLj34WFhfWr8K+dt//XXIkwnKfzYef6PY/sCKDwjOH14USy9sep3reO1BZzyaTN3EcQyY5TJIuzZs3C9OnTLVfh1XPmlVO6O9acYbpyQFdaczlzXx+yLUd635r7XpRDW8wlhzZzH0HkkJwiWXz99deRmJho/LuwsBBBQUF1r7BBc/PKPbbasf57ztoBfD645nLmvj5kW470vjX3vSiHtphLDm3mPoLIIcnyApfacnd3h4+PT6VbvYT0NFx9CIWJAgrD+qYhPev3PPbGWdvtKIz9Z4qM+s8Z34tyaLMcYiSiWnOKZNHiXJSGaSoA3L5TvPH3wNmOdwJ3te2+wRHb7Sgq9d+tZPa+dcbPoBzaLIcYiajWZJksXr16FQcPHsTBgwcBAFlZWTh48CCys7NtF0S7eGDEF4CPf+X7fQIM99vLfGeWZqrdANC8HdB2iO1jIvO1izeM7NxKju9bZ/wMyqHNcoiRiGpFIYSoZo4D+7R161bcc889t90/duxYJCcn17h9YWEhNBoNCgoK6n9IWg4rKVjDze3W64C1LxpObH/oU6DDSKmjI1Mu/Al8GAUolMColUDZNfm/b53xMyiHNsshRidk0e8/chqyvMClX79+sJsc10XpOCfQ18at7S7IATbPADa8Dtx5H5fwslfHfjD8DOsLtJFwrlJLcsbPoBzaLIcYicgssjwMTXao5ySgWVug+CKw8U2poyFTKpLF8PuljYOIiGSDySJZhqsbMOR9w+8HVwBZ26WNh25XmAv8vdfwexsmi0REZB4mi2Q5wd2AzuMNv69LAMpLJA2HbpH5k+Fni85VX6BERERUBSaLZFkDpgEN/IBLJ4Ed86SOhm5WcQi6rRmTJhMREd3AZJEsS60B4m7Ms/brfCD/mLTxkMH1K/+cGhDO6Y2IiMh8TBbJ8to9ALQeCOjLgR8SAL1e6ojoRCqg1wJN2wBN75A6GiIikhEmi2R5CgUw6D+AygvI3gWkfSF1RHRsneEnD0ETEVEtMVkk62gYBNz7b8PvG6cCReekjceZlZcAJ34x/M4pc4iIqJaYLJL1dHsG8L8bKC0A1k+ROhrndWorUH4N8A4AAjpJHQ0REckMk0WyHhelYe5FhQtwJMVw3hzZ3s0TcSsU0sZCRESyw2SRrCvgbqD784bff0g0rEVMtqPXAZk/G37n+YpERFQHTBbJ+vq9DmiCgYJsYOssqaNxLjl7gOILgLohEBItdTRERCRDTBbJ+twbAPffmKB71yIg9w9p43EmR28cgm49EFCqpI2FiIhkicki2UbrGOCuhwChA9a9ZDg8StYlBKfMISKiemOySLYz8F3AXQOcTQN+/1TqaBzfucPAlWzAVQ20ulfqaIiISKaYLJLteDcH7ksy/L75baDgb0nDcXgVh6Bb9QfcvKSNhYiIZIvJItlWp3FAUHeg7Crw0/8ZDpWSdRz70fCTh6CJiKgemCySbbm4GOZedFEBmT8BR9dJHZFjunwaOJcOKJSGi1uIiIjqiMki2Z5vONArwfD7T/8HlBRIGo5DqjgEHdIT8GwsbSxERCRrTBZJGr1fBRq3Aq7mAZvekjoax1NxCDqch6CJiKh+XKUOgJyUSg0MWQB8PgTYuxSIGAbotcDVc0CD5oYRMRel1FFWptcBf+20/xiP/QRk7zT8zUPQRERUT7IdWVy0aBFCQ0OhVqsRFRWFHTt2SB0S1VZoH+DuRwEI4PPBhtvq8YafCyKAjLVSR/iPjLWGmOQQ47eP/XNfcpx9xUhERLIjy2Txm2++QUJCAt544w2kpaWhd+/eiIuLQ3Z2ttShUW217GX4qddWvr8wF/j2cftIdDLWGmIpPFv5fsZIREROQCGE/OYu6datGzp16oSPP/7YeF/btm3x4IMPYtasmtceLiwshEajQUFBAXx8fKwZKlVHrzOMhN2a4BgpAJ8AICFdusO9jJGIHAi//6guZHfOYllZGfbv348pU6ZUuj8mJgY7d+6scpvS0lKUlpYa/y4sLLRqjGSmv3ZWk+AAgAAKzwAfRwNqjc3CqqSkwHFi/GsnENrbZmEREZFjkF2yeOHCBeh0OjRv3rzS/c2bN0deXl6V28yaNQvTp0+3RXhUG1fPmVfu/FHrxmEJcojR3NebiIjoJrJLFisoFIpKfwshbruvwuuvv47ExETj34WFhQgKCrJqfGSGBs1rLgMA9/wL8G1n3VhMyc8AtrxTczk5xGju601ERHQT2SWLTZs2hVKpvG0UMT8//7bRxgru7u5wd3e3RXhUGyE9DefSFeYCqOrU2Rvn2vV+Vbpz7doMAvYnO0aMIT1tHBgRETkC2V0N7ebmhqioKKSmpla6PzU1FT178stQVlyUwMB3b/xx66jwjb8Hzpb2ogzGSERETk52ySIAJCYm4rPPPsOyZctw9OhRvPzyy8jOzsazzz4rdWhUW+3igRFfAD7+le/3CTDc3y5emrhuxhiJiMiJyXLqHMAwKfecOXOQm5uLiIgIzJ8/H3369DFrW04dYIfksjoKYyQiGeP3H9WFbJPF+uCHhYiInBG//6guZHeBiyVU5Mecb5GIiJxJxfeeE44TUT04ZbJYVFQEAJw+h4iInFJRURE0GokWEiDZccrD0Hq9HmfPnoW3t7fJuRlro2LexpycHNkP67Mt9oltsU9si31iW0wTQqCoqAgBAQFwcZHlNa4kAaccWXRxcUFgYKDF6/Xx8ZH9jqkC22Kf2Bb7xLbYJ7alahxRpNrivxVEREREZBKTRSIiIiIyicmiBbi7u2PatGkOsaQg22Kf2Bb7xLbYJ7aFyLKc8gIXIiIiIjIPRxaJiIiIyCQmi0RERERkEpNFIiIiIjKJySIRERERmcRkkYiIiIhMYrJopkWLFiE0NBRqtRpRUVHYsWNHteW3bduGqKgoqNVqhIWF4ZNPPrFRpKbNmjULXbp0gbe3N3x9ffHggw8iMzOz2m22bt0KhUJx2+3YsWM2irpqSUlJt8Xk5+dX7Tb22CcA0LJlyypf44kTJ1ZZ3p76ZPv27RgyZAgCAgKgUCjw/fffV3pcCIGkpCQEBATAw8MD/fr1w5EjR2qsd/Xq1WjXrh3c3d3Rrl07rFmzxkot+Ed1bSkvL8fkyZMRGRkJLy8vBAQE4PHHH8fZs2errTM5ObnKviopKZGsLQAwbty422Lq3r17jfXaW78AqPL1VSgUmDt3rsk6peoXc/bBcvrMkPNgsmiGb775BgkJCXjjjTeQlpaG3r17Iy4uDtnZ2VWWz8rKwqBBg9C7d2+kpaXhX//6FyZNmoTVq1fbOPLKtm3bhokTJ2L37t1ITU2FVqtFTEwMrl27VuO2mZmZyM3NNd7uvPNOG0RcvbvuuqtSTOnp6SbL2mufAMDevXsrtSM1NRUAMHz48Gq3s4c+uXbtGjp06IAPP/ywysfnzJmD9957Dx9++CH27t0LPz8/3HfffSgqKjJZ565duzBy5EiMGTMGf/zxB8aMGYMRI0Zgz5491moGgOrbUlxcjAMHDuDNN9/EgQMHkJKSguPHjyM+Pr7Gen18fCr1U25uLtRqtTWaYFRTvwDAwIEDK8X0008/VVunPfYLgNte22XLlkGhUODhhx+utl4p+sWcfbCcPjPkRATVqGvXruLZZ5+tdF94eLiYMmVKleVfe+01ER4eXum+Z555RnTv3t1qMdZFfn6+ACC2bdtmssyWLVsEAHH58mXbBWaGadOmiQ4dOphdXi59IoQQL730kmjVqpXQ6/VVPm6vfQJArFmzxvi3Xq8Xfn5+Yvbs2cb7SkpKhEajEZ988onJekaMGCEGDhxY6b7Y2FgxatQoi8dsyq1tqcrvv/8uAIi//vrLZJnly5cLjUZj2eBqqaq2jB07VjzwwAO1qkcu/fLAAw+Ie++9t9oy9tAvQty+D5bzZ4YcG0cWa1BWVob9+/cjJiam0v0xMTHYuXNnldvs2rXrtvKxsbHYt28fysvLrRZrbRUUFAAAGjduXGPZjh07wt/fH/3798eWLVusHZpZTpw4gYCAAISGhmLUqFE4deqUybJy6ZOysjKsWLECTz75JBQKRbVl7bFPbpaVlYW8vLxKr7u7uzv69u1r8rMDmO6r6raRQkFBARQKBRo2bFhtuatXryIkJASBgYEYPHgw0tLSbBNgDbZu3QpfX1+0bt0aEyZMQH5+frXl5dAv586dw48//ojx48fXWNYe+uXWfbCjf2ZIvpgs1uDChQvQ6XRo3rx5pfubN2+OvLy8KrfJy8ursrxWq8WFCxesFmttCCGQmJiIXr16ISIiwmQ5f39/fPrpp1i9ejVSUlLQpk0b9O/fH9u3b7dhtLfr1q0bvvjiC2zYsAFLlixBXl4eevbsiYsXL1ZZXg59AgDff/89rly5gnHjxpksY699cquKz0dtPjsV29V2G1srKSnBlClTMHr0aPj4+JgsFx4ejuTkZKxduxZfffUV1Go1oqOjceLECRtGe7u4uDisXLkSmzdvxrx587B3717ce++9KC0tNbmNHPrl888/h7e3N4YOHVptOXvol6r2wY78mSF5c5U6ALm4dZRHCFHtyE9V5au6XyovvPACDh06hF9//bXacm3atEGbNm2Mf/fo0QM5OTn4z3/+gz59+lg7TJPi4uKMv0dGRqJHjx5o1aoVPv/8cyQmJla5jb33CQAsXboUcXFxCAgIMFnGXvvElNp+duq6ja2Ul5dj1KhR0Ov1WLRoUbVlu3fvXunCkejoaHTq1AkLFy7EBx98YO1QTRo5cqTx94iICHTu3BkhISH48ccfq0207LlfAGDZsmV49NFHazz30B76pbp9sKN9Zkj+OLJYg6ZNm0KpVN72H1p+fv5t/8lV8PPzq7K8q6srmjRpYrVYzfXiiy9i7dq12LJlCwIDA2u9fffu3SUfGbmVl5cXIiMjTcZl730CAH/99Rd++eUXPPXUU7Xe1h77pOLq9Np8diq2q+02tlJeXo4RI0YgKysLqamp1Y4qVsXFxQVdunSxu77y9/dHSEhItXHZc78AwI4dO5CZmVmnz4+t+8XUPtgRPzPkGJgs1sDNzQ1RUVHGK1QrpKamomfPnlVu06NHj9vKb9y4EZ07d4ZKpbJarDURQuCFF15ASkoKNm/ejNDQ0DrVk5aWBn9/fwtHVz+lpaU4evSoybjstU9utnz5cvj6+uL++++v9bb22CehoaHw8/Or9LqXlZVh27ZtJj87gOm+qm4bW6hIFE+cOIFffvmlTv9kCCFw8OBBu+urixcvIicnp9q47LVfKixduhRRUVHo0KFDrbe1Vb/UtA92tM8MORAprqqRm6+//lqoVCqxdOlSkZGRIRISEoSXl5c4ffq0EEKIKVOmiDFjxhjLnzp1Snh6eoqXX35ZZGRkiKVLlwqVSiVWrVolVROEEEI899xzQqPRiK1bt4rc3Fzjrbi42Fjm1rbMnz9frFmzRhw/flwcPnxYTJkyRQAQq1evlqIJRq+88orYunWrOHXqlNi9e7cYPHiw8Pb2ll2fVNDpdCI4OFhMnjz5tsfsuU+KiopEWlqaSEtLEwDEe++9J9LS0oxXCM+ePVtoNBqRkpIi0tPTxSOPPCL8/f1FYWGhsY4xY8ZUmlngt99+E0qlUsyePVscPXpUzJ49W7i6uordu3dL1pby8nIRHx8vAgMDxcGDByt9fkpLS022JSkpSaxfv16cPHlSpKWliSeeeEK4urqKPXv2SNaWoqIi8corr4idO3eKrKwssWXLFtGjRw/RokUL2fVLhYKCAuHp6Sk+/vjjKuuwl34xZx8sp88MOQ8mi2b66KOPREhIiHBzcxOdOnWqNN3M2LFjRd++fSuV37p1q+jYsaNwc3MTLVu2NLkTsyUAVd6WL19uLHNrW959913RqlUroVarRaNGjUSvXr3Ejz/+aPvgbzFy5Ejh7+8vVCqVCAgIEEOHDhVHjhwxPi6XPqmwYcMGAUBkZmbe9pg990nFND633saOHSuEMEwFMm3aNOHn5yfc3d1Fnz59RHp6eqU6+vbtayxf4bvvvhNt2rQRKpVKhIeH2yQRrq4tWVlZJj8/W7ZsMdmWhIQEERwcLNzc3ESzZs1ETEyM2Llzp6RtKS4uFjExMaJZs2ZCpVKJ4OBgMXbsWJGdnV2pDjn0S4XFixcLDw8PceXKlSrrsJd+MWcfLKfPDDkPhRA3zvInIiIiIroFz1kkIiIiIpOYLBIRERGRSUwWiYiIiMgkJotEREREZBKTRSIiIiIyickiEREREZnEZJGIiIiITGKySEREREQmMVkkIiIiIpOYLBIRERGRSa5SByAFvV6Ps2fPwtvbGwqFQupwiIiIbEIIgaKiIgQEBMDFheNFZB7ZJYuzZs1CSkoKjh07Bg8PD/Ts2RPvvvsu2rRpY3YdZ8+eRVBQkBWjJCIisl85OTkIDAyUOgySCdkli9u2bcPEiRPRpUsXaLVavPHGG4iJiUFGRga8vLzMqsPb2xuA4cPi4+NjzXCJbOPvtcDOMaYf7/klEBhvu3iIyC4VFhYiKCjI+D1IZA6FEEJIHUR9nD9/Hr6+vti2bRv69Olj1jaFhYXQaDQoKChgskjyp9cBa1sCxX+bKKAAPAOB+CzARWnLyIjIzvD7j+pCdiOLtyooKAAANG7c2GSZ0tJSlJaWGv8uLCy0elxENnN+RzWJIgAIoDjHUK55P1tFRUREDkLWZ7cKIZCYmIhevXohIiLCZLlZs2ZBo9EYbzxfkRzK9VzLliMiIrqJrJPFF154AYcOHcJXX31VbbnXX38dBQUFxltOTo6NIiSyAQ9/y5YjIiK6iWwPQ7/44otYu3Yttm/fXuMVXe7u7nB3d7dRZEQ21qy34ZzE4jMAqjoF+cY5i8162zoyIiJyALIbWRRC4IUXXkBKSgo2b96M0NBQqUMikpaLEoh6HwAgcOu8oTf+jlrAi1uIiKhOZJcsTpw4EStWrMB///tfeHt7Iy8vD3l5ebh+/brUoRFJJ2go0HsVhEeLyvd7BgK9VxkeJyIiqgPZTZ1jasWV5cuXY9y4cWbVwakDyFEVl5biyTnvwVd1GXMei4U6oB9HFMkh6HTAjh1Abi7g7w/07g0o+dauNX7/UV3I7pxFmeW2RLalUGL3tfYAgNm+fZkokkNISQFeegn4+6YZogIDgfffB4Zy0JzI6mR3GJqIiJxHSgowbFjlRBEAzpwx3J+SIk1cRM6EySIREdklnc4woljVAaWK+xISDOWIyHqYLBIRkV3aseP2EcWbCQHk5BjKEZH1MFkkIiK7lGvmokPmliOiumGySEREdsnfzEWHzC1HRHXDZJGIiOxS796Gq55NzJgGhQIICjKUIyLrYbJIRER2Sak0TI8DAApF5atcKhLIBQs43yKRtTFZJCIiuzV0KLBqFRAQUPn+wEDD/Zxnkcj6ZDcpNzkovQ44vwO4ngt4+APNenNCaSICYEgI74vT4c7xe6G7qsaXL7bHffcqOaJIZCNMFkl6OSnA/peA4pvmyPAMBKLe55rGRATAcKhZHXwJANCnbyQTRQnodDqUl5dLHQZZiEqlgtLMDxKTRZJWTgqwYxiAW2bdLT5juL/3KiaMREQSEkIgLy8PV65ckToUsrCGDRvCz88PClNXkd3AZJGko9cZRhRvTRSBG/cpgP0JQIsHeEiaiEgiFYmir68vPD09a0wsyP4JIVBcXIz8/HwAgH8N808xWSTpnN9R+dDzbQRQnGMo17yfraIiIqIbdDqdMVFs0qSJ1OGQBXl4eAAA8vPz4evrW+0haV4NTdK5buayC+aWIyIii6o4R9HT01PiSMgaKvq1pnNRmSySdDzMXHbB3HJERGQVPPTsmMztVx6GJuk062246rn4DKo+b1FheLwZl2cgIpIznQ7YscOwjre/v2HVHV7RLh8cWSTpuCgN0+MAELj1v5sbf0ctsJuLW3Q6YOtW4KuvDD91OqkjIiKyfykpQMuWwD33AKNHG362bGm4n+SBySJJK2go0HsVhEeLyvd7BtrVtDnc2RER1V5KCjBsGPD3LdcynjljuF9O+9Dk5GQ0bNjQavUnJSXh7rvvtlr99cFkkaQXNBQlg/7EqJPvYFL2/6Gkzy9AfJZdJYqOsrMjshm9Dji3FTj9leGnnkPxzkanA156CRBVnGVUcV9CgnyO0owcORLHjx+3Wv2vvvoqNm3aVOftU1JSEBsbi6ZNm0KhUODgwYMWi43JItkHhRK7r7XH2it9offta1eHnh1pZ0dkEzkpwNqWwKZ7gJ2jDT/XtjTcT05jx47b/8m+mRBATo6hnL0rLy+Hh4cHfH19rfYcDRo0qNf0RNeuXUN0dDRmz55twagMmCwSVcORdnZENlGxKtOtc6hWrMrEhNFp5Jo565m55Wpj1apViIyMhIeHB5o0aYIBAwbg2rVrxseXL1+Otm3bQq1WIzw8HIsWLTI+dvr0aSgUCnz77bfo168f1Go1VqxYUeVh6HXr1iEqKgpqtRphYWGYPn06tFqt8fGkpCQEBwfD3d0dAQEBmDRpksmY63sYesyYMZg6dSoGDBhQ5zpMYbJIVA0pd3ZEslPjqkwwrMrEQ9JOoYZFQWpdzly5ubl45JFH8OSTT+Lo0aPYunUrhg4dCnHjcNCSJUvwxhtvYObMmTh69CjeeecdvPnmm/j8888r1TN58mRMmjQJR48eRWxs7G3Ps2HDBjz22GOYNGkSMjIysHjxYiQnJ2PmzJkADAnr/PnzsXjxYpw4cQLff/89IiMjzW7HypUr0aBBg2pvK1eurMcrZT5OnUNUDal2dkSyxFWZ6Ca9ewOBgYbzu6s6lUehMDze28Kzo+Xm5kKr1WLo0KEICQkBgEpJ2owZMzBv3jwMHWo4Lz40NNSY7I0dO9ZYLiEhwVimKjNnzsSUKVOM24SFhWHGjBl47bXXMG3aNGRnZ8PPzw8DBgyASqVCcHAwunbtanY74uPj0a1bt2rLNG/e3Oz66oPJIlE1pNrZEckSV2WimyiVwPvvGy4EVCgEhPhnirSKuaAXLLD8fIsdOnRA//79ERkZidjYWMTExGDYsGFo1KgRzp8/j5ycHIwfPx4TJkwwbqPVaqHRaCrV07lz52qfZ//+/di7d69xJBEwLI9YUlKC4uJiDB8+HAsWLEBYWBgGDhyIQYMGYciQIXB1NS/18vb2hre3dy1abj08DE1UjYqdHWDY2d3Mmjs7Ilniqkx0i6FDgVWrgICAyvcHBhrur2bgrs6USiVSU1Px888/o127dli4cCHatGmDrKws6PV6AIZD0QcPHjTeDh8+jN27d1eqx8vLq9rn0ev1mD59eqV60tPTceLECajVagQFBSEzMxMfffQRPDw88Pzzz6NPnz41Lq1XQdaHod99911MnjwZhw4dQtu2baFSqawRF5HdqNjZTZpkGGGsEBhoSBStsbMjkiWuyuTwhBAQ16/XapsHBwL9/9Ai4vmD0F5VY9lzd2FAXwWUSkBfbH49Cg8Ps5enUygUiI6ORnR0NKZOnYqQkBCsWbMGiYmJaNGiBU6dOoVHH320Vu24VadOnZCZmYk77rjDZBkPDw/Ex8cjPj4eEydORHh4ONLT09GpU6ca65f1YehevXoBMFy1c/ToUahUKkRERCAyMhKRkZHo0qWLzYInspWhQ4H74nS4c/xe6K6q8eWL7XHfvUqOKBLdrGJVph3DIKCAolLCaH+rMlHtievXkdkpqk7bbqz45VXgzzps3+bAfig8PWsst2fPHmzatAkxMTHw9fXFnj17cP78ebRt2xaAIX+ZNGkSfHx8EBcXh9LSUuzbtw+XL19GYmKi2fFMnToVgwcPRlBQEIYPHw4XFxccOnQI6enpePvtt5GcnAydTodu3brB09MTX375JTw8PIznUdaktoehL126hOzsbJw9exYAkJmZCQDw8/ODn5+f2fVUpdaHoaOjowEYJn88evQodu/ejYSEBPj6+uKXX37BoEGD8Oabb9YrKCJ7pFQC6uBL8Gp3Fn36CiaKRFWRyapM5Lh8fHywfft2DBo0CK1bt8a///1vzJs3D3FxcQCAp556Cp999hmSk5MRGRmJvn37Ijk5GaGhobV6ntjYWPzwww9ITU1Fly5d0L17d7z33nvGZLBhw4ZYsmQJoqOj0b59e2zatAnr1q2r11yK1Vm7di06duyI+++/HwAwatQodOzYEZ988km961YIUdVp+/UTFRWF/fv3W7paiyksLIRGo0FBQQF8fHzqVZeuXIf0zTtQfDEXnk38EXlvbyhVjp9FWHpR+KLrWkRO3wAAWDwgFgPucbW7ZKy4tBRPznkPvqrLmPNYLNQB/exuhMTc19HS/Scpvc5wde31XMO5cM16212/WJwM2mzpz0txmRbtphre2xlvxcLTzfGvz7TG90ttv/9KSkqQlZWF0NBQqNXqOh2GtpTaHIYm89zav6ZY5dN260mijmr3dykIPv8S7m54Y6qIi8DZJYHIbvY+ug933P+eU1IMq5rcPFl1YKDhQpC6nL+XkgK8lAgoRxn+jhsEtGhe9/qsIicF6n0v4etWNxq9fa5hpCTqfbsZKTH3dbR0/0kqJ8Uwr9/N07XYWb9YnFzafGNVJgCYbUerMsmFvX6/KBQKsw4Fk2Ox2NXQubm5KC0tBQCbXPSyaNEiYyYcFRWFHTZeQmP3dynoWjYMfprKc4r5+ZxB17Jh2P2dY65SYOl1kivqO3PWMvVZxY0VKRTX7XdFCnNfR4da59oZVwpxxjY7IWf9fiH7ZbGRxTFjxuDkyZN4+OGH8Z///MdS1Vbpm2++QUJCAhYtWoTo6GgsXrwYcXFxyMjIQHBwsFWfGzAcGgg+/xKgEVAA0JdVHhbX6xUIuvASSs/1gtLVcf6b1umAqS8BDRUwnqt+M4UCmJoA3N/LvEOalepz0cFda/hno6HLJUChrHV9ViF0UPw2CSjDbafrA4CAAtj5EkT/XoBCmiDNfR0H9rBs/0nqln657WE76BeLk1mbtWX/vBe1Fy9B61a/mLRlWviUFN2o7yK0DnoYWqfVIfjsJOh9ABeFAvob9ytUAi4uwvD9cj4BuvIHnOKUJ7IPFj9n8dixYwgPD7dklbfp1q0bOnXqhI8//th4X9u2bfHggw9i1qxZt5UvLS01jnoChnM2goKC6nzO4sENW3H3xXsAGBLFzNlt6tAKIiIi87SZkgkXt3++rg822YK7Y/vVup76nrNIjsXc/q31Yej3b8xQnJmZaZzc8mbWThTLysqwf/9+xMTEVLo/JiYGO3furHKbWbNmQaPRGG9BQUH1iqH4IlcfICIi6fB7iGyp1uP4ERERAICXX34Zf/75Jxo0aIC77roLERERiIiIMF6ybS0XLlyATqe7bS7H5s2bIy8vr8ptXn/99UpzJ1WMLNaVZxN/4KLhd4VKoM2UzCrLHfJeg4h+ver8PPbm19+AoQ/d+EOlQ7MXfgUAnP+wF1D+z+GQlDVAr+ha1lcNc+uzigu/wmVXzUHqe6wBmkrT1+a+jjNmAMZZrSzQf5KSQb9YnMzaXFymRa93twAAfp18T72vXi4u06HLPMN7du8rveBZz8Pa9urw1l/Rvuj2flaoKh8E9GzCVXDIdmr96e3fvz8A4KeffgJgSLwOHz6Mw4cPIzU11erJYoVbL58XQpi8pN7d3R3u7u4We+7Ie3vj7JJA+PmcgYuLgMKt8odYr1cgtzAQkY8OdqhzSnrfD3gF3FjFRK+Fj6vhNb2ibwyhd/1nneT7zTvn7eb6ql132cz6rKLRYOBoixpXpHC5c7BkV3ua+zpOeBWY/bHl+k9SMugXi5NZm13LtChUGyYUdm3SBK71TBZdy7QovfGedW3SuN712avI+MHIW9LC+P1yK+P3y3CugkO2Y/ZhaJ1Oh8WLF+Pll1/GvHnzsHnzZly8eBE+Pj7o2bMnnn76aSxYsMCKoRo0bdoUSqXytlHE/Px8m60co1Qpkd3sfUBh+ODeTK9XAAogp9kCh0oUgVvXSa78WF3WSbZ0fVZRsSIFgNuvCrGPFSnMfR3d3GTwepvrpn4RdtovFueMbXZCDvv9otcB57YCp78y/NTrpI6IasHsZPHFF1/Em2++ifz8fLz++usYNGgQfH19ERwcjPj4eGvGWImbmxuioqKQmppa6f7U1FT07NnTZnF0Hz4Uv7utQl5h5VUKcgsD8bvbKoedZ9G4KPwtR0Dquih8RX0tblnswZqLzNfajRUp4Gm/K1KY+zpauv8k5YwrhThjm52Qw32/5KQAa1sCm+4Bdo42/FzbklM9yYjZV0P7+fnh888/R2xsLLy9vbFr1y5s27YN06dPx8iRI7Fw4UJrx2r0zTffYMyYMfjkk0/Qo0cPfPrpp1iyZAmOHDli1pqLXMGl/iy94oosVhSRwaoZ5r6Oclgxx1xyWFnH0uTQZkuvuMIVXOxjBZdaq5gb9LZTJ26Mmsron5zk5GQkJCTgypUrVqk/KSkJ33//PQ4ePGiV+qti8RVcrl69inbt2gEwTLqtVCoxceJElJWVGRettpWRI0fi4sWLeOutt5Cbm4uIiAj89NNPZi/ObUlKlbJO0xfI3c2JRe8+9U/slEqgX7/61WF1LkqgeT+po6iWua+jpftPUs64Uogztlno0N3rEHxVl+GSrwbsMEG2NNl/v+h1htWGqjzHVgBQAPsTgBYPyKIvR44ciUGDBlmt/ldffRUvvvhinbYtLy/Hv//9b/z00084deoUNBoNBgwYgNmzZyMgIKDesZl9GDosLMyYFLZo0QJnzpwBAAwZMgQrVqyodyC19fzzz+P06dMoLS3F/v370adPH5vHQERENpCTAvVPd+DrVv/CB8Fzod4+gIcx5eD8jttXG6pEAMU5hnJ2rry8HB4eHvD19bXaczRo0ABNmjSp07bFxcU4cOAA3nzzTRw4cAApKSk4fvy4xU4TNDtZHD58ONavXw8A6NevH5YtWwYAyMjIwHWJFhUnItLpgJLsxriWEYDt2xTQ8bx5xyKD5TbJhOtmzgVpbrlaWLVqFSIjI+Hh4YEmTZpgwIABuHbtmvHx5cuXo23btlCr1QgPD8eiRYuMj50+fRoKhQLffvst+vXrB7VajRUrViA5ORkNGzas9Dzr1q1DVFQU1Go1wsLCMH36dGi1WuPjSUlJCA4Ohru7OwICAjBp0iSTMSclJeHuu++uU3s1Gg1SU1MxYsQItGnTBt27d8fChQuxf/9+ZGdn16nOm5l9GPpN4wRtwP/93/+ha9euaNasGQoLCzF+/Ph6B0JEVFspKcCkSUqcO9MDABC3znCxzvvvy+xiHaraTYcxq1rgUG6HMZ2Oh5lzQZpbzky5ubl45JFHMGfOHDz00EMoKirCjh07UHGJxpIlSzBt2jR8+OGH6NixI9LS0jBhwgR4eXlh7NixxnomT56MefPmYfny5XB3d8fGjRsrPc+GDRvw2GOP4YMPPkDv3r1x8uRJPP300wCAadOmYdWqVZg/fz6+/vpr3HXXXcjLy8Mff/xhdjtWrlyJZ555ptoyixcvxqOPPlrlYwUFBVAoFLcluHVRpzOEg4ODceTIEfz0009o3LixzeZWJCKqkJICDBt2+9ySZ84Y7pfd1d10u9ocxrTz84mdUrPehiv1a5gbFM0sO2dkbm4utFothg4daryWITIy0vj4jBkzMG/ePAy9sYMIDQ1FRkYGFi9eXClZTEhIMJapysyZMzFlyhTjNmFhYZgxYwZee+01TJs2DdnZ2fDz88OAAQOgUqkQHByMrl27mt2O+Ph4dOvWrdoypqYMLCkpwZQpUzB69Oh6X8gL1DFZBIAmTZpgzJgx9Q6AiKi2dDrgpZcqEsVbJ+g3zBuZkAA88IDML94xoeLQu+6qGtu3KXDfvfbXTovEKOFhTLKAirlBdwyDgAKKSgmj9eYG7dChA/r374/IyEjExsYiJiYGw4YNQ6NGjXD+/Hnk5ORg/PjxmDBhgnEbrVYLjUZTqZ7OnTtX+zz79+/H3r17MXPmTON9Op0OJSUlKC4uxvDhw7FgwQKEhYVh4MCBGDRoEIYMGQJXV/NSL29vb3h7e9ei5Qbl5eUYNWoU9Hp9pcPr9VHrtaGJiKS2YwfwdzUDTkIAOTmGco4mJQVoe6cS577qgQvrOiIuRomWLQ332wuLxSjRYUyyIAnmBlUqlUhNTcXPP/+Mdu3aYeHChWjTpg2ysrKg1+sBGA5FHzx40Hg7fPgwdu/eXakeLy+vap9Hr9dj+vTplepJT0/HiRMnoFarERQUhMzMTHz00Ufw8PDA888/jz59+qC8vNysdqxcuRINGjSo9rZy5cpK25SXl2PEiBHIyspCamqqRUYVgXqMLBIRSSXXzIEkc8vJhRwOvVs0RokOY5JpQgiI2l7U2mQgivsdwnPvL0Qz1WXMHDUA7n59AIUSKC42uxqFh4fJZX1vK6tQIDo6GtHR0Zg6dSpCQkKwZs0aJCYmokWLFjh16pTJc/3M1alTJ2RmZuKOO+4wWcbDwwPx8fGIj4/HxIkTER4ejvT0dHTq1KnG+mt7GLoiUTxx4gS2bNlS5yurq8JkkYhkx9/MgSRzy8mBHA69WzzGmw5jGuqzzWFMMk1cv47MTlF12vZfN36e/nJ9nbZvc2A/FJ6eNZbbs2cPNm3ahJiYGPj6+mLPnj04f/482rZtC8Bw1fGkSZPg4+ODuLg4lJaWYt++fbh8+TISExPNjmfq1KkYPHgwgoKCMHz4cLi4uODQoUNIT0/H22+/jeTkZOh0OnTr1g2enp748ssv4eHhYfac0LU5DK3VajFs2DAcOHAAP/zwA3Q6nXFZ5MaNG8PNzc3sdlWFh6GJSHZ69zZc9WxqkEGhAIKCDOUchRwOvVslRhkst0n2xcfHB9u3b8egQYPQunVr/Pvf/8a8efMQFxcHAHjqqafw2WefITk5GZGRkejbty+Sk5MRGhpaq+eJjY3FDz/8gNTUVHTp0gXdu3fHe++9Z0wGGzZsiCVLliA6Ohrt27fHpk2bsG7dOouO+FX4+++/sXbtWvz999+4++674e/vb7zt3Lmz3vWbvdyfI7Hkcn/OSg7LjZFpjrB0WsXhTqDyIc+KBNIeDsla0ldfAaNH11zuv/8FHnnE+vFUxaoxymC5TTmo73J/dToMbSG1OQxN5rH4cn9ERjkpUO97CV+3ujGEsH2u4b/8qPf5Xz7ZzNChhoTwpZcqj2YFBgILFjhWogjI49C7VWOUwXKbzkChUJh1KJgcC5NFqp2K1RRuPdm8YjUFHhYiGxo61HD+244dhotZ/P0Nh57tbRoZS6g49H7mzO0XjwCGEdXAQGkPvcshRiKqPZ6zSOarcTUFGFZT0HO9NbIdpRLo189wWLNfP8dMFAFDu95/3/D7rUfiKv5esEDa9sshRiKqPSaLZD4HWhSeSI4qDr23uOVaj8BA+zlHUw4xElHt8DA0mY+rKRBJTg6H3uUQI9VOxWTW5FjM7Vcmi2Q+rqZAZBcqDr3bMznESDVzc3ODi4sLzp49i2bNmsHNzY1XJDsAIQTKyspw/vx5uLi41DgPI5NFMh9XUyAiciouLi4IDQ1Fbm4uzp49K3U4ZGGenp4IDg6Gi0v1ZyUyWSTzcTUFIiKn4+bmhuDgYGi1Wuh0vIDRUSiVSri6upo1UsxkkWqnYjWF/S9VvtjFM9CQKHLaHCIih6NQKKBSqaBSqaQOhSTAZJFqL2go0OIBrqZARETkBJgsUt1wNQUiIiKnwHkWiYiIiMgkJotEREREZBKTRSIiIiIyickiEREREZnEZJGIiIiITGKySEREREQmMVkkIiIiIpOYLBIRERGRSZyUm8gZCR26ex2Cr+oyXPLVQEA/rsBDRERVktXI4unTpzF+/HiEhobCw8MDrVq1wrRp01BWViZ1aETykZMC9U934OtW/8IHwXOh3j4AWNsSyEmROjIiIrJDshpZPHbsGPR6PRYvXow77rgDhw8fxoQJE3Dt2jX85z//kTo8IvuXkwLsGAYFROX7i88AO4YBvVcZ1v4mIiK6QSGEEDUXs19z587Fxx9/jFOnTpm9TWFhITQaDQoKCuDj42PF6IjsiF5nGEEs/ttEAQXgGQjEZ/GQNJGD4vcf1YWsRharUlBQgMaNG1dbprS0FKWlpca/CwsLrR0Wkf05v6OaRBEABFCcYyjXvJ+toiIiIjsnq3MWb3Xy5EksXLgQzz77bLXlZs2aBY1GY7wFBQXZKEIiO3I917LliIjIKdhFspiUlASFQlHtbd++fZW2OXv2LAYOHIjhw4fjqaeeqrb+119/HQUFBcZbTk6ONZtDZJ88/C1bjoiInIJdnLN44cIFXLhwodoyLVu2hFqtBmBIFO+55x5069YNycnJcHGpXc7LczbIKRnPWTwD3HqBCwCes0jk+Pj9R3VhF+csNm3aFE2bNjWr7JkzZ3DPPfcgKioKy5cvr3WiCAAV+THPXSSn03oWsHOMiQcFcPc7wNVrNg2JiGyn4nvPDsaJSEbsYmTRXGfPnkXfvn0RHByML774AkrlP6Mffn5+Ztfz999/87xFIiJyWjk5OQgMDJQ6DJIJWSWLycnJeOKJJ6p8rDbN0Ov1OHv2LLy9vaFQKOodV2FhIYKCgpCTkyP7YX22xT6xLfaJbbFPbItpQggUFRUhICCgTkfmyDnZxWFoc40bNw7jxo2rdz0uLi5W+Y/Kx8dH9jumCmyLfWJb7BPbYp/YlqppNBqL1EPOg/9WEBEREZFJTBaJiIiIyCQmixbg7u6OadOmwd3dXepQ6o1tsU9si31iW+wT20JkWbK6wIWIiIiIbIsji0RERERkEpNFIiIiIjKJySIRERERmcRkkYiIiIhMYrJopkWLFiE0NBRqtRpRUVHYsWNHteW3bduGqKgoqNVqhIWF4ZNPPrFRpKbNmjULXbp0gbe3N3x9ffHggw8iMzOz2m22bt0KhUJx2+3YsWM2irpqSUlJt8VU05KP9tgnANCyZcsqX+OJEydWWd6e+mT79u0YMmQIAgICoFAo8P3331d6XAiBpKQkBAQEwMPDA/369cORI0dqrHf16tVo164d3N3d0a5dO6xZs8ZKLfhHdW0pLy/H5MmTERkZCS8vLwQEBODxxx/H2bNnq60zOTm5yr4qKSmRrC2AYYGDW2Pq3r17jfXaW78AqPL1VSgUmDt3rsk6peoXc/bBcvrMkPNgsmiGb775BgkJCXjjjTeQlpaG3r17Iy4uDtnZ2VWWz8rKwqBBg9C7d2+kpaXhX//6FyZNmoTVq1fbOPLKtm3bhokTJ2L37t1ITU2FVqtFTEwMrl27VuO2mZmZyM3NNd7uvPNOG0RcvbvuuqtSTOnp6SbL2mufAMDevXsrtSM1NRUAMHz48Gq3s4c+uXbtGjp06IAPP/ywysfnzJmD9957Dx9++CH27t0LPz8/3HfffSgqKjJZ565duzBy5EiMGTMGf/zxB8aMGYMRI0Zgz5491moGgOrbUlxcjAMHDuDNN9/EgQMHkJKSguPHjyM+Pr7Gen18fCr1U25uLtRqtTWaYFRTvwDAwIEDK8X0008/VVunPfYLgNte22XLlkGhUODhhx+utl4p+sWcfbCcPjPkRATVqGvXruLZZ5+tdF94eLiYMmVKleVfe+01ER4eXum+Z555RnTv3t1qMdZFfn6+ACC2bdtmssyWLVsEAHH58mXbBWaGadOmiQ4dOphdXi59IoQQL730kmjVqpXQ6/VVPm6vfQJArFmzxvi3Xq8Xfn5+Yvbs2cb7SkpKhEajEZ988onJekaMGCEGDhxY6b7Y2FgxatQoi8dsyq1tqcrvv/8uAIi//vrLZJnly5cLjUZj2eBqqaq2jB07VjzwwAO1qkcu/fLAAw+Ie++9t9oy9tAvQty+D5bzZ4YcG0cWa1BWVob9+/cjJiam0v0xMTHYuXNnldvs2rXrtvKxsbHYt28fysvLrRZrbRUUFAAAGjduXGPZjh07wt/fH/3798eWLVusHZpZTpw4gYCAAISGhmLUqFE4deqUybJy6ZOysjKsWLECTz75JBQKRbVl7bFPbpaVlYW8vLxKr7u7uzv69u1r8rMDmO6r6raRQkFBARQKBRo2bFhtuatXryIkJASBgYEYPHgw0tLSbBNgDbZu3QpfX1+0bt0aEyZMQH5+frXl5dAv586dw48//ojx48fXWNYe+uXWfbCjf2ZIvpgs1uDChQvQ6XRo3rx5pfubN2+OvLy8KrfJy8ursrxWq8WFCxesFmttCCGQmJiIXr16ISIiwmQ5f39/fPrpp1i9ejVSUlLQpk0b9O/fH9u3b7dhtLfr1q0bvvjiC2zYsAFLlixBXl4eevbsiYsXL1ZZXg59AgDff/89rly5gnHjxpksY699cquKz0dtPjsV29V2G1srKSnBlClTMHr0aPj4+JgsFx4ejuTkZKxduxZfffUV1Go1oqOjceLECRtGe7u4uDisXLkSmzdvxrx587B3717ce++9KC0tNbmNHPrl888/h7e3N4YOHVptOXvol6r2wY78mSF5c5U6ALm4dZRHCFHtyE9V5au6XyovvPACDh06hF9//bXacm3atEGbNm2Mf/fo0QM5OTn4z3/+gz59+lg7TJPi4uKMv0dGRqJHjx5o1aoVPv/8cyQmJla5jb33CQAsXboUcXFxCAgIMFnGXvvElNp+duq6ja2Ul5dj1KhR0Ov1WLRoUbVlu3fvXunCkejoaHTq1AkLFy7EBx98YO1QTRo5cqTx94iICHTu3BkhISH48ccfq0207LlfAGDZsmV49NFHazz30B76pbp9sKN9Zkj+OLJYg6ZNm0KpVN72H1p+fv5t/8lV8PPzq7K8q6srmjRpYrVYzfXiiy9i7dq12LJlCwIDA2u9fffu3SUfGbmVl5cXIiMjTcZl730CAH/99Rd++eUXPPXUU7Xe1h77pOLq9Np8diq2q+02tlJeXo4RI0YgKysLqamp1Y4qVsXFxQVdunSxu77y9/dHSEhItXHZc78AwI4dO5CZmVmnz4+t+8XUPtgRPzPkGJgs1sDNzQ1RUVHGK1QrpKamomfPnlVu06NHj9vKb9y4EZ07d4ZKpbJarDURQuCFF15ASkoKNm/ejNDQ0DrVk5aWBn9/fwtHVz+lpaU4evSoybjstU9utnz5cvj6+uL++++v9bb22CehoaHw8/Or9LqXlZVh27ZtJj87gOm+qm4bW6hIFE+cOIFffvmlTv9kCCFw8OBBu+urixcvIicnp9q47LVfKixduhRRUVHo0KFDrbe1Vb/UtA92tM8MORAprqqRm6+//lqoVCqxdOlSkZGRIRISEoSXl5c4ffq0EEKIKVOmiDFjxhjLnzp1Snh6eoqXX35ZZGRkiKVLlwqVSiVWrVolVROEEEI899xzQqPRiK1bt4rc3Fzjrbi42Fjm1rbMnz9frFmzRhw/flwcPnxYTJkyRQAQq1evlqIJRq+88orYunWrOHXqlNi9e7cYPHiw8Pb2ll2fVNDpdCI4OFhMnjz5tsfsuU+KiopEWlqaSEtLEwDEe++9J9LS0oxXCM+ePVtoNBqRkpIi0tPTxSOPPCL8/f1FYWGhsY4xY8ZUmlngt99+E0qlUsyePVscPXpUzJ49W7i6uordu3dL1pby8nIRHx8vAgMDxcGDByt9fkpLS022JSkpSaxfv16cPHlSpKWliSeeeEK4urqKPXv2SNaWoqIi8corr4idO3eKrKwssWXLFtGjRw/RokUL2fVLhYKCAuHp6Sk+/vjjKuuwl34xZx8sp88MOQ8mi2b66KOPREhIiHBzcxOdOnWqNN3M2LFjRd++fSuV37p1q+jYsaNwc3MTLVu2NLkTsyUAVd6WL19uLHNrW959913RqlUroVarRaNGjUSvXr3Ejz/+aPvgbzFy5Ejh7+8vVCqVCAgIEEOHDhVHjhwxPi6XPqmwYcMGAUBkZmbe9pg990nFND633saOHSuEMEwFMm3aNOHn5yfc3d1Fnz59RHp6eqU6+vbtayxf4bvvvhNt2rQRKpVKhIeH2yQRrq4tWVlZJj8/W7ZsMdmWhIQEERwcLNzc3ESzZs1ETEyM2Llzp6RtKS4uFjExMaJZs2ZCpVKJ4OBgMXbsWJGdnV2pDjn0S4XFixcLDw8PceXKlSrrsJd+MWcfLKfPDDkPhRA3zvInIiIiIroFz1kkIiIiIpOYLBIRERGRSUwWiYiIiMgkJotEREREZBKTRSIiIiIyickiEREREZnEZJGIiIiITGKySEREREQmMVkkIiIiIpOYLBIRERGRSa5SByAFvV6Ps2fPwtvbGwqFQupwiIiIbEIIgaKiIgQEBMDFheNFZB6nTBbPnj2LoKAgqcMgIiKSRE5ODgIDA6UOg2RCdsnirFmzkJKSgmPHjsHDwwM9e/bEu+++izZt2phdh7e3NwDDh8XHx8daoRIREdmVwsJCBAUFGb8Hicwhu2Rx27ZtmDhxIrp06QKtVos33ngDMTExyMjIgJeXl1l1VBx69vHxYbJIREROh6dgUW0ohBBC6iDq4/z58/D19cW2bdvQp0+fKsuUlpaitLTU+HfFf1YFBQVMFsmxlF8HfnsfuHqu5rLu3kCPF4AGvtaPi4jsQmFhITQaDb//qFZkN7J4q4KCAgBA48aNTZaZNWsWpk+fbpXn12q1WL9+PQBg4MCBcHWV/UtqFku3Ww6voyxizNyI9VcjAERg4IEn4aovrbqcizvWd1oGbP3dbttiLjn0i6XJoc3OuI+wNGdsM9knWb/zhBBITExEr169EBERYbLc66+/jsTEROPfFSOLRA6n7BqAG6MFvV8FFCYOHJzcYrOQiIhI3mSdLL7wwgs4dOgQfv3112rLubu7w93d3UZREUlI6P/5vXciYGokouQqoLVNSEREJG+yTRZffPFFrF27Ftu3b+fl/0QV9DrzyvHkdiIiMpPskkUhBF588UWsWbMGW7duRWhoqNQhEdmPm0cWq6PgZLxERGQe2SWLEydOxH//+1/873//g7e3N/Ly8gAAGo0GHh4eEkdHJDUmi0REZFmy+8b4+OOPUVBQgH79+sHf3994++abb6QOjUh65s6EpVBaNw4iInIYshtZtLdpIRUKBXx9fY2/OwtLt1sOr6MsYhQ6+F5JA7z9qo1RoVAYyjUMttu2mEsO/WJpcmizM+4jLM0Z20z2SfaTctcFJyUlh7X7Y2D9FCDiYWDYMtPltswCts0GOo8HBr9nu/iISFL8/qO6kN1haCKqRsUFLjWdk1jxuLkXxBARkdNiskjkSIzJYg3nJLowWSQiIvPI7pxFe6PVapGamgoAuO+++5xmOSZLt1sOr6MsYtTpkdpxKeDiivu0WpMxaqE0lFNUX04O5NAvliaHNjvjPsLSnLHNZJ/4zrMAnc7MiZAdjKXbLYfX0e5jFHrolOqayylczCsnE3bfL1YghzY74z7C0pyxzWR/eBiayJFwUm4iIrIwfmMQORLOs0hERBbGZJHIkQhz14bmR5+IiMzDbwwiR8LD0EREZGG8wIXIkZh9GJrJIhHVjk6nQ3l5udRhkIWoVCooleadksRksZ4UCgUaN25s/N1ZWLrdcngdZRGj0KFxUQbg3aL65f5cFIZynk3tti3mkkO/WJoc2uyM+whLs5c2CyGQl5eHK1euSBYDWUfDhg3h51f98rAAl/vjckfkWFKnAb8tALpPBAa+Y7rcvmXADy8D4YOBUSttFh4RSasu33+5ubm4cuUKfH194enp6TTJuiMTQqC4uBj5+flo2LAh/P39qy3PkUUiR2JcwaWGnTmX+yMiM+h0OmOi2KRJE6nDIQvy8PAAAOTn58PX17faQ9I8cYnIkVQkfy41nIdSMXUOk0UiqkbFOYqenp4SR0LWUNGvNZ2LypHFetJqtdi8eTMA4N5773Wa5Zgs3W45vI6yiFGvwOYOHwOlHri3uuX+hIuhnNKt2nJyIId+sTQ5tNkZ9xGWZk9t5qFnx2Ruvzr+p80GysrKpA5BEpZutxxeR7uPUehQpjLjPCSFC8pU3taPx0bsvl+sQA5tdsZ9hKU5Y5vJ/vAwNJEj4TyLRERWkZycjIYNG1qt/qSkJNx9991Wq78++I1B5EjMTRZrOqeRiIgqGTlyJI4fP261+l999VVs2rSpztunpKQgNjYWTZsapkQ7ePCgxWJjskjkUMwdWeT5R0RE5iovL4eHhwd8fX2t9hwNGjSo1xXn165dQ3R0NGbPnm3BqAyYLBI5Ej0PQxMRrVq1CpGRkfDw8ECTJk0wYMAAXLt2zfj48uXL0bZtW6jVaoSHh2PRokXGx06fPg2FQoFvv/0W/fr1g1qtxooVK6o8DL1u3TpERUVBrVYjLCwM06dPh1arNT6elJSE4OBguLu7IyAgAJMmTTIZc30PQ48ZMwZTp07FgAED6lyHKbzAhciR8JxFIrIiIQSul+skeW4PldKsq3dzc3PxyCOPYM6cOXjooYdQVFSEHTt2oGINkiVLlmDatGn48MMP0bFjR6SlpWHChAnw8vLC2LFjjfVMnjwZ8+bNw/Lly+Hu7o6NGzdWep4NGzbgsccewwcffIDevXvj5MmTePrppwEA06ZNw6pVqzB//nx8/fXXuOuuu5CXl4c//vjD7PauXLkSzzzzTLVlFi9ejEcffdTsOuuKyWI9KRQKaDQa4+/OwtLtlsPrKIsYhR6aaycBnxqW+1MoDeXcGthtW8wlh36xNDm02Rn3EZZmj22+Xq5Du6kbJHnujLdi4elWc9qSm5sLrVaLoUOHIiQkBAAQGRlpfHzGjBmYN28ehg4dCgAIDQ1FRkYGFi9eXClZTEhIMJapysyZMzFlyhTjNmFhYZgxYwZee+01TJs2DdnZ2fDz88OAAQOgUqkQHByMrl27mt3e+Ph4dOvWrdoyzZs3N7u++mCyWE9KpRK9e/eWOgybs3S75fA6yiJGoUXvo1OB+94CqpmNX6l0MZQL7AIox9swQsuTQ79Ymhza7Iz7CEtzxjZbQocOHdC/f39ERkYiNjYWMTExGDZsGBo1aoTz588jJycH48ePx4QJE4zbaLVaY2JeoXPnztU+z/79+7F3717MnDnTeJ9Op0NJSQmKi4sxfPhwLFiwAGFhYRg4cCAGDRqEIUOGmD1fpre3N7y97WOKMyaLRLVwvqgUOZeLpQ7DpJCr19EEwN9XSpGffdlkOZ8LxbgDwLWSMmRWU04uGri74k5f+Y+SEtk7D5USGW/FSvbc5lAqlUhNTcXOnTuxceNGLFy4EG+88Qb27NljXLFkyZIlt43a3brcnZeXV7XPo9frMX369CpHH9VqNYKCgpCZmYnU1FT88ssveP755zF37lxs27YNKpWqxnbwMDSRDF0pLkPvOZtRUm6/S+QtUJ3Dg0pg+c6/sHTHTpPl+rkcQ7IbcDK/EEMXmS4nJ3Mebo8RXYKkDoPIoSkUCrMOBUtNoVAgOjoa0dHRmDp1KkJCQrBmzRokJiaiRYsWOHXqVL2TrE6dOiEzMxN33HGHyTIeHh6Ij49HfHw8Jk6ciPDwcKSnp6NTp0411s/D0A5Ep9Nh69atAIB+/fpVuxC3I7F0u+XwOuZcuorXOxhO7P48yxNaYX+jWF5l7tjUdgEiFA0QdtLDZIzN9F7YdOcClMMVYadMl5ODqyWleKmdFq7njkCnC7DL946lyeHz4oz7CEtzxjZbwp49e7Bp0ybExMTA19cXe/bswfnz59G2bVsAhquOJ02aBB8fH8TFxaG0tBT79u3D5cuXkZiYaPbzTJ06FYMHD0ZQUBCGDx8OFxcXHDp0COnp6Xj77beRnJwMnU6Hbt26wdPTE19++SU8PDyM51HWpLaHoS9duoTs7GycPXsWAJCZmQkA8PPzg5+fn9n1VMXqyaJSqYROJ82VU7YghMD169eNvzsLS7dbDq9jqVaPJmrD7xtf7mOXa9Nqv1uB9e7NoEL1MWpP6LE+0/B622tbzDVv/VE00Z4EoLXb946lyeHzInTaf2KcHwmI+i1bJxRuuB4+x/C7nbbZ0uTQz/bIx8cH27dvx4IFC1BYWIiQkBDMmzcPcXFxAICnnnoKnp6emDt3Ll577TV4eXkhMjISCQkJtXqe2NhY/PDDD3jrrbcwZ84cqFQqhIeH46mnngIANGzYELNnz0ZiYiJ0Oh0iIyOxbt26es2lWJ21a9fiiSeeMP49atQoAIYrs5OSkupVt9W/IfgGJ0dRrrXfw89GejP/MXOgqXPcXV0Abc3lyMYun/7n9+ILgL60fvW5uNdve3Iabdu2xfr166stM3r0aIwePbrKx1q2bFll7jJu3DiMGzeu0n2xsbGIja36HM4HH3wQDz74oFkxA4YRz/okdVXFZylWTxZ5wjk5ilI5JItOOM+im6vjtMWhlN90IdhTmwBlPfup6Dxw7Krhdw5CENlUvZPFGTNmwNXVFR06dECHDh3QokWLasuXlpbC3Z3/IZL8lMkiWTTzS1ThOOc+MVm0U9qbRhJ92wL1PdXB6wJwbLfhd105YMbVpERkGXXey7744ou4cOEC3nzzTYwePRparRbLli0zzl5uSs+ePW+7z5oLcxNZSplWBufeOuXIouMkvg5Fe92y9Sk9/vldV2LZuomoWnX+V69fv36IjY3FQw89hFdeecV4abgpP/zwA44dO4Zr167h7NmzCAgIMD42fPjwWi2BQySFMq0eblIHURMnTBbdlTzVxS6VWzihc73p02fpuomoWnX+xnj44YexZ88eaDQadO/eHZ999lm1F7Pcdddd8PT0RH5+PkaNGoWwsDD06dMHo0aNkvV0AAqFAg0aNECDBs41IbCl2y2H17FUJ5B7DSjUutptjAqhQ4Prf6OBq66G5f5cDOVKz9ltW8zl5qpE7jXgSrl568Y6Ajl8XhTaEsN7THvJMvsIFxc0uH4GDa7/DYXWOZJFOfQzOYd6nUTi6uqKoUOHQqPRIDExEfPnz8fcuXMxaNCg28oqlUo8//zziIiIQJ8+fQAAZ86cQVZWFiIiIuoThqSUSiX69esndRg2Z+l2y+F1LNMBM9MUuL+9L0bb6T84SlGOfkcmAw99Wv1yf65KQzlNsOyX+3N3c8XMNAU6BvvgMTvtF0uTw+dFqSsxvMda9QeUj9e/PqUS/U7NAq5fBvr2q3+AMiCHfibnUOdkMS4uDhkZGQgKCkLXrl2xcOFCtG7dGh999BFSU1Mxf/78SuUff/xx/PXXXwgKCkL79u2Nt44dO9a4pA6RPSi9cc6ie32v6rSmiqlzajrMXPG4kMF5mDWouMCl1I5X1nFKFecsqjyqL1cbKk9Dsmjp8yGJqFq1/tbLzs4GALzzzjs4deoUfv31V7z33nt45JFHEBUVhWXLluHnn3++bbutW7ciKysLDz30EHJycvDnn3/i3//+Nxo2bIg2bdrUOvBFixYhNDQUarUaUVFR2LFjR63rIKqNiquh3VV2nCxWnLNY0yErY7Io/wTL/cYFLmU6+bfFoVScV+iqtlydFXXxnEUim6r1yKKpEcLIyEjjCOFPP/1kcvsvvvgCaWlpxr83btyI//73v7WK4ZtvvkFCQgIWLVqE6OhoLF682DjSGRwcXNsm1YtOpzMmqr1795b1+Ze1Yel2y+F1LCvX4o2OAg1VudDp7rLLGHVCiR13vQucdUPvdjqTMer0MJRzcUVvnelycqByEXijo4BKeRU6mbfFXHL4vOjKSwzvMXeNRd5jOp0OO4JeBPxL0busGPbXYsuTQz+Tc6h1slixTuX8+fOxdetWeHp64ttvv8XWrVvRqlUrZGZmIiwszOT2arUamZmZxtHEmJgY/Otf/6pVDO+99x7Gjx9vXFJnwYIF2LBhAz7++GPMmjWrtk2qFyEErl41TBS7/+dlcIFzTBarhwJX4Q8A+OGPM4BLPXdieh2Udv46NvjrCvxD2gEohzi8ph6Xh1mPuJaPq80CgfLqV08SCiWuegQafrfTtpir6cVi+Hs1BiDs9r1jaTd//uy1/0TuH7jqNdzwu4WWBL2qagaogJN7fkBh5sl612nvbu7nooJLaNi4mcQRObfk5GQkJCTgypUrVqk/KSkJ33//PQ4ePGiV+uujzucs1nWE8LPPPsPw4cNxzz33oH379jhy5EitnresrAz79+/HlClTKt0fExODnTt3VrlNaWkpSkv/mSC2sLCwVs9prk4H3oBrfZe0kgmtizvWd1oGAPi/VYdQpq/flXpuLgLv3ZiC015fx7td3LE+xNBm/O/5+i9fZg0u7kBLM8q53jQxvr22xUyBLu44fOO9aK/vHUu7+fNnt/3n4g50Gm6VqsNOrXS6fr5y/m8mixIbOXJklRfwWsqrr76KF198sU7blpeX49///jd++uknnDp1ChqNBgMGDMDs2bMrTVVYV3VOFus6QnjXXXfh999/x5o1a3DkyBEEBQVVeY6jKRcuXIBOp0Pz5s0r3d+8eXPk5eVVuc2sWbMwffp0s5+jrjLc2sNFlFv9eexBsPYv4+/dWzWBTtQvWVQqBICLAIDjqrawx7PP9IqbVowI6QX7XJDYzI90o5YAMg2/221bzCNuarOzfAZL9Tf1s932n/VWk/1DGQl3lT222bJu3ue4qRtIGAmVl5fDw8MDHh4WvGDrFhXTJNVFcXExDhw4gDfffBMdOnTA5cuXkZCQgPj4eOzbt6/esdX501yfEUK1Wo1HHnmkrk8N4PY1p4UQJuehev3115GYmGj8u7CwEEFBQfV6/qq0e3U9XOu7pJVMaJc/YPx92dgu9W63tqwM6zduBAC0nvgdXH1861WfNWi1WmRXLE7/6Lf1X77MGrRaoCLG6tz8WbHXtphJcVObneUzuD79DPDXjSM79tp/5r4X6yBv8Oe4v0OgVeq2Jzfvc3xbhEocjbysWrUK06dPx59//glPT0907NgR//vf/4zXVixfvhxz5sxBVlYWWrZsiUmTJuH5558HAJw+fRqhoaH45ptvsGjRIuzevRsff/wxFArFbYeh161bh6SkJBw5cgQBAQEYO3Ys3njjDeN+KCkpCcuWLcO5c+fQpEkTDBs2DB988EGVMdfnMLRGo0Fqamql+xYuXIiuXbsiOzu73tdz1HkPY+4IoSXOVblZ06ZNoVQqbxtFzM/Pv220sYK7uzvXo7a0+p6jeCv9TaMElq6byMG4KV3scizRVlw4QbV0hADKi6V5bpVnzTM9AMjNzcUjjzyCOXPm4KGHHkJRURF27NhhzEeWLFmCadOm4cMPP0THjh2RlpaGCRMmwMvLC2PHjjXWM3nyZMybNw/Lly+Hu7s7Nt4Y0KiwYcMGPPbYY/jggw/Qu3dvnDx50rjk8bRp07Bq1SrMnz8fX3/9Ne666y7k5eXVarW6lStX4plnnqm2zOLFi/Hoo49W+VhBQQEUCgUaNmxo9nOaUq9/R80ZIdTrLXtA0c3NDVFRUUhNTcVDDz1kvD81NRUPPPBANVuSRSlVNZepjUrJooXrJnIwKnue69MGlC5MFiVTXgy8U/9z4OrkX2cBt5rnZc7NzYVWq8XQoUMREhICAIiMjDQ+PmPGDMybNw9Dhw4FAISGhiIjIwOLFy+ulCwmJCQYy1Rl5syZmDJlinGbsLAwzJgxA6+99hqmTZuG7Oxs+Pn5YcCAAVCpVAgODkbXrl3Nbm58fDy6detWbRlTg2QlJSWYMmUKRo8eDR8fH7Of0xQ7PHZRs8TERIwZMwadO3dGjx498OmnnyI7OxvPPvuszWNRKBTGcxicaTkmhYsSHqXnAc/GllnKS2gN9QFQKO3zbSmHvjY3Rjm0xVyO1BZzqZQKXCwBXJUKu22zpftFoVCgqFyBMp1AoJPkys743raEDh06oH///oiMjERsbCxiYmIwbNgwNGrUCOfPn0dOTg7Gjx+PCRMmGLfRarXQaDSV6uncuXO1z7N//37s3bsXM2fONN6n0+lQUlKC4uJiDB8+HAsWLEBYWBgGDhyIQYMGYciQIWafKuPt7Q1vb+9atNygvLwco0aNgl6vx6JFi2q9fVXs81u5BiNHjsTFixfx1ltvITc3FxEREfjpp5+M/0HYklKpRP/+/W3+vFJTuijQPz0BGDy/2mXlzK4PwlAfADw0ttqyUpFDX5sboxzaYi5Haou53NxUmLZPgaDGHnjUTufes3S/KJVKfJntg4zcQiRHyvKrq9bs8r2t8jSM8En13GZQKpVITU3Fzp07sXHjRixcuBBvvPEG9uzZA09PQx1Lliy5bdTu1nksa1pdTq/XY/r06VWOPqrVagQFBSEzMxOpqan45Zdf8Pzzz2Pu3LnYtm0bVKqaj6DV5TB0eXk5RowYgaysLGzevNkio4qATJNFAHj++eeNJ6OSBFxuvHV0FjpzyngYWgG4OMmwAVEdud04DF2udfw5JW+mv3HOGQ9DS0ihMOtQsNQUCgWio6MRHR2NqVOnIiQkBGvWrEFiYiJatGiBU6dOmTzXz1ydOnVCZmYm7rjjDpNlPDw8EB8fj/j4eEycOBHh4eFIT09Hp06daqy/toehKxLFEydOYMuWLWjSpIn5jamBbJNFkljFOYt6SyWL5ZXrJSKTKs5ZLHeyJQ51eiaLVLM9e/Zg06ZNiImJga+vL/bs2YPz58+jbdu2AAxXHU+aNAk+Pj6Ii4tDaWkp9u3bh8uXL1eaOaUmU6dOxeDBgxEUFIThw4fDxcUFhw4dQnp6Ot5++20kJydDp9OhW7du8PT0xJdffgkPDw+zj4LW5jC0VqvFsGHDcODAAfzwww/Q6XTGC4EbN24MNzc3s9tVFSaL9aTT6YyTgffs2dNplmPSubhhZ9u3gCu+6GmJpbzKywz1KZQWqc8a5NDX5sYoh7aYy5HaYi6lQo//6yCgdCmz2yUOLd0vOp0OI4OuodRfwMUB1jQ3hzO+ty3Bx8cH27dvx4IFC1BYWIiQkBDMmzcPcXFxAICnnnoKnp6emDt3Ll577TV4eXkhMjISCQkJtXqe2NhY/PDDD3jrrbcwZ84cqFQqhIeHG1eXa9iwIWbPno3ExETodDpERkZi3bp1Fh3xq/D3339j7dq1AIC777670mNbtmxBv3796lU/k8V6EkKgoKDA+LuzEC4qFHi1AnQWWspLpzXUB/t9HeXQ1+bGKIe2mMuR2mIuVxcXhHgDgLDbNlu6X4QQ8PMwJIkuTjKy6IzvbUto27Yt1tcwx+fo0aMxevToKh9r2bJlla/3uHHjMG7cuEr3xcbGIjY2tsp6HnzwQTz44INmxQwYRjyTkpLMLn8zUzFbCk8Oo7qx5jyLRFQtN1fn3nU7S7JIZC+ce49Ddaew8KC03vGXaCOyFJXyn2Sp4jw+Z6IEk0UiW2KySHVjzUm5iahaNyeLWie7yAXgBS5EtsZkkerG0oehdTrL1kfkwFQ3TS9VbuFVsuSAy/0R2RaTRaobSy/JJziySGQu15uW+yvXOeFhaI4s2hwvsHFM5vYrr4a2gPrOXyRLLq5wKy8EVGrL1KcvN9Rn6RFLC5NDX5sboxzaYi5Haos5lC4KXC0HBIBSrR56OzxvUa8XFu+Xa+WAHoAzLY0t9Xu7YqWR4uJi49KD5DiKi4sBoMYVZRTCCf9dKCwshEajQUFBgcWWwnE62+YAW2YCUeOAIe/Xv76Tm4EvHwKaRwLP/Vr/+ogcXPibP6Ok3L4PQXu6KfHxY1Ho27qZReqLTNqAohItNr/SF2HNGlikTmdTl++/3NxcXLlyBb6+vvD09OQ61Q5ACIHi4mLk5+ejYcOG8Pf3r7Y8RxapbipGAC11YUrFsoF2PrJIZC96hDXBlszzUodRreIyHX49cd5iyaKeK7hIws/PDwCQn58vcSRkaQ0bNjT2b3WYLFLdVJyzaOm1obncH5FZlo3rgsvF9jvl1Edb/sTSX7Msek6l7saBMF7gYlsKhQL+/v7w9fVFebn9vueodlQqldmrAjFZrCedToc9e/YAALp16+Y0yzHpFCrsafMG4NoM3Syx3J+23FCfu8Yi9VmDHPra3Bjl0BZzOVJbzCWHNmvclXgpUsAXZ6HThVtkub9nw3UQABROtNyfPfWzUqmUPAaSBpPFehJC4NKlS8bfnYVQqnDJu53hd0ss5aXXWrQ+a5BDX5sboxzaYi5Haou55NBmlasCYRoAKLXYcn93aAy/O8sKLnLoZ3IOTBapbhQ3XY749z6gvjvvCycARNSvDiKyG25WvGSZh6GJbIvJItWN8qbpHL6IB/Sl9avPxR3otKx+dRCR3XC1cLJ48/RAzjR1DpE9YLJIdRPWF8j9w/B7o5aAqOdJzwpe2ELkSFQWzuh0Nx2GVXJkkcimmCxS3Xj7A7iRLD63E3Ct51tJqwXWr693WERkH1QWHv27eWTRhSOLRDbFjxwREVmcm6uFD0NzZJFIMhxZtABnnUrA0u2Ww+voSDHKoS3mcqS2mMve2+yqdEGJrv7XvlXQCYFSneF3Z7rAxd77mZwDl/vjcn9ERBaXmnEOE77Yh7uDGuL7idH1rq+guBwd3toIADgxM87i50Q6C37/UV3w00ZERBanUhpG/8p1lplAmxe4EEmHySIREVlcxTyLFksWK13gwmSRyJZ4zmI96XQ67N+/HwAQFRXlNOeXWLrdcngdHSlGObTFXI7UFnPJoc2uLsCz7QQ83Yqhs8ASnlqdFs+2E1AoYJH65EAO/UzOgcliPQkhkJ+fb/zdWVi63XJ4HR0pRjm0xVyO1BZzyaHNri4KRDQGAJ1FYtTqxI367LfNliaHfibnwMPQRERkcRXnLFrKzYehici2OLJIREQWd/Pa0InfHISunmMTJWXleKhZfaMiorpgskhERBbXyOuf9ePXH8lDmb5+I41uLoLJIpFEmCwSEZHFNb4pWfz34HaAop5nPQk9cCmjnlERUV0wWSQiIqsa3TUYrvVcP16r1WL9eiaLRFJwymSx4qqywsLCetel1WpRXFxsrK++O0S5sHS75fA6OlKMcmiLuRypLeaSQ5udcR9hadZoc8X3Hq+uptpw/E9bFYqKigAAQUFBEkdCRERke0VFRdBoNFKHQTLhlGtD6/V6nD17Ft7e3lBYYNmowsJCBAUFIScnR/ZrbbIt9oltsU9si31iW0wTQqCoqAgBAQFwceHseWQepxxZdHFxQWBgoMXr9fHxkf2OqQLbYp/YFvvEttgntqVqHFGk2uK/FURERERkEpNFIiIiIjKJyaIFuLu7Y9q0aXB3d5c6lHpjW+wT22Kf2Bb7xLYQWZZTXuBCRERERObhyCIRERERmcRkkYiIiIhMYrJIRERERCYxWSQiIiIik5gsEhEREZFJTBbNtGjRIoSGhkKtViMqKgo7duyotvy2bdsQFRUFtVqNsLAwfPLJJzaK1LRZs2ahS5cu8Pb2hq+vLx588EFkZmZWu83WrVuhUChuux07dsxGUVctKSnptpj8/Pyq3cYe+wQAWrZsWeVrPHHixCrL21OfbN++HUOGDEFAQAAUCgW+//77So8LIZCUlISAgAB4eHigX79+OHLkSI31rl69Gu3atYO7uzvatWuHNWvWWKkF/6iuLeXl5Zg8eTIiIyPh5eWFgIAAPP744zh79my1dSYnJ1fZVyUlJZK1BQDGjRt3W0zdu3evsV576xcAVb6+CoUCc+fONVmnVP1izj5YTp8Zch5MFs3wzTffICEhAW+88QbS0tLQu3dvxMXFITs7u8ryWVlZGDRoEHr37o20tDT861//wqRJk7B69WobR17Ztm3bMHHiROzevRupqanQarWIiYnBtWvXatw2MzMTubm5xtudd95pg4ird9ddd1WKKT093WRZe+0TANi7d2+ldqSmpgIAhg8fXu129tAn165dQ4cOHfDhhx9W+ficOXPw3nvv4cMPP8TevXvh5+eH++67D0VFRSbr3LVrF0aOHIkxY8bgjz/+wJgxYzBixAjs2bPHWs0AUH1biouLceDAAbz55ps4cOAAUlJScPz4ccTHx9dYr4+PT6V+ys3NhVqttkYTjGrqFwAYOHBgpZh++umnauu0x34BcNtru2zZMigUCjz88MPV1itFv5izD5bTZ4aciKAade3aVTz77LOV7gsPDxdTpkypsvxrr70mwsPDK933zDPPiO7du1stxrrIz88XAMS2bdtMltmyZYsAIC5fvmy7wMwwbdo00aFDB7PLy6VPhBDipZdeEq1atRJ6vb7Kx+21TwCINWvWGP/W6/XCz89PzJ4923hfSUmJ0Gg04pNPPjFZz4gRI8TAgQMr3RcbGytGjRpl8ZhNubUtVfn9998FAPHXX3+ZLLN8+XKh0WgsG1wtVdWWsWPHigceeKBW9cilXx544AFx7733VlvGHvpFiNv3wXL+zJBj48hiDcrKyrB//37ExMRUuj8mJgY7d+6scptdu3bdVj42Nhb79u1DeXm51WKtrYKCAgBA48aNayzbsWNH+Pv7o3///tiyZYu1QzPLiRMnEBAQgNDQUIwaNQqnTp0yWVYufVJWVoYVK1bgySefhEKhqLasPfbJzbKyspCXl1fpdXd3d0ffvn1NfnYA031V3TZSKCgogEKhQMOGDastd/XqVYSEhCAwMBCDBw9GWlqabQKswdatW+Hr64vWrVtjwoQJyM/Pr7a8HPrl3Llz+PHHHzF+/Pgay9pDv9y6D3b0zwzJF5PFGly4cAE6nQ7NmzevdH/z5s2Rl5dX5TZ5eXlVltdqtbhw4YLVYq0NIQQSExPRq1cvREREmCzn7++PTz/9FKtXr0ZKSgratGmD/v37Y/v27TaM9nbdunXDF198gQ0bNmDJkiXIy8tDz549cfHixSrLy6FPAOD777/HlStXMG7cOJNl7LVPblXx+ajNZ6diu9puY2slJSWYMmUKRo8eDR8fH5PlwsPDkZycjLVr1+Krr76CWq1GdHQ0Tpw4YcNobxcXF4eVK1di8+bNmDdvHvbu3Yt7770XpaWlJreRQ798/vnn8Pb2xtChQ6stZw/9UtU+2JE/MyRvrlIHIBe3jvIIIaod+amqfFX3S+WFF17AoUOH8Ouvv1Zbrk2bNmjTpo3x7x49eiAnJwf/+c9/0KdPH2uHaVJcXJzx98jISPTo0QOtWrXC559/jsTExCq3sfc+AYClS5ciLi4OAQEBJsvYa5+YUtvPTl23sZXy8nKMGjUKer0eixYtqrZs9+7dK104Eh0djU6dOmHhwoX44IMPrB2qSSNHjjT+HhERgc6dOyMkJAQ//vhjtYmWPfcLACxbtgyPPvpojece2kO/VLcPdrTPDMkfRxZr0LRpUyiVytv+Q8vPz7/tP7kKfn5+VZZ3dXVFkyZNrBaruV588UWsXbsWW7ZsQWBgYK237969u+QjI7fy8vJCZGSkybjsvU8A4K+//sIvv/yCp556qtbb2mOfVFydXpvPTsV2td3GVsrLyzFixAhkZWUhNTW12lHFqri4uKBLly5211f+/v4ICQmpNi577hcA2LFjBzIzM+v0+bF1v5jaBzviZ4YcA5PFGri5uSEqKsp4hWqF1NRU9OzZs8ptevTocVv5jRs3onPnzlCpVFaLtSZCCLzwwgtISUnB5s2bERoaWqd60tLS4O/vb+Ho6qe0tBRHjx41GZe99snNli9fDl9fX9x///213tYe+yQ0NBR+fn6VXveysjJs27bN5GcHMN1X1W1jCxWJ4okTJ/DLL7/U6Z8MIQQOHjxod3118eJF5OTkVBuXvfZLhaVLlyIqKgodOnSo9ba26pea9sGO9pkhByLFVTVy8/XXXwuVSiWWLl0qMjIyREJCgvDy8hKnT58WQggxZcoUMWbMGGP5U6dOCU9PT/Hyyy+LjIwMsXTpUqFSqcSqVaukaoIQQojnnntOaDQasXXrVpGbm2u8FRcXG8vc2pb58+eLNWvWiOPHj4vDhw+LKVOmCABi9erVUjTB6JVXXhFbt24Vp06dErt37xaDBw8W3t7esuuTCjqdTgQHB4vJkyff9pg990lRUZFIS0sTaWlpAoB47733RFpamvEK4dmzZwuNRiNSUlJEenq6eOSRR4S/v78oLCw01jFmzJhKMwv89ttvQqlUitmzZ4ujR4+K2bNnC1dXV7F7927J2lJeXi7i4+NFYGCgOHjwYKXPT2lpqcm2JCUlifXr14uTJ0+KtLQ08cQTTwhXV1exZ88eydpSVFQkXnnlFbFz506RlZUltmzZInr06CFatGghu36pUFBQIDw9PcXHH39cZR320i/m7IPl9Jkh58Fk0UwfffSRCAkJEW5ubqJTp06VppsZO3as6Nu3b6XyW7duFR07dhRubm6iZcuWJnditgSgytvy5cuNZW5ty7vvvitatWol1Gq1aNSokejVq5f48ccfbR/8LUaOHCn8/f2FSqUSAQEBYujQoeLIkSPGx+XSJxU2bNggAIjMzMzbHrPnPqmYxufW29ixY4UQhqlApk2bJvz8/IS7u7vo06ePSE9Pr1RH3759jeUrfPfdd6JNmzZCpVKJ8PBwmyTC1bUlKyvL5Odny5YtJtuSkJAggoODhZubm2jWrJmIiYkRO3fulLQtxcXFIiYmRjRr1kyoVCoRHBwsxo4dK7KzsyvVIYd+qbB48WLh4eEhrly5UmUd9tIv5uyD5fSZIeehEOLGWf5ERERERLfgOYtEREREZBKTRSIiIiIyickiEREREZnEZJGIiIiI/r+9+wdpq4vDOP7cRhNiDalJiilSaBUHRU0UrFQCzVLS1aFbijpkcBAiuEgkggoBh9JMooilo+AgFJGOUhHaQium2jU4GCzoIKSUlKbv8iKIOaF/bIP6/UzJ4Z7f+Z07PZzkco0IiwAAADAiLAIAAMCIsAgAAAAjwiIAAACMCIsAAAAwIiwC+GXhcFjxeLzSbQAA/gFe9wegrHA4rGAwqGfPnp2MHR0dqbq6Wi6X65/3E4/Hlc1mtbKy8s/XBoCriJNFAL/M4/FUJChK0rt373Tv3r2KrA0AVxFhEYDRwMCA1tfXlU6nZVmWLMtSNps98zN0OBzW8PCw4vG46urqVF9fr/n5eeXzeQ0ODsrlcqmpqUlra2snc378+KGZmRk1NjbK6XQqEAhoeXnZ2Mu3b99kt9u1ubmpRCIhy7LU09PzN7cPABBhEUAZ6XRa9+/fVywWUy6XUy6X0+3bt0te++LFC/l8Pr19+1bDw8MaGhrS48eP1dvbq/fv3ysSiejJkyf68uWLJGl8fFzPnz/X7OysdnZ2NDIyomg0qvX19ZL1bTabNjY2JElbW1vK5XJ69erV39k4AOAEYRGAkdvtlt1uV01Njfx+v/x+v2w2W8lrA4GAxsfH1dzcrLGxMTmdTvl8PsViMTU3NyuZTOrw8FDb29vK5/N6+vSpFhcXFYlE1NjYqIGBAUWjUc3NzZWsf+3aNe3v78vr9SoQCMjv9+vGjRuSpI2NDaVSqb91GwDgSquqdAMALoeOjo6TzzabTV6vV+3t7Sdj9fX1kqTPnz9rd3dXX79+1cOHD0/VKBQK6uzsNK7x4cMHBQKBM+OhUEihUOhPtwAAKIGwCOBcVFdXn/puWdapMcuyJEnFYlHFYlGStLq6qoaGhlPzHA6HcY2tra2SYbGvr08TExMKBoO/2z4AwICwCKAsu92u79+/n2vN1tZWORwO7e3t6cGDBz89L5PJqK+v78z4p0+f1NLScp4tAgD+R1gEUNadO3f05s0bZbNZ1dbWyuPx/HFNl8ul0dFRjYyMqFgsKhQK6fj4WJubm6qtrVV/f3/JecViUdvb29rf39f169fldruVz+dVVVVV9kQSAPD7eMAFQFmjo6Oy2WxqbW3VzZs3tbe3dy51p6amlEwmlUql1NLSokgkopcvX+ru3bvGOdPT01paWlJDQ4MmJyclSR8/flRbW9u59AQAOIs3uAC40BYWFnRwcKBEIlHpVgDgUuJkEcCFlslkTj11DQA4X4RFABfa69ev1d3dXek2AODSIiwCuJAKhYK6urr06NEj3bp1q9LtAMClxX8WAQAAYMTJIgAAAIwIiwAAADAiLAIAAMCIsAgAAAAjwiIAAACMCIsAAAAwIiwCAADAiLAIAAAAI8IiAAAAjAiLAAAAMCIsAgAAwOg/WpwksHfbCEEAAAAASUVORK5CYII=",
"text/plain": [
"<Figure size 640x480 with 3 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from numpy import array,linspace\n",
"from matplotlib import pyplot as plt\n",
"\n",
"steps_times=array([\n",
" [0,0,0],\n",
" [3,0,-1],\n",
" [4,0,1],\n",
" [6,3/5,0],\n",
" [8,0,2],\n",
" [8.5,0,-2],\n",
" [11,-3/5,0],\n",
" [14,-1.5/2,0],\n",
" [16,-1.5/2,0.5],\n",
" [17,1.5,-0.5],\n",
" [20,0,0]\n",
"])\n",
"times=steps_times[:,0]\n",
"steps=steps_times[:,[1,2]]\n",
"\n",
"lt=array([[times[i]-times[j] if times[i]-times[j]>=0 else 0 for j in range(times.shape[0])] for i in range(times.shape[0])])\n",
"synth_ramps=array([2,1])+lt.dot(steps)\n",
"\n",
"times_extended=linspace(times.min(),times.max(),1000)\n",
"sigma=array([[1 if times_extended[i]-times[j]>=0 else 0 for j in range(times.shape[0])] for i in range(times_extended.shape[0])])\n",
"slopes=sigma.dot(steps)\n",
"\n",
"print('lower triangular (t-tj)sigma(t-tj)) \\n',lt)\n",
"fig,ax_list=plt.subplots(nrows=3,ncols=1,height_ratios=[1,1/2,1/2],constrained_layout=True)\n",
"ax_list[0].plot(times,synth_ramps,'-o')\n",
"for j in range(steps.shape[1]):\n",
" ax_list[1].stem(times,steps[:,j],markerfmt=['blue','orange'][j],label=f'series i={j+1}')\n",
" ax_list[2].plot(times_extended,slopes[:,j],label=f'series i={j+1}')\n",
"for j in range(len(times)):\n",
" ax_list[2].axvline(times[j],linestyle='--',color=[0.7,0.7,0.7])\n",
"ax_list[-1].set_xlabel('time $t_j$')\n",
"ax_list[0].set_ylabel('$y_i$')\n",
"ax_list[1].set_ylabel('$a_{i,j}$')\n",
"ax_list[2].set_ylabel(r'$\\frac{dy_i}{dt}$')\n",
"ax_list[0].legend([f'series i={j+1}' for j in range(steps.shape[1])],loc='center left',bbox_to_anchor=[1.01,0.5]);\n",
"ax_list[1].legend(loc='center left',bbox_to_anchor=[1.01,0.5])\n",
"ax_list[2].legend(loc='center left',bbox_to_anchor=[1.01,0.5]);"
]
}
],
"metadata": {
"jupytext": {
"cell_metadata_filter": "-all",
"formats": "ipynb,py:percent"
},
"kernelspec": {
"display_name": "main",
"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.11.4"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Raw
synthetic_ramp_method.py
# ---
# jupyter:
# jupytext:
# cell_metadata_filter: -all
# formats: ipynb,py:percent
# text_representation:
# extension: .py
# format_name: percent
# format_version: '1.3'
# jupytext_version: 1.16.7
# kernelspec:
# display_name: main
# language: python
# name: python3
# ---
# %% [markdown]
# synthetic ramp method
#
# calculate synthetic ramps using convolution for a given series of time points:
# * define matrix $A$ with the coefficients $a_{i,j}$ of the unit steps in $n$ different series $i$ for each time $t_j$.
# * define matrix $\sigma(t-t_j)*\sigma(t)=(t-t_j)\cdot \sigma(t-t_j)$ which is lower triangular with permutations of $(t-t_j)\cdot \sigma(t-t_j)$
# * example with $t_0<t_1$: $(t_0-t_1)\cdot \sigma(t_0-t_1)=0$, and $(t_1-t_0)\cdot \sigma(t_1-t_0)=(t_1-t_0)$ due to step definition
# * the matrix $\sigma(t-t_j)\cdot A$ has then the slopes $\frac{dy_i}{dt}$ at times $t_j$ and $(t-t_j)\cdot \sigma(t-t_j)\cdot A$ has the synthetic ramps by convolution.
# $$
# \underbrace{\left[\begin{array}{cccc}
# 1 & 0 & 0 & ...\\
# 1 & 1 & 0 & ...\\
# 1 & 1 & 1 & ...\\
# ... & ... & ... & ...
# \end{array}\right]}_{\sigma(t-t_j)}\cdot
# \underbrace{\left[\begin{array}{ccc}
# a_{1,0} & a_{2,0} & ...\\
# a_{1,1} & a_{2,1} & ...\\
# a_{1,2} & a_{2,2} & ...\\
# ... & ... & ...\\
# \end{array}\right]}_{A} =
# \underbrace{\left[\begin{array}{ccc}
# a_{1,0} & a_{2,0} & ...\\
# a_{1,0}+a_{1,1} & a_{2,0}+a_{2,1} & ...\\
# a_{1,0}+a_{1,1}+a_{1,2} & a_{2,0}+a_{2,1}a_{2,2} & ...\\
# ... & ... & ...\\
# \end{array}\right]}_{\mathrm{slopes}}
# $$
#
# $$
# \underbrace{\left[\begin{array}{cccc}
# t_0-t_0 & 0 & 0 & ...\\
# t_1-t_0 & t_1-t_1 & 0 & ...\\
# t_2-t_0 & t_2-t_1 & t_2-t_2 & ...\\
# ... & ... & ... & ...
# \end{array}\right]}_{(t-t_j)\cdot \sigma(t-t_j)}\cdot
# \underbrace{\left[\begin{array}{ccc}
# a_{1,0} & a_{2,0} & ...\\
# a_{1,1} & a_{2,1} & ...\\
# a_{1,2} & a_{2,2} & ...\\
# ... & ... & ...\\
# \end{array}\right]}_{A} =
# \underbrace{\left[\begin{array}{ccc}
# a_{1,0}\cdot(t_0-t_0) & a_{2,0}\cdot(t_0-t_0) & ...\\
# a_{1,0}\cdot(t_1-t_0)+a_{1,1}\cdot(t_1-t_1) & a_{2,0}\cdot(t_1-t_0)+a_{2,1}\cdot(t_1-t_1) & ...\\
# a_{1,0}\cdot(t_2-t_0)+a_{1,1}\cdot(t_2-t_1)+a_{1,2}\cdot(t_2-t_2) & a_{2,0}\cdot(t_2-t_0)+a_{2,1}\cdot(t_2-t_1)+a_{2,2}\cdot(t_2-t_2) & ...\\
# ... & ... & ...\\
# \end{array}\right]}_{\mathrm{synth\ curves}}
# $$
#
# example below with two series $i=\{1,2\}$
# %%
from numpy import array,linspace
from matplotlib import pyplot as plt
steps_times=array([
[0,0,0],
[3,0,-1],
[4,0,1],
[6,3/5,0],
[8,0,2],
[8.5,0,-2],
[11,-3/5,0],
[14,-1.5/2,0],
[16,-1.5/2,0.5],
[17,1.5,-0.5],
[20,0,0]
])
times=steps_times[:,0]
steps=steps_times[:,[1,2]]
lt=array([[times[i]-times[j] if times[i]-times[j]>=0 else 0 for j in range(times.shape[0])] for i in range(times.shape[0])])
synth_ramps=array([2,1])+lt.dot(steps)
times_extended=linspace(times.min(),times.max(),1000)
sigma=array([[1 if times_extended[i]-times[j]>=0 else 0 for j in range(times.shape[0])] for i in range(times_extended.shape[0])])
slopes=sigma.dot(steps)
print('lower triangular (t-tj)sigma(t-tj)) \n',lt)
fig,ax_list=plt.subplots(nrows=3,ncols=1,height_ratios=[1,1/2,1/2],constrained_layout=True)
ax_list[0].plot(times,synth_ramps,'-o')
for j in range(steps.shape[1]):
ax_list[1].stem(times,steps[:,j],markerfmt=['blue','orange'][j],label=f'series i={j+1}')
ax_list[2].plot(times_extended,slopes[:,j],label=f'series i={j+1}')
for j in range(len(times)):
ax_list[2].axvline(times[j],linestyle='--',color=[0.7,0.7,0.7])
ax_list[-1].set_xlabel('time $t_j$')
ax_list[0].set_ylabel('$y_i$')
ax_list[1].set_ylabel('$a_{i,j}$')
ax_list[2].set_ylabel(r'$\frac{dy_i}{dt}$')
ax_list[0].legend([f'series i={j+1}' for j in range(steps.shape[1])],loc='center left',bbox_to_anchor=[1.01,0.5]);
ax_list[1].legend(loc='center left',bbox_to_anchor=[1.01,0.5])
ax_list[2].legend(loc='center left',bbox_to_anchor=[1.01,0.5]);
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment