saccadic · January 22, 2021 10:11
diff --git a/polyfit.ipynb b/polyfit.ipynb
 {
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "polyfit",
      "provenance": [],
      "collapsed_sections": [],
      "authorship_tag": "ABX9TyNlTrSpvUDzQrcjq22dL/Ny",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/gist/saccadic/8a369b5dff3ea2f2de03e7401a04d9a3/polyfit.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4r0Bcs_mKyHe"
      },
      "source": [
        "import cv2\r\n",
        "import numpy as np\r\n",
        "import pandas as pd\r\n",
        "import matplotlib.pyplot as plt"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "oDqu0XEwK3Ku"
      },
      "source": [
        "def powerXArray(a, order):\r\n",
        "    data = []\r\n",
        "    for i in range(len(a)):\r\n",
        "      data.append(np.power(a[i], order))\r\n",
        "    return np.array(data)\r\n",
        "\r\n",
        "def powerXYArray(a, b, order):\r\n",
        "    data = []\r\n",
        "    for i in range(len(a)):\r\n",
        "        data.append(np.power(a[i],order) * b[i])\r\n",
        "    return np.array(data)\r\n",
        "\r\n",
        "def MakeLeastSquaresMatrix(x, y, order):\r\n",
        "    x_order = []\r\n",
        "    xy_order = []\r\n",
        "    \r\n",
        "    for i in range(order*2+1):\r\n",
        "        data_x = powerXArray(x,i)\r\n",
        "        data_xy = powerXYArray(x,y,i)\r\n",
        "        x_order.append(data_x.sum())\r\n",
        "        xy_order.append(data_xy.sum())       \r\n",
        "    \r\n",
        "    # AX = B\r\n",
        "    # X = A^-1 * B\r\n",
        "\r\n",
        "    matA = []\r\n",
        "    for i in range(order+1):\r\n",
        "        mat_x =[]\r\n",
        "        for j in range(order+1):\r\n",
        "            if (i == 0) and (j == 0):\r\n",
        "                mat_x.append(1)\r\n",
        "            else:\r\n",
        "                mat_x.append(x_order[j+i])\r\n",
        "        matA.append(mat_x)\r\n",
        "        \r\n",
        "    matB = []\r\n",
        "    for i in range(order+1):\r\n",
        "        matB.append(xy_order[i])\r\n",
        "    \r\n",
        "    return np.array(matA), np.array(matB).T\r\n",
        "\r\n",
        "def model(x, param):\r\n",
        "    data = []\r\n",
        "    for i in range(len(x)):\r\n",
        "        y = param[0]\r\n",
        "        for k in range(1,len(param)):\r\n",
        "            y = y + param[k] * np.power(x[i],k)\r\n",
        "        data.append(y)\r\n",
        "    return np.array(data)\r\n",
        "\r\n",
        "def polyfit(x, y, order):\r\n",
        "  A, B = MakeLeastSquaresMatrix(x, y, order)\r\n",
        "  # print(A)\r\n",
        "  # print(B)\r\n",
        "  \r\n",
        "  # result, param = cv2.solve(A, B,cv2.DECOMP_QR)\r\n",
        "  # if result == False:\r\n",
        "  #   return -1\r\n",
        "\r\n",
        "  A_inv = np.linalg.inv(A)\r\n",
        "  param = np.dot(A_inv, B)\r\n",
        "\r\n",
        "  return param"
      ],
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 686
        },
        "id": "Op9ydXllKzaK",
        "outputId": "6b1a9fb9-193d-4744-9d28-f7561f215a1b"
      },
      "source": [
        "pt = np.array([\r\n",
        "    [0.0,0.0],\r\n",
        "    [0.1,0.2],\r\n",
        "    [0.2,0.4],\r\n",
        "    [0.3,0.6],\r\n",
        "    [0.4,0.7],\r\n",
        "    [0.41,0.71],\r\n",
        "    [0.42,0.72],\r\n",
        "    [0.46,0.73],\r\n",
        "    [0.5,0.8],\r\n",
        "    [0.51,0.82],\r\n",
        "    [0.52,0.83],\r\n",
        "    [0.53,0.87],\r\n",
        "    [0.6,0.9],\r\n",
        "    [0.7,0.91],\r\n",
        "    [0.8,0.93],\r\n",
        "    [0.9,0.96],\r\n",
        "    [1.0,1.0]\r\n",
        "],dtype=np.double)\r\n",
        "\r\n",
        "data = pd.DataFrame.from_dict({\r\n",
        "    'x': pt.T[0],\r\n",
        "    'y': pt.T[1]\r\n",
        "})\r\n",
        "\r\n",
        "print(data)\r\n",
        "\r\n",
        "print(\"-----\")\r\n",
        "\r\n",
        "order = 5\r\n",
        "param = polyfit(data[\"x\"],data[\"y\"],order)\r\n",
        "param = param.T\r\n",
        "#param_np = np.polyfit(data[\"x\"], data[\"y\"], deg=order)\r\n",
        "\r\n",
        "print(\"param\")\r\n",
        "print(param.T)\r\n",
        "# print(\"param_np\")\r\n",
        "# print(param_np)\r\n",
        "\r\n",
        "print(\"-----\")\r\n",
        "\r\n",
        "x = np.array(range(1000)) / 1000\r\n",
        "y1 = model(x, param)\r\n",
        "# y2 = model(x, param_np)\r\n",
        "\r\n",
        "plt.clf()\r\n",
        "plt.plot(data[\"x\"],data[\"y\"], c=\"k\")\r\n",
        "plt.plot(x,y1, c=\"r\")\r\n",
        "# plt.plot(x,y2, c=\"g\")\r\n",
        "plt.show()"
      ],
      "execution_count": 3,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "       x     y\n",
            "0   0.00  0.00\n",
            "1   0.10  0.20\n",
            "2   0.20  0.40\n",
            "3   0.30  0.60\n",
            "4   0.40  0.70\n",
            "5   0.41  0.71\n",
            "6   0.42  0.72\n",
            "7   0.46  0.73\n",
            "8   0.50  0.80\n",
            "9   0.51  0.82\n",
            "10  0.52  0.83\n",
            "11  0.53  0.87\n",
            "12  0.60  0.90\n",
            "13  0.70  0.91\n",
            "14  0.80  0.93\n",
            "15  0.90  0.96\n",
            "16  1.00  1.00\n",
            "-----\n",
            "param\n",
            "[ 2.35978781e-04  2.21484544e+00 -1.13490566e+00  8.40614535e-01\n",
            " -2.89088057e+00  1.97100168e+00]\n",
            "-----\n"
          ],
          "name": "stdout"
        },
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hUZd7G8e8vCRA6kaJIjyBK0QUDyIKKijSlKAgEwQLKqgu2XUXWgmJb1sX2WlZgWVcEKbYFpQhBwAJCQEQIgpgAElrC0iGF5Hn/OJENiBJgMiczuT/XNVcyMyeZ+5jk9uGZc55jzjlERCT0RfgdQEREAkOFLiISJlToIiJhQoUuIhImVOgiImEiyq8XrlKliqtbt65fLy8iEpKWL1+e7pyreqLnfCv0unXrkpiY6NfLi4iEJDPb9GvPacpFRCRMqNBFRMKECl1EJEyo0EVEwoQKXUQkTJy00M1svJntNLPVv/K8mdkrZrbBzFaZWfPAxxQRkZMpyAj9LaDTbzzfGWiQdxsMvHHmsURE5FSdtNCdc4uA//7GJt2Bt51nCVDJzKoHKqCISDjIyclhzcSJLL3sMr779ttCeY1AnFhUA/gp3/0teY9tO35DMxuMN4qndu3aAXhpEZGiyTnH2rVrSZg3j92TJ3P5smW0O3KEfcD0adNoevHFAX/NoJ4p6pwbA4wBiIuL05U1RCSsbNy4kYSEBObPn8/X8+Zx5c6d3AM0BXaXKcPKG2+k+ogR9G/YsFBePxCFngrUyne/Zt5jIiJhbceOHcyfP5/58+eTkJBASkoKLYB7Spfmn9nZRAOZF14Iw4cT06cPMSVLFmqeQBT6dGCImU0GWgF7nXO/mG4REQl1e/bsYeHChUcLfM2aNQA0L1eOp2rWpFP16lTetg1nht18M9xxB6VatQKzoOQ7aaGb2btAO6CKmW0BRgAlAJxz/wBmAl2ADcAh4LbCCisiEkyHDh3iyy+/PFrgy5cvJzc3l9LR0dzUvDmvdexIi5QUyqxfD99/D23bwogRWHw8VKgQ9LwnLXTnXPxJnnfAHwOWSETEJ9nZ2SxbtoyEhAQSEhJYvHgxWVlZREVFcWVcHO/06kW7w4c5Z+VK7KuvvC9q3RpefBF69YKaNX3N79vyuSIifsvNzeXbb789Og++aNEiDhw4gJlxZdOmvNmlC1dERlJ782YiExNhyRIoXx7at4dHH4XOnaFWrZO/UJCo0EWk2NizZw8rV65kxYoVLF68mM8++4xdu3ZRHuheuzbjL76YliVKUGPrVqJWrYJVq6BkSWjRAv78Z+jYEX7/e++xIkiFLiJhKS0tjRUrVhxzS09OpiHQELi8QgUeionhgshIyu/cCZs3e7ezz4ZLLoFbb/XmxFu0gOhon/emYFToIhLSnHOkpqayYsUKVi1dypbFi9mzahVl0tOpDdQB/hwdzflATP4vPHjQK++4OGjW7H+36qF7orsKXUSKNudg717YsQO3bRs7v/uOrStXsuf778nctAlLSyMmK4sWQLfjv9QMV706EfXqwfnnwwUXQMOG3i02tshOnZwuFbqIBF9ODuzaBTt3HnvbsePozW3fzpEtW4hITyfyyBEADDg773YE2B0VxeEKFaBaNSLq1yfrd7+jZP36UKcO1KmD1aiBhVlp/xYVuogEhnOwZw9s2QKpqd5tx45jy/rnz9PTITf3F98iNyKC/dHR7AA2ZmSwNTeXHcCuyEii69ThrEaNqHnJJTRo25YL2rShaunSQd/NokyFLiIFk5EBmzZBcrJ327z5f8X9c4kfPvzLrytfHqpVw519Npk1a7LnvPNIMyM1O5uUgwf5ZutWvk5JYUt2NrtzcykbEUGzZs1o3rw5zZs3p0uzZlxwwQWUKFEi+PscYlToIvI/mZnwww+QlATr1v2vvJOTvcJ2+dbUK1kSzj3XO5nmkkugWzcyKldme2Qkm3Jy+OHQIZJ27WLd5s2kpKSw8ZtvOHxc4Z9zzjk0btyYjt27Hy3w+vXrExGhi6mdDhW6SHF05Ih3qvq333rlnZQEa9fChg3e/PbPzj0XzjsPrr4aYmM5Urs220qXZkNODuv27iVl40ZSUlJI3rCBlLlz+e9/j710QoUKFahXrx4NGzakc+fO1KtX7+itbt26lClTJsg7Ht5U6CLhLisL1qyBFSu82/LlXpFnZHjPR0VB/frQuDGuVy92V6/OxjJlWJuTw4bUVK+wk5NJSUggNTUVl2+UXrJkSerUqUNsbCwtWrQ4prBjY2OJiYnBgrQwlajQRcLP7t3w1Ve4L74gZ9EiIpYvJyIzE4Ds0qVJr12bra1akVKpEuvLlmW9c2z/739JWb2aTZ98QmbetgBmxrnnnktsbCxXXXXVMWVdr149zj33XE2PFCEqdJEiLDs7m71797J79252797Nnj17jvm4e/dusrdvp9b69TT46Scu3LWL8/JG3keAb4Avga+B5cCPhw/j1q2DdeswMypWrEhMTAyVK1fmoosuonv37scUdp06dShVqpR//wHklKjQRXyWlpbG3LlzmTt3LsnJyceU94EDB36xfSmgDXAN0NuMZs4RARyMjGRtTAxTGzRga2wse88/n7LVqhETE0PvmBgGV6pETEwMlfI+VqhQQaPrMKNCFwmy7OxslixZwpw5c5gzZw7Lly/HOUflypVp3Ljx0bnn/OVbPSKCBt9/T61vvqHS8uVEZGTgoqKgdWvsmmugfXvKtmhBXFQUcX7voPhGhS4SBBs3bjxa4AkJCezbt4/IyEguvfRSnnzySTp16kTz5s2JjIz83xd9/z385z/w7rvesq3OeWdA3n47dOyIXXGFd4y3SB4VukghOHjwIAsXLmTOnDnMnj2b9evXA1C7dm369OlDx44dufrqq6lUqdKxX/jjjzB5slfieZc345JL4IknoHt3uOiioF3OTEKPCl0kAJxzrF69+ugofNGiRWRlZVG6dGmuuOIK7r77bjp27EjDhg1/eRjf1q0wdapX4kuXeo+1bQv/93/Qo4fvV8GR0KFCFzlNu3btYt68eUdLfOvWrQA0btyYIUOG0KlTJy677DKiT7SWdlYWzJgB48bBp59665o0awZ/+xv06QO1awd5byQcqNBFTtH+/fu54YYbSEhIwDlHTEwM7du3p1OnTnTo0IGavzWiTkqCf/4TJkyAtDSoUQOGD4f+/b2lXUXOgApd5BRNnDiRefPmMWzYMHr06EGLFi2OfTPzeFlZ8P778Oqr8NVX3pmZ3brBoEHeJc1+62tFToHlP403mOLi4lxiYqIvry1yJlq2bElGRgbffvvtb5/WvmMHjBkDb7wB27Z5p9f/4Q9w881QrVrwAktYMbPlzrkTHp2qEbrIKfjuu+9YtmwZL7300q+X+YoV8PLL3tEqWVneKHzcOOjUCXQijxQiFbrIKRg/fjwlSpTgpptuOvYJ52DBAnjuOZg7F8qVg8GDYcgQ73JnIkGgQhcpoMzMTCZMmECPHj2oUqWK92BuLkyf7hX50qXeVMpzz8Fdd0HFiv4GlmJHhS5SQNOnT2fXrl0MHDjQWzN88mR49lnvyJV69eD11+HWW0GXRROfqNBFCmj8+PHUrlmTa/bv987YTEqCJk1g4kTo3ds7ekXER3qHRuQknHNMnTKFqNmzWXzkCJG9e3tTLVOmeBeK6NdPZS5Fgn4LRX7D1q1beb1PH6774gtmADmlSsHbb3slruPHpYhRoYv8ijUffMDGPn14+sgR9lWsSM6oUUQOHAi6+rwUUZpyETleejoZgwdzfs+eXJ6TQ/oDD1Bh2zYi//AHlbkUaQUqdDPrZGbrzGyDmT18gudrm9lnZvaNma0ysy6BjypSyDIz4fnncfXrU2LcOP4ZEcG6jz+myujROnJFQsJJC93MIoHXgM5AIyDezBodt9mjwFTnXDOgL/B6oIOKFKpPP/WOXHnoITZUr04T58h+6SXiumhsIqGjICP0lsAG51yycy4LmAx0P24bB1TI+7wisDVwEUUK0ebN0KuXd3p+bi6JTz3FBevX87u+fRkyZIjf6UROSUHeFK0B/JTv/hag1XHbPAF8amZDgbJA+xN9IzMbDAwG78otIr7JzIQXXoCnn/ZO23/mGb656iratW9P06ZNGTNmzG8vvCVSBAXqTdF44C3nXE2gCzDBzH7xvZ1zY5xzcc65uKpVqwbopUVO0eefw8UXw1/+4o3M165lY79+dLn+es466yxmzpxJeV2rU0JQQQo9FaiV737NvMfyGwRMBXDOLQaigSqBCCgSMPv2eWusXH65N0KfORM++IDdFSrQpUsXDh8+zKxZszj33HP9TipyWgpS6MuABmZWz8xK4r3pOf24bTYDVwOY2YV4hZ4WyKAiZ2TGDGjUyFuf/P77YfVq6NyZjIwMevTowY8//shHH31E48aN/U4qctpOWujOuSPAEGAOsBbvaJY1ZjbSzLrlbfYn4A4z+xZ4F7jV+XXlDJH8du6Evn29KwTFxMDixd7cedmy5Obmcsstt7Bo0SLeeust2rVr53dakTNSoDNFnXMzgZnHPfZ4vs+TgDaBjSZyhqZN86ZY9u+HkSNh2DAoWfLo08OGDWPq1KmMGjWK+Ph4H4OKBIZO/Zfws3s3DB3qrYLYogW89ZY33ZLPq6++yt///nfuvvtuHnzwQX9yigSYTv2X8DJ3LjRt6q1V/uST3kWZjyvzDz/8kHvuuYfu3bvzyiuv6PBECRsqdAkPBw96l3vr0AHKl4clS+Dxx3+xrO3ixYvp168fLVu2ZNKkSURqxUQJI5pykdC3YgXEx8P69XDffd5VhE6w9soPP/xA165dqVGjBjNmzKBMmTI+hBUpPBqhS+hyDl5+GS691BuhJyTAiy+esMx37txJ586dMTNmz56NTmyTcKQRuoSm9HS47Tb4+GPvkMTx46Fy5RNueujQIbp27crWrVuZP38+9evXD3JYkeBQoUvoWbjQu2JQejq88oo3d/4rb2xmZWURHx/PsmXL+OCDD7j00kuDHFYkeFToEjpycrzFtEaOhPr1vdF5s2YAZGRksG7dOtauXUtSUtLRjz/88APZ2dm8+uqr9OjRw+cdEClcKnQJDWlp3qh83jzSr72WOV27smryZNaOGEFSUhIpKSnk5uYCEBERwXnnnUejRo3o1q0bbdq04brrrvN5B0QKn/l1hn5cXJxLTEz05bUlxCxdCr16kbtjB4/GxPDcjh0AlChRgoYNG9KoUSMuvPDCox8bNGhAdHS0z6FFCoeZLXfOxZ3oOY3QpehyDt58E3fPPeyKjqZDVha7S5dm4sSJxMXFERsbS1SUfoVFfqa/BimaDh3y1mF5+23mR0fTe/9+br7vPp5++mnKli3rdzqRIkmFLkXPjz9ypHt3Itas4Qng/bp1+Xj8eFq3bu13MpEiTScWSZHi5s0j6+KL2Z+URNeICHjsMVasXKkyFykAjdClaHCOPc8+S/nHHmOdczzapAnPvvMOF198sd/JREKGCl185zIzWdehAxcsWsSMiAh+fPJJ3h8+XG94ipwi/cWIrzYuW8a+Dh24aM8eJtSqxaVz59K1YUO/Y4mEJM2hiy9ycnJ456GHsFataLBnD/MGDuSmjRtpoDIXOW0aoUvQJSUlMe6GGxi5bh0ZpUqx7733aK8zOUXOmApdgiYrK4tRo0ax88kneSknhz1161L5iy+wGjX8jiYSFjTlIkGxbNkyWlxyCdGPP87/5eRwpH17Kn/3ncpcJIBU6FKoDh06xEMPPUS7Vq14dsMGHgT44x8pNWsWlCvndzyRsKIpFyk0Cxcu5Pbbb2fPhg2sqlaN2LQ0eOEF7zJxujCzSMBphC4Bt2/fPu666y7atWtH7YwMNlevznn79mHvvQf3368yFykkKnQJqE8++YTGjRszZswYXurXj3mHD1P6yBFYsABuuMHveCJhTYUuAZGenk7//v257rrrqFSpEmtefJF7//MfrGJF+OoraNXK74giYU9z6HJGnHNMmTKFoUOHsnfvXp544gkeiY0latAgaNQIZs+Gc87xO6ZIsaARupy21NRUunfvTnx8PLGxsaxYsYIRMTFE3Xwz/P733sWcVeYiQaNCl1PmnGPs2LE0atSIefPmMXr0aL768kuaTJ4M994LPXp4I/OKFf2OKlKsaMpFTsmPP/7IHXfcwWeffcaVV17J2LFjOa9uXe/qQmPHwu23wxtvgFZKFAm6Ao3QzayTma0zsw1m9vCvbNPbzJLMbI2ZTQpsTPFbTk4OL7zwAk2bNmX58uWMHTuWhIQEzqtVC/r08cr8kUdgzBiVuYhPTvqXZ2aRwGvANcAWYJmZTXfOJeXbpgEwHGjjnNttZtUKK7AE3+rVqxk0aBBLly6la9euvPHGG9SoUQMOH4aePWHWLHjxRe+EIRHxTUFG6C2BDc65ZOdcFjAZ6H7cNncArznndgM453YGNqb4ITc3l5EjR9K8eXNSUlKYPHky//nPf7wy378funTx5srHjlWZixQBBSn0GsBP+e5vyXssv/OB883sSzNbYmadTvSNzGywmSWaWWJaWtrpJZag+etf/8qIESPo1asXSUlJ9OnTBzODPXugQwf4/HOYMMGbNxcR3wVqsjMKaAC0A2oCi8ysqXNuT/6NnHNjgDEAcXFxLkCvLYUgISGBxx57jPj4eCZOnOgVOUBamlfma9bAtGlw/fX+BhWRowoyQk8FauW7XzPvsfy2ANOdc9nOuRRgPV7BSwhKTU0lPj6ehg0bMmbMmP+V+dat0K4dfP89TJ+uMhcpYgpS6MuABmZWz8xKAn2B6cdt8xHe6Bwzq4I3BZMcwJwSJNnZ2fTu3ZtDhw7x/vvvU+7nJW43bYLLL/c+zpoFnU44qyYiPjrplItz7oiZDQHmAJHAeOfcGjMbCSQ656bnPdfBzJKAHOBB59yuwgwuheOhhx7iq6++YvLkyVx44YXegxs3eiPzPXtg3jy49FI/I4rIrzDn/JnKjouLc4mJib68tpzYtGnT6N27N/fccw8vv/yy9+DPZb53r1fml1ziZ0SRYs/Mljvn4k70nE79FwDWrVvHwIEDad26Nc8//7z3oMpcJKSo0IWDBw/Ss2dPoqOjmTp1KiVLllSZi4QgnaNdzDnnuPPOO0lKSmLOnDnUrFlTZS4SolToxdybb77JO++8w8iRI7nmmmtU5iIhTFMuxdiyZcu499576dy5M4888ojKXCTEqdCLqV27dtGrVy/OOeccJkyYQMSWLSpzkRCnKZdiKDc3lwEDBrB9+3a++OILKmdlwdVXe8eZJySozEVClAq9GHrmmWeYNWsWr7/+Oi3q1vVG5tu2wdy5KnOREKZCL2bmzp3LiBEj6N+/P3f27euNzJOTYeZMaN3a73gicgZU6MXITz/9RL9+/WjUqBH/eP55rEsXWL3aW2jryiv9jiciZ0iFXkxkZWXRu3dvMjMz+eCddygbHw/LlnlL4GqhLZGwoEIvJh588EGWLFnC+5Mmcf7DD8PChfDOO1oCVySMqNCLgSlTpvDKK6/wwNCh3DB1KsyZA+PGQb9+fkcTkQDScehhbu3atQwaNIg2rVvzfHo6fPQRvPIKDBrkdzQRCTAVehg7cOAAPXv2pEzp0sy+4AIi3n0XnnkGhg71O5qIFAIVephyzjF48GDWrVvH4uuuo9y//gUPPADDh/sdTUQKiQo9TL3++uu8++67fHLttZz31ltw883w/PPw8/VBRSTsqNDD0Ndff83999/PqGbN6Pjxx9C1q/cmaIR+3CLhTH/hYSY9PZ0bb7yR+LPO4sHvvsPatoUpU6BECb+jiUgh02GLYSQnJ4ebbrqJOtu2MT4qCmvcGGbMgNKl/Y4mIkGgQg8jTz31FFs+/ZTlZcoQWb26d7x5xYp+xxKRIFGhh4nZs2fz7yefZEXp0pSqWNFbOfHss/2OJSJBpDn0MLB582aGxMczv1QpKpUqhc2ZA/Xq+R1LRIJMhR7iMjMz6d+zJ+/s20ed3Fxs+nRo2tTvWCLiA025hLg/338/9yUm0soMmzgRLrvM70gi4hMVegibNHEi57/xBjcAvPAC3Hij35FExEeacglRSUlJrLntNoYCOffdB/fd53ckEfGZRughaP/+/Yxv356/Z2dzuGtXSo8e7XckESkCNEIPMc45XuzalWe3bWPPRRdReupUndIvIoAKPeRMGj6cexcuZO/ZZ1NpwQKIjvY7kogUESr0ELL8o4+4YtQojkRHU/nrryEmxu9IIlKEFKjQzayTma0zsw1m9vBvbNfTzJyZxQUuogCkJydTpndvKppRYu5cIurU8TuSiBQxJy10M4sEXgM6A42AeDNrdILtygP3Al8HOmRxl5OZSXLLljTIzmbHa69RoW1bvyOJSBFUkBF6S2CDcy7ZOZcFTAa6n2C7p4BRQEYA84lzJLZtS8tdu1g8YAD177rL70QiUkQVpNBrAD/lu78l77GjzKw5UMs598lvfSMzG2xmiWaWmJaWdsphi6OkO++kVWIiM5s0oe2//+13HBEpws74TVEziwBeAP50sm2dc2Occ3HOubiqVaue6UuHvR3jxnHBmDHMq1CBdosXY7p8nIj8hoIUeipQK9/9mnmP/aw80ARYYGYbgUuB6Xpj9Mxkff015f/wB1ZGRFDviy8oU66c35FEpIgrSKEvAxqYWT0zKwn0Bab//KRzbq9zropzrq5zri6wBOjmnEsslMTFQWoqh666ivTcXHaMHct5Wj1RRArgpIXunDsCDAHmAGuBqc65NWY20sy6FXbAYufAAf7bpg0Rhw7x3q230nngQL8TiUiIKNBaLs65mcDM4x57/Fe2bXfmsYqpnBz2de1KxU2bGN6kCc+OHet3IhEJIVqcqwjJvOceKixYwMPly3P/p58SFaUfj4gUnE79LyLcP/5Bqddf52UzOs+YQfXq1f2OJCIhRkPAomDBAnL/+EdmA1nPPssVV1zhdyIRCUEqdL8lJ5Pdowc/5OYysUsXJg4b5nciEQlRKnQ/7d/PkWuv5cC+fQytXZv3J07UyUMicto0h+6X3FzcTTdh69bRLyqKF6ZPp1KlSn6nEpEQphG6Xx59FJsxg3uB3m++ycUXX+x3IhEJcRqh+2HiRHjuOd4EMgYO5LbbbvM7kYiEAY3Qg23pUnIHDmRxZCTjGjdm0auv+p1IRMKECj2YUlNx3buzDRhQpgzzPvyQ0qVL+51KRMKECj1YDh+GHj3I3LWLTtnZvDRtGrGxsX6nEpEwokIPBufg9ttxy5fT2zmuHTaMbt20rpmIBJYKPRhGj4ZJk3giKor9bdvy9NNP+51IRMKQCr2wzZ2LGzaM2eXKMbZcOb6ZPFmLbolIoVCzFKaUFFzfvmwuV46+Bw7w8cyZnH322X6nEpEwpUIvLIcOwfXXk3n4MFcfPszjf/87l112md+pRCSMqdALg3Nwxx24VavoZcbFN9zAAw884HcqEQlzKvTC8OKLMGkSfy1fnnVnn03i+PFadEtECp0KPdDmz8c9+CCfV6nCUwcOsOT996lYsaLfqUSkGFChB9KmTdC7NzvPOotr09N54623uOiii/xOJSLFhAo9UA4fhuuvJzsjg7YHD9Jv8GBuueUWv1OJSDGi1RYDwTkYPBi3ciU3OUeF5s15+eWX/U4lIsWMRuiB8Mor8M47vF69OvMyMlj+3ntER0f7nUpEihkV+pn6/HP405/4pk4dhm7axIyPP6ZevXp+pxKRYkhTLmdi+3bo04e9VarQbtMmhv/lL1x77bV+pxKRYkoj9NN15Aj07Uvu7t1cnZtL3FVXMXLkSL9TiUgxpkI/XY88AgsX8mC1amyLimLmu+8SGRnpdyoRKcZU6Kfjww/hb39jVr16vPLTTyxYsIBq1ar5nUpEijkV+qn64Qe49Va21apFj5QUnn/xRdq0aeN3KhERvSl6Sg4dgp49yTbj96mpdL/xRu69916/U4mIAAUsdDPrZGbrzGyDmT18gucfMLMkM1tlZglmVifwUX3mHNx1F271am6OiKBk/fqMGzdOi26JSJFx0kI3s0jgNaAz0AiIN7NGx232DRDnnLsIeA/4W6CD+m7sWHj7bf5VuzbTMzN5//33qVChgt+pRESOKsgIvSWwwTmX7JzLAiYD3fNv4Jz7zDl3KO/uEqBmYGP6LDERhg7l+9hYbt+0iTfffJMmTZr4nUpE5BgFKfQawE/57m/Je+zXDAJmnegJMxtsZolmlpiWllbwlH7atQt69eJQxYq0SU7mD3feSf/+/f1OJSLyCwF9U9TM+gNxwPMnet45N8Y5F+eci6tatWogX7pw5ObCgAG4rVvpcugQsXFxvPTSS36nEhE5oYIctpgK1Mp3v2beY8cws/bAI8AVzrnMwMTz2ahRMGsWT9esyXeHDrF82jRKlSrldyoRkRMqSKEvAxqYWT28Iu8L9Mu/gZk1A94EOjnndgY8pR8+/xwefZQlsbGMSE7mk5kzqVu3rt+pRER+1UmnXJxzR4AhwBxgLTDVObfGzEaaWbe8zZ4HygHTzGylmU0vtMTBkJYGffuyt2pVrklO5tHHHqNz585+pxIR+U3mnPPlhePi4lxiYqIvr/2bcnOhSxdyP/uMVs4R064ds2bN0jotIlIkmNly51zciZ7Tqf/HGzUK5szhsSpV2B4dzcyJE1XmIhISVOj5ff457tFHWVSjBn/bvp1Fn39OSByNIyKCCv1/8ubNd591FtelpjL65Zdp3bq136lERApMhQ5HjzfPTUuj/ZEjdOndm6FDh/qdSkTklKjQ4ei8+bDy5Tlco4YW3RKRkKRCz5s3T6hShdcPHmTpe+9Rvnx5v1OJiJyy4l3oefPm6RUrckN6OmMnTqRx48Z+pxIROS3Ft9Dz5s1z0tK4Jjubm//4R/r163fyrxMRKaKK7xWLRo+GOXP4c2QkpVq2ZPTo0X4nEhE5I8VzhL50Ke4vf2FuhQpMiIpihRbdEpEwUPwKfd8+XHw8u0qVou++fbw7eza1a9f2O5WIyBkrXoX+83VBU1Lo7hz3jBhBx44d/U4lIhIQxavQJ0yASZN4KjKScldfzWOPPeZ3IhGRgCk+hb5+Pe7uu/m6VCn+VbUqiVp0S0TCTB8kWtcAAAa4SURBVPEo9KwsXHw8+7OyiHeOKe+9R5UqVfxOJSISUMWj0IcPx1asYADw51dfpVWrVn4nEhEJuPA/Dn3WLHjhBV4DysbHc/fdd/udSESkUIT3CH37dnIGDGBdZCTj6tfn8zFjtOiWiISt8B2h5+aSO2AA2bt3c3PJkkz68EPKlSvndyoRkUITviP00aOJmDePocCD//oXF154od+JREQKVXgW+rJl5A4fzvtA6SFD6NOnj9+JREQKXfgV+oEDZPXuzY7cXMbGxfGxFt0SkWIi7ObQs4cOJWrjRoZUqMA/P/iAkiVL+h1JRCQowqrQ3UcfUeKtt3geGDJtGrVq1fI7kohI0ITPlMv27WQMGMBa4Mjjj3PNNdf4nUhEJKjCo9CdY+8NN1DywAHevOwy3hgxwu9EIiJBFxZTLgdHjaLi4sU8ExPDsx9+SEREWOyWiMgpCfnmy/3uO0o88ggzzeg+ezaVK1f2O5KIiC9Cu9AzM9nZoQN7cnPZ+dxztGjZ0u9EIiK+CelC3zhgAOds387bV1zBLQ895HccERFfhWyhp02ZQu1p05gSE8Ndn3yiRbdEpNgrUKGbWSczW2dmG8zs4RM8X8rMpuQ9/7WZ1Q100Pyytm/H3XwzG8xoNn8+ZcuWLcyXExEJCSctdDOLBF4DOgONgHgza3TcZoOA3c65+sCLwKhABz3KOdZcdhkxWVlseu45zv/d7wrtpUREQklBRugtgQ3OuWTnXBYwGeh+3DbdgX/nff4ecLUV0hzI10OG0GzDBua0acM1w4YVxkuIiISkghR6DeCnfPe35D12wm2cc0eAvcAvjh80s8FmlmhmiWlpaacVOKpBAxbUrEnHefNO6+tFRMJVUN8Udc6Ncc7FOefiqlatelrf45L77qPdTz9RIjo6wOlEREJbQQo9Fci/ylXNvMdOuI2ZRQEVgV2BCCgiIgVTkEJfBjQws3pmVhLoC0w/bpvpwC15n/cC5jvnXOBiiojIyZx0cS7n3BEzGwLMASKB8c65NWY2Ekh0zk0H/glMMLMNwH/xSl9ERIKoQKstOudmAjOPe+zxfJ9nADcGNpqIiJyKkD1TVEREjqVCFxEJEyp0EZEwoUIXEQkT5tfRhWaWBmw6zS+vAqQHME4o0D4XD9rn4uFM9rmOc+6EZ2b6VuhnwswSnXNxfucIJu1z8aB9Lh4Ka5815SIiEiZU6CIiYSJUC32M3wF8oH0uHrTPxUOh7HNIzqGLiMgvheoIXUREjqNCFxEJE0W60IvaxamDoQD7/ICZJZnZKjNLMLM6fuQMpJPtc77tepqZM7OQP8StIPtsZr3zftZrzGxSsDMGWgF+t2ub2Wdm9k3e73cXP3IGipmNN7OdZrb6V543M3sl77/HKjNrfsYv6pwrkje8pXp/BGKBksC3QKPjtrkb+Efe532BKX7nDsI+XwmUyfv8ruKwz3nblQcWAUuAOL9zB+Hn3AD4BojJu1/N79xB2OcxwF15nzcCNvqd+wz3+XKgObD6V57vAswCDLgU+PpMX7Moj9CL1MWpg+Sk++yc+8w5dyjv7hK8K0iFsoL8nAGeAkYBGcEMV0gKss93AK8553YDOOd2BjljoBVknx1QIe/zisDWIOYLOOfcIrzrQ/ya7sDbzrMEqGRm1c/kNYtyoQfs4tQhpCD7nN8gvP/Dh7KT7nPeP0VrOec+CWawQlSQn/P5wPlm9qWZLTGzTkFLVzgKss9PAP3NbAve9ReGBieab0717/2kCnSBCyl6zKw/EAdc4XeWwmRmEcALwK0+Rwm2KLxpl3Z4/wpbZGZNnXN7fE1VuOKBt5xzo82sNd5V0Jo453L9DhYqivIIvThenLog+4yZtQceAbo55zKDlK2wnGyfywNNgAVmthFvrnF6iL8xWpCf8xZgunMu2zmXAqzHK/hQVZB9HgRMBXDOLQai8RaxClcF+ns/FUW50IvjxalPus9m1gx4E6/MQ31eFU6yz865vc65Ks65us65unjvG3RzziX6EzcgCvK7/RHe6Bwzq4I3BZMczJABVpB93gxcDWBmF+IVelpQUwbXdODmvKNdLgX2Oue2ndF39Pud4JO8S9wFb2TyI/BI3mMj8f6gwfuBTwM2AEuBWL8zB2Gf5wE7gJV5t+l+Zy7sfT5u2wWE+FEuBfw5G95UUxLwHdDX78xB2OdGwJd4R8CsBDr4nfkM9/ddYBuQjfcvrkHAncCd+X7Gr+X99/guEL/XOvVfRCRMFOUpFxEROQUqdBGRMKFCFxEJEyp0EZEwoUIXEQkTKnQRkTChQhcRCRP/DzPzVjYkukt0AAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "tags": [],
            "needs_background": "light"
          }
        }
      ]
    }
  ]
 }
	{
	"nbformat": 4,
	"nbformat_minor": 0,
	"metadata": {
	"colab": {
	"name": "polyfit",
	"provenance": [],
	"collapsed_sections": [],
	"authorship_tag": "ABX9TyNlTrSpvUDzQrcjq22dL/Ny",
	"include_colab_link": true
	},
	"kernelspec": {
	"name": "python3",
	"display_name": "Python 3"
	}
	},
	"cells": [
	{
	"cell_type": "markdown",
	"metadata": {
	"id": "view-in-github",
	"colab_type": "text"
	},
	"source": [
	"<a href=\"https://colab.research.google.com/gist/saccadic/8a369b5dff3ea2f2de03e7401a04d9a3/polyfit.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
	]
	},
	{
	"cell_type": "code",
	"metadata": {
	"id": "4r0Bcs_mKyHe"
	},
	"source": [
	"import cv2\r\n",
	"import numpy as np\r\n",
	"import pandas as pd\r\n",
	"import matplotlib.pyplot as plt"
	],
	"execution_count": 1,
	"outputs": []
	},
	{
	"cell_type": "code",
	"metadata": {
	"id": "oDqu0XEwK3Ku"
	},
	"source": [
	"def powerXArray(a, order):\r\n",
	" data = []\r\n",
	" for i in range(len(a)):\r\n",
	" data.append(np.power(a[i], order))\r\n",
	" return np.array(data)\r\n",
	"\r\n",
	"def powerXYArray(a, b, order):\r\n",
	" data = []\r\n",
	" for i in range(len(a)):\r\n",
	" data.append(np.power(a[i],order) * b[i])\r\n",
	" return np.array(data)\r\n",
	"\r\n",
	"def MakeLeastSquaresMatrix(x, y, order):\r\n",
	" x_order = []\r\n",
	" xy_order = []\r\n",
	" \r\n",
	" for i in range(order*2+1):\r\n",
	" data_x = powerXArray(x,i)\r\n",
	" data_xy = powerXYArray(x,y,i)\r\n",
	" x_order.append(data_x.sum())\r\n",
	" xy_order.append(data_xy.sum()) \r\n",
	" \r\n",
	" # AX = B\r\n",
	" # X = A^-1 * B\r\n",
	"\r\n",
	" matA = []\r\n",
	" for i in range(order+1):\r\n",
	" mat_x =[]\r\n",
	" for j in range(order+1):\r\n",
	" if (i == 0) and (j == 0):\r\n",
	" mat_x.append(1)\r\n",
	" else:\r\n",
	" mat_x.append(x_order[j+i])\r\n",
	" matA.append(mat_x)\r\n",
	" \r\n",
	" matB = []\r\n",
	" for i in range(order+1):\r\n",
	" matB.append(xy_order[i])\r\n",
	" \r\n",
	" return np.array(matA), np.array(matB).T\r\n",
	"\r\n",
	"def model(x, param):\r\n",
	" data = []\r\n",
	" for i in range(len(x)):\r\n",
	" y = param[0]\r\n",
	" for k in range(1,len(param)):\r\n",
	" y = y + param[k] * np.power(x[i],k)\r\n",
	" data.append(y)\r\n",
	" return np.array(data)\r\n",
	"\r\n",
	"def polyfit(x, y, order):\r\n",
	" A, B = MakeLeastSquaresMatrix(x, y, order)\r\n",
	" # print(A)\r\n",
	" # print(B)\r\n",
	" \r\n",
	" # result, param = cv2.solve(A, B,cv2.DECOMP_QR)\r\n",
	" # if result == False:\r\n",
	" # return -1\r\n",
	"\r\n",
	" A_inv = np.linalg.inv(A)\r\n",
	" param = np.dot(A_inv, B)\r\n",
	"\r\n",
	" return param"
	],
	"execution_count": 2,
	"outputs": []
	},
	{
	"cell_type": "code",
	"metadata": {
	"colab": {
	"base_uri": "https://localhost:8080/",
	"height": 686
	},
	"id": "Op9ydXllKzaK",
	"outputId": "6b1a9fb9-193d-4744-9d28-f7561f215a1b"
	},
	"source": [
	"pt = np.array([\r\n",
	" [0.0,0.0],\r\n",
	" [0.1,0.2],\r\n",
	" [0.2,0.4],\r\n",
	" [0.3,0.6],\r\n",
	" [0.4,0.7],\r\n",
	" [0.41,0.71],\r\n",
	" [0.42,0.72],\r\n",
	" [0.46,0.73],\r\n",
	" [0.5,0.8],\r\n",
	" [0.51,0.82],\r\n",
	" [0.52,0.83],\r\n",
	" [0.53,0.87],\r\n",
	" [0.6,0.9],\r\n",
	" [0.7,0.91],\r\n",
	" [0.8,0.93],\r\n",
	" [0.9,0.96],\r\n",
	" [1.0,1.0]\r\n",
	"],dtype=np.double)\r\n",
	"\r\n",
	"data = pd.DataFrame.from_dict({\r\n",
	" 'x': pt.T[0],\r\n",
	" 'y': pt.T[1]\r\n",
	"})\r\n",
	"\r\n",
	"print(data)\r\n",
	"\r\n",
	"print(\"-----\")\r\n",
	"\r\n",
	"order = 5\r\n",
	"param = polyfit(data[\"x\"],data[\"y\"],order)\r\n",
	"param = param.T\r\n",
	"#param_np = np.polyfit(data[\"x\"], data[\"y\"], deg=order)\r\n",
	"\r\n",
	"print(\"param\")\r\n",
	"print(param.T)\r\n",
	"# print(\"param_np\")\r\n",
	"# print(param_np)\r\n",
	"\r\n",
	"print(\"-----\")\r\n",
	"\r\n",
	"x = np.array(range(1000)) / 1000\r\n",
	"y1 = model(x, param)\r\n",
	"# y2 = model(x, param_np)\r\n",
	"\r\n",
	"plt.clf()\r\n",
	"plt.plot(data[\"x\"],data[\"y\"], c=\"k\")\r\n",
	"plt.plot(x,y1, c=\"r\")\r\n",
	"# plt.plot(x,y2, c=\"g\")\r\n",
	"plt.show()"
	],
	"execution_count": 3,
	"outputs": [
	{
	"output_type": "stream",
	"text": [
	" x y\n",
	"0 0.00 0.00\n",
	"1 0.10 0.20\n",
	"2 0.20 0.40\n",
	"3 0.30 0.60\n",
	"4 0.40 0.70\n",
	"5 0.41 0.71\n",
	"6 0.42 0.72\n",
	"7 0.46 0.73\n",
	"8 0.50 0.80\n",
	"9 0.51 0.82\n",
	"10 0.52 0.83\n",
	"11 0.53 0.87\n",
	"12 0.60 0.90\n",
	"13 0.70 0.91\n",
	"14 0.80 0.93\n",
	"15 0.90 0.96\n",
	"16 1.00 1.00\n",
	"-----\n",
	"param\n",
	"[ 2.35978781e-04 2.21484544e+00 -1.13490566e+00 8.40614535e-01\n",
	" -2.89088057e+00 1.97100168e+00]\n",
	"-----\n"
	],
	"name": "stdout"
	},
	{
	"output_type": "display_data",
	"data": {
	"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3hUZd7G8e8vCRA6kaJIjyBK0QUDyIKKijSlKAgEwQLKqgu2XUXWgmJb1sX2WlZgWVcEKbYFpQhBwAJCQEQIgpgAElrC0iGF5Hn/OJENiBJgMiczuT/XNVcyMyeZ+5jk9uGZc55jzjlERCT0RfgdQEREAkOFLiISJlToIiJhQoUuIhImVOgiImEiyq8XrlKliqtbt65fLy8iEpKWL1+e7pyreqLnfCv0unXrkpiY6NfLi4iEJDPb9GvPacpFRCRMqNBFRMKECl1EJEyo0EVEwoQKXUQkTJy00M1svJntNLPVv/K8mdkrZrbBzFaZWfPAxxQRkZMpyAj9LaDTbzzfGWiQdxsMvHHmsURE5FSdtNCdc4uA//7GJt2Bt51nCVDJzKoHKqCISDjIyclhzcSJLL3sMr779ttCeY1AnFhUA/gp3/0teY9tO35DMxuMN4qndu3aAXhpEZGiyTnH2rVrSZg3j92TJ3P5smW0O3KEfcD0adNoevHFAX/NoJ4p6pwbA4wBiIuL05U1RCSsbNy4kYSEBObPn8/X8+Zx5c6d3AM0BXaXKcPKG2+k+ogR9G/YsFBePxCFngrUyne/Zt5jIiJhbceOHcyfP5/58+eTkJBASkoKLYB7Spfmn9nZRAOZF14Iw4cT06cPMSVLFmqeQBT6dGCImU0GWgF7nXO/mG4REQl1e/bsYeHChUcLfM2aNQA0L1eOp2rWpFP16lTetg1nht18M9xxB6VatQKzoOQ7aaGb2btAO6CKmW0BRgAlAJxz/wBmAl2ADcAh4LbCCisiEkyHDh3iyy+/PFrgy5cvJzc3l9LR0dzUvDmvdexIi5QUyqxfD99/D23bwogRWHw8VKgQ9LwnLXTnXPxJnnfAHwOWSETEJ9nZ2SxbtoyEhAQSEhJYvHgxWVlZREVFcWVcHO/06kW7w4c5Z+VK7KuvvC9q3RpefBF69YKaNX3N79vyuSIifsvNzeXbb789Og++aNEiDhw4gJlxZdOmvNmlC1dERlJ782YiExNhyRIoXx7at4dHH4XOnaFWrZO/UJCo0EWk2NizZw8rV65kxYoVLF68mM8++4xdu3ZRHuheuzbjL76YliVKUGPrVqJWrYJVq6BkSWjRAv78Z+jYEX7/e++xIkiFLiJhKS0tjRUrVhxzS09OpiHQELi8QgUeionhgshIyu/cCZs3e7ezz4ZLLoFbb/XmxFu0gOhon/emYFToIhLSnHOkpqayYsUKVi1dypbFi9mzahVl0tOpDdQB/hwdzflATP4vPHjQK++4OGjW7H+36qF7orsKXUSKNudg717YsQO3bRs7v/uOrStXsuf778nctAlLSyMmK4sWQLfjv9QMV706EfXqwfnnwwUXQMOG3i02tshOnZwuFbqIBF9ODuzaBTt3HnvbsePozW3fzpEtW4hITyfyyBEADDg773YE2B0VxeEKFaBaNSLq1yfrd7+jZP36UKcO1KmD1aiBhVlp/xYVuogEhnOwZw9s2QKpqd5tx45jy/rnz9PTITf3F98iNyKC/dHR7AA2ZmSwNTeXHcCuyEii69ThrEaNqHnJJTRo25YL2rShaunSQd/NokyFLiIFk5EBmzZBcrJ327z5f8X9c4kfPvzLrytfHqpVw519Npk1a7LnvPNIMyM1O5uUgwf5ZutWvk5JYUt2NrtzcykbEUGzZs1o3rw5zZs3p0uzZlxwwQWUKFEi+PscYlToIvI/mZnwww+QlATr1v2vvJOTvcJ2+dbUK1kSzj3XO5nmkkugWzcyKldme2Qkm3Jy+OHQIZJ27WLd5s2kpKSw8ZtvOHxc4Z9zzjk0btyYjt27Hy3w+vXrExGhi6mdDhW6SHF05Ih3qvq333rlnZQEa9fChg3e/PbPzj0XzjsPrr4aYmM5Urs220qXZkNODuv27iVl40ZSUlJI3rCBlLlz+e9/j710QoUKFahXrx4NGzakc+fO1KtX7+itbt26lClTJsg7Ht5U6CLhLisL1qyBFSu82/LlXpFnZHjPR0VB/frQuDGuVy92V6/OxjJlWJuTw4bUVK+wk5NJSUggNTUVl2+UXrJkSerUqUNsbCwtWrQ4prBjY2OJiYnBgrQwlajQRcLP7t3w1Ve4L74gZ9EiIpYvJyIzE4Ds0qVJr12bra1akVKpEuvLlmW9c2z/739JWb2aTZ98QmbetgBmxrnnnktsbCxXXXXVMWVdr149zj33XE2PFCEqdJEiLDs7m71797J79252797Nnj17jvm4e/dusrdvp9b69TT46Scu3LWL8/JG3keAb4Avga+B5cCPhw/j1q2DdeswMypWrEhMTAyVK1fmoosuonv37scUdp06dShVqpR//wHklKjQRXyWlpbG3LlzmTt3LsnJyceU94EDB36xfSmgDXAN0NuMZs4RARyMjGRtTAxTGzRga2wse88/n7LVqhETE0PvmBgGV6pETEwMlfI+VqhQQaPrMKNCFwmy7OxslixZwpw5c5gzZw7Lly/HOUflypVp3Ljx0bnn/OVbPSKCBt9/T61vvqHS8uVEZGTgoqKgdWvsmmugfXvKtmhBXFQUcX7voPhGhS4SBBs3bjxa4AkJCezbt4/IyEguvfRSnnzySTp16kTz5s2JjIz83xd9/z385z/w7rvesq3OeWdA3n47dOyIXXGFd4y3SB4VukghOHjwIAsXLmTOnDnMnj2b9evXA1C7dm369OlDx44dufrqq6lUqdKxX/jjjzB5slfieZc345JL4IknoHt3uOiioF3OTEKPCl0kAJxzrF69+ugofNGiRWRlZVG6dGmuuOIK7r77bjp27EjDhg1/eRjf1q0wdapX4kuXeo+1bQv/93/Qo4fvV8GR0KFCFzlNu3btYt68eUdLfOvWrQA0btyYIUOG0KlTJy677DKiT7SWdlYWzJgB48bBp59665o0awZ/+xv06QO1awd5byQcqNBFTtH+/fu54YYbSEhIwDlHTEwM7du3p1OnTnTo0IGavzWiTkqCf/4TJkyAtDSoUQOGD4f+/b2lXUXOgApd5BRNnDiRefPmMWzYMHr06EGLFi2OfTPzeFlZ8P778Oqr8NVX3pmZ3brBoEHeJc1+62tFToHlP403mOLi4lxiYqIvry1yJlq2bElGRgbffvvtb5/WvmMHjBkDb7wB27Z5p9f/4Q9w881QrVrwAktYMbPlzrkTHp2qEbrIKfjuu+9YtmwZL7300q+X+YoV8PLL3tEqWVneKHzcOOjUCXQijxQiFbrIKRg/fjwlSpTgpptuOvYJ52DBAnjuOZg7F8qVg8GDYcgQ73JnIkGgQhcpoMzMTCZMmECPHj2oUqWK92BuLkyf7hX50qXeVMpzz8Fdd0HFiv4GlmJHhS5SQNOnT2fXrl0MHDjQWzN88mR49lnvyJV69eD11+HWW0GXRROfqNBFCmj8+PHUrlmTa/bv987YTEqCJk1g4kTo3ds7ekXER3qHRuQknHNMnTKFqNmzWXzkCJG9e3tTLVOmeBeK6NdPZS5Fgn4LRX7D1q1beb1PH6774gtmADmlSsHbb3slruPHpYhRoYv8ijUffMDGPn14+sgR9lWsSM6oUUQOHAi6+rwUUZpyETleejoZgwdzfs+eXJ6TQ/oDD1Bh2zYi//AHlbkUaQUqdDPrZGbrzGyDmT18gudrm9lnZvaNma0ysy6BjypSyDIz4fnncfXrU2LcOP4ZEcG6jz+myujROnJFQsJJC93MIoHXgM5AIyDezBodt9mjwFTnXDOgL/B6oIOKFKpPP/WOXHnoITZUr04T58h+6SXiumhsIqGjICP0lsAG51yycy4LmAx0P24bB1TI+7wisDVwEUUK0ebN0KuXd3p+bi6JTz3FBevX87u+fRkyZIjf6UROSUHeFK0B/JTv/hag1XHbPAF8amZDgbJA+xN9IzMbDAwG78otIr7JzIQXXoCnn/ZO23/mGb656iratW9P06ZNGTNmzG8vvCVSBAXqTdF44C3nXE2gCzDBzH7xvZ1zY5xzcc65uKpVqwbopUVO0eefw8UXw1/+4o3M165lY79+dLn+es466yxmzpxJeV2rU0JQQQo9FaiV737NvMfyGwRMBXDOLQaigSqBCCgSMPv2eWusXH65N0KfORM++IDdFSrQpUsXDh8+zKxZszj33HP9TipyWgpS6MuABmZWz8xK4r3pOf24bTYDVwOY2YV4hZ4WyKAiZ2TGDGjUyFuf/P77YfVq6NyZjIwMevTowY8//shHH31E48aN/U4qctpOWujOuSPAEGAOsBbvaJY1ZjbSzLrlbfYn4A4z+xZ4F7jV+XXlDJH8du6Evn29KwTFxMDixd7cedmy5Obmcsstt7Bo0SLeeust2rVr53dakTNSoDNFnXMzgZnHPfZ4vs+TgDaBjSZyhqZN86ZY9u+HkSNh2DAoWfLo08OGDWPq1KmMGjWK+Ph4H4OKBIZO/Zfws3s3DB3qrYLYogW89ZY33ZLPq6++yt///nfuvvtuHnzwQX9yigSYTv2X8DJ3LjRt6q1V/uST3kWZjyvzDz/8kHvuuYfu3bvzyiuv6PBECRsqdAkPBw96l3vr0AHKl4clS+Dxx3+xrO3ixYvp168fLVu2ZNKkSURqxUQJI5pykdC3YgXEx8P69XDffd5VhE6w9soPP/xA165dqVGjBjNmzKBMmTI+hBUpPBqhS+hyDl5+GS691BuhJyTAiy+esMx37txJ586dMTNmz56NTmyTcKQRuoSm9HS47Tb4+GPvkMTx46Fy5RNueujQIbp27crWrVuZP38+9evXD3JYkeBQoUvoWbjQu2JQejq88oo3d/4rb2xmZWURHx/PsmXL+OCDD7j00kuDHFYkeFToEjpycrzFtEaOhPr1vdF5s2YAZGRksG7dOtauXUtSUtLRjz/88APZ2dm8+uqr9OjRw+cdEClcKnQJDWlp3qh83jzSr72WOV27smryZNaOGEFSUhIpKSnk5uYCEBERwXnnnUejRo3o1q0bbdq04brrrvN5B0QKn/l1hn5cXJxLTEz05bUlxCxdCr16kbtjB4/GxPDcjh0AlChRgoYNG9KoUSMuvPDCox8bNGhAdHS0z6FFCoeZLXfOxZ3oOY3QpehyDt58E3fPPeyKjqZDVha7S5dm4sSJxMXFERsbS1SUfoVFfqa/BimaDh3y1mF5+23mR0fTe/9+br7vPp5++mnKli3rdzqRIkmFLkXPjz9ypHt3Itas4Qng/bp1+Xj8eFq3bu13MpEiTScWSZHi5s0j6+KL2Z+URNeICHjsMVasXKkyFykAjdClaHCOPc8+S/nHHmOdczzapAnPvvMOF198sd/JREKGCl185zIzWdehAxcsWsSMiAh+fPJJ3h8+XG94ipwi/cWIrzYuW8a+Dh24aM8eJtSqxaVz59K1YUO/Y4mEJM2hiy9ycnJ456GHsFataLBnD/MGDuSmjRtpoDIXOW0aoUvQJSUlMe6GGxi5bh0ZpUqx7733aK8zOUXOmApdgiYrK4tRo0ax88kneSknhz1161L5iy+wGjX8jiYSFjTlIkGxbNkyWlxyCdGPP87/5eRwpH17Kn/3ncpcJIBU6FKoDh06xEMPPUS7Vq14dsMGHgT44x8pNWsWlCvndzyRsKIpFyk0Cxcu5Pbbb2fPhg2sqlaN2LQ0eOEF7zJxujCzSMBphC4Bt2/fPu666y7atWtH7YwMNlevznn79mHvvQf3368yFykkKnQJqE8++YTGjRszZswYXurXj3mHD1P6yBFYsABuuMHveCJhTYUuAZGenk7//v257rrrqFSpEmtefJF7//MfrGJF+OoraNXK74giYU9z6HJGnHNMmTKFoUOHsnfvXp544gkeiY0latAgaNQIZs+Gc87xO6ZIsaARupy21NRUunfvTnx8PLGxsaxYsYIRMTFE3Xwz/P733sWcVeYiQaNCl1PmnGPs2LE0atSIefPmMXr0aL768kuaTJ4M994LPXp4I/OKFf2OKlKsaMpFTsmPP/7IHXfcwWeffcaVV17J2LFjOa9uXe/qQmPHwu23wxtvgFZKFAm6Ao3QzayTma0zsw1m9vCvbNPbzJLMbI2ZTQpsTPFbTk4OL7zwAk2bNmX58uWMHTuWhIQEzqtVC/r08cr8kUdgzBiVuYhPTvqXZ2aRwGvANcAWYJmZTXfOJeXbpgEwHGjjnNttZtUKK7AE3+rVqxk0aBBLly6la9euvPHGG9SoUQMOH4aePWHWLHjxRe+EIRHxTUFG6C2BDc65ZOdcFjAZ6H7cNncArznndgM453YGNqb4ITc3l5EjR9K8eXNSUlKYPHky//nPf7wy378funTx5srHjlWZixQBBSn0GsBP+e5vyXssv/OB883sSzNbYmadTvSNzGywmSWaWWJaWtrpJZag+etf/8qIESPo1asXSUlJ9OnTBzODPXugQwf4/HOYMMGbNxcR3wVqsjMKaAC0A2oCi8ysqXNuT/6NnHNjgDEAcXFxLkCvLYUgISGBxx57jPj4eCZOnOgVOUBamlfma9bAtGlw/fX+BhWRowoyQk8FauW7XzPvsfy2ANOdc9nOuRRgPV7BSwhKTU0lPj6ehg0bMmbMmP+V+dat0K4dfP89TJ+uMhcpYgpS6MuABmZWz8xKAn2B6cdt8xHe6Bwzq4I3BZMcwJwSJNnZ2fTu3ZtDhw7x/vvvU+7nJW43bYLLL/c+zpoFnU44qyYiPjrplItz7oiZDQHmAJHAeOfcGjMbCSQ656bnPdfBzJKAHOBB59yuwgwuheOhhx7iq6++YvLkyVx44YXegxs3eiPzPXtg3jy49FI/I4rIrzDn/JnKjouLc4mJib68tpzYtGnT6N27N/fccw8vv/yy9+DPZb53r1fml1ziZ0SRYs/Mljvn4k70nE79FwDWrVvHwIEDad26Nc8//7z3oMpcJKSo0IWDBw/Ss2dPoqOjmTp1KiVLllSZi4QgnaNdzDnnuPPOO0lKSmLOnDnUrFlTZS4SolToxdybb77JO++8w8iRI7nmmmtU5iIhTFMuxdiyZcu499576dy5M4888ojKXCTEqdCLqV27dtGrVy/OOeccJkyYQMSWLSpzkRCnKZdiKDc3lwEDBrB9+3a++OILKmdlwdVXe8eZJySozEVClAq9GHrmmWeYNWsWr7/+Oi3q1vVG5tu2wdy5KnOREKZCL2bmzp3LiBEj6N+/P3f27euNzJOTYeZMaN3a73gicgZU6MXITz/9RL9+/WjUqBH/eP55rEsXWL3aW2jryiv9jiciZ0iFXkxkZWXRu3dvMjMz+eCddygbHw/LlnlL4GqhLZGwoEIvJh588EGWLFnC+5Mmcf7DD8PChfDOO1oCVySMqNCLgSlTpvDKK6/wwNCh3DB1KsyZA+PGQb9+fkcTkQDScehhbu3atQwaNIg2rVvzfHo6fPQRvPIKDBrkdzQRCTAVehg7cOAAPXv2pEzp0sy+4AIi3n0XnnkGhg71O5qIFAIVephyzjF48GDWrVvH4uuuo9y//gUPPADDh/sdTUQKiQo9TL3++uu8++67fHLttZz31ltw883w/PPw8/VBRSTsqNDD0Ndff83999/PqGbN6Pjxx9C1q/cmaIR+3CLhTH/hYSY9PZ0bb7yR+LPO4sHvvsPatoUpU6BECb+jiUgh02GLYSQnJ4ebbrqJOtu2MT4qCmvcGGbMgNKl/Y4mIkGgQg8jTz31FFs+/ZTlZcoQWb26d7x5xYp+xxKRIFGhh4nZs2fz7yefZEXp0pSqWNFbOfHss/2OJSJBpDn0MLB582aGxMczv1QpKpUqhc2ZA/Xq+R1LRIJMhR7iMjMz6d+zJ+/s20ed3Fxs+nRo2tTvWCLiA025hLg/338/9yUm0soMmzgRLrvM70gi4hMVegibNHEi57/xBjcAvPAC3Hij35FExEeacglRSUlJrLntNoYCOffdB/fd53ckEfGZRughaP/+/Yxv356/Z2dzuGtXSo8e7XckESkCNEIPMc45XuzalWe3bWPPRRdReupUndIvIoAKPeRMGj6cexcuZO/ZZ1NpwQKIjvY7kogUESr0ELL8o4+4YtQojkRHU/nrryEmxu9IIlKEFKjQzayTma0zsw1m9vBvbNfTzJyZxQUuogCkJydTpndvKppRYu5cIurU8TuSiBQxJy10M4sEXgM6A42AeDNrdILtygP3Al8HOmRxl5OZSXLLljTIzmbHa69RoW1bvyOJSBFUkBF6S2CDcy7ZOZcFTAa6n2C7p4BRQEYA84lzJLZtS8tdu1g8YAD177rL70QiUkQVpNBrAD/lu78l77GjzKw5UMs598lvfSMzG2xmiWaWmJaWdsphi6OkO++kVWIiM5s0oe2//+13HBEpws74TVEziwBeAP50sm2dc2Occ3HOubiqVaue6UuHvR3jxnHBmDHMq1CBdosXY7p8nIj8hoIUeipQK9/9mnmP/aw80ARYYGYbgUuB6Xpj9Mxkff015f/wB1ZGRFDviy8oU66c35FEpIgrSKEvAxqYWT0zKwn0Bab//KRzbq9zropzrq5zri6wBOjmnEsslMTFQWoqh666ivTcXHaMHct5Wj1RRArgpIXunDsCDAHmAGuBqc65NWY20sy6FXbAYufAAf7bpg0Rhw7x3q230nngQL8TiUiIKNBaLs65mcDM4x57/Fe2bXfmsYqpnBz2de1KxU2bGN6kCc+OHet3IhEJIVqcqwjJvOceKixYwMPly3P/p58SFaUfj4gUnE79LyLcP/5Bqddf52UzOs+YQfXq1f2OJCIhRkPAomDBAnL/+EdmA1nPPssVV1zhdyIRCUEqdL8lJ5Pdowc/5OYysUsXJg4b5nciEQlRKnQ/7d/PkWuv5cC+fQytXZv3J07UyUMicto0h+6X3FzcTTdh69bRLyqKF6ZPp1KlSn6nEpEQphG6Xx59FJsxg3uB3m++ycUXX+x3IhEJcRqh+2HiRHjuOd4EMgYO5LbbbvM7kYiEAY3Qg23pUnIHDmRxZCTjGjdm0auv+p1IRMKECj2YUlNx3buzDRhQpgzzPvyQ0qVL+51KRMKECj1YDh+GHj3I3LWLTtnZvDRtGrGxsX6nEpEwokIPBufg9ttxy5fT2zmuHTaMbt20rpmIBJYKPRhGj4ZJk3giKor9bdvy9NNP+51IRMKQCr2wzZ2LGzaM2eXKMbZcOb6ZPFmLbolIoVCzFKaUFFzfvmwuV46+Bw7w8cyZnH322X6nEpEwpUIvLIcOwfXXk3n4MFcfPszjf/87l112md+pRCSMqdALg3Nwxx24VavoZcbFN9zAAw884HcqEQlzKvTC8OKLMGkSfy1fnnVnn03i+PFadEtECp0KPdDmz8c9+CCfV6nCUwcOsOT996lYsaLfqUSkGFChB9KmTdC7NzvPOotr09N54623uOiii/xOJSLFhAo9UA4fhuuvJzsjg7YHD9Jv8GBuueUWv1OJSDGi1RYDwTkYPBi3ciU3OUeF5s15+eWX/U4lIsWMRuiB8Mor8M47vF69OvMyMlj+3ntER0f7nUpEihkV+pn6/HP405/4pk4dhm7axIyPP6ZevXp+pxKRYkhTLmdi+3bo04e9VarQbtMmhv/lL1x77bV+pxKRYkoj9NN15Aj07Uvu7t1cnZtL3FVXMXLkSL9TiUgxpkI/XY88AgsX8mC1amyLimLmu+8SGRnpdyoRKcZU6Kfjww/hb39jVr16vPLTTyxYsIBq1ar5nUpEijkV+qn64Qe49Va21apFj5QUnn/xRdq0aeN3KhERvSl6Sg4dgp49yTbj96mpdL/xRu69916/U4mIAAUsdDPrZGbrzGyDmT18gucfMLMkM1tlZglmVifwUX3mHNx1F271am6OiKBk/fqMGzdOi26JSJFx0kI3s0jgNaAz0AiIN7NGx232DRDnnLsIeA/4W6CD+m7sWHj7bf5VuzbTMzN5//33qVChgt+pRESOKsgIvSWwwTmX7JzLAiYD3fNv4Jz7zDl3KO/uEqBmYGP6LDERhg7l+9hYbt+0iTfffJMmTZr4nUpE5BgFKfQawE/57m/Je+zXDAJmnegJMxtsZolmlpiWllbwlH7atQt69eJQxYq0SU7mD3feSf/+/f1OJSLyCwF9U9TM+gNxwPMnet45N8Y5F+eci6tatWogX7pw5ObCgAG4rVvpcugQsXFxvPTSS36nEhE5oYIctpgK1Mp3v2beY8cws/bAI8AVzrnMwMTz2ahRMGsWT9esyXeHDrF82jRKlSrldyoRkRMqSKEvAxqYWT28Iu8L9Mu/gZk1A94EOjnndgY8pR8+/xwefZQlsbGMSE7mk5kzqVu3rt+pRER+1UmnXJxzR4AhwBxgLTDVObfGzEaaWbe8zZ4HygHTzGylmU0vtMTBkJYGffuyt2pVrklO5tHHHqNz585+pxIR+U3mnPPlhePi4lxiYqIvr/2bcnOhSxdyP/uMVs4R064ds2bN0jotIlIkmNly51zciZ7Tqf/HGzUK5szhsSpV2B4dzcyJE1XmIhISVOj5ff457tFHWVSjBn/bvp1Fn39OSByNIyKCCv1/8ubNd591FtelpjL65Zdp3bq136lERApMhQ5HjzfPTUuj/ZEjdOndm6FDh/qdSkTklKjQ4ei8+bDy5Tlco4YW3RKRkKRCz5s3T6hShdcPHmTpe+9Rvnx5v1OJiJyy4l3oefPm6RUrckN6OmMnTqRx48Z+pxIROS3Ft9Dz5s1z0tK4Jjubm//4R/r163fyrxMRKaKK7xWLRo+GOXP4c2QkpVq2ZPTo0X4nEhE5I8VzhL50Ke4vf2FuhQpMiIpihRbdEpEwUPwKfd8+XHw8u0qVou++fbw7eza1a9f2O5WIyBkrXoX+83VBU1Lo7hz3jBhBx44d/U4lIhIQxavQJ0yASZN4KjKScldfzWOPPeZ3IhGRgCk+hb5+Pe7uu/m6VCn+VbUqiVp0S0TCTB8kWtcAAAa4SURBVPEo9KwsXHw8+7OyiHeOKe+9R5UqVfxOJSISUMWj0IcPx1asYADw51dfpVWrVn4nEhEJuPA/Dn3WLHjhBV4DysbHc/fdd/udSESkUIT3CH37dnIGDGBdZCTj6tfn8zFjtOiWiISt8B2h5+aSO2AA2bt3c3PJkkz68EPKlSvndyoRkUITviP00aOJmDePocCD//oXF154od+JREQKVXgW+rJl5A4fzvtA6SFD6NOnj9+JREQKXfgV+oEDZPXuzY7cXMbGxfGxFt0SkWIi7ObQs4cOJWrjRoZUqMA/P/iAkiVL+h1JRCQowqrQ3UcfUeKtt3geGDJtGrVq1fI7kohI0ITPlMv27WQMGMBa4Mjjj3PNNdf4nUhEJKjCo9CdY+8NN1DywAHevOwy3hgxwu9EIiJBFxZTLgdHjaLi4sU8ExPDsx9+SEREWOyWiMgpCfnmy/3uO0o88ggzzeg+ezaVK1f2O5KIiC9Cu9AzM9nZoQN7cnPZ+dxztGjZ0u9EIiK+CelC3zhgAOds387bV1zBLQ895HccERFfhWyhp02ZQu1p05gSE8Ndn3yiRbdEpNgrUKGbWSczW2dmG8zs4RM8X8rMpuQ9/7WZ1Q100Pyytm/H3XwzG8xoNn8+ZcuWLcyXExEJCSctdDOLBF4DOgONgHgza3TcZoOA3c65+sCLwKhABz3KOdZcdhkxWVlseu45zv/d7wrtpUREQklBRugtgQ3OuWTnXBYwGeh+3DbdgX/nff4ecLUV0hzI10OG0GzDBua0acM1w4YVxkuIiISkghR6DeCnfPe35D12wm2cc0eAvcAvjh80s8FmlmhmiWlpaacVOKpBAxbUrEnHefNO6+tFRMJVUN8Udc6Ncc7FOefiqlatelrf45L77qPdTz9RIjo6wOlEREJbQQo9Fci/ylXNvMdOuI2ZRQEVgV2BCCgiIgVTkEJfBjQws3pmVhLoC0w/bpvpwC15n/cC5jvnXOBiiojIyZx0cS7n3BEzGwLMASKB8c65NWY2Ekh0zk0H/glMMLMNwH/xSl9ERIKoQKstOudmAjOPe+zxfJ9nADcGNpqIiJyKkD1TVEREjqVCFxEJEyp0EZEwoUIXEQkT5tfRhWaWBmw6zS+vAqQHME4o0D4XD9rn4uFM9rmOc+6EZ2b6VuhnwswSnXNxfucIJu1z8aB9Lh4Ka5815SIiEiZU6CIiYSJUC32M3wF8oH0uHrTPxUOh7HNIzqGLiMgvheoIXUREjqNCFxEJE0W60IvaxamDoQD7/ICZJZnZKjNLMLM6fuQMpJPtc77tepqZM7OQP8StIPtsZr3zftZrzGxSsDMGWgF+t2ub2Wdm9k3e73cXP3IGipmNN7OdZrb6V543M3sl77/HKjNrfsYv6pwrkje8pXp/BGKBksC3QKPjtrkb+Efe532BKX7nDsI+XwmUyfv8ruKwz3nblQcWAUuAOL9zB+Hn3AD4BojJu1/N79xB2OcxwF15nzcCNvqd+wz3+XKgObD6V57vAswCDLgU+PpMX7Moj9CL1MWpg+Sk++yc+8w5dyjv7hK8K0iFsoL8nAGeAkYBGcEMV0gKss93AK8553YDOOd2BjljoBVknx1QIe/zisDWIOYLOOfcIrzrQ/ya7sDbzrMEqGRm1c/kNYtyoQfs4tQhpCD7nN8gvP/Dh7KT7nPeP0VrOec+CWawQlSQn/P5wPlm9qWZLTGzTkFLVzgKss9PAP3NbAve9ReGBieab0717/2kCnSBCyl6zKw/EAdc4XeWwmRmEcALwK0+Rwm2KLxpl3Z4/wpbZGZNnXN7fE1VuOKBt5xzo82sNd5V0Jo453L9DhYqivIIvThenLog+4yZtQceAbo55zKDlK2wnGyfywNNgAVmthFvrnF6iL8xWpCf8xZgunMu2zmXAqzHK/hQVZB9HgRMBXDOLQai8RaxClcF+ns/FUW50IvjxalPus9m1gx4E6/MQ31eFU6yz865vc65Ks65us65unjvG3RzziX6EzcgCvK7/RHe6Bwzq4I3BZMczJABVpB93gxcDWBmF+IVelpQUwbXdODmvKNdLgX2Oue2ndF39Pud4JO8S9wFb2TyI/BI3mMj8f6gwfuBTwM2AEuBWL8zB2Gf5wE7gJV5t+l+Zy7sfT5u2wWE+FEuBfw5G95UUxLwHdDX78xB2OdGwJd4R8CsBDr4nfkM9/ddYBuQjfcvrkHAncCd+X7Gr+X99/guEL/XOvVfRCRMFOUpFxEROQUqdBGRMKFCFxEJEyp0EZEwoUIXEQkTKnQRkTChQhcRCRP/DzPzVjYkukt0AAAAAElFTkSuQmCC\n",
	"text/plain": [
	"<Figure size 432x288 with 1 Axes>"
	]
	},
	"metadata": {
	"tags": [],
	"needs_background": "light"
	}
	}
	]
	}
	]
	}