Basic GAM example with formula after merge

ex_gam_mpg_basic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generalized Additive Models (GAM) in statsmodels\n",
    "\n",
    "This notebook illustrates the basic usage of the new GAM support.\n",
    "\n",
    "We are using formulas for the linear part and define additive smooth basis for the nonlinear terms.\n",
    "The following illustrates a Gaussian and a Poisson regression where\n",
    "categorical variables are treated as a linear term and the effect of\n",
    "two explanatory variables is captured by penalized B-splines.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports and Loading data\n",
    "\n",
    "The classes and functions can be either imported through the `api` or directly from the specific modules, similar to other model classes.\n",
    "\n",
    "The data is from the automobile dataset https://archive.ics.uci.edu/ml/datasets/automobile which has miles per gallon as dependent variable and various characteristics of car models as explanatory variables. We can load a dataframe with selected columns from the unit test module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import statsmodels.api as sm\n",
    "from statsmodels.gam.api import GLMGam, BSplines\n",
    "\n",
    "# import data\n",
    "from statsmodels.gam.tests.test_penalized import df_autos"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spline Basis and penalization weight alpha\n",
    "\n",
    "We need to provide an instance of a spline basis class to the model. Currently there are two verified spline basis, B-Splines and CyclicCubicSplines.\n",
    "\n",
    "The model class also takes the penalization weight ``alpha`` as argument. This is currently a model attribute and changing it on an existing instance might cause the results instance to produce inconsistent results if it uses the changed attribute.\n",
    "\n",
    "We specify here alpha values from the unit tests which were obtained by the R mgcv package. At the end we will illustrate the methods to select an optimal alpha."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# create spline basis for weight and hp\n",
    "x_spline = df_autos[['weight', 'hp']]\n",
    "bs = BSplines(x_spline, df=[12, 10], degree=[3, 3])\n",
    "\n",
    "# penalization weight\n",
    "alpha = np.array([21833888.8, 6460.38479])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specifying and Estimating the Model\n",
    "\n",
    "The model class for generalized additive models is ``GLMGam``. As other models, it can be used either with numpy arrays, pandas DataFrames or using the formula interface. The formula interface currently only supports the linear part. Specifying penalized splines inside the formulas is not yet supported. (patsy supports unpenalized splines in the formulas.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:               city_mpg   No. Observations:                  203\n",
      "Model:                         GLMGam   Df Residuals:                   189.13\n",
      "Model Family:                Gaussian   Df Model:                        12.87\n",
      "Link Function:               identity   Scale:                          4.8825\n",
      "Method:                         PIRLS   Log-Likelihood:                -441.81\n",
      "Date:                Tue, 29 Jan 2019   Deviance:                       923.45\n",
      "Time:                        19:46:42   Pearson chi2:                     923.\n",
      "No. Iterations:                     3   Covariance Type:             nonrobust\n",
      "================================================================================\n",
      "                   coef    std err          z      P>|z|      [0.025      0.975]\n",
      "--------------------------------------------------------------------------------\n",
      "Intercept       51.9923      1.997     26.034      0.000      48.078      55.906\n",
      "fuel[T.gas]     -5.8099      0.727     -7.989      0.000      -7.235      -4.385\n",
      "drive[T.fwd]     1.3910      0.819      1.699      0.089      -0.213       2.995\n",
      "drive[T.rwd]     1.0638      0.842      1.263      0.207      -0.587       2.715\n",
      "weight_s0       -3.5556      0.959     -3.707      0.000      -5.436      -1.676\n",
      "weight_s1       -9.0876      1.750     -5.193      0.000     -12.518      -5.658\n",
      "weight_s2      -13.0303      1.827     -7.132      0.000     -16.611      -9.450\n",
      "weight_s3      -14.2641      1.854     -7.695      0.000     -17.897     -10.631\n",
      "weight_s4      -15.1805      1.892     -8.024      0.000     -18.889     -11.472\n",
      "weight_s5      -15.9557      1.963     -8.128      0.000     -19.803     -12.108\n",
      "weight_s6      -16.6297      2.038     -8.161      0.000     -20.624     -12.636\n",
      "weight_s7      -16.9928      2.045     -8.308      0.000     -21.002     -12.984\n",
      "weight_s8      -19.3480      2.367     -8.174      0.000     -23.987     -14.709\n",
      "weight_s9      -20.7978      2.455     -8.472      0.000     -25.609     -15.986\n",
      "weight_s10     -20.8062      2.443     -8.517      0.000     -25.594     -16.018\n",
      "hp_s0           -1.4473      0.558     -2.592      0.010      -2.542      -0.353\n",
      "hp_s1           -3.4228      1.012     -3.381      0.001      -5.407      -1.438\n",
      "hp_s2           -5.9026      1.251     -4.717      0.000      -8.355      -3.450\n",
      "hp_s3           -7.2389      1.352     -5.354      0.000      -9.889      -4.589\n",
      "hp_s4           -9.1052      1.384     -6.581      0.000     -11.817      -6.393\n",
      "hp_s5           -9.9865      1.525     -6.547      0.000     -12.976      -6.997\n",
      "hp_s6          -13.3639      2.228     -5.998      0.000     -17.731      -8.997\n",
      "hp_s7          -13.8902      3.194     -4.349      0.000     -20.150      -7.630\n",
      "hp_s8          -11.9752      2.556     -4.685      0.000     -16.985      -6.965\n",
      "================================================================================\n"
     ]
    }
   ],
   "source": [
    "gam_bs = GLMGam.from_formula('city_mpg ~ fuel + drive', data=df_autos,\n",
    "                             smoother=bs, alpha=alpha)\n",
    "res_bs = gam_bs.fit()\n",
    "print(res_bs.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot the smooth components\n",
    "\n",
    "The results classes have a ``plot_partial`` method that plots the partial linear prediction of a smooth component. The partial residual or component plus residual can be added as scatter point with `cpr=True`.\n",
    "\n",
    "The plots show that both smooth explanatory variables, weight and hp have a downward sloping effect, i.e. cars with larger weight and horsepower have higher fuel consumption and a smaller miles per gallon. However, the effect levels of at higher values, and a linear function would not have been appropriate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg0AAAF5CAYAAAAcQxneAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xd8VFX+//HXCYHQAgqCNKXIioCIhKKAgtiVtbLqRiyr\nYt1VFxv23lbXXta2uipudl3Frl+7gIIoAREUGygISkBUqrTk/P74ZH5TMpNMS2YmeT8fj/uYmXvv\n3DlzU+7nnvI5znuPiIiISE3yMl0AERERyQ0KGkRERCQuChpEREQkLgoaREREJC4KGkRERCQuChpE\nREQkLgoaREREJC4KGkRERCQuChpEREQkLgoaREREJC4JBw3OuU7OuSedcz8559Y75+Y454oi9rnW\nOfdD5fY3nXM901dkERERyYSEggbn3FbAB8BG4ACgN3A+8EvIPhOAvwCnA0OAdcDrzrkmaSqziIiI\nZIBLZMIq59zNwFDv/chq9vkBuNV7f0fl61ZAGXCi9/7pFMsrIiIiGZJo88QhwEzn3NPOuTLn3Czn\n3LjARudcd6AD8HZgnfd+NTADGJqOAouIiEhmJBo09ADOBL4E9gceAO52zh1Xub0D4LGahVBlldtE\nREQkR+UnuH8e8JH3/orK13Occ32xQGJiNe9zWDBRdYNzbbH+Ed8BGxIsj4iISEPWFOgGvO69X1nb\nH5Zo0PAjMD9i3XzgyMrny7AAYVvCaxvaA7NjHPMA4KkEyyEiIiJBY4F/1/aHJBo0fAD0iljXC1gE\n4L3/1jm3DNgH+BT+f0fI3YD7YhzzO4CJEyfSu3fvBItT/4wfP5477rgj08XIOJ0Ho/MQpHNhdB6C\ndC5g/vz5HHfccVB5La1tiQYNdwAfOOcuAZ7GgoFxwKkh+9wJXO6c+wb7EtcBS4AXYhxzA0Dv3r0p\nKiqKsUvD0bp1a50HdB4CdB6CdC6MzkOQzkWYOmneTyho8N7PdM4dAdwMXAF8C5zrvf9PyD63OOea\nAw8CWwFTgYO895vSV2wRERGpa4nWNOC9fxV4tYZ9rgauTq5IIiIiko0094SIiIjERUFDlikuLs50\nEbKCzoPReQjSuTA6D0E6F3UvoTTStVIAm+yqtLS0VB1aREREEjBr1iwGDhwIMNB7P6u2P081DSIi\nIhIXBQ0iIiISFwUNIiIiEhcFDSIiIhIXBQ0iIiISFwUNIiIiEhcFDSIiIhIXBQ0iIiISFwUNIiIi\nEhcFDSIiIhIXBQ0iIiISFwUNIiIiEhcFDSIiIhIXBQ0iIiISFwUNIiIiEhcFDSIiIhIXBQ0iIiIS\nFwUNIiIiEhcFDSIiIhIXBQ0iIiISFwUNIiIiEhcFDSIiIhIXBQ0iIiISFwUNIiIiEhcFDSIiIhIX\nBQ0iIiISFwUNWcB7WLcu06UQERGpXn6mCyAwdixs2QJPP53pkoiIiMSmmoYssN9+8L//wcyZmS6J\niIhIbAoassDxx0Pv3nDppZkuiYiISGwJBQ3OuauccxURy+ch29+L2FbunLs//cWuX/Lz4YYb4M03\n4e23M10aERGR6JLp0zAP2Adwla+3hGzzwEPAFSHb1yddugbk8MNhyBC45BKYMQOcq/k9IiIidSmZ\n5okt3vsV3vvllcvPEdvXR2xfm46C1nfOwc03w8cfw3PPZbo0IiIiVSUTNPzOObfUObfAOTfRObdd\nxPaxzrkVzrm5zrkbnXPN0lHQhmDUKNh/f7jsMhtNISIikk0SDRo+BP4EHACcAXQHpjrnWlRufwo4\nDtgLuBE4HngyHQVtKG68Eb74Ap54ItMlERERCZdQnwbv/eshL+c55z4CFgFHA4957x8J2f6Zc24Z\n8JZzrrv3/tvUi1v/DRwIRx8NV10Fxx4LTZtmukQiIiImpeRO3vtVzrmvgJ4xdpmBdYjsCVQbNIwf\nP57WrVuHrSsuLqa4uDiVIuak666DPn3g/vvhvPMyXRoREckGJSUllJSUhK1btWpVnZbBee+Tf7Nz\nLbGahqu89/dG2T4cmAL0997Pi3GMIqC0tLSUoqKipMtS35x+Ojz7LCxYABGxlIiICACzZs1i4MCB\nAAO997Nq+/MSzdNwq3NuhHOuq3NuGPAcNuSyxDnXwzl3uXOuqHL7ocDjwORYAYPEduWVNh/Fbbdl\nuiQiIiIm0Y6QXYB/A18A/wFWALt771cCm4B9gdeB+cCtwP+AQ9NW2gakc2c45xy4/XYoK8t0aURE\nRBLvCBmzg4H3fgk2akLSZMIEePBByxZ5992ZLo2IiDR0mnsii7VpY4HDAw/Atxp7IiIiGaagIcud\ncw60bWtDMEVERDJJQUOWa9HCOkVOnAhz52a6NCIi0pApaMgB48ZBjx6WXlpERCRTFDTkgMaNLeHT\nSy/BBx9kujQiItJQKWjIEcccA7vuChdfDCnk4xIREUmagoYckZcHN90E778Pr72W6dKIiEhDpKAh\nhxxwAIwcCZdcAhUVmS6NiIg0NAoacohzVtvw6acQMWeJiIhIrVPQkGOGDoXDDoMrroBNmzJdGhER\naUgUNOSgG26A776Dhx/OdElERKQhUdCQg/r2hRNOsGGYa9dmujQiItJQKGjIUVdfDb/8AnfdlemS\niIhIQ6GgIUd16wZnngm33AIrV2a6NCIi0hAoaMhhl15qQy9vvjnTJRERkYZAQUMOa98ezj8f7rkH\nlizJdGlERKS+U9CQ4847DwoL4ZprMl0SERGp7xQ05LhWrWz2y0cfhS++yHRpRESkPlPQUA+ccQZ0\n6WIJn0RERGqLgoYs8NtvNhFVspo2teaJZ56Bjz9OX7lERERCKWjIsLfegn794MADYdGi5I9z/PHQ\np4+NqBAREakNChoybNIkWLAA1q2zvAveJ3ecRo0svfRbb9kiIiKSbgoaMuymm6BTJ3v+2mvw9NPJ\nH+uww2C33Wzq7GSDDxERkVgUNGRY69Zw773B1+ecAz//nNyxnLNETzNnWg2GiIhIOiloyAJHHAGH\nH27Ply+Hiy5K/lh77QUHHGDDMLdsSUvxREREAAUNWePeey1JE8A//wnvvZf8sW68Eb78Eh5/PC1F\nExERARQ0ZI3OncPnkDjtNNiwIbljFRXBMcfYTJi//ZaW4omIiChoyLSyMthjD9hhB3jqKRg0yNZ/\n/bWNhkjWddfBsmVw//3pKaeIiIiChgwbMwY++AAWLoRp06C8HBo3tm033wzz5iV33N/9DsaNs6aK\nVavSV14REWm4FDRk2I8/hr9etQomTLDnW7bAqafa9NfJuOIKa5649dbUyigiIgIKGjKuY8eqry+7\nDHbc0V5/+CE88EByx+7UCc49F+64w5oqREREUqGgIcMmTYLhw6FHD3ucNMnmknjooeA+F18MS5Yk\nd/yLLoImTeD669NTXhERabgUNGRY+/Y2WdWCBfbYvr2tHzkSTjnFnq9ZA2efndzxt97ago4HH7R+\nEyIiIslS0JDFbrklGEQ8/zw891xyxzn7bGjXDq68Mn1lExGRhiehoME5d5VzriJi+Txke4Fz7j7n\n3E/OuTXOuWecc+3TX+yGoU0buPvu4Os//zm5kRDNm8NVV8G//w1z5qSvfCIi0rAkU9MwD9gW6FC5\n7BGy7U5gNDAGGAF0Ap5NsYwN2tFHw8EH2/Mff7TJqJJx8smWC+Kyy9JXNhERaViSCRq2eO9XeO+X\nVy4/AzjnWgEnA+O995O997OBk4DhzrkhaSxzg+KcJWhq0cJe/+MfltchUY0bW2fIV16BqVPTW0YR\nEWkYkgkafuecW+qcW+Ccm+ic265y/UAgH3g7sKP3/ktgMTA09aI2XF27ho9+OO002Lgx8eMcdRQM\nGKCps0VEJDmJBg0fAn8CDgDOALoDU5xzLbCmik3e+9UR7ymr3CYpOPvsYIrpzz+3TpKJysuDm26y\nmopXXklv+dIhNKX2HnvYjJ8iIpI9nE/hltM51xpYBIwHNgCPeu+bRezzEfCW9/7SGMcoAkpHjBhB\n69atw7YVFxdTXFycdPnqm08+scChvNxyL3z6KfTqldgxvIe994affrLjNWpUO2VNxh57hDe9DB9u\nw1BFRARKSkooKSkJW7dq1SqmTJkCMNB7P6u2y5BS0AD/Pyh4E3irctk6tLbBOfcdcIf3/q4Y7y8C\nSktLSykqKkqpLA3BhAnBWoYRI+Ddd60GIREffghDh8KTT8Jxx6W/jMnaYYfwXBI9elj+ChERiW7W\nrFkMHDgQ6ihoSClPg3OuJbAD8ANQCmwB9gnZviOwPTA9lc+RoKuusospwJQp8OijiR9j993h8MNt\nbopNm9JbvlRES6ktIiLZI9E8Dbc650Y457o654YBz2GBwn8qaxf+CdzunNvLOTcQeAz4wHv/UdpL\n3kA1bx4+F8WFFyY3r8T118PixcnPa1EboqXUFhGR7JFoTUMX4N/AF8B/gBXA7t77lZXbxwMvA88A\n72E1EGPSUlL5//bbD44/3p7/+qtNSpWovn0td8OVV1oHxGwQK6W2iIhkh4SCBu99sfe+i/e+mfd+\ne+/9sd77b0O2b/Ten+2938Z7X+i9P8p7rz7wteC226BtW3v+9NPw8suJH+PmmyE/H84/P71lExGR\n+klzT+Sodu1syuuAs86yia0S0batdap86il45530lk9EROofBQ057LjjYN997fn331vHxkT96U82\n1PGss5JLGCUiIg2HgoYc5px1ZGxWmRnj7rvhowS7nOblWWrqBQvg739PfxlFRKT+UNCQ43bYAa6+\n2p57D6eeCps3J3aMnXeG8eNtREVongQREZFQChrqgfHjoX9/e/7pp3D77Ykf46qrbLTCX/6ieSlE\nRCQ6BQ31QOPG8PDDwcyQV1+deCbFFi2seeO115QfQUREolPQUE8MHgznnGPPN2yA009PvMbgsMPg\nkEMs70OiIzFERKT+U9BQj1x3HWy/vT1/+22bWyJR99wDv/xizRUiIiKhFDTUIy1bwv33B1+fdx6s\nWJHYMbp2tSyRd98Nc+akt3wiIpLbFDTUM6NHwzHH2POVKy1wSNT48Tbl9hlnQEVFessnIiK5S0FD\nFvjnP8NrCFJ1552w1Vb2fOJEeOONxN7fpInlbvjwQ3jkkfSVS0REcpuChizw+efW+XDq1PQcr0OH\n8ERNZ5wB69cndowRIyxb5MUXw3LNHiIiIihoyAp/+5tNBX3UUfDDD+k55sknw8iR9vzbb4MJoBJx\nyy32eNFF6SmTiIjkNgUNWSA/H/77X3s86ijYtCn1YzoHDz4IBQX2+vbbYfbsxI7Rrp0FNI8/DpMn\np14mERHJbQoassS228Izz8DHH6dvqupeveDyy+15ebmlmC4vT+wYp5wCQ4fCmWemJ5gREZHcpaAh\ni+y+uw11vPfe5HIsRHPRRdCnjz0vLbXjJyIvzybF+uqr5NJTi4hI/aGgIcucfrp1QDztNPjkk9SP\n16SJpZh2zl5ffjl8911ix9hlF+uoee21ib9XRETqDwUNWcY5G37Zpw8ceST8/HPqxxw2zJoXwEZR\nnHVW4immr74a2raFs8/WhFYiIg2VgoYs1KwZPPssrFoFY8cm3g8hmhtvhE6d7Plrr1nHy0QUFsJd\nd8HLL8MLL6ReHhERyT0KGrJUt25QUgKvvw7XXJP68Vq3tr4SAeeem3gtxhFHwMEH28RYa9emXiYR\nEcktChqy2P77w/XX20RUL72U+vGOOMIWsIRNF16Y2Puds8BjxYrk8j6IiEhuU9CQ5S6+2KasPv54\n+Prr1I93zz3W1ADw6KPw7ruJvb97dwsY7rjD0kyLiEjDoaAhy+XlWXKlbbe1jpHr1qV2vM6d4eab\ng69PPx1++y2xY5x/PgwaBCedBBs2pFYeERHJHQoackDr1vDcc5YOety41EcvnHGGjagAq71ItKkh\nPx8eewwWLrRptBNVVgZ77AE77GCPmttCRCQ3KGjIEX362IX6P/+xUQypyMuz3A1Nmtjrv/8dZsxI\nvDzXXgu33Qbvv5/Ye8eMgQ8+sKDjgw+sBkVERLKfgoYcctRRcMEFtqQ6F0SfPsEahoqK5JoaLrjA\nUkwffzysXh3/+378sfrXIiKSnRQ05JibbrJpq48+GpYuTe1YF14IAwfa8/nzEx/a2aiRpbteudKG\ncMarY8fqX4uISHZS0JBj8vOtiaJJE/jDH2DjxtSO9a9/QePG9vqWW2zCrER0724jMv71L5twKx6T\nJtlU4D162OOkSYl9poiIZIaChhzUvr1doGfNgvHjUzvWzjvDVVfZ84oKm/ci0UDkhBMsgDn99Phq\nP9q3t34QCxbYY/v2CRdbREQyQEFDjtptN7vD/8c/bEhmokJHMLz6KvTrZ+s//9w6OCbCOZsJs6DA\n+kZUVMT/2Ro9ISKSOxQ05LBTT4WTT7YhlLNnJ/be0BEM06ZZ/4RAM8Xf/gYzZyZ2vLZtrYnizTct\nmIn3szV6QkQkdyhoyGHOwX33Qd++duFduTL+90aOWFi9Gq64wp6Xl1uNQaLNFPvvbx0iJ0yAefPi\n/2yNnhARyQ0KGnJc06Y2I+aaNXDssfHPiBltBMPFF8OAAfZ63jyb9yJRN90EPXva7Jyxgg6NnhAR\nyU0pBQ3OuUuccxXOudtD1r1XuS6wlDvn7k+9qBJL1642ouKtt4KdGmsSbQRD48bWxJCfb/vcdJN1\ntkxEs2bw1FPwxRfBmot4PltERLJf0kGDc24wcCowJ2KTBx4CtgU6AB2Bi5L9HInPvvvCDTfY8sIL\nNe8fawTDLruEN1P86U+waVNiZenf32op/v536+MQ72eLiEh2SypocM61BCYC44Bfo+yy3nu/wnu/\nvHJZm0ohJT4TJtjU1yecAF99lfxxLrkEdt3Vns+dm1wzxfnnw377WTPFDz8kXxYREckeydY03Ae8\n5L1/J8b2sc65Fc65uc65G51zzZL8HEmAc9a80KmTBQ9rkwzVGje2eS4CzRQ33ph4M0VenmWLzM+3\nvhZbtiRXFhERyR4JBw3OuT8CA4BLYuzyFHAcsBdwI3A88GSS5ZMEtWplfQQWL4ZTTkl+Rsxdd4XL\nLrPngdEUiTZTtG8PJSUwdWriuR9ERCT7OJ/AVcU51wWYCeznvZ9bue5dYLb3/rwY7xkFvAX09N5/\nG2V7EVA6YsQIWrduHbatuLiY4uLiuMsnQc88YxNc/f3v1lSQjE2bYPBg+PRTe33VVYlPow3Wz+KK\nK+CNN6zvhYiIJK6kpISSkpKwdatWrWLKlCkAA733CdYJJy7RoOEwYBJQDrjK1Y2wzo/lQIGPOKBz\nrjmwFjjAe1+lW1wgaCgtLaWoqCipLyHRTZhgU1e/9RbstVdyx5g92wKH8nJravj442B/h3hVVMCB\nB8Inn1gzR5cuyZVFRETCzZo1i4E282CdBA2JNk+8BfQDdgX6Vy4zsU6R/SMDhkoDsKBCKXzq2A03\nwMiRNiPmkiXJHWPAALj0Unu+ZYuNpti8ObFj5OXZMMymTS0bZCqTbImISOYkFDR479d57z8PXYB1\nwErv/XznXA/n3OXOuSLnXFfn3KHA48Bk7301OQKlNgRmxGzaNLUZMS+/PDg3xZw5lr8hUe3aWV+L\nOXPg7LOTK4eIiGRWOjJChtYubAL2BV4H5gO3Av8DDk3D50gS2rWzjJGzZ1uK52Q0aWKjMho1stfX\nXWcX/0QNGmQTbD38sGWN1IRVIiK5JeWgwXu/d6ATpPd+ifd+L+99O+99c+99L+/9JcrTkFmDB9sc\nFQ8+aEMpk1FUZPkbwJophgyxbI6JXvBPOgk6dLDETpqwSkQkt2juiQZi3DhbzjwTSkuTO8bll0Pz\n5vZ80yabHTOZC36ziKwdmrBKRCQ3KGhoQO65x/omHHkk/PRT4u8vKICttw5f9+GHiTcxdOoU/nrb\nbRMvi4iI1D0FDQ1IYEbM9euhuDj+GTFDdesW/rq8PPEmhsCEVZ06WRbLbt2ST0IlIiJ1R0FDA7P9\n9vDf/8I778SehbI6kybB0KFV1yfSxBCYsGrpUutgWVJiSajSoazMaj7UyVJEJP0UNDRAe+9twyZv\nugmefz6x97Zvb30Z+vcPX9+yZXJlOeEES1c9YULiZYlmzBir+VAnSxGR9FPQ0EBdeKFNanXyyfD9\n94m//403YLvtgq9//RVWrUquLNdeaxf7sWMTnxgrUmSNhzpZioikj4KGBso5eOQRGw1xwgmJ929o\n397u5ocNs9eLF8NZZyXXNyEvDx5/HPr2hUMOST57JUDHjtW/FhGR5CloaMDatLHpqydPTq5PQX6+\npYcOzDP273/b8ZLRvDm88IId84AD4OefkztOoJNljx72OGlScscREZGqFDQ0cKNGwUUXWQ6GmTMT\nf3+3bpY0KuDPf4ZvvkmuLB07WrPH8uUwejSsW5fY+8vKrA/Djz/asSZNshoRERFJDwUNwrXXWsfG\nY49N/EINcMwxlukRYO1aG865aVNyZenVC157DebNs/kyEpkcS50gRURql4IGoUkTa2ZYuhT++tfk\njnH33bDjjvZ85szkhnMGDBoEzz0Hb79twUhFRXzvUydIEZHapaBBALvDv+su6xyZTD+Ali0t30Lj\nxvb6llvgzTeTL8+++8LEidZP4rzz4utgqU6QIiK1S0GD/H+nnGJV+qeearUOiSoqCp82+4QTYMWK\n5Mtz9NE20dZdd1kTSk3UCVJEpHblZ7oAkj2cg4cesv4NJ5xgNQV5CYaV48dbZ8Y33oBly+CPf4QN\nG+x5Mp0TzzzTckBceqnVYlx6aex9A5kmRUSkdihokDBt28ITT1jzwG23WRKoRARyLvTvb6Mg3nkn\nuG3hQsvDMGNGYse85BLrEHnZZTYk86KLEnu/iIikh5onpIq997Zg4bLLksvQ2KGDdax0ruq2jz9O\nbj6IK6+0zpUTJlgTiCa4apg0t4hIZilokKiuu86m0S4uTm4Y5r77wgUXVF3vvdU2JOOaa+Dqq62J\nYsIEBQ4NkYbVimSWggaJqkkTG7mwZImNXkjGu+9GXz9nTnLHcw6uugruvBNuvRVOOy256b0ld2lY\nrUhmKWiQmHr1sgv0Qw9Z3oREzZ2b/jIBnHuu9Zt47DGrCdm4sXY+R7KPhtWKZJaCBqnWuHE2G+a4\ncfDDD+k5ZqdOqR/jhBPg2WfhxRfh0EOTa0KR3KNhtSKZpaBBquUcPPwwNG1qF+p4szOCjaCIZulS\n6xAZS7yd3Q47zFJOT5sG++0HK1fGXzbJTYFhtQsW2KPmFhGpWwoapEaBYZhvvw23317z/oGLflkZ\nFBZCly722LKlbd+0yWovli2L/v5EOruNGmXDOr/+GnbfHb74IvHvJyIi8VHQIHHZZx8bDXHppTB7\ndvX7Bi76ixbBmjWwapU9rl0b3GfpUpurYsmSqu+P7Nz2/ffV1zwMHgwffQQFBRY4pJK+WkREYlPQ\nIHG74QbYeWebDXP9+tj7RV70f/st+n5r1tjkVJEiO7f98kvNNQ/du1szxbBhcNBBln4605RTQETq\nGwUNErfAMMxFi+D882PvF3nRb9Ys9r5lZXDvveHrIju7tWkTvj3WMLtWreCll+Dss+Evf7EU1Jkc\nWaGcAiJS3yhokITstBPccQc88AC88EL4tsCd9ZIl1oeha1e76E+fHgwCCgurHvPcc+Hll4OvQzu7\nPfts1Q6O22wTu3yNGln5HnoIHn0URo605o1MUE4BEalvFDRIwk47zUYunHJK+DDMyL4MXbrYxb9v\n32AQMG2aBQ75+VZzATYi45hjoLS06meNGRPeFwLiywR56qkwdaqVr6gIXnkl+e+bLOUUEJH6RkGD\nJMw5eOQRu+ifeGJwGGY8d9ZnnGEBxZYtNooiUGuwfj38/vcWcNR0jHiHVg4ZYnNn7LabHXv8+Lpt\nrlBOARGpbxQ0SFK22SY4DPOWW2xd5J30ihVVO/9FBgGFhdakATYEc/RomwobrLljxYqqn53IHfs2\n21g/hzvvhPvvt46Sn30W//sjJdK5UTkFRKS+UdAgSdt3XxuCefnl1hQwaVJ4n4U1a6p2/ou84Hfq\nBM8/D7/7nb3+7DN7z4YN1jSxZk1wX+es9iDRO3bnrN/E9OlWozFgAFx/vU23nSh1bhSRhkxBg6Tk\n6qut6r24GPLyoF278O2RNQvRquzbtoVXXw02Vbz7rvVxiExbvd120LgxDB2a3BDGoiLLMXHBBVbu\nIUPgk08SO4Y6N4pIQ6agQVKSn2/DMDdutDTTbduGb1++PPziHqvKvmdPG0HRooW9fvHFqh0g48nX\nUJOmTeHGG2HGDJshc/BguPLK+Ps6qHOjiDRkChokZZ07w8SJNg9E5J332rXxX9y7dbMlYMUK6NAh\n8XwN8Rg4EGbOtKaVm26yJou33675fercKCINWUpBg3PuEudchXPu9pB1Bc65+5xzPznn1jjnnnHO\nqQtYPXfAAXDJJfGlhY5lzJiqnRSXLbOmivfftyGcoVK9y2/SBK66yoZ6tmljfTSOOgoWL479HnVu\nFJGGLOmgwTk3GDgVmBOx6U5gNDAGGAF0Ap5N9nMkd1x7bfTkTfFe3GMFFzfdBH/7W+3d5e+yi3Xk\nnDjRmj122smCidBOmCIikmTQ4JxrCUwExgG/hqxvBZwMjPfeT/bezwZOAoY754akobySxfLzYcoU\ny8qYl2cTSNU02iF0CGNZWfi20JqFiy+2+SRC7/K9T9/cDs7B2LHw5ZeWhvpvf7N+FvfdZ/kkREQk\n+ZqG+4CXvPfvRKwfBOQD/7912Hv/JbAYGJrkZ0kO2XVX69BYUWE1DzNmVF+Ff8ghwc6N69aFb4sc\nPXHttXDzzcHXtTH8sbDQAoavvoKDD7YAok8f+O9/g0msREQaqoSDBufcH4EBwCVRNm8LbPLer45Y\nXwZ0SLx4kosOPNBqBi691NJGV+fTT2Nvi3aRvuQSy7Hgfezhj/Pm2eRVjRvbYzLJnLbfHh57DObM\nseaKP/7Rak1eeim+NNYiIvVRQkGDc64L1mdhrPc+kdQ4DtC/2gbkuutg992tE2O0tM+BZolk0jpf\ncYWlhO4QEYZ27GjH3XXXYKrqNWssr0Oy+vWzmpP33oPmzeHQQ23kxfPPq+ZBRBoe5xO4bXLOHQZM\nAsqxQAAHV+AhAAAgAElEQVSgERYQlAMHAm8BW4XWNjjnvgPu8N7fFeWYRUDpiBEjaN26ddi24uJi\niouLE/k+kkWWLLELeFGRJW/Kzw9u22MPa1KIxbnwO/rCQgsUrr02uO6II2x0RVmZBQyTJlkTReRx\n8/OTy/4YzXvvWRnefdc6UF5xhZWjUaP0HF9EJJaSkhJKSkrC1q1atYopU6YADPTez6r1Qnjv416A\nFkCfiOUj4HGgN9AK2AgcEfKeHYEKYEiMYxYBvrS01EtuW7bM++HDve/Rwx7Lyrx/6y3v8/K8v/DC\n8H179PDewoLoy+DBVY/lvfePPGLHC+y3777er15d/XELC9P/XadMsc8G73v29P6++7xfty79nxNN\ntPMsIg1TaWmpx27ci3wC1/Nkl9QPAO8Ct4e8vh/4FtgLGAh8AEyt5v0KGuqJ4cPDL9bDh9v6226z\n1yUlsfcNLPn5NV8IX3zR+6ZNg+8ZODC4f+RxGzXyft682vvOM2Z4f/TRFshsvbX3553n/Zdf1t7n\neR/7PNc1BS8imVfXQUM6MkJGtm+MB14GngHeA37AcjZIPRerY+L48Tac8eSTrWMhBHMuFBSEv2e3\n3WpOmnTIIZa9ceut7XVpqc1euWBB1VwOP/wAffum5/tFM2SIjaz4+ms45RT417+gVy9LFPXss+lr\nFgmVLfNfaPIukYYn5aDBe7+39/68kNcbvfdne++38d4Xeu+P8t6nMIJeckWseRmcg4cesovp4Ydb\nx8hAZsXFi5NL2DRsWHiWyAULLOCYOzd9GRsTmQa7Rw+49VZYutSmDF+/Hv7wB+ja1fo9fPllYser\nTrbMf5EtwYuI1KG6qM6obkHNE/VGWVn11dXffef9NttYX4DNm6s/VrxV34sXe9+nT7CqPi/PmkMq\nKlL/Pqk2A3zyifdnnOF9q1b2/ubNw483eHBy5arpPNeVbGkmEWnI6rp5IqHRE7UhMHqitLSUoqKi\njJZFat+778J++8Ff/gJ33hl7v8jRFcOHW81BNL/+as0fr74aXDd2rNVuNG+efFl32MGq3gN69LAa\njET99hu88ooNPw0dppmfbzUTuTp/xfLl1iTx44/B0Su5+l1EctWsWbMYOHAg1NHoCc1yKXVq1Ci4\n+2646y5L0RxNWZnNQBmquqrvrbayqbQvuyy47qmnLPBI5iIfkK5mgGbNrKmicePw9Vu2WK6JPfeE\n228PD1BygSbvEml4FDRInTvrLPjrX+Gcc8JrBwLGjKma9KmmC3ajRpYp8plnoEULWzd7tk15/Z//\nJFfOdE+Q1b9/+OsBA+CRR6xD56WXWs1Gnz5w7rlWM7F2bWqfJyKSbmqekIwoL7eq7XfesbvU0Atq\nZLNAQYF1mIz3Tvazzyzh0tdfB9cdd5zVbrRpk57yJ6O66vy1a+H114PL4sVWMzF8OOy/v9XQDBxY\ntbZCRBq2um6eUNAgGbN2LYwcaRfTGTOgUydbn0h/hljWrIEzz7RmioAOHayfwyGHpF722uS9BTxv\nvAFvvmmB1dq11j9j2DAYMcKW3XaDpk0zXVoRyST1aZAGo2VLmwAK4Pe/h9WVicfT0SxQWAhPPmnD\nHwPZyZcts7kjjj8efv45uG9gKGTXrjbBVbdusYdEpmvYZHWcgx13tM6iL7xgZZ0xA66+2oKE226D\nvfay7zViBFx+uQUYqyOniRMRSTPVNEjGffqpXfwGDIDXXkv/3fPSpXDaaeH9Jzp0gAcftCAi1jwY\n0Wo40lELkqrycpvJc/JkmDLFlhUrgsHGoEG2DBxo57Rly9ovU1mZ9UXJlpEU2VYekdqi5glpECL/\nqV98MRx9tA3HfPbZ8Mmt0sF7q3U491xYtSq4fswY+Phj60MQTSBjZf/+VisydGh6hmGmk/fwxRfw\n0Uc26mTmTPjkE9iwwQKJ3r2DQcSgQTaJWCpDUaPJhmAqm8sjUlvqOmhI879mkfgEUhCDXYRvvtmC\nhUMPhXHj4NFHIS/FxrNod5v77gunn26jE8A+s7rPCYzi+Ogja0JZsSJ8e6ayMYYKBAa9e8OJJ9q6\nzZth/vxgEFFaaqNINm2y79unjwUPgff17m1NLsl2tMy27JDZVh6R+kJBg2REtH/qBx1ktQFjx9ow\nxNtvtwtisiIDkyOPtLvNl16CiRPhggusT0Ig4VKjRvY8VuXbnDl20Q1o2TL1YZi1pXFjm7p7l11s\nzg+wsn/2WTCQmDsXXn7ZkmMF3tOzZ3gg0auXNXm0alX953XsGF4Dk6lgKhAoLl0avj4bgjuR+kBB\ng2RErItMcTH88gv8+c/Qtq118ktWrLtN56wz5KGHWufCe+6xfgLl5dUfL3Lyqc2bc6udvEkT6+Mw\nYACceqqt894utPPnhy+PPWaTfQV06GDBQ+TSo4c14UyaVHU4aSaEBopgZRs0KHuDO5Fco6BBMqK6\ni8xZZ9mIgSuusH/6F16Y3GfUdPfbujXccYfNTnn22fDee+HbGzWyxTnr0zBnTtWkU6mqqcNebXfo\nc84Cgg4dLBdEqNWr4auvwpfZs21WzzVrbJ+8PBtt0quXXZwDwcSGDVZrk2oTU6IiA8XOndWXQSSd\nFDRIRgRSEMdy2WV2gb7oIrsbvuiixD8j3rvfnXe2XAgvv2wdMj//3NaXl9tnn3SSBTBHH219GwIi\nMzzGKzQQWLEieAFeuBC23z54Z9y+fdUmlt//3moM6mJUQKtWwZEYoQK1E4FA4ssv7fH11+H++4M1\nMk2bwu9+F14z0auXpf0+9dTa+Q7Z0kwiUl9p9IRkndCL6ubN8P331lFywoTa/+wtW6xfxZVXhreL\nN25s03q/9prdRTdrBtOnQ9++iX9GrCGeoQK9/aNlxwyt7ci2UQFbtsCiRVVrKL76KvYIle23t/4r\nu+4K3bunVjuhSbSkodHoCWnwItult9vOagC8t8falJ9vHQf/+EebWOvmm22I5ubN8L//Bfdbs8Zq\nIJK564+nJ39gn8g752SOVZfy8y3Q2WEH69gaav16+OYbS4tdVhZcv2SJTegFlpSrf38LIHbdFXbf\n3TpkxgokQgPMtm2tueWnnxQwiNQWZYSUrBN5IWzcGK66Ci65xDou1kXlWPPmFqAsXGidMQsLq+5T\nWmrBzcKF9njkkba+pqyRkVXmhYXBfBCR+0Rmx4xsEsml6vfmzW00R8+e4euHDrVsna+/bud6u+3g\n7bctIdfOO0O7dtZp9dZb4cMPw0ewBALMhQst38ZHHwV/HtmeLlwkF6mmQbJOtHbpQArlSy6xfgB3\n322dFGtbmzZw3XUwfrxdsJcsCW4LDNUMCAQ7kf0QIvspxOprEW1dZN+PaNXvuSba92/f3mog9t8/\nuN+6dZY+e+pUW66+2mormjWzGogRI+Dbb2N/zpw5tf5VRBoc9WmQrFNdu/TDD9tEVAcfDCUlwWmw\n66pchx5qHf/WrKk6RLNVK0tKdeGF0S9mmex/UB/SKm/ebKM33n/fgojJk214biwFBdb/RKQ+Uxpp\nkRq89hocdZRlNQxMeFWbF8R582x2yd9+swvRjjvaxaq83KrVI/M3NG5cdR1UTTldlxfy+phWubwc\n+vWzvBLR9Otn85oE1IfASSSSOkKK1OCgg2ySptGjrT18661hVuWfSmjmx3QZNiw4LHLLFrvbDRg6\n1NJe33qrzf8A0QMGsAtVdcMt013uUPUxrXKjRtHzZgSyiM6da4HloYfCYYdZBtBp02xbbZ9vkfpK\nHSElJxUVWae4Zs1scqZQ6b4g/vZb7G1Ll1qTxKZNNvwytE0+VPPmNu/FIYcEO+4FAoaA2ryQR3aY\nrI0OlHUxbXikaN/De9htN3j+eev78OijFvgFAoaA+hA4idQ11TRIzura1e4Uu3Wz7IUB6b4gNmtW\n9QIf8Msv4fkHvvjC9o8MNNavh2uuqf5zanMkRG2meQ7UnsycGbzzr6s7+UmTYNttq65fvtxqFw47\nzJoxWrWyn0GoX3+1dNmjR6uZQiReqmmQnLb11tbnIHDhaN/eJqNKp+nTq47UKCiwfgFt2oSvLy+v\nGjA0aRL72Pn5VvahQ2t3JERgFMaCBfaYzotkYLRIZFNBrDv5WDUSydRUtG8ffThsaADWqFH4MM2A\nHXawFOIdOtjn3XqrdXIVkdgUNEjO224765D4r39ZjcBhh8HXX8feP9GLU9++NnlTaL6ExYvt4tul\nS3xljJWcaMsWK89nn9moixdfrHpHnO1iBQexak5CcyuE5reItb4m06cHR9E4ZxNyRQZgzZqFvy4s\ntJwOy5bBI4/ANttYLpCddrLloousDDVNYibS0ChokHrjxBNtXP9vv8HAgeEZHEMdemj4xSmeJECR\nd+reW8ARmrchIC/PAotAwqZNm6rmdGjZ0ibMCli92tJXH3aYZTY85BB46KHwmSazVWRwEKiFiVVz\nEqtTZrKdNfv2hbVr7WdSUWGdYiNrUqZPt0AhP98ep0+39e3bWwbQ55+3TJIvvGA/18cft8eOHW37\nCy/kXjAnUhsUNEi90q+fta0feKBNMHXiiZYGOlRk0p+akgBFq5kI3BUvWlR1/06dLLDo3Dn2Mfv3\nt+O+/DKccIK1uQds2GDrTz/djjFokM2FMXVq9Gr2dEilE2No1srBg+27/fij1RREO06sTpm12Vmz\nb18LzDZvtsdttqn6fZs3t4DykUcsWJs2zQKG6dNt3pG2bW37P/8ZngZbpEHx3md0AYoAX1pa6kXS\npaLC+3/9y/vCQu+32877l14Kbiso8N7uS20pKKj+WIMHh+8/eLD3PXqErwtdhg+39w0fHnufAQO8\nLysLfsaGDVbG007zvmPH2O9r0cL70aO9v+MO7+fOte+ZDpFlDXyH2jhOWZmt79HDHgPnIdb62pDo\n9/3yS+9vvdX7PfbwPi/Pe+e8HzjQ+65dve/SpfbLKxJLaWmpBzxQ5Oviml0XH1JtARQ0SC1Ztsz+\nsTdrZr/po0d7v2RJ9CCgOtGCjMiLTmFh7Itg167eN2oUO7iIVF7u/ccfe3/llRZcxAogwPsOHbw/\n7jgLkBYvTv5cRQZBPXrE975ly8Iv9F27Jnecupbs9/Xe++XLvX/sMe/btAk/RseO3r/8svdr19Za\nsUWqqOugQRkhpd6KzILYuLHNX3HxxfDKK9YJLp7MgE2bho8MKCiwjpCJTMEcOcU1VM0QGcvSpfDW\nW8Fl2bLY+3bubCMxdt/dHouKrPw1STZjZOT7CgvDh6fWdJy6zNIYK7FWPOWMJvJnmp9vHVubNIE9\n97QmsgMPtKaRQMIpkXSr64yQqmmQeivybrKgwPtWrex5//7ez5oV33GGDAk/zpAhtj7yLru66ulo\nTRWJNgEEPq9LF++7d/d+n32suaK6mojGjb3fbTfvzz3X+//8x/tFi6I3aSTbNBB5jrt2Tew46WoW\niUc8tUOpHG/4cO+/+ML7u+7y/qCDvG/a1NZ37uz9ySd7//TT3v/8c/q/lzRsap4QSZPq+hQ0b25t\n0+PGef/DD9HfH7hId+1qF5jABTFwgUnkgldWZs0gBQW2DBlS/YUqWkAS7fM2bvR+yhRrythnH+9b\ntqw+iADvO3Xy/sgjrY1+6lTvV69O9gynftFPpZkgUen+rJoCrfXrvX/9de/PO8/7Pn3sM/PyvB82\nzPtrrvF+xgzvt2xJrQwiChpE0iT0n3pkv4Tu3b2/805rl27e3PurrvJ+zZrw99d0QazNC160z478\nvPz8qherLVu8nzPH+wcf9P5Pf/K+V6+ag4hA2Y84wvurr/b+uee8X7gwvk6WqXZeTCToSKRmJ9XP\nqg2LFnn/0EPejxkTrPFq29b74mLrk/Ljj3VbHqkfsjpoAM4A5gCrKpdpwIEh298DKkKWcuD+Go6p\noEFqXawLxi+/eH/RRRZUdOhgF9vNm21btKAg9MJVWFh7F7xonx2r5qSmi9/Kld6/+qr3V1zh/b77\nVi13rKVVK2vaOPpoq2XZdlvvi4qsM2m6JBJ0pHrRr8vRGTXZtMlqeS67zPtBg4LfqV8/7885xwK3\nbGjKSDVQk9qX7UHDaOBAoGflcj2wEehduf1d4AGgHdC+cmlZwzEVNEitq+mC8d133o8da38RvXtb\n+/OwYVUvUsm2iyd6wYu2f+A75OdXDSjiEbgAdO9uozL+/ncb4rnbblbbEk8gEVh22MH7ww/3/uKL\n7S55+vTav8jVZVNGXSsr8/7JJ612KDACxTn7OZ13ng3H/fXXui9XpmtnpGY5N3rCObcSuMB7/5hz\n7l1gtvf+vATer9ETkjVKS+HSS+GNN6B3b/tXuWlTsGf/0KHhPebjGQFRVmaTa4WOwKjpfcuXxx6d\nkc6RDt98Y8ctL7fyzJljy6ef2hIteVV12rWDXr2qLj162OiVVCT7vXPRd9/Bu+8GlyVLLNPowIGw\n114wapSdj2jzbqRT5AiReEf8SN3JmdETWDbJPwK/Ab18sKahDFgBzAVuBJrVcBzVNEjWef9961gI\n3u+6q/fPP29t/MnceUVrVigsTL6qN10jHQLlmDs39vG6dUusBiLWkp/v/Y47en/IId5fcIH3Dz9s\nHTiXL6/97x2vZKri66L6vqLC+6+/tnN27LHet28fPK/Nmnl/1FE2YuODD7xfty61z4r8PpEjh1TT\nkH2yunnC20V+Z2ANsBn4mfA+DeOA/YC+QDHwPfBMDcdT0CBZa/Jk7/fay/5Sioq8f+QRa7aI5yIR\n+Acc2ZyQqX/A1Y0miVWuWO8ZNsz7b7/1/rXXvL/9du9PP93OU3XZLGMt227r/X77WTDxxBMWxAT6\nldSGWBf6dASEdfEzjfzMli29b9LEnufleb/zzt6feKL399zj/bRpiQUSkccePFh9GrJdLgQN+UCP\nyov9DcByYKcY+47COkN2r+Z4Chok673zjvd7721/MTvu6P0//2nDHatT00U6Wpt8bd65lpXF1wky\ntFxz59p7GjWyJZ6UyatWWUbLiROt8+VRR3m/yy7BzJzxLM2a2ee2bu39Tjt5P39++s5VrAt9Mn0m\nMtHPItpnbtxoeUceesgCuIEDLUcH2M+tXz/vTzrJ+3vv9f7DD204aLZ8H0lNXQcN+Uk0Z2wBAq1c\ns5xzQ4BzgTOj7D4DcFinyW+rO+748eNpHTrtH1BcXExxcXGiRRQB0pttcNQoWz76CG66CU45xaZS\nPv98OPXU4NTMoWqapTHahEyBibDA2pKPPLLmdvt4v2f79taHoWfP8GyIkZYuhSFDLIvhnDnhfTG6\ndq1anmifP2iQLaEqKuD77+HLL4PL559b34mVK8P3/e234PNVq6x/Sd++lmkxsBQXJ36uIPZsmh07\nhrffxzNhVjLvSUVZmWWzjCxDkyY2JfiAAfb7CPZzmzfP+unMnGmPEyfapF2NGsHOO1sfiUGD7HGX\nXer++9QF7+tPRs6SkhJKSkrC1q2KnJGvtqUadQBvA4/G2DYcq2nYuZr3q6ZBakVtVh1/9pn3J5xg\nd3Ft21qynpUrw/eJbA8eMKDmO+Pq7vTSVa0eb41DvLUjqZ7nigrvly61po6bbrKaiVhNOqFLsqNI\nYpU3mT4TdT2MM9ronUQ+c8MGqwX6xz8ssdmuuwbPY36+9337WnNR27Y2imjatJpr1LLFxo2Wo+SJ\nJ6ypa7/97LvMnZvpktWurK5pcM7dALyG9VUoBMYCI4H9nXM9gGOBV4GVQH/gdmCy935eqsGNSKJi\n3VGmQ58+8PjjcO218Pe/W+3D3/5m01yffbZtt5g4KD+/5jvh6u70YtVCJPo9AzUORx5pvfJ//tlG\niITWKFRXPgivXVi6NLHPj+ScTSfeqZPN1QCw225WqxPQooVNGV5eHly3ZUv4cX76Ce68E37/e6tN\niWXSpKojU8DOS6KjMZJ5Tyoiz227donVnhUUVK0F2rDBantCayS++MJqf4YNs1qJrl3tnP7ud/bY\nsyd06wbbbx8+rXtdWbYsOMInMNpn/nyrRQEr2y67wGmnZaZ89VmizRPbAk8AHbHkTp8C+3vv33HO\ndQH2xZoqWmCBxf+wfg8ida4uqlq7doV77oErroAHHoB//MMe99vPJrUKFVkFH02sCxpUvWB8+KH9\nQ9y0KXx9PN8z8mIXbUKtUIGLTaA8oQFMpHSc58iAq08feOcdmD4dpk615b33wvdZvRrGj7elXz8r\n45gxVSeMSuZCX5cTa1WnNn6nmza15qghQ4LrKiosGPz6awswA4/vvQf//KcFGgEtWsBWW9W8tG4d\nfV2TJuHl8d4Cwg0bbPnhh2BgEHhcvjz42f362QRtp59ugUK/fnZcqR2a5VLqrepyHdSWTZvgf/+D\nu+8Ov1OG1PMKROYpiBR6YU/0e0Yeu0ULu9j+9FP0cxcZZBQU2Ayb6TrP8eQH6N7d8hnUZMcd7fdg\nzBhru0+mfTtbckRk4nc6UkWFXcgXL7blhx+s38mvv1ZdAuura3Zv3hxatrRagkCgEO2y1KOHBQW7\n7AL9+9tjjx6Wv6Ihq+s8DQl3hBTJFXVddQx21zR2rC2vvWbNFT/9ZP/Ytt8epk2zBFHJXLgCtRAz\nZlStmge7aCf7faPVcFR3MYq84x00KL3nOp476s6dw4OGAQPgj3+0ss+YEVz/1Vdw8822bL89HHWU\n7ZdIAFGbTV2JyMTvdKS8POjSxZZhw+J7T3m5db6NDCYCy5o19rfTtGlwKSiwx3btrNOmmhmyg2oa\nRGpJoEr7++/tdUWF9SHYaSc4+WQ49li78CUqVo1DXd791vYdbzzHr26fJUvguefg2WetKaOioupn\n9OgBxxxjyy67VB9AZEtNg0ikuq5pUNAgUksiLzTDhsE111ib8HPPWVPGXntZx72nn7ahdPFcgAMX\ny0AnxjZt7K4vU+3skbKl/T9g+XJ44QULIN5+O3otTa9ewQCiTx9bF/o92ra1oCJWc01DlG0/54ZK\nQYNIPVFdu/yvv1rg8NRTdiELNXSoNWNEypV/0tl8V75ypZ23//7X5nSIVgPRr58FD5MmwayQf8GF\nhVZVns3nvi5l88+5IVGfBpF6orp2+a22gpNOsqVr1/CRFh9+aH0iRo2yWonly+29mzcHO1cmksyo\nrmVL+380bdta8qNTT7Ug7NlnLYCYOjXY+W7uXFsirVljSzaf+7qUzT9nqT0NvN+pSO2ZNMnuvnr0\nsMfQ4ZOhttuu6uvPPrML24cf2kXqgw9suFmobP0nHdlpMVuzCm67LZx1FkyebP1O7rjDhu7FI1vP\nfV3KlZ+zpJeCBpFaEujpvmCBPcaqzo4MLj7+GD75pGowEZl8KVv/SccbLGWTzp3hr3+1PBDffQe3\n3GKdI2P56Se47z6rrWiocvHnLKlTnwaRLBXZZtytmw1V++UXe921K4weDQccYE0ZhYUZKWZSovXP\n8D77+mx88411Un3qKZsnI1JennVmPeYYOPzwzJdXGh51hBQRIPaQwrVrLTPf//2f5YJYuBAaN7aq\n9ZEjYcQIG6kRbRKtbBGtEx1kd8e6+fOt/8N//2tpliPl5dn3OvJIOOIIywkhUtsUNIjUoVwZkVCd\nb76xAOK992DKFBu6mZ9vyYsCQcQee2RXat1oI0ug5iyQ2cB76ygZCCBilXHgwGAA0bt36p9bH35X\nJf0UNIjUofo2bMx7uwuePNkCiMmTLc2vc7DrrhZAjBxpU0tvs03myhl53gsL7cIamno7F34W3tuw\nzGeftYv4l19G32+nnYIBRK6nspbsoqBBpA7FM8dBLvPevl9oEBFIvdy3rwUPgcmKdtrJZjSsC8uX\n20yJa9YE1w0ebKmE47mTzta77vnzrSyROR5CbbcdHHKI9UcZNQqaNYvv2PX9d1WSo6BBpA41xLu3\nxYstL8Hkyfbd58+34KJlS7sLHjzYgojBg62zZTJ3xfFI5SKYCz+3776D55+3AOL996NPwtSsGeyz\nj2UFHT3aMnvGkgvfWeqeggaROpQNswZm2po1UFpqQz0/+siWQLKp1q0tQ2K/fsFph3feOT39I5K9\nCJaVWTATOgQ12++6y8rgxRetGeOddyxRVzT9+1vwMHo07LZbeM2PflclGgUNIpJxZWUwcyZ8+ql1\n+vv0U+srUV5u2zt3tj4IffrYY+B5u3bxf0ayF8FoE3bl0l33mjXw1lvw8svw6quwbFn0/dq2hYMO\nslqIAw6wLKIikRQ0iEhW2rjRAoe5c61JY/58y13wzTfBYKJt26qBRO/eVu2ermaOyGaNggKrGUnn\nXXdZmfU7+PRTe92/P7z0Uvrv7CsqYPZsCyBeecVqe2Jp2tSyWD78sDVp5Ck1n6CgIaNlEZHEbdoE\nX38dHkjMn28BRqAJobDQOlpG1k706JF458u6aNvPVG3GsmWWe+OVV+D11y0nRzTbbGOBQ2BIbe/e\nCiIaKgUNIlIvlJdbZ8DQQCLwPDBqoqAAdtwxPJAYMMBGVsSqmaiLtv3I2gyo+34TmzZZh9Ujj4TV\nq6vft00bC3T23NOWoiJL+CX1n2a5FJF6oVEju/jusIO1ywd4b7kjIgOJd96xxFRg7feDBtnd/YgR\nlu2yeXPbFpjTozZFzlAaWFeXmjSx2oR+/cJrPdq0gS1bwgOJn3+2jpYvvmivmze3cxYIInbfPTsz\nhGbr0FmJTTUN0uDoH1X2WrHC8hsERnJ88IFdEPPzbQjonntaEDF8eO12DFy+3AKd2u7TEG9ZImtW\ntt7a+kJMnRpcfv459jHy8632IRBE7LGH9T/JNA0jTZ2aJ0RqWeQ/qoICu6tV8FC9TARbFRVWCzFl\nSnD58cdghssDDoCDD4ahQ+3C2FBVVFiNTWgQ8f331b+nT59gELHnnpmZK0MJq1KnoEGklkVrrwbr\nrPfNNwocYsmGu8JAhsspU6w54//+z6ap3mor2G8/CyAOPjj4M2zItUqLFoUHEfPnV7//9ttbLc6e\ne1oQ1rt37Qdi2fA7lesUNIjUsmg94wP0Tyu2bLwrrKiwxFSvvmrLxx9bLcSee1qw8MQTlm8ioCH/\nfDg4ooIAABYjSURBVFessO8eCCJmzw4OlY2moMD6UwwYEFz69Utv3wglrEqdggaRWhb4RzVzZnhW\nQciOC2G2yoW7wuXLrTPgpEmWQCky86J+vkFr18L06cEg4sMPYcOG6t+Tl2ejXQYMsAyhO+9szRzb\nb9+wm4cySUGDSB2JNmlSNl4Is0Um7wqTaWb49VebQ+Prr4PrttoK7r4bDj/cmqMkaNMmq7WZOtU6\no86ebecunktEo0aW2nuHHSwwCzwGnrdqVfvlj9RQmqYUNIjUIVWP5oZkazkCP9+lS+3Cts02MGOG\nTRR12GFwwgmw//51N7tnrlm7FubMgU8+sSBi9myYN88CjES0bVs1kAg8dupUO+c/F2rG0kF5GkTq\nUF2M+ZfU/fhj9a9jifbzXbQISkpg4kTrNNmpkwUPJ51kVe+1IVfvelu2tIvt8OHBdZs3W7bPefMs\npfhXX1mTz4IF4bV2oVautOWjj6pua9LEaik6drRzErm0aWN5J5o3t2Av8Lx5c0tgFSsJWLK/M1I9\nBQ0ikvUiky2lkmipa1e4+GKYMMGq4R97DB58EG6+GYYNs+Dh6KPTW6U+ZkzwrnfhQqv9yNVg9eef\n4cwzqwZA3ltgsHChLQsWhD8uWRK9qSOQhjy0GSleeXnhQURoUBGZt6Kuk3PVV2qeEJGsV9vNSBs2\nWAfKxx6DN96wkQN/+AP86U+w116pz+uQjSNPkpVstf/GjZZWPDSQCDz/7rvY82ykKi/PhpDmSu1O\notQ8ISISobabkZo2tdqFo4+2/g9PPGEBxJNPQrducOKJtnTvntzxI2tKVqywQCgXL2LJVvsXFECv\nXrZEs3598LyELr/8Ar/9ZtsDj9U9X78+fChply65W6uTjRQ0iIiE6NwZLrnEmjCmTbPg4bbb4Jpr\nYNQoa74YMyY4F0Y8Jk0KH6mzZk3uNlGks6koVPPm1nTUtWvqx9q8ORhAbNmS+vEkSJOpiohE4ZxV\nvT/yiE1Z/fjj1iZ/wgnQoQOMG2fV9PG08LZvD+3aha/L1Y55kybZeenRwx4nTcp0iapq3Bhat7aA\nZrvtMl2a+kVBg4hIDVq0sGDh3XetDX78eEsetccesNNOcOON1i5fncg78lzrmFdWZt936FB7PX26\n1ZTkYhOLJE9Bg4hIAnr0sKaKhQstcBg8GK6/3vo77LEH3H+/zYcRKRfu0KsTGAGycKE9Hnlkpksk\nmZBQ0OCcO8M5N8c5t6pymeacOzBke4Fz7j7n3E/OuTXOuWecc4pDRaTeycuDffaxfA/Ll9tjq1Zw\nzjlWizB6NPz737Bune0f6My5YEFu3qEr74FA4jUN3wMTgIGVyzvAC8653pXb7wRGA2OAEUAn4Nn0\nFFVEJDu1bAljx9qkWT/+CHfeaWmsx4614GDsWHjllapzYeSSXG9ekfRIKGjw3r/ivf8/7/03lcvl\nwFpgd+dcK+BkYLz3frL3fjZwEjDcOTck/UUXEck+7drBn/8crMq/7DJLw/z739uF9qyzbFrvXOvV\nn+vNK5IeSfdpcM7lOef+CDQHpmM1D/nA24F9vPdfAouBoSmWU0Qk53TvDpdeaimXP/kETjnFahz2\n2cdGYJxyitVORM62mo1yrXkl0HFzhx3scfnyTJeofkg4aHDO7eycWwNsBO4HjvDefwF0ADZ571dH\nvKWscpuISIPkHPTvD3/7m42y+PhjOPVUm1Fy9GibSKt9ewskhg3TBS4d1HGzdiST3OkLoD+wFdZ3\n4Qnn3Ihq9ndAjSOZx48fT+vWrcPWFRcXU1xcnEQRRUSyk3MwaJAtN95otRCjR8P339v2sjLo0wfu\nvdcm1MrEtNLZINVJvupjx82SkhJKSkrC1q1atapOy5Dy3BPOuTeBb4CngbeArUNrG5xz3wF3eO/v\nivF+zT0hIg1a5NwUTZrYRE4FBbDffnbxPPRQm/GxoUh1amtNjV070pGnIQ8oAEqBLcA+gQ3OuR2B\n7bE+DyIiEkXkSITBg60Z46abbO6Fk0+2u+xRo+COO3J3sqtEpFpToI6btSOh5gnn3A3Aa9jQy0Jg\nLDAS2N97v9o590/gdufcL8Aa4G7gA+99lFnURUQE7IIWbRbP8eNt+fFHeOEFeOklmxfjvPOgd284\n4ABbRoxIbC6MXJDqHBe1PclZQ5Von4ZtgSeAjsAq4FMsYHincvt4oBx4Bqt9+D/gz+kpqohI/VTT\nBa5jRzjjDFvWroU334SXX4b//c9yQhQUwJ57BoOInXe2vhO5LFogJZmXcp+GlAugPg0iIknxHubP\nhzfegNdfh8mTbZrojh1h//0tgNh336qTZUn9Udd9GjQ1tohIjnLORlr06QN//Sts2GA1Fq+/boHE\n44/bPkVFFkDsv79NONWkSaZLLrlKE1aJiNQTTZtazcKtt8KcOfDDD/DYY9CrFzz0EOy1F7RtayMx\n7rsPvvkmvqm9RQJU0yAiUk917AgnnmhLRYVlpXz9dVv++ldLZd29e7AWYu+9ISJdjkgYBQ0iIg1A\nXp41UxQV2QiMNWvgvfeC/SEeeAAaNbLmi733tjwHQ4faZFwiAQoaREQaoMJCOOQQWwC+/TYYQNx3\nH1x7rQURu+5qAcSee9rjtttmttySWQoaRESE7t3h9NNt8R6++MI6VU6dCi++CHdV5vT93e8seAgE\nEj175v7wTomfggYREQnjnCWP6t3bJtYCWLrUgohAIPGvf1lwsc02lsEysAwaZBNvSf2koEFERGrU\nuTMcc4wtAKtWwbRpMH26zdp5332wcqVt69IlPJAYOBC23jpzZZf0UdAgIiIJa90aDjrIFrBah0WL\nLICYOdMeb74ZVldOX9izZ3htRFERtGiRufJLchQ0iIhIypyDbt1sOeooW1dRAV9/bQFEYHnuOUtC\nlZdnSakCQcTgwbDLLpYSW7KXggYREakVeXmWWKpXLzjuOFu3ZQt89ll4jcSTT9r6xo0tcAht2ujd\nG/J1pcoa+lGIiEidyc+H/v1tGTfO1m3YYBksA0HElCnw4IPW5FFYCCNH2rTge+9tQUWechlnjIIG\nERHJqKZNYbfdbAlYuxZmzbLRGu++C5ddZsFFmzbBAGLUKNhpJw35rEsKGkREJOu0bAkjRthy6aWw\ncSN8+CG8844tf/0rbN5swzv33ju4dO+e6ZLXbwoaREQk6xUUWDPFyJFwzTWwbh188EEwiPjPf6zj\nZdeuwQBi1CgbKirpo6BBRERyTosWNsnW/vvb619/taRTgSDiscds/eTJVlsh6aGgQUREct5WW4XP\npbFihU3INWhQRotV7yhoEBGReqddu2C+CEkfDVwRERGRuChoEBERkbgoaBAREZG4KGgQERGRuCho\nEBERkbgoaBAREZG4KGgQERGRuChoEBERkbgoaBAREZG4KGgQERGRuChoEBERkbgoaBAREZG4KGgQ\nERGRuChoEBERkbgoaBAREZG4JBQ0OOcucc595Jxb7Zwrc84955zbMWKf95xzFSFLuXPu/vQWu/4q\nKSnJdBGygs6D0XkI0rkwOg9BOhd1L9Gahj2Be4DdgH2BxsAbzrlmIft44CFgW6AD0BG4KPWiNgz6\nIzA6D0bnIUjnwug8BOlc1L38RHb23h8c+to59ydgOTAQeD9k03rv/YqUSyciIiJZI9U+DVthNQs/\nR6wf65xb4Zyb65y7MaImQkRERHJQQjUNoZxzDrgTeN97/3nIpqeARcAPwC7ALcCOwB9SKKeIiIhk\nWNJBA3A/0AcYHrrSe/9IyMvPnHPLgLecc929999GOU5TgPnz56dQlPpj1apVzJo1K9PFyDidB6Pz\nEKRzYXQegnQuwq6dTevi85z3PvE3OXcvcAiwp/d+cQ37NgfWAgd479+Msv1YrHZCREREkjPWe//v\n2v6QhGsaKgOGw4CRNQUMlQZg/R5+jLH9dWAs8B2wIdHyiIiINGBNgW7YtbTWJVTTUJlvoRg4FPgq\nZNMq7/0G51wP4FjgVWAl0B+4HVjsvd87baUWERGROpdo0FCB1RpEOsl7/4RzrgswEegLtAC+ByYB\nN3jv16ahvCIiIpIhSfVpEBERkYZHc0+IiIhIXBQ0iIiISFzSEjQ45/Z0zr3onFtaOUnVoVH2udY5\n94Nzbr1z7k3nXM+I7Vs7555yzq1yzv3inHvEOdciYp9dnHNTnHO/OecWOecuTEf506mmc+Gceyxi\nQq8K59yrEfvk/LmIc3KzAufcfc65n5xza5xzzzjn2kfss51z7hXn3Drn3DLn3C3OubyIffZyzpU6\n5zY4575yzp1YF98xHuma5K0enIcznHNzKn+nVznnpjnnDgzZXu9/FwLiOBf1/vchmsq/lQrn3O0h\n6xrM70VAjPOQPb8T3vuUF+BA4FrgcKAcODRi+wQs1fQhwM7A88ACoEnIPq8Bs4BBwDBsdMbEkO2F\n2LDNx4HewNHAOmBcOr5DupY4zsVjwCtAO6B95dI6Yp+cPxfYCJrjK8vXD3gZG1bbLGSff1SuG4kN\nzZ0GTA3ZngfMxYYS9QMOwOY6uT5kn25YHpBbgF7An4HNwH6ZPgcJnId3gQcifida1rPzMLryb6Nn\n5XI9sBHo3VB+FxI4F/X+9yHKORkMLARmA7eHrG8wvxc1nIes+Z2ojS9dQdUL5Q/A+JDXrYDfgKMr\nX/eufN+AkH0OALYAHSpf/7/27jbGjqqO4/j3V5BiqW2BPqABammlJlJKKU8NpS2tqdEoaHyBLxQk\nJhp9IbxAEoNREmxI8SGQxickGtBo1Eg0JqXVUreNtkJoFYrYWtJqjdBHa1ss0FL+vjjnsrPTe7ez\nZbf37tzfJ9ns3pkzszP/+9+9/3vmzD2fBfYCpxfa3As81+4neoCx+CHwaD/bvLumsRifz2tuIQde\nBT5SaDM9t7kqP35/TurxhTafAfY3zh1YCjxT+l0/BZa3+5yrxCEv+33xH0STbWoXh3x8+4BbuzUX\nmsWiG/MBGA1sARYWz73b8qJVHDotJ4Z8TIOkKaQpsh9vLIuIg8ATwJy86Bpgf0T8ubDpKtLtnVcX\n2qyNiNcKbVYC0yWNHaLDHyoLclf1ZknflnROYd0c6hmL8uRms0kfLlbMiy3ADvrmxaaI2FvYz0pg\nLOm23kabVaXftbKwj05zMpO81SoOkkZI+hgwClhP9+ZCORbrCqu6Jh+AbwG/iYjVpeVX0F150SoO\nDR2RE29m7omqziP9k9xVWr4rr2u02V1cGRHHJP2n1GZbk3001h0YrAMeYo8BvwS2A1NJPQTLJc2J\nVPrVLhZS08nNzgOO5AKyqJwXzfKmse7pftqMkTQyIl4dhFMYFC3iACee5K0WcZB0CalIOBM4RHoH\nuVnSLLovF5rFYkte3RX5AJALplmkwrFsEl2SFyeIA3RQTpyKoqEV0fyDogbSRvn7sPmwiYj4eeHh\nXyVtIo3vWEDqgmplOMeiMbnZ3Aptq+QFJ2jTqbEYyCRvj6v1JG99Nu9nXafFYTPpU2LHAR8FHpE0\nr5/2dc6FprGIiM3dkg9KHwZ4P+ma+tGBbEqN8qJKHDopJ07FLZc7SQc2qbR8Ir1Vz878+A2STgPO\nzusabZrtA46vnoaN/ITvJQ2IgprFQmmukg8ACyLihcKqncAZksaUNinnRfk8JxXWtWozETgYEUfe\nzLEPplIcWs3D0vBE/l7MiWEfh4h4LSK2RcTGiLiL9O7nNrosF6DfWDRTy3wgvaueAGyQdFTSUdKA\nx9skHSE99yO7IC/6jUPuoSxrW04MedGQXxR3Aosay3ISXE3vNbz1wLjcTdmwiFRsPFloMy+/gDYs\nBrZEREd1xw9ErjLPpXdCr9rEQr2Tm10fx09utoE0uLOYFxcDF9I3L2ZIGl/YbjHp8svfCm0W0dfi\nvLwjnCAOzZQneatFHJoYAYyki3KhH41YNFPXfFhFGul/GanXZSbwFGkqgsbPR6l/XvQbh3zZuqx9\nOTFIoz7Pyid6GWlk6+358QV5/Z2k0cEfysH5FbCVvrdcLs+BupLUfbsF+FFh/RjS9ZyHSV28N5Fu\nH/nUYJzDYH31F4u87j5SwTQ5P4FP5Sf1LXWKBakrfj9wHam6bXydWWqznXRpZjbwR46/nepp0jiQ\nS0l3kewC7im0eWc+96WkkdWfA44A7213DKrEAbgI+BJwec6JG4DngdU1i8MS0uWpyaTbru8lFQoL\nuyUXqsSiW/Khn9iU7xromrxoFYdOy4nBOsH5pBfIY6WvHxTa3E16oTtMGrE5rbSPcaTK6gDpn+z3\ngVGlNjOANXkfO4A72v3kDiQWpEFPK0g9L6+QBjN+B5hQt1i0iMEx4OZCm5HAMtLlmUPAL4CJpf1c\nQPpsg5fyH8FSYESTmG8g3ca7FfhEu8+/ahyA84EeYE9+LreQXkRGl/Yz3OPwUM73l3P+/5ZcMHRL\nLlSJRbfkQz+xWU3foqFr8qJVHDotJzxhlZmZmVXiuSfMzMysEhcNZmZmVomLBjMzM6vERYOZmZlV\n4qLBzMzMKnHRYGZmZpW4aDAzM7NKXDSYmZlZJS4azMzMrBIXDWbWlKTtkj4/gPaTJb0u6dKhPC4z\nax8XDWbWyhXAgwPcpt/PpZd0i6T9J39IZtZOp7f7AMysM0XEvpPYTBXWe8Ibs2HKPQ1mNSHpg8V3\n8ZJm5ssFSwrLHpL0cP55rqS1kg5L+qekBySNKrTtc3lC0nRJf5D0sqRnJS3K+7+hdChTJa2W9D9J\nf5F0Td5+Pmm217F5u2OSvjxE4TCzIeCiwaw+1gKjJc3Kj+eTptNdUGgzD+iRdBHwGGmq4UuAm4Br\nSdMQH0eSgF+Tpie+Evg0sITmvQZfBe4DZgJ/B34iaQSwDrgdOAhMAt4OfP3kTtXM2sFFg1lNRMRB\n4Bl6i4QFwDeByyWNkvQOYCqwBvgi8OOIWBYR2yLiT6QX9FskndFk9+8DpgA3R8SzEbEOuIvmlyO+\nFhErIuJ54CvAZGBaRBwFDqRDjT0RsTsiDg/O2ZvZqeCiwaxeeugtGq4DHgU2k3oR5gMvRMQ2Ui/A\nJyUdanwBK/J2U5rs92LgXxGxp7DsyRbHsKnw84ukwmLiwE/FzDqNB0Ka1csa4FZJM4EjEbFV0hrg\neuAcUlEBMBr4HvAAx/cW7Giy34EMYDxa+Lmxjd+gmNWAiwazelkLjCFdaujJy3qAO4GzgW/kZRuB\n90TE9or73QxcKGlCobfhqibtTlRYHAFOq/g7zazDuPo3q5GI+C/p8sDH6S0a1gCzSZcYGsuWAnMk\nLct3WUyTdKOkpgMhgd8B24BHJM2QdC1pwGPQt1A40S2X/yAN1lwo6VxJbx3QCZpZW7loMKufHtLf\ndg9AROwHngNezIMTiYhNpDEO7yL1TmwE7gb+XdjPG8VARLwO3AicRRrL8CBwD6lIeKXZNi32sx74\nLvAzYDfwhZM8RzNrA0X4c1bMbOByb8Na0p0RVS9zmNkw5qLBzCqR9GHgJWArqYfifmBfRMxv64GZ\n2SnjgZBmVtXbSB/adD6wlzTO4Y62HpGZnVLuaTAzM7NKPBDSzMzMKnHRYGZmZpW4aDAzM7NKXDSY\nmZlZJS4azMzMrBIXDWZmZlaJiwYzMzOrxEWDmZmZVfJ/YdMqfqAKjPkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xbc66710>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot smooth components\n",
    "res_bs.plot_partial(0, cpr=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgkAAAF5CAYAAAAVqLmkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XmcU9X9//HXmRkYQHYQQUEUEKRYF0YURHCtVK0rWsUF\ntbWC1bZqf99at2rrUreqtYpba13QcQP9qrXiLi4gAoqigBuiyA7Koqwz5/fHZ/LNMjczk0ySm0ze\nz8cjj0xys5xcQu7nnvM5n+O894iIiIgkKgm7ASIiIpKfFCSIiIhIIAUJIiIiEkhBgoiIiARSkCAi\nIiKBFCSIiIhIIAUJIiIiEkhBgoiIiARSkCAiIiKBFCSIiIhIoJSDBOfcts65B51zK5xzPzjnZjnn\nBiY85i/OuUU12190zvXJXJNFREQkF1IKEpxz7YG3gI3ACKA/8Hvg25jHXAicC4wB9gK+ByY555pn\nqM0iIiKSAy6VBZ6cc9cCQ7z3+9XxmEXADd77m2tutwWWAqd57x9rZHtFREQkR1IdbjgCmO6ce8w5\nt9Q5N9M5d2Zko3NuR6Ar8HLkPu/9GuAdYEgmGiwiIiK5kWqQ0As4G5gHHALcCdzqnDulZntXwGM9\nB7GW1mwTERGRAlGW4uNLgGne+8tqbs9yzg3AAofxdTzPYcFD7Q3OdcLyG74ENqTYHhERkWLWAtgB\nmOS9X5npF081SFgMzEm4bw5wbM3fS7CAYBviexO6AO8lec0RwEMptkNERESiTgYezvSLphokvAX0\nS7ivH7AAwHs/3zm3BDgI+AD+L3Fxb+D2JK/5JcD48ePp379/is0pbueffz4333xz2M0oKNpn6dF+\nS532WXq031IzZ84cTjnlFKg5lmZaqkHCzcBbzrmLgMewg/+ZwK9iHnMLcKlz7jOs0VcCC4H/TfKa\nGwD69+/PwIEDkzxEgrRr1077LEXaZ+nRfkud9ll6tN/SlpXh+pSCBO/9dOfcMcC1wGXAfOB33vtH\nYh5zvXOuFXAX0B54AzjUe78pc80WERGRbEu1JwHv/XPAc/U85grgivSaJCIiIvlAazeIiIhIIAUJ\nBWzUqFFhN6HgaJ+lR/stddpn6dF+yy8plWXOSgNscagZM2bMULKKiIhICmbOnElFRQVAhfd+ZqZf\nXz0JIiIiEkhBgoiIiARSkCAiIiKBFCSIiIhIIAUJIiIiEkhBgoiIiARSkCAiIiKBFCSIiIhIIAUJ\nIiIiEkhBgoiIiARSkCAiIiKBFCSIiIhIIAUJIiIiEkhBgoiIiARSkCAiIiKBFCSIiIhIIAUJIiIi\nEkhBgoiIiARSkCAiIiKBFCSIiIhIIAUJIiIiEkhBgoiIiARSkCAiIiKBFCSIiIhIIAUJUq+lS2Hf\nfaF3b7tetizsFomISC4oSJB6jRwJb70FX3xh18ceG3aLREQkFxQkSL0WL677toiINE0KEqRe3brV\nfVtERJqmsrAbIPlv4kQbYli82AKEiRPDbpGIiOSCggSpV5cu8OabYbdCRERyTcMNIiIiEkhBgoiI\niARSkFBANmwIuwUiIlJMFCQUiE8+gV694O23w26JiIgUCwUJBaJ3b9hhBxg9Gr7/PuzWiIhIMVCQ\nUCBKS+H++20a4h/+EHZrRESkGChIKCA77QTXXw/jxsGLL4bdGhERaeoUJBSYs8+Ggw+GM86A774L\nuzUiItKUpRQkOOcud85VJ1w+jtn+WsK2KufcuMw3u3iVlMC998LatfDb34bdGhERacrS6UmYDWwD\ndK257BuzzQN3x2zvBmgEPcN69IBbb4UHH4Qnnwy7NSIi0lSlU5Z5i/d+eR3bf6hnu2TA6NEWIIwZ\nA0OHWulkERGRTEqnJ2En59w3zrnPnXPjnXM9Eraf7Jxb7pz70Dl3jXOuZSYaKvGcg7vvBu8tUPA+\n7BaJiEhTk2qQMBU4HRgBjAV2BN5wzm1Vs/0h4BRgf+Aa4FTgwUw0VGrr0gXuugueegoeeCDs1oiI\nSFPjfCNOQZ1z7YAFwPne+38HbD8AeAno472fn+Q1BgIzhg8fTrt27eK2jRo1ilGjRqXdvmJx+ukw\nYQLMnGnTJEVEpOmprKyksrIy7r7Vq1czefJkgArv/cxMv2ejggQA59w04EXv/SUB21oB64AR3vvA\nmf2RIGHGjBkMHDiwUW0pVmvXwsCB0K6dlW1u3jzsFomISC7MnDmTiooKyFKQ0Kg6Cc651kBvYHGS\nh+yBzXhItl0yoE0beOQR+OADuKRWqCYiIpKeVOsk3OCcG+6c6+mc2wd4EtgCVDrnejnnLnXODazZ\nfiRwP/C69352FtouMSoq4K9/hRtvhEmTwm6NiIg0Ban2JHQHHgbmAo8Ay4HB3vuVwCbgYGASMAe4\nAXgcODJjrZU6nX8+jBhh0yOXLg27NSIiUuhSqpPgvU+aRei9X4jNapCQlJTYIlC77gqnnQbPPWf3\niYiIpEOHkCZmm21sOuSkSTb0ICIiki4FCU3QiBHwxz/CxRfDW2+F3RoRESlUChKaqCuvhCFD4IQT\nYNmysFsjIiKFSEFCE1VWBo8+Clu2wIkn2rWIiEgqFCQ0Ydtua4HC5Mlw6aVht0ZERAqNgoQmbr/9\n4Npr4brrbI0HERGRhlKQUAR+/3sYOdKmRX76aditERGRQqEgoQg4B/feC926wbHHwvffh90iEREp\nBAoSikTbtrZS5BdfwJgx0Mh1vUREpAgoSCgiAwbAv/4FDz0E48aF3RoREcl3KZVllsJ34okwdSqc\nd54FDfvvH3aLREQkX6knoQjdeKPNejjuOBt+EBERCaIgoQiVlcFjj0GHDnDkkbBmTdgtEhGRfKQg\noUh17AhPPw1ffw0nnwxVVWG3SERE8o2ChCLWvz888ogtKX3JJWG3RkRE8o2ChCJ36KFw/fVWkXH8\n+LBbIyIi+USzG4QLLoDZs+HMM2GnnWDvvcNukYiI5AP1JAjOwZ13QkUFHH00LFgQdotERCQfKEgQ\nAMrLYeJEaNnShiBWrQq7RSIiEjYFCfJ/ttkGnn8eli2Do46CDRvs/qVLYd99oXdvu162LNx2iohI\nbihIkDh9+8Izz8D06XDKKTY1cuRIeOstK7z01lu2SJSIiDR9ChKkliFDbGrkk09aUuPixfHbE2+L\niEjTpCBBAh11FNx+O9x6K1RXx2/r1i2cNomISG4pSCgQ3sO//w1ffpm79xw7Fi6+2N6zb1/o1QsG\nDYLNm5WfICJSDBQkFIDly22NhV/8As44o/aZfTZddRWMHg3z58M990Dz5jBtmvITRESKgYKEAlBe\nDh9+aH+/9hrcdlvu3ts5Cw723x+OOaZ2DQXlJ4iINF0KEgpA27Y21BBx4YUwb17u3v/bb2HtWpsS\n+c038duUnyAi0nQpSCgQBxwAv/2t/b1hA5x2GmzZkpv3HjkSpk6FTZssN8I56NEDhg61AkwiItI0\nKUgoIH/9qyUQArzzji3MlAuJQwqlpXaprIQuXXLTBhERyT0FCQWkVSu4/34oqflXu+IKmDUr+++b\nOKSw++52feCBsGhR9t9fRETCoSChwAwebDkJYFMRR4+2YYBsmjjRhhYiUyBLSqwS41dfwX77Wdlm\nERFpehQkFKDLL4ddd7W/P/gA/vzn7L5fly7w5pvw+efRKZBff23ByYIFNvMhMaFRREQKn4KEAhG7\nyNJBB8Ett0CzZrbt2mstRyEXEvMTunaF77+H4cPhs89y0wYREckNBQkFInGRpcsusx4FsOJKw4bB\nPvtkvwJiYn7C9tvDG29YwLLPPtbLICIiTYOChAIRtMjShRdC69Z2e/NmmDIl+xUQY/MTIlMge/a0\nwGWnnWyq5rPPZrcNIiKSGwoSCkTiGXy3blBWBh06xN//xRfZbUdsfsKbb0anQHbqBC+9BIccYotD\n3XNPdtshIiLZpyChQASdwYN198datQrWrct9+wBatoQnnoCzz4azzoJLL83tOhMiIpJZChIKRLIz\n+IkTLRegRQu7vXEj/M//hNfO0lL4xz+s0NM119h6D2vWhNceERFJn4KEAteli+UDzJ5txZYA7rzT\nhiPCWsrZOQtUnnnGFqQaMkQzH0RECpGChCaid2+44Ybo7SVLwl/K+fDDbc2HzZthzz1tKEJERApH\nSkGCc+5y51x1wuXjmO3lzrnbnXMrnHNrnXNPOOdU3T9Hxo6NDjtEZGIp59gaDan2TvTvD+++awmN\nxx8P55xjC1SJiEj+S6cnYTawDdC15rJvzLZbgMOBkcBwYFtgQiPbKHWIPYAPHw79+sVvLy9v/Hsk\n1mhItXeiXTt49FEbBvnXv6y09CefNL5dIiKSXekECVu898u998tqLqsAnHNtgV8A53vvX/fevwec\nAQx1zu2VwTZLjMQDePPm0KdPdPuKFbByZePeI6hGQ6qcgzFjrDLkhg0wcCCMH9+4domISHalEyTs\n5Jz7xjn3uXNuvHOuR839FUAZ8HLkgd77ecBXwJDGN1WCJB6wV660s/TDD7fby5fDuec27j2CajSk\na7fdYPp064049VT4xS+srLOIiOSfVIOEqcDpwAhgLLAjMNk5txU29LDJe5844W1pzTbJgqADuHNW\nzChSaOmRR+Dxx9N/j2Q1GtLVujU88ADcd58NQwwaZLMzREQkv6QUJHjvJ3nvJ3jvZ3vvXwQOAzoA\nP6/jaQ7wjWij1CHZAbxbN7jttujjzj47/SWdk9VoaKzTTrNehdJSCxTuuQe8vikiInnD+Ub+Kjvn\npgEvAi/VXDrE9iY4574Ebvbe/z3J8wcCM4YPH067du3ito0aNYpRo0Y1qn3FzHubUTChJnX06KMt\niHAu3HYlWr8ezjsP7r7b2nvXXbXLTYuIFLvKykoqKyvj7lu9ejWTJ08GqPDez8z0ezYqSHDOtQYW\nAH8CHgSWAyd675+s2d4XmAsM9t4Hrg8YCRJmzJjBwIED026LBFu+HAYMsGuwbv5TTw23Tck8/riV\nc27bFh56yGZtNFVLl1rS6eLF1uszcWLmemhEpHjMnDmTiooKyFKQkGqdhBucc8Odcz2dc/sATwJb\ngEdqeg/+BdzknNvfOVcB/Bt4K1mAINm39dY29TDiN7+Bb74Jrz11Of54mDXL1qPYbz+44grYsiXs\nVmVHY6eViojkQqqJi92Bh7HegUewnoPB3vvIJLvzgWeBJ4DXgEVYzQQJ0bHHwskn29+rV8OZZ+bv\n2P/228Orr8Lll8OVV8L++8OCBWG3KvMyMa1URCTbUk1cHOW97+69b+m93957f5L3fn7M9o3e+994\n7zt779t474/33oeweoAk+sc/ojMhnn8e/vnPcNtTl7Iy+NOfYPJkWLjQpk1OaGIluTI5rVREJFu0\ndkOR6NAhPjC44AL48svQmtMgQ4fC++9bSefjjrN6D02lpHOmp5WKiGSDgoQicthh8Mtf2t/r1lkh\no+rqcNtUn/btrZbCuHEW5OyzD3z6aditarxsTSsVEckkBQlF5qabbNwfbOy/T5/0Fm7KJeeszsPU\nqRbcVFRYgSgREckuBQlFpm1buPfe6O3586MZ9kcckfx5jVkJMlO6dYNOnaz3Y9QoGD3aaiyIiEh2\nKEgoQgcdZEs2J5o1K/lz8mHK3siR1psQWeth/HjYe2+YOzf3bRERKQYKEorUddc1/LFLl1r55Fhh\nTNlLfM/ttoPNm2HPPeHBB3PfHhGRpk5BQpHaaivo3z/+vtglpmONHAkbN8bfF8aUvcT37NnTgpfj\njrOhhzPO0IqSIiKZpCChiN1xR/w6DmvXBo/xJ57Bl5eHM2UvaNrgVlvZapL33QePPQZ77QUffZT7\ntomINEUKEorYEUfEV1786ivo2rV2YmLiGXzz5tlpT33JkXVNGzztNHj3XSgpsRUl//Wv/K0qKSJS\nKBQkNFENmY0Q1GuwZk3txMSJE6F16+jttWvrngmRrsYmR/7oR/DOO1aC+swzbSGr777LfDtFRIqF\ngoQmqiEH3JYtkz//66+jf3fpYgmCseqaCZGuTKxn0KoV3HOPrSL59NOWd/HYY+pVEBFJh4KEJqoh\nB9wpU6BNG1srobQ0ftvq1VBVFb2deJDduNFqLmRy/D+T6xmcdBJ8/DEMGQInnACHH241IUREpOEU\nJDRRDTngDhhgwwubN8OiRZb0FwkWVq+Gq66KPrasrPbz166FwYMz1+ZMr2fQvbu9xlNPwYcf2ue9\n4YbavSIiIhJMQUIBWbGi4Y9N9YDbpYuN57/6qiX/Afz5z3YbYMuW4OetW9fwNtUnW+sZHHWU9SqM\nHQt//KPVVZg8OTOvLSLSlClIKBBz5sAOO8BzzzXs8ekecIcNg7/8xf723rrtZ8+GTZuCHx87hTKf\ntWlj61a8+65N4dxvP1vwavJk5SuIiCSjIKFA7LwzDB9uWfurVmXmNZPNgLjoIvjJT+zvJUsscEhm\nt90a/rr5YOBAK+38yCOwYIEFC4MHwxNPxOdgiIiIgoSC4Zwtlbx+PfzmN5l5zWQzIEpKrMxx1652\nO2gaYfPmNowxaVLDXzdflJRYMuOHH8J//mMFmY4/Hvr2hdtvV9VGEZEIBQkFZNtt4bbb4OGHYcKE\nxr9eXTMgttnGphEGDSe0aWNTJJMNY2RiKmMulJTYkMMrr9gwxF57we9+Bz16wFlnwYsvJs/FEBEp\nBgoSCsxJJ9mZ+dixje/Gr28GxIEHwmWXRW+XltqB9LPP6s5xyORUxlzZc0+orLTPNmYMvPwyHHKI\ntf2ss+CllxQwiEjxUZBQYJyLrrkwZkzjku4aMgPiT3+Cffaxv6uq4NNPobq68a+br3bYAf76VwsW\npk+HX/7SAoSf/MQChjFjbIjlhx/CbqmISPY5H3Jqt3NuIDBjxowZDBw4MNS2FJKJE23s/8EH4ZRT\n6n7s0qX22MWL7UA3cWJq0wv32su64yN69oQvv0yr2QXJe5g50yo3Pv64FWVq3tyCp4MOssugQcG1\nJEREsmnmzJlUVFQAVHjvZ2b69dWTUKCOPdbWKDj3XFi4sO7HNjaRcOnS+NsLFsD48am9RiFzDioq\n4LrrbErphx/C9ddbbsb111uw0LGjrWdx881Wsrq+3hYRkUKgIKGA/eMflpl/5pl1DzskJg5On57a\n9MSgKZenn261G4qNc7DLLpbg+PTTtm+mTIELL7RZERddBLvvbj01xx1nsyXmzFEtBhEpTAoSCliH\nDjYtctIkWxo5mcTEwY0bU+tVCCqkVFVllQyLfZXFsjKrs3DJJTZL4rvv7PrXv7YaE+edZ6tTbrut\nJZ3ec4/1RihoEJFCoCChwB16qJ3V/8//1B4WiIhNJCwvj9/WkOmJydY6+PRTGDVKRYhitWgBBxxg\nVSvffNOChkmT4LTTLDgYOxb69LG8jtNPt2GbfCo2JSISS0FCE3DDDTY98YILgrfHlmjec8/4bQ2Z\nnhiUkBdZCOr5562LXYJttZVNpbz2WlsbY9UqeOYZG4p4/3049VSrSVFRARdfDK+/nrwEtohIrilI\naAI6d4a//c2KLL3wgt2XrDRyOtMTE3sfwAKPSPBwww3xwx2zZ9sy0s2aZX456ULXrh387Ge2jsT7\n79uQxIMPQv/+NnS0//7QqZMN5YwbZ8NCIiJh0RTIJsJ7K3701Vd2kP7JTyznIGLoUOtNiNXQqZE9\netSeQTF0qM2u+PWv7XZpqZU4HjHCAoO1a6OPbdPGlqSWulVXW+AwaZJd3nrLCjj16WP7dcQIG8po\n3TrslopIYyxbBvfdZ4Xa2rdv3GtpCqQ0iHNw5512ML/qqoaVRm7o1MjVq2vft3ChlW0+6yy7XVVl\nz58ypXahoXwrPJTJBagy+VolJbYA1UUXwWuvwcqV8NRTFvD9979w5JE21fLAA2065vvvKwFSpFBU\nV1tS8wknQPfuVqjunXfCblX9FCQ0If362bj29ddbt3asoNyDhq6xEHTmumCBBRazZ0eDix9+sLUQ\nEhMZq6rya0XITC5Alc3FrNq2jQ47fP65JYrefLPlOVx5Jeyxh/27jh5tAdvy5Zl7bxHJjOXLbUi2\nXz8rvPbBBxbkf/ON9Q7mO9WIa2L++EfLTSgrsyI/S5ZEhxISdesWP+adLImxrgP7kiUWHR9+uK13\nkGxK5Bdf2OXYY2sPe+RaJhegyuViVn362OWcc2wa69tvR4cmHnzQHlNRER2aGDLE8kJEJLe8t97A\nu++2317nLFn5X/+CYcOCF87LV+pJaGLKy+Guu6yM8i9/aWegyVZrzMQaC9262Xs+9RTsvXf9j8/W\nQTSVbv9MLkAV1mJW5eWWn3DttfDee7Zf77/fzlbuvhv2288SII8+2tb6UAKkSPatWGFJ5DvvbMOC\n771na8F8841Ndx4+vLACBFDiYpN16qk2jj1vnh0sGiMxEbG01Ob5JyY7rlplB6fZs6OPLSmJL1Ec\nlECZCYnrSwwaBNOmBT922TLr0Uh3LYtsvVamVFfbj9Pzz1svw5QplgC5007RXob991cCpEgmeA+T\nJ9vJ2YQJdt/IkbYYXC6CgmwnLipIaKKWLLFo9uc/tzPLxnj9dYuKIwf7XXaxoYWgg+HixXYmHzlz\n7dfP8iNWrMjuQbRFC+uCjygvhw0bMv8+hWjNGhsSigxNzJ9vwxDDhkWDhl13LbwzHJEwrVxpvXd3\n320nY337WiL3aafZtPRc0ewGSUvXrjbL4Z//TH5G3VDnnx/fGxCbrJioWzd48UXL3gX7z7N2rSX1\nJRv2iMjkTAGJats2Ouzw+efwySfWJdqyJfz5z7bWRNeulnV95532b6ZZEyK1eQ9vvGEr7263neWA\n7b67BeFz58Lvf5/bACEX1JPQhG3ZEl3CeOrUaJXEVJWU1D5o9OplB5xkPv/cxsy//tpu9+tnvQ/b\nbZf8OfvuW39th2T23js+GNprr8KYXhS2jRttn7/8Mrz6qg3ZbNlia00ccED0suOO6mmQ4rVqFTzw\ngPUazJljCcRnnWWl1bfeOty2ZbsnQbMbmrCyMluFcOhQ61EYMya91wmKI+tL0Ovd24YpDjjApkvO\nm2ezLSZNsmGQWJGiTokH9VSSHJ95pnZugNSvvNyGkg480G6vW2eB2SuvWNBQWWm9SNtvb4+JBA09\neoTbbpFs894C6Lvugscft/8HxxwDt91mOT0lRdIPr56EIvCLX9jsg08+Sa8rrE0bO3hElJTYwbgh\nuQULFtjBJZKj0KkTPPecnelHJPYgRGQryVEa7rvvrHs1EjTMmmX39+4dHzR07RpuO0Uy5dtvo70G\nH39s3/VIr0HYSclBlJMgjXbttRYVp7sQ09SpFiiUldn1Bx80/D9Lz54WAOy+u91eudIOLpMmRR8T\nGZKIKCtLf0qmZFb79nDEEVbE6f33LQF1wgT46U/t3/Wkk6znpn9/OPtsK+r01Vdht1okNZFeg9NO\ns6G2//f/YMAAy6/65BP4wx/yM0DIBfUkFIlx46wIz9SpDatnkGmrV1vy3Guv2e2yMssMPukkrfVQ\nyJYutX/TV1+1aWBz5tj9PXrY7Il997XLgAHF0z0rheO776wQ2d13W0J2r17wq1/BGWfY6qyFQFMg\nJSOqqiyJsaTExv7TTWJs6KJQQTZssEWhYnsIbrnFzlIXLIje17MnfPlleu2TcK1YYWdkb7xhQ0Uz\nZlgiZPv21js0dKgFqRUVtUuHi+SC93aydNdd8NhjsHmzlT8fM8bKJhdaMJvXww3OuYucc9XOuZti\n7nut5r7Ipco5N67xTZXGKC21JMYZM+Cee9J/nSOPjF+r4IgjGv7cFi3sP2VsAuV559lBJFZk+mSh\n0NTNqM6d7Qf3xhvth3j1astnuOAC2LTJqs8ddBB06GBDFKedZt/LadPi61yIZNp331nS4W67WRL1\n66/DZZfZcOcTT9hCaoUWIORC2j0JzrlBwKPAauBV7/0FNfe/CswDLgMik6Z+8N6vS/I66knIoV/+\nEp580mYbpDN1JxNFi7yHK66Av/wlel/HjjbssN12+VG1MBWNmbpZbKqr7bs3bZpd3n3Xch02b7YC\nT7vvbj1ee+1ll3799MMt6fPevmd33QWPPGKB6pFH2olKUwkK8nIKpHOuNTAeOBMLBhL94L3XmnR5\n6Npr7SB80UU2LTIMzlkRn223tTyJqiqbh7xuXfhzjtORy0WeCl1JifUgRHoRwILOWbMsYJg2zXoe\nxtX0PbZpA3vuGQ0aBg2ynibVbJBE3ltu04oVdpk+3YKDDz6wIcxLLrGZXrlaX6WpSLdOwu3AM977\nV5xzQUHCyc65U4ElwDPAld779ek2UjJn663hmmvg17+GM8+EwYMb9rxILkJix9Ouu6bfljFjrEhP\nZHnpTZvsQLHjjjalLl/WQqhPQ1fTlGDl5dEg4Jxz7L7Vq21oLNLb8NBDtrwu2HcjEjBErjt0CK/9\nkh0bN0YP+JHL8uV13960Kfr80lIbDr3uOus1SDcPq9ilPNzgnDsRuBjr2thcM7zwXsxww5nAAmAR\nsCtwPfCO9/64JK+n4YYcq6qyH1fv7Qe4If95ErvUy8vtDC8TB/EePWDhwuBtgwZB8+b5tYBSonxc\n5KkpWrw42tsQCR4iS5PvsIOtKTJgQPR6552t9LSEr7ra6g/Ud5CPvR074ymiWTPLe9l6a7uOXBJv\nd+5sBcA6dsz9Z821vJrd4JzrDkwHfuK9/7DmvrggIeA5BwAvAX289/MDtitICME771gvwu23W69C\nfXr3jj9brq8sc31iZ0ksXx78gwAWjMTmQDR0vL8xszDSFcZ7FjPv4bPPLGCYNcumsH30UbROQ0mJ\nfW8Tg4e+fS3wlPR4D99/37Az+8jtVavi13+J6NixYQf8yH1t2mioKVG+BQlHAROBKqJJiaWAr7mv\n3Ce8oHOuFbAOGOG9fzHgNQcCM4YPH067hDlRo0aNYtSoUQ3/NJKSX/3Ksnrnzav/YJbp5LzEpZ1L\nSuySONMhcanp8nJLbqzvIJzK0tGZogTG/LBmjVXKiwQNketIrkhZmdXe335768Xq0cPyHGKv27QJ\n9zPk0qZNVuQslbP8oJkoW21V9wE+8XaHDvZvIQ1XWVlJZWVl3H2rV69m8uTJkCdBwlZAz4S77wPm\nANd67+dmgeedAAAgAElEQVQEPGcoMBnYzXs/O2C7ehJCsmKFnVUdfTTce2/dj810l3riLIlYzZvH\njy0mU9dBONksjNmzbfrT+vXWFT1lip1dZkKme1sks1atigYNc+fa1LeFC+166dL4fJu2beODhqBA\nonXr8D5LMtXVNgSTyll+UOGysrLUDvidOmloJyx5NbvBe/898HHsfc6574GV3vs5zrlewEnAc8BK\nYDfgJuD1oABBwtW5s81bHzvWehWGDAl+XGw3eqdONl1tyJDsdal362azME4/ve658+nMIthnn+jQ\nxtq19jkyVd2xqScwFvpwSseOVgVy2LDa2zZtgkWL4gOHyPX779sCYkuXxj+nXbvg4CHyd9eu1iOR\nbsKc9/DDDw1P2lu+3HoEgrr1O3SIP8APGFB3N3/bturWF9PoiovOuVeA9733F9TkLIwHBgBbAV9j\nwxNXq05Cfqqqsgp4VVU2ZSjoBy3ZAkyQ/pLMiUs7x4qUZZ461XoLgn70oO6ehMRFqVq3tqCgWbP4\nIY2yMgt6MqEQExhTOfAX+3DKxo3JA4nI30GFtFq2tO9j0KV1a7uurg4+6AfVIGnVKrVx/I4d1a3f\nlOVVT0IQ7/2BMX8vBPZv7GtK7pSW2pz0wYPhzjujU9BiJZt5ANFVAVMVu7TzwoXxwwubNlnXfbdu\nlqH+8cfxz+3Uye6vawGoTp3ig4ROney6Zcv4JMlMdpF26VJ4B82RI6MH/i++sH+TZJ+h2OtBlJfb\n9Nwdd0z+mA0booHE0qX2HVy7Nv4SuW/FCpg/3/4uKYke1Hv3rrtbv1Wr3H1mEcWXwl57Wc2ESy6B\n44+vfSa5alXm3zP2gJp4hrpxox2wvvjCkqESrVxp9RnqmhvfvXv8ehCRUs9TptgQQ2xOQjFL5cDf\n1IdTMqFFC8tF6dUr7JaIZEYTKEopmXDNNdarcOGFtbfVNde4McWUIiZOtK7rXr3sbC1WspyEO+6w\n5YqTBTCxrxm77PSAATaUsXmzXccmLRbjGgyJB/q6DvzJ9qmINF0KEgSIJjHed58tfBIrccGlNm2i\nB4pnn238e0d6FT7/3Ao0xUoMGnr2jM5xf+UVy22YOzczB/hI13tk8apjj03v8xSSVA78sf9Ob76Z\n//kWItJ4Wipa/k91NQwfbgfYDz6wrlPITkJesoS5xCmKVVWW4R2x1VbwwgtwzDHRQKBdOwtkPvoo\n+rg2beJzDxqSZKcpjCJSaPJ6qWhpWkpKbBnpBQvgyiuj92fjDPKII+LP2n/2M7t/7Fg7uG/ZYtex\nAQLY8MM++9jMiMhQx+rV8QECWJARqyFJdql0vReiYhxOEZHGUZAgcfr3h0svheuvT3/mQkN88EHw\n7foO5pGOr549Lbg46qjgxyXOWmjIAb+xY+6zZ9v88mbN7DoxcAlb0HCKAgcRqYuCBKnlwguhXz+b\n8VBVldv3TjyYN2sWvH3pUktc/PDD2jkTLVvCww+nfsBPp8ck9iC7++7xvSDJilPlUmz7pk+P37Z4\nMRx5ZHzgcMQR4bRTRPKTggSppXlz+Oc/baneW2/Nznvstlvw7cSz+T32iH9cz5qi4LFnxQsXWlAT\nKZO7fj2MGgW/+132k+xi25EYUEWGPMI8W49tX+JMkW7davcWZbP3SEQKj+okSKDBg+E3v7Ghh6OP\nrruATDpiiylFEhehdkGioKRJqD0sMb9mfVHnbEhi3Tr4+c/h7LNt6CRbdfbrGh6JDHmkUrAo0xLb\nl7hA1vbb56YdIlKY1JMgSV11lU2NHDs2fvGbTEj2eoln3RA8BJA4LLFpk10SX/eOOyzB8bXXMtr8\n/5PYjtJSK4Hbpk20UFOYlQoT27fnnvH7MrHORSbqXohI06EgQZJq08ZKNb/wAowfn9nXTlaToKG1\nCuoqwAQW3ETO5OfPhwMOsJ6RdYEriKQvth2DBkFFhZ2d77qrldKFcGdNBCVjxgZiYBU3M1n3QkSa\nEO99qBdgIOBnzJjhJT+ddJL3HTt6v3Sp90uWeD90qPe9etn10qXpvWavXt7beb9devWq+/66DB0a\n/xzwftAg7z/91Pthw+Lv33FH7197LfX2NuRzDxpUuw3e22Mzsc8yJXF/lZZ6X1bmfZs23s+e3fDX\nydR3QUTSN2PGDA94YKDPwjFaPQlSr1tusbH+88/PXFXCZGfX6Zx1T5xYe42HLVugTx8bZrjllvhe\nhf33h1//2taAaKiGfO5k0zrzrVJh4nBHVVXDZmQkDgUlzowohAqVmvIpkhoFCVKvrbeGm2+2aYWJ\nFQjTHV9P1g2+aZMNH5SXWzd4Q6cuxi7/DNGVI0tKbJbDBx/AsGHR7XfcATvtZHkX331X/3skyyuI\nPejErmSZz+oKvBKLUMVKDJTefTd+eyGsClmMpbdFGkNBgjTIKafAIYfUXlAp3fH1oLPrkSPtwLNx\no12aNcvcWXdsr0Kk1+Hbb+Gyy2xa5cUXw/LlyZ+frIcj9qCTmDSZOM0zX8QGaKWl8dvqWjo7MQhI\n/LyFUKGy2Je7FkmVggRpEOesdkJ5uSUFZmMlwMb8gDckSz/Sq/DJJ3DGGdED5Jo1trhVz55w3nlW\ndyFRsmqMQVMMI4955pmGtz+XYgO0WbMsQTVxRkaQuoKA8vLk34Wwu/hj3z8xECyEwEYkVNlIdEjl\nghIXC8pDD1my28MPZ/61ExPqhg5t+HMbmhwYm2xXUeH9qad637x5/Ps2a+b9CSd4/+yz3m/alJ02\nF2LSX+w+btOm4Z+7Mf+umZD4/m3aFNZ+F6lLthMXFSRIykaN8r5dO++/+iqzr5vuLICGHHAjjykv\nr33AWrjQ+/PO875ly/ht4P3WW3v/2996/+673ldXZ67NYR84GyuVz53OjJVMCvv9RbIp20GCloqW\nlH37rXXn77QTvPSSdePnWuxS08uX178s9L77Rqsexiovh6++si74ZcssZ+Gee2DFitqP3XlnOPVU\nOPnkaHnodBXTstSJ+74hy3Y3pfcXySYtFS15p0MHuP9+ePVVO6iGITZhMDZAgOBchmT5DRs3RjPc\nu3SBa66BRYvg6afh+OPjCzXNnQuXXAI77GD5DKWllgQ5e3bq7W9Ky1LXl3PQ2NU1Gyvs9xcpZFq7\nQdJy4IFwwQVw0UVw8MG5L+dbV1Jj0AG3W7f4M/e6XqtZM1sN8YgjbHrkE0/Agw/C5MnRx1RX2/UP\nP9hnHzvW6gYccEBwBchEEycGr0lRiOpbmyJxPY5cC/v9RQqZehIkbVdfbV3wJ54I33+f2/fu3Dn+\ndqtW0foKmzfXfTbbpk38trrO4tu3tyWzX3/dSi4H8d7qLhx6KHTqZAfNBx+sPV00Vr4VWGoMTSsU\naboUJEjaWrSARx+FBQvg3HNz+96JqTQlJdH6CtOm1S6SE3tQ/uyz1LufZ8+GmQ0Y7fv+e3u90aPt\nPQ84AP7+d/jyywZ/tIKT70MnYU/BFClkGm6QRtl5ZzuLPu00K+l71lm5ed/EksobNsTfrutsNp3u\n58GDg1eubNMGXn4Zvv7a6iI8+2w06bGqygo4vfaa1V/YbTc46ii77LGH1Z4oZJHk0YULbT907Ajd\nu+ff0EmYS3WLFDr1JEijjR4NZ59tvQlTp+bmPRPPVhMrBWb6bDZoOGX2bCvENGiQHXj+/W9YssRy\nF37/++gqixGzZsFf/mLDFi1awC9/aQFGYknpQhE5+C5YYMmj3bvn59CJhkNE0qcgQTLillvsYDly\npB0os+3OO+MrBT77bMOHEDLV/TxmTO3XWrnS1oi48Ub49FMLJK6+2tahiLVpE9x7ryV9du1qeQ/P\nP1846z9Aww++YXf35/twiEhey0bxhVQuqJhSQYstZDRokPdduni/777eb9yY3fdtTDGidJ6bWGQJ\nvO/ePbXX2n774NeJvbRv7/3o0d4//bT369c3/DOFoaGfPezCUR9+aFUW01kOWyTfaaloyWux9Qre\nfRe22Qbeece627Mp1S7kyNlsz57w9tupPReC8wcWL06tHV27xt/u3dtqMbRqFb3vu+/ggQdsOuXW\nW8NJJ8GECTbVMt80tP5A2N39Y8facEhkOewxY3L7/iKFTEGCNEriD/7338Ott8Jtt9kYfSY1ZqGe\nSDDz1VfprV44YEDw/al0ZSe+b8eO8Nhj9lkmTrRKjrHTM9etg8pKOO44CxiOOw4eeSS+eFSYXfkN\nncYZdnd/2EGKSCFTkCCNEnQAGDPGxtjHjLHM/kxJrLLYpk3DpzHWdWC45pr63zvoTL5Vq9q5EXfd\nlfw1EmdkRG63agXHHAPjx1vA8OyzcPrpVtky9v0nTIBRoyxgOOoo63E46qjoPnnrrdpTP/NB2BUP\nww5SRAqZpkBKowRVDnQOxo2D+fPt4Ddlik2VbKzEA/3WWzd8vYO6Ki7+7Gc2SyGV55eW2ucaMyZ6\nZh/pyk42vS7xNYIOVuXlcPjhdtm82UpfT5gATz4Z7T3ZuNHKRj/9dO3n59NZcuz6Gt262f4KY+ZD\nU6puKZJr6kmQRknW5dysmZUz3m47OOywzHSDN+aMMHI2G7Qw0/r1DX9+5Gx40SIbgkilKzvVM+pm\nzeCQQ6x3YtEiCxjOOafuz71ypdWtyMUMk/rE9vyE2cuRz9Utw575IVIfrQIpWbVggRUi6tHDagIk\nlkROxbJltc8I0/nBb9s2fly/TZv6exKSCWOFwepqOys/7LDk7XbO2nbccbbPunfPXnsSewwi/y7F\ntNJlurRCpTSWVoGUgtazJ/znPzBvHhx9dO3KiKnI1BnhlCnxeQRTpqTfpjDG20tK7L1+/OP4+2MX\nlvIe3ngDfvc7C9CGDLHaDfPnZ749yXoMlAtQPyVVSr5TkCBZN3CgJeO9/bYtBrV5c/bfs65u3AED\n7Ax882a7TjZzoSFiA5cJE+wAmauu48QAZcECeO89W846MQdk6lT4n/+xx+65J/z1r/DJJ5lpR7ID\nXSpJncVKgZTkvWwUX0jlgoopFY3nnrOCNqec4n1VVXbfK4wCPmEXDUr00Ufe//nP3u+6a/LiTT/+\nsT3mww+9r65O732Sfe582x/5aOnSaDGyoUPttkgqVExJmoxDD7Vpfg89BL/9bfCCSZmS7Ow2m4li\n+dZ1/KMfwZ/+ZGtGzJtnUz0Tl7v+8EO4/HIbuujd24YnXn45td6eZEMu+bY/8lE+J1WKgIYbJMdO\nOMG6nW+/3Q5g2ZKsGzebGfdt29Z9O0x9+8JFF8H06fbZb7zREkpjzZ9vhbAOPhg6d7ahoYcfhm+/\nrfu1kx3o1JUuUvgUJEja0j0r/9Wv4IYb4Kqr7GCVDcnObhcujH9c4u3G+Oyzum/nix13tLLZU6bY\nEteRwKAspmrKmjXw6KNWBXLrreGAA+Cmm1L7TGEXURKRxtMUSElbY6dvXXqprZB4001w/vmZb1+Q\nTE5/TNSsWfyyz2VluUnSzJTVq2HSJCvS9NxzyXsQdt7Z1pY48kjrjSgtzW07RSRKUyAlbyWOMb/9\nNnz0UcOff+WV8Mc/wgUXZK9HIVHHjnXfboyWLeu+ne/atYOf/9zyRpYts5LaF1wAffrEP27uXLj+\negsSu3a1EtITJsQHX2DLZLdta8FT27apfTdEJD80Kkhwzl3knKt2zt0Uc1+5c+5259wK59xa59wT\nzjml4zRBiWPM3tt8/IZyzpLpLr3Upudde21m2xcksahQJosMPfNM9Ky6tNSmfRaqsjLYbz/4299s\nquScOXDddRYYlMT8aqxYAfffb0WbOnWC/fe3f9Pp02GffeJXX0zlu9FQqlgokmXpTosABgFfAO8B\nN8XcfwfwJbAfsAfwNvBGHa+jKZAFaunS2lPqyspSf53qau8vv9yef9VVGW9mnGxOOUtlyt+SJYU7\n9W35cu/vv9/7kSO9b906+fTKTHw36qNpllLssj0FMq0FnpxzrYHxwJnAZTH3twV+AZzovX+95r4z\ngDnOub2899PSeT/JT1262Jh+bDdzOl3szsEVV9gZ6qWXQlVV9mY+RDLxsyGVKX+RWRZgsw2OPbZw\nyvF27gyjR9tl40YblnjmGctnqCuxsbrack+GD7ez/q23bnxbNM1SJLvSXQXyduAZ7/0rzrnLYu7f\ns+Y1X47c4b2f55z7ChgCKEhoYqZMsW7k9estQGhMieM//cm6uS+5xAKFK66wAKJQNGSVx4hszrLI\npfJyGDHCLmCf/4UX7DJpUvwS29XVcMstdgFLgBw2LHrp2TP1f+9U9rlIPti82eqTTJ1qQ3G//W3Y\nLapbykGCc+5EbBihImDzNsAm731ivvhSoGvqzZN8FylxnCkXX2zj+X/8owUe111XOIFCKksSr1pV\n9+1C1asXjB1rly1bYNo0CxZeeMH+rq6OPnbuXLvcc4/d7t7dgoV997XrAQPi8x+CaBloyXcLF1pA\n8M47dj19enQNm27d4De/ye/fuJSmQDrnugPTgZ947z+sue9V4D3v/QXOuVHAvd77lgnPmwa85L2/\nOOA1BwIzhg8fTrt27eK2jRo1ilGjRqX6maQJuPVWq/531lkwblzTm2a3ww621kJEz57w5ZdhtSY3\nvv3WFp2KXGbMiJ8ymqhDB5tWG+lpqKiA5s1z195cSLaCphSmH36w73VsUPDNN3U/56uvbBG2hqis\nrKSysjLuvtWrVzN58mTI0hTIVIOEo4CJQBUQiX1KsaSJKuCnwEtA+9jeBOfcl8DN3vu/B7ym6iRI\noPvug1/+0qblPfCATaVrKrREMHz/vf2QvvGGffYpU+y+ZFq0gL33jgYNQ4Y0bunxfKDvQeHy3mb+\nRIKBqVPhgw9sqLQuvXpZfZHBg+37vMcejftty3adhFSHG14CEhao5T5gDnAt8A2wGTgIeBLAOdcX\n2B5oxGi1FKPTT7f59SeeaMmRjz9eeLUHklE3OWy1FRx4oF3Axmrffz/a0/DmmzbFMmLDBnj9dbuA\n9S4NGGCP2bLFzsb+8x/YZpvcf5Z0KfGycKxaZUNmkaDgnXfqL1nepg3stVd8UJCJhN1cSilI8N5/\nD3wce59z7ntgpfd+Ts3tfwE3Oee+BdYCtwJvaWaDpOPYY63ewNFHw2GHWTXAQj97hOzOsihUzZrB\noEF2ueACO1ObOzc+aIgdkqmqsjO3iGXLYLvtrL7DwIE2PDFwoBWDqi+3ISxKvMy8TAzhbNkSTS6M\nBAXz5tX9HOdgl10sEIgEBTvvXPhDpY0uy+ycewV433t/Qc3tcuBGYBRQDjwPnOO9DyxzouEGaYi3\n3rIgoV8/+O9/rXCPFJ+FC+PzGmbPrv85bdpYl25s4NCvX378eC9bVrtHSTkJjZPOEM6iRdEhg6lT\nLa8gdmZOkC5dor0DgwfDnnuGs6hbtocbtHaDFIz33rOpdp06wfPPW7KfFLfBg+1ML6JZs4atl9Gy\npa2M2a9f/HXfvtC+ffbaK9nXu3d870yvXrZCacT69RYExOYS1DcFuXlzCzRjhw122CE/ZiXkW06C\nSGj22MPWhxgxwpLWnnsOdt897FZJmJ5+uvaZeFWVBZQzZsDMmXb99dfxz1u/HmbNskuiLl2CA4je\nvZve7IqmKHEIp317ePDBaFAwa1bds2rAAoBIQDB4sP3OlJdntdl5Sz0JUnCWLoXDD7fM4iefhIMO\nCrtFku+WL48PHGbNsgNJfZnosUpKbJnt2F6HyN/bbZcfZ5WFassWS0zduDH4uq5tiderVtkJxJo1\n1qtUX0DQurXlwcT2EhRS8qt6EkQSbLONlQI+7jg49FCbKnnSSWG3SvLZ1lvDIYfYJWLzZpg/3xLS\nPvkkev3JJ8GzDKqrrdv688/tIBSrVavawxZ9+liBqG22yd/pu1VVqR+EM/WY2MfGFtnKJufgRz+K\nTy780Y/yIz8lXylIkILUurWtF/CrX8HJJ9uP+u9/n/t2qBhOcrNn20qQsSW7BwwIu1VRzZpFD+iJ\n1qyBTz+NDxwif69bV/vxP/xg0zfff7/2NucsUNh2W+tx2G676N+R69at0z/ANuYxqfSkFKLOneOT\nCwcNsiXRpeEUJEjBatYM/v1v+5H9f//PpsfdfLOt/5ArhbxQU7ZFloqG6FLRmSzhnU1t29pMiIqE\n4vPeW0CYGDjMm5d8+MJ7WLLELjMz3hlcOJo3t4JY5eXB13VtS+cx7dpZ7QwNAzWOggQpaM7B1VfD\n9tvDuefCnDnw2GPQsWNu3j+VYjjF1uuwfn3dtwuRc3b2v+22sP/+8ds2b7ZAIRJAzJ9vJXkXLbLr\nJUvCOXNv1iz1A22mD9zNm+dvrYowFNJvgYIEaRLGjLHx4OOOswpnTz9tY43ZlkoxnGLrdWjZsvHL\niBeSZs3sO9ivX/D2qiqri/DNN/HBwzffWABV3wE3nQN2ebkOzvmokH4LFCRIk7H//lY29aijbPzx\n4YfhZz/L7numUl652ErwZnIZ8UKS7CyxtNRud+tmhXekeBXSb4FiTGlSevWyWgoHHABHHmlLT2/a\nlL33i5RX/vxzu66ryzCxl6Gpl+CNLCO+ebNd51PSYjZFzhK/+MKujz027BZJvimk3wIFCdLktGlj\n9ROuvhpuuMES6Oqru54LEydaidhevey6mBZ1WrrUyuX27m3XywKLtDeNdhTSWaKEo5B+CxQkSJNU\nUgIXXWRd3GvXWr3+u++2TPOwpNLr0NTky9l1LtpRSGeJEo5C+i1QkCBN2p572rSzU06x5MZjjrHq\ne5Jb+XJ2nYt2FNJZokh9FCRIk7fVVnDXXfDUUxa177qr/R1yRfKikrhqZ1ireObiLL+QzhJF6qMg\nQYrGUUfZGvEVFdajcNhhxV3cJlWNGc9PLGgTVoEbneWLpEZBghSVbt2snPOTT9q4dCRgCFoNUOI1\nZjx/xYq6b+eKzvJFUqMgQYqOc3D00fDRR3D//da7sPvucPzxtt6ABGvMeL6S+UQKk4IEKVplZTB6\ntJVyvvdemD7d8hVOOMHOMpWzEK8xB3p184sUJgUJUvSaNYMzzrBaCnfdZcHCsGFW1vlvf9NsiIjG\nHOjVzS9SmBQkiNRo3tyWnv70U3jpJRuCuPhiW2Xy5z+HF17I3br3+UgH+szLlyJTIskoSBBJUFIC\nBx0ElZW2CM/111v+wogRdhb9l7/EL+okkq58KTIlkoyCBJE6dOoE551nCY1TpljwcMMNdua3zz5w\n++0ajpD05UuRKZFkFCSINIBztrLkv/5lXcSVldEAYtttbbXJykr44YewWyqFRLM+JN8pSBBJUatW\ncOKJVm9h8WL4+99h1So46SQbpx89GiZNgi1bwm6p5DvN+pB8pyBBpBE6d4Zf/9qWp/78c/jjH2Ha\nNPjpT6F7d+tpePddTaeUYEoGlXynIEEkQ3r1gksvtboL06dbz8Kjj8Jee8HOO1vC4+efh91KEZGG\nU5AgkmHOWbnnm26ChQtt6uSQIZbw2KeP/X3bbUp4bAxNHRTJDQUJIllUWgo/+Qncd58d2B55xIYo\nzj/fktSOOAImTICNG8NuaWHR1EGR3FCQIJIjrVpZyedIwuOtt1pvwnHH2QyJc8+1YQrlL9RPUwdF\nckNBgkgIIgmPU6fCxx9bpceJE2HQIPjxj+HGG3Xgq4umDorkhoIEkZD17w/XXgtffQX//S/ssosl\nQPboAYcfDo8/Dhs2hN3K/KKpgyK5URZ2A0TElJXZ1Mmf/hS+/dZmRtx/v60b0aEDjBoFp58Oe+5p\nyZHFLDJ1UESySz0JInmoQwcYO9ZKQc+ZA2PGwFNP2XTKXXaxmRIajhCRbFOQIJLndt4Z/vpXG454\n/nnYdVf405+sWNNhh8Fjj2k4QkSyQ0GCSIEoLbWVKCsrrRfhjjvgu+9sxkS3bpYIOW2aZkeISOYo\nSBApQO3bw1lnWTnouXPh7LPh6adh773hRz+C666zZa5FRBpDQYJIgevXD665BhYssIWlBg6EK66w\n2RGHHmoJkBqOEJF0KEgQaSJKS+GQQ+Chh2DJErjzTlizxlas7NbNehumTtVwhIg0nIIEkSaoXTsr\n0PTWWzBvHpxzDjz7rK0bEanL8M03YbdSRPKdggSRJq5vX7jqKvjyS3jxRauz8Oc/w/bbW02GykpY\nvz7sVopIPlKQIFIkSkvh4INh/HgbjrjrLli3zpa03mYbK9T0wguwZUvYLRWRfJFSkOCcG+ucm+Wc\nW11zeds599OY7a8556pjLlXOuXGZb7aINEa7dnDmmVa18JNP4Pe/t5kSI0ZY/YXf/U7TKUUk9Z6E\nr4ELgYqayyvA/zrn+tds98DdwDZAV6Ab8IfMNFVEsmGnneDyyy134d13rWfh8cdtOmXfvtFtIlJ8\nUgoSvPf/8d4/773/rOZyKbAOGBzzsB+898u998tqLusy2mIRyQrnLF/hppvg66/hpZdg2DC45Rar\n+hjZpvoLIsUj7ZwE51yJc+5EoBXwdsymk51zy51zHzrnrnHOtWx0K0Ukp0pL4aCD4N57YelSeOIJ\n6NkTLrrIhiMi2777LuyWikg2pRwkOOd2cc6tBTYC44BjvPeRzsiHgFOA/YFrgFOBBzPTVBEJQ4sW\nMHIkTJhgAcM//2n3n3kmdO0Kxx5rtRkUMIg0PeksFT0X2A1oD4wEHnDODffez/Xe/zPmcR8555YA\nLznndvTez89Ae0UkRO3bwy9+YZdFi6yaY2UlnHIKNGsGBx4IxxwDRx1lAYSIFDbnG5m+7Jx7EfjM\ne392wLZWWM7CCO/9i0mePxCYMXz4cNq1axe3bdSoUYwaNapR7ROR7Fu40JayfvJJeP11qK62wk3H\nHmtBQ69eYbdQpPBVVlZSWVkZd9/q1auZPHkyQIX3fmam3zMTQcLLwALv/S8Ctg0FJgO7ee9nJ3n+\nQGDGjBkzGDhwYKPaIiLhW7kSnnnGAoYXXrB1I3bd1XoXDj/cEiBLS8NupUjTMHPmTCoqKiBLQUKq\ndRKuds7t65zrWZOb8FdgP2C8c66Xc+5S59zAmu1HAvcDrycLEESk6enUyQoz/e//wvLlNp1yl13g\ntoSfc9kAAA2hSURBVNtg8GDo3NlyHO64Az79VLUYRPJZqjkJ2wAPYPUPVgMfAId4719xznUHDgZ+\nB2yF1VR4HLg6c80VkULSujUcd5xdtmyBd96xqZUvvQS//a3dt/32VgnyJz+xnIYuXcJutYhEpBQk\neO/PrGPbQmxWg4hILWVlMHSoXS6/HNauhcmTo0HDvffa43bbzYKGgw+2Og1bbRVuu0WKWTqzG0RE\nGq1NG8tROPxwu714Mbz8sgUMjzwCf/sbNG8O++wTDRoqKizYEJHc0AJPIpIXunWzqZT33WcVH+fM\nsUChXTu4/vpoPsMxx8Dtt1upaOUziGSXYnIRyTvOWSnonXeGc8+13IV3340OTZx/PmzeDD16WA/D\nfvtZj0OfPvZcEckMBQkikvfKyqzuwpAhcNlltsT1G29YwPDii/Dvf9vjOne2HofIYwcNsuRJEUmP\nggQRKTitW8Ohh9oFrCT0O+/YctdTpsB118GaNVaP4cc/tl6GSODQq5d6G0QaSkGCiBS89u1hxAi7\nAFRVWU7DlCkWOLz8MowbZ9u23tp6GwYOhN13t0vPngocRIIoSBCRJqe01Ao47bIL/OpXdt/Kldbb\nMGWKXd9+O6xYYdvat48GDHvsYdf9+9t6FCLFTEGCiBSFTp3gsMPsAjYzYtEieP99u7z3npWTvuUW\n2968OQwYEA0adt/daji0bRveZxDJNQUJIlKUnIPttrNLpFYDWC7DrFnxwcODD9psCoDeveN7HHbf\nHbbdVsMV0jQpSBARidG2rVV6HDYset+mTTB3rgUMkeDhxhstYRIszyESMEQCiL59tZCVFD4FCSIi\n9Wje3Fay3HVXOO00u897WLAgvsfh0Ufhhhtse3m5BQo772z5DZHrvn2hVavwPotIKhQkiIikwTnY\nYQe7HH109P5Vq2y44oMPrPdh7lxbo2Lp0ujzevaMDx623daSJJs1s4Ak8nfQJXF7aamGOiR7FCSI\niGRQx45wwAF2ifXtt9GgYe5cm6L57LPw979DdXXj3rOuIKKh2xq7vbGvHXaws2ULrF9vlw0bon9n\n4nayx/znP1YxNJ8pSBARyYEOHaIFnWJt3Gi9D5s3J79s2pT+9vq2ff99+u+bSc5lNgApKUntgL5l\nS2rtbdHCLi1bxl9i7+vcue7H9O2b2X2YDQoSRERCVF5ui1sVGu/twJrroCbZtnXr4rdXV8cflNu2\nDT6QJ7uvrseUl1sQUgwUJIiISMpiz/yl6SqSWEhERERSpSBBREREAilIEBERkUAKEkRERCSQggQR\nEREJpCBBREREAilIEBERkUAKEkRERCSQggQREREJpCBBREREAilIEBERkUAKEkRERCSQggQREREJ\npCBBREREAilIEBERkUAKEkRERCSQggQREREJpCBBREREAilIEBERkUAKEkRERCSQggQREREJpCBB\nREREAilIEBERkUAKEkRERCSQggQREREJlFKQ4Jwb65yb5ZxbXXN52zn305jt5c65251zK5xza51z\nTzjnumS+2QJQWVkZdhMKjvZZerTfUqd9lh7tt/ySak/C18CFQEXN5RXgf51z/Wu23wIcDowEhgPb\nAhMy01RJpP9MqdM+S4/2W+q0z9Kj/ZZfylJ5sPf+Pwl3XeqcOxsY7Jz7BvgFcKL3/nUA59wZwBzn\n3F7e+2kZabGIiIjkRNo5Cc65EufciUArYArWs1AGvBx5jPd+HvAVMKSR7RQREZEcS6knAcA5twsW\nFLQA1gLHeO/nOuf2ADZ579ckPGUp0LXRLRUREZGcSjlIAOYCuwHtsdyDB5xzw+t4vAN8HdtbAMyZ\nMyeNphS31atXM3PmzLCbUVC0z9Kj/ZY67bP0aL+lJubY2SIbr++8r+v43YAXcO5F4DPgMeAloENs\nb4Jz7kvgZu/935M8/yTgoUY1QkREpLid7L1/ONMvmk5PQqISoByYAWwBDgKeBHDO9QW2x4YnkpkE\nnAx8CWzIQHtERESKRQtgB+xYmnEp9SQ4564G/otNhWyDHdz/BzjEe/+Kc24ccChwBpavcCtQ7b0f\nlumGi4iISHal2pOwDfAA0A1YDXxATYBQs/18oAp4AutdeB44JzNNFRERkVxqdE6CiIiINE1au0FE\nREQCKUgQERGRQKEHCc65c5xz851z651zU51zg8JuU75wzl3unKtOuHwcs10LagHOuWHOuaedc9/U\n7KMjAx7zF+fcIufcD865F51zfRK2d3DOPVSzcNm3zrl/Oue2yt2nyK369plz7t8B373nEh5TbPvs\nIufcNOfcGufcUufckzUzuGIfU+//SedcD+fcf5xz3zvnljjnrnfOhf5bnC0N3G+vJXzXqmoS4WMf\nUzT7LROLKWZqf4W6g51zJwB/Ay4H9gBmAZOcc53DbFeemY0ljHatuewbs00LapmtgPexJNlaSTbO\nuQuBc4ExwF7A99j3rHnMwx4G+mNTeA/H9udd2W12qOrcZzX+S/x3b1TC9mLbZ8OAfwB7AwcDzYAX\nnHMtYx5T5//Jmh/p57Ck8cHAacDpwF+y3/zQNGS/eeBuot+3bsAfIhuLcL81ajHFjO4v731oF2Aq\n8PeY2w5YCPwhzHblywULnmYm2dYW2IiVxY7c1w+oBvYKu+0h7rNq4MiE+xYB5yfsu/XAz2tu9695\n3h4xjxmB1f3oGvZnCmmf/RuYWMdzdi7mfVbzeTvX7IN9Y75Xdf6fxKaIbwY6xzxmDPAtUBb2Zwpj\nv9Xc9ypwUx3P0X6DlVh5gZx+z0LrSXDONcMipNgFoTxWtVELQkXtVNMl/LlzbrxzrkfN/VpQqwGc\ncztiZyax+2kN8A7R/TQY+NZ7/17MU1/Czm72zlFT89H+Nd3Dc51z45xzHWO2DUH7rD32eVfV3G7I\n/8nBwIfe+xUxrzMJaAcMyHaD80Tifos42Tm33Dn3oXPumoSehqLdby69xRQztr/CHG7oDJRiC0DF\n0oJQUVOxLqIRwFhgR2ByzbhvV7SgVkN0xX6Q6vqedQWWxW703ldhP2LFui//C4wGDsS6ffcDnnPO\nuZrtRb3PavbDLcCb3vtInlBD/k92Jfi7CMW738BK858C7A9cA5wKPBizvej2m3NuF+fcWqzXYBw1\niymS4+9ZJsoyZ1p9C0IVDe99bJnN2c65acAC4OckL2Gt/dcwDdlPRbsvvfePxdz8yDn3IfA59iP+\nah1PLZZ9Ng74EfE5Qsk0dJ8U034bGnun9/6fMTc/cs4tAV52zu3ovZ9fz2s21f2W6cUUI1LaX2H2\nJKzAqjNuk3B/F2pHQAJ471cDnwB9gCVAc+dc24SHaf/FW4L956nre7ak5vb/cc6VAh3QvgSg5od6\nBfbdgyLeZ86524DDgP2994tiNjXk/+QSan8XI7eLab8trufh79Rcx37fimq/ee+3eO+/8N7P9N5f\ngiX2/44cf89CCxK895uxRaEOitxX0xV1EPB2WO3KZ8651kBvLBEvdkGtyPaGLKhVVGoObkuI309t\nsXHzyPdsCtDeObdHzFMPwoKLdxCcc92BTkDkx70o91nNge4o4ADv/VcJm+v6Pxn7XftxwgyuQ7Ay\n97Hd701KPfstyB7YGW/s963o9luCoMUUgSx/z0LO1vw5lmU+GsuWvgvL4Nw67EzSfLgAN2DTW3oC\n+wAvYlFgp5rt44D5WBdwBfAW8EbY7Q5hP22FdcvtjmX4nldzu0fN9j/UfK+OAH4MPAV8CjSPeY3n\ngOnAIKwrdB7wYNifLYx9VrPteiyQ6lnzYzQdmAM0K+J9Ng7LDh+GnZVFLi0SHpP0/yT2Qz8Ly/nY\nFcs3WgpcGfbnC2u/Ab2AS4GBNd+3I4HPgFeKdb8BV2NDWT2BXYC/YoHBgbn+nuXDzvg1tkz0eiz6\n2TPsNuXLBajEpoSuxzJXHwZ2jNlejs0/XoGtuvk40CXsdoewn/arOdBVJVzujXnMFVgPzA9Ylm+f\nhNdoD4zHIu1vgXuAVmF/tjD2Gbb07PNYD8wG4AvgDhKC9yLcZ0H7qwoYHfOYev9PYoHYs8C6mh/u\n64CSsD9fWPsN6A68Biyv+f85r+ag2LpY9xvwz5r/d+tr/h++QE2AkOvvmRZ4EhERkUBNsqSliIiI\nNJ6CBBEREQmkIEFEREQCKUgQERGRQAoSREREJJCCBBEREQmkIEFEREQCKUgQERGRQAoSREREJJCC\nBBHBOfeqc+6msNshIvlFQYKIiIgEUpAgIiIigRQkiEhEiXPuOufcSufcYufc5ZENzrlq59xY59xz\nzrkfnHOfO+dGhtlYEck+BQkiEnEatqzsXsAfgD855w6K2f4XbEnaXYGHgEecc/1y3koRyRktFS0i\nOOdexdaa3y/mvneAl733FzvnqoFx3vtzY7ZPAWbE3iciTYt6EkQk4oOE24uBLjG3pyZsnwL0z2qL\nRCRUChJEJGJzwm1P/b8R6ooUacIUJIhIQw0OuD03jIaISG6Uhd0AESkYxzvnZgBvAqcAg4Azwm2S\niGSTggQRgeTDBrH3Xw6cCNyO5Suc6L2fl+2GiUh4NLtBROpVM7vhaO/902G3RURyRzkJIiIiEkhB\ngog0hLocRYqQhhtEREQkkHoSREREJJCCBBEREQmkIEFEREQCKUgQERGRQAoSREREJJCCBBEREQmk\nIEFEREQCKUgQkf+/UTAKRsEowAoAZ3h0d3UcEiEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0xc0e07b8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "res_bs.plot_partial(1, cpr=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Poisson\n",
    "\n",
    "GLMGam is a subclass of GLM and supports the same families. However, currently only Gaussian and Poisson have been tested and verified.\n",
    "\n",
    "Miles per gallon can take on only positive values, so we can use Poisson with a log-link to impose that the prediction is non-negative. Except for adding the family, everything works the same as in the Gaussian case. One difference is that the non-linear link implies that predict_partial and plot_partial are not in terms of the response variable. Those are still for the linear prediction, i.e. the values before applying the inverse link, ``exp`` in this case. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                 Generalized Linear Model Regression Results                  \n",
      "==============================================================================\n",
      "Dep. Variable:               city_mpg   No. Observations:                  203\n",
      "Model:                         GLMGam   Df Residuals:                   194.75\n",
      "Model Family:                 Poisson   Df Model:                         7.25\n",
      "Link Function:                    log   Scale:                          1.0000\n",
      "Method:                         PIRLS   Log-Likelihood:                -530.38\n",
      "Date:                Tue, 29 Jan 2019   Deviance:                       37.569\n",
      "Time:                        19:46:43   Pearson chi2:                     37.4\n",
      "No. Iterations:                     6   Covariance Type:             nonrobust\n",
      "================================================================================\n",
      "                   coef    std err          z      P>|z|      [0.025      0.975]\n",
      "--------------------------------------------------------------------------------\n",
      "Intercept        3.9960      0.130     30.844      0.000       3.742       4.250\n",
      "fuel[T.gas]     -0.2398      0.057     -4.222      0.000      -0.351      -0.128\n",
      "drive[T.fwd]     0.0386      0.075      0.513      0.608      -0.109       0.186\n",
      "drive[T.rwd]     0.0309      0.078      0.395      0.693      -0.122       0.184\n",
      "weight_s0       -0.0811      0.030     -2.689      0.007      -0.140      -0.022\n",
      "weight_s1       -0.1938      0.063     -3.067      0.002      -0.318      -0.070\n",
      "weight_s2       -0.3160      0.082     -3.864      0.000      -0.476      -0.156\n",
      "weight_s3       -0.3735      0.090     -4.160      0.000      -0.549      -0.198\n",
      "weight_s4       -0.4187      0.096     -4.360      0.000      -0.607      -0.230\n",
      "weight_s5       -0.4645      0.103     -4.495      0.000      -0.667      -0.262\n",
      "weight_s6       -0.5092      0.112     -4.555      0.000      -0.728      -0.290\n",
      "weight_s7       -0.5469      0.119     -4.598      0.000      -0.780      -0.314\n",
      "weight_s8       -0.6211      0.137     -4.528      0.000      -0.890      -0.352\n",
      "weight_s9       -0.6866      0.153     -4.486      0.000      -0.987      -0.387\n",
      "weight_s10      -0.7370      0.174     -4.228      0.000      -1.079      -0.395\n",
      "hp_s0           -0.0247      0.010     -2.378      0.017      -0.045      -0.004\n",
      "hp_s1           -0.0557      0.022     -2.479      0.013      -0.100      -0.012\n",
      "hp_s2           -0.1046      0.038     -2.719      0.007      -0.180      -0.029\n",
      "hp_s3           -0.1438      0.050     -2.857      0.004      -0.242      -0.045\n",
      "hp_s4           -0.1919      0.063     -3.047      0.002      -0.315      -0.068\n",
      "hp_s5           -0.2567      0.079     -3.231      0.001      -0.412      -0.101\n",
      "hp_s6           -0.4152      0.120     -3.455      0.001      -0.651      -0.180\n",
      "hp_s7           -0.4889      0.152     -3.214      0.001      -0.787      -0.191\n",
      "hp_s8           -0.5470      0.195     -2.810      0.005      -0.928      -0.166\n",
      "================================================================================\n"
     ]
    }
   ],
   "source": [
    "alpha = np.array([8283989284.5829611, 14628207.58927821])\n",
    "gam_bs = GLMGam.from_formula('city_mpg ~ fuel + drive', data=df_autos,\n",
    "                             smoother=bs, alpha=alpha,\n",
    "                             family=sm.families.Poisson())\n",
    "res_bs = gam_bs.fit()\n",
    "print(res_bs.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Selection of penalty weight alpha\n",
    "\n",
    "The model class has two methods to select the penalization weight alpha that optimizes some goodness of fit criteria.\n",
    "\n",
    "`select_penweight` uses a one step approximation to the leave-one-out change in parameters. This only requires one fit for each alpha. This uses scipy optimizer, by default `basinhopping` to minimize either an information criterion 'aic', 'bic' or generalized cross-validation criterion 'gcv'. The default is 'aic'.\n",
    "\n",
    "`select_penweight_kfold` uses k-fold cross validation on a fixed grid to find the alpha with the lowest cost function. The cost function can be specified by the user, the default is mean squared prediction error. There are currently no premade cost functions for non-Gaussian models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([  8.28404434e+09,   1.46278733e+07]), 5.243746995925903)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import time\n",
    "t0 = time.time()\n",
    "gam_bs.select_penweight()[0], time.time() - t0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((1995262.3149688789, 125.89254117941663), 3.335373878479004)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t1 = time.time()\n",
    "gam_bs.select_penweight_kfold()[0], time.time() - t1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}

Hello, thanks for this valuable resource on how to use GAM in Python. I would like to know how I should specify knot locations.
I tried bs.knot_kwds={'knots':(num1,num2)}
It runs without error. But, is this the correct way?

The docstring says list of dict, because splines could have several variables and we need several splines. Try brackets [ ] around the dict.

You can inspect the bsplines instance to see how the knots where set.
Note: in Bsplines everything is in list of splines, e.g. from above bs = BSplines(x_spline, df=[12, 10], degree=[3, 3])

AFAICS, [s.knots for s in bs.smoothers] should show the knots for each univariate bspline.