tspspi/leastsquares.ipynb

## leastsquares.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cdef5c8e",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "In this notebook there will be a short implementation of the least squares minimization approach using QR decomposition as has been described in the blog article at https://www.tspi.at/2021/12/18/leastsquaresqr.html. Basically the idea is to solve an overdetermined system of linear equations with the best possible solution - a solution that minimizes the residuum (or error).\n",
    "\n",
    "The set of linear equations that should be satisfied is described by $A \\vec{x} = \\vec{b}$. The matrix $A$ contains the coefficients for the variables $\\vec{x}$, the constants are all adsorbed in a constant vector $\\vec{b}$. The basic idea is to minimize the residuum $\\mid A \\vec{x} - \\vec{b} \\mid^2_2 \\to \\text{min}$. This will be done by applying QR decomposition to the coefficient matrix $A$. Details can be found in the blog post linked above.\n",
    "\n",
    "The basic idea is that we want to solve\n",
    "\n",
    "$$\n",
    "A \\vec{x} = \\vec{b}\n",
    "$$\n",
    "\n",
    "In the optimal case this would imply\n",
    "\n",
    "$$\n",
    "A \\vec{x} - \\vec{b} = 0\n",
    "$$\n",
    "\n",
    "Since we only search for the best approximation (also in case there is no exact solution) we want to minimize the residuum:\n",
    "\n",
    "$$\n",
    "\\mid A \\vec{x} - \\vec{b} \\mid^2_2 \\to \\text{min}\n",
    "$$\n",
    "\n",
    "Now one can apply QR decomposition to $A$:\n",
    "\n",
    "$$\n",
    "\\mid QR \\vec{x} - \\vec{b} \\mid^2_2 \\to \\text{min}\n",
    "$$\n",
    "\n",
    "Since Q is an orthogonal matrix - this means it doesn't change the length of the vector but only it's orientation - and we're only interested in minimizing it's length one can simply apply the orthogonal transformation $Q^T$:\n",
    "\n",
    "$$\n",
    "\\mid Q^T (QR \\vec{x} - \\vec{b}) \\mid^2_2 \\to \\text{min} \\\\ \n",
    "\\mid \\underbrace{Q^T Q}_{1} R \\vec{x} - Q^T \\vec{b} \\mid^2_2 \\to \\text{min} \\\\ \n",
    "\\mid R \\vec{x} - Q^T \\vec{b} \\mid^2_2 \\to \\text{min} \\\\ \n",
    "\\to R \\vec{x} - Q^T \\vec{b} = 0 \\\\\n",
    "\\to R \\vec{x} = Q^T \\vec{b}\n",
    "$$\n",
    "\n",
    "As one can see this reduces the problem of least squares to a simple set of linear equations with an upper triangular matrix $R$ and a simple to calculate matrix $Q^T$. The matrix vector product $Q^T \\vec{b}$ is easy to calculate.\n",
    "\n",
    "The matrix $R$ might contain rows that are fully set to zero. These do not contribute to the solution of the problem but can be used to determine the quality of the fit. The structure allows one to simply perform backsubstitution to solve the remaining equations:\n",
    "\n",
    "$$\n",
    "\\vec{c} = Q^T \\vec{b} \\\\ \n",
    "x_{\\alpha} = \\frac{c_{\\alpha} - \\sum_{i = \\alpha+1}^{n} x_i * r_{\\alpha,i}}{r_{\\alpha,\\alpha}}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "583110e0",
   "metadata": {},
   "source": [
    "### Simple matrix class\n",
    "\n",
    "Again a simple matrix class will be used. In addition the [QR decomposition function using Givens rotations](/2021/12/08/qrgivens.html) will be included to allow simple QR decomposition though in many cases an implementation based on Householder reflections might be more appropriate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "84186f6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3248bc96",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Matrix:\n",
    "    def __init__(self, rows, cols, *, rounding = 8, initList = None, name = None, description = None):\n",
    "        self.name = name\n",
    "        self.description = description\n",
    "        self.rows = rows\n",
    "        self.cols = cols\n",
    "        self.data = [ 0 ] * rows\n",
    "        self.rounding = rounding\n",
    "        for row in range(rows):\n",
    "            self.data[row] = [ 0.0 ] * cols\n",
    "        # The following would NOT work since we create the same copy of the\n",
    "        # lists and reference them\n",
    "        # self.data = [ [ 0.0 ] * cols ] * rows\n",
    "        if initList:\n",
    "            if not isinstance(initList, list):\n",
    "                raise ValueError(\"List initializer has to be a list ...\")\n",
    "            # Two ways: Either one list with rows * cols entries just ordered in row major\n",
    "            # order (first row, second row, etc.) or nested list\n",
    "            if len(initList) == (rows * cols):\n",
    "                for row in range(rows):\n",
    "                    for col in range(cols):\n",
    "                        try:\n",
    "                            k = float(initList[col + row * cols])\n",
    "                            self.data[row][col] = round(k, self.rounding)\n",
    "                        except ValueError:\n",
    "                            raise ValueError(\"Initializer list has to contain numbers\")\n",
    "            else:\n",
    "                if len(initList) != rows:\n",
    "                    raise ValueError(\"Row mismatch in initializer list\")\n",
    "                for row in range(rows):\n",
    "                    if not isinstance(initList[row], list):\n",
    "                        raise ValueError(\"Columns have to be lists in initializer\")\n",
    "                    if len(initList[row]) != cols:\n",
    "                        raise ValueError(\"Column mismatch in initializer list\")\n",
    "                    for col in range(cols):\n",
    "                        try:\n",
    "                            k = float(initList[row][col])\n",
    "                            self.data[row][col] = round(k, self.rounding)\n",
    "                        except ValueError:\n",
    "                            raise ValueError(\"Initializer list has to contain numbers\")\n",
    "    def get(self, row, col):\n",
    "        if (row < 0) or (row >= self.rows):\n",
    "            raise IndexError(\"Row {} invalid (matrix is {} x {})\".format(row, self.rows, self.cols))\n",
    "        if (col < 0) or (col >= self.cols):\n",
    "            raise IndexError(\"Column {} invalid (matrix is {} x {})\".format(col, self.rows, self.cols))\n",
    "        return self.data[row][col]\n",
    "    def set(self, row, col, value):\n",
    "        if (row < 0) or (row >= self.rows):\n",
    "            raise IndexError(\"Row {} invalid (matrix is {} x {})\".format(row, self.rows, self.cols))\n",
    "        if (col < 0) or (col >= self.cols):\n",
    "            raise IndexError(\"Column {} invalid (matrix is {} x {})\".format(col, self.rows, self.cols))\n",
    "        self.data[row][col] = round(value, self.rounding)\n",
    "    def copy(self, *, name = None, description = None):\n",
    "        result = Matrix(self.rows, self.cols, name = name, description = description)\n",
    "        for row in range(self.rows):\n",
    "            for col in range(self.cols):\n",
    "                result.set(row, col, self.get(row, col))\n",
    "        return result\n",
    "    def getName(self):\n",
    "        return self.name\n",
    "    def getDescription(self):\n",
    "        return self.description\n",
    "    def __str__(self):\n",
    "        resString = \"\"\n",
    "        if self.name:\n",
    "            resString = resString + \"{} = \".format(self.name)\n",
    "        for row in range(self.rows):\n",
    "            for col in range(self.cols):\n",
    "                resString = resString + \"{}\".format(self.data[row][col]) + \"\\t\"\n",
    "            resString = resString + \"\\n\"\n",
    "        return resString\n",
    "    def transpose(self, *, name = None, description = None):\n",
    "        result = Matrix(self.cols, self.rows, name = name, description = description)\n",
    "        for col in range(self.cols):\n",
    "            for row in range(self.rows):\n",
    "                result.set(col, row, self.get(row, col))\n",
    "        return result\n",
    "    def add(self, other, *, name = None, description = None):\n",
    "        if not isinstance(other, Matrix):\n",
    "            raise ValueError(\"Adding matrix with non matrix is not supported\")\n",
    "        if (other.cols != self.cols) or (other.rows != self.rows):\n",
    "            raise ValueError(\"Addition requires same dimensions\")\n",
    "        result = Matrix(self.rows, self.cols, name = name, description = description)\n",
    "        for col in range(self.cols):\n",
    "            for row in range(self.rows):\n",
    "                result.set(row, col, round(self.get(row, col) + other.get(row, col), self.rounding))\n",
    "        return result\n",
    "    def sub(self, other, *, name = None, description = None):\n",
    "        if not isinstance(other, Matrix):\n",
    "            raise ValueError(\"Adding matrix with non matrix is not supported\")\n",
    "        if (other.cols != self.cols) or (other.rows != self.rows):\n",
    "            raise ValueError(\"Addition requires same dimensions\")\n",
    "        result = Matrix(self.rows, self.cols, name = name, description = description)\n",
    "        for col in range(self.cols):\n",
    "            for row in range(self.rows):\n",
    "                result.set(row, col, round(self.get(row, col) - other.get(row, col), self.rounding))\n",
    "        return result\n",
    "    def multiply(self, other, *, name = None, description = None):\n",
    "        # Calculate self * other\n",
    "        if isinstance(other, Matrix):\n",
    "            if not self.cols == other.rows:\n",
    "                raise ValueError(\"Colums of left side {} have to match rows of right side {}\".format(self.cols, other.rows))\n",
    "            result = Matrix(self.rows, other.cols, name = name, description = description)\n",
    "            for row in range(self.rows):\n",
    "                for col in range(other.cols):\n",
    "                    t = 0\n",
    "                    for k in range(self.cols):\n",
    "                        t = t + self.get(row, k) * other.get(k, col)\n",
    "                    t = round(t, 4)\n",
    "                    result.set(row, col, round(t, self.rounding))\n",
    "            return result\n",
    "        else:\n",
    "            try:\n",
    "                otherFloat = float(other)\n",
    "                result = Matrix(self.rows, self.cols, name = name, description = description)\n",
    "                for r in range(self.rows):\n",
    "                    for c in range(self.cols):\n",
    "                        result.set(r,c, round(self.get(r,c,)*otherFloat, self.rounding))\n",
    "                return result\n",
    "            except ValueError:\n",
    "                raise ValueError(\"Matrix has to be multiplied by matrix or scalar\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "76ffd930",
   "metadata": {},
   "outputs": [],
   "source": [
    "class IdentityMatrix(Matrix):\n",
    "    def __init__(self, size, *, name = None, description = None):\n",
    "        super().__init__(size, size, name = name, description = description)\n",
    "        for k in range(size):\n",
    "            self.set(k,k,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "571a598d",
   "metadata": {},
   "outputs": [],
   "source": [
    "class GivensRotationMatrix(IdentityMatrix):\n",
    "    def __init__(self, size, i, j, x1, x2, *, name = None, description = None):\n",
    "        super().__init__(size, name = name, description = description)\n",
    "        if (i >= size) or (i < 0) or (j >= size) or (j < 0):\n",
    "            raise ValueError(\"Invalid first index i or j\")\n",
    "        # Calculate s and c\n",
    "        s = 1.0 / math.sqrt(x1 * x1 / (x2 * x2) + 1)\n",
    "        c = s * x1 / x2\n",
    "        self.set(i,i,round(c, self.rounding))\n",
    "        self.set(j,i,round(-1 * s, self.rounding))\n",
    "        self.set(i,j,round(s, self.rounding))\n",
    "        self.set(j,j,round(c, self.rounding))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "bfeb5286",
   "metadata": {},
   "outputs": [],
   "source": [
    "def QRDecomposeWithGivens(matIn):\n",
    "    matR = matIn.copy()\n",
    "    matQ = IdentityMatrix(matIn.rows)\n",
    "    for col in range(matIn.cols):\n",
    "        for row in range(matIn.rows-1, col, -1):\n",
    "            if matR.get(row,col) == 0:\n",
    "                continue\n",
    "            rotationMatrix = GivensRotationMatrix(matIn.rows, row-1, row, matR.get(row-1, col), matR.get(row,col))\n",
    "            matR = rotationMatrix.multiply(matR)\n",
    "            matQ = matQ.multiply(rotationMatrix.transpose())\n",
    "\n",
    "    return (matQ, matR)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f634341",
   "metadata": {},
   "source": [
    "## Implementing least squares"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81798434",
   "metadata": {},
   "source": [
    "The whole procedure will be wrapped into a single ```leastSquares``` function that just accepts the coefficient matrix and the constant vector. In the end it will return a variable vector $\\vec{x}$. Note that this naive implementation assumes that the system is not underdetermined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3500b858",
   "metadata": {},
   "outputs": [],
   "source": [
    "def leastSquares_CheckRowEmpty(matA, row, *, epsilon = 1e-3):\n",
    "    for k in range(matA.cols):\n",
    "        v = matA.get(row, k)\n",
    "        if v < (-1.0 * epsilon) or v > epsilon:\n",
    "            return False\n",
    "    return True\n",
    "\n",
    "\n",
    "def leastSquares(matA, vecB):\n",
    "    if not isinstance(matA, Matrix):\n",
    "        raise ValueError(\"Coefficient matrix has to be a matrix type\")\n",
    "    if not isinstance(vecB, Matrix):\n",
    "        raise ValueError(\"Constant matrix has to be a vector (matrix) type\")\n",
    "    if vecB.cols != 1:\n",
    "        raise ValueError(\"Constant vector has to have only one row\")\n",
    "    if vecB.rows != matA.rows:\n",
    "        raise ValueError(\"Columns in coefficient matrix have to match the number of constants in constant vector\")\n",
    "\n",
    "    q,r = QRDecomposeWithGivens(matA)\n",
    "    bTransformed = q.transpose().multiply(vecB)\n",
    "\n",
    "    # Check how many rows we have to drop (i.e. how many variables we're going to ignore)\n",
    "    populatedRows = matA.rows\n",
    "    while leastSquares_CheckRowEmpty(r, populatedRows - 1) and (populatedRows > 0):\n",
    "        populatedRows = populatedRows - 1\n",
    "    if populatedRows == 0:\n",
    "        return None\n",
    "\n",
    "    # Initialize our result vector that will contain our best fit for our unknowns\n",
    "    xResult = [ 0.0 ] * matA.cols\n",
    "\n",
    "    # Now perform backsubstitution ...\n",
    "    for row in range(populatedRows-1, -1, -1):\n",
    "        t = bTransformed.get(row, 0)\n",
    "        for i in range(row+1, populatedRows):\n",
    "            t = t - xResult[i] * r.get(row, i)\n",
    "        t = t / r.get(row, row)\n",
    "        xResult[row] = t\n",
    "\n",
    "    return xResult"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e60c59eb",
   "metadata": {},
   "source": [
    "### Two mini examples\n",
    "\n",
    "In the following section I've ran my implementation against two sets of equations that I knew the solutions for by having them calculated by hand and having them compared to solutions someone else had calculated. This is of course no formal proof that the implementation works but it's more likely it's correct in this case."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "0955a832",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A:\n",
      "2.0\t0.0\t\n",
      "0.0\t1.0\t\n",
      "1.0\t2.0\t\n",
      "\n",
      "b:\n",
      "1.0\t\n",
      "0.0\t\n",
      "3.0\t\n",
      "\n",
      "approxX:\n",
      "[0.6190413105929794, 0.9523275104908754]\n",
      "\n",
      "\n",
      "A:\n",
      "2.0\t3.0\t\n",
      "-3.0\t-4.0\t\n",
      "\n",
      "b:\n",
      "-6.0\t\n",
      "7.0\t\n",
      "\n",
      "approxX:\n",
      "[3.000360460272565, -4.000360490266763]\n"
     ]
    }
   ],
   "source": [
    "# Example 1\n",
    "matA = Matrix(3, 2, initList = [ [ 2, 0 ], [ 0, 1 ], [ 1, 2 ] ] )\n",
    "matB = Matrix(3, 1, initList = [ 1, 0, 3 ] )\n",
    "print(\"A:\")\n",
    "print(matA)\n",
    "print(\"b:\")\n",
    "print(matB)\n",
    "print(\"approxX:\")\n",
    "print(leastSquares(matA, matB))\n",
    "\n",
    "# Example 2\n",
    "matA = Matrix(2, 2, initList = [ [ 2, 3 ], [ -3, -4 ] ] )\n",
    "matB = Matrix(2, 1, initList = [ -6, 7 ] )\n",
    "print(\"\\n\\nA:\")\n",
    "print(matA)\n",
    "print(\"b:\")\n",
    "print(matB)\n",
    "print(\"approxX:\")\n",
    "print(leastSquares(matA, matB))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46bfa3c7",
   "metadata": {},
   "source": [
    "## Using least squares to fit linear functions\n",
    "\n",
    "This is a pretty simple application of least squares. The idea is to have a set of data points $(x,y)$ that should in optimal case fulfill a linear equation $f(x) = k * x + d$. Usually measured data of course doesn't exactly fit one's model due to noise or modelling errors. So one can use least squares to perform linear regression - one can try to find the best approximation for a given dataset with ones model function.\n",
    "\n",
    "In this case the coefficient matrix $A$ is formed by the prefactors of $k$ and the constant $1$ for the unknowns $k$ and $d$ in the equation $k*x + d*1$. The constants are the $y$ values.\n",
    "\n",
    "In the following example a measured dataset was given by\n",
    "\n",
    "| y   | x   |\n",
    "| --- | --- |\n",
    "|  4  | 2   |\n",
    "|  5  | 3   |\n",
    "|  7  | 5   |\n",
    "| 10  | 7   |\n",
    "| 15  | 9   |\n",
    "\n",
    "Least squares will be used to find the best solution for $k$ and $d$ for the following equations:\n",
    "\n",
    "$$\n",
    " 4 = 2 * k + d \\\\\n",
    " 5 = 3 * k + d \\\\\n",
    " 7 = 5 * k + d \\\\\n",
    "10 = 7 * k + d \\\\\n",
    "15 = 9 * k + d\n",
    "$$\n",
    "\n",
    "The coefficient matrix is given by\n",
    "\n",
    "$$\n",
    "A = \\begin{pmatrix}\n",
    "2 & 1 \\\\ \n",
    "3 & 1 \\\\ \n",
    "5 & 1 \\\\ \n",
    "7 & 1 \\\\ \n",
    "9 & 1\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "The constant vector is formed by\n",
    "\n",
    "$$\n",
    "b = \\begin{pmatrix} 4 \\\\ 5 \\\\ 7 \\\\ 10 \\\\ 15 \\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "088289d2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "A:\n",
      "2.0\t1.0\t\n",
      "3.0\t1.0\t\n",
      "5.0\t1.0\t\n",
      "7.0\t1.0\t\n",
      "9.0\t1.0\t\n",
      "\n",
      "b:\n",
      "4.0\t\n",
      "5.0\t\n",
      "7.0\t\n",
      "10.0\t\n",
      "15.0\t\n",
      "\n",
      "approxX:\n",
      "[1.517810419839562, 0.3084008097165992]\n"
     ]
    }
   ],
   "source": [
    "matA = Matrix(5, 2, initList = [ [ 2, 1 ], [ 3, 1 ], [ 5, 1 ], [ 7, 1 ], [ 9, 1 ] ] )\n",
    "matB = Matrix(5, 1, initList = [ 4, 5, 7, 10, 15 ] )\n",
    "print(\"\\n\\nA:\")\n",
    "print(matA)\n",
    "print(\"b:\")\n",
    "print(matB)\n",
    "print(\"approxX:\")\n",
    "fitResult = leastSquares(matA, matB)\n",
    "print(fitResult)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73356a70",
   "metadata": {},
   "source": [
    "Now lets visualize the result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "947fba63",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "badb68d4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEICAYAAABRSj9aAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAArrElEQVR4nO3dd3iUdbr/8fdNL4EAAqGEKk0IsRDBbrAC6iLqupa1u6hnPe6e465gZ7Gsuuuqu1bs/nRlWZqKKCoSsCtYEnoNEFqogZBA2v37Y4azERMJkwlT8nldV66Zeer9zYTPPDzzzD3m7oiISPyqE+kCRESkZinoRUTinIJeRCTOKehFROKcgl5EJM4p6EVE4pyCXiLGzE42syWRriNWmdkCM0uPdB0S/UzX0UtNM7Ns4Hp3/yjStYjURjqil1rHzOqFY5mD2J+Zmf6tScToj08ixszSzSyn3ONsM/uDmWWaWZ6Z/cvMGpWbf66ZfW9mO8zsczNLLTdvtJmtMLNdZrbQzEaUm3e1mX1mZo+Z2TZgTAW1jDGziWb2upntBK42s0Qze9HMNpjZOjO738zqBpeva2aPmtkWM1tlZjebme97gTCzDDN7wMw+AwqA7mbWx8w+NLNtZrbEzC4ut/9hwbp3Bff1h+D01mY2LTjmbWb2yb4XjeDv64zg/YZm9riZrQ/+PG5mDcv/ns3sVjPLDY7nmrA8iRITFPQSbS4GhgDdgFTgagAzOwZ4CbgBOAx4Dnh7X5gBK4CTgUTgT8DrZta+3HYHASuBtsADlex7ODARaAG8AbwKlAA9gKOBs4Drg8v+BhgKHAUcA5xfwfauAEYCzYDNwIfAP4M1XAo8bWb9gsu+CNzg7s2AFODj4PRbgRygDZAE3AFUdL71TuC4YD1HAgOBu8rNb0fgd9MRuA54ysxaVvJ7kDijoJdo83d3X+/u24B3CAQXBIL1OXf/yt1L3f1VYC+BcMPd/x1cr8zd/wUsIxB2+6x393+4e4m7F1ay7y/cfaq7lwHNCQT57919t7vnAo8BlwSXvRh4wt1z3H078FAF23vF3Re4ewmBF69sd385WMO3wCTgouCyxUBfM2vu7tuD8/dNbw90cfdid//EK35j7XJgrLvnuvtmAi92V5SbXxycX+zu04F8oHclvweJMwp6iTYby90vABKC97sAtwZPYewwsx1AJ6ADgJldWe60zg4CR8Wty21rbRX2XX6ZLkB9YEO5bT5H4Gic4H7XVrJuZdsbtF/9lxM40ga4EBgGrDaz2WZ2fHD6X4DlwAdmttLMRldSewdgdbnHq4PT9tkafMHZp/zvVuJc2N5wEqlha4EH3P0np13MrAvwPHA6gaPyUjP7HrByi1Xl8rLyy6wl8D+G1vsF5D4bgORyjztVYXuz3f3MCnfs/g0w3MzqAzcDE4BO7r6LwOmbW4OneWaZ2TfuPnO/Tawn8GKyIPi4c3CaiI7o5ZCpb2aNyv0c7EHG88CNZjYoeBVLUzM7x8yaAU0JhOpmgOAbjSnVKdbdNwAfAI+aWXMzq2Nmh5vZqcFFJgC/M7OOZtYCGHWATU4DepnZFWZWP/hzrJkdYWYNzOxyM0t092JgJ1AaHMu5ZtbDzKzc9NIKtv8mcJeZtTGz1sA9wOvV+R1I/FDQy6EyHSgs9zPmYFZ297kEztM/CWwncDrj6uC8hcCjwBfAJqA/8FkYar4SaAAsDO5zIoHz5RB44fkAyAS+IzC+EioOYYJH5mcROMe/nsApqoeBfW8mXwFkB6/4uRH4dXB6T+AjAufUvwCedveMCnZxPzA3WE8W8G1wmog+MCUSDmY2FHjW3btEuhaR/emIXiQEZtY4eO17PTPrCNwLTIl0XSIV0RG9SAjMrAkwG+hD4FTUu8Dv3H1nRAsTqYCCXkQkzunUjYhInIvK6+hbt27tXbt2DWnd3bt307Rp0/AWFCU0ttgVz+PT2KLDvHnztrh7m4rmVaWL30vAuUCuu6cEp/2L/3x8ugWww92PqmDdbGAXgUvOStw9rSoFd+3alblz51Zl0Z/IyMggPT09pHWjncYWu+J5fBpbdDCz1ZXNq8oR/SsErl1+bd8Ed/9VuY0/CuT9zPqD3X1LFfYjIiI14IBB7+5zzKxrRfOCn9a7GDgtzHWJiEiYVOmqm2DQT9t36qbc9FOAv1V2SsbMVhH4RKET6Dw47mf2MZJAS1eSkpIGjB8/vqpj+JH8/HwSEuKzV5PGFrvieXwaW3QYPHjwvEpPj7v7AX+ArsD8CqY/A9z6M+t1CN62BX4ATqnK/gYMGOChmjVrVsjrRjuNLXbF8/g0tugAzPVKMjXkyyuDTakuAP5V2TLuvj54m0vgU4MDK1tWRERqRnWuoz8DWOzuORXNDHYXbLbvPoGGTvOrsT8REQnBAYPezN4k0DWvd/B7J68LzrqEQGvU8st2MLPpwYdJwKdm9gPwNfCuu78fvtJFRKQqqnLVzaWVTL+6gmnrCXxLDu6+ksB3V4qIyAF8k72Neau3c+Oph4d921H5yVgRkdoif28Jj7y/mNe+WE3nVk248vguNGkQ3mhW0IuIRMjspZu5Y3IW6/MKuebErvzhrN5hD3lQ0IuIHHLbdxdx37sLmfztOnq0TWDijScwoEvLGtufgl5E5BBxd96bv5F73prPjoJi/vu0Htx8Wg8a1qtbo/tV0IuIHAK5O/dw91vzmbFgE/07JvLatYPo26H5Idm3gl5EpAa5O/+el8P90xayt6SM0UP7cP1J3ahX99B9HYiCXkSkhqzdVsDtk7P4dPkWBnZrxUMX9Kd7m0PfO0dBLyISZqVlzqufZ/OXGUuoW8e4//wULhvYmTp1LCL1KOhFRMJo2aZdjJqUybdrdpDeuw0PjuhPhxaNI1qTgl5EJAyKS8t4NmMF//h4OU0b1uXxXx3F8KM6EPjajshS0IuIVFNWTh5/nPgDizfu4tzU9oz5RT9aJzSMdFn/R0EvIhKiPcWlPPbRUp6fs5I2zRoy7ooBnNWvXaTL+gkFvYhICL5cuZXRkzLJ3lrApQM7MXroESQ2rh/psiqkoBcROQi79hTz0HuLeeOrNXRu1YR/Xj+IE3q0jnRZP0tBLyJSRbMW53LHlCw27dzD9Sd143/P6lUjTcjCLforFBGJsG27ixj7zgKmfr+enm0TePqmEzi6c801IQs3Bb2ISCXcnXd+WM+Ytxewc08xvzu9J/81+PAab0IWbgp6EZEKbMzbw9+/28t3ud9xZHIiD180iD7tDk0TsnBT0IuIlOPujP9mLQ++u4i9xaXcOewIrj2pG3Uj1L4gHBT0IiJBq7fuZvSkLL5YuZXjurdiRMdCfnVK90iXVW0KehGp9UrLnJc/W8VfP1hC/Tp1eHBEfy45thNz5syOdGlhccCGyGb2kpnlmtn8ctPGmNk6M/s++DOsknWHmNkSM1tuZqPDWbiI1DKZE+CxFBjTInCbOSEsm12ycRcXPPM597+7iBMPb80H/3sKlw2KXKfJmlCVI/pXgCeB1/ab/pi7/7WylcysLvAUcCaQA3xjZm+7+8IQaxWR2ipzArxzCxQXBh7nrQ08Bki9OKRNFpWU8XTGcp6atZxmjerz90uP5rzU9lHRhCzcDhj07j7HzLqGsO2BwHJ3XwlgZuOB4YCCXkQOzsyx/wn5fYoLA9NDCPrv1+5g1MRMlmzaxfCjOnDvef1o1bRBmIqNPubuB14oEPTT3D0l+HgMcDWwE5gL3Oru2/db5yJgiLtfH3x8BTDI3W+uZB8jgZEASUlJA8aPHx/SgPLz80lIOPTf4HIoaGyxK57HdyjGdmrG+Rg/zSrHmJ0+tcrb2VvqTFlWxIzsElo0NK7q14Cj2lZ+vBtLz9vgwYPnuXtaRfNCfTP2GeA+wIO3jwLX7rdMRf//qfRVxd3HAeMA0tLSPD09PaTCMjIyCHXdaKexxa54Ht8hGdt3yYHTNfuxxOQq7/vzFVu4Z1IWa7aVcNmgzowe2ofmjX6+CVm8PG8hfTutu29y91J3LwOeJ3CaZn85QKdyj5OB9aHsT0RqudPvgfr7fUtT/caB6Qewc08xt0/O4rLnv6KOwZu/OY4HR/Q/YMjHk5CO6M2svbtvCD4cAcyvYLFvgJ5m1g1YB1wCXBZSlSJSu+07Dz9zLOTlQGJyIOQPcH7+o4WbuHNqFpt37eWGU7rz+zN60bhBbLUvCIcDBr2ZvQmkA63NLAe4F0g3s6MInIrJBm4ILtsBeMHdh7l7iZndDMwA6gIvufuCmhiEiNQCqRdX+Y3Xrfl7GfPOQt75YT192jXj+SvTSE1uUbP1RbGqXHVzaQWTX6xk2fXAsHKPpwPTQ65OROQguDtvB5uQ5e8t4X/P7MWNpx5Og3ohnaWOG/pkrIjEhfU7Crlr6nw+XpzLUZ1a8MhFqfRKahbpsqKCgl5EYlpZmfPmN2v48/TFlJY5d5/bl6tP6BrTTcjCTUEvIjFr1ZbdjJ6UyVertnFij8P484hUOh/WJNJlRR0FvYjEnJLSMl78dBV/+3ApDerV4ZELU/llWnJcti8IBwW9iMSURRt2MmpSJpk5eZzZN4n7z08hqXmjSJcV1RT0IhIT9paU8tTHy3k6YwUtmtTnqcuOYVj/djqKrwIFvYhEvW/XbGfUxEyW5eYz4uiO3HNuX1rGcROycFPQi0jUKigq4a8zlvLy56to37wRL19zLIN7t410WTFHQS8iUemz5VsYPTmTtdsKueK4Ltw2pDfNalF/mnBS0ItIVMkrLObBdxfxr7lr6da6KRNuOJ6B3VpFuqyYpqAXkagxY8FG7p46n627i7gp/XB+d3pPGtWvfU3Iwk1BLyIRt3nXXsa8vYB3szZwRPvmvHjVsfRPTox0WXFDQS8iEePuTPluHWOnLaRgbyl/PLs3I0/pTv26tbsJWbgp6EUkItbtKOTOKVlkLNnMMZ0DTch6tFUTspqgoBeRQ6qszHnjq9U89N5iHBhzXl+uOF5NyGqSgl5EDpkVm/O5fVIWX2dv4+SerXlwRH86tVITspqmoBeRGldSWsa4T1by+EfLaFSvDn+5KJWLBqgJ2aGioBeRGrVgfR6jJmUyf91Ozu6XxH3DU2irJmSHlIJeRGrEnuJS/vHxMp6dvZKWTRrwzOXHMLR/+0iXVSsp6EUk7Oat3sZtEzNZsXk3Fx6TzN3nHkGLJmpCFikHDHozewk4F8h195TgtL8A5wFFwArgGnffUcG62cAuoBQocfe0sFUuIlFn994S/jJjCa9+kU2HxMa8eu1ATu3VJtJl1XpV+VTCK8CQ/aZ9CKS4eyqwFLj9Z9Yf7O5HKeRF4tucpZs567E5vPpFNlce14UZ/3OKQj5KHPCI3t3nmFnX/aZ9UO7hl8BFYa5LRGJEfpHzh3//wMR5OXRvE2hCdmxXNSGLJuE4R38t8K9K5jnwgZk58Jy7jwvD/kQkSrw/fwN3fFrI7pJ1/Ff64dyiJmRRydz9wAsFjuin7TtHX276nUAacIFXsCEz6+Du682sLYHTPf/t7nMq2cdIYCRAUlLSgPHjxx/sWADIz88nISEhpHWjncYWu+JtfDv2lvH6wiLmbioluanzmyMb06V5/AV8LD1vgwcPnlfZKfKQj+jN7CoCb9KeXlHIA7j7+uBtrplNAQYCFQZ98Gh/HEBaWpqnp6eHVFdGRgahrhvtNLbYFS/jc3cmzsvh/tmLKCx2/nh2b3r7Ws44bXCkS6sR8fK8hRT0ZjYEGAWc6u4FlSzTFKjj7ruC988CxoZcqYhE1NptBdwxJYtPlm0hrUtLHrowlR5tE8jIyIl0aXIAVbm88k0gHWhtZjnAvQSusmkIfBj8CPOX7n6jmXUAXnD3YUASMCU4vx7wT3d/v0ZGISI1pqzMee2LbB6ZsQQDxg7vx68HdaGOmpDFjKpcdXNpBZNfrGTZ9cCw4P2VwJHVqk5EImp57i5GTcpi3urtnNKrDQ+OSCG5pZqQxRp9MlZEfqK4tIxxc1byxEfLaNygLo/+8kguOKajmpDFKAW9iPzI/HV53DYxk4UbdnJO//aM+UU/2jRrGOmypBoU9CICBJqQPTFzGePmrKRV0wY8++sBDElpF+myJAwU9CLCN9nbGDUxk5VbdnNxWjJ3DutLYpP6kS5LwkRBL1KL5e8t4ZH3F/PaF6tJbtmY168bxEk9W0e6LAkzBb1ILTVrSS53Ts5iw849XHtiN249qxdNGyoS4pGeVZFaZvvuIu6btpDJ362jR9sEJt54AgO6tIx0WVKDFPQitYS7Mz1rI/e+PZ8dBcXccloPfntaDxrWi78eNfJjCnqRWiB35x7umjqfDxZuon/HRF67dhB9OzSPdFlyiCjoReKYu/PvuTnc9+5CikrKuH1oH647qRv16lblO4ckXijoReLUmq0F3D4lk8+Wb2Vgt1Y8fGEq3Vo3jXRZEgEKepE4U1rmvPJ5Nn+dsYS6dYz7z0/hsoGd1YSsFlPQi8SRZZt2cdukTL5bs4PBvdvwwIj+dGjRONJlSYQp6EXiQFFJGc/OXsGTHy+nacO6PP6roxh+VAc1IRNAQS8S8zJzdnDbxEwWb9zFeUd2YMx5fTksQU3I5D8U9CIxqrColMc/Wsrzn6ykTbOGPH9lGmf2TYp0WRKFFPQiMejLlVsZPSmT7K0FXDqwE7cPO4LmjdSETCqmoBeJIbv2FPPQe4t546s1dG7VhH9eP4gTeqgJmfw8Bb1IjPh48SbunDKfTTv3cP1J3bj1rN40bqD2BXJgCnqRKLdtdxFj31nA1O/X0yspgacvP4GjO6sJmVSdgl4kSrk772RuYMzbC9i1p5jfn9GT/0rvQYN6al8gB+eAfzFm9pKZ5ZrZ/HLTWpnZh2a2LHhb4eGFmQ0xsyVmttzMRoezcJF4tjFvD795bS63vPkdnVo1Ydp/n8zvz+j18yGfOQEeS4ExLQK3mRMOWb0S3apyaPAKMGS/aaOBme7eE5gZfPwjZlYXeAoYCvQFLjWzvtWqViTOuTtvfr2GM/82m0+Xb+Guc45g8k0n0Ltds59fMXMCvHML5K0FPHD7zi0KewGqEPTuPgfYtt/k4cCrwfuvAudXsOpAYLm7r3T3ImB8cD0RqcDqrbu57PmvuH1yFikdE5nx+1O4/uTu1K1Kj5qZY6G48MfTigsD06XWM3c/8EJmXYFp7p4SfLzD3VuUm7/d3Vvut85FwBB3vz74+ApgkLvfXMk+RgIjAZKSkgaMHz8+pAHl5+eTkJAQ0rrRTmOLXT83vjJ3PsguYfKyIurWgV/1bsCpyfUOqn3BqRnnY/z037JjzE6fGmrZVRLPz10sjW3w4MHz3D2tonk1+WZsRX+llb6quPs4YBxAWlqap6enh7TTjIwMQl032mlssauy8S3ZGGhC9sPaAs44oi33n9+fdomNDn4H3yUHT9v8mCUm1/jvNZ6fu3gZW6hBv8nM2rv7BjNrD+RWsEwO0Knc42RgfYj7E4krRSVlPDVrOU9nLKd5o/r849KjOTe1fehNyE6/J3BOvvzpm/qNA9Ol1gs16N8GrgIeCt6+VcEy3wA9zawbsA64BLgsxP2JxI3v1+7gtok/sHRTPucf1YF7zutHq6YNqrfR1IsDtzPHQl4OJCYHQn7fdKnVDhj0ZvYmkA60NrMc4F4CAT/BzK4D1gC/DC7bAXjB3Ye5e4mZ3QzMAOoCL7n7gpoZhkj0Kywq5dEPlvDSZ6tIat6Il65O47Q+YWxClnqxgl0qdMCgd/dLK5l1egXLrgeGlXs8HZgecnUicWLR1lLueXwOa7YVcPmgzowe2odmakImh4g+GStSg3buKebP0xfx5jd76HpYE8aPPI7juh8W6bKkllHQi9SQDxdu4q6pWWzetZeh3erz2LWn0Ki+mpDJoaegFwmzLfl7GfP2AqZlbqBPu2Y8f2Ua25Z/r5CXiFHQi4SJu/PW9+v50zsL2L23lFvP7MUNpx5Og3p1yFge6eqkNlPQi4TB+h2F3DV1Ph8vzuXozi145MJUeiYdoD+NyCGioBephrIy559fr+Gh9xZTWubcc25frjqha9X604gcIgp6kRCt2rKb0ZMy+WrVNk7q0Zo/X9CfTq2aRLoskZ9Q0IscpJLSMl78dBV/+3ApDerV4ZELU/llWnLo7QtEapiCXuQgLFy/k1GTMslal8dZfZO47/wUkpqH0IRM5BBS0ItUwd6SUp78eDnPZKygRZP6PHXZMQzr305H8RITFPQiBzBv9XZGTcpkeW4+FxzTkbvP6UvL6jYhEzmEFPQilSgoKuEvM5bwyufZtG/eiJevOZbBvdtGuiyRg6agF6nAp8u2MHpyJjnbC7ny+C7cNqQPCQ31z0Vik/5yRcrJKyzmgXcXMmFuDt1bN2XCDcczsFurSJclUi0KepGgGQs2cvfU+WzdXcRN6Yfzu9N7qj+NxAUFvdR6m3cFmpC9m7WBvu2b89LVx5LSMTHSZYmEjYJeai13Z/K36xg7bSGFRaX88ezejDylO/Xr1ol0aSJhpaCXWmndjkLumJzF7KWbGdClJQ9fmEqPtgmRLkukRijopVYpK3Pe+Go1D723GAfGnNeXK4/vSh01IZM4pqCXWmPF5nxGT8rkm+ztnNyzNQ+OUBMyqR0U9BL3SkrLGPfJSh7/aBmN69flr788kguP6aj2BVJrhBz0ZtYb+Fe5Sd2Be9z98XLLpANvAauCkya7+9hQ9ylysBasz2PUpEzmr9vJ0JR2/Gl4P9o2UxMyqV1CDnp3XwIcBWBmdYF1wJQKFv3E3c8NdT8iodhTXMo/Pl7Gs7NX0rJJA565/BiG9m8f6bJEIiJcp25OB1a4++owbU8kZHOzt3HbpExWbt7NRQOSueucI2jRRE3IpPYyd6/+RsxeAr519yf3m54OTAJygPXAH9x9QSXbGAmMBEhKShowfvz4kGrJz88nISE+L5PT2H7enhJn4tIiZq4poVUj45qUBqS0jo63ofTcxaZYGtvgwYPnuXtahTPdvVo/QANgC5BUwbzmQELw/jBgWVW2OWDAAA/VrFmzQl432mlslZu9JNdP+PNM7zp6mt/71nzP31McnsLCRM9dbIqlsQFzvZJMDcfhzlACR/ObKngR2Vnu/nQze9rMWrv7ljDsV4QdBUXc/+4iJs7L4fA2Tfn3DceT1lVNyETKC0fQXwq8WdEMM2sHbHJ3N7OBQB1gaxj2KcL78zdw19QFbC8o4ubBPbj5tB4HbkKWOQFmjoW8HEhMhtPvgdSLD03BIhFSraA3sybAmcAN5abdCODuzwIXATeZWQlQCFwS/C+GSMhyd+3h3rcW8N78jfTr0JxXrz2Wfh2q0IQscwK8cwsUFwYe560NPAaFvcS1agW9uxcAh+037dly958Entx/PZFQuDsT5+Vw37SF7CkpY9SQPvzm5G7Uq2oTsplj/xPy+xQXBqYr6CWORcclCSIHsHZbAXdMyeKTZVsY2LUVf76wP4e3OcirIfJyDm66SJxQ0EtUKytzXvsim0dmLMGA+4b34/JBXUJrQpaYHDhdU9F0kTimoJeotTx3F6MmZTFv9XZO7dWGBy/oT8cWjUPf4On3/PgcPUD9xoHpInFMQS9Rp7i0jHFzVvLER8to0rAuf7v4SEYcHYYmZPvOw+uqG6llFPQSVeavy+OPEzNZtGEn56S2Z8x5/WjTrGH4dpB6sYJdah0FvUSFolLnofcW8/wnKzmsaQOeu2IAZ/drF+myROKCgl4i7utV27jns0I2FqzgV2mduOOcI0hsXD/SZYnEDQW9REz+3hIefm8x/+/L1bRpbLxx/SBO7NE60mWJxB0FvUTErCW53Dk5iw0793Dtid0Y1HiTQl6khijo5ZDavruI+6YtZPJ36+jZNoFJN53AMZ1bkpGRG+nSROKWgl4OCXfn3awN3PvWAvIKi7nltB789rQeNKx3gCZkIlJtCnqpcZt27uHuqfP5YOEmUpMTef36QRzRvnmkyxKpNRT0UmPcnQlz13L/u4soKinjjmF9uPbEg2hCJiJhoaCXGrFmawG3T8nks+VbGdStFQ9fmErX1k0jXZZIraSgl7AqLXNe+Tybv85YQt06xgMjUrj02M6hNSETkbBQ0EvYLN20i9smZvL92h2c1qctD4xIoX1iNZqQiUhYKOil2opKynh29gr+8fEyEhrW44lLjuIXR3aofhMyEQkLBb1Uyw9rdzBqUiaLN+7iF0d24N7z+nJYQhibkIlItSnoJSSFRaU89tFSXvhkJW2bNeKFK9M4o29SpMsSkQoo6OWgfbFiK7dPziR7awGXDuzM7cP60LyRmpCJRKtqBb2ZZQO7gFKgxN3T9ptvwBPAMKAAuNrdv63OPuUQyJxQ4Zdz7NxTzEPvLeafX62hy2FN+OdvBnHC4epPIxLtwnFEP9jdt1QybyjQM/gzCHgmeCvRKnPCj79uL28tvHMLC9bncd233cndtYffnNyN/z2zN40bqH2BSCyo6VM3w4HX3N2BL82shZm1d/cNNbxfCdXMsT/+TlWA4kISP3+IxJYv8+wVAziqU4uIlCYiobFABoe4stkqYDvgwHPuPm6/+dOAh9z90+DjmcAod59bwbZGAiMBkpKSBowfPz6kmvLz80lISAhp3Wh3KMZ2asb5GD/9m3CMmadMoV4NffApnp83iO/xaWzRYfDgwfP2P32+T3WP6E909/Vm1hb40MwWu/uccvMrSoUKX1mCLxLjANLS0jw9PT2kgjIyMgh13Wh3SMb2XXLgdM1+LDGZM04bXGO7jefnDeJ7fBpb9KtWdyl3Xx+8zQWmAAP3WyQH6FTucTKwvjr7lJpTVuZ81uW3FHqDH8+o3zjwhqyIxKSQg97MmppZs333gbOA+fst9jZwpQUcB+Tp/Hx0yt6ym8te+JLLv+7MCy3/h5JmHQGDxE5w3t8h9eJIlygiIarOqZskYErwY+71gH+6+/tmdiOAuz8LTCdwaeVyApdXXlO9ciXcSsuclz5dxaMfLqF+nTo8dEF/fnXsMMzuiHRpIhImIQe9u68Ejqxg+rPl7jvw21D3ITVrycZd3DbxB37IyeOMI5K4//wU2iU2inRZIhJm+mRsLbS3pJSnZ63g6YzlNG9Un39cejTnprZXEzKROKWgr2W+W7OdUZMyWbopnxFHd+Tuc/vSqmmDA68oIjFLQV9LFBSV8OgHS3nps1W0a96Il65O47Q+akImUhso6GNBsPfMqXk5gevcg71nqurz5VsYPTmLNdsK+PVxnRk1pA/N1IRMpNZQ0Ee7cr1nDP6v9wxwwLDPKyzmz9MXMf6btXQ9rAnjRx7Hcd0Pq/GSRSS6KOijXSW9Z5g59meD/sOFm7hrahabd+3lhlO78z9n9KJRfTUhE6mNFPTRLi/noKZvyd/LmLcXMC1zA33aNeP5K9NITW5Rc/WJSNRT0Ee7xIp7z5CY/KOH7s7U79fxp3cWUrC3lFvP7MWN6YdTv261ulyISBxQ0Ee70+/5cX94+EnvmfU7CrlzShazlmzm6M4teOTCVHomNYtAsSISjRT00W7fefiZY/G8HKzcNz6VlTlvfL2Gh99bTGmZc8+5fbnqhK7UraFWwiISmxT0sSD1Yki9mNnlWqau2rKbUZMy+XrVNk7q0Zo/X9CfTq2aRLZOEYlKCvoYU1JaxgufruKxD5fSsF4dHrkolV8OSFb7AhGplII+hqzZWcqIpz8na10eZ/dL4r7hKbRtriZkIvLzFPQxYG9JKU9+vJynv9hDy6ZlPH35MQxNaaejeBGpEgV9lJu3OtCEbHluPid0qMdT151KSzUhE5GDoKCPUrv3lvDXD5bwyufZdEhszCvXHAsbFirkReSgKeij0CfLNnP75Cxythdy5fFduG1IHxIa1iNjw8JIlyYiMUhBH0XyCop5YPpCJszNoXvrpky44XgGdmsV6bJEJMYp6KPE+/M3cvdb89m2u4j/Sj+cW07vqSZkIhIWCvoIy921hzFvL2B61kb6tm/Oy1cfS0rHxEiXJSJxREEfIe7O5G/XMXbaQgqLS/nj2b0ZeUp3NSETkbALOejNrBPwGtAOKAPGufsT+y2TDrwFrApOmuzuY0PdZ7zI2V7AHVPmM2fpZgZ0acnDF6bSo21CpMsSkThVnSP6EuBWd//WzJoB88zsQ3ff/9KQT9z93GrsJ26UlTmvf7Wah99bjAN/+kU/rjiuC3XUhExEalDIQe/uG4ANwfu7zGwR0BHQNYAVWLE5n9GTMvkmezun9GrDgyNSSG6pJmQiUvPM3au/EbOuwBwgxd13lpueDkwCcoD1wB/cfUEl2xgJjARISkoaMH78+JBqyc/PJyEhek6DlJQ572cXM3V5MQ3rwqV9GnBih3ohtS+ItrGFUzyPDeJ7fBpbdBg8ePA8d0+rcKa7V+sHSADmARdUMK85kBC8PwxYVpVtDhgwwEM1a9askNcNt6ycHT7siTneZdQ0v+n1ub5pZ2G1thdNYwu3eB6be3yPT2OLDsBcryRTq3XVjZnVJ3DE/oa7T67gRWRnufvTzexpM2vt7luqs99ot6e4lL/PXMZzc1bSskkDnv31MQxJaR/pskSklqrOVTcGvAgscve/VbJMO2CTu7uZDQTqAFtD3WcsmJu9jdsmZbJy825+OSCZu87pS2KT+pEuS0Rqseoc0Z8IXAFkmdn3wWl3AJ0B3P1Z4CLgJjMrAQqBS4L/xYg7+XtL+Mv7i3nty9V0bNGY164dyCm92kS6LBGRal118ynws+8ouvuTwJOh7iNWzF66mTsmZ7E+r5Crju/KH8/uTdOG+iyaiEQHpVE17CgoYuy0hUz+dh2Ht2nKxBuPZ0AXNSETkeiioA/Re1kbuPutBewoKOLmwT24+bQeakImIlFJQX+Qcnfu4Z63FvD+go2kdGzOq9ceS78OakImItFLQV9F7s6/5+Vw/7SF7CkpY9SQPvzm5G7UUxMyEYlyCvoqWLutgDumZPHJsi0M7NqKhy7sT/c2sfFpORERBf3PKC1zXvsim0feX0Idg/vOT+HygZ3VhExEYoqCvhLLc3dx28RMvl2zg/TebXhgRH86tmgc6bJERA6agn4/xaVlPDd7BX+fuZwmDevy2K+O5PyjOobUhExEJBoo6MvJysnjjxN/YPHGXZyT2p4//aIfrRMaRrosEZFqUdATaEL22EdLeeGTVRzWtAHPXTGAs/u1i3RZIiJhUeuD/quVWxk9OYtVW3ZzybGduH3YESQ2VhMyEYkftTbod+0p5uH3F/P6l2vo1Koxb1w/iBN7tI50WSIiYVcrg37W4lzunJLFhp17uO6kbtx6Vi+aNKiVvwoRqQVqVbpt213EfdMWMuW7dfRsm8Ckm07gmM4tI12WiEiNqhVB7+5My9zAmLcXkFdYzC2n9+S3gw+nYT01IROR+Bf3Qb9p5x7unDKfjxZtIjU5kdevH8QR7ZtHuiwRkUMmboPe3fnXN2t5YPoiikrKuGNYH649UU3IRKT2icugX7O1gNGTM/l8xVYGdWvFwxem0rV100iXJSISEXEV9KVlzozsYqbMnE29OnV4cER/Ljm2k5qQiUitFjdBn1dQzFUvf833a4s4rU9bHhiRQvtENSETEanWCWszG2JmS8xsuZmNrmC+mdnfg/MzzeyY6uzv5zRvXI8uhzXhhtSGvHhVmkJeRCQo5KA3s7rAU8BQoC9wqZn13W+xoUDP4M9I4JlQ91eFenjikqM5vkM9dZoUESmnOkf0A4Hl7r7S3YuA8cDw/ZYZDrzmAV8CLcysfTX2KSIiB6k6Qd8RWFvucU5w2sEuIyIiNag6b8ZWdH7EQ1gmsKDZSAKnd0hKSiIjIyOkovLz80NeN9ppbLErnsensUW/6gR9DtCp3ONkYH0IywDg7uOAcQBpaWmenp4eUlEZGRmEum6009hiVzyPT2OLftU5dfMN0NPMuplZA+AS4O39lnkbuDJ49c1xQJ67b6jGPkVE5CCFfETv7iVmdjMwA6gLvOTuC8zsxuD8Z4HpwDBgOVAAXFP9kkVE5GBU6wNT7j6dQJiXn/ZsufsO/LY6+xARkepRhy8RkThngYPu6GJmm4HVIa7eGtgSxnKiicYWu+J5fBpbdOji7m0qmhGVQV8dZjbX3dMiXUdN0NhiVzyPT2OLfjp1IyIS5xT0IiJxLh6DflykC6hBGlvsiufxaWxRLu7O0YuIyI/F4xG9iIiUo6AXEYlzcRP0B/q2q1hmZp3MbJaZLTKzBWb2u0jXFG5mVtfMvjOzaZGuJZzMrIWZTTSzxcHn7/hI1xROZvY/wb/J+Wb2ppk1inRNoTKzl8ws18zml5vWysw+NLNlwduWkawxVHER9FX8tqtYVgLc6u5HAMcBv42z8QH8DlgU6SJqwBPA++7eBziSOBqjmXUEbgHS3D2FQM+rSyJbVbW8AgzZb9poYKa79wRmBh/HnLgIeqr2bVcxy903uPu3wfu7CIRF3HyBi5klA+cAL0S6lnAys+bAKcCLAO5e5O47IlpU+NUDGptZPaAJlbQhjwXuPgfYtt/k4cCrwfuvAucfyprCJV6CvtZ8k5WZdQWOBr6KcCnh9DhwG1AW4TrCrTuwGXg5eFrqBTNrGumiwsXd1wF/BdYAGwi0If8gslWFXdK+1urB27YRrick8RL0Vf4mq1hmZgnAJOD37r4z0vWEg5mdC+S6+7xI11ID6gHHAM+4+9HAbmL0v/4VCZ6vHg50AzoATc3s15GtSioSL0Ff5W+yilVmVp9AyL/h7pMjXU8YnQj8wsyyCZxyO83MXo9sSWGTA+S4+77/fU0kEPzx4gxglbtvdvdiYDJwQoRrCrdNZtYeIHibG+F6QhIvQV+Vb7uKWWZmBM7zLnL3v0W6nnBy99vdPdnduxJ43j5297g4KnT3jcBaM+sdnHQ6sDCCJYXbGuA4M2sS/Bs9nTh6sznobeCq4P2rgLciWEvIqvXFI9Gism+7inBZ4XQicAWQZWbfB6fdEfziF4lu/w28ETwAWUkcfcuau39lZhOBbwlcGfYdMdwywMzeBNKB1maWA9wLPARMMLPrCLyw/TJyFYZOLRBEROJcvJy6ERGRSijoRUTinIJeRCTOKehFROKcgl5EJM4p6EVE4pyCXkQkzv1/W3v88TZE1vUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# First plot our fitted model\n",
    "x = range(0,12,1)\n",
    "y = [ 0.0 ] * len(x)\n",
    "for i in range(len(x)):\n",
    "    y[i] = x[i] * fitResult[0] + fitResult[1]\n",
    "\n",
    "# Then plot our data\n",
    "x2 = [ 2, 3, 5, 7, 9 ]\n",
    "y2 = [ 4, 5, 7, 10, 15 ]\n",
    "\n",
    "plt.plot(x,y)\n",
    "plt.plot(x2, y2, linestyle=\"\", marker=\"o\")\n",
    "plt.title(\"Linear regression\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "142e83ae",
   "metadata": {},
   "source": [
    "## Simple quadratic regression\n",
    "\n",
    "One can of course use the linear least squares method to fit arbitrary polynomial functions. One can for example also apply this to quadratic functions. Let's assume our model function is\n",
    "\n",
    "$$\n",
    "y = ax^2 + bx + c\n",
    "$$\n",
    "\n",
    "And we supply a set of points that might reflect our measured values:\n",
    "\n",
    "$$\n",
    "(0,1) \\\\\n",
    "(1,3) \\\\ \n",
    "(4,7) \\\\ \n",
    "(5,8) \\\\ \n",
    "(7,20) \\\\ \n",
    "(10,25)\n",
    "$$\n",
    "\n",
    "This would lead to a (not exactly solvable) set of linear equations by simply substituting the points into the equations:\n",
    "\n",
    "$$\n",
    "1 = a*0 + b*0 + c \\\\ \n",
    "3 = a*1^2 + b*1 + c \\\\ \n",
    "7 = a*4^2 + b*4 + c \\\\ \n",
    "8 = a*5^2 + b*5 + c \\\\ \n",
    "20 = a*7^2 + b*7 + c \\\\ \n",
    "25 = a*10^2 + b*10 + c\n",
    "$$\n",
    "\n",
    "As one can see this is a simple set of linear equations as before:\n",
    "\n",
    "$$\n",
    "1 = c \\\\ \n",
    "3 = a*1 + b*1 + c \\\\ \n",
    "7 = a*16 + b*4 + c \\\\ \n",
    "8 = a*25 + b*5 + c \\\\ \n",
    "20 = a*49 + b*7 + c \\\\ \n",
    "25 = a*100 + b*10 + c\n",
    "$$\n",
    "\n",
    "The following codeblock just builds the coefficient matrix $A$ by calculating the terms $x^2$, $x$ and the $1$ as well as the constant vector $\\vec{b}$ by inserting $y$ coordinates. The function accepts an arbitrary number of points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "01875f1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def quadraticRegression(points, trace = False, plot = False, plotRange = [-10, 10], plotSteps = 1):\n",
    "    matA2 = Matrix(len(points), 3)\n",
    "    matB2 = Matrix(len(points), 1)\n",
    "    for row in range(len(points)):\n",
    "        matA2.set(row, 0, points[row][0] * points[row][0])\n",
    "        matA2.set(row, 1, points[row][0])\n",
    "        matA2.set(row, 2, 1.0)\n",
    "        matB2.set(row, 0, points[row][1])\n",
    "\n",
    "    if trace:\n",
    "        print(\"A:\")\n",
    "        print(matA2)\n",
    "        print(\"b:\")\n",
    "        print(matB2)\n",
    "\n",
    "    fitResult = leastSquares(matA2, matB2)\n",
    "\n",
    "    if trace:\n",
    "        print(\"Result:\")\n",
    "        print(fitResult)\n",
    "\n",
    "    if plot:\n",
    "        x2 = [ 0.0 ] * len(points)\n",
    "        y2 = [ 0.0 ] * len(points)\n",
    "        for row in range(len(points)):\n",
    "            x2[row] = points[row][0]\n",
    "            y2[row] = points[row][1]\n",
    "\n",
    "        x = range(plotRange[0], plotRange[1], plotSteps)\n",
    "        y = [ 0.0 ] * len(x)\n",
    "        for i in range(len(x)):\n",
    "            y[i] = x[i] * x[i] * fitResult[0] + x[i] * fitResult[1] + fitResult[2]\n",
    "\n",
    "        plt.plot(x,y)\n",
    "        plt.plot(x2, y2, linestyle=\"\", marker=\"o\")\n",
    "        plt.title(\"Quadratic regression\")\n",
    "        plt.grid()\n",
    "\n",
    "    return fitResult"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "198fdede",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameters: [0.17054737581685814, 0.7278125943022667, 2.1137956163462075]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAEICAYAAABCnX+uAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvUElEQVR4nO3dd3xUVfrH8c+T3iAhBEJCAqG30KugQkQRRQV3FbGwyKpYsPwsu+u6rt3VdVddy9pFscFiQ8SKSFBwkU4IvUMgBAIkJIGElPP7Yy7uEBOSTCa5U57365VXZu7c8hyGfHNz5p57xBiDUkop7xVgdwFKKaXqR4NcKaW8nAa5Ukp5OQ1ypZTychrkSinl5TTIlVLKy2mQK1uJSIqIGBEJcuM+7xORN9y1P08lIm1EpFBEAu2uRdlLg1z9iohcKyJrReSYiOwXkZdEJNruuqoiIiNEJMt5mTHmb8aY6+2qqbEYY3YbY6KMMeV216LspUGuTiEidwN/B/4ARANDgBTgWxEJbuRaRERs+T/q7mO78y8OpSrTIFe/EJGmwMPAbcaYr40xpcaYncB4oB1wlbXe2yLymNN2p5wVi8i9IrJNRApEZL2IXOr0WqCI/FNEckVkOzCmUg3pIvK4iCwGjgHtRWSyiGyw9rddRG601o0EvgISrS6GQhFJFJGHROQ9p32eKSI/iUieiOwRkWuraX9Vx+4qIvNE5LCIbBKR8U7rNxeRz0XkqIgsE5HHRGSR0+tGRKaKyBZgi7XsIhFZbdXyk4j0clr/TyKy12rnJhEZaS0fJCLLrePkiMgz1vJTuqWsts+xat0qIjc47fshEZklIu9Y+18nIgNO899BeRNjjH7pF8YYgNFAGRBUxWvTgfetx28Djzm9NgLIcnp+OZCI40ThCqAISLBeuwnYCCQDscACwJw8JpAO7AZ6AEFAMI6w7wAIMBxHyPar6tjWsoeA96zHbYAC4EprX82BPtW0v/Kxo4E9wGTreT8gF+hhrT/T+ooAulvrLnLanwHmWe0Mt7Y/AAwGAoFJwE4gFOhibZ9obZsCdLAe/xeYaD2OAoY4reP8b7cQeAkIA/oAB4GRTv8mxcCF1rGfAJbY/X9Ov9zzpWfkylkckGuMKavitWygRW12Yoz50BizzxhTYYz5D46z0UHWy+OBfxlj9hhjDuMIlMreNsasM8aUGcdfBV8YY7YZh4XAt8BZtWzT1cB3xpgZ1r4OGWNWn2b9X46N4xfbTmPMW1YtK4GPgcusDxh/CzxojDlmjFmP45ddZU8YYw4bY44DNwCvGmN+NsaUG2OmAyU4uq/KcQR6dxEJNsbsNMZss/ZRCnQUkThjTKExZknlg4hIMnAm8CdjTLHVxjeAiU6rLTLGfGkcfervAr1r8w+oPJ8GuXKWC8RV05+bgOMMr0Yi8jun7oM8IBXHLwlwnKnvcVp9VxW7cH4dEblARJZYXQZ5OM4q46rYrirJwLYa16r62G2BwSfbYR37aqAVjl9qQZXWP6XuavZ3d6X9JeM4C98K/B+OM+cDIjJTRBKt7a4DOgMbrS6ci6o4TiJw2BhT4LRsF9Da6fl+p8fHgDDtu/cNGuTK2X9xnCH+xnmh1Rd9AY4/3cHRVRLhtEorp3XbAq8DtwLNjTExQCaObhFwnNknO23bpoo6frklp4iE4jgL/icQb+3vS6f91XT7zj04umVqy3l/e4CFxpgYp68oY8zNOH6plQFJTus7t6u6/T1eaX8RxpgZAMaYD4wxZ+IIfIPjQ2eMMVuMMVcCLa1lH1nvibN9QKyINHFa1gbYW4e2Ky+lQa5+YYzJx/Fh5wsiMlpEgkUkBfgQx9n6+9aqq4ELRSRWRFrhOJM8KRJHCB0EEJHJOM7IT5oF3C4iSSLSDLi3hrJCcHQ5HATKROQCYJTT6zlAc6n+8sj3gXNFZLyIBFkfUPap4ZgnzQU6i8hE698iWEQGikg3q3viE+AhEYkQka7A72rY3+vATSIyWBwiRWSMiDQRkS4ico71i6sYOI6juwURuUZEWhhjKoA8a1+nXHJojNkD/AQ8ISJh1oeo1/G/90z5MA1ydQpjzFPAfTjOgAuAHTjOvs81xhRZq70LrMHxQd23wH+ctl8PPI3j7D4H6AksdjrE68A31vYrcYTh6eopAG7H8QvgCI4rZ+Y4vb4RmAFst7orEittvxtHV8zdwGEcv4Rq1TdsHXsUMAHHGe9+HGfEodYqt+L4QHQ/jn+TGTj+oqluf8tx9JO/aLVlK3Ct9XIo8CSOX5j7cZx932e9NhpYJyKFwHPABGNMcRWHuBLHB6D7gE9x9N/Pq01blXcTY3RiCVU9Efk9jrP0YVYoqmqIyN+BVsaYSXbXovyLftChTssYM01ESoGhOC7NUxarOyUEWAsMxNGV4fMjSpXn0TNypVwkIgNxdKck4rg+/FXgSaM/VKqRaZArpZSX0w87lVLKyzVqH3lcXJxJSUlxaduioiIiIytfOuvdfK1NvtYe8L02+Vp7wPfaVFV7VqxYkWuMqXZkdaMGeUpKCsuXL3dp2/T0dEaMGOHegmzma23ytfaA77XJ19oDvtemqtojIlWNgP6Fdq0opZSX0yBXSikvp0GulFJeToNcKaW8nAa5Ukp5OQ1ypZTychrkSinl5TTIlVKqAR0/Uc7Dn69jz+FjDXYMDXKllGpAH6/M4q3FO8nOr+oW8u6hQa6UUg2kosIwbdEOeiVFMzClWYMdR4NcKaUayIJNB9ieW8R1Z7ZDRGrewEUa5Eop1UDe+HEHCdFhXNgzoUGPo0GulFINIHNvPv/dfohrh6YQHNiwUatBrpRSDWDaoh1EhAQyYVCbBj+WBrlSSrnZ/vxi5qzZx/gByUSHBzf48TTIlVLKzd75707KjeH3w9o1yvE0yJVSyo2OnSjj/Z93c373VrRpHtEox9QgV0opN/p4RRb5x0u5/qzGORsHDXKllHKbigrDm4t20Ds5hv5tG24AUGUa5Eop5SbzNx5g56FjXN/AA4Aq0yBXSik3eePH7bSOCeeC1FaNelwNcqWUcoO1Wfn8vOMw1w5NIaiBBwBVpkGulFJu8Oai7USGBHLFoORGP7YGuVJK1VN2/nHmZmRzxcA2NA1r+AFAlWmQK6VUPU3/aRcVxjB5WIotx9cgV0qpeigqKeODn3cxOrUVybGNMwCoMg1ypZSqh49WZHG0uIzrzmxvWw0a5Eop5aLyCsO0xTvo26ZxBwBVpkGulFIu+m5DDrsOHeN6G8/GQYNcKaVc9uaPO2gdE875PeJtraPWQS4igSKySkTmWs9jRWSeiGyxvtv3d4VSSjWyNXvyWLrzMJOHNf4AoMrqcvQ7gA1Oz+8F5htjOgHzredKKeUX3li0g6jQIK4Y2PgDgCqrVZCLSBIwBnjDafFYYLr1eDowzq2VKaWUh9qZW8QXGfu4anAbmtgwAKiy2p6R/wv4I1DhtCzeGJMNYH1v6d7SlFLKM736wzaCAgO4/szGu+f46Ygx5vQriFwEXGiMuUVERgD3GGMuEpE8Y0yM03pHjDG/6icXkSnAFID4+Pj+M2fOdKnQwsJCoqKiXNrWU/lam3ytPeB7bfK19kDjt+lIcQX3LDzO8KQgftcj1O37r6o9aWlpK4wxA6rdyBhz2i/gCSAL2AnsB44B7wGbgARrnQRgU0376t+/v3HVggULXN7WU/lam3ytPcb4Xpt8rT3GNH6bHv18nWn/5y/M7kNFDbL/qtoDLDenydYau1aMMX82xiQZY1KACcD3xphrgDnAJGu1ScBntf6Vo5RSXuhI0Qne/3k3l/ROtG04flXqc83Mk8B5IrIFOM96rpRSPuutn3ZyvLScm0d0sLuUUwTVZWVjTDqQbj0+BIx0f0lKKeV5CkvKeHvxDkZ1j6dzfBO7yzmFjuxUSqlaeH/JLo4Wl3FLWke7S/kVDXKllKpBcWk5byzawZkd4+iTHGN3Ob+iQa6UUjX4aEUWBwtKuMXD+sZP0iBXSqnTKCuv4JWF2+iTHMMZHZrbXU6VNMiVUuo0Ps/YR9aR40xN64iI2F1OlTTIlVKqGhUVhpcWbKNLfBNGdvXcu5BokCulVDXmbchhy4FCbknrQECAZ56Ngwa5UkpVyRjDS+nbaBMbwZieCXaXc1oa5EopVYWfth1izZ48bhze3vaJI2ri2dUppZRN/r1gKy2bhPLbfkl2l1IjDXKllKpk5e4j/LTtEDec1Z6w4EC7y6mRBrlSSlXy0oJtRIcHc9XgNnaXUisa5Eop5WTT/gK+25DD5GEpRIbW6b6CttEgV0opJy+nbyUiJJBrh6bYXUqtaZArpZRl16Ei5qzZx9WD2xATEWJ3ObWmQa6UUpbn5m8hODCAG85qb3cpdaJBrpRSwNYDhcxetZeJQ9rSsmmY3eXUiQa5UkoBz8/fQlhwIDd56K1qT0eDXCnl9zbtL+DzjH1MGppCXFSo3eXUmQa5Usrv/eu7zUSGBDHFy/rGT9IgV0r5tXX78vkqcz+/H5ZCs0jvuVLFmQa5UsqvPTtvC03DgrjOS8/GQYNcKeXH1uzJ47sNOdxwVnuiw4PtLsdlGuRKKb/1zLzNxEQEM/nMdnaXUi8a5Eopv7Ri12EWbj7IjWd3IMpL7qlSHQ1ypZRfembeZuKiQpg0tK3dpdSbBrlSyu8s2X6IxVsPcdPwDkSEePfZOGiQK6X8jDGGZ77dTMsmoVwzxPvPxkGDXCnlZxZvPcTSnYeZmtbRK2b/qQ0NcqWU3zDG8PS8TSRGhzFhULLd5biNBrlSym+kbzrIqt153HpOJ0KDfONsHDTIlVJ+whjDM/M2kxwbzuUDkuwux600yJVSfmHe+hzW7s3n9nM6ERzoW9HnW61RSqkqVFQ4zsbbxUVyad/WdpfjdhrkSimf91XmfjbuL+COkZ0I8rGzcdAgV0r5uNLyCp6et4mOLaO4uHei3eU0iBqDXETCRGSpiKwRkXUi8rC1PFZE5onIFut7s4YvVyml6mbmsj1sP1jEn0Z3JTBA7C6nQdTmjLwEOMcY0xvoA4wWkSHAvcB8Y0wnYL71XCmlPEZhSRnPfbeZQe1iObdbS7vLaTA1BrlxKLSeBltfBhgLTLeWTwfGNUSBSinlqlcXbiO38AR/ubAbIr55Ng4gxpiaVxIJBFYAHYF/G2P+JCJ5xpgYp3WOGGN+1b0iIlOAKQDx8fH9Z86cWeciK4xhZ24R7VtE1XlbT1ZYWEhUlO+0ydfaA77XJl9rD1TfpiPFFfzph+P0bRnIzX3CbKjMNVW1Jy0tbYUxZkC1Gxljav0FxAALgFQgr9JrR2ravn///sYV98xabXo/8IUpKil1aXtPtWDBArtLcCtfa48xvtcmX2uPMdW36Z5Zq02n+740uw8VNW5B9VRVe4Dl5jTZWqerVowxeUA6MBrIEZEEAOv7gbrsqy4mDGpDXonh1YXbG+oQSikfsiH7KB+tzGLS0LYkx0bYXU6Dq81VKy1EJMZ6HA6cC2wE5gCTrNUmAZ81UI30b9uMQa0CefWHbezPL26owyilfMQTX22kaVgwt6Z1sruURlGbM/IEYIGIZADLgHnGmLnAk8B5IrIFOM963mAu7xxChYF/fLOpIQ+jlPJyP245yA+bD3LbOR2JjvDeCZXrosapMYwxGUDfKpYfAkY2RFFVaRERwO+HteOVhdu4dmgKPZOiG+vQSikvUV5h+NuXG0mODWfiGb4xaURteNXIzlvSOtA8MoTHvlh/8gNWpZT6xaer9rIh+yh/OL+rT92mtiZeFeRNw4K587zO/LzjMN+uz7G7HKWUBykuLefpbzfROymai3sl2F1Oo/KqIAeYMDCZTi2jeOLLDZwoq7C7HKWUh3hz0Q6y84u5z8cH/1TF64I8KDCAv4zpxs5Dx3h3yS67y1FKeYBDhSW8nL6Nc7vFM7h9c7vLaXReF+QAI7q05OzOLXjuu80cKTphdzlKKZs9P38Lx0vLufeCrnaXYguvDHKA+8d0c9wQZ/4Wu0tRStlof1EF7/+8mwkDk+nY0rduP1BbXhvkneObcOWgNry3ZBfbDhbWvIFSyid9uPkEoUEB/N+5ne0uxTZeG+QAd57XmbDgQJ74cqPdpSilbLBs52FW5JRz4/AOtGgSanc5tvHqII+LCmVqWke+25DDT1tz7S5HKdWIjDH87csNxIQK15/Vzu5ybOXVQQ4weVgKrWPCeeyLDZRX6CAhpfzFnDX7WLU7j990CiYipMZB6j7N64M8LDiQey/oyvrso3y8MsvucpRSjaCguJTHv9hAr6Rozmzt3yEOPhDkABf1SqBvmxj+8c0mikrK7C5HKdXAnp+/hYOFJTwyNpUAPxv8UxWfCHIR4a8XdedgQQmvLtxmdzlKqQa0OaeAtxbv5IoByfRJjrG7HI/gE0EO0K9NMy7unchrP25nX95xu8tRSjUAYwwPfJZJZGgQ97fJhGdTGZ4+Dp5NhYxZdpdnG58JcoA/je6CMfD4FxvsLkUp1QA+z8hmyfbDvNhzK1Hf3gX5exAM5O+Bz2/32zD3qSBPahbBrWkd+WJtNgs3H7S7HKWUGxWWlPH4F+tJbd2UM3e9BKWV/vIuPQ7zH7GnOJv5VJADTBnenvZxkTz4WSbFpeV2l6OUcpMX5m8h56jjA07Jr+YKteqW+zifC/LQoEAeHZfKzkPHeDldP/hUyhdsPVDAm4t2MH5AEv3aNIPopKpXrG65j/O5IAcY1jGOS3on8nL6NnbkFtldjlKqHhwfcK4jIiSQP4227m448gEIDj91xeBwx3I/5JNBDnD/Rd0IDQrggc8ydVo4pbzYF2uz+WnbIe45vwvNo6z7qfQaDxc/D9HJGASikx3Pe423t1ib+GyQt2wSxj3nd+HHLbl8sTbb7nKUUi4oKinjsbkb6J7QlKsHV5pMudd4uDOThSNmw52Zfhvi4MNBDnDNkLaktm7KI5+vp6C41O5ylFJ19ML3W9l/tJhHx/UgMEBHcFbHp4M8MEB4fFxPDhaW8My8zXaXo5Sqg60HCnlz0XYu659E/7axdpfj0Xw6yAF6J8dw9eA2TP9pJ5l78+0uRylVC8YYHpqz7peb4qnT8/kgB/jD+V2JjQzh/tmZVOitbpXyeF9l7mfR1lzuPq8zcVH+O2FEbflFkEeHB/OXMd1YvSePGct2212OUuo0jp0o47G56+naqgnXDGlb8wbKP4IcYFyf1pzRvjlPfb2J3MISu8tRSlXjhe+3si+/mEfHpRIU6DcRVS9+868kIjw6rgfHTpTpHJ9KeajMvfm89oPjA86BKfoBZ235TZADdGzZhBvOas/HK7P4efshu8tRSjkpLa/gDx9lOD7PGtPN7nK8il8FOcBt53QiqVk498/O5ERZhd3lKKUsr6RvY0P2UR4bl0pMRIjd5XgVvwvy8JBAHr6kB1sOFPLmoh12l6OUwjHrzwvfb2VMrwTO79HK7nK8jt8FOcDIbvGM6h7Pc/M36021lLJZeYXhjx9lEBnqOMlSdeeXQQ7w6LhUQgID+MOHayjXa8uVss1bi3ewek8eD13SQ68Zd5HfBnl80zAeHtuD5buO8NZi7WJRyg47cov4xzebOLdbSy7pnWh3OV7Lb4McHNeWn9c9nn98s4mtBwrtLkcpv1JRYfjTxxmEBAXw2LieiOhNsVxVY5CLSLKILBCRDSKyTkTusJbHisg8EdlifW/W8OW6l4jw+KWphIcEco92sShVNxmzHLPXPxTj0iz27y/dzdIdh7l/TDdaRYc1TI1+ojZn5GXA3caYbsAQYKqIdAfuBeYbYzoB863nXqdlkzAeGZvK6j15vP7jdrvLUco7ZMxyzFqfvwdcmMU+68gxnvxyA2d2jGP8gOSGrdUP1BjkxphsY8xK63EBsAFoDYwFplurTQfGNVCNDe7iXglckNqKZ77dzJacArvLUcrzzX/E5VnsjTH8+ZO1GOCJ32iXijtIXaZBE5EU4AcgFdhtjIlxeu2IMeZX3SsiMgWYAhAfH99/5syZLhVaWFhIVFSUS9vWxtESw18WHSMuPID7h4Q1yk3sG7pNjc3X2gO+1yZ3tWd4+jiEX2eHQRwz9pzGj1mlvJl5gmu6hXBu2+B61+IP71FaWtoKY8yA6rYJqu3ORSQK+Bj4P2PM0dr+FjXGvAa8BjBgwAAzYsSI2h7yFOnp6bi6bW0FJmQz9YOVbJRkpo7o2KDHgsZpU2PytfaA77XJbe1ZlWR1q5xKopNOu/8DR4u5PX0hA1Oa8cjEMwhwwwmTvke1vGpFRIJxhPj7xphPrMU5IpJgvZ4AHKjTkT3QmF4JXNQrgX99t5mN+4/aXY5SnsuFWeyNMfxldiYlZRX8/be93BLiyqE2V60I8CawwRjzjNNLc4BJ1uNJwGfuL6/xPTI2lejwYO6etYbScr0Xi1JVcprFnlrOYj83I5t563O467zOtG/hO10hnqA2XSvDgInAWhFZbS27D3gSmCUi1wG7gcsbpMJGFhsZwmPjenLTeyt4acE27ji3k90lKeWZeo2v9cz1uYUlPDhnHb2TornuzHYNXJj/qTHIjTGLgOr+Bhrp3nI8w+jUVozrk8gL32/h3O4t6ZEYbXdJSnktYwx/+HANhSVlPHVZb50sogHov2g1HrqkB80iQ7jnwwy93a1S9fDW4p0s2HSQv1zYjS6tmthdjk/SIK9GTEQIT1zakw3ZR3lxwVa7y1HKK63bl8+TX23k3G4t+d0ZOv9mQ9EgP41zu8fz235J/HvBVjKy8uwuRymvcuxEGbfPWEVMRDBPXdZbB/40IA3yGjxwcXfim4Ry24xVFBSX2l2OUl7j0bnr2Z5bxLNX9CE2Umf8aUga5DWIDg/m+Sv7knXkuGNYcR1Gwirlr75cm82MpXu48ewODOsYZ3c5Pk+DvBYGpMRy96jOzM1w/OdUSlVvb95x7v04g95J0dw9qrPd5fgFDfJauunsDpzduQUPf76ODdk66lOpqpRXGO6cuZryCsPzV/YlWC81bBT6r1xLAQHCM+N7Ex0ezNQPVlJUUmZ3SUp5nBe/38rSnYd5dFwqbZtH2l2O39Agr4O4qFCem9CXnblF/PWzTLvLUcqjLN95mOfmb2Zcn0R+0y/J7nL8igZ5HZ3RoTm3j+zEJyv38tGKLLvLUcoj5B8v5Y6Zq0lqFsGj41LtLsfvaJC74LZzOnFG++b8dXamTkSh/J4xhvs+XUvO0WKem9CHJmH1v8e4qhsNchcEBgjPTehDREggt36wiuMnyu0uSSnbfLg8iy8ysrlrVGf6tvG6qXt9gga5i1o2DePZK/qw+UABj8xdZ3c5Stli28FCHpyzjqEdmnPT2R3sLsdvaZDXw9mdW3DLiA7MWLqHz1bvtbscpRpVUUkZU99fSVhwAM9e0UcnirCRBnk93XluZwa0bcZ9n6xlR26R3eUo1SiMMfzhozVszinguQl9iW8aZndJfk2DvJ6CAgMcAx+CApj6/kqKS7W/XPm+l9K38eXa/dx7QVfO7tzC7nL8nga5GyTGhPP05b1Zn32Ux7/YYHc5SjWo7zfm8M9vNzG2TyI3nNXe7nIUGuRuM7JbPFPObs+7S3Yxc+luu8tRqkFsO1jIHTNW0z2hKU/+ppfemtZDaJC70R/P78LZnVtw/+xMlmw/ZHc5SrlVQXEpU95ZTnBQAK9O7E94SKDdJSmLBrkbBQUG8OJVfWnbPIKb31vB7kPH7C5JKbeoqDDc+Z/V7Dx0jJeu7kdSswi7S1JONMjdrGlYMG9OGogBrpu+jKM6GYXyAf/6bjPfbTjAAxd1Z0j75naXoyrRIG8AKXGRvHR1P3bkFnH7jFWUV+hkFMp7fZ2ZzfPfb2X8gCSdd9NDaZA3kKEd4nh4bA/SNx3kb1/qlSzKO23aX8Bds9bQJzmGR8am6oebHirI7gJ82dWD27Ilp5A3F+2gc3wUVwxsY3dJStVa3rET3PDOciJDg3h1Yn/CgvXDTU+lZ+QN7P4x3TirUxz3z87kZ72SRXmJsvIKbpuxiuz847xyTX8duenhNMgbmONKln4kx0Zwk17JorzEP77ZxI9bcnl0bCr92+odDT2dBnkjiA53XMlSYRxXshTolSzKg81YuptXf9jOxCFtmTBIuwO9gQZ5I2kXF8nLV/dju17JojzY15nZ/OXTtYzo0oIHLu5udzmqljTIG9HQjnE8dEkPFmw6yJNf6ZUsyrMs2X6I22eupndyDC9d3Y/gQI0Hb6FXrTSyiUPasiWngNd/3EFJ9xBG2F2QUsC6ffncMH05bWIjmDZpIBEhGg3eRH/l2uCBi7ozsmtL3l1/gk9X6QTOyl67DhUxadoymoQF8c7vB9EsMsTuklQdaZDbICgwgH9f3Y+usQHc82EG36zbb3dJyk8dKCjmd9OWUlZRwTvXDSIxJtzukpQLNMhtEhYcyO39wujZOprbPljFj1sO2l2S8jNHi0u5dtoyDhwt4a1rB9KxZRO7S1Iu0iC3UXiQ8PbkgbRvEcmUd1awfOdhu0tSfqK4tJwp7yxnc04Br0zsT982eq24N9Mgt1lMRAjvXjeYVtFhTH5rGZl78+0uSfm48grD/81czZLth3l6fG+G61RtXq/GIBeRaSJyQEQynZbFisg8Edlifddf53WRMQueTWV4+jh4NpUWOz7jvesH0zQ8mN9NW8rWAwUu75OHYhzfM2a5u2rlA4wx3D87k6/X7efBi7sztk9ru0tSblCbM/K3gdGVlt0LzDfGdALmW89VbWTMgs9vh/w9CAby98Dnt9N69+e8d/1gAkS4+o2f2XO4DkP5nfaJ0z41zFVln24tZcbS3UxN68DkYe3sLke5SY1Bboz5AajceTsWmG49ng6Mc29ZPmz+I1B6/NRlpcdh/iO0i4vkvesHUVxawVVvLGF/fnG996nUSW8u2sGcbaVMGJjMPaO62F2OciMxpuah4iKSAsw1xqRaz/OMMTFOrx8xxlTZvSIiU4ApAPHx8f1nzpzpUqGFhYVERUW5tK0nGZ4+znEmXolBWDhiNgDb88p5alkxsWHCnweH0yTk9PeArs0+G4OvvEfOfKVNc7ef4KPNpfRpbritfySBAb5zX3FfeY9Oqqo9aWlpK4wxA6rbpsGHbxljXgNeAxgwYIAZMWKES/tJT0/H1W09yqokqwvkVBKd9Ev7RgDdeh7i2reW8uqmID64YQhNw4Lrtc/G4DPvkRNvb5Mxhqe/3cxHm7cyrk8iF7fMY+Q5aXaX5Vbe/h5V5kp7XL1qJUdEEgCs7wdc3I//GfkABFcadBEc7lju5IwOzXnlmv5szC7g2mlLyT9+mjsm1nKfyr8YY3hk7npeXLCVKwcl8/T4Pj51Jq7+x9UgnwNMsh5PAj5zTzl+oNd4uPh5iE7GIBCd7Hjea/yvVk3r2pIXruzL2r35XPHqfzlwtJo+c6d9UsM+lX8orzD8+ZO1vLV4J9ed2Y6/XdpTQ9yH1di1IiIzcPy1HyciWcCDwJPALBG5DtgNXN6QRfqcXuOh13gW1uJPqAt6JjAtLIgb313Bb17+iXevG0y7uMhq96lUaXkFd89aw5w1+7j9nI7ceV5nnWvTx9XmqpUrjTEJxphgY0ySMeZNY8whY8xIY0wn67sOSWxAZ3VqwYwbhlBUUsblr/ykg4ZUtYpLy7nl/ZXMWbOPey/oyl2jumiI+wEd2ekleifH8OFNQwkNCmTCa0v4aVuu3SUpD3PsRBk3vLOceetzeHRsD24a3uF/L1YahKZjDHyLBrkX6dgyio9uPoPEmDCunbaMrzOz7S5JeYiC4lImTVvK4q25/PPy3kw8I+V/L1YzCE3D3HdokHuZhOhwZt14Bqmtm3LL+yv54OfddpekbHak6ARXv/Ezq3bn8cKV/bisf9KpK+iAMZ+nQe6FYiJCeP/6IQzv3IL7Pl3Li99voTYDu5Tv2Zt3nAmvLWHj/gJendifMb0Sfr1SfjWTl1S3XHkdDXIvFR4SyGu/G8ClfVvzz2838/Dn66nQCZ39ytIdh3npub/xVt5kNgVdycivR1bdXRKd9Otlp1uuvI5OzOfFggMDePry3sRGhvDmoh0cLjrBPy/vTUiQ/n72de8t2cWKua/yRNAbhEmJY+HJvm849VLUkQ84ljt3r+iAMZ+iQe7lAgKE+8d0Iy4qlL9/vZGsI8f499X9SIjWKbt80YmyCh6cs44ZS3ezPPJDwspLTl3hZN+3c5CffDz/EUx+FhKd5AhxHXfgM/TUzQeICDeP6MC/r+rHpv0FXPT8IhZv1csTfc2BgmKuen0JM5bu5pYRHWheXs30gFX1ffcaD3dmOm6idmemhriP0SD3IWN6JfDZrWcSGxnCxDd/5sXvt2i/uY/IyMrjkhcWs27fUV68qi9/HN3VcWZdFe379jsa5D6mY8soZk8dxsW9E/nnt5u5/p3l5B07YXdZqh4+WZnFZa/8l8AA4eObh3JRr0THC3qzNGXRIPdBkaFB/OuKPjw6tgc/bjnIRS8sYm2WDuv3NmXlFTw2dz13zVpD/zbN+Py2M+me2PR/K+jN0pRFP+z0USLCxDNSSG0dzdT3V/Lbl3/ioUt6cOWgZL33hhfIO3aCWz9YxaKtuVw7NIW/jOlGcGAV5116szSFnpH7vL5tmjH39rMY3D6W+z5dyz0fZnD8RLndZanTWLDxAOf/6weW7jjMU5f14qFLelQd4kpZ9H+HH4iNDOHtyYO4Y2QnPlmVxaUvLWZHbpHdZalK8o+Xcs+Ha5j89jKiw4P5+OahjB+QbHdZygtokPuJwADhzvM689a1A9l/tJgLn/uR13/YTll5hd2lKeD7jTmMenYhn67ay61pHfn8tjPpmRRtd1nKS2iQ+5kRXVry5e1nMbRDcx7/cgOXvLiYNXvy7C7Lb+UfK+XuWWv4/dvLiQkPYfYtw7jn/C6EBgXaXZryIhrkfigxJpw3Jg3g5av7kVtYwriXFvPQnHUUFJ9mXlDldt9vzGHUvxYye/VebjunI3NuG6Zn4coletWKnxIRLuiZwLBOcfzj601M/+9Ovs7cz8Nje3B+j1Z2l+fT8o+V8vDcdXyyci9dWzXhjd8N1ABX9aJn5H6uaVgwj45L5eObhxITEcyN767ghneWsy/veM0bqzr7bn0O5z27kM9WO+bTnHOr9oWr+tMgVwD0swac3HtBV37ccpDznlnItEU7KNch/m6xek8e17zxM9e/s5zYyBA+mzqMu0Z10TtVKrfQrhX1i+DAAG4a3oExPRO4f3Ymj8xdz+zVe/nzBd0Y0j5WBxK5YOP+ozz97Wbmrc+heWQIf72oOxOHtNUAV26lQa5+JTk2grcnD2RuRjaPzl3Pla8vYUDbZkw9pyMjOrfQQK+FnblFPPvdZuas2UdUaBD3jOrM5GHtiAzVHznlfvq/SlVJRLi4dyLndY/nw+V7eGXhdia/tYzU1k25Na0To7rHExCggU7GLMf9v/OzIDqJI2fcy1P7ejFreRYhgQHcPLwDN57dgeiIYLsrVT5Mg1ydVlhwIBPPSOGKgW2YvWovL6Vv5ab3VtA5PoqpaR0Z0zOBIH8dPn5ydvqTM+/k7yH0qzspKZ/CxCFXcUtaB1o2CbO3RuUX/PQnUNVVSFAA4wcm891dw3luQh8A7pi5mnOfWcisZXs4UeaHI0SrmJ0+Qk7wVLPZPHRJDw1x1Wg0yFWdBAUGMLZPa76+42xeuaY/UWFB/PHjDEb8YwFf7ygl52ix3SU2uLxjJ5j+004qqpmFPqhgbyNXpPyddq0olwQECKNTW3F+j3jSNx/k399vZeamI/znifkMTInl4l4JXNAzgbioULtLdYuy8gp+3JLLhyv28N36A5wor2B0RBzxFVVMt6Yz9KhGpkGu6kVESOvSkrQuLflg7vccCEtibkY2f/1sHQ/OWccZHZpzUa9ERvdoRbPIELvLrbN9hRU8+dVGPlmZxYGCEmIjQ7h6SBsu759MfO7fdHZ65RE0yJXbJEYFcNWIztwxshObcgqYuyabuRn7+PMna/nr7EyGdYzjol4JjOrRiuhwz7yKI/94KSt3HWHZzsMs3naINXuOExiwnbQuLbisfxLndI3/3zXgif+bnf7kVSs6O72ygwa5cjsRoWurpnRt1ZS7R3Vm3b6jfJ6xj7lrsvnDRxnc9+lauiU0pUdiND1bR5Pauimd45sQFtz4d/zbl3ecZTsPs3ynI7w35RRgDAQFCD1aR3NFlxDuvuys6j+41Bl6lAfQIFcNSkRIbR1Nauto7h3dldV78vhmXQ5r9+bxRcY+ZizdDTiCs1N8E1ITm9IzKZoeidF0T2hKeEj9w90Yw9HiMnILSzhYUMKWA4Ust8J7r3VPmciQQPq1bcaFPRMYkNKMPskxRIQEkZ6erlefKI+nQa4ajYjQt00z+rZpBjgCNuvIcTL35rN2bz6Z+44yf+MBPlzhuBokQCAhOpwmYUFEhgYR5fQVGRpEVFgQUaGBRIUGEx4SQP6xUnILT3CwoITcwpJfgju36MSvLo9s2SSUgSmx3HBWOwakxNK1VZNTr4e3BvoMz8+CVdplojybBrmyjYiQHBtBcmwEF/RMABzhnp1fTObefDL35pN15DiFJWUUnSgj79gJso4cczwvKaewpOxX+wwQaB4VSlxUKC2ahNKhZRQtrMdx1vI2sREkx4ZXf6sBp4E+ApC/x/EcNMyVR9IgVx5FREiMCScxJpxRNdwXvaLCcKy0nKKSMopKymgaHkyziBAC63vrgCoG+lB63LFcg1x5IA1y5bUCAuSXrha3qmagT7XLlbJZvUZ2ishoEdkkIltF5F53FaWUraob0KMDfZSHcjnIRSQQ+DdwAdAduFJEururMKVsM/IBx8AeZzrQR3mw+pyRDwK2GmO2G2NOADOBse4pSykb9RoPFz8P0ckYBKKTHc+1f1x5KDHGtam8ROQyYLQx5nrr+URgsDHm1krrTQGmAMTHx/efOXOmS8crLCwkKirKpW09la+1ydfaA77XJl9rD/hem6pqT1pa2gpjzIDqtqnPp0RVXRrwq98KxpjXgNcABgwYYEaMGOHSwdLT03F1W0/la23ytfaA77XJ19oDvtcmV9pTn66VLCDZ6XkSsK8e+1NKKeWC+gT5MqCTiLQTkRBgAjDHPWUppZSqLZe7VowxZSJyK/ANEAhMM8asc1tlSimlaqVeIymMMV8CX7qpFqWUUi5w+aoVlw4mchDY5eLmcUCuG8vxBL7WJl9rD/hem3ytPeB7baqqPW2NMS2q26BRg7w+RGT56S6/8Ua+1iZfaw/4Xpt8rT3ge21ypT06+bJSSnk5DXKllPJy3hTkr9ldQAPwtTb5WnvA99rka+0B32tTndvjNX3kSimlquZNZ+RKKaWqoEGulFJezuODXEQuF5F1IlIhIgOclqeIyHERWW19vWJnnbVVXXus1/5sTdKxSUTOt6vG+hCRh0Rkr9P7cqHdNbnCFydNEZGdIrLWel+W211PXYnINBE5ICKZTstiRWSeiGyxvjezs8a6qqZNdf4Z8vggBzKB3wA/VPHaNmNMH+vrpkauy1VVtsealGMC0AMYDbxkTd7hjZ51el+8buSvj0+akma9L9543fXbOH42nN0LzDfGdALmW8+9ydv8uk1Qx58hjw9yY8wGY8wmu+twl9O0Zyww0xhTYozZAWzFMXmHanw6aYoHMsb8AByutHgsMN16PB0Y15g11Vc1baozjw/yGrQTkVUislBEzrK7mHpqDexxep5lLfNGt4pIhvVno1f9qWvxpffCmQG+FZEV1oQvviDeGJMNYH1vaXM97lKnnyGPCHIR+U5EMqv4Ot1ZUDbQxhjTF7gL+EBEmjZOxafnYntqNVGHJ6ihfS8DHYA+ON6jp+2s1UVe817U0TBjTD8cXUZTReRsuwtSVarzz1C97n7oLsaYc13YpgQosR6vEJFtQGfA9g9xXGkPXjRRR23bJyKvA3MbuJyG4DXvRV0YY/ZZ3w+IyKc4upCq+uzJm+SISIIxJltEEoADdhdUX8aYnJOPa/sz5BFn5K4QkRYnPwwUkfZAJ2C7vVXVyxxggoiEikg7HO1ZanNNdWb9MJ10KY4Pd72Nz02aIiKRItLk5GNgFN753lQ2B5hkPZ4EfGZjLW7hys+QR5yRn46IXAq8ALQAvhCR1caY84GzgUdEpAwoB24yxtT7Q4OGVl17jDHrRGQWsB4oA6YaY8rtrNVFT4lIHxxdETuBG22txgU+OmlKPPCpiIDj5/4DY8zX9pZUNyIyAxgBxIlIFvAg8CQwS0SuA3YDl9tXYd1V06YRdf0Z0iH6Sinl5by2a0UppZSDBrlSSnk5DXKllPJyGuRKKeXlNMiVUsrLaZArpZSX0yBXSikv9//5HA9/mVAlPQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "res = quadraticRegression(\n",
    "    [\n",
    "        [ -10, 10 ],\n",
    "        [ -8, 10 ],\n",
    "        [ 0, 1 ],\n",
    "        [ 1, 3 ],\n",
    "        [ 4, 7 ],\n",
    "        [ 5, 8 ],\n",
    "        [ 7, 20 ],\n",
    "        [ 10, 25 ]\n",
    "    ],\n",
    "    plot = True,\n",
    "    plotRange = [ - 15, 15 ]\n",
    ")\n",
    "\n",
    "print(\"Parameters: {}\".format(res))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "412741a4",
   "metadata": {},
   "source": [
    "## Simple cubic fit\n",
    "\n",
    "As mentioned before this can be applied to any problem that reduces to a set of linear equations - for example also a cubic fitting problem - with minor modifications:\n",
    "\n",
    "$$\n",
    "f(x) = a*x^3 + b*x^2 + c*x + d\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "2c2ea734",
   "metadata": {},
   "outputs": [],
   "source": [
    "def cubicFit(points, trace = False, plot = False, plotRange = [-10, 10], plotSteps = 1):\n",
    "    matA2 = Matrix(len(points), 4)\n",
    "    matB2 = Matrix(len(points), 1)\n",
    "    for row in range(len(points)):\n",
    "        matA2.set(row, 0, points[row][0] * points[row][0] * points[row][0])\n",
    "        matA2.set(row, 1, points[row][0] * points[row][0])\n",
    "        matA2.set(row, 2, points[row][0])\n",
    "        matA2.set(row, 3, 1.0)\n",
    "        matB2.set(row, 0, points[row][1])\n",
    "\n",
    "    if trace:\n",
    "        print(\"A:\")\n",
    "        print(matA2)\n",
    "        print(\"b:\")\n",
    "        print(matB2)\n",
    "\n",
    "    fitResult = leastSquares(matA2, matB2)\n",
    "\n",
    "    if trace:\n",
    "        print(\"Result:\")\n",
    "        print(fitResult)\n",
    "\n",
    "    if plot:\n",
    "        x2 = [ 0.0 ] * len(points)\n",
    "        y2 = [ 0.0 ] * len(points)\n",
    "        for row in range(len(points)):\n",
    "            x2[row] = points[row][0]\n",
    "            y2[row] = points[row][1]\n",
    "\n",
    "        x = range(plotRange[0], plotRange[1], plotSteps)\n",
    "        y = [ 0.0 ] * len(x)\n",
    "        for i in range(len(x)):\n",
    "            y[i] = x[i] * x[i] * x[i] * fitResult[0] + x[i] * x[i] * fitResult[1] + x[i] * fitResult[2] + fitResult[3]\n",
    "\n",
    "        plt.plot(x,y)\n",
    "        plt.plot(x2, y2, linestyle=\"\", marker=\"o\")\n",
    "        plt.title(\"Cubic regression\")\n",
    "        plt.grid()\n",
    "\n",
    "    return fitResult"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "c66ae3c4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parameters: [0.17054737581685814, 0.7278125943022667, 2.1137956163462075]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEICAYAAAC3Y/QeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAArBUlEQVR4nO3deXxU9b3/8dcnC2EJOxjCjggqIKJEtC5XcLlYlwu2lWJbSqsU26u3Wm2t2l7bWulte2v1tlVb+nNrrUXqSq1alRpcKrIoZRHByBp2AgQC2efz+2MOGuKEBDKTk5l5Px+PPObM92yf7wQ+8833nPP9mrsjIiLpJSPsAEREpOUp+YuIpCElfxGRNKTkLyKShpT8RUTSkJK/iEgaUvKXpGZm68zsggbWnWNmq1o6pjCY2QozGxt2HJI8lPwldGb2BTNbZGZlZrbFzF4ws7Obe1x3f93dj49HjK2duw9398Kw45DkoeQvoTKzG4F7gJ8AeUB/4D5gQohhHcLMsuJ4LDMz/b+T0OkfoYTGzDoDdwDXuvtT7r7f3avd/a/u/p1gm4fN7M46+4w1s+J6hzrNzN4zs91m9pCZtY21rZn1M7OnzGyHmZWY2W8aiOuHZvaEmT1qZnuBr5hZZzN7IPjLZJOZ3WlmmcH2mWZ2l5ntNLO1ZnadmfnBLw0zKzSzGWb2JnAAONbMTjCzl81sl5mtMrNJdc5/cVCffcG5vh2U9zCz58xsT7Df6we/SOp2f5lZjpndY2abg597zCyn7mdiZjeZ2fagPl9tzu9RkpOSv4TpU0Bb4OlmHueLwHhgMDAU+H79DYJE/RywHhgI9AFmHeaYE4AngC7An4BHgBrgOOAU4N+BacG2XwM+DYwCTgUmxjjeFGA60BHYAbwMPAYcA1wJ3Gdmw4NtHwCucfeOwAjgH0H5TUAx0JPoX0m3AbHGZ/kecEYQz8nAGA79THoBnYPP4GrgXjPrepjPQlKQkr+EqTuw091rmnmc37j7RnffBcwgmkzrGwP0Br4T/IVR4e5vHOaYb7n7M+4eAToRTe43BPtuB+4GJgfbTgL+z92L3X038NMYx3vY3VcEdb0IWOfuD7l7jbu/AzwJfC7YthoYZmad3H13sP5geT4wIPgL6XWPPTjXF4E73H27u+8AfkT0y4c6x7kjOMbzQBmQFtdG5GNK/hKmEqBHHPrUN9ZZXk80ydfXD1h/BF80dY85AMgGtgRdLnuA3xFttROcb2MD+zZ0vNMPHis43heJtsgBPgtcDKw3s3lm9qmg/H+BIuAlM1tjZrc0EHtvop/DQfU/k5J6n8MBILeBY0mKUvKXML0FVBC7m+Sg/UD7Ou97xdimX53l/sDmGNtsBPofwRdN3Rb1RqAS6OHuXYKfTu5+sJtmC9C3gXgaOt68Osfq4u657v4NAHdf6O4TiH65PAPMDsr3uftN7n4scBlwo5mdH+Ncm4l+wRzU0GciaUzJX0Lj7qXA7UT7nCeaWXszyzazT5vZz4PNlgAXm1k3M+sF3BDjUNeaWV8z60a0H/zxGNssIJqkf2pmHcysrZmd1cQ4twAvAXeZWSczyzCzwWZ2brDJbOB6M+tjZl2A7zZyyOeAoWY2JahvtpmdZmYnmlkbM/uimXV292pgL1ALYGaXmtlxZmZ1ymtjHP/PwPfNrKeZ9SD6GT/alLpK+lDyl1C5+y+BG4lekNxBtFV8HdEWL8AfgX8B64gm4FiJ/bFg3Zrg5876G7h7LdHW8nHABqIXTj9/BKF+GWgDvAfsJnoxOD9Y9/vg/EuBd4HniV4cjpWYcfd9RC8YTybaIt8K/AzICTaZAqwL7jT6OvCloHwI8ArRPvq3gPsauLf/TmBREM8y4B1ifCaS3kyTuYjEl5l9Gvituw9odGORkKjlL9JMZtYuuDc/y8z6AD+g+beviiSUWv4izWRm7YF5wAlAOfA34Hp33xtqYCKHoeQvIpKG1O0jIpKG4jZgVaL16NHDBw4cGHYYR2z//v106NAh7DBanOqdXlTv1mvx4sU73b1n/fKkSf4DBw5k0aJFYYdxxAoLCxk7dmzYYbQ41Tu9qN6tl5mtj1Xe7G6f4GGZBWb2L4tOKPGjoLxbMGrhB8Fr1zr73GpmRcFohuObG4OIiByZePT5VwLnufvJREcRvMjMzgBuAea6+xBgbvAeMxtG9OGW4UQHuLrv4NC4IiLSMpqd/D2qLHibHfw40SFxHwnKH+Hj8VsmALPcvdLd1xIdqGpMc+MQEZGmi0uff9ByX0z00fl73f1tM8sLxkTB3beY2cEREPsA8+vsXhyUxTrudKJjoJOXl0dhYWE8wm1RZWVlSRl3c6ne6UX1Tj5xSf7BuCmjgkGtnjazEYfZ3GIdooHjzgRmAhQUFHhrv7ASSzJcEEoE1Tu9qN7JJ673+bv7HqCQaF/+NjPLBwhetwebFXPokLd90XCzIiItKh53+/QMWvyYWTvgAuB9YA4wNdhsKvBssDwHmBzMMzqI6EiFC5obh4iINF08un3ygUeCfv8MYLa7P2dmbwGzzexqokPoXgHg7ivMbDbRoXFriE7eHXPoWxGRdLa0eA+vrd7B1DMH0rFtdlyP3ezk7+5LiU5oXb+8BIg1yxDuPoPoXKsiItKAh95cx8vvbeOrZw2K+7E1to+ISCtUUlbJ35Zu4bOn9qFDTvwHY1DyFxFphf6yuJiq2ghfOiMxcwIp+YuItDK1EedPb6/n9EHdGJLXMSHnUPIXEWllXlu9g427ypnyqcTNBKrkLyLSyjw6fz09cnP492G9EnYOJX8RkVZk464D/GPVdq4c0482WYlL0Ur+IiKtyGMLNmDAlWP6J/Q8Sv4iIq1EZU0tsxdu5PwT8+jdpV1Cz6XkLyLSSry4fCsl+6uYkqDbO+tS8hcRaSX++NZ6BnRvz9nH9Uj4uZT8RURagZVb9rJo/W6+dPoAMjJijXwfX0r+IiKtwKPz15OTlcHnRvdtkfMp+YuIhGxfRTVPv7uJS0f2pmuHNi1yTiV/EZGQPfPuJg5U1Sb0id76lPxFREIUiTgP/3MdJ/XpzMl9O7fYeZX8RURCNO+DHXy4Yz9XnT0Qs8Rf6D1IyV9EJEQPvrGWYzrmcMlJvVv0vPGYw7efmb1qZivNbIWZXR+U/9DMNpnZkuDn4jr73GpmRWa2yszGNzcGEZFk9P7Wvbz+wU6mnjkwoeP4xBKP6WFqgJvc/R0z6wgsNrOXg3V3u/sv6m5sZsOAycBwoDfwipkN1Ty+IpJuHnxjLW2zM/hCgsfxiaXZXzXuvsXd3wmW9wErgT6H2WUCMMvdK919LVAEjGluHCIiyWRnWSXPLNnMZ0/t22K3d9YV14khzWwg0cnc3wbOAq4zsy8Di4j+dbCb6BfD/Dq7FdPAl4WZTQemA+Tl5VFYWBjPcFtEWVlZUsbdXKp3elG9j9wzRVVU1UQYkb0jlM8ubsnfzHKBJ4Eb3H2vmd0P/Bjw4PUu4Cog1uVsj3VMd58JzAQoKCjwsWPHxivcFlNYWEgyxt1cqnd6Ub2PTEV1LTe9/g/GHd+TKy8Np+MjLlcYzCybaOL/k7s/BeDu29y91t0jwO/5uGunGOhXZ/e+wOZ4xCEikgzmLNlMyf4qpp1zbGgxxONuHwMeAFa6+y/rlOfX2exyYHmwPAeYbGY5ZjYIGAIsaG4cIiLJwN158M21nNCrI2cO7h5aHPHo9jkLmAIsM7MlQdltwJVmNopol8464BoAd19hZrOB94jeKXSt7vQRkXTxZlEJ72/dx88/N7JFH+qqr9nJ393fIHY//vOH2WcGMKO55xYRSTYPvLGGHrk5TBjVsg911acnfEVEWkjR9n28umoHU84YQE5WZqixKPmLiLSQB99cR5usDL54Rss/1FWfkr+ISAvYWVbJk4uLuXxUH3rk5oQdjpK/iEhLePCNtVTVRrjm3PBu76xLyV9EJMFKy6v541vrufikfI7tmRt2OICSv4hIwj06fz37Kmv4xrmDww7lI0r+IiIJVF5Vy4NvrGXs8T0Z0aflZupqjJK/iEgCPb5wAyX7q7h23HFhh3IIJX8RkQSpqokw87U1nDawK6cN7BZ2OIdQ8hcRSZBnl2xic2kF/9nKWv2g5C8ikhC1Eef+eR8yLL8TY4f2DDucT1DyFxFJgJdWbGXNjv3857jBoQ7g1hAlfxGROHN37i0sYlCPDnx6RH7jO4RAyV9EJM5e+2Anyzft5RvnDiYzo/W1+kHJX0Qk7u57tYj8zm2ZeErM6clbBSV/EZE4entNCW+v3cW0c46lTVbrTbGtNzIRkSTj7tz10mqO6ZjDF08Pf9jmw4nHHL79zOxVM1tpZivM7PqgvJuZvWxmHwSvXevsc6uZFZnZKjMb39wYRERagzeKdrJg3S6uO+842maHO1lLY+LR8q8BbnL3E4EzgGvNbBhwCzDX3YcAc4P3BOsmA8OBi4D7zKx1f0oiIo042Orv3bktnz+tX9jhNKrZyd/dt7j7O8HyPmAl0AeYADwSbPYIMDFYngDMcvdKd18LFAFjmhuHiEiYXl21nSUb9/Bf5w8JfYrGpjB3j9/BzAYCrwEjgA3u3qXOut3u3tXMfgPMd/dHg/IHgBfc/YkYx5sOTAfIy8sbPWvWrLjF2lLKysrIzW0d43e3JNU7vaR7vd2dH75VwYFq53/OaUdWK7q9c9y4cYvdvaB+eVa8TmBmucCTwA3uvvcwT7TFWhHzG8jdZwIzAQoKCnzs2LFxiLRlFRYWkoxxN5fqnV7Svd4vLt/K+r2L+cUVJ3PB6L5hh9Ukcbnbx8yyiSb+P7n7U0HxNjPLD9bnA9uD8mKgbodYX2BzPOIQEWlpkYhz98urObZHByaO6h12OE0Wj7t9DHgAWOnuv6yzag4wNVieCjxbp3yymeWY2SBgCLCguXGIiIThuWVbWLVtH9dfMISszOS5ez4e3T5nAVOAZWa2JCi7DfgpMNvMrgY2AFcAuPsKM5sNvEf0TqFr3b02DnGIiLSo2ohzzyurGZqXy2Ujk6fVD3FI/u7+BrH78QHOb2CfGcCM5p5bRCRM87fUsGZHFb/90qlktKKLvE2RPH+jiIi0ItW1EZ4pqmZ4706MH94r7HCOmJK/iMhReHzhRnaUOzdeOLRVjtffGCV/EZEjVFZZE+3r75rBeSccE3Y4R0XJX0TkCP1u3ofsLKti8vFtkrLVD0r+IiJHZGtpBb9/fQ2XndybY7u0/mEcGqLkLyJyBO56aRWRCNw8/viwQ2kWJX8RkSZ6b/NenninmKlnDqBft/Zhh9MsSv4iIk30Py+spFPbbK4bNyTsUJpNyV9EpAnmrd7B6x/s5JvnD6Fz++yww2k2JX8RkUbURpyf/G0lA7q3Z8oZA8IOJy6U/EVEGvHE4o2s2raPm8ef0KonZT8SqVELEZEEOVBVw10vreaU/l24+KTkG8ahIUr+IiKH8bt5a9i+r5LvX3Ji0j7QFYuSv4hIAzaUHOD+eR9y6ch8Rg/oFnY4caXkLyLSgB/9dQXZGcb3LxkWdihxp+QvIhLDK+9tY+7727n+giH06tw27HDiTslfRKSeiupafvTcCoYck8tXzxoUdjgJEa8J3B80s+1mtrxO2Q/NbJOZLQl+Lq6z7lYzKzKzVWY2Ph4xiIjEy/2FH7JxVzk/mjCc7CSal/dIxKtWDwMXxSi/291HBT/PA5jZMGAyMDzY5z4zS96h8UQkpRy8yHvZyb05c3CPsMNJmLgkf3d/DdjVxM0nALPcvdLd1wJFwJh4xCEi0lwHL/J+7+ITww4loRL998x1ZrY06BbqGpT1ATbW2aY4KBMRCdXBi7w3XDA0JS/y1mXuHp8DmQ0EnnP3EcH7PGAn4MCPgXx3v8rM7gXecvdHg+0eAJ539ydjHHM6MB0gLy9v9KxZs+ISa0sqKysjNzc37DBanOqdXlKh3lW1zvfeKCc7E+44sx1ZGY0/0JUM9R43btxidy+oX56VqBO6+7aDy2b2e+C54G0x0K/Opn2BzQ0cYyYwE6CgoMDHjh2bkFgTqbCwkGSMu7lU7/SSCvW+++XV7Cj/gMe+dnqT+/qTud4J6/Yxs/w6by8HDt4JNAeYbGY5ZjYIGAIsSFQcIiKNWb1tH/cVFjFhVGpf5K0rLi1/M/szMBboYWbFwA+AsWY2imi3zzrgGgB3X2Fms4H3gBrgWnevjUccIiJHqjbi3PzEUjq2zeb2S1PvSd6GxCX5u/uVMYofOMz2M4AZ8Ti3iEhzPPTmWpZs3MOvrjyF7rk5YYfTYlLz6QURkSZYt3M/v3hpFRecmMdlI/Mb3yGFKPmLSFqKRJxbnlpKdmYGMy4fkVLDNTeFkr+IpKXHFmxg/ppdfP+SE8nrlNr39Mei5C8iaWfznnJ++sL7nH1cDyYV9Gt8hxSk5C8iacXdue3pZdRGnP/5zElp191zkJK/iKSVp9/dROGqHdx80fH069Y+7HBCo+QvImljS2k5P/rre4we0JWpnxoYdjihUvIXkbQQiTg3Pv4vqmsj/OKKk8lowtg9qSxhY/uIiLQmM19fw1trSvj5Z0cyqEeHsMMJnVr+IpLylhWXctdLq7j4pF5cUdA37HBaBSV/EUlpB6pquH7Wu3TvkMNPLk/fu3vqU7ePiKS0Hz+3krUl+/nTtNPp0r5N2OG0Gmr5i0jKenH5Vv68YAPX/NvgtBmquamU/EUkJW0treCWp5ZyUp/O3Hjh0LDDaXWU/EUk5UQizk1/WUJldYR7Jo+iTZZSXX36REQk5fz6H0W8WVTC7ZcNY3DP1j3HbliU/EUkpby6ajv3zF3NZ07tw+TT0nPQtqaIS/I3swfNbLuZLa9T1s3MXjazD4LXrnXW3WpmRWa2yszGxyMGEZGNuw5ww6wlnNCrEzMm6rbOw4lXy/9h4KJ6ZbcAc919CDA3eI+ZDQMmA8ODfe4zs8w4xSEi6WLpbLh7BPywC9w9gqp3Z/H1Rxfj7vz2S6fSro3SyuHEJfm7+2vArnrFE4BHguVHgIl1yme5e6W7rwWKgDHxiENE0sTS2fDXb0LpRsCjr3O+yeCtz3PP5FEM6K7hGxqTyD7/PHffAhC8HhOU9wE21tmuOCgTEWmauXdAdfkhRW28kjs7PsV5J+SFFFRyCeMJ31idcB5zQ7PpwHSAvLw8CgsLExhWYpSVlSVl3M2leqeXlq73uaXFMRNJx8ptLRpHMv++E5n8t5lZvrtvMbN8YHtQXgzUvQTfF9gc6wDuPhOYCVBQUOBjx45NYLiJUVhYSDLG3Vyqd3pp8Xq/2zfo8jmUde7bonEk8+87kd0+c4CpwfJU4Nk65ZPNLMfMBgFDgAUJjENEUs35t+NZ7Q4ty24H598eTjxJKF63ev4ZeAs43syKzexq4KfAhWb2AXBh8B53XwHMBt4DXgSudffaeMQhIukhMuIKHulxI8WRHjgGnfvBZb+CkZPCDi1pxKXbx92vbGDV+Q1sPwOYEY9zi0j6uevlVdy7bjiVn36Ja84dHHY4SUlP+IpIUpm9cCP3vvohV47pz/R/OzbscJKWkr+IJI03PtjJbU8v45whPbhjwnA9wdsMSv4ikhRWb9vHNx5dzOCeudz7xVPJzlT6ag59eiLS6u3YV8lXH1pI2zaZPPjV0+jUNjvskJKekr+ItGql5dVMfXABu/ZX8eDU0+jTpV3jO0mjlPxFpNXaX1nDVx5awAfb93H/l07lpL6dww4pZWgCdxFplSqqa5n2yCKWFpdy7xdOZezxxzS+kzSZWv4i0upU1UT4xqOLmb+2hF9cMZKLRvQKO6SUo+QvIq1KTW2Ebz2+hFdX7WDGxJO4/JS+YYeUkpT8RaTViESc7z65jL8t28L3LzmRL5zeP+yQUpaSv4i0CpGI89/PLufJd4r51gVDmXaOnt5NJF3wFZHQ1dRGuPnJpTz1zia+fu5gvnn+cWGHlPKU/EUkVJU1tVz/5yW8uGIrN104lOvOO07DNrQAJX8RCc2Bqhqu+eNiXv9gJ7dfOoyrzh4UdkhpQ8lfREKxt6Kaqx5ayDsbdvPzz45k0mn9Gt9J4kbJX0RaXElZJVMfWsCqrfv49ZWncsnI/LBDSjtK/iLSojbuOsBXHlpA8e5yZn65gHF6cjcUSv4i0mIWrdvF9D8upqY2wiNXjeGMY7uHHVLaSnjyN7N1wD6gFqhx9wIz6wY8DgwE1gGT3H13omMRkfA8/W4x331iGX26tuOBqQUc2zM37JDSWks95DXO3Ue5e0Hw/hZgrrsPAeYG70UkBUUizl0vreJbj/+LUwd04en/PFOJvxUIq9tnAjA2WH4EKAS+G1IsIpIg5VW1fPsv/+Jvy7YwqaAvd048iTZZGligNTB3T+wJzNYCuwEHfufuM81sj7t3qbPNbnfvGmPf6cB0gLy8vNGzZs1KaKyJUFZWRm5u+rVyVO/0EqveJeURfrOkknWlESYd34aLBmal3MNbyfD7Hjdu3OI6vS4faYmW/1nuvtnMjgFeNrP3m7qju88EZgIUFBT42LFjExRi4hQWFpKMcTeX6p1e6tf71fe3c+fsJVTVGL+dMprxw1NzSOZk/n0nPPm7++bgdbuZPQ2MAbaZWb67bzGzfGB7ouMQkcSrro3wi5dW8bt5azgxvxP3fuEU9e+3UgntfDOzDmbW8eAy8O/AcmAOMDXYbCrwbCLjEJHE27ynnMkz5/O7eWv4wun9dWG3lUt0yz8PeDro58sCHnP3F81sITDbzK4GNgBXJDgOEUmgJdtruOFXr1NdE+FXV57Cf5zcO+yQpBEJTf7uvgY4OUZ5CXB+Is8tIom3v7KGn7/4Po+8U6luniSjJ3xF5Ki8WbST7z65lE17yrlwQBa/nnYmbbMzww5LmkjJX0SOyL6Kan7y/Pv8ecEGBvXowOxrPsX+dUuV+JOMkr+INNm81Tu49cmlbN1bwdfOGcSNFx5PuzaZFK4LOzI5Ukr+ItKoraUV/OzF93n63U0M7tmBJ75xJqf2/8RzmZJElPxFpEEV1bU88MZa7n21iJpa57pxx3Hdecd93MWzdDbMvYNzS4vh3b5w/u0wclK4QUuTKPmLyCe4Oy8u38qM51dSvLuci4b34raLT6R/9/Yfb7R0Nvz1m1BdjgGUboy+B30BJAElf5F0FLTYKS2Gzoe22FdsLuXHz73H/DW7OKFXRx772umcObjHJ48x9w6oLj+0rLo8Wq7k3+op+YukmzotduCjFnvx7nLu3DCCF1dspWv7bO6cOILJp/UjK7OBgQBKi4+sXFoVJX+RdHOYFvubdh/fPH8IV581iM7tsw9/nM59o18cscql1VPyF0k3DbTM+2SU8MbN5zWe9A86//ZD/4IAyG4XLZdWT7MqiKSRmtoI5e3zY66zzn2bnvgh2q9/2a+gcz8cg879ou/V358U1PIXSQO791cxa+FGHp2/ntF7J/KzNg/QjsqPNzjaFvvISTByEvOSeFz7dKXkL5Ki3J3F63cze9FGnl2ymcqaCGcO7s4ll11PTs0p8I/Yd/tIelDyF0kWDd2eWa9825ib+dP+03l6ySY27iqnXXYmnx3dl6mfGsjxvToGB5sEJyvZpzMlf5Fk0MDtmWyYD/967JDyji/dxIaaaQw8diLfumAo44f3okOO/qvLofQvQiSRDvMw1RFp4PbMyKKHyCBySHF7q+KubnPIvPp/mhG4pDolf5F4iDXGDcRurcORfwE0cHumeYTo2AqHyty36ciOL2kntFs9zewiM1tlZkVmdktYcYg028EumdKNGP5xkn/huw0Pf9AE2/dWMOdfm/ne08vYZjGGVwDIaGAMfT1oJY0IpeVvZpnAvcCFQDGw0MzmuPt7YcQj0iwNPTFbv+ygGK34krJKlm/ey/JNpazYXMqyTaVs3BXdPzcni77dr2bannvIjlR8vFN2O+zkLxza5x+U60EraUxY3T5jgKJgjl/MbBYwAVDyl+RzhGPZVHbozTMLN7Bmx34+3FHGis172VL6cVLv3609J/XpzJfPGMjpx3ZjWH4nsjLHw9LBsa8f9D8jPtcVJK2Yu7f8Sc0+B1zk7tOC91OA0939unrbTQemA+Tl5Y2eNWtWi8faXGVlZeTmpt+E1ulU7zPemkbbyh2fKC+zXLK9ihyqPio74G24pXoacyJnk2VwTAejf8cMBnTKZGCnDPp3yqBDdoxO/FYunX7fdSVDvceNG7fY3Qvql4fV8o/1r/sT30LuPhOYCVBQUODJ+ARhYZo++ZhM9XZ3qmudippaKqprqayOUFlTS0V1hP2VNZQFP/sqguWKGnYfqGL3gSpKyqr4p03hW/yGdvWS/G1VXyYrM4Obsx7nGN9JafYxLBl6PZePvIJv98ilT9d2ZGYkX6KPJZl+3/GUzPUOK/kXA/3qvO8LbA4pFomT2ohTvvjP5My7k3PLNlOxMJ8PRtzImvyLOVBVy/7KGsqraqmsiVBVG6GqJkJlTYThO//Oxdtn0rVmO7uyevJUl6v5Z4fzqIk4tREn4k7Eo0nanUPfw0dlB18PxlK3rNadSARqIhFqap3q2gg1EY8uRyIcyR/AmRlGp7ZZdOvQhu4dctjQ51LmVOVy0faZdKrcRkX7fEpO/y4/GH0l3Tq0wWwGAF2BcfH/2EWOSljJfyEwxMwGAZuAycAXQopFDiMScXaUVbK1tIJteyvYtq+S7Xujy9v3VbJrfxV7DlSz50AVY6sK+Wn2/yPboi3gtvs3M3j+bfy+eg1zImd/dEwzaJOZQU5WBpfZm3zWf/vRODPda7YzZecvKa2o5s3255FpRkaGkWFglkFGBmRYtLWcYYZZ8ApY8N6IJujoftF9o69GdqaRmWFkZ2aQlWFkZWaQnWnkZGXQNjuTnKwMcrIzP1rOzcmK/rTNomNOFh3bZtM2OwOz+i320cBNH7UE+yHSuoWS/N29xsyuA/4OZAIPuvuKMGKRaAt6694KVm3dx/qSA6wvOcCGXfuD1wNU1hz6EFGGQc+OORzTsS3dc9swuGcundtlc+OKp2hfWXXIttEHjp7ltq/9gPY5mbTLziQrwz5OnnffBKWVh+zTlkq+nfk43/7P/05ovUXSWWgPebn788DzYZ0/XVVU17Jyy17e37qP97fsZeXWfazauo/S8uqPtmmbncGAbh0Y2KMD5w7tSf/u7enduR15ndqS1ymH7rk5sfuq39kW85zZZZvp1blt7IA0G5RIKPSEb4rbWlrB4vW7Wbx+N+9s2M2KzaVU10Y7uDu0yeT4Xh25ZGQ+J/TqyNC8jhzbowM9O+bE6NZogqOZ2UmzQYmEQsk/xewsq2TtPx5m8NK76FKznZpId/5eM4m/Z/wbJ/frwrRzjmVUvy4My+9Eny7tyIjn3SZHM7OTZoMSCYWSf5Krro3wzvrdvPbBDuat3sGxW57np9n/j/bBRde+GTu5p91DRC4bRdaoTyc2mIMPFs29Ay8txprywFGdffSQkkjLUfJPQtW1Ed4o2snflm7hpRVb2VtRQ2aGMbp/V+7s+MmLrhm15WS8+mMY9fnEB3c0MzsF+4hIy1HyTxI1tRH++WEJf1u6hRdXbKW0vJqOOVlcODyPfx+Wx5nH9aBT22z4YeyLrrqAKiJ1Kfm3cmt2lDFr4UaeXFxMyf4qOrTJ5MJheVw6sjfnDO1BTla9UR11AVVEmkDJP1Fije/exK6NqpoIf1+xlcfe3sBba0rIzDAuPDGPy0/tw7lDe9I2u4FhfEEXUEWkSZT8E6HOlHsGTZ7EY9Oecv7w1jqeWBRt5fft2o7vjD+eK0b35ZhODdwnX58uoIpIEyj5J0JD47vPvSNmEi7aXsZv533IM+9uwoELTjyGK8f059+G9Dy6WzF1AVVEGqHknwhNfGr1Xxv3cH/hh/z9va3kZGXwpTMGMO2cQfTt2r4FghSRdKbknwiNXHRduG4X97yymjeLSujUNovrxh3HV84cSPfcnBYOVETSlZJ/IjRw0XXraTfz339YxMvvbaNnxxxuu/gErhzTn45ts8OLVUTSkpJ/ItR70jXSsQ9PdLmK257vQbvsEr4z/niuOmsQ7doc5q4dEZEEUvJPlJGT2H/8Z7jtD//g5Y0RqkoiTDljAP913nHq3hGR0Cn5J8iLy7fygznL2ba3mktOyuc7449nYI8OYYclIgIo+cfd5j3l/GDOCl5+bxsn9OrI1040pl1+athhiYgcQsk/TmojziP/XMddL62i1p1bP30CV509iDdffy3s0EREPiFhyd/Mfgh8DdgRFN0WzN6Fmd0KXA3UAt90978nKo6WsGJzKbc+tYylxaWcO7Qnd04cQb9uuldfRFqvRLf873b3X9QtMLNhRCdsHw70Bl4xs6HuXpvgWOKuNuL8dt6H3P3yarq0z+ZXV57CZSPzj24WLBGRFhRGt88EYJa7VwJrzawIGAO8FUIsR23znnK+9fgS3l67i0tG5jNj4gi6tG8TdlgiIk1i7p6YA0e7fb4C7AUWATe5+24z+w0w390fDbZ7AHjB3Z+IcYzpwHSAvLy80bNmzUpIrEdq4dYaHl5RSU0EvnRiG87uk9Vga7+srIzc3NwWjjB8qnd6Ub1br3Hjxi1294L65c1q+ZvZK0CvGKu+B9wP/Bjw4PUu4CogVpaM+Q3k7jOBmQAFBQXe5JmhEmR/ZQ0/+usKZi8p5uS+nfm/yac0evtm4ZHMaJVCVO/0ononn2Ylf3e/oCnbmdnvgeeCt8VAvzqr+wKbmxNHS3h/616+8eg7rCvZz7XjBnPDBUPJzswIOywRkaOSsOxlZvl13l4OLA+W5wCTzSzHzAYBQ4AFiYojHl5YtoXP3PdP9lfW8Ni0M/jO+BOU+EUkqSXygu/PzWwU0S6ddcA1AO6+wsxmA+8BNcC1rfVOn0jEufuV1fz6H0WM6teF300ZTV5TJ1UREWnFEpb83X3KYdbNAGYk6tzxsK+imm89voRXVm7nitF9+fHEEYefPlFEJInoCd8YPtxRxvQ/LGJdyQF+9B/D+fKnBujefRFJKUr+dS2dTcWLP2DQgS38ke6UjfseQ8+8OOyoRETiTsn/oKWzqXn2v2hbWwFAb3bCgtsgL1fz4YpIytEtK4Gy528nK0j8Hzk46bqISIpJ++Tv7vzsxfdpX74l9gYNTcYuIpLE0jr519RGuPmJpdxf+CGlbfJibxRMui4ikkrSNvmXV9VyzR8X85fFxXzz/CF0uexOyG536EbZ7aKTsYuIpJi0vOB7oKqGrz60kAXrdvHjiSOYcsYAYGh05dw7ol09nftGE78u9opICkq75F9eVcvVDy9i4bpd3PP5UUwY1efjlSMnKdmLSFpIq26fiupapv1hIW+vLeGXk+olfhGRNJI2Lf+K6lq+9odF/PPDEu664mQmnqLELyLpKy1a/gcT/xtFO/nfz53MZ07VHTwikt5SPvlX1tTy9UcX8/oHO/nZZ0byudFK/CIiKZ38q2sjfOPRdyhctYOffuYkJp3Wr/GdRETSQEr3+WdlGMf26MAFl5/E5DH9ww5HRKTVSOnkb2Z8/9JhYYchItLqpHS3j4iIxNas5G9mV5jZCjOLmFlBvXW3mlmRma0ys/F1ykeb2bJg3a9Ms6SIiLS45rb8lwOfAV6rW2hmw4DJwHDgIuA+Mzs4B+L9wHSiE7cPCdaLiEgLalbyd/eV7r4qxqoJwCx3r3T3tUARMMbM8oFO7v6WuzvwB2Bic2IQEZEjl6gLvn2A+XXeFwdl1cFy/fKYzGw60b8SyMvLo7CwMO6BJlpZWVlSxt1cqnd6Ub2TT6PJ38xeAXrFWPU9d3+2od1ilPlhymNy95nATICCggIfO3bs4YNthQoLC0nGuJtL9U4vqnfyaTT5u/sFR3HcYqDuE1V9gc1Bed8Y5SIi0oISdavnHGCymeWY2SCiF3YXuPsWYJ+ZnRHc5fNloKG/HkREJEEset31KHc2uxz4NdAT2AMscffxwbrvAVcBNcAN7v5CUF4APAy0A14A/subEISZ7QDWH3Ww4ekB7Aw7iBCo3ulF9W69Brh7z/qFzUr+0jgzW+TuBY1vmVpU7/SieicfPeErIpKGlPxFRNKQkn/izQw7gJCo3ulF9U4y6vMXEUlDavmLiKQhJX8RkTSk5J8gZva/Zva+mS01s6fNrEuddTGHu04FRzPMd6ows4uCuhWZ2S1hx5MoZvagmW03s+V1yrqZ2ctm9kHw2jXMGBPBzPqZ2atmtjL4N359UJ6UdVfyT5yXgRHuPhJYDdwKjQ53nQqOZpjvpBfU5V7g08Aw4MqgzqnoYT45FPstwFx3HwLMDd6nmhrgJnc/ETgDuDb4HSdl3ZX8E8TdX3L3muDtfD4e0yjmcNdhxJgIRzrMd8tGl1BjgCJ3X+PuVcAsonVOOe7+GrCrXvEE4JFg+RFScKh2d9/i7u8Ey/uAlURHJU7Kuiv5t4yriA5lAdF/LBvrrDvssNYpJNXrner1a0xeMHYXwesxIceTUGY2EDgFeJskrXtKT+CeaE0Z7joY46gG+NPB3WJsn1T328Z5mO9Uker1k4CZ5QJPEh2zbG+yzkSr5N8MjQ13bWZTgUuB8+sMXtfQcNdJI87DfKeKVK9fY7aZWb67bwlm7NsedkCJYGbZRBP/n9z9qaA4Keuubp8EMbOLgO8C/+HuB+qsijncdRgxtrBUr/dCYIiZDTKzNkQvbs8JOaaWNAeYGixPJQWHag+GoX8AWOnuv6yzKinrrid8E8TMioAcoCQomu/uXw/WxRzuOhUczTDfqcLMLgbuATKBB919RrgRJYaZ/RkYS3Q4423AD4BngNlAf2ADcIW7178onNTM7GzgdWAZEAmKbyPa7590dVfyFxFJQ+r2ERFJQ0r+IiJpSMlfRCQNKfmLiKQhJX8RkTSk5C8ikoaU/EVE0tD/B1xe5Zm+nffeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ret = cubicFit(\n",
    "    [\n",
    "        [ -23, -70 ],\n",
    "        [ -20, -30 ],\n",
    "        [ -15, -20 ],\n",
    "        [ -10, -10 ],\n",
    "        [ -8, -10 ],\n",
    "        [ 0, 1 ],\n",
    "        [ 1, 3 ],\n",
    "        [ 4, 7 ],\n",
    "        [ 5, 8 ],\n",
    "        [ 7, 20 ],\n",
    "        [ 10, 25 ],\n",
    "        [ 15, 100 ]\n",
    "    ],\n",
    "    trace = False,\n",
    "    plot = True,\n",
    "    plotRange = [ - 25, 25 ]\n",
    ")\n",
    "\n",
    "print(\"Parameters: {}\".format(res))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}