{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Testing your code\n",
    "_building confidence that it is doing what it is supposed to!_ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from matplotlib import rcParams\n",
    "rcParams[\"font.size\"] = 14"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mathematical identities & analytic solutions\n",
    "_We know the answer_\n",
    "\n",
    "For example, for a moving average calculation, we know that if we compute the moving average along a straight line, the average over each segment should be the same as if we evaluate the line at the centre of that segment. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from agu_oss import moving_average"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def line(x, slope=1, intercept=0):\n",
    "    return slope * x + intercept"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "window_width = 5\n",
    "n = 20\n",
    "x = np.linspace(0, 100, n)\n",
    "y = line(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "if np.mod(window_width, 2) == 1:\n",
    "    moving_average_x = x  \n",
    "elif np.mod(window_width, 2) == 0: \n",
    "    # average assigned at the centers\n",
    "    moving_average_x = x - np.diff(x[:2])/2 # assumes evenly spaced x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "moving_average_y = moving_average(y, window_width)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1180c9eb0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAArwklEQVR4nO3deXxU1fn48c+TBbKwhbKEBBCQLSwq2UgsaEUUKiCKUvpFBfxVKVJFWgREsIZqFcENFVCsIu60oLaiggtQlEIgCIIsAVkUAoSwJGA2ksz5/TFDCJMEEnIzN5l53q9XXsnce+be50ySZ87cexYxxqCUUsp3+NkdgFJKKc/SxK+UUj5GE79SSvkYTfxKKeVjNPErpZSPCbA7gIpo0qSJadOmjd1hKKVUrbJx48Zjxpim7ttrReJv06YNKSkpdoehlFK1ioj8VNZ2vdSjlFI+RhO/Ukr5GE38SinlYzTxK6WUj9HEr5RSPqZCvXpE5BrgISAGiADuNsa8WWK/AI8Bo4EwIBn4kzFmW4kydYFngP8DgoGvgbHGmINVrcSpU6c4evQoBQUFVT2UslFoaCgtW7bEz0/bI8q3fbwpjVnLUzmUmUtEo2Am9uvELT0iLTt+Rbtz1gN+AN5yfbmbBEwARgGpwF+BL0WkkzHmtKvMC8BgnIn/OPAcsFREYowxRZdagVOnTpGenk5kZCTBwcE434NUbeNwOEhLS+PYsWM0a9bM7nCUss3Hm9KY8uFWcgucaTEtM5cpH24FsCz5V6hpZYz5zBjziDFmMeAouc/V2h8PzDDGLDHG/ACMBOoDw11lGgJ/ACYaY740xnwH3AVcAfStSgWOHj1KZGQkISEhmvRrMT8/P5o3b05WVpbdoShlq1nLU4uT/lm5BUXMWp5q2Tms+EzdFggHvji7wRiTC6wGrnZtigEC3cocAHaUKHMeERktIikikpKRkVHuyQsKCggODq5qHVQNEBgYSGFhod1hKGWrQ5m5ldp+KaxI/OGu7+lu29NL7AsHioBjFyhzHmPMfGNMrDEmtmnTUiOOz6Mtfe+gv0fl6/Zm/EKgf9lpOaKRdQ1cK6dscF/KS8rY5q4iZZRSyqsVFjl47Zt9PP/VLvzEEOgvFBSdS43Bgf5M7NfJsvNZ0eI/4vru3nJvxrlPAUcAf6DJBcoopZTP2XYoi1vmruHpZTvp06kZqyf3YdbtVxLZKBgBIhsF89SQ7pb26rEi8e/DmdhvOLtBRIKA3sD/XJs2AgVuZVoCUSXKqIt48803qVevXqWek5SURLdu3SyPRURYvHix5cdVylfkFRQxa/lObn55DUey8pl3RzSv3BVDs/pB3NIjkjUP92HfjAGsebiPpUkfKt6Pvx7Q3vXQD2gtIlcBJ4wxP4vIC8BUEdkJ7AKmAb8A7wEYY7JE5HVglogc5Vx3zi3AV9ZVx7sNGzaMm266qVLPeeihh3jggQeqKSKl1KXY+NMJJi3ewp6MbG6LbsmjA6NoFFLHY+ev6DX+WGBlicfTXV8Lcfbdn4lzUNYczg3gurFEH36APwOFwCLODeAaUZU+/Fap7sESVjjbe6myPZjq1atX6U8JSqnqkZ1fyKzlqSxcu5+IhsEs/H/xXNvxwp1XqkNF+/GvMsZIGV+jXPuNMSbJGNPCGBNkjLnW1Z+/5DHyjDEPGGN+ZYwJMcYMcnXptNXZwRJpmbkYzg2W+HhTWrWeNz8/n/Hjx9O8eXOCgoJISEjg22+/BWDVqlWICJ999hnx8fHUqVOH5cuXl3mp56mnnqJ58+bUq1ePESNGMH36dEouWuN+qWfUqFEMHDiQ2bNnExkZSVhYGHfffTc5OTnFZZYtW0bv3r0JCwujcePG9OvXjx07dlTr66GUt1u9K4Mbn1/NwrX7GZnYhuV/vsaWpA+1ZCGWypj+yTa2HzpV4fKbfs7kTNF5Y9LILShi0uItvL/+5wodo0tEAx4b1LVScU6aNIl//vOfvPHGG7Rr147nnnuO/v37s3v37uIykydP5tlnn6V9+/bUr1+fTz/99LxjfPDBB0yfPp2XX36Za665hiVLljBjxgzCwsIueO5vvvmGFi1a8NVXX3HgwAF+97vf0bFjR6ZMmQJAdnY248eP54orriA3N5cnnniCQYMGsX37durU8dzHUaVqq5JXEcIbBtEqLJj1+0/Srmko//pjIrFtGtsan9cl/spyT/oX226F7Oxs5s2bxz/+8Q8GDBgAwCuvvMKKFSuYM2cOffs6BzMnJSVx4403lnuc2bNnM2rUKO655x4ApkyZwsqVK9m1a9cFz9+gQQPmzZtHQEAAUVFRDB06lK+//ro48d92223nlV+wYAENGjRg/fr19OrV65LrrZQvcJ9y4XBWHoez8rghqhkvDY8mKNDf5gi9MPFXtuX96xkrSCtjRFxko2AW/THRqrDOs2fPHgoKCvj1r39dvM3f35/ExES2b99enPhjY2MveJydO3dy7733nretZ8+eF038Xbp0ISDg3K8+IiKC5OTk8+J79NFHSU5OJiMjA4fDgcPh4OefK/YJSClfVtaUCwDbD5+uEUkfdFpmJvbrRLDbL8PqwRLujHEOzChrpGrJbaGhoRc91qWMdg0MDCx1DIfj3CecQYMGkZGRwauvvkpycjKbNm0iICCAM2fOVPpcSvkSY0yZDUmwdsqFqvL5xH9Lj0ieGtK9WgdLuGvfvj116tQpvpkLUFRUxNq1a+nSpUuFj9O5c2fWr19/3jb3x5V1/PhxduzYwSOPPELfvn2Jiori9OnTOoeOUhdx4EQOI94o///PyikXqsrrLvVcilt6RHq0+2ZoaCj33XcfDz/8ME2aNKFt27Y8//zzpKenM3bsWFJTKzYL34MPPsjdd99NXFwcvXv35qOPPiI5OfmiN3cvJCwsjCZNmvDaa6/RqlUr0tLSmDhx4nmXhpRS5zgchrfW7mfm8lQEuD06kk+3Hia34Nyn6Oq+ilBZ+t9sk6effhqAu+++m8zMTHr06MGyZcto0aJFhRP/73//e/bu3cvDDz9MTk4OQ4YMYcyYMfz73/++5Lj8/PxYtGgR48aNo1u3brRv355nn3221A1fpRT8ePQ0k5dsZeNPJ7m2Y1P+fms3WoaF0KtD0xo9NkjOXm+uyWJjY01KSkqZ+3bs2EFUVJSHI6q5br31VgoLC/nkk0/sDuWS6O9T1QYFRQ7mr97L7K92E1LXn78O7MKtPSJr3AyzIrLRGFOql4i2+GuxnJwc5s2bR//+/QkICGDJkiX8+9//ZsmSJXaHppTX+iEti0mLt7D98CkGXNGCpEFdaVq/rt1hVYom/lpMRPj888958sknyc3NpUOHDrz99tvceuutdoemlNfJKyhi9te7mb96L41D6/DqXTH061rmciI1nib+Wiw4OJivvtI57pSy1KwOkH201OYcGjEvby7DYlvxyE1RNAwJLOPJtYMmfqWUKqmMpA/QmEze+UNPenVwX1ak9vH5fvxKKVVR3pD0QRO/Ukr5HE38SikFGIeDjZ+9YXcYHqHX+JVSPi/j0H4OvjOWmJw1dofiEdriV0r5LONwsH7JC9Sdn0hU9nrWXf4gJrScxVFCm3k2uGqkLX7lUYsXL2bo0KHUhhHjyrul7d3ByUVjiM/fzPY63an/u3kktO8O/M3u0KqdJn6llE8pKixkwz+f4orUl2iIH8ldpxF321/w868Zc+V7gib+cgZrENoMJu4uvb2WKygoKDUfv1K+Yv+OFPKXjCWhMJXvQ+JpPnwePVu1tzssj9Nr/OUM1ih3u4UutKh5YmIiEyZMOK/8qVOnCA4O5qOPPgLgzJkzTJ48mZYtWxIaGkpcXBzLly8vLl/eou179uxh8ODBhIeHExoaSnR0NEuXLj3vXOnp6dx8880EBwdz2WWXsWDBArp160ZSUlJxmaysLEaPHk2zZs2oX78+1157Le6T6b311ltcdtllhISEMHDgQNLT0618CZUqbVYHSGpY6uvM38KJ+OBGmhUeIiX6aa6YuJxwH0z64I0t/s8fhiNbrTnWggEVKxfeHX47o9KHv9Ci5nfeeSdPPvkks2bNws/P+f68ZMkSgoODi9fpvfvuu9mzZw/vvfceLVu25LPPPmPQoEFs2LCBK6+8svg87ou2Hzp0iN/+9rc88cQTBAcHs2jRIoYMGcKWLVvo3LkzACNHjuTw4cOsWLGC4OBgJkyYwE8//VR8TGMMAwYMoGHDhixdupTGjRuzcOFC+vTpQ2pqKi1atCA5OZlRo0bx+OOPM3ToUFauXMkjjzxS6ddJqUopp9FWx5HLxgZ9aHPnS8Q2b+nhoGoW75uWubKJ/6dvy993WQUXFr/ExO8uOzubBg0a8N///peoqChatGjB559/zvXXXw9A3759ufzyy3n11VfZs2cPHTp0YP/+/bRu3br4GLfccgsRERHMnTuXVatWcd1117F48eKLzqefkJDAwIEDmTZtGqmpqXTu3Jm1a9eSkJAAwIEDB2jTpg2PPvooSUlJrFixgptvvpmMjAyCg8+tLHTVVVcxfPhwJk2axPDhw8nIyODLL78s3n/PPffw+uuvl3tzV6dlVlWW1PAC+7I8F0cN4DvTMlc2AV/oj+TuT6sWy0VcaFHzXr160a9fP959912uv/56Dh8+zMqVK3nssccA+O677zDGlFqqMT8/nz59+py3zX3R9uzsbKZPn87SpUs5fPgwBQUF5OXlccUVVwDORdz9/PzOe16rVq2IiIgofrxx40ZycnJo2vT8rm95eXns2bMHcCbxQYMGnbc/MTGR119//VJeLqWURbwv8dcigwYNIjIykldffZXIyEgCAgLo0qVL8aLmd955J6NHj2bu3Lm8//77tGrVil69nJ9CHA4HIsKGDRtK3awt2QKH0ou2P/TQQyxbtoxnnnmGDh06EBISwogRI4rPW5FPgQ6Hg+bNm/PNN9+U2tegQYMKH0cpK53KPE4Du4OoBTTxhzYrv1dPNTq7qPmcOXO47rrrAGcrvuSi5oMHD2b06NEsXbqUd999lzvuuKN4hZ8ePXpgjOHIkSPFz6+ob7/9lhEjRhRf/jnbSu/YsSMAUVFROBwONm7cSM+ePQE4ePAghw4dKj5GdHQ06enp+Pn50a5duzLP06VLF9atW3feNvfHSlnl+xUf0GL1FE38FaCJ36YumxVZ1DwoKIghQ4bwxBNP8P333/POO+8U7+vYsSN33HEHo0aN4tlnnyU6OpoTJ06watUq2rVrx5AhQ8o9d8eOHfnoo48YPHgwgYGBTJ8+nby8vOL9nTp1ol+/fowZM4Z58+YRFBTExIkTCQkJKX7j6du3L7/+9a8ZPHgwM2fOpHPnzhw5coRly5bRt29fevfuzbhx47j66qt56qmnuP3221m1alVxjySlrHLiaBp73xlH7Kmv2O/XmsIABwFnMksX9KKRt1VmjKnxXzExMaY827dvL3dfTff111+brl27mrp165quXbuaZcuWmdDQULNgwYLzygAmOjq61PPPnDljHnvsMdO2bVsTGBhomjdvbgYNGmRSUlKMMcasXLnSACYjI+O85+3fv99cf/31JiQkxERGRppZs2aZAQMGmJEjRxaXOXz4sBk4cKCpW7euadWqlVmwYIFp166dmTFjRnGZU6dOmXHjxpnIyEgTGBhoWrZsaYYNG2Z+/PHH4jJvvPGGadWqlQkKCjL9+/c3L730knH+2ZWtNv8+lWc5iorMhqXzzYnHWpr8v4aZ/73+kMnPy7U7rBoFSDFl5FRLevWIiD+QBNwJtAAOA+8CScaYQlcZAR4DRgNhQDLwJ2PMtosdXxdbt9+xY8eIiIjg/fffv2gPoarQ36eqiKNp+zj07n1clbOWXQEdCRwyl7Zd4uwOq8ap7l49k4E/ASOBrcAVwEIgH3jcVWYSMAEYBaQCfwW+FJFOxpjTFsWhLLJixQpOnz5N9+7dOXr0KFOnTqVJkyb079/f7tCUDzMOBxs+fIGorTPpRBHrOv6FuGFT8Q/Qq9aVYdWrdTXwiTHmE9fj/SLyH6AnFLf2xwMzjDFLXNtGAkeB4cCrFsWhLFJQUMC0adPYu3cvISEh9OzZk9WrV5fqIaSUp6Tt3UbmB/cRf+Z7ttW9kka/n0dCu652h1UrWZX4vwXGikhnY8xOEekC9AGecu1vC4QDX5x9gjEmV0RW43zT0MRfw/Tr149+/frZHYbyReXMnxVhoAHBrO/+GHFDxiN+OuPMpbIq8T8N1Ae2i0iR67h/N8bMde0Pd313n6glHYgs64AiMhrn/YDzRqYqpbxcOVMuiEDuvWuJj2zr4YC8j1VvmcOAETgv20S7fh4rIn9wK+d+J1nK2OYsaMx8Y0ysMSbWfXRoGWUvKWhVs+jvUV1MM036lrCqxT8LeMYY84Hr8VYRuQyYArwOHHFtDwcOlHheM0p/CqiUwMBAcnNzCQkJqcphVA1QUFBw3jgGpVT1sKrFHwIUuW0rKnH8fTiT/w1nd4pIENAb+F9VTtysWTPS0tLIycnRFmMt5nA4SE9Pp2HDC8ydpLxabvZp1s0bY3cYPsGq5tUnwMMisg/YBvQA/gK8Bc7ROiLyAjBVRHYCu4BpwC/Ae1U58dl5YQ4dOkRBQUFVDqVsFhoaSpMmTewOQ9ngh2//Q9jXD5FgdL0GT7Aq8T+As7/+XJyXbw4Dr3H+4pUzgWBgDucGcN1oRR/+Bg0aFL8BKKVqj6yTx0h960HiTy7loLRg243v03XNg7bMn+VLav18/Eqp2mnzl+8RsWYqvzInWR9xBz3uepqgkHp2h+VVfGc+fqVUjXY8/SD737mfmNMr2efXhqxBC0nscY3dYfkUTfxKKY8wDgcbl87n8u+eoLvJZW2bMcQMn06dukF2h+ZzNPErpaxVzsjbQgKIpZDUgM4E3TaXxKgYG4JToIlfKWW1ckbeBlLIuo4Tifvdwzqpms301VdKeUzC8Gl2h6CwbgCXUkqpWkITv1LKMnu26prKtYFe6lFKVVl+Xg7fvTOV2AMLnVMvqhpNW/xKqSrZueErjsyMJ/HgG2xu1BdHcDnTbujI2xpDW/xKqUuSfTqTrW9PJD79XxyVJmy59nXirrvd7rBUBWjiV0pV2tbVH/GrlZNIMEdJbjqEriOeI7xBmN1hqQrSxK+UqrCsExmkvjWO+MzPOCARbO+/iJ4J/e0OS1WSJn6lVIV8t/xtWq+dRrQ5xdrIEfS4awatgkPtDktdAk38SqnzlTPlQjSwx78dmTe/S+KVvTwfl7KMJn6l1PnKmXIBoPXkdQTWqevBYFR10O6cSqkK06TvHTTxK6WKOYrcl85W3kgTv1IKgJ93bSZ1Rm+7w1AeoNf4lfJxBWfySXn/b0TvfZU8qWN3OMoDtMWvlA/78fs1/PR0Ion7XmZbvUQKxqwrf2oFnXLBa2iLXykflJebzaZ3HiHu4FtkSgO+S5hNdP9Rzp0Td9sam6p+mviV8jE7k78gZPl4Eh1prA+7iU4jXiS6cVO7w1IepIlfKR/xy6mTbHtrAnEZH5IuTdh63QLirx1id1jKBpr4lfI25Yy8DUaIM7C+2e10H/EMLeo38nxsqkbQxK+Utyln5K0/hp0DFpMQf4OHA1I1jfbqUcqHdNakr9DEr5RSPkcTv1JewjgcrP/oRbvDULWAXuNXygsc2p/K8ffHEJ//nd2hqFrAsha/iLQQkYUikiEieSKyXUSuLbFfRCRJRA6JSK6IrBKRrladXylfVFRYyLr3/06jBb1pl7ed5KhHMDryVl2EJS1+EWkErAG+BQYAGUA7oGT3gknABGAUkAr8FfhSRDoZY05bEYdSvuSnnd+Ru3gsCYU7+D44jubD59GzdQdgst2hqRrOqks9k4DDxpgRJbbtO/uDiAgwHphhjFni2jYS5xvDcOBVi+JQyusVnMkn5b3HiNn3GjkSREr0DGIG/hHx01t2qmKs+ku5BUgWkUUiclRENovI/a6ED9AWCAe+OPsEY0wusBq4uqwDishoEUkRkZSMjAyLwlSqdtu9+RsOzIgncf88ttbvRdF9ycTefJ8mfVUpVrX42wFjgeeBGcBVwEuufS/jTPoA6W7PSwciyzqgMWY+MB8gNjbWWBSnUjVfOSNvz/gF0bboDCelIZuunkPMjXfaEJzyBlYlfj8gxRgzxfV4k4h0AP6EM/Gf5Z7ApYxtSvm2ckbe1nHksb7xADqNeJEeYU08HJTyJlZ9PjwMbHfbtgNo7fr5iOt7uFuZZpT+FKCUKkf8g+/RUJO+qiKrEv8aoJPbto7AT66f9+FM/sXjxUUkCOgN/M+iGJRSSlWAVYn/eSBBRKaKSHsRGQqMA+YAGGMM8ALwsIgMEZFuwJvAL8B7FsWgVK13MuOw3SEoH2DJNX5jzAYRuQV4EngU+Nn1fW6JYjOBYJxvBmFAMnCj9uFXyjndwnfLFtB2fZLdoSgfYNmUDcaYT4FPL7DfAEmuL6WUS8ah/Rx85z5icv7H7oAONPD3JyD/ZOmCOvJWWUTn6lHKJsbhYMNHs+m8dSZRpoB1HcYTO2wqAYF17A5NeTlN/ErZIG3vDk4uGkN8/ma21elOw9+9QkL7bnaHpXyEJn6lPKiosJANi57kyl0v0RB/krs9StyQP+Pn7293aMqHaOJXykP270ghf8lYEgpT+T6kJ+F3zKNny8vtDkv5IE38SlmtnCkXLjOQKfVJiZlJzIB7dX4dZRtN/EpZrZwpF0TAjE0mtlmZ01Mp5THa5FDKgxpr0lc1gCZ+pZTyMZr4lbLIqczjJL804uIFlbKZXuNXygKbv/6AiG+mEGtOOicbV6oG0xa/UlVw4mgaKc8O4apv/ki2X332DP64/KkVdMoFVUNoi1+pS2AcDjZ+9g8uT3mcK0w2ay8bTcwdj1OnbhBE77Y7PKUuSBO/UpWUfnAPh98dS2zuOnYFdCRzyFwSu8TZHZZSFaaJX6kKchQVseHDF+jywyw6UcS6jhOIG/YI/gH6b6RqF/2LVcpdOSNviwigJ4Vsq3sljX4/j4R2XW0ITqmq08SvlLtyRt4GUsj67knE3fqgTregajVN/EpVQvxtf7Y7BKWqTJstSinlYzTxK1VCasoKu0NQqtrppR6lgJxfstjy9iTijyzSkbfK62mLX/m8H779DyefjSMh/QM2NBmMI6RJ2QV15K3yEtriVz4r6+QxUt96kPiTSzkoLdh24/v0vPomu8NSqtpp4lc+afOX7xGxZiox5iRrI+6kx11P0zKknt1hKeURmviVTzmefpD979xPzOmV7PNrQ9aghST2uMbusJTyKE38yicYh4ONS+dz+XdP0N3ksrbNGGKGT3dOqqaUj9HEr7xPGVMuCBALpAZ0Jui2uSRGxdgSmlI1gSZ+5X3KmXIBoP3Da3RSNeXztDun8ima9JWqpsQvIo+IiBGRl0tsExFJEpFDIpIrIqtERKc3VJYqLDhjdwhK1XiWJ34RSQDuBba47ZoETAAeAOKAo8CXIlLf6hiUb9qzdR37ZiTaHYZSNZ6liV9EGgLvAn8ATpbYLsB4YIYxZokx5gdgJFAfGG5lDMr35OflsPYff6b14ptoXJRhdzhK1XhWt/jnA4uNMe4zXbUFwoEvzm4wxuQCq4GrLY5B+ZCdG77iyMx4Eg++weZGffG/f70udq7URVh2p0tE7gXaA3eVsTvc9T3dbXs6EFnO8UYDowFat25tUZTKW2SfzmTr2xOJT/8XR6UJW659nbjrbnfunKiLnSt1IZYkfhHpBDwJ9DbGXOjumnF/ahnbnAWNmY/zEwSxsbFlllG+aevqj/jVykkkmKMkNx1C1xHPEd4gzO6wlKo1rGrxJwJNgB+cl/MB8AeuEZExwNneO+HAgRLPa0bpTwFKlSnrRAapb40jPvMzDkgE2/svomdCf7vDUqrWsSrxfwykuG1bAOzG+UlgF3AEuAHYACAiQUBvYKJFMShvUc5i5/URoo2wNnIkPe56ilbBoTYEp1TtZ0niN8ZkApklt4lINnDC1YMHEXkBmCoiO3G+EUwDfgHesyIG5UXKGXnrh2HvkKUkXtnLwwEp5V08OYxxJhAMzAHCgGTgRmPMaQ/GoGq59pr0laqyakv8xpjfuD02QJLrSymllE10rh5VYziKikheNMPuMJTyejpjlaoRft61mex/jaVnwTa7Q1HK62mLX9mq4Ew+axc+QvN3+xJRsJ8NV/0doyNvlapW2uJXtvnx+zXwnwdILNrDd/WvofWdc4gLbw3cb3doSnk1TfzK4/Jys9n09hTi0t4mSxqwKfFFovuNtDsspXyGJn7lUTuSlxO67M8kmjTWh91EpxEv0qNxU7vDUsqnaOJX1itn5O0ZqUuUyeeQNGNrnwXEXzPEhuCUUpr4lfXKGXlbx+SzrulQuo94hoj6jTwbk1KqmCZ+5VEJf/qH3SEo5fO0O6dSSvkYTfzKUscO/WR3CEqpi9DEryxhHA7Wf/gidebrYudK1XR6jV9V2aF9Ozn+wRji8zexPbAbnfwP4Z93onRBHXmrVI2giV9dsqLCQjb882muSJ1NQ4TkrlOJu20Cfv7+doemlLoATfzqkuzf+R15i8eSULiD74PjaD58Hj1bd7A7LKVUBWjiV5VScCaflHcfI2b/a+RIECnRM4gZ+EfET28XKVVbaOJXFbZ78zf4/+d+Eh372Vj/N7S5aw6xzVvaHZZSqpI08avSyplyob2BYxLGpqvnEHPjnTYEppSygiZ+VVo5Uy6IQJ0HU+gR1sTDASmlrKQXZlWlNNSkr1Stp4lfKaV8jCZ+VSzz2BE2PD/U7jCUUtVMr/ErjMPBd8sW0HZ9EleZbBC7I1JKVSdt8fu4jEP72fzMAGLW/4XjAc05MPTz8qdW0CkXlPIK2uL3UcbhYMNHL9J569NEmQLWdRhP7LCpBATWgW677Q5PKVWNNPH7oLS9Ozi5aAzx+ZvZVqc7DX/3Cgntu9kdllLKQzTx+5CiwkI2LHqSK3e9REP8Se72V+KGjNdJ1ZTyMZr4vVE5I28N/iRQxPchPQm/Yx49W15uQ3BKKbtp4vdG5Yy8DaCIlNhZxNx0j06qppQPs+S/X0SmiMgGETklIhki8omIdHMrIyKSJCKHRCRXRFaJSFcrzq8qLnbgaE36Svk4qzLAb4C5wNVAH6AQ+EpEGpcoMwmYADwAxAFHgS9FpL5FMSillKoASy71GGP6lXwsIncBWcCvgU9ERIDxwAxjzBJXmZE4k/9w4FUr4lCwbc2n6McopdSFVNdn/vquY590PW4LhANfnC1gjMkFVuP8lFCKiIwWkRQRScnIyKimML3HqczjJL94F12/HG53KEqpGq66Ev9sYDOw1vU43PU93a1ceol95zHGzDfGxBpjYps2bVotQXqLzV+9T94LscQe/4R14XdgQsp5vXTkrVKKaujVIyLPAb2AXsaYIrfdxr14GdtUBZ04msbetx8g9vTX7PNrQ+bAN0mIvtbusJRSNZyliV9Engd+D1xnjNlbYtcR1/dw4ECJ7c0o/SlAXYRxONj46WtcvvFxrjA5rL3sj8Tc8Tfq1A2yOzSlVC1gWeIXkdk4k/5vjDE73Xbvw5n8bwA2uMoHAb2BiVbF4AvSD+7h8Ltjic1dR2pAJ+reNpfEqFi7w1JK1SKWJH4RmQPcBdwCnBSRs9ftfzHG/GKMMSLyAjBVRHYCu4BpwC/Ae1bE4O0cRUVs+PB5uv7wDJ0oYl3HCcQNewT/AB2Dp5SqHKuyxljX96/dtk8Hklw/zwSCgTlAGJAM3GiMOW1RDN6jjCkX/ICewA9BVxE27BUS2kXZEppSqvazqh//RZfuMMYYnG8CSVac06uVM+UCQNfJK3XkrVKqSjSD1DKa9JVSVaVZpIbJz8uxOwSllJfTxF+D7Ez5msMz4+0OQynl5bRLSA2Q80sWW96aSHz6P8k4b147pZSynrb4bfbDN/8m89k4Eo4uYkOTwYT8OUUXO1dKVStt8dsk6+QxUt8aR/zJTzkgEWzv9wE9E3/r3DlRFztXSlUfTfw22PTFO7T83zSiTRZrI0bQ466naBVSz+6wlFI+QhO/Bx07coCf3rmfmF9Wsce/LZmD3ibxqt52h6WU8jGa+KtDOYudNwYamADWtr2P2OHTCaxT1/OxKaV8nib+6lDOyFs/4PD/fUli52jPxqOUUiVorx4Pu0yTvlLKZpr4lVLKx2jit1BhkYN5q/bYHYZSSl2QXuO3yLZDWUxesoUf0k5xny6EpZSqwbTFX0V5BUXMWr6Tm19ew5GsfObdEa0jb5VSNZq2+Ktg408nmLR4C3sysrktuiWPDoyiUUgd6K4jb5VSNZcm/kuQnV/IrOWpLFy7n4iGwSz8f/Fc27Gp3WEppVSFaOKvpNW7Mpjy4VbSMnMZmXgZE/t3pl5dfRmVUrWHZqwKysop4PFPt7N440HaNQ3lX2MSiWujUygrpWofTfzl+HhTGrOWp3IoM5ew0EAKCh3kFDgY+5vLGXd9B4IC/e0OUSmlLokm/jJ8vCmNKR9uJbegCIAT2QUI8JcbO/JAnw72BqeUUlWk3TnLMGv5zuKkf5YBPlh/wJ6AlFLKQpr43Rw4kUNaZl6Z+w5l5no4GqWUsp5e6nFxOAxvrd3PzOWpCM4WvruIRsGeDksppSyniR/48egvPLxkCyk/neSajk35TccmzFq+67zLPcGB/kzs18nGKJVSyho+nfgLihzMX72X2V/tJriOP88MvZLboiMRERqH1i3u1RPRKJiJ/TpxS49Iu0NWSqkq89nE/0NaFpMWb2H74VPc1D2cpJu70qz+udnVbukRqYleKeWVfC7x5xUUMfvr3cxfvZfGoXV45c5o+ndrYXdYSinlMR5P/CIyFpgItAC2AeONMd944twb9p9g8uIt7D2WzdCYlkwb0IWGIYGeOLVSStUYHk38IjIMmA2MBb51ff9cRLoYY3628lwlR96GNwzi8qahfPvjcVqGBfP2H+Lp3UEnVVNK+SZPt/j/ArxpjHnN9fgBEekP3AdMseok7iNvD2flcTgrj94dmvDKnTGE6qRqSikf5rEBXCJSB4gBvnDb9QVwtZXnmrU8tdTIW4C9Gdma9JVSPs+TI3ebAP5Autv2dCDcvbCIjBaRFBFJycjIqNSJyhthqyNvlVLKnikb3AfFljlQ1hgz3xgTa4yJbdq0ctfjyxthqyNvlVLKs4n/GFBE6dZ9M0p/CqiSif06Eew2bbKOvFVKKSePJX5jzBlgI3CD264bgP9Zea5bekTy1JDuRDYKRoDIRsE8NaS7DshSSik836vnOeBtEVkPrAHGABHAK1afSEfeKqVU2Tya+I0xi0TkV8A0nAO4fgBuMsb85Mk4lFLKl3m8b6MxZi4w19PnVUop5aQLsSillI/RxK+UUj5GE79SSvkYMaasRQZrFhHJAC71BnATnGMIfInW2Tdonb1fVet7mTGm1AjYWpH4q0JEUowxsXbH4UlaZ9+gdfZ+1VVfvdSjlFI+RhO/Ukr5GF9I/PPtDsAGWmffoHX2ftVSX6+/xq+UUup8vtDiV0opVYImfqWU8jGa+JVSysd4beIXkbEisk9E8kRko4j0tjsmq4jIFBHZICKnRCRDRD4RkW5uZUREkkTkkIjkisgqEelqV8xWE5FHRMSIyMsltnldnUWkhYgsdP2e80Rku4hcW2K/V9VZRPxF5PES/7v7ROQJEQkoUaZW11lErhGR/4hImutveJTb/ovWT0TqishLInJMRLJdx2tZ0Ri8MvGLyDBgNvAk0APnQi+fi0hrWwOzzm9wznB6NdAHKAS+EpHGJcpMAiYADwBxwFHgSxGp79lQrSciCcC9wBa3XV5VZxFphHPdCgEGAFE463a0RDGvqjMwGfgTMA7oDDzoejylRJnaXud6OKekfxAoayHwitTvBeA24P+A3kADYKmInL/0YHmMMV73BSQDr7lt2w08ZXds1VTfejiXtRzkeizAYWBqiTLBwGngj3bHW8W6NgT24HzDWwW87K11xtlwWXOB/d5Y56XAQrdtC4Gl3lhn4BdgVGV+p67/gTPAHSXKtAIcQL+KnNfrWvwiUgeIAb5w2/UFzhayN6qP89PbSdfjtjjXNi5+DYwxucBqav9rMB9YbIxZ4bbdG+t8C5AsIotE5KiIbBaR+0VEXPu9sc7fAteJSGcAEemC803+M9d+b6xzSRWpXwwQ6FbmALCDCr4GHl+IxQOaAP6UXsA9Hejr+XA8YjawGVjrenx2QfuyXoNaux6liNwLtAfuKmO3N9a5HTAWeB6YAVwFvOTa9zLeWeencTZktotIEc4c9XfjXMAJvLPOJVWkfuE4P+G7T96WXuL5F+SNif8s95FpUsa2Wk9EngN6Ab2MMUVuu73mNRCRTjgvffQ2xpy5QFGvqTPOT3Epxpiz17c3iUgHnNe8Xy5RzpvqPAwYAQwHtuF8s5stIvuMMa+XKOdNdS7LpdSvwq+B113qwfkuWETpd75mlH4XrdVE5HmcN3f6GGP2lth1xPXdm16DRJyf5n4QkUIRKQSuBca6fj7uKudNdT4MbHfbtgM420nBG3/Ps4BnjDEfGGO2GmPeBp7j3M1db6xzSRWp3xGcVzWaXKDMBXld4ne1BjcCN7jtugFn7x6vICKzcbaK+hhjdrrt3ofzj+OGEuWDcN79r62vwcdAd5wtwLNfKcAHrp934X11XgN0ctvWkXNrU3jj7zkEZ8OtpCLO5SpvrHNJFanfRqDArUxLnL2+KvYa2H1Xu5rulA/Dedf7HteLMRvn3fPL7I7NovrNAU7hvOkVXuKrXokyk11lhgDdcCbIQ0B9u+O38HVYhatXjzfWGWdXvgJgKs57G0OBLOBPXlznN4GDOLuvtgFuBTKAZ72lzjh74V3l+soB/ur6uXVF6wfMA9Jw3rfsAazEeZ/Pv0Ix2P0iVOOLOxbYD+TjfIe8xu6YLKybKecrqUQZAZJwXi7IA/4LdLM7dotfB/fE73V1diXA71312YWzf7t4a51x3th9AeenmlxgL857O0HeUmec43DK+v99s6L1A4Jw3ug/7nrz+ARoVdEYdHZOpZTyMV53jV8ppdSFaeJXSikfo4lfKaV8jCZ+pZTyMZr4lVLKx2jiV0opH6OJXymlfIwmfqWU8jH/H2TFTf9pyjlqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(x, y, '-o', label=\"original\")\n",
    "ax.plot(moving_average_x, moving_average_y, '-s', label=\"averaged\")\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_test = line(moving_average_x)\n",
    "\n",
    "start_ind = window_width // 2\n",
    "end_ind = len(y) - (window_width - start_ind)\n",
    "passed = np.allclose(\n",
    "    moving_average_y[start_ind:end_ind],\n",
    "    y_test[start_ind:end_ind]\n",
    ")\n",
    "assert(passed)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Code comparisons\n",
    "_Compare against an independent implementation_\n",
    "\n",
    "We can also implement a moving average using pandas. Lets compare with that. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame({\n",
    "    \"x\": x,\n",
    "    \"y\": y\n",
    "})\n",
    "df.set_index(\"x\", inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "pd_averaged = df[\"y\"].rolling(window=window_width).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pd_averaged"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "the pandas approach puts the value at the end of the window, whereas our implementation puts it at the center. So we need to account for that in our comparison "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "passed_comparison = np.allclose(\n",
    "    pd_averaged.values[window_width-1:-1],\n",
    "    moving_average_y[start_ind:end_ind]\n",
    ")\n",
    "assert(passed_comparison)\n",
    "passed_comparison"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checking Behaviour: Convergence \n",
    "_We expect certain behaviours of the code_\n",
    "\n",
    "Even if we don't know the real solution, we know certain things about how the solution to converge.\n",
    "\n",
    "For example: computing derivatives. Consider the Taylor expansion of a function \n",
    "\n",
    "$$\n",
    "f(x + h) \\simeq f(x) + \\frac{\\partial f(x)}{\\partial x} h + \\mathcal{O}(h^2)\n",
    "$$\n",
    "\n",
    "If we approximate the solution to $f(x + h)$ by $f(x)$, the we expect that as we decrease $h$, then the error $\\varepsilon$ should decrease as $\\mathcal{O}(h)$:\n",
    "\n",
    "$$\n",
    "\\varepsilon(h) = |f(x + h) - f(x)|\n",
    "$$\n",
    "\n",
    "e.g. for small enough $h$, if we decrease $h$ by a factor of 2, $\\varepsilon$ should decrease by a factor of 2. \n",
    "\n",
    "If instead, we include gradient information in our approximation of $f(x + h)$, namely using $f(x) + \\frac{\\partial f(x)}{\\partial x} h$, then as we decrease h, the error $\\hat{\\varepsilon}$ should decrease as $\\mathcal{O}(h)$:\n",
    "\n",
    "$$\n",
    "\\hat{\\varepsilon}(h) = \\left|f(x + h) - \\left(f(x) + \\frac{\\partial f(x)}{\\partial x} h)\\right)\\right|\n",
    "$$\n",
    "\n",
    "e.g. for small enough $h$, if we decrease $h$ by a factor of 2, $\\hat{\\varepsilon}$ should decrease by a factor of 4. \n",
    "\n",
    "#### Why do we care about this behaviour?\n",
    "\n",
    "If we are using gradient-based optimization (e.g. in an inversion or a machine learning application), we need to provide methods for computing derivatives at any point. This test checks the behaviour of the derivative. \n",
    "\n",
    "**For reference:** https://doi.org/10.1137/1.9781611973808"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f(x):\n",
    "    \"\"\"\n",
    "    A decaying sinusoid\n",
    "    \"\"\"\n",
    "    return np.exp(-x**2/10) * np.sin(2*x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def f_deriv(x):\n",
    "    \"\"\"\n",
    "    The derivative of our decaying sinusoid\n",
    "    \"\"\"\n",
    "    return -x/5 * np.exp(-x**2/10) * np.sin(2*x) + 2*np.exp(-x**2/10) * np.cos(2*x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x11810f160>"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD8CAYAAAB9y7/cAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6XklEQVR4nO3deXzM1/rA8c/JZF+RTYIkdmLfgpRSpVpcWtUFXXRvXd1+3XVzu96ly+16W92rvaU7Sle0ai2KCLGTBCFIJLJv5/fHEdcuiZn5zkye9+s1r5HJd77nmUieOXO+5zxHaa0RQgjhebysDkAIIYRjSIIXQggPJQleCCE8lCR4IYTwUJLghRDCQ3lbHcCxIiIidEJCgtVhCCGEW1m1atUBrXXkiY+7VIJPSEhg5cqVVochhBBuRSmVfqrHZYhGCCE8lCR4IYTwUJLghRDCQ0mCF0IIDyUJXgghPJTdZ9EopR4BRgNtgVJgGfCI1jrV3m0JIdxPfn4+2dnZlJeXWx2KW/Dx8SEqKorQ0NBaP9cR0yQHAm8CKwAFPAX8opRK1FrnOKA9IYSbyM/PZ9++fTRp0oSAgACUUlaH5NK01hQXF7N7926AWid5uyd4rfXQY79WSl0L5AHnAbPt3R4Aa2dAeRG0GwHBJ831F0K4iOzsbJo0aUJgYKDVobgFpRSBgYE0adKEPXv21DrBO2MMPuRIO7mn+qZS6lal1Eql1Mr9+/fXrYXUr+C7e+DFNvDhCFj/TZ2DFUI4Tnl5OQEBAVaH4XYCAgLqNKTljAT/CrAGWHqqb2qtp2qte2qte0ZG1rH3PW4G3L4Y+t8Ph7Pgiwmw5LW6xiuEcCAZlqm9uv7MHFqqQCn1EtAP6Ke1rnRgQ9C4o7kNeBC+uhl+egwqy6D/fQ5rVgghXJnDErxS6mXgauACrfV2R7VzEpsPXP6euZ/3FFRVmqQvhBD1jEOGaJRSrwDjgEFa642OaOOMbN5w2dvQ+SpY8CykL3F6CEIIz1FVVcVtt91GeHg4Sil+/fVXcnNziY6OZtu2bTU6R2lpKXFxcU4tqGj3BK+UegO4ARgL5CqlGh+5Bdu7rTPyssGIlyGsGXx3L1SUObV5IYTnmDt3Lh988AGzZ88mKyuL5ORknnvuOYYNG0bLli1rdA4/Pz8eeOABHnroIQdH+z+O6MFPxMycmQdkHXO73wFtnZlvEAx7AfZvhCWvOr15IYRn2Lp1KzExMSQnJ9O4cWMqKip49913uemmm2p1nvHjx7No0SLWr1/voEiPZ/cEr7VWp7lNsXdbNdL2YkgcBQv/BTnOuxQghPAMEyZM4N577yUjIwOlFAkJCcydOxcvLy/OO++8o8c9/fTTNG7cmOzs7KOPjR07lu7du1NWZkYQGjVqxHnnncdnn33mlNhdasMPh7n4H7B1Psy5D6752sy6EUK4hL/NXs+GPflObTMxNpQn/9KhRse+8sorxMfH8/7777NixQpsNhvPPPMMPXr0OG764uTJk/npp5+48cYb+e677/j444+ZOXMmf/75J76+vkePS0pK4rfffrP7azqV+lFsLDQGLngEts2HjGVWRyOEcCNhYWGEhIRgs9lo3LgxkZGRpKenExMTc9xxNpuNTz75hEWLFvHggw8yadIkXnzxRdq1a3fccbGxsezcudMpsdePHjxAjxtg4QtmLD6+r9XRCCGOqGlP2pUUFxcTHR190uPx8fG88sorTJgwgeHDh3PHHXecdExAQADFxcXOCLOe9OABfAMh6RbYNBcObLE6GiGEG4uIiCA395TVV1i4cCE2m42MjAxKS0tP+n5OTg51XrVfS/UnwQP0ugVsflLGQAhxTrp168aGDRtOevzrr7/m008/Zf78+eTn5/PII4+cdExqairdu3d3Rpj1LMEHR0LXsbB2OhRkn/14IYQ4haFDh5KWlsbBgwePPrZnzx5uueUWnnvuOc4//3w++eQTXnvtNX7++efjnvv7779z8cUXOyXO+pXgAfpOMjVq/phqdSRCCDfVqVMnkpKSmD59OmDqtl9//fV069aNe++9F4B+/frx8MMPM2HChKNvBEuXLiUvL48xY8Y4JU6ltXZKQzXRs2dP7ZRlvJ+Ng4wlcO8GMzYvhHCKtLQ02rdvb3UYdvHDDz9w9913s2HDBmw2W42ec8UVV9CtWzcmT55c6/bO9LNTSq3SWvc88fH614MH6H0bFOeaC65CCFEHF198MX/961/ZtWtXjY4vLS2lS5cuR3v4zlA/E3xCfwhtAikzrI5ECOHG7rrrLuLj42t0rJ+fH4899phTNzypnwneyws6Xwlb58nFViGEx6qfCR6g89WgK2Hdl1ZHIoQQDlF/E3xUO4jpCinTrY5ECCEcov4meIAuV0PWWshOszoSIYSwu/qd4DuOAWUzC5+EEMLD1O8EHxwJrQbDui+gqsrqaIQQwq7qd4IHM5smfzdkLLU6EiGEsCtJ8K0vApuvLHoSQpzWqTbdPtGWLVuIjo4mLy+vRufMzs4mMjKyxgul6kISvH8oND8fNs4BFyrbIIRwHSduuj116lSmTJly3DGTJ09m4sSJhIWF1eicUVFRXHfddTz55JMOiNiQBA/Qdhjk7jCbcwshxAlO3HT72C34ADIzM/n222+54YYbanXeG264gU8//ZScnBx7hnuUJHgwCR5ML14IIY5xqk23TzRjxgw6depEXFzc0cduuukmOnTocHT3psrKSvr168eIESOOHtOxY0diY2P5+uuvHRJ7/dmy70xCY6BJD5Pgz7/f6miEqF++fxj2rnNum407wSV/r9Ghp9p0+4EHHjjumN9//52ePY8v5vjqq6/SrVs37r//ft544w2effZZtm7dSkpKynHHVW/CffPNN5/bazoFSfDV2g6D+U9DfpZJ+EIIwcmbbgN8+OGHxx2Tnp5O165dj3ssKCiI//73vyQnJxMeHs7zzz/PrFmziIqKOu642NhYVqxY4ZDYJcFXazfcJPhNc6HXTVZHI0T9UcOetCsrLi7G39//pMd79uzJo48+ypQpU5g4cSKXXHLJScc4chNuGYOvFtkOGrWQcXghRK2dbhNurTWLFi3CZrOxbds2TrXBkiM34ZYEX00pM0yzYyGU5FsdjRDCjZxuE+6XXnqJP//8k4ULF7Js2TJee+21k45x5CbckuCP1XYYVJXD9gVWRyKEcCNDhw5l2bJlVFRUHH1s7dq1PProo0ydOpXk5GT+85//8NBDD5Gamnr0mKKiIlatWuWwTbglwR+rWRL4hsA2SfBCiJobNmwYAQEB/PjjjwCUlJQwfvx4xo0bx+WXXw7A2LFjGTNmDOPHj6e0tBSAmTNnEhcXR//+/R0SV/3cdPtMPhsH+9bB3Slm2EYIYTeetOn2id566y2++OIL5s2bV+PnJCUlcc899zBu3LizHiubbttDywvgUAbkbLc6EiGEG7nlllu48MILa1WLZsyYMYwdO9ZhMUmCP1HLQeZ+23xr4xBCuBWbzcbkyZNrVYvmwQcfRDlwpEAS/IkatYAG8ZLghRBuTxL8iZSCVhea6ZKV5VZHI4QQdSYJ/lRaDoKyAtjlmOXDQtRnrjSxw13U9WcmCf5UEvqbvVplmEYIu/Lx8XHYsnxPVlxcjI+PT62fJwn+VAIaQNOekuCFsLOoqCh2795NUVGR9ORrQGtNUVERu3fvPqlIWU1IsbHTaTkIfv07FOVAYCOroxHCI4SGhgKwZ88eysvlGldN+Pj4EB0dffRnVxuS4E+n5SD49XnY8Rt0uMzqaITwGKGhoXVKVqL2ZIjmdGK7gW8w7FxsdSRCCFEnkuBPx+YDcX1g5+9WRyKEEHUiCf5MEvqZjbgL9lsdiRBC1Jok+DNJOFLhLX2RtXEIIUQdSII/k5guR8bhJcELIdyPQxK8Uup8pdQspdRupZRWSk1wRDsOd3QcXhK8EML9OKoHHwykAncD7r1sTcbhhRBuyiEJXms9V2s9WWv9JVDliDacRsbhhRBuSsbgzyamq4zDCyHckuUJXil1q1JqpVJq5f79LjgMYvOGuL6wQ+bDCyHci+UJXms9VWvdU2vdMzIy0upwTi2hHxzYBAXZVkcihBA1JrVoaqJ6HH7nIug4+qRva61J3Z3PzDW72ba/gJzCMg4WlhEW4EOnJmF0bBJG/9YRxIcHOTlwIUR9Jgm+JmK6gE8QZCw9LsGXV1YxbWk6n/2RwZbsAny9vWgTHUyjID+aRwRxsLCM71P3Mn1FJl4KLukUw8SBLekQW7M9G4UQ4lw4JMErpYKBVke+9ALilFJdgRytdYYj2nQom7epD5+x9OhDqbvzePDLFDZk5dMjviHPXdaJ4Z1iCAs8vii/1pqMnCKmr8hk2tJ05qRkcXGHxjw3uhONgnyd/UqEEPWIckTRfaXUQGDBKb71kdZ6wume17NnT71y5Uq7x2MXC56Hhf+k8oEd/HvRPt78dRuNgnx55tKODO3QuEanyCsu56MlO3l9/lYaBvnwytXd6NMi3MGBCwDysyB9MWSngW8QBDSE4Chofj74hVgdnRDnRCm1Smvd88THHdKD11r/CihHnNsycX1AV/HWp9N5bVszxvRoyuPDE0/qsZ9JWIAPd13YmkHtorjrs9WMe2cZ9wxuw52DWqGUZ/24XEJlOax4F/54B3K2HXlQAcd0arwDoN0w6HQltL4IvCyfdyCE3cgYfA0VRXfDDy+q0pcyedgQbj2/ZZ3P1bFJGLPv7Mej36zjpZ83k1tUxhMjEiXJ24vWsOVn+HEyHNwC8f2g540QnwyNO0NlGZQcgpwdsP5rSP0aUr+CZn1gxMsQnWj1KxDCLiTB10BhaQXXTdvAk1UJXBW1i6hzSO7Vgvy8efmqrjQM8uWDxTspq6ji6VEd8fKSJH9OKivgu3tg9TQIbwXjPjc982PfPG3e4BsIobGQcB4MfR5SZsDPT8Db/aHvJBj4MPgEWPYyhLAHSfBnUVWluWfGGlZn5BLSoT9RO2dARRl4n/sFUqUUT4xIxNfbi7d/205Fpebvl3eSnnxdlRXBlzfC5u+h/30w4OGa/T95+0L3a6HtMPjlCVj8bzNef/VnEOyiazOEqAEZcDyLf/y4kZ837OPxEYk0734hVJRA1lq7nV8pxcMXt2PSBa2YsTKTV+dttdu565XiXJh2GWz+AYa9ABc+Ufs34aBwGPUGXDkN9qbCu4PMRVkh3JQk+DP4YmUmb/+2nfG945iQnGDGaOG46ZL2oJTivovaMLp7E17+ZTPfpeyx6/k9XkUZfDYW9vwJV3wISbec2/kSR8INc6CiFN67SPblFW5LEvxppOw6xORv1tGvVQRTRnYwwyYh0dCohd0TPJgk//zoTvSMb8h9n69lTeYhu7fhsX58xPyfXPof6HCpfc7ZpAfcPA9CGsN/r4Rdq+xzXiGcSBL8KRSXVXLPjDVEBPvxxrju+NiO+THFJUPGMqiyfxVkP28bb1/bg8gQP275eCX7D5favQ2P8+c0MxUy+S7oNMa+527QDK6bCYHh8MloM2wjhBuRBH8Kz81NY/v+Ql68osvJ89zj+kBxjpl+5wDhwX68e31P8ovLuf+LtVRV2X8hmsfYtQrm/B+0uAAGT3FMG6GxcP0s8AmEaZfCAblGItyHJPgTLNiUzbRl6dzUrznJrSJOPiCur7lPX+KwGNo1DuWxEYn8tnk/7y/e4bB23Fp5MXx9sxlCGfM+eNkc11bDBNOT11VmuKY413FtCWFHkuCPkVtYxoNfptA2OoQHhrY99UHhLc1H9sw/HBrLNb3juCgxmn/8sJF1u/Ic2pZbWvgvyNkOI1+HwEaOby+yDVz1KRzKgC9uMPPthXBxkuCP8Y8fNpJTWMZLV3XB3+c0PUKloFlvyFzu0FiUUvxzTGcigv24a/pqisokoRy1bz0sfgW6jIMWA5zXbnxfGPESbF8APz/uvHaFqCNJ8Ef8mZHL9BWZ3HhewtnL+TZLMrVNCg86NKYGgb68fFVXdhwo5IUfNzu0LbdRVQWz7wb/MLjoGee33/066H0HLHsTVn/q/PaFqAVJ8EBllebxb1OJDvXj7sFtzv6EZr3N/S7HDtMA9GkRzrV94vlgyQ5WpcvYLyvfg10rYOhzZmGSFS56BpoPgDn3wb4N1sQgRA1Iggc+WZbO+j35PD4ikWC/GlRviO0GXt4OH6ap9uDFbYkJ9eehr1Iorah0SpsuqSQfFjxrSvx2vsq6OGzecPm7pszwFxOgrNC6WIQ4g3qf4PcfLuWFnzbRr1UEwzvF1OxJPgFmlycHX2itFuLvw7OjO7E1u4DX59fjaXpLXzczWIY8dXzxMCsER8Hl78CBzTD3AWtjEeI06n2C//cvmykuq+RvozrUrshX0yTYvcrUHHeCC9pGMbpbE/7z6zY27s13SpsupfAALH0DEkeZT1CuoMVAOP8BWPMprPnM6miEOEm9TvDb9hcwfUUm43vH0TIyuHZPbpZkCo/tTXFMcKfw+IhEgv29eeLb9ThiJy6XtuhlKC+CCx61OpLjDXwY4s+DufdD7k6roxHiOPU6wf/zh434e3tx54Wta//k6gutThqmAWgY5MtDF7fjj505fLtmt9PatVzeLrMrU5dxEHma9QlW8bLBZW+B8oKvb4OqenyNRLicepvgV6Xn8OP6fdw2oCURwX61P0FYEwht6tQED3BVz2Z0adaAZ+dsJL/EOcNDlvvtn4CGgQ9ZHcmpNYiDYf+CzGWmlrwQLqJeJnitNc/P3UhkiB83929e9xM1S3J6gvfyUjw9qgMHC0v598+OqYfjUvKzYM1/zfzzBnFWR3N6na+CxEthwXN23S9AiHNRLxP8L2nZrEzP5Z7BrQn0PYdNrZr1hvxdZgjBiTo3bcC4pDg+WrrT8y+4/vE26EqzjZ4rU8rs5xoUCV/fCuUlVkckRP1L8FprXv55MwnhgVzVs9m5naxZkrl3ci8e4IGhbQn28+bp7zZ47gXX0sOw8n1oPxIancMnLWcJbAQjX4P9G+G3v1sdjRD1L8H/vGEfG7LyuXNQa7xt5/jyG3cC7wBLEnyDQF/uHdyaxVsPMi8t2+ntO8Wf06AkD5LvtDqSmms9BLpda2rl7FppdTSinqtXCV5rzb9/2UJCeCCjusae+wltPmZOthNKFpzK+D7xtIwM4tm5aZRV2H8DEktVVph6L3HJ0LSn1dHUztDnICQWvr3DlDUWwiL1KsH/ZM/ee7VmvSArxZIxVx+bF48NT2THgUI+XrrT6e071IZvIS/TvXrv1fxDYdTrZpXrfAsKoglxRL1J8FprXrFn771a0ySoKoesNfY7Zy0MbBvJ+W0ieXXeFnIKyyyJwe60NmUJwltBm4utjqZuWl4APW80q28znFOzSIgT1ZsE75DeO/zvQuuuFfY7Zy0opXhseHsKSit4dZ6HTJvc8yfsWQ29bwcvN/4VHfIUhDWFmX+VoRphCTf+66k5rTVvLNhKXCM7997BFJ1qEG/JhdZqbaJDuDopjk+WpbPjgAdUNlz5gdkDtfOVVkdybvxCzKyag1vM/HghnKxeJPgl2w6SsiuP2wa0sG/vvVqzJNODt3C64j2DW+Pr7cU/vt9oWQx2UZIHqV9BpzFmUw931/IC6DHBDDllWvMpT9Rf9SLBv/nrViJD/Li8e1PHNNC0FxzOcvqCp2NFhfhz+4CW/LB+Lyt35lgWxzlL+dwUFetxg9WR2M+Qp82smpkTZQGUcCqPT/BrMw+xeOtBbu7X/PT7rJ6rpr3MvUXTJavd3L85USF+PDs3zT0XP2ltFjbFdIEm3a2Oxn78Q2Hkq2ZWzYJnrY5G1CMen+D/8+s2Qv29Gd8n3nGNHF3wZO1H8EBfb+6/qC2rMw4xd91eS2Opk8w/IHuDmX3iaVpdeMxQjbUdAVF/eHSC35pdwI8b9nJ9ckLNtuKrq6MLnqwfY728R1PaRofwrx83ut/ip5Xvg28IdBxjdSSOMeRpCG0iC6CE03h0gp+6cBt+3l5MSE5wfGPNepkqghaPsdq8FA8Pa8fOg0V89keGpbHUSkmeWdzU+Qrwq+XmK+7CP/TIrJqtZ1wApbWmpLySvKJyissqqapyw+E24RIc2K21VnZ+Cd+u3sPVSc0Ir0u999pqmgRVr5gkH9fb8e2dwcA2kfRtEc4r87YwunsTQvx9LI2nRtZ/Y3bI6nqN1ZE4VssLoOdN6KVvsCvqApZUtGHDnnwycorIzC1mX14JhWUVnJjT/X28iAj2o3GoP43D/GkRGUyb6GDaRIfQIiLIMbPDhNvz2AT/4ZKdVFRVcVM/J1UhPPZCq8UJXinFI8PaMfL1xUxduJ37LnKxXZBOZe10iGjjWRdXT5BbWMb8jdn8dvAy7tffwTe381TZ3/HyC6FZo0BaRgbRr1UEwX7eBPrZ8PO2UVZRRUl5JcXllWTnl7A3v4R1u/OYuy7r6JtAoK+NLk0b0C2uAX1ahNMroREBvg6aUCDcikcm+MLSCj5Zls7FHRsTHx7knEZDos2GFC5yAa1z0waM7BLLO79v55o+8USH+lsd0unlbIeMpXDhk6auugeprNIs2JjNp8vTWbjlAJVVmsah/sxt9SS3bZvE0u7zCLnizdpt+A6UlFeybX8Bm/cdZk3GIVZnHmLqwu28+es2fG1e9IhvyIC2kQxuH0XLyOBan194Bo9M8DNWZJJfUsEt/Vs4t+GmSZC+2Ez3c4E/qAeGtuX71Cxe+mkz/xjT2epwTm/tDECZXZE8RFFZBdOWpvPRkp3sySshKsSPW89vwSUdG9OpSZhJuL9sIXTRy7B5JLS9pFbn9/ex0SE2jA6xYVzWrenRNlfszGXx1gP8vuUAf/9+I3//fiNxjQK5KDGaSzo1pluzhnh5Wf+7KZzD4xJ8RWUV7y3aQVJCI7rFNXRu482SIPVLs+CpwTluJmKPcBoFcl3fBD5YvIOb+jenTXSI1SGdrKoK1n4GLQaYfW7dXHFZJdOW7eTt37ZzsLCM5JbhPD4ikcGJ0ficOE4+cDJs+QVm3Ql3LIXgyHNqO9DXmwFtIhnQxpxnz6Fi5m/M5pe0fXy0dCfvLtpBdKgfl3SM4S9dYiTZ1wMel+C/T93L7kPFTBnZwfmNHzsO7wIJHmDSBa34fGUmz89N44MbkqwO52SZy+BQOlww2epIzonWmtkpWTw3J429+SX0bx3BPYPb0CP+DJ0Mb18YPRWmDjRTJ8d9btfiarENArimTzzX9Iknv6Sc+WnZzF2XxX//yODDJTtp0iCAEZ1jGNk1lsSYUBnG8UAeleC11rzz+3ZaRARxYbso5wdw7IKnjpc7v/1TaBjky6QLWvH89xtZsvUAya0irA7peGv+C77B0P4vVkdSZ5v3HeaJmaks255DxyahvDauG70SGtXsydGJMPRZmHs/LH8L+k50SIyh/j5c2q0Jl3ZrwuGScn5J28esNXt4b9EO3l64nVZRwYzqEsvIrrHOu24lHM6jEvyKnbmk7Mrj2cs6WvPR0+Idnk7n+uQEPl6azrNz05g9qZ/rfCwvL4YNMyFxFPi6X1KprNJMXbidl37eRKCvN89c2pGxSXHYavvz7XUzbFsAvzwJCeeZUg0OFOLvw2XdmnJZt6bkFpYxNzWLmav38OLPm3nx5810bWYu0I/oHEOUK1+cF2flUQn+3d+30zDQh9HdHFRUrCaa9YKlb5oFTz6u8cfh72Pj/qFtuHfGWmau3X30opzltvwEpfnQ6QqrI6m1nQcKue+LtaxKz+WSjo155tKOdV9voZTZAeo/58GXN8GtvzptsVfDIF/G945nfO949hwqZvbaPXy7Zg9PfbeBp+dsoE/zcP7SJZahHaKds57E2cpLoHC/uRUdNAvuygqgrBAqy6Cq0tyUAuUFXt7g7Qc+AaaktW+wqXrqHwYBDSEw3GX+7gGUKxWl6tmzp165sm4bFe88UMgFL/7KpAtaWTvvO+07mDEebvwR4vpYF8cJqqo0I99YRE5BGfPuG+ga86Q/vw7Sl8J9G8HLBeKpoVlr9/DIVynYvBRPjerIqK6x9hm/3vE7fDwSEi+FMe9bOhNra/ZhZq3NYvbaPew4UIjNS9GnRSOGdYphSGI0USGuk8TOqqoKcneYOkfZG2H/RnPdJzcdCh2wYb1PEASFQ1AUBEebPSNCGh+5xZhbaKx5M7DT/7FSapXW+qTNix3Wg1dKTQQeAGKA9cA9WuvfHdXeB4t34OPlxbV9HVhUrCaqd3jK/MOlEryXl+Kx4YlcPXUZ7y3azqRBra0NqPQwbP4Rul/nNsm9tKKS5+ak8dHSdHrEN+T1cd2ICQuwXwPN+8Ogx2He38zvUZ877HfuWmoVFcL/DQnh3sGt2ZCVz9x1Wcxdt5dHv0nlsW9T6R7XkCGJ0VzYLopWUS42z76sEDKWQfoS2L0Sdq+G0rz/fb9BHDRsDm2Gms16QqIhKBICI0w5Cd9g8A0Em5/psVf/flZVmFtFqRleLC8yv8el+VB8CEoOmU8BRTnmE0FBtnkjyVwORQdOjtPmB6ExppR0aCwMftLEZkcOSfBKqauAV4CJwKIj998rpRK11nYvkJJXVM7nK3cxsmus9T2L6h2eXGwcHqBPi3CGdojmzV+3cWXPZtaOr26ca0oTuElhsez8Em6dtoo1mYe4uV9zHrqk3cnTHu2h372wayX89Ji5nmNxJ0EpdXS+/f0XtWXTvsP8tH4fP23Ye3SefdOGAVzQNop+rSPo0yKcsAAnl8bQ2vTKN30PW38xnauqclA2iO4AHUebFdLRHSCibd2Hv7xsgF/drhdVlJlPC/lZcHjP8ff5e8wbkbL/75NDhmiUUsuBFK31Lcc8tgX4Umv9yOmeV9chmv/8uo1//LCR7+/uT/uY0DrFbFdf3Ww+bt+30SUWPB1r54FChrz8G6O7NbV28dOnV5qPzHenuPy+q6m787j5o5Xkl5Tz4hVduKRTjGMbLD5kpk6WF8Ntv5mP9i5oz6Fift20n/kbs1m89QDF5ZV4KejUtAF9mjeiZ0IjeiU0pEGgr/0b1xr2psC6LyFtFuTuNI837gwtBkLzAebN0eLCdeWVVRwuqSC/uJzDJRUcLi2noKSCgtIKCksrKCyrpOjI/cSBLet8ncNpQzRKKV+gB/DCCd/6CUg+xfG3ArcCxMXV7ePJ2sxD9G8d4RrJHcyK1nVfQF6m3T9ynauEiCCu75vAe4t3cH1yAomxFvzMinJg2zzoM9Hlk/sPqXu5d8YaGgb68OXtyc75eQU0gKs+gfcugs+uhglzzZCBi4ltEMC43nGM6x1HWUUVqzPMKtrF2w7y/mIz/RKgeUQQnZuG0blpAxJjQmnXOISGQXVM+nm7zdTalBlmr1svb5PQk++CNhc7bLGc1pri8koOFpSRU1hGTlEZuYXm37lFZeQWlXOoqIxDReUcKionr9h8XVhWedZzKwVBvt6MTYqz+4Vsu/fglVKxwG5ggNZ64TGPPwGM11qf9grouVxkLSitcGzN99rYs9r0wC5/z+wt6mLyissZ+K8FtG0cwme39HH++OmqD2H23XDbQodPCTwXHyzewVPfbaBL0wZMva6H84f/Nn0Pn42F9iPgio9d/s3wWCXllazNPMTK9FzWZh4iZVcee/P/V0o7KsSPlpHBJEQEkhAeRNOGgcQ28Ce2QQARwX7HTzWtrIDNP5jfm23zQFdBfD/zt9V+pLmgWQvVyTqv2CTi3MLyI0m6Ommbrw8WlpFTWEpOgfl36Wn2V7B5KRoG+hAW4EPDQF8aBPoQFuBLWIB5LDTAm1B/H0IDfAj28ybE35tgP2+C/My9v4/XOf8NOv0iK3DiO4c6xWN24zLJHSC6o1nwtGuFSyb4sAAf7ruoLY99m8qcdVmM6Bzr3ABSv4LwVubjtAvSWvOPHzbx1m/buCgxmlfHdnPcdo9n0vYSswjqx8kw/ykYPMX5MdSRv4+N3i3C6d3if8k3O7+EjXsPs2nvYTbuPcyOAwX8uH4fOYVlxz1XKWgQ4EPzwFJGM4+Li+cQUZlNnncka6KuJTVyBAXBcXgdBK/FBwBTxK2ySlNeqSmtqKT0SBXOorJKCksrKCqr5HCJGSbJLymnvPL0qSjI10bDIF/Cg3yJDPajTXQI4UG+NAryO3Lve/T7DQN9CfH3dp21JSdwRFY8AFQCJw4cRgH7HNCe67H5mIs6LlJZ8lTGJsXx2R8ZPDsnjUHtogj0ddIb5OF9sHMRnP+Ay12fADNm+tBXKXz9527G947jqVEda79wyZ76TDQbhCx62ewGlXTL2Z/joqJC/YkK9ef8NsfX3MkrLmd3bjFZecXsySuhPHsLHdKn0S3ne3x1KSk+nXnV50Z+qexO8X5FaVYV5ZXb0RqqtEYD3l4KL6Xw9lL4+9jw8/bCz8dGkJ+NQF9vwoN9SYgIIsTf9KAbHNPDbhhket7VvW9L3swdxO5/1VrrMqXUKmAI8MUx3xoCfGXv9lxWs96w5FUoK3LJ8VObl+JvIzsw5q2lvLlgG/cPddLagbRZ5iN2h9HOaa8WSsormfTf1fySto//G9KGOwe1sn76n1Jwyb/MG+Pc+80Mjq7jrI3JzqoTbaLeCikvmbUkNh/oehX0mUjn6EQ6A09ZHagbclS37SVgmlLqD2AxcDsQC7zloPZcT7PeZs7sntVm+bkL6pnQiEu7xjJ14Xau6NnUOTVINsw0U9Wi2jm+rVooLK3g1mkrWbz1IE+P6sC1fROsDul/bN5m4dNnV8PMv5oVlB0utToq+9m5CH5/EbbNNytC+/8fJN1m5qeLc+KQqzZa6xnAPcBjwBqgHzBMa53uiPZcUnVlycxl1sZxFo8Ma4+PTTFl1nocvqq5YL+pl584yrHt1FJecTnXvrecpdsO8uIVXVwruVfz8YerPzUdh69ugg2zrI7o3GhtEvr7l8CHw2FvKgz+G9yTChc+IcndThx2WV5r/abWOkFr7ae17nHsjJp6ISgcwlu79Dg8QHSoP/cOacOCTfuZu26vYxvbOPvI8Myljm2nFg4VlXHNu8tZtzuPN8d35/IeLlKn51R8g2DcDIjtDl9cb2aVuButYcvP8N4QmHaZmb9+yT/hnhTod49ZSSrsxn3mXbmjuN5mmbIL1fs5lQnJCXRsEsqU2evJKy53XEMbZprZM1GJjmujFnIKyxj7znI27T3M29f24OKODl7AZA/+YXDdTGh5oZlquvAFl//9AkyMm36AdwbBp2Pg8F4Y/hLcvQZ632aKdwm7kwTvSM16Q3EuHNhidSRn5G3z4vnLOnOwoJR//rDRMY0UHjSrexNHucTsmQMFpYyduozt+wt45/qeDGrnRkMCvoEw9jPodCXMfxpmTjKrXl1RVZV5Y3+7P3x2lanJ8pdX4c4/oddNpjKjcBgXmjzugZr1NveZyyGyjbWxnEWnpmFMSG7O+4t3MLp7E3rE13DDipra+B3oSlMp0WL7D5cy7p1lZOYW8f6EXpznapug1ITNBy57GxrGw8J/wd61cOXH0MjJ+xCfTmW5We+w6N+wPw0atYRRb0LnK03swimkB+9I4a1NjejM5VZHUiP3XdSGJg0CeOirdZSUn32Jda1s+NZU8Gvcyb7nraXs/BKunrqUXbnFfHhDknsm92peXjDoMbPV36FMeHsgpHxh7ZBN6WFY9ha82g2+uc08NvpdmLQCuo2X5O5kkuAdycvL1KVx8Qut1YL8vHludCe2Zhfwwo+b7HfiohzYsdDy4Zl9+SVc/c4ysvJK+PCGXvRpUbsl7i6rzVBTlCyiNXx9s7l4mbPduTHk7IAfHoGXEuGHhyCsqXnjmbgUOl/hNiWhPY0M0ThasyTY8qNJcoF2HvZwgAFtIrmmTxzvLd7B4MRo+yTBTd+bNQGJI8/9XHWUlVfMuHeWk51fwkc3JtV8z1R30TABbvoJVr4P856CN/tC79uh719NCWtHqCiFjXNg9TSz5aCXzQzB9bkDmp5UFkVYQHrwjlZdz3vXCmvjqIXJw9oT1yiQ+79YS0FpxbmfMG02hDUz0/sssPtQMVe9vYz9h0v5+CYPTO7VvGymlMFf/zCbmC9+Bf7dCebcb78L/RVlZprjrDvhxbbw5Q3m3AMegnvWwZj3JLm7EOnBO1psd7PxQOZy81HaDQT6evPiFV244u2lPD17w7nVjS89bBa09LrJkuGZzJwixr27jEOF5Uy7KYlucQ2dHoPThcbA5e/CgIdh8b/NfPkV75gieImXQuvBZqpqTWawVJSZuv3pi2HnYtj5u9nByDfE/D53HWfK9coQjEuSBO9ovoEQ0xky3ONCa7WeCY24Y0BL3vx1G31aNqr7Rt1bfobKUtOjdLIdBwoZ/84yCkor+OTm3nRp1sDpMVgqopXZzPuCR81F7vXfwIJnzM3LB6LaQ3hLMxHAv4G5AFpaYDadLsg29dZzdpjZT2CGgRJHQbsRJqm70ObS4tQkwTtDXF8zNlpR6lbzfv9vSBtWpucy+etUOsSG0SY6pPYnSZtt9rusnjLqJJv3HWb8u8uprNJ8dmsfOsSGObV9lxIaY8bF+9xhtofLXA5Za2HPGti77n/7iVZVmA2j/YLNhtBRiabHH9Xe/A47aDMN4TiS4J0hri8se9P8UVVvyu0GvG1evD62G8NeXcTtn6xi1qR+tau7X14CW34yNfGd+BE+dXce1763HB+bF5/f1odWUXV4Y/JUobHQ4TJzO5bWpoyEDLV4FLnI6gxxfc19+hJr46iDqFB/Xh3blZ0HCnn4q5TaFSTbvsB83G/vvNkzS7Ye4Oqpywj09ebz2/pKcq8ppSS5eyBJ8M4QHGlqsGQstTqSOkluGcEDQ9vxXUoWL/28ueZPTJttaqck9HdccMeYuy6LCR+soEmDAL66I5mECCeUPxbChckQjbPE9TUJr6rKrfbWrHb7gBakHyzktflbiQkzmy2fUWU5bJoLbS4B7zpusFwLHy3ZyZTZ6+kR15D3ru9FWKCsmBTC/TKNu4pPNhey9qdZHUmdKKV4+tKODGwbyeMzU1mwMfvMT0hfbAqtOXj2TEVlFU/OTOXJWeu5sF00027qLcldiCMkwTtL9Ti8mw7TAPjYvHhjXHfax4Qw8dM/WbTlwOkPTpttdh5qOchh8eSXlHPjRyv5aGk6t/RvztvX9iDAV8aRhagmCd5ZGiZASAyku2+CB1Ov5sMbkogPD+TGD1fwy4ZT7KNeVWWWsLca7LD9aNOy8rn09cUs2XqA50d34tHhidZuji2EC5IE7yxKmbIFGUvdY4OGM4gI9mP6rX1oHxPC7Z+sYvbaPccfsHsVHM5yyPCM1prpf2Rw6RuLKSit4NObezM26SzXA4SopyTBO1NcMuTvhkMZVkdyzhoE+vLJzb3pHteQu6av5qWfN1NZdeSNK22WWSnZ+iK7tnmwoJS7pq/h4a/X0SuhEXPu6k9vT6kIKYQDSIJ3pvjqcXjX3oi7pkL8ffjoxiRGd2vKq/O2cN37yzlwuMRs7tFiAAQ0sEs7Wmu+Wb2LwS/9xg+pWdx/URs+ujGJyBD3WRUshBUkwTtTVCL4hUGG+y14Op0AXxsvXtmFf17emZU7c7nzlU8hZztVbUfY5fyrM3K55r3l3DtjLQkRQcy5qz+TBrWW8XYhakDmwTuTl81sxO2GK1rP5spezejUNIw/P36IqnLFhKUR3BmZU6fSvFpr1mQe4rX5W5m/MZtGQb48NaoD43vHS2IXohYkwTtb/HmmPsvhfRDiRhs910D7mFDahaZwIKA7mwsCueKtpSTGhDK6exNGdo0lKuT01Qe11mTkFPFdShbfrN7N1uwCwgJ8eGBoWyYkJxBUmxo4QghAErzzNT+ybH/n76YIlyfJ2YHal0rk0OdY0GMgn6/M5Os/d/HMnDSenZtGs4aBtIkOpmVUML42L8oqqygtr2JrdgGpe/I4VFQOQK+Ehjx3WSf+0iWGEH9ZtCREXUmCd7bGXcAv1DMTfNosc99uBAG+Nq5PTuD65AS2ZhfwQ2oWaXsPs2XfYX7bvJ+KKo2PzQtfmxfx4YFc3KExHZuEcX7rSOLCHTN3Xoj6RhK8s9m8TdmCHb9bHYn9bZgFMV2gYfxxD7eKCmbSoNZHv9ZaoyzcfFuI+kJm0VghoT/kbDObL3iKvN2we2WNSgNLchfCOSTBW6F6HN6TevEbvzP3iaOsjUMIcZQkeCtEdzJ7YO5caHUk9rNhFkS2h4jWZz9WCOEUkuCt4OUFCf08pwdfsN8s3kp03s5NQoizkwRvlYT+cCjdI+rSsPE7s5+nE7fmE0KcnSR4qxydD7/I2jjsIW0WNGwO0R2sjkQIcQxJ8FaJbA8Bjdx/mKY4F3YsNMMzMjtGCJciCd4qR8fhF7p3ffhN30NVBbSX2TNCuBpJ8FZqeQHk74IDW6yOpO7WfwthzaBJd6sjEUKcQBK8lVpeaO63zbM2jroqzoVt86HDpTI8I4QLkgRvpYbxEN4Ktrppgt84F6rKocNlVkcihDgFSfBWa3mhmUlTXmJ1JLW3/htoEAexMjwjhCuSBG+1VhdCRbHZjNudFOXA9gWm9y7DM0K4JEnwVkvoBzZf9xuH3zjHzJ6R4RkhXJYkeKv5BkFcH9g63+pIamf9N9AwAWK6Wh2JEOI0JMG7gpYXQvZ69ykfXJQD23+V4RkhXJwkeFfQqnq6pJv04tNmga6U4RkhXJzdE7xS6lal1AKl1CGllFZKJdi7DY8T3RGCo91nuuS6L830zsadrY5ECHEGjujBBwI/AVMccG7PpJQZptm+ACorrI7mzPJ2HdlP9koZnhHCxdk9wWut/621fh7wgDKJTtTmIrMyNHO51ZGc2bovzX3nK6yNQwhxVjIG7ypaDTbTJTfNtTqSM0v5HJr2gkYtrI5ECHEWlif4I2P2K5VSK/fv3291ONbxC4EWA49snuGi1SX3pprZPp2vsjoSIUQN1CjBK6WeOXLB9Ey3gXUJQGs9VWvdU2vdMzIysi6n8BzthkPuTsjeYHUkp7buc1A2mT0jhJvwruFx/wY+OcsxHrD3nMXaXALcY1aJutruSFVVZvy91WAIirA6GiFEDdQowWutDwAHHByLCImGZklmmGbAg1ZHc7z0xZC/G4Y8ZXUkQogacsQ8+MZKqa5AmyMPJSqluiqlGtm7LY/UbjhkrYVDmVZHcry108E3GNpeYnUkQogacsRF1tuB1cCnR76ec+TrkQ5oy/O0G2HuXWk2TUk+rP8aOo42tXOEEG7BEfPgp2it1SluH9q7LY8U3hIi25lhGleR+hWUF0H3662ORAhRC5ZPkxSn0G447FwMhS5y2ePPjyCqAzTpYXUkQohakATvijqMNsW8Ur+2OhLISoE9q6H7dVKaQAg3IwneFTXuCNGdIGW61ZHA6mlg84POV1odiRCiliTBu6ouV8HuVXBgi3UxlBdDygxIHAmBMglKCHcjCd5VdboClJdJsFbZMAtK8szwjBDC7UiCd1UhjU1tmrUzzCpSK6x41xQVS+hvTftCiHMiCd6VdRkLeRmQsdT5bWf+Abv+gN53yMVVIdyUJHhX1m44+ARZc7F16RvgHwZdxzm/bSGEXUiCd2W+QeYC5/qZ5oKns+TuNPuu9rgB/IKd164Qwq4kwbu6ruOhNM9stOEsy982F3iTbnVem0IIu5ME7+oS+kHjTrDsTedsBFKSB39+DB0vh7Amjm9PCOEwkuBdnVLQdxLs3wjb5jm+vT8/hrIC6DPR8W0JIRxKErw76DAaghubC5+OVFoAi1810yJjuzq2LSGEw0mCdwfevtD7Vtg2H/Y5cDu/Zf+Bwmy48EnHtSGEcBpJ8O6ixw3gHQDLHNSLLzwIi18x9eib9XJMG0IIp5IE7y4CG5k56SmfQ36W/c+/6CUoL4RBj9v/3EIIS0iCdyfJk8xMmvlP2/e8hzLhj6nQZRxEtbPvuYUQlpEE704atYC+E2HNp7Brlf3OO/8ZQMHAh+13TiGE5STBu5vzH4DgaPj+QfsUIdv0vSmFkDwJGjQ79/MJIVyGJHh34xcCg6fA7pXnXqOm8ADMuhOiO8KAh+wSnhDCdUiCd0edrzb7o/4yBUry63YOreG7e8zK1cveBm8/e0YohHABkuDdkZcXDPuX6YF/dRNUVtT+HGunQ9psuOBRs0WgEMLjSIJ3V016wPAXYMtP8OMjtXvu1nmm9x7XF5LvdEh4QgjreVsdgDgHPW+Eg9tg6evQqCX0uf3sz9n8E8y4BiJaw1WfgJfN8XEKISwhCd7dDXnK1G//4WEoOwzJd51+PD3tO/jyBohsB9fNlI20hfBwMkTj7rxsMHoqJI4y89n/cx5s//V/39faDMl8OAJmjIfoDnD9LEnuQtQD0oP3BL5BcOVHsOVnmHs/fDwKbL7gFwo2HzicBSGxMPQ5U9PGN9DqiIUQTiAJ3pO0HgITl5mVrocyoTTflABO6AddrpapkELUM5LgPY1PAPS62eoohBAuQMbghRDCQ0mCF0IIDyUJXgghPJQkeCGE8FCS4IUQwkNJghdCCA8lCV4IITyUJHghhPBQSmttdQxHKaX2A+l1fHoEcMCO4Tibu8cP7v8a3D1+cP/X4O7xgzWvIV5rHXnigy6V4M+FUmql1rqn1XHUlbvHD+7/Gtw9fnD/1+Du8YNrvQYZohFCCA8lCV4IITyUJyX4qVYHcI7cPX5w/9fg7vGD+78Gd48fXOg1eMwYvBBCiON5Ug9eCCHEMSTBCyGEh5IEL4QQHsrtE7xSaqJSaodSqkQptUop1d/qmGpDKXW+UmqWUmq3UkorpSZYHVNNKaUeUUqtUErlK6X2K6VmK6U6Wh1XbSil/qqUSjnyGvKVUkuVUsOtjquulFKTj/wevW51LDWllJpyJOZjb3utjqs2lFIxSqmPjvwdlCilNiilBlgdl1sneKXUVcArwHNAN2AJ8L1SKs7SwGonGEgF7gaKLY6ltgYCbwLJwCCgAvhFKdXIyqBqaRfwENAd6AnMB75VSnW2NKo6UEr1AW4BUqyOpQ42ATHH3DpZG07NKaUaAIsBBQwH2gN3AtkWhgW4+SwapdRyIEVrfcsxj20BvtRaP2JdZHWjlCoAJmmtP7Q6lrpQSgUDecClWuvZVsdTV0qpHOARrfXbVsdSU0qpMOBPTIJ/AkjVWk+yNqqaUUpNAcZord3q0181pdRzwACt9XlWx3Iit+3BK6V8gR7ATyd86ydMj1I4XwjmdyrX6kDqQillU0pdjflUtcTqeGppKqZjM9/qQOqoxZFhyh1KqelKqRZWB1QLlwLLlVIzlFLZSqk1SqlJSilldWBum+AxBX1swL4THt8HNHZ+OAIzXLYGWGpxHLWilOp05NNTKfAWcJnWep3FYdWYUuoWoBXwuNWx1NFyYAJwCeYTSGNgiVIq3MqgaqEFMBHYDgzF/B38HfirlUEBeFsdgB2cOMakTvGYcDCl1EtAP6Cf1rrS6nhqaRPQFWgAXA58pJQaqLVOtTKomlBKtcVcg+qvtS6zOp660Fp/f+zXSqllmGR5PfCSJUHVjhew8phh4dVKqdaYBG/pxW537sEfACo5ubcexcm9euFASqmXgbHAIK31dqvjqS2tdZnWeqvWuvqPdA1wr8Vh1VRfzKfZVKVUhVKqAhgATDzytZ+14dWe1roAWA+0tjqWGsoCNpzwWBpg+WQPt03wR3orq4AhJ3xrCO43fuq2lFKvAOMwyX2j1fHYiRfgLonxW8yMk67H3FYC04/82+169Uopf6AdJnG6g8VA2xMea0Pd97awG3cfonkJmKaU+gPzQ74diMWMo7qFIzNPWh350guIU0p1BXK01hmWBVYDSqk3gGsxF5lylVLVn6YKjvTCXJ5S6u/AHCATc5F4HGb6p1vMhddaHwIOHfuYUqoQ8/vj8kNMAEqpF4DZQAbmE/jjQBDwkZVx1cLLmGsGjwIzMFO27wImWxoVgNbarW+Yixs7MRfIVgHnWx1TLeMfiLlmcOLtQ6tjq0Hsp4pbA1Osjq0Wr+FDTE+rFDNv+RdgqNVxneNr+hV43eo4ahHvdGAP5tPGbuArINHquGr5GoYDa4ESYDMmwSur43LrefBCCCFOz23H4IUQQpyZJHghhPBQkuCFEMJDSYIXQggPJQleCCE8lCR4IYTwUJLghRDCQ0mCF0IID/X/omdqey8oOQIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "\n",
    "x = np.linspace(0, 2*np.pi, 100)\n",
    "ax.plot(x, f(x), label=\"f(x)\")\n",
    "ax.plot(x, f_deriv(x), label=\"f'(x)\")\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks okay (?), how do we test it? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.00e-01  7.10e-01  nan  2.51e-01  nan\n",
      "2.50e-01  3.11e-01  1.14  8.15e-02  1.54\n",
      "1.25e-01  1.37e-01  1.13  2.25e-02  1.81\n",
      "6.25e-02  6.32e-02  1.09  5.87e-03  1.92\n",
      "3.12e-02  3.02e-02  1.05  1.50e-03  1.96\n",
      "1.56e-02  1.47e-02  1.03  3.78e-04  1.98\n",
      "7.81e-03  7.26e-03  1.01  9.48e-05  1.99\n",
      "3.91e-03  3.61e-03  1.01  2.38e-05  2.00\n"
     ]
    }
   ],
   "source": [
    "x = 1\n",
    "h0 = 0.5\n",
    "reduce_by = 2\n",
    "n_iter = 8\n",
    "\n",
    "# initiate with nans \n",
    "err1 = np.nan*np.zeros(n_iter)\n",
    "err2 = np.nan*np.zeros(n_iter)\n",
    "order1 = np.nan*np.zeros(n_iter)\n",
    "order2 = np.nan*np.zeros(n_iter)\n",
    "\n",
    "# compute function and derivative at x\n",
    "f_x = f(x)\n",
    "f_deriv_x = f_deriv(x)\n",
    "\n",
    "# create our vector of h-values \n",
    "h = h0 * (1/reduce_by)**(np.arange(n_iter))\n",
    "\n",
    "for i, hi in enumerate(h):\n",
    "    f_true = f(x+hi)\n",
    "    err1[i] = np.abs(f_true - f_x)\n",
    "    err2[i] = np.abs(f_true - (f_x + hi*f_deriv_x))\n",
    "    \n",
    "    if i > 0: \n",
    "        order1[i] = (err1[i-1] / err1[i])/reduce_by\n",
    "        order2[i] = (err2[i-1] / err2[i])/reduce_by\n",
    "    \n",
    "    print(f\"{hi:1.2e}  {err1[i]:1.2e}  {order1[i]:1.2f}  {err2[i]:1.2e}  {order2[i]:1.2f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'error')"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Substituting symbol O from STIXNonUnicode\n",
      "Substituting symbol O from STIXNonUnicode\n",
      "Substituting symbol O from STIXNonUnicode\n",
      "Substituting symbol O from STIXNonUnicode\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAGCCAYAAADddtWoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABGLUlEQVR4nO3dd3yNd//H8dc3QyRCrCTErMTexBYSW9Ss0VJ7tdVaraq7vW93x90WLdWiSqnWLLXVJonYNUKJ2js2Rezx/f1xhR8aJHLOuc74PB+P82hynetc1/vEaT65ru9SWmuEEEIIS3MzO4AQQgjnJAVGCCGEVUiBEUIIYRVSYIQQQliFFBghhBBWIQVGCCGEVUiBEUIIYRVSYIQQQliFh9kBrEEpFQmMwCigI7TWY573muzZs+v8+fNbO5qwkWvXrpEhQwazYwgHJp+hlNm6det5rbV/cs85XYFRSnkAI4FawAVgi1Jqrtb61LNelz9/frZs2WKLiMIGoqOjCQ8PNzuGcGDyGUoZpdTRpz3njLfIKgLxWuvjWuvrwFzgZZMzCSGEy7G7AqOUqqGUWqCUOqmU0kqpTsns85ZS6rBS6qZSaqtSKuyRp4OA4498fwLIZeXYQgghnmB3BQbwBXYBfYAbTz6plGqDcQvsc6AssB5YopTK+2CXZI4pM3oKIYSN2V0bjNZ6MbAYQCk1KZld+gOTtNbjk75/RynVAHgTGAScBPI8sn9u4Ehy51JK9QB6AAQGBhIdHZ32NyDsQmJiovx7ijSRz1Da2V2BeRalVDqgPPDVE08tB6omfb0ZKK6UygOcB5oDdZI7ntZ6HDAOIDQ0VEuDnvOQBlqRVvIZSjuHKjBAdsAdOPPE9jMkFRGt9V2lVD9gFcYtwJFa6wSbphRC2NSVK1c4e/Ysd+7csdgx/fz82LNnj8WO54g8PT0JCAggU6ZML/R6RyswDzzZpqIe3aa1XggsTMmBlFKNgcYhISGWSyeEsJkrV65w5swZcuXKhbe3N0ol1wybelevXiVjxowWOZYj0lpz48YNTp48CfBCRcYeG/mf5TxwD8jxxPYA/nlVkyJa64Va6x5+fn5pzSaEMMHZs2fJlSsXPj4+FisuApRS+Pj4kCtXLs6ePftCx3CoAqO1vg1sBeo+8VRdjN5kQggXc+fOHby9vc2O4bS8vb1f+Naj3d0iU0r5Ag/uV7kBeZVSZYCLWutjwHBgslJqM7AOeANj7MvYFzyf3CITwsHJlYv1pOVna49XMKHA9qSHN/Bx0tefAGitfwX6Ah8BcUB1IFJr/dTpCp5FbpEJIYR12N0VjNY6muQHSz66zxjguRNY2sK9+5p79zXpPOyxVgshhHnkt2Iazd1+knojYlj85ym0lgkDhBDiAZcvMEqpxkqpcZcvX36h1+f0S086DzfemrqNFt+vZ/PhixZOKIQQjsnlC0xa22CqhWRnSZ8aDH2lFAl/36D1Dxvo/ssWDpxNtHBSIYQru3TpEoGBgRw8ePDhtoEDB1K37pOdav9fy5YtGT58uC3iJcvlC4wluLspWlfIQ/R7EQyoX5gNBy9Q/5s1/Gvun5y9etPseEIIO3f9+nWmTp3K3bt3n7rP559/TmRkJMHBwQ+3xcXFUbp06ae+ZvDgwXz22We86B2atJICY0He6dzpFRFCzIBw2lfOx8w/jhM+LJoRK/aReOvpHxwhhOuKi4ujVatW1KxZk927d9OpUyeCg4MpVqwYmzZtAowC9OOPP9K1a9fHXrtjx45nFpiSJUtSoEABpkyZYtX38DR214vMGWTz9eK/TYrTqWp+hi3by8hV+5m66Rh96xSkTYU8eLpLXRfCmj5euJv4hCtpOsa9e/dwd3dP8f7FgjIxuHHxVJ3jyJEjREZGEhUVxeXLlxk1ahQTJkzg9u3b1K5dm3bt2rF//34WL16Mm5sb1apVe/ja06dPc+bMGdKlS0dkZCQxMTHkyJGDH3/8kYiIiIf7NWnShOnTp9OrV69UZbMEl/9Nl9ZG/mfJnz0Do9uVY+5bVSmQPQMfzdtF/W/WsGz3aelxJoSgb9++hIeHU7hwYaZNm8a5c+fo3bs33t7etG3bloMHD5KQkEBsbCzly5d/bNDj9u3bARg9ejT9+vVjx44dlChRgv79+z92jooVK7J582Zu3PjH8lpW5/JXMA8mxgwNDe1urXOUzZuFX3tWZtWes3y59C96Tt5K+XxZ+FdkEcrny2qt0wrhslJ7JZEca092efz4cRYsWMDixYsBqFGjBqNGjeKTTz4B4P79+wD4+Phw9OhRcubM+djr4+Li8PPzY+bMmeTIYUzP2LJlSwYNGvTYfkFBQdy5c4eEhITH2m9sweULjK0opahTLJDwwv7M2nqC4Sv28cr3G2hQPAfvNyhMAX9fsyMKIWwoKioKrTXFihUDoH79+o81xu/cuZNSpUqRJUsWbty4QWBg4GOvj4uLo3Hjxg+LC8CBAwd4ctqrB/O0mXEF4/K3yGzNw92N1yrmJWZAOP3rFiJ2/znqjljDv+ft4tzVW2bHE0LYSEKCsUzVk4UDjNU058yZQ9++fQHInj07ly5demyfuLg4qlSp8ti27du3U6ZMmce2XbxojM3z9/e3UPKUkwJjEp90HvSuXZDoARG0rZiXaZuPET4sim9X7ef6belxJoSzCwoKAoxbZU/q168flSpVomPHjgCULVuW+Pj4h89fv36dAwcOULZs2cdel1yB2bVrF0FBQckWMmtz+QJjzUb+lPDP6MWnzUqwvF8Nwgr6M3zFPmoOi2b65mPcvXfflExCCOtr1KgRgYGBfPrpp9y7dw8w2n169uzJ6dOnmTVrFm5uxq/o+vXrs2fPHi5cuAAY3ZMBSpUq9fB4Fy5c4MSJE/8oMLGxsTRo0MAG7+ifXL7A2MtsysH+voxtX57Zb1Yhb1YfBs35kwYjY1kZf0Z6nAnhhLJly8bKlSs5deoU5cqVo0WLFnTs2JFatWqxcOFCfH3/v122ZMmSVKxYkRkzZgBGgSlYsCAZMmR4uM/27dvx9PR82KYDcPPmTebOnUv37lbrw/RMSn55GUJDQ/WWLVvMjgEYS5Uu232GoUv/4tD5a1R8KSuDGhahbN4sZkdzGNHR0YSHh5sdQ9jAnj17KFq0qMWPa29LJi9dupQ+ffoQHx+f4vE5o0ePZv78+SxfvjxN537Wz1gptVVrHZrccy5/BWOPlFI0KJGDZf1q8GmzEhw6l0jzMevpNXUbR85fMzueEMIEDRo0oFevXpw4cSLFr/H09OS7776zYqpnk27KdszT3Y32lfPRvGwuxq05xPg1h1gef5p2lfLxTq0Qsvl6mR1RCGFDvXv3TtX+PXr0sFKSlJErGAfg6+VB/7qFiBkQTqvQPEzeeJSaw6IZHXWAG7fvmR1PCCGS5fIFxuxeZKkRkCk9nzcvybK+YVQJzsawZXsJ/yqKmX8c5959aUsTQtgXly8w9tKLLDVCAjIyvkMoM3tWIaefN+/P3knkyFii/jorPc6EEHbD5QuMI6v4UlbmvlWVMe3KcevuPTpP+oO24zex88TfZkcTQggpMI5OKUVkyZws71eTj5sUZ++ZqzQZtY7e07dz/OJ1s+MJIVyYFBgnkc7DjY5V8xMzIJy3I0JYHn+aWl9H88nCeC5du212PCGEC5IC42QypvfkvfqFiX4vghZlczNp/WFqDIvi++iD3LwjPc6EELYjBcZJ5fBLz5CWpVjatwYV82dlyNK/iPgqmt+2npAeZ0IIm5AC4+QKBWZkQqcKTOteCf+MXrw3aweNvo0laq/0OBPC2R0/fpzw8HCKFStG6dKlmTNnjk3PLwXGRVQNzs68t6rx3WtluX77Hp1/Mnqc7Tj+t9nRhBBW4uHhwTfffEN8fDwrVqygT58+XL9uu84/Ll9gHGmgZVq5uSkalw5iZf+a/LdxMfaeuUrT0evoNU3mOBPC3l26dInAwEAOHjz4cNvAgQOpW7fuU1+TM2fOh9P3BwQEkCVLFs6fP//w+ZYtWzJ8+HCrZXb5AuOIAy3TKp2HG52qvUTMgHB61wph9Z6z1Bkew3/m7+J8oqyqKYStXb9+nalTp3L37tMXG/z888+JjIwkODj44ba4uDhKly6donNs2bKFO3fukCdPnofbBg8ezGeffYa1/sB2+QLjyjKm96R/vcLEvB9Omwp5mLrpGDWHRvHNyn1cuyWragphC3FxcbRq1YqaNWuye/duOnXqRHBwMMWKFWPTpk2AUYB+/PFHunbt+thrd+zYkaICc+HCBTp06MCECRNQSj3cXrJkSQoUKMCUKVMs+6aSSIFJq/v3IfGc2SnSJCBjev7XvCQr+tWgRiF/vlm5n5rDopi84Qh3ZFVNIazmyJEjREZGMnz4cC5fvsyoUaOYMGECu3btInPmzLRr1w6tNYsXL8bNzY1q1ao9fO3p06c5c+YM6dKlIzIykgwZMhAcHExUVNRj57h16xbNmzdn0KBBVK1a9R8ZmjRpwvTp063y/qTApNXG0TCmEuxfaXaSNCvg78v3r5dnzltVKeDvy7/n76bu8Bh+33lKepwJYQV9+/YlPDycwoULM23aNM6dO0fv3r3x9vambdu2HDx4kISEBGJjYylfvvxjVx/bt28HjEXF+vXrx44dOyhRogT9+/d/uI/Wmk6dOlGrVi3at2+fbIaKFSuyefNmbty4YfH3J+vBpFXB+hA3Haa+AtX6Qq2PwN3T7FRpUi5vFn7tUZnVf51lyNK/6DVtG6XzZOaDBkWoEpzN7HhCPN+SD+D0n2k6hPe9u+Ceil+ROUpCwy9TvPvx48dZsGABixcvBqBGjRqMGjWKTz75BID79427Bz4+Phw9epScOXM+9vq4uDj8/PyYOXMmOXLkAIxG+0GDBj3cZ926dfz666+UKlWKefPmATB58mRKliz5cJ+goCDu3LlDQkLCY+07liAFJq38C0H3VbD0A1j3DRzbAC0ngl9us5OliVKK2kUDCS8cwOxtJxixYh+vjd9IRGF/BjYsQpEcmcyOKIRDi4qKQmtNsWLFAKhfv/5jje07d+6kVKlSZMmShRs3bhAYGPjY6+Pi4mjcuPHD4gJw4MABQkJCHn5fvXr1h4Xqaby9vQHkCsZueXpD45GQPwwW9oGx1aHZWCjcwOxkaebupmgdmocmpYOYtP4IY6IO0HBkLC3K5qZ/vULkyuxtdkQh/ikVVxJPc+PqVTJmzGiBMMlLSEgA+EfhAEhMTGTOnDl8/fXXAGTPnp1Lly49tk9cXBx9+vR5bNv27dsfdktOqYsXLwLg7++fqtelhLTBWFLJltBzDfjlgeltYNmHcNc5JppM7+nOGzWDWfN+BN3DCrBwZwIRX0Xz+eI9/H3dOd6jELYUFBQEGLfKntSvXz8qVapEx44dAShbtizx8fEPn79+/ToHDhygbNmyj73uRQrMrl27CAoKSrbQpZUUGEvLFgxdV0CF7rBhFPzUEC4dNTuVxWT2Sce/IosS9V44jUsFMT72EDWGRjE2RibTFCI1GjVqRGBgIJ9++in37hn/71y9epWePXty+vRpZs2ahZub8Su6fv367NmzhwsXLgBG92SAUqVKPTzehQsXOHHiRKoLTGxsLA0aWOduixQYa/BMD42+glY/w/l98EMY7FlkdiqLypXZm69bl2Zx7zDK5cvCl0uMyTRnbpHlm4VIiWzZsrFy5UpOnTpFuXLlaNGiBR07dqRWrVosXLgQX1/fh/uWLFmSihUrMmPGDMAoMAULFiRDhgwP99m+fTuenp4P23RS4ubNm8ydO5fu3btb7o09Qrl691OlVGOgcUhISPf9+/db/gQXD8NvnSFhO1R6A+p+Ah5elj+PydYfPM+QJX+x48RlCgdmZGDDwkQUDnisW6UtRUdHEx4ebsq5hW3t2bOHokWLWvy4V63cBpNaS5cupU+fPsTHx+Pu7m6RY44ePZr58+ezfPnyZ+73rJ+xUmqr1jo0uedc/grG6lPFZH0JuiyDym/BprEwoR5cPGSdc5moanB25vWqxui2xvLNXSZtoc24jWw/dun5LxZCPFeDBg3o1asXJ06csNgxPT09+e677yx2vCe5fIGxCQ8vaPAFvDoNLh2GH2rC7rlmp7I4pRSNSuVkRf+afNq0OIfOJdJ8zHrenLKVQ+cSzY4nhMPr3bs3+fLls9jxevToQeHChS12vCdJgbGlIo3gjbXgXxhmdYJF/eHOTbNTWZynuxvtq+QnekAEfWoXJGbfOeqOWMOHc//k7FXne79CiORJgbG1zHmh8xKo+g5smQAT6sCFg89/nQPy9fKgX91CxAyIoG3FvPz6x3HCh0UzfMU+EmUyTSGcnhQYM7h7Qr3PoO1MuHwSfqgBf/5mdiqr8c/oxafNSrCif00iCgfw7ar91Bwaxc/rj3D7rkymKYSzkgJjpkL1jVtmgSVgdldY8A7csfx0DfbipewZGN2uHPN6VaNgoC+DF+ym7ogYFu5I4L50bRbC6UiBMZtfLuj0O1TvD9t+gfG14Nw+s1NZVZk8mZnevTI/da6At6c770zfTrMx61h/4PzzXyxEMlx9uIU1peVnKwXGHrh7QJ3B8PpsSDwL42oaMzQ7MaUUEYUD+L13GF+3Ks35q7do++MmOk7cTHzCFbPjCQfi6elplYkaheHGjRt4er7YDPFSYOxJSB3jlllQOZj3Bsx7C25fMzuVVbm7KV4pn5vV74Xzr8gixB3/m0bfxdL/1zhOXLpudjzhAAICAjh58iTXr1+XKxkL0lpz/fp1Tp48SUBAwAsdQ2ZTtjeZckKH+bBmKMQMhZNbodUkCLD8SGV7kt7TnR41gmkTmpcxMQf4ad0RFu08Rfsq+XgzPJjsvs43+4GwjEyZjKUjEhISuHPnjsWOe/PmTdKnT2+x4zkiT09PAgMDH/6MU8vlp4p5IDQ0VG/ZssXsGI87FA2zu8OtqxA5FMq2B5OmXrG1hL9vMGLFPmZvO0F6T3c6Vc1PjxoFyOyTLkWvl6liRFrJZyhlZKoYR1Ug3Lhllqei0cNsTg+j2LiAoMzeDGtVmuX9alK7aCDfxxwkbEgUI1bs48pNy/2VKoSwHikw9i5jILSfCxEfwq7fYFx4mpeCdSQhAb5891pZlvQJo2pINkau2k/YkChGRx3gmgzWFMKuSYFxBG7uUPN96LjQaPQfXxu2TAQXur1ZJEcmfmgfyqJ3qlM+XxaGLdtLjaFRjF9zSNahEcJOOW2BUUotUEpdUko5zxD5/NWNW2b5q8OifvBbF7jpWl16S+TyY2KnCsx+sypFc2bif4v3UCNpVoBbd6XQCGFPnLbAACOADmaHsLgM2aHdb1B7MMTPN6aZSYgzO5XNlc+XhSndKvFrj8rkz5aBwQt2EzEsmumbj3Hnnkw/I4Q9cNoCo7WOApyzRdzNDcL6Q+fFcO82TKgLm8a51C2zByoVyMavPSszuWtFAjKlZ9CcP6n9dQxrT97hrhQaIUxl8wKjlKqRdPvqpFJKK6U6JbPPW0qpw0qpm0qprUqpMFvndAh5Kxu3zApEwJIBMLM93Pjb7FQ2p5QirKA/c9+qysROoWRM78GPf96m3og1zI87KfOcCWESM65gfIFdQB/gH/M7KKXaACOBz4GywHpgiVIq7yP77HrKI49t3oId8ckKr80wZmfeuwR+CIMTW81OZQqlFLWKBLLoneq8U9YLT3c3+syIo+HIWJbuOiWjvIWwMVMHWiqlEoG3tdaTHtm2Cdipte7+yLb9wG9a60GpPH540vFbPuX5HkAPgMDAwPIzZsxI7VuwK5ku76VY/Feku32RQwU6cCJ3E5cZmPmkxMREfDJkYPPpe8w7cJvT1zT5MrnRPMST0v7uKBf9uYiUS0xMxNfX1+wYdi8iIuKpAy3tqsAopdIB14HXtNazHtlvNFBCa10zlccP5xkF5lF2OZL/Rdy4BPPfhr8WQeFIaDrauMpxMY+Owr577z7z4hIYuWofxy/eoEyezLxbrxDVQ7JLoRFPJSP5U8aRRvJnB9yBM09sPwPkSM2BlFIrgVlApFLqhFKqimUi2jnvLNBmCjQYAvtXGL3Mjm82O5WpPNzdaFk+N6vfDeeLFiU5e+Um7Sdsps0PG9l06ILZ8YRwWvZWYB548rJKJbPt2QfQuo7W2l9r7aO1zq213pDcfkqpxkqpcZcvX37RrPZHKaj8BnRdbgzS/KkhrBsJ9127V5WnuxuvVcxL1IBwPm5SnCMXrtFm3EZe/3ET245dMjueEE7H3grMeeAe/7xaCeCfVzUWobVeqLXu4efnZ43DmytXOei5Boo0ghX/gelt4Jr8xe7l4U7HqvlZ834EHzUqSvypK7QYs57OP23mzxNO9IeGECazqwKjtb4NbAXqPvFUXYzeZCK10vtBq5+h0ddwKAbGVnf5W2YPpPd0p1tYAWLfj2BA/cJsO/Y3jUetpefkLfx12rVmSBDCGswYB+OrlCqjlCqTdP68Sd8/6IY8HOiklOqmlCqqlBoJBAFjbZ3VaSgFFbpBt5Xg4QWTXoY/nWcGnbTK4OVBr4gQYgdG0LdOQdYfuEDDkbG8PW0bB84mmh1PCIdlxhVMKLA96eENfJz09ScAWutfgb7AR0AcUB2I1FoftUYYp2yDeZqcpaD7asgdCrO7GguaydiQhzKl96RvnULEDozgzZrBrP7rLPVGxNB/ZhxHLzj3yqJCWIMsOJbEabopp8TdW7CwL+yYBqXaQJPvjCsbJ2KJLqbnE2/xQ8xBftlwlLv3Na3K5+ad2gXJldnbMiGFXZNuyinjSN2UhS14eEGzMVDrI9j5K/zSVBr/k5Hd14sPGxUj9v0IXq+UlznbThIxLJr/zN/FmSs3zY4nhN2TAuOqlIIaA6DlT3ByG/xYG87vNzuVXQrIlJ6Pm5YgakA4r5TPzbRNx6gxNIpPF8VzPvGW2fGEsFsuX2Bcqg0mOSVaQKffjaWYf6wNh9eYnchu5crszRctSrL63XBeLhXET+sOEzYkiiFL/+LStdtmxxPC7rh8gXHqcTAplacCdF8FGXPC5OawfYrZiexa3mw+fN26NCv616RusUDGxhwkbGgUw1fs4/KNO2bHE8JuuHyBEUmy5DdG/ucPg/m9YOV/XX7k//ME+/vy7WtlWdqnBtVDsvPtqv1U/3I1Xy/fK1c0QiAFRjwqvR+0mwXlO8PaETCrI9y+bnYqu1c4R0bGti/P772rU71gdr5bfYDqQ1bz5ZK/pI1GuDSXLzAu3wbzJHdPeHkE1P8c9iyESY3gqlVm6XE6xYP8+P718izrW4NaRQP5Yc1Bqg9ZzaeL4jkrvc6EC3L5AiNtMMlQCqr0glenwbm/jMb/M7vNTuUwCufIyHevlWVFv5pElsjJpPVHqD40isHzd5Hw9z/W2BPCabl8gRHPUCQSOi+B+3dhQn1j+n+RYiEBvgxvU4bV79akeZlcTN10jJrDohg050+OX5Rbj8L5SYERzxZUxpheJutLMK01bB5vdiKHky9bBoa0LEXUe+G0Ds3D7K0niPgqmgGzdnDkvExBI5yXFBjxfJmCjCuZgvVh8XuwZCDcv2d2KoeTJ6sP/2tekpj3w3m9cj4W7Eig1tfR9Ps1TibVFE7J5QuMNPKnkJcvvDoVKveCTWNh+mvG4EyRajn9vPlvk+LEDoyga/WXWLrrNHVHxNBr2jZZJkA4FZcvMNLInwpu7tDgc2NtmQMrYWJDuHzC7FQOKyBjej5sVIy1AyN4o2Yw0X+dpcE3sfScvIVdJ+UPHuH4XL7AiBdQoRu0mwl/H4XxtSFhu9mJHFo2Xy8GNijCug9q0btWCOsPXuDl79bSZdIfbJelnIUDkwIjXkxIHeiyDNzTGVcyexaancjhZfZJR/96hVk7sBbv1i3EtmOXaD5mPe0nbOKPIxfNjidEqkmBES8usJgxh1lgcfi1PawbKQuYWYCftyfv1C7I2oG1GNigCPEJV2g1dgOvjdvI+oPnkTWchKOQAiPSxjcAOi2CYk1hxX9gYW+4JxM+WoKvlwdvhgcTOzCCjxoV5cC5RNqO30SrsRtYs++cFBph96TAiLTz9DbWlQl7D7b9AlNegRt/m53Kafik86BbWAFi34/g4ybFOfn3DTpM3EyzMetZteeMFBpht1y+wEg3ZQtxc4Pa/4Zm38PR9TChLlw8bHYqp5Le052OVfMTPSCcz5uX5ELiLbr+vIWXv1vL0l2nuX9fCo2wLy5fYKSbsoWVaQsd5kHiWWMOs2MbzU7kdLw83GlbKS9R74UztGUprt26yxtTttJwZCwLdyRwTwqNsBMuX2CEFeSvDt1WGdP//9wYds4yO5FT8nR3o3VoHlb2r8k3bcpw9/593pm+nXojYpi7/QR378l6PsJcUmCEdWQPMYpM7gowpxtEfyk9zKzEw92NZmVzsbxfTUa1LYunuxv9ft1BneExzNxynDtSaIRJpMAI6/HJCu3nQum2EP0FzOkBd2UBLmtxd1O8XCqIxb3DGPt6eTJ4efD+bzuJ+CqaqZuOcuuuzB8nbEsKjLAuDy9oNgZqfQR/zoRfmsK1C2ancmpubooGJXKw6J3qTOwUSjZfLz6cu4vwYdH8vP4IN+9IoRG2IQVGWJ9SUGOA0ZX55Db4sRac22d2KqenlKJWkUDmvVWVX7pUJHcWbwYv2E3Y0Ch+jD3EjdtSaIR1SYERtlOiBXT6HW4lwoQ6cCjG7EQuQSlFjUL+zOxZhendKxPi78tnv++h+pDVjI46wOUbMjBWWIfLFxgZB2NjeSoY08tkzAlTWhgDM4VNKKWoEpyN6T0qM+uNKpTM7cewZXup/uVqhiz9i3NXpX1MWJbLFxgZB2OCLPmh63LIHwYL3oEVg+G+9HSypQr5szKpc0UWvVOdGoX9GRtzkOpDVjN4/i5OXJLlnIVluHyBESZJ7wftZkH5zrDuG5jVAW7LLzZbK5HLj9Fty7Gqf02alcnFtM3HCB8Wzbszd3DgrCwoJ9JGCowwj7snvDwC6n8OexbBpEi4etrsVC6pgL8vQ1qWImZABO2r5OP3PxOoO2INb0zeys4Tf5sdTzgoKTDCXEpBlV7w6jQ4t9dYwOz0LrNTuaygzN4MblycdQNr8XZECOsOnqfJqHW0n7CJDQcvyMSaIlWkwAj7UCQSOi8BfQ8m1od9y81O5NKy+Xrxbr3CrP+gFh80LMKeU1d5bfxGXvl+PSvjZQZnkTJSYIT9CCoD3VdD1gIwvQ1sGmd2IpeXMb0nb9QMZu3ACD5tWpyzV2/R7ZctNBwZy/y4kzLfmXgmKTDCvmQKMq5kCtaHJQNg8ftwXwYEmi29pzvtq+Qn6r1whrcuzd37mj4z4qj1dQzTNh2TaWhEsqTACPvj5QuvToXKb8HmH2BmB7hzw+xUAmMG5xblcrO8bw1+aF+eLD6e/Gvun4QNiWL8mkNcu3XX7IjCjkiBEfbJzR0afAENvoS/fjfmMLt+0exUIombm6J+8RzM61WNKV0rERLgy/8W76HakNWMWLGPS9dumx1R2AEpMMK+VX4TWv0ECdthQj24dNTsROIRSimqF8zOtO6VmftWVSrkz8rIVfupNmQ1//s9njNXbpodUZhICoywf8WbG9P+XztrLMV8aqfZiUQyyubNwvgOoSzrW4N6xQKZsPYwYUOiGDTnT45euGZ2PGECly8wMheZg8hfHbosAzcP+CkSDq42O5F4isI5MvLNq2WJfi+CVqG5mb31BBFfRdN7+nb2nLpidjxhQy5fYGQuMgcSUBS6roDMeWBqK9gxw+xE4hnyZvPhf81LsnZgBN3DCrBqzxkajoyl66Q/2Hr0ktnxhA24fIERDsYvl9GNOW8VmNsT1o6QpZjtXECm9AyKLMq6D2rRr04hth67xCvfr+fVcRtYs++cDNp0YlJghOPxzgyvz4biLWDlf2GJjJVxBJl90tGnTkHWDazFR42Kcvj8NTpM3EyTUetYuusU9+9LoXE2HmYHEOKFeHjBKxOMgZkbRsGVBHjlR/D0NjuZeI4MXh50CytA+yr5mLvtJN/HHOSNKdsI9s/Am+EhNC0ThKe7/O3rDORfUTguNzeo/z+o/0XSWJlmMlbGgXh5uPNqxbys6l+Tb18ri6e7G+/N2kH4sGh+Xn+Em3fkqtTRSYERjq/KW9ByIiRsMybK/PuY2YlEKni4u9GkdBBL+oQxsVMoOfzSM3jBbqoPWc2Y6ANcuSlLOjsqKTDCOZRoYYyVuXoGfqyL79VDZicSqaSUolaRQH57owq/9qhMsSA/hi7dS7UvVjNs2V+cT5QlnR2NFBjhPPJXhy5Lwc2dMnH/goNRZicSL0ApRaUC2fili7Gkc1ih7IyJPki1L1fz73m7OH5RVj51FFJghHMJLAZdV3AzfQBMbQk7fjU7kUiDErn8GNOuPCuTlnSe8ccxwpMGbcYnyKBNeycFRjgfv1zElfk8aaxMDxkr4wSCk5Z0jn2/Fl2rv8SqPWeI/DaWTj9tZuMhWWnTXkmBEU7prqevjJVxQjn80vOvyKKs/6A279UrxJ8nLvPquI20+H49y3aflrE0dkbGwQjnJWNlnJafjydv1ypIt7ACzNpynHGxh+g5eSvB/hl4o2YwTcvkIp2H/P1sNvkXEM5Nxso4tYcrbb4bzshXy+Dp7saA33ZSc1gUP8bKAmhmc8oCo5TKo5SKVkrFK6V2KKVamJ1JmEzGyjg1D3c3mpbJxZI+YUzqXIG8WX347Pc9VP1yNcOX7+WCdHE2hVMWGOAu0FdrXQyoC4xUSvmYnEmY7YmxMpz+0+xEwsKUUoQXDuDXnlWY81ZVKr2UlW9XH6DakNUMni9dnG3NKQuM1vqU1jou6euzwCUgu6mhhH14MFZGucHEhnAo2uxEwkrK5c3CuA6hrOxfg8algpi22eji3HfGdv46LV2cbcHmBUYpVUMptUApdVIppZVSnZLZ5y2l1GGl1E2l1FalVFgazhcKeALH0xBbOJPAYtBtpbGuzJSWsHOm2YmEFYUEZGRYq9KseT+CzlXzszz+DA2+iaXzT5vZfPiidHG2IjOuYHyBXUAf4MaTTyql2gAjgc+BssB6YIlSKu8j++x6yiPPE8fKBvwCdNXyKRKPerCuTJ5KMKc7rP1Gxso4uZx+3nz0cjHWf1CLd+sWYseJy7T+YQMtx25gRfwZ6eJsBTbvpqy1XgwsBlBKTUpml/7AJK31+KTv31FKNQDeBAYlHaPE886jlPIC5gJfaK3XWyC6cDbemaH9HGPhspWDjW7MDb4AN3ezkwkryuyTjndqG12cZ245zrg1h+j+yxYKBvjSs2awLBdgQcrMP+yVUonA21rrSUnfpwOuA69prWc9st9ooITWumYKj6uAacBerfV/n7FfD6AHQGBgYPkZM2QJXmeRmJiIr69vynbW9wk+OIk8J+ZzLnsV9hTtx313L+sGFHbj7n3N5tP3WHzoNicSNVnTKxrk96R8lltk80vhZ8iFRUREbNVahyb3nL0VmCDgJFBTa73mkf3+A7TTWhdO4XGrA2uAnY9sbq+1fmq3odDQUL1ly5bUvwlhl6KjowkPD0/dizaMhmX/MqaYeXUa+GS1SjZhn7TWRO89x/fRB9l85CIZPKFbjYJ0rJqfrBnSmR3Pbimlnlpg7HUk/5NVTyWz7ekv1notTtpDTlhRlV6QMQfMfQMmNoDXf4PMeZ//OuEUlFJEFAkgokgAW49e5LPZmxm5aj/j1hyiTYU8dAt7idxZZLRDatjbL+HzwD0gxxPbA4Az1jihUqqxUmrc5cuXrXF44WhKvAKvz4Grp2WsjAsrny8rfcqlZ0W/GkSWzMmUjUcJHxZN/1/j2Hv6qtnxHIZdFRit9W1gK8bgyEfVxehNZo1zLtRa9/Dz87PG4YUjeilMxsoIAAoGZuTr1qWJeT+CDlXys2TXaep/s4auk/5gyxGZcuh5zBgH46uUKqOUKpN0/rxJ3z+4FzEc6KSU6qaUKqqUGgkEAWNtnVW4MBkrIx6RK7M3/2lsdHHuW6cg245douXYDbT8fj2r9kgX56cx4womFNie9PAGPk76+hMArfWvQF/gIyAOqA5Eaq2PWiOM3CITT/XkWJl1I2WsjIvLkiEdfesUYt0HtRjcuBgJf9+g689baDgyljnbTnDn3n2zI9oVU3uR2RPpReZcXqgX2dPcvWWMldk9Fyr2lLEyLiIln6E79+6zcEcCY2MOsu9MIrkye9Mt7CXaVMiDTzp77UNlWY7Yi0wI++HhBa9MhIxBsHE0XD0FLcaDZ3qzkwmTebq70aJcbpqVyUXU3rN8H32QjxfG8+2q/XSokt/luzhLgREiJdzcoMHnxuJlyz+Eyefh1akyVkYA4OamqF00kNpFA/njyEV+iDnIyFX7+WHNQdqE5qFbWAHyZHW9Ls4uX2CUUo2BxiEhIWZHEY6g6tvGWJl5byaNlZltdAQQIkmF/FmpkD8r+85cZdyaQ0zddIwpm47xcqmc9KwRTLGgTGZHtJkUNfIrpTyVUqeVUsWtHcjWpJuySLWSLf9/rMwEGSsjklcoMCNftSpN7MAIulTLz8r4M0R+G0uHiZtZf+C8S8zinKICo7W+A9whFaPphXBqL4VBlyWAkrEy4ply+nnzYaNirP+gNgPqFyY+4TJtf9xE09Hr+H3nKe45cRfn1HRT/g4YpJRy+dtqQgAQWBy6rQC/3EljZWY9/zXCZfn5eNIrIoS1A2vxv+YluHLjDr2mbaPW19FM2XiUm3fumR3R4lJTYMKApsBJpdSqpEXDHj6slM/qZByMSBO/3Mao/zyVYE43WPet2YmEnUvv6U67SvlY9W44Y9qVw8/bk4/m7aL6kNWMWr2fy9fvmB3RYlJTYM4DszHWcjkGXHji4ZCkDUak2YN1ZYo1gxX/hmUfwn0ZcCeezd1NEVkyJ/N7VWNa90oUD/Ljq+X7qPrlKj5bFE/C3/9Yj9HhpPh2l9a6szWDCOHQPLyg5URY4g8bRsG1c9B0NLh7mp1M2DmlFFWDs1M1ODvxCVcYt+YgP60/wqT1R2haJhc9axagUGBGs2O+kFS3pyilCgDFMBr892itD1k8lRCOyM0dIodBxkBY/RlcOw+tfwEvWbRKpEyxoEx882pZ3q1XmAlrD/PrH8eZve0EtYsE8EZ4MKH5smCsp+gYUnyLTCmVSSk1CzgAzAPmA/uVUjOVUo5ZXoWwNKWgxgBo8h0cioKfGxuFRohUyJPVh/82Kc76D2rRr04hth//m1ZjN/DK9+tZvvu0w0yumZo2mJFAKSACY5JKb6B20rZvLJ7MRqSRX1hFuQ7QZiqcjYcJ9eDSEbMTCQeUJUM6+tQpyLqBtfikaXHOXr1Fj8lbqTsihpl/HOfWXfvueZaaAtME6Ka1jtFa30l6RGOsad/MGuFsQRr5hdUUiYQO8+H6BaPIyIBM8YK807nToUp+ot8L59vXyuLl4c77s3dSY2gUP8Qc5OpN++x5lpoC403yvcUuAjLrnxDJyVvZ6Mbs5gE/RcLhWLMTCQfm4e5Gk9JB/N67Or90qUhIgC9fLPmLql+s5sslf3H2yk2zIz4mNQVmHfCpUurhjG1KqQwY67lYZbVJIZxCQFHouhwy5oQpLSB+vtmJhINTSlGjkD9Tu1VmwdvVqFHYn3FrDlJ9SBQfzN7JwXOJZkcEUldg+gOVMAZaxiilooETSdv6Wj6aEE7kwYDMnGVgZkf440ezEwknUSp3Zka3Lcfqd8NpXSE3c7efpM7wGHpO3sL2Y5dMzZaacTB/KqUKAq8DRQAFTAGmaq0df0SQENbmk9Vok/mtM/z+LiSehfBBRs8zIdIof/YMfNasJH3rFOLn9Uf4ZcNRlu0+Q8WXsvJmzWDCC/vbvItzigqMUsoTOA7U1lqPt24kIZxYOh+jd9nCPhAzBBLPQOTX4C5T/AnLyO7rxbv1CvNGzWBm/HGcCbGH6DzpDwoHZqRnzQI0Lh2Ep3tqbl69OJefTVm6KQubc/eApqOgen/YOglmdYQ7chNAWFYGLw+6Vn+JmPcjGN66NAD9Z+6g5tAoJqw9zLVbd62eweVnU5ZuysIUSkGdwdBwKPz1O0xuATf+NjuVcEIPlnVe2jeMnzpVIHdWHz5dFE/VL1fz9fK9nE+8ZbVzp6ZYhAE1MRr5dwHXHn1Sa93EksGEcAmVekIGf5jTA35qaKyQmSnI7FTCCSmliCgSQESRALYdu8QPMQcZFXWAcWsO8XnzkrxSPrfFz5maAvNgNmUhhCWVaGF0AJjxujEg8/U54F/I7FTCiZXLm4Uf2ody8Fwi49ccstoyzilt5PcA/gDmaa0TrJJECFdWIBw6LYKpLWFiPWg7C/JUMDuVcHLB/r58+Uopqx0/pY38d4FhgMw9LoS1BJUxBmSmz2xMkrlvudmJhEiT1DTybwTKWyuIEALIWsAoMv6FYPqrEDfN7ERCvLDUtMGMB75SSuUFtvLPRv5tlgwmhMvyDYBOv8OMdjDvTWNAZrU+MiBTOJzUFJgHf0oNT+Y5DbinPY7tKaUaA41DQkLMjiLE//PKCO1mwdw3YOVgY0Bmvf+Bm20GyAlhCakpMC9ZLYWJtNYLgYWhoaHdzc4ixGM8vOCVCcYVzcYxxpVMs+/BI53ZyYRIkRT/OaS1PoqxVPJoYAlwP2lbXUD+/BfCGtzcoMGXUHsw7PoNprWGW1fNTiVEiqRmyeR2wExgP8bVzIMeZe7A+5aPJoQAjLaXsP7QdAwcXgOTXobEc2anEuK5UnND932gu9a6H/DoJDYbgTKWDCWESEbZdvDadDi31xgrc/Gw2YmEeKbUFJiCwIZkticC1hkGKoR4XKH60HEB3LhkjPo/tdPsREI8VWoKTAKQ3PwVNYCDlokjhHiuPBWhyzJwT5e0DPMasxMJkazUFJhxwLdKqWpJ3+dRSnUEhgLfWzyZEOLp/AsbAzL9csOUV2D3XLMTCfEPqVnRcqhSyg9YAaQHooBbwFda69FWyieEeBq/XNBlCUx7FWZ1Nhr+K/UwO5UQD6Vq1JbW+kMgO1ARqAz4a63/bY1gQogU8M4CHeZB4YawZACs+hS0060LKBxUqocFa62va623aK03a60TrRFKCJEKnt7QejKU6wCxX8GCd+Ce9VcrFOJ5nGp1SiFclrsHNP4WfANhzTC4fsGYBSCdj9nJhAtz+YmNlFKNlVLjLl++bHYUIdJGKaj1EUR+BXuXwOTmcP2i2amEC3P5AqO1Xqi17uHn52d2FCEso2J3aDUJErYZyzBfPml2IuGiXL7ACOGUijeD12fDlQRjQOa5vWYnEi5ICowQzuqlGsa6MvfvwMT6cHyz2YmEi5ECI4Qzy1nKGJDpnQV+bgJ7l5qdSLgQKTBCOLss+aHLcggoAjPawvYpZicSLkIKjBCuwNcfOi40bpvN7wWxw2VAprA6KTBCuAqvjNB2JpRsBas+hqUfwP37ZqcSTkwGWgrhSjzSQfNxkCEANo6Ga+eg2VhZhllYhRQYIVyNmxvU/x/4BsDKwcbaMq0ng5ev2cmEk5FbZEK4IqWgel9oOhoOxcDPjeHaebNTCScjBUYIV1b2dXh1KpyNN8bK/H3M7ETCiUiBEcLVFW4I7ecZ7TET6sGZeLMTCSchBUYIAfmqQOclxtc/NYBjG83NI5yCFBghhCGwOHRZBhn84ZemMupfpJlTFhil1Hql1A6l1C6l1H/MziOEw8iSzygyAUWNUf9x08xOJByYUxYYoIHWujRQGmiolCpjch4hHEeG7Emj/sNg3puwbqTZiYSDcsoCo7W+kvRluqSHECI1vDJC21lQvAWs+A8s/0hG/YtUs2mBUUrVUEotUEqdVEpppVSnZPZ5Syl1WCl1Uym1VSkV9oLn2gScBVZqrePSllwIF+SRzlh2uWIPWP8dzH8L7t0xO5VwILYeye8L7AJ+SXo8RinVBhgJvAWsTfrvEqVUMa31saR9dj3l2A211scffKO1rqSUygTMUkqV0Fo/7XVCiKdxc4OGQ42pZaI+M5ZgbjUJ0vmYnUw4AKVNmlFVKZUIvK21nvTItk3ATq1190e27Qd+01oPesHzDATuaa2/Sua5HkAPgMDAwPIzZsx4kVMIO5SYmIivr0x9Ykk5E5ZRaN9YrmQqyJ8l/81dz4xmR7Iq+QylTERExFatdWhyz9nNXGRKqXRAeeDJQrAcqJqK42QGPLTW55VS6YF6wIjk9tVajwPGAYSGhurw8PDUBxd2KTo6Gvn3tLRwiK+C3+xuVN/3Gbw+B/xymR3KauQzlHb21MifHXAHzjyx/QyQIxXHyQosU0rtBLYAMVrrRZaJKISLK9YEXp8NVxKMUf/n9pmdSNgxeyowDzx5z04ls+3pL9b6kNa6vNa6lNa6hNb6k2ftr5RqrJQad/ny5RfJKoTreSkMOv0O924b85ed2GJ2ImGn7KnAnAfu8c+rlQD+eVVjMVrrhVrrHn5+ftY6hRDOJ2cp6LoM0mcyZmI+sNLsRMIO2U2B0VrfBrYCdZ94qi6w3vaJhBDPlLUAdFkO2YJhWhvYOcvsRMLO2HocjK9SqkzSyHo3IG/S93mTdhkOdFJKdVNKFVVKjQSCgLFWzCS3yIR4URkDjdtleavAnG6w8XuzEwk7YusrmFBge9LDG/g46etPALTWvwJ9gY+AOKA6EKm1PmqtQHKLTIg0Su8H7X6Doo1h6Qew8mMwafiDsC827aastY7GaLR/1j5jgDE2CSSEsAzP9NDqZ/i9P6wdbqwt8/I34G43IyGECeRfXwhhGW7uRlHJEABrhhqj/ltOAE9vs5MJk9hNI79ZpA1GCAtSCmp9CA2Hwd7FMLkF3Pjb7FTCJC5fYKQNRggrqNTDuHo58QdMagRXT5udSJjA5QuMEMJKSrwC7WbCxcMwoS5cOGh2ImFjUmCEENYTXAs6LYLb14ypZRK2m51I2JDLFxhpgxHCynKVMwZkevrApJfhULTZiYSNuHyBkTYYIWwgewh0XQ6Z88LUVrB7rtmJhA24fIERQthIppzQeTHkKg+zOsPm8WYnElYmBUYIYTveWaD9XCjUABa/B1FfyKh/JyYFRghhW57e0GYKlHkdYr40Rv/fv2d2KmEFMpJfCGF77h7QdBRkyA7rvoHrF6DFePDwMjuZsCCXLzBKqcZA45CQELOjCOFalIK6H4NvACz7lzG1zKvTjDVmhFNw+Vtk0otMCJNV6QXNx8GxDcao/8SzZicSFuLyBUYIYQdKt4HXZsCFA8aAzIuHzU4kLEAKjBDCPhSsCx0WwM2/YWJ9OP2n2YlEGkmBEULYjzwVoMsycPOAnyLhyDqzE4k0kAIjhLAv/oWNUf8Zc8Lk5rBnkdmJxAty+QIjc5EJYYf8ckOXpZCzFMxsD1t/NjuReAEuX2CkF5kQdsonK3SYD8G1YWFvWPOVjPp3MC5fYIQQdixdBnhtOpRqA6s/haWD4P59s1OJFHL5gZZCCDvn7gnNxoJPNtg4Bq6fh6ZjwCOd2cnEc0iBEULYPzc3qP85ZPCHVR/DjUvQ+hfjCkfYLblFJoRwDEpBWH9o/C0cXA0/NzGmlxF2SwqMEMKxlO8IrScbAzEnNoDLJ8xOJJ5CCowQwvEUfRnaz4Grp2BCfTi3z+xEIhkuX2BkHIwQDip/dej0O9y7bUwtc2Kr2YnEE1y+wMg4GCEcWM5S0HWZMcX/z43hwCqzE4lHuHyBEUI4uKwFjPnLshaAaW3gz9/MTiSSSIERQji+jDmg0yLIUxFmd4NN48xOJJACI4RwFt6Z4fXZUDgSlgyAqM9lahmTSYERQjgPT29jAGbZ1yFmCPzeH+7fMzuVy5KR/EII5+LuAU1GGaP+146A6xegxXjw8DI7mcuRAiOEcD5KQZ3/gk92WP6hMbXMq9PAK6PZyVyK3CITQjivqm9D8x+MlTEnvQyJ58xO5FKkwAghnFvpV40p/8/tNQZkXjpqdiKXIQVGCOH8CtU3Fi+7fh4m1IMzu81O5BKkwAghXEPeStB5qdE+81NDOLbR7EROz+ULjMxFJoQLCSwGXZcbPcx+aQp7l5qdyKm5fIGRuciEcDGZ8xpTywQUhRltIW662YmclssXGCGEC8qQHTouNGZknvcGrP/O7EROSQqMEMI1eWWEdrOgWDNY/hGs+I9MLWNhMtBSCOG6PLyg5URYnA3WjYRrF6DxSGM2AJFm8lMUQrg2N3do9LXR8B/zJdy4aBQdkWZyi0wIIZSCiEEQ+RXsXQKTW+BxJ9HsVA5PCowQQjxQsTu0nAAn/qBM3Idw9bTZiRyaFBghhHhUiVeg3Uy8b5w2Rv1fOGh2IoclBUYIIZ4UXIu4Mp/BravG/GWndpidyCFJgRFCiGRczVTQGJDpkR5+agSHY82O5HCkwAghxNP4FzKKjF8umNIC4heYncihSIERQohn8csFnZdAzjIwqyNsnWR2IochBUYIIZ7HJyt0mAfBtWFhH1jzlYz6TwEpMEIIkRLpMhgLl5VqA6s/haWD4P59s1PZNacdya+UcgM2AUe11i3NziOEcALuntBsLPhkg41jjAXMmo4Bj3RmJ7NLTltggDeBgzj3exRC2JqbG9T/3JhaZtXHcOMStP7FuMIRj3HKW2RKqQCgBTDO7CxCCCekFIT1h8bfwsHV8HMTuH7R7FR2x6YFRilVQym1QCl1UimllVKdktnnLaXUYaXUTaXUVqVU2Aucahjwb0BukAohrKd8R2g9GU7/CRMbwOUTZieyK7a+gvEFdgF9gBtPPqmUagOMBD4HygLrgSVKqbyP7LPrKY88Sc/XALTWer0N3o8QwtUVfRnaz4Grp2BCfTi3z+xEdkNpk7raKaUSgbe11pMe2bYJ2Km17v7Itv3Ab1rrQSk87gfAO8AdID2QEZihte6azL49gB4AgYGB5WfMmPHib0jYlcTERHx9fc2OIRxYaj9DvlcPUWrnxyh9j52l/sPVTIWsmM5+REREbNVahyb3nN0UGKVUOuA68JrWetYj+40GSmita77AOcKTzvHcXmShoaF6y5YtqT2FsFPR0dGEh4ebHUM4sBf6DF08BJObQ+I5eHUKBNeySjZ7opR6aoGxp0b+7IA7cOaJ7WeAHLaPI4QQqZS1AHRZbvx3amvYNdvsRKaypwLzwJOXVCqZbSk7kNbRz7t6UUo1VkqNu3z58oucQgghHpcxEDotgtwV4LeusHm82YlMY08F5jxwj39erQTwz6sai9FaL9Ra9/Dz87PWKYQQrsY7s9HwX7ghLH4Por5wyall7KbAaK1vA1uBuk88VRejN5kQQjgOT2+jC3OZ1yHmS6PQ3L9ndiqbsukod6WULxCS9K0bkFcpVQa4qLU+BgwHJiulNgPrgDeAIGCsFTM1BhqHhIQ8d18hhEgVdw9oOgoyZIN1I+H6BWj+A3h4mZ3MJmx9BRMKbE96eAMfJ339CYDW+legL/AREAdUByK11ketFUhukQkhrEopqPsJ1P0Uds+Faa3hVqLZqWzCplcwWutojEb7Z+0zBhhjk0BCCGEr1Xobk2QueAd+bgztfjOubJyY3bTBCCGE0yvbDl6dCmfjYWJ9+Pu42YmsyuULjHRTFkLYVOGG0H4uJJ41iszZv8xOZDUuX2CkDUYIYXP5qkLnxXD/LvzUAI7/YXYiq3D5AiOEEKbIUQK6LIP0meGXJnBgpdmJLE4KjBBCmCXrS9B1OWQLhmlt4M/fzE5kUVJghBDCTL4B0Ol3yFMZZneDTT+YnchiXL7ASCO/EMJ06f3g9dlQpBEseR9W/88pppZx+QIjjfxCCLvgmR5a/Qxl28OaofB7f4efWsamAy2FEEI8g7sHNPkOMmSHtSOMqWVajHfYqWWkwAghhD1RCur8F3yyw/IP4cYleHUaeGU0O1mqufwtMiGEsEtV34ZmY+HIOpj0Mlw7b3aiVHP5AiON/EIIu1XmNXhtOpzbmzS1zDGzE6WKyxcYaeQXQti1QvWhwzy4dg4m1IOze8xOlGIuX2CEEMLu5a0MnZcYXZcnNoBjm8xOlCJSYIQQwhEEFjdG/ftkhV+awv4VZid6LikwQgjhKLLkgy7Lwb8QTH8Vds40O9EzSYERQghH4usPHRdB3iowpzts/N7sRE/l8gVGepEJIRxO+kzGiphFG8PSD2DVJ3Y5tYzLFxjpRSaEcEgPppYp1xFiv4aFfexuahkZyS+EEI7KzR0aj4QM/hD7Fdy4CC1+NIqPHXD5KxghhHBoSkHtf0ODL2HPQpjaEm5eMTsVIAVGCCGcQ+U3jYkxj22An1+GxHNmJ5ICI4QQTqNUa3htBpzbBxPrwaUjpsaRAiOEEM6kYF3ouACuX4QJ9eHMbtOiSIERQghnk6cidFlqtM/81BCObjAlhhQYIYRwRgFFjallMvjD5Gawd6nNI7h8gZGBlkIIp5U5L3RZBv5FYEZbiJtu09O7fIGRgZZCCKeWITt0WgT5q8O8N2D9KJud2uULjBBCOD2vjNBuFhRraizDvGKwTaaWkQIjhBCuwMMLWv4E5TvDum9gwdtw7651T2nVowshhLAfbu7w8gij4X/NULjxN7wywWpTy8gVjBBCuBKloNaH0HAo/LUIprwCN63TyUkKjBBCuKJKPY2rl+MbYdtkq5xCbpEJIYSrKtkSsheCwBJWObwUGCGEcGU5S1nt0HKLTAghhFVIgRFCCGEVLl9gZKoYIYSwDpcvMDJVjBBCWIfLFxghhBDWIQVGCCGEVUiBEUIIYRVSYIQQQliFFBghhBBWIQVGCCGEVUiBEUIIYRVSYIQQQliF0jZYNtMRKKXOAUcf2eQHpHR4f2r2tcTrrHUcezmPJWQHzj/luWe9j5S8R0vtk5b90/o6ax3H7HNY0tM+Q897H2l9PrX7pXZfS7zuUfm01v7JPqO1lkcyD2CcNfa1xOusdRx7OY+Fsm55kfeRkvdoqX0s8bN1pM+QI31+kvIm+xl63vtI6/Mv8vMy+/PztIfcInu6hVba1xKvs9Zx7OU81vas95GS92ipfdKyf1pfZ63jmH0OW3je+0jr86ndL7X7WuJ1KSK3yIRTUkpt0VqHmp1DOC75DKWdXMEIZzXO7ADC4clnKI3kCkYIIYRVyBWMEEIIq5ACI4QQwiqkwAghhLAKKTBCCCGsQgqMEIBSaoFS6pJS6jezswjHopTKo5SKVkrFK6V2KKVamJ3JXkgvMiEApVQE4At01Fq3NDuPcBxKqZxAoNY6TikVAGwFCmutr5sczXRyBSMEoLWOAq6anUM4Hq31Ka11XNLXZ4FLGPOYuTwpMMKuKaVqJN2+OqmU0kqpTsns85ZS6rBS6qZSaqtSKsyEqMIO2frzo5QKBTyB42mI7TQ8zA4gxHP4AruAX5Iej1FKtQFGAm8Ba5P+u0QpVUxrfSxpn11POXZDrbX8InBuNvv8KKWyJZ2jq5a2B0DaYIQDUUolAm9rrSc9sm0TsFNr3f2RbfuB37TWg1J5/PCk40sbjBOy5udHKeUFrADGa60nWy61Y5NbZMJhKaXSAeWB5U88tRyoavtEwpFY6vOjlFLAJGC1FJfHSYERjiw74A6ceWL7GSBHag6klFoJzAIilVInlFJVLBNR2DFLfX6qAW2AZkqpuKRHSQtldGjSBiOcwZP3eVUy2559AK3rWC6OcDBp+vxordcif6wnS34owpGdB+7xz782A/jnX6VCPEk+P1YmBUY4LK31bYxBbXWfeKousN72iYQjkc+P9cktMmHXlFK+QEjSt25AXqVUGeBiUjfS4cBkpdRmYB3wBhAEjDUhrrAz8vkxl3RTFnYtqetwVDJP/ay17pS0z1vA+0BOjDEP/bTWa2wUUdgx+fyYSwqMEEIIq5A2GCGEEFYhBUYIIYRVSIERQghhFVJghBBCWIUUGCGEEFYhBUYIIYRVSIERQghhFVJghLBDSqlopdQos3MIkRZSYIQQQliFFBghhBBWIQVGCPvlppT6XCl1Xil1Vin1lVJK/p8VDkM+rELYr3bAXYzle98G+mKsnCiEQ5DJLoWwQ0qpaMBLa13lkW0rgKNa626mBRMiFeQKRgj7tfOJ7xMwVlsUwiFIgRHCft154nuN/D8rHIh8WIUQQliFFBghhBBWIQVGCCGEVUgvMiGEEFYhVzBCCCGsQgqMEEIIq5ACI4QQwiqkwAghhLAKKTBCCCGsQgqMEEIIq5ACI4QQwiqkwAghhLCK/wPGgoEIW4zwngAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(6, 6))\n",
    "\n",
    "\n",
    "ax.loglog(h, err1, label=\"$\\mathcal{O}(h)$\");\n",
    "ax.loglog(h, err2, label=\"$\\mathcal{O}(h^2)$\");\n",
    "ax.invert_xaxis()\n",
    "ax.grid(\"both\")\n",
    "ax.legend()\n",
    "# ax.set_aspect(\"equal\")\n",
    "\n",
    "ax.set_xlabel(\"h\")\n",
    "ax.set_ylabel(\"error\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Related: error handling \n",
    "_sometimes the coude should \"break\"_ \n",
    "\n",
    "Then it is important to give users (including yourself!) useful error messages to diagnose the problem. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def moving_average(data, window_size):\n",
    "    \n",
    "    # checking types\n",
    "    if not isinstance(data, (np.ndarray, list)):\n",
    "        raise Exception(\n",
    "            f\"data must be a numpy.ndarray or a list. the input is type: {type(data)}\"\n",
    "        )\n",
    "    \n",
    "    if not isinstance(window_size, int):\n",
    "        raise Exception(\n",
    "            f\"window_size must be an int. The provided input type is {type(window_size)}\"\n",
    "        )\n",
    "    \n",
    "    # checking valid inputs\n",
    "    if window_size <= 0:\n",
    "        raise Exception(\n",
    "            f\"window_size must be strictly positive. The input value is {window_size}\"\n",
    "        )\n",
    "    \n",
    "    data = data.squeeze() # remove single-dimensional entries \n",
    "    if data.ndim > 1:\n",
    "        raise Exception(\n",
    "            f\"The current implementation only supports 1-dimensional arrays. The provided dimension is {data.ndim}\"\n",
    "        )\n",
    "        \n",
    "    average = np.full(data.size, np.nan)\n",
    "    half_window = window_size // 2\n",
    "    for i in range(half_window, data.size - half_window):\n",
    "        average[i] = np.mean(data[i - half_window: i + window_size-half_window])\n",
    "    return average"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Other ideas for testing \n",
    "\n",
    "## Examples as tests & consistency over time\n",
    "https://arxiv.org/pdf/1508.07231.pdf\n",
    "\n",
    "- save some result, check that you can still reproduce them if you make changes\n",
    "- visual checks, save and re-create figures\n",
    "\n",
    "## Realm of reasonable solutions\n",
    "- should the solution or inputs be strictly positive?\n",
    "    - throw random numbers at it, are the results physical\n",
    "    - what should your code do with non-sense\n",
    "- conserved quantities (e.g. conservation of energy)\n",
    "    - if not, then how does it behave"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}