{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spline Regression\n",
    "\n",
    "In the first steps, I followed the example by Christopher Krapu (https://ckrapu.github.io/2018/07/09/Spline-Regression.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import theano\n",
    "import theano.tensor as tt\n",
    "import pymc3 as pm\n",
    "\n",
    "import pandas as pd\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "theano.config.floatX = 'float64'\n",
    "theano.config.compute_test_value = 'ignore'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_B_spline_deg_zero_degree_basis_fns(breaks, x):\n",
    "    \"\"\"Build B spline 0 order basis coefficients with knots at 'breaks'. \n",
    "    N_{i,0}(x) = { 1 if u_i <= x < u_{i+1}, 0 otherwise }\n",
    "    \"\"\"\n",
    "    expr = []\n",
    "    expr.append(tt.switch(x<breaks[1], 1, 0))\n",
    "    for i in range(1, len(breaks)-2):\n",
    "        l_break = breaks[i]\n",
    "        u_break = breaks[i+1]\n",
    "        expr.append(\n",
    "            tt.switch((x>=l_break)&(x<u_break), 1, 0) )\n",
    "    expr.append( tt.switch(x>=breaks[-2], 1, 0) )\n",
    "    return expr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_B_spline_higher_degree_basis_fns(\n",
    "        breaks, prev_degree_coefs, degree, x):\n",
    "    \"\"\"Build the higer order B spline basis coefficients\n",
    "    N_{i,p}(x) = ((x-u_i)/(u_{i+p}-u_i))N_{i,p-1}(x) \\\n",
    "               + ((u_{i+p+1}-x)/(u_{i+p+1}-u_{i+1}))N_{i+1,p-1}(x)\n",
    "    \"\"\"\n",
    "    assert degree > 0\n",
    "    coefs = []\n",
    "    for i in range(len(prev_degree_coefs)-1):\n",
    "        alpha1 = (x-breaks[i])/(breaks[i+degree]-breaks[i]+1e-12)\n",
    "        alpha2 = (breaks[i+degree+1]-x)/(breaks[i+degree+1]-breaks[i+1]+1e-12)\n",
    "        coef = alpha1*prev_degree_coefs[i] + alpha2*prev_degree_coefs[i+1]\n",
    "        coefs.append(coef)\n",
    "    return coefs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_B_spline_basis_fns(breaks, max_degree, x):\n",
    "    curr_basis_coefs = build_B_spline_deg_zero_degree_basis_fns(breaks, x)\n",
    "    for degree in range(1, max_degree+1):\n",
    "        curr_basis_coefs = build_B_spline_higher_degree_basis_fns(\n",
    "            breaks, curr_basis_coefs, degree, x)\n",
    "    return curr_basis_coefs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "def spline_fn_expr(breaks, intercepts, degree, x):\n",
    "    basis_fns = build_B_spline_basis_fns(breaks, degree, x)\n",
    "    spline = 0\n",
    "    for i, basis in enumerate(basis_fns):\n",
    "        spline += intercepts[i]*basis\n",
    "    return spline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compile_spline(data,n_bins,degree,intercepts):\n",
    "    breaks = np.histogram(data, n_bins)[1][1:-1]\n",
    "    for i in range(degree+1):\n",
    "        breaks = np.insert(breaks, 0, data.min()-1e-6)\n",
    "        breaks = np.append(breaks, data.max()+1e-6)\n",
    "    xs = tt.vector(dtype=theano.config.floatX)\n",
    "    f = theano.function([intercepts, xs],spline_fn_expr(breaks, intercepts, degree, xs))\n",
    "    return f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "domain = np.linspace(-4,4,100)\n",
    "n_bins = 2\n",
    "degree = 4\n",
    "coefficients = tt.vector(dtype=theano.config.floatX)\n",
    "spline = compile_spline(domain,n_bins,degree,coefficients)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.asarray([-3.69,  1.60, -2.05, -1.63, -3.07,\n",
    "       -1.06,  3.64, -1.61, -2.58, -3.57,\n",
    "       -1.55, -1.93, -3.94,  3.51, -0.17,\n",
    "       -3.92,  0.52,  3.06 ,  3.40, -0.21])\n",
    "true_coef = np.asarray([ 0.82, -0.34,  1.75, -0.78,  0.25,\n",
    "       -0.76,  0.59,  1.13 ,  1.32])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'f(x)')"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcEAAAEfCAYAAAA5j323AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XlcTP37P/BXC5UlWdJGmzWRLBVRWcqayoci620pIW5JiEhkXyNLkj1bpGQLCaWFbhJZChVRUVRItP3+8Ot8jZk2zdI01/Px6HHf856zXPM2M9ecc97nfYnl5uaWgRBCCBFB4oIOgBBCCBEUSoKEEEJEFiVBQgghIouSICGEEJFFSZAQQojIoiRICCFEZFESJIQQIrIoCRJCCBFZlAS5LDk5WdAh1DvUp9xHfcob1K/cx+s+pSRICCFEZFESJIQQIrIoCRJCCBFZlAQJIYSILEqChBBCRBYlQUIIISKLkiAPpKWlwc7ODubm5rCzs0NaWpqgQyKEEMKBpKADqG/evXsHJycnpKSkMG1xcXEICgqCmpqaACMjhBDyJzoS5LJ9+/axJEAASElJgaenp4AiIoQQUhFKglz28eNHju2ZmZl8joQQQkhVKAlymby8PMd2RUVFPkdCCCGkKpQEuczBwQEaGhosbRoaGnBzcxNQRIQQQipCA2O4TEVFBUFBQfD09ERmZiYUFRXh5uZGg2IIIaQOoiTIA2pqavD19RV0GIQQQqpAp0MJIYSILEqChBBCRBYlQUIIISKLkiAhhBCRRUmQEEKIyKIkSAghRGRREiSEECKyKAkSQggRWZQECSGEiCxKgoQQQkQWJUFCCCEii5IgIYQQkUVJkBBCiMiiJEgIIURkURIkhBAisigJEkIIEVmUBAkhhIgsoUqCd+/exfjx46GlpQU5OTn4+/tXuU5iYiJGjBgBRUVFaGlpYePGjSgrK+NDtIQQQuo6oUqC3759Q5cuXbBhwwbIyMhUuXx+fj5Gjx6N1q1b4+bNm9iwYQN27doFb29vPkRLCCGkrpMUdAA1MWTIEAwZMgQAMGfOnCqXDwgIwPfv37F3717IyMigS5cuSEpKwp49e+Do6AgxMTFeh0wIIaQOE6ojwZq6d+8e+vbty3LUOHjwYGRkZCAtLU2AkRFCCKkLhOpIsKY+fPgAZWVlljZ5eXnmOXV1dY7rJScn12q/tV2fsKM+5T7qU96gfuW+2vRphw4dKn2+XidBAGynPMsHxVR2KrSqTqtMcnJyrdYn7KhPuY/6lDeoX7mP131ar0+Htm7dGh8+fGBpy87OBvB/R4SEEEJEV71Ogvr6+oiOjkZhYSHTFh4eDiUlJaipqQkwMkIIIXWBUCXBr1+/IiEhAQkJCSgtLUV6ejoSEhLw9u1bAICHhwcsLCyY5ceOHQsZGRnMmTMHT58+xYULF7Bjxw7MmTOHRoYSQggRriT48OFDGBsbw9jYGN+/f8f69ethbGyMdevWAQAyMzORkpLCLN+sWTOcP38eGRkZGDhwIFxcXDB37lw4OjoK6iUQQgipQ4RqYIyRkRFyc3MrfH7v3r1sbdra2rhy5QovwyKEECKkhOpIkBBCCOEmSoKEEEJEFiVBQgghIouSICGEEJFFSZAQQojIoiRICCFEZFESJIQQIrIoCZJ6LS0tDXZ2djA3N4ednR2V0CKEsBCqm+UJqYm0tDRYWVmxzCIUFxeHoKAgmjuWEAKAjgRJPebp6cmSAAEgJSUFnp6eAoqIEFLXUBIk9VZGRgbH9szMTD5HQgipqygJcllGRgZu374t6DAIACUlJY7tioqKfI6EEFJXURLkgsLCQpw6dQoWFhawsLCAg4MDSkpKBB2WyHNzc4OGhgZLm4aGBtzc3AQUESGkrqEkyAXFxcVwdnbGnTt3ANDRYF2hpqaGoKAgWFtbw8jICNbW1jQohhDCgkaHckGTJk1gaWmJEydOMG0nTpzAoEGDBBgVAX4lQl9fX0GHQQipo+hIkEtsbW1ZHl+8eBF5eXkCioYQQkh1UBLkkn79+kFVVZV5XFhYiKCgIAFGRAghpCqUBLlEXFwc48ePZ2n7/fQoIYSQuoeSIBf9eUo0NjYWr169ElA0wuvPqc7evXsn6JAIIfUUJUEu0tDQQI8ePVjaTp48KaBohFP5VGcBAQGIjIxEQEAAHB0dac5PQghPUBLkMnNzc5bHp06donsGa4DTVGfp6ek01RkhhCcoCXLZoEGDICMjwzxOT0/HtWvXBBiRcOHVVGdUTYIQwgklQS5r0qQJRo8ezdJG96lVHy+mOuN0itXKyooSISGEkiAv2Nvbszy+efMmkpOTBRSNcOE01VmbNm1qNdUZVZMghFSEkiAP6OrqQk9Pj6XtwIEDAopGuHCa6szb27tWU51RNQlCSEUoCfKInZ0dy+OTJ0/i69evAopGuJRPdRYSEgJfX1+oqKjUantUTYIQUhFKgjxiaWkJeXl55nF+fj7OnDkjwIhEF1WTIIRUhJIgj0hJSWHq1Kksbb6+vigrKxNQRKKLqkkQQipCVSR4aNq0adi+fTtzn+CzZ88QEREBY2NjAUcmeqiaBCGEEzoS5CEVFRW2m+d37twpoGgIIYT8iZIgj82aNYvl8Y0bNxAfHy+gaAghhPyOkiCP9e3bF3369GFp27p1q4CiIYQQ8jtKgjwmJiYGZ2dnlraQkBA8f/5cQBERQggpR0mQD0xNTdG9e3eWtu3btwsoGkIIIeUoCfKBmJgYFi5cyNJ29uxZpKamCiYgQgghACgJ8s2oUaPQsWNH5nFJSQm8vLwEGBEhhBBKgnwiLi4OJycnlrZjx47h9evXAoqI1BaVZyJE+FES5KOxY8eyzFJSXFwMDw8PlmXoi1U4UHkmQuoHSoJ81KBBAyxfvpylLTg4GPfu3QNAX6zChMozEVI/UBLks7Fjx7KNFF2xYgXKysroi1WIUHkmQuoHSoJ8Ji4ujjVr1rC0xcbGIiQkhL5YhQiVZyKkfhC6JHjgwAHo6OhAQUEBJiYmiIqKqnDZiIgIyMnJsf0lJSXxMWJ2xsbGGDp0KEubh4cHFBQUOC5PX6x1D5VnIqR+EKokGBgYiKVLl8LZ2Rl37tyBvr4+rK2t8fbt20rXi4mJwYsXL5i/du3a8Sniinl4eEBc/P+6/9WrV2jbti19sQoJKs9ESP0gVKWUdu/ejQkTJjB1+jZv3oywsDAcPHgQ7u7uFa4nLy+Pli1b8ivMauncuTOmTJmCw4cPM2379u3D2bNnceTIEWRmZkJRURFubm70xVpHfPv2DdnZ2fj+/TsKCgpQUFCAKVOmQEJCApKSkvj8+TNKS0vRvHlzyMrKsvzIIYTUTUKTBH/+/In4+HjMmzePpX3QoEGIjY2tdN0BAwbg58+f6NSpExYtWlRn6vktX74cwcHB+Pz5MwCgsLAQmzdvRlBQEMTExAQcnegpLS3Fq1ev8PTpU7x+/RqvXr3C69evkZGRgQ8fPuDbt2/V3paYmBhat24NVVVVqKqqQk1NDVpaWtDW1kaHDh3QoEEDHr4SQkh1CU0SzMnJQUlJCeTl5Vna5eXl8eHDB47rKCoqYtu2bejZsyd+/vyJ06dPw9LSEhcvXkS/fv14Eue7d++wadMmZGRkQElJqdIjOXl5eXh6emLu3LlM2+3bt3HixAlMnDiRJ/GR/5OVlYWYmBhER0cjPj4ejx8/rlGiq0xZWRmysrKQlZWF+/fvszzXoEEDdO3aFQYGBkyVkYquBxNCeEssNze3TNBBVEdGRga0tLRw+fJlGBoaMu0bNmzAuXPn2L5oKmJtbQ0JCQmcOnWqwmWSk5P/KsZ3797B0dER6enpTFubNm3g7e0NFRUVjuuUlZVh7ty5LPHLysrizJkzde4UrrD7+vUr4uLiEBMTg3v37lV5LZmf2rdvD0NDQxgaGqJ79+6QlBSa36eE1GkdOnSo9Hmh+aS1bNkSEhISbEd92dnZbEeHlenVqxcCAwMrXaaqTqvIpk2bWBIgAKSnp8Pf3x++vr4Vrrd//34YGhri+/fvAID8/Hx4e3vj2LFjdFoUv36U/O2/ydu3b3Hp0iVcvHgRMTExKC4u/us4GjZsCHl5eTRp0gQyMjKQkZGBuLg4SktLUVxcjMLCQuTm5iI3Nxdfvnyp0bZfvnyJly9f4ujRo2jevDnMzc0xevRoGBkZ8eTUaW36lFSM+pX7eN2nQpMEGzZsCF1dXYSHh8PKyoppDw8Ph4WFRbW38/jxY56devrb+/w0NDSwbNkyrFixgmm7ePEifH19YW9vz9UYRUF6ejoCAwMRGBiI+Pj4Gq0rJycHHR0ddOrUCe3atYOmpibU1NSgoKCAZs2aVftHyc+fP/H+/Xu8efMGaWlpePnyJRITE5GYmFjh+6Tc58+fcezYMRw7dgwtWrSAtbU1Jk+ejK5du9botRBCqiY0SRAA5s6di1mzZqFXr14wMDDAwYMHkZmZiWnTpgEAZs2aBQDw8fEBAOzZsweqqqrQ0tLCz58/cebMGVy6dAlHjx7lSXy1uYF69uzZCAwMxMOHD5m25cuXQ09PDz169OBajPVVXl4ezp8/jzNnzlR67+jvJCUl0b17d/Tt2xf6+vrQ1dVF27ZtuXL03bBhQ6irq0NdXZ3tuQ8fPiA2NhaxsbGIjo7Gw4cPUVpaynE7nz59go+PD3x8fKCrq4upU6fCxsYGjRs3rnWMhBAhuiZY7sCBA/Dy8kJWVha0tLSwbt06ZpDLyJEjAQCXLl0CAHh5eeHw4cPIyMiAtLQ0tLS04OTkhCFDhvAktrS0NIwcOZLllKiGhka17x9LSUmBiYkJ8vPzmTY1NTXcvn0bcnJyPIlZGFR0OqS0tBQRERHw9/fHhQsXUFhYWOW2OnTogEGDBsHU1BSGhoZ1Ipl8/vwZN2/exPXr1xEaGsqMFq5Is2bNMGXKFMycOfOvb5+h03a8Qf3KfbzuU6FLgnXdrVu34O/vX637/NLS0uDp6ckykjQ+Pp65D7Kcubm5SF8f/PND8PHjR/j7++PIkSNsc61y0rRpU2hoaGDt2rUwMjLiZai1VlRUhIiICAQGBiIkJAR5eXkVLisuLg4rKyssWLAAOjo6NdoPfVnzBvUr91ESFDLV/Qcrrxjx+5d4+VGjt7c320CaBQsWYNWqVdwOVygkJyejffv2iImJga+vL0JCQlBUVFTpOl27dsX79+/x6dMnpq0mR+V1QWFhIS5fvozjx48jPDwcZWUVf1TNzMywcOFCKCsrs/2w4vR66cuaN6hfuY+SoJCp7j+YnZ0dAgIC2Nqtra3h7e2NIUOG4NGjRyzPbdiwAQ4ODlyLVRgUFBTA29sbwcHBSExMrHTZtm3bwtbWFuPHj8f69esr7N/KRurWVW/evMHhw4dx+PBhlsT+JxkZGWaUMVBx4qcva96gfuU+XvcpzeskIJWNJJWSksKRI0fQqlUrludcXV1x/vx5foQncG/evIG7uzu6dOmCdevWVZgApaSkYG1tjeDgYDx69AjLli2DpqZmvavIoaqqipUrVyIxMRG7du2ClpYWx+V+T4AAleIipCqUBAWkqpGk6urqCAgIYBm4UVZWhlmzZuHmzZt8iZHfysrKEBUVhSlTpkBXVxdeXl7Izc3luGz79u3h6emJZ8+ewdfXFyYmJixzddbXUkcyMjKYPHky7t69i5MnT0JfX7/KdVJTU3kfGCFCipKggFSnFE+PHj1w5MgRltlDfv78iXHjxtWrI8LCwkL4+/vDxMQEI0aMwIULFzjeMiAuLo6RI0ciKCgI9+/fh6OjI1q0aMFxm/W91JG4uDiGDx+O0NDQKqcB/O+//+Du7s4y6pgQ8gtdE+Sympy/Lh8dWtVI0pMnT2L27NksbWJiYtiyZQtmzJjBlbgFISMjA35+fjh8+DCys7MrXK5Zs2aYNm0apk+fDlVV1Wpvv7r9Wx+UlZUhICAA8+bNw48fPzguIy8vDzc3N0yaNAmvX7+ma1c8QNcEuY8GxggZXv2D7du3D0uXLmVrX7RoEVxdXSEhIcH1ffJCWVkZoqOj4efnh+Dg4EqnMevSpQscHBzQo0cPdOvWjY9RCq/U1FQ4OjriwYMHKCgo4LiMtrY25s+fj3HjxvE5uvqPkiD30cAYAgBwcHDAvn372JLdli1bYGlpiffv3wsosur58uULDh06hP79+2PEiBE4d+4cxwQoJiaG4cOH48KFC7h79y6mTJkCaWlpAUQsnNTV1XHx4kW8ffsWu3btgrKyMtsyiYmJmDVrFuzs7Kqcwo2Q+q5GR4Lx8fGIjo5GUlIScnJyICYmhpYtW6Jjx44wMDCg6b3A+18toaGh+Oeff9hGAbZo0QJ79uzBsGHDeLbvmiorK8ODBw9w+PBhBAYGVlqmSFZWFhMnToS9vT3btTz6df33vn37hp07d2Lnzp1s7xkAaNKkCZYsWQIHBweqccgF9F7lPoGfDv348SN8fX1x8uRJvHv3DmVlZWjQoAGaN2+OsrIy5ObmoqioCGJiYlBWVoatrS3s7OzQunVrngVdl/HjQxAbG4uJEydyvI42atQorF69mi2R8NO7d+8QEBCA06dP49mzZ5Uu26FDB9jb28PW1hZNmjThuAx9sdReeno6PDw8ON47CQBaWlrYunUrS5kyUnP0XuU+gSZBDw8P7N+/H40bN4alpSVMTEzQs2dPtlMs7969w4MHDxAeHo6LFy/i27dvsLe3h7u7O88Cr6v49SHIzMzErFmzcPv2bbbnGjZsiFmzZmHBggV8q0mYlZWFixcvIjg4GBEREZXOblI+stHOzg4mJiZVTgdHXyzcEx0dDRcXFzx58oTj87a2tlizZg3bPaqkeui9yn0CTYKDBg3Cv//+i1GjRrHcg1WZ0tJShISEwMvLq97ez1YZfn4ISkpK4OXlhbVr16KkpITteSkpKVhZWcHc3BzBwcHIzMysstp9dZWWliIhIQHh4eG4du0aYmJiKk18AKCsrIyJEydi6tSpaNOmTbX3RV8s3FVSUoKNGzfCx8eH49ykzZs3h4eHByZNmlTtzz35hd6r3Cfw06GkZgTxIYiLi4OLiwtLGabKqKmp4cKFCzVKhN++fUN8fDz+++8/xMXF4e7du8jJyalyPQkJCZiZmWHq1KkwMzP7q4rp9MXCfcnJyZCTk4O7uztOnDjBcZm+ffti+/bt6Ny5M5+jE170XuU+SoJCRlAfgtLSUgQEBGD16tV49+5dlcs3btwYhoaGaNeuHZo3b45GjRqhUaNGKC0tRX5+Pr58+YLs7GykpKQgJSWlxqNPdXV1MX78eIwZMwby8vJ/+7IA0BcLL/zep1FRUVi0aBGePn3KtlyDBg3w77//wtnZGTIyMvwOU+jQe5X76lQS3LRpExYtWlThKZLPnz/DyckJhw8f5lZ8QkfQH4KCggL4+fnhwIEDSEtL4+u+u3fvDktLS1hYWKB9+/Zc266g+7Q++rNPi4qKsGfPHmzYsIHjKFJNTU1s27YNAwYM4GOUwofeq9xXp+4TXL9+PczMzJCcnMz23JUrV9CnTx/cuHGDa8GRmmvUqBHmzZuHhw8fIiAgoMI5NLlBVlYW5ubm2LZtGx49eoTbt29j4cKFXE2AhD/Kj/hiYmIwdOhQtudfv34NKysr2NvbVzq7DyHCpkZJMCAgAO/fv4eJiQn27t0L4NdN0HPmzMHEiROhrq6OO3fu8CRQUjPi4uIwMzPD1atXuTJVmLi4ONq1awcbGxts3LgRx48fh5mZGXJzcxEdHc2FiCuXlpYGOzs7mJubw87Oju9HuaJCTU0Np06dwpEjRzhONn7mzBno6enh2LFjVQ6EIkQY1PiaYG5uLhYtWoTAwEAYGBjg7du3+PjxI1xdXTF//nyRH01WF0+H/DmH5sKFC1FUVISXL1/izZs3+Pr1KwoKClBQUAAxMTHIysoyf6qqqtDU1ISqqiqkpKSY7VVUEJgXc3PeunULTk5OfNufKKjO+zQvLw9r1qyBn58fx4TXr18/bN++HR07duRVmEKnLn7+hV2duiZY7vv377CwsEBcXBzExMSwZs0azJ07lxfxCR1R+BBUVhCYFwVrx48fj6tXr/Jtf6KgJu/T+/fv499//61w4MyCBQvg7OxM09tBND7//FanrgkCv4bjGxsbIz4+Ho6OjujWrRtWrFgBFxeXCifsJfULvwvWfvz4ka/7I6z09PRw+/ZteHh4sI0QLSoqwubNm2FoaIjw8HABRUjI36tRElyzZg2GDx8OMTExhIaGYs2aNQgLC8PChQtx+PBhGBkZITY2llexkjqC3wVrK7rFQtgL5AqT8oEz0dHRMDU1ZXv+9evXGD16NGbOnImsrCwBREjI36lREty+fTtmzJiBO3fuoGfPngAASUlJuLm54dq1a5CUlMTIkSN5EiipO/hdsNbBwaFeF8gVJurq6ggICMChQ4egoKDA9vzZs2ehp6cHHx+fSstkEVJXSCxdunRVdRfu27cvZsyYwXHWDyUlJUyZMgUFBQUYNGgQN2MUKp8+feLbfJ2CIicnh+HDhyMnJwctW7aEgYEB9uzZw7NBKkVFRZg4cSLf9icKavM+FRMTg5aWFiZPnowvX74gPj6e5fkfP37gxo0buHLlCrp27QoVFRVuhCwUROHzz2+87lOaMYbL6MI491Gfch83+zQuLg5OTk54/Pgxx+cnTpyIVatW1XrmIGFA71XuE+jAmNzc3L/ecG3WJYQIj969eyM8PBzr169H06ZN2Z739/dHr169sHfvXjpFSuqcSpNgt27d4OHhUaMbk1NTU7FixQro6OjUOjhCiHCQlJTE7Nmzce/ePVhbW7M9n5+fD1dXVxgbG3Ms/0WIoFSaBPfu3YvQ0FD06NEDZmZm8PT0xKVLl/DkyROkp6fj7du3ePz4MS5evAhPT08MHjwYPXv2RFhYGDOjDCFEdCgpKcHX1xcXLlxAp06d2J5/+vQpLC0tYWtri1evXgkgQkJYVXlNsKysDNeuXYO/vz+uX7+OwsJCtiKoZWVlkJaWxuDBgzF58mQMGTKkykKp9RVdE+A+6lPu40efFhUVYf/+/diwYQO+fPnC9nyDBg1gb2+PRYsWoXnz5jyNhV/ovcp9Ap0x5u7du+jUqRNTZbqoqAgPHz5EUlISPn36BABo0aIFOnXqBF1dXTRo0IBngQoL+hBwH/Up9/GzT7OysrBq1SqcPHmS4/NycnJYtGgR7OzsmKn5hBW9V7lPoANjRo0axTILRO/evZGdnY1JkyZh/vz5mD9/PiZNmgQ9PT1KgIQQjhQUFLB3716EhYVBX1+f7fnc3Fy4ublBX18fZ86cQWlpqQCiJKKq0iTYuHFjfPv2jXn85s0blseEEFJdvXr1QmhoKA4cOIA2bdqwPZ+WlgZ7e3sYGRkhNDSUqlQQvmC/6/03Xbt2hZeXF378+AFZWVkAQHR0dJXDnG1tbbkXISGk3hATE8PYsWMxcuRI7Nu3D9u3b0d+fj7LMomJiRg3bhwMDAywbNkyGBsbi+wYA8J7lV4TjI+Px7Rp05CamvprYTGxKn+diYmJMdcLRRFdE+A+6lPu40WflpfsysjIgJKSEtzc3Kqc1Sc7OxsbN27EoUOHKvxxbWhoCFdXV/Tv37/OJ0N6r3KfwEsplZaWIj09HR8/foSpqSlcXV2rnBatd+/eXA1SmNCHgPuoT7mP231a2xqTKSkpWLduHccSXeX69u2LRYsWYdCgQXU2GdJ7lfsEngR/N2fOHEyfPl2kk1xV6EPAfdSn3MftPuVWjcmEhASsXbsWoaGhFS7To0cPLFiwAObm5pCQkPireHmF3qvcV6fqCe7Zs4cSICGEDbdqTOro6OD06dO4ceMGBg8ezHGZhw8fYurUqdDT04Ofnx++f/9e43gJKVfjorqEEPInbteY7N27N86dO4erV69yrF8I/Kph6OzsDG1tbaxevRrv3r37q30R0UZJkBBSa7yqMdmnTx+cPXsW4eHhFdYq/fTpE7Zt2wYdHR1Mnz4dd+/epdsrSLVREiSE1JqamhqCgoJgbW0NIyMjWFtbV3tQTHX06NED/v7+iI6OxsSJEzlOzlFSUoLAwECMHDkSffv2hY+PD1WzIVWieoJcRhfGuY/6lPuEvU/fv38PHx8fHD58GHl5eRUuJy0tjVGjRmHSpEkwMjKCuDhvf/cLe7/WRXVqYAwhhNQFysrK8PDwwNOnT7Ft2zaOFSsAoLCwEAEBAbC0tISuri48PT2RlJTE52hJXUZJkBAitBo3bozp06cjJiYGQUFBsLCwqPC2iTdv3mDLli3Q19fHgAEDsGvXLrx9+5bPEZO6RuiS4IEDB6CjowMFBQWYmJggKiqq0uUjIyNhYmICBQUFdO/eHQcPHuRTpIQQfhETE8OAAQNw9OhRPHnyBMuXL4eqqmqFy8fHx2PFihXo1q0bzMzMsHv37hoVDyf1h1AlwcDAQCxduhTOzs64c+cO9PX1YW1tXeGvudTUVNjY2EBfXx937tzBwoULsXjxYgQHB/M5ckL4Jy0tDXZ2djA3N4ednZ3IfbkrKSnBxcUF8fHxCA4Oho2NDaSlpStc/v79+1i+fDm6d+8OIyMjbNiwAY8ePaIRpiJCqAbGDB48GNra2ti5cyfT1rNnT1haWsLd3Z1teXd3d4SEhODBgwdM27x58/D8+XNcv36dJzHShXHuoz6tvupOXyZqfZqfn4+QkBCcOXMGd+7cqVaCU1ZWxtChQzFkyBAYGRmhSZMmVa4jav3KDzQw5v/7+fMn4uPj2eYtHTRoEGJjYzmuc+/ePbblBw8ejIcPH6KoqIhnsRIiKJ6eniwJEPg1L6enp6eAIqobZGVlMXHiRAQHByMxMRHr1q3jWNvwd+/fv8ehQ4dga2sLTU1NWFpaYufOnXj8+DHVPKxHKi2lVJfk5OSgpKQE8vLyLO3y8vL48OEDx3U+fPiAAQMGsC1fXFyMnJycCmezSE5OrlWstV2fsKM+rZ7Xr19X2P5nH4pyn5qZmcHMzAyZmZm4efMmbt26hUePHlWY3H7+/Inbt2/j9u3bAIAWLVpAT08P+vr60NPTY5kxR5T7lVdq06dVHUUKTRIs9+fs8WVlZZXOKM9peU7tv6vNoTedDuFoqrivAAAgAElEQVQ+6tPq09TUxH///cex/fc+pD79pUOHDjAyMoK7uzs+fvyIK1eu4OrVq7h16xYKCgoqXO/Tp08IDQ1lJvpWV1eHsbExOnToAGtr67+eLo6w4/V7VWiSYMuWLSEhIcF21Jednc12dFiudevWHJeXlJREixYteBYrIYLi5uaGuLg4tmuCtZ2+TBTIy8tjypQpmDJlCgoLCxEZGYnQ0FCEhYVVeIRdLjU1lam7umLFCia5GhkZoV+/fmjdujUfXgH5G0KTBBs2bAhdXV2Eh4fDysqKaQ8PD4eFhQXHdfT19XHp0iWWtvDwcPTo0YPjtEuECLvy6cs8PT2RmZkJRUXFahW3JaykpaVhamrKTN79+vVr3LhxA+Hh4YiMjMSXL18qXT85ORnJycnMLVmdO3dmkmL//v3pR3gdIlSjQwMDAzFr1ixs3boVBgYGOHjwII4fP47o6Gioqqpi1qxZAAAfHx8Av36dGRoaYsqUKZg2bRpiY2Ph7OyMAwcOwNLSkicx0mkm7qM+5T7q079XVFSEBw8e4NatW7hz5w7u37+Pnz9/Vnt9MTEx6OjoYMCAARgwYAD69OkDGRkZHkYs3OpUUd264MCBA/Dy8kJWVha0tLSwbt069OvXDwCYWeZ/P/qLjIzEsmXL8Pz5cygqKmLBggWYPn06z+KjLxfuoz7lPupT7ikoKEBMTAwiIiJw/fp1PHv2DCUlJdVeX0pKCoaGhhg0aBBMTU3RuXPnSscsiBpKgkKGvly4j/qU+6hPeSM5ORkKCgqIiYlBZGQkIiIiKh11yknbtm0xbNgwDB06FP3796/0Rn9RQANjCCFEiMjKymLIkCEYMmQIACAvLw9RUVGIiIjAnTt38OTJk0rXf/v2LXx9feHr64smTZpg6NChsLCwgKmpKRo3bsyPlyBSKAkSQggPNWvWDMOHD8fw4cMB/BqhHhERgVu3buHmzZuVTuL99etXnDt3DufOnYOMjAyGDx8Oa2trDB48GA0bNuTXS6jXKAkSQggftWrVCqNHj8bo0aNRVlaGV69e4caNGwgLC0NERAQKCws5rvf9+3cEBgYiMDAQLVq0wJgxYzB16lR07dqVz6+gfhGaadMIIaS+ERMTQ/v27eHg4ICAgAC8fv0aJ0+exD///AMFBYUK1/v06RN8fX3Rv39/mJqa4vjx4/j+/TsfI68/KAkSQkgd0ahRIwwfPhw7duzAs2fPcPXqVcyePRtt2rSpcJ24uDg4Ojqia9euWLduXYXTSBLOKAkSQkgdJC4ujj59+mD9+vVISEjA5cuXMW3aNMjJyXFcPicnB5s2bULXrl0xf/58ZgYbUjlKgoQQUseJi4vD0NAQ27dvx4sXL3DkyBG24gDlfv78iaNHj6J3796YP3++yNWTrClKgoQQIkSkpKRgaWmJoKAgPHz4EPPnz4esrCzbcsXFxTh69Ch69eqFJUuW4PPnzwKItu6jJEgIIUJKQ0MDq1evRmJiIjZs2MBxjtji4mL4+PigR48e2Lt3L9VS/QMlQUIIEXJNmzaFg4MD/vvvP+zduxcaGhpsy+Tm5sLV1RX9+vVDdHS0AKKsmygJEkJEWlpaGuzs7GBubg47OzuhvoYmKSkJW1tb3L9/H3v27EHbtm3ZlklKSsLw4cOxaNEi5OfnCyDKuoWSICFEZKWlpcHKygoBAQGIjIxEQEAArKyshDoRAr+S4YQJE3D//n24u7ujadOmbMscOHAAffr0wc2bNwUQYd1BSZAQIrI8PT1ZChADQEpKCjw9PQUUEXdJS0vDyckJ//33HyZNmsT2/Pv37/G///0PK1asqFE5qPqEkiAhRGRlZGRwbM/MzORzJLzVunVreHt748KFC9DU1GR7fteuXRgwYABevnwpgOgEi+YOrYFv376huLi40mWkpaWRl5fHp4hEA/Up9/GqTxs3bgxJSeH5WlFSUuLYrqioyOdI+MPY2Bh3797FsmXLcOjQIZbnnj59ChMTExw4cICZ7FsUCM+7VcB+/PgB4NeM8JWRkpIS+fpf3EZ9yn286NOysjLk5uaiadOmQpMI3dzcEBcXx3JKVENDA25ubgKMirdkZGTw9etXjs99+/YNEyZMgJubGxYuXCgSxX3pdGg1FRYWolGjRoIOg5A6S0xMDHJycvj27ZugQ6k2NTU1BAUFwdraGkZGRrC2tkZQUBDH++3qk4pOAwO/fsysWbMG06dPR0FBAR+jEgzh+LlWR4jCryJCakMYPyNqamrw9fUVdBh8VdFp4N+dP38e79+/x6lTp9C8eXM+RCUYdCRICCEixs3Nje2GekVFRbbLPbGxsRg5cmSlR47CjpIgIYSIGE6ngUNDQ3H79m106dKFZdmnT59i6NChePXqlYCi5S1KgoSvIiIi0LdvX7Rq1QpjxowRaCxJSUmQk5PD06dPBRpHdRw8eJDj0HZC/lb5aeCQkBD4+vpCTU0N6urquHz5Mvr27cuy7Js3bzB8+HAkJycLKFreoSRYT8nJyVX6N3v2bIHEtXjxYvTu3RuPHj2Cn58f3/ZramrKNuKvXbt2ePHiBTp27Mi3OAip6+Tk5BAYGMh2m8SHDx9gYWHBNrmAsKMkWE+9ePGC+du5cydb24YNGziux+sZ5l+/fg0TExOoqKhUWByUXyQkJKCgoCA0w/kJ4RcZGRkcO3YMEyZMYGnPyMjAqFGj8ObNGwFFxn2UBOspBQUF5q/8YvefbeWnA4OCgjBixAgoKCjg5MmTHE+93bhxA3Jyciz3F929exfDhg2DoqIitLW14eLiUuH9R+X7+vHjB2bOnAk5OTmcO3eO43b/PE158+ZNyMnJISIiAgMGDICSkhIGDx6MxMREln1ERUVhxIgRUFJSgqqqKqysrJCdnY3p06cjLi4O3t7ezJFwVlYWx9Oht2/fxsCBA9G6dWt06tQJK1euZPlhYGpqCldXV6xYsQLq6uro2LEjVq9ejbKyMo6v+9OnT5CXl0d4eDhL++XLl6GgoIDc3FwAwLJly9CzZ08oKipCR0cHa9asqXQaq1WrVrEVVeX073bhwgUYGRlBQUEB3bt3x/r166mUDqkWSUlJeHt7Y8qUKSzt6enpsLCwwPv37wUUGXdREqwFTqcZFRUVqzwV+bd/vLJq1SrMmTMHsbGxMDMzq9Y68fHxsLa2hpWVFaKionDo0CHcv38fCxcu5Lh8+alHCQkJpjr2yJEjaxTnmjVrsHbtWty6dQsyMjKwt7dnnnvw4AGsrKygra2Na9euITQ0FObm5iguLsb27dvRvXt3zJgxgzkSlpeXZ9t+WloabGxs0Lt3b0RGRmLr1q3w9/dnO2r29/eHrKwswsLC4OnpiR07duDixYscY27RogUGDRqEM2fOsLQHBARgyJAhzL+rrKws9u3bh9jYWGzcuBH+/v7MEfzfunz5MhwdHZl/Wy8vL5w+fRobN26s1XaJ6BAXF8eOHTswbtw4lvbU1FTY2Njgy5cvAoqMeygJEsydOxfm5uZQV1ev1v1DALBjxw5MmDABDg4O0NTUhL6+PjZv3owzZ85wLM9SfuoR+PWFr6CgUOMZS1auXIl+/fqhU6dOcHFxQWJiInJycgAA27dvh56eHjZv3oxu3bpBS0sLM2fOZIZ9N2jQADIyMsyRsLg4+1t///790NDQwKZNm9CxY0eYm5tj+fLl2LNnD8vRk46ODlxcXNCuXTvY2NjAwMAAd+7cqTDucePG4eLFi/j+/TsA4MuXL7h69SpsbGyYZZYuXQp9fX2oqalh+PDhmD9/Ps6dO1ej/vnTli1bsGjRItja2kJdXR0DBgzAihUrcODAgVptl4gWcXFx7N69G//73/9Y2p88eYJp06ZVOZVkXUdJkKBHjx41Xic+Ph5Hjx6FiooK82dpaQkAPLtw3rVrV+b/y+d2/PjxIwAgISEBJiYmtdp+UlIS9PX1WW747tOnD75//85SWkdbW5tlPUVFRSYOTkaMGAExMTFcvXoVABASEgIpKSkMHTqUWebs2bMYMmQIOnbsCBUVFXh4eCA9Pf2vX0tZWRkSEhKwbt06ln8jR0dH5ObmMqdhCakOSUlJ+Pj4YNiwYSztN27cwKJFiyq8HCAMaEQAYZsOTlxcnO1N/ed1pNLSUsycORMzZ85k256Kikq1911+RPb7/iq6ZvX7AJbyRFVaWsq2/t8qKyurcMaT39sbNGjA9lxJSUmF25WWloa5uTnOnDmD0aNHMzXrGjZsCACIjIzErFmz4ObmhgEDBkBWVhZBQUHYtGlThdus6t+orKwMpaWlWL58OcfTzpzqyxFSmQYNGsDPzw8jR45EfHw803748GHIycnh3bt3yMjIgJKSEtzc3IRm6jlKgrXA6dd0YWGh0E/23LJlS+Tl5bG8lsePH7Ms0717dzx//rzW9661bNkSAJCVlcV8Mf+5r+ro3r07bt++jcWLF3N8vkGDBpUmKgDo1KkTwsLCWJJhTEwMZGRkoKqqWuOYfjdu3DiMHTsWz58/x507dxASEsI8FxMTAw0NDTg5OTFtVY2+a9myJT58+MDS9nu/iYuLo1u3bnj16hXHf6PCwsK/fSlEhDVu3BinT5+Gqakp3r59y7Tv2LGDZbm4uDihmYOVTocSNgYGBmjYsCE8PDzw+vVrBAYG4ujRoyzLODs7IzIyEosXL0ZCQgJevXqFy5cvY9GiRTXaV+fOnaGgoIB169bh1atXuH79OtsHqjoWLFiAe/fuwcXFBU+ePEFSUhIOHTrE1IVTVVVFXFwc3rx5g5ycHI5Hjvb29khJScGSJUuQlJSES5cuYe3atZgzZw7b0V9NGRkZQV5eHjNnzoSysjLLzcjt27dHWloazp8/j5SUFOzbtw8XLlyodHvGxsbIyMjAzp07kZKSgoMHD+LKlSssyyxZsgTHjx/Hxo0b8ezZM7x48QLnz5/H6tWra/VaiGhTUFDAmTNnICsrW+EywlSYmJIgYdO6dWvs27cPV69ehaGhIU6dOgVXV1eWZXR1dXHp0iUkJSVh+PDhMDY2hqenJ1q3bl2jfUlJScHPzw/Pnz9Hv379sHXrVqxYsaLGMffq1QuBgYFISEjA4MGDYWZmhpCQEOYUqpOTE0pKSmBgYIB27dqxHUUBv2bQOHPmDO7du4f+/fvDyckJEydOxNKlS2scz5/ExcUxduxYPHnyBDY2NiynVy0tLWFvbw9nZ2cYGRkhJiYGS5YsqXR73bp1w8aNG+Hj44P+/fsjJiYGCxYsYFlmxIgROHHiBMLCwjBw4ECYmZlh165daNu2ba1fDxFtWlpaOHLkCMcBZuWEpTCxWG5urvBe0eSjvLy8KmsJAvXjdGhdQ33Kfbzs0+p+Vuqj5ORkdOjQQdBh8M3WrVuxZs0ajs9ZW1tzpToHr/uUjgQJIYT8FScnJ46jsuXl5YWmMDElQUIIIX9FXFwcR44cYTvFXlhYCAkJCQFFVTOUBAkhhPw1OTk5nDhxAlJSUkzbly9fMGfOHOYWprqMkiAhhJBa6datG1atWsXSdufOHezdu1cwAdUAJUFCCCG1NmvWLLbrg6tXr67z9TopCRJCCKk1cXFx7Nmzh2Vk8I8fPzBr1qw6Pb8oJUFCCCFcoaKigm3btrG0PX78uE6fFqUkSAghhGvGjBnDVnFi/fr1LJPQ1yWUBAkhhHDVxo0bWWqgFhQUwMXFpU5WmxCaJPjjxw+4uLhAU1MTysrKGD9+PN69e1fpOuvXr2crTNuxY0c+RVy/zJ49m62wprBYv349y1ydfz7mFTk5OQQHB/N8P7V16dIl9OzZEy1btsTs2bMFHQ6pB+Tl5dnmqL127RqCgoIEFFHFhCYJurq6IiQkBH5+frh8+TK+fPmCcePGVVkZoEOHDkw18RcvXiAqKopPEQve7NmzIScnh82bN7O0R0REQE5OjilIWx0bNmyAj48Pt0MUiHnz5uHSpUtc215FPxBevHjBVn+tLpo/fz4sLCzw+PFjbNiwQdDhkHpi0qRJbD82lyxZUudqWQpFEszLy8OxY8ewevVqDBw4ELq6uvDx8UFiYiJu3bpV6bqSkpJMNXEFBQW0atWKP0HXEdLS0ti5cyeys7NrtZ1mzZqxnN7gtZ8/f/Js202aNEGLFi14tv1yCgoKLDcQ10W5ubnIycnBoEGDoKysLLJzfhLuExcXx44dO1gqsHz48AEbN24UYFTshCIJxsfHo6ioCIMGDWLa2rRpg06dOiE2NrbSdVNTU6GlpQUdHR1Mnz4dqampPI62bjEyMkLbtm0rLdAKAHfv3sXgwYOhoKCADh06wNXVlSUR/Xm0c/fuXZiamkJFRQWqqqoYPHgwnj59im/fvqFt27ZspwHDw8PRqlUrjtUbft/+jh070KVLF3Tp0gXAr2S4Zs0adOnSBcrKyhg4cCDCwsKY9UpKSuDo6AgdHR0oKiqiZ8+e8PLyqnSmit9Ph6alpbGdMpeTk0O3bt2qtf3169fj5MmTCA0NZdaNiIgAwH46NDExEZaWllBUVIS6ujpmz56NvLw8tj7Yu3cvtLS0oKamhjlz5qCgoKDKfq9Ibm4uHBwcoKamBkVFRVhaWuLZs2cAfp0RUFdXBwBYWFiwxE4IN3Tq1Imtuomvry9evnwpoIjYCUVR3Q8fPkBCQoIpwFpOXl6+wi9VAOjduzf27NmDDh06IDs7G5s3b8aQIUMQExNT6ZFAcnIyW5u0tHS1f9XXlYKlJSUlKCsrw7JlyzBt2jRMnz4d6urqTHIrLCxEYWEhMjIyMHbsWFhbW2P79u1ITU2Fs7MzSktL4eHhwWyrpKQEhYWFKC4uxoQJE2Brawtvb28UFRXh8ePHKC4uhoSEBKysrHDkyBEMHTqUieXIkSMwMzODrKwsx/4pKSlBZGQkGjduDH9/fya+OXPmIDU1Fbt374aysjJu3LiB8ePH4+rVq9DW1kZRURHk5eXh4+ODli1b4uHDh3BxcUHTpk0xYcIEAEBxcTFKS0uZ/f7+uFWrVkhISGDi+Pr1K2xsbNC3b18UFhZWuX17e3s8e/YMubm58Pb2BvAr+ZXv6+fPnygsLERBQQHGjBkDXV1dXLlyBZ8/f8aiRYswZ84c+Pn5MX0QFRWFVq1a4fTp03j//j3s7e2hrq6O+fPnV9rvFb3nZs2ahVevXjHVv9evX48xY8bg7t27TCFiExMT+Pn5QU9PjyX2v5Wfn1/p57K+4/T9IcosLCxw7NgxprRScXExnJ2dsWXLlmpvozZ9WlUFCoEmQU9Pzyo74vcK3H/6vQI4J2ZmZiyPe/fuDV1dXZw4cQKOjo4Vrsep0/Ly8qpVeqYulf2RkJCAhIQEzM3NYWBggE2bNuHgwYNo2LAhgF+JXVpaGsePH4eioiJ27NgBcXFx6Ojo4MuXL3BycoK7uzsaNWrEbEtaWhqfP39GXl4ezM3N0blzZwBgjpwAYPr06TA1NcWnT5+grKyM3NxcXL16FYcPH66wb8q3vXfvXubHRkpKCs6fP4/79++jffv2AH79srx79y5OnDiBrVu3QlpaGu7u7sx2OnbsiGfPniE4OBjTp08H8OuUuLi4OLPvPx83btwYAFBaWgpbW1soKipi586dTP9Utn1paWk0btwYP3784Fh9vmHDhpCWlsbp06dRUFAAX19fNG3aFADg5eWFUaNG4f3799DU1ISEhARkZWXh5eUFSUlJ6OjowMrKCnfv3sXixYur7Pc/vXr1CqGhobh06RL69esH4Nev8G7duiEwMBAzZsxAmzZtAPyqIckp/r8hKysrsjULRa2UUnV5enpi5syZzOPbt28jIyMDxsbGVa7L6z4VaBKcPXs2bGxsKl2mTZs2uH//PkpKSpCTk8NyTS87OxuGhobV3l+TJk3QuXNnvH79+q9jFlarV6+Gqakp5s2bx/bcixcvoKenx1Igs2/fvvj58ydev36Nrl27sizfvHlzTJgwAWPGjIGJiQmMjY1hZWXFfKH26NEDXbp0wcmTJ+Hs7IyAgADIycmx/Sj5k5aWFsvR9qNHj1BWVgZjY2OWHzs/fvxg+fAcPHgQR48exdu3b5mjt7/5EnZ3d0diYiLCwsJYkjU3tv/ixQtoa2szCRAADAwMIC4ujufPn0NTUxPAryRfXggYABQVFREXFweg6n7ntE9xcXHo6+szbc2aNUOXLl2QlJRUo/gJqY0xY8bAx8cH9+/fZ9qWLVuG27dvC7zahECvCbZs2RIdO3as9K9Ro0bQ1dVFgwYNEB4ezqz77t07vHjxAgYGBtXeX2FhIZKTk6GgoMCLl1On9ezZExYWFixHNeUqO6KuqH3Pnj24ceMGDA0NceXKFfTu3ZvlWt2UKVOY05rHjx/HhAkTqnyzlx+RlSstLYWYmBiuXr2KiIgI5u/evXvMqcfAwEC4urpiwoQJOHfuHCIiIjBjxowaD6w5ceIEDh06hJMnT7K8P7i1/cruj/q9j38fRFD+3O/rVtXv1d0nIfwkJiaGdevWsbQ9efIEJ06cEFBE/0coBsY0a9YMkydPxsqVK3Hr1i08evQIs2bNgra2NgYMGMAsp6enh/379zOP3dzcEBkZidTUVMTFxWHq1KkoKCiAra2tAF6F4K1cuRLR0dFsX5qdO3fG/fv3WQaTREdHo2HDhtDQ0Khwe926dcOCBQtw6dIl9O/fHydPnmSes7GxQUZGBvbv349Hjx5h4sSJNY5XR0cHZWVl+PDhAzQ1NVn+lJWVmTh79eoFe3t76OrqQlNTEykpKTXaT2xsLJydneHj48N2erE622/YsGGVt+p07twZiYmJ+PLlC8t+S0tL0alTpxrFW1m//7nP0tJS3Lt3j2nLz8/H06dPa7xPQmpLT08PY8eOZWlbu3Ytvn//LqCIfhGKJAgA69atg7m5OaZNm4Zhw4ahcePGOHXqFMvRRXJyMsu9b+/fv8fMmTOhp6eHyZMno2HDhrh+/TrXrn0IG01NTfzzzz/Yt28fS/uMGTOQmZkJZ2dnvHjxAqGhofDw8ICdnR0aNWrEtp3U1FSsWrUKsbGxePPmDe7cuYPExESWL9ZmzZrB0tISbm5uMDQ0RLt27Wocb/v27WFjY4N///0XwcHBSE1NxcOHD7Fr1y5cuHCBWSYhIQHXr1/Hq1evsGnTphrdC5qVlYVJkyZhxowZ6N27N7KyspCVlcXcUlKd7auqquLZs2fM+6+oqIhtP9bW1mjUqBEcHByQmJiIu3fvwsnJCaNGjWJOhValOv3+u3bt2mHEiBFwcnJCVFQUEhMTYW9vj6ZNm2L06NHV7iNCuMXd3Z3lUkNmZiYOHjwowIiEZHQo8GsQx+bNm9lu/P7dnzdhCrpz66LFixezHTkoKysjICAAK1euhJGREZo1a4axY8di5cqVHLfRqFEjvHz5Ev/88w9ycnLQunVrWFtbsw2Fnjx5Mk6dOoXJkyf/dby7d+/Ghg0bsHLlSrx//x7NmzdHz549YWRkBACYNm0aHj9+jJkzZ6KsrAwWFhaYO3cujh8/Xq3tJyUl4ePHj/D29mZOsQJA27Zt8fjx42ptf+rUqYiMjMTAgQPx9etXhISEMPGVa9SoEc6dOwdXV1cMHjwYUlJSGDFiRI1uTq9uv/9uz549WLp0KWxtbfHjxw8YGBjg7NmzkJGRqfZ+CeGWtm3bws7ODrt27WLatm/fjqlTp6JJkyYCiUksNzeXLhxUQ15eXrVuJK5Lo0MFLTAwEAsWLMDz5885HlFWF/Up9/GyT6v7WamPRGV0aFpaGjw9PZGRkQElJSW4ublBTU2tWutmZ2dDV1cXX79+Zdrc3d3h5OTEcXle96nQnA4lwqOgoADPnz/H1q1bMXXq1FolQEJI3ZKWlgYrKysEBAQgMjISAQEBsLKyqnaViFatWsHBwYGlbefOncjPz+dFuFWiJMhHaWlpsLOzg7m5Oezs7OpsaZHa8vLyQv/+/dG8eXO4uLgIOhxCCBd5enqyDQ5LSUmBp6dntbfh6OgIWVlZ5vHnz58FVnOQkiCf1PbXkzBxdXVFdnY2Ll68yPJGJ4QIv4yMDI7t5TPCVIecnBzmzp3L0rZ79258/vy5VrH9DUqCfMKNX0+EECJoSkpKHNsVFRVrtJ3Zs2ejefPmzOP8/HyWW9z4hZIgn3Dj1xO/+fv7Q0VFpcLHdVn5xNgPHz7k+b5cXFwwcuRInu+HkLrAzc2N7f5hDQ0NuLm51Wg7srKy+Pfff1na9u/fzzJhPD9QEuQTbv16EqT//e9/iI+PF3QYhBABUlNTQ1BQEKytrWFkZARra2sEBQVVe3To76ZPn85yySQnJ6fatzdxCyVBPuHWrydBkpGRgby8vKDDEAmlpaVVzkJDiKCoqanB19cXISEh8PX1/asECPw6GpwxYwZLm7e3N4qLi7kRZrVQEuQTbv56qq7Kas+Vn9q8cuUKevXqBQUFBZibm1dab/HP06HldfnOnTsHXV1dtGnTBhMmTGCrWH/8+HEYGBhAQUEBvXr1wu7duyut92dmZobly5eztOXn50NRUZGpKnL69GkMHDgQbdq0Qfv27TF16lS8f/++wm1GRERATk6OJTZOp0yfP38OGxsbZrszZsxAVlYW83xJSQlzT5SamhqWLl1aZbIq3/fVq1fRv39/KCgowMTEhOWourxvr127hr59+0JeXh4vXrxAaWkpNm3aBG1tbbRu3RqGhoa4dOkSy/YzMjJgZ2cHDQ0NKCkpoX///rhz5w7z/JUrV2BiYgIFBQXo6OhgzZo1LHOfXrhwAYaGhkydwxEjRjClkNLT02Frawt1dXUoKSlBT08P586dq/T1ElITDg4OLBPnv3nzBufPn+fb/ikJ8hG3fj1VR3ntuT59+iAyMhI3btyAg4MDyzRzP378wMaNG7F7925cu3YNJSUlmDhxYo0mXn7z5g0CAwNx/PhxBAYGIiEhAWvWrGGeP3LkCNasWYNly5YhNjYWnp6e8PLywoEDByrcpuAKnLEAAA/BSURBVI2NDQIDA1kS5YULFyAtLc3UKPz58ydcXV0RGRmJ06dPIycnh+0XZU1lZmZixIgR0NLSQlhYGIKCgvD161fY2toysXh7e+Po0aPYsWMHrl+/jpKSEgQEBFRr+ytWrICHhwfCw8Ohrq4OGxsblusfhYWF2LJlC7Zv347Y2Fi0bdsWe/fuxa5du7Bq1SpERUVh5MiRmDx5MlMD8du3bxg5ciTevHmD48ePIyoqCosXL2a2GRYWBnt7e9jZ2SEmJgbe3t4IDg5mJjPOysrCjBkzYGtri9jYWFy+fBnjx49n1nd2dsb3798REhKC6OhorF+/XmRvhCe8oaCgwDafs5eXF98mgBeaadNIzXz58gV5eXkYNmwYcxq2Y8eOLMsUFxdjw4YN6NOnDwDAx8cHurq6uH37NsvE5JUpLi7Gnj17mC/Gf/75h6keAQCbN2+Gh4cHLC0tAQDq6upISUmBn58f7O3tOW5zzJgxWLZsGSIiImBiYgIAzC0l5bUQf5+KTV1dHdu2bYO+vj7evXv314N3/Pz80LVrV6aQMPCrT9TV1fHw4UP06tULe/fuxfz585m5Nzdu3IibN29Wa/suLi4YPHgwgF/Dwbt06YKzZ89iypQpAH4dZW7atAm6urrMOt7e3nB0dIS1tTUAYPny5YiKioK3tzf279+Ps2fP4sOHD7h+/TpTdPr30+5btmzBvHnzMGnSJOa5VatWwd7eHuvXr0dGRgaKiopgaWnJzKnbpUsXZv23b9/CwsKCmVi8vBI9Idw0b948HDlyhEl8T548QVhYGExNTXm+bzoSrKd+rz1nY2MDb29vpKensywjLi6OXr16MY9VVVWhpKSE58+fV3s/bdu2ZTkyUFRUZCafzs7ORnp6OpycnKCiosL8eXh4MLeL/PkcALRo0QKDBg3CmTNnAPw6WomIiGCpPRkfHw9bW1t07doVbdq0wcCBAwGA7TXWxKNHjxAVFcUSj7a2NoBft7Pk5eUhMzMTenp6zDp/9mFlfq/r16RJE2hra7P0taSkJEsVi/z8fGRkZDA/Usr17duXWS8hIQHa2tpMAuT0mrZu3crymuzs7FBQUICsrCx069YNAwYMgKGhISZPngw/Pz/m3w/4dapqy5YtMDMzg6enJw2MIjzRrl07WFhYsLR5eXnxZd90JFiP7dmzB7Nnz0ZYWBiuXLkCT09P+Pv7M0cj3MCp/l35qcPy/27btq3Cuo/Lli3jWOh33LhxWLBgAbZu3Yrz589DRUUFffv2BfDrFOCYMWMwYMAA+Pj4QF5eHjk5ORg+fHiFdf7KCwb/forlz4vvpaWlGDJkCMd7N+Xl5Su9jskNUlJS1S4wWl6DsKpTRqWlpViyZAmsrKxY2n/8+IFWrVpBQkIC58+fx/3793Hz5k0cO3YMHh4euHTpErp164YpU6Zg8ODBuH79Om7duoUhQ4bAyckJrq6uf/ciCanAggULEBwczDyOiIjAs2fPWIpM8wIdCdZzldWeKy0txYMHD5jHb9++RUZGBtdqzbVu3RrKyspISUlhqwdYXj5IXl6erQ0ARowYAQAIDQ1FYGAgbGxsmC/+8pJFK1asQL9+/dCxY0d8/Pix0lhatWoFgPW+zMePH7Ms0717dzx//hxt27Zli7Vp06Zo1qwZS6V34FcS+r0PK/N7Ve1v375VWddPVlYWSkpKiImJYWmPjo5m1uvevTsSExPZBiP9/pqSkpLYXo+Ghgbz5SImJgZ9fX0sXboU4eHhUFJSYhmYoKKign/++QeHDx/GsmXLcOTIkWq9XkJqokePHjA0NGRp8/X15fl+KQnWU9WpPScpKQlXV1fcu3cPCQkJmD17Njp37lzt64HVsXTpUuzcuRO7d+9GcnIynj59ipMnT2Lbtm2VrictLQ1zc3Ns3rwZCQkJLKdC27RpAykpKfj6+iI1NRWhoaFsVav/pKmpiTZt2mDDhg14+fIlbt68yVaWa+bMmcjPz8e0adMQFxeH1NRU3Lp1C//++y9TDNfBwQFeXl4IDg5GcnIyli5dyjJ6tDJbtmxBeHg4nj17BkdHRzRs2JCtyOif5s2bB29vb5w9exYvX77E2rVrER0dDUdHRwDA2LFj0apVK0ycOBFRUVFITU3F5cuXmdGhixcvxtmzZ7F27Vo8ffoUSUlJCA4OxurVqwH8SsybN2/GgwcP8PbtW1y+fBnv3r1j3idLlizBjRs3kJqaioSEBNy4cYMK8hKe+XOcwKlTp1gKUfMCnQ6tp6pTe05KSgrOzs5wcHBAeno6evfujePHjzNHXNwwZcoUNGrUCDt37sTq1ashLS0NLS0t2NnZVbnuuHHjcOLECejo6LB88bZq1Qp79+7F6tWrceDAAWhra2Pt2rUYM2ZMhdtq0KAB/Pz84OzsjP79+6Nbt25YuXIlxo0bxyyjpKTEFBQeM2YMfvz4wVxvLB/C7ejoiKysLOYU7rhx42BtbY0X/6+9ewtp8v/jAP6eWJ77zQvLM4VZHkpo0oJKVk6C0FIMMSk64EWBN+XZskTRm2aF4oWGCUajtBDKLItOlKzEKIWSSgtKPJI625zhab+LaP+/Rf2s9vS0Pe/X3b4sn7cP4Xv77nn2efXqP3+fgoICHD16FN3d3QgJCUFdXR3c3Nx++G8OHjwIo9GIgoICDA0NITg4GOfOnUNERAQAwM3NDU1NTcjPz8fOnTsxNTWF5cuXW14UqNVq1NfXQ6PRoKKiAo6OjggKCrJcaLNo0SK0trbizJkzGBsbg5+fH7KysiznZXZ2FtnZ2ejt7YW7uztUKhW/6o8EExsbC19fX8vtTiaTCY2NjVAoFIIdk/ME58ne5glqtVrLH7e/na2c0+95+PAhtm3bhjdv3nz3ApY/jfMEhSGVeYLz9StzBzUaDUpKSiyPAwIC0NHRYflc39q4HUpERFb3q5Nz9u7dO+eCu56eHty5c0ewnCxBIiKyul+dnLN48WLLfbhfCHmBDEtQonbt2mUTW6H2ICoqCnq9/q/ZCiX6E35ncs7XF8i8fv0a4+PjVsn1NV4YQ0REVvc7k3MiIyOhUCjwzz//IC4uDvv27Zv3PbQ/iyVIRERWl5+fjydPnszZEp3v5ByZTIampia4uLigq6tLsAIEWIJERCSAL5NziouLMTAwAG9v73ldHfqFi4uLwAk/YwnOk4ODAyYnJy1f4ExEc5nNZphMJsG/5opsx5fJOX8z/m+dJ3d3dxiNRkxMTPzweR8/fpwzKZl+H8+p9Ql1Tp2dnefMhiP627EE50kmk8HDw+M/nzc0NISAgIA/kEg6eE6tj+eU6DPeIkFERJLFEiQiIsliCRIRkWSxBImISLI4RYKIiCSL7wSJiEiyWIJERCRZLEEiIpIsliAREUkWS5CIiCSLJSgws9mMHTt2QC6X48qVK2LHsVmjo6PIysrC2rVr4e3tjfDwcKSnp2NkZETsaDanuroaERERWLJkCVQqFXQ6ndiRbNapU6ewefNmBAQEICgoCMnJyejs7BQ7lt05efIk5HI5srKyrP6zWYICq6ioEHQWllT09/ejv78fhYWF0Ol0qKqqgk6nQ2pqqtjRbEpDQwNyc3ORkZGBBw8eQKlUIikpCT09PWJHs0ktLS1ITU3FzZs3cfXqVTg6OiIhIQGjo6NiR7MbbW1tqK2tRXh4uCA/n/cJCujZs2fYvXs37t+/j+DgYNTW1iI+Pl7sWHbj1q1bSE5Oxrt37zhlYp7UajXCw8NRXl5uWVMoFIiPj0dBQYGIyeyD0WhEYGAgtFottm7dKnYcmzc2NgaVSoWysjKcOHECYWFh0Gg0Vj0G3wkKxGAwIDU1FadPn4aXl5fYceySwWCAk5MTXF1dxY5iEyYnJ9He3o7o6Og569HR0WhtbRUplX0xGo2YnZ2FXC4XO4pdOHToEOLj46FSqQQ7BkcpCSQ9PR1qtRpbtmwRO4pd0uv1KCkpwZ49ezjEdZ6Gh4cxMzPzzYsyLy8vDA0NiZTKvuTm5mL16tVQKpViR7F5tbW1ePv2LaqqqgQ9Dv96/ITi4mKUlpb+8DmNjY3o7e3F8+fPce/evT+UzHbN95xGRUVZHo+PjyMlJQU+Pj4oKioSOqLdkclkcx6bzeZv1ujnHTlyBI8fP0ZzczOvA/hNXV1dKCoqwo0bN7Bw4UJBj8XPBH/C8PAwhoeHf/gcf39/ZGRk4OLFi3Bw+N9u88zMDBwcHKBUKtHc3Cx0VJsx33P6ZcvTaDQiKSkJAHDp0iW4u7sLntFeTE5OwsfHB2fPnkVCQoJlPTMzE52dnbh+/bqI6WxbXl4eGhoa0NjYiBUrVogdx+ZptVqkpaXNeTExMzMDmUwGBwcH9PX1wcnJySrHYgkKoK+vD3q9fs7a+vXrUVJSgtjYWCxdulScYDbOYDAgKSkJZrMZly9fhoeHh9iRbI5arcaqVatQVlZmWYuMjMT27dt5YcwvysnJQUNDA65du4aVK1eKHccu6PV69PX1zVlLS0tDUFAQ0tPTERoaarXdC26HCsDX1xe+vr7frPv7+7MAf5HBYEBiYiIMBgO0Wi1MJhNMJhMAwNPTU/AtE3uRlpaGAwcOIDIyEuvWrUNNTQ0GBgawf/9+saPZpMzMTNTV1eH8+fOQy+UYHBwEALi5uXGX4jfI5fJvLi5ydXWFp6cnwsLCrHosliDZhPb2drS1tQH4/M7l/339mSF9X2JiIkZGRqDRaDA4OIjQ0FDU19cjMDBQ7Gg2qbq6GgC+ufUpJycHeXl5YkSin8TtUCIikizeJ0hERJLFEiQiIsliCRIRkWSxBImISLJYgkREJFksQSIikiyWIBERSRZLkIiIJIslSEREksUSJCIiyWIJEknAxMQElEolFAoFxsfHLevj4+NYs2YNlEolPn36JGJCInGwBIkkwMXFBZWVlXj//j2OHz9uWT927Bh6enpQWVkJZ2dnERMSiYNTJIgkQqFQ4PDhw9BoNIiNjQUA1NTUIDs7GwqFQuR0ROLgFAkiCZmamkJMTAw+fPgAs9kMLy8v3L59GwsWLBA7GpEoWIJEEvPixQts2LABjo6OaGlpQUhIiNiRiETDzwSJJObu3bsAgOnpabx69UrkNETi4jtBIgl5+fIlVCoV4uLi0Nvbi+7ubjx69AheXl5iRyMSBUuQSCKmp6cRExODwcFB6HQ66PV6bNy4EZs2bYJWqxU7HpEouB1KJBGlpaVob29HWVkZPD09sWzZMhQWFqKpqQkXLlwQOx6RKPhOkEgCOjo6EBMTg5SUFJSXl1vWzWYzEhMT8fTpU+h0Ovj5+YmYkujPYwkSEZFkcTuUiIgkiyVIRESSxRIkIiLJYgkSEZFksQSJiEiyWIJERCRZLEEiIpIsliAREUkWS5CIiCSLJUhERJL1L8McEa21od4kAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "true_mean   = spline(true_coef,domain)\n",
    "y_noiseless = spline(true_coef,x)\n",
    "\n",
    "y_noisy = y_noiseless + np.random.randn(x.shape[0])*0.2\n",
    "\n",
    "plt.style.use('fivethirtyeight')\n",
    "plt.scatter(x,y_noisy,label = 'Noisy realizations of \\nspline-valued process',color='k')\n",
    "plt.plot(domain,true_mean,color='k',label = 'True function value')\n",
    "plt.legend(loc = 'lower left')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('f(x)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Estimating spline parameters\n",
    "\n",
    "Now, I would like to use the same technique on my data in order to estimate spline parameters for creating a representative function. My data shows a linear trend first and then changes into a quadratic dependency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fa2e6f0e630>]"
      ]
     },
     "execution_count": 78,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEJCAYAAADbzlMFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGzZJREFUeJzt3X9wVNX9//HXEhKy8sONYdmoCfhVIgQaBdFAcRQMCKUWUq2pUMdpqQoFKupEB2JrENoxYMRKFKgK6GihrQt0RPvD0TH4Awj8UZwwonFHSkwUkiY0SjAxMbnfPzT7YZPdkE32x9nk+ZjZmebcs3vfZ4/ltffu2Xtt9fX1lgAAiLIB0S4AAACJQAIAGIJAAgAYgUACABiBQAIAGIFAAgAYgUACABiBQAIAGKHPBpLH44l2Cb3GGKIv1uuXYn8MsV6/xBi6q88GEgAgthBIAAAjEEgAACMQSAAAIxBIAAAjDIx2AQAAsw2eMkWTPvrI+3fr2LE6U1oa8v1whAQACGjwlCmK++gj2STvI+6jjzR4ypSQ74tAAgAE1B5GZ2sPpVAjkAAARiCQAABGIJAAAAG1jh0rq0Ob9V17qBFIAICAzpSWekOp/RGuVXYs+wYAdOlMaak8Ho/S09PDuh+OkAAARiCQAABGIJAAAEYgkAAARiCQAABGIJAAAEYgkAAARiCQAABGIJAAAEYgkAAARiCQAABGIJAAAEYgkAAARiCQAABGIJAAAEYgkAAARiCQAABGIJAAAEYgkAAARiCQAABGIJAAAEboViDt27dP8+fPV0ZGhhwOh7Zv3+6zfcmSJXI4HD6PmTNn+vT5+uuv9eCDD+rSSy/VRRddpPnz5+uzzz4L3UgAADGtW4F05swZjRs3TmvXrpXdbvfbZ/r06SovL/c+3G63z/b8/Hy9+uqr2rp1q/7xj3/o9OnTuu2229Ta2tr7UQAAYt7A7nSaNWuWZs2aJUlaunSp3z6DBg2Sy+Xyu+2LL77QSy+9pI0bN+qGG26QJD3zzDPKzMzU3r17NWPGjJ7UDgDoQ0L2HdKBAwc0evRoTZo0ScuXL9d///tf77b3339fLS0tys7O9ralpqZqzJgxOnjwYKhKAADEsG4dIZ3LzJkzNXfuXI0aNUqffvqpfv/732vevHnau3evBg0apJqaGsXFxSk5OdnneU6nUzU1NQFf1+Px9Kqu3j7fBIwh+mK9fin2xxDr9UuMQZLS09O73B6SQPrJT37i/d/jx4/XhAkTlJmZqddff13z5s0L+DzLsmSz2QJuP1fxXfF4PL16vgkYQ/TFev1S7I8h1uuXGEN3hWXZ94UXXqiLLrpIx44dkySNGDFCra2tqqur8+lXW1srp9MZjhIAADEmLIFUV1enEydOeBc5TJgwQfHx8SopKfH2+eyzz1ReXq7JkyeHowQAQIzp1im7hoYG79FOW1ubqqqqVFZWpqSkJCUlJWnt2rWaN2+eXC6XPv30U61Zs0ZOp1M/+tGPJEnnn3++7rjjDhUUFMjpdCopKUm/+c1vNH78eE2fPj1sgwMAxI5uBdLhw4c1d+5c79+FhYUqLCzUggUL9MQTT+jo0aP6y1/+oi+++EIul0vXXXednn/+eQ0dOtT7nEcffVRxcXFauHChmpqadP311+uPf/yj4uLiQj8qAEDM6VYgXXfddaqvrw+4fffu3ed8jcTERBUVFamoqKj71QEA+g2uZQcAMAKBBAAwAoEEADACgQQAMAKBBAAwAoEEADACgQQAMAKBBAAwAoEEADACgQQAMAKBBAAwAoEEADACgQQAMAKBBAAwAoEEADACgQQAMAKBBAAwAoEEADACgQQAMAKBBAAwAoEEADACgQQAMAKBBAAwAoEEADACgQQAMAKBBAAwAoEEADACgQQAMAKBBAAwAoEEADACgQQAMAKBBAAwAoEEADACgQQAMAKBBAAwAoEEADACgQQAMAKBBAAwAoEEADACgQQAMAKBBAAwAoEEADACgQQAMAKBBABhFO92K3PuXA1LStLQzEzFu93RLslYA6NdAAD0VfFut+zLl8vW2ChJslVWyr58uSSpJTc3mqUZiSMkAAiTxDVrvGHUztbYqMQ1a6JUkdm6FUj79u3T/PnzlZGRIYfDoe3bt/tstyxLhYWFGjt2rFJSUnTTTTfpww8/9OlTX1+vRYsWaeTIkRo5cqQWLVqk+vr60I0EAAxjq6oKqr2/61YgnTlzRuPGjdPatWtlt9s7bd+wYYM2btyodevW6a233pLT6dTNN9+s06dPe/vcddddKisrk9vt1s6dO1VWVqbFixeHbiQAYBgrNTWo9v6uW4E0a9YsFRQUKCcnRwMG+D7Fsixt3rxZ9913n3JycjRu3Dht3rxZDQ0N2rlzpySpvLxcb775pp588klNnjxZWVlZ+sMf/qDXX39dHo8n9KMCAAM0FRTI6vAh3rLb1VRQEKWKzNbr75AqKipUXV2t7Oxsb5vdbtfUqVN18OBBSdKhQ4c0ZMgQTZ482dtnypQpGjx4sLcPAPQ1Lbm5aiwu1tcpKbJsNrWlpamxuJgFDQH0epVddXW1JMnpdPq0O51OnThxQpJUU1Oj5ORk2Ww273abzabhw4erpqYm4Gv39uipLxx9MYboi/X6pdgfQ0zXP2GC9Oqrvm0xOp7ezkN6enqX20O27PvssJG+PZXXMYA66tino3MV3xWPx9Or55uAMURfrNcvxf4YYr1+iTF0V69P2blcLknqdKRTW1vrPWoaMWKEamtrZVmWd7tlWaqrq+t0ZAUA6J96HUijRo2Sy+VSSUmJt62pqUkHDhzwfmeUlZWlhoYGHTp0yNvn0KFDOnPmjM/3SgCA/qtbp+waGhp07NgxSVJbW5uqqqpUVlampKQkpaWlacmSJVq/fr3S09M1evRoPf744xo8eLBuvfVWSdKYMWM0c+ZM3X///dqwYYMsy9L999+v2bNnx/xhLAAgNLoVSIcPH9bcuXO9fxcWFqqwsFALFizQ5s2bde+996qxsVEPPvig6uvrNWnSJO3evVtDhw71Pue5557TihUrdMstt0iS5syZo8ceeyzEwwEAxKpuBdJ1113X5VUVbDab8vPzlZ+fH7BPUlKSnn322eArBAD0C1zLDgBgBAIJAGAEAgkAYAQCCQBgBAIJAGAEAgkAYAQCCQBgBAIJAGAEAgkAYAQCCQBgBAIJAGAEAgkAYAQCCQBgBAIJAGAEAgkAYAQCCQBgBAIJAGAEAgkAYAQCCQBgBAIJAGAEAgkAYAQCCQBgBAIJAGAEAgkAYAQCCQBgBAIJAGAEAgkAYAQCCQBgBAIJAGAEAgkAYAQCCQBgBAIJQJ8S73ZraGamhiUlaWhmpuLd7miXhG4aGO0CACBU4t1u2Zcvl62xUZJkq6yUfflySVJLbm40S0M3cIQEoM9IXLPGG0btbI2NSlyzJkoVIRgEEoA+w1ZVFVQ7zEIgAegzrNTUoNphFgIJQJ/RVFAgy273abPsdjUVFESpIgSDQALQZ7Tk5qqxuFhtaWmybDa1paWpsbiYBQ0xglV2APqUltxcAihGcYQEADACgQQAMAKBBAAwAoEEADACgQQAMAKBBAAwQkgCqbCwUA6Hw+dx+eWXe7dblqXCwkKNHTtWKSkpuummm/Thhx+GYtcAgD4iZEdI6enpKi8v9z7279/v3bZhwwZt3LhR69at01tvvSWn06mbb75Zp0+fDtXuAQAxLmSBNHDgQLlcLu9j+PDhkr49Otq8ebPuu+8+5eTkaNy4cdq8ebMaGhq0c+fOUO0eABDjQhZIx48fV0ZGhq644gr98pe/1PHjxyVJFRUVqq6uVnZ2trev3W7X1KlTdfDgwVDtHgAQ42z19fVWb1/kjTfeUENDg9LT01VbW6uioiJ5PB6VlpbK4/Fo9uzZOnLkiNLS0rzPWbZsmU6cOKHdu3cHfF2Px9Pb0gAAhkhPT+9ye0iuZXfjjTf6/H311VdrwoQJ2rFjh6655hpJks1m8+ljWVanto7OVXxXPB5Pr55vAsYQfbFevxT7Y4j1+iXG0F1hWfY9ZMgQjR07VseOHZPL5ZIk1dTU+PSpra2V0+kMx+4BADEoLIHU1NQkj8cjl8ulUaNGyeVyqaSkxGf7gQMHNHny5HDsHgAQg0Jyyu63v/2tfvCDHyg1NdX7HdJXX32lBQsWyGazacmSJVq/fr3S09M1evRoPf744xo8eLBuvfXWUOweANAHhCSQPv/8c911112qq6vT8OHDdfXVV+uNN97QyJEjJUn33nuvGhsb9eCDD6q+vl6TJk3S7t27NXTo0FDsHgDQB4QkkLZt29bldpvNpvz8fOXn54didwCAPohr2QEwSmJenoYlJ2uYw6FhyclKW7s22iUhQggkAMZIzMtTwtatsrW2yibJ1tqqEbt2KTEvL9qlIQIIJADGSHjhBXX8daLtu3b0fQQSAHO0tgbXjj6FQAJghHi3O/DGuLjIFYKoIZAAGCFxzZpOp+skyZLU/ItfRLgaRAOBBMAItqqqgNua1q+PYCWIFgIJgBGs1FS/7c0pKRGuBNFCIAEwQlNBgSy73afNstv12dKlUaoIkUYgATBCS26uGouL1ZaWJstmU1tamhqLi3Vqzpxol4YICcmlgwAgFFpyc9WSm+vbyI06+w2OkAAARiCQAABGIJAAAEYgkAAARiCQAABGIJAAAEYgkAAARiCQAABGIJAAAEYgkAAARiCQAABGIJAAAEYgkAAARiCQAABGIJAAAEYgkAAARiCQAABGIJAAAEYgkAAARiCQAHjFu90ampmpYUlJGpqZqXi3O9oloR8ZGO0CAJgh3u2Wffly2RobJUm2ykrZly+XJLXk5kazNPQTHCEBkCQlrlnjDaN2tsZGJa5ZE6WK0N8QSAAkSbaqqqDagVAjkABIkqzU1KDagVAjkABIkpoKCmTZ7T5tlt2upoKCKFWE/oZAAiDp24ULjcXFaktLk2WzqS0tTY3FxSxoQMSwyg7oY+Ld7m8XKFRVyUpN1QV33y2lp3fruS25uQQQooYjJKAPaV+6PaCyUjbL0oDKSo169FF+T4SYQCABfYi/pdtxTU0s3UZMIJCAPoSl24hlBBLQh7B0G7GMQAL6EH9Lt1sTE1m6jZjAKjugD2lfIXf2KruKu+9WMivnEAMIJKCP6bh0+5THo+Qo1gN0V8RP2W3ZskVXXHGFXC6Xpk2bpv3790e6BACAgSIaSLt379bKlSuVl5end955R1lZWcrNzVVlZWXY9hmp+7sk5uVpWHKyhjkcGpacrMS8vLDsBwD6qogG0saNG/Wzn/1MP//5zzVmzBgVFRXJ5XJp27ZtIdtHewBNysrS0EsvlX3ZMp8fCdoXLdIwhyOk4ZSYl6eErVtla22VTZKttVUJW7eGJJS4YRqA/iJigdTc3Kz3339f2dnZPu3Z2dk6ePBgSPbR6Vfqp07J1tzs08dmWbJJ34bT8uU9/gf+7KBI2LpVtg7bbZISXnihR6/d7oJ//rPTr+57UzMAmCxigVRXV6fW1lY5nU6fdqfTqZqampDsw9+v1LvS05uPdQy+jmHk1doa9Guf7eJNm7hhGoB+I+Kr7Gw233++Lcvq1NbO4/EE9dqTevBrdFtVVdD7yXz44e4F34ABQb/22SZVV/tt70nN0RRLtfoT6fov+Oc/dfGmTUqorlazy6XPli7VqTlzevWazEH0MQYp/RwX+Y1YICUnJysuLq7T0VBtbW2no6Z25yq+Iys1VbYgF0hYqalB7ychQFD4vK6k5oULg37tszW7XBp08mTn1+5BzYGcl5OjgW+/7f37m2nT9NUrr4TktaVv/wMOVa3REOn6491u2QsLvR94Bp08qf9XWKiUCy/s8VW4mYPoYwzdE7FTdgkJCZowYYJKSkp82ktKSjR58uSQ7MPvDcbi49V2wQWyJFkdj856ePOxgJdnaX/Exan5zjvVtH590K99ts+WLg3rDdPaw8gmeR8D335b5+XksJgiSvydduY0LfqLiK6yW7ZsmXbs2KEXX3xR5eXlWrFihU6ePKmFCxeG5PX93mBs0yadPnZMX9bXq/HZZ0Ny87FAd9ZsfO45fVlfry/r6nodRpJ0as6csN4wrT2MztYeSr1dTHFeTo6GORyadM01GuZw6LycnJDU3FcECnwujor+LKLfId1yyy06deqUioqKVF1drYyMDL388ssaOXJkyPbR/it1f4eXobr5mL/LszQVFITlxmbRumFaoE/p3anl7COvdu1HXqE8HRir2hfFtL/Htu8CXwp82pmLo6I/iPiVGu666y4dOXJENTU1evvtt3XttddGuoSQaMnN1ekjR/Tl//6n00eO9Iu7bHb3U3pXR16xKpSnMLs6LRfo6JuLo6I/4Grf/dmgQUF1N/VTeqCwCFWI+LsLa29+D9bVaTm/p51DeJoWMBmB1I81Pv20rAG+/wlYAwao+c47Y+ZTeqCwSMzLC1mIhHqhwbnuWdQfj74BiUDq11pyc9X4zDO+n8afeUZN69f36lP6N9OmyerQZn3XHmqBwiLhhRdCFiKhXmjAaTnAP24/0c8FWjTRm8UUX73ySth/39QuYCgEuEpGd0Ik3u32LljJdLlkJSXJdupUp349PYUZyUUxQCwhkBAW7eET7h/TBfwxdFyc31A6V4h0XAE36ORJWfHxshISfK6L2NsjmmitngRMxik7xLRAp7+af/GLHp0W83sKsKVF1pAhLDQAwowjJMS0rk5/tU6ZEvRpsYDfF/3vf/ry2LGQ1w/g/xBIiHmh/B6MH6YC0cMpO+AsrIADoodAAs7S8YepX6ek8H0RECGcsgM6OPtUX1+4bQAQKzhCAgAYgUACABiBQAIAGIFAAgAYwVZfX9/xOpgAAEQcR0gAACMQSAAAIxBIAAAjEEgAACMQSAAAI8RsIG3ZskVXXHGFXC6Xpk2bpv3793fZ/7333tO0adPkcrl05ZVXatu2bRGqNLBgxvDuu+/K4XB0enz88ccRrPj/7Nu3T/Pnz1dGRoYcDoe2b99+zud88MEH+uEPf6iUlBRlZGRo3bp1sqzoLfIMdgwVFRV+5+DNN9+MUMW+nnjiCd1www1KS0vTZZddpttuu01Hjx495/NMmYee1G/aHDz33HOaOnWq0tLSlJaWphtvvFGvv/56l88x5f1vF+wYwjkHMXktu927d2vlypVav369pkyZoi1btig3N1elpaVKS0vr1P/48eP66U9/qttvv13PPvusSktLlZeXp+TkZOXk5ERhBMGPoV1paamSkpK8fw8fPjwS5XZy5swZjRs3TgsWLNCvfvWrc/b/8ssvdfPNN2vq1Kl666235PF4tGzZMp133nm65557IlBxZ8GOod2uXbv0ve99z/v32fMRSe+9957uvPNOXXXVVbIsS48++qh+/OMf6+DBgwFrMmkeelJ/O1Pm4KKLLtLq1at12WWXqa2tTX/+8591++23a+/evT71tTPp/W8X7BjahWMOYvJ3SDNmzND48eNVXFzsbbvqqquUk5OjVatWdeq/atUqvfrqq/r3v//tbbvnnnv00Ucf6Y033ohIzR0FO4Z3331Xc+fO1SeffKLk5ORIlnpOF198sR577DHdfvvtAfts3bpVjzzyiD7++GPZv7u9Q1FRkbZt26ajR4/KZrNFqly/ujOGiooKXXnllSopKdHEiRMjWF33NDQ0aOTIkdq+fbvmzJnjt4/J89Cd+k2fA0m65JJLtGrVKi1cuLDTNpPf/7N1NYZwzkHMnbJrbm7W+++/r+zsbJ/27OxsHTx40O9zDh061Kn/jBkzdPjwYbW0tISt1kB6MoZ206dP15gxYzRv3jy988474SwzpA4dOqTvf//73v8TSt/OwYkTJ1RRURHFyoJ3xx13aPTo0Zo9e7ZeeeWVaJfj1dDQoLa2NjkcjoB9TJ6H7tTfzsQ5aG1t1a5du3TmzBllZWX57WPy+y91bwztwjEHMRdIdXV1am1tldPp9Gl3Op2qqanx+5yamhq//b/55hvV1dWFrdZAejKGlJQUPfHEE3rppZf00ksvKT09XTk5Odq3b18kSu61QHPQvi0WDBkyRL/73e/0/PPPy+126/rrr9fChQv117/+NdqlSZJWrlypzMzMLv8hMXkeulO/iXPwwQcf6OKLL9aIESN0//33609/+pPGjx/vt6+p738wYwjnHMTkd0iSOh3aWpbV5eGuv/7+2iMpmDGkp6f73JcnKytLn376qZ566ilde+21Ya0zVEycg2AkJyf7nOefOHGiTp06pQ0bNui2226LYmXSQw89pNLSUv3rX/9SXFxcl31NnIfu1m/iHKSnp+vdd9/VF198oT179mjJkiV67bXXNG7cOL/9TXz/gxlDOOcg5o6QkpOTFRcX1+nTRG1tbadPHu1GjBjht//AgQN1wQUXhK3WQHoyBn8mTZqkY8eOhbq8sAg0B5KCGrNpTJiD/Px87dq1S3v27NEll1zSZV8T5yGY+v2J9hwkJCTo0ksv1cSJE7Vq1SplZmZq06ZNfvua+P5LwY3Bn1DNQcwFUkJCgiZMmKCSkhKf9pKSEk2ePNnvc7KysrR3795O/SdOnKj4+PhwlRpQT8bgz5EjR+RyuUJdXlhkZWXpwIEDampq8raVlJTowgsv1KhRo6JYWe9Eew5WrFihnTt3as+ePbr88svP2d+0eQi2fn+iPQcdtbW1qbm52e82097/QLoagz+hmoOYCyRJWrZsmXbs2KEXX3xR5eXlWrFihU6ePOldEbJ48WItXrzY23/hwoX6/PPPtXLlSpWXl+vFF1/Ujh079Otf/zpaQwh6DJs2bdJrr72mTz75RB9++KFWr16tv//977r77rujUn9DQ4PKyspUVlamtrY2VVVVqaysTJWVlZKk1atXa968ed7+t956q+x2u5YuXaqjR49qz549evLJJ7V06dKonaoIdgw7duyQ2+1WeXm5PB6PnnrqKW3ZskWLFi2KSv0PPPCAduzYoS1btsjhcKi6ulrV1dVqaGjw9jF5HnpSv2lz8Mgjj2j//v2qqKjQBx98oNWrV+u9995Tbm6u3/pNev97OoZwzkFMfod0yy236NSpUyoqKlJ1dbUyMjL08ssva+TIkZKkqqoqn/6XXHKJXn75ZT300EPatm2bUlJStG7duqj9BkkKfgwtLS16+OGHdeLECSUmJnr7z5o1Kxrl6/Dhw5o7d67378LCQhUWFmrBggXavHmzTp48qf/85z/e7eeff77+9re/6YEHHtANN9wgh8OhZcuWRfVDQbBjkKTHH39clZWViouL02WXXaann346at9dbNmyRZI6/Xe8YsUK5efnS5LR89CT+iWz5qC6ulqLFi1STU2Nhg0bpvHjx2vnzp2aMWOGJLPf/3bBjkEK3xzE5O+QAAB9T0yesgMA9D0EEgDACAQSAMAIBBIAwAgEEgDACAQSAMAIBBIAwAgEEgDACAQSAMAI/x/UiHV10ne4nAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.array([0.        , 0.03364583, 0.06396991, 0.26822917, 0.30600694,\n",
    "       0.30619213, 0.53136574, 0.7740625 , 0.77430556, 0.7750463 ,\n",
    "       0.77590278, 0.8390625 , 1.00467593, 1.0059375 , 1.00693287,\n",
    "       1.2787963 , 1.34362269, 1.48939815, 1.53609954, 1.74604167,\n",
    "       1.81266204, 1.97753472, 2.14038194, 2.46284722, 2.47528935,\n",
    "       2.47613426, 2.63649306, 3.01434028, 3.51010417, 3.51114583,\n",
    "       3.51288194, 3.5137037 ])\n",
    "y_noisy = np.array([-1.79115574e-01,  6.25958906e-01, -6.47491364e-01, -1.29351477e+00,\n",
    "        2.64397486e-01,  3.63929988e-01, -9.37538928e-02, -3.98027000e+00,\n",
    "       -4.06193312e+00, -4.17972094e+00, -4.84284664e+00, -4.53880276e+00,\n",
    "       -5.52900380e+00, -5.52900380e+00, -5.52900380e+00, -8.22763816e+00,\n",
    "       -8.05190714e+00, -1.06254842e+01, -6.94957746e+00, -5.19495296e+00,\n",
    "       -3.17329480e+00,  8.30268444e+00,  1.48384712e+01,  3.85102710e+01,\n",
    "        4.15606490e+01,  4.15606490e+01,  6.35088300e+01,  1.03625013e+02,\n",
    "        1.77265790e+02,  1.77265790e+02,  1.77265790e+02,  1.77265790e+02])\n",
    "\n",
    "plt.plot(x, y_noisy, 'ro')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_bins = 2\n",
    "degree = 4\n",
    "num_coef = n_bins + degree + 1 # Restraint that must be obeyed by parameter sets of B-splines\n",
    "domain   = np.linspace(x.min(),x.max(),20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.75685185]\n"
     ]
    }
   ],
   "source": [
    "breaks = np.histogram(domain, n_bins)[1][1:-1]\n",
    "print(breaks)\n",
    "for i in range(degree+1):\n",
    "    breaks = np.insert(breaks, 0, domain.min()-1e-6)\n",
    "    breaks = np.append(breaks, domain.max()+1e-6)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This break seems quite a good estimate, since the function is linear for smaller values and quadratic afterwards. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [],
   "source": [
    "with pm.Model() as model:   \n",
    "    coef = pm.Flat('coef',shape = num_coef,testval = np.zeros(num_coef))\n",
    "    x_as_tensor = tt.as_tensor(x)\n",
    "    s = spline_fn_expr(breaks, coef, degree, x_as_tensor)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (2 chains in 4 jobs)\n",
      "NUTS: [sigma, coef]\n",
      "Sampling 2 chains: 100%|██████████| 5200/5200 [01:13<00:00, 70.29draws/s] \n",
      "The chain reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.\n",
      "The chain reached the maximum tree depth. Increase max_treedepth, increase target_accept or reparameterize.\n",
      "The gelman-rubin statistic is larger than 1.05 for some parameters. This indicates slight problems during sampling.\n",
      "The estimated number of effective samples is smaller than 200 for some parameters.\n"
     ]
    }
   ],
   "source": [
    "SAMPLES = 2000\n",
    "BURN = 600 \n",
    "\n",
    "with model:\n",
    "    sigma  = pm.HalfCauchy('sigma',beta=2.0)\n",
    "    y_hat  = pm.Deterministic('y_hat',s)\n",
    "    y      = pm.Normal('y',mu = y_hat,sd = sigma,observed = y_noisy)\n",
    "    #step = pm.Metropolis()\n",
    "    trace  = pm.sample(SAMPLES, chains=2, tune=BURN)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apparently, the sampling process faced some problems. I already tried to use a greater sample size, but I didn't change much in the outcoming fit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-2.17213140e+00,  7.46408966e+00, -2.54286865e+01, -1.72642752e+01,\n",
       "        1.06338317e+02,  1.76751366e+02, -2.55739008e+99])"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "T = len(domain)\n",
    "n_samples = 500\n",
    "function_samples= np.zeros([T,n_samples])\n",
    "for i in range(n_samples):\n",
    "    coef_sample = trace['coef'][i,:]\n",
    "    function_samples[:,i] = spline(coef_sample[:],domain)\n",
    "    \n",
    "ci = np.percentile(function_samples,[95,50,5],axis = 1)\n",
    "#ci = np.percentile(function_samples,[65,45,5],axis = 1)\n",
    "trace['coef'][i,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fa2e686cf28>]"
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaQAAAEJCAYAAADbzlMFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XdYVNe+PvCXooBiQBEBFVBglGpDsV2OCMYaEEWueowFRUSKguUqkUQ5YrBEjZ1EgomJmmM/ip4kFjxGBQ0qwWMhxIKNInjB0ARhfn+cH3ODMyBlyh54P8/jk7DXGubLsryz9157LY2CggIxiIiIVExT1QUQEREBDCQiIhIIBhIREQkCA4mIiASBgURERILAQCIiIkFgIBERkSAwkIiISBBaRCBlZGSouoQmYf2qpc71q3PtAOtXNWXX3yICiYiIhI+BREREgsBAIiIiQWAgERGRIDCQiIhIEBhIREQkCAwkIiISBAYSERFJiY2NxcCBA2FoaAgjIyPExsYq/D21Ff4ORESkVmJjY7F8+XLJ15WVlZKvAwMDFfa+PEMiIqIaVqxY0aDj8sJAIiKiGiorKxt0XF4YSEREVIOWllaDjssLA4mIiGpYs2ZNg47LCwOJiIhqCAwMxNq1a6Gp+Z+I0NLSwtq1axU6oQHgLDsiIpIhMDAQHh4eEIlESntPniEREZEgMJCIiEgQGEhERCQIDCQiIhIEBhIREQkCA4mIiASBgURERILAQCIiIkFgIBERkSAwkIiISBAYSEREJAgMJCIiEgQGEhERCQIDiYiIBIGBREREgsBAIiIiQWAgERGRIDCQiIhIEBhIREQkCAwkIiISBAYSEREJQr0C6fLly5gyZQrs7OxgaGiIffv21WifP38+DA0Na/waMWJEjT6vX7/G0qVLYWVlhc6dO2PKlCl49uyZ/H4SIiJSa/UKpOLiYtjb22Pt2rXQ09OT2cfNzQ3p6emSX4cOHarRHhERgZMnT+Krr77C6dOn8ccff2Dy5MmorKxs+k9BRERqT7s+nUaOHImRI0cCAIKCgmT20dHRgYmJicy2wsJCfPvtt9ixYweGDx8OAPjiiy/g5OSECxcuwMPDozG1ExFRMyK3e0hJSUmwsbGBs7MzFixYgBcvXkjaUlNTUVFRAXd3d8mxrl27omfPnrh69aq8SiAiIjVWrzOkdxkxYgQ8PT1haWmJx48fIzo6Gl5eXrhw4QJ0dHSQm5sLLS0tGBkZ1XidsbExcnNza/2+GRkZ8ihP7t9LFVi/aqlz/epcO8D6VU2e9YtEojrb5RJIPj4+kv93cHBAnz594OTkhB9//BFeXl61vk4sFkNDQ6PW9ncVX18ZGRly+16qwPpVS53rV+faAdavSlVVVcjIyEDHjh2lTiYURSHTvs3MzNC5c2c8ePAAANCpUydUVlYiPz+/Rr+8vDwYGxsrogQiImqk9PR0uLq64sSJE7h//77S3lchgZSfn4+srCzJJIc+ffqgVatWSExMlPR59uwZ0tPTMXDgQEWUQEREDSQWi/Hdd9/Bzc0Nt2/fxvbt2/HHH38o7f3rdcmuqKhIcrZTVVWFp0+fIi0tDe3bt0f79u2xdu1aeHl5wcTEBI8fP8bf/vY3GBsb44MPPgAAGBgYYPr06fjkk09gbGyM9u3bY8WKFXBwcICbm5vCfjgiIqqfV69eISwsDEePHpUcKywsxLp16zB8+HBoaip+HYV6BdLNmzfh6ekp+TomJgYxMTGYOnUqNm3ahDt37uD7779HYWEhTExM4Orqij179qBdu3aS13z66afQ0tKCn58fysrK8Je//AWxsbHQ0tKS/09FRET1duPGDcyaNQuPHz+Wamvbti1ev35d6zOo8lSvQHJ1dUVBQUGt7X9O1Nro6upiw4YN2LBhQ/2rIyIihamqqsL27dsRFRUltUiBjo4OPvzwQ8yePVspYQTIaZYdERGplxcvXiAgIKDGvf1q3bt3x5IlS9CtWzfY29srrSYGEhFRC3PhwgX4+/sjLy9Pqm3cuHHw9/dH37598eLFizofzZE3BhIRUQtRUVGBNWvWYMuWLRCLxTXa9PX1sWDBAowcORIODg5o1apVjRV3lIGBRETUAmRmZsLPzw83btyQarO3t8fixYvh7OwMS0tLpZ4V/RkDiYiomTt+/DhCQkJQVFRU47iGhgYmT56MadOmoVevXjAwMFBRhf/BQCIiaqZKS0uxbNky7N27V6rNyMgIixYtwtChQ+Ho6AhtbdXHgeorICIiubt79y5mzJghc3HUAQMGYMGCBXByclLpJbq3MZCIiJoRsViMr7/+GsuXL8fr169rtGlra2PWrFmYMGECHB0dVX6J7m0MJCKiZqKgoADBwcE4deqUVFvnzp2xdOlS9OvXDw4ODoK4RPc24VVEREQNlpSUhNmzZyMrK0uqbfjw4Zg3bx7s7e1hbm4umEt0b2MgERGpscrKSqxfvx4bNmxAVVVVjTZdXV3Mnz9f8mzRe++9p6Iq64eBRESkpp49ewY/Pz9cu3ZNqs3a2hpLliyBg4MD7O3tBXmJ7m3Cr5CIiKScPHkSwcHBePXqlVTbhAkTMH36dNjY2Aj6Et3bGEhERGqktLQUy5cvxzfffCPVZmBggLCwMAwePBiOjo7Q19dXQYWNx0AiIlITd+7cwYwZM/D7779LtfXt2xdhYWEQiUTo2bOnWu41x0AiIhI4sViMr776CitWrJB6tkhLSwszZ86Et7c3evToAVNTUxVV2XQMJCIiAfvf//1fzJs3Dz/99JNUm5mZGZYsWYLevXvD0dERurq6KqhQfhhIREQCdfnyZcyePRs5OTlSbdXPFllbW8PGxkZtJi7UhYFERCQwlZWV+PTTT7F582apZ4v09PQQFBQEd3d32NnZoUOHDiqqUv4YSEREAvLkyRPMmjUL169fl2oTiURYsmQJevbsCXt7e7Ru3VoFFSoOA4mISCCOHTuGBQsW4I8//pBq8/HxwV//+lfY2NjAwsKiWVyiexsDiYhIxYqLi7F48WJ8//33Um2GhoYIDw/HwIED4eDggHbt2qmgQuVgIBERqVBqaipmzpyJzMxMqTZnZ2csWLAANjY2sLW1VctnixqCgUREpAJVVVXYvHkzYmJi8ObNmxpt1fsWeXp6okePHjAzM1NRlcrFQCIiUrLnz59j9uzZSE5Olmrr0qULlixZAkdHRzg6OkJPT08FFaoGA4mISIn+8Y9/IDQ0VOaiqKNGjcKcOXPQrVs32NjYQFNTUwUVqg4DiYhICUpKSrB48WIcOHBAqk1fXx+hoaFwdXWFra1ts3q2qCEYSERECpaamopZs2bh0aNHUm29e/dGWFgYrK2t1WbfIkVpuT85EZGCVU9cWLNmjdTEBS0tLcyYMQPjx49Hjx490LlzZxVVKRwMJCIiBcjKysL8+fNx48YNqbbOnTtj8eLFzWZRVHlhIBERydnJkycREhKCwsJCqbaRI0dizpw5sLGxgZWVVbNccaGxGEhERHJSUlKCpUuXYt++fVJt+vr6CAkJwbBhw2Bvbw8DAwMVVChsDCQiIjlIS0vDzJkz8fDhQ6k2JycnhIWFwdbWtkWsuNBYDCQioiaoqqrCli1b8Omnn6KioqJGm5aWFqZNmwYfHx/Y2trCxMRERVWqBwYSEVEjPXnyBHPmzMG1a9ek2szMzBAQEIBhw4bBwcEBOjo6KqhQvTCQiIga4cCBA1iyZAmKi4ul2kaMGAF/f3+0a9cOffv25cSFemIgERE1QEFBAUJCQpCQkCDVpq+vj6CgIHh4eMDBwQFZWVkMowZgIBER1VNiYiICAgLw4sULqbbevXtj4cKFcHBwgEgkanHr0MkDA4mI6B1KS0vx8ccfIy4uTqqtVatWmDVrFsaPHw87OzsYGRmpoMLmgYFERFSHtLQ0zJo1Cw8ePJBqs7KyQnh4OPr06QM7O7sWvQ6dPHD0iIhkqKysxObNm7F27Vqpdeg0NDTg4+ODadOmwdbWtsVsoKdoDCQiordkZmZizpw5SElJkWrr1KkTFi1aBBcXlzqnc2dmZuLjjz9GUVERzMzMEBkZCUtLS0WXrtYYSERE/59YLMa+ffvwP//zPygpKZFqHzFiBObOnQs7OzuYm5vXOoMuMzMT3t7eNVZtSElJwfHjxxlKdWAgEREBePnyJYKDg/HPf/5Tqu29996rMZ27TZs2dX6v6OhoqSWEHj58iOjoaOzevVuudTcn9ZqXePnyZUyZMgV2dnYwNDSUWjhQLBYjJiYGtra2MDU1xbhx43D37t0afQoKChAQEAALCwtYWFggICAABQUF8vtJiIga6cyZM3BxcZEZRs7Ozti6dSt8fX3Rv3//d4YR8J+tJ2TJzs5ucq3NWb0Cqbi4GPb29li7di309PSk2rds2YIdO3Zg3bp1OH/+PIyNjTFhwgT88ccfkj7+/v5IS0vDoUOHcPjwYaSlpWHevHny+0mIiBqouLgYYWFh8PX1RV5eXo02HR0dBAYGYs2aNXB3d4e1tXW9H3KtbZKDqalpk2tuzup1yW7kyJEYOXIkACAoKKhGm1gsxq5duxAWFobx48cDAHbt2gWRSITDhw/Dz88P6enpOHv2LH744QcMHDgQALB582aMGTMGGRkZEIlE8vyZiIjeKSkpCfPmzcPjx4+l2kQiEcLDw+Hs7AwbG5sGr84dGRmJlJSUGpftunfvjsjIyCbX3Zw1+VHizMxM5OTkwN3dXXJMT08PQ4YMwdWrVwEA165dg76+viSMAGDQoEFo27atpA8RkTKUlZUhIiICY8eOlQojTU1NTJkyBRs3bsTo0aPRs2fPRm0VYWlpiePHj2P06NFwdXWFr68vJzTUQ5MnNeTk5AAAjI2Naxw3NjaWXEfNzc2FkZFRjdNdDQ0NdOzYEbm5ubV+74yMjKaWp5DvpQqsX7XUuX51rh2Qb/23b9/GJ598IvOsyMTEBAEBAejTpw86deqE/Px85OfnN+n9Vq9eLfn/8vJytfy9kGfN77oaJrdZdm9fWxWLxVIB9La3+7xNXpfy1P2yIOtXLXWuX51rB+RXf3l5OWJiYrB161ZUVlZKtY8bNw5+fn5wdHSU630ejn/DNDmQqjecys3NRdeuXSXH8/LyJGdNnTp1Ql5eXo0AEovFyM/PlzqzIiKSp7S0NAQEBODevXtSbZ06dUJISAhcXV1hb2/PPYtUrMn3kCwtLWFiYoLExETJsbKyMiQlJUnuGbm4uKCoqKjGJlbXrl1DcXFxjftKRETyUlFRgXXr1mH48OEyw2jkyJHYvn07Jk2ahD59+jCMBKBeZ0hFRUWShQWrqqrw9OlTpKWloX379jA3N8f8+fOxceNGiEQi2NjY4LPPPkPbtm0xadIkAEDPnj0xYsQIhIeHY8uWLRCLxQgPD8eoUaPU+nSWiITp7t27CAgIwK1bt6TaOnTogJCQELi7u8POzo5BJCD1CqSbN2/C09NT8nVMTAxiYmIwdepU7Nq1CwsXLkRpaSmWLl2KgoICODs74+jRo2jXrp3kNbt378ayZcswceJEAMCYMWOwfv16Of84RNSSVVZWYtu2bVizZg0qKiqk2ocPH4558+ahV69eMDMz4+Z5AlOvQHJ1da1zVQUNDQ1EREQgIiKi1j7t27fHl19+2fAKiYjq4ffff8e8efNw/fp1qTYDAwMEBQVh1KhRPCsSMK5lR0RqraqqCl988QWioqJQVlYm1T506FAEBwejb9++MDU15VmRgDGQiEhtPXr0CPPnz0dSUpJUW7t27TBv3jx4enrCzs4OrVu3VkGF1BAMJCJSO1VVVdizZw8iIyNRWloq1e7i4oLQ0FAMGDCA68epEQYSEamVBw8eIDg4WOZZUdu2beHv7w8fHx/Y2tryrEjNMJCISC1UVlZi586dWLNmjcx7Rf369cPChQsxcOBAnhWpKQYSEQne3bt3ERwcjBs3bki16erqYvbs2ZgyZQrPitQcA4mIBKuiogLr16/H+vXr8ebNG6n2fv36ITQ0FEOGDJEsY0bqi4FERIKUmpoKf39//P7771Jt+vr6mDNnDiZNmgQ7Oztoa/OfsuaAv4tEJCilpaVYu3Yttm/fLnNl7sGDByMkJASDBg2CkZGRCiokRWEgEZFgJCUlISQkBPfv35dqMzQ0REBAACZOnAiRSNSojfNI2BhIRKRyRUVFiIqKQlxcHMRisVS7m5sbQkJCMHDgwBprZFLzwkAiIpVKTExEaGgonj59KtXWvn17hIaGwtvbG927d+eyP80cA4mIVKKgoAArVqzAvn37ZLaPHj0akydPxujRo6Gnp6fk6kgVGEhEpHQJCQlYvHgxcnJypNpMTU2xYMECeHl5oaSkhGHUgjCQiEhpsrKysGzZMpw4cUKqTUNDA15eXggKCpLs4JqRkaGCKklVGEhEpHCVlZWIj49HVFQUioqKpNotLCwQFhYGT09PGBsbq6BCEgIGEhEp1K1bt7Bw4UKZy/5oaWnBx8cHwcHBcHBw4AOuLRx/94lIIYqLi7Fu3Trs2LFD5gOuIpEI4eHhGDt2LAwNDVVQIQkNA4mI5O6nn37C4sWL8eTJE6k2PT09TJ8+Hf7+/rCxsYGmpqYKKiQhYiARkdxkZWUhIiICx48fl9k+ePBghIeHw9XVlbPnSAoDiYiarLKyEnv27EFUVBT++OMPqfaOHTsiODgY06ZNg7GxMR9wJZkYSETUJLdu3UJ4eDhSUlKk2jQ1NeHp6Ynw8HA4OTlx/TmqEwOJiBrlXZMWbGxssHjxYnh5eaFt27YqqJDUDQOJiBrszJkzWLx4MR4/fizVpqenhxkzZiAkJARdu3bl5TmqNwYSEdVbdnY2IiIicOzYMZntgwYNQkREBIYOHcpniqjB+CeGiN6poqICu3fvRkxMTK2TFkJDQzF79mxuD0GNxkAiojr961//wrJly3Dv3j2pNk1NTXh5eSEiIgI9evTg5TlqEgYSEcn0+PFjREZGylwIFQCsra3x0UcfYfz48bw8R3LBP0VEVENpaSm2bt2Kzz//HKWlpVLtbdq0wcyZM7FkyRIYGRmpoEJqrhhIRAQAEIvFSEhIwIoVK2TOngOA999/Hx999BH69OnDy3MkdwwkIsJvv/2G5cuX4/z58zLbRSIRli9fzstzpFD8k0XUgr169Qrr169HbGws3rx5I9VuYGCAuXPnIiwsDPr6+iqosP4yMzMRHR2NrKwsmJmZITIyEpaWlqouixqAgUTUAlVVVeH777/HqlWrkJubK9VePXvuk08+gZWVlQoqbJjMzEx4e3vj4cOHkmMpKSk4fvw4Q0mNcN13ohYmNTUVo0ePRlBQkMww6tWrFw4ePIj4+Hi1CCMAiI6OrhFGAPDw4UNER0erqCJqDJ4hEbUQeXl5WL16Nfbu3QuxWCzV3rFjR4SFhSEgIACtW7dWQYWNl5WVJfN4dna2kiuhpmAgETVz5eXliI+PR0xMDAoLC6XaW7VqhcmTJ2PlypUwNjZWQYVNZ2ZmJvO4qampkiuhpmAgETVTYrEYiYmJ+OKLL3D//n2ZfQYPHozVq1fD2dlZradxR0ZGIiUlpcZlu+7duyMyMlKFVVFDMZCImqEbN25gxYoVSEpKktnetWtXREREYOrUqc1iC3FLS0scP34c0dHRyM7OhqmpKWfZqSEGElEzkpmZidWrV+Pw4cMy2/X09DB79mxEREQIfhp3Q1laWmL37t2qLoOagIFE1AwUFBRg8+bNiI2NxevXr6XaNTU1MXbsWERFRcHa2loFFRK9GwOJSI1VVFRgz549WLt2LV6+fCmzj4uLC/72t79h0KBBSq6OqGEYSERqSCwW49SpU1i5cmWtExasrKwwd+5cBAYGqvWEBWo5GEhEaubmzZtYsWIFrly5IrO9Q4cOWLBgAYKCgpCZmckwIrXBQCJSE48fP0Z0dDQOHjwos11XVxczZszAihUrYGBgoOTqiJpOLvM9Y2JiYGhoWONXjx49JO1isRgxMTGwtbWFqakpxo0bh7t378rjrYmavcLCQqxatQoDBgyQGUaamprw9PTE1atXsX79eoYRqS25nSGJRCIkJCRIvtbS0pL8/5YtW7Bjxw7s2LEDIpEI69evx4QJE/DLL7+gXbt28iqBqFkpKytDfHw8Nm7ciPz8fJl9XFxcsGbNGgwYMEDJ1RHJn9wCSVtbGyYmJlLHxWIxdu3ahbCwMIwfPx4AsGvXLohEIhw+fBh+fn7yKoGoWSgvL8e+ffuwYcMGPH/+XGaf7t2745NPPoG3tzfvEVGzIbdHtB89egQ7Ozv06tULs2fPxqNHjwD850G9nJwcuLu7S/rq6elhyJAhuHr1qrzenkjtVVZW4sCBAxgwYADCw8NlhpGRkRFWr16NX375BRMmTGAYUbOiUVBQIL3sbwOdOXMGRUVFEIlEyMvLw4YNG5CRkYHk5GRkZGRg1KhRuHXrFszNzSWvCQ4ORlZWFo4ePVrr983IyGhqaUSCV1VVhfPnz+OLL76QfJB7m46ODnx8fBAQEIC2bdsqt0AiORGJRHW2y+WS3fvvv1/j6/79+6NPnz7Yv3+/5Nr225/kxGLxOz/dvav4+srIyJDb91IF1q9aiqpfLBbjxx9/xJo1a3Dr1i2ZfVq3bo1Jkybh448/rnVF67pw7FWL9TeMQlZV1NfXh62tLR48eCC5r/T2RmB5eXlqu9Q9UVP961//wsiRIzFlyhSZYaStrY2JEyfi6tWr2LlzZ6PCiEjdKCSQysrKkJGRARMTE1haWsLExASJiYk12pOSkjBw4EBFvD2RYF29ehWenp4YP348fvnlF6n26jXnrly5gvj4eHTv3l0FVRKphlwu2UVGRmL06NHo2rWr5B5SSUkJpk6dCg0NDcyfPx8bN26ESCSCjY0NPvvsM7Rt2xaTJk2Sx9sTCV5qairWrFmDM2fO1NrH3d0dK1euRO/evZVYGZFwyCWQnj9/Dn9/f+Tn56Njx47o378/zpw5AwsLCwDAwoULUVpaiqVLl6KgoADOzs44evQon0GiZu/u3buIiYnBiRMnau0zZMgQrFq1Ci4uLkqsjEh45BJI8fHxdbZraGggIiICERER8ng7IsG7d+8eNm3ahEOHDkEslj2RtX///oiMjISbm5tyiyMSKPXfKpJIQFJTUzF9+nQMHjwYBw8elBlGjo6OOHDgAM6cOcMwesvly5fRq1cvWFhYoFevXrh+/bqqSyIl4uKqRHKQlJSEjRs34uzZs7X26dGjB5YtW4YJEyY0i23D5e3y5csYP3483rx5AwB49eoVgoODYW5ujqFDh6q4OlIGBhJRI4nFYiQmJuKzzz6rdSsIAOjWrRsWLVqEDz/8kEFUh/nz50vCqFplZSXmz5+PtLQ0FVVFysRAImqgqqoqnDp1Cps2bcLNmzdr7WdjY4Pg4GBMnz4d2tr8q/Yute14W1hYqORKSFX4t4Sont68eYOjR49i8+bNdW6f4uDggAULFsDX15dnRPWUmZmJ0tJSmW3cTqPlYCARvUN5eTm++eYbfP7553j48GGt/ZydnbFo0SKMHTuWi542UHR0NCorK2W27dq1S8nVkKowkIhqUVJSgm+++QabN2+WWvrqz4YOHYolS5bAzc2NQdRIWVlZMo/b2NhwQkMLwkAiekthYSG++uor7Ny5E3l5eTL7aGpqYvjw4Vi2bBkfaJWD2tbqs7GxUXIlpEoMJKL/79GjR4iNjcV3332HoqIimX20tbUxatQoLF++HE5OTkqusPmKjIxESkpKjUui3bt3R2BgoAqrImVjIFGLJhaLceXKFezcuROnT5+udVWF1q1bw8vLC8uWLVPr7QSEytLSEsePH0d0dDSys7NhamqKyMhIlJeXq7o0UiIGErVI5eXlOHbsGHbu3Ilff/211n56enoYM2YMoqKiamwwSfJnaWmJ3bt31zjGTTpbFgYStSgvX77Enj17sHv3bmRnZ9faz9DQEJMnT8aiRYvw6tUrhhGREjCQqEVIT09HbGwsvv/++1qfdwEAKysrzJo1C3PmzJFsFf7q1StllUnUojGQqNmqXtpn586dda4xp6GhgUGDBmH+/Pn44IMP+DArkYowkKjZKS0txaFDh7Br1646V1TQ1dXF2LFjsXDhQm6KRyQADCRqNp49e4Y9e/Zgz549yM/Pr7WfsbExJk+ejJCQEJiamiqxQiKqCwOJ1NqbN29w5swZfP311zhz5gyqqqpq7Wtraws/Pz/MnDkTurq6SqySiOqDgURq6fHjx/j222/x3Xff1brsDPCfFRVcXV0RFBSEkSNHcmkfIgFjIJHaqKiowA8//IC9e/fi7NmztT7ECgBt27aFp6cnFi5cCDs7OyVWSUSNxUAiwXv06JHkbCgnJ6fOvlZWVvDx8UFgYCCMjIyUVCERyQMDiQSpvLwc//znP/H1118jMTGxzr66urpwd3eHn58f3N3doaWlpaQqiUieGEgkKA8ePMA333yD/fv348WLF3X2tbGxwaRJk+Dn5wcTExMlVUhEisJAIpUrKSnB6dOnsXfvXly8eLHOvnp6ehgxYgRmzZqF4cOH8yFWomaEgUQq8ebNG1y8eBEHDx5EQkJCrds9VOvRowd8fX3h5+eHjh07KqlKIlImBhIpjVgsRmpqKg4ePIgjR47UuQsrALRp0wbvv/8+Zs+eDVdXV54NETVzDCRSuKdPn+LYsWM4dOhQvbYTsLe3x6RJkzBz5kzOlCNqQRhIpBB5eXmSELp27do7+3fo0AEeHh6YMWMGhg4dyrMhohaIgURyUz054dChQzh37hzevHlTZ/82bdrA1dUVEydOhJeXF/T09JRUKREJEQOJmqR6csLf//53nDp16p2TE7S1tTFgwAB4enpi8uTJvCRHRBIMJGqw0tJSJCYmIiEhAT/88ANevnz5ztc4ODhg7Nix+PDDD2FhYcE15QQqMzMT0dHRyMrKgpmZGSIjI2FpaanqsqiFYCBRvRQUFOCnn35CQkICzp07h+Li4ne+xsLCAqNGjYKbmxvGjBnD+0ICl5mZCW9vbzx8+FByLCUlBcePH2cokVIwkKhW2dnZOH36NBISEnDx4sV33hMCACMjI3h4eMDX1xdubm5o1aoVMjLWRFaGAAAQVElEQVQyGEZqIDo6ukYYAcDDhw8RHR2N3bt3q6gqakkYSFTD/fv3kZCQgFOnTuGXX36pc0Xtau3atcOQIUPg7e0NLy8vtG3bVgmVkrzVto1Hdna2kiuhloqB1MKJxWL8+uuvkhCqa8vvPzMyMsJ//dd/YfTo0Rg7diwMDAwUXCkpmpmZmczj3FWXlIWB1AKVlJTgypUrOHv2LE6dOoUnT57U63Vdu3aFq6srxo0bBw8PD07TbmYiIyORkpJS47Jd9+7dERkZqcKqqCVhILUAYrEY//73v5GYmIjz588jKSkJr1+/rtdrRSIRhg0bBk9PTwwZMgStWrVScLWkKpaWljh+/Diio6ORnZ0NU1NTzrIjpWIgNVO5ubmSAEpMTHznunHVNDU14eTkBDc3N3h7e8PJyQna2vxjoi7enrY9bdo0iESier/e0tKSExhIZfgvTTNRVlaGq1ev4ty5czh//jz+/e9/1/u1rVu3Rr9+/eDu7g5vb29YW1tzkzs1JGvadlJSEk6dOsWzHFILDCQ1JRaLkZ6ejnPnziExMRGXL19GaWlpvV/fqVMn9O/fH66urvD09ETnzp05NVvNyZq2/fTpU07bJrXBQFITYrEYGRkZSE5OxpUrV3Dx4kU8f/683q/X1dVF7969MXjwYIwYMQL9+vVDmzZtFFgxKRunbZO6YyAJVEVFBdLS0pCUlISzZ8/i1q1byM/Pr/frNTQ0YGNjg4EDB8LV1RUeHh4wMjLikj3NGKdtk7pjIAlEUVERUlJSkJSUhKSkJKSkpKCkpKRB36Njx47o378/hgwZgtGjR8PKyooTEloQWdO2u3btymnbpDb4r5WK5ObmIjk5GUlJSUhOTkZaWhoqKysb9D10dHTg5OSEwYMHw8PDAy4uLrwM14LJmrY9bdo0TmggtcFAUoKysjLcvXsXqampuH79OpKTk/H77783+Pvo6enBzs4Offr0gYuLC9zd3WFsbMzLcCTx9rTt+uzQSyQUSg+kuLg4bN26FTk5ObC1tUVMTAyGDBmi7DIUpqSkBLdv30Zqaip+/fVXpKam4t69e/VamPRt7du3h6OjI2xsbDBq1CgMGjQIBgYGDCAiapaUGkhHjx7F8uXLsXHjRgwaNAhxcXHw9fVFcnIyzM3NFfreitjnpaioCLdu3ZIET1paGu7du4eqqqpGfb8uXbrAyckJ/fr1w7Bhw9C7d2/o6uoiIyOjQQ83EhGpI6UG0o4dO/DXv/4VM2fOBABs2LAB586dQ3x8PFauXCn396sOoTt37uDRo0c19vA5ffo07OzsJGt1vSucCgsLkZaWhl9//VXyKyMjo16rYcuiqakJGxsb9OrVCwMGDICbmxusra1rnYTAjdOIqLlTWiCVl5cjNTUVoaGhNY67u7vj6tWrcn8/WU+t/1lxcTFSUlIkv6o3IcvPz8e9e/fw22+/Sf6bnp7eoGd+ZOnUqRNEIhHs7Ozg7OwMNzc3mJiY1Oth1GfPniE8PJwbpxFRs6a0QMrPz0dlZSWMjY1rHDc2Nq73OmsNIeup9do8fPgQHh4eAIC8vDy516Kvr487d+40egp2bGwsN04jomZP6ZMa3r4hLxaLa71J35QZQg8ePGhQ/8YEkYaGBnR0dFBWVlZnv3bt2tU7HGV58eKFzOMPHjxQm1lU6lJnbZRV/7NnzxAbG4sXL17A2NgYgYGB6NKlS5O+J8detVj//3nXvXClBZKRkRG0tLSkzoby8vKkzpqqNeVGvpWVFa5fv97o179NS0sLXbp0Qc+ePSWX3QYOHAh/f39cunSp1tdpa2sjLi6uST9LbeNjZWUll8kOR44cQWhoKF6/fg0dHR1s27YNPj4+Tf6+1dR9Uoay6s/MzJS6NJuent6kS7Mce9Vi/Q2jtEBq3bo1+vTpg8TERHh7e0uOJyYmwsvLS+7vJ+up9fpo1aoVzM3NYWlpCWtra9ja2sLJyQn29vbQ19eXOpurbbkWbW1tdO7cGbt27cLQoUMb/XMAQGBgINLT0xWycdqRI0cwZ84cydclJSWSr/v378+JFAoia5KKrMvMvDRLLYlSL9kFBwdj3rx5krOL+Ph4ZGdnw8/PT+7v9een1n/66ScUFhbWaNfS0kLr1q1haGgIDw8PDB48GE5OTujZsyd0dHTq/T617bIpzwkHXbp0UdjGaUFBQTKPBwYGomvXrk2aSFF95lVWVgZdXV25n3mpK1kTblJSUmBkZCSzPxdHpZZCqYE0ceJEvHz5Ehs2bEBOTg7s7Oxw8OBBWFhYKOT9qp9aj46ORnp6OkQiERwdHdG7d290795dLtstKGuXTUVtnFZeXi7zeEVFRZM+rdd15qWOoSTPafe1nQnVtnQUF0ellkLpkxr8/f3h7++v1PecPHmyQq+DqvMumxoaGg16lqq+n9bfnt7/5+PKDCR5BIm8p93Xtk1Ep06doKWlpZBLs0TqgDuytXADBgyQeby2y0f1/bT++vXrBh1viszMTMydOxcffPAB5s6di8zMTMlxb29vHDp0CJcuXcKhQ4fg7e0taa+vuqbdN0Zt9x2rL/X6+vrC1dUVvr6+fNaMWhQurtrCffnllxgxYkSNqeXGxsb4+uuvERIS0uhP6zo6OjK3z2jI/bn6qO1+TPVlVHlMEqht2n1j7+3Udt+x+uxNXc+2iZqKZ0gtnKWlJc6ePVvjU/nZs2cxdOjQJn1a37ZtW4OON1ZdodOUHVT/fNZV2yodjb23U33fkWdCRDXxDIlq/VTelE/r1feJFD3Lrq7QaewOqrLOurS1tWus2N7Uezs8EyKSxkAihfHx8YGPj49CH66rK3TqujRWF1lnXW/evIGFhQUsLS0VNpOSqKVjIJFae9f9mMZMya/trMvS0hInT56Ua/1E9H8YSKTW3hU6jbk01thLfUTUNAwkUnvyvh/T2Et9RNQ0DCSit7x91tW2bVusW7eO94yIFIyBRCTDn8+6MjIyGEZESsDnkIiISBAYSEREJAgMJCIiEgQGEhERCQIDiYiIBIGBREREgsBAIiIiQdAoKCio/3ahRERECsIzJCIiEgQGEhERCQIDiYiIBIGBREREgsBAIiIiQWgWgRQXF4devXrBxMQEw4YNw5UrV+rsf+nSJQwbNgwmJibo3bs34uPjlVSpbA2p/+eff4ahoaHUr99++02JFf/H5cuXMWXKFNjZ2cHQ0BD79u1752tu376NsWPHwtTUFHZ2dli3bh3EYtVM9Gxo/ZmZmTLH/uzZs0qq+P9s2rQJw4cPh7m5OaytrTF58mTcuXPnna8Tyvg3pn4hjf/u3bsxZMgQmJubw9zcHO+//z5+/PHHOl8jlLEHGl6/ssZe7befOHr0KJYvX46NGzdi0KBBiIuLg6+vL5KTk2Fubi7V/9GjR/jv//5vTJs2DV9++SWSk5OxePFiGBkZYfz48YKvv1pycjLat28v+bpjx47KKLeG4uJi2NvbY+rUqQgMDHxn/1evXmHChAkYMmQIzp8/j4yMDAQHB6NNmzYIDQ1VQsU1NbT+akeOHIGjo6Pk6z//PijLpUuXMGfOHPTr1w9isRiffvopvL29cfXq1VrrEdL4N6b+akIY/86dOyMqKgrW1taoqqrCgQMHMG3aNFy4cKFGbdWENPaNqb+aosde7Z9D8vDwgIODA7Zu3So51q9fP4wfPx4rV66U6r9y5UqcPHkSN27ckBwLDQ3FvXv3cObMGaXU/GcNrf/nn3+Gp6cn7t+/DyMjI2WWWqcuXbpg/fr1mDZtWq19vvrqK6xatQq//fYb9PT0AAAbNmxAfHw87ty5Aw0NDWWVK6U+9WdmZqJ3795ITExE3759lVjduxUVFcHCwgL79u3DmDFjZPYR8vjXp34hjz8AdOvWDStXroSfn59Um5DHvlpd9Str7NX6kl15eTlSU1Ph7u5e47i7uzuuXr0q8zXXrl2T6u/h4YGbN2+ioqJCYbXK0pj6q7m5uaFnz57w8vLCxYsXFVmm3Fy7dg2DBw+W/IUE/jP2WVlZyMzMVGFlDTN9+nTY2Nhg1KhR+Mc//qHqcgD85x/0qqoqGBoa1tpHyONfn/qrCW38KysrceTIERQXF8PFxUVmHyGPfX3qr6bosVfrQMrPz0dlZSWMjY1rHDc2NkZubq7M1+Tm5srs/+bNG+Tn5yusVlkaU7+pqSk2bdqEb7/9Ft9++y1EIhHGjx+Py5cvK6PkJqlt7KvbhE5fXx+rV6/Gnj17cOjQIfzlL3+Bn58f/v73v6u6NCxfvhxOTk51/oMi5PGvT/1CG//bt2+jS5cu6NSpE8LDw/Hdd9/BwcFBZl8hjn1D6lfW2Kv9PSQAUqe7YrG4zlNgWf1lHVeWhtQvEokgEokkX7u4uODx48fYtm0bhg4dqtA65UFoY98QRkZGNa739+3bFy9fvsSWLVswefJkldX10UcfITk5GT/88AO0tLTq7CvE8a9v/UIbf5FIhJ9//hmFhYU4ceIE5s+fj4SEBNjb28vsL7Sxb0j9yhp7tT5DMjIygpaWltQnjLy8PKlPI9U6deoks7+2tjY6dOigsFplaUz9sjg7O+PBgwfyLk/uaht7AA36eYVE1WMfERGBI0eO4MSJE+jWrVudfYU4/g2pXxZVjn/r1q1hZWWFvn37YuXKlXBycsLOnTtl9hXi2DekflkUMfZqHUitW7dGnz59kJiYWON4YmIiBg4cKPM1Li4uuHDhglT/vn37olWrVooqVabG1C/LrVu3YGJiIu/y5M7FxQVJSUkoKyuTHEtMTISZmRksLS1VWFnjqXLsly1bhsOHD+PEiRPo0aPHO/sLbfwbWr8sQvqzX1VVhfLycpltQht7WeqqXxZFjL1aBxIABAcHY//+/di7dy/S09OxbNkyZGdnS2aKzJs3D/PmzZP09/Pzw/Pnz7F8+XKkp6dj79692L9/P0JCQtSi/p07dyIhIQH379/H3bt3ERUVhVOnTmHu3LlKr72oqAhpaWlIS0tDVVUVnj59irS0NDx58gQAEBUVBS8vL0n/SZMmQU9PD0FBQbhz5w5OnDiBzz//HEFBQSq5bNHQ+vfv349Dhw4hPT0dGRkZ2LZtG+Li4hAQEKD02pcsWYL9+/cjLi4OhoaGyMnJQU5ODoqKiiR9hDz+jalfSOO/atUqXLlyBZmZmbh9+zaioqJw6dIl+Pr6yqxdSGPfmPqVNfZqfw9p4sSJePnyJTZs2ICcnBzY2dnh4MGDsLCwAAA8ffq0Rv9u3brh4MGD+OijjxAfHw9TU1OsW7dOJc8gAQ2vv6KiAh9//DGysrKgq6sr6T9y5Eil137z5k14enpKvo6JiUFMTAymTp2KXbt2ITs7Gw8fPpS0GxgY4NixY1iyZAmGDx8OQ0NDBAcHq+zDQEPrB4DPPvsMT548gZaWFqytrbF9+3aV3L+Ii4sDAKk/t8uWLUNERAQACHr8G1M/IJzxz8nJQUBAAHJzc/Hee+/BwcEBhw8fhoeHh8zahTT2jakfUM7Yq/1zSERE1Dyo/SU7IiJqHhhIREQkCAwkIiISBAYSEREJAgOJiIgEgYFERESCwEAiIiJBYCAREZEgMJCIiEgQ/h/s0T3e3YPovwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.fill_between(domain,ci[0],ci[2],color='0.5',alpha = 0.5,label = '90% CI')\n",
    "plt.scatter(x,y_noisy,label = 'Noisy realizations of \\nspline-valued process',color='k')\n",
    "plt.plot(domain, ci[1],color='k',label = 'True function value')\n",
    "#plt.legend(loc = 'upper left')\n",
    "#plt.xlabel('x')\n",
    "#plt.ylabel('f(x)')"
   ]
  },
  {
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}