Lecture_08.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import norm\n",
    "from scipy.spatial import distance_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 align=\"center\"><b><u>Lecture 8 - Hilbert Spaces</u></b></h1>\n",
    "\n",
    "<h2>Reproducing Kernel Hilbert Space</h2>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lecture 8 introduces the idea of function approximation through the use of Reproducing Kernel Hilbert Space (RKHS).  This example reflects the regularised empirical error minimisation problem seen at the end of the lecture.  The derivations are shown here for a Radial Basis Function (RBF) kernel, and then the example is evaluated using a Gaussian kernel function.<br/>\n",
    "\n",
    "The dataset is made from pairs of 1D inputs and outputs $(x_1,y_1), (x_2,y_2), ... ,(x_N,y_N) \\in \\mathcal{X} \\times \\mathbb{R}$<br/> \n",
    "\n",
    "The regularised empirical error is given by:\n",
    "$$\n",
    "\\frac{1}{N} \\sum_{i=1}^{N} | f(x_i) - y_i |^2 + \\lambda \\left\\Vert f \\right\\Vert^2\n",
    "$$\n",
    "\n",
    "and is minimised by functions of the form:\n",
    "\n",
    "$$\n",
    "\\hat{f}(\\cdot) = \\sum_{n=1}^{N} \\alpha_n k(\\cdot, x_n)\n",
    "$$\n",
    "\n",
    "Where $k(\\cdot,\\cdot)$ is a general reproducing kernel function.  The coefficients can be found by evaluating:\n",
    "\n",
    "$$\n",
    "\\underset{\\alpha \\in \\mathbb{R}^N}{\\mathrm{argmin}} \\frac{1}{N} \\sum_{i=1}^{N} | f(x_i) - y_i |^2 + \\lambda \\left\\Vert f \\right\\Vert^2\n",
    "$$\n",
    "\n",
    "Defining the following vectors and matrix allows the error to be defined in vector form:\n",
    "\n",
    "$$\n",
    "\\alpha = \n",
    "\\begin{bmatrix}\n",
    "\\alpha_1\\\\\n",
    "\\alpha_2\\\\\n",
    "\\vdots\\\\\n",
    "\\alpha_N\\\\\n",
    "\\end{bmatrix}\n",
    "\\text{    ,    }\n",
    "y = \n",
    "\\begin{bmatrix}\n",
    "y_1\\\\\n",
    "y_2\\\\\n",
    "\\vdots\\\\\n",
    "y_N\\\\\n",
    "\\end{bmatrix}\n",
    "\\text{    ,    }\n",
    "K_{ij} = k(x_i,x_j)\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\alpha^* = \\underset{\\alpha \\in \\mathbb{R}^N}{\\mathrm{argmin}}\\frac{1}{N}  \n",
    "(K\\alpha-y)^T(K\\alpha-y) + \\lambda \\left\\Vert f \\right\\Vert^2\n",
    "$$\n",
    "\n",
    "The error here is convex, and so the solution for alpha is found by taking the derivative and setting to zero. First the $L_2$ norm of $f$ needs to be rewritten.<br/>\n",
    "\n",
    "$$\n",
    "\\left\\Vert f \\right\\Vert^2 = \\langle f,f \\rangle = \\sum_{i}\\sum_{j} \\alpha_i\\alpha_j\n",
    "\\langle k(\\cdot,x_i),k(\\cdot,x_j) \\rangle = \\sum_{i}\\sum_{j} \\alpha_i\\alpha_j k(x_i,x_j) = \n",
    "\\alpha^T K \\alpha\n",
    "$$\n",
    "\n",
    "The third equality here is the significant step, and uses the reproducing property of the kernel function. i.e. $f = \\langle f, k \\rangle$.\n",
    "\n",
    "$$\n",
    "\\frac{\\partial}{\\partial \\alpha} \\frac{1}{N}  \n",
    "(K\\alpha-y)^T(K\\alpha-y) + \\lambda \\alpha^T K \\alpha = 0 \n",
    "$$\n",
    "$$\n",
    "\\frac{1}{N}K (K\\alpha - y) + \\lambda K \\alpha = 0\n",
    "$$\n",
    "\n",
    "The linear system that must then be solved is given below, this step is legitimate since $K$ must be positive semi-definite (positive definite if there are no repeated $x_i$ in the dataset).\n",
    "\n",
    "$$\n",
    "(K + \\lambda N I )\\alpha = y\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dataset for the example will be generated using a generating function, with some small amount of noise so that the final approximating function can be compared to the generating function to see how good the approximation is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "###--- Data generation ---###\n",
    "sigma = 0.1    # The small added noise\n",
    "N = 10         # Number of training points\n",
    "\n",
    "# Generate data with some small error - uncomment the desired function\n",
    "g = lambda x: 2*x + np.sin(10*x) + 2      # Sine wave - l=0.1 works well, smaller length scale\n",
    "#g = lambda x: x                          # Ramp function - l=1 works well, larger length scale"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The generating function is defined as:\n",
    "\n",
    "$$\n",
    "g(x) = \\sin\\left(10x\\right) + 2x + 2\n",
    "$$\n",
    "\n",
    "and is evaluated at the training points $\\{x_i\\}_{i=1}^{N}$, giving $\\{y_i\\}_{i=1}^{N}$, with added noise to generate the data.  For plotting and further inspection of the function, it is also evaluated  at $\\{d_j\\}_{j=1}^{J}$, the evaluation points, with $J>>N$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Training points\n",
    "x = np.linspace(0,1,N)\n",
    "y = g(x) + norm.rvs(size=N, loc=0,scale=sigma)\n",
    "\n",
    "# Evaluation points\n",
    "d = np.linspace(0,1,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEGCAYAAAB1iW6ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3hUZfbA8e9LE0GKC8iiQIJdeglNXZpYVhALqCsoigVQRF1FBXEFVFYB0QXpiooaBEWw8ENXkaorSEA6CqgBQcRQRDSUQM7vjxNiCCmTZGbuzJ3zeZ48mXIz99yUk3fecl4nIhhjjIl+xbwOwBhjTHBYQjfGGJ+whG6MMT5hCd0YY3zCEroxxvhECa9OXLlyZYmPj/fq9MYYE5WWL1++S0Sq5PScZwk9Pj6epKQkr05vjDFRyTm3JbfnrMvFGGN8whK6Mcb4hCV0Y4zxCc/60HOSlpbGtm3bOHjwoNehmAhTunRpqlevTsmSJb0OxZiIFVEJfdu2bZQrV474+Hicc16HYyKEiLB79262bdtGrVq1vA7HmIgVUV0uBw8epFKlSpbMzXGcc1SqVMneuRmTj4hK6IAlc5Mj+70wJn8R1eVijPE/EfjhB0hJgX374MABOPNMOO88KFXK6+iiW8S10L22c+dOunbtyplnnkmTJk1o2bIls2bN8iyeBQsW8L///S/z/oQJE3j99deD8tqLFy+mTp06NGzYkAMHDgTlNQH+/e9/H3f/wgsvDNprmyiSmAjx8VCsGOlxtfj4kXn06aMPnXUWtGgBl18O11wD9etD2bJQrx48+SRs3+518FFKRDz5aNKkiWS3fv36Ex4Lp/T0dGnRooWMHz8+87Hk5GQZPXp0SM+blpaW63ODBg2SESNGhOS8vXr1kldeeSXor1u2bNmgv6aI978fpgDefFOkTBlJB/mIy6UhKwREyp50WK65RmTCBJHZs0UWLxZZtkxk6lSRxx4TadtWBESKFxe5+mqRVau8vpDIAyRJLnnVEnoWc+fOlVatWuX6/JEjR6Rfv36SkJAg9erVkwkTJoiIyPz586V169bSuXNnOe+886Rr166Snp4uIiJJSUnSqlUrady4sVx22WXy008/iYhI69atZcCAAdKqVSt57rnn5IMPPpBmzZpJw4YN5ZJLLpGff/5ZfvjhB6lataqcfvrp0qBBA1m0aNFxCb5169byyCOPSNOmTeWcc86RRYsWiYjIH3/8Iddff73Uq1dPbrjhBmnWrJksW7bsuGt56aWX5NRTT5X4+Hjp2rWrzJ8/Xzp06JD5fJ8+feTVV18VEZG4uDh54oknpFGjRlK3bl3ZsGGDiIjs379fbrvtNqlbt67Uq1dPZsyYIY8++qgUK1ZMGjRoIF27dhWRPxN8enq69OvXT+rUqSN169aVadOm5fv9y8rr3w9TAHFx8h21pC2fCYjU4jt5g25yoOa5+X7p5s0i/fuLVKqkib1fP5H9+8MQc5TIK6FHbB/6Aw/AypXBfc2GDeE//8n9+XXr1tG4ceNcn588eTIVKlRg2bJlHDp0iIsuuojLLrsMgK+//pp169Zx+umnc9FFF/HFF1/QvHlz+vbty/vvv0+VKlWYPn06AwcO5JVXXgHg119/ZeHChQDs3buXJUuW4Jzj5ZdfZvjw4YwcOZLevXtzyimn0K9fPwA+++yz42I6cuQIX331FXPmzGHIkCHMnTuXcePGceqpp7J69WrWrl1Lw4YNT7iWO++8k88//5yOHTvSpUsXFixYkOf3rnLlyqxYsYJx48bx3HPP8fLLL/PUU09RoUIF1qxZk3kNnTt3ZsyYMazM4Yc3c+ZMVq5cyapVq9i1axdNmzalVatWuX7/Lr744jxjMpHrwy31uYXXcQij6UsvJlKKNPgx/8Hts86CZ56Bhx+G/v3huedg+nTtwfnb38IQfBQLOKE754oDScB2EemY7bnbgBHAsZ6vMSLycrCC9EqfPn34/PPPKVWqFMuWLeOTTz5h9erVzJgxA4B9+/axadMmSpUqRbNmzahevToADRs2JDk5mYoVK7J27VouvfRSAI4ePUq1atUyX//GG2/MvL1t2zZuvPFGduzYweHDhwOeb33dddcB0KRJE5KTkwH4/PPPuf/++wGoW7cu9evXL9o3Itt5Zs6cCcDcuXOZNm1a5jGnnnpqnq/x+eefc9NNN1G8eHGqVq1K69atWbZsGeXLl8/x+2cJPfocPQpPPAH/5gMas5wZdKEWyX8eULNmwK/1l7/ApElw663QowdccgmMHw933BH8uP2iIC30+4ENQPlcnp8uIvcWPSSVV0s6VOrUqcO7776beX/s2LHs2rWLhIQEQLunXnzxRS6//PLjvm7BggWcdNJJmfeLFy/OkSNHEBHq1KnDl19+meP5ypYtm3m7b9++PPjgg3Tq1IkFCxYwePDggGI+dt5j5zwWZ0GVKFGC9PT0zPvZ53zndp6CTCfMK66cvn8muqSnw513wmuvwV1tNzF6yaWUPrD3zwPKlIGhQwv8uhddBEuXwo036uuvWaOt9hIR27/gnYBmuTjnqgMdgKhvdeelXbt2HDx4kPHjx2c+lpqamnn78ssvZ/z48aSlpQGwceNG/vjjj1xf77zzziMlJSUzoaelpbFu3bocj923bx9nnHEGAFOmTMl8vFy5cuzfv79A13HxxRfz9ttvA7B+/frMLpG8xMXFsX79eg4dOsS+fftO6NrJyWWXXcaYMWMy7+/dq3+8JUuWzPweZdWqVSumT5/O0aNHSUlJYdGiRTRr1izQyzIRLD0d7rpLk/mQITBp3jmUfulFiIsD5/TzpEnQrVuhXv/UU2HOHO2KHTUKunfXdwPmeIFOW/wP8AiQnscxnZ1zq51zM5xzNXI6wDnX0zmX5JxLSklJKWisIeec47333mPhwoXUqlWLZs2aceuttzJs2DBA+51r165N48aNqVu3Lr169cqzJVmqVClmzJjBo48+SoMGDWjYsOFxUxCzGjx4MNdffz1/+9vfqFy5cubjV111FbNmzaJhw4YsXrw4oOu45557SElJoX79+gwbNoz69etToUKFPL+mRo0a3HDDDdSvX59u3brRqFGjfM/z+OOPs3fvXurWrUuDBg2YP38+AD179sx8nayuvfZa6tevT4MGDWjXrh3Dhw/nr3/9a0DXZCJXejr06gWvvAKDBmmXC6DJOzlZD0hOLnQyP6ZECXjhBXj2WXjrLe2GsaSeTW6jpcc+gI7AuIzbbYDZORxTCTgp43ZvYF5+rxuJs1z84siRI3LgwAEREdm8ebPExcXJoUOHPI6q6Oz3IzI98YTOl3v8cZEcJieFxNNP6zl79BA5ejQ854wUFHGWy0VAJ+fclUBpoLxz7k0RuTnLP4XdWY5/CRhW1H80pvBSU1Np27YtaWlpiAjjx4+nlC3BMyEwc6YuBOrRQz+Hq0LDwIGQlqbdO5UqwYgR4TlvpMs3oYvIAGAAgHOuDdAvazLPeLyaiOzIuNsJHTw1HilXrpxt72dCbu1a7ctu3lxnn4S73M6gQbBrlw6Q1q6t/1RiXaHHiZ1zT6JN/w+A+5xznYAjwB7gtuCEZ4yJRHv3wtVXQ/ny2krPMkkpbJzT2XDffqt9+GefbfPUC5TQRWQBsCDj9hNZHs9sxRtj/O/ee2HrVli0CE4/3bs4SpSAt9+Gli3h2mth2TKI5ZL5VpzLGFMgM2bA1Kk6m6VlS6+j0SmNH34IR47oXPXDh72OyDuW0I0xAfv5Z+jdG5o2hQER9J78nHN02uSyZTpgGqssoWcT6eVzAxUfH8+uXbvyPOa1117j3nvzXtxb2PMb/xGBnj3hjz/g9dcjb6XmddfBPffoIOlHH3kdjTcsoWchIlxzzTW0atWK77//nuXLlzNt2jS2bdsW0vPmtTjJ64Tq9flN5HjrLe3aeOYZOP98r6PJ2ciRWlu9e3f46Sevowm/6E7oWQroEx+v94tg3rx5lCpVit69e2c+FhcXR9++fQEtrvXwww/TtGlT6tevz8SJEwFNem3atKFLly6cf/75dOvWLbNuyfLly2ndujVNmjTh8ssvZ8cOnd3Zpk0bHnvsMVq3bs2oUaP48MMPad68OY0aNaJ9+/bs3LmT5ORkJkyYwAsvvJC5UjQlJYXOnTvTtGlTmjZtyhdffAHA7t27ueyyy2jUqBG9evXKtW7Kq6++yrnnnkvr1q0zvxYI+Pw5HWf877ff4KGHtKvlvvu8jiZ3pUtrZcbUVK37UoiyRtEttxVHof4o8krRjAL6oj8z/ShTRh8vpFGjRskDDzyQ6/MTJ06Up556SkREDh48KE2aNJHvv/9e5s+fL+XLl5cff/xRjh49Ki1atJDFixfL4cOHpWXLlvLLL7+IiMi0adOkR48eIqK1zO++++7M196zZ09mDfCXXnpJHnzwQRE5cYOLm266SRYvXiwiIlu2bJHzzz9fRET69u0rQ4YMERGR2bNnCyApKSnHxf/TTz9JjRo15JdffpFDhw7JhRdeKH369CnQ+XM7Lhxspah3HnxQxDmRr77yOpLAjB6tKSEE+7d4jmish56vgQP133BWqan6eBFrRhwTieVz586dy/r16zPv//bbb+zfv59FixZllrXt0KFDjqVsly5dSps2bahSpUrm+Tdu3Fig8xe2zK+JXuvWaUGsO+/UFno06NNHZ+M88ABceilk/Gn6XvR2uWzdWrDHA1CnTh1WrFiReX/s2LF89tlnHCskJhnlc1euXMnKlSv54YcfMje4yKt87rHj16xZwyeffJJ5XPbyuffeey9r1qxh4sSJJ5SvPSY9PZ0vv/wy8zW3b99OuXLlAAIqZZvbMYGeP9DjjD+I6JzzChUg21axEa1YMZ31cuSIVoGMla6X6E3ouRXKL0AB/eyioXxu9pK1x3YGatWqFYkZYwgfffRRZinbrJo3b86CBQvYvXs3aWlpvPPOOwU+f27HGX96911YsEDLmGcpAhoVzjoLnu2yjI8/hinFegRlnC3SRW9CHzpUC+ZnVcgC+sdEQ/nc0aNHk5SURP369alduzYTJkwAYNCgQSxatIjGjRvzySefUDOHf2zVqlVj8ODBtGzZkvbt2x+33V6g58/tOOM/aWnw2GNQp462cqNOYiJ93mnLhXxBP0awe8t+nXfp56SeW+d6qD+CUj73zTdF4uJ0tCYurkgDoiby2aBoeE2cqAOL77/vdSSFFBcnArKaulKcNLmTSXpBcXFeR1Yk5DEoGr0tdAh6AX1jjEpNhcGDdfu3q67yOppCyhhPq8daHuR5XuYuvuDCIo2zRbroTujGmJAYPRp27NDdgcJdFjdosnQ7DmIINdlCbyaQVuNMD4MKrYhL6BIrw9GmQOz3Inz27oVhw6BjR7j4Yq+jKYIs42xlSWU097GWeoy+eLrHgYVORCX00qVLs3v3bvvjNccREXbv3k3p0qW9DiUmvPAC/PprkeYXRIZu3XRj6oyNqq+OW8WVDbbz5Owm+HWBs/MqeSYkJEj2XXXS0tLYtm2bzW02JyhdujTVq1enZMmSXofia/v2af5r314X5vjNxo1Qty7ccgtMnux1NIXjnFsuIgk5PRdRK0VLlixpKw+N8dCYMZrU/VqC9txz4f77tYjX3XdDQo5pMXpFVJeLMcY7v/+u3S0dOkCjRl5HEzqPPw5VqmiRMb/17lpCN8YAMGEC7N6tCc/PKlTQEsBffqk7L/lJRPWhG2O8ceCA7sVZrx58+qnX0YReeroWGtu1C775Bk4+2euIApdXH7q10I0xvPoq7Nzp/9b5McWK6c5GW7fqnHu/sIRuTIw7elT7zps1g1atvI4mfNq21VWw//43ZBRUjXqW0I2JcR9+CJs3645EUbsqtJCGDdM9Up980utIgsMSujExbuRInXt+3XVeRxJ+F1ygBRgnTIBvv/U6mqKzhG5MDPvqK/j8c93Zp0RErUoJn8GDdVD0sce8jqToLKEbE8NGjtRpfHfc4XUk3jntNHj4YZg5E5Ys8Tqaogk4oTvnijvnvnbOzc7huZOcc9Odc5udc0udc/HBDNIYE3zJybq8v2dPyNjFMGb9859QtSo88kh0LzYqSAv9fmBDLs/dAewVkbOBF4BhRQ3MGBNa48bpIGjfvl5H4r1TToFBg2DxYpgzx+toCi+ghO6cqw50AF7O5ZCrgWMbTM4ALnGB7FhsjPHEgQNanOraa6FGDa+jiQx33glnnw39++tUzmgUaAv9P8AjQHouz58B/AggIkeAfUCl7Ac553o655Kcc0kpfpn4aUwUmjYN9uyBPn28jiRylCypJYPXroU33vA6msLJN6E75zoCv4jI8rwOy+GxE3qiRGSSiCSISEKVKlUKEKYxJlhEtKpinTrQurXX0USW66/XCoyDBsGhQ15HU3CBtNAvAjo555KBaUA759yb2Y7ZBtQAcM6VACoAe4IYpzEmSJYuhRUrtHVuHaPHc05Xjm7dqnPTo02+CV1EBohIdRGJB/4BzBORm7Md9gFwa8btLhnHRPFYsTH+NWYMlC+vmzyYE7VvD+3aaffL/v1eR1MwhZ6H7px70jnXKePuZKCSc24z8CDQPxjBGWOC65df4J134LbbdGaHOZFzWl43JUVr3EQTK59rTAx59lkYMAA2bIDzz/c6msh23XUwdy58/z1Urux1NH+y8rnGGNLT4aWXoE0bS+aBGDpUC3c984zXkQTOEroxMWLePG1t9uzpdSTR4YILoHt3GDsWfvzR62gCYwndmBgxaRJUqqSLiUxgBg/WaZ7RUl7XEroxMWDnTpg1C269FUqX9jqa6BEXB717645OGzd6HU3+LKEbEwOmTIEjR+Cuu7yOJPoMHKj/BP/1L68jyZ8ldGN8TkQHQ1u1ssHQwjjtNK3G+PbbuiArkllCN8bnFizQLeZsMLTw+vWDv/wl8jfBsIRujM9NngwVK0Lnzl5HEr0qVNBk/t//wvz5XkeTO0voxvjYr7/Cu+9C1642GFpUffpA9epaXjdSC5tYQjfGx6ZPh4MHoUcPryOJfqVLw5Ahug/re+95HU3ObOm/MT7WvLluZrFqlVVWDIYjR6BePb29Zo03G2vb0n9jYtC6ddqa7NHDknmwlCih5XW/+UangkYaS+jG+NSrr2oCujl7sWtTJNdcAy1a6CYYqaleR3M8S+jG+FBamm6jdtVVYJuDBZdzMHw4bN8Oo0Z5Hc3xLKEb40Nz5mjt89tv9zoSf/rb3/Sf5bPPwq5dXkfzJ0voxvjQlClQtSpccYXXkfjXs8/C779rmd1IYQndGJ/ZPeEdZs86TLedz1Pi7HhITPQ6JF+qXVvfAY0dCz/84HU0yhK6MX6SmMi0+/5HGqXozhTYskXX/FtSD4khQ3TgecAAryNRltCN8ZOBA5mSdhMNWEkDVutjqalaMtAE3emnw8MP6wKuL7/0OhpL6Mb4yoYtZVhGM7rz+vFPbN3qTUAx4JFHoFo1rciYnp7DAYmJEB8PxYrp5xC+W7KEboyPvF7+XopzhK5MPf6JmjW9CSgGlC2ri42WLtWW+nESE7XLa8sWLQAT4i4wW/pvjE8cPQrxp/1B/V8X83/pf//ziTJldP+5bt28C87n0tMhIUGnMH77LZx8csYT8fGaxLOLi4Pk5EKdy5b+GxMDFiyAbXvK0v2eUzRhOKefLZmHXLFi8Pzzupn0yJFZnsitqytEXWCW0I3xiTfegPLlodPwi7X1l56uny2Zh0WbNnDdddr9kpmvc+vqClEXmCV0Y3wgNVXrnl9/fZa3+ybsnn9eu8r79ct4YOhQ7fLKqkyZkK1GsoRujA988IGuWrRCXN6Ki9Odjd55Bz77DH13NGlS2LrA8h0Udc6VBhYBJwElgBkiMijbMbcBI4DtGQ+NEZGX83pdGxQ1Jng6dtSa51u2aH+u8c7Bg1Cnjm6IsXIllCwZ3Ncv6qDoIaCdiDQAGgJXOOda5HDcdBFpmPGRZzI3xgRPSgp8/LE2+iyZe690afjPf2D9ehg9OrznzvfHL+r3jLslMz4idEc9Y2LP9Ok6ZdG6WyJHx4768cQThZ6dWCgB/T93zhV3zq0EfgE+FZGlORzW2Tm32jk3wzlXI6hRGmNy9eab0KAB1K3rdSTmGOe0aFexYnDPPeHbVDqghC4iR0WkIVAdaOacy/6r8yEQLyL1gblAjpszOed6OueSnHNJKSkpRYnbGANs3KgrFK11Hnlq1tTJLB99lMMK0hApUI+biPwKLACuyPb4bhE5lHH3JaBJLl8/SUQSRCShim2jEjxhrBVhIsvUqdoavOkmryMxOenTB5o2hfvvhz17Qn++fBO6c66Kc65ixu2TgfbAN9mOqZblbidgQzCDNHkIc60IEzlE9Mfcti2ccYbX0ZicFC8OL70Eu3fDo4+G/nyBtNCrAfOdc6uBZWgf+mzn3JPOuU4Zx9znnFvnnFsF3AfcFppwzQkGDjxxp1orlxoTli2DzZttIWika9AAxo/XVnqoWXGuaFesWM4jLs7lUsvT+MX998OECbBzJ1Ss6HU0JlysOJefhblWhIkMR47oQFvHjpbMzZ8soUe7MNeKMJFh3jxtmVt3i8nKEnq0C3OtCBMZEhOhQgW48kqvIzGRpITXAZgg6NbNEngMOXAAZs6EG27QZebGHGMtdGOizIcfamVF+x9usrMWeoTYvh2WLIE//tC6HAC1a0OjRlCqlLexmcgydaruNt+6tdeRmEhjCd1Da9fCmDFaN3nz5pyPOekkaNJEl3bfeuuJ458mtuzdC3PmQN++umjFmKysy8UDX38NnTtDvXpaWKl2bd2H8KuvNLH/8IN+njFD/3BTU7XAT40a8Pjj4VlCbCLTu+9CWhp07ep1JCYS2cKiMDp8WBPyiBE6Q+G++3RxSKVKeX+dCHzxhW5v9d57ULUqTJwInTrl/XXGf9q10+65b77RSU0m9tjCogiweTNcdJEm8549tUbyk0/mn8xB/3AvvlhnNixfrgn96qt1UGzv3pCHbiLE9u2wYIG2zi2Zm5xYQg+DefN0cHPzZn3LPHFi4Vf3NWqkXTNDhui+hRdeqF00xv+mTdN3a1ZZ0eTGEnqIffQRdOigVW1XrYLrriv6a5YqpTuhzJ2rqwVbtNCa2Mbfpk6FhAQ491yvIzGRyhJ6CM2cqV0jtWvrW+Vgl1dp1Qq+/BJOOQXatNF/HsafvvkGVqywwVCTN0voIfLxx7qSLyFBpyUG0ldeGOedp/PXa9fW1v+CBaE5j/HWW29pv/mNN3odiYlkltBDYN06TeZ168J//xv6anhVquh5zjwTrrrKul/8RkS7W9q21QVFxuTGEnqQpaRoUi1bVpdolysXnvNWrqx96lWrwhVXwJo14TmvCb2kJB1Qt+4Wkx9L6EF0+LAuGNqxA95/XxcChVO1aprUy5TROeq7doX3/CY0pk7VgfDOnb2OxEQ6S+hB9MQTsHgxvPIKNGvmTQzx8TBrlv5TueEGXVVootfRozpdsUMH28jC5M8SepAsWADDh8Ndd3k/T7hZMy2JPn8+PPSQt7GYopk/H37+2bpbTGCsOFcQ7N0L3bvD2WfDCy94HY3q3l3nvT//vBb3uvVWryMyhTF1qo7DdOjgdSQmGlgLPQj69IGfftJdZMqW9TqaPw0bpvPT+/SBjRu9jsYU1MGDurL4uuvg5JO9jsZEA0voRTRzps4RHjwYmjb1OprjlSgBb7yhJXi7dtVBWxM95syB336zjSxM4CyhF8H+/VoxsUED6N/f62hyVr06vPyyFvV6/HGvozEFkZio01DbtvU6EhMtLKEXweDBWgFv/HhtDUeqa6+FXr200uNnn3kdjQnEr7/C7Nnwj39E9u+WiSyW0Atp9WoYNUpntbRs6XU0+Xv+eS0TcMcd+s7CRLYZM7SLzLpbTEFYQi+E9HTo3RtOPRWefdbraAJTpozOj9+6FQYM8Doak5/ERK2qmJDjNgbG5MwSeiEkJmqVwxEj4C8fJepqnmLF9HNiotfh5erCC7XPf+xYXQBlItO2bbBwobbObSMLUxD5JnTnXGnn3FfOuVXOuXXOuSE5HHOSc266c26zc26pcy4+FMFGgoMHdXCxSRPoXjxRtx/askUrKG3ZovcjOKkPHQq1asHtt+tepSbyvPWW/jrZYiJTUIG00A8B7USkAdAQuMI51yLbMXcAe0XkbOAFYFhww4wcY8Zot8Xw4VDsXwNPzIqpqTBwoDfBBaBsWZ31snmzboFnIs+bb0Lz5rpQzZiCyDehi/o9427JjI/sO0tfDUzJuD0DuMQ5/71Z3LNHW7h//7tu1svWrTkfmNvjEaJdO+jRA0aOhPXrvY7GZLV2rQ6433yz15GYaBRQH7pzrrhzbiXwC/CpiGSvuH0G8COAiBwB9gEnbOngnOvpnEtyziWlpKQULXIPPPMM7NuXZSA0ty2Igr01UQgMG6ZLyu+5R9/em8jw5ptQvLgWVjOmoAJK6CJyVEQaAtWBZs65utkOyak1fkKaEJFJIpIgIglVqlQpeLQe2roVRo/Wmij162c8OHSoTh/JqkwZfTzCVami/6AWLtR6IcZ76ek6/HL55XDaaV5HY6JRgWa5iMivwALgimxPbQNqADjnSgAVgD1BiC9i/Pvf2pIdknVIuFs3LWsYF6fTEeLi9H6UTB6+806tzPjQQ7qQxXhr4UKd4XLLLV5HYqJVILNcqjjnKmbcPhloD3yT7bAPgGP1/LoA80T880b+xx91Dvcdd+TQm9KtGyQna/MqOTlqkjnoW/vx43WXpUGDvI7GvPGGdoNdfbXXkZhoFUgLvRow3zm3GliG9qHPds496ZzrlHHMZKCSc24z8CAQoZVNCudYn7kfF+Q0bqyrXceOhQ0bvI4mdqWm6urQLl2ssqIpPOdVQzohIUGSkpI8OXdBbN+umy/fdhtMnOh1NKGRkgLnnKMlDD76yOtoYtO0aboxyrx5VozL5M05t1xEclxDbCtF8zFsmPam+LF1fkyVKrp93scfa8lWE35vvKF70LZu7XUkJppZQs/Djh06xnnrrbqq38/uvVdrh/zzn1Y3Pdx27oT//leHX4rZX6QpAvv1ycOoUbrJsp9b58eUKqUVGTduhHHjvI4mtkybpptB25LJ9cMAABMNSURBVGIiU1TWh56L/fv1LfBll8Hbb3sdTXiI6BzopCT47jutJmlCr3FjnfW6fLnXkZhoYH3ohfDyy7oq9OGHvY4kfJzTCpK//qrz7k3orV4NX3+tg+7GFJUl9BykpcELL0CrVpG3T2ioNWgA3bvrqtjkZK+j8b8pU6BkSZ3hYkxRWULPwdtv62KiWGqdZ/X00zo4F8FFI30hLU2X+nfsCJUrex2N8QNL6NmIwHPPwfnnw5VXeh2NN6pX19kuU6dqf7oJjf/+V2e43Hpr/scaEwhL6NnMnw8rV0K/frE9hax/f201PvKIVWMMlSlTdA1ArDYcTPDFcMrK2ejRmsiiqCRLSJQvD//6l/6D++QTr6Pxnz174IMPdFeikiW9jsb4hSX0LJKT4cMPtbZJ6dJeR+O9Xr10QVX//rpa1gTPtGm6gMtmt5hgsoSexbhxOnXv7ru9jiQynHSSDpCuXAnTp3sdjb9Mnqwziho29DoS4yeW0DOkpurc82uu0QVFRt10kyaexx+3kgDB8vXXsGKF1qM3JpgsoWeYOhX27oW+fb2OJLIUK6blg7//XuvamKKbPFnf/cT6OI0JPlv6j87iaNRI+4lXrdJuF/MnEd1YesMG2LwZTjnF64ii14EDcPrputG4bf1nCsOW/ufjiy80kffta8k8J87p/qM7d2rBMlN4s2ZpaYU77vA6EuNHltCBCRN0ml7Xrl5HErlatIBOnWD4cNi92+tootfkyVCrlm1iYUIj5hP67t269dctt0DZsl5HE9mGDtUqlMOGeR1JdPruO92R6PbbY3vRmgmdmP+1mjIFDh3SOdcmb3Xras3uF1/UrflMDhITdfJ+sWL6OTEx86lXXtGHbe65CZWYTugiOnOjZUuoV8/raKLDkCG6GcOTT3odSQRKTISePWHLFv3l2rJF7ycmkpam3S1XXqm1cowJhZhO6AsXwrffWuu8IGrV0u/X5MmwaZPX0USYgQN1QUNWqakwcCDvv6+Dyr17exOaiQ0xndAnToSKFeGGG7yOJLo8/rjOo37iCa8jiTBbt+b6+IQJULMmXHFFeEMysSVmE3pKCrz7rm7mcPLJXkcTXapWhQce0HokK1d6HU0EqVkzx4c3VmvNZ59p70vx4mGOycSUmE3or7+uGwz07Ol1JNHp4Yd1z1HbBCOLoUOhTJnjHytThkkNx1KihM5uMSaUYjKhi+iMgxYtoE4dr6OJThUrwqOPwpw5sHix19FEiG7ddJQ9Lk5XY8XFcXDMy7y6pDbXXAPVqnkdoPG7mEzoX30F69fbar2i6ttXk9SAAbYJRqZu3bQOc3o6JCfzTomb2LPHBkNNeOSb0J1zNZxz851zG5xz65xz9+dwTBvn3D7n3MqMj4geLps8Wd8Z22Bo0ZQpA4MGaemE2bO9jiYyjRkD555rK0NNeATSQj8CPCQiFwAtgD7Oudo5HLdYRBpmfETsLOU//tDBvBtu0OX+pmhuvx3OOUdb6UePeh1NZFm6VN8N9u1rK0NNeOT7ayYiO0RkRcbt/cAG4IxQBxYqM2bo8nUboAqOkiV1LHDdOnjzTa+jiSyjRmmjwTaBNuFSoHaDcy4eaAQszeHpls65Vc65j5xzOQ41Oud6OueSnHNJKSkpBQ42GF55RVuUF1/syel9qUsXSEjQeemHDnkdTWT46Sd45x1tOJQr53U0JlYEnNCdc6cA7wIPiMhv2Z5eAcSJSAPgReC9nF5DRCaJSIKIJFSpUqWwMRfapk2waJH+kVmZ3OBxTjfB2LpVt/EzWsHz6FG4916vIzGxJKCE7pwriSbzRBGZmf15EflNRH7PuD0HKOmcqxzUSINgyhTty+ze3etI/OeSS+DSS3UP0l9/9Toabx06pAm9Y0c46yyvozGxJJBZLg6YDGwQkedzOeavGcfhnGuW8boRVTU7PR3eeEOTzumnex2NPw0frtv4Pfus15F4a9o0XYl8331eR2JiTSAt9IuAW4B2WaYlXumc6+2cOza7tguw1jm3ChgN/EO82tsuFwsWaJeADVCFSGIiDa+J5xaZwn+GHWLrqFleR+QJERg5UhesXXKJ19GYWFMivwNE5HMgzx5nERkDjAlWUKEwZYrOOLjmGq8j8aFjZWNTU3mKfzGdG/nXQ6lMqZwYczshf/wxrFkDr71m4zQm/GJiduzvv2shrhtusEJcIZGlbGxNfuR+RvHG0ZtY+XBiPl/oP8OHa73zm27yOhITi2Iiob/7ri4osu6WEMlWNnYAz3Aqe+m348GYKgnw1VfatffPf0KpUl5HY2JRTCT0KVN0tsFFF3kdiU9lKxtbkX0MZjCf0T6mSgIMH65Fy+66y+tITKzyfULfsgXmz9epitanGSI5lI3tffLrnH/6Ph56CA4f9iiuMNq4EWbOhHvusYVExju+T+jH9ui95RZv4/C1HMrGlnxpHCNfqsCmTTB2rNcBht6IEdrN0rev15GYWOa8ml2YkJAgSUlJIT2HCNSuDVWq6ApRE14i8Pe/a5GqTZugcsQtNQuOH37Qioq9e8OLL3odjfE759xyEUnI6Tlft9BXrIBvvoGbb/Y6ktjknM7J3r/f3/uPPv20bi03YIDXkZhY5+uE/sYb+jb4+uu9jiR21amj/coTJug/WL/57jsddO/Vy1YgG+/5NqEfOQJvvaX1NE491etoYtuTT2q3V58+WoLBT556SksI9+/vdSTG+Dihf/op/PKLDYZGgooVdUrfkiW6gtIvNm3Sd4G9e9t+oSYy+Dahv/mmtsz//nevIzGg/1gvukg3lt6zx+togmPwYO3Se/RRryMxRvkyoe/fD7NmwY03wkkneR2NAS1bPHasJvPHHvM6mqJbtgymTtVVoX/9q9fRGKN8mdBnzYIDB2KuLlTEa9AA7r8fJk6M7mmkIvDQQ3DaadZ3biKLLxN6YiLEx9tS/0j01FNQq5Yujz940OtoCue992DxYhgyxDYaN5HFdwn9559h7lzo2tWW+keismV1UenGjTr7JdocPgyPPKIL1u680+tojDme7xL69Ok6Nc66WyJX+/Zw220682XlSq+jKZjx42HzZl3qXyLf3QSMCS/fLf1v1gzS0uDrr4P+0iaI9uzRVu5f/6qlAaJh8Hr7drjgArjwQvjoI3sHaLwRM0v/N23S2QfWOo98f/kLvPwyrFoFjz/udTSB6dtXF6yNG2fJ3EQmXyX0qVP1D812i4kOHTvqopyRI2HePK+jydt77+nsqUGD4MwzvY7GmJz5pstFBM47T7f/ivTkYP6UmgqNG+uOUqtXR2aZht9+0+6hSpUgKUmX+hvjlZjocklK0i4X626JLmXK6DTTn3/WWSORuGVd//7w0086O8eSuYlkvknoU6fqMuzOnb2OxBRUkybw7LO648/w4V5Hc7z33tOZLf/8JzRv7nU0xuTNF10uR49qV0vLlpoUTPQR0bGPd96Bjz+GSy/1OiL48Udd3XrmmfC//9nGzyYy+L7LZf58fcvetavXkZjCcg4mT9a+6n/8Q3cB8tLRo9p9l5amZZgtmZto4IuEPnWqbszboYPXkZiiKFtWZ5IcPQpXXeVtVcZBg3R5/7hxcM453sVhTEFEfUI/eBDefReuuw5OPtnraExRnX22dptt2qTTGlNTC/DFx4r4FCumn4/tEF5Ar70GQ4fC7bdbPX0TXaI+oc+Zo9PKbHaLf7Rrp7l4yRLdPjAtLYAvSkyEnj1hyxbtkN+yRe8XMKnPnauFw9q3123zjIkm+SZ051wN59x859wG59w659z9ORzjnHOjnXObnXOrnXONQxPuiaZOhapVoW3bcJ3RhEOXLppQ58zRTb4PH87nCwYOPLE5n5qqjwdozRqdJXXBBTBjhk1RNNEnkPJCR4CHRGSFc64csNw596mIrM9yzN+BczI+mgPjMz6H1L59MHu2btBrhZL8p2dPfff18MOwd692rZUrl8vBW7cW7PFsli2DK6+EU07RfyIVKhQuZmO8lG8LXUR2iMiKjNv7gQ3AGdkOuxp4XdQSoKJzLuS7LM6cCYcO2ewWP+vXD159VVf/tmun+8TmqGbNgj2exaef6ju8cuVg4UKdAmtMNCpQH7pzLh5oBCzN9tQZwI9Z7m/jxKSPc66ncy7JOZeUkpJSsEhz0KqVLkhp1qzIL2Ui2G236QKfdet0EdL8+TkcNHSoLjvNqkwZfTwXIjoA2qEDnHUWfPGFDsoaE60CTujOuVOAd4EHROS37E/n8CUnrFgSkUkikiAiCVWqVClYpDk46yzdoNcq3/lfx446jbBMGbjkEt1k4tChLAd066Zr8+Pi9BciLk7v5zJavnOn9pf36KE7Wy1cCNVC/p7SmNAKKKE750qiyTxRRHJai7kNqJHlfnXgp6KHZ8yfmjSBFSt0zGTECKhbV1vYmbNgunWD5GTd4SQ5OcdkfvCgLmCqW1f7ykeM0JktFSuG8UKMCZFAZrk4YDKwQUSez+WwD4DuGbNdWgD7RGRHEOM0BtDFR+PH6wYTp5yiLezzztPEvHp1zsW9RHRe++DB2nC/805dzr9ihfbRFy8e9sswJiTyreXinLsYWAysAdIzHn4MqAkgIhMykv4Y4AogFeghInkWagnVjkUmdojA//0fPP207noEOoW1Xj3tmilTRrtWVqzQGVGgXTcPPKADrNZVZ6JRXrVcfFGcy5gff4TPPtMZK99/DwcO6Ef58tpVk5AAbdrYoKeJfpbQjTHGJ3xfbdEYY0w0J/QgFWIyxhi/iM4F88cKMR2r3XGsEBNYlS5jTMyKzhZ6EAoxGWOM30RnQi9iISZjjPGj6EzoRSjEZIwxfhWdCb0QhZiMMcbvojOhF7AQkzHGxILonOUCmrwtgRtjTKbobKEbY4w5gSV0Y4zxCUvoxhjjE5bQjTHGJyyhG2OMT1hCN8YYn7CEbowxPmEJ3RhjfMKzHYuccynAliC8VGVgVxBeJ1rY9fpXLF0r2PUWVpyIVMnpCc8SerA455Jy247Jj+x6/SuWrhXsekPBulyMMcYnLKEbY4xP+CGhT/I6gDCz6/WvWLpWsOsNuqjvQzfGGKP80EI3xhiDJXRjjPGNqEnozrkrnHPfOuc2O+f65/D8Sc656RnPL3XOxYc/yuAJ4HofdM6td86tds595pyL8yLOYMnverMc18U5J865qJ3uFsi1OuduyPj5rnPOTQ13jMEUwO9yTefcfOfc1xm/z1d6EWcwOOdecc794pxbm8vzzjk3OuN7sdo51zioAYhIxH8AxYHvgDOBUsAqoHa2Y+4BJmTc/gcw3eu4Q3y9bYEyGbfv9vv1ZhxXDlgELAESvI47hD/bc4CvgVMz7p/mddwhvt5JwN0Zt2sDyV7HXYTrbQU0Btbm8vyVwEeAA1oAS4N5/mhpoTcDNovI9yJyGJgGXJ3tmKuBKRm3ZwCXOOdcGGMMpnyvV0Tmi0hqxt0lQPUwxxhMgfx8AZ4ChgMHwxlckAVyrXcBY0VkL4CI/BLmGIMpkOsVoHzG7QrAT2GML6hEZBGwJ49DrgZeF7UEqOicqxas80dLQj8D+DHL/W0Zj+V4jIgcAfYBlcISXfAFcr1Z3YH+149W+V6vc64RUENEZoczsBAI5Gd7LnCuc+4L59wS59wVYYsu+AK53sHAzc65bcAcoG94QvNEQf+2CyRaNonOqaWdfb5lIMdEi4CvxTl3M5AAtA5pRKGV5/U654oBLwC3hSugEArkZ1sC7XZpg77zWuycqysiv4Y4tlAI5HpvAl4TkZHOuZbAGxnXmx768MIupHkqWlro24AaWe5X58S3ZZnHOOdKoG/d8nrrE8kCuV6cc+2BgUAnETkUpthCIb/rLQfUBRY455LRvscPonRgNNDf5fdFJE1EfgC+RRN8NArkeu8A3gYQkS+B0mghKz8K6G+7sKIloS8DznHO1XLOlUIHPT/IdswHwK0Zt7sA8yRjFCIK5Xu9GV0QE9FkHs19rJDP9YrIPhGpLCLxIhKPjhl0EpEkb8ItkkB+l99DB71xzlVGu2C+D2uUwRPI9W4FLgFwzl2AJvSUsEYZPh8A3TNmu7QA9onIjqC9utejwgUYPb4S2IiOmA/MeOxJ9A8b9JfgHWAz8BVwptcxh/h65wI7gZUZHx94HXMorzfbsQuI0lkuAf5sHfA8sB5YA/zD65hDfL21gS/QGTArgcu8jrkI1/oWsANIQ1vjdwC9gd5ZfrZjM74Xa4L9e2xL/40xxieipcvFGGNMPiyhG2OMT1hCN8YYn7CEbowxPmEJ3RhjfMISujHG+IQldGOM8QlL6MZkyKjJfWnG7aedc6O9jsmYgoiW4lzGhMMg4Enn3GlAI6CTx/EYUyC2UtSYLJxzC4FTgDYist/reIwpCOtyMSaDc64eUA04ZMncRCNL6MYAGbvGJKI7yvzhnLvc45CMKTBL6CbmOefKADOBh0RkA7rV3WBPgzKmEKwP3RhjfMJa6MYY4xOW0I0xxicsoRtjjE9YQjfGGJ+whG6MMT5hCd0YY3zCEroxxvjE/wME0Isrw9gZrAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Ploting\n",
    "plt.plot(d,g(d),color='b',label=\"Generating function\")\n",
    "plt.scatter(x,y,color='r',label=\"Generated data\")\n",
    "plt.xlabel(r\"$x$\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The RBF kernel must now be defined, parameterised with a length scale parameter $l$.  The form of the RBF kernel is defined below:\n",
    "\n",
    "$$\n",
    "k(x,y) = \\exp\\left(-\\frac{|y-x|^2}{2l^2}\\right)\n",
    "$$\n",
    "\n",
    "For the implementation, the $\\texttt{scipy.spatial.distance-matrix}$ function is used to generate a matrix from the training point locations $x$, this is then scaled by $\\frac{-1}{2l^2}$ an fed through the $\\texttt{np.exp}$ function which performs the element wise exponential, finally giving $K_{xx}$ - the kernel matrix generated by the training point locations $x$.  This is an efficient vectorised procedure for evaluating the kernel matrix.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def RBFKernel(x_data, l):\n",
    "    \"\"\" Generate RBF kernel matrix efficiently using scipy's distance matrix function\"\"\"\n",
    "    x = x_data.reshape((len(x_data), 1))\n",
    "    Y = distance_matrix(x,x)\n",
    "    return np.exp(-pow(Y,2)/(2*pow(l,2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to solve the linear systems equation $(K\\alpha + \\lambda N I ) = y$, which  is done using $\\texttt{np.linalg.lstsq}$, a least sqaures approach for evaluating the inverse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve(K, x_data, y_data, reg):\n",
    "    return np.linalg.lstsq((K + (reg*len(x_data)*np.eye(len(x_data)))), y_data, rcond=None)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last function approximates the function using the dataset and the evaluation points, so that we can compare the approximating function to the original generating function over the space.\n",
    "\n",
    "Each evaluation point gives an approximation:\n",
    "\n",
    "$$\n",
    "\\hat{f}(d_j) = \\sum_{n=1}^{N} \\alpha_n k(d_j, x_n)\n",
    "$$\n",
    "\n",
    "or vectorised:\n",
    "\n",
    "$$\n",
    "\\hat{f}_{d} =   K_{dx} \\alpha\n",
    "$$\n",
    "\n",
    "where $K_{dx}$ is generated from the vectors $d$ and $x$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def approximate(x_new, l, alpha, x_data):\n",
    "    # get input data into (N,1) shape\n",
    "    x_data = x_data.reshape((len(x_data), 1))\n",
    "    x_new = x_new.reshape((len(x_new), 1))\n",
    "    \n",
    "    # Get kernel matrix\n",
    "    Y = distance_matrix(x_data, x_new)\n",
    "    Kxd = np.exp(-pow(Y.T, 2) / (2 * pow(l, 2)))\n",
    "    \n",
    "    # Perform approximation (equivalent to vectorised version of above markdown equation)\n",
    "    y_new = Kxd @ alpha\n",
    "\n",
    "    return y_new"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "###--- Parameters ---###\n",
    "lam = 0.01   # Regularisation constant\n",
    "l = 0.1     # Gaussian kernel length scale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = RBFKernel(x,l)                   # Generate RBF kernel matirx\n",
    "\n",
    "alpha = solve(K,x,y,lam)                  # Solve equations for the coefficients\n",
    "\n",
    "y_approx = approximate(d,l,alpha,x)       # Evaluate the approximating function given the coefficients"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEGCAYAAAB1iW6ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3hUxdfA8e+kQAgEpHeSIJ1AAkS6gHSlCAoWOqIIIohKEVSqKAqIiArSpIUiILGhP16kKwKhSu8BpPeQEEg57x8TYgi7ECDZu9nM53n2SXLv7L1nyeZwd+7MGSUiGIZhGOmfm9UBGIZhGKnDJHTDMAwXYRK6YRiGizAJ3TAMw0WYhG4YhuEiPKw6cZ48ecTPz8+q0xuGYaRLW7ZsuSAieW3tsyyh+/n5ERYWZtXpDcMw0iWlVLi9fabLxTAMw0WYhG4YhuEiTEI3DMNwEZb1odsSExPDyZMniY6OtjoUw0V5eXlRpEgRPD09rQ7FMFKdUyX0kydP4uPjg5+fH0opq8MxXIyIcPHiRU6ePIm/v7/V4RhGqnOqLpfo6Ghy585tkrmRJpRS5M6d23wCNFyWUyV0wCRzI02Z95fhypyqy8UwDNcnAocPw4ULcPUq3LgB/v5Qpgxkzmx1dOmb012hW83d3Z2goCDKly9PYGAgn3/+OfHx8fd9Xv/+/Slfvjz9+/dn8uTJzJ49+57tw8LC6NOnT6rEfP78eapVq0alSpVYt25dqhwTYObMmZw6dSrx51dffZU9e/ak2vENFxcSAn5+4OZGnG9xfn53Nd27Q9GiULIk1KgBTZtC69YQFARZs0K5cjBkCITbnTpj3JOIWPKoUqWKJLdnz567tjla1qxZE78/e/asNGjQQIYMGXLf5/n4+Eh0dHRahmbX/PnzpVOnTql+3Lp168rmzZtT/bhWc4b3mcubO1fE21viQZbQWsqyW0DEx+uWPPecyLffivz6q8j69SJhYSILFoh88IFIw4YiSulHkyYiGzda/UKcDxAmdvKqSejJJE3oIiKHDx+WXLlySXx8vMTGxkq/fv0kODhYKlSoIJMnTxYRkRYtWoibm5sEBgbKggULZOjQoTJmzBgR0UlxwIAB8sQTT0jJkiVl7dq1IiKyatUqadasmYiIDB06VLp27Sp169YVf39/mTBhQuL5R4wYIaVLl5aGDRvKSy+9lHjc27Zt2yZFixaVPHnySGBgoERFRd3xGhYtWiSdO3cWEZHOnTtL7969pUaNGuLv7y+LFi1KbPfpp59KQECAVKxYUQYOHCiLFi2SrFmzSqlSpRKPmzTBz5s3TwICAqR8+fIyYMCAO/79Bg8eLBUrVpRq1arJmTNnHun3kRac4X3m8nx9ZRfl5Ak2CoiUZq8spK3cLFbivk8NDxcZOlQkf36d2F9/XeTixbQPOb24V0J32j70vn1h+/bUPWZQEHzxxYM9p3jx4sTHx3Pu3Dl+/PFHcuTIwebNm7l58ya1atWicePG/PTTT2TLlo3tCQEPGzbsjmPExsayadMmli1bxvDhw1mxYsVd59m3bx+rVq0iIiKC0qVL07NnT3bs2MGSJUvYtm0bsbGxVK5cmSpVqiR7TUGMGDGCsLAwvvrqq/u+ntOnT7N+/Xr27dtHy5YtadOmDb/99huhoaFs3LgRb29vLl26RK5cufjqq68YO3YswcHBdxzj1KlTDBw4kC1btpAzZ04aN25MaGgorVq1IjIykurVqzNq1CgGDBjA1KlT+eCDDx7sH91I92aH16UHk/Ahghl0pSNz8CAOTtz/pnSxYjBsGLzzDgwdChMnwpIlMHMmNGuW5qGnaynuQ1dKuSultimlfrGxr4tS6rxSanvC49XUDdNakrDu6vLly5k9ezZBQUFUq1aNixcvcvDgwfs+/7nnngOgSpUqHDt2zGabZs2akTlzZvLkyUO+fPk4e/Ys69ev59lnnyVLliz4+PjQokWLR34trVq1ws3NjXLlynH27FkAVqxYQdeuXfH29gYgV65c9zzG5s2bqVevHnnz5sXDw4P27duzdu1aADJlykTz5s3v+3oN13TzJrz6KnRmFlXZxHaC6MpMncxBZ+sUyp4dxo+HLVugSBFo0QI++UTfVDVse5Ar9LeAvUB2O/sXisibjx6S9qBX0mnlyJEjuLu7ky9fPkSEiRMn0qRJkwc6RuaEW/fu7u7Exsbes03SdvKQ79ykQ/OSj7lOep7bxxeRBxrOd6+4PD09E491r9druJ7oaHj+eVi2DAa33MXw/2uJx42I/xp4e8OoUQ983MBA+PNP6NYNBg+GHTtgxgx9OONOKbpCV0oVAZoB09I2HOdy/vx5evTowZtvvolSiiZNmjBp0iRiYmIAOHDgAJGRkWl2/tq1a/Pzzz8THR3N9evX+fXXX1P0vPz587N3717i4+NZunTpfds3btyYGTNmEBUVBcClS5cA8PHxISIi4q721apVY82aNVy4cIG4uDjmz59P3bp1H+CVGa4mOlqPVlm2DL79Fkb9GIDH1Eng6wtK6a9TpkD79g91fG9vmDcPPv0Uvv8emjeHhLerkURKr9C/AAYAPvdo87xSqg5wAHhbRE4kb6CU6g50Byj2AB+9HOnGjRsEBQURExODh4cHHTt25J133gH0sL1jx45RuXJlRIS8efMSGhqaZrE88cQTtGzZksDAQHx9fQkODiZHjhz3fd7o0aNp3rw5RYsWJSAggOvXr9+zfdOmTdm+fTvBwcFkypSJZ555ho8//pguXbrQo0cPsmTJwoYNGxLbFyxYkE8++YSnnnoKEeGZZ57h2WeffeTXa6RP0dHQqhX8738wdarucgF08n7IBG6LUjBgABQqBJ0766T+yy/mSv0O9u6W3n4AzYFvEr6vB/xio01uIHPC9z2Alfc7rrOOcnE2ERERIiISGRkpVapUkS1btlgcUfpn3mepJz5epGNHPV5u+nTHnXfOHBE3N5GnnhKJjHTceZ0BjzjKpRbQUin1DOAFZFdKzRWRDkn+U7iYpP1U4NNH/Y/G0Lp3786ePXuIjo6mc+fOVK5c2eqQDCPR2LEwZw6MGAGvvOK483ZIyD6dO+sPAUuWgJuZJnn/hC4ig4BBAEqpekC/pMk8YXtBETmd8GNL9M1TIxXMmzfP6hAMw6Zly2DgQGjbFqwYmdqhA1y6BG+9Be+9B5995vgYnM1Dj0NXSo1AX/r/BPRRSrUEYoFLQJfUCc8wDGd04AC8/LKe2/Hdd7p/2wq9e+tYxoyBUqWS9N9nUA+U0EVkNbA64fshSbYnXsUbhuHaYmL01bG7O4SG6hosVlFKD3E+fBh69oQSJaBePevisZrpdTIM44F8/DFs3qyHJzrDYDUPD1i4UCfzl1+GhPlyGZJJ6IZhpNjGjTBypL5Cb9vW6mj+kz07LFqky/F26ABxcVZHZA2T0G1YunQpSin27dtnaRxDhgyxWfflQV25coVvvvkm8edTp07Rpk2bRz4uwM2bN2nYsCFBQUEsXLgwVY4JEBoaekep3tT6tzAeXmQkdOwIhQtDCsoGOVxAgI5rxQr9KSJDsjeeMa0fzjwOvW3btlK7dm0ZOnRoqh0zJiYm1Y71oI4ePSrly5dPk2Nv2LBB6tSpk+rH7dy58x3VIFOTs7zP0psBA/R485UrrY7Evtvj4t3cRFavtjqatIEpn5tyERERUqhQIdm/f7+ULl06cfuqVavkySeflFatWknZsmXl9ddfl7i4OBHRJWPfeecdqVSpktSvX1/OnTsnIrp07qBBg6ROnToyduxYOXbsmNSvX18qVKgg9evXl/DwcBERadmypcyaNUtERCZPnizt2rUTkTuTmq+vrwwaNEiqV6+eOMGocePGUrx4cZk0aVJi7PXr15dKlSpJQECAhIaGiojIiy++KF5eXhIYGCj9+vW7I8F/99130rp1a2nSpImUKFFC+vfvn/iap02bJiVLlpS6devKq6++Kr169brj3+rs2bPy+OOPS/bs2SUwMFAOHTokvr6+cv78eRER2bx5s9StW1dE7l0ieNasWVKhQgWpWLGidOjQQf7880/JmTOn+Pn5JR436b/FihUrJCgoSAICAqRr166Jdeh9fX1lyJAhia9/7969Nn/HzvA+S2927RLx8BDp1s3qSO4vIkKkZEkRX1+Ra9esjib13SuhO235XLb0hcupXD83ZxBUuXfVr9DQUJo2bUqpUqXIlSsXW7duTZzMs2nTJvbs2YOvry9Nmzblhx9+oE2bNkRGRlK5cmXGjRvHiBEjGD58eGIp2ytXrrBmzRoAWrRoQadOnejcuTMzZsygT58+hIaGMmXKFGrVqoW/vz/jxo3j77//thlb0aJF2bBhA2+//TZdunThzz//JDo6mvLly9OjRw+8vLxYunQp2bNn58KFC1SvXp2WLVsyevRodu3alVjeN3kFxO3bt7Nt2zYyZ85M6dKl6d27N+7u7owcOZKtW7fi4+ND/fr1CQwMvON5+fLlY9q0aYwdO5ZffrmrCOddbJUIPnDgAKNGjeLPP/8kT548iaV7W7ZsSfPmze/qGoqOjqZLly788ccflCpVik6dOjFp0iT69u0LQJ48edi6dSvffPMNY8eOZdq0DFV+KE2IwBtv6H7q0aMf4UA3L8HVPXBtL1w/Am6e4JEVPLPDY4GQqzK4ez1yvNmy6VK7Tz4J776rS8hkFM6b0C0yf/78xOTw0ksvMX/+/MSEXrVqVYoXLw7Ayy+/zPr162nTpg1ubm68+OKLAHTo0CGxXC6QuB1gw4YN/PDDDwB07NiRAQMGALqY1ogRI3jqqadYunSp3fK1LVu2BKBChQpcv34dHx8ffHx88PLy4sqVK2TNmpXBgwezdu1a3Nzc+PfffxNL5N5LgwYNEmvElCtXjvDwcC5cuEDdunUTY2nbti0HDhxI4b+ibbdLBGfOnDmxRPDKlStp06YNefLkAe5funf//v34+/tTqlQpADp37szXX3+d+DtLWqr49r+18WhCQmDtWp0YE35NKRcXDSdC4cgMOLMCSKjUqTxAklXidPOEnFWg2PPg1wGyFHjomGvWhP79dTGv1q3h6acf+lDpivMm9PtcSaeFixcvsnLlSnbt2oVSiri4OJRSfJYwBS15iVl7JWeTbs96j0G6Sdv9888/5M6d+441PJO7XfrWzc3tjjK4bm5uxMbGEhISwvnz59myZQuenp74+fndVT73XseFRy/d6+HhkbgG671K9yY9T2qV7k16DlO6N3VcuaKvcqtV0+VrUyw+Dg59Czs/hFuXIKsvBHwAuatDjrL6ZxTE3dD7L4bBhQ1wdhVs6w/b34OCTaDMO5C//kPNXBo+XBfv6tYNdu2C+1wruAQzyiWJxYsX06lTJ8LDwzl27BgnTpzA39+f9evXA7rL5ejRo8THx7Nw4UJq164NQHx8PIsXLwb0VP3b25OrWbMmCxYsACAkJCSx3aZNm/jtt9/Ytm0bY8eO5ejRow8V/9WrV8mXLx+enp6sWrWK8ISVdu2Vwb2XqlWrsmbNGi5fvkxsbCxLlixJ0fP8/PzYsmULQIqe06BBA77//nsuXtTlgO5XurdMmTIcO3aMQ4cOATBnzhxTujcNjRoF58/DpEkPUCvl8nb4v5oQ1gtyVoL6/wctj0DFEVD4GcjmD8pNJ2kPb/AuAkVbQaVPoekmaL4Pyg6Ay9tgZUNY8SScXv7AK1tkzgyzX1rGudOx9Mv9nV6wOiTkgf8N0hOT0JOYP38+rVu3vmPb888/n1hPpUaNGrz33nsEBATg7++f2DZr1qzs3r2bKlWqsHLlSoYMGXLXsQG+/PJLvvvuOypWrMicOXOYMGECN2/e5LXXXmPGjBkUKlSIcePG8corrzzUFXL79u0JCwsjODiYkJAQypQpA0Du3LmpVasWAQEB9O/fP0XHKly4MIMHD6ZatWo0bNiQcuXKpah079ChQ3nrrbd48skncXd3v2/78uXL8/7771O3bl0CAwMTSxW/9NJLjBkzhkqVKnH48OHE9l5eXnz33Xe0bduWChUq4ObmRo8ePVL0mowHc/y4Xv6tUyeoVCmFTzrwDfweDJHHoGaITuYFGuoEnlLZS0PQx/o/geCvITIcVjWBVY3h6gOUiQoJofInbenHWL6jK6vD/aB7d9dO6vbulqb1w1lHudiTdFHn5JIvLO0qbpfujYmJkebNm8sPP/xgcUSpw5nfZ86kUyeRzJn1os33FR8nsuVdkRBEVrcQuXkp9QKJjRbZN0Hk+8dE5nmIhL0tcuvq/Z/n6ysCEkkW8eewlGavRJNJb0/HuMcoF3OFbtg1bNgwgoKCEj+RtGrVyuqQDAfZsUOXxe3TJwXT++OiYf0LsG8clHoTnlwKmXKmXjDumaF0H2hxAIp3hf1fwK8BcOr3ez/v+HEAvLnBJHqynzJ8wqDE7a5IyUPe/HpUwcHBEhYWdse2vXv3UrZsWUviMTIO8z67v6ef1tP8Dx+GnPfKzfFx8OcLcOIHqPw5lO6b9qUXL/wNf3eFa/vAv7MeQJHpsbvb+flBwn0kgHaEsITn2VHoGcr8+0faxpiGlFJbRCTY1j5zhW4Yxh1WrYLff4f3379PMheBLX0Skvl4KPO2Y+ro5qkOT2+D8oPh2FxYVhHOrrm73ahRd6xPN5638SaK3rnmPuj91XTDJHTDMO4wbJhet7NXr/s03P0xHPxGj0gp09cRof3H3QsCR0GjP8EtM/zxlB7qGHfrvzbt2+vB8wkLVef3zcKITodYsasgKVg7PV0yCd0wjERr1uhJRAMHgte9Jm2Gfw87PwC/jhD0icPiu0ueavpq/fFusOdTPcTx+rH/9rdvD8eOQXw8HDtGz+lPUKECvP02REVZFXTaMQndMIxEw4dDgQLw2mv3aBRxGDa+CnlqQLVpDzYkMS14ZoNqU6H2It2v/lslOGH7EtzDQw/FPH5czyJ1NSahJ3P27FnatWtH8eLFqVKlCjVq1GCphZ/PVq9ezV9//fXAz/Pz8+PChQv3bDNz5kzefPPNNDm/kf6sW6f7zwcOhCxZ7DSKuwl/vghuHlBrPrhncmiM91Ssjb5a9ykB656DLe9AfMxdzerWhZde0gn9yBEL4kxDJqEnISK0atWKOnXqcOTIEbZs2cKCBQs4efJkmp73XlPUrU6oVp/fcJyRIyF/fj33xq7tA+HSFqj+XcL0fSeTrTg0Wq+HT+4fDysbwY276xmNGaOX0Esop+Q67A1QT+tHqkwsmjtXTxJQSn+dO/fBnp/MihUr7lnbOzY2Vvr16yfBwcFSoUIFmTx5sojoSUd169aV559/XkqXLi3t2rWT+Ph4EREJCwuTOnXqSOXKlaVx48Zy6tQpEbm7tO5PP/0kVatWlaCgIGnQoIGcOXNGjh49Kvnz55dChQpJYGCgrF27Vs6dOyfPPfecBAcHS3BwsKxfv15ERC5cuCCNGjWSoKAg6d69uxQrViyxjG1SM2bMkJIlS0qdOnXuKImb0vPbapfemIlFd9uwQRfTHjv2Ho3+XaYnDoW95bC4HsmROSILsoj8UEjk/N937R4+XL/mNWssiO0R4JL10OfOFfH21i/h9sPb+5GS+oQJE6Rv375293/77bcycuRIERGJjo6WKlWqyJEjR2TVqlWSPXt2OXHihMTFxUn16tVl3bp1cuvWLalRo0ZiffQFCxZI165dRUQn9J49eyYe+9KlS4n/CUydOlXeeecdEdF1xMeMGZPY7uWXX5Z169aJiEh4eLiUKVNGRER69+4tw4cPFxGRX375RYC7EvqpU6ekaNGicu7cObl586bUrFkzMaGn9Pz22qUnJqHfrVUrkVy5RK5ft9PgVoTI0mIiP5fVMzfTi0vbRUL9ReZnFjlyZ26IjBQpUkSkShWRhKUN0oV7JXTnrbZ4P++/f/dt6qgovb19+1Q5Ra9evVi/fj2ZMmVi8+bNLF++nJ07dyYW4rp69SoHDx4kU6ZMVK1alSJFigAQFBTEsWPHeOyxx9i1axeNGjUCIC4ujoIFCyYeP2lp3ZMnT/Liiy9y+vRpbt26hb+/v82YVqxYccfSbNeuXSMiIoK1a9cmlott1qwZOW0MIN64cSP16tUjb968iee/XRI3pedPaTsj/ThwAH78ET74AOwWB935AUSd0N0Z7pntNHJCOQOhySZY3wY2dICru/RwR+WGtzd88oleVm/uXF2zJr1Lv33o9qbvPsK03vLly7N169bEn7/++mv++OMPzp8/D+hPMxMnTmT79u1s376do0eP0rhxY8B+adjy5csntv/nn39Yvnx5YrukpXV79+7Nm2++yT///MO3335rt+xtfHw8GzZsSDzmv//+i4+PD2C/nG9S9tqk9PwpbWekH+PGQaZMYPf++IWNsP9LKNkT8tZ0aGypwisPPLUcSnSHPaN1mYI4/b5t1w6eeAIGDdJrpqZ36Teh2yswcd/CE/bVr1+f6OhoJk2alLgtKsmngCZNmjBp0iRiYvSd8wMHDhB5j3dB6dKlOX/+PBs2bAAgJiaG3bt322x79epVChcuDMCsWbMStycvI9u4cePE1ZCAxFWI6tSpQ0hCFbnffvuNy5cv33WOatWqsXr1ai5evEhMTAyLFi164PPba2ekT2fPwqxZ0Lkz5Mtno0F8DGx6DbIUsna8+aNyzwRPTNblCU4sgZWN4eYl3Nxg/Hg4dQrGjrU6yEeXfhN6smm9gP551KiHPqRSitDQUNasWYO/vz9Vq1alc+fOfJowYPXVV1+lXLlyVK5cmYCAAF5//fV7jlDJlCkTixcvZuDAgQQGBhIUFGR3xMiwYcNo27YtTz75ZOLqPaCXrVu6dClBQUGsW7eOL7/8krCwMCpWrEi5cuWYPHkyoMvWrl27lsqVK7N8+XKK2fiPrWDBggwbNowaNWrQsGHDxJWYHuT89toZ6dNXX8GtW3oRC5v2T4Qr/8AT3+il4tIzpXR5gloL4eJG+L/aEHmcWrWgTRs98uXMGauDfET2OteTPwB3YBvwi419mYGFwCFgI+B3v+M54ygXI2MwN0W169f1jdBWrew0uHFe5PscIiubOjQuhzizWr+2UF+RiMNy4IBeBLtHD6sDuz9SqXzuW4C96vLdgMsiUgIYDzhmDlayab2pdTPUMDKCWbPg0iW99qZN/wyD2OtQeZwjw3KM/HWhwUqIiYAVdSmZ/wA9e8LUqbBvn9XBPbwUJXSlVBGgGWBvCfVngdsdqouBBupBFoo0DMOh4uP1FPgnntALKt/l6h44NBlKvA45yjk8PofIVRkarNKzX1fUZdjbe8maFd57z+rAHl5Kr9C/AAYA8Xb2FwZOAIhILHAVyJ28kVKqu1IqTCkVdnvkSHLiqnUtDadg3l/aihX6SrRPHzsNtvYDj2xQYbhD43K4nBWh4RpAyLW9ER+/H86PP+oyCOnRfRO6Uqo5cE5EttyrmY1td/3liMgUEQkWkeDbY6GT8vLy4uLFi+aPzkgTIsLFixfxumcZwYxh4kQ9qqVtWxs7z6yA079BwId6yJ+ry1FWD2uMvU7PMo2pUOo8AwY88JrUTiElE4tqAS2VUs8AXkB2pdRcEemQpM1JoChwUinlAeQALj1oMEWKFOHkyZPYu3o3jEfl5eWVOAEsozp8GH79FT78EDInnyMkAjveB+9iuh5KRpGzItT9BbdVjVj5YVOKd19FaGh2kq0Z7/Tum9BFZBAwCEApVQ/olyyZA/wEdAY2AG2AlfIQl9menp5m5qFhpLGvv9aFqV5/3cbOf3+Bi5t0Wdz0NCM0NeSrDbUXk3ttK35+72V6Dv6JFi3c8UhH8+kfehy6UmqEUqplwo/TgdxKqUPAO0A6vq1gGK7r+nWYPl13tRQqlGynxMPODyFbCfB3gXnwD6NwM1Twl9QtsYwOFT9kxgyrA3owD/R/j4isBlYnfD8kyfZowFZvnGEYTmTuXLh2DXr3trHzxBK4sgNqzAU3T4fH5jRK9EAubWPws5/w+pwgItu/YL/GjZNRVt2ADA4OlrCwMEvObRgZkQhUqqQnTG7dmmw95/g4WBagVx96eie4uVsWp1OIu8W10Kdwv7qdeRc28NqAilZHlEgptUVEgm3tS79T/w3DeCAbN8KOHdCjR7JkDhC+QC/fVmG4SeYA7pnI/swSbkoOnnR/iYvnblgdUYqYhG4YGcTkyZAtm64weAeJhz0fQ44AKPqcJbE5pSwFuFZ+FmUK7mXPvIFWR5MiJqEbRgZw+TIsXAgdOkBCteX/nAzVM0PLD7Z+wWcn41e9Ef93/C2ezDeRs9t/tzqc+zK/PcPIAGbPhuhoG0MVRWDXKD2ypdgLlsTm7Mq+9Am7T5Yn87auEH3vhdetZhK6Ybg4Ed3dUq0aBAUl23n6f3B5K5R/z/Sd21HELwv/ux6Ct/tFLq/sa3U492QSumG4uLVrdd2WHj1s7Nw9CryLgF9Hh8eVnnTpG8iE5QPJeSUEzq62Ohy7TEI3DBc3bRrkyAEvJO9RObcOzq+Hsv31ij6GXblygVvFQRw950fkml56JScnZBK6YbiwK1dg8WI9siX5Al/s+Qwy54HHX7UktvTmjd7ejFj2JVlj9yD7vrA6HJtMQjcMFzZvnr4Z2q1bsh1X98KpX6BkL/BInukNW7JkgTovt+CnLS2I2z4cIk9YHdJdTEI3DBc2fbq+EZpk+Vht3zhw94JSvSyJK73q1Am+/nsCsTFxxG8bZHU4dzEJ3TBc1LZteop/t27JZobeOANH50DxruB197oEhn3u7tB7kD8Tfu+DCp8HV3ZZHdIdTEI3DBc1fbqud37XUrsHJuqbeqXftiSu9K5ZM/jz0gAion2I2fKh1eHcwSR0w3BBN27oyorPPw85cybZEXMdDk6Coq0he0nL4kvPlIIho3Iz5pd+eJ4NhQubrA4pkUnohuGCfvgBrl61cTP0yHdw6zKU6WdJXK4iOBhOZ+/L+Yg8RG963+pwEpmEbhguaOZM8PODevWSbJR42D8BcleHvDWsCcyFfDjCh89+GYzXlRVwZqXV4QAmoRuGyzk+YSl/rIin87HhuBX3g5AQvePfX+D6YShj+s5Tg68vZA7oyb+XCnHt71FWhwOYhG4YriUkhDn9dyC40ZmZEB4O3Vy77+gAACAASURBVLvrpL7/C/AuakrkpqL+73kxdV1fsketRC5utTock9ANw5XI4PeZGdOeeqzCn2N6Y1QUfNUfzq6CUr3BLR2teuzkcuSAovW6c+2GDyf+b6zV4ZiEbhiu5K/jRThESbow884dgafB3RtKmGn+qa3zazn4YWd3CsV8T/Sl8LsbhIToGxpubvrr7S6wNGASumG4kJnZepGV6zzPkv825gBqoicSZcpp76nGQ/LwgBLN3kJEsfP7CXfuDAnRXV7h4bqOcdIusDRgErphuIioKFgY24a27kvJRuR/O5p6ggdQuo9lsbm62o2L8ufJlyibaSpnjl/5b8f77+tfTFJRUXp7GjAJ3TBcxNKlEBHtSZf3CughGEpB8WLQwhsKPQPZS1kdoksr3rwfPl7XWT19yn8bjx+33dje9kdkErphuIjZs3UX7ZMjGsGxYxAfDys+BnUVSpmr87RWrGIgRyLr8cRj37JubXzCxmJ2GtvZ/ohMQjcMF3DqFKxYoReBdkv6V33gS8heGgo2siy2jKRQne48nv8Is8esJCYGGDXq7kL03t56exowCd0wXMC8efqCvGPSleQubISLm/RQRWX+1B3Bq0RrbpKbRv5T+OordGW0KVP+6wLz9dU/31UxLXXc97eslPJSSm1SSu1QSu1WSg230aaLUuq8Ump7wsOMjTIMBxGBWbOgenUolbSbfP+X4Jkd/DtZFluG4+5FptKdea7qUr4ae5ZTp9DJ+3YX2LFjaZbMIWVX6DeB+iISCAQBTZVS1W20WygiQQmPaakapWEYdu3YAbt26cUXEkWdguPfQ/FXwNPHstgyIlXiNTzcYnnxiVn07evYc983oYt2PeFHz4SHpGlUhmGk2Jw54OmZbBHoQ9+CxJkViayQowzkq0O/1lNZvDieH3903KlT1LGmlHJXSm0HzgH/JyIbbTR7Xim1Uym1WClV1M5xuiulwpRSYefPn3+EsA3DAIiN1XNUmjeH3LkTNsbd0gm90DPgU8LS+DKsx7uTy/MQ3Zqt5o03dCljR0hRQheROBEJAooAVZVSAcma/Az4iUhFYAUwy85xpohIsIgE581rlr5KNQ6cWmw4l//7Pzh7Nll3y/FFEH1W3ww1rFHseciUk9Hdp3PmDAwY4JjTPtCtbxG5AqwGmibbflFEbib8OBWokirRGffn4KnFhnOZM0evSPTMM0k2HpgIPqXMUEUruXtBsbbkjgplYL9IpkyB1avT/rQpGeWSVyn1WML3WYCGwL5kbQom+bElsDc1gzTuwcFTiw3nEREBoaHw4ouQKVPCxoub4eJGKPWmGapoNb/2EBfF0Fd+5PHHYeHCtD9lSupoFgRmKaXc0f8BfC8ivyilRgBhIvIT0Ecp1RKIBS4BXdIqYCMZB08tNpzH0qV67dA7xp7vnwge2aB4Z8viMhLkrQ3eRcl8eh7r17cjf/60P+V9E7qI7AQq2dg+JMn3g4BBqRuakSLFiuluFlvbDZc2Zw74+0ON26vJ3TgLxxdCie56/LlhLeUGvi/Dvs8pUP0CqDxpfkrzmSy9c/DUYsM5nDoFf/yhp/orlbDx0BSIvwUlzVBFp+HXDiQWTixyyOlMQk/vHDy12HAO8+fre+AdOiRsiI+BQ5OhQGM9DtpwDo9VhBzl4ZhjBimYhO4KHDi12HAOc+dC1apJpvqf+AFunILSZqiiU1FKX6Wf/xOuH0vz05mEbhjpzK5dsH17spuhByZCtuJQ8GnL4jLs8H1Zfw2fn+anMgndMNKZuXPB3T3JVP9L2/QVYMle4OZuaWyGDdn8wfclh9yoNst/W2z7dli5EtauhY0bITIS4uL0vrJldQW96tWhRQu9wriRscXH61K5TZtCvnwJGw9M1AtAP/6KpbEZ91Ar7a/OwSR0S4jA77/Dxx/D+vV6W4kS0KQJ5Mqlr75iY3UVvVmz4OuvIVs26NoV+vTRbY2Mad06OHECPvssYUP0BTg2L2EB6Mcsjc2wnknoDrZtG7z6KmzdCkWLwoQJ0KYNFCpku31cHISF6aQ+eTJ89RX06AGjR0N2M9Q4w5k7V//n3rJlwobD0yD+pp4ZamR4SsSaSrjBwcESFhZmybmtEB8Pn38OgwdD3rzw0Ud6MErilO3kbl3R07gvbYbLO8EzGxFxhfl1VVH6ff4MyrsQ336brIaH4dKio6FAAXj2Wf3JjfgY+Kk4+JSGBiusDs9wEKXUFhEJtrXPXKE7wKVL+gbWH39A69YwdWqSUqfJRZ6A3aPgyAz9BwuQ1Q/iovGJPstL/sILEz1Yvuc5Pn63F8uXP8m4cQp3cy/M5f36qy7Dmjj2/GQoRJ2E4G8sjctwHiahp7EzZ6BRIzhwQCfybt2SzOxLKuY67Bis61gjULybLsGZqwpkyqnbxMdAxEHcDs+gied0mpb7nh82t6Zdm+lMm50TH7MwjUsLCdFX6PXrJ2zYP0EPVSxkPqYZmknoaSg8HBo2hNOnYdkyaNDATsNrB2Hdc3Btj07kAe9DVt+727l5Qo5yUHksquII2P8lz8qHVPEP4s0X5zFqSi2KFEnTl2RY5NIlfYXeq5e+ac6lLXqoYuXxZqiikciMQ08j4eFQuzZcuKAXIbCbzP/9Bf73hJ7lV+93qDbFdjJPzsMbyr+He5M/yZvPg+kv12VSP11M33A9ixfDrVtJulv2f5lQVbGrpXEZzsUk9DRw9apeEiwiAlatSlINL7mjIbCmpf7Y3HTLwy1IkKcq3s9tJSJbI0Y2f40v+87mwoVHCt9wQnPnQpkyUKkSuqpi+AIo3gUymckJxn9MQk9lsbF6wYF9+2DJEggKstPw5I/wd2fIXw8arYdsfg9/0kw5yNniB6541Wdk866MeXMhly8//OEM5xIersefd+yYcP/l0Le6qqIZqmgkYxJ6KhLRE3/+9z+YNOke3SxnVsD6FyBXMNT5UXefPCqPLORq9SNXM9ViVLP2jHl7GTExj35Yw3rz5umv7doBcTfh4Df6Rmj20pbGZTgfk9BT0axZOpEPGKAnD9l0eQeseVb/MdZbBp6pODTFIyu5Wv/KVbcK9KvZgZHvHUu9YxuWENELWTz5pF7/m/AFegHo0n2tDs1wQiahp5KDB+HNN+Gpp/SUfptuXYF1z+sp2k8th8y5Uj8QTx9yt1xMlizxtHisLd9Ounn/5xhOa/t22Ls3oSKyCOz/QtfXLtDQ6tAMJ2QSeiq4dUt/HM6cGWbPxvYkHxH4uytEhkPt7yFLgbQLyOdxMtWZyROPhxG76V2HrDZupI25c8HTE9q2Bc6tgcvb9dW5zckMRkZnEnoqGDJE11uZNg3748D3jdMz+yp9BnlrpXlM7r6tuFn8XXo1+pqQ0Ys5fz7NT2mksrg43X/erJku2sb+LyBzHr2avGHYYBL6I1q3Tle+e/11Pa3fpgubYPt7ULSNQ/s+M1f9hKgsTzCq1Ru8+dpF4uMddmojFaxcqWcad+gARByCkz9BiR7gkcXq0AwnZRL6I7h1SyfyYsVg3Dg7jeJuwsaukKUgVJvm2I/Kbp5415tKbp/LNM3fjy++cNypjUc3Z46ugd+sGXqav5sHlOxpdViGEzMJ/RGMGaNvWH3zDWTNaqfRrpFwdQ9UnWLNJJCcgbiVH0DXujNZPnsFGajAZboWGQk//KD7zr3UJTg8A3zbgbedOsuGgUnoD+3QIRg5Uv/B2S1he2kr7BkN/p2hkHVrPaoKHxKXtRTfduvO692iiI62LBQjhUJDdVLv2BE9kSguCsq8Y3VYhpMzCf0hiEDPnnpUi91ujLhbelSLVz6oMt6h8d3F3Qv36lPwzX2U58t8xPDh1oZj3N+cOeDrC7Vr3NR1Wwo0hpwVrQ7LcHImoT+EJUtgxQo93rzQqhA948PNTX8NCdGNDkyEKzvhiUn/lb+1Uv664NeB/s0/Z8H0cDZtsjogw57Tp3VBtw4dwO3EfIg+A2XftTosIx24b0JXSnkppTYppXYopXYrpe66vlNKZVZKLVRKHVJKbVRK+aVFsM7g1i0YOBACAqBH9nnQvbsutiGiv3bvDiFfwz/DoVAzKPKs1SH/J/BjPDwV4zsPoksXTNeLk5o/X69w1aG9wN5x8FgFKPAQhduMDCclV+g3gfoiEggEAU2VUtWTtekGXBaREsB44NPUDdN5fP01HDkCY8eC+4eDISrqzgZRUbB1IMRH61rVziRrUVTZfrSqNJ/sMX8zcqTVARm2zJkDwcFQJsdyuLoLyrxrJhIZKXLfhC7a9YQfPRMeyRcifRaYlfD9YqCBUq73Drx0Sd8IbdJEPzh+/O5GxYEqkVD6bche0tEh3l+5geBVgDl932HMGGHPHqsDMpLatUtP9+/YEdj7mR7u6vuy1WEZ6USK+tCVUu5Kqe3AOeD/RGRjsiaFgRMAIhILXAXuWjVTKdVdKRWmlAo7nw6nLo4apWudjxmTsKFYsTsbKKAzcM0dAj5wcHQp5JkNAj+iZM4NdHhyET176t4iwznMmaNLR3R4JgzOrtQXBu72VhI3jDulKKGLSJyIBAFFgKpKqYBkTWxdjd+VJkRkiogEi0hw3rx5HzxaCx09ChMnQteuUKFCwsZRo8A7SenbGkAJwOfV1K2imNr8u0CO8nzedQh/ro9l5kyrAzJAT/WfOxeefhpynf4UPHNAydetDstIRx5olIuIXAFWA02T7ToJFAVQSnkAOYBLqRCf0xg5Ul85jRiRZGP79jBlih5f5gG87AHiC+2cfBV2N3eoMJzH3PYzrMs8+vfHrHLkBP74A06dgl6dDsKJJXpWqGd2q8My0pGUjHLJq5R6LOH7LEBDYF+yZj+hOxsA2gArRVzng/yhQ7qKYo8eUCj5RL327eHYMdg1EXLFwlOTQaWD0aBFW0POIAY8PZzI6zEMGmR1QMasWfDYY9Cw8FhwywSl37I6JCOdSUnmKQisUkrtBDaj+9B/UUqNUEq1TGgzHcitlDoEvAO8lzbhWuOjj3QJ0wED7DSIuQ67R0K+elCwiSNDe3jKDSqMINOtI8waMovp02HLFquDyriuXYOlS6F7xzN4HJ8FxTunbYllwyV53K+BiOwEKtnYPiTJ99FA29QNzTkcPKhvVL31FhQsaKfRvs8h+hzU+Sl9DS8r3BxyV6VNlpEULtCR3r0z8+ef6esluIrFi+HGDejTZAJcvQVl+lkdkpEOpYO+AWt99JGe4j9woJ0G0Rdg71go+hzkqebQ2B6ZUlBhBG43jvP9J9PZsOG/ia6GY82eDVUqXKZQ5NdQ7AXnHPJqOD2T0O/h4EE96uCNNyB/fjuN9o2F2OtQ8SOHxpZqCjaGPDWonuNTqleLYcAAiIiwOqiM5ehRWLMGxnWfiIqNgPKDrQ7JSKdMQr+HMWN033n//nYaRF+AA1+B70uQo6xDY0s1SkH591FRx5k7MoTTp++xJqqRJubMgayZr1M73wQo3MIU4TIemknodpw+rUcddO16r6vzzyE2CgI+dGhsqa7QM5AziMejP6FTxzjGj9cDd4y0Fx8PM2fC2Ncn4x57Ccq/b3VIRjpmErodEyZAbCz0s3dv6uZFXVGx2Avp9+r8NqX0x/yIA4zruwQ3NxhsPvU7xNq1cOpkNJ2rjoP8DdLffRjDqZiEbsO1azBpErRpA48/bqfRvvEQG5n+r85vK/IcZC9NnjMf8+67wvz5sDF5gQcj1c2YAb2azCCLOgMB5urceDQmodvw7bc6qdsdd37rsl50oFgbeKy8Q2NLM27uUG4QXNnB4C7LKFAA3nnH1HlJS9euwU+hN3m/9SeQt5aex2AYj8Ak9GRu3oTx46FBA6hSxU6j/V9BbITrXJ3f5tcOsvqS5choRo6Ev/7S46ONtPH99/By1Rnk8joJFYaZCQDGIzMJPZn58/UNUbtX57FRcOBLKNRcLzzgStw89bqV59fTtcVfVKgAgwbpRT2M1Dd31k2GtPkYyVNL958bxiMyCT0JEX0zNCAAGtlbIObwDLh5Acq7VHWD/zzeDTLlwn3/GD79FA4f1l1QRuratw/KZp5BwewnURWHmatzI1WYhJ7EunV6cYE+fez8fcXH6IlEeWvrPk9X5JEVSvWCkz/StOY+6tfXFSavXrU6MNcye+ZNBj/7MbdymKtzI/WYhJ7EhAmQK5cuoGhT+EKIDNer/riyUm+Ce2bUvnF89pkurfvZZ1YH5Tpu3YK4A9MpmvskmaoMM1fnRqoxCT3BsWMQGqrXeE66ZkUiiYc9oyFHgJ6I48q88kHxrnB0NlXKnaZdO/j8czh50urAXMOyn6J4q8FHXHKvba7OjVRlEnqCr7/WF0pvvGGnwanf4OpuKDcgfdQ7f1Rl3gWJhf0T+OgjPaNx6FCrg3IN5//6hkI5T5Oj7sfm6txIVRkgM91fZCRMmwbPPQdFi9pptO9z8C6i67ZkBD6PQ9Hn4eC3+BeJoFcvPUXdLCr9aMIPXaN16dEcjGyCe4EnrQ7HcDEmoaNLxl65om+G2nR5u16wt1QfPbQvoyjzLsRcgcPTGTwYsmUzJQEe1eFfx5PH5yLZaqTT6pyGU8vwCV1ET/OvWBFq2Ru4svdzPfqjxGsOjc1yeapB3idh33jy5Ipl4ED48Uf480+rA0ufYiMv8oTPONaHP0fBgGCrwzFcUIZP6Bs36qGKPXva6c6M+hfC50PxbpDpMYfHZ7my/SDqOBxfnLhq08CBpiTAwwj/7VOyZrrOjZIj7t/YMB5Chk/okybprgS7QxUPfAXEQ5kMumBv4ebgUwr2jiGrtzBsmL5C/+knqwNzUiEh4OcHbm766+0loCKPU/TGlywK60S9Z12k/o/hdDJ0Qr90CRYuhA4dwMfHRoOY63BwMhRpDdmKOzw+p6DcoOy7cHkrnFvDK69AqVK6Lz0uzurgnExIiB73Gh6uP8KEh+ufQ0K4/teHxMfBiZwj8MxAt2EMx8rQCX3mTF2Mq2dPOw2OzdE3Bcu87ciwnI9fR8icF/aOwcNDr2i0Z49eB9NI4v33ISrqzm1RUTCxP97n5vDl8rd48ZVi1sRmZAgZNqHHx8PkyVCzpr4hehcRXSI3VzDkqenw+JyKRxYo1RtOLYMru3nuOahaFYYM0SvVGwmOH7e5Ob7uGa5GPcbO2PfsD4s1jFSQYRP6qlV6EWi7V+dnVsC1fVDaXmGXDKbUG+DuDfvGohR8+qmeOfrVV1YH5kSK2bj6LgdugcJHoe/T+bWcjo/JyFAybEKfOhVy5tSrEtm0f4KeAl/sBYfG5bQy54bHX4FjIRD1L/XqwdNP6+6Xy5etDs5JjBp1Z90IBXRUnL5SlN+P9KKBmeVvpLEMmdAvXIClS6FjR/DystHg2kE49SuU6AHumR0en9Mq8w5InO6KAj75RFdhHD3a4ricRfv2MGUK+PrqT3Wtc0Ex4a3ZY+n2mhduGfKvzXCkDPkWmztXV7zr1s1Og4Nf6xmhJXs4NC6nl80firaBQ5Mh5hqBgTqHTZgAJ05YHZyTaN9eV3q7eQXae3Loai1+3tGWLl2sDszICO6b0JVSRZVSq5RSe5VSu5VSdw3IVkrVU0pdVUptT3gMSZtwH52IrttStaqdm6ExEXoRi2IvQJaCDo/P6ZXtDzHX4NBUAD76SP+bDnHa37hFdn8C0WfpPOELOnZU5MpldUBGRpCSK/RY4F0RKQtUB3oppcrZaLdORIISHk47FW7jRti9G1591U6Do7P1eqGlejs0rnQjdzDkr6+LlcXdxNcXeveGWbPgn3+sDs5JXD8K+8az83on/tofTN++VgdkZBT3TegiclpEtiZ8HwHsBQqndWBpZdo0yJoVXrJVNFEEDnydMFSxmsNjSzfKvQc3TsGxuYCeZJQjB7znoqvyPbBt/RDlTtfxH9O4MZSzdfljGGnggfrQlVJ+QCVgo43dNZRSO5RSvymlbM5tVkp1V0qFKaXCzp8//8DBPqqICFiwAF580c7M0HNr4NpeKGmvKLoBQIGGkLMy7PkM4uPIlUsvJr1smR4OmqGdXg4nfmBH7Ads3VfYXJ0bDpXihK6UygYsAfqKyLVku7cCviISCEwEQm0dQ0SmiEiwiATnzZv3YWN+aAsX6trndrtbDnwNmXJlnJrnD0spvUh2xAE4qX/VvXvrWvL9++tJWxlS3E0I6434lOSNie9Qpgw0aWJ1UEZGkqKErpTyRCfzEBH5Ifl+EbkmItcTvl8GeCql8qRqpKngu++gTBmoXt3Gzqh/4eRSPdbaI4vDY0t3ijwH2UroZflEyJJFD8PesgXmzbM6OIvsGw8RB9iTZSIbNmbmrbcwQxUNh0rJKBcFTAf2isjndtoUSGiHUqpqwnEvpmagj+rAAfjrL+ja1c7Ez0NT9bqhJcxQxRRxc9fL8V0Kg7N/QEgI7T/wpwphDO7yLze+W2B1hI4VeQJ2jYQirfjg6ybkzAmdOlkdlJHRpOT6oRbQEaifZFjiM0qpHkqp29mvDbBLKbUD+BJ4ScS5KmbPnKmvljp0sLEzPgYOT4GCTfXSa0bK+HfSQzv/6A3du+N2/BjjeJcTcYX54vW9/5WOzQi29gXiOZR9PKGhugvK5mLjhpGGlFV5Nzg4WMLCwhxyrrg4PXkvMBB+/dVGg+OLYP0LUPdnXf/bSLl9X8DWt2EksE9vasVSVlKfQ0WeIt+JLZaG5xAnf4K1z0LgJ3QZ8x6LFunKuXmcrtPRcAVKqS0iYnPJqwzRw7diBfz7r+5usengZMjqCwWfdmhcLqFEd7gCtPpv02cM4AZZGHKyu2VhOUxMBIT1gscqcDzru4SEwGuvmWRuWCNDJPTvvoNcuaBFCxs7r+3XC0CX6K77hY0H4+ENf+WECkBJvakUB+nF10zlVXbssDS6tLfzQ31DveoUxo3XK1e8+67FMRkZlssn9MuXITQU2rWDzLbqbB2aAspDrxlqPJynx8A1oPV/m4ZmGUPObDG89ZYLrz96MQwOTISSPTkv1Zk6Vd+jMTXPDau4fEJfuFCvSmSzOFLsDTgyE4o+B1nyOzgyF9K+G/i8BIHA44CvLzmnfsZHY7xYswYWL7Y6wDQQdws2vgJe+SHwYyZMgOhoGDDA6sCMjMzlb4rWrKlniO7caWO44tE5sKETNFgJ+Z9K81hcWsx1+MlPl0146ndA34yuXBmuXIF9+yCLKw3v3zlED1Os+zMXMjfH31/Xh//+e6sDM1xdhr0pevAgbNigxwPbHHt+cLJe0T5fPUeH5no8s0G5QXD6f3B2NQDu7rq07vHj8Nln1oaXqi5thd0f62GbhZszZoyegTxsmNWBGRmdSyf02bP12PP27W3svLwTLvyla56bJeZSR8k3IEth2D4oseO8Xj144QW9GMbhw9aGlyribsHfXfRqVlW+4MwZmDhR36MxRbgMq7lsQo+PhzlzoGFDKFTIRoNDU8AtM/h3dnhsLssjC1QYBhf/hn9/Ttz8+efg6akn26T7G6S7RsKVf6DqFMiUk9Gj9WIpQ4daHZhhuHBCX7dOT+7obCtfx0bBsTlQrC1kNisPpKriXXQ31o7BEB8HQOHCMHIk/PabXvov3Tq3HvZ8rF9j4eacPAmTJ+v3WMmSVgdnGC6c0GfPhmzZoFUrGzuPf69X3SmRASa+OJqbBwR+BFd3J9ZLB3jzTT1T96234Pp1C+N7WLeuwoaOkNUPqug1VUeO1J8EP/zQ2tAM4zaXTOhRUbBoEbRta6eexqEpkL0M5K3t8NgyhKLP69EuO96H2EgAPDxg0iQ4eTKddk+EvQlRJ6DGXPD04Z9/9GIpPXuCn5/VwRmG5pIJPTRUD1Xs2NHGziu74MIGePw1czM0rSg3qDIBbvyr19ZMUKMG9OgBX3yhlwJMN47N1582AoZA3hqIwDvv6FWa0uV/TobLcsmEPncuFCsGdeva2HloKrhl0kPOjLSTtyb4tYe9Y+H6kcTNn36q+9S7dtUTvpzetf2wqTvkrQXlBwO6wNuKFXqYoln82XAmLpfQz56F5cv1UMW7FheIvaEXgS76PHiZ6klpLuhT3ae+9b/iJtmzw7ffwt69ug/aqcVGwbo24O4FtRaAmwe3bulaLaVL6+4Ww3AmLpfQFyzQMxRt1j0/sRhirkCJ1xweV4bkXRjKv6+XqTuzInHz00/rkSGjR8O2bRbGdy8isLmnvrlbMwS8iwDwzTd6sZRx4/RQTMNwJi439T84WP8tbrFVhntFXYg6BS0OmP5zR4mLhl8DAAXP7NDVGYFLl6B8ecidGzZvdsKyAIem6q6WgKFQcRigZ7yWLw+1a+sFsc1byLBChpn6v3evTuQ2r86vHYBza6HEq+Yv0ZHcvaDaVLh+SNc/SZArly5rvHu3Exa0OrdW1zgv0BgC9JhEEd3FEh+vR+uYt5DhjFwqoc+dq/vNX37Zxs7D00G5m5mhVsj/lF6rdf94uPB34uamTaFvX/jqKzsrSVnh+lFY9zxk9YfaCxJr5C9YoK/KR40ywxQN5+UyXS7x8eDvD2XLwu+/J98ZA6FFIE9NqJOepyqmYzERuuvFwxue3qav3NEjXapV0ytK7dwJBQtaGeM1WF5TL1jRZCNkLwXAhQv6fVW8uF5o3N2sg2JYKEN0uaxfr/s4bY49//dniD4Hj7/q8LiMBJ4+uv7JtX16lZ8EmTPDvHm6WmH79hAba1F8cbf0urLX9sGTixKTuYie5Xrlip5IZJK54cxcJqHPnQtZs8Kzz9rYeWiargJYsInD4zKSKNREd73sHasXVk5QrpyuibJqFfTvb0FcEq8rKJ7+HzwxGQo0TNw1bZpeJGX4cKhQwYLYDOMBuERCv3lTT/Vv3VrXb7lD5Ak4/TsU76rHRBvWqjIeclbWC4skmXDUqZOu8/LFF7oOj8OIwJa3IHw+BI3WN80T/PMP9OkDjRrBe+85MCbDeEgukdCXLdMfiW2ObjnyHSDw+CuODsuwxd0LnlysywOsa6OHNSYYMwaeegq64aCb6wAAC3VJREFUd9dDGdOciO7+OfAVlO0HZf8bbhMZqeu458ihyzDfNUnNMJyQS7xN586F/PmhQYNkOyReJ/T8DSCbvyWxGTZk84cac+DyNtj4mv49oSfqfP89FCgAzZvrFafSjAhs6w+7R+l7K0GfJY5FjI+HV1+F/fshJES/twwjPUj3Cf3yZfjlF3jpJV3R7w5nV0HkMXi8mxWhGfdSuBlUHKmLXm15O3Hlizx59Cil+Hi9OMnJkw9wzJAQPabQzU1/DQmx3U7iYfMbsG8clHoTqn57x8Dy997TwxRHjbJxkWAYzkxE7vkAigKrgL3AbuAtG20U8CVwCNgJVL7fcatUqSKpYcoUERDZvNnGzvUviyzKKRJ7I1XOZaSy+HiRsL4iIYjsGHrHrrAwER8fkTJlRM6fT8Gx5s4V8fbWb4bbD29vvT2pmEiRdS/qc24bqGNIYsIE/dSePe/aZRhOAQgTe/na3g75L1kXvJ2gAR/gAFAuWZtngN8SEnt1YOP9jptaCb1uXZHSpW388d28JDI/s8jmN1PlPEYaiY8T2dBFJ9hdn9zxi1y9WsTLS6RCBZGTJ+9zHF/fO5P57Yev739troeLLKskEqJEdo++602zYIGIUiKtWonExqbaKzSMVHWvhH7fLhcROS0iWxO+j0i4Ui+crNmzwOyE8/0NPKaUSvMpIsePw5o1+mboXVOxj4ZA/E3T3eLslBtUnQrFXoQdg2DTa3pMOLr88S+/wLFjupb67t33OM7x4/fefnY1/B4MEYeg7k9QbuAdb5pJk/QM45o19bh4M97cSI8eqA9dKeUHVAKSL09QGDiR5OeT3J30UUp1V0qFKaXCzp8//2CR2rBkif7arp2NnUemQ85KkDPokc9jpDE3D6g1T1dmPDwdVjWC6AuA7sNeuxZiYnRRrNWr7RyjWDHb20sUgc294I+nIFNOPQO0cPPE3SJ6Cbk33oBnnoH//c8JC4UZRgqlOKErpbIBS4C+InIt+W4bT7mrpoCITBGRYBEJzps374NFakOfPvD333pK9h0ubYPL283VeXqi3PRapDVD4MJG+LUcHP4OJJ6gINiwQY9+qV9f37S8a3GMUaPuXG9QATUzwdAoODjp/9u7+9iq7jqO4+9veWjLUyXCEATadQEHwYQVgviQbYLbKjEQFhwoy9wgJnuIMdP9YWjUZTITH7aZxSUTJ4ILU6Yj2ExNYAiy4QosDMdkmYFChcEYTsXx/PT1j99Z13WlPYVzz+Gcfl7JTc6999fb77fn3G9/9/c793fgY1+Hz2+DqnGtTd5+O/TKFy+GBQvCla76908lW5GSiFXQzawPoZivcPdVHTTZT5g8fddI4MClh9e5Xr3COiAf0LwUysqhpqOuu1zWar4MNzXBgKtg8wJY82k4tIGa6vNs2QILF4arHk2ZAtu2tfm5+fNhyRKoGQ2fAn7UB+45DVVD4YYXYNJPoPd71Xr16rAU7jPPhP8FTzzRwVlSIjnTZUE3MwN+Abzm7g9foFkjcJsFU4Ej7n4wwTjjO3cS9q6AUTeHj9iSP4Mnwo2bYOov4VhzGC5prGVgcwM/f3Aja1a/waFD55k0Cb548wm2b2qBlqehdi08dBLuAa4eE851n7EjXA6PMLyybl04x3327NDj37oVFi3ScrhSDF2utmhmnwGeB3YA56OHFwGjAdz98ajo/xSoB44Dd7h7p0spluoCF7SshE3zYNra963JITl19hjsWw17n4Q317Z+CcnLyjl9tg/lZUdbmx4/N5h3+tdTefUt9K6eSUVlGYcPh578tm1hsnPnThg6FO69F+67T1cdkvzpbLXFLj9kuvsLdDxG3raNE/pF2du9FPqNhmHTso5EktC7P1w5P9xOHApzI8f2YEebKT9/hlM2jPVNQ1n13DiWPzuF02cufEhPmgTLlsHcuVBRkV4KImkp1qjhsX2hFzfhO2GSTYqlchhUvn/FzHKgvg7q74Yf/y+cEdPcDCdOhFtVFdTVwcSJYVukyIpV0PcsBxxqb886EsnAoEFhfFykpypON7Z1Ia5pMKAm62hERFKX34LefiGmp74d1teu1TK5ItIz5bOgr1gRFs1uaQnnorW0wIs/AK+EUbOzjk5EJBP5LOgNDXD8+Hv3K4G6c7C1V7gIsYhID5TPgt5+IaZPABXAH4521FpEpEfIZ0FvvxDTdcAbwJkLLNAkItID5LOgt12IaTgwFvhrH3jw+1lGJSKSqXwW9HcXYqquhmuBc8CcR8LjIiI9VD4LOoTi3bwb5o6AUV+A+ZfHygMiIlnJb0EHeHMNnDgAV92RdSQiIpnLd0Ev6wsjZsAIfd9bRCTfa7l8ZHq4iYhIznvoIiLSSgVdRKQgVNBFRApCBV1EpCBU0EVECkIFXUSkIFTQRUQKQgVdRKQgzN2z+cVmh4GWBF5qCPCvBF4nL5RvcfWkXEH5Xqxqdx/a0ROZFfSkmNlL7j456zjSonyLqyflCsq3FDTkIiJSECroIiIFUYSCviTrAFKmfIurJ+UKyjdxuR9DFxGRoAg9dBERQQVdRKQwclPQzazezF43s11m9q0Oni83s5XR85vNrCb9KJMTI99vmNlOM3vFzNaZWXUWcSalq3zbtJtjZm5muT3dLU6uZnZLtH//bmZPpR1jkmIcy6PNbL2ZvRwdzzOyiDMJZrbUzN4ys1cv8LyZ2aPR3+IVM6tLNAB3v+xvQC9gN1AL9AX+Boxv1+Zu4PFoex6wMuu4S5zvZ4F+0fZdRc83ajcQ2Ag0AZOzjruE+3YM8DIwOLp/RdZxlzjfJcBd0fZ4YG/WcV9CvtcCdcCrF3h+BvAnwICpwOYkf39eeuhTgF3u3uzup4HfALPatZkFLI+2fwdMNzNLMcYkdZmvu6939+PR3SZgZMoxJinO/gX4HvBD4GSawSUsTq5fBR5z9/8AuPtbKceYpDj5OjAo2q4CDqQYX6LcfSPw706azAJ+5UET8CEzG57U789LQf8osK/N/f3RYx22cfezwBHgw6lEl7w4+ba1kPBfP6+6zNfMrgFGufuzaQZWAnH27VhgrJltMrMmM6tPLbrkxcn3fuBWM9sP/BH4WjqhZaK77+1uyctFojvqabc/3zJOm7yInYuZ3QpMBq4raUSl1Wm+ZlYGPALcnlZAJRRn3/YmDLtcT/jk9byZTXD3/5Y4tlKIk++XgGXu/pCZfRJ4Msr3fOnDS11J61Reeuj7gVFt7o/kgx/LWtuYWW/CR7fOPvpczuLki5l9DmgAZrr7qZRiK4Wu8h0ITAA2mNlewthjY04nRuMey7939zPuvgd4nVDg8yhOvguBpwHc/UWggrCQVRHFem9frLwU9K3AGDO70sz6EiY9G9u1aQS+Em3PAf7s0SxEDnWZbzQE8TNCMc/zGCt0ka+7H3H3Ie5e4+41hDmDme7+UjbhXpI4x/JqwqQ3ZjaEMATTnGqUyYmT7z+B6QBmNo5Q0A+nGmV6GoHborNdpgJH3P1gYq+e9axwN2aPZwD/IMyYN0SPPUB4Y0M4CH4L7AK2ALVZx1zifJ8DDgHbo1tj1jGXMt92bTeQ07NcYu5bAx4GdgI7gHlZx1zifMcDmwhnwGwHbsw65kvI9dfAQeAMoTe+ELgTuLPNvn0s+lvsSPo41lf/RUQKIi9DLiIi0gUVdBGRglBBFxEpCBV0EZGCUEEXESkIFXQRkYJQQRcRKQgVdJFItCb3DdH2YjN7NOuYRLojL4tziaThu8ADZnYFcA0wM+N4RLpF3xQVacPM/gIMAK5393eyjkekOzTkIhIxs48Dw4FTKuaSRyroIkB01ZgVhCvKHDOzmzIOSaTbVNClxzOzfsAq4Jvu/hrhUnf3ZxqUyEXQGLqISEGohy4iUhAq6CIiBaGCLiJSECroIiIFoYIuIlIQKugiIgWhgi4iUhD/B9mPY6PyZCXDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Plotting\n",
    "plt.plot(d,g(d),color='b',label=\"Defining function\")\n",
    "plt.scatter(x,y,color='r',label=\"Generated data\")\n",
    "plt.plot(d,y_approx,color='orange',label=\"Approximating function\")\n",
    "plt.xlabel(r\"$x$\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now adjust both the generating function $g$, and the parameters ($lam$ and $l$).\n",
    "\n",
    "Try:\n",
    "<ul>\n",
    "    <li>Changing the generating function to the ramp function $g(x) = x$ with $l=0.1$.</li>\n",
    "    <li>Now increase $l$ to 1, does the approximation improve? why?</li>\n",
    "    <li>Now go back to the original generating function, with $l=0.1$, $lam=0.01$.  Adjust $lam$ between $0.001 - 10$ and observe the fit.  what is $\\lambda$ doing?</li>\n",
    "    <li>Adjust $l$ between 0.001 and 1, when is the approximation the best?  how does this relate to $g(x)$?</li>\n",
    "</ul>\n",
    "\n"
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}