Skip to content

Instantly share code, notes, and snippets.

@ruoyu0088
Last active December 14, 2023 23:24
Show Gist options
  • Star 19 You must be signed in to star a gist
  • Fork 5 You must be signed in to fork a gist
  • Save ruoyu0088/70effade57483355bbd18b31dc370f2a to your computer and use it in GitHub Desktop.
Save ruoyu0088/70effade57483355bbd18b31dc370f2a to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"import pylab as pl\n",
"\n",
"def make_test_data(seg_count, point_count):\n",
" x = np.random.uniform(2, 10, seg_count)\n",
" x = np.cumsum(x)\n",
" x *= 10 / x.max()\n",
" y = np.cumsum(np.random.uniform(-1, 1, seg_count))\n",
" X = np.random.uniform(0, 10, point_count)\n",
" Y = np.interp(X, x, y) + np.random.normal(0, 0.05, point_count)\n",
" return X, Y"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"from scipy import optimize\n",
"\n",
"def segments_fit(X, Y, count):\n",
" xmin = X.min()\n",
" xmax = X.max()\n",
"\n",
" seg = np.full(count - 1, (xmax - xmin) / count)\n",
"\n",
" px_init = np.r_[np.r_[xmin, seg].cumsum(), xmax]\n",
" py_init = np.array([Y[np.abs(X - x) < (xmax - xmin) * 0.01].mean() for x in px_init])\n",
"\n",
" def func(p):\n",
" seg = p[:count - 1]\n",
" py = p[count - 1:]\n",
" px = np.r_[np.r_[xmin, seg].cumsum(), xmax]\n",
" return px, py\n",
"\n",
" def err(p):\n",
" px, py = func(p)\n",
" Y2 = np.interp(X, px, py)\n",
" return np.mean((Y - Y2)**2)\n",
"\n",
" r = optimize.minimize(err, x0=np.r_[seg, py_init], method='Nelder-Mead')\n",
" return func(r.x)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4lFX2wPHvmUkhgQQiCS2hhSbFAgQIFhQVFVZBRUF0\nxY7u6lp2FVFXUGzouqvuDxXBioINFZAVFBXBQiIZLDRpkYQAQghjEggpk7m/P6Ywk0zIAAlJJufz\nPHnMvO9937mzC2cu9z33XDHGoJRSKrRY6roDSimlap4Gd6WUCkEa3JVSKgRpcFdKqRCkwV0ppUKQ\nBnellApBGtyVUioEaXBXSqkQpMFdKaVCUFhdvXF8fLzp1KlTXb29Uko1SDabba8xJqG6dnUW3Dt1\n6kRGRkZdvb1SSjVIIpIVTDudllFKqRCkwV0ppUKQBnellApBGtyVUioEaXBXSqkQpMFdKaVCUNDB\nXUSsIvKjiCwKcC5SRN4TkS0iki4inWqyk0opVR/Ysuy8sGwLtix7XXelWkeS534nsAGIDXDuRsBu\njOkqIlcCTwFja6B/SilVL9iy7IybuZLScoNV4OYzk4mJCic1uSX9O8bVdfcqCSq4i0gS8CfgceDv\nAZqMAh52/z4PmC4iYnSDVqVUA2TLspOWmUdqcksA8ma+zkkvPMWv+bnsjI3n6SHjmeGObmEWYeqo\nPlw1qEMd9riyYEfuzwETgZgqzicC2wGMMQ4RyQdaAnuPuYdKKXUc2bLsjJuVRqnDiUXg4rXLeHLJ\ndKIdJQAkFeQybcl0ABb2HorDaZi8YC092sTUqxF8tXPuInIRsMcYYzvWNxORCSKSISIZubm5x3o7\npZSqcR+tzqHU4QTAaeDeFbO9gd0j2lHCxBWzva/LnYa0zLzj2s/qBPNA9XRgpIhsA94FzhGRtyu0\n2QG0BxCRMKA5UOmTGmNmGmNSjDEpCQnV1r1RSqnjypZlZ+2OfL9j7QoCT0D4HjfA5t2Ftdm1I1bt\ntIwx5n7gfgARORu4xxjz5wrNFgLXAiuBy4GvdL5dKdVQ2LLsfLg6h3m2Q6N2gJ57MnGKYAkQzvZF\n++eWzP9pJ/sOlDL7xkG13t9gHHWeu4hMFZGR7pevAi1FZAuuB66TaqJzSilV26Z9uoErZnzP3PRs\nv8B+6dqv+Oite9kfEUWxNdzvGidCfFE+k5a9Rli5w3t8xea93PXuj/UiXVLqaoCdkpJitOSvUqqu\n2LLsvLx8K5+v3+13PLy8jIe+fIXxP/6PlR1O4m8jJ3Latp+YuGI27Qr2sjM2nmdPv4pTf9/MNT9+\nyo9te/C3kfeS06KN330sAmd0jccAw/u0rbFsGhGxGWNSqm2nwV0p1djYsuxc/UoaxWVOv+NtCvby\n4oIn6bdzIzMGXsa/zrqWcou1yvsM//VbnlryfwBMHH4HS3qcXmXbJy49qUYCfLDBXcsPKKUaFVuW\nnee+2ERJhcA+OOsXFr15J933ZnPrJfczbegNhw3sAItPPIMR1z1P5gmJzJj/JI9+/iKRjtKAbd9b\nlV1jnyEYGtyVUo2GZ8T+zea9eOcsjGFC+oe8/d4/sUfFMmr8f6ocgVuAAZ38c9lzWrThiquf4uWB\nl3HNj5/y8Vv/IDkvp9K1v+TkMzf9+AV4De5KqUZj2uINflMxzUqKeGn+kzzw9ess7n4al1zzb7a2\nbF/l9RaLMGl4T87v1drveJk1nCeH3sB1l0+hTWEen7x5F5et/dKvjQEe+HjNcQvwOueulAppngen\nP2bbyd1/aMqk695sXv74CTrad/Lk2dfz6oBLQKTS9VaLYIzBIofKDNiy7Ix5eSXlzsrxs3XhXp7/\n5BlSt6/lwz7n8NCwv1AUEeU9bxG4cmAHRvdLOqoVrfpAVSnV6LmC8PeU+0+v86cN3/D04ucpimjC\nbaMm8UP7PgGvj7AKD4/sg72otFKBsLnp2Ty0YG3AAG9xlnPH9+9yx3fv8tsJidw+aiIbWiX7tWkS\nbmHOTalHHOCDDe5HUhVSKaUaBE/hr5+3/+EX2MPKHUz6+nVuylhARmJP/jpqEntiWga8hwAPj6y6\nINhVgzrQo00MH63O4euNe9jxR7H3nNNi5bkzria9/Uk8t+gZ5s/+B4+eezMFEVF+KZXrnZPg8UC1\nGI+djtyVUiHF89C01OHEd1CdsN/O9AXTGJSzjtf7X8wTQ2+grMLiJF8C3HNBD24b2jWo9/SUA66o\n5YE/+Pf/nuXs32w4xEKYOfRtUx4VhXXWLLj66qA/n6ZCKqUapbTMPErK/AN7/5z1LHrzTk7avYU7\nLr6HR867JWBgt/hMuYdbxVvytzr9O8bxzoTBXDWoAwM7xRHfLMJ7Lq9pC66/Ygr5kU39AjuA9eBB\nePDBI/uAQdJpGaVUSJibns17q7LJP1jml+Z4ne0THlz2KjnNWzF+zFQ2JnSq8h4Wi3Bej1bEx0Qe\n8QPP/h3j6N8xzvsvB19GLMSUFAW+MLt2smc0uCulGixblp2PVuewaXchq7b513KJKi1m2pL/Y9SG\n5XzeLZV//OluCiObVrpHm9hIfi9wlfQ1TsMp7VsENRVTlbTMPL8aNR47Y+NJKghQ6rxD7WzyocFd\nKdVgVNwhybOpRkWd9+3gpY+foFvedp4eMp6XUi/HSOBZ6EtOTeSNldsoczgJD7MEPRVTldTklkSE\nWShz9yshJpLurWN4esh4pvls+gFAdDQ8/vgxvV9VNLgrpRoE3welEWEWhnRLCBjYh21O49+L/kOZ\nNYzxY6byXadTq7ynBYiJCmfOTaneL41j3U2pf8e4Sve75tV0FvYeCuDNlrF07OAK7EfwMPVIaHBX\nSjUInukOp4FSh7NSNUeLs5x/fPM2t6V9wE9tu/HXS+5nZ2yrw94zzP3Q1DNfXlMq3m94n7Z8s3kv\nC3sPZWHvoZzfqzW3nNWlVrfl0+CulGoQPNMdFTNhAE4oyuf5hf/izKyfmHvKhTxy3gRKwiL82giu\nbBhPtqIAV6S0Py77nnpy5d9blc26XQV8sWE3KzbnHtUipmBpcFdKNQgbfy+kZdMIv8VCAKfs3MiL\n86cRX/QH9w6/kw9OHuZ33iIw4cxkYqLCiYuOYOqidd759cv6JR23/l81qAP2olLW7MjHaaDM4SQt\nM6/ugruINAFWAJHu9vOMMVMqtLkO+BeuvVQBphtjXqnZriqlGiNblp0Zy7eytMI0DMYw7ufPePiL\nGexp1pLL/vwv1rXxz3JJbNGE/47r5xdAe7SJqbH59SPl+7C1Jh7eHk4wI/cS4BxjzH4RCQe+FZHF\nxpi0Cu3eM8bcXvNdVEo1VnPTs5m8YC2OCvMwkWUlPLr0Jcas+YKvO/fnzovvIT8qxq+NVagU2KHy\nfPjxFOhha20JZoNsA+x3vwx3/+jm10qpWmXLsvPQ/DVUXNGf9MfvzJj/JH12b+W508fx/OnjKqU5\nWgUeveSkOgvih3O8vlyCmnMXEStgA7oCLxhj0gM0Gy0iQ4BNwN3GmO01102lVGNiy7Iz9ZN1lQL7\n2VszeG7RM4gxXH/5FJZ1GVDlPexFgXdEaiyCqi1jjCk3xpwKJAEDRaRifcxPgE7GmJOBpcCbge4j\nIhNEJENEMnJzA6zUUko1ep4iXD/n5HuPiXFy57dzeW3eI+yMTeCi656vFNgFV4leq1Dr89kNwRFl\nyxhj/hCRZcCFwFqf43k+zV4Bnq7i+pnATHBVhTzi3iqlQt6M5Vv9qis2P1jIs4v+zTmZGczrcy7/\nPP8vFIc3qXTd45eeVKcPS+ubYLJlEoAyd2CPAoYBT1Vo09YYs8v9ciSwocZ7qpQKWZ4aMauz7Gz4\nvdB7vPfurbz08RO0KczjgQtuY+4pFwbcLckirmmYunxYWt8EM3JvC7zpnne3AO8bYxaJyFQgwxiz\nELhDREYCDmAfcF1tdVgpFVpsWfaANWJGr/mSxz9/gX1RsYy5+il+atejyntYJPjyvI1FMNkyvwB9\nAxyf7PP7/cD9Nds1pVQo890tyTewRzjKmPzlTP7802K+63gyd1w8kbymLaq8T5jFtbepjtj96QpV\npdRx5ykCVlLm9MurbluQy0vzn+TUXZt4MfVy/n3mNZRbrFXe55Sk5ky+uLcG9gA0uCulalXFMr1p\nmXns/ONgpcB+2raf+L+FTxNRXsYtlz7AZ91Pq/befRKba2CvggZ3pVSt8S3TG2a1gDE4nIYw3/3s\njOHW9A+5d8VstrRM4tZLH+S3ExIr3evWIckUljh4L2M75eXmuNeGaWg0uCulak3FMr0enlTHmJID\n/OvT57hw00o+OfFM7ht+B0URUZXuc+uQZCaN6AnAZf2SNN0xCBrclVK1JjW5JVaL4Ky41BTonruN\nGR8/Qfv83Txy7s283n9kwDRHcG2o4aHpjsHR4K6UqlWOAIF95PrlTFvyX/ZHRHPVlY+zqn3FRe+H\nWAVNczwKGtyVUrUmLTPP76FpWLmDB5a9xg22hfyQ1IvbRk0it9kJVV5vtQiPaprjUdHgrpSqcZ4M\nmcKDZd5jCfv38cKCaQzMWc8rKaOYdvb1OKyBQ5AAZ3SL567zumtgP0oa3JVSNWpuejYPzV/jtxXe\ngO1reWHBUzQrLeJvF9/LJ73OOuw9rBbRwH6MNLgrpY6ZZ6QeFx3Bgx+vOTQVYww3ZCzkgWWvkhXX\nlqvHPsbmhI4B72ERMMYV2HXF6bHT4K6UOiaeXPbiMv/aMNGlB3lq8X+5+NdvWNJ9MPeMuJv9kdHe\n8xFWweE0CK4SvZMv6o29qFRTHGuIBnel1FHzbKpRMbAn5+Uw4+Mn6LIvhyfPvo6XB46ulObocBrC\nrBYu75/E6H5JGtBrmAZ3pdRR8WyqUVpuGLluGRNXzKZdwV7sUTFElx7kQGQ0fx77KCs7nhLweqeB\n8nIniS2iNLDXAg3uSqmjkpaZR5k7sE9bMp1oRwkALQ8W4ESYNnhslYEdXHPsumNS7dHgrpQ6KqnJ\nLbFYYOKK2d7A7mHBcPOqj3kzZaTf8aQWTRjSoxV92jXX+fVapsFdKXVU+neMo3PLprQr2BvwfMXj\nEWEWnh/XT4P5cVLtBtki0kREfhCRn0VknYg8EqBNpIi8JyJbRCRdRDrVRmeVUvXH3PRstuYeYGds\nfMDzvscFuLy/PjQ9nqoN7kAJcI4x5hTgVOBCEUmt0OZGwG6M6Qo8S4U9VpVSocWWZWfygrUYYEn3\nwZXOF4VF8vSQ8d7XVoswWsvzHlfVBnfjst/9Mtz9U7ES0CjgTffv84BzRaoo76aUatBsWXae+2IT\nDqchqrSY4RtXsiMmnh2xCTgRcmITmHTh7SzsPRTQbfDqSlBz7u7NsW1AV+AFY0x6hSaJwHYAY4xD\nRPKBlkDgyTilVINUccHSHd+/S2JhLqOvfhpbUi+/tmFWYUxKe81hryNBBXdjTDlwqoi0AD4WkT7G\nmLVH+mYiMgGYANChQ4cjvVwpVcc8m28AdNm7nZtWfcwHfc7zBvbze7UmISYSAxrU69gRZcsYY/4Q\nkWXAhYBvcN8BtAdyRCQMaA7kBbh+JjATICUlpXKRZ6VUvRYXHYExgDE8uvQlisKbMO3s67znz+7R\niqsG6cCtPggmWybBPWJHRKKAYcCvFZotBK51/3458JUxRoO3UiHEU+3RACM3rOC07F/411nXkte0\nBeDKiLEXldZpH9UhwYzc2wJvuufdLcD7xphFIjIVyDDGLAReBd4SkS3APuDKWuuxUuq4m5uezYPz\n12AMNCsp4sFlr/Jzm27MPeUCwBUYIsJ1tWl9Um1wN8b8AvQNcHyyz+/FwBU12zWlVH1gy7LzT3dg\nB7j72zkk7Ldz82X/xGmxIgJXDuygc+z1TDB57kqpRuzD1TnejTdO3PMb19o+4Z1TL+CXtt1dBw1a\n/Kse0uCulDosz4IVMU4e/fwl8ps04+kh13rPR+p0TL2ktWWUUod1oMQBwOi1XzFgx3ruHX4H+VEx\nAFw1SKdj6isN7kqpKs1Nz2b+TzuJLd7PpK9fx9buROaddB4AAzvF8cSlJ9VxD1VVdFpGKVWlxWt3\nAXDvitnEHSzkofP/ihELAtw3vGfddk4dlgZ3pVRAtiw7m38v5KRdm7n6x8W82e8i1rdOBiClU5xO\nxdRzOi2jlKrElmVn7MyVOMsczFj6InubtuDZM68GIMwCk3TUXu/pyF0pVUlaZh6OcsOVv3zOqbs2\n89g5N1IY2RSAqaNO0lF7A6DBXSlVSVx0BCcU5TNx+Zt83+FkFvY8y3tOSww0DDoto5QCXFMxaZl5\nxEVH8N6qbO77+g2alh5k8rBbwb09g0XQnPYGQoO7Uspbp72kzIkB+uVsYOyapcwYNJot8YeqPI48\npZ1OyTQQOi2jlPLWaTeA1VnOY0tfZGdMPP89zb8GYLfWMXXTQXXENLgrpUhNbonFXWfgmtX/o9ee\n35h67s0URUR52+iUTMOiwV0pBUC5ExL27+Pv37zN8s79WNL9NL/zj12iWTINiQZ3pRRPLd6AAR5c\n9iqR5aVMOe8W70NUgK4JTXWHpQZGg7tSjZwty876XQUMzvqFS9YvZ8agy9l2QqJfm84Jzeqod+po\nBbPNXnsRWSYi60VknYjcGaDN2SKSLyI/uX8mB7qXUqp+8WbJFBUzdelLZDdvzYup/vvuWARuPatL\nHfVQHa1gUiEdwD+MMatFJAawichSY8z6Cu2+McZcVPNdVErVlrTMPIrLnNySsYBuedu5YfRkSsIj\nvefDLMLUUX10rr0BqnbkbozZZYxZ7f69ENgAJB7+KqVUQxAXHUHbglzu/O4dlnYdxFddB/qdnzqq\nj861N1BHNOcuIp1w7aeaHuD0YBH5WUQWi0jvGuibUqqWvZ22jYe+nIUYeOS8CZXOa6mBhivo4C4i\nzYAPgbuMMQUVTq8GOhpjTgH+D5hfxT0miEiGiGTk5uYebZ+VUjXgrnd/JP67rxmx6XumDx5DTvPW\nfucjrKJ57Q1YUMFdRMJxBfY5xpiPKp43xhQYY/a7f/8UCBeR+ADtZhpjUowxKQkJCcfYdaXU0bJl\n2VmcsY1HvpjB1hMSmTXwMr/zCc0ieGfCYJ1rb8CCyZYR4FVggzHmP1W0aeNuh4gMdN83ryY7qpSq\nGbYsO899sYkJ6R/S2b6LKefdSmlYuF+bYb3baGBv4ILJljkduAZYIyI/uY89AHQAMMbMAC4H/iIi\nDuAgcKUxxtRCf5VSx8CT+hifu5NZaR+w6MQz+bZzX782EVZhdL+kOuqhqinVBndjzLeAVNNmOjC9\npjqllKp5nhF7SZmTh794GYfFyqPn3OjXZliv1tx6VhcdtYcALfmrVCPgW9L33M3pnLd1FY8NvYHd\nMYcejQ3pFs+s8Sl12EtVk7T8gFKNQFpmHiVlTiLLinn4i5fZGN+BN/qP9Guj86ihRUfuSjUCcdER\nGOC2lR+QVLCHMVdNw2H1/+s/vE/buumcqhUa3JVqBNbuzKfzvh1M+OFDPuw9lB/a9/GeE+CWIcm6\nEjXEaHBXKoTZsux8uDqHd9OzeGPpDErCInly6A1+bQxQWOKomw6qWqPBXakQ5XmIWlzmZMSv3zFk\n249MPu8W9jatnAmj8+2hRx+oKhWiPA9Rm5YU8dCXs1jbugtv9x1RqZ3mtYcmHbkrFYJsWXZ2/HEQ\ngDu+f5e2+/P46yX347RYAdcmS7ecmUxMVDipyS01rz0EaXBXKsT4Tsd0y83ihowFvHPy+fyYeKK3\nzbCerZk0omcd9lLVNp2WUSrEfLg6h+IyJxjDY0tfYn9ENE+fda1fm/iYyCquVqFCg7tSIcSWZWee\nLQeAS9ctY9D2tTx11rXYo5t721gFnWNvBDS4KxVC0jLzcJQ7iS3ezwPLXuPHtj1475Tz/do8eslJ\nOsfeCGhwVypEeB6iisDfv3mbEw4W8M/z/4KRQ3/NE1s00cVKjYQ+UFUqBPgWBuv1+xau+fFT3u47\nnHVtuvq1i4rQv/KNhY7clQoBaZl5lDqcYJw89vlL7IuK5d9nXlOp3Q2nd66D3qm6oF/jSoWAuOgI\nnAbG/rKUvrs2cvef/k5Bk2aAq3bMyUnNGTugg07JNCLVBncRaQ/MBlrjWqU80xjzfIU2AjwPjACK\ngOuMMatrvrtKqYqmfbqBGSsyaXGwgPuWv0l6+z583Huo9/x5ugFHoxTMtIwD+IcxpheQCtwmIr0q\ntBkOdHP/TABeqtFeKqUCmpuezYwVmQBMXP4mscX7eWjYra4lqG5L1+/m6lfSsGXZ66qbqg5UG9yN\nMbs8o3BjTCGwAUis0GwUMNu4pAEtRESLQytVy95blQ3AqTs3cuXPn/Nayig2JXSq1K60zElapu5Z\n35gc0QNVEekE9AXSK5xKBLb7vM6h8heAUqqGRYRZsDjLefTzF9nTLI7nTx8XsJ3FIqQmtzzOvVN1\nKegHqiLSDPgQuMsYU3A0byYiE3BN29Chgz7YUepo2bLsvLx8Kxnb7Pz5p8WctHsrt428jwOR0ZXa\nhlmEqaP66Jx7IxNUcBeRcFyBfY4x5qMATXYA7X1eJ7mP+THGzARmAqSkpGgJaaWOgi3LzrhZaZQ6\nnMQfsHPvirf4puOp/O/EM/zaCXBGt3juOq+7BvZGqNppGXcmzKvABmPMf6pothAYLy6pQL4xZlcN\n9lMphSuwT/1knSunHbj/69dpUlbClAoPUQGsFtHA3ogFM3I/HbgGWCMiP7mPPQB0ADDGzAA+xZUG\nuQVXKuT1Nd9VpRo33xE7wIDtaxm99iumDx5DZkv/QmA6FaOqDe7GmG9x/QvvcG0McFtNdUopdYgt\ny05aZh47/jjoDexh5Q4e/fwlcmJbMX3wGMD1lzTcKlyR0p7L+iVpYG/kdIWqUvWYp2ZMqcPpN+ty\nre0TTtybxc2X/ZPi8CaAaxXq5It7a1BXgAZ3peo1T80Yp8G7i3Xrwr3c/d1cvuwygKVdB3nbjh3Q\nQQO78tLCYUrVY6nJLQmz+M+K/vOrVwlzlvPwebf4PURdvHaXrkJVXhrclarH+neM44qUQ1nGp237\niYt//YYXU69ge4s2fm2/27JXywwoLw3uStVzvdu5tsiLcJTx6NIZbGvRlpcHja7UzmmgzKFlBpSL\nBnelDsOWZeeFZVvqdDT89cY9ANy06mO67MthyrBbKQmL8J4Ps7jKEFgFwsMsWmZAAfpAVakq+Waq\nRIRZmHNT6nF9YGnLsvPU4g38sM1OUv5u/vb9eyzufhrLk/t725zfqzW3nNUFcD18TU1uqQ9VFaDB\nXakqfbg6h5IyJwYoKXNy37yfGZTc0ptD7sk/r42AasuyM+bl7yl3pbUz+ctZGIGp597sbWMRuMWn\nTrsGdeVLg7tSAdiy7Myz5XiyDzHAltwDbMk9wLurtnPxyW1Z9Msuyp2GyPCaG9V7vjAW/bzTG9iH\nbl3F+ZvTePLs69gVm+DXPi0zT4O6CkiDu1IVzZlDt79P5Nc9u9gZG8/TQ8az0Gdno3KnYf5PO72v\ni8ucfLg655hH876bXHu+VCLLSnhk6Qw2t2zPaymjvG0F1zy7zq+rqmhwV8pH5vMzaX/fncSWFAOQ\nVJDLtCXTAfwCfEUfZGynT7vmTF20jpIyJ1Z3bZcj2bPUs2DJt1zqX9Pm0SF/N1eOe4IyazgAYVZh\nTEp7RmuJAXUYmi2jlJsty07k5IcIdwd2j2hHCRNXzD7stWXlhte++8076nY4DZMXrD2iLJuKC5Y6\n2ndya/o85vc6i7QOJwPQJjaS9yYM5olLT9LArg5Lg7tSbh+tzqFtQW7Ac+2qOO5ry579fqNuh9N4\nc86DTalMiHXVicEYHln6MiXWcB4feqP3/B3naglfFRydllGN3tz0bBav3UV2XhF/RMVwwsHKG40J\n8M479/N5t1Q+7zaYHc1bBXXv5Rv3UHiwjFe+/Q2nMZVSKj1z9HHRETy0YI33IeoFm1Zy9m82Hjn3\nZnKbnQDAJae2O6JpHtW4aXBXjdrc9Gwe+HgNAIOy1xBTXEi5CFZzaAxebA3n68796fTHLqZ8OYsp\nX85ibesufNYtlc+7D2ZjfMdKG2V4/LDNzg/bDo3WS90rSPt3jGNuejaTF6zFaQwi4g3sUaXFTP5y\nFutbdWZ2v4sAaBEdznNX9q2l/xVUKNLgrhq191ZlA9B1bzYzP3qMrLhEXhkwittXvk+7gr2VsmU6\n2ndy/qY0Lti8kru/ncs/vp3DthZt+bxbKp91H8yP7XrgtFgP+55x0RHMTc/mn/PXuKo9Avh8mdzx\n/bskFubyt5ETKXffa39xGbYsu07JqKCJMYffylREXgMuAvYYY/oEOH82sAD4zX3oI2PM1OreOCUl\nxWRkZBxxh5WqCbYsOx+uzuHdH7KJL8jjo7fvIaLcwWXXPENO89ZB3SNhv51zt6RzweaVnJb1M5Hl\nDnKjW/BF14F81n0w33c8ldKw8ErXifvHGeCeXfZuZ8nrt/Nx76FMHHGX97hV4O/n9+C2oV2P7gOr\nkCEiNmNMSnXtghm5vwFMBw6XLvCNMeaiIPumVJ3wzG8XHixj1re/Ue40NC0p4vV5jxB3sJAxV00L\nOrAD5DaL491TL+TdUy+kWUkRZ2dmcP7mNC769RvG/fI5+yOi+Do5hc+7pbKsSwqFkU0B14KogEMq\nY3h06UsciIhi2tmHdqoUtGaMOnLBbLO3QkQ61X5XlKo9gRYIhZU7eHHBNHrkbuPGy6ewrk3Vo+Ke\nbWJo1iSMXX8cZMcfxRhcy//DrRZKHE72R0azqOcQFvUcQoSjjMHZv3DBppUM2+IK9qWWMFZ2PJnP\nug9maddB3oekACPXLWPiitkkFuQiwPt9zmNftLsSZJiFy/snaU67OmI1Nec+WER+BnYC9xhj1tXQ\nfZWqEWmZeRSX+UyEGMPjn73AWb+tZuKFd/gV4/KV0CyCu4f18MtSsWXZ+Wh1Dh9kbKfEUXlypTQs\nnOXJ/Vme3J8HzV/pu2Mj529eyQWbV/LEZy/w2Gcv8mO7HnzePRWHWPjHN3OIdpR4r7/o12/4ttMp\nFF9xpV/tGKWORE0E99VAR2PMfhEZAcwHugVqKCITgAkAHTpoSpc6fuKiI/xe3/ndO4xds5TnTxvH\n+6ecX6lF+mJGAAAbf0lEQVR9dISVS/smBtxoun/HONIy83A4D02uNIu0cqCkvNJ0ixELq5N6sjqp\nJ9POvp5ue7O5YPNKzt+cxv1fvxGwr9GOEu7/5i12Tn9QA7s6ase8iMkYU2CM2e/+/VMgXETiq2g7\n0xiTYoxJSUhICNREqVqxdme+9/crflnK3d/N5YM+5/HsGVcFbF/iKOedH7Kr3NkoNbmlt4Z6k3AL\nb94wqPocdBE2J3Rk+mlXMvLa5zjtL68FnnsHWufn6q5K6pgc88hdRNoAu40xRkQG4vrC0K1gVJ3y\nZMMIEBMZxpK1uwAYkmnjic+ms6JTX+6/8PYq89M9OedlPnnpvvp3jGPOTamVioS9t2q734j+cHbG\ntmJHbAJJAVa/7oyNr/K9lQpGtcFdRN4BzgbiRSQHmAKEAxhjZgCXA38REQdwELjSVJdfqVQtsmXZ\nGTdzJaXl/n8Me+/eyosLprE5vgN/veR+HFb/P/4igDmUyVJdlkr/jnF+gbd/xzimjurD5AVrKwV4\ni/veFWfonx4ynmlLpvvNuReFRfLMWeM1Q0Ydk2CyZcZVc346rlRJpeqFtMw8yioE9nYFe3ht3iPk\nRzbj+sunsD8y2u9814Sm3HBGMlMXraPM4arqeEVK+4Bz7odz1aAO9GgTQ1pmHj9v/4Ol63e7Uh+r\nGO54FkdNXDHbu2jq+xv+TreRl3ON7qqkjoGuUFUhJzW5Je6BMgCxxft54/2HiSorYfTVT7M7pvIj\noWz7QXq0iQk41XKkPCN6W5adFZtzXV8WVgsYg6PcVBq9L+w91BvkwyzCe7cMZowGdXWMNLirkOFJ\nUfx2815vAI1wlDHzo8foZN/JtWOmsjmhY8Bry8td89u3De1aY6PlivPygLdI2Ncb93hH9R6eGvA6\nWlc1QYO7Cgm2LDvjZrk2s/YQ4+SZT58ldfta7rj4HlZ2PLnK62trfjvQvDy4pm88X0Z7CktoFRN5\nxFNASh2OBncVEjy7GPm6b/mbjNywgmlnXcfCXmcHvM4qcOXADnUSWCsGfqVqkgZ3FRLSM/2zb69Z\nvYhb0z/krb4jmDFotN+5nm1i6NcxDgO6rF+FLA3uqkHzTG2s2LzXe2zY5jQe/mImS7sOYsp5t/jl\nslsEHtMt6lQjoMFdNVi2LDtjXv7eu+AI4NSdG/nvwn+xpk1X/jby3kq11U9KbK6BXTUKGtxVg+Mp\n3fv+qu1+gb2jfSevznuEPc3iuHH0ZIrDm1S6duwArWmkGgcN7qpBCVS6F+CEonzefH8KAlx7xSPk\nNW1R6VqLQI82Mcetr0rVpWMuHKbU8ZSWmVcpsDcpK+bVeVNpsz+Pm0Y/xLYTEgNfbFzXK9UYaHBX\nDUpcdIRfYLc4y/nvJ89wyq5N3HnxPaxO7Am46sI0Cbdw65BkwiyCBYgI11otqvHQaRnVoLy1ctuh\nF8Yw5cuZnL85jSnn3cJn3U8DIMLqXxdmWO82x1xSQKmGRoO7ajDmpmez4fdC7+sJP3zEtav/x8sD\nL+PN/hdjAcYNqrwgSRcLqcZIg3uI8WSShMIo1fezLF33OzNWZHrPXbx+OQ98/TqfnHgm086+DoAJ\nQ5KZNKJnHfVWqfpFg3sDFSiIezJJSh1OIsIszLkptcEG+Lnp2UxesJZyp8Ei4FvBd1D2Gp759FnS\n2/fhnj/djRHXo6OVmXm8sGxLSHyxKXWsNLg3QL6BLzL8UBD31FdxGihtwLv42LLsfhte+Ab2rnuz\nmfnRY2S3aMuESx+kJOzQ3qjrduazZkd+g/9iU6omVJstIyKvicgeEVlbxXkRkf+KyBYR+UVE+tV8\nNxW4gt6DH6/hoflrcDgNBiguc3LfvJ+xZdmJi47AswGQ01TeFDrQ/V5YtqXe7dOZlpmHM8DuFq0K\n83jjgymUhEVw3RWPkB/lylnvmtCUYb1a4zSuz+3Znk6pxiyYkfsbuHZaml3F+eFAN/fPIOAl939V\nDbJl2Zl715Pc/dUbPOresefpIeNZ2HsoW3IPMHbmSob2aOVtb8G1KXRV0xT1eQonLjoCi4hfgG9a\nUsTr8x4h7mAhY66axo7mrYhtEsak4T295XO/cW+ModvTKRXcNnsrRKTTYZqMAma7901NE5EWItLW\nGLOrhvqogK3Pvcyji/7r3WszqSCXaUtcuxsu7D0UR7nhyw27/a75IGO7d7u5mEgr5/ZsTbfWMaQm\nt+Sj1TnexUAlZU6mfrKOyRf3PmyAPx4Pa21ZdqYuWue3B2lYuYMXF0yjR+42brx8CuvadHUdt1q4\napCrnEBVG1Yr1VjVxJx7IrDd53WO+5gG9xpiy7Jz+mv/8dtEGSDaUcLEFbO9W7T57snsBJw+k9WF\nJeXM/2kn4FqG79vWAD/n5DP25e9575bTAgbG4zXS96xAPdQ5wxOfTees31Yz8cI7WJ7c33uqa0JT\nv2s15VGpQ47rClURmSAiGSKSkZubezzfukH7cHUObQv2BjyXWJDL0K2riC3eH/T9nFVs1uxwwsvL\nt3pfe+bk56Zn89wXm7wPa2tzTjs1uaXfCtQ7v3uHMWu+4PnTxvH+Ked7j1sE7huuaY9KVaUmRu47\ngPY+r5PcxyoxxswEZgKkpKRUEWJURQLsjI0nqSDwF+Lr8x4BYGN8B2yJvViV1IuMpF5sb97ar5Z5\nML7csNv7gNUzUvf9MrBI4C3pjnXKxnN94cEy77ErflnK3d/N5YM+5/HsGVf5tQ+zauUMpQ6nJoL7\nQuB2EXkX14PUfJ1vr1m92zXn32f8mf98+iy+obooLJKHht3CjhZt6J+zgQE567lowwqu+nkJALub\nncCqxF7YknqyKqk3G1p1prxCffOKnAZmLN/K5t2FlQp0AbRtEcV/r+zrF8CPdcrG93rPF8mQTBtP\nLvk/VnTqy/0X3l7pS8qzobVOwygVWLXBXUTeAc4G4kUkB5gChAMYY2YAnwIjgC1AEXB9bXW2sVq3\nM589MS0RIC8qlriDhX7ZMgBpHVybP1uc5XTfm03Kjg2k5KwjJWcDF238FoAD4U34sV0P7+j+x3Y9\nOBAZ7fdeBli63v/BrK8d9oPMWL4VAeJjIhndL4kZy7dS7J4n952yCXYk75ufD9B791ZeXDCNTQkd\n+esl9+Ow+v8xFWpvQ2ulQkUw2TLjqjlvgNtqrEfKjy3Lznursnl443ccCG/CaX95nZLwyCrbOy1W\nfm3VmV9bdebtviMAaFOwl5Qd60nJWU/Kjg3cvvI9rMZJuVjY0KozGYk9yUjqRUZiL36Pja+2T77B\n/930bHy3pXYa2Ly7kGeXbsJpTFAjed/8/MT8Pbw27xHyI5tx3eUPs9/95SPun7AKRcGUUoHpCtV6\nLi0zD+Mo54JNK1nWZcBhA3tVfo+NZ1HsEBb1HAJAs5IiTt25kQE560nZsY4xa5Zy3epFAOTEtiIj\nqScZSb3JSOzJpvgOlbaq8+Ws8NqANysHKq+UrTg370l9BIgt3s8bH0whqqyE0Vc/zZ6YQyPzW4Yk\nExMVrmmOSgVJg3s9l5rckoE715NQ9AeL3SVtAe/ce3iYhRF92vDpml2Ulgf3jHp/ZDTfdu7Lt537\nAq488p57fvOO7gdnr+GS9csBKIhsyup2J7IqyTV3/1Pb7gG3r6uKiFB4sIxR07+l1OFk0579GJ8R\nvSf1McJRxsyPHqPDH7u4dsxUNid09LtPTFQ4tw3tGvT7KtXYaXBvAC7NTKc4LIJlXVK8x8YN6kBi\niyjvSPaawZ0YM+N7bx0WAVI6xbE1dz/7DpQFvrGbwxrGmrbdWNO2G6+njAJjSMrfzYCc9QzIWU//\nHeu595u3ACizWFnXuosrI8f9sHZv06pH0uVO41fN0aO4zMlHq3O4rF8SFpw88+mzpG5fyx0X3+t9\nfuBhFXR+XakjpMG9HrNl2bnq5e9YvmYFXyf3pygiCgCrRRgdoGb5o5ecxOQFa71z3ZPceeDjZrky\nUYImQk6LNuS0aMPHfc4BoPnBQvrt/JUBOevon7OB8av/x82r5gPwW1xbMhJ7u6dzerH1hKSgUjDn\npGezbkc+9379JiM3rGDaWdexsNdZldo9eslJOhWj1BHS4F6Pfbg6h97Zv9Jm/z4+7XG69/jYAe0D\nBrurBnWgR5uYSlkq79ycynNfbOLbzXsrpTYGKz8qhmVdBrCsywAAIhxl9Nm9xZWCuWM95279gSvW\nfgHAvqhYbIk9yUjqyarE3qxt05XSsHDvvUauW8bEFbNpV7CX/CZNiSvez+y+f2LGoNGV3veJS0/y\nlhhQSgVPg3s9trewhBEbv6XEGsZXXQYCEBFmYXS/pCqvCbQEv3/HOO46rzurtu2jzOHEarVwVvcE\nBNhdUMzPOflH3LfSsHBWJ/ZkdWJPZnEZGEPyvh3ujBzX3P2wLekAlFjD+bltN2yJvbA4yxn/4/+I\ncpQCEFe8n3KxsLpdj0qj/VuHJGtgV+ooaXCvx7bnHWDypu/5plNfb0rgKUnNj2qKoqrCWrYsO1fM\n+L7KkgRBEyGzZRKZLZO8ZQJaHviDlB3rvQusblw1nwino9KlVuPknm/eYr57CgigZ5sY3VVJqWOg\nwb2esmXZafLzapIKcnn2jD97jx/R3HkFVY3qJ5yZ7PfQs2ebGDbvKeQY3gqAvKYt+Kz7ad6Nq5uU\nFbP+P5cHLGjUzqd2jtUiPHbpScf25ko1chrc6yFblp2xL3/PPb9+R5nFytJuh8rjjx1Q89MUk0b0\npEPLpixeu4vhfdp666N7RvngKii2bmc+O/4oPur3KQ5vws7YhIA1cnb6LJ6q6pmCUip4Wn2pHvpw\ndQ6OcsPwjd/xfcdTKGjSDIBLTm1Xa3PQVw3qwFs3DvKrj37b0K7e0f7M8SlcNahjNXep3tNDxlMU\n5r8QqygskqeHjAcgwiqHfaaglAqOjtzrmbnp2by3aju99vxGxz9+58XUK7znurWOqcOeuXLNm4Rb\nKC1zIgLn9mzNLWd1AVzFxr7csLvauXtPLRxPtoxvjZxhvVpz61lddNSuVA3Q4F6PeDaGLne6Ru0O\nsfB5t1Tv+bpeyHO43Y5mjU/BlmXnw9U57C0s4csNu6m4YDbcKpzdoxX/swz1BnmP2CZhzBqfglKq\nZmhwr0fSMvNc28sZw4iN35HeoQ/26OaAa0qmPoxoD7fbke+5Bz9ew9z0bAyuGvCnd43nrvO6e+vJ\njH81nQOl5d5rE1tEHY/uK9Vo6Jx7PRIXHQFAt73ZdNmXw+IeZ9As0sqtQ5J57sq+ddy7I3NZvyQi\nwy1YxZWb7wns4PoSuKRvol/7fvXgi0upUKIj93rEXuRa2DNi43c4ET7rNphebWMbZL734fLq0zLz\n6N2uORFhFsocTsLDLFymD1GVqlEa3OsRz8h9+MbvWJXUi9xmcbQ71mTzOlRxCqfijk0PX9wbe1Gp\nlvFVqhYENS0jIheKyEYR2SIikwKcv05EckXkJ/fPTTXf1dBnLyolOS+HE/dmsdhdS6Y28trriu+O\nS2UOJ/aiUm+6pVKqZgWzzZ4VeAEYBuQAq0RkoTFmfYWm7xljbq+FPjYKtiw7P23/gxGbvwfgs+6n\n1Wpee11ITW7pNxVT19k/SoWyYKZlBgJbjDGZAO6NsEcBFYO7Okq2LDvjZq6ktNxw56/fYWt3Ir/H\nxrNk3e/YsuwhM7I9XCqlUqpmBTMtkwhs93md4z5W0WgR+UVE5olI+xrpXSORlplHWbmhg30XfXZv\nZXGP0zAc2qIulPiufFVK1Z6aSoX8BOhkjDkZWAq8GaiRiEwQkQwRycjNrVxfpLGKi47AAMM3fQfA\nEvd8u1D3C5eUUg1TMMF9B+A7Ek9yH/MyxuQZY0rcL18B+ge6kTFmpjEmxRiTkpCQcDT9DUkf/5gD\nwPCN3/Nzm27kNG8NHH7BkFJKHU4wwX0V0E1EOotIBHAlsNC3gYi09Xk5EthQc10Mfb/+Xkhi/h5O\n3bWJJT0ObYLdwp0aqZRSR6raB6rGGIeI3A58BliB14wx60RkKpBhjFkI3CEiIwEHsA+4rhb7HFJs\nWXb2FzsYs8mVJbO4+6HgHh8TWdVlSil1WEEtYjLGfAp8WuHYZJ/f7wfur9mu1R+eVZVx0RE1vujG\n88D0wo3fs75VZ7ad4HpWbRW09K1S6qjpCtVq+K6qdBrXQ87IcAtzbkr1FsE6ltS+1OSWJBXlMWDH\nep4507XjklXg0UtO0vl2pdRRC5ngbsuy89HqHAzQp13zGhth+66qBDBASZmTGcu30iomkg8ytuNw\nGiLCDgX8iv2qqr6KJxNmTI4NgMXdT0eAKwd2CKnFS0qp46/BBnffALnx90IectdB97C4qxFOvqhy\n/RJP3XHBVb2wYv2TtMw8Cg+WsW5XAS2bRiAiYA7d2wBL1+/2609pmZPnvtjE8D5tve8H+NVSmXOT\nqzb7uJkrKSs3WC2ue83J+IpNLTuwNd6VlNS7XfNa+F9MKdWYNMjgPjc9m8kL1uI0hjCrBUe5s9IO\nQE4DxWVOHlqwFmMMFhGmjuoDwD/nr/G2fz9jO1ektKdPu+as3ZnPPFsOZQ4ngTYUEgh4HMAJfLN5\nL99sdm30bAGSWzWjpMzpHe1P/WQdrWKbUOrexcLhhPgDdgbkrGf64LHee3mqQyql1NFqcMHds1uR\nwx2dS6upmugZzTuN4Z8fr8HgH6DLyg1z07OBwwdvqjlXkRPYsme/37U/5+QD+X7tzt+chtU4WexO\ngYzQmitKqRrQsIL7nDl0+/tENu3Z5bf3ZrCqK557JMH7WI1ct4yJK2aTWJBLmcVK99xt7OncnVnj\nB+iDVKXUMWs4wX3OHJgwgdiiIgCSCnKZtmQ6wBEF+LogxkmYsxyrs5wwp5OLNixnypeziHK4pl/C\nneVMWzKdpSe2pn/H8+u4t0qpUCDGHM/x6iEpKSkmIyMj+As6dYKsrEqH90XF8tg5N2J1lhPuDaCH\nAqnfa1OO1Xko0IY7HX6vD7WrcJ3TdZ2r/aHXVlO5je97hZW7jluC/DdBQat2xO7eUX1DpVSjJSI2\nY0y1u8k3nOBusfhlrBypcrHgsFgot1hxiNX1X4uVcosFh8WKwxLmPmY5dE782xy6xkqZ72upcJ3v\nvX3ey2F1vfc/l72KBOijEUGcDXfnJaVU7Qs2uDecaZkOHQKO3H9vdgJXXP105UDqE2jLLRaM1J+9\nwK+3LSSpoHJVzLy41sTXQX+UUqGn/kS86jz+OGWRTfwOFYVF8sTZ17O9RRt2xrZiT0xL8pq2ID8q\nhgOR0ZSER+KwhtWrwA7w9JDxFIX5140pCotk4x0hW8FBKXWc1a+odzhXX80P9z9FTmwCToSc2AQm\nXXh7vX+YGsjC3kOZdOHtfp/luSvu4fQpd9R115RSIaLhTMsAPw0ZwZ8Pdqn0eFLElaNecSFTTbFI\n4Hs3jwoj/6Aj4DU928QQ0ySMH7bZK507Nak5nzDU74vpiUtPqrH+KqVUgwruqcktCbeKd4WnhwUY\nO7ADewtL+GLD7iMO8uf3as3S9burzGmpKrh3atmUtTvy8e3OwE5x3De8pzdXfdqnG3h5RSYG1xfQ\nLUOSmTSiJ7YsOy8v38rugmLGDtBaMkqpmtVwsmXcPHVh9haW8PXGPZQ7DeFh/lUaP1qdw9od+fyS\nk+8NquIToAUIswq92sYydkAHerSJ4epX0ihzOBERnMb4BXOrwDk9W/NFhS+AYb1a86XPl8mwXq2Z\nNb7yQ+xjrRyplFIeoZcKGcDhgubc9Gwe+HiN9/WtQ5KJiQqvsia7b832qYvWeStBWoCI8ENFv3xH\n275fCuFVVIVUSqmaVKOpkCJyIfA8rp2YXjHGTKtwPhKYjWvv1DxgrDFm25F2+kgdbo9Re1Gpt1aM\nBYiJCue2oV2DulePNjFVbs4xs8LIfM5NqToqV0rVO9UGdxGxAi8Aw4AcYJWILDTGrPdpdiNgN8Z0\nFZErgaeAsZXvdvykJrckMtziHVUfSTGuI9mYWjexVkrVR8GM3AcCW4wxmQAi8i4wCvAN7qOAh92/\nzwOmi4iYuprzwRV0dVStlGqsggnuicB2n9c5wKCq2rg31M4HWgJ7a6KTR0tH1Uqpxuq4LmISkQki\nkiEiGbm5lZffK6WUqhnBBPcdQHuf10nuYwHbiEgY0BzXg1U/xpiZxpgUY0xKQkLC0fVYKaVUtYIJ\n7quAbiLSWUQigCuBhRXaLASudf9+OfBVXc63K6VUY1ftnLt7Dv124DNcqZCvGWPWichUIMMYsxB4\nFXhLRLYA+3B9ASillKojQeW5G2M+BT6tcGyyz+/FwBU12zWllFJHq+FUhVRKKRW0Ois/ICK5QOXd\nN4ITTx2nWdYB/cyNg37mxuFYPnNHY0y1GSl1FtyPhYhkBFNbIZToZ24c9DM3DsfjM+u0jFJKhSAN\n7kopFYIaanCfWdcdqAP6mRsH/cyNQ61/5gY5566UUurwGurIXSml1GE0uOAuIheKyEYR2SIik+q6\nP7VNRNqLyDIRWS8i60Tkzrru0/EgIlYR+VFEFtV1X44XEWkhIvNE5FcR2SAig+u6T7VJRO52/5le\nKyLviEiTuu5TbRCR10Rkj4is9Tl2gogsFZHN7v/WePnaBhXcfTYOGQ70AsaJSK+67VWtcwD/MMb0\nAlKB2xrBZwa4E9hQ1504zp4HlhhjTgROIYQ/v4gkAncAKcaYPrhKm4Rq2ZI3gAsrHJsEfGmM6QZ8\n6X5doxpUcMdn4xBjTCng2TgkZBljdhljVrt/L8T1Fz6xbntVu0QkCfgT8Epd9+V4EZHmwBBcdZow\nxpQaY/6o217VujAgyl1JNhrYWcf9qRXGmBW4am75GgW86f79TeCSmn7fhhbcA20cEtKBzpeIdAL6\nAul125Na9xwwEXDWdUeOo85ALvC6ezrqFRFpWtedqi3GmB3AM0A2sAvIN8Z8Xre9Oq5aG2N2uX//\nHWhd02/Q0IJ7oyUizYAPgbuMMQV13Z/aIiIXAXuMMba67stxFgb0A14yxvQFDlAL/1SvL9xzzKNw\nfam1A5qKyJ/rtld1w10evcbTFhtacA9m45CQIyLhuAL7HGPMR3Xdn1p2OjBSRLbhmnY7R0Tertsu\nHRc5QI4xxvOvsnm4gn2oOg/4zRiTa4wpAz4CTqvjPh1Pu0WkLYD7v3tq+g0aWnAPZuOQkCIigmse\ndoMx5j913Z/aZoy53xiTZIzphOv/36+MMSE/ojPG/A5sF5Ee7kPn4r8JfajJBlJFJNr9Z/xcQvgB\ncgC+GxxdCyyo6TcIqp57fVHVxiF13K3adjpwDbBGRH5yH3vAXWNfhZa/AXPcA5dM4Po67k+tMcak\ni8g8YDWujLAfCdGVqiLyDnA2EC8iOcAUYBrwvojciKs67pgaf19doaqUUqGnoU3LKKWUCoIGd6WU\nCkEa3JVSKgRpcFdKqRCkwV0ppUKQBnellApBGtyVUioEaXBXSqkQ9P+mwUKtbcZcAAAAAABJRU5E\nrkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x1d425eeec18>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"X, Y = make_test_data(10, 2000)\n",
"px, py = segments_fit(X, Y, 8)\n",
"\n",
"pl.plot(X, Y, \".\")\n",
"pl.plot(px, py, \"-or\");"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXl4k1d69/85kmzZBhuEDcZgLGJW4w0wiyEsYbHJQkLY\nl3Qy7UxCpu2017xtf53MpOXNlTTpTOft+07bmTYhmXYmLYuBACGZhH0LE2ywHbxhzOIgYxZjjDAG\ng21J5/eHpMeSbQiLFy3nc12+kJ7nSHqwz3PrPve57+8tpJQoFAqFIrjQ9fQFKBQKhaL7UcZfoVAo\nghBl/BUKhSIIUcZfoVAoghBl/BUKhSIIUcZfoVAoghBl/BUKhSIIUcZfoVAoghBl/BUKhSIIMfT0\nBdyLmJgYOXTo0J6+DIVCofArCgoKrkkp+3/bOJ81/kOHDiU/P7+nL0OhUCj8CiGE5UHGqbCPQqFQ\nBCHK+CsUCkUQooy/QqFQBCHK+CsUCkUQooy/QqFQBCHK+CsUCkUQooy/QqHwOQosVn594CwFFmtP\nX0rA4rN5/gqFIjgpsFh56cNcmm0OQg061r2SSYbZ1NOXFXAoz1+hUPgUuZV1NNscOCS02BzkVtb1\n9CUFJMr4KxQKnyIzMZpQgw69gBCDjszE6J6+pIBEhX0UCoVPkWE2se6VTHIr68hMjFYhny5CGX+F\nQuFzZJhNyuh3MSrso1AofAKV4dO9KM+/hyiwWO+5rHWfM0WEYm1s1mKeahms8Hc85z3g9did4aMT\ngrcWpLBqckJPXmrAo4x/D+BOZWtqcaDXeU90zzQ3hwSBc9MLKbE5pEp9U/gtnnPboNfhkBKbXaLX\nCeaMHqDNeYeUrPmklFEDI9U870KU8e8BcivraGpxIAGbo3WiA/xy72nutji0sRJotjkQHo/dqW9q\nJaDwJ9qmcErXcbtDsvdkDUInQErtWG5lndfcLrBY2VpYjQQWj49X8/4xUca/mymwWLl04w46AXbX\n7HdIyceF1WwtrKbJw/B74r5RHBIa7rSwcu1RWuySEL1gw+op6kZQ+BxtQ5vuFM4WmwOEwO6Q2lgH\nMDGhL8fPO+P9EjBFhHq918oPnKsGgC35F9S8f0yU8e9k2k74tjFO97JXrxPgkDgk6ITgWkMTzR7e\n0L0QwNHKOppd3xzNducXh7oJFL7Evap03SmcpohQ/v6TEuwevo4pIlRb4eoAa2Ozdi63ss75peGi\n2S755d7TPJMSp+2LqXvg4VDGvxNpO+HXzE/mrc/KtOeLxsdr4R5pl0wcaiLfYnUue8trcHyb5cdp\n/Fvs3quDaw1NXfL/USgeFc8Qz90WB+8dOscHL0/QUjgLLFaeiOnN2au3tNdUWxsxhuicq18BZ2oa\n+PWBs2QmRpOZGI1eL7DZW2+SL89c48sz19AJ1F7YI9Apxl8I8Z/AfOCqlDKlg/MC+BfgWaAR+GMp\nZWFnfLYv0Tam+UXpZe15s83B4YqrmmcvgWPnW1Pa5AMYfnAuj09dbvA6FhNp7JTrVyg6i2eK9/Pi\nv/8dcTevcSkqhn+a8TI/i+lFZHgIDXdaWPtlZTtn5+TlBuL7hlF94y5SwvYTl7SEh5kj+zM0upfX\nl4UbTxkIZfwfnM7y/H8L/Ar46B7nnwFGuH4mA//h+jeg8Ixphhh0JMdFcfRcHVI6wzvVN+52yue0\n3RWIMho0D0lNfkWPs24d5p/8CP2dOwDE36zlZzt/xU8EvDdm1n1f2vYecSc57DlZc8/X6JQMxCMh\n5IO6nN/2RkIMBT67h+f/PnBQSrnB9bwCeEpKefle7zdhwgSZn5/fKdfWnXjm6LtDPpIH9+wfBb1O\nIKVKA1X4CEOHgsXS7vCV3v3I/PN7+YePhk7A3KRYXps5TM17F0KIAinlhG8b110VvoOBCx7Pq13H\nvBBCrBZC5Ash8mtra7vp0jqXDLOJP581nLJL9TS1OEM+XWn4wZkWN7/0APv+9buMfyLaefOtW9e1\nH6pQ3Iuqqg4PD7x1nT0f/ik/OfCfZFYVY7DbHvujHBL2n7r62O8TjPjUhq+Uci2wFpyefw9fziNT\nYLHS+Nvf8eXBjxjkEfPckXz/Je+j8kLZAX6281dE2FwbvxYLrF7tfPzSS13ymQrFPUlI6NDzvxHW\nm8uRMfxJ/g5eO7aVm8ZeHB46jgPDJnIwMYO6Xn0f6eNsDslbn5axfGKCyvx5CLrL+F8Ehng8j3cd\nCyjcIZ/QnA2883mrMXbHPIGH+gIQ0kGorQWjvYVQ94/N/dimnVuz/8NWw++msRHeeEMZf0X382d/\nBj/+sdehRoORNXNfY0fyLHo1NTLNcoJZ5/KZVZnP/IojOBAUx41g37CJ7B82kbLYYSDEA39kUXU9\nRdUlgDPzZ8OrKvz5bXSX8d8B/FAIsRHnRm/9/eL9/ohnmueXG/6tnTGOsDXx852/YkH5Yc2AGz0M\nurGdcW8hxGF/vIu6x/JboehSLl/GIXRc6d2PgQ117Va+t40R7Bo5lV0jpyKkgzE1lcw+d5zZ5/L5\nX0fW89dH1lHTux8HEidwYNgEjpjHctsY8cAf32xzsFXVvnwrnZXquQF4CogRQlQD/xsIAZBSvgd8\njjPN8yzOVM8/6YzP9SU80zzjbl7rcEyYrYn+t60060No1odwIySMZoPzcZPe4Dzueq79GEJo1hk6\nPu56/C87/okBjTfaf2CCEsZSdC85h07z3Nr/5HTmXJbM/MtvrV2RQkfZwOGUDRzOvz25kujbN5j5\nTQGzzx7n2VNHWFG8m2adgbwhKRwYNpF9wydiMQ361us4XdPwrWOCnU4x/lLKld9yXgJ/3hmf5au4\n0zybbQ4uRcUQf7P9hvXFqP688N1fdvpn/8Ps73vH/AF7WDj6d97p9M9SKO7F+rwqjv/ifZY33uQX\n5qceqGixLXW9+rI1ZQ5bU+ZgsNuYcPEks87lM/vccdbs/4A1+z/gXL/BHEicwP5hEzk+JJkWfUi7\n98m3WPnZ5+VEhoeoPYB74FMbvv5MhtnEmvnJ/Ou+02xLfoq/OLrZ63yjwcg/zXi5w9f2DTdw486j\nZz64l9N/e/gjBt2sRQdsHjGN86bxvP7I76pQPBw5x6t4o2gnlaZBHE1Ifez3s+kN5CakkZuQxj/O\n+h5DblxxhYeO852vP+eV/E9oCA3nyNBx7B82gYOJE5liOeG6D65x6d9j+MWMl/m3sXNUCnQHKOPf\nSRRYrLy5o5Rmu2TM1fPcDAmnIawXcR3EPNtSf8eGXuctdPWw7Eie5Xx/Kfniv/6C0dfO8/rhSgBe\nfzbpkd9XoXgQCixWmoqKmVR9knee+h4Ioen0PCh6p9zVPV9zoe9AfpfxPL/LeJ7w5rs8aSli9rnj\nzDp3nGdOfwWAA4HO9Q7xN2v5x52/QgjIrRyhjH8blPHvJN4/dI5mu2T4tSrmnDvOP097iX978r7R\nMA0J9Ak3cP12y+NfiBDkpGXz5r61JF2tZO2XkJU8UE18RZeytbCaZV/voklvYEvqHAZGGontE0bJ\nxfpvDf8IYO6YWH4wcxgVVxr4+09KvRwhlwaiF3dCw9g7YjJ7R0wGKUmq/YaN639Cn6bbXuMibE38\nf4c+4t//5LsUWKzqPvBAGf9OYH1elVZ+/uqxbdwxGPmfcc8+1Ht0iuF3sS15Fj85+F8sK97DW3Nf\nU5onii7BU7HWWlfP4tJ97Bo5FWtEH2ho4koHgoMTh5o4V3vLa76vnJzAuwudYaIMs4lRAyO9Otk1\n3GnhgyPf3HtlLATlAxKJbGrs8HTczWtsOFbFx4XVKvzjgerh+5isz6vip9tKkED/W9d58eQBNqfO\ndd4Aj8jw/r0I0bfPcdYLMHRwvC314ZHsHDmVRaX7CW1p8tJFVyg6A3dq8z/vruClD3MZc2QXfZpu\ns37s0/d8TahesHBcPA13vfe3UgZ53ytu7X9rYzOmiFB+e/Q8Ujo7fhn0Ar1wvtdLkxPoG97qv16K\niunwcy9FxXiJvymcKM//MXAbfjffLfyMELud30xc8K2vDQ/R0Wx3IB3eQm2hesHkxGgqr3kvXw2u\ndo+jBkbycWE1Occv3HePYGN6NgvKD/H06a84WJGgeVMq80HRGXh2o2tqcTB571bO9RtM7pB7b/R+\n78knsDY2e8kyC7x1+8G7ZkYnBA6XMKIeyfKJCQzqG67N4+RBfbR78J9mvNwu661JH6IlWijxN2+U\n8X9ECixW/n57q+GPaL7DH339ObtGTrlnHrJw7YBJoMnVx3TJxHhSBvWh7FK91p4O4OPCalpsDvR6\nHUsy4rXjuZV1LB4fT5TRwPuHK++5OZabkMr5vnGsKN7NyuRZ7CmvQaB0zxWdgykiVJt7w2stTLx4\nkn+Y9b37VuUeraxjzfPJhLhSogFC9KKdQfasmUFKdDqBQBLi6onhHgNo7U+hbdbbNaQQnI5OYEfy\nLIb378XPl6Sree+BMv6PSG5lndcm1NKSvfS9e4u1kxbd8zV6AcmD+2ibYHa7g8F9w7Xm7Z64Ox55\ndgTzbH6NlPfNpJBCx6a0LP728EcMvX6R8/0Ge/UAVjeB4lEpsFj5orS1QH9V0U6a9AY+Tplz39fF\nRoWRYTax4dXM+/bibSuNvmZ+sqbZA3jdBwmmcK/XallvwN8e+i2r87bS/9Z1vtEJKq40qHnvgTL+\nj4C7D69eBzYH6B12Xjm+neODx/D14NH3fJ3d4TT+FTUN2sR+0GVo20YxD8KWlDn81Zf/w7KSPfzT\nzD8GWnsAKxSPgtsJcfeaDmu5y6LS/Xwx6sn77nOF6AWvzRwGoHXzuhee7R7bhil/feCsV4Oks7W3\n7/k+m1Oz+LPcLSwq28/7k5fwd66VekfOVjCijP9D4umBu3m64iuG1Nfw9uxX7vvaEL1g8XhnCOd+\n8feO+p96ekN6l+dvd0iEzru13WBTOJdv3MEh4WpkNAeGTWRJyT7+77Q/wqZ3/rnfO1xJQnQvdRMo\nHhq3E+Kecc+d+gN9mm6zIb11o7dfRAj9ehtJjOnFU6MGPJLS5r2+INz3gXu/4X58028wx+LHsKx4\nL+9PWowDwd9/UsqogZFqBYAy/g+NVzwSQEpePb6VStMg9g6f5DU2Pb4PyycmeMXz3ZPufpOvrZef\nW1nHn88a7uUNucdlJkZTcaWBL0ovkxwXxc0mG1sKqrUvp43p2WSdzWNWZT57RmRqn/FF6WVl/BUP\nTVvju+rEF5zrF0/ekNYeTtcbW2hssfPzxWmdbmTdq4KthdVszr+gOUB2e8dh0M2pWfzii39h/MVT\nFMYnYXdIJfrmQhn/h8QUEeoV659UXcbYy2d4I/vPcOj0XmNjo8IeycC2jXm6jX1bb8jzi2TUwEiP\nDInW9zqYOIGa3v1YXrTLy/g/kxL30NelULiN799tK8FRUkLGpVO8Pev77TZ6m1u6bm/JfR8s8lhB\nd1QcBvD70dP43/vWsqxkD4Xxzkp3v20U0sko4/+QtE1Le/XYVurCo9jSwWbXo06y+8U874XnasFz\n/tt1ejanzuVPc7cw8OY1rkTFEGbQeWVJKBQPSoHFys+/KKf8SgNvnthJkz6Ej1Pbz30HdHl9iacz\n5Fkcdqamge0nLgHQGBrOZ6OnM//Ul7w151Waw8IRrv9HsHv/qsjrIfHcLB127QJZZ4/x0fj5NIUY\n240dENn+2IPibgf5oBM0MzEag67jNLtNqVnopYMlpXsBuGtz8NKHuRRYrI98fYrgY31eFcve+4pj\n563Ojd6yA3w+6kluhEe1G6sT7R2lrsZ9z/xyxTgmDW29bzalZtG7+Q7PnTqC3QEbjlWp+Y8y/g/F\n+rwq1n5ZqT1/5fg27hpC+e/xzwHOghU3Br3QcpK7gwyziaUThmjXoMO58QZQZYrjiDmd5cV7ENK5\nF6CqHRUPg7uuxZ1bMP/UEaKabrPhHhW9Br3vFFQVDh7NuX7xLC3Zg8S5Mm5qUfNfGf8HpMBiZc0n\npVpIpf8tK4vK9rMlZQ7XPVLc9DpB1phYclZP6fZl5aLx8RhDdM7y9xAdyya0ds7MSctmSH0NUy3F\nAOh07YtrFIp7kVtZh0dSGatOfMGZ6CEci08GnMKEL01OwKBzOUGy5yLrBRart1cvBJvS5jKp+iRP\nXHd2j5V0fVjK11HG/wFxFnW1TujvaFIOL2rHJICUjB3St0fiie69gr/KHsW6VzKJDA/RVgK7R07B\nGhbJiqJdAJj7RQR9zFPx4GQmRmtzafTVbxh/qcKZ3una6J01agCD+oZrksx2h+wxz7ptASbA1uTZ\n2ISOpSV7tGPdHZbyNZTxf0DcGTgCCG++y3e+/pw9IybzTb/BXuP0Pbzc9dwryEyMRu/aB2gyhLIt\neRbZZ45iaqynX6/g9noUD0fFlQZtT2llkWujN2W2dr7udrN2j+hFz+roZCZGa8KIOgF6HdT27seB\nYRNYXLofvas3drAXOyrj/4C4O3XpdIKlJXsw3W3g/UmL241bktG+XL0nmT16gPY4Jz0bo93GwrKD\nnLhwI+g3vBTfToHFyuqP8vnpthJaHJLw5rssLD3A70dPoz68NWPsmZS4divPHr0PXM1kwFlZD7Ap\nLZvYW9eZ8U0hAGsPVwb1PaBSPb+F9XlV5ByvwmjQcenGHaTNxivHt1MwaLSWN+xJW3nansKzSlgI\nZwi2ov9Qvo4bxfLiXfzXhBf4WBW7KO6Dew7dbWmtZp9/6kuimhu9pJsFrQJr3ybd0B3kVtZhszuL\n0Dy3Hg4kTqA2oi/Li3dzYNhEHBDU94Dy/O+DW7K5qLqeY+etVN+4y7zTR0mor2HtpEX06xXilVKm\nw3fiiN7KiK3HN6ZnM+paFeMunWLjsSpWf5Qf1N6P4t6455Anq07s5Ez0EPIHj9GOCfCpzBnPEK0n\nNr2BrSmzmXP2GNG3bwBwrYOGM8GCMv73Ied4lfcBKXnt2FbO941jz4jJPJ0Sx4+fSSLMI8PGVzJo\nPOOvxhAdP5iRyPQRMZyZ9Sy3QsNZXrQbh4TdJ2tY9v5RfrqtRH0JKLzITIx2Ksi6SLpaybjLFU6v\n37XRK/CteQ+tiQ8rJye0a360OXUuIQ47L5YdACDmMWpx/B0V9rkPsVFhQL32fMLFk4y9fJq/y/pT\nHDo9KYP6PFI1bndwr+v66bYSPh09nQXlh3h7zqvcMkZgd0g25FWxVbW5U7QhIlSvef8rT+yiSR/C\n1mTnRu/Q6AiWThjiU/PejTv8JIB1ea1O3NmYBAoHjWJ58R5+M/FFrU9GMKI8//vw2sxhXjo5r+Vt\n5Xp4FFtc5ew5x6u0MvGHqcbtLjq6rsXj48lJn0dESxPzyw9rxyWq8EvRSoHFyoq1R7nR6MyICW++\ny4tlB/jMY6P3fF0jpohQn5v3niwaH4+hjZXblJrFyLoqxtWc4edflGsr3gKLlV8fOBs0K2Bl/O9D\nhtnE3KRYABLrqsk6m8d/j3uOuyFhABRV17Ny7VG/miwZZhPXktI4FWNmRfFur3OqzZ3CTW5lHS0e\nVV3Plx92bfQ+4zXOs6mLL5JhNrF8ore44mdJM2gMMbKkaDfHzltZ/B9f8bPPy716EvvTPf2oKON/\nH9bnVbH3ZA0ArxzfTpM+hI9cUg5uWuw9V8zysLg9m8jwUHLSsxl7+TSjr34DQJ8IA2vmJ/u0F6fo\nPjxrRMDZret0dAIFg70z3PxBHXbR+HivFfwtYwSfj5rG8ycPE9ZyF3D2uLjb4giqRu/K+N8Dt5aJ\nA4i+fYPFpfv4OGUOdb36eo3TCfzCW3an7f1iVwXlVxrYljyLJr2B5S7vv77RxpuflgWFx6N4MNx7\npWNqKhl7+bS20Stw9qp4d2GqX/SEyDCbWD090evYprQsopobeabiq3bje7pQs7tQxv8eeGqZvFz4\ne0LsNj70kHIAZ6bD2y+m+oW3nFtZp7XeA7gRHsWukVNZWHYAo82ZntocJB6P4v4UWKz8cu9pLeyz\nsmgndw2hbHVV9E4cauKTH07zC8Pv5vVnk3h3YSrDB/QG4Fh8Mt+Y4ljmIffgxtcKNbsKZfzvgSki\nFL1w9ij9zte/Z++IyVRGOzMDBGDQCd7xE88H2i/jATamZdP37i3mnW71foK95D3YKbBYWflBLl+e\nuYYEIprv8GKZs6L3ZpjTcPbxU0G0VZMT2PtXM/nBjESETrA5NYspVSUkWL33LXylULOrUcbfhedO\n//q8KtZ8UopdwtKSffS7c5ON05cyaaiJdxem8jfzRpHz2hS/MfzgXPq+tSAFg6617P2oOQ1L34Gs\nKGrd+F37ZXCXvAc7WwurvQq7ni8/TGTzHdalt270Pk6fCl8gMjwEpFPszS50LCnZ63W+9FL9PV4Z\nWKg8f7ylEAx6HXa7A7sEncPO949v5+u4UeyPGUnYxXp+/EyS3y4JV01OYNTASH6597TTsxM6ctKy\n+dvDH2G2XsJiGoRDwt9vL+G5tEE+mb+t6FraCjGvLNpJRUwChYNHa8eS/dwzbrjTggSuRMVw+Ilx\nLCndxy+nrdLasOYcr/Lqtx2oKM+f9g3T3XKwWWfyGHrjMmsnLQQhtL6k/kyG2cSP5o4k1JX8vCVl\nDnahY1lxa+zz5OUGfrEreFLeFK14hjySa84x9vIZZ3qnR0Wvr0iYPCpll29qj3PSshnUcI1p509o\nx+wO+PGWooCf+8r4Qzsp2hC9MzSy+thWLH0HsmvkFCBwGqBkmE1seDWTSUNNXI2MZv+wCSwt2atJ\n3boJhC87xcPhadhXnfiCuy4pcDf6ALgHPNNT9w2fRF14lJfzA3C29nbAOz/K+OPdBOV7U4eSEN2L\nF25VknHpFB9OfFFbDr4y7YmAWQpmmE1s+sFUssfEkpM2jwG3rcw6l+81RvhJGqui8zBFhCJwbvQu\nOHmIz0ZP52ZYby3J4a0FKX5/D6yanECSS4W0RR/C9uRZZJ/JxdToHesP9Hz/TjH+QoinhRAVQoiz\nQojXOzj/x0KIWiHECdfPK53xuZ1JhtmEKSKU9w5XcvbqLZ7dvR5rWCRbUuYCzl9UZHhIz15kF/Da\nzGEcHDaBmt79WF68y+vcwL7hPXRVip6gwGLlrc/KkMALJw/Ru/kO68c+TUSojlWTE/wuyeF+hHpo\nPmxKyyLUYWPByUNeYwJhlXM/Htv4CyH0wK+BZ4AxwEohxJgOhuZIKce6fj583M/tCtwqnk9cv0jW\nmTz+e9yz3Al1SjkYAlT6IMNs4pWnRrAlZQ6zz+UT23BNO3fReoeVHwT20lfRiqeE86qinZT3H0rh\noNE0NjvYnH+hh6+uc/GUfKjoP5SigSOcBY8eDQCS4qJ64tK6jc7w/CcBZ6WUlVLKZmAjsKAT3rdb\nKbBYKal2Lvu+f3w7LXo9H2XM187PHNmf3Mq6gDSErz+bhP6V76OXDpaU7PM612xzsLWwuoeuTNEd\nFFisvLGthEMVV3FISLlylrQrZ9mQPk/b6PUnGZMHYdXkBN5dmMr0ETFkj4llc1oWSbXnSa45p40p\nqq4PaOenM4z/YMDTLah2HWvLYiFEsRBiixBiSCd8bqdRYLHy1qdlOIB+jfUsKd3H1uTZXOvVGts8\nUHE1oEWfqvoN4g/mNJYX70ZI7wYebdP/FIGDu6hrXV4Vx8475/WqEzu5YzCyPaV1ozdEH3ghkFWT\nE/jv70/mqVED2JE0g7uG0HYVv4Hs/HTXhu+nwFApZRqwB/hdR4OEEKuFEPlCiPza2touv6gCi5Wf\nbith5Qe5FLu8/pcLPyPM1syHkxZ6jbXbZUCLPl1raCInbR4J9TVMsRRrx0P1Iqg1zwOdtt26ejU1\n8kK5a6PX2Bu9TpA1JpYNq6f4/UbvvbA2NnMzrDdfjJzKi2UHMbZ4d/cqvVgfkA5fZxj/i4CnJx/v\nOqYhpayTUrp/ox8CGR29kZRyrZRygpRyQv/+/Tvh0u6Nu7BrQ14VzTZnv8+wlrt8p/D37Bk+iXPR\n3osTvQ4tFTTQPCBwdjTaNXIKN8J6a1LPOgHfe/KJgA13KZzZXJ6qHy+UH9Y2egGklIwd0jdgDT+4\nOpbpnBu/fZpuM+9Mrtf54ur6gFzxd4bxPw6MEEI8IYQIBVYAOzwHCCE8dV9fAMo74XMfC7fH4w5p\nCGBx6X6i79zkg0mLvMbqcG4Q/VX2qIDtdLV4fDz2UGdO97zTX9H3zk2khA+PfBPQ4a5gpsBi5ePC\nandYH3Dm9pf3H8rXg0YBzvsiEJ2dtuh0OnITUrnQJ5albXL+JYEpevjYxl9KaQN+COzCadQ3SSnL\nhBBvCSFecA37SyFEmRCiCPhL4I8f93MfF8/CLoNeENfbwPePb+dE3AiOxScTEarDoBdab95F4+N9\nsltXZ+HU/knlq5kLMNptLCw7gABsDme4KxAnfzDjufK1u6I+KVfOklpzzqtH75yk2ICd825yK+uw\n2R1IoWNz6lyetBQRX1/jNUYnAm/Po1Ni/lLKz6WUI6WUw6SU77iOrZFS7nA9/omUMllKmS6lnCWl\nPNUZn/s4uAu7VkxKQAekFH5JovUSH0xcBELQ2OwAKVkxKSFgvX1P3Dnee0MGciJuJCuKduHwSHtz\nSGcBkCIwaLvyBafX3xhiZLurojdEL3ht5rCeucBuxNMR3JrmrOtZ3CbrLZAKPN0EdYVvhtnEoL7h\ntNglq49t5UKfWHaOmqqdtztgUN/wgPujd4SnMdiYls2oa1WMu1ThNcbfNV0UrTjj3K3xnt5NjSw4\neYhPR8+gwdgLAHGvFwcYno5gWOJQjgwdy9KSvVrWW6AWeAat8XdLOJ+paWBcdTkTLpbz4cQXsbuk\nHCB44p3Q6v3ogE+TZnA7JEzr8uVGaf0HGB7B/hfKD9Gr5S4bXBu9AHZHYOX2fxsfF1ZzrvY2m9Ky\niL95Vct6EzqBKSI04Bq8B6Wkc4HFysq1R2mxSyTwH8e3cSOsN5tT53qN0+uDxfdp9X5yK+vIq6zj\n06QZPF9+mLdnv8JtYwQAH3xZSVbywKBYCQU67jg3AFKy6sROTg54ghNxIwGnVxiomW0d4bny3TMi\nkxthvVlevIevho7F7pC8+WkZSInNIQk16AIiFByUnv/HhdU0uwy/2XqJeaeP8j/jnqUx1FvLxhFk\nnk+G2cQNz8hsAAAgAElEQVSfzxrOkH4R5KRl06vlLvNPfamdt0t4S/X5DQg8O7ulXTlDSs05TbpZ\nAE+OiAkIA/eguMNgArCHhrIj+SmePv0VUXdvAc76npYAq/UJSuPv6c9///gntOj1/G58q5SDp7xz\nsHg+bgosVjbnX+DrQaOoiEnw6vIFzpL3Ze99xfq8qh66QkVnofXoPbGTxhAjO8bMRCfAGKLjR3NH\nBo3h13D3LBCCnNQsjPYWXnCJvel0ghBX9l+g2IWgDPu4OxGZGutZWrKX7WNmUdu7HwDmfhH83+Vj\nya2sC8pOVrmVddgcEoQgJ20ea/Z/wKja81T0H6qNsUtY80kpowZGBt3vJxAosFj57m/yAOdG7wvl\nh9mRNJPMcYmkD+kbvPPe7gz7tNglZbHDKBuQyLKSPfzP+OdASt5ckIq1sTlgfj9B6fm7s1a+8/Xn\nhNua+MBDysHmcGjhj0D4Az8s2savgG3JT9GkN7C8jfcP4JDBFRILFNbnVbH0P77iVrOzcc+LJw86\nN3rT5xETaQwYw/aweKZ7uvfBN6VlkXblLElXK5E47UYg2YWgNP6ZidFE2Jt4ufAz9g2byNmYVnnX\n8BD9fV4Z+GSYTayZn4xOCKwRfdg9YgoLyw5gtHmneYYGyNI3mCiwWPn7T0rRlHykZNWJLygbkEhR\n3Eg2HqsK2kpuz4ZOr01PBGD7GKfzs7R4LwIounCDn24rCZjfT9AZ/wKLla2F1SwsPUBMY72X1w/w\nvWmJPXRlvoO1sRm7q5HxhvR5mO42MO/0Ue28ANbMTw4YDyhY2FpYrf1dAdIvn2bM1W+0it5A2sx8\nFNwr/tefTSJ7TCz14ZGa86NvaWH3yRrW51Wxcu3RgPgCCCrjX2Cxsuz9r1ife57vHdtG8cDh5A5J\npX/vUNLj+/DuwtSA6VT0OHhmghw1p1HVJ9ary5cEyi7V3+PVCl/laoO3WuWqEzu5HRLGJ2OeAgJb\nuPBheW3mMAw6waa0LEx3G5h7Nk871xwgvQ2Cyvj//Ity7A6Ye/YYw65f5IOJC0EIam81U375JqNc\nfT2DHafOTwoAUujIScvmSUsxCdbL2pjaNoZE4bu4i5M8s9wim27z/KnD7EiawS1XHcecpNigSu+8\nHxlmE69Me4I/mNO5GNm/XYP3QJA6CSrjf7bWmbP76rGtVEcN4PPR07Rzgdap6HFZNTmB7DGxAGxJ\nnYNd6LwaXeyvuBoQS99Axy3g9s+7Kzh4uhaDq3BxQdlBIlqanLn9Lu602JXh9yAyPASp07MldQ4z\nvikk7mZrj5E3A6DeJaiMf6QxhHEXTzGp+iS/mbjAS8pBJ4JHyuFBeW3mMEINOmoiYziQmMHSkr3o\nHc4sEbv6svQL3JWrDgl2u4OYXqEgJS+d+ILS2GGUDByujX0mJe4+7xR8ZCZGE6IXbEmdiw7JotL9\n2rlmm4P3D527z6t9n6Ax/uvzqrBcb+SV49uoN/ZiU2qWdk4Ab7+YqryeNmSYTWx4NZPpI2LYlD6P\n2FvXeaoyH3DG/RvutASU1kkg4pnCiIArN5sYe/k0SbXntY3ewX3D1H5XB2SYTYwd0pcLfQfyB3Ma\ny0r2eLU43XfKv1e/QWH8CyxW1h4+R4L1Mk+fPsq6cc9oejXgDHGoid8xGWYTP5o7kiMjJ3G1l4kV\nRa0bv2u/rFSNXnwczxTGPuHOOPWqE19wOySMHUkznc8nm9X8vwfXG51ihptSszDfuMLkC6XaOYdD\n8su9p/127ge88V+fV8Xy949yvq6R7+dvxy50/Hb8815j3BW/io7JMJt4MimOzalzmX0un9iGa4BT\n4z/Y0wP9gQyziczEaG40Njs3esu/5JMxM7lljECvwp33JTHGKW+9c+RUbhp7eXX5ksCXZ66x/P2j\nfil3EtDGv8BiZc0npdgckr53brKseC/bk5/iaqT3ZFc69d/OgEgjm9Ky0EsHiz1inwB6XeB1OQo0\nthZW45DwYtkBwm2tG70ZZpMKd96H12YOw6AXNIUY2ZE0g2crviKy6bbXGJtDsuaTUr9bAQS08c+t\nrNOKWv7ILeUw0buoS230PhiLxsdTHT2IrxLSWF682yv2+dSoAcqA+DAFFisbj1dp0s0lscModW30\nFlZZ/c5odScZZhM5q6eQHt+HTWlZhNuamF/+Zbtx/ih3EtDGPzMxGr1eYLQ1892CzziQmMGZ/mav\nMcHQo7QzyDCbeHtBKlvGzcN84wqZVSU9fUmKB6DAYuXHHxdjd8C4SxWujd7W9E67A78zWt1NhtlE\nbFQYxQNHcCrG3C7nH8Dgh6vfgDb+GWYTyXFRvFh2gP6NN1g7aVG7MT8Igh6lnYG7x+/vh0/hRlhv\nL6nnfaeu+mXMM9Bx5/ifveqsb1l1Yie3QsPZkTRDGxOi9z+j1d0UWKwcrLgKQrApLZtxlysYUWvR\nzgtg6YQhfudEBrTxB5gy1MTqY9sojR3G0YQ0r3M/mJHod3+wnsKdL95kCGVb8iyePv0H+t65CTjb\n/f3d9hL1BeBj5FbW0dTiDM9F3b3F/FNfsiNpJreNESQNjGTV5AQ2rJ6i7oFvQZM5B7YnP0WzzsDS\nkr3aeYNeUHqx3u/mf0Ab//V5VZz9bQ7DrlezdtJCr56lEJhNmbsKz3zxreOexmi38WLZQe28w6Xx\nr+LHvsOZmgbcMm7ujd51Y58mRC/4h4WpvLtQ1bY8CJ79ra9H9GHv8EksKttPiL2FMXGRtNglRdX1\n/HSbfzlAAWv8f/Z5OT/dVsKrx7ZRHdWfz0dN8zqvUtweDs988ZWvzKcoboQz51+2qkT646ZXoFJg\nsfJJ0SXnE9dGb/HA4ZQNHB507UkfF/fc/+t5o3h3YSpn5y8jprGe1Q3lXLnprXGVc1wZ/x5lfV4V\n7x2uJP1SBZMvlPJfExZg07c2LVMVvY+GW/K29FI9G9LmMfqahbGXT2vnDXqlCOkrbC2s1r6Xx188\nxehrFtanP41AKXc+Cu65v2pyAk/+2Spqekczbu82rt/2ThMfEBXWQ1f48ASk8f+i1Kk++erx7dw0\n9mJjWrZ2TifgHVXK/lhca2ji06QZ3A4JY7mr4leA1ypA0bN4/iVWFTk3ej8dM5NVkxOUcudjkltV\nz5aU2cyqLGBAQ+sKSi/8K4EkII3/MylxDLlxhWcq/sD6sU97STkY9Dol3fyYxEQauW2M4LPR03mh\n/DC9mhq13qd/vemEX8U9A5XF4+PR0brR+8mYmdwODQdQhv8xMUWEagWPi8oOAE7nZ8WkBL/63Qak\n8V81OYGfVOzELnT8V4a3lIPdrqQIHpfF4+PRC8hJz6ZXy12eO3UEcHqb5+sa/W7jK9AosFh5/9A5\nHMDCsgOE2Zq13P4zNQ09e3EBgLWxGYtpEHlDUlhaskdb8dY2NPlVwkPgGf9162DIEJ45uIUWvYHJ\nHsVIKt7ZOWSYTbz9YionBo/mdHQCKzy6fLlxh94U3UuBxcrKD3LZfbJG69F7Im4EZbHOcIRbqEzx\n6GQmRqMTTrG3YdcvMuHiSSSw+2QNKz/wH5HDwDL+69bB6tVQXY0AerXc5Wc7f8ULrqVZ1hjVqaiz\nWDU5gbi+4eSkZzP+UgUja897nVfa8D2Dux4DIONiOaOuVbEh/Wnt/BMuoTLFo5NhNrF6eiKfj3qS\nW6HhXhW/zTYHWwure/DqHpzAMv5vvAGNjV6HImxN/O3hjxBA+pC+yvB3EgUWK5fr77I1eRZNeoNX\nxa+7gEjR/Xi2F1xVtJOG0HA+9ajo9acNSV/m9WeTEL178eno6Tx36gi9mlrtzsbjVX7h/QeW8a/q\nOM486OY1JeDWyeRW1iElWCP6sHvEFBaWHcBoc6a9DekX8S2vVnQVpZfqAedG73OnjrA9eRaNro3e\n7DFKx6qzKLBYaWqxszkti14td3m24oh2zu6Aj/3A+w8s45/Qsbd5KSoGu4Q9ZVe6+YICl8zEaIwh\nzumzMX0eprsNZJ8+CsCF642qw1cPca3BWXS0uHQfYbZmr5DPa8rr7zRyK+twSCgcNJqz/eJZVrzX\n63zO8Qs+n/QQWMb/nXcgwtvrbDQY+acZLwOwUxn/TiPDbGLN/GQMOsFX5jQu9IllebEz9FN+pYFf\n7Krwq82vQKF/pBGkZOWJXZyIG8nJ2ETAVYeh6DQyE6PR6wQIQU5aNhMvniSxrtXbt/uBxn9gGf+X\nXoK1a8FsRgpBdVR/Xn/6h+xIngXA08kDe/gCAwtrYzMOKZFCR05aFtMsRQy50foFGwhNrv2N5EF9\nmHDxJCPrqljv4fWDkm7uTDLMJt5akIJeJ9iWMgub0LGsxFvq2e7jMhqdYvyFEE8LISqEEGeFEK93\ncN4ohMhxnc8TQgztjM/tkJdegvPnEQ4HNcUVXH5uEQOjjPxgRiKvP5vUZR8bjLgFrwSwJWUudqFr\np3W+92QNb2wr8WkPKGBYt45nns9k87of40Dg8DhlUNLNnc6ogZHoBVzrZWL/8EksLt2HwW7Tzku8\nN+B9jcc2/kIIPfBr4BlgDLBSCDGmzbDvA1Yp5XDg/wE/f9zPvR8FFiu/PnCWiisNPDVqAL9+KUMZ\n/i7ALXi1anICV6NiOJiYwdKSPegddm2MA6fWkmry3sW40pxNtZcRgA7JW3vf19Kc/VFv3tfJrayj\nxe4s8NqUmkX/2zeY+U2Bdl7g2y1iDd8+5FuZBJyVUlYCCCE2AguAkx5jFgBvuh5vAX4lhBBSdr4Y\nTIHFyoq1R7U/CkBYiE7l93cR7h6wkUYDOaezWbvtHZ6qzGff8MnaGElrk3f1N+gi7pPmvDNtNovH\nx/fQhQUu7ri/zSE5mJhBba++LCveo839gPf8gcHABY/n1a5jHY6RUtqAeqBL1qDvHzrnZfgBmluU\npENX09BkY/+widT26uuV8w+qsro7kPdJc1Y9lruGDLOJV6Y9AYBNb+Dj5NnMPnecmNutK1xf9vx9\nasNXCLFaCJEvhMivra19pPeouXm33TGdH/bX9Dckzhtgc+pcZp07rqkdGvSCtPg+rJmfrAxQF7E+\nr4orffp3eO5SVIwzA0jRJUSGh6BzpVJtTssixGHnRVeoDQLf878IDPF4Hu861uEYIYQB6AO0c8Wl\nlGullBOklBP69+94Mn8byyd65/rrBLy1IEUZni6kwGKl1pVfvik1C4N0sKR0HwBSSkou1vPWZ2Uq\n5t8FrM+r4qfbSvjHad/BJrxv50aDkf8z82UWqZBPl5GZGI3BZf3PRQ8hf3CSM+ffFdFe80mpz+b7\nd4bxPw6MEEI8IYQIBVYAO9qM2QF81/V4CbC/K+L94NSceXdhKunxfcgeE8vmH0xVUgNdiLtJ+N6T\nNQCc7zeYowmpLCveg5AO7A5ni8emFv/RPPEn3AJ6+4dPxi4Et0PCcOBMc377hR/xnX97Qzk+XUiG\n2cTSCa2+76bULEbWVTHuUgUANh/O939s4++K4f8Q2AWUA5uklGVCiLeEEC+4hv0GiBZCnAX+CmiX\nDtqZrJqcwCc/nMbalyeoid/FuIXEPL/JN6ZlM/TGZTKrSrVjEticf8EnbwJ/JjkuCoDnyw9jdNh5\nacU7JP74U6b96X+h+6OX1PzvBhaNj9e8/9+PnkZjiNEp9ezCV/P9OyXmL6X8XEo5Uko5TEr5juvY\nGinlDtfju1LKpVLK4VLKSe7MIIX/49nY3eCaTTtHTqXe2IvlbaSebT56E/grBRYrvz16HoDlxbs4\nFWPmRNxIAEINOhXu6SY8N35vGyP4/ajpPF9+mPBm5/6jr2b9+NSGb2exPq+K7/wmz2djbYGEZ2P3\nnNemMmNEDE0hRrYlz+KZiq/oc6e1eYhD+uZN4K+4V11JVysZe/kMG9PngRBkj4llw6sqtbk7iQwP\n0R5vSptLZPMdnjn9B+3YwYqrPXFZ9yXgjL97A+zLM9dUR6luwt3cOsNs4rSrU1ROejZGewsLPTIf\nAMpcqpOKxyczMRohBMuLdtOkD2GbS8YkJtKoDH8345lNeDw+mUrTIK9q973lNT4X8gw449+2g5Tq\nKNW9uOWcywckciJuBCuKdnk1dj+t2gg+Np4V7OZwZ6vGnSOnUh/u7E3tVvZUdB8ZZhNZY2KdT4Rg\nc1oWmRdKMVsvAc5Vr6+FPAPO+LftIKU6SnUvC8e1xplz0uYx+pqF9MuntWNW1UbwsXBnV/1iVwU/\n3VZC6rH99Gm6zcb0bG1MjMrr7xF+MHMYoXqBALYmz8YudCwtaZV6brjjW3M/4Iy/O9Vz+ogY3l2Y\nqtI8u5ECi5W3PivT5IM/TZpBY4hRk3oGCNEpceHHIbeyjqaWVsm2FcW7sfQdSG5CqnYsZVCfnri0\noCfDbGLD6in8zbxRzJ07jkOJGSwp2YvOpXX13uFKfvZ5eQ9fZSsBZ/zB+QXw39+frAx/N9M27fOW\nMYLPRk/nhfLDRDTfAaCipsHnYp/+hDPO73w89PpFplSVkJOWjfQo8PJlSYFAx73/tWh8PNvGZTPw\n1nVmfPO1dv79w5U+M/8D0vgreobMxGgMeu8ptTFtHr2b7/DcqS8BZ+xTFXs9OhlmE9OGxwCwrGQP\nNqFjS8oc7bxqV+obZJhN/Mm7f0FdeJRXzr/Ed2L/yvgrOo0Ms4klGfFeXaMKB4/mTPQQL7E3Vez1\neNTfacFgt7G0ZC8Hhk3kamSrsV89PVFl+vQwBRYrb2wr4R/3nWNb8iyyzuRhamzNcvOVL2dl/BWd\nyuLx8VpvXwCEYGNaNhmXTjGi1gL4bsWjv2A06Jh97jj9b9/w2ugFqLx2u4euSgFOw7/yg1zW5VVx\n/LyVTWlZhDpsvHjyIOBb7TSV8Vd0Ku6ir+kjYrRj21Jm06wzsMK18atXKquPTIHFSmGVleXFu6np\n3Y+DiRO8znekaqvoPnIr62ixtW7In+4/lBNxI505/66UZ19xfJTxV3Q6GWYTP5o7klCX3sP1iD7s\nHjmFRaX7MdpaVFepR6TAYuWtT8uIuXGNpyoL2Jw6F7tO7zWmraqtonvJTIwmxOBtVjelZZFUe57U\nK2eRwKGKqz4R9lTGX9ElZJhNvPl8sqZ1vjEtG9PdBrLOHOWL0is+lfLmD7jz+4ur61lSuhe9dLAp\nNUs7L4AfzEhUGW49TIbZxIZXM3lpcgIThzodnE+TZnDXEKo1eD923sry97/q8S8AZfwVXUbppXoc\nrrzPPwxNpzpqAMuLdnP9djPvHa7k5d/k9ewF+hFafr90sLx4D0fM6VSZWgsYBd76MoqeI8Ns4p2F\nqTw1agAADcZefD7qSRacPISxxVl9bXPAxz2c9aaMv6LL8NzckkJHTloW0y0niL9xBYDDZ64p7aUH\nxBQRigSmWooZUl9D/pxF/GBGIgadQAeEhqg2mb5GZmI0eneXr9QsoppuM+/MUe182cX6HvX+lfFX\ndBmLxsdrcX+ALalzsQudl+CV0l56MEpdgngri3ZhDYvk/X6pZCUPJOe1Kfz1vFGse0WpePoiOlfc\nMzchBUvfgV5zv6i6nuXvH+0xB0gZf0WX0TbufzmqP4eeGM/Skr3oXSXv7mYkivsjAFNjPdmnj7It\neRZ39KHkVtZ5KaoqfIvcyjrsrrinFDo2p85lmqWI+PoabUxPdvpSxl/RpbSVGshJzybuVh0zKwsA\naGiy9cRl+R3Jg/qwqOwAoQ4bOa7c/jNKIdWncTc6cvNxyhwcCJYW7/Ua55A9U/eijL+iS3HfAO74\n/75hk6jt1VfL+c/Jv8BPt5X0eOaDr2O93cTyot18HTeKiv5DAdh+4pLaM/Fh3DUvw/v3Apwr3y+f\nGOcl9gbOrms9sV+jjL+iS3HfAH8zbxRj4/tg0xvYkjKX2WeP0f/WdWx2yfq8KlauPaq+AO7DiMpS\nRtZVsSF9ntdxtWfi22SYTUzyMOybUrMY3FDLVEsxAIP7hrFmfnKPhO2U8Vd0Oe64tFvtMyctC4N0\nsKR0nzam2S57PPXNF3HrxNT/239wKzScz5Kme51Xeya+z+Lx8YS60n72jMjEGhap5fxfvHGXtz4r\nUzF/RWATGxUGwPl+g8kdksLyot1eXb58SffEF3AXdu04fIrnyg/z6ejpNIaGe425qfZMfJ4Ms4k3\nX0hBJ6DZEML25KeYd/qo1t+6xeZQMX9FYPPazGG4FZ83ps9j6I3LZF4oAcCgFywaH3+fVwcfWwur\nudviYH75YSJamshpE/IB9YXpL1gbm7WCx01pWRjtLbxQfgjoOa0rZfwV3UaG2cSm16Yy2BTOFyOn\nctPYy+n9A2+9kKLSFT0osFjJyb8AwPLiXZyKMXMibmS7ccmqa5dfkJkYraU8lw9IpCR2GMtdOf8t\ndsmesivdfk3K+Cu6nSv1d2kKMbIt+SmerfgDfe40aEVMCidbC6ux2SVJVysZe/mMM71TePv5OlTX\nLn8hw2xi9fRE7fmmtCxSas6RXHMOibPFY3dnbinjr+hWcivrcLjWvzlp8zDaW3jx5EG+tlj59YGz\nKuMHp9e/2eX1LyveQ5M+hK3Js9uNU5IO/sXrzybx7sJU+kWEsCNpJk36EJZ4NHjv7swtZfwV3Ypb\nowbgZGwiRQNHsKJoF+WXb/KLXRW89GFu0H8B5FbW0WKXGFuaWFS6n50jp1IfHuk1ZviA3krSwQ9Z\nNTmBCUP7UR8eya6RU1hYdgCjzbl6a2qxd+vcV8Zf0a10VPGbVHuetCtngJ7LfPAl3PHhp09/RZ+m\n2+26del18PPFacrw+yEFFqvWcCcnLZu+d28x94xT3fb4eWu3Oj/K+Cu6lczEaAy61tj1jqSZNIYY\ntR6/ep3g0o07Qe/9C51gRfFuLH0HkpuQ6nXu7QWpyvD7Ie7U3ZKLzv2tr8xpVEf113L+Jd3r/Cjj\nr+hWMswm3lqQou1d3jJG8PtR03mh/BDPJ0YiwVnx+0Hwhn8+LqxmSG01U6pKyEnLRgrv21Rt8von\nuZV1NNscOKQzRVcKHVtS5zL9m68ZdPMqAHp99+3jKOOv6HZGDYxE7+H9b0rPpnfzHRL2f06LXSKB\nZpuDrUFY8VtgsbKloJqlJXuwCx1bUuZ4ndcJ1Cavn+LWudILCHH9uyV1Ljoki0tc1e4eRY9djTL+\nim7HM+MHIH9wEmf7xTPryA6vcVcbmrr70nqc3Mo6aG5mSck+9g+bwNVIb0OfOriPCvn4KW6dq7/K\nHsWGV52b9dV9YjliTmdpyV6EdNBi7z6FT2X8Fd2OpwekFyCFYGN6NhMuljOi1qKNGxBp7MGr7Bka\n7rQw81w+A25b2dhBRe8U5fX7NW6dK0ALa25KyyKhvobMqlIkzoy47kAZf0W34+kBub3YrSlzaNYZ\nWO6SehaCoJN7KLBY+fDIN6wo2kVN734cTJzgdV716Q0ccivrsLsWv7tGTOGmsRdLXRu/Byqudss1\nKOOv6BHaKn1ej+jD7hGZzoYlthZem54YdOGN9w+dI+ZGLU9VFrA5dS52nd7rvFEVdQUMmYnR2r5X\nU4iR7WOc1e6RTbfZc7KGn31e3uXXoIy/oscosFgprLqhPc9Jz6bfnZv8aX0JWckDe/DKup/1eVXs\nPlnDktK96KWDTalZXuezx8Sqoq4AQ9C677UpLYswWzPPlx8Gukfu4bGMvxCinxBijxDijOvfDmem\nEMIuhDjh+tnR0RhF8OHZ4xTgyNCxVEcNIGPvtqCr9M05XoWQDpYX7+EP5jSqTHHaOaNBx9qXJyjD\nH0DkVtbhMfUpjR1Gef+hLHOFPaHr5R4e1/N/HdgnpRwB7HM974g7Usqxrp8XHvMzFQFC24IvKXRs\nSstixvmviam9FDSVvgUWK2WXbzLVUsyQ+hpy0rw3eptsjqD6IgwG2rY3RQg2pWUx9vIZRtWeB+CZ\nlLh7vbxTeFzjvwD4nevx74AXH/P9FEFE24IvgM2pc51Nrkv2cqamgV8fOMv6vKqAFn3LrazDZpes\nKNqFNcyp+dLRGEXg4E56WDU5Qevy9cmYp2jRG/jRhSO8uzCVVZMTuvQaDI/5+lgppXttcgWIvce4\nMCFEPmADfial3N7RICHEamA1QEJC1/7HFb7BqIGR6AXYXEvgy1H9OZQ4nmXFe5j+9UocOr1WEWkM\n0QVk3NsUEYqpsZ7sM0dZN/ZZmgztU/26K/1P0X1kmE1kmE0kD+pDzvEqYqNiuZo4mnl7cxB7cmga\nPBjjz38GL73UJZ//rZ6/EGKvEKK0g58FnuOklBK4V3maWUo5AVgF/FIIMayjQVLKtVLKCVLKCf37\n93/Y/4vCD3HG/b2PbUybR9ytOqZXFmpx0e7WPekuCixW3txRyqKyAxjtNqdufxsEStIhUHH//Yuq\n6wnbvJGYygp0UiKQGC9WY3/1VVi3rks++1uNv5RyrpQypYOfT4AaIUQcgOvfDhNUpZQXXf9WAgeB\ncZ32P1D4NZ4Sz272DZ9EbURfVhTt8jouAlDa4OPCapptDpYX7ebruFFU9B/aboySdAhc3PLdAH97\n+COM9hav8/o7d+CNN7rksx835r8D+K7r8XeBT9oOEEKYhBBG1+MY4Eng5GN+riJAsDY2o2vTiNam\nN/Bx6hzmnD1G/1utcX4ZYB1r3U1bxl86xci6qnbSzW7mJMUGXKhL4SQzMZoQV8x/0M1rHQ+q6pqU\nz8c1/j8DsoQQZ4C5rucIISYIIT50jUkC8oUQRcABnDF/ZfwVgLfUg6dpz0nLxiAdLC7dpx2Tju7T\nPekOthZW02KXLC/aze2QMD4bPV07p9c5fx+hesFrMzuMkioChKUThpA1JhZr9D22TLto//OxNnyl\nlHXAnA6O5wOvuB5/BaS2HaNQQGvWQ25lHZmJ0fz8i3KOnbfyTb/B5A1JYXnxLt6bvBiEwGAIrApX\nCfRuauT5U4f5JGkmt40R2rnZo2MZO6QvmYnRyusPUNz6/s02B6EGHVuX/Bkv/eYdImytgob28HD0\n77zTJZ+vKnwVPY5b6iHDbOLHzyQRFuLMf96Yls0T1stMvlDqHNiNcrfdQcqgPjxffpiIliZy2oi4\nDUR/M4sAABFLSURBVIg0ar8TRWDiqe/fYnPwWcosXn/6h1RH9ceBoDqqP/v+1z90WbbP46Z6KhSd\ninslsLWwmm22qdzc+z7Li3eTl5BKs13ycWF1wBhEa2Mzy4t3cSrGzIm4kdpxvU4EnahdMOIOebbY\nHIQYdExJjGZtyix2JM8CnCG/Davb13x0FsrzV/gcGWYTg/qG02gI0wSvou7eAmBLQbXfF3sVWKz8\n+sBZLh3KZezlM870To9Kt5EDegfMF5zi3niq266Zn8xvj55HSqfMefaYWDasntKl80AZf4VP4o7t\nb0yfR5itmRfLDgD+n+vvjvP+YlcFwz/bRJPewDaXp+fG3eBbEfi4Q57WxmaabQ4t7Tl9SN8udwCU\n8Vf4JBVXGgA4GZtI8cDhrCzaBVJ2a7OLrmBrYTV3WxwYW5pYVLqfXSOnciM8ymvM8AG9e+jqFD1F\n2xaP3ZHYoIy/wucosFhZ80mp9jwnLZuk2vOkXjkLwMFuanbR2RRYrOQcd+ZszztzlD5Nt9nQQbeu\nHz+T1N2XpuhhPENA3SVhooy/wudoK/W8Y8xM7hiMrCh2VvzuLa/xy7h/bmUdNpeUxcqiXVj6DiQ3\nwTsLetJQk4r3BymeWW/dgTL+Cp8jMzEaY0jr1Gww9uL3o6fxwslDhDffxSHhvUPn/E7p0x2uGnr9\nIlOqSshJy0YK71tweGxkT1yaIghRxl/hc7iXwC95SNpuTM8msvkOz1UcAWD/qav88+4Kv2j6UmCx\nsvqjfP5132kAlpXswS50bEnxro8M1QsWqxRPRTehjL/CJ3Gne7rJHzyGc/3iNbE3u0PikNDs49k/\nBRYrK9YeZffJGq7cbMJgt7GkZB/7h03gamTrpl56fJ8uT+1TKDxRxl/hszibXLueCMHGtGwmXCxn\n+LVWoSuH9O3sH0/VRoBZlfkMuG1t161rQFSYMvyKbkUZf4XPkmE2sem1qZj7OTVvtqbMpllnYLlH\nn1OAskv1PXF5D0TbL6blRbuo6d2PA8MmeB0PLL1ShT+gjL/Cp8kwmzRVy7pefdkzYjKLSvcTamvV\nPd+QV8X6vK6RvX1U3FW8BzzSUgfevMasygI2p87FrtN7jY+JNHb3JSqCHKXto/B53Lnx4Mz5f67i\nD8w9m8fno6cB4ADe2FZC2aV6Fo2P7/HwSYHFysq1R73CPQBLSveilw42pWZ5HdcL1EavottRnr/C\n54mNCtMeHxk6luqo/u26fElgXV6VT2T/fFxYTbPdWY3sNv9COlhevIc/mNOoMsV5jX91emKPf2Ep\ngg9l/BU+z2szh2ndjhw6PZtTs5h2/gTx9TXtxvqC9k9H8fuplmKG1Ne02+gVQGR4SLdcl0LhiTL+\nCp8nw2xi4+oppMf3AWBz2lwAlhbvbTe2u3RR7sei8fHo2/SmXFG0C2tYJLtGekv0+rtWkcJ/UcZf\n4RdkmE1a+OdS1AAOPzGepSV70Dns2hgBrJmf3OMhlAyziVenPaE9NzXWk33mKNuSZ9Fk8Db0Aqeu\nv0LR3Sjjr/ALCixW9ntkzmxIn8eghmvM+KZQOyaB/7PrFK9+lN8jcf8Ci5VXP8pnwa+OcO7abe34\norIDGO02p26/BwIwhvT8SkURnKhsH4VfkFtZh80je2b/8InURvRlRfFuDg6bqB2/3tjCnpM17C+v\nYdMPpnbbKqDAYmX5+19pwm3gqj2QzgbtX8eNoqL/UMBp9PV6wbIJQ1jsA9lJiuBEef4Kv8BZ7dsa\nR2/Rh/BxymzmnD1G/1vtvXy7dGbduPPtu3ol4KnY6cn4S6cYWVfFRg+vXwLSIRncN1wZfkWPoYy/\nwi/IMJt4e0GK9gWg1wmKsxcT4rCzqGxfh6+51tDESx/mdosAXGZiNIYO7qblRbu5HRLGZ6OnO69b\n0K0NOxSKe6GMv8JvWDU5gbcXpGDQCRwOyRctfciLT2Z50W6Q3gVVQsANV2s8h+z6FNAMs4m3FqQS\nEdpaudu7qZHnTx3m06QZ3DZGoBPw9oup3dqwQ6G4F8r4K/wKa2MzDtlaQJWTnk2i9RKTqsu8xkkJ\nx85bkdIZY5dAw52WDt7xwblfCKnAYuWtz8pobG7NPppffpiIliY2urp1rZ6eyKrJCd3asEOhuBdq\nw1fhVzjDK0KTTvh81JO8uXctK4p2cWxISrvx7vWAlPDe4UrAWVSVmRj9UAbY3Xi92eYg1KDTPPcC\ni5Xcyjou3rhDU4t30H9F8S5OxZg5ETcSgIYm2yP8jxWKrkEZf4X/IQTSZdbvhoSxfcxTLCvZw5v/\nf3t3HxxVdcZx/PvsLokEJWwAo9lgDIxawGI1qSK2+AIqWivWQavQYp1WnY5W69RxqH+otcNUbWvr\ntI4Vqa22KHYUWsY3sIDVzghCUN6l2KTE3QQTSSAImpfdp3/c3ZBklwTY7N5k7/OZyWT3srv3uUPm\nd8+ee+4502+j5bjeFz9/6u1qROgMcHAu1vZ1MlhTvaezC6m1PcaSDWGAzjl8/H7BJ86FZoDxDdV8\npX4nP5t2i9MHxaETkTEDgYW/GVScIZ/dW9gvnnUZc99/lZnb3uIv51zV6/sV51tAe0eMP/zrv6z+\nsIGYarfWfCqTx44k4PfR1hFDgcXv1bLzk/20xdO+I6qUFRUQ3vs50Zhy/aY3afUHWDrx4s7POLOk\nMK1jN6Y/WZ+/GVQmjx1JXsDX+YcrwNbicWwuHseNG5cnXfjtKTHaxu8TVn3YQMcRrghWURZkVsWh\nmTej8WsKXe1qOkg0puS3t3LtllUsP30Ke4cO7/z3h17Z6vqkc8YkWMvfDCqJ9X3XVO8hWJBH88E2\nggV5bKqfxZznHuG9J+Yy6sBe6oaP4tGpc1nWpeUNTst/2vhiGlq+YGO4+yIwh5tjJ9Gvf2ZJIf4u\nXTuHc/nOdylsPcDiSd3v6E2MOLKLvWYgsPA3g05FWTA5QDeNQ5+DEw84LevSlkYefuP3AN1OADGF\nlduTZwNVdVrmZ5x0AhVlQZ5fW8vrW+oZOSyPZR/UEQP8PigalkfjZ73PxXPDxhXUFhbzbtmkzm2C\nje03A4uFv8kN8+cnTaVc0NHKvW8/l9T6j8Vb7j6BE0/IZ3dLK4rTMl+yIcwjr29P6tIBiMboM/jL\nmuuYUruJR6fOJRDwE4spfr+PWRWlNpWDGVAs/E1O0NralPPoh1oaeWjFk9QUlVATDFFdFCJSeCJR\nnx9V2N3SeujFIry4rjblNA1H6tubVhAVH0u+PI3rKscQGjH0qIeVGpMNFv4mJ+wffTLDG+qStrf7\nAlyz7S2Gtx6aZbPNF2BX8GRqikJUB0uoKQo5P8EQjcNGdA7N7I3fJ0Rj3Tv/A9EOZm1eyapxlTSN\nGG0tfTOgWfibnPDpfQ8QuOdOCjoOteQPBvKZN+MOlk24iKLPWyhvijC2KUJ5c/x3U4QLq6vIjx66\n87clryB+MnC+KdQUOd8WaoIlHMgv6Hxdz+C/eutq7l+5kFGf7+PsyA5+fmAjFWVXZP7AjTlGon0M\njev1zSLXAQ8C44FzVXX9YV43A3gc8AMLVfXhvj67srJS169P+XHGpFT9+AKC8x+ksHH3YUf79OSL\nRSnZ/ynl8ZNBeXMd4/aEKW+uI7SvAV+XW7M+Ob6ImmBJ/GQQ6jxJTKrbwfwVT3Y78USHDsX/9NMw\nZ07GjteYVESkSlUr+3xdmuE/HogBTwH3pAp/EfED/wEuBcLAOuBGVd3W22db+Jtj9fzaWu5bujnt\nz8nvaOOU5nrGNkcob6pjbFOY8qY6ypsjjDp4aJioknrd3tZQKfnhj9Ouw5ijcaThn1a3j6puj++s\nt5edC3ykqtXx1y4GZgK9hr8xx2r2eacA8OK62qSx/EejNZDHztFl7BxdlvRvw7/4rPObwm9e+XXK\n9w+JRKja1Wz9/mZAysYdviGga/MnHN9mTMbMPu8Uzs/gmPqW445nY8kZ/H3ixUSGj075mvrhozI6\njbQx6egz/EXknyKyJcXPzP4uRkRuFZH1IrK+sbGxvz/eeEjVrmYWvFOdlX09OnUuBwP53bYdDOTz\n2CXfs5u6zIDVZ7ePqk5Pcx8RYEyX56Xxban2tQBYAE6ff5r7NR62pnpP0jQ/fh9oLD6SU5KHah6r\nZRMvJuATflG1mPy6CK0lIf59093MvvVm6/IxA1Y2hnquA04TkXKc0L8BmJ2F/RoPmzx2JPlDfLS2\nxxBxFlK5dOJJvLwhzNbIPjalcS0gwe8TbvlaubM+wA+nkF/2SwDygct6f6sxrksr/EXkW8DvgNHA\nqyLygapeLiIlOEM6r1TVDhG5A1iOM9TzGVXd2svHGpO2rhPAJe6wrdrVzJINYVrbY0lz6ydW+zoS\nXz01yOnFJ3Ct3cRlBrF0R/ssBZam2F4HXNnl+WvAa+nsy5ij1XMCuMSCLJ1DM6XPGaCT5AV8zLti\nvIW+GfRsPn/jGYm1APwC+UN8TB9f3Dk+/0jPAdFoZheCNyZbbHoH4xk9u4IA3tnZmLIbKBWf2LTM\nJndY+BtP6dkVtOgHk7l9UVX32T178Avc8vWxx7TwuzEDlYW/8bSKsiCTSkewe1vyAi8CTCot5P5v\nTkwK/MTqXnYyMIOVhb/xvNsuHMfqHQ20RxW/T5ylGmPKkIDvsME/Z+Ea2jpifS78bsxAZeFvPK+i\nLMjiW8/vdi2gt1Z9YtRQTG1dXjN4WfgbQ/K1gN7CPDFqqL0jZheAzaBl4W/MUUp1A5kxg42FvzHH\noOc3BWMGG7vJyxhjPMjC3xhjPMjC3xhjPMjC3xhjPMjC3xhjPMjC3xhjPEj0aCc0zxIRaQR2HePb\nRwGf9mM5g4EdszfYMXtDOsdcpqqj+3rRgA3/dIjIelWtdLuObLJj9gY7Zm/IxjFbt48xxniQhb8x\nxnhQrob/ArcLcIEdszfYMXtDxo85J/v8jTHG9C5XW/7GGGN6kXPhLyIzRGSHiHwkIvPcrifTRGSM\niKwWkW0islVE7nK7pmwREb+IvC8ir7hdSzaIyAgReUlEPhSR7SJyvts1ZZqI3B3/u94iIi+IyHFu\n19TfROQZEWkQkS1dthWJyJsisjP+u9+nkM2p8BcRP/AEcAUwAbhRRCa4W1XGdQA/UdUJwGTgdg8c\nc8JdwHa3i8iix4E3VPVLwFnk+LGLSAi4E6hU1TMBP3CDu1VlxJ+BGT22zQNWquppwMr4836VU+EP\nnAt8pKrVqtoGLAZmulxTRqlqvapuiD/ejxMIIXeryjwRKQW+ASx0u5ZsEJFCYCrwRwBVbVPVve5W\nlRUBYKiIBIACoM7levqdqr4NNPXYPBN4Nv74WeCa/t5vroV/CPi4y/MwHgjCBBE5FTgbWOtuJVnx\nW+BeIOZ2IVlSDjQCf4p3dS0UkWFuF5VJqhoBfgXUAvXAPlVd4W5VWVOsqvXxx7uB4v7eQa6Fv2eJ\nyPHAy8CPVbXF7XoySUSuAhpUtcrtWrIoAJwDPKmqZwMHyEBXwEAS7+eeiXPiKwGGich33K0q+9QZ\nktnvwzJzLfwjwJguz0vj23KaiAzBCf5FqrrE7Xqy4ALgahH5H07X3iUi8ld3S8q4MBBW1cS3updw\nTga5bDpQo6qNqtoOLAGmuFxTtnwiIicDxH839PcOci381wGniUi5iOThXBxa5nJNGSUigtMPvF1V\nH3O7nmxQ1Z+qaqmqnorzf7xKVXO6Raiqu4GPReSM+KZpwDYXS8qGWmCyiBTE/86nkeMXubtYBtwU\nf3wT8I/+3kFOLeCuqh0icgewHGdkwDOqutXlsjLtAuC7wGYR+SC+7T5Vfc3Fmkxm/AhYFG/YVAM3\nu1xPRqnqWhF5CdiAM6rtfXLwbl8ReQG4CBglImHgAeBh4G8i8n2c2Y2v7/f92h2+xhjjPbnW7WOM\nMeYIWPgbY4wHWfgbY4wHWfgbY4wHWfgbY4wHWfgbY4wHWfgbY4wHWfgbY4wH/R/JsVZGqz49aQAA\nAABJRU5ErkJggg==\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x1d4278a3358>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"X = np.random.uniform(0, 10, 2000)\n",
"Y = np.sin(X) + np.random.normal(0, 0.05, X.shape)\n",
"px, py = segments_fit(X, Y, 8)\n",
"pl.plot(X, Y, \".\")\n",
"pl.plot(px, py, \"-or\");"
]
}
],
"metadata": {
"anaconda-cloud": {},
"kernelspec": {
"display_name": "Python [default]",
"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.5.2"
}
},
"nbformat": 4,
"nbformat_minor": 1
}
@cerlymarco
Copy link

This can be done automatically by Linear Trees... A sklearn compatible implementation is available here

@dankoc
Copy link

dankoc commented Jul 3, 2021

Thanks for posting this!

Here are a few suggested edits to pick the number of segments automatically by optimizing AIC and/ or BIC. The implementation is similar to the heuristic strategy presented in this paper: https://discovery.ucl.ac.uk/id/eprint/10070516/1/AIC_BIC_Paper.pdf

def segments_fit(X, Y, maxcount):
    xmin = X.min()
    xmax = X.max()
    
    n = len(X)
    
    AIC_ = float('inf')
    BIC_ = float('inf')
    r_   = None
    
    for count in range(1, maxcount+1):
        
        seg = np.full(count - 1, (xmax - xmin) / count)

        px_init = np.r_[np.r_[xmin, seg].cumsum(), xmax]
        py_init = np.array([Y[np.abs(X - x) < (xmax - xmin) * 0.1].mean() for x in px_init])

        def func(p):
            seg = p[:count - 1]
            py = p[count - 1:]
            px = np.r_[np.r_[xmin, seg].cumsum(), xmax]
            return px, py

        def err(p): # This is RSS / n
            px, py = func(p)
            Y2 = np.interp(X, px, py)
            return np.mean((Y - Y2)**2)

        r = optimize.minimize(err, x0=np.r_[seg, py_init], method='Nelder-Mead')
    
        # Compute AIC/ BIC. 
        AIC = n * np.log10(err(r.x)) + 4 * count
        BIC = n * np.log10(err(r.x)) + 2 * count * np.log(n)
        
        if((BIC < BIC_) & (AIC < AIC_)): # Continue adding complexity.
            r_ = r
            AIC_ = AIC
            BIC_ = BIC
        else: # Stop.
            count = count - 1
            break
        
    return func(r_.x) ## Return the last (n-1)

@jmox0351
Copy link

Great stuff, thanks! I used this to model draining out a tank to capture the formation of vorticies.

@mylo19
Copy link

mylo19 commented Nov 29, 2021

That's really helpful, thanks. I want to use a part from your implementation in a project I am working on. Is that ok? Could you add a license, or specify whether you want some reference etc?

@dankoc
Copy link

dankoc commented Nov 29, 2021

No objections from me - feel free to use anything I've contributed under MIT no attribution license: https://opensource.org/licenses/MIT-0.

Be sure to credit ruoyu0088 for their contributions as they wish!

@KateZi
Copy link

KateZi commented Jun 13, 2022

Hi! Sometimes, depending on the number of segments and the signals I pass, I get None in the return value. I believe it stems from the py_init initializing some indices to NaN. Why could it be? It leads me to a question: why do you multiply by 0.01?
Thank you!

@dushyant-fire
Copy link

dushyant-fire commented Sep 20, 2022

@KateZi the multiplication factor of 0.01 ensures that only a small set of y-axis data (here 1% around the initial x-values) is selected to initialize the y-values corresponding to the initial x-values.

I had a data set which had a relatively larger difference between x-intervals than 1%, wherein I had to increase this parameter value to 0.02 (2%). That ensured that there are some y-values to initialize the optimization parameters.

@jef07
Copy link

jef07 commented Feb 12, 2023

Hello, thanks for the material, really useful! I have a question: is there any such code available to accomodate for 3 dimensions?

@QINQINKONG
Copy link

QINQINKONG commented Feb 20, 2023

Thanks for posting this!

Here are a few suggested edits to pick the number of segments automatically by optimizing AIC and/ or BIC. The implementation is similar to the heuristic strategy presented in this paper: https://discovery.ucl.ac.uk/id/eprint/10070516/1/AIC_BIC_Paper.pdf

def segments_fit(X, Y, maxcount):
    xmin = X.min()
    xmax = X.max()
    
    n = len(X)
    
    AIC_ = float('inf')
    BIC_ = float('inf')
    r_   = None
    
    for count in range(1, maxcount+1):
        
        seg = np.full(count - 1, (xmax - xmin) / count)

        px_init = np.r_[np.r_[xmin, seg].cumsum(), xmax]
        py_init = np.array([Y[np.abs(X - x) < (xmax - xmin) * 0.1].mean() for x in px_init])

        def func(p):
            seg = p[:count - 1]
            py = p[count - 1:]
            px = np.r_[np.r_[xmin, seg].cumsum(), xmax]
            return px, py

        def err(p): # This is RSS / n
            px, py = func(p)
            Y2 = np.interp(X, px, py)
            return np.mean((Y - Y2)**2)

        r = optimize.minimize(err, x0=np.r_[seg, py_init], method='Nelder-Mead')
    
        # Compute AIC/ BIC. 
        AIC = n * np.log10(err(r.x)) + 4 * count
        BIC = n * np.log10(err(r.x)) + 2 * count * np.log(n)
        
        if((BIC < BIC_) & (AIC < AIC_)): # Continue adding complexity.
            r_ = r
            AIC_ = AIC
            BIC_ = BIC
        else: # Stop.
            count = count - 1
            break
        
    return func(r_.x) ## Return the last (n-1)

Hi dankoc! Thanks for sharing the code which really helps me. One little clarification, shouldn't it be natural log in the AIC/BIC formula?

@dankoc
Copy link

dankoc commented Feb 20, 2023 via email

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment