segments_fit.ipynb

segments_fit2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "abc13f67-34d3-4207-ae5d-b40d1d9dda9d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import polars as pl\n",
    "import io\n",
    "from matplotlib import pyplot as plt\n",
    "import shapely"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "11709392-f1a2-40c0-9d9b-574c12faf3ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = \"\"\"3.0;6.37\n",
    "4;6.38\n",
    "5;6.38\n",
    "11;6.42\n",
    "12;6.37\n",
    "14;6.33\n",
    "15;6.34\n",
    "16;6.35\n",
    "17;6.35\n",
    "19;6.37\n",
    "20;6.4\n",
    "21;6.41\n",
    "22;6.39\n",
    "23;6.37\n",
    "24;6.39\n",
    "25;6.4\n",
    "26;6.35\n",
    "27;6.32\n",
    "28;6.31\n",
    "30;6.31\n",
    "31;6.3\n",
    "32;6.3\n",
    "33;6.29\n",
    "34;6.29\n",
    "35;6.27\n",
    "36;6.28\n",
    "37;6.28\n",
    "38;6.29\n",
    "40;6.34\n",
    "41;6.33\n",
    "42;6.33\n",
    "43;6.35\n",
    "45;6.36\n",
    "46;6.35\n",
    "47;6.36\n",
    "49;6.33\n",
    "50;6.35\n",
    "51;6.35\n",
    "52;6.35\n",
    "54;6.35\n",
    "55;6.34\n",
    "56;6.34\n",
    "57;6.35\n",
    "58;6.35\n",
    "59;6.34\n",
    "60;6.33\n",
    "61;6.32\n",
    "62;6.33\n",
    "63;6.32\n",
    "64;6.31\n",
    "65;6.31\n",
    "66;6.31\n",
    "67;6.31\n",
    "69;6.32\n",
    "70;6.3\n",
    "72;6.29\n",
    "73;6.29\n",
    "74;6.29\n",
    "75;6.3\n",
    "76;6.28\n",
    "77;6.28\n",
    "78;6.29\n",
    "79;6.29\n",
    "80;6.3\n",
    "81;6.31\n",
    "82;6.32\n",
    "84;6.32\n",
    "86;6.33\n",
    "87;6.33\n",
    "88;6.35\n",
    "90;6.33\n",
    "92;6.34\n",
    "93;6.33\n",
    "94;6.32\n",
    "96;6.33\n",
    "97;6.33\n",
    "98;6.34\n",
    "99;6.33\n",
    "101;6.31\n",
    "102;6.32\n",
    "103;6.32\n",
    "105;6.32\n",
    "106;6.33\n",
    "108;6.33\n",
    "109;6.31\n",
    "110;6.31\n",
    "111;6.34\n",
    "112;6.35\n",
    "114;6.35\n",
    "116;6.34\n",
    "117;6.35\n",
    "118;6.35\n",
    "119;6.34\n",
    "120;6.34\n",
    "121;6.33\n",
    "123;6.32\n",
    "124;6.32\n",
    "125;6.32\n",
    "126;6.33\n",
    "128;6.32\n",
    "129;6.33\n",
    "130;6.33\n",
    "131;6.32\n",
    "132;6.33\n",
    "133;6.33\n",
    "134;6.32\n",
    "135;6.32\n",
    "136;6.32\n",
    "137;6.32\n",
    "138;6.32\n",
    "139;6.32\n",
    "141;6.32\n",
    "143;6.33\n",
    "144;6.33\n",
    "145;6.32\n",
    "146;6.32\n",
    "147;6.32\n",
    "148;6.3\n",
    "149;6.28\n",
    "150;6.29\n",
    "151;6.28\n",
    "152;6.28\n",
    "153;6.27\n",
    "154;6.26\n",
    "155;6.26\n",
    "156;6.28\n",
    "157;6.25\n",
    "158;6.25\n",
    "159;6.24\n",
    "160;6.24\n",
    "161;6.22\n",
    "162;6.21\n",
    "163;6.2\n",
    "164;6.18\n",
    "165;6.18\n",
    "166;6.18\n",
    "168;6.17\n",
    "169;6.17\n",
    "170;6.16\n",
    "171;6.15\n",
    "172;6.15\n",
    "173;6.18\n",
    "174;6.16\n",
    "175;6.16\n",
    "176;6.16\n",
    "177;6.17\n",
    "178;6.15\n",
    "179;6.15\n",
    "180;6.16\n",
    "181;6.16\n",
    "182;6.16\n",
    "183;6.16\n",
    "184;6.16\n",
    "185;6.16\n",
    "186;6.17\n",
    "187;6.16\n",
    "188;6.16\n",
    "189;6.16\n",
    "190;6.16\n",
    "191;6.16\n",
    "192;6.16\n",
    "193;6.17\n",
    "194;6.16\n",
    "195;6.15\n",
    "196;6.16\n",
    "197;6.15\n",
    "198;6.16\n",
    "199;6.14\n",
    "200;6.14\n",
    "201;6.13\n",
    "202;6.13\n",
    "203;6.13\n",
    "204;6.13\n",
    "205;6.13\n",
    "206;6.13\n",
    "207;6.14\n",
    "208;6.14\n",
    "209;6.14\n",
    "210;6.14\n",
    "211;6.14\n",
    "212;6.14\n",
    "213;6.15\n",
    "214;6.15\n",
    "215;6.14\n",
    "216;6.12\n",
    "217;6.12\n",
    "218;6.12\n",
    "219;6.12\n",
    "220;6.11\n",
    "221;6.12\n",
    "222;6.11\n",
    "223;6.11\n",
    "224;6.09\n",
    "225;6.11\n",
    "226;6.1\n",
    "227;6.09\n",
    "228;6.1\n",
    "229;6.1\n",
    "230;6.1\n",
    "231;6.1\n",
    "232;6.09\n",
    "233;6.09\n",
    "235;6.1\n",
    "237;6.08\n",
    "238;6.08\n",
    "239;6.08\n",
    "240;6.09\n",
    "241;6.08\n",
    "242;6.1\n",
    "243;6.11\n",
    "244;6.1\n",
    "245;6.12\n",
    "246;6.13\n",
    "247;6.14\n",
    "248;6.15\n",
    "250;6.13\n",
    "251;6.12\n",
    "252;6.11\n",
    "253;6.12\n",
    "254;6.14\n",
    "255;6.14\n",
    "256;6.15\n",
    "257;6.13\n",
    "258;6.12\n",
    "259;6.1\n",
    "260;6.12\n",
    "261;6.12\n",
    "262;6.12\n",
    "263;6.11\n",
    "264;6.12\n",
    "265;6.13\n",
    "266;6.13\n",
    "267;6.15\n",
    "268;6.15\n",
    "269;6.14\n",
    "270;6.16\n",
    "271;6.12\n",
    "272;6.12\n",
    "273;6.11\n",
    "274;6.11\n",
    "275;6.11\n",
    "276;6.1\n",
    "277;6.09\n",
    "278;6.08\n",
    "279;6.06\n",
    "280;6.07\n",
    "281;6.07\n",
    "282;6.08\n",
    "283;6.08\n",
    "284;6.06\n",
    "285;6.08\n",
    "286;6.08\n",
    "287;6.07\n",
    "288;6.08\n",
    "289;6.07\n",
    "290;6.08\n",
    "291;6.08\n",
    "292;6.08\n",
    "293;6.09\n",
    "294;6.08\n",
    "295;6.08\n",
    "296;6.06\n",
    "297;6.05\n",
    "298;6.05\n",
    "299;6.06\n",
    "300;6.05\n",
    "301;6.06\n",
    "302;6.05\n",
    "303;6.04\n",
    "304;6.05\n",
    "305;6.05\n",
    "306;6.06\n",
    "307;6.05\n",
    "308;6.05\n",
    "309;6.04\n",
    "310;6.05\n",
    "311;6.04\n",
    "312;6.05\n",
    "313;6.05\n",
    "314;6.07\n",
    "315;6.06\n",
    "316;6.06\n",
    "317;6.06\n",
    "318;6.05\n",
    "319;6.05\n",
    "320;6.05\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "fcd6e4d3-c163-4b30-abc6-4344d3cd2d5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import optimize\n",
    "import numpy as np\n",
    "\n",
    "def segments_fit(X, Y, count):\n",
    "    xmin = X.min()\n",
    "    xmax = X.max()\n",
    "    ymin = Y.min()\n",
    "    ymax = Y.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(\n",
    "        err, x0=np.r_[seg, py_init], method='SLSQP', \n",
    "        bounds=[(0, None)] * len(seg) + [(ymin, ymax)] * len(py_init), \n",
    "        constraints=[\n",
    "            optimize.LinearConstraint([1] * len(seg) + [0] * len(py_init), 0, xmax - xmin)\n",
    "        ])\n",
    "    \n",
    "    return func(r.x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "2b27e07d-91b9-4349-a451-f062a7eb2854",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiwAAAGdCAYAAAAxCSikAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/GU6VOAAAACXBIWXMAAA9hAAAPYQGoP6dpAABx5UlEQVR4nO3deVxU9f7H8dewigugqICKuCbgBoqmUqaWW7abmrlki6VZltWtrH6V1XXpZtlKaV3NVNrMsuu+ayouBO47rgi5Ai6ILPP748AgIsIoMAO8n4/HPJw553vO+ZzD3Dufzvl+P1+T2Ww2IyIiImLHHGwdgIiIiEhBlLCIiIiI3VPCIiIiInZPCYuIiIjYPSUsIiIiYveUsIiIiIjdU8IiIiIidk8Ji4iIiNg9J1sHUFQyMzM5fvw4VapUwWQy2TocERERKQSz2cy5c+eoVasWDg7530cpMwnL8ePH8fPzs3UYIiIicgOOHj1KnTp18l1fZhKWKlWqAMYJu7u72zgaERERKYzk5GT8/Pwsv+P5KTMJS/ZjIHd3dyUsIiIipUxB3TnU6VZERETsnhIWERERsXtKWERERMTuKWERERERu6eERUREROyeEhYRERGxe0pYRERExO4pYRERERG7p4RFRERE7J4SFhEREbF7SlhERETE7ilhEREREbunhKWIzYg8TNj45cyIPGzrUERERMoMJSxFLHzlAeISUwhfecDWoYiIiJQZSliK2PBODant6cbwTg1tHYqIiEiZYTKbzWZbB1EUkpOT8fDwICkpCXd3d1uHIyIiIoVQ2N9v3WERERERu6eERUREROyeEhYRERGxe0pYRERExO4pYRERERG7p4TFTqjgnIiISP6UsNgJFZwTERHJnxIWO6GCcyIiIvlT4TgRERGxGRWOExERkTJDCYuIiIjYPasTlri4OAYOHIiXlxcVK1YkODiYqKioQm27du1anJycCA4OzrNu9uzZBAUF4erqSlBQEHPmzLE2NBERESmjrEpYzp49S1hYGM7OzixYsICdO3cyceJEPD09C9w2KSmJwYMHc+edd+ZZt379evr168egQYPYsmULgwYNom/fvmzYsMGa8ERERKSMsqrT7euvv87atWtZs2aN1Qd65JFHaNy4MY6Ojvz+++/ExMRY1vXr14/k5GQWLFhgWdajRw+qVq1KREREofavTrciIiKlT7F0up07dy6hoaH06dOHmjVrEhISwpQpUwrcburUqRw4cIB33nnnmuvXr19Pt27dci3r3r0769aty3efqampJCcn53qVRgUVjFNBORERESsTltjYWMLDw2ncuDGLFi1i2LBhjBw5kunTp+e7zb59+3j99deZOXMmTk5O12yTkJCAt7d3rmXe3t4kJCTku99x48bh4eFhefn5+VlzKnajoIJxKignIiJiZcKSmZlJq1atGDt2LCEhITzzzDMMHTqU8PDwa7bPyMjg0UcfZcyYMdxyyy3X3bfJZMr12Ww251l2pdGjR5OUlGR5HT161JpTsRsFFYxTQTkRERG49i2PfPj6+hIUFJRrWWBgILNnz75m+3PnzrF582aio6N57rnnACPpMZvNODk5sXjxYrp06YKPj0+euyknTpzIc9flSq6urri6uloTvl0a2M6fge38b3i9iIhIeWDVHZawsDD27NmTa9nevXvx97/2D6q7uzvbtm0jJibG8ho2bBhNmjQhJiaGW2+9FYD27duzZMmSXNsuXryYDh06WBOeiIiIlFFW3WEZNWoUHTp0YOzYsfTt25eNGzcyefJkJk+ebGkzevRo4uLimD59Og4ODjRr1izXPmrWrEmFChVyLX/hhRfo2LEjEyZM4P777+ePP/5g6dKl/PXXXzd5eiIiIlIWWHWHpU2bNsyZM4eIiAiaNWvG+++/z6RJkxgwYIClTXx8PEeOHLEqiA4dOvDjjz8ydepUWrRowbRp0/jpp58sd2BERESkfNPkhyIiImIzmvxQREREygyr+rBI8ZoReZjwlQcsQ5iz32uUkIiIlHdKWOzI1UXist8rYRERkfJOj4TsyJVF4lQwTkREJIc63YqIiIjNqNOtiIiIlBlKWERERMTuKWERERERu6eERUREROyeEpaCZGTAzp22jkJERKRcUx2W60lP57uRtxO3L4qDPfqw6HxFTDjRNbA2YQ29cXZ0xsnBCWcHZ5wdnfP9t6A2Tg5OODs642hyxGQyGcdOTYXYWAgIgOxlIiIi5ZQSluv57DO+vRxJ5G3A+VmWxd/vMl7FwdnBGWeTE84XL+GbbGZ+7VepP3pC8RxMRESklFAdluu5dIkPn23JgZN7SXOEi86OrG/Qkpr1auLr6Ux6ZjppmWmkZaTl+bcw6wrj4x11GPXz0aI5HxERETtT2N9v3WG5ngoVePXLGLjvPli6FMgAz1hYNgVatbqpXZvNZjLMGbkTmt3bSevbh/TTJ/mkPXx+K+xOTwCzWY+FRESkXFOn24K4ucEff0CnTsbnxES46y7YsuWmdmsymXBycMLN2Q13V3e8YuPx6dEHvwMnqZ8I7Y4Z7XZ7pENc3E0dS0REpLRTwlIYFSvCn3/C7bcbn8+ehTvvhG3bimb/MTFGQnTypPG5dWsCOtwLwK4awO7dRXMcERGRUkoJS2FVrgzz5kGHDsbn06eNpOVmhzxHRUGXLsb+ANq2haVLadKqOwAnK8HpnVE3dwwREZFSTgmLNapUgQULjKQCjDsiXbrAnj03tr+NG42k5+xZ43P79rB4MXh6UimoJXUTjcW7YzfcdOgiIiKlmRIWa7m7w6JF0Lq18fmff6BzZ9i3z7r9rF8PXbtCUpLx+fbbjf16eBifAwMJPGW83X2ymMZQi4iIlBJKWAowI/IwYeOXMyPycM5CT0/jTkhwsPE5Pt640xIbm/82V1qzBrp1g+RkABJatTPu3FSpkrP9lC00OucKwO4UDWsWEZHyTQlLAcJXHiAuMYXwlQdyr6hWzRjq3Ly58fnYMeNOy6FD+W8DsHIl9OgB588DsMY/mEfv/z+oVCnPMT0veQGwy+1Czp0YERGRckgJSwGGd2pIbU83hndqmHellxcsWwZNmxqfjxyBzp15OaDCtbdZuhTuvhsuXgQgrt0dvP3kWJ7o1vSaxwysHgjA7urceD8ZERGRMkCVbovCP/8Yw5Kzhx83bAirVkHt2jltFi6EBx4w5ggCuOce+OUXqFAh/91+NAafC+/ikAkX6k2mwuNDi+0UREREbKGwv9+6w1IUvL1h+XJo3Nj4fOCA8XgoPt74/L//wf335yQrDzwAs2dfN1kBqBnUhqopkOkA+/ZFFl/8IiIidk4JS1Hx9YUVK4y7K2CMGurSBb77Dh56CC5fNpY//DD8/DO4uBS4S1NgIAFZI4V2Hd9aTIGLiIjYPyUsRal2beNOS716xufdu+GppyAta6LDRx6BiAhwdi7c/vz9CTjraOzq3MGij1dERKSUUMJS1OrWNe601K2be/nAgfDDD+BkxXyTDg4EOtQEYLfDmZzER0REpJxRwlIc6tUzkpb69Y3PTz0F06ZZl6xkCfAwHjHt8jLD/v1FF6OIiEgpYv0vqBROgwbGUORjx3ISlxsQ6BcC/MWe6pC5aycOgYFFF+M1zIg8TPjKA5Yh2dnvB7bzL9bjioiIXI/usBQnZ+ebSlYA6jVph0s6pDjDkV3FP1LoyqJ31y2AJyIiUoKUsNg5p6BmND5jvN995O9iP96VhfKuWzRPRESkBOmRkL1r3JjAk7CjJuw6u5cexXy4e+L+ovr+KXRt9hKO99ynR0EiImIXdIfF3rm5EZDhCcDu9AQozsLEqam8NnMIPf1W8fKX9xsF7g7nM4GjiIhICbI6YYmLi2PgwIF4eXlRsWJFgoODiYqKyrf9X3/9RVhYGF5eXri5uREQEMAnn3ySq820adMwmUx5XpcuXbL+jMqgwErGXY7dHulw/HjxHWj7dlbWTgfg03Ywf+cfEBgI48blFL4TERGxAaseCZ09e5awsDA6d+7MggULqFmzJgcOHMDT0zPfbSpVqsRzzz1HixYtqFSpEn/99RfPPPMMlSpV4umnn7a0c3d3Z89VE/xVKKB0fXkR4N0M2MKuGsCuXbnnKCpCyRvXsLd6zuchD8CW8BR833gDpk+HL780qveKiIiUMKsSlgkTJuDn58fUqVMty+plV3XNR0hICCEhIbna//bbb6xZsyZXwmIymfDx8bEmnHKjSeP2ED+Tk5Xg9M4ovO66q1iO8/eOpVADaidDde/6bOEgjz0IC2eAw+7dcOed0L8/TJxoTEUgIiJSQqx6JDR37lxCQ0Pp06cPNWvWJCQkhClTplh1wOjoaNatW8cdd9yRa/n58+fx9/enTp063HPPPURHR193P6mpqSQnJ+d62bsZkYcJHrOY4DGLmRFZ+L4hf5m9qZtovN8du6FY4gobv5w1RzcD0CYOKqS8iotDBZY0hI/7XVG1NyKCSw0b82HP4bR+ewEzIg9bts8+p+zPIyOicy23V1fHbw9xlPZrKiJS1KxKWGJjYwkPD6dx48YsWrSIYcOGMXLkSKZPn17gtnXq1MHV1ZXQ0FBGjBjBU089ZVkXEBDAtGnTmDt3LhEREVSoUIGwsDD27duX7/7GjRuHh4eH5eXn52fNqdhE+MoDJKakkZiSZlVtk0nHHC2TIO4+uatY4jp1KontLicAaJjsSUKKH3UchgPwRlA8UV+9RZJbFQAqpFzg1YVf88M3I1j9/dw89VqyP8/berxU1HGxl3oz16uBU9quqYhIUbMqYcnMzKRVq1aMHTuWkJAQnnnmGYYOHUp4eHiB265Zs4bNmzfz9ddfM2nSJCIiIizr2rVrx8CBA2nZsiW33347P//8M7fccguff/55vvsbPXo0SUlJltfRo0etORWbGN6pIZ5uzni6OVtV22RAr1bUTzRmd96dUvTnObxTQ25LiSfa1xiB1NSjKbU93Xi3y0geCnyItMw0+qf/xNxfFjC3zd2W7YJOHGTy1yP5ft1kgpwvW84pu35Lrxa1SkUdF3upN3O9Gjil7ZqKiBQ1k9lc+HGy/v7+dO3alW+//dayLDw8nA8++IC4uLhCH/SDDz7ghx9+yNPJ9kpDhw7l2LFjLFiwoFD7TE5OxsPDg6SkJNzd3QsdS2nx9SONGB54gF574X9fJYKHR5HuP+mrT/A8+RIAJz3HUv2F0QCcSTlDy69bciz5GEOChzD1/qmwfj08+yzExOTsoFo1GD8ennwSHDRaXkRECqewv99W/bKEhYXlSTL27t2Lv791xcXMZjOpqanXXR8TE4OvOnZaBFQPAGB3dYw5iorY3zuWAOCfCNXb5PQvquZWjZkPzcTB5MC0mGn8uP1HaN8eNm2CTz+F7C/XmTPw9NPQoQPExhZ5fCIiUr5ZlbCMGjWKyMhIxo4dy/79+5k1axaTJ09mxIgRljajR49m8ODBls9ffvklf/75J/v27WPfvn1MnTqVjz76iIEDB1rajBkzhkWLFhEbG0tMTAxPPvkkMTExDBs2rAhOsWwIrNcGgIOecGnHliLf/+YTxj5DjwPBwbnWdfTvyJu3vwnAM/97hoNnDxozT48cCbt3w6OP5jTesAF69y7eAnciIlLuWJWwtGnThjlz5hAREUGzZs14//33mTRpEgMGDLC0iY+P58iRI5bPmZmZjB49muDgYEJDQ/n8888ZP3487733nqVNYmIiTz/9NIGBgXTr1o24uDhWr15N27Zti+AUy4aaQW3wTIFMB9i3r4gnQUxJIcoUD0DrDG+oWDFPk7fveJsOfh1ITk1mwG8DSM80Cszh6wszZ8KyZZA9xD0mBubPL9oYRUSkXLOqD4s9K+t9WDh4kA7/bsB6P/jpYCh9p20qun1v2ECjWe04UA0WJ3Sja/iiazY7lHiI4K+DSUpN4q3b3+L9Lu/nbvDHH0Y5fzAeG61dCyZT0cUphp9/hg8+gH794M03bR2NiMhNKZY+LGJD/v4EnHUEYPe5g0W667Ob1nCgmvG+dWD+lWzredbjm3u+AeDfa/7NqkOrcje4915o2tR4v349rF5dpHEKxt2sRx6Bbdvgrbfgq69sHZGISIlQwlICiqQwmYMDAQ41AdjtcAbS0oosnr93LgWg/lmo1vaO621Kv2b9eDz4ccyYGThnIGdSzuTs88OVTL+jf07jsWPtpijbjSgo9qI8t+vty1LY74MvYPDgXP2D0p9/nqWfzbzp44uI2DslLCWgqAqTBXoYtTd2eZlh//4iiycqu8NtvAlatixw+896fsYtXrdwLPkYT819CrPZbNnnB5Wbc8TD22i4eDErZswvtYXOCvq7FWXBuevtK3zlAZpvWEr7d16AzEwA9noZ1YedMjNp+9ow2Lv3pmMQEbFnSlhKQFEVJgvwM+Zk2lMdMnftLJp4UlLY7PgPAK0zvcHNrcDtK7tUJqJ3BM4OzszZPYfJUZMt++wR7MePnXPusryzY65dFzrL3LuHHWNGcGDlHC6mXcy1rqC/W1EWnLvevv7tcIAv5n6IU1aysu/+/jw16lu2teoIgPul83DPPXD27E3HISJir9TpthRJ/3EWlXYM4LITHHR6hXpv/ufmdxoZScOI9sRWg6X/9ODOrwpXqA9g4rqJvLLkFdyc3Nj89GaCagQZKy5dggYNIN4YecSOHRAUdPOxFjWzmc8eqs0LwfGWRVWcK+NTxRffKr74VPbBt/JV/2Ytr16xOg6mEsj3582DBx/MeQT4+OPw7bdGcb5z5yAszOjPAsbklAsWgLNz8cclIlJECvv7bdVszWJbTkHNaLwGdtSE3Uf+pl4R7PPMxlXEZnW4bRV0p1Xbjmo/isWxi1l8YDGP/PoIG4dupIJTBahQAV5+GV55xWg4YQJ8/30RRFvEtm9nqreRrDhnQJojnEs7z7kz+9h3Jv95rAAcTY54V/a+blKTvczNueC7Vte0aBE89FBOsjJoEEyZklNJuEoVmDsX2raFkyeNoeUvvghffnljxxMRsWO6w1KapKTw8JCKzA6Cj3f6MeqnIwVvU4Clw7vT1WcxDc/A/gEbjB8/KyScT6Dl1y05ceEEz7V5js/vzpr/6fx58Pc3KuA6Ohp9brLrtNiJA++OpJHpcxwz4Z//gHMmJFQ2XvEhjUgY+CDxHg4knE8g/ny88e+5eE5dPIWZwv/Pxt3VPXcyUylvUuNT2Qevil45d22WLTMe81y6ZHx+5BGYMcO4lldbuxa6dIHLl43PX3wBVxRzFBGxZ7rDUha5uRGY7gkksjstwRgtcpN1Tjaf3AI+0DreBC1aWL29T2Ufpt0/jbtn3c0Xm76gW8Nu3NvkXqhcGV54Ad55BzIy4MMP7W4I7uxtP0ML6HQIvFZvgo8/xj0igltOA4f3w9yJ8Nxz8N6nueZuSstI4+TFk8Sfi8+TzCRcSMi1/FL6JZJTk0lOTWbP6etPqeDk4IR3JW98qYxP9D58u2bicx58Gofg+3xvfI5vtCQ6FZwq5GwYFgaTJ8OQIcbnF16AW26Brl2L/JqJiNiK7rCUMjMHBzOw4RY6HoJVHxyD2rVvfGcXL9Ln8Ur8GgQfbq/Fv34p/ASWVxu1cBSTNkzCy82LrcO3UqtKLePuir+/cbfF1RUOHjQq49qDPXu49aMANtaB8B31GfZz1vxHy5YZdyeunK/JxwcmToT+/a1KEM1mM8mpyXmTmis/Z/176uIpq8L3rOCZ6+6MT2UffFf/jc//VuF7HnxMVfD9bTHVmt+KScX7RMSOFfb3WwlLKfP3vwbSuvJMalyAE92WwF133fjO1q2j/s9hHKoKy070pMuXN15OPzU9lXbftSMmIYY769/J4kGLjccbr74K/8nqHPyvfxl3WuzAkX+/in/6fzCZ4Xilt/H515iclZcvGwnK++9DSkrO8s6djf4hgYFFHk9aRhr/rJ5PwtD+xDumGI+lWt9Cwj2diE85kSvZSc3If+LQqzk7OONd2TvffjbZn70re+e+ayMiUkKUsNiJGZGHCV95gOGdGjKwnXWzWl/L+clfUCX+eQBOVR2P18jXbjiuQxPe5/+CvwNgyYU3uevDD67bvqDz2H1qN60nt+Zi2kXG3zme1257DeLjyahXD8fLl0mrWIk5v63l0+jTRXY9btS/7/flrVYJ3H4YatSYyYMPhjGwnX/u8/SFI4Ofpu7qxTkbOjmxvf9TPOXXk0suFXilexOAQv2Nr9x3nm02bjQe4SQnG427d4fffzc6MF/BbDaTlJpkSV5+id7Ggp27CPLLpJr7ReITj5GwbT0Jzqmczjsl1HVVcvKgrmetfEdI+ZgrUmvHEaq26Qi+vpbzae1flajDZ23+NxWR0kl9WOzElQXBiuL/zCsHBeO3G456wO7YjYTdRFxdzm0GoNFpiMioxfXu1RTmPAKqB/Bpj08Z+udQ3lrxFp3rd6Zt7bb8EdKDhzbMxfniBZI/+oS41n2K7HrckIMHWeiVAED7Y9X5qaYHCVnx5DrP17vQv+dr3FKzPR8sm0zts/GQnk6zH77mV/fZvH/nUMJdHMFkKtTf+OricJZtXE4bCUp2snLXXTBnTp5kBcBkMuFZwRPPCp4E1gjkrZ/MmM83JDHejXmPZU2rcPgwtGlD6pmTnKgEs9rdTuBnrzDq1xWcTjlBpukM6aZEHJ0SuWw+w2XzGTClcyE9iV2nkth1atd1L99733nwf7PiLOeTkJRChhnb/k1FpMxT4bhiVpTFxQAIDCTwpPF298mbKx6X4mKUgW8Vb+LOftd/tFTY83gy5EkeDnqY9Mx0Hp39KOdSz2F69V+kZw3FHbTxDxq5YdNCcgm/TmOtUSiWJjV75Tqvq89zeKeG7G3dkdW/rYD/+z9wcQGgTvJJvpkzlp///IDx6TvpmnyIF1t45Cqbf7Ur9539fnTtVCNBSUw0GnXqZEwiWYgCfteKFzD6Dc2Zg5ODC37J8NriNdy39Chj7nyRppWfZEjQWEIq/pvwbgv5rvtGWpj+R1N+YXzYEpYOWsqMB2fwn7s+5GWfh3j0lC9dYiHoBHhmPR37ulES5ilTLMfu1aKWXRcHFJGyQY+ESqEXHnLjs5aXeHlLJT767fyN7eTCBXo/WZnfAuE/22vxyk10uL3a2ZSzBH8TzJGkIwxqMYjpD0435sD54QejwcSJ8NJLRXY8a4X3rc+zTQ9x6zGI/NceY0RNYe3bZ4wcWrz42uvd3Izh21e+6tfPeV+9ek7H3e3bjX4xp7I63N52m1H4rXLlGz213KZPh8ceM947Ohr7Lmjk0IULxnaffQa7d+dadamKG17Pp3DRBWLmeNNy4xFLAicicqM0W3MZFuhijAza7XYBkpJubCcxMWzOGrAT6h1SRJEZqrpVZeZDM3EwOfDD1h+YuXUmvP56ToOPPoLUwnccLVLHjvGr2yEAep/1sS5ZAWjcGBYuhF9+ufYIrZQU2LXLSA7Cw+G116BvX6O+Tc2aRjLStCn06mXUTslOVtq3h/nziy5ZASNJfC2rj1NGBvTpkycJsTh82OgUXacOPPts7nb+/vDRR1Q4cpy7LhjzRM1z/wciIoouVhGRAihhKYUCqgcAsLs6uYffWuHkxpUc8TTet2p6EyON8nFb3dt4u+PbAAyfN5xY3wpGiXkwSvbbqPLtqdk/sKqe8b538743thOTCR5+2Lj2s2cblXyHD4eePY0RRNd7nHPxIuzcaSQnJ7Oe7bVtayQ4VarcWDzXM3Ys3Hef8T4pCe691xhuDsbjqzVrjHNp0MBIJLMfTQF07Gic3/79RuViT0/u7jAYgPmNMc47a34jEZHipkdCpVDCR+/ie2EMDplwod5kKjw+1Op9LHzmTnrWWs4tp2DPkM3QunWRx5memU7n7zvz15G/aFu7LX81n4Rzuw7GygYNjB98p5Lt9/1d/yY8FbCXkHj4+7lt0KxZ0R/EbIYTJ+DQobyvgweNf7PvMLVpY5Tgr1q16OPIdu6c8bhp61bjc5cuxt2XTz+F6OjcbV1c4NFHjeJzwcF5dnUk6Qj+k/xxyIQT/wGviN/h/vuLL3YRKfP0SKgM8w5qi2cKZDrAvn2RN7SPqFPGhHmhCabi+dHGqNw686GZeFbwZGPcRt45/2dOH4rYWPj552I5br7++YfZjnsB6H3Cy3g0UxxMJvD2hltvhX79jMcy4eHGXZTdu427LPHxsGUL/PVX8SYrkDPnUM2axufly42quFcmKz4+8N57cPQoTJ16zWQFoK5HXZq71SPTARY3BMaNu25HYxGRoqKEpRQyBQYSkNX1YVf8Nut3cO4cm52NxxGtTbWNKrTFpK5HXabcOwWA8X+NZ/nwHjkrx40r0UcKib/NYmkD433vwN43Pa3BDXNwMBKEFi1KrtNq1sihPMcLDTU6Qx8+bIyCyk5qruPuEONR2rxbgA0bYNWqYghYRCQ31WEpJXIVNGvrT+BZRyL9MtidHGv9zmJiiKplvC3qDrfX8nDQwzwV8hTfRn9Lz61j2dKqBQF/b4Xt21n50Xd0ejXvI62iLLiXva/7/vqStDbGEN2AJ5+5oX0UJp6rC6qVVGG1GZGH+WiR0aep4y01ch1zRuRhwldf4t/vf0GnHz7jkE89xjXqxu2D72Vg+3q5tk9Nz8DVyZGOt9Rg9d6TefbnnB4KwMJGkGECx/HjmVGhvuWcs7d5pXsT1WURkSKjOyylRK6iYw4OBDgY/yW82+EMpKVZta8TG1dw1ANMZghpVjIT5E3qMQk3U10uc5p77rxsmeu4xucfX/ORwtVF1m5G+MoDnI8/wdoqxr66HawEIdYlatbEk9123tbjuf4tinMp6LiJKWkkpqTlOWZ2TG9mNoRt2xhw5ygWeTQgfFVsnu1T0jIt+7jW/pbFeOJgrsTpirCxNrBoEUtnLbSca/Y2xX2+IlK+KGEpJa4uEBbgYfy7y8tsjOKwQtSeFQA0OQVV2txorVzrVHKpxJvtwjHhzIFKu5l4tzE8tumx3UafiqsUZcG94Z0acn9CJEuydtXRu7vVj4OsiefqgmolVVhteKeGeLo54+nmnOeY1yqId3VM2du7OTtY9nGt/T3b+RaqOrYBskYLAe/u+p/lXLO3USE5ESlKGiVUSu0bM5Jb+By3NDjf6lccHupd6G3f712dt1ucZsB2B2bMTCnR4l+TIicxatEoKpic2fRlGs1OYIxaWbasWI/786BW9GsUTaPTsLf/Okzt2xfr8cq66Vum89jvj9HqhBNRX6Ub/XL27IFGjWwdmoiUMholVMbVb9IOl3RIcYYjuzcUfsPkZKJcTgMQaqpd4pVKX7j1BXo26sklcxqPPOpCihPGHZbIGxvtVCjJyfyatgWA3kcrY7r11uI7VjnRo1EPTJj4u2Y68ZUxOk9/9JGtwxKRMkwJSynlFNSMxln1v3Yfjir8htHRbM7qcNvap1XRB1YAk8nEtAem4V3Jmx2el3mlW9aKceOK7Zgpc39jfkNjNFLvencbdwPkptSsVJM2tY3HQguaZ03SOHWqMVxbRKQY6P+5S6vGjXOGNifuK/RmCZtWEOee1eG2ecl0uL1azUo1jfmFgK/awh9NMOqEbLuBIdqFsGjZN1xwgbqJEHr/8GI5Rnl0d6O7AZh3Z9ZMkpcvG8XoRESKgRKW0srNjYB0TwB2pyUUunhXdofbgFNQuYQ63F5Lt4bdeLn9ywA8cT/EVQHGjy/6A124wOzzmwDofcgN0+23F/0xyqlet/QCYInbcS5XcDYWfvVV7vL+IiJFRAlLKRZYyahxsdsjDY4fL9Q2UWd2ABCa4FB8lV4LaeydY2nlHcyZijDoIcj4KQIOFO1Q2NR5c5nbKAOA3nW6GbMWS5Fo5dsK70renEs7z19PZT3bO3fOqOorIlLElLCUMjMiDxM2fjkzIg/j4lgfgF01yH8W3islJbHZ1ehw29qxDjg7F2OkBft5UzwX/3kelwwnVtSHD9ub+f3hZ5kRebhI9j8j8jAffTue5Argew7a31v8j4Ou/PsUpm3wmMUEj1lcZOdckhxMDvRs3BOA+bf75PQNmjTJmLVaRKQIKWEpZa4sYLbzfB0ATlaC0zs2F7xxdDRRvsbbUN+in+zQWuErD5CS4k2NS08B8H9dwPfkYjZP+alI9v/dkl3sqWTcUXow1hWHzl2KZL/XY22BudJeZC27H8v8k+uMWZ/BmPhx2jTbBSUiZZISllLmyoJf7e66A78kY/nu2I0Fbhu/cTnH3cEhE4JLqMLt9WQXKvOq9AAdzgeQ4QD9H4bXFk6Es2dvev9vuR5m/i3G46CHvTuXyB0lawvMlfYia10bdsXR5MiuU7s4+PzAnBX/+Q+kp9suMBEpc1Q4rjQ7fZpuo6qzpCF8uzeAJ2fuum7z/z3VkXv91tD0BGx/Zosx+Z6dSLx4huD3a3O4wiUGbIUZrv1h1qyb2ueyZ7pxV60lVL8A8R3/wOme+4ooWrlSp2mdWHV4FV/0/IIR/zcXFi82VsyaBf372zY4EbF7KhxXHnh5EXjeqIGxK+Vogc03n90OQOg/DhAUVKyhWcuzYjVm9fkRx0yY2QJ+2BEBP/544zu8fJlfTxmzCD8Q64xT1+5FFKlc7e7GWcOb982D11/PWTF+fKFHr4mIFMTqhCUuLo6BAwfi5eVFxYoVCQ4OJioq/8Jlf/31F2FhYXh5eeHm5kZAQACffPJJnnazZ88mKCgIV1dXgoKCmDNnjrWhlUsBLrUB2O12AZKS8m+YmEiUq/GYpbVjXXCyv4m6O7S6n3dqGP0gnu0F+19/GuLibmhfGcuXMqfBZQB6Vw0DV9cii1Ny69XYGN684tAKLoa1hexKwlu3wsKFNoxMRMoSqxKWs2fPEhYWhrOzMwsWLGDnzp1MnDgRT0/PfLepVKkSzz33HKtXr2bXrl289dZbvPXWW0yePNnSZv369fTr149BgwaxZcsWBg0aRN++fdmwwYqS8+VUYPUAAHZXx5jLJT9//22pcGsPHW7z88azP9LxQg3Ou0L/bue4/OSQG/qv9HV/hvNPZfBMgS49VCyuOAXVCKKuR10upV9ixaGVue+yFGMFYxEpX6xKWCZMmICfnx9Tp06lbdu21KtXjzvvvJOGDfPvMBgSEkL//v1p2rQp9erVY+DAgXTv3p01a9ZY2kyaNImuXbsyevRoAgICGD16NHfeeSeTJk264RMrLwLqhQJw0BMu7dyab7vjm5aTUMXocNuyhe073ObH0cGRGSOWUfWSic214f8ylhrFyKyRns7s40sBuO+AIy497ymGSCWbyWSy3GWZv28+3HcfBBiJNGvWwNq1NoxORMoKqxKWuXPnEhoaSp8+fahZsyYhISFMmTLFqgNGR0ezbt067rjjDsuy9evX061bt1ztunfvzrp166zad3nkHdQWj0uQ6QD79q7Pt93mvSsBaHoSKrbpUELR3Rg//+Z82+JNAD68DZZ+8dL17x5dJXP1KmbXvwRA7yq3QsWKxRKn5MjuxzJ//3zMJhO89lrOygkTbBSViJQlViUssbGxhIeH07hxYxYtWsSwYcMYOXIk06dPL3DbOnXq4OrqSmhoKCNGjOCpp56yrEtISMDb2ztXe29vbxISEvLdX2pqKsnJyble5dHvl9wJPGm8nx21Nt8CZFFndwLQ+oQjBAaWVHg37KH+7/N0itExeFCvy+zu3zvPMNn8irRtmhvOMQ+onAoV6/UusZjLs871OuPq6MqhxEPsOrULHn2UCzWziv78+Sds327bAEWk1LMqYcnMzKRVq1aMHTuWkJAQnnnmGYYOHUp4IUpxr1mzhs2bN/P1118zadIkIiIicq03mUy5PpvN5jzLrjRu3Dg8PDwsLz8/P2tOpcyYuPsSjc8Yf8ZEh/hrFyA7e5bNbkaH21A77XB7LZ+MXkWjM04kVIFXAndg/ve/c62/ZpG2zExmHzY6evbcZ2KC+ZaSDLncquRSic71OwNZj4VcXPi2zYM5DT780EaRiUhZYVXC4uvrS9BVw2EDAwM5cuRIgdvWr1+f5s2bM3ToUEaNGsW7775rWefj45PnbsqJEyfy3HW50ujRo0lKSrK8jh4teFhvWTSsS2N8UjwB+KfSOZ69rW6eNubNm4nK6nDbunZoCUZ3cyp6VOfFgP/gkg7zboEvFo6BzTkVfa9VpM28bh2z614AoHlSYx7r3rzE4y6vLLM375sHgPdLI0hyq2KsnDULDpe+6QdExH5YlbCEhYWx56q+BHv37sXf39+qg5rNZlJTUy2f27dvz5IlS3K1Wbx4MR065N/XwtXVFXd391yv8mhgO39uq20kkburmxlQPW910bioFfxTGRwzoWWLbnnW27MRw1/kPw5GDZV/3WVm63N9LPPUDGznz9rXuzCwXc73L+b3r4mtBm5p8NKjr+RaJ8Ure/bmv478RdKlJB7pEoTHqy8ZKzMyYOJEG0YnIqWdVQnLqFGjiIyMZOzYsezfv59Zs2YxefJkRowYYWkzevRoBg8ebPn85Zdf8ueff7Jv3z727dvH1KlT+eijjxg4MKeM9wsvvMDixYuZMGECu3fvZsKECSxdupQXX3zx5s+wHAjwCwZgT3XI3Lkjz3pLh9sT4NamfQlGVjSef3Muvf7xINUJHgk9xMXXX752Q7OZ2Qf+BKDHAROV7n+4BKOUBlUb0MSrCemZ6SyJzfoPkOefz+n0/O23cPKk7QIUkVLNqoSlTZs2zJkzh4iICJo1a8b777/PpEmTGDBggKVNfHx8rkdEmZmZjB49muDgYEJDQ/n8888ZP3487733nqVNhw4d+PHHH5k6dSotWrRg2rRp/PTTT9yaXYBKrqtBk/Y4Z0CKMxzZnbd2TVSiUbI/9IRjznDTUsTk4sLUp+fjc96Ymfql2HC46o4cAJs3M7u20fm6t2NzqFq1hCOVXMObAby84OmnjfcpKfD55zaKTERKO80lVBZs3UrTb1qysyYsON6FHt8sy1l3+jQ9X6jOwsbw5a4GPPtj6ZwVGGDpxBF0O/cVZhPMXlKNh/7cnysp2Tn6KZpW+A6XdDjh9ykeT4+0YbTl07LYZdz1w114V/Lm+MvHcTA5wNGj0KCBMcqralWIj1flYRGx0FxC5UnjxgSeMt7uStyXa5U5KsrS4Ta0dpsSDqxo3TXqc/51zOiT8tRtZzj64uM5K81mft1jTOfQNRY8HtSke7Zwu//tVHapzD8X/iE6PtpY6OcHvbOGl589CzExNotPREovJSxlgZsbAemeAGxLjSf43UWW2iRHNy/jZCVwyoA/j/nlW6fF3s2IPEzYhytxavk5IfEOnHWDQQ5/kPFT1gSJ27Yx2/sMACEn/Jhx4KINoy2/XBxd6NrAqKSc/VhoRuRh/nPJJ6dRZKQtQhORUk4JSxkRUMm483CgajquJxMstUmi9q8GoNkJWFel4bXrtJQC2TVXfjzuQLPTj1M5FVbVg3H/fRyOH2f/r5PZ6mOMhLrofGepPc+yINfszRh/u5VVG+Q0UMIiIjdACUsZEejdDDA6pbY8n2CpTbI50ahw2+qEI5cbNMxVs6Q0ya650qtFLbY27cOwLU0AeLfdJdaPfIjZO34BoPNB2BrYrdSeZ1mQnbBsjNvIyQsnGd6pIecaBZLuWsFooIRFRG5A6Sh5KgVq0rgdxM/kZCUY1xq82vnDqVNEVTJGzbRxrc93b9jvpIcFGdjOP1dNFfPptRx/oQ6zGl+iv98G3LPK+jyc2pBnJqr/ii3VqlKLYJ9gYhJiWLh/IYPaDTL+dovaGJMhHjoE//wD1ykMKSJyNd1hKSMqBwXjl2S83x27ETAq3G4uIx1ur2by8iL8kRnUPwuHPWGbN5jM8EDbwQVuK8XPMrx5//yche3a5bzfkHf4vYjI9ShhKSsCAwnIGim0+6TxGOhI1DJOVwTnDGgeXLoq3BaG+z29mWV+EMdM4/Pth8H74cdsG5QAOY+FFu5fSHpmVvXlKxMWPRYSESspYSkrvLwIOG/0EdiVYsyrtDmrw23zf8A1tF2+m5Zm7cb+wEcxNXHKgOdSW4CV00RI8bi19q1Uc6tG4qVEIo9lJSdKWETkJihhKUMCXWoDsNvtAiQnE5W0G4DWJ53gljI6a3GlSrw4fS8pTWfR54sVto5Gsjg6ONKjUQ/giqq3tWpBnTrG+02bjPmFREQKSQlLGRJQ3Si7v7s6sGYNmysbHW5DXeuDQxn+U3t44NSvP1SrZutI5ApXz94M5NxlOX8edu60QVQiUlqV4V+x8iewXigABz0hZdb3RPkayzPTGxE2fnmpLRp3PTMiDxM2fjkjI6LL7DmWVj0a9cCEia3/bOWzlZHMiDzMxOQr5ne64rFQ9t9Rfz8RyY8SljLEO6gtHpcg0wGWbPuDMxXBJR22X25CXGJKmSymll1Qbt7W42X2HEsrr4peVDYFAvDp2p8JX3mAdTUa5zS4ImHJ/jvq7yci+VHCUoaYrhgpNCPgMgDNT4D3HXdQ29OtTBZTu7KgXFk9x9Kse0PjsZBnte0M79SQY/UCSHNwNFZeMbQ5+++ov5+I5EezNZclmZk83seFaS0ycE2HVCd4ZoszX/96qWz3YRG7FR0fTavJrajkXInTr57G1ckVQkMhKgpMJmMyRA8PW4cpIjak2ZrLIwcHAhxqAkayAtC6QhnvcCt2LdgnGN/KvlxIu8Dqw8Ywe0vHW7PZGC0kIlII+iUrYwI9ct9SD/W71UaRiIDJZMozGaLqsYjIjVDCUsYE+AVb3rukQ9MyWOFWSpfshMVSj+XWK5JolegXkUJSwlLGNGjSHueselwt/wGXNmWzwq2UHnc1uAtnB2f2ndnHvtP7oFGjnJo5kZHGoyERkQIoYSljnIKa0fi08T70lAs01KgLsS13V3du978dyLrLYjLlPBY6dQpiY20YnYiUFkpYyprGjekQb/S4rZ1clxkbjtg4IJG8szdvqR2Qs1L9WESkEJSwlDVubnz06FQm/NWaZfVeUiEusQvZ/VhWHlrJ+cvnmZbunbNS/VhEpBCUsJRBHg8PpNb7s0lt1ESFuMQuNPFqQn3P+lzOuMzyg8tp179HzkrdYRGRQlDCUkYNbOfP2te7MLCdv61DEcFkMlkeC83bO49+XVtAoFG2n+hoSEmxYXQiUhooYRGREmEZ3rx/PmazOafjbXq6kbSIiFyHEhYRKRGd6nXCzcmNY8nH2H5ie+56LHosJCIFUMIiIiXCzdmNLvW7AFlVb6+seKuOtyJSACUsIlJiLMOb982Hpk2hUiVjhe6wiEgBlLCISInp2bgnAOuOruNs2jlo08ZYceQIHD9uw8hExN4pYRGRElPPsx5BNYLIMGfwfwtn8t+0mjkr9VhIRK5DCYuIlKjsx0I/b59LZM3GOSuUsIjIdShhEZESlT28+ZJjFAcaBOasUD8WEbkOJSwiUqLC/MJwd3XnXNoZxo1qDP5ZxQ03bTJqsoiIXIMSFhEpUc6OznRr2A0wqt5a6rFcvAjbt9swMhGxZ0pYRKTE5Zq9WfVYRKQQrE5Y4uLiGDhwIF5eXlSsWJHg4GCioqLybf/bb7/RtWtXatSogbu7O+3bt2fRokW52kybNg2TyZTndenSJevPSETsXo9GxuSHm49vJiHkio636sciIvmwKmE5e/YsYWFhODs7s2DBAnbu3MnEiRPx9PTMd5vVq1fTtWtX5s+fT1RUFJ07d+bee+8l+qq5Q9zd3YmPj8/1qlChwg2dlIjYN5/KPrT2bQ3Awkrx4OxsrFDCIiL5cLKm8YQJE/Dz82Pq1KmWZfXq1bvuNpMmTcr1eezYsfzxxx/8+eefhISEWJabTCZ8fHysCUdESrFejXsRFR/FvEOLGRIcbHS63b0bzp6FqlVtHZ6I2Bmr7rDMnTuX0NBQ+vTpQ82aNQkJCWHKlClWHTAzM5Nz585RrVq1XMvPnz+Pv78/derU4Z577slzB+ZqqampJCcn53qJSOmRPbx58YHFpLVrY1n+1YcRhI1fzozIw7YKTUTskFUJS2xsLOHh4TRu3JhFixYxbNgwRo4cyfTp0wu9j4kTJ3LhwgX69u1rWRYQEMC0adOYO3cuERERVKhQgbCwMPbt25fvfsaNG4eHh4fl5efnZ82piIiNtandhhoVa5CcmszaFjl3VFL/WkdcYgrhKw/YMDoRsTcms9lsLmxjFxcXQkNDWbdunWXZyJEj2bRpE+vXry9w+4iICJ566in++OMP7rrrrnzbZWZm0qpVKzp27Mhnn312zTapqamkpqZaPicnJ+Pn50dSUhLu7u6FPSURsaHBcwbzw9Yf+FfQUD7sa9yt3dEyjKcfeY/hnRoysJ2/jSMUkeKWnJyMh4dHgb/fVt1h8fX1JSgoKNeywMBAjhw5UuC2P/30E08++SQ///zzdZMVAAcHB9q0aXPdOyyurq64u7vneolI6WIZ3nxyHVSvDkDTo7tY+1pnJSsikotVCUtYWBh79uzJtWzv3r34+1///1giIiIYMmQIs2bNolevXgUex2w2ExMTg6+vrzXhiUgp061hNxxMDuw4uYPDHVsYC8+cgf37bRuYiNgdqxKWUaNGERkZydixY9m/fz+zZs1i8uTJjBgxwtJm9OjRDB482PI5IiKCwYMHM3HiRNq1a0dCQgIJCQkkJSVZ2owZM4ZFixYRGxtLTEwMTz75JDExMQwbNqwITlFE7FVVt6p08OsAwPzgSjkrNLxZRK5iVcLSpk0b5syZQ0REBM2aNeP9999n0qRJDBgwwNImPj4+1yOib775hvT0dEaMGIGvr6/l9cILL1jaJCYm8vTTTxMYGEi3bt2Ii4tj9erVtG3btghOUUTsWfZjoXkeJ3IWKmERkatY1enWnhW2046I2Jet/2yl5dctcXNy4/SYFNzSgFat4DoVtEWk7CiWTrciIkWtec3m1HGvQ0p6Cis71jUWbt1qTIYoIpJFCYuI2JTJZOLuRkYRufnBlY2F6ekMf2kKIyOiVURORAAlLCJiB3rdkjW82fMk2c+o6+zZyrytx1VETkQAJSwiYge61O+Ci6MLsRkn2WOUY6HDqX30alGL2p5uDO/U0LYBiojNWTX5oYhIcajsUpk7/O9gSewS5jd1IWDVZTonHqRz/5CCNxaRckF3WETELliGN2fXYzl2zHiJiKCERUTsRPbszWs8k0h2zVq4YYPtAhIRu6KERUTsQmOvxjSu1pg0UyZLG2QtVAE5EcmihEVE7Eb2XZb5jbMW6A6LiGRRwiIidsMye3OAozG8efNmSEuzaUwiYh+UsIiI3ejo35GKzhWJr5hBjA+QkgLbttk6LBGxA0pYRMRuuDq5cleDu4ArHgupH4uIoIRFROyMZXjzLVkLlLCICEpYRMTO9GzUE4DIOnCqIup4KyKAEhYRsTN+Hn608G6B2QSLGgJ798Lp07YOS0RsTAmLiNgdy+zN2f1YNm60XTAiYheUsIiI3cmevXlhI8gwoX4sIqKERUTsT7s67fB0cedMRdhQB/VjERElLCJif5wcnOje2Oh8O78xJK9ay4x1B20clYjYkhIWEbFLluHNjcH90nnm/brKxhGJiC0pYRERu9SjUQ9MmIjxhbgq8HylU7YOSURsSAmLiNilGpVq0LZKAAALGkOHUwdsHJGI2JISFhGxW3c3exDIGt6skUIi5ZoSFhGxW72aPQTAkgaQumMLnD9v44hExFaUsIiI3QrxDcE7w43zrvBXHTOsUsdbkfJKCYuI2C0HkwM9q7cDsh4LLVpk24BExGaUsIiIXevVYQiQNXvzwoU2jUVEbEcJi4jYta7N7scp08Se6nDg1D44qAJyIuWREhYRsWseFTxocM4b0GMhkfJMCYuI2L3b/HsASlhEyjMlLCJi915+5GUAVtSHC6uWQlqajSMSkZKmhEVE7F5gzab4p1Ui1QlW1DgP69fbOiQRKWFKWETE7plMJu6uditgTIaox0Ii5Y/VCUtcXBwDBw7Ey8uLihUrEhwcTFRUVL7tf/vtN7p27UqNGjVwd3enffv2LLrG/9nMnj2boKAgXF1dCQoKYs6cOdaGJiJlWK+wxwGjH4t5kYY3i5Q3ViUsZ8+eJSwsDGdnZxYsWMDOnTuZOHEinp6e+W6zevVqunbtyvz584mKiqJz587ce++9REdHW9qsX7+efv36MWjQILZs2cKgQYPo27cvGzZsuOETE5GypXOrh6iQYeKIJ+w8+jecPGnrkESkBJnMZrO5sI1ff/111q5dy5o1a27qoE2bNqVfv368/fbbAPTr14/k5GQWLFhgadOjRw+qVq1KREREofaZnJyMh4cHSUlJuLu731R8ImKf7n67IQscY5mwBF59diY8+qitQxKRm1TY32+r7rDMnTuX0NBQ+vTpQ82aNQkJCWHKlClWBZaZmcm5c+eoVq2aZdn69evp1q1brnbdu3dn3bp1+e4nNTWV5OTkXC8RKdsaVLkdyBrerKq3IuWKVQlLbGws4eHhNG7cmEWLFjFs2DBGjhzJ9OnTC72PiRMncuHCBfr27WtZlpCQgLe3d6523t7eJCQk5LufcePG4eHhYXn5+flZcyoiUgptunQHAH/VhcSVCyEz08YRiUhJsSphyczMpFWrVowdO5aQkBCeeeYZhg4dSnh4eKG2j4iI4N133+Wnn36iZs2audaZTKZcn81mc55lVxo9ejRJSUmW19GjR605FREphZ7v2gX/pApkOMCSKidh61ZbhyQiJcSqhMXX15egoKBcywIDAzly5EiB2/700088+eST/Pzzz9x111251vn4+OS5m3LixIk8d12u5Orqiru7e66XiJRtA9v509u7PZA1GaKGN4uUG1YlLGFhYezZsyfXsr179+Lv73/d7SIiIhgyZAizZs2iV69eeda3b9+eJUuW5Fq2ePFiOnToYE14IlIOZM/evKARZC5ccP3GIlJmWJWwjBo1isjISMaOHcv+/fuZNWsWkydPZsSIEZY2o0ePZvDgwZbPERERDB48mIkTJ9KuXTsSEhJISEggKSnJ0uaFF15g8eLFTJgwgd27dzNhwgSWLl3Kiy++ePNnKCJlym3t+1HlsokTleHv2L/g/HlbhyQiJcCqhKVNmzbMmTOHiIgImjVrxvvvv8+kSZMYMGCApU18fHyuR0TffPMN6enpjBgxAl9fX8vrhRdesLTp0KEDP/74I1OnTqVFixZMmzaNn376iVtvvbUITlFEyhIXJ1e6ZtYDYF79DFixwrYBiUiJsKoOiz1THRaR8uO7/z7HU0e/pO0x2OA6Ar74wtYhicgNKpY6LCIi9qBnj+cB2FQbTqycZ+NoRKQkKGERkVKnVq0mhJyrjNkEC50OwYEDtg5JRIqZEhYRKZXu9ggFNHuzSHmhhEVESqVe7R8DYFEjSF+k4c0iZZ0SFhEpldp2GoBXiomkCrBu7zK4fNnWIYlIMVLCIiKlkqOTMz3SjKKV8+ukwPr1No5IRIqTEhYRKbXubtQT0OzNIuWBEhYRKbW69xqJQyZs84Yjq/+0dTgiUoyUsIhIqeVVN4B2ZysBsCB1B/zzj40jEpHiooRFREq1u91bA1mzN181iaqIlB1KWESkVJkReZiw8csZGRFN8JjFbD7RCoBl9eHSomtXvc3eZkbk4ZIMVUSKkBIWESlVwlceIC4xhXlbj5OYksY2j474noOLLrBq5wLIzMx3m/CVqogrUlopYRGRUmV4p4bU9nSjV4taeLo5U7lyJe5IrgXA/JpJEBOT7zbDOzUs4WhFpKhotmYRKfXmfPI0DyVPodFp2Oc7FkaPtnVIIlJImq1ZRMqNu3qOwDkD9nvB3tVzbB2OiBQDJSwiUupVCWhJxxNuAMxPjoLkZBtHJCJFTQmLiJQJd1cxRgvNa5gJK1bYOBoRKWpKWESkTOjVdgAAq+rB+cX/s20wIlLklLCISJlwS4+BNDgLaY6wdIfK9IuUNUpYRKRMmLnjDKHxXgDMd/8H9u8vcJsbKShX2G1UrE6kaClhEZEyIXzlASqmhQDG7M3mBQsKtY21BeUKu42K1YkULSUsIlImDO/UkHM178UtDeLcYeuaXwu1jbUF5Qq7jYrViRQtFY4TkbIjM5N7n3Tjf/UuM3a1M6MXnAcXF1tHJSLXocJxIlL+ODhwd6VgAObVS4O1a20bj4gUGSUsIlKm3N3mUQDW14EzC0u46m1GhtHZ9/ff4aOPYPXqkj2+SBnmZOsARESKkv/dj9L03RfZURMW7fiD/nxW9AfJzIQjR2D7dtixw3ht3w67dsGlSzntHB1hwwZo3broYxApZ5SwiEjZUqMGvZK92VHzH+a7HKF/QgL4+NzYvsxmiIvLSUiyk5MdO+DChYK3z8iA116DJUvAZLqxGEQEUMIiImVQQKX2wO8saAwZixbi+NgQwKiN8tGiPaSmZ+Dq5Mgr3ZswsJ2/ZbsZkYf5aOFugo7vY9ihNTRbu4hqF5MKdcwMkwOnfPzY7lGHAzX96b9/DVWOH4Vly1j26Q/Et7uD8JUHLKOGst9feXyxrRmRh/V3sWNKWESkzIk034bHpd85XRE2rY6gXVbCEr7yAIkpaQCkpGUSvvJAzg9TQgJn3x/HjxsXEnDq+sXe4qr6stPLj33V67Kvhj97vPw54FWHdGcXMrLGXR6r3YD3Ij4AoO6EMbz7/NfEJV+21GXJrtGiH0b7cWXtHP1d7I8SFhEpczoOup+TX7zGnMAM5iesoV1mJjg4MLxTw1x3WEaE+cHs2TBtGixYwPMZGbn2k+Lkyt/+zbgUEMQqp5rsr1mPe/vfSYZbRT5atMc41i01OL73JG5Z71fvPQnALc89wemYuXjt2krjhFjGX9rGq57Bee6wiP0Y3qmh/i52THVYRKRMmvZkKx6vG02r4xA1dBOEhhorzGaIjjaSlFmz4PTpvBt36ACPPw59+oCHx40HsWIFdOlivPfzg717oUKFG9+fSBmkOiwiUq71bPUIAH/XgviFv8KJE/DJJ9CypTFq5/PPcycrderAG2/Anj1G/Zannrq5ZAWgc2e4+27j/dGjxjFF5IboDouIlE0HD9JmXAM214bvFrvxxMY0SE/P3aZCBXjwQRgyBO680xiGXNS2bTOSJLMZPD3hwAGoVq3ojyNSShXbHZa4uDgGDhyIl5cXFStWJDg4mKioqHzbx8fH8+ijj9KkSRMcHBx48cUX87SZNm0aJpMpz+vSlfUMRESsUb8+vc4YicHoDik88kA642+DhY0goWMr+PpriI83Hgt161Y8yQpA8+ZGQgSQmAhjxxbPcUTKOKs63Z49e5awsDA6d+7MggULqFmzJgcOHMDT0zPfbVJTU6lRowZvvvkmn3zySb7t3N3d2bNnT65lFfSsV0RuQr8mDzM+fTInKsNPzYyX4W+8z71DyLzfCfYOJtjHeDWq1ghHh2JIXN57DyIijKJyn38Ozz0H9eoV/XFEyjCrEpYJEybg5+fH1KlTLcvqFfA/unr16vHpp58C8N///jffdiaTCZ8bLe4kInINgf/3KXGzWrC5UiIxtRyIObGVmIQY9pzawz8X/mHh/oUs3L/Q0r6ic0VaeLcg2DuYEN8Qgn2CaVazGRWdK95cIHXqwIsvwvjxcPky/N//wQ8/3Nw+RcoZqxKWuXPn0r17d/r06cOqVauoXbs2zz77LEOHDr3pQM6fP4+/vz8ZGRkEBwfz/vvvExISkm/71NRUUlNTLZ+Tk5NvOgYRKTuyi4C19u9A1MGzDPdvyGu9RwNw4fIFtp/YTnRCNL9sWcPGuGgucZCLaReJPBZJ5LFIy34cTA408WpiuQuT/apZqaZVcYzs0p8Hvvoa1+REzDNnYnrpJbjO/8fdyPlmD7W+uiBeaWJN8babLfRWVq5ZeWFVp9vsRzQvvfQSffr0YePGjbz44ot88803DB48uMDtO3XqRHBwMJMmTcq1PDIykv3799O8eXOSk5P59NNPmT9/Plu2bKFx48bX3Ne7777LmDFj8ixXp1sRAQgbv5y4xBQcTZBhhtqebqx9vUu+7RxMGaRynMqVj/HgrelEJ0QTHR/NyYsnr7n/WlVqGcnLFY+UGlZriIPJ4Zr7r+3pRt+/ZvPCvK+MFV27wuLFRX6+kP+5lgZXXq+CzsGattfbHkr3NSvtCtvp1qqExcXFhdDQUNatW2dZNnLkSDZt2sT69esL3D6/hOVqmZmZtGrVio4dO/LZZ9eeuOxad1j8/PyUsIgIcOUdlqpEHT6b73+FX6+d2Wwm4XwCMQkxxusf4999p/dhJu//dVZ2qUxL75a57sRsPViF79bEMbxTQxwup9LpgTuodTbB2GDRIqPDbxGdb1m4W6A7LOVPsSQs/v7+dO3alW+//dayLDw8nA8++IC4uLgCty9swgIwdOhQjh07xoIFCwoVm4Y1i0hJOZd6jm0ntuUkMgkxbDuxjUvpeUc2OpocCageYElgQnadpeXIsVS/CAQHQ1QUOKgklpRfhf39tqoPS1hYWJ6RPHv37sXfv2izUrPZTExMDM2bNy/S/YqIFIUqrlXo4NeBDn4dLMvSM9PZe3ovMQkxRMdHE/OP8e/plNPsOLmDHSd3MHPbTKPxq1AnCYITYgj+4iGC7xpEiG8I9Tzr5XmkJCIGqxKWUaNG0aFDB8aOHUvfvn3ZuHEjkydPZvLkyZY2o0ePJi4ujunTp1uWxcTEAEbH2pMnTxITE4OLiwtBQUEAjBkzhnbt2tG4cWOSk5P57LPPiImJ4csvvyyCUxQRKX5ODk4E1QgiqEYQjzZ/FDD+4+v4ueN5HintP7OfYx5wzAP+d/YP+OUPANxd3fM8UmpaoymuTq62PDURu2B1pdv//e9/jB49mn379lG/fn1eeumlXKOEhgwZwqFDh1i5cmXOQUymPPvx9/fn0KFDgJEI/fbbbyQkJODh4UFISAjvvvsu7du3L3RceiQkIqVFcmoyWwd1I+bwBmJ8ICa0DttMJ7iccTlPWycHJwKrBxrDrLM6+Lb0aUk1N1XLlbKhWPqw2DMlLCJSqmzdavRhMZuhalXS9u5mT+YJ43HSFXdjzqScuebmdT3qGn1ifEIsd2P8PfwxAVzjPxJF7JUSFhERezdkCHz/vfH+X/+CDz/MtdpsNnMs+ZjlkVJ0gpHMHEw8eM3deVCB4OOZPHeiHg//EAWVKxfzCZRyly4ZUzI4O9s6knJNszWLiNixGZGH6VH1LlIds34sP/sMjhzJ1cZkMuHn4UfS2RYsXB9Gncz/wz3xa5oxm4fqTKEWw6lm7or/ZR+cMyCJS6yqdZlHWuxl9ONPMyPysA3O7MbMiDxM2PjlBcac3W5kRHSh2ucn46/VvNLfi5d7V8G8d2+hjy+2o4RFRMQGwlceYLdrVaaG3mcsSE01Svbn0zYuMYV5W4+TmJLGuRRXYvb50HNrddZ9uYdDYxM4PxZiwqHrAchwgI0Vf+Pr5ftK8IxuTvY5hq88UKh287YeL1T7a4qM5F8T7mJi8EU+bp3Klm8/KPTxxXaUsIiI2MDwTg3xdHNmZqf+pFbxMBb+8ANs2XLNtrU93ejVohaebs50OrmXxXP+j29/e58mp4y7Mi4Z4NXoTj4/2gynDFjeIJX7XJaU5CndlOxzHN6pYaHa9WpRq1Dt89i0ia9e7cwnoWmWRfMOLCj08cV21IdFRMTWPv4YXn7ZeN+9OyxceO12O3fCG2/AH3/kXt6+PUyYALffDnPm8MLkh/isHbS4WIW/x50tnhmoS6OoKBY8dQf33HeBTAcIjYPNtaH9UVg3ej80VLJiC+rDIiJSWowYAdkFOBctgqVLc68/dgyeegqaN8+drAQEwJw5sHatkawA3Hsvb++vhWcKbK14jmlL/lMy52DvYmLY0r8zfXsaycqQeB9+8xkJQGQdOPXnTzYOUAqihEVExNZcXeHf/875/OqrkJkJZ8/C669D48bw3XfGMoBatWDKFNi2DR54IPcwZicnvJ54jv9bbXx8a90HnL98vsROxS5t3crx+7twzz3nOO8Knc948M2EHfj1eYoWCWA2wcJNEbaOUgqghEVExB707w8hIcb76GgYNMh4RDFhgjH8FsDDA8aNg337jDsuTvkUK3/qKUbEONPwDCQ4XODDlf++drvyYPt2LnTvwr09znLMA5pccGP2W1tx8agGzZrRK6EKAPPSdsL5cp7Y2TklLCIi9sDBwUhOss2aZdxhAXBxMfq4HDhg3HGpWPH6+6pRA9c+/fkwq8/tR5EfczTpaPHEbc927SLjri4M6HSav2tB9ctOzH9uPVVr1DXWm0z08usCwMIGmaQvXWzDYKUgSlhEROxF167QrVvOZ5MJHnsM9u6Fjz4CL6/C7+u553hwF9x+GFLMl3lz+ZtFH68927MHunThX8En+SMAXDNM/DFwHg3qtszVrF3Xx6l2ERLdYP3SqTYKVgpDCYuIiI3NiDxM8JjFBI9ZzNv3vshfjdqwoFknxo37keD6jxI8bY/1Bc3atMHUti0fLzI+/rD1BzYf31z0wReB/Iq2WVtMztJu/37o0oWv/BL4JGtKumm9JtMhsFuuaz0j8jCOd3alx0Hjp3Be3EpjqgSxS0pYRERsLHzlARJT0khMSWNmXCYDe7/D8F6v8G1iJcvyGypo9txzhB6HQVmlXV5a9BL2WMkiv6Jt1haTC195AGJjoXNnFlQ8zvM9jfUftHuDR259ytI21zWtWJFers0BmOd7/pp1cMQ+KGEREbGx7CJynm7OluJwV7+/oYJmfftCjRr8exm4pcGaI2uYs3tO0Z/ATcqvaJu1xeReaeICnTuzJe0YfftgDF8OeIQ3un2Qq+3V17RH20dxyITt3nDkfzOL/gSlSKhwnIhIWfbWW/Dvf/N2Z3j/DmhYtSE7nt2Bq5OrrSMrWkeOwB13cPz0IW59Co55QOdaYSx8Yjkuji7X3/bQIW57vz5r68JXO+sz/KfYkolZABWOExERgGeeAUdHXl0LPhccOHD2AF9u+tLWURWtY8egc2cuxB3i3v5GstLEsxGzB/5ZcLICUK8evRJrAjDP6SCcOlXMAcuNUMIiIlKW+fnBAw9Q+TJ8sNQoPPf+6vc5ffG0jQMrIsePQ5cuZByMZUBvjOHLFaoxf/AiqrpVLfRuejUyOrwsrw8p8+cWV7RyE/RISESkrFu5Ejp3JsMErV6qyNYqFxnZdiSf9vzU1pHldukSJCYar6Ska7+/+vO+fXDyJC91h0/ag6ujK8sfW04Hvw5WHdq8ciV153bmmAfMO9qRu79dVbTnJvkq7O+3EhYRkbLObDbmIdqxg6UNoOtgcHJwYvvw7TSp3sS2sf3zD9x/P8TEQGrqDe3iqzYwopfxPqJ3BI80e8T6naSlMeyRynzT4jLPbnHhy58v5F9JWIqU+rCIiIjBZILnngPgrli4J8WP9Mx0Xl36qo0DAz78EDZsuOFkZUEzV56/23j/QecPbixZAXB2ppdnWwDm+V/GvG7dje1Hio3SRxGRMmhG5GHCVx6gtX9VVu89iVuqL8vdKuOWcp5/f5/AvGccmLtnLssPLqdL/S42iTFi5W56hU/BHcDZGdq3B0/PnJeHB1FJZhbFXSIstBFplaowa28y990RxAOdm7L10iH6Tu9E5uVUbq/1MPPWtsPfxSge99GiPQB0vKUGUYfPWoYwh688wPBODRnYzj9PPF1uH4zrgb847AnffvEx09al59tWSp4eCYmIlEFh45cTl5iCowkysv5f/p1lk3l8s9GhtP1TgUTW2UWwTzCbh27G0cGxxGP898Ov8ubs/xgfBg+G77/P0yb7PGp7ugFY3v8yIoBbv72VY8nH6FyvMynHXyM+KT1XO8By/ldvv/b1ayRpJ0/SY1RNFjWGV9Z58kvIjPzbSpHRIyERkXIsu5jalcXnjvV/3LJ+8rwzVHSqQkxCDNO3TLdJjM/svGKywWefvWabK4vHZb9//DZf7o24l2PJx2ji1YTZfWczonNArnZXFt+7evt8C9HVqEGnc7UB2FgrkVYk31jBPikWusMiIlKe9OgBi4wJhiZ++wSvHPsvvpV92fv8Xiq7VC65ODZtgrZGnxFatYLNm42+NgXIyMyg98+9+WPPH1SvWJ3IJyNpWK3okoqD771EA/MnOGbCyZofUXXEy0W2b7k23WEREZG8sjrfAjz361EaVG1A/Pl4/rP2PyUbR3h4zvvhwwuVrAD8a8m/+GPPH7g6uvLHI38UabICUP+egQSehAwHWLxBZfrtiRIWEZHypGdPqF8fANeFS5gQ9AIA/1n3H+KS40omhjNnICLCeO/hAf37F2qz8E3hfBL5CQDTHphmda2VQgkJodfxSgDMS9kKKSlFfwy5IUpYRETKE0fHXP1Fev/vAGF+YaSkp/Dm8jdLJobvvzeKxAEMGQKVKhW4yYJ9C3hugXF36KaGLxfEZKKX7x3GMetnkLF8afEcR6ymhEVEpLx54gmoUAEA09RpfHy7MZvx91u+J+p4VPEeOzMz7+OgAmz9Zyt9f+1LpjmTIcFDeOP2N4oxQAi7cwgel+BUJdi0eFqxHksKTwmLiEh5U60aDBhgvE9Opu2y3Qxobnx+efHLFOtYjGXLjHL6AF26QJPrV9o9fu44vWb14vzl83Su15lv7vkGUyH7u9wo56496B5rHGPekWVGpWCxOSUsIiLlyIzIw4SNX8682x+0LNv/9nhCPYfj7ODKqsOrGDV3SvEFcMXdlVV39blu0wuXL1iGL7uZ/OjXYFLhZl++SvY5z4g8XLgNqlShl1MgAPO8kxj70WyCxywmeMziwu9DipwSFhGRciR85QHiElMYG+/GFv9mADQ6eZj9P27A2/QwAN9seY/LGZeL/uDHjsEffwCQULka/5fZIN+mGZkZDPhtAH/H/40THlRNeZvpa0/d0GGzzzl85YFCb9OzdX9MZoj2hQtLfyExJY3ElDSr9iFFSwmLiEg5cmXxtOQnn7Ys/9f+pbzdaTTOVOWSOY6vNn1V9AefPNnowwLMu/Uenr4z/8dBry551TJ8eXTb/1LPo8ENF3ErsGDcNdS4tx9tswZNVXKKtBSiUyE521HhOBGR8uryZfD3h4QEY/TQwYN8e3IRQ/8cStUKVdk/cj/V3KoVzbHS0qBu3ZxjHT4MtWtfs2n4pnCenW+MZLrh2ZeLwPsPefF2yzM8sBvmfHUGqla1SRxlnQrHiYjI9bm4wDPPGO8zMuCbb3g8+HGa12zO2UtneX/V+0V3rN9/N5IVgAceyDdZKbHhy4XQq0EPAJY0gNSF82wWhxisTlji4uIYOHAgXl5eVKxYkeDgYKKi8h8GFx8fz6OPPkqTJk1wcHDgxRdfvGa72bNnExQUhKurK0FBQcyZM8fa0ERExFpPPw1OTsb7yZNxTEtnYreJAHyx6Qv2nt5bNMf56opHTPkMZS7p4csFCek+BN9zcMEFVq+cZtNYxMqE5ezZs4SFheHs7MyCBQvYuXMnEydOxNPTM99tUlNTqVGjBm+++SYtW7a8Zpv169fTr18/Bg0axJYtWxg0aBB9+/Zlw4YNVp2MiIhYqVYteNjobMvJk/DLL3Rt2JW7G99NemY6ry197eaPsWsXrFxpvG/SxBjOfJX4c/HcM+ueEh2+XBDTHXdw90EjmZt3Yq1xF0psxqo+LK+//jpr165lzZo1N3SwTp06ERwczKRJk3It79evH8nJySxYsMCyrEePHlStWpWI7PLNBVAfFhGRG7R2Ldx2m/H+1lshMpKdJ3fSIrwFGeYMVjy2gk71Ot34/keOhM8/N95PmgQvvJBr9YXLF7hj2h1ExUfRxKsJ659cT1U3++gvMmdIOx6qv4FGp2Hfo+uhXTtbh1TmFEsflrlz5xIaGkqfPn2oWbMmISEhTJly8+P1169fT7du3XIt6969O+vWrct3m9TUVJKTk3O9RETkBnToAMHBxvsNG3hj4Bj+PlCJTnUeBaDvj8NpOWYhwWMWMzIiOk9Nk+vWOTl/3ijFD6S7VqD7mfrMiDxs2Wb6+lgG/DaAqPgoqleszrxH591wsmJ1vZVCuCtsEM4ZsN8L9s6bXizHkMKxKmGJjY0lPDycxo0bs2jRIoYNG8bIkSOZPn36TQWRkJCAt7d3rmXe3t4kZHfQuoZx48bh4eFhefn5+d1UDCIi5ZbJBKNGWT6+9fM4ls5aSPI/D2AyV+Rk6m7iLi8mMSWNeVuP56lpct06J7NmQdZ/UC5o3pk9qY6Erzxg2eb1Za8V2ezLN1JvpSBV7nmIOw4Z7+ftmlssx5DCsSphyczMpFWrVowdO5aQkBCeeeYZhg4dSviV80LcoKufVZrN5us+vxw9ejRJSUmW19GjR286BhGRcmvQIHjUuKNSMS2VL398l9ea1aGu00AAklymU8Utg14tauWpaZJvnROzOVdlW8fnRljaDe/UEMcqS4jP+BUomtmXb6TeSoF8femVUgeAeZXieKlp5aI/hhSKkzWNfX19CQoKyrUsMDCQ2bNn31QQPj4+ee6mnDhxIs9dlyu5urri6up6U8cVEZEsJhN89x0cPAjr11PpZAK93xtBr2WLCJq6hIOJB3n4jmje6XRfnk0HtvNnYDv/vPuMjISYGON927bc/Vgv7s5atXD/Qg5nfAEU3fDlfOO4Sb0C72MUX7HaH+bEb6T3688V+TGkYFbdYQkLC2PPnj25lu3duxd//5v7grRv354lS5bkWrZ48WI6dLi5bFtERKxQoQLMmWMUeAPYvJkKTz7DhDvHAfDhug85fu544fd35VDmZ5+1vN36z1b6/tKXDHMGj7V8zObDlwvSuNdgGp+GNEdYsu4HW4dTblmVsIwaNYrIyEjGjh3L/v37mTVrFpMnT2bEiBGWNqNHj2bw4MG5touJiSEmJobz589z8uRJYmJi2Llzp2X9Cy+8wOLFi5kwYQK7d+9mwoQJLF26NN+aLSIiUky8veF//4PKlY3Pv/7Kwz9to4NfBy6mXeSt5W8Vbj+nTsHPPxvvq1WDvn2BnOHL5y6fo1O9Tky+d7LNhy8XqE0beh11A2Deub8hNdXGAZVTZiv9+eef5mbNmpldXV3NAQEB5smTJ+da/9hjj5nvuOOOXMuAPC9/f/9cbX755RdzkyZNzM7OzuaAgADz7NmzrYorKSnJDJiTkpKsPSUREbna//5nNjs4mM1GTxTz+m/fNfMuZtO7JvPfx/8uePsJEyzbml9+2Ww2m83nU8+bW3/T2sy7mJt83sR85uKZYj6JorPkmbvMvIvZ52XMGYsW2jqcMqWwv9+aS0hERK7tk0/gpZeM966uPPppRyISltC5XmeWDV6W/52RzExo1MjoDwOwbx8ZDerT++fe/LHnD6pXrE7kk5E3NSKopF3+aRZeWwZw3hU2Jz9C64mFqxEmBdNcQiIicnNefBGGDjXep6Yy7sO/cXVwYcWhFfy598/8t1u0KCdZ6d4dGjXKNfvyzQ5ftgWX7nfTNeuU5h1cbNtgyiklLCIikqcg2ozIw4RNWMHMwa9aSun7x57mpZ3GfwH3+fFZWoyZd+0Cald0tl3ZpTcN//0iH0d+DBTN8GWb8PSkF7cAMM/rDI+8/H2xFI9TYbr8KWEREZE8BdGyP3/11xH45Rdo3BiA1+ecouolJy4Tx7HLc/MWUDt0COZlzWzs58dzqSeJTTPK8tt69uWbdXeI0XF4U23w37GkWIrHqTBd/pSwiIhInqJruT5Xq2aMHKpaFfdUGL8kHYBkl1kM6uCVe0eTJxtdbYGtT9/PAdM4MGVyW63edj98uSC+9/Sn1XEwm6BC5l/FUjyuWIrflRHqdCsiIoWzfDl07056Zjohz8B2bxjVbhQfdzce95CaCn5+cPIk8Z6O3PqWN0fPH6dTvU4sGrgIF0cX28Z/s8xm3n6oKu8HJ9Fnp4mfv02CKlVsHVWpp063IiJStLp0ga++wikTJmb1O/1iw+fsP7Pf+DB7Npw8yQVnuPcZd46eP04Tryb81ve30p+sAJhM9PK/C4BFDcykLV5g44DKFyUsIiJSeEOHwqhRdDsAPfZBmjmd1+ZkFQ/96isyTDCgN0S5nb3p2ZftUZtuj1PjAiRXgLXLptk6nHJFCYuIiFjnP/+BXr34aDE4ZMJvxxaz+vdPYe1aXu0KfwRQaocvF8Shcxd6xjoCMC9+lVFzRkqEEhYREbGOoyNERNDUpzlPRxmLXlo4ii/bwMdZI5ZL7fDlgri50atiSwDm1rlI8pD+8NlnsGKFMR2BFBt1uhURkRtz+DAn7gil0aOnOOeas/iDDm/xZtf3bRdXMUv86mO841/mshO4pUGfHfBkNNx+GEy+vtC8ObRoYfzbvDkEBhoTS8o1Ffb3WwmLiIjcuPXrmTD6dl7vnAHAg2fqk+AzhdB61Vi99yQAr3RvAhg1RoZ3asjAdv42C7dInD/PtL6t+LDBPnbVyFnc6DQ8HgOPxUDtcznL0x0cOOXrz45q/ni1b01iwyb8J7ka/e8JBeCjRXtITc/A1ckx17Vq7V+VqMNny8Y1uw4lLCIiUiIuzfyehxY+jke6Ixm1vmajsw+OJsjI+nWp7WnMdByXmEJtTzfWvt7FhtEWjbDxyzn7zykaZ27Au84mfjwfyTmHNMDo19NjPzwRDffuBZeMvNtfcK7Ae4+8yV9BHYhLTLEsv/JaZV/DsnLN8qOERURESs7Ro+DqyozYFMvdgTJ7hwWjhP6V53Ph8gV+3fkr41Z/xZ6zGy3tPC+78MBuN15Yf4Hg+PRc+0hzq8iiybN564CD7rAoYRERESlZ+07vY2rMVKbFTCP+fLxledtqzXnCuS2PzDuMx/+WGgv9/WHTJqhRI5+9lX1KWERERGwoPTOdxQcW8130d8zdM5f0TOMOi5uTG70PufHk0jN0PAwOt90OS5eCSxkorncDlLCIiIjYiZMXTjJj6wy+i/6OHSd3WJY3OGP0dXksqD91wmeCyWTDKG1DCYuIiIidMZvNbDq+if9G/5eILTNITr8AGB11u7kE8uSD73HvLffi6uRawJ7KDiUsIiIiduxi2kVmf/cy3234mlX1cpZ7uXkxsMVAngh5ghbeLWwWX0lRwiIiIlIavPEG+78Zx7RgmBZiIq5Kzs9yaK1Qngh+gv7N++NZwdNmIRYnJSwiIiIlIHuIc0HDkK8cCg1XDPNu6wcPPghz55JhgsUda/Pfoa35I3YBaZlGbZcKThXoHdibJ0KeoFO9TjiYinZmnauHaZckJSwiIiIlIGz88kIVestud81CeufOQYcOsH270bhbN0798j0zdv7Id9Hfsf3Edst+6nvW5/Hgx3ks+DHqetQt0nOwRZG6wv5+a/JDERGRmzC8U0Nqe7rRq0Utanu6We6g5NdueKeGud4DUKUKzJ0L1asbnxcvpvo7E3ix3YtsHbaVTUM3MTx0OB6uHhxMPMjbK9+m3qR6dJ/RnZ93/ExqemqRnEN+sdsD3WERERGxF6tXw513QnpWVdxvv4Unn7Ssvph2kTm75vBd9HesOLTCsryaWzUGNB/AkyFP0tKnZUlHfVP0SEhERKQ0mjIFnn7aeO/sDMuXw2235WkWezaWaTHTmBozlWPJxyzLW/m24ongJ3i0+aNUdataUlHfMCUsIiIipdXIkfD558b7GjVg40aoV++aTTMyM1gau5T/xvyX33f/zuWMywC4OrryUOBDPBHyBF3qdynyjrpFRQmLiIhIaZWeDj17GiX7AVq0gLVroXLl6252+uJpZm6byXfR37H1n62W5f4e/jwe/DhDgofg72lfEykqYRERESnNzp6FW2+FffuMzw88ALNng0PBd0rMZjPRCdF89/d3zNo+i8RLiQCYMHFXg7t4IuQJHgh4gApOFYov/kJSwiIiIlLa7d4N7dpBUpLx+a234P33rdpFSloKv+/+ne+iv2PZwWWW5Z4VPC0ddUN8Q4oyaqsoYRERESlm1yu4VhTF2GZEHiZq8o9M/P4NHDIzAXj7kTe55YWhAFYf+1DiIUtH3SNJRyzLg32Cae75AH/vbooT7nS8pcZ1i+AVJSUsIiIixex6BdeKohhb9j5e2j6PkfPCAbjk5MLkbk8wp919HLyQeUPHzjRnsix2Gf+N+S9zds0hNSOrjovZiYoZ7XHP7IpLRkvqeFYu9kJyKhwnIiJSzK5XcK0oirFl76Pam6/CE08AUCH9MiPnf828z4bwws6FjGhf2+pjO5gc6NqwKxG9Izj+8nE+7/k5/lWagimdi05rSHB5mwS3odSs8wcHzx40NkpNhay7PLagOywiIiKlQWoqDBsG06blXl6nDrz5ppHQuLjc1CGi46P5b/R/mbltJmcvnbUsvzO1Nk+sSOTBcb/j1umumzrG1XSHRUREpCxxdYWpU2HrVujdO2f5sWMwfDg0bmxUxk1Lu+FDhPiG8Pndn3P83pX8mPEgXeMqYDLDMtc4BvS4wB9zxhXBidwYqxOWuLg4Bg4ciJeXFxUrViQ4OJioqKjrbrNq1Spat25NhQoVaNCgAV9//XWu9dOmTcNkMuV5Xbp0ydrwREREyrbmzeHXXyE6Gu67L2f5kSMwdCgEBMD33+eU9y+sEyfgs8+gTRsqNG1Jv/fnsHjKJQ5OgjErIOQfEw9kNC7SU7GGVQnL2bNnCQsLw9nZmQULFrBz504mTpyIp6dnvtscPHiQu+++m9tvv53o6GjeeOMNRo4cyezZs3O1c3d3Jz4+PterQgXbjw8XERGxS8HB8McfsGkT3H13zvLYWBgyBIKCYOZMyMjIfx8XL0JEBPTqBbVqwQsvwObNOesdHPBv1523n/qBvz9MpsIXX+e/r2JmVR+W119/nbVr17JmzZpCH+C1115j7ty57Nq1y7Js2LBhbNmyhfXr1wPGHZYXX3yRxMTEwkd+FfVhERGRci0yEt55BxYvzr08IADefRf69DGKzmVkwIoVMGOGUYju/Pm8+2rVCgYOhEceAV/fYg27WPqwzJ07l9DQUPr06UPNmjUJCQlhypQp191m/fr1dOvWLdey7t27s3nzZtKueM52/vx5/P39qVOnDvfccw/R0dHWhCYiIlK+tWsHixYZMz537pyzfPduI/Fo2dKYo6huXeja1XhsdGWyUrcujB4NO3ZAVBSMGlXsyYo1nKxpHBsbS3h4OC+99BJvvPEGGzduZOTIkbi6ujJ48OBrbpOQkIC3t3euZd7e3qSnp3Pq1Cl8fX0JCAhg2rRpNG/enOTkZD799FPCwsLYsmULjRtf+3lZamoqqampls/JycnWnIqIiEipkF0ErrV/VUsxN8gpGpf9Pnt9a/+qRHV7i3d7D6XllE+ouWWTsaPt243XFc65VmJJ09up/OQQug3rk6vs/5XHXb33JACvdG9S7IXk8mPVIyEXFxdCQ0NZt26dZdnIkSPZtGmT5fHO1W655RYef/xxRo8ebVm2du1abrvtNuLj4/Hx8cmzTWZmJq1ataJjx4589tln19zvu+++y5gxY/Is1yMhEREpS7KLwDmaIMMMtT3dACyF4bLfZ6/P1c5spt6WSF5e8wOtju8B4LKDE5EBbVkUfBe/+gaT6uRy3eJz2fsDbqoIXn4K+0jIqjssvr6+BAUF5VoWGBiYpwPtlXx8fEhISMi17MSJEzg5OeHl5XXNbRwcHGjTpg37sid8uobRo0fz0ksvWT4nJyfj5+dXmNMQEREpNYZ3amj9HZYr2n10OYMnAkKZ4HsOt9MnGZfiw4BerQgE3BbtwS3rGNc7bvYdlpspgnezrEpYwsLC2LNnT65le/fuxd8//9tD7du3588//8y1bPHixYSGhuLs7HzNbcxmMzExMTRv3jzf/bq6uuLq6mpF9CIiIqXPwHb+13wMc+Wy6z2muXpdx+usK8xxbcWqTrejRo0iMjKSsWPHsn//fmbNmsXkyZMZMWKEpc3o0aNz9WcZNmwYhw8f5qWXXmLXrl3897//5bvvvuOVV16xtBkzZgyLFi0iNjaWmJgYnnzySWJiYhg2bFgRnKKIiIiUdlbdYWnTpg1z5sxh9OjRvPfee9SvX59JkyYxYMAAS5v4+HiOHMmZAbJ+/frMnz+fUaNG8eWXX1KrVi0+++wzel9RpS8xMZGnn36ahIQEPDw8CAkJYfXq1bRt27YITlFERERKO80lJCIiIjajuYRERESkzFDCIiIiUsrNiDxM2PjlzIg8XGCbkRHRBba1R0pYRERESrnwlQeIS0whfOWBAtvM23q8wLb2SAmLiIhIKTe8U0Nqe7pdt05KdpteLWoV2NYeqdOtiIiI2Iw63YqIiEiZoYRFRERE7J4SFhEREbF7SlhERETE7ilhEREREbunhEVERETsnhIWERERsXtKWERERMTuKWERERERu6eERUREROyeEhYRERGxe0pYRERExO4pYRERERG752TrAIpK9qTTycnJNo5ERERECiv7dzv7dzw/ZSZhOXfuHAB+fn42jkRERESsde7cOTw8PPJdbzIXlNKUEpmZmRw/fpwqVapgMplueD/Jycn4+flx9OhR3N3dizDC0k/XJn+6NvnTtcmfrs316frkryxdG7PZzLlz56hVqxYODvn3VCkzd1gcHByoU6dOke3P3d291H8JiouuTf50bfKna5M/XZvr0/XJX1m5Nte7s5JNnW5FRETE7ilhEREREbunhOUqrq6uvPPOO7i6uto6FLuja5M/XZv86drkT9fm+nR98lcer02Z6XQrIiIiZZfusIiIiIjdU8IiIiIidk8Ji4iIiNg9JSwiIiJi95SwXOGrr76ifv36VKhQgdatW7NmzRpbh1Ti3n33XUwmU66Xj4+PZb3ZbObdd9+lVq1auLm50alTJ3bs2GHDiIvP6tWruffee6lVqxYmk4nff/891/rCXIvU1FSef/55qlevTqVKlbjvvvs4duxYCZ5F8Sno+gwZMiTPd6ldu3a52pTF6zNu3DjatGlDlSpVqFmzJg888AB79uzJ1aa8fncKc23K6/cmPDycFi1aWArBtW/fngULFljWl9fvzJWUsGT56aefePHFF3nzzTeJjo7m9ttvp2fPnhw5csTWoZW4pk2bEh8fb3lt27bNsu7DDz/k448/5osvvmDTpk34+PjQtWtXy1xOZcmFCxdo2bIlX3zxxTXXF+ZavPjii8yZM4cff/yRv/76i/Pnz3PPPfeQkZFRUqdRbAq6PgA9evTI9V2aP39+rvVl8fqsWrWKESNGEBkZyZIlS0hPT6dbt25cuHDB0qa8fncKc22gfH5v6tSpw/jx49m8eTObN2+mS5cu3H///ZakpLx+Z3Ixi9lsNpvbtm1rHjZsWK5lAQEB5tdff91GEdnGO++8Y27ZsuU112VmZpp9fHzM48ePtyy7dOmS2cPDw/z111+XUIS2AZjnzJlj+VyYa5GYmGh2dnY2//jjj5Y2cXFxZgcHB/PChQtLLPaScPX1MZvN5scee8x8//3357tNebk+J06cMAPmVatWmc1mfXeudPW1MZv1vblS1apVzd9++62+M1l0hwW4fPkyUVFRdOvWLdfybt26sW7dOhtFZTv79u2jVq1a1K9fn0ceeYTY2FgADh48SEJCQq7r5Orqyh133FHurlNhrkVUVBRpaWm52tSqVYtmzZqVm+u1cuVKatasyS233MLQoUM5ceKEZV15uT5JSUkAVKtWDdB350pXX5ts5f17k5GRwY8//siFCxdo3769vjNZlLAAp06dIiMjA29v71zLvb29SUhIsFFUtnHrrbcyffp0Fi1axJQpU0hISKBDhw6cPn3aci10nSjUtUhISMDFxYWqVavm26Ys69mzJzNnzmT58uVMnDiRTZs20aVLF1JTU4HycX3MZjMvvfQSt912G82aNQP03cl2rWsD5ft7s23bNipXroyrqyvDhg1jzpw5BAUF6TuTpczM1lwUTCZTrs9msznPsrKuZ8+elvfNmzenffv2NGzYkO+//97S8U3XKceNXIvycr369etned+sWTNCQ0Px9/dn3rx5PPTQQ/luV5auz3PPPcfWrVv566+/8qwr79+d/K5Nef7eNGnShJiYGBITE5k9ezaPPfYYq1atsqwv798Z3WEBqlevjqOjY54s9MSJE3ky2vKmUqVKNG/enH379llGC+k6Uahr4ePjw+XLlzl79my+bcoTX19f/P392bdvH1D2r8/zzz/P3LlzWbFiBXXq1LEs13cn/2tzLeXpe+Pi4kKjRo0IDQ1l3LhxtGzZkk8//VTfmSxKWDC+JK1bt2bJkiW5li9ZsoQOHTrYKCr7kJqayq5du/D19aV+/fr4+Pjkuk6XL19m1apV5e46FeZatG7dGmdn51xt4uPj2b59e7m7XgCnT5/m6NGj+Pr6AmX3+pjNZp577jl+++03li9fTv369XOtL8/fnYKuzbWUl+/NtZjNZlJTU8v1dyYXG3T0tUs//vij2dnZ2fzdd9+Zd+7caX7xxRfNlSpVMh86dMjWoZWol19+2bxy5UpzbGysOTIy0nzPPfeYq1SpYrkO48ePN3t4eJh/++0387Zt28z9+/c3+/r6mpOTk20cedE7d+6cOTo62hwdHW0GzB9//LE5OjrafPjwYbPZXLhrMWzYMHOdOnXMS5cuNf/999/mLl26mFu2bGlOT0+31WkVmetdn3Pnzplffvll87p168wHDx40r1ixwty+fXtz7dq1y/z1GT58uNnDw8O8cuVKc3x8vOV18eJFS5vy+t0p6NqU5+/N6NGjzatXrzYfPHjQvHXrVvMbb7xhdnBwMC9evNhsNpff78yVlLBc4csvvzT7+/ubXVxczK1atco11K686Nevn9nX19fs7OxsrlWrlvmhhx4y79ixw7I+MzPT/M4775h9fHzMrq6u5o4dO5q3bdtmw4iLz4oVK8xAntdjjz1mNpsLdy1SUlLMzz33nLlatWpmNzc38z333GM+cuSIDc6m6F3v+ly8eNHcrVs3c40aNczOzs7munXrmh977LE8514Wr8+1rglgnjp1qqVNef3uFHRtyvP35oknnrD8/tSoUcN85513WpIVs7n8fmeuZDKbzeaSu58jIiIiYj31YRERERG7p4RFRERE7J4SFhEREbF7SlhERETE7ilhEREREbunhEVERETsnhIWERERsXtKWERERMTuKWERERERu6eERUREROyeEhYRERGxe0pYRERExO79PwIzn6qj/DsNAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "df = pl.read_csv(io.StringIO(data), has_header=False, separator=';', new_columns=['x', 'y'])\n",
    "x, y = df.to_numpy().T\n",
    "x2, y2 = segments_fit(x, y, 30)\n",
    "\n",
    "scale_y = 1000\n",
    "\n",
    "line = shapely.LineString(np.c_[x2, y2 * scale_y])\n",
    "line = line.simplify(20)\n",
    "x3, y3 = shapely.get_coordinates(line).T\n",
    "y3 /= scale_y\n",
    "plt.scatter(x, y, s=1)\n",
    "# plt.scatter(x2, y2, color='r');\n",
    "plt.plot(x2, y2, 'r', lw=2)\n",
    "plt.plot(x3, y3, 'g');"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}

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

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)

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

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?

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!

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!

@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.

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

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?

Yes, it should use the natural log.

…

On Mon, Feb 20, 2023 at 4:58 PM Qinkong ***@***.***> wrote: ***@***.**** commented on this gist. ------------------------------ 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) Thanks 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? — Reply to this email directly, view it on GitHub <https://gist.github.com/70effade57483355bbd18b31dc370f2a#gistcomment-4477503> or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAYUH7OJLOUYQZRZ5DFMF6LWYPSH7BFKMF2HI4TJMJ2XIZLTSKBKK5TBNR2WLJDHNFZXJJDOMFWWLK3UNBZGKYLEL52HS4DFQKSXMYLMOVS2I5DSOVS2I3TBNVS3W5DIOJSWCZC7OBQXE5DJMNUXAYLOORPWCY3UNF3GS5DZVRZXKYTKMVRXIX3UPFYGLK2HNFZXIQ3PNVWWK3TUUZ2G64DJMNZZDAVEOR4XAZNEM5UXG5FFOZQWY5LFVA2DONZZGY4TGOFHORZGSZ3HMVZKMY3SMVQXIZI> . You are receiving this email because you commented on the thread. Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub> .

-- Associate Professor Baker Institute for Animal Health College of Veterinary Medicine Cornell University Website: http://www.dankolab.org E-mail: ***@***.***

ruoyu0088 Thanks for the original version - it was exactly what I needed for a piecewise linear approximation of a CDF I wanted.

First off, thanks a bunch for this, @ruoyu0088 ! One of the fastest solvers I've found outthere. One small hiccup though, as you can see, the segment 2 is going back in time, which is a big no-no in the timeseries world and I could not figure out how to fix that

@josesydor

Please try the following code:

from scipy import optimize

def segments_fit(X, Y, count):
    xmin = X.min()
    xmax = X.max()
    ymin = Y.min()
    ymax = Y.max()    

    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.01].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):
        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='SLSQP', 
        bounds=[(0, None)] * len(seg) + [(ymin, ymax)] * len(py_init), 
        constraints=[
            optimize.LinearConstraint([1] * len(seg) + [0] * len(py_init), 0, xmax - xmin)
        ])
    
    return func(r.x)

@ruoyu0088 Thanks, but unfortunatelly it didnt work:

These are my X, Y:
3;6.37
4;6.38
5;6.38
11;6.42
12;6.37
14;6.33
15;6.34
16;6.35
17;6.35
19;6.37
20;6.4
21;6.41
22;6.39
23;6.37
24;6.39
25;6.4
26;6.35
27;6.32
28;6.31
30;6.31
31;6.3
32;6.3
33;6.29
34;6.29
35;6.27
36;6.28
37;6.28
38;6.29
40;6.34
41;6.33
42;6.33
43;6.35
45;6.36
46;6.35
47;6.36
49;6.33
50;6.35
51;6.35
52;6.35
54;6.35
55;6.34
56;6.34
57;6.35
58;6.35
59;6.34
60;6.33
61;6.32
62;6.33
63;6.32
64;6.31
65;6.31
66;6.31
67;6.31
69;6.32
70;6.3
72;6.29
73;6.29
74;6.29
75;6.3
76;6.28
77;6.28
78;6.29
79;6.29
80;6.3
81;6.31
82;6.32
84;6.32
86;6.33
87;6.33
88;6.35
90;6.33
92;6.34
93;6.33
94;6.32
96;6.33
97;6.33
98;6.34
99;6.33
101;6.31
102;6.32
103;6.32
105;6.32
106;6.33
108;6.33
109;6.31
110;6.31
111;6.34
112;6.35
114;6.35
116;6.34
117;6.35
118;6.35
119;6.34
120;6.34
121;6.33
123;6.32
124;6.32
125;6.32
126;6.33
128;6.32
129;6.33
130;6.33
131;6.32
132;6.33
133;6.33
134;6.32
135;6.32
136;6.32
137;6.32
138;6.32
139;6.32
141;6.32
143;6.33
144;6.33
145;6.32
146;6.32
147;6.32
148;6.3
149;6.28
150;6.29
151;6.28
152;6.28
153;6.27
154;6.26
155;6.26
156;6.28
157;6.25
158;6.25
159;6.24
160;6.24
161;6.22
162;6.21
163;6.2
164;6.18
165;6.18
166;6.18
168;6.17
169;6.17
170;6.16
171;6.15
172;6.15
173;6.18
174;6.16
175;6.16
176;6.16
177;6.17
178;6.15
179;6.15
180;6.16
181;6.16
182;6.16
183;6.16
184;6.16
185;6.16
186;6.17
187;6.16
188;6.16
189;6.16
190;6.16
191;6.16
192;6.16
193;6.17
194;6.16
195;6.15
196;6.16
197;6.15
198;6.16
199;6.14
200;6.14
201;6.13
202;6.13
203;6.13
204;6.13
205;6.13
206;6.13
207;6.14
208;6.14
209;6.14
210;6.14
211;6.14
212;6.14
213;6.15
214;6.15
215;6.14
216;6.12
217;6.12
218;6.12
219;6.12
220;6.11
221;6.12
222;6.11
223;6.11
224;6.09
225;6.11
226;6.1
227;6.09
228;6.1
229;6.1
230;6.1
231;6.1
232;6.09
233;6.09
235;6.1
237;6.08
238;6.08
239;6.08
240;6.09
241;6.08
242;6.1
243;6.11
244;6.1
245;6.12
246;6.13
247;6.14
248;6.15
250;6.13
251;6.12
252;6.11
253;6.12
254;6.14
255;6.14
256;6.15
257;6.13
258;6.12
259;6.1
260;6.12
261;6.12
262;6.12
263;6.11
264;6.12
265;6.13
266;6.13
267;6.15
268;6.15
269;6.14
270;6.16
271;6.12
272;6.12
273;6.11
274;6.11
275;6.11
276;6.1
277;6.09
278;6.08
279;6.06
280;6.07
281;6.07
282;6.08
283;6.08
284;6.06
285;6.08
286;6.08
287;6.07
288;6.08
289;6.07
290;6.08
291;6.08
292;6.08
293;6.09
294;6.08
295;6.08
296;6.06
297;6.05
298;6.05
299;6.06
300;6.05
301;6.06
302;6.05
303;6.04
304;6.05
305;6.05
306;6.06
307;6.05
308;6.05
309;6.04
310;6.05
311;6.04
312;6.05
313;6.05
314;6.07
315;6.06
316;6.06
317;6.06
318;6.05
319;6.05
320;6.05

This is the function original version. Again, segment 2 is going backwards:

@josesydor I uploaded another notebook to process your data, it looks ok.