{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9ff541fc-9f04-4d11-a6d9-ce34de342377",
   "metadata": {},
   "source": [
    "# Introduction to `lamatrix`: Getting started."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ecc5fe4-c734-4fc3-a375-01a86e3c5872",
   "metadata": {},
   "source": [
    "`lamatrix` is a tool to help you fit linear models to data. This tool is designed to make it easier to build models and fit them. \n",
    "\n",
    "These are the basic principles behind `lamatrix`. \n",
    "\n",
    "Firstly, let's create some data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1998681d-dc6b-41df-9072-cf6286e7a1f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from astropy.io import fits\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "85427fad-2cdc-4c4b-a3ba-e638ff542060",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(-1, 1, 0.01)\n",
    "y = x * 2.3459 + 1.435 + np.random.normal(0, 0.3, size=len(x))\n",
    "ye = np.ones(len(x)) * 0.3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a3653ae6-79c2-475b-a192-9478fcd2373d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjYAAAHHCAYAAACskBIUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAwg0lEQVR4nO3dCXgT1drA8VNAyl5BZC8ogiKyCQiCCwgoIuKK+5VFBUFEURSpXsENi8t1uYqAXhHvvSIKAi4XEVEQXEAEkU1QEKTKJqItixYt8z3vuXfyJWnaJmmSmTnz/z1PKE0nkzMzSebNe95zJs2yLEsBAAAYoIzTDQAAAEgUAhsAAGAMAhsAAGAMAhsAAGAMAhsAAGAMAhsAAGAMAhsAAGAMAhsAAGAMAhsAAGAMAhvAYwYMGKCOOeYY3z03ijd16lSVlpamtm7d6nRTAEcR2AAlnCjsW7ly5VT9+vX1yf3HH390unnGOHz4sPrnP/+pzj77bFWzZk11xBFHqFq1aqlzzjlHPf/88yo/P1+Z4uDBg+q+++5TixYtcqwN8vzBr+tKlSqphg0bqj59+qiXXnqpVPt77ty5ev2Ak8o5+uyABzzwwAPq2GOPVb///rtaunSpDng+/vhjtXbtWlWhQgWnm+dpv/32m7r44ovVe++9pzp37qzuuOMOVbt2bbV371710UcfqZtuukktW7ZMvfjii8qUwOb+++/X/+/ataujbZk4caKqUqWKDmQkUJdjcN1116mnnnpKvfPOOyozMzOuwGbChAkEN3AUgQ1Qgl69eqn27dvr/99www06q/DII4+ot956S11++eVON8/TbrvtNn1ClZPprbfeGvK3kSNHqm+//Va9//77yq3+/PNPnXEqX7688pq+ffvq17JtzJgx6pVXXlH9+vVTl112mQ7iAS+iKwqI0RlnnKF/bt68OXDfoUOH9ImhXbt2KiMjQ1WuXFkvt3DhwpDHSv2DpP8ff/xx3c1y3HHHqfT0dHXKKaeo5cuXF3quOXPmqBYtWujMkPycPXt2xDYdOHBABwLyLVvWd8IJJ+jnsCwrZDl57ptvvlnNmDFDNW/eXFWsWFF16tRJrVmzRv998uTJqkmTJvr5JKNQXL2GrFvqbS688MJCf5PsluyHG2+8scjH5+TkqH/84x/q3HPPLRTU2Jo2baqzNsEkkJBA6KSTTtLtlAyPPM8vv/wSspy07fzzz9fZtQ4dOuhlGzdurLu9wv36669qxIgRgf0n+0CCV3muSMdOnt8+duvXr4/q+Mvjjz76aP1/ydrYXUHB2Y0NGzbogKNGjRq6vRJQSwAdbt26dapbt276+DVo0EA99NBDIW2N1zXXXKODd8mSBQeUS5Ys0cGOdFnJNst+kqBUMm426aKVbI0I7uqyyX6TrNxRRx2l2y37aubMmaVuMxCOjA0QI/tkX7169cB9eXl5+iR91VVXqUGDBql9+/bp7pOePXuqzz//XLVp0yZkHdOmTdPLyAlZPvwfffRRdckll6jvvvtO15iI+fPnq0svvVQHINnZ2ernn39WAwcO1Cey8ADjggsu0CfR66+/Xj+XZEHuvPNO3cXw5JNPhiwvJyk5WQ4bNkz/LuuWAGDUqFHqueee04GEBAnSJuma+PDDDyPuB2n3X/7yF72cdB3Jydj29ttv630ify/Ku+++qwoKCopdJhLZZ9IdKPvilltuUVu2bFHPPvus+vLLL9Unn3wS2H9i06ZNOlCQ/dK/f381ZcoUfQKWk6oERnb3UJcuXfS+knXLyfvTTz9VWVlZaseOHTqICSZ1KBK4DR48WJ/kZbujOf4S1Ej3z9ChQ3X3mxxv0apVq0Cwctppp+k6rtGjR+vg6PXXX1cXXXSReuONN/RjxM6dO9VZZ52ls0X2chIkS7CQCNdee61en7z+pO5JSCAs+0naLoGJbNMzzzyjfvjhB/03+7hs375dB0T/+te/Cq336aef1q9TCZ4kEJw+fboOlqTbq3fv3glpO6BZACJ66aWXJN1hLViwwPrpp5+snJwca+bMmdbRRx9tpaen699tf/75p5Wfnx/y+F9++cWqXbu2dd111wXu27Jli17nUUcdZe3duzdw/5tvvqnvf/vttwP3tWnTxqpbt67166+/Bu6bP3++Xq5Ro0aB++bMmaPve+ihh0Kev2/fvlZaWpq1adOmwH2ynLRd2mGbPHmyvr9OnTpWXl5e4P6srCx9f/Cy/fv3D3nujRs36mUmTpwY8twXXHCBdcwxx1iHDx8ucv/edttt+rGrVq0KuV/2o+xv+7Znz57A35YsWaIf88orr4Q8Zt68eYXul3bKfYsXLw7ct3v3br39I0eODNz34IMPWpUrV7a++eabkHWOHj3aKlu2rLVt27aQY1etWjW9nmDRHn/ZHlnH2LFjC+2P7t27Wy1btrR+//33wH2y/zp37mw1bdo0cN+IESP0OpYtWxayXRkZGYWOVyTy3LKctCUSabf8/eKLLw7cd/DgwULLZWdn69fX999/H7hv2LBh+rGRhK/j0KFDVosWLaxu3boV214gVnRFASXo0aOH/rYt6Xf59i/fkCXjEZw5KVu2bKDOQroEJIMh36ilK2HlypWF1nnFFVeEZHzs7i3J2AjJFKxatUpnGaRrwybfoCWDE16wKc8v2Ytg0jUlsYxkRoJ17949ZMh2x44d9U/JDlWtWrXQ/XabIjn++OP1clKbYZNtl+eUb+bBXRHhJMshpIA1fHtkf9u3Ro0aBf4m2QHZH7If9uzZE7hJBkbWE971J/vK3rdC1ifddMHbJOuUZeR4BK9TjrtklBYvXhyyTtlPdpdSvMc/nCwvmTGp2ZJsj90GydJJ1kdqjeyReLJ/Tj31VN29Frxdsr8TwT4e0g5bcDZIuj2lbdKtJK8vyZRFI3gdkhHMzc3V+z2a/QPEgq4ooARSNyAncPkglq4MOdFJF0S4l19+Wf3tb3/TdRJ//PFH4H4ZURVOujuC2UGOXSfy/fffB2pMwsmJOfhkIMvWq1cvJCgRJ554Ysi6inpuO3AKHwVj3x9euxJOik2lbkeeR4IQCRRk+6VLozh2e/fv3x9yv3TH2PUdjz32mO5esskJXo6DDAePZPfu3cVuq72vg7dJ1rl69epCwUpR64x0PGM9/uGky0yChHvvvVffimqHdFPJfraDzvDXRSLYxyP49bRt2zZdQyQBffjrQY5HNKTLSWqBJGAPHlJeXPALxIPABiiBfDO2R0VJvcPpp5+urr76arVx48bAt9t///vfunZD/i61LXLilW/xUr8SXGRsk79FEl7smwxFPXe8bbryyit1Ialkbe6++269L2R/lXSibdasmf4pw+Zbt24duF8CDMmWCFlXMMmGyL4NzhAFi5RJKWmbZJ2SAZIao0gkqA0WqZYl1uMfzi78leHukqGJRAqaU0GOR/DzSdZK9o9kle666y593CRrKRkk2eZoipalrkvqa84880xdx1W3bl1dCyX1SlJvBiQSgQ0QA/tkJcWbUrAqxZtCRnfIiJtZs2aFfAMdO3ZsXM9jd79INiGcBFThyy5YsEB3HQR/y5bMQfC6kkWKZ6X4U4IN6Q6RDEt4wW1Rw+hlf9qPi4aMRJJtlaxOooplZZ2SpbCDqXhEe/yLyk7IY4Wc7EtqhxzPaF4X8bILf+0AS0bMffPNNzojJdk5W6Rh+EVtnxQ/yygvKWoPznZKYAMkGjU2QIxkGLRkceTkLaNjgjMDwZkAGTL72WefxfUc8o1WRtLIySQ41S8nExleHOy8887T36ol0Aomo6HkRCMBRLJJt5O0S7IVsi8ki1MS6SaSUVdSjxPe9qKyRVKDItv64IMPFlpWalpk2HasZJ1ynOSkG07WJ+stSbTHX2b5tdcbTDI88rqS4fZSXxXup59+CjneMseMjEwK/ntRWaxYSPZERnfJFABSi1XUtsn/ZZRTOMnkRNo+WYe8FuXYBY8ulOkMgEQjYwPEQU7gMlRVhh0PGTJED5eWb+syJFeyFzIEedKkSbp4NbyGJFqSGZJ1SdeXBADSFSBDbGWYcvA6ZSp8ySDdc889+mQh3ToyVPfNN9/Uc7NIRiLZpJ0yDFjqaySQKqoGJpwEh7Kvhg8frof/yrbIY6U4VTI/Mmw8uEtLhmXLsGLZN1KrIZddkCyHZDDkueVkKwXesR5LqR2RY2gPBZcCWclUSCZG9mnwRHaRRHv8Jcsk97322mu6i0uyXTI/kdyklkuOdcuWLfWQccni7Nq1SwdHMqz6q6++0uuQLjPJqtjz/9jDvSWTI7VC0ZJtk65UGXptzzws+1xeP/YQbiFdT/Iakm4yWa5atWo6AxOp9kr2nZBCdsn42EGu7JMnnnhCt1m6caVeSLZXurtiaTMQlZjHUQE+G+69fPnyQn8rKCiwjjvuOH2Tob4yLPfhhx/WQ4xlOPHJJ59svfPOO4WGR9tDhh977LFC64w0DPiNN96wTjzxRL3O5s2bW7NmzSq0TrFv3z49fLpevXrWEUccoYcHy3OED7eW55AhucGKatPChQv1/TNmzAjcF+m5bTfddJNeftq0aVYsZP/JvpZhvzVq1LDKlStn1axZUw9/njRpkvXbb78Veszzzz9vtWvXzqpYsaJVtWpVPUx61KhR1vbt2wPLSDt79+5d6LFdunTRt/D9J8PbmzRpYpUvX14/vwyzfvzxx/Ww5OL2k4j2+ItPP/1Ut12eJ/yYb9682erXr58eei/HsX79+tb555+vpxkItnr1ar0NFSpU0MvIkPUXX3wxpuHe9k3W0aBBA/08U6ZMCRlublu/fr3Vo0cPq0qVKnrfDBo0yPrqq6/04+XYBR/L4cOH6ykRZCh48ClG2ievS9k/zZo104+z2wIkUpr8E10IBABFkwJimZROJpCzu1wAINWosQFQalJrJCODZI4XghoATqLGBkDcpFZCRilJvYZMJlfUNZ8AIFUIbADETUZCyVBtKfj9+9//XuiaWACQatTYAAAAY1BjAwAAjEFgAwAAjOG7Ghu5rsn27dv11PNcfA0AAG+Qyhm5dIxc9LdMmaLzMr4LbCSoCb+KMQAA8IacnBzVoEGDIv/uu8DGvkig7BiZGhwAALhfXl6eTkwEX+w3Et8FNnb3kwQ1BDYAAHhLSWUkFA8DAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAICYHThwQKWlpemb/N8tCGwAAIAxCGwAAIAxCGwAAIAxCGwAAIAxCGwAAIAxCGwAAIAxCGwAAIAxCGwAAECp565xy7w2BDYAAMAYBDYAAMAYBDYAAMAYBDYAAMAYBDYAAMAYBDYAAMAYBDYAAMAYBDYAAMAY5ZxuAAAA8LYqVaootyBjAwAAjEFgAwAAInLLZRJiQWADAACMQWADAACMQWADAICHeLF7KJUIbAAAgDEIbAAAgDEIbAAAgDE8HdiMHz9e9zGOGDHC6aYAAAAX8Gxgs3z5cjV58mTVqlUrp5sCAABcwpOBzf79+9U111yjXnjhBVW9enWnmwMAAFzCk4HNsGHDVO/evVWPHj1KXDY/P1/l5eWF3AAAMNGBOIaCmzZ83HOBzfTp09XKlStVdnZ2VMvLchkZGYFbZmZm0tsIAIBfAgm38VRgk5OTo2699Vb1yiuvqAoVKkT1mKysLJWbmxu4yToAAPDDFbfT/hdA+SmIKqc8ZMWKFWr37t2qbdu2gfsKCgrU4sWL1bPPPqu7ncqWLRvymPT0dH0DAAClC5S8wFOBTffu3dWaNWtC7hs4cKBq1qyZuuuuuwoFNQAAIDkqV66sLMvSmSA3BT2eCmyqVq2qWrRoUWjHHnXUUYXuBwAA/uOpGhsAAEzjxmLiA/9rk5syMUZmbCJZtGiR000AAAAuQcYGAAAYg8AGAAAfDf/evXu367q+EonABgAAGMPzNTYAACA+VTxYHFwSMjYAAHh09JIbL1JtWZaeisUpZGwAAEBMwYuTgUtJyNgAAFAEsjTeQ2ADAACMQWADAACMQWADAACMQWADAACMQWADAACMQWADAACMQWADAIBHh4nb139ycih65cqV9aR8Tk/MZyOwAQDAZZg/J34ENgAAuJidlSku0HFD5sYtCGwAAECxl1CQm1cQ2AAAAGMQ2AAAAGMQ2AAA4CH2KCQvdQ+lEoENAAAwBoENAAAwBoENAMAxzNeCRCOwAQAAxijndAMAAPArmVgPiUXGBgAAGIPABgAAnw8F379/v2suYllaBDYAABg8D81+j10SobSosQEAIIyM0Ep0/UvwOpMdaFT+X/AU/NyJXJ+bkbEBACCJuPJ2ahHYAACMwbw48TNl3xHYAAB8wc0nbpOKd51GjQ0AAB6qn0mFyh6qqQlHxgYAAJfUzdhZJSbuix8ZGwAAEsiNQUnlsAyM/D8ZI7/cgIwNAAAwBhkbAIDnmZp98NJcN25BxgYAABiDwAYAABdeCgHxIbABAADGoMYGAIBiSIaGifO8g4wNAMCoWXzhbwQ2AADAGAQ2AADAGNTYAACMxhw3/uKpjM3EiRNVq1atVLVq1fStU6dO6t1333W6WQAAwCU8Fdg0aNBAjR8/Xq1YsUJ98cUXqlu3burCCy9U69atc7ppAAADeOkilMyrY0BXVJ8+fUJ+HzdunM7iLF26VJ100kmOtQsA4D52cLJr1y5Hnt8OOLwQJJnEU4FNsIKCAjVjxgwdXUuXVFHy8/P1zZaXl5eiFgIAkPz5c8Kv+eR3nuqKEmvWrNHRb3p6uhoyZIiaPXu2at68eZHLZ2dnq4yMjMAtMzMzpe0FgGQwdR4ZU7cLqeO5wOaEE05Qq1atUsuWLVNDhw5V/fv3V+vXry9y+aysLJWbmxu45eTkpLS9AID4yJdYAhwY3xVVvnx51aRJE/3/du3aqeXLl6unn35aTZ48OeLyktmRGwAguiHRXEIAXua5jE24w4cPh9TQAAC82c0jgZXXCm3JKrmPpzI20q3Uq1cv1bBhQ7Vv3z41bdo0tWjRIvXee+853TQAAOACngpsdu/erfr166d27NihC4Flsj4Jas4++2ynmwYAcACzCsPTgc2LL77odBMAAICLeb7GBgDg3fodZs9FohHYAAB8XbTsFDuokxuj0BKHwAYAABiDwAYAABiDwAYAPMQLV59O9NwudGXB2FFRAAB4FRerTA0CGwAAIuDSEt5EYAMAQJJ4KUtT2UNtLQ41NgBgqHhrU7x+/SPJtBQ1L47cz/BqsxHYAECCUOTqX0w06B50RQEA4BKmdAc5icAGAHygqOHhXEQy/uAj0r4jMHEeXVEAAMAYBDYAAMAYBDYAACMw4gmCwAYAABiDwAYA4CkMq0dxCGwAAL5CYGQ2AhsAgG8DAq/PsozCmMcGAJDQuVwAJxHYAAB8jwDNHAQ2AACjhc8GTLeT2QhsAACeQVYFJaF4GAB8XAgcTaAQy7KxPDdFu0gGMjYAAITNYMzsxd5FxgYAXC7WDEeiMyxerqvhEgv+Q2ADAACMQVcUAHiUqV0msl3CzxknxI+MDQAUgSJX59CVhHiRsQEARMyUMGldbPPjwB3I2ACAj7JLErw4mQWRQClSsESGBolCxgYAXCY4U2JnURAZWROEI2MDAACMQcYGAJAU1OjACWRsAAAJId1mdJ3BaQQ2AICkI+hBqhDYAEAcmOMGcCcCGwAAYAyKhwHAQ9w+vDm4fYnKZNnrpBgZ0SCwAQAXzy3DiRyIDV1RAODh2hx7fQRAwH+RsQEAg8V6pexkdCUBqUTGBgAAGIOMDQDEqKjsh9sLewE/8FTGJjs7W51yyimqatWqqlatWuqiiy5SGzdudLpZABzAPDLR48rZ8BNPBTYfffSRGjZsmFq6dKl6//331R9//KHOOeccPtQAoJR1OAQ9MIWnuqLmzZsX8vvUqVN15mbFihXqzDPPdKxdAADAHTyVsQmXm5urf9aoUcPppgBAXLU6dlcaXWuADzM2wQ4fPqxGjBihTjvtNNWiRYsil8vPz9c3W15eXopaCAAAUs2zGRuptVm7dq2aPn16iQXHGRkZgVtmZmbK2ggAAFLLk4HNzTffrN555x21cOFC1aBBg2KXzcrK0l1W9i0nJydl7QQAuK/g2R4llspLYyB1PNUVJS/E4cOHq9mzZ6tFixapY489tsTHpKen6xsAADBfOa91P02bNk29+eabei6bnTt36vuli6lixYpONw8AjCdZjuAsiZcnJeTyEWbyVFfUxIkTdXdS165dVd26dQO31157zemmAfCY4kYhcWFJ/2ESQ3N4KmPj1W8FAFAcAijAp4ENAMC/vNzthdTxVFcUAP9h4jozL7vAZRyQLGRsAMCHyH7AVGRsAACAMcjYAPAU6Y7yY7Ft+DDr4pCNgZ+RsQHgK04M5aZOCEgdMjYAkCR+zS4BTiJjAyClUpG9SORzxDN6x5QJ/pi0Dl5ExgYA/ocMS3JR+4NUIGMDwPOKy9BIoOK12pZ4rz7N3DAAgQ0AADAIXVEAPINuIm+h6wlOIGMDAACMQWAD+JjT86s4/fwAzENXFABfdFclY8RTLLMBA0gNMjaAociGAPCjmAKbnJyc5LUEgO8RjAFIaWDTrFkzNWbMGHXw4MFSPzEAeEm8c8sAcHFg8/7776v33ntPNW3aVE2dOjV5rQKAFHLrJRCiuaQBlz0AShHYdO7cWS1btkxlZ2ere++9V7Vr104tWbIkllUAgK+R+QFcWDzcr18/tXHjRtW7d2/Vq1cv1bdvX7Vly5bEtw4AACBVo6LOOeccdcMNN6jZs2er5s2bq1GjRvEtBDAEhbyR0fUDGDSPzaRJk9Ty5cv17euvv1ZlypRRLVq0UEOGDFGtW7dW06dP1wHOrFmzVPv27ZPXagAImkeGwAtAXIHNuHHjVMeOHXVX1KmnnqprbCpWrBj4++DBg9XDDz+sBgwYoNauXRvLqgHAk4In/mPCPsBjgU0089hcf/31urAYAADA8zMP16pVS3344YeJXi3gatSjIBEk40PtDuCywEY+2Lt06ZLo1QLwOLfOFeMUhn0DycFFMAEkVTIuPgkAReEimAAAwBgENgAAwBgENgAAwBgENoAPSI0LI7YA+AGBDQDfKm54NUOvAW9iVBSAEjGzbvKHfQNIDDI2AFyN+V4AxILABgAAGIPABkDCcGkJAE6jxgaA51DrA6AoBDaAj7i5CDiVl11wYvvduM8BE9EVBTjAjV02iZ7rhrlzADiBjA0AJBhDuAHnkLGBcdyYDUHqJ7+T52CIOOA/BDZAkpgYYIV3L5m4jcw4DHgbgQ3gMGpRkjeJn30/gQrgHwQ2gMsVlxUxMWNSFIIUAEYGNosXL1Z9+vRR9erV0x/mc+bMcbpJQEKQuQEAHwY28qHfunVrNWHCBKebAgAAXMZzw7179eqlbwBK9wUheLK+VE+QF969FNweAPBVYAMgOn4MGoqbP4a5ZQB/MD6wyc/P1zdbXl6eo+0B/Hi5hFh4vf0AnOW5GptYZWdnq4yMjMAtMzPT6SYBnpfsQmc/jfYCkFjGBzZZWVkqNzc3cMvJyXG6SQBchqHkgDmM74pKT0/XNyAZ/FS/AgBeUM6L/e+bNm0K/L5lyxa1atUqVaNGDdWwYUNH2wYAAJzlucDmiy++UGeddVbg99tvv13/7N+/v5o6daqDLQOSz5QCYQBIFs8FNl27dmXIJjzLS11XsQ4X98p2ATCb5wIbwC+8FAT5GfPjAO5i/KgowG9DjLnmFAA/I2MDJEii6l7I1ABA/MjYwJNMya74MfhjrhgAyURgAwAAjEFXFICkorgWQCqRsQEQguJjAF5GYAMAAIxBYAMAAIxBYAP4pMZFRiQBgOkIbOBLDBcHADMR2MBoFMKWHhkfAF5CYAMAAIxBYAPfI6uTWJLZIbsDwCkENkAMqM0BAHcjsIEvEJD8f62MG6/V5Oa2AfAWAhvAIW7qsqFAGIApuFYUgITj+lAAnEJgA2NIF5MUAgs782D/DgDwB7qiAA/XDKUycKO7CoAXENgASQw6JAgwuSCWYAeA2xDYwHck4Igl05HMEVXRBgZuKjQGADcjsIEjGH6duOHOBD0A8P8IbAAAgDEIbAAAgDEIbJCQ7iS6lbzV7cZMvwBMRWADGIDRSQDwX0zQB8SJyf+Sh5mLAcSLwAZATAg6ALgZXVFwbX2JG2tTAADuRsYGKCWpa6EAFwDcgYwNjKh1SXRWJxnrLC1GMgFAyQhsEEDXT/IRnABActEVBbh45JNXurkoKAbgFgQ2gEEIMAD4HYENEBYUSDeck5kaghMAiB81NvAkZtoFAERCYAMAAIxBVxRSyuluHgCA2cjYJADDpGEKhqMD8DoCGwAAYAy6ouA6dFcBAOJFxgYAABiDwAYAABiDrigYK3iiu0hF3fYcOHR7AYA5PJmxmTBhgjrmmGNUhQoVVMeOHdXnn3/udJMAAIALeC6wee2119Ttt9+uxo4dq1auXKlat26tevbsqXbv3q38jmHnAAC/81xg88QTT6hBgwapgQMHqubNm6tJkyapSpUqqSlTpjjdNFcj6ImMeVsAwCyeCmwOHTqkVqxYoXr06BG4r0yZMvr3zz77LOJj8vPzVV5eXsgNiBfXqAIAd/NUYLNnzx5VUFCgateuHXK//L5z586Ij8nOzlYZGRmBW2ZmZopai3BSpOtUoS4BCQD4g6cCm3hkZWWp3NzcwC0nJ8fpJqGE7jJGKQEAfDHcu2bNmqps2bJq165dIffL73Xq1In4mPT0dH1D8TP8JiqTEbw+alYAAKnmqYxN+fLlVbt27dQHH3wQuO/w4cP6906dOjnaNq9zspsoUfxUCOynbQUAYzM2QoZ69+/fX7Vv31516NBBPfXUUzrzIKOk4E1cGwoA4NvA5oorrlA//fSTGjNmjC4YbtOmjZo3b16hgmI/dSEV942doCExYulaC57xWDC8HgBSx3OBjbj55pv1De4KomJlry+VJ/7woAMAYBZP1dggMaOJZDkm6gMAmIjABgAAGIPABnFndYJH5jD5HQDADTxZYwN4CXU9AJA6ZGyQFNTxAACcQGADT3Db1cmZIA8A3InABgAAGIPAxkFuy0K48ara8jiyIgCAaBHYwNPBGF1CAIBgBDZAgrNMAADnMNwbIeyTOdeXYpg2AHgRGRsAAGAMAhsAAGAMAhsAAGAMAhuPj0aSWhjqYQAA+C8CGwAAYAxGRSFqicoMycgr5pwBACQDgU0CMCwYAAB3ILCBK5HVAQDEgxobAABgDAIbl3LbNZkAAPACAhsXIZgBAKB0CGwAAIAxCGwAAIAxGBVlGEYTAQD8jIyNz4IeuQEAYCoCmwSjABgAAOfQFZWCSxDQPRT97MzsKwBAaZCxQcKDl2R3d0nASEYMABAJgY2Pu6jsQCT4RrYEAOBldEXBcVxEFACQKGRsDO3ysTNEdp0PAAB+QGDjAxIAxdvNlKq6GQAAEoHABgAAGIMamxQI7w6SbqJEdxFFO0yaehYAgMkIbDyAOhkAAKJDYINCyOoAALyKGhsAAGAMMjZIODI+AACnkLGBJzDsHAAQDQIbAABgDAIbD0+eBwAAQlFj4wLJmNcmmaKdMwcAgFQjYwMAAIxBYAMAAIzhqcBm3LhxqnPnzqpSpUrqyCOPdLo5AADAZTwV2Bw6dEhddtllaujQoU43BQAAuJCniofvv/9+/XPq1KlONwUAALiQpwKbeOTn5+ubLS8vT3lV+GgkZvgFAMDDXVHxyM7OVhkZGYFbZmam003yJDuIYs4dAICbOR7YjB49WqWlpRV727BhQ9zrz8rKUrm5uYFbTk5OQtsPAADcw/GuqJEjR6oBAwYUu0zjxo3jXn96erq+AQAA8zke2Bx99NH6BgAA4PnAJhbbtm1Te/fu1T8LCgrUqlWr9P1NmjTx1CUJAABAcngqsBkzZox6+eWXA7+ffPLJ+ufChQtV165dlQmCRz5Joa7XriMFAICvi4djIfPX2CNzgm+mBDUAAMBHgQ0AAEBxCGwAAIAxCGwAAIAxCGwAAIAxCGySeOkBuckoJwAAkBqeGu4NcOFPAEBxyNikkGRvuIgkAADJQ2ADAACMQWADAACMQY2NoahFAQD4ERkbAABgDAIbAABgDLqiHOgSCr9PruANAABKj8DG5aiVAQAgenRFAQAAYxDYAAAAYxDYAAAAYxDYAAAAYxDYAAAAYxDYAAAAYzDc2wUY0g0AQGKQsQEAAMYgsAEAAMYgsAEAAMYgsAEAAMYgsAEAAMYgsAEAAMYgsAEAAMYgsAEAAMYgsAEAAMYgsAEAAMYgsAEAAMYgsAEAAMYgsAEAAMYgsAEAAMYgsAEAAMYop3zGsiz9My8vz+mmAACAKNnnbfs8XhTfBTb79u3TPzMzM51uCgAAiOM8npGRUeTf06ySQh/DHD58WG3fvl1VrVpVpaWlJTSSlGApJydHVatWTZmIbfQ+07dPsI3eZ/r2+WEb85KwfRKuSFBTr149VaZM0ZU0vsvYyM5o0KBB0tYvB9DEF2kwttH7TN8+wTZ6n+nb54dtrJbg7SsuU2OjeBgAABiDwAYAABiDwCZB0tPT1dixY/VPU7GN3mf69gm20ftM3z4/bGO6g9vnu+JhAABgLjI2AADAGAQ2AADAGAQ2AADAGAQ2AADAGAQ2MRg3bpzq3LmzqlSpkjryyCOjeozUZo8ZM0bVrVtXVaxYUfXo0UN9++23Icvs3btXXXPNNXoSI1nv9ddfr/bv369SLdZ2bN26Vc/eHOk2Y8aMwHKR/j59+nTlhHj2ddeuXQu1f8iQISHLbNu2TfXu3Vu/NmrVqqXuvPNO9eeffyovbKMsP3z4cHXCCSfo12jDhg3VLbfconJzc0OWc+o4TpgwQR1zzDGqQoUKqmPHjurzzz8vdnl57TVr1kwv37JlSzV37tyY35OpFss2vvDCC+qMM85Q1atX1zdpf/jyAwYMKHSszj33XOWVbZw6dWqh9svj3HwcY9m+SJ8pcpPPELcew8WLF6s+ffroWX+lLXPmzCnxMYsWLVJt27bVI6OaNGmij2tp399RkVFRiM6YMWOsJ554wrr99tutjIyMqB4zfvx4veycOXOsr776yrrgggusY4891vrtt98Cy5x77rlW69atraVLl1pLliyxmjRpYl111VVWqsXajj///NPasWNHyO3++++3qlSpYu3bty+wnLzMXnrppZDlgrc/leLZ1126dLEGDRoU0v7c3NyQ/dCiRQurR48e1pdffmnNnTvXqlmzppWVlWV5YRvXrFljXXLJJdZbb71lbdq0yfrggw+spk2bWpdeemnIck4cx+nTp1vly5e3pkyZYq1bt04fhyOPPNLatWtXxOU/+eQTq2zZstajjz5qrV+/3vrrX/9qHXHEEXobY3lPplKs23j11VdbEyZM0K+1r7/+2howYIDenh9++CGwTP/+/fXrIPhY7d2713JKrNsor7Nq1aqFtH/nzp0hy7jpOMa6fT///HPItq1du1a/bmW73XoM586da91zzz3WrFmz9GfB7Nmzi13+u+++sypVqqTPl/JefOaZZ/Q2zps3L+79Fi0CmzjIiy+awObw4cNWnTp1rMceeyxw36+//mqlp6dbr776qv5dDri8SJYvXx5Y5t1337XS0tKsH3/80UqVRLWjTZs21nXXXRdyXzRvAjdvowQ2t956a7Fv+DJlyoR88E6cOFF/MOfn51uplKjj+Prrr+sPnD/++MPR49ihQwdr2LBhgd8LCgqsevXqWdnZ2RGXv/zyy63evXuH3NexY0frxhtvjPo9mWqxbmM4CayrVq1qvfzyyyEnxQsvvNByi1i3saTPWLcdx9IewyeffFIfw/3797v2GAaL5rNg1KhR1kknnRRy3xVXXGH17NkzYfutKHRFJdGWLVvUzp07dYo0+DoXkm777LPP9O/yU7oL2rdvH1hGlpdrWi1btixlbU1EO1asWKFWrVqluz7CDRs2TNWsWVN16NBBTZkypcTLzrttG1955RXd/hYtWqisrCx18ODBkPVKl0ft2rUD9/Xs2VNfBG7dunUqlRL1epJuKOnKKleunGPH8dChQ/o1Ffz+ke2Q3+33Tzi5P3h5+1jYy0fznkyleLYxnLwW//jjD1WjRo1C3QDSLSpdjEOHDlU///yzckK82yjdp40aNdIXUrzwwgtD3ktuOo6JOIYvvviiuvLKK1XlypVdeQzjUdJ7MRH7rSi+uwhmKskbTwSf8Ozf7b/JT3nhBpOTiXxI2cukqq2lbYe8OU888URdhxTsgQceUN26ddP1J/Pnz1c33XST/tCSOo5Uincbr776av0BK33Lq1evVnfddZfauHGjmjVrVmC9kY6x/bdUSsRx3LNnj3rwwQfV4MGDHT2O0o6CgoKI+3bDhg0RH1PUsQh+v9n3FbVMKsWzjeHk9SivzeAThNRiXHLJJerYY49VmzdvVnfffbfq1auXPmGULVtWuX0b5UQugXOrVq10kP3444/rzxUJbuQixm46jqU9hlJTsnbtWv35GcxNxzAeRb0X5Qvfb7/9pn755ZdSv/aL4vvAZvTo0eqRRx4pdpmvv/5aFyOavH2lJS/UadOmqXvvvbfQ34LvO/nkk9WBAwfUY489lrATYrK3MfgEL5kZKVbs3r27/rA57rjjlEnHUT50pICxefPm6r777kvpcUTsxo8frwu45Zt9cHGtfPsPfs1KgCCvVVlOXrtu16lTJ32zSVAjX5omT56sg26TSEAjx0iyoMG8fgyd5PvAZuTIkbr6vDiNGzeOa9116tTRP3ft2qVPhjb5vU2bNoFldu/eHfI4GU0jI1Xsx6di+0rbjpkzZ+qUeL9+/UpcVtLF8uGUn5+fkOuIpGobg9svNm3apD9o5LHhlfxyjEUijmGqtnHfvn36W2LVqlXV7Nmz1RFHHJHS4xhOurzkm6m9L23ye1HbIvcXt3w078lUimcbbZLFkMBmwYIF+qRX0mtDnktes6k+KZZmG23yWpRgWtrvtuNYmu2TLwcSmEo2tCROHsN4FPVelC5uGcUm+6y0r4silapCx6diLR5+/PHHA/fJaJpIxcNffPFFYJn33nvPseLheNshBbbho2iK8tBDD1nVq1e3Ui1R+/rjjz/W65GRGMHFw8GV/JMnT9bFw7///rvlhW2U1+Wpp56qj+OBAwdccxyluPDmm28OKS6sX79+scXD559/fsh9nTp1KlQ8XNx7MtVi3UbxyCOP6NfXZ599FtVz5OTk6NfAm2++aXllG8MLpE844QTrtttuc+VxjHf75Fwibd6zZ4/rj2E8xcMyWjSYjM4MLx4uzeuiKAQ2Mfj+++/1EEt7SLP8X27BQ5vlzSfD4YKHJMrwNXkxrl69Wle5RxruffLJJ1vLli3TJ00ZauvUcO/i2iHDSWX75O/Bvv32W/2Gk9E34WQI8QsvvKCH28pyzz33nB4CKEPnnRDrNsrw5wceeEAHClu2bNHHsXHjxtaZZ55ZaLj3OeecY61atUoPZzz66KMdHe4dyzbKCUFGDrVs2VJvb/DwUtk2J4+jDAeVD/6pU6fqoG3w4MH6/WSPQLv22mut0aNHhwz3LleunD7hyVDosWPHRhzuXdJ7MpVi3UZpv4xYmzlzZsixsj+H5Ocdd9yhgx55zS5YsMBq27atfh2kOtCOdxvlM1YC8s2bN1srVqywrrzySqtChQp6SLAbj2Os22c7/fTT9UihcG48hvv27Quc8ySwkalP5P9yXhSyfbKd4cO977zzTv1elCkKIg33Lm6/xYvAJgYy/E4OaPht4cKFheb6sMk3i3vvvdeqXbu2PoDdu3e3Nm7cWGhOAznxSLAk38IGDhwYEiylSkntkDdY+PYKOYFnZmbqaDucBDsyBFzWWblyZT2/yqRJkyIu68Zt3LZtmw5iatSooY+fzAkjb9TgeWzE1q1brV69elkVK1bUc9iMHDkyZKi0m7dRfkZ6XctNlnX6OMr8Fw0bNtQnc/mGJ/Pz2CTDJO/L8KHqxx9/vF5ehpv+5z//Cfl7NO/JVItlGxs1ahTxWEkQJw4ePKiDbAmuJaiT5WV+kNKeLFK5jSNGjAgsK8fpvPPOs1auXOnq4xjr63TDhg36uM2fP7/Qutx4DBcW8Tlhb5f8lO0Mf4x8bsg+kS+EwefGaPZbvNLkn9J1ZgEAALgD89gAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgA8LRXX31VVaxYUe3YsSNw38CBA1WrVq1Ubm6uo20DkHpcKwqAp8lHWJs2bdSZZ56pnnnmGTV27Fg1ZcoUtXTpUlW/fn2nmwcgxcql+gkBIJHS0tLUuHHjVN++fVWdOnV0cLNkyRKCGsCnyNgAMELbtm3VunXr1Pz581WXLl2cbg4Ah1BjA8Dz5s2bpzZs2KAKCgpU7dq1nW4OAAeRsQHgaStXrlRdu3ZVkydPVlOnTlXVqlVTM2bMcLpZABxCjQ0Az9q6davq3bu3uvvuu9VVV12lGjdurDp16qSDHemaAuA/ZGwAeNLevXtV586ddbZm0qRJgfsl0JEuKemeAuA/BDYAAMAYFA8DAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABjENgAAABliv8DwgZ5A2NHyJsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.errorbar(x, y, ye, ls='', c='k')\n",
    "ax.set(xlabel='$x$', ylabel='$y$', title='Randomly Generated Data');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d98dfab-fbb7-4157-b7be-6e04c43833e3",
   "metadata": {},
   "source": [
    "I have generated some synthetic data with the relation \n",
    "\n",
    "$y = mx + c$\n",
    "\n",
    "where I have fixed $m$ and $c$. This is a very simple model, where I know the answer for the two variables I would like to fit. We can fit this with `lamatrix`.\n",
    "\n",
    "Our first question is what model is reasonable to fit this data. `lamatrix` has the following model options\n",
    "\n",
    "* constant\n",
    "* polynomial\n",
    "* step\n",
    "* sinusoid\n",
    "* Gaussian\n",
    "* spline\n",
    "\n",
    "From this list, we know that our data is drawn from a polynomial. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cba0696-bfe9-41be-b901-a796aacd18f1",
   "metadata": {},
   "source": [
    "## `Model` Class\n",
    "\n",
    "`lamatrix` provides a `Model` class to work with. `Model`s are classes that are (usually) initialized with the name(s) of any variables that the model is a function of. These names are strings. `Model`s may have additional keyword arguments that are required on initialization. `Model`s have the following class methods that you will usually interact with\n",
    "\n",
    "* `design_matrix`: The `design_matrix` method will return the design matrix that describes the model for all the input variables.\n",
    "* `fit`: The `fit` method takes in data and optionally errors to fit the functional form of the model to data\n",
    "* `evaluate` The `evaluate` method will return the best fit of the model to the data\n",
    "\n",
    "Let's look at an example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b5377e92-fa80-4bf0-808e-0a9f7aa99111",
   "metadata": {},
   "outputs": [],
   "source": [
    "from lamatrix import Polynomial"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "278c17ae-c7a6-42a0-a439-e52d954d2417",
   "metadata": {},
   "source": [
    "The `Polynomial` class is a type of `Model` in `lamatrix`. We will initialize it with the name of the variable of the polynomial, and the order of the polynomial. Let's make a model using this class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "eec9bfe8-d474-4230-aa8c-31b89762d306",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Polynomial(x)[n, 1]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Polynomial('x', order=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5e65ddf-182d-481a-9596-3dd214c8d4a8",
   "metadata": {},
   "source": [
    "This is an `n` by `1` polynomial model. Here `n` indicates the number of points in the variable (`x`) and `1` indicates how many vectors are in the design matrix. If we increase the order, we will see this number increase"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "03de4547-51be-4b0a-ab88-06d29c911219",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Polynomial(x)[n, 4]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Polynomial('x', order=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64bbdb0a-44bc-4334-8272-c35dbc5c5a19",
   "metadata": {},
   "source": [
    "We can check the equation of our model object using the `equation` property"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9aefefa0-01aa-46b2-9f8e-0f9cd6c6e76b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\\[f(\\mathbf{x}) = w_{0} \\mathbf{x}^{1}\\]</div>"
      ],
      "text/plain": [
       "'\\\\[f(\\\\mathbf{x}) = w_{0} \\\\mathbf{x}^{1}\\\\]'"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Polynomial('x', order=1).equation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1219af15-b049-4533-9718-4f85125e9245",
   "metadata": {},
   "source": [
    "This latex equation shows the model. Let's display the latex here:\n",
    "$$\n",
    "f(\\mathbf{x}) = w_{0} \\mathbf{x}^{1}\n",
    "$$\n",
    "\n",
    "This is a simple model of some unknown weight $w_0$ multiplied by $x$. \n",
    "\n",
    "To make this model fit our data we will need to include a constant term."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "cb1914bb-f909-41db-9892-b768b7e0a152",
   "metadata": {},
   "outputs": [],
   "source": [
    "from lamatrix import Constant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "6bebd5a7-0869-4031-a6ff-afba6e04738e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Constant()[n, 1]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Constant()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d4270d1-7508-4681-bbc2-706df00fa0a1",
   "metadata": {},
   "source": [
    "This `Model` does not require any input, there is no variable. If we print the equation we see it contains only an unknown weight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "d00a2270-c407-4d1a-a511-89a5f64b9127",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\\[f() = w_{0} \\]</div>"
      ],
      "text/plain": [
       "'\\\\[f() = w_{0} \\\\]'"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Constant().equation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ea51358-1f80-4e88-ad42-cf7ad34135e3",
   "metadata": {},
   "source": [
    "Let's combine these `Model`s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "dadfd9ce-ad04-4f92-8b01-38be332b19b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = Polynomial('x', order=1) + Constant()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "f54f1166-8177-4ba5-820a-9af752389f2e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "JointModel\n",
       "\tPolynomial(x)[n, 1]\n",
       "\tConstant()[n, 1]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5d8e08d-b0c3-4c4b-a1bd-141c286ef7c5",
   "metadata": {},
   "source": [
    "Now we have combined the model we have a `JointModel`. Let's look at the equation now"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "5a4720db-a0b8-47e3-82d4-db9e51dcb8d5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\\[f(\\mathbf{x}) = w_{0} \\mathbf{x}^{1} + w_{1} \\]</div>"
      ],
      "text/plain": [
       "'\\\\[f(\\\\mathbf{x}) = w_{0} \\\\mathbf{x}^{1} + w_{1} \\\\]'"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.equation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3673ccdd-8457-4739-9199-adb1da14aba0",
   "metadata": {},
   "source": [
    "This is now our first order polynomial, including the offset term. We have two unknown weights: $w_0$ and $w_1$ representing the slope and intercept of our model."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "578b2355-1231-468b-89fb-c138b7cf9c4d",
   "metadata": {},
   "source": [
    "Our model contains a design matrix, which is the matrix representation of the equation above. When dotted with a vector of weights ($\\mathbf{w}$) this gives us a realization of our model. We can access the design matrix inside the `Model` object. The `design_matrix` method will return the design matrix. You must input examples of all variables where you want the design matrix to be evaluated. In this case, that is `x`.\n",
    "\n",
    "This can be any `x` that could be generated. Here I've generated an short example `x` vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "34c8958c-c215-42ba-87e0-a45ace7eaf5d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.  ,  1.  ],\n",
       "       [-0.75,  1.  ],\n",
       "       [-0.5 ,  1.  ],\n",
       "       [-0.25,  1.  ],\n",
       "       [ 0.  ,  1.  ],\n",
       "       [ 0.25,  1.  ],\n",
       "       [ 0.5 ,  1.  ],\n",
       "       [ 0.75,  1.  ]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_example = np.arange(-1, 1, 0.25)\n",
    "model.design_matrix(x=x_example)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e685a230-f2f8-4766-9f5c-1a2d737c8a80",
   "metadata": {},
   "source": [
    "This design matrix represents `x` and a constant. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9799b55e-38a0-4c62-bd8d-176624eee839",
   "metadata": {},
   "source": [
    "Now we have the model, let's fit it to the data. We will input the variables (`x`) and the data we want to fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "bab19e73-1912-42b1-a694-95e13e452179",
   "metadata": {},
   "outputs": [],
   "source": [
    "model.fit(x=x, data=y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "559b490c-644c-46f5-a1b3-4acbc1505ba2",
   "metadata": {},
   "source": [
    "We can plot the model over the data to see the fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "635592e5-04f7-4ec1-8430-20e8e20e2868",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjYAAAHHCAYAAACskBIUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAABJkUlEQVR4nO3dC5xV4/7H8e801cx0VUhX126SSuVSqBBJcj8ijsrtoMidUCGdJJJLRx1Hcs5BjpDLQemIQiUlt5RbkUsl6aLLVDPr//o97T3/PdNMzWXPXnut/Xm/Xlv2mrX3ftZee/b6ze/5Pc+T5nmeJwAAgBCo4HcDAAAA4oXABgAAhAaBDQAACA0CGwAAEBoENgAAIDQIbAAAQGgQ2AAAgNAgsAEAAKFBYAMAAEKDwAYImL59+2r//fdPudfGrk2cOFFpaWlatmyZ300BfEVgA+zmQhG9VaxYUQ0aNHAX959++snv5oVGbm6u/vnPf+rEE0/UXnvtpUqVKqlOnTo66aST9Pe//13Z2dkKi02bNunOO+/UO++841sb7PVjP9dVqlTRvvvuq549e+rJJ58s0/v9+uuvu+cH/FTR11cHAuDuu+/WAQccoC1btmjOnDku4Hnvvff0+eefKzMz0+/mBdrmzZt15plnaurUqerYsaNuvPFG7bPPPlqzZo3effddXXXVVZo7d66eeOIJhSWwueuuu9z/d+nSxde2PPbYY6pWrZoLZCxQt3Nw8cUXa8yYMXrttdfUqFGjUgU2Y8eOJbiBrwhsgN3o3r272rdv7/7/0ksvdVmFkSNH6pVXXtG5557rd/MC7brrrnMXVLuYDhw4MN/PbrjhBn399dd66623lKy2b9/uMk6VK1dW0Jxzzjnusxw1ZMgQPf3007rooov0pz/9yQXxQBDRFQWU0LHHHuv+/fbbb/O2bd261V0Y2rVrp5o1a6pq1apuvxkzZuR7rNU/WPr//vvvd90sBx10kDIyMnT44Ydr3rx5O73WlClT1LJlS5cZsn9feumlQtu0ceNGFwjYX9n2fM2aNXOv4Xlevv3stQcMGKDnn39eLVq0UFZWljp06KDPPvvM/Xz8+PFq3Lixez3LKOyqXsOe2+ptTj/99J1+Ztktex/+8pe/FPn45cuX6x//+IdOPvnknYKaqCZNmrisTSwLJCwQOuSQQ1w7LcNjr/P777/n28/aduqpp7rs2hFHHOH2PfDAA123V0Fr167Vtddem/f+2Xtgwau9VmHnzl4/eu4WLVpUrPNvj997773d/1vWJtoVFJvdWLx4sQs4ateu7dprAbUF0AV98cUXOv744935a9iwoe655558bS2tCy64wAXvliWLDShnzZrlgh3rsrJjtvfJglLLuEVZF61la0xsV1eUvW+Wldtzzz1du+29mjx5cpnbDBRExgYooejFvlatWnnb1q9f7y7S559/vi677DJt2LDBdZ9069ZNH374odq0aZPvOZ555hm3j12Q7cv/vvvu01lnnaXvvvvO1ZiYadOm6eyzz3YByIgRI/Tbb7+pX79+7kJWMMA47bTT3EX0kksuca9lWZCbbrrJdTE8+OCD+fa3i5RdLPv37+/u23NbAHDzzTfrb3/7mwskLEiwNlnXxNtvv13o+2DtvvDCC91+1nVkF+OoV1991b0n9vOivPHGG8rJydnlPoWx98y6A+29uOaaa7R06VI9+uij+vjjj/X+++/nvX/mm2++cYGCvS99+vTRhAkT3AXYLqoWGEW7hzp37uzeK3tuu3h/8MEHGjRokH755RcXxMSyOhQL3C6//HJ3kbfjLs75t6DGun+uvPJK1/1m59u0atUqL1g5+uijXR3Xrbfe6oKj//znPzrjjDP0wgsvuMeYFStW6LjjjnPZouh+FiRbsBAPf/7zn93z2efP6p6MBcL2PlnbLTCxY3rkkUf0448/up9Fz8vPP//sAqJ//etfOz3vQw895D6nFjxZIDhp0iQXLFm3V48ePeLSdsDxABTqySeftHSHN336dO/XX3/1li9f7k2ePNnbe++9vYyMDHc/avv27V52dna+x//+++/ePvvs41188cV525YuXeqec8899/TWrFmTt/3ll19221999dW8bW3atPHq1avnrV27Nm/btGnT3H777bdf3rYpU6a4bffcc0++1z/nnHO8tLQ075tvvsnbZvtZ260dUePHj3fb69at661fvz5v+6BBg9z22H379OmT77WXLFni9nnsscfyvfZpp53m7b///l5ubm6R7+91113nHrtw4cJ82+19tPc7elu9enXez2bNmuUe8/TTT+d7zJtvvrnTdmunbZs5c2betlWrVrnjv+GGG/K2DRs2zKtatar31Vdf5XvOW2+91UtPT/d++OGHfOeuRo0a7nliFff82/HYcwwdOnSn9+OEE07wDj30UG/Lli152+z969ixo9ekSZO8bddee617jrlz5+Y7rpo1a+50vgpjr237WVsKY+22n5955pl52zZt2rTTfiNGjHCfr++//z5vW//+/d1jC1PwObZu3eq1bNnSO/7443fZXqCk6IoCdqNr167ur21Lv9tf//YXsmU8YjMn6enpeXUW1iVgGQz7i9q6EhYsWLDTc/bq1StfxifavWUZG2OZgoULF7osg3VtRNlf0JbBKViwaa9v2YtY1jVlsYxlRmKdcMIJ+YZsH3nkke5fyw5Vr159p+3RNhWmadOmbj+rzYiyY7fXtL/MY7siCrIsh7EC1oLHY+939Lbffvvl/cyyA/Z+2PuwevXqvJtlYOx5Cnb92XsVfW+NPZ9108Uekz2n7WPnI/Y57bxbRmnmzJn5ntPep2iXUmnPf0G2v2XGrGbLsj3RNliWzrI+VmsUHYln789RRx3lutdij8ve73iIng9rR1RsNsi6Pa1t1q1kny/LlBVH7HNYRnDdunXufS/O+wOUBF1RwG5Y3YBdwO2L2Loy7EJnXRAFPfXUU3rggQdcncS2bdvyttuIqoKsuyNWNMiJ1ol8//33eTUmBdmFOfZiYPvWr18/X1BiDj744HzPVdRrRwOngqNgotsL1q4UZMWmVrdjr2NBiAUKdvzWpbEr0fb+8ccf+bZbd0y0vmPUqFGueynKLvB2Hmw4eGFWrVq1y2ONvtexx2TP+emnn+4UrBT1nIWdz5Ke/4Ksy8yChMGDB7tbUe2wbip7n6NBZ8HPRTxEz0fs5+mHH35wNUQW0Bf8PNj5KA7rcrJaIAvYY4eU7yr4BUqDwAbYDfvLODoqyuodjjnmGPXu3VtLlizJ++v23//+t6vdsJ9bbYtdeO2veKtfiS0yjrKfFaZgsW95KOq1S9um8847zxWSWtbmtttuc++FvV+7u9A2b97c/WvD5lu3bp233QIMy5YYe65Ylg2x9zY2QxSrsEzK7o7JntMyQFZjVBgLamMVVstS0vNfULTw14a7W4amMFbQnAh2PmJfz7JW9v5YVumWW25x582ylpZBsmMuTtGy1XVZfU2nTp1cHVe9evVcLZTVK1m9GRBPBDZACUQvVla8aQWrVrxpbHSHjbh58cUX8/0FOnTo0FK9TrT7xbIJBVlAVXDf6dOnu66D2L+yLXMQ+1zlxYpnrfjTgg3rDrEMS8GC26KG0dv7GX1ccdhIJDtWy+rEq1jWntOyFNFgqjSKe/6Lyk7YY41d7HfXDjufxflclFa08DcaYNmIua+++splpCw7F1XYMPyijs+Kn22UlxW1x2Y7LbAB4o0aG6CEbBi0ZXHs4m2jY2IzA7GZABsyO3v27FK9hv1FayNp7GISm+q3i4kNL451yimnuL+qLdCKZaOh7EJjAUR5s24na5dlK+y9sCzO7lg3kY26snqcgm0vKltkNSh2rMOGDdtpX6tpsWHbJWXPaefJLroF2fPZ8+5Occ+/zfIbfd5YluGxz5UNt7f6qoJ+/fXXfOfb5pixkUmxPy8qi1USlj2x0V02BYDVYhV1bPb/NsqpIMvkFHZ89hz2WbRzFzu60KYzAOKNjA1QCnYBt6GqNuz4iiuucMOl7a91G5Jr2Qsbgjxu3DhXvFqwhqS4LDNkz2VdXxYAWFeADbG1Ycqxz2lT4VsG6fbbb3cXC+vWsaG6L7/8spubxTIS5c3aacOArb7GAqmiamAKsuDQ3qurr77aDf+1Y7HHWnGqZX5s2Hhsl5YNy7ZhxfbeWK2GLbtgWQ7LYNhr28XWCrxLei6tdsTOYXQouBXIWqbCMjH2nsZOZFeY4p5/yzLZtueee851cVm2y+YnspvVctm5PvTQQ92QccvirFy50gVHNqz6k08+cc9hXWaWVYnO/xMd7m2ZHKsVKi47NutKtaHX0ZmH7T23z090CLexrif7DFk3me1Xo0YNl4EprPbK3jtjheyW8YkGufaejB492rXZunGtXsiO17q7StJmoFhKPI4KSLHh3vPmzdvpZzk5Od5BBx3kbjbU14bl/vWvf3VDjG048WGHHea99tprOw2Pjg4ZHjVq1E7PWdgw4BdeeME7+OCD3XO2aNHCe/HFF3d6TrNhwwY3fLp+/fpepUqV3PBge42Cw63tNWxIbqyi2jRjxgy3/fnnn8/bVthrR1111VVu/2eeecYrCXv/7L22Yb+1a9f2Klas6O21115u+PO4ceO8zZs37/SYv//97167du28rKwsr3r16m6Y9M033+z9/PPPeftYO3v06LHTYzt37uxuBd8/G97euHFjr3Llyu71bZj1/fff74Yl7+p9MsU9/+aDDz5wbbfXKXjOv/32W++iiy5yQ+/tPDZo0MA79dRT3TQDsT799FN3DJmZmW4fG7L+xBNPlGi4d/Rmz9GwYUP3OhMmTMg33Dxq0aJFXteuXb1q1aq59+ayyy7zPvnkE/d4O3ex5/Lqq692UyLYUPDYS4y1zz6X9v40b97cPS7aFiCe0uw/xQuBAKBoVkBsk9LZBHLRLhcASDRqbACUmdUa2cggm+OFoAaAn6ixAVBqVitho5SsXsMmkytqzScASBQCGwClZiOhbKi2Ffw+/PDDO62JBQCJRo0NAAAIDWpsAABAaBDYAACA0Ei5Ghtb1+Tnn392U8+z+BoAAMFglTO2dIwt+luhQtF5mZQLbCyoKbiKMQAACIbly5erYcOGRf485QKb6CKB9sbY1OAAACD5rV+/3iUmYhf7LUzKBTbR7icLaghsAAAIlt2VkVA8DAAAQoPABgAAhAaBDQAACA0CGwAAEBoENgAAIDQIbAAAQGgQ2AAAgNAgsAEAAKFBYAMAAEKDwAYAAIQGgQ0AAAgNAhsAABAaBDYAACA0CGwAAEBoENgAAIAS27hxo9LS0tzN/j9ZENgAAIDQILABAAChQWADAABCg8AGAACEBoENAAAIDQIbAAAQGgQ2AAAgNAhsAABAmeeuSZZ5bQhsAABAaBDYAACA0CCwAQAAoUFgAwAAQoPABgAAhAaBDQAACA0CGwAAEBoENgAAIDQq+t0AAAAQbNWqVVOyIGMDAABCg8AGAAAUKlmWSSgJAhsAABAaBDYAACA0CGwAAAiQIHYPJRKBDQAACA0CGwAAEBoENgAAIDQCHdjce++9ro/x2muv9bspAAAgCQQ2sJk3b57Gjx+vVq1a+d0UAACQJAIZ2Pzxxx+64IIL9Pjjj6tWrVp+NwcAACSJQAY2/fv3V48ePdS1a9fd7pudna3169fnuwEAEEYbSzEUPGzDxwMX2EyaNEkLFizQiBEjirW/7VezZs28W6NGjcq9jQAApEogkWwCFdgsX75cAwcO1NNPP63MzMxiPWbQoEFat25d3s2eAwCAVFhxOy0SQKVSEFVRATJ//nytWrVKbdu2zduWk5OjmTNn6tFHH3XdTunp6fkek5GR4W4AAKBsgVIQBCqwOeGEE/TZZ5/l29avXz81b95ct9xyy05BDQAAKB9Vq1aV53kuE5RMQU+gApvq1aurZcuWO72xe+65507bAQBA6glUjQ0AAGGTjMXEGyNtSqZMTCgzNoV55513/G4CAABIEmRsAABAaBDYAACQQsO/V61alXRdX/FEYAMAAEIj8DU2AACgdKqVQ3HwvvIXGRsAAAI6eimZ7C9pW8+e+j4zU1VXr/atHQQ2AACg2P744w83MZ/NI+ds3Khhkr60bqBXX5W2bZNmzJBfCGwAACgCWZpd8Dzp2WeVddhhukOSreCY06WLtHCh1Lev/EJgAwAASmbBAunYY6XevVXh55+1VNKZkrZYxsbnlQAIbAAAQLHsLanygAFS+/bS++9LVapo65AhOljSFNshLU1+I7ABAAC7HUI9UNJXkipNnLijG6p3b2nJEm27+WZlK3kw3BsAABTpREljJLWI3M9p00bpjz4qHX30jg1JVntExgYAAOzkwEj30rRIUPOrpMusjubdd/8/qElCBDYAAOD//fGHKg0dqkWSTre5aSQ9KKmJpH/Yz9PTlcwIbAAACOgw8ej6T/F4jTSrT3n2WalpU1V+4AFlSJoqqZWk6yWtK+JxNp+NzWuTb24bHxHYAACQ4vPntJf0vqSMyy6TfvlFuQceqNMknSxpsYKFwAYAgCQWzcrsKtApbeamjqQnJM2T1MHm3LOMy4gR2jxvnl5VMDEqCgCAFFNJ0jWShkiqEdn2T0nnLFyoKo0b5xvpZEsolNeCmeWBjA0AACmk8vTp+kzS/ZGgJpqt6WMZm3r1FHQENgAApIDGkute2uOCC9RM0gpJ/SQdKWmOwoPABgCAAImOQop2Ee1OdUkjJX0h6VTLylSsqFGSmkqaaPcVLgQ2AACEUFqke8mWQbjZuqAkvS5pzbvvuvsbFE4ENgAAhMwRkmZHMjJ1I8HNKZJ62JIIVhwcYgQ2AICUma8l7CqsXOmCmbmR2pn1km6U1FLSG0oNDPcGACDgrJvpWkm1O3Rw3U/mSUmDJK1UaiGwAQDAJ/GYG8YKgkdH1nLSxo1uhNM1kWHcqYiuKAAAAihtyRLXvfRqJKj5xbqeHnlEHVM4qDEENgAABMm6ddL11yvryCPdWk7Zku6NDN9O79dPuSUYCh5l+yfLIpZlRVcUAAA+zENjxdIl6YpKi0yoV6V1a2n1anf/lcjK29/u4nF/BGxJhLIiYwMAQBGjteIZDJRlBJgtefBhZMHKtNWrpWbNtGXKFJ1eRFBTNRI8xSsLE+/nK08ENgAAlKPSrrxt6kv6l6QPJLW3Xijrerr3Xumzz5TTtWu5tDfoCGwAAKERmnlxtmxxQ7WXSLpQUq6kxyNFwtsHDJAq2frc8bUxJO8dgQ0AICUk84U7r3i3ShVpyhRlHX64/mrZHknvSzpc0uWSfvW7oQFA8TAAAGUUWwhc0hFJeRYtkgYOlKZPd1mHnyJrPD2jxKsaqakJIjI2AAAksG6moD1s5uCbbpJatXJBzRZJwyU18ymoCToCGwAA4hz0FGs0VU6O616yBSorPfaYu7/91FPVQtIdlgWKY5uqxoxqqlOnjvu31JmlJEdgAwBAglV4/31lHnusxkva24qDmzeX3npL2ZMmaanfjQs4amwAAIFX0snu/NJQ0ihJWd26ufu/SxoqacTs2aq6xx5urSe/6mg2JllBdWmRsQEAoJxl2ozBo0e74dvnSfLS0rTt4ovdMgiP2A7lMHw7VRHYAABQggxHSWtTzpb0pdXejBypKpJm2jQ177+vrQ8/rNXl1trURWADAEA5aCnpf5ImS9rfaoUbNFAvSZ2tpsZGQKFcENgAABDHla9rRbqXFko6XtJmSXdJ+u299/Sfcm8tCGwAAKGaxdcv6ZKulPS1pAGR+5atOVjSnbaDzSqMcseoKAAAyqjCzJlaICnawfSZpGskveNzu1IRGRsAAEop7YcfpD/9SVmnnOKCmjWS+ks6jKDGN2RsAAChVh5z3GRJusX+bdvWrcTtVaigv+XmakgkuIF/ApWxeeyxx9SqVSvVqFHD3Tp06KA33njD72YBAFKF5yl98mQtjkysl7Zli9SlizZ/8IGrqyGo8V+gApuGDRvq3nvv1fz58/XRRx/p+OOP1+mnn64vvvjC76YBAEJUFF1ohmfhQhfEZPbtq30lfW/z0fzrX9Lbb8traYO7gzGvTtgFKrDp2bOnTjnlFDVp0kRNmzbV8OHD3Ydvzpw5fjcNABDWFbhXr5auuEJq106aOVNeVpbrcmpuc9OceaaUllbowyzgIOhIvMDW2OTk5Oj55593H1jrkipKdna2u0WtX78+QS0EAIRh+HaVNm2ktWt3bOzVS5vvvFPDDrZB3OXLgqLizJ1TcM2nVBeojI357LPPXBSekZGhK664Qi+99JJatLBF3gs3YsQI1axZM+/WqFGjhLYXAMpDWOeRSZbjOj4ywZ5NtJdmQU3r1tK770qTJsnjOpLUAhfYNGvWTAsXLtTcuXN15ZVXqk+fPlq0aFGR+w8aNEjr1q3Luy1fvjyh7QUA+NyVVAJpy5bphchSCFY1Y2s5ZY8ZI82fL3XqlLB2IIUCm8qVK6tx48Zq166dy8a0bt1aDz30UJH7W2YnOooqegMAJGemxDd2zHfcoax27XSWpO2S7MrSxP7/0kuldOuYQhAEtsYmKjc3N18NDQBg9/O5FLd+I5HiPddMsVhtyrPPSjffLP30k6wMeLqkgZKK7gvYuc0rV64s96YihIGNdSt1795d++67rzZs2KBnnnlG77zzjqZOnep30wAAQbNggXTNNdL77++4f8AB2jJ8uE7s3dvvliFVAptVq1bpoosu0i+//OIKgW2yPgtqTjzxRL+bBgAIyKzCe1tZw4AB0lNP7cjY2OKUt90m3XCDG3GLYAtUYPPEE0/43QQAQIAveAMiMwZXmjhxx0bLzowcaTPA7rifivVFIRO44mEAQHiKjxM1e+5Jkj6V9KCkPWwuNJub5r33pKef/v+gBqFAYAMACO2Iq4Mk1fzzn2WVmDal3ipJl9pSCDYnzdFHy0/RoM5uyVbIHWQENgCA0LGqm79KspUEM6ZN0zZJoyU1tbIG24Hh26FFYAMACA0brn2hpCU2ktaCGptg77jj1ErSDZLW+d1AlDsCGwAIy+rTIZ0xuLhdWe0l2cDtf0mqL+kbWzzZgplnn9XiuLQEQRCoUVEAABS0T6Tb6eLIfStDHiZpjKSttqGI1bcTjcUqE4PABgAQTFu3quJDD+krSdHFcv4p6VZJv8Th6ZNxdmbsHoENACB4Xn9duu46ZXz1lauj+VDSNZLmKrkEKUtTNUBt3RVqbAAgpEo7zNqPVbWLK+3rr6UePXbcvvpKuXXqqJ+ko2KCGsu0FDUvjm1neHW4EdgAQJyk/ArZ5ai6pPskZR1xxI5sTaVK0o03avPChbI5hP3OMyRqokHsHl1RAICkZWW/mZMmuTqaurZh2zbplFOkBx+UmjYN3RIIYekO8hOBDQCkgKKGh5dmEclEOULSw1YYPHCgKw624GbfF15Q5llnJUXwUdh7R2DiP7qiAABJpZ6kpyI1M0dKyq1WTTdKamlrPHXr5nfzkOQIbAAASaGypJsjswZfFNk2QdKa2bP1gPVC+dw+BANdUQAAf3me0l9/XZ9LahLZNDsyfPsjSSvr1CnxvDMUb6cuMjYAAP8sXix1767Mc891Qc0vkWzN0ZGgBigpAhsAQOKtWyddf7106KHS1KnyKlfWvZHVt22tp12V3zKsHrtCYAMASJycHOkf/5CaNNkxZHv7dqlnT22eN8+txp2IWWAIjMKNGhsAQLHEDm9euXJliR/f0eak6dxZWrhwx4bmzaUxY6Ru3eT5FGBEj4d1ocKDjA0AoEwz7VpwsKu5cNJ+/ln/lvS+pHQLamrUkEaPlj791AU1QDwR2AAAyseWLdLw4cpq00YX2Hw0NmS7b1/J1nu67rodyyIkiWj3VLJOVojioysKABB3Z9i6Tu3bS8uWuWUR3o8M35756KOqlOAun4KzAVNXE25kbAAAcdNC0luSXrILzLJlUoMG2jJhgo6RtCAOz7+7bi+AwAYAUlBJul6Ks+8eksZI+kRSV+uFkrT1ppvcPDU5555b6PMxKgnlgcAGAFB6OTm6XNLXkgZG6hteimRutg0daikWBY0VRFvXFaOkgonABgCSXEkzHAkrhJ05U5nHHqvxkvaS3JIIlq2xtbeXyl/RuhoClNRD8TAAoEQaSsro00d64QWlS/pd0hBJj1kCx+/GIeWRsQGAgEp0l0mmpMGR1bcrvvCClJambRdf7JZBeDSOQY0dl92A0iCwAYAiUOQa4Xk6W9KXku6WVMWCmI4dpQULtPXhh7W6HF6SriSUFoENAKDQTIkFc4empenD6tU1WdL+kn6Q1MtGPU2dKrVpo1RG8JWcCGwAIIWyS8Xtvqptc8bceqtsVafjJG2WdJct7yTpP7ZDmk27F795aAgSEC8ENgCQZHztAtu+XVdJ+sq6nJ580hUHPy/pYEl3RgKcZEJAhIIIbAAAToV331XW0UdrrKQ9LcY5+GCXrbHp9b73u3FAMRHYAECK2zfSvZTVo4cqfPGFfpNc1mbN9Ol6pwzPy8KS8AOBDQCkqKxI3cxiSX+ywU8VKmjb5Ze74ds2J40qlmyqM4ZpIxkQ2ABAqvE8pU+e7AKaIZEAZ4bVz3zwgbaOHq015fCSBD1IFAIbAEihOW4qfPqp1KWLMvv2dV1QyySdI+l4i3datvS7eUCZsaQCAKQAKwa+x2YPPuYYKTdXXlaWhm7erFGRlbiBsCBjAwABUtLhzfbX69WR1bevsOlncnOlXr20ecECDSuHoKY8hl9Hn5OuLBQHGRsAKAbrboqO7knkBbYsI4pOkPSQpEMi922yvWZvvqmsbt3kBaj7DCgJMjYAEODanMKGVFdYtkwvSpoeCWpsLae/SGonKde6ooAQI7ABgJCoGqmj2bNTJ51pE+xJ2nbllcpavlx/t6CmOM/BTL4IOAIbAAiB3pH5aG63OprsbL0lqbWkraNGSbVq+d08IGEIbAAgwAs5Vvj4Y70n6WlJDSV9J2ntk0/qJEmLEtICILkEKrAZMWKEDj/8cFWvXl116tTRGWecoSVLlvjdLAA+COo8MnGzapV02WXK7NRJR9v7Iek2SS0sS3PKKfl2pXsJqSRQgc27776r/v37a86cOXrrrbe0bds2nXTSSan5pQYgNW3bJj34oNS0qfSPfyjN81y2ppn98ScpuxRPaaO8CHoQFoEa7v3mm2/muz9x4kSXuZk/f746derkW7sAIBGseynryCOlr77asaFtW20eOVIXnnii300DkkagMjYFrVu3zv1bu3Ztv5sCACVmdTrRrrRdda0dJOllSVPtS9uCmr33lh5/XPrwQ+V26OBb+4FkFKiMTazc3Fxde+21Ovroo9VyF+ubZGdnu1vU+vXrE9RCACibapFRTtdJyrBeKNs4YIAqDRsm7bGH380DklJgMzZWa/P5559r0qRJuy04rlmzZt6tUaNGCWsjAJRKbq7+LMmGRtwaCWqsI76VFQbfey9BDRC2wGbAgAF67bXXNGPGDDVsaAMcizZo0CDXZRW9LV++PGHtBICSai+pVo8e+qek+pK+kdRTUvfIPDUoe8Eza0+FW6C6ouyDePXVV+ull17SO++8owMOOGC3j8nIyHA3AEhmaStXaoKkfnZnwQJtiMwiPMayNH43DgiQCkHrfvr3v/+tZ555xs1ls2LFCnfbvHmz300DgFKpZLcxY5TVps2OoEbS5nPPdcO370vCoKZgpiTIc+QEue0ISWDz2GOPue6kLl26qF69enm35557zu+mAQiYXY1CKmxhyfJg3UufS6p8xx1K27BBH0o6StKGRx7RL+X6yiiIICc8AtcVBQBB10TSg5J6RO6vsHpASU/Z95zPbQOCLlCBDQAE2vr1qjRkyI4sTaSbaUyklsZqalC8rAoQmq4oAKknFGtC5eZKTz7plkGo/NBDLqh5TZLNwHVLQIKaeC+7wDIOKC9kbACgPM2ZI11zjTRvnrub26SJTv36a73hc7PIfiCsyNgAQHn4+WfpooskW/LAgprq1aVRo7R57lzfgxogzMjYAAgU644q79FKZWJLuNjq2/fcY43dsa1fP+mvf5Xq1v3/baXouilutw3ZGKQyAhsAKaXcAiMLJF59Vbr+eunbb3dsO+oo6eGHtbFFi7zXXLlyZfxfG0AeAhsAKKO0xYul226Tpk3bsaFePWnkSG084wxVq1HD7+YBKYXABoBvGZOSdK/49RrRxxc2IqumpKGSsiwzs327VLmyy9hsHDhQ1Sy4CTi6tBBEBDYAUMJuKht1Ycsf/FVSHdtgQc1pp0kPPCA1blzqOpqwI1BCIjAqCkCo57qxQCWec+B0lNzSB/+IBDVfStoyZYr08ss7gpo4KO3q08wNAxDYAECxpP30k/4t6X1J7SStlXStpFaScrp29bt5ACIIbAAEhmVfEj7Ue8sWafhwZR12mC6wCfYk/V1SU0kPWS9UYlsTKCwsCT9QYwMARTjDCoPbt5eWLVOapPckDZS0wO+GASgSGRsghfm9DpPfr1+UFpLekvSSfUkuWyY1aKAtEyboWIIaIOkR2ABIie6qaBC1y66s33933UufSLKqmS22AvdNN0lLlijn3HN32p1iXSD5ENgAIZWs2ZCklJMjjRunKq1b65pIH/2LkczNtqFDrVjE7xYCKI8am+XLl6tRo0YleQgAJNXkfQVZ91LmMcdIn33m6mg+j9TRvF3urwzA94xN8+bNNWTIEG3atKlcGgMAiWJ/ok2SNFNS+mefSXvsoexRo9SmiKCmtHPLAEjiwOatt97S1KlT1aRJE02cOLH8WgUA5SRT0hBJiyX1sl4oSZv69NFea9cq86ab3P0gDZdmSDVQhsCmY8eOmjt3rkaMGKHBgwerXbt2mjVrVkmeAgD84Xk6OxLQ3CWpiqR3I5Pt/XHfffotQc0g8wMkYfHwRRddpCVLlqhHjx7q3r27zjnnHC1dujT+rQOAODhU0h5nnaXJkvaT9IMkG+PUJTICCkB4lGlU1EknnaRLL71UL730klq0aKGbb76Zv0KAkAjDqKrakh6V9LGkyh98oM2RbE1zSc+X8jnp+gFCNCpq3Lhxmjdvnrt9+eWXqlChglq2bKkrrrhCrVu31qRJk1yA8+KLL6q9zdYJAOUoOnKqYOCVLukvkoZFghuz5bTT1PyVV1y2BkB4lShjM3z4cK1bt851Rc2YMUNr167V/PnzNXbsWF1++eV6++23XZDTt2/f8msxAOxCl8jswGMjQc2nkW3rH3+8XIKaMGS2gJSex2Z3LrnkEldYDACJZLUz90s6J3LfioEHRxasTKaRTgACtghmnTp1XOYGSCV+TCyHiE2bVGnYMH1pC1ZGgpjHIkO6f1ew8NkBkjCwsXRs586d4/20AEIU/MVL+vPPS4MHq/KPP7r7b0dmDbbZg5NdtAi5PN4XIJWxVhSAclWsxSdLqE1kDprMfv2kH39U7r77ujlqTghIUAOg/BDYAAiMPSPdTPMldbI597KypLvv1ub5892ilQAQ964oACiPL6orI3PQ1Ipse1bS6R9/rCrNmllayOcWAkgWBDYAklqFGTO0UNIhkfs22d41kt6zYtuGDX1uHYBkQ1cUkAKsviVo86yk2TItZ56prJ49XVDzq6TLJbWPBDUAUBgyNgCSig12HmRDt2328uxseenpeignx3VDrU3QzMWxPwMQLGRsACRHxsfzlP7cc1oi6XbL2GRnS127avPs2bquHIKaZMHaU0B8kbEB4Lu2NnT7xBOVPmeOGkj6TlL9Z59VZq9eqpKWxnwvAIqNjA0A36T9+qselzTPJtubM0delSq6TVILm0G4Z0+b8dPvJgIIGAIbAAlfELKS5LqX9uzQQZdGvoi29+qlzQsXaoSk7IS2GkCYENgASKhukRW3R9sX0IYN+khSRwtmnnhCXv36xXoOK+ylJgVAYQhsgBTi57DvgyS9LOlNSc0lrZK0fvRoHSFpdqRtiaqh8SMwIhgDEoPABkjiLpsgBz3R59u0cqXulbRI0mmStkl6QFITSVsuuEBeXF4NAHZgVBSAcmFlvxdKqt2hg26JbHsjUltjQ7rDLDqEG0DikbFB6CRjNiTVHC7pA0n/tNFOq1bpG0mnSjolEtQkolvGXsNuAFILgQ1QTsIYYBXsrip4jGkrV2qCpA8lHSVpgwUYd9zhlkT4r4KBWhgg2AhsAJ8FcR2nwoZvVxozRllt2qhfZNtTkppK2nT11dqaoK6fghkaZvUFUg+BDRDgzE8yZIWse+lzSZXvuENpGza4bM2RkvpKWhHH1yFIARDKwGbmzJnq2bOn6tev777Mp0yZ4neTgNTM3CxZooyzznJdTJaZya1TR9njxrkuKAtuAMAPgQts7Eu/devWGjt2rN9NAVJSdcvO3Hab1LKlKk6b5rqZ7pPcrMHbL7yQ4dsAfBW44d7du3d3NwClF7ugZLQuZXeT46VFupdsyYNKDz/stm0/+WS1fPNNfS2pf40ape5eYoFLACmbsQFQtoLa0rDupbmSG/G0j3U7NWkivf66sidPdkFNsthVHQ41OkBqCFzGpqSys7PdLWr9+vW+tgcoL7EZmHhduOtJbtbgiyL37bfnLkl3z52rqrVqWepH8RbP9gNIPaHP2IwYMUI1a9bMuzVq1MjvJgHJLztblR54QF/FBDVPRJZBsMUrq9WuXa6Fzskw2gtAMIU+sBk0aJDWrVuXd1u+fLnfTQKSWvp//ysdcogqDx2qapEFKm0m4UsjC1eGEd1UQHiEvisqIyPD3YDyEKaiV1txe4ykzF693P3cunV10YoVekZipBOAwAhcYGP97998YyvP7LB06VItXLhQtWvX1r777utr24BAWrtWlUeM0KeRGYS9ypWVdv312jxwoJ6uZ1U2ABAcgQtsPvroIx133HF596+//nr3b58+fTRx4kQfWwYEq0DY+qEvllSlTRulrV7ttr0s6cR581SlVatyKQwGgPIWuMCmS5curh8cCKJk6bqq9OGHmieprd1ZvVq5zZrp5CVL9JYFTQcdVKo5ZpLhuAAg9MXDQFBFRwbFM2BoIOlpSbV69nRBzVobADVypDbPmeOCGpQchcdAciGwAUI2xLiwNaesfP52W95JUm+ro0lL098jazxt799fqmTVNQAQfAQ2QJxY3Us8/mqPa6bG85T+yiv6UtI9ll2Q9J6k36dN018k/Vr2VwCApEJgg0AKS3alPKUtWiSdeKIye/fWAZJ+lHS+pGMtS2PFwQEO/gAgNMXDAHZtj8iyB1kdOkg5OfIyMnRPdrZbGmGT340DgHJGxgYI0S+zdS/ZopTXWMYmJ0c66yxtnj9fQ3wMaiiuBZBIZGyAELDupYcltYnc/1zSQa++qqxTT5VXwq668lhMEwAShYwNEGBpy5drkqSZkaDmd0lXR/4/N2YiSwBIFQQ2QBBt3izddZey2raVreyUI+mxyOrbj0buA0AqoisKCJhzrDC4bVtp+XKlSXo3UlNjaz0VpqQzCANAkJGxQUoK4nDxQyW9Lel5+8Vdvlxq1EhbnnpKXXYR1ABAqiGwQagVNgtv0NSWNFbSx5KsamazpK2DBkmLFyvn7LPL/fWjGR8rJgaAZEdgAySpdElXRYZvXxW5/x9JzSVtu/12qUoVv5sIAEmHwAYpLxmzOsdFMjRjIxmbT2xle8kVCv+g5GaZHbI7APxC8TBQArEFuOVx8d5PUo1LLnG1NOY3SXdIepyRTgBQLGRskBKSvVi4SmQZBFusMvO111wQ80hk+Pa4OAU1yTwDcDK3DUCwENgAfnbZbNjgupcWS27ZgywrDD76aDfB3jWRCfcSgQJhAGFBYAP4pMInnyizWzc3c3AjScsk2RintS+84JZECDIyMAD8QmCD0HU3WQ1MMhew7hXpXso85hilf/CBW5xysKSDJb1oO6TZtHsAgNIgsAESZds2VRw7Vl9FVuFO8zxtP+ccNZN0j6QtpQziEoXuKgBBQGADlHPmyHXHfPCB1KaNMm65RbUiQ7k3T52q7IkT9aOCi2AHQLIhsEHKsYCjJJmOsoyoSlu6VDrjDOmkk6RFi+Ttuacul9TeVt8++uhiBwbJ3LUGAMmEwAa+SPbh12Vl5bLDbZRTu3bSyy9L6enSwIHatHChm5MmN47FtgQ9APD/mKAPiCfPU/pzz2mJpAZ2f+tW6cQTpTFjpBYtLKLzu4UAEGpkbIA4qfDxx9IxxyjzkktcUPOtFQRPmiRNnbojqAEAlDsCG8SlOynM3Uq7s3dkyYPMTp2kDz6QV7WqbpN0iM0YfOqp2rhpU9K9P8wzAyCsCGyA0tq6VRUfecStvn1pZPi2LrxQmz/+WCMkZSewKYxOAoAdqLEBSqGbpOW1a6u5pAxJH1mGZvp0ZZ1wgrwkycoEWTRQA4CSImMDlED6d9/pFUlvSi6oWSnpEklH2Eino45SKqAbC0AyI7BB0g7rTqYh4Tbrzb2SanfqpJ42ibCkByQ1lTTBBkP52joAQBRdUcCu5Oaq4tNPu2UQ6tn9bdv0hqTrJDek21hdC5kLAEgOZGwQeDaLcLyzOvacR6SlKeeoo5Txl7+4oMaKhNf+6186JSaoSSS6gABg9whskJRdP37aJ9K99KHV1MybJ69aNd0sqaUNhLKlEcqA4AQAyhddUUDU1q26UdJgSTUimyZKGvTHH1rhU5OC0s3FKCYAyYLABpBc91Ltzp01KnJ/rqRrIlmbICHAAJDq6IpCSrNRTf+N3Cp+951Ut66yx49XBx+DGrqrAKD0CGwQSGWdade6miw783kkW7PVbtddJy1Zou0XXMDwbQAIKAIbpJQ0Sf0kN3zb6mkqSXotUhi8bdgwqUa0ugYAEETU2CChbLSVDaX2g80L/LCkwyP3bcj2tZFZhAEA4UBgE+eLdVBGsaQSm4OmRv/+mh25v17SXZIeicwgjP9H8TGAoKMrCuG1ZYtujXQ7ZU6erFxJT0hqImk0QQ0AhBKBDZJ2osCydFml//e/yjr8cI2IrPO0rX17t1DlpZJWxbW1AIBkQmCDUGkeqZnJ7NVLFZYu1c+SLpT0+2uvab7fjQMAlDsCG4RCTcvM3HGHPpPUzVbbrlxZW2+80c1T87TtkGbjoQAAYUdgg2DLyXHdS7ZAZZXHH3fV8FMkbf7oI1UeNUp/7GKiOyv0Lu08OACA5BTIwGbs2LHaf//9lZmZqSOPPFIffhi0ie8RF++9p8xOnfS4pL0lbW/aVLZE5ZmWsTnwQL9bBwDwQeACm+eee07XX3+9hg4dqgULFqh169bq1q2bVq2iJDRlVuf+8Uepd2/p2GOV/sknWhuZj2bN22/rLb/bBgDwVeACm9GjR+uyyy5Tv3791KJFC40bN05VqlTRhAkT/G5aUgtF0LNli3TPPVKzZtKzz7q6mW19+7rh2w/ZzyvZPMIlw7pMABAugQpstm7dqvnz56tr16552ypUqODuz54dnX4tv+zsbK1fvz7fDcGT/vLL0sEHS4MHS5s2SUcfLX30kbY++qhWB2iNKgBA+QpUYLN69Wrl5ORon332ybfd7q9YsaLQx4wYMUI1a9bMuzVq1ChBrUVBNi9NSeemOUTSdBu+fcEF0rJlUoMG0jPPSLNmSW3bFvt5CEgAIDUEKrApjUGDBmndunV5t+XLl/vdJBRjYr5akXWdFko6wYqBMzKkO+5wq2/r/PMZvg0ACP5aUXvttZfS09O1cuXKfNvtft26dQt9TEZGhrth92tcxUNZ18yySPtyScPsfEe2vSCp+/z5qnKI5W8AAAhJxqZy5cpq166d/ve//+Vty83Ndfc7dOjga9tSsZso3irMmqUFkh6LBDWfR7I151jGZv/9d/v4VCoETqVjBYDQZmyMDfXu06eP2tvaP0ccoTFjxrjMg42SQjBtXrxY/23fXudKam3DtiUNkTTO5t/zu3EAgEAJXGDTq1cv/frrrxoyZIgrGG7Tpo3efPPNnQqKU6kLaVd/scfum2yyJN0kac9jjnFBjQUx4yNBzW9KLiXpWotmU6ICO7weAAIocIGNGTBggLshuYKokrDupecaNVIFK+beskXvSrpG0qcqXwWDDgBAuASqxga7H01UHLafXxP1HSpphqTn7cNnQU2jRtry1FPqkoCgBgAQfgQ2SIjatsaXpI8lF8RstgkXBw2SFi9Wztln+908AEBIENig1Fmd2JE5RU5+t327+kdW377KZhCW9B9JzSVtu/12qUqV8j0QAEBKIbBB+Xn7bWV17KhHIxmbTyLZml6SflDqYGg2ACROIIuHkdz2k/RmtWo6OxI52winOyQ9zvBtAEA5I2OD+Nm4UXdLWiy5oGa7pE0XX+xW3y7rnDTJtjo5WRgASE4ENig7Gz797LPKattWg23BSkk2N3QbGyI+YoR+97t9AICUQWDjo2TLQpTKxx9LnTpJvXurwk8/aamksyR1lfRFHFbVtseRFQEAFBeBDUoVjO2dlqZtF18stWsnvfeeG920dfBgtZD0UgLbQ5cQACAWgQ1KZts2DYwM36705JM7uqHOP9/NR7Ptllu0RcFX1iwTAMA/BDbIxy7mRV3QK9jw7aOO0hhJe1gxcKtW0syZ0jPPuBmEw4ZsEAAED8O9sVsHSnrAFq087TR3/1dJt0t6cNYsVa1Rw+/mAQCQh4wNimQ5iqrDh2uRpDNs8FN6urZddZWaRuakUbrNIwwAQPIgsMHOPE8XSFpigc3DDytD0jRb32nOHG297z6t9bt9AAAUgcAm4EPDbU2n4q7qXRwVFixQZteu+rekBlZHs99+Ol1SN4t3Dj44bq8DAEB5ILCBU0fSP2xyvc6dlT53rqx8eJAthzBzpl7xu3EAABQTgU2KqyTpeklfSbpEUprnaft556mZpHtth0ybRzi+2SEm3QMAlBdGRcVxWHDQVP7f//SppOaR+x9JOmT6dOUedZR+njTJ59YBAFByZGxSUGNJr9pcNL17u6BmpaSLJR0huaAmGZDVAQCUBoFNCqkW6V6yNZxOtWLgihV1v+SGbz9p9/1uIAAAZURgkwILZKZJuihSR3OLdUFJekPSmnfe0U2S1set1QAA+IvAJuSrfVeYN0+zJT0lqV5kjSfL1pxiQ7mbNInLawAAkCwIbEIqbcUKqW9fZR13nI6UtEHSzZJaSvqv340DAKCcMCoqZKybyVbfzmrTxipw3baJkTlpVvjdOAAAyhkZmxDpYcseHHSQ7rOMjQU1RxyhzTNmqB9BDQAgRRDYhICNanpd0mt2Qr/9VqpbV5o4UZo9W7mHH55vCLXdAAAIKwKbABQAF6WG5IZrfy6pu6StdrvuOmnJEqlPH6kCpxcAkFq48pUjW36gXAKc3Fw3oZ4N374hsiyCTbh3iKRtw4ZJNSzkSc7ZmXc36R4T8wEAyoLAJmA6SKp18sl6QtI+khZLOlnSaZK+SZLgpby7u8otYAQABB6BTUC6qNJ++UX/lPSBZWg++UTrItmaVpKmxiGLUpxsCgAAyY7AJsllWCBz//1u+PafrRfKRj717u0Khkdb15OCr7jdVAAA7A6BTbLyPNe9ZOs6Vb7zTqVt3OiyNW6hyr//XSt30+UTzRBZtw0AAKmCwCYJHWwrb/fqpZclHWSBTL162vLEEzpa0vxSPF9ZCnITVTcDAEA8ENgkk7Vr9aCkTy1L8+67ypb0V+t6+vhj5fTq5XfrAABIegQ2CWDdQbFdQgW7iewkVJwwQVVat9a1kXUusk8+WS0k3b7jCeKWlaGeBQAQZgQ2PrPupXlWJHzNNUr77TctknSipHVPPaXvigiMAABA4QhsfNJQ0jOS3pPU1mqFa9ZU9siRai1pus9tI6sDAAgqApsEy5RUZfRoN7He+ZHh2+MlbVq4UNv799d2vxsIAECAEdgk0JmS62qqNnKkLA8yS1I7SVfYD/feW2FBxgcA4BcCmwRoGeleelHSAZJy6tfXeZI6SVrod+MCgmHnAIDiILApR7UkPRwJXk6QtEXS3ZJ+e+89Ped34wAACCECm3J6U6+IrL59taR0SZMlNZc01HaI6Z5hNWsAAOLHpkxBHFWYNUsLJDe6yXwmaaCkGbt4jM1rE6Th3BaMEYgBAJIRGZt4+f576dxzldW9uwtq1kgaIOmw3QQ1AAAgfsjYxMP27VLnzi648SpU0GO5uRocCW4AAEDiBCpjM3z4cHXs2FFVqlTRHnvsoaRRsaJ0xx0uuNny/vvqT1ADAIAvAhXYbN26VX/605905ZVXKulcfLE0Y4ZyDz3U75YAAJCyAtUVddddd7l/J06cqKRTIVAxIgAAoRSowKY0srOz3S1q/fr1CqqCo5Gik9YBAIAdQp9mGDFihGrWrJl3a9Sokd9NCiSWSQAABIHvgc2tt96qtLS0Xd4WL7YlI0tn0KBBWrduXd5t+fLlcW0/AABIHr53Rd1www3q27fvLvc58MADS/38GRkZ7gYAAMLP98Bm7733djcAAIDABzYl8cMPP2jNmjXu35ycHC1cuGNt7MaNGwdqSQIAAFA+AhXYDBkyRE899VTe/cMOswULbPqYGerSpYvCIHbkkxXqBm0dKQAAUrp4uCRs/proyJzYW1iCGgAAkEKBDQAAwK4Q2AAAgNAgsAEAAKFBYAMAAEKDwKYclx6wm41yAgAAiRGo4d4AC38CAHaFjE0CWfaGRSQBACg/BDYAACA0CGwAAEBoUGMTUtSiAABSERkbAAAQGgQ2AAAgNOiK8qFLqOA2W8EbAACUHYFNkqNWBgCA4qMrCgAAhAaBDQAACA0CGwAAEBoENgAAIDQIbAAAQGgQ2AAAgNBguHcSYEg3AADxQcYGAACEBoENAAAIDQIbAAAQGgQ2AAAgNAhsAABAaBDYAACA0CCwAQAAoUFgAwAAQoPABgAAhAaBDQAACA0CGwAAEBoENgAAIDQIbAAAQGgQ2AAAgNAgsAEAAKFRUSnG8zz37/r16/1uCgAAKKbodTt6HS9KygU2GzZscP82atTI76YAAIBSXMdr1qxZ5M/TvN2FPiGTm5urn3/+WdWrV1daWlpcI0kLlpYvX64aNWoojDjG4Av78RmOMfjCfnypcIzry+H4LFyxoKZ+/fqqUKHoSpqUy9jYm9GwYcNye347gWH8kMbiGIMv7MdnOMbgC/vxpcIx1ojz8e0qUxNF8TAAAAgNAhsAABAaBDZxkpGRoaFDh7p/w4pjDL6wH5/hGIMv7MeXCseY4ePxpVzxMAAACC8yNgAAIDQIbAAAQGgQ2AAAgNAgsAEAAKFBYFMCw4cPV8eOHVWlShXtsccexXqM1WYPGTJE9erVU1ZWlrp27aqvv/463z5r1qzRBRdc4CYxsue95JJL9McffyjRStqOZcuWudmbC7s9//zzefsV9vNJkybJD6V5r7t06bJT+6+44op8+/zwww/q0aOH+2zUqVNHN910k7Zv364gHKPtf/XVV6tZs2buM7rvvvvqmmuu0bp16/Lt59d5HDt2rPbff39lZmbqyCOP1IcffrjL/e2z17x5c7f/oYceqtdff73Ev5OJVpJjfPzxx3XssceqVq1a7mbtL7h/3759dzpXJ598soJyjBMnTtyp/fa4ZD6PJTm+wr5T7GbfIcl6DmfOnKmePXu6WX+tLVOmTNntY9555x21bdvWjYxq3LixO69l/f0uFhsVheIZMmSIN3r0aO/666/3atasWazH3HvvvW7fKVOmeJ988ol32mmneQcccIC3efPmvH1OPvlkr3Xr1t6cOXO8WbNmeY0bN/bOP/98L9FK2o7t27d7v/zyS77bXXfd5VWrVs3bsGFD3n72MXvyySfz7Rd7/IlUmve6c+fO3mWXXZav/evWrcv3PrRs2dLr2rWr9/HHH3uvv/66t9dee3mDBg3ygnCMn332mXfWWWd5r7zyivfNN994//vf/7wmTZp4Z599dr79/DiPkyZN8ipXruxNmDDB++KLL9x52GOPPbyVK1cWuv/777/vpaene/fdd5+3aNEi74477vAqVarkjrEkv5OJVNJj7N27tzd27Fj3Wfvyyy+9vn37uuP58ccf8/bp06eP+xzEnqs1a9Z4finpMdrnrEaNGvnav2LFinz7JNN5LOnx/fbbb/mO7fPPP3efWzvuZD2Hr7/+unf77bd7L774ovsueOmll3a5/3fffedVqVLFXS/td/GRRx5xx/jmm2+W+n0rLgKbUrAPX3ECm9zcXK9u3breqFGj8ratXbvWy8jI8J599ll33064fUjmzZuXt88bb7zhpaWleT/99JOXKPFqR5s2bbyLL74437bi/BIk8zFaYDNw4MBd/sJXqFAh3xfvY4895r6Ys7OzvUSK13n8z3/+475wtm3b5ut5POKII7z+/fvn3c/JyfHq16/vjRgxotD9zz33XK9Hjx75th155JHeX/7yl2L/TiZaSY+xIAusq1ev7j311FP5Loqnn366lyxKeoy7+45NtvNY1nP44IMPunP4xx9/JO05jFWc74Kbb77ZO+SQQ/Jt69Wrl9etW7e4vW9FoSuqHC1dulQrVqxwKdLYdS4s3TZ79mx33/617oL27dvn7WP725pWc+fOTVhb49GO+fPna+HCha7ro6D+/ftrr7320hFHHKEJEybsdtn5ZDvGp59+2rW/ZcuWGjRokDZt2pTvea3LY5999snb1q1bN7cI3BdffKFEitfnybqhrCurYsWKvp3HrVu3us9U7O+PHYfdj/7+FGTbY/ePnovo/sX5nUyk0hxjQfZZ3LZtm2rXrr1TN4B1i1oX45VXXqnffvtNfijtMVr36X777ecWUjz99NPz/S4l03mMxzl84okndN5556lq1apJeQ5LY3e/i/F434qScotgJpL94pnYC170fvRn9q99cGPZxcS+pKL7JKqtZW2H/XIefPDBrg4p1t13363jjz/e1Z9MmzZNV111lfvSsjqORCrtMfbu3dt9wVrf8qeffqpbbrlFS5Ys0Ysvvpj3vIWd4+jPEike53H16tUaNmyYLr/8cl/Po7UjJyen0Pd28eLFhT6mqHMR+/sW3VbUPolUmmMsyD6P9tmMvUBYLcZZZ52lAw44QN9++61uu+02de/e3V0w0tPTlezHaBdyC5xbtWrlguz777/ffa9YcGOLGCfTeSzrObSaks8//9x9f8ZKpnNYGkX9LtoffJs3b9bvv/9e5s9+UVI+sLn11ls1cuTIXe7z5ZdfumLEMB9fWdkH9ZlnntHgwYN3+lnstsMOO0wbN27UqFGj4nZBLO9jjL3AW2bGihVPOOEE92Vz0EEHKUzn0b50rICxRYsWuvPOOxN6HlFy9957ryvgtr/sY4tr7a//2M+sBQj2WbX97LOb7Dp06OBuURbU2B9N48ePd0F3mFhAY+fIsqCxgn4O/ZTygc0NN9zgqs935cADDyzVc9etW9f9u3LlSncxjLL7bdq0ydtn1apV+R5no2lspEr08Yk4vrK2Y/LkyS4lftFFF+12X0sX25dTdnZ2XNYRSdQxxrbffPPNN+6Lxh5bsJLfzrGJxzlM1DFu2LDB/ZVYvXp1vfTSS6pUqVJCz2NB1uVlf5lG38sou1/Usdj2Xe1fnN/JRCrNMUZZFsMCm+nTp7uL3u4+G/Za9plN9EWxLMcYZZ9FC6at/cl2HstyfPbHgQWmlg3dHT/PYWkU9btoXdw2is3es7J+LopUpgqdFFXS4uH7778/b5uNpimsePijjz7K22fq1Km+FQ+Xth1WYFtwFE1R7rnnHq9WrVpeosXrvX7vvffc89hIjNji4dhK/vHjx7vi4S1btnhBOEb7XB511FHuPG7cuDFpzqMVFw4YMCBfcWGDBg12WTx86qmn5tvWoUOHnYqHd/U7mWglPUYzcuRI9/maPXt2sV5j+fLl7jPw8ssve0E5xoIF0s2aNfOuu+66pDyPpT0+u5ZYm1evXp3057A0xcM2WjSWjc4sWDxcls9FUQhsSuD77793QyyjQ5rt/+0WO7TZfvlsOFzskEQbvmYfxk8//dRVuRc23Puwww7z5s6d6y6aNtTWr+Heu2qHDSe147Ofx/r666/dL5yNvinIhhA//vjjbrit7fe3v/3NDQG0ofN+KOkx2vDnu+++2wUKS5cudefxwAMP9Dp16rTTcO+TTjrJW7hwoRvOuPfee/s63Lskx2gXBBs5dOihh7rjjR1easfm53m04aD2xT9x4kQXtF1++eXu9yk6Au3Pf/6zd+utt+Yb7l2xYkV3wbOh0EOHDi10uPfuficTqaTHaO23EWuTJ0/Od66i30P274033uiCHvvMTp8+3Wvbtq37HCQ60C7tMdp3rAXk3377rTd//nzvvPPO8zIzM92Q4GQ8jyU9vqhjjjnGjRQqKBnP4YYNG/KueRbY2NQn9v92XTR2fHacBYd733TTTe530aYoKGy4967et9IisCkBG35nJ7TgbcaMGTvN9RFlf1kMHjzY22effdwJPOGEE7wlS5bsNKeBXXgsWLK/wvr165cvWEqU3bXDfsEKHq+xC3ijRo1ctF2QBTs2BNyes2rVqm5+lXHjxhW6bzIe4w8//OCCmNq1a7vzZ3PC2C9q7Dw2ZtmyZV737t29rKwsN4fNDTfckG+odDIfo/1b2Ofabrav3+fR5r/Yd9993cXc/sKz+XmiLMNkv5cFh6o3bdrU7W/DTf/73//m+3lxficTrSTHuN9++xV6riyIM5s2bXJBtgXXFtTZ/jY/SFkvFok8xmuvvTZvXztPp5xyirdgwYKkPo8l/ZwuXrzYnbdp06bt9FzJeA5nFPE9ET0u+9eOs+Bj7HvD3hP7gzD22lic96200uw/ZevMAgAASA7MYwMAAEKDwAYAAIQGgQ0AAAgNAhsAABAaBDYAACA0CGwAAEBoENgAAIDQILABAAChQWADAABCg8AGAACEBoENgEB79tlnlZWVpV9++SVvW79+/dSqVSutW7fO17YBSDzWigIQaPYV1qZNG3Xq1EmPPPKIhg4dqgkTJmjOnDlq0KCB380DkGAVE/2CABBPaWlpGj58uM455xzVrVvXBTezZs0iqAFSFBkbAKHQtm1bffHFF5o2bZo6d+7sd3MA+IQaGwCB9+abb2rx4sXKycnRPvvs43dzAPiIjA2AQFuwYIG6dOmi8ePHa+LEiapRo4aef/55v5sFwCfU2AAIrGXLlqlHjx667bbbdP755+vAAw9Uhw4dXLBjXVMAUg8ZGwCBtGbNGnXs2NFla8aNG5e33QId65Ky7ikAqYfABgAAhAbFwwAAIDQIbAAAQGgQ2AAAgNAgsAEAAKFBYAMAAEKDwAYAAIQGgQ0AAAgNAhsAABAaBDYAACA0CGwAAEBoENgAAIDQILABAAAKi/8D8mKwNkgY94sAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.errorbar(x, y, ye, ls='', c='k')\n",
    "ax.plot(x, model.evaluate(x=x), c='r')\n",
    "ax.set(xlabel='$x$', ylabel='$y$', title='Randomly Generated Data');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f32954cb-f7b4-4ee0-a523-f9ff57b897b1",
   "metadata": {},
   "source": [
    "This is a good fit. To see the results of the fit we must first introduce a new class: `Distribution`s."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c8471f8-50a5-411c-b32b-995a169721aa",
   "metadata": {},
   "source": [
    "There are two occasions inside `lamatrix` where we require some concept of a distribution.\n",
    "\n",
    "1. The **prior** distribution: this holds our prior belief of the values of best fitting weights. By default, all priors are infintely wide Gaussians (i.e. we have no prior understanding).\n",
    "2. The **posterior** distribution: this holds the best fitting values after we have fit our model to data, informed by the prior.\n",
    "\n",
    "Both of these are Gaussian distributions in `lamatrix`. We use the following Python objects to hold these distributions.\n",
    "\n",
    "## `Distribution`\n",
    "\n",
    "In `lamatrix` all distributions are by definition Gaussian. These are specified by two parameters: the **mean** and the **standard deviation**. These are specified as `Tuple` objects with `(mean, standard deviation)`, e.g.\n",
    "```\n",
    "(0, 1)\n",
    "```\n",
    "\n",
    "Would indicate a Gaussian with mean 0 and standard deviation 1.\n",
    "\n",
    "Distributions have special properties compared with regular tuples:\n",
    "\n",
    "* `mean` and `std` properties return the mean and standard deviation as floats respectively\n",
    "* freezing/thawing: In `lamatrix` you may want to \"freeze\" a certain variable to its mean value, i.e. make the standard deviation 0. `Distribution` classes have the `freeze` and `thaw` methods to enable this.\n",
    "* you can draw a sample from the distribution using the `sample` method\n",
    "\n",
    "## `DistributionContainer`\n",
    "\n",
    "One `Tuple` `Distribution` object is valid for each weight we fit in the model. We must hold many weights to account for all of the components in the model. These are held in a `DistributionContainer`. This is a special case of a `List` which has the following methods\n",
    "\n",
    "* `mean` and `std` properties return the mean and standard deviation as `numpy.NDArray`'s respectively\n",
    "* `freeze` and `thaw` methods as above\n",
    "* `sample` method as above"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc8a983e-5dc9-4cb1-b6aa-e93a20ca67a8",
   "metadata": {},
   "source": [
    "Now that we understand distibutions, we can return to the model. To find the best fit of our model, we need to access the `posteriors`. These are the results of the fit. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "100c6159-79ad-472b-b22d-4e84387c0360",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DistributionContainer\n",
       "\t[(2.36, 0.12), (1.445, 0.071)]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.posteriors"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6379d0f8-7237-4e6b-9032-222a500c70a3",
   "metadata": {},
   "source": [
    "These posteriors represent the mean and standard deviation of the best fit weights for our model. We can look at the means"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "15dcd6fa-7f06-4f8a-99d1-17cfbf3ee713",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.36184521, 1.44454686])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.posteriors.mean"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab842980-1711-49cd-b0f5-f9fe64dead14",
   "metadata": {},
   "source": [
    "These values match the input $m$ and $c$ very closely! These represent our best fitting weights. \n",
    "\n",
    "Let's improve our fit by including errors on the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "f2551827-bf62-4f1a-8b18-bcfa85867d61",
   "metadata": {},
   "outputs": [],
   "source": [
    "model.fit(x=x, data=y, errors=ye)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "454fef0b-406b-4cd4-9240-caab690e9508",
   "metadata": {},
   "source": [
    "Now we have fit again, let's look at the mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "a11c3d06-22f6-4a61-bfd6-871776c61e00",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.36184521, 1.44454686])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.posteriors.mean"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cb0bf81-e27b-419c-b680-0008f121eeaa",
   "metadata": {},
   "source": [
    "The posterior standard deviation now includes a more reasonable estimate of the error on these best fit parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "acd25b3d-48b6-4d3b-846b-c2bc2f58fa42",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.03674281, 0.021214  ])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model.posteriors.std"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab2a7014-cc7b-4d49-a386-94dcf5a5bdae",
   "metadata": {},
   "source": [
    "We can plot this new model fit over the data. Here we draw 100 samples of the model, which draws weights from the posterior distribution, rather than using the best fit mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "7ab78f54-627b-4d7a-968c-05bd7159ffaa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjYAAAHHCAYAAACskBIUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAABjXElEQVR4nO3dB3hb5dk+8Ft2lrOcvfdehCwyySAJZLJTNiRAgZbRUiiU9CtQCjS08OdrP8oIFAgtBEqAEKAQMsiC7L23s/dyQuw4jn3+1/3qPUaW5cRb0tH9uy5lSLJ8jmRbj9/3GT7HcRyIiIiIeEBcuA9AREREpKgosBERERHPUGAjIiIinqHARkRERDxDgY2IiIh4hgIbERER8QwFNiIiIuIZCmxERETEMxTYiIiIiGcosBGJMmPGjEGTJk1i7nPL+U2YMAE+nw87duwI96GIhJUCG5ELvFG4l1KlSqF+/frmzX3v3r3hPjzPyMzMxL/+9S9cfvnlqFGjBkqXLo1atWrhiiuuwJtvvom0tDR4RUpKCv74xz9i9uzZYTsGfv7Ar+vy5cujUaNGuPLKK/Huu+8W6vn++uuvzeOLhFOpsH52kSjwpz/9CU2bNsWZM2ewcOFCE/B8//33WLt2LcqVKxfuw4tqqampuPbaa/Htt9+id+/e+O1vf4vatWvj2LFjmDNnDu6//34sWrQIb7/9NrwS2DzzzDPm3wMGDAjrsbz++uuoWLGiCWQYqPM1uOuuu/C3v/0NX331FRo2bFigwObVV19VcCNhpcBG5AKGDRuGbt26mX///Oc/N6sKf/nLX/DFF1/ghhtuCPfhRbXf/OY35g2Vb6a//vWvs9326KOPYsuWLZg+fToi1blz58yKU5kyZRBtRo0aZb6WXU899RQ++OAD3HHHHfjZz35mgniRaKStKJF86tu3r/l727ZtWdedPXvWvDF07doViYmJqFChgrnfrFmzsn0s8x+4/P/SSy+ZbZbmzZujbNmyuOSSS7BkyZIcn+vzzz9Hhw4dzMoQ/548eXLIYzp9+rQJBPhbNh+vdevW5nM4jpPtfvzcDz74ICZNmoR27dohISEBvXr1wpo1a8zt48ePR4sWLczn44rC+fI1+NjMt7n66qtz3MbVLT4P9913X64fv3v3bvzzn//E0KFDcwQ1rpYtW5pVm0AMJBgItW/f3hwnV3j4eY4fP57tfjy2kSNHmtW17t27m/s2a9bMbHsFO3HiBB5++OGs54/PAYNXfq5Qrx0/v/varV+/Pk+vPz++Zs2a5t9ctXG3ggJXNzZu3GgCjmrVqpnjZUDNADrYunXrMHDgQPP6NWjQAM8991y2Yy2oW2+91QTvXCULDCjnzZtngh1uWfGc+TwxKOWKm4tbtFytocCtLhefN67KVa9e3Rw3n6tPPvmk0McsEkwrNiL55L7ZV61aNeu6kydPmjfpm2++Gffccw9OnTpltk+GDBmCxYsXo1OnTtkeY+LEieY+fEPmD/+//vWvuO6667B9+3aTY0LTpk3D9ddfbwKQcePG4ejRo7jzzjvNG1lwgHHVVVeZN9G7777bfC6ugjz22GNmi+F///d/s92fb1J8s3zggQfM//nYDAAef/xxvPbaayaQYJDAY+LWxHfffRfyeeBx33bbbeZ+3Drim7Hryy+/NM8Jb8/NN998g4yMjPPeJxQ+Z9wO5HPxq1/9CklJSfjHP/6BFStW4Icffsh6/mjr1q0mUODzMnr0aLzzzjvmDZhvqgyM3O2h/v37m+eKj8037/nz52Ps2LHYv3+/CWICMQ+Fgdu9995r3uR53nl5/RnUcPvnl7/8pdl+4+tNHTt2zApW+vTpY/K4nnjiCRMcffzxx7jmmmvw6aefmo+hAwcO4LLLLjOrRe79GCQzWCgKt99+u3k8fv0x74kYCPN54rEzMOE5vfLKK9izZ4+5zX1d9u3bZwKif//73zke9+9//7v5OmXwxEDwo48+MsESt71GjBhRJMcuYjgiEtK7777L5Q5nxowZzuHDh53du3c7n3zyiVOzZk2nbNmy5v+uc+fOOWlpadk+/vjx407t2rWdu+66K+u6pKQk85jVq1d3jh07lnX9lClTzPVffvll1nWdOnVy6tat65w4cSLrumnTppn7NW7cOOu6zz//3Fz33HPPZfv8o0aNcnw+n7N169as63g/HjuPwzV+/HhzfZ06dZyTJ09mXT927FhzfeB9R48ene1zb9q0ydzn9ddfz/a5r7rqKqdJkyZOZmZmrs/vb37zG/OxK1euzHY9n0c+3+7lyJEjWbfNmzfPfMwHH3yQ7WOmTp2a43oeJ6+bO3du1nWHDh0y5//oo49mXffss886FSpUcDZv3pztMZ944gknPj7e2bVrV7bXrnLlyuZxAuX19ef58DGefvrpHM/HoEGDnIsuusg5c+ZM1nV8/nr37u20bNky67qHH37YPMaiRYuynVdiYmKO1ysUfm7ej8cSCo+bt1977bVZ16WkpOS437hx48zX186dO7Oue+CBB8zHhhL8GGfPnnU6dOjgDBw48LzHK5Jf2ooSuYDBgweb37a5/M7f/vkbMlc8AldO4uPjs/IsuCXAFQz+Rs2thOXLl+d4zBtvvDHbio+7vcUVG+JKwcqVK80qA7c2XPwNmis4wQmb/PxcvQjErSnGMlwZCTRo0KBsJds9evQwf3N1qFKlSjmud48plFatWpn7MTfDxXPn5+Rv5oFbEcG4ykFMYA0+Hz7f7qVx48ZZt3F1gM8Hn4cjR45kXbgCw8cJ3vrjc+U+t8TH4zZd4DnxMXkfvh6Bj8nXnStKc+fOzfaYfJ7cLaWCvv7BeH+ujDFni6s97jFwlY6rPsw1civx+Pz07NnTbK8Fnhef76Lgvh48DlfgahC3PXls3Fbi1xdXyvIi8DG4IpicnGye97w8PyL5oa0okQtg3gDfwPmDmFsZfKPjFkSw9957D//v//0/kyeRnp6edT0rqoJxuyOQG+S4eSI7d+7MyjEJxjfmwDcD3rdevXrZghJq27ZttsfK7XO7gVNwFYx7fXDuSjAmmzJvh5+HQQgDBZ4/tzTOxz3eH3/8Mdv13I5x8ztefPFFs73k4hs8XweWg4dy6NCh856r+1wHnhMfc/Xq1TmCldweM9Trmd/XPxi3zBgkPPnkk+aS23Fwm4rPsxt0Bn9dFAX39Qj8etq1a5fJIWJAH/z1wNcjL7jlxFwgBuyBJeXnC35FCkKBjcgF8DdjtyqK+Q6XXnopbrnlFmzatCnrt9v333/f5G7wdua28I2Xv8UzfyUwydjF20IJTvYtDrl97oIe00033WQSSblq8/vf/948F3y+LvRG26ZNG/M3y+YvvvjirOsZYHC1hPhYgbgawuc2cIUoUKiVlAudEx+TK0DMMQqFQW2gULks+X39g7mJvyx35wpNKExoLgl8PQI/H1et+PxwVel3v/uded24askVJJ5zXpKWmdfF/Jp+/fqZPK66deuaXCjmKzHfTKQoKbARyQf3zYrJm0xYZfImsbqDFTefffZZtt9An3766QJ9Hnf7hasJwRhQBd93xowZZusg8LdsrhwEPlZxYfIskz8ZbHA7hCsswQm3uZXR8/l0Py4vWInEc+WqTlEly/IxuUrhBlMFkdfXP7fVCX4s8c3+QsfB1zMvXxcF5Sb+ugEWK+Y2b95sVqS4OucKVYaf2/kx+ZlVXkxqD1ztZGAjUtSUYyOSTyyD5ioO37xZHRO4MhC4EsCS2QULFhToc/A3WlbS8M0kcKmfbyYsLw40fPhw81s1A61ArIbiGw0DiOLGbSceF1cr+FxwFedCuE3Eqivm4wQfe26rRcxB4bk+++yzOe7LnBaWbecXH5OvE990g/Hx+LgXktfXn11+3ccNxBUefl2x3J75VcEOHz6c7fVmjxlWJgXentsqVn5w9YTVXWwBwFys3M6N/2aVUzCu5IQ6Pz4Gvxb52gVWF7KdgUhR04qNSAHwDZylqiw7/sUvfmHKpfnbOktyuXrBEuQ33njDJK8G55DkFVeG+Fjc+mIAwK0AltiyTDnwMdkKnytI//M//2PeLLitw1LdKVOmmN4sXJEobjxOlgEzv4aBVG45MMEYHPK5euihh0z5L8+FH8vkVK78sGw8cEuLZdksK+Zzw1wNjl3gKgdXMPi5+WbLBO/8vpbMHeFr6JaCM0GWKxVcieFzGtjILpS8vv5cZeJ1//nPf8wWF1e72J+IF+Zy8bW+6KKLTMk4V3EOHjxogiOWVa9atco8BrfMuKri9v9xy725ksNcobziuXErlaXXbudhPuf8+nFLuIlbT/wa4jYZ71e5cmWzAhMq94rPHTGRnSs+bpDL5+Tll182x8xtXOYL8Xy53ZWfYxbJk3zXUYnEWLn3kiVLctyWkZHhNG/e3FxY6suy3D//+c+mxJjlxJ07d3a++uqrHOXRbsnwiy++mOMxQ5UBf/rpp07btm3NY7Zr18757LPPcjwmnTp1ypRP16tXzyldurQpD+bnCC635udgSW6g3I5p1qxZ5vpJkyZlXRfqc7vuv/9+c/+JEyc6+cHnj881y36rVavmlCpVyqlRo4Ypf37jjTec1NTUHB/z5ptvOl27dnUSEhKcSpUqmTLpxx9/3Nm3b1/WfXicI0aMyPGx/fv3N5fg54/l7S1atHDKlCljPj/LrF966SVTlny+54ny+vrT/PnzzbHz8wS/5tu2bXPuuOMOU3rP17F+/frOyJEjTZuBQKtXrzbnUK5cOXMflqy//fbb+Sr3di98jAYNGpjP884772QrN3etX7/eGTx4sFOxYkXz3Nxzzz3OqlWrzMfztQt8LR966CHTEoGl4IFvMTw+fl3y+WnTpo35OPdYRIqSj3/kLQQSEckdE4jZlI4N5NwtFxGRkqYcGxEpNOYasTKIPV4U1IhIOCnHRkQKjLkSrFJivgabyeU280lEpKQosBGRAmMlFEu1mfD7f//3fzlmYomIlDTl2IiIiIhnKMdGREREPEOBjYiIiHhGzOXYcK7Jvn37TOt5DV8TERGJDsyc4egYDv2Ni8t9XSbmAhsGNcFTjEVERCQ67N69Gw0aNMj19pgLbNwhgXxi2BpcREREIt/JkyfNwkTgsN9QYi6wcbefGNQosBEREYkuF0ojUfKwiIiIeIYCGxEREfEMBTYiIiLiGQpsRERExDMU2IiIiIhnKLARERERz1BgIyIiIp6hwEZEREQ8Q4GNiIiIeIYCGxEREfEMBTYiIiLiGQpsRERExDMU2IiIiIhnKLARERERz1BgIyIiIvl2+vRp+Hw+c+G/I4UCGxEREfEMBTYiIiLiGQpsRERExDMU2IiIiIhnKLARERERz1BgIyIiIp6hwEZEREQ8Q4GNiIiIFLp3TaT0tVFgIyIiIp6hwEZEREQ8Q4GNiIiIeIYCGxEREfEMBTYiIiLiGQpsRERExDMU2IiIiIhnKLARERERzygV7gMQERGR6FaxYkVECq3YiIiIiGcosBEREZGQImVMQn4osBERERHPUGAjIiIinqHARkREJIpE4/ZQSVJgIyIiIp6hwEZEREQ8Q4GNiIiIeEZUBzYvvPCC2WN8+OGHw30oIiIiEgGiNrBZsmQJxo8fj44dO4b7UERERCRCRGVg8+OPP+LWW2/FW2+9hapVq4b7cERERCRCRGVg88ADD2DEiBEYPHjwBe+blpaGkydPZruIiIh40ekClIJ7rXw86gKbjz76CMuXL8e4cePydH/eLzExMevSsGHDYj9GERGRWAkkIk1UBTa7d+/Gr3/9a3zwwQcoV65cnj5m7NixSE5OzrrwMURERGJh4rbPBlCxFESVQhRZtmwZDh06hC5dumRdl5GRgblz5+If//iH2XaKj4/P9jFly5Y1FxERESlcoBQNoiqwGTRoENasWZPtujvvvBNt2rTB7373uxxBjYiIiBSPChUqwHEcsxIUSUFPVAU2lSpVQocOHXI8sdWrV89xvYiIiMSeqMqxERER8ZpITCY+bY+pICsx4U7+iKoVm1Bmz54d7kMQERGJeWUBtAJgao9TUrilEpbjiPrARkRERMIoIwONALQHwKSQFrxu+3agZs2wHI4CGxERkRhQ0W4rHTx4ELVr187q5F9gjgMcOQLfkiUYCoD1yo3crah9+xAuCmxEREQkX0wnuVWrgM2b4Vu0CMMBcMBRLQZQvC2MXf4V2IiIiMSoivlMDo6zOTTcdor78ktg0yZgzx40tkFNVuvcH34ARo9GOCiwERERyYPAfi3cwmG7kVhSFUBHm0fTiVfMmwccPw4cOID6NqAo5/MhgT3lGnFTKjwU2IiIiEiuuArT1gY0DGweGDgQCceOAVu3+recMjJQBkAm7xwXB5QqBfh8CBcFNiIiInlYpWHSbSyJt8nA3HbqaoMbk3LM7afkZCAtDcjMNIHNOeYS87aMDH9ScWGSkgtJgY2IiIj8hIHJoUPoa7ecOtgAh1tRdXj74cP+gIYX3jc+HvEZGTBrNFyp4YpNairCRYGNiIiI+LGx3vr1iFuwANcAaAqgpk0MrgSgNO+Tnu6/L7ed7IqNz12xKVcOqFMHaM91nvBQYCMiIhLj4gFT2eT7+mtg7Vpz6WyDmuo2oMkxg4nbTjaXJt1xcBxAg9atgYsuArpy8yo8FNiIiIjEsGo2KfgSBjZffOFvrnfgAJrbxGEGCgxf0u2qTBluPwVuO1WogG0nTmAzgFEDBwIMbqozHAoPBTYiIiIxqDyANjag6WC3nbgNxfwaJv9WsCs5TAzOCAhwDJZ0s9y9Rg2gYUN8N2cOdjOw6dMHaNdOgY2IiIiUkPR0sxrT2Y5BaGVzaEzFU1KSP2/mHMMZ/ypNnA0W4tySbubRMKBp3NgfxDRogEVz5mAbb7/4YqBePSAhIWynp8BGREQkSpv55fdz1OAb/6xZuAFAOwD17HW8mFDk7Fl/pZPddooPWKVhqHMMQP1WrfwBTdu2QJMmSGjZEp/fdRdQpYp/FSeMPWxIgY2IiIjHA6Pyth8Nt53KfvklLg1YpXG3nLLCEQYmjpN1HQMaTn46BJhVmfqDBwMMbnjhqg1Xb3isrJKKAApsREREomieU6hGgbk1ESwFoIkNaLjtxC2ouNWrTW4NP4Idg+NtDg0DmNJuYrBdsWGbvSPcoWJfPvv31X37Am3aALVqAZUr+xOII0hkHY2IiIgUms9uLzGY6QGgtd124ioNduwwlVDu+kpgHo3/g31A2bLYf+YM9gHodccdaNmiBd5+6il/Hk2nTv48mjIMiyKPAhsREREPqWCrnHraPBp38jZ70pgNrTNnsqqdAgMBE+iULu3PlalZE9+vX48tDGz69YPTogWWAqZXDW+L1KCGFNiIiIh4gKlc2rgR19lRCNx2qmsDmkS75eSu0px1uwgH5NEkc5WHOTM2f8YNbJxeveDUqIH9iA6RkekjIiIiecJEYsdxTFKxG5gweBnKxOD//AfDbE5NG9ubhjOeygQ02YMNatwxCKcA7AGwmjcMGgQMHQpcdhnmAVjD+zRs6K92ihJasREREYlSFWzX4F5226nU4sVmG6q6Ld8ulUseDYOaNFu+zaBmPYCtAAayc3DbtshMTMQOt29NhFQ75ZUCGxERkWjDfjNr1+ImG9g0tdtO2LgR9QO2nc7ZIMYt3TYl3fHxOJKRYRKDGcxsAMyW03YGMpzxxBWa9HR/UBOFFNiIiIjnG9l5hQlMdu8GVq5E3A8/YLBNDuZWVBXelpJitp0y7CoNt5wyA/JoTnM2VKNGWJaUlBXQMLjhOITD/HhO5mZisDvBOwopsBEREYkClWxScNy//w1s3w5s3mxWa6oGbjtlZobcdmJAcxTALgB9+vfHt0lJJqDZyd43jIfgHQpsREREIqT5XihcdWkJoL/No8Hcuf4J3AcPmt40cQErMvEBwyp5OROQR7PBNtnrOXAgvp8wwWxFMXH4p5Z83qDARkREJAIxMKkDoB+ArnbLqbE7gZsVUWfPmjdxtx+Nm+Lrs9tQHIOw3243rQWw2ebRnOvSxfw7WnNoLkSBjYiISKQ5cQIDADPTqblNDmaTvWr2NjOB2wYwpQIGVrqrQKk1a2J2UhLW2YqnJLsNxTwap2bNHEGNl/KbFNiIiIiEoQ9NYOJ0lrQ0U9nkmzoV1wc02WP5dkU3eLFBDQKqn9y5TpXr1jUTt1my/d+kJLNa8/Xq1UDduqjIjsExQIGNiIhIkJBBRxE+pttczxVnAxjfRx+ZrSbf5s3obqdvV7a3u+XbgW/c8TaPhqMOmDNzCQdUdu+OhLZt8c9HH/WXblesiNMpKUUSjEUDBTYiIiLFKLfJ2y4GLr1tkz3f9OlAUhKwZ4/ZfiptgxcnIDk4S1wcjmRm4pDNo1nJwObqq81qDRo1AqpWjbrmekVBgY2IiHhGVPXFOXPGTN9mLk0LW/mERYuA48dNYnCCDWgybEJwXGBgU748UL06Vu3ebUYh8GI6BXfvDjRtaprwFea5yy0IiwYKbEREJCZEStBj1lCSkuD77jvcCqAZgCFNmyLhyBGAF9sc71xAMMO/M20/mopsosctpvbtMXnCBNOPZpfdigJzbAoQ1HiJAhsREZFizJ/J4jhm24n9aOL++U9g2zaz/cRqJ+zd6x+TECBwDALzaE7Y8u0uffoAvXsDzZubfjR7AaQW8flUiKKcmmAKbERERPIpvys/ZRnMzJ+PX9rSbXz/vcmjaW5vc4OaVLv15JZvu/1oDtoeNEsY2IwaBbRpg8yaNc1qjWSnwEZERKQIBVZTcQupEYArGMB8/jn62P+zpBupqWY6dzB32+mk7TuTZBODV9t//7FHD6BxY/PxBV2BcXIrN/cABTYiIiLFgEMpBwLoBqAJ33B/+AEXAUjkjaeZLZM9MdiVYuc67bXBzDIb0Ox082hq1YrJaqe8UmAjIiJRL5JWH7i1xB40fe0IhNaAmenErSd2Dg5M7XVLuQ0GK5UrY8uJE6Zb8FIu7NiAhkFO9gycwqsQtIrD59ALFNiIiIgUgTi3wsmWbre0206cvm1CrrS0rCZ7ZQI+hnk08QzKWNHUtCk+mTbNzHbaY4Mab4QbJUeBjYiISGFGIXDV4+BB3BAwrLKp7Rpc0a7guEMQfEErNqdsYnCLjh2Bnj2Biy/Gd9OmmfJtTuWW/FNgIyIiUlDcvpk3D3GzZmEwgDZ2NEIVuyoTZ/NoEBDYGKVKmc7AGw4fxnIGNnfcYfrSZNaqZRKFpeAU2IiIiJxHqJLueBvE+N58E9i508x36m0DmgoBqzLubCd3xSbdznWq3by5Kdn+YMoUbAZwL/vStG6d1ZxPCk5p1SIikm/cjvH5fObilaTTPOG20+7duA3AGAY2c+cCM2YAq1aZBOHKNqjJDGqyZ1J0ExJMdRN70eCmm4Dbb8dcAIv5/2bNgDJu5o0UhlZsRERE8oIdhWfORNy8eWa+Uxt3ttOZM6bBXmkb0Ljdgt2Ahn2IjzD3hiMQli7FCgDDBw8GmjTBpnCfkwcpsBERETkfbg8tXAjMmgVs3QqsXm22nVjthJNso+eXabdBHPt3mg1o9tteNGPGjMG3S5eaVRt06gT4sjJupAgpsBEREU8rTI+b+lx5+de/zFwnbNhgcmlw7Bjq5PIGylAl3c51YrXTBsBsN3H0we09ephtKLNNpaCm2ERVjs3rr7+Ojh07onLlyubSq1cvfPPNN+E+LBER8ZoTJ3Alk3oZg8ycCUyb5t92OnjQjDIoFVTtBLtSw8TgXXaFZiKA921gs4gBTaNGWbk3UnyiasWmQYMGeOGFF9CyZUvTS+C9997D1VdfjRUrVqB9+/bhPjwREYlyp48fx+Bq1cwE7qEAzDvLvHn+su6A6dvudlOW0qWxJz3ddAj+3nYNZnO9HTbHpkT76sS4qApsrryS8fNPnn/+ebOKs3DhQgU2IiKSjftmf5CrLHmpdtq8Gb4vv8TPALQCzORtjkDAca7D5GQ2k+LjgSpVgNq1MW39ehPU7LCTuHedOmW2nBR0lKyoCmwCZWRkYNKkSSZS5ZZUbtLS0szFdTIg0UtERARHjvi3mlasgG/DBrNaw67BlbgQk8uHcBuK7ybVmzY13YLRty8mPfywSQw+UER5NKH65+Rl5lOsi7rAZs2aNSaQOXPmjImCJ0+ejHbt2uV6/3HjxuGZZ54p0WMUESlugdsPeX0DjAYlel4s02Yfmh9+8CcHr15tetRwHEK5wP4zwRISsDs11TTWG3r77UCvXsisWdP0o1F7vfCLusCmdevWWLlyJZKTk/HJJ59g9OjRmDNnTq7BzdixY/HII49kW7Fp2LBhCR6xiIgURHEFOGYdZc0aYP58f/k2/71li0kY5qDKBBvQZAT0ozFKlwaqVwcaNsSkJUtMIHPF0KFA27ZmMreCmsgQdYFNmTJl0KJFC/Pvrl27YsmSJfj73/+O8ePHh7x/2bJlzUVERGJvBShYdQYjTPydOBFISvIHNdyKCuie7PajiXODGm4pVa3qn77dowdw6aX4cskSk0cD/lLN5y2Wui9HuKgLbIJlZmZmy6EREZHoDF6KM8mWqzCXAugOoBGvYE7N/v3+gCQjuHDbL8NO367apAnQpQtwxRVmdSazTh0zqJIBUMVKlfKeoCwlIqoCG24rDRs2DI0aNcKpU6cwceJEzJ49G99++224D01ERCIQZzXFrVpl+tE0seXb/Nvk1JzjiMrQuP7CUGULgCFjxpg8GjPPiakM6enqRxPBoiqwOXToEO644w7s378fiYmJplkfg5rLL7883IcmIiJhcL4eLnUZlDCF4dNPTaUTS7hr2tWbXIOa+Hjszcgw/Whm2cZ6Vwwf7s+jcVe4NIE7okVVYPP222+H+xBERCTCMfxgINONQ7OZ8zt9Oi6x5du5zs9mHo3tR/Ptxo0mqNlue9I4bdr8FNRIxIuqwEZERLyVv1OU3XP5htYRMJO3G9t/82/s3Gka7QWWb2frMFO+vD8xuGNHYNAgvP/gg6YfzVHNdYpKCmxERCSqk5Z9NiF4MICWANqyNYjtGlyed8jICB3UsGKW5dtckWFicI8eph/NihLqR6PGesVDgY2IiESvw4dxFYDOdtupo82jqZhL12AGNayjLcuAhl2DBw0C+vcHWPnUuLEZraAMmuimwEZERKJPaiqwcCHi5szB5TawaWDzaNi5zF2hCdxEYrpwMoB9DICYEHz11QD7ojGgSUz0bzmpH03UU2AjIhJFomGSc1FvO2XbykpORoVdu/yjEDZuBJYuNZVPVey2U27ZMCds+fZaAFMBvHn33f4tqBo1/IMsxTMU2IiISFSowcDl00+BTZuAdeuAVavYBwR1bL8aBK3S8O8UBjyc43T4MKYDWM5cYt7voouAamZ2d4lRTk3JUGAjIiKR7dQpU77dg8HKN9/4g5o9ezj8z9wcuN7iBjXMozkGYA+Ai2+8Ef/3j3+YgIbbUKm8Qx5G7URid2a5MI7CEBERiTilbafgsh99ZBKETSvWGTP8KzY2qHE5AXk0LNPeBOArAC+ywmnUKCyxfWlMUBOGVRpeIj1IqhBFx3o+CmxERDyKuSk+n89c+O+8Yj5Lfj+mqN+Y2DV4FIAbuCIzcyaG2oonHD+e62ynHwHTf6b6yJHo+MoreBfAfPaiadnSzHxy7CpMtL9xy/lpK0pEJEaGTEaDynbLqQuA5vbv+AULTMVTqPJt2FUYrtJsBMDJgX/99a/hNGqE9QDO8g6liv+trigbDUrhKLAREZHwS001QUxvAC0AdLJdg6vytjNnQgc1LM9OTMTaEycwA8Bcmxj8F3YQrlDBH9REGSUYF54CGxGRGJDbKkLYVxg4UHLnTvhmzMA1AC4C0C6ga3Cuwwy4GsYxCD174i/vv49tNjGY21HFPdfpfKszCkzCT4GNiIiUvMxM4MgRYMECYMUK+FaswJUAatsme7l2lilTxgyqRJcuwFVXmQGVC99/3zTeC515I7FGgY2IiJSsU6dMMIPFi4H1602TPSQl+Sdx5/IhZ2z5dr3u3U1Agz59TMdgJyHBXC/iUmAjIiIlogz/YCCzfDmwYYP/b/4/ORk4dy5kUMPy7VKVKmHtqVMmj+Z3TzwBtGsH1Krl33IKqNwKTNgOV0WXhJ8CGxERKV7nzpktJlY7xX3yiT+YYdfggweBs6FTfDO5sGPzZpqPGIE/f/SRyaN5vHNnoE4dIE7dSiQ0fWWIiEjxYBLt8ePwff89bgJMLxpMmQLMmgXs3p1rUMMxCLsBfAfgJcZF996LxbbBHqugTqemFqg/j8QGrdiIiEiRS2BF0+rVwJo18C1fjuEAWvKGzZv9AU8IDHOYL7PZdg2eb1dsMlq0KNI8GvUb8jYFNiIiku+A4CC3kUJgngyHUl7CwOazz/xznVauRFcb7OQW1DBfZvXp02aVZobtR3PUrt7kZa5TQSnA8R4FNiIiUuheLnG290xH21yPM57ASdwHDpgGe+xJE1JCQlY/mucnTsRWAIdtfo26wUhBKLAREZGCcxxUtN2CuUrTxo5BYOk2duzI/eNKlwZq1gSYDHzddXA6dMCSiRPD1o8m7I0KpcgosBERkYI5cwbYuhXD7ArNxbZrcA0A5XL7GI5BqFoVaNUKuPZa4LLLgPr14ZQrV2z9aIK7ASvh2NsU2IiISP7fOPbtM0FN3MKFGGmDGubWVMyl3JarMCc5+6l1a2DoUOCaa4AmTYAaNYDy5YEUk01zQVpVkQtRYCMiEoPys/Xi3pdzmxIBtGLwwrJtNtlbtAiD7FTuXN9QypbFtrQ0fA/gzj/+ESlt2qBup04mj+YUk3a5iiNSRBTYiIjI+TmOSf5taPNnTGIwK55Yun3smJnAHTI0iY/3r8h07IgXpk/HCgBjevQwW1FcvYlUqpCKbgpsREQ8UGad2/0LLS3NNNPrZyueLrYTuLFkiekojBBBTTqA4wBqdeoE/OxncPr0wezp0811JmG4BGjKduxSYCMiIjmY6dqHDpmgJm7+fFxjy7ibAP7SbRvUBGIYwaqmXQCmA3jkz38GWrc2gyoPhOEcJDYpsBERiVLFtWVSyW47+aZO9c91mj/fjEOoep7p26wzYvDyA4DJALYA+E2HDv65Tqmp+T4vUqKwFIQCGxGRXMRa6/0EW9nUFkBnBjaTJ5uRCFy5qZlbHk1cHPZlZpr8mc8ALAdwBMAJ3paYWOBhlYFbSSrPlvxQYCMiEuNK267B27/6Cr4VK/DBk0+abSfMmYOUs2dNQBMyqGHg0rQp/r5yJWbbjsHMo0lDbFAeT2RSYCMiEkMJxu7KE+8fZ8u0m9jEYDPbaflyDLHbTpy+HRzQ8G2cZdqVW7QAhg+HM2QIPh8xwjTXy8u6Sm7bSwoSpKgosBERifAtsCLHAOLUKTMGgdtOnezWE9ib5uhR0zk4lFQ7bXshgFuefRbo2hVOhQrYE8a5TgqIJJgCGxGRWMJE3iNH4Fu1yiQEd7a9aRrwtqOcp53TWTtpexXnWgJYBuDmnj3NKASu6iiskEiiwEZEJEZ+2FfhdhOTgTdtgm/+fIwC0NJWQeWW4suAhhO3pwCYGZBHY/rRcJDlWYY9oWmwpISDAhsREQ/z2cClgZ28bfJoVq401U7MqykT4mO4AsMNsIpNmuD9HTvwOYD9Nsg5Xx6NyrQlEhSsDk9ERCIeG+k1AtADwHAAP+OVn35q+tIgOTlkUMMxCDsBfMHh3X/4A94CsBHA7jwmB58v6CmWfCGRIApsREQKgNssPp/PXCKuz0pamulH0wHA5QCuB8y200Detn8/kJmZ40M4BgHVqwO9e+N5AP+PE7l79TJdhNlNWHk0Ei20FSUi4qExCCzfjt+8GZcCuMQmBre25dtmTEKQDFu+zTyabvfeC2f4cHzXt68p33Zq1EDOwQkikU2BjYhIFAlZ3swVmJMn0QxAU/5g//pr3GCHVdbNJY+GuDG01851Yh7N9Ntug1O9usmnKYrjK6qVLPcxlYwseaHARkQkEnrL5OKCb+QpKcDhw/Bt2GC2mroxkPn4YwyyIxJC4piDmjXxzcGD+Njm0HAMgtOwYYFHIIhECgU2IiLROH+KZdbHjyNlyxbc1bevSRC+GUAr3rZnT8ighttKpapUATic8uqr8cJjj5nybc51MltOCmrEAxTYiIhEk4wM4NgxfxLw5s2mH80tANoDqHWeH+pMAN7GnJuf/xwYNQqZVati+2OP5ZjrpOGTEu0U2IiIRAmu8/h27jTTtrFihb9se/Fi9AVQLsT9GZ742ESvTh18vHs3JgH4dswYoHFjc3usDKuU2KLARkQkn0p6kGM5W9VUn4HK3LlmUCUWLAC2bTNbUqGCGgYthwA0vOwyODfeiP+7+25/Hk2jRjwBLscU+XGKRIKo2lAdN24cLrnkElSqVAm1atXCNddcg02bNoX7sEQkDCK6j0wR/uZZHUBzAF0BXMMr//1v4JNPgA0bQo4zMLkylStjKYC/MZD5059QfuRIbD17FiccBxUqsQ+xiHdFVWAzZ84cPPDAA1i4cCGmT5+O9PR0XHHFFZ79oSYiMSozE4kAmgBm7MEQALcBuJa3LVpkSrtzfAgAXruBqzW33ILHbQm307w5UKuWf65TLpjEzJWmAicyi0SQqNqKmjp1arb/T5gwwazcLFu2DP369QvbcYmIFAluY3Elau9eM9epfUCTveZ2REIoKQAO8Jc/AJMBfHj33dj8xhv+EQgJuRZ9i3hSVAU2wZKTmecPVKtWLdyHIiKSb9n64qSloVn16qah3qw33sDVNqhp507lzu1BatTAd0eO4FMAy+z07cz69Qs110kkmkVtYJOZmYmHH34Yffr0QQf2ZMhFWlqaubhOhljCFREJFzPm4OhR+A4fNttOzKUp/eGHZmBl7fP9kE5MBFq0AEaOxB+fecYkCmf1o4kPNTxBJDZEbWDDXJu1a9fi+++/v2DC8TPPPFNixyUikhdcgalok4N969fDt3q1abDXmXHJnDmmAipk+Xa5cgArm4YOBa6/Hpl16mD7M8/gTBjOQSQSRWVg8+CDD+Krr77C3Llz0aBBg/Ped+zYsXjkkUeyrdg0ZNtwEZEwSbDbS/Vs7oyPVU6LFmEEgErnyaM5yFlQw4YBt90GdOpkxiLA54vJoKYwXZs1e8rboiqw4RfiQw89hMmTJ2P27Nlo2pTj3s6vbNmy5iIiEnYszz50yAyqbGRzaLrz+o8/Bk6dMpO5g6XbvJnVAD4DMP4PfwCaNQNYts0tJ1WFikRvYMPtp4kTJ2LKlCmml82BAwfsVnMiEpT5LyKRPAaB+X0cVrllCy4D0AtAJ9t0j0FNqPLtY3YMwtesCgXM1G2HQQ3nPUXISklxNSUsCRof4U1RFdi8/vrr5u8BAwZku/7dd9/FGLYJFxEpgmGURbZFwTdNVjxxttO+fcC6daZz8O0AWtgtqZDVTmXKYNPZs5htV2m2AzgO+LecztOPRgoumgM0ieLARl90IhI1UlOBEyeAgweBLVvMTCfMmwesWYOLcgloUu0YhAYDB+LZqVOxggVTDLzCcPgi0SqqAhsRkYiXns4mW2bbCbt3A0uX+odVLlxoyroRIqhhHs1R2zX4CwDPjR2L2VOnmk7C3JISP62qSF4osBGRqN0yiiiZmf48GgYve/f6ZzlxYCUDmqQk/7ZUEMduMSUBmA7gSwC7AfypWTPTkyaSFPVzH9GvpUQ1BTYiIoXFPJrjx4H9+4GtW/2Tt3/4AVi/HjgTuhib20t7AcwHTNfgDYF5NGXKFPsha/VDvEqBjYhIAbGRBEu0fTt3mjJus+3EVZoVK/z5NaH4fNjlOFhlE4MXATiiPBqRIqPARkSiSiQ0VePAAk7frgqY2U6+WbOAJUv8KzWsfspN1apAq1Z4bdEis/XE8u285tHkZ+tGqzESyxTYiEhMKUxg5I5BYBeZOrZrMJvs4Z13/Hk0TBwOIbV8eaxIScG848fx87Fj8eE115gcmtD3FpHCUGAjIpIXKSlmKGUNjjWwM516AqZ0G5s3h/yQ1NKlsSU9HYtTUkweDbef7uzUyXQSFpHiocBGRDxX5VSUn+PHY8dQ4dw5pOzejbYA2gHoA6CLnfVUKpd+NLsArE9PxxQA82x/GpNHE0Vd0rWlJdFIgY2ISIiAKM5NDN6zxyQC+1atwnV2laaZ7RocEscdNGuGCcuX47+2fFv9aPwUKElJUGAjIlEveIUmUEFWbiraoKYWAxtWOK1cCcyejett0nDIMQjlywONGwP9+sEZORJvX3llgfNoCjp9Wr1hRBTYiIj8JDXVBDPVATQE0J7XvfsusGqV6VMTavRkCmOaRo2Azp2BG24A+vSBU7688mhEwkSBjYhEjWIr82Y105Ej8O3dizZ2QGVPe+FKTShnARwAsAXAwHvvBa69FmjYkAdpEo1FW08SHgpsRCRm+dw8GvaeSU6Gb/VqjADQGzCJwuVDfEymnevE3JlZACYDuGz0aKBuXSCeHW5EJJwU2IjEsHDPYQrn569ggxpuO/mYQ8OS7RkzcAuAarnk0ZwCsBMwU7f/A2AN4J/pxMZ7CmpEIoICGxGJie0qN4jiGISatsleI8BsPeHDD4HFi81YBAY6OcTFmRWZOXv3mhWa+bZ8+7CSdUUiDisaRcSD+Ebu8/nMhf+OeefOmZWY2gE5NCPZMI+3ffWVf9ZTEFPRVK0a0LMnnDvvxNO8q121Ya8aEYnywGb3bu4qi4hETzDm5tHEHTyIxgC6ARgM4A4AN9qGe8EcO5hyHf99/fXAc8/BuecebLLbUSLikcCmTZs2eOqpp5CijH8RiQJM/q1lOwTHrV2LvgBuBjAGQC+7HRWM4RQDmBkA/szA5re/NSXcFRo1Qobj5OiTIyJRHNhMnz4d3377LVq2bIkJEyYU31GJiBRCGTvTqbatbmKVU6nPP8fP7WpN3RA//E7YLaa5AMYDeA7ATN5Qvz5Qho8YvnJpXnLL5cnLfURiSb4Cm969e2PRokUYN24cnnzySXTt2hXz5nEKiohIZPxAq2JXaTj2oAeAywHcBaD0F1+Y3BomDwc6xz8qVzaVTh8CeAbARADbiymPxg1EtPIjEkHJw3fccQc2bdqEESNGYNiwYRg1ahSSkpKK/uhERPLCccwYBK7QNADQFUB/ALfaXJrONnk424fYVZr1vGnkSBPQ/MPm1SiPRiRGq6KuuOIK/PznP8fkyZPRrl07PP744/otRMQjoqWqqhwThA8fRh0AHezk7WsBs+00wO1TE/QxPJtttsHe3wCceeQRLABwDEDGBT6ftn5EPNTH5o033sCSJUvMZcOGDYiLi0OHDh3wi1/8AhdffDE++ugjE+B89tln6NaNtQciIsX3w+vEjh2okJmJlE2bTPk2V2b6AWhuA57cxiBstfkzLN3eA+DvDRua20QkxgKb559/Hj169DBbUT179jQ5NgkJCVm333vvvfjzn/+MMWPGYO3atcVxvCIS49zybW49+XbsAPbuNRO4R9sVG07fDjUGIa5CBaw8fRoLAXxqK5+45eREefdmESlEYJOXPjZ33323SSwWESlqDBkq2QRhrsr4vv8e4GX+fFO+Hby37tjghT+52gwdipc+/RSL8rjlJCLRqcg7D9eqVQvfffddUT+sSESLlnyUaOWOQahly7eZRzOcN7z6KjB1KnDyZI4fZqm2sol1m2/w/48/jmkcgxDBQQ1XfJS7IxJhs6L4g71/f9YjiIiE3rLJM1YyHT1qgpqGAJrYTsED3RlP+/eHHIPA6dubAXxvp2+zZvPPTZpEVB6Nm4RcoOdFRHKlIZgiUqwK8sZtqpiSk4HMTPj27MElAJoCuMwmCFcNUenEPJpku0qzDMBnAFbb6wqbRyMi0UOBjYhE3BgE5tH4DhwwKzK+tWtxJYBLAVPSHfKHVlwcNmRmYiOAL+z2UyRvOYlI8VFgIyKRIS3N5MpwFAKHVeKHHwBWV06bZvrSBHcMzlKvHtC6Nd6fNcuUbzNROJK2nESkZCmwEZGwirMl2myyxxWaXnbbyffuu/7AJj095BgEVjbV6toV6NcPzvDheHvWLKg9qIgosBGJAZHaZ8X344+m0ok5M1i2DL7t2zEKQHf+fwWnN2XHXJmTAHbY/JnbHnkEGDwYTkKCghoRMRTYiEiJK2eb7PkOHjTl2xxY6ZsyBZgzxwytDDVLOwXAXjvLabq93Dp8OFClCjOUC3QcbqAXqkQ/0oJAEckbBTYiUmIrPqXsthP7lTcCEL9woZnndAVvnDTJDLMMFdRw7AETgznPaYqd85TGG0qXRrRzy75FpGgosBGREsmjYaUTQyL2pGlpe9KUnjjRDKtkJ2EGNTlUrIiE1q3RskcP1L/8ctx+7bWmfFtEJDcKbESk+DBYsXk0le2WExODewDoyxWbpUv9+TVBmEdTuVUrgMnBQ4YAl18Op1IlBTUickEKbESkyGQbCHnkCCpkZJhqp9Y2oLkIMDk0Dc9Tvr3D5tGMuO024IYbgKZNgTJlCpxHIyKxRYGNiBSpeLd8+8gR4PhxYM0a9AMw2I5DCJWhw/LtUkwCbtMGny1caMYgDLv3XqB27ZCfQ4m9IpIbBTYiMaRYy74zM80YBG47VeP/lywB9u6F7+uvcZ/Nowk1dfcEgK0A1pw4ga8WLsyavs38muISjsBIwZhIyVBgIxLuLZsIecMr1PFwm+jUKbPt1N5uO71/++1mWGV9N9AJwr4zFRs1wsxdu/AdYCZv79NcJxEpJAU2IlK4MQgcVpmaCuzcaToFc1DlMDsWgWXdwTju4CCA9RxqOWoUxr38slmx8dIYBJVwi4SPAhvxnEhcDfGcjIyfApr9+4FNm4DNm+H76itTvl0xxPTtDLvtxH40fe67D62uuw6n27fH+pdfLpZD5GtP+Z0sLiLRTYGNSDHxYoDF82E/mgNbt6L8mTNIWbUKj996q5m8fXWlSmY7ircHOwVgO/NobMfgXo895q92YmAUYbzyWonEKgU2ImEWLcGPOwbBbC+tXw/s3g3fDz/gIQB1eN0phi/ZneHH1a+P7/buxfcA/munb79bpw4QFyqVuHBbP4HBZOD1FGpsgoh4T9H9ZBGRYsE3ZJ/PZy7Bb87nu62ocGhBdZsAzG7BfbjNNH068MorwMSJZjRC8BiEdAD7AfzA3JlrrwU3m96xPWq4JVUQbpDCSyQHgCISXlEX2MydOxdXXnkl6tWrZ36Yf/755+E+JJEiwZWG4gxQ8i0jw/SjqQGgbkC34Nt425tvAhuZLZPTcQ7mBvAVgFeZX3z//VgCIPI2nUTEi6IusOEP/YsvvhivvsofmSJSXGMQWLpdw3YLZkBzLWC2nfq41VDBSpUCmjUzZdvvAngJwLe8vjrXe0RESkbU5dgMGzbMXESk6BKbif83Yw4OHQJ8PlPlxKTgTrZrcINcxiCwTPsAp3VfdhnQowdeeu45bMlj+XZuuTEiIjGzYiMieeMGDW7gcl7nzpkcGl58u3YB8+fDN2OG6Rh8mx1eGXK2U+XKWAjgMy70PPggMn/xCzPnKVw9ac6Xh6McHZHYEHUrNvmVlpZmLq6TJzk3WMR7ClJdZXrNJCfDl5ICTmVqzuvmzAFmzjQjES6ys58CZfI3Ig6lbNkS6NQJ//zgA8wF8Cuu2BRBpVOkV4eJSGTz/IrNuHHjkJiYmHVp2JBzhUWkPGDmOplhlWvWmDya63gDq53mzQPOnMkR1DABOAnAJ2fP4v5163D6vvvA9P1DvLEIy7dLotpLRLzJ84HN2LFjkZycnHXZvZtdNERiF0uzmRhc1W4x+ebNg++//8WDbLLHOxw9muNjzto5TlyZ+RjAPwBM4OpN69ZmMne00zaViHd4fiuqbNmy5iJSHKIp6TXeNthjoz3fvn3oCaA1bxg/Htiwwb8NFfQx7DnDzdtNAFbarsFZ07dFRCJQ1AU23H/fupUj8/ySkpKwcuVKVKtWDY0asVWYiGTjOGbMQQUb2HCVptTcubgdQDfevoaDDnIGNSzf3nTunOlJw5Wa2XbVpqAN9kRESkLUBTZLly7FZUxStB555BHz9+jRozFhAhfHRbwr3wnCKSnwHTpkqp04bbsFgN5cyXzjDfTPrdKJOPKgUyd8OHWq6UXD8u0QnWtERCJO1AU2AwYMyJr9IhJtSmzr6uxZlgD6h0zu3IleADoDGGL70eDIkRxBzTnbNbhm9+5Anz5IGDIEL334IZ4pVQoVOeDyAqJlS05EvC3qAhuRWFGQIMhUAxw/7m+wd/gwsGULfEuXmn407ex2VMiKgdKlsTY9HcsA3HXPPcDQoUDdukB8PA+kiM7ImwIHbYpI+Hm+KkokVkqMK9jy7U4NG2JIgwZI/fpr/0yn8ePNak2lEN/wnL6N+vWB4cPxb1vt5Fx3HdCggT+oERGJMlqxESkiRdVYLr8rNdxSSrR/N7LJwQPcfjQHOOwg5/Rt2ETgtQCuuP12OP37470pU5BiHlBVhCISvRTYiCdmHcVk75H0dNNzhiMmawJoBaArgMttgOMGNcFdg5MBVLn4Yny2apWpdBrMraeaNf1BTTGL2ddKREqMAhuRKOOOQTD/PnDAbDO154BYG9ywpDvUJhI313YAWA5g1O234/VVq7CHN9TmMAUREW9QYCMSZWMQmCvjY3fg/fvhW7YMdwPowlWYXLacaBdg+tFwYOUcANdcdx2SfvvbEjlmJdeKSElSYCMSBcoE5NHUs2MQMHu2mcLNfjQJuX1gYiLQti0+XrjQbDtxtcas9ZxnOyjmt/hEJKopsBGJgjEIDFyq2wZ7PXjD3/4GcO7ZmTM5ghqujZzix7VpA1x6KZxevfDWwoXYb3vViIh4mQIbkQgVPAbhooA8GvanCYXdgXfbbaef3XUXMGQInLp1zXUiIrFAgY1IpElJAdN5E2yX4DYAhrM/DWAmcocqxuZKzGEAqwHMB/A9gOtvvtn0qKng85kcl2ga2CkiUlAKbCQmRWK5eGmbRxN37Bjq25WZwYAZh+AGOjkGVdpJ21vs1G0GNEttkIOqVf0diEVEYogCG/G0SAteQokLyKPhsMq45cvxMwCDbD+airm1CC9TxqzIfJmUZFZpFthtqLRiqmrSio+IRAMFNiJhVNFeGHI1sYMqy736Km615duhvkHPMqbh9O1u3YAePfDPJ5/EJgZvYTh+EZFIo1lREvO4ClHS86HK2rlO1WxAwxEIDwBmWKVv5UrUCBHUZAA4ZFdmcMstwC9+gcxbbjGJwpEU1HB1jBcRkXDQio1IIXJz8i093QQz3HZi8NKas5psCXcd24APIcq3TwLYaYOaHwD0veMOgOXc586ZMQkiIuKnwEZiQqEDksLKzDRjEHyHD5tE4KZ2lWaArXyqlMsYhFQ7rHKJrXZaDGAbb2jRwj+s8ty5AnUAjrTp5cHdiSPt+EQkeiiwESlmXIVhQIOMDGDTJowAMBRA7ypVgBMnQn5Muq1sWmMDGo5CWG8roIpjOIEShEXEKxTYiBSXtDQzddt0Bt6+HVi3Dr4vvsC9dhsqVFDDbaW4ChWw8vRpU7rNoGatXbVhsBMtNB9KRMJFgY14erspLKsP3B46eTJr26kdV2xefx1YscIMrmSPmlDdZU7bcu3mQ4fizU8/NUnBnMadUvJnICIStVQVJVJUuEJx8iRw6BBw8CCwcqUp236Ut02fDiQlmdlOvhDl25y+PRPA21yZueMOfANgw3mCGgZxrOQqycDNXYVRxZOIRDKt2IgUAdMVmAFNZiZSNm7EM1deaaqdfl21KhKOHwdOcSxlCJUrY/HJk5hr82g2Aniqa1eTSxMNlJsjIpFGgY3EnPy+AQe+aR/kSkygs2dNvowp0961C1i1yuTR3AWY/BowqMmlfDuxTRsk9O2LywcMQO/WrVGjWzec4R1Kc7hCdmHdWhMRiSIKbCQsInFWU773cI8fh+/06Z/yaF55xQQ22LvXlHCH3OeNi8O2zEys4pDK++4DevUCWrZEhapVkVrAUmwFPSIiP1FgI5IfjpM1BsFsPW3ejNFsmMfbZs40vWqQW1BTrRrQvj3emzfP9KW5dvhwoFkzoJS+DUVEiop+oorkUTnbj4ZbTPX470mTgO++wygO0uYdbFCTA1ejWrYE+vUDevbE5/PmYQ+vr19fQY2ISBHTT1UpsFAJo9G4rXRB6emobvNo4vbtw+UAhvD6994Djh0zM59CrdBwg6hi8+ZI7doVP//4Y6xduRIzHnzQ3zk4zNRnRkS8SuXeIhcag3DokAleejGv99VX8Ut2Debt+/ebJnyhvonYj2YO/3H33XDuusuUcrPRnpOYWCyHqlJsERE/rdiIhMIEXvakYZn28uUYA+BSBjazZ5s5T2Vy+7jKlYHWrfH+kiVYBGDYiBFmyymolkouQCtKIlJQCmxEAjBgSXRnOx04AMyZA9+0abgOMNtRDHZyFmP7y7crt2oF9O0L9OmDz5YsMdO40bw5vEZBh4hEMgU2ErFl3SVZEs7J2r7jx03pdjX++8svTVCDtWuBI0dQN2AMgi9o+vYBwIw/uO7++4GuXZFZr54ZWKm3fhGRkqfARmIagxS3fDtu/36z3TSIN7z1ln80gs1ZCc6j4UDKI3bswcDf/x7t+vcHWrTwVzqdO6egRkQkTBTYSNTLtStwHsYgVLZBTSt+M/zrX3gAQEM2FN6+3T9pO8THsZfwNjt5mxO4L7v6aqBtW6BSpZ+GYBbzFlB+m/iJiMQKBTbimW7AeVXa5tHwTBvZaqd+zK/54gu0BVDWruQED6vkQMqERo0wddcufAdgM4Dt3HIKDGouQPkpIiLFS+XeElNf7FUAk0fDgGYYgIcA3A7gEt7h6FHThC8uKKg5G1C+fXb0aEwAzNBK5tUcNQ9cfN9GDDAZCEV6kOkGbNFwrCLibVqxkZjg+/FHE9AwsOnC7SMAF9sOwly5cSudAgOaDLvttAnAPDt9u9/gwVj67LMmaTgSaUVIRGKdAhvxtHI2j8Z34IAJaDjTqY/No6lsy7vjgvJuGLScArDDBjPzbE4N/5/ZunWxBzUKTkRECk6BjUQl980/1FgHIz3ddA3mKg0b6pX+5BOTGNzClnOXDbHlBHdbqXZtzNi/33QL5mpNEhOT7QoOfDk+QkREIogCG/EUX8AKDfbuxUjbMbjMpEno6s57ChHQsI6pVNWqZvo2Lr0U777wgmmwx9yaM2E6FxERyT8FNlKicl1hKQK+06dRB0AN/nvGDGDePIy2eTTYs8fk0gQHNNzwSbZbTV2vvx4YMgROnTpY8sIL5noREYkuCmyKQKyUSUf8GIR9+9ADMBe88w6wY4fZenLnOgUHNal2RWY5YMq3x994o+lH41SqFLNBjfJ7RCTaKbCR6MVGeMeOmRUZBjClJ03C3bbZHtatA86ezepJE9w1+JCdtj0LwGob4DidOwPVq/sHYIqISFRSYCPRt12VmekfdXD8OHwbN5o8mt5cmZk82eTRmI88ezZk+fYJAFsCyrcZ0Oy0VVAoxxoqERGJZgpsJLqkpAAnTgD79wPLlsE3ezZuA8yQSuza5U8cDvoQbqz4SpXCxnPnsNSu0uywl8MaViki4ikKbCQ6cAWGQQ0HU65eDcyfDyxaZPJoWtrybQruAZwGYD/71gwYgH/OmGGGVu4CsNduSYmIiLcosJGIxkDFrMLs3m2CGCxZ4g9qNm8GjhwBzpwxTfhCdZdJqFYNCW3bosqAAUjt1Anfz5hhtp3887r9id5UXFVaIiJS8qJyVtSrr76KJk2aoFy5cujRowcWL14c7kOSYlDRbjF1YODy5ZfA228DH30E8PXes8cENQgKajiR+yT/wcGUo0YBDzwAXHklMrt2xfqAoEZERLwp6lZs/vOf/+CRRx7BG2+8YYKav/3tbxgyZAg2bdqEWrVqIZZ5puw8NRW17NiDjgC68boPPvDn1Zw6BWSYHsA5xcVha2Ym1gC4/t57gU6dgCZNgHr1TCdi5dKIiHhf1K3YvPzyy7jnnntw5513ol27dibAKV++PN5h3xI5b9Dj8/nMhf+OSByDcOQIfFu2mCGV1wG4BcDlvG3rVn/ScIighl2D93G1pmdPvA/gPf770kuBbt38gU0Zt5NNTppKLSLiLVEV2Jw9exbLli3D4MGDs66Li4sz/1+wYEHIj0lLS8PJkyezXSTCsHybQcu2bcDcufBNmYIbAVxjJ3DXcYOeUCpVMg32PuZrfeedptHeCl7fpg2TZ4r8UN1AyM3PERGRyBJVgc2RI0eQkZGB2rU52vAn/P8BzgYKYdy4cUhMTMy6NGzIDQ4JB26TBSfqco3Et2uXP2/miy+A998HPv3UdA+uHzDbKYfSpYEWLYBhw/BPAF8xPurSBasAHA0xrFIBiYhIbIi6HJv8Gjt2rMnJcXHFRsFN+JW1YxCa8T/TpwOrViF16VIcXLvWXM+Lk1tQU6cOwC7BAwaYlZnFH39sGu1lNmjgn8AtIiIxK6oCmxo1aiA+Ph4HDx7Mdj3/X4dvdiGULVvWXOTCycZF4ULJy/E2aGkAoL1b8cRKJ5ZzHzliBlgit4AmMdFf7dSvH9CjB9C4MTJr1cLmIjlyERHxgqgKbMqUKYOuXbti5syZuOaaa2x6Rqb5/4MPPhjuw4tqxd3LhRtDleyk7TYA2nGaNoDWvHHt2qz5THEhyrcZclVmzkyvXkD//kDLlgBX3VjtZEu+Y22AYyydq4iIZwMb4rbS6NGj0a1bN3Tv3t2Ue3PlgVVSEpkSANRkaxkbyHQCcJHtUWPWdE6fxukQQU2qrXZi3sx1Y8b4q5waNAAaNQIS+KgiIiJRHtjceOONOHz4MJ566imTMNypUydMnTo1R0Kxl+WnX80FB0oW9xiEfftMH5o2duupq121qcIVuIA5TYFbT+l2htMGO6iS852u5fYTV22qVg3LqeSnL1DwakrElteLiHhQ1AU2xG0nbT1FbtM/E6QcPWqCGt/y5bjCNtlrCqAah2jblRnm2yTYlRnKsFO2mTOzzAY1O+1sJ6dDB1PaXVjawhER8baoDGykcKsy7n2Dk7CLAh+ZNWe+mTOBpCT4Fi7EVUz8tttOpWxQEzzbKcUGMCsBzAewyQ6q3MOFH94hLqo6E4iISJgosJEiUc7mzDApuDkDl//+F9iwwazaNGLbmfOVb9eqhTmHDmGBDWz22yCHqzciIiL5ocBGCryqY7Z1mEfD+6xfj9Q1azD+t79FZ944dy6QnGzybOJDBDQcg1CKW0ss3+7WDe+/9prpRbPL5tdos0hERApCgY0UfAzCsWPAxo0/XRYuNHOdqvP2wwxPcn6RsXybqbRJHHA5YoS/yV7jxlj+2mtm28mLDfaU1yMiUnIU2Ej+8A2aE7Y3b/ZvNXE45bJl/sDm0CHTeC/eBjCpgR9mc2Xc8u3FDGxuuglo3hyZNWqYJGEREZHCUmAjeZea6u8QvGaNf2DlypVmCwr79/uHWJ47Z7acGNgEbzsdteXbLN1eYhODMzt29E/fTmHqcNFXZxUnrcKIiEQmBTZyQQxUOLDCN2eOSQbG6tVmthP27PGXdTPPJpdtp7iyZbEyLc1M3GZQs8Wu2uy1ScPBwypFREQKQ4FNGEXaKkQwU5Z9+DB6u5VO8+b5V2u2bzfbTmYFJ8SqhWO3oZgz02roUPTt0wfd2rTB61ddZZKDL7w+8xP3eVGTOxERyQsFNhISW+G1ZJn2zJnow3wYXvn110g9dAhnk5PNKo7bYA8hugavAbCc09WZR9OiBTJr1sTGYjhObQmJiEggBTaSDeegN7aznDj6IH72bIy0/2bDPTePJmQ/mooVseLHH7PyaHYA+F2nTiawQVoaooUbLIV1HIWIiBSI2rlK1hcCOwYf+/BDrHzxRbNKw4Cm1IIFaGU7BzOoQUD34CxlyviDl6FD8RGALwCTU8PqJzOFu1R0xs9ugMNLpG0TiohIaNH5jiNFqobdamrCgGXhQmDFClwGwIwVPXgwx3YTcfPHx8TfOnUAVjd17w506ID5n3xi8mhOlPxpiIiIKLCJacnJ6A6gNYD6dvo2pk41wUxTdznPcUKWbx8HULNrV6B3b+CSS/z9aKpXN7k1yngREZFwUWATo3k0LZgYPGMGetqAppXdijIVT+npWdO3c35wWaxNS8NqAHcwMZjBTf36/i2njAwFNSIiElbKsYkSTGTl1g8vgaXPTG7Na4JrnE0MHgqgLwOXhQsxDMClNripyTuls67JP7TSzaPJ6iLM4GXgQEwG8AlXZvr0AZgc3LIlUI5jMEVERMJLKzYxorbNo2Fgw7EHHRi8zJuHiwGUt6szTqg8GqbZAFgPYNCVV5ptp1nffGN61Djt2wMcZCkiIhIhFNh43bFjpsFeS9s9uJ3dhmJODccjuGs9wdtOXLfhulCV9u3xzbp1WAZg3Guvma7BDGpMfVRcwRb8IrEZoYiIeIMCG682ieN21fr1iFu1Ct0AtLUBTWM7fbsM73PunFmRyRGelC6NpPR0rARw9TXX4LN168yqDaudToXlZERERPJGgY0HX1CWbcd9841/MOXatbjcXlfbJg4zmHGCVmrc6dtlOb+pbVt8MWeOabQ3sm9fU+nEbsIlSas6IiJSEApsPMJnuwNf7HYJ5uTtzZuBLVvQ1ebRuCszcUEN9tLtSgzzaC4dMcJUOn09Zw72A8ho377EgxoREZGCUmATJQMyz6eaTQZuavNo2vDK2bOBvXuB48dNHo0TFNi4frSjD9YyFmJgc+21yKxb16zSnHa7CouIiEQJBTYRGswcPMislvOrYIOY1rZ7cGs7hZt5NMyvwVluLvlXZ9xVGncbipObytWpg+8OHDArNdx24gpNZocOQM2a/qBGREQkyiiwiUZnzphkYFY41bL5M61tQFOXAQvvY4Oa4K2ndNs1eBOAXsOGYcq772I7YKqdeAFzbERERKKUGvRFkXjbHbjUrFlmSCW3nwbYBnvd7FZUheABlfb/Z20ezToA3wKYwiBn+HBTxs0uwruZTxOm8xIRESkqWrGJBpmZwOHD6G9XZeLXrUMvu1LDS9WAF5LBT3Af4hS7zbTGXphLs5fV3u3bY3MYTkdERKS4KLCJcJW54rJoEXzbtpkVGubQlJ471wQ2NW35dq4ddNhAr0oVfH/smAlgHnzzTVzTvDlaDRpkAhtUZ0cbERER79BWVAnNdCrAA5lAhnOdfCzdXr4c/ey2k2/NGlPSzaAGAUnB2aJUJiFfdBEwbBi+AvAFF366dUNmx47+cQiFPE8REZFIpBWbYhRYrp3XZnNl7LiDUrNnm8nbpifNrFnAzp3obldwkJKSlUcTWPHEMQelWJ7NYZUMajp3Ns32Fn3wgcmhQSvO8I7s7sxqzCciIoWhwCZSnDvHGm+TR8ME4fiNG00gwxwaLF1qZj5VCbi7z+bT+GzSL5OD2UivkTtxu2tXoHlzZFarZnrUOCUYvASWrUdKwCgiIrFBgU2YGu65b8gMTBiw+BYuhG/7dlPGzTyaUt9/b1ZsWL6NQ4fMfTMDKp7cF44BzUkAWwBsAHAXJ3B36QLUqwc0amQCJiefqyiF2kITEREJIwU24cKA4uRJdOIOEQOWVavMdlNfW/nE4ZWNgpKg4gP+z340qQB22lEIrHZaBWAMV2y45VSlyk8rQREuIoeIiohIVFJgEw4pKSaI8S1fbuY4cQwCfvjBXHcJgET+n0nIQVtI8XbbKb5MGVSuVQulW7fGtzNnYrntQ8OkYObUoFKlYt8OEhERiUQKbEoQE4PZ19c3dy5w9Ch8e/aYxnrMqcHixeY6zn0KFLhiwzEIyXwMJgUzj6ZtW3wzc6bpUbPXruKYEu8ggVtf+d1mKqm8GRERkaKgwKYkpKeboKWHTQb2bdgA7NplJm9nVT4dOJCVRxPyxSlfHhtSUsy20y2cwH3JJVmJwewoLCIiIgpsilWcLc8eULUqWgBob1dnMufNw4Jp07I6B8fl8qK41U7lmQTcpg2mTZuGFQBu7t8faN3aBDtuUJPXCiHls4iIiJcpsCkmFW0ScEvbIZhbUC1sYBO3eLGZyh1Yrl0qYLvKzGyKi8PhzEwzrHLDrl1YsWsXttpBlU7Hjv7kYFUviYiIZKPApqilpJhhlK3sFhOHFvD/jWwZd23e58QJ88RnBPSjCcyj+RFAtTZtMGP9ejOgkoMr99sEYc59QunSxXoKWtUREZFopcCmqKSlmQZ7vtWrzdiDKnbFhltNzWw34XIBd8+wW1BxgYFNuXJIOnMGGwEMHzIEX65fj112leZIeM5KREQkqmhWVFFhH5o5c+Bbv95sOTFRmInBveyKTYIt3XaTg0vbqNJxV2Hq1DEJwdMBM9spY+BA05dmbRQGNe6KDy/qDCwiIiVJgU1R2b8f2LzZlG0zoOkDoJ3NtXH7zyDEttNxG7xg0CBg+HB8DWAJ++p16GC2ngKrpGKZGywxSVpERCQ32ooqKuvXAwsWABs2oItNAkZAHo277eSORDhpt5g22sCm+5VXwmnY0HQQPsY7lC8fvnMRERGJUgpsisp//wusXAmcOZMV1LiznUoHTN9m+XaFOnWw6MABM9vp4X/+E7eyW3CDBjidmOgPakRERKRAtBVVVDZt8pdfZ2SYVRrHbju59UtnABy1M53ODRxotpy+YfDTtStOt24NX+PGqOjOd4pw3A5S/oyIiEQiBTZFpVy5bMtgbo+adDsGgT1o5tlgJn3IEMy3QQ5atsz2sSIiIlJw2ooqKvE/pQX77DYUp28fBrDZTuBeCWAHgMcvvtj0phEREZEYXrF5/vnn0bt3b5QvXx5VIm3b5iTTgX/adjpkA5kZTL8BMBUw4xC4SpNZt25YD1VERMSroiqwOXv2LH72s5/hl7/8JSJOmTKAzwckJmILgLkAvgVMLs1CG9Bst8nDIiIiUjyiaivqmWeeMX9PmDABEeeSS4B160zOzAw7rHKHLek+GO5jExERiRFRFdgURFpamrm4TgZsGRWpkSNNyTbatsU306aZ2U57bIl3UQme4K2ZTiIiIjEW2IwbNy5rpadYde0KNG6MzBo1zOBKr83dVhAlIiLRIOw5Nk888QR8Pt95Lxs3sj9vwYwdOxbJyclZl927OaigGLDJXufOQJs2ngtqREREokXYV2weffRRjBkz5rz3adaM87ELpmzZsuZS7LhFxAub9ImIiEhsBjY1a9Y0FxEREZGoD2zyY9euXTh27Jj5OyMjAys5mwlAixYtULEi52iLiIhILIuqwOapp57Ce++9l/X/zsxpATBr1iwMGDAAXhBY+cRk3dOnTytoExERiZbk4fxg/xq+2QdfvBLUiIiISAwFNiIiIiLno8BGREREPEOBjYiIiHiGAhsRERHxDAU2xTR6wL2wyklERERKRlSVe4toZpWIiJyPVmxKEFdv+KYcOKFbREREio4CGxEREfEMBTYiIiLiGcqx8SjlooiISCzSio2IiIh4hgIbERER8QxtRYVhSyj4Ok7wFhERkcJTYBPhlCsjIiKSd9qKEhEREc9QYCMiIiKeocBGREREPEOBjYiIiHiGAhsRERHxDAU2IiIi4hkq944AKukWEREpGlqxEREREc9QYCMiIiKeocBGREREPEOBjYiIiHiGAhsRERHxDAU2IiIi4hkKbERERMQzFNiIiIiIZyiwEREREc9QYCMiIiKeocBGREREPEOBjYiIiHiGAhsRERHxDAU2IiIi4hkKbERERMQzSiHGOI5j/j558mS4D0VERETyyH3fdt/HcxNzgc2pU6fM3w0bNgz3oYiIiEgB3scTExNzvd3nXCj08ZjMzEzs27cPlSpVgs/nK9JIksHS7t27UblyZXiRzjH6ef38SOcY/bx+frFwjieL4fwYrjCoqVevHuLics+kibkVGz4ZDRo0KLbH5wvoxS/SQDrH6Of18yOdY/Tz+vnFwjlWLuLzO99KjUvJwyIiIuIZCmxERETEMxTYFJGyZcvi6aefNn97lc4x+nn9/EjnGP28fn6xcI5lw3h+MZc8LCIiIt6lFRsRERHxDAU2IiIi4hkKbERERMQzFNiIiIiIZyiwyYfnn38evXv3Rvny5VGlSpU8fQxzs5966inUrVsXCQkJGDx4MLZs2ZLtPseOHcOtt95qmhjxce+++278+OOPKGn5PY4dO3aY7s2hLpMmTcq6X6jbP/roI4RDQZ7rAQMG5Dj+X/ziF9nus2vXLowYMcJ8bdSqVQuPPfYYzp07h2g4R97/oYceQuvWrc3XaKNGjfCrX/0KycnJ2e4Xrtfx1VdfRZMmTVCuXDn06NEDixcvPu/9+bXXpk0bc/+LLroIX3/9db6/J0tafs7xrbfeQt++fVG1alVz4fEH33/MmDE5XquhQ4ciWs5xwoQJOY6fHxfJr2N+zi/UzxRe+DMkUl/DuXPn4sorrzRdf3ksn3/++QU/Zvbs2ejSpYupjGrRooV5XQv7/Z0nrIqSvHnqqaecl19+2XnkkUecxMTEPH3MCy+8YO77+eefO6tWrXKuuuoqp2nTpk5qamrWfYYOHepcfPHFzsKFC5158+Y5LVq0cG6++WanpOX3OM6dO+fs378/2+WZZ55xKlas6Jw6dSrrfvwye/fdd7PdL/D8S1JBnuv+/fs799xzT7bjT05OzvY8dOjQwRk8eLCzYsUK5+uvv3Zq1KjhjB071omGc1yzZo1z3XXXOV988YWzdetWZ+bMmU7Lli2d66+/Ptv9wvE6fvTRR06ZMmWcd955x1m3bp15HapUqeIcPHgw5P1/+OEHJz4+3vnrX//qrF+/3vnDH/7glC5d2pxjfr4nS1J+z/GWW25xXn31VfO1tmHDBmfMmDHmfPbs2ZN1n9GjR5uvg8DX6tixY0645Pcc+XVWuXLlbMd/4MCBbPeJpNcxv+d39OjRbOe2du1a83XL847U1/Drr792/ud//sf57LPPzM+CyZMnn/f+27dvd8qXL2/eL/m9+Morr5hznDp1aoGft7xSYFMA/OLLS2CTmZnp1KlTx3nxxRezrjtx4oRTtmxZ58MPPzT/5wvOL5IlS5Zk3eebb75xfD6fs3fvXqekFNVxdOrUybnrrruyXZeXb4JIPkcGNr/+9a/P+w0fFxeX7Qfv66+/bn4wp6WlOSWpqF7Hjz/+2PzASU9PD+vr2L17d+eBBx7I+n9GRoZTr149Z9y4cSHvf8MNNzgjRozIdl2PHj2c++67L8/fkyUtv+cYjIF1pUqVnPfeey/bm+LVV1/tRIr8nuOFfsZG2utY2Nfwf//3f81r+OOPP0bsaxgoLz8LHn/8cad9+/bZrrvxxhudIUOGFNnzlhttRRWjpKQkHDhwwCyRBs654HLbggULzP/5N7cLunXrlnUf3p8zrRYtWlRix1oUx7Fs2TKsXLnSbH0Ee+CBB1CjRg10794d77zzzgXHzkfaOX7wwQfm+Dt06ICxY8ciJSUl2+Nyy6N27dpZ1w0ZMsQMgVu3bh1KUlF9PXEbiltZpUqVCtvrePbsWfM1Ffj9w/Pg/93vn2C8PvD+7mvh3j8v35MlqSDnGIxfi+np6ahWrVqObQBui3KL8Ze//CWOHj2KcCjoOXL7tHHjxmaQ4tVXX53teymSXseieA3ffvtt3HTTTahQoUJEvoYFcaHvxaJ43nITc0MwSxK/8SjwDc/9v3sb/+YXbiC+mfCHlHufkjrWwh4Hvznbtm1r8pAC/elPf8LAgQNN/sm0adNw//33mx9azOMoSQU9x1tuucX8gOXe8urVq/G73/0OmzZtwmeffZb1uKFeY/e2klQUr+ORI0fw7LPP4t577w3r68jjyMjICPncbty4MeTH5PZaBH6/udfldp+SVJBzDMavR35tBr5BMBfjuuuuQ9OmTbFt2zb8/ve/x7Bhw8wbRnx8PCL9HPlGzsC5Y8eOJsh+6aWXzM8VBjccYhxJr2NhX0PmlKxdu9b8/AwUSa9hQeT2vchf+FJTU3H8+PFCf+3nJuYDmyeeeAJ/+ctfznufDRs2mGREL59fYfELdeLEiXjyySdz3BZ4XefOnXH69Gm8+OKLRfaGWNznGPgGz5UZJisOGjTI/LBp3rw5vPQ68ocOExjbtWuHP/7xjyX6Okr+vfDCCyaBm7/ZBybX8rf/wK9ZBgj8WuX9+LUb6Xr16mUuLgY1/KVp/PjxJuj2EgY0fI24Choo2l/DcIr5wObRRx812efn06xZswI9dp06dczfBw8eNG+GLv6/U6dOWfc5dOhQto9jNQ0rVdyPL4nzK+xxfPLJJ2ZJ/I477rjgfblczB9OaWlpRTJHpKTOMfD4aevWreYHDT82OJOfrzEVxWtYUud46tQp81tipUqVMHnyZJQuXbpEX8dg3PLib6buc+ni/3M7F15/vvvn5XuyJBXkHF1cxWBgM2PGDPOmd6GvDX4ufs2W9JtiYc7Rxa9FBtM8/kh7HQtzfvzlgIEpV0MvJJyvYUHk9r3ILW5WsfE5K+zXRa4KlaETo/KbPPzSSy9lXcdqmlDJw0uXLs26z7fffhu25OGCHgcTbIOraHLz3HPPOVWrVnVKWlE9199//715HFZiBCYPB2byjx8/3iQPnzlzxomGc+TXZc+ePc3rePr06Yh5HZlc+OCDD2ZLLqxfv/55k4dHjhyZ7bpevXrlSB4+3/dkScvvOdJf/vIX8/W1YMGCPH2O3bt3m6+BKVOmONFyjsEJ0q1bt3Z+85vfROTrWNDz43sJj/nIkSMR/xoWJHmY1aKBWJ0ZnDxcmK+L3CiwyYedO3eaEku3pJn/5iWwtJnffCyHCyxJZPkavxhXr15tstxDlXt37tzZWbRokXnTZKltuMq9z3ccLCfl+fH2QFu2bDHfcKy+CcYS4rfeesuU2/J+r732mikBZOl8OOT3HFn+/Kc//ckECklJSeZ1bNasmdOvX78c5d5XXHGFs3LlSlPOWLNmzbCWe+fnHPmGwMqhiy66yJxvYHkpzy2cryPLQfmDf8KECSZou/fee833k1uBdvvttztPPPFEtnLvUqVKmTc8lkI//fTTIcu9L/Q9WZLye448flasffLJJ9leK/fnEP/+7W9/a4Iefs3OmDHD6dKli/k6KOlAu6DnyJ+xDMi3bdvmLFu2zLnpppuccuXKmZLgSHwd83t+rksvvdRUCgWLxNfw1KlTWe95DGzY+oT/5vsi8fx4nsHl3o899pj5XmSLglDl3ud73gpKgU0+sPyOL2jwZdasWTl6fbj4m8WTTz7p1K5d27yAgwYNcjZt2pSjpwHfeBgs8bewO++8M1uwVFIudBz8Bgs+X+IbeMOGDU20HYzBDkvA+ZgVKlQw/VXeeOONkPeNxHPctWuXCWKqVatmXj/2hOE3amAfG9qxY4czbNgwJyEhwfSwefTRR7OVSkfyOfLvUF/XvPC+4X4d2f+iUaNG5s2cv+GxP4+LK0z8vgwuVW/VqpW5P8tN//vf/2a7PS/fkyUtP+fYuHHjkK8VgzhKSUkxQTaDawZ1vD/7gxT2zaIkz/Hhhx/Oui9fp+HDhzvLly+P6Ncxv1+nGzduNK/btGnTcjxWJL6Gs3L5OeGeF//meQZ/DH9u8DnhL4SB7415ed4Kysc/CreZJSIiIhIZ1MdGREREPEOBjYiIiHiGAhsRERHxDAU2IiIi4hkKbERERMQzFNiIiIiIZyiwEREREc9QYCMiIiKeocBGREREPEOBjYiIiHiGAhsRiWoffvghEhISsH///qzr7rzzTnTs2BHJyclhPTYRKXmaFSUiUY0/wjp16oR+/frhlVdewdNPP4133nkHCxcuRP369cN9eCJSwkqV9CcUESlKPp8Pzz//PEaNGoU6deqY4GbevHkKakRilFZsRMQTunTpgnXr1mHatGno379/uA9HRMJEOTYiEvWmTp2KjRs3IiMjA7Vr1w734YhIGGnFRkSi2vLlyzFgwACMHz8eEyZMQOXKlTFp0qRwH5aIhIlybEQkau3YsQMjRozA73//e9x8881o1qwZevXqZYIdbk2JSOzRio2IRKVjx46hd+/eZrXmjTfeyLqegQ63pLg9JSKxR4GNiIiIeIaSh0VERMQzFNiIiIiIZyiwEREREc9QYCMiIiKeocBGREREPEOBjYiIiHiGAhsRERHxDAU2IiIi4hkKbERERMQzFNiIiIiIZyiwEREREc9QYCMiIiLwiv8PqNNiCPGaCrYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.errorbar(x, y, ye, ls='', c='k')\n",
    "ax.plot(x, np.asarray([model.sample(x=x) for count in range(100)]).T, c='r', alpha=0.1);\n",
    "ax.set(xlabel='$x$', ylabel='$y$', title='Randomly Generated Data');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1299efe-5b11-423d-8e6f-7451c02e11db",
   "metadata": {},
   "source": [
    "Often we want to be able to save a model to share it with colleagues. `Model` objects are able to save and load models using `json` files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "1b179373-c6b5-4a08-a06c-aaa379b45433",
   "metadata": {},
   "outputs": [],
   "source": [
    "model.save('output.json')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3074fd5-1ac6-4a9b-9cfd-f690989fec5d",
   "metadata": {},
   "source": [
    "This saves a juman readable model. Below we show the first 20 lines of this file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "e9e07094-ba16-49d4-8e59-9add41b55fb3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{\n",
      "    \"object_type\": \"JointModel\",\n",
      "    \"initializing_kwargs\": {},\n",
      "    \"priors\": {\n",
      "        \"mean\": [\n",
      "            0,\n",
      "            0\n",
      "        ],\n",
      "        \"std\": [\n",
      "            Infinity,\n",
      "            Infinity\n",
      "        ]\n",
      "    },\n",
      "    \"posteriors\": {\n",
      "        \"mean\": [\n",
      "            2.361845212485081,\n",
      "            1.4445468615716888\n",
      "        ],\n",
      "        \"std\": [\n",
      "            0.03674280542968608,\n"
     ]
    }
   ],
   "source": [
    "!head -20 output.json"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f5107f6-992f-4101-9ced-2f601c471a8b",
   "metadata": {},
   "source": [
    "The model json can then be loaded back in to lamatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "db81b733-cf38-430b-8a42-8e094ef42858",
   "metadata": {},
   "outputs": [],
   "source": [
    "from lamatrix import load"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "627657e9-be6d-433d-b4c6-cd46e61a0889",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "JointModel\n",
       "\tPolynomial(x)[n, 1]\n",
       "\tConstant()[n, 1]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "load('output.json')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a25c695-fb79-4833-9f39-cbb83917b07e",
   "metadata": {},
   "source": [
    "Finally you may want to include the results of your model as a latex table in your work. You can do this by using the `to_latex` method, which will also return the equation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "1beb24eb-ff19-4353-bea5-f2aa98ec5821",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\\[f(\\mathbf{x}) = w_{0} \\mathbf{x}^{1} + w_{1} \\]\n",
      "\\begin{table}[h!]\n",
      "\\centering\n",
      "\\begin{tabular}{|c|c|c|}\n",
      "\\hline\n",
      "Coefficient & Posterior & Prior \\\\\\hline\n",
      "w & $2.36 \\pm 0.04$  & $0 \\pm \\infty$ \\\\\\hline\n",
      "w & $1.44 \\pm 0.02$  & $0 \\pm \\infty$ \\\\\\hline\n",
      "\\end{tabular}\n",
      "\\end{table}\n"
     ]
    }
   ],
   "source": [
    "print(model.to_latex())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "870ae1a5-af9e-4093-8573-36b1757da720",
   "metadata": {},
   "source": [
    "That concludes the basics of `lamatrix`. To see how to use the package for more advanced cases check out the next tutorial."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}