SNN_python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# **Semiparametric Neural Networks and the \"OLS Trick\" -- in Keras, via Python**\n",
    "\n",
    "This notebook has a couple of purposes.  The first: it is a chance for me to teach myself how to use jupyter notebooks.  (They seem pretty simple so far).  I am new to python.  Please forgive any amateurish code in here (or better yet, send me some suggestions).\n",
    "\n",
    "The main purpose is to give python examples of a semiparametric neural network, and the \"OLS trick.\"  I've gotten a good amount of interest in the methods since developing them for [this paper](https://iopscience.iop.org/article/10.1088/1748-9326/aae159).  I implemented that paper using my hand-coded R package `panelNNET` ([here](https://github.com/cranedroesch/panelNNET)).  Developing that package was a fantastic way to really learn how neural nets work.  But I've since learned that all of my ideas can be implemented in Keras, with a bit of hacking.  Keras is somewhat faster, but its main advantages are that it is more general and that it is maintained by an active community.  \n",
    "\n",
    "### The basic idea behind semiparametric neural nets\n",
    "\n",
    "A feed-forward neural network with a single continuous output and one hidden layer is fully described by the following equation:\n",
    "\n",
    "$$y = a(\\mathbf{X}\\Gamma_2)\\Gamma_1 + \\epsilon $$\n",
    "\n",
    "where $\\mathbf{X}$ is $N\\times p_2$, $\\Gamma_2$ is $p_2 \\times p_1$, and $\\Gamma_1$ is $p_1 \\times 1$.  (We're ignoring intercepts and \"bias\" terms for simplicity.)\n",
    "We can re-write the contents of the outer activation $a(.) \\equiv \\mathbf{V}$, giving us\n",
    "\n",
    "$$y = \\mathbf{V}\\Gamma_1 + \\epsilon$$\n",
    "\n",
    "This is a linear regression!  Basically, a neural network is nothing more than an algorithm for selecting a good set of regressors $\\mathbf{V}$ from your input data $\\mathbf{X}$.  This is most obvious in simple feedforward neural networks, but it's also true for other more exotic forms.  \n",
    "\n",
    "Why would you want to fit a neural network instead of an linear regression?  Typically it's because you don't know the function mapping $\\mathbf{X}$ to $y$ , except that it isn't linear.\n",
    "\n",
    "But what if you know *something* about the functional form mapping $\\mathbf{X}$ to $y$.  Maybe you know that there is a strong linear relationship between some $z = f(x), x \\in \\mathbf{X}$.  For example, say that one of your $\\mathbf{X}$'s is temperature, and you know that there are declines in crop yield ($y$) for each hour spent above 30$^\\circ$C.  If you didn't know about neural networks, you could just regress $z$ on $y$.  But, you wouldn't be using the information in the rest of $\\mathbf{X}$.  Plus, what if the true relationship between $z$ and $y$ isn't linear, even through there is a significant OLS slope coeffcient?  \n",
    "\n",
    "The basic idea of a semiparametric neural net (SNN) is to start with what you suppose to be a reasonable linear regression specification, and let the neural network (1) fix any mis-specification, and (2) use all of the data that you're not quite sure what to do with.  It's really quite simple:\n",
    "\n",
    "$$y = \\mathbf{Z}\\beta + a(\\mathbf{X}\\Gamma_2)\\Gamma_1 + \\epsilon $$\n",
    "\n",
    "(in the one-layer case), where $\\mathbf{Z}$ are variables that you believe might have a linear relationship with $y$. In practice, $\\mathbf{Z}$ and $\\mathbf{X}$ can be anything, including the same thing. \n",
    "\n",
    "### The \"OLS trick\"\n",
    "\n",
    "The \"OLS trick\" short-circuits backpropagation by replacing the top-level weights with OLS (or ridge) estimates, which have a closed-form solution.  You simply replace $\\Gamma_1$ and $\\beta$ with the relevant parts of $\\mathbf{(W}^T\\mathbf{W}+\\lambda I)^{-1}\\mathbf{W}^T\\mathbf{y}$, where $\\mathbf{W} \\equiv [\\mathbf{Z, V}]$.  Why would you want to do this?  \n",
    " - It can speed model fitting\n",
    " - You might want some of the coefficients to have a causal interpretation, which you think can be done by controlling for other variables (in this case you'd think carefully about what you're penalizing). (I won't focus on causal inference in this notebook, though in principal this class of models could be used for certain problems that involve causal inference, under certain assumptions.  See [here](https://arxiv.org/abs/1702.06512) for more on this.)\n",
    " - You want to try to do approximate inference on marginal effects from parametric terms\n",
    " \n",
    "The rest of the notebook demonstrates an implementation of all of this in Python, using Keras.\n",
    "\n",
    "### An example\n",
    "\n",
    "We'll do a simplified version of the yield model for corn, fitting the following model:\n",
    "\n",
    "$$\n",
    "yield_{it} = \\alpha_i + TimeTrend_{t}\\gamma + \\displaystyle\\sum_h \\beta_h TempBin(h)_{it} + f(Weather_{it}, Soil_i) + \\epsilon\n",
    "$$\n",
    "\n",
    "Basically, county-level yield is a linear function of a time trend and the proportion of time spent in each 1$^\\circ$C temperature bin (from -1 to 40+), and some unknown nonlinear function of weather and soil.  I account for unobserved cross-sectional heterogeneity using a county-level fixed effect.\n",
    "\n",
    "(This specification is fairly arbitrary.  I could improve efficiency with random rather than fixed effects.  I could use a different specification for time and temperature, and put precip in the parametric part of the model.  This specification isn't neccesarially what I'd use for any particular purpose, other than demonstrating the method.)\n",
    "\n",
    "[Here](https://www.dropbox.com/s/cbvgghuoyd9kpve/IA_corn.csv?dl=0) is a dataset of yield and weather over time for the state of Iowa.  I'll be using that for this demo.\n",
    "\n",
    "We'll get started by loading the packages that we'll be using, load the data, and then partition the covariates into those that we'll be using parametrically, nonparametrically, as fixed effects, and as the outcome"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using TensorFlow backend.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3754, 43)\n",
      "(3754, 1754)\n",
      "99\n",
      "[1979 1980 1981 1982 1983 1984 1985 1986 1987 1989 1990 1991 1992 1993\n",
      " 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007\n",
      " 2008 2009 2010 2011 2012 2013 2014 2015 2016 2017]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.stats import multivariate_normal\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "import pandas as pd\n",
    "import re\n",
    "import warnings\n",
    "from patsy import dmatrices\n",
    "from keras.layers import Input, Dense, concatenate, Dropout, BatchNormalization\n",
    "from keras.models import Model\n",
    "from keras.callbacks import EarlyStopping, ModelCheckpoint\n",
    "from keras import regularizers, optimizers\n",
    "import os\n",
    "from keras import backend as K\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "\n",
    "datadir = \"/home/andrew/Dropbox/USDA/random/\" # set your own directory here\n",
    "dat = pd.read_csv(datadir + \"IA_corn.csv\")\n",
    "dat = dat.drop(columns = 'Unnamed: 0')\n",
    "# partition the covariates\n",
    "y = dat[[\"yield\"]]\n",
    "fe = dat[[\"fips\"]]\n",
    "parametric = pd.concat([dat[[\"year\"]], dat.loc[:,\"tbinneg1\":\"tbin40\"]], axis = 1)\n",
    "pattern = 'jday|soil' # all of the daily weather variables have the string \"jday\" in them, likewise \"soil\" for soil\n",
    "idx = [i for i in range(len(dat.columns)) if re.search(pattern, dat.columns[i])]\n",
    "nonparametric = dat.iloc[:,idx]\n",
    "\n",
    "# look at the dimensions\n",
    "print(parametric.shape)\n",
    "print(nonparametric.shape)\n",
    "print(len(dat.fips.unique()))\n",
    "print(dat.year.unique())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The weather data is very wide -- it's daily measurements of min and max temp and relative humidity, precip, insolation, and average win velocity, pushed to the county level.  The source is [here](https://climate.northwestknowledge.net/MACA/data_portal.php).\n",
    "\n",
    "Let's start by seeing how far a simple fixed-effects model will take us.  Turns out it'll take us most of the way, but let's pretend that we don't know that in advance.  It won't be true in every application.\n",
    "\n",
    "First we'll define functions to take group-level means of data, run OLS (or L2-penalized OLS), and another to get back estimates of fixed effects.  For those who haven't taken (or don't remember) intro econometrics, you can estimate a model of the form \n",
    "\n",
    "$$\n",
    "y_{it} = \\alpha_i + \\mathbf{X}_{it}\\beta +\\epsilon_{it}\n",
    "$$\n",
    "by subtracting off group level means\n",
    "$$\n",
    "y_{it} - \\bar{y}_i = \\alpha_i - \\bar{\\alpha}_i + (\\mathbf{X}_{it} - \\bar{\\mathbf{X}}_{i})\\beta +\\epsilon_{it} - \\bar{\\epsilon}_{it} \n",
    "$$\n",
    "leaving\n",
    "$$\n",
    "ydm_{it} = \\mathbf{Xdm}_{it}\\beta +\\epsilon dm_{it} \n",
    "$$\n",
    "$\\beta$ is invariant to this transformation because $f(\\mathbf{X})\\equiv \\mathbf{X}\\beta$ is a linear function.  A big motivation for using neural nets for semiparametric predictive panel models is that they yield a linear representation of the data, in a way that tree-based methods (for example) don't.  \n",
    "\n",
    "To recover the $\\alpha$'s, you can simply compute\n",
    "$$\n",
    "\\hat\\alpha_i = \\bar{y}_i - \\bar{\\mathbf{X}}_{it}\\hat\\beta\n",
    "$$\n",
    "\n",
    "I implement all of that below.  First I'll define functions to do it all:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def groupmeans(G, X):\n",
    "    gm = pd.concat([G, X], axis = 1).groupby(G.columns[0], axis = 0).mean()\n",
    "    GM = np.empty(shape = X.shape)\n",
    "    for j in np.unique(G.values):\n",
    "        idx = np.where(j == G.values)[0]\n",
    "        toput = np.repeat(gm.loc[j,].values.reshape(1,GM.shape[1]), repeats = len(idx), axis = 0)\n",
    "        GM[idx,:] = toput\n",
    "    return GM\n",
    "\n",
    "# ols function\n",
    "def ols(X, y, lam = 0, parapen = None):\n",
    "    # default value is unpenalized\n",
    "    if parapen is None:\n",
    "        parapen = np.zeros(X.shape[1])\n",
    "    # sanity check\n",
    "    assert (len(parapen) == X.shape[1]), \"wrong length for parapen\"\n",
    "    penmat = np.diag(parapen*lam) \n",
    "    b = np.dot(np.linalg.inv(np.dot(X.T, X)+ penmat), np.dot(X.T, y))\n",
    "    return(b)\n",
    "\n",
    "# function to get FE estimates\n",
    "def getfe(b, GMx, GMy):\n",
    "    fe = GMy - np.dot(GMx, b)\n",
    "    return fe\n",
    "\n",
    "# mse function\n",
    "def mse(x, y):\n",
    "    return ((x-y)**2).mean()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we'll partition our data into a training and test set.  If we were doing this for real, we'd get efficiency gains using bootstrap aggregation.  It'll suffice here to use the period up to 2012 to predict 2012 onward."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MSE is 443.1393808958281\n"
     ]
    }
   ],
   "source": [
    "# generate training and testing indices\n",
    "tr_idx = np.where(dat.year < 2012)[0]\n",
    "te_idx = np.where(dat.year >= 2012)[0]\n",
    "\n",
    "# fit parametric model\n",
    "xbar = groupmeans(fe.loc[tr_idx,:], parametric.loc[tr_idx,:])\n",
    "ybar = groupmeans(fe.loc[tr_idx,:], y.loc[tr_idx,:])\n",
    "b = ols(X = np.array(parametric.loc[tr_idx,:]) - xbar, \n",
    "        y = np.array(y.loc[tr_idx,:]) - ybar)\n",
    "alpha = getfe(b, xbar, ybar)\n",
    "\n",
    "# get error for test set\n",
    "Xb = parametric.loc[te_idx,:] @ b\n",
    "mergeleft = pd.concat([fe.loc[te_idx,:], y.loc[te_idx,:],Xb], axis = 1)\n",
    "mergeleft.columns = [\"fips\", \"yield\", \"Xb\"]\n",
    "mergeright = pd.concat([fe.loc[tr_idx,:], pd.DataFrame(alpha)], axis = 1).groupby(by = \"fips\").mean()\n",
    "mergeright.columns = [\"alpha\"]\n",
    "predframe = mergeleft.merge(mergeright, left_on = \"fips\", right_on = mergeright.index)\n",
    "predframe['pred'] = predframe.Xb + predframe.alpha\n",
    "\n",
    "print(\"MSE is \" + str(mse(predframe.pred, predframe['yield'])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, we can predict corn yields plus or minus about 21 bushels/acre with this model.\n",
    "\n",
    "Now we're ready to define a semiparametric neural net.  I'll define a constructor function to define a model based on data shapes, hyperparameters, and model architecture.  The point of the constructor function is to facilitate experimentation with different sorts of networks.  \n",
    "\n",
    "If you haven't seen a basic neural net implemented in Keras before, I'm likely to lose you from here downwards.  Here's a [fairly good basic tutorial](https://machinelearningmastery.com/keras-functional-api-deep-learning/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to build a SNN for a fixed effects model\n",
    "def SNN_constructor(P, NP, FE, lam, nlayers = 10, nnodes = None, drate = 0,\n",
    "                    NP_out_dim = 1, use_BN = True, LR = .001):\n",
    "    # default to 25 nodes per layer if not specified\n",
    "    if nnodes is None:\n",
    "        nnodes = np.repeat(25, nlayers)\n",
    "    if nnodes[-1] != NP_out_dim:\n",
    "        nnodes[-1] = nnodes[-1]**0 * NP_out_dim\n",
    "    # nonparametric part\n",
    "    input_np = Input(shape = (NP.shape[1],), name = \"nonparametric_input\")\n",
    "    # first layer\n",
    "    lay = Dense(nnodes[0], activation = \"relu\", \n",
    "                kernel_regularizer=regularizers.l2(lam))(input_np)\n",
    "    if use_BN == True:\n",
    "        lay = BatchNormalization()(lay)\n",
    "    # subsequent layers\n",
    "    for i in range(1,nlayers):\n",
    "        lay = Dropout(rate = drate)(lay)\n",
    "        if i == max(range(nlayers)):\n",
    "            layname = \"np_top\"\n",
    "        else:\n",
    "            layname = \"layer_\" + str(i)\n",
    "        lay = Dense(nnodes[i], activation = \"relu\", name = layname,\n",
    "                    kernel_regularizer=regularizers.l2(lam))(lay)\n",
    "        if use_BN == True:\n",
    "            lay = BatchNormalization()(lay)\n",
    "    # parametric part\n",
    "    input_p = Input(shape = (P.shape[1],), name = \"parametric_input\")\n",
    "    # FE part\n",
    "    input_fe = Input(shape = (FE.shape[1],), name = \"FE_input\")\n",
    "    fe_lay = Dense(1, use_bias = False, name = \"FE_output\")(input_fe)\n",
    "    # concatenate them\n",
    "    conc = concatenate([fe_lay, input_p, lay], name = \"conc\")\n",
    "    # output\n",
    "    out = Dense(1, name = \"top\", use_bias = False,\n",
    "                kernel_regularizer=regularizers.l2(lam))(conc)\n",
    "    # compile\n",
    "    model = Model(inputs = [input_fe, input_p, input_np], outputs = out)\n",
    "    optimizer = optimizers.Nadam(lr = LR)\n",
    "    model.compile(optimizer=optimizer, loss='mean_squared_error')\n",
    "    return model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A few things to note here.  First, I define the nonparametric, parametric, and FE as separate inputs, and concatenate them.  I haven't told you yet what shape `FE` is yet -- they'll be represented by dummy coding, but we won't just rely on the neural net to determine their value (more later).  Next, the fixed effects feed into a layer with output shape `(,1)` -- I don't add them directly into the top layer as I do with the parametric covariates.  Finally, the model will expect three tensors (all 2D matrices) as inputs.\n",
    "\n",
    "Here are the main functions for implementing the OLS trick in the context of a fixed effects regression:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# extraction of top layer function\n",
    "def get_top_layer(model, trdat, \n",
    "                  conclayername = \"conc\", toplayername = \"top\"):\n",
    "    # define model to extract from \n",
    "    extr = Model(inputs = model.input,\n",
    "                 outputs = model.get_layer(conclayername).output)\n",
    "    toplayer = extr.predict(trdat)\n",
    "    return(toplayer)\n",
    "\n",
    "# put weights back into keras model \n",
    "def replace_top_weights_fe(model, beta, alpha, toplayername = \"top\", felayername = \"FE_output\"):\n",
    "    # beta\n",
    "    old = model.get_layer(toplayername).get_weights()[0]\n",
    "    for b in range(len(beta)):\n",
    "        if b == 0: # set the FE layer to one\n",
    "            old[b] =(old[b]**0).reshape(old[b].shape) \n",
    "        else:\n",
    "            old[b] = (old[b]**0)*beta[b-1].reshape(old[b].shape)\n",
    "    model.get_layer(toplayername).set_weights([old])\n",
    "    # alpha\n",
    "    old = model.get_layer(felayername).get_weights()[0]\n",
    "    for a in range(len(alpha)):\n",
    "        old[a] = (old[a]**0)*alpha.values[a].reshape(old[a].shape)\n",
    "    model.get_layer(felayername).set_weights([old])\n",
    "\n",
    "# olstrick wrapper function\n",
    "def fe_olstrick(model, trdat, y, fe, lam, parapen,\n",
    "                conclayername = \"conc\", toplayername = \"top\"):\n",
    "    toplayer = get_top_layer(model, trdat, conclayername, toplayername)[:,1:] # remove FE row\n",
    "    xbar = groupmeans(fe, pd.DataFrame(toplayer))\n",
    "    ybar = groupmeans(fe, y)\n",
    "    newb = ols(X = toplayer - xbar, \n",
    "            y = np.array(y) - ybar,\n",
    "            lam = lam, parapen = parapen)\n",
    "    alpha = getfe(newb, xbar, ybar)\n",
    "    alpha = pd.concat([fe, pd.DataFrame(alpha)], axis = 1).groupby(by = \"fips\").mean()\n",
    "    replace_top_weights_fe(model, newb, alpha)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first function takes a compiled SNN model and extracts it's top layer -- the concatenation $[\\alpha, \\mathbf{Z, V}]$ following our notation earlier. \n",
    "\n",
    "The second one replaces the top level weights $\\beta$ and the fixed effects estimates $\\alpha$ with ones you got in the next function. \n",
    "\n",
    "The final one actually does the OLS trick (in the FE context).  \n",
    "\n",
    "These functions all depend on the layers of the SNN being identifiable via names that I've applied to them.  Also, this code would probably break if the fixed effects were out of order.  NB: this code isn't seemless!  There are plenty of corner cases that I'd identify and protect against if this were to go into any sort of production setting.\n",
    "\n",
    "Next, we'll prep the data for Keras:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prep the data for Keras\n",
    "FE = np.asarray(dmatrices(\"1~ C(fips)-1\", fe)[1]) # this is the dummy expansion (or \"one-hot encoding\" if you prefer)\n",
    "P = np.array(parametric)\n",
    "NP = np.array(nonparametric)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is the function that puts it all together:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_with_FE_olstrick(y, P, NP, FE, # input data\n",
    "                         tr_idx, te_idx, # train and test indices\n",
    "                         lam, drate, # tunable hyperparameters\n",
    "                         batchsize, epochs, LR, # training parameter\n",
    "                         nnodes, NP_out_dim, nlayers, # architecture\n",
    "                         weight_location): # path to weights\n",
    "    trdat = [FE[tr_idx,:],\n",
    "             P[tr_idx,:],\n",
    "             NP[tr_idx,:]]\n",
    "\n",
    "    tedat = [FE[te_idx,:],\n",
    "             P[te_idx,:],\n",
    "             NP[te_idx,:]]\n",
    "    print(\"starting constructor\")\n",
    "    model = SNN_constructor(P = P,\n",
    "                        NP = NP,\n",
    "                        FE = FE,\n",
    "                        lam = lam,\n",
    "                        nlayers = nlayers,\n",
    "                        nnodes = np.repeat(nnodes, nlayers),\n",
    "                        NP_out_dim = NP_out_dim,\n",
    "                        drate= drate,\n",
    "                        LR = LR)\n",
    "    global bestweights\n",
    "    if bestweights is not None:\n",
    "        model.set_weights(bestweights)\n",
    "    best_mse = mse(model.predict(tedat), y.loc[te_idx,:])\n",
    "    callbacks = [EarlyStopping(monitor='val_loss', patience=4),\n",
    "                 ModelCheckpoint(filepath=weight_location+'temp_weights.h5', monitor='val_loss', save_best_only=True)]\n",
    "    print(\"fitting model\")\n",
    "    model.fit(x = trdat, \n",
    "              y = y.loc[tr_idx,:],\n",
    "              validation_data = (tedat, y.loc[te_idx,:]),\n",
    "              callbacks = callbacks,\n",
    "              batch_size  = batchsize,\n",
    "              verbose = 0,\n",
    "              epochs = epochs)\n",
    "    print(\"loading weights\")\n",
    "    if \"temp_weights.h5\" in os.listdir(weight_location):\n",
    "        model.load_weights(weight_location+\"temp_weights.h5\")    \n",
    "    print(\"doing OLStrick\")\n",
    "    parapen = np.repeat(1, P.shape[1]+NP_out_dim) # we're not leaving some terms unpenalized here.  we could if we wanted.\n",
    "    fe_olstrick(model, trdat, y.loc[tr_idx,:], fe.loc[tr_idx,:], lam, parapen)\n",
    "    current_mse = mse(model.predict(tedat), y.loc[te_idx,:])\n",
    "    improvement = best_mse - current_mse\n",
    "    if current_mse.values < best_mse.values:\n",
    "        bestweights = model.get_weights()\n",
    "    # output\n",
    "    out = pd.DataFrame(data = {\"lam\":lam,\n",
    "                        \"drate\":drate,\n",
    "                        \"batchsize\":batchsize,\n",
    "                        \"LR\":LR,\n",
    "                        \"improvement\":improvement.values[0]}, index = [0])\n",
    "    return {\"model\": model, \"out\":out}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what it does:\n",
    "1.  It defines lists of the inputs for the training and testing sets.\n",
    "2.  It constructs the SNN with hyperparameters of your choice\n",
    "3.  It looks in the global environment for weights from previous runs, and loads them if they are there.  \n",
    "4.  It sets the step size (\"learning rate\") of the optimizer\n",
    "5.  It establishes the baseline mean squared error\n",
    "6.  Specifies the callbacks to enable early stopping\n",
    "7.  Fits the model \n",
    "8.  Loads the best weights from the fitting run\n",
    "9.  Does the OLS trick\n",
    "10.  Assesses fit, and improvement.  If the model improves, save the best weights (to the global env)\n",
    "11.  Output the hyperparameters and the improvement, along with the actual model\n",
    "\n",
    "I set it up this way so that I can do hyperparameter optimization along the lines of the approach I took in the ERL paper.  Rather than seeking a single optimal hyperparameter, I continually try new hyperparameters, saving the resultant weights when I get improvement.  I start with 10 random tries, and then I start trying to pick the best ones using a random forest.  Here's how I implement that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    74658.953148\n",
      "dtype: float64\n",
      "        lam  drate  batchsize     LR   improvement\n",
      "0  0.000215    0.7         64  0.001 -19840.017848\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    377.814817\n",
      "dtype: float64\n",
      "        lam  drate  batchsize    LR  improvement\n",
      "0  0.004642    0.1       2048  0.01  7353.305644\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    9564.323161\n",
      "dtype: float64\n",
      "     lam  drate  batchsize        LR  improvement\n",
      "0  0.001    0.3        256  0.004642 -9186.508343\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    569.029702\n",
      "dtype: float64\n",
      "        lam  drate  batchsize        LR  improvement\n",
      "0  2.154435    0.8         64  0.000215  -191.214884\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    769.079445\n",
      "dtype: float64\n",
      "     lam  drate  batchsize    LR  improvement\n",
      "0  0.001    0.7        256  0.01  -391.264628\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    441.823516\n",
      "dtype: float64\n",
      "        lam  drate  batchsize       LR  improvement\n",
      "0  0.004642    0.5        512  0.00001   -64.008698\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    363.352305\n",
      "dtype: float64\n",
      "       lam  drate  batchsize        LR  improvement\n",
      "0  0.00001    0.3       2048  0.000464    14.462512\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    311.307842\n",
      "dtype: float64\n",
      "        lam  drate  batchsize        LR  improvement\n",
      "0  0.000046    0.3       1024  0.000022    52.044463\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    8.281185e+09\n",
      "dtype: float64\n",
      "       lam  drate  batchsize        LR   improvement\n",
      "0  0.00001    0.2        128  0.002154 -8.281185e+09\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    1071.18126\n",
      "dtype: float64\n",
      "        lam  drate  batchsize        LR  improvement\n",
      "0  2.154435    0.5         64  0.000046  -759.873418\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    349.864933\n",
      "dtype: float64\n",
      "        lam  drate  batchsize        LR  improvement\n",
      "0  0.021544    0.4        128  0.000046   -38.557091\n",
      "max pred was -82812233.94618163 improvement\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    530.323985\n",
      "dtype: float64\n",
      "        lam  drate  batchsize        LR  improvement\n",
      "0  0.464159    0.7         64  0.000022  -219.016143\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    103612.512855\n",
      "dtype: float64\n",
      "        lam  drate  batchsize    LR    improvement\n",
      "0  0.000046    0.7       1024  0.01 -103301.205013\n",
      "max pred was -9569.156333652349 improvement\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    415.415301\n",
      "dtype: float64\n",
      "        lam  drate  batchsize       LR  improvement\n",
      "0  0.464159    0.5         64  0.00001  -104.107459\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    478.472291\n",
      "dtype: float64\n",
      "        lam  drate  batchsize        LR  improvement\n",
      "0  0.000215    0.0        256  0.002154  -167.164449\n",
      "max pred was -2211.0308382128846 improvement\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    380.583579\n",
      "dtype: float64\n",
      "   lam  drate  batchsize       LR  improvement\n",
      "0  0.1    0.5         64  0.00001   -69.275737\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    308.21337\n",
      "dtype: float64\n",
      "        lam  drate  batchsize      LR  improvement\n",
      "0  0.004642    0.2        128  0.0001     3.094473\n",
      "max pred was -141.71343913086844 improvement\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    417.30847\n",
      "dtype: float64\n",
      "   lam  drate  batchsize       LR  improvement\n",
      "0  0.1    0.4        128  0.00001  -109.095101\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    376.567235\n",
      "dtype: float64\n",
      "        lam  drate  batchsize      LR  improvement\n",
      "0  0.021544    0.1       1024  0.0001   -68.353865\n",
      "max pred was -70.39143291757762 improvement\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    456.410466\n",
      "dtype: float64\n",
      "   lam  drate  batchsize      LR  improvement\n",
      "0  0.1    0.4        128  0.0001  -148.197097\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    591.99682\n",
      "dtype: float64\n",
      "        lam  drate  batchsize        LR  improvement\n",
      "0  2.154435    0.3         64  0.000022   -283.78345\n",
      "max pred was -103.20522295167693 improvement\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    420.229896\n",
      "dtype: float64\n",
      "   lam  drate  batchsize       LR  improvement\n",
      "0  0.1    0.5        128  0.00001  -112.016526\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    78373.892845\n",
      "dtype: float64\n",
      "     lam  drate  batchsize    LR   improvement\n",
      "0  0.001    0.8         64  0.01 -78065.679476\n",
      "max pred was -8608.666799437768 improvement\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    392.497087\n",
      "dtype: float64\n",
      "   lam  drate  batchsize       LR  improvement\n",
      "0  0.1    0.5        512  0.00001   -84.283717\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    447.318214\n",
      "dtype: float64\n",
      "        lam  drate  batchsize    LR  improvement\n",
      "0  0.004642    0.4        128  0.01  -139.104844\n",
      "max pred was -1645.250578387694 improvement\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    382.106347\n",
      "dtype: float64\n",
      "        lam  drate  batchsize        LR  improvement\n",
      "0  0.021544    0.4        128  0.000046   -73.892978\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    344.094513\n",
      "dtype: float64\n",
      "     lam  drate  batchsize        LR  improvement\n",
      "0  0.001    0.4        128  0.000215   -35.881143\n",
      "max pred was -31.010815576731456 improvement\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    4450.638411\n",
      "dtype: float64\n",
      "   lam  drate  batchsize     LR  improvement\n",
      "0  0.1    0.3       2048  0.001 -4142.425042\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    299.803468\n",
      "dtype: float64\n",
      "        lam  drate  batchsize        LR  improvement\n",
      "0  0.000215    0.3        256  0.000022     8.409902\n",
      "max pred was -961.1045285127497 improvement\n",
      "starting constructor\n",
      "fitting model\n",
      "loading weights\n",
      "doing OLStrick\n",
      "yield    2584.961142\n",
      "dtype: float64\n",
      "        lam  drate  batchsize       LR  improvement\n",
      "0  0.004642    0.5       2048  0.00001 -2285.157675\n"
     ]
    }
   ],
   "source": [
    "# grid of hyperperameters\n",
    "hpgrid = np.array([(lam, drate, batchsize, LR) \\\n",
    "             for lam in 10**np.linspace(-5, 1, num = 10) \\\n",
    "             for drate in np.linspace(0, .8, num = 9) \\\n",
    "             for batchsize in 2**np.linspace(6, 11, num = 6) \\\n",
    "             for LR in 10**np.linspace(-5, -2, num = 10) \\\n",
    "             ])\n",
    "hpgrid = pd.DataFrame(hpgrid, columns = [\"lam\", \"drate\", \"batchsize\", \"LR\"])\n",
    "\n",
    "trdat = [FE[tr_idx,:],\n",
    "         P[tr_idx,:],\n",
    "         NP[tr_idx,:]]\n",
    "\n",
    "tedat = [FE[te_idx,:],\n",
    "         P[te_idx,:],\n",
    "         NP[te_idx,:]]\n",
    "\n",
    "global bestweights\n",
    "bestweights = None\n",
    "# 10 initial random experiments\n",
    "BHOframe = pd.DataFrame(columns=[\"lam\", \"drate\", \"batchsize\", \"LR\", \"improvement\"])\n",
    "for i in range(10):\n",
    "    ridx = np.random.choice(range(hpgrid.shape[0]), size = 1)\n",
    "    output = fit_with_FE_olstrick(y, P, NP, FE, # input data\n",
    "                     tr_idx, te_idx, # train and test indices\n",
    "                     lam = hpgrid.lam[ridx].values[0], drate = hpgrid.drate[ridx].values[0], # tunable hyperparameters\n",
    "                     batchsize = int(hpgrid.batchsize[ridx]), epochs = 100, LR = hpgrid.LR[ridx].values[0], # training parameter\n",
    "                     nnodes=25, NP_out_dim=25, nlayers=10, # architecture\n",
    "                     weight_location = '/home/andrew/Dropbox/USDA/random/')\n",
    "    print(mse(output['model'].predict(tedat), y.loc[te_idx,:]))\n",
    "    print(output[\"out\"])\n",
    "    K.clear_session()\n",
    "    BHOframe = BHOframe.append(output[\"out\"])\n",
    "    \n",
    "# every other iteration, alternate random choices and ones predicted (with a RF) to be optimal\n",
    "rf = RandomForestRegressor(n_estimators = 100, \n",
    "                           random_state = 1, \n",
    "                           max_features = .33,\n",
    "                           n_jobs = -1)\n",
    "for i in range(20):\n",
    "    if i % 2 == 1:\n",
    "        X = BHOframe.iloc[1:, :-1]\n",
    "        X[\"iter\"] = np.array(range(X.shape[0]))\n",
    "        \n",
    "        improvement = BHOframe.improvement[1:]\n",
    "        \n",
    "        # Train the model on training data\n",
    "        rf.fit(X, improvement)\n",
    "        pgrid = hpgrid.copy()\n",
    "        pgrid[\"iter\"] = X.shape[0]\n",
    "        pred = rf.predict(pgrid)\n",
    "        ridx = np.where(pred == max(pred))[0]\n",
    "        if (len(ridx)>1):\n",
    "            ridx = np.random.choice(ridx, size = 1)\n",
    "        print(\"max pred was \" + str(max(pred)) + \" improvement\")\n",
    "    else:\n",
    "        ridx = np.random.choice(range(hpgrid.shape[0]), size = 1)\n",
    "    \n",
    "    output = fit_with_FE_olstrick(y, P, NP, FE, # input data\n",
    "                     tr_idx, te_idx, # train and test indices\n",
    "                     lam = hpgrid.lam[ridx].values[0], drate = hpgrid.drate[ridx].values[0], # tunable hyperparameters\n",
    "                     batchsize = int(hpgrid.batchsize[ridx]), epochs = 100, LR = hpgrid.LR[ridx].values[0], # training parameter\n",
    "                     nnodes=25, NP_out_dim=25, nlayers=10, # architecture\n",
    "                     weight_location = '/home/andrew/Dropbox/USDA/random/')\n",
    "    print(mse(output['model'].predict(tedat), y.loc[te_idx,:]))\n",
    "    print(output[\"out\"])\n",
    "    BHOframe = BHOframe.append(output[\"out\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Note about Keras:  I've found that repeatedly fitting different models under the same kernel leads to slowdowns after the first few models, even when using `K.clear_session()`.  One way around this can be to put you model-fitting script in something like `run_models.py`, and then repeatedly call this script from another script.  That way the OS cleans up whatever cruft is slowing the system down.  It is likely that there is a better way.  If so, I'd love to hear about it.)\n",
    "\n",
    "In this instance, the SNN only gives a modest improvement in performance over the OLS regression.  If we wanted to improve performance further, we could:\n",
    " - implement random, rather than fixed effects.  This would involve a different approach to the OLS trick, and would rely on Keras to estimate the coefficients $\\alpha$\n",
    " - implement bootstrap aggregation -- fitting the model multiple times to different training and test sets, averaging the results\n",
    " - try different architectures\n",
    " - improve the parametric component of the model\n",
    " - reduce the dimensionality of the inputs by taking their principal components\n",
    " - etc etc\n",
    "\n",
    "It's also worth mentioning that my python/Keras implementation of hyperparameter optimization didn't work as well as did my R implementation using `panelNNET`.  After a few early initial improvements, subsequent trials were entirely negative.  I'm not entirely sure why this is -- and I'd really appreciate suggestions!\n",
    "\n",
    "But, the goal here is simply to show, in broad strokes, how the method works.  We'll finish it by finalizing the model and plotting predictions against observations, and computing MSE.  We got it down to about 17.3 bushels/acre with the SNN."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "yield    17.314834\n",
      "dtype: float64\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnX+cXGV56L/PTiYyi5oNJdfCQAjy0aARkpWtpI3tFeolCgorIJHi1dp+SlttKyldGyrXBIslbaqo11u9tFJ/gJjww20QNVhBbdMbaGI2xCipIBBYUVCyqMlAJrPP/WPO2ZyZPT/eMzNnfuw+389nP5l558w5z5yZvM/7Pj9FVTEMwzCMevo6LYBhGIbRnZiCMAzDMEIxBWEYhmGEYgrCMAzDCMUUhGEYhhGKKQjDMAwjFFMQhmEYRiimIAzDMIxQTEEYhmEYoczptADNcOyxx+qiRYs6LYZhGEZPsWPHjp+q6oKk43paQSxatIjt27d3WgzDMIyeQkQecznOTEyGYRhGKKYgDMMwjFBMQRiGYRihmIIwDMMwQjEFYRiGYYTS01FMhmEY3cToznE2bNnLjyZKHD9QYGTlYoYHi50Wq2FMQRiGYbSA0Z3jXHXHbkrlCgDjEyWuumM3QM8qCTMxGYZhtIANW/ZOKQefUrnClZt2MbpzvENSNYcpCMMwjBbwo4lS6HhFldUbx1i05i5WrL+np5SFmZgMwzAiiPMpXD26m1vue5yKKjkR+ufmOHCoEnoe9f7tNbOTKQjDMAyPoEKYV8hz4NBhypXq9B6c3Lc/9gw3bds39b6KKgcOVcj1CZVJDT23T6lcYcOWvT2hIMzEZBiGwREn8/hECQUmSuUp5eDjT+633Pd46DkmJ6u7iSTGJ0o9YW4yBWEYhkG4kzmMH02UqGj4LkGBD1+ylEI+l3gef0fSzUrCFIRhGAbRTuZ64g1IVa678DSKA4XE40rlCtfcucfpup3AFIRhGAZwvMOE7sI1d+5heLDI1jVn89FVy8jn4k1O+w+Wu3YXYQrCMAwDGFm5mGTvQTL7D5anHl9z555pfowwNmzZ24Irtx6LYjIMY8aStvRFXKhqGlasv4dxR5MVuJu32o0pCMMwZiRpSl/UH9ssaZQDwLxCviXXbTVmYjIMY0YSVfoizJzjGsGUFQcOHe5KP4QpCMMwZiRRZpuw8bQr/lZTrmhX+iFMQRiGMSOJikoKG3dJbsua8YkSV4/uTjxudOc4K9bfw8ltqO1kCsIwjBnJyMrF0xLWCvkcIysXTzs2KvGt3dy0bV+skqjP9s462c4UhGEYMwp/hb164xgvmNPH/P48AhQHClx34WmhUUwuSW3t4gv37Yt8LY1fpRVkpiBE5EQRuVdEvicie0Tkvd74MSLydRH5gffvfG9cROTjIvKQiDwgIq/OSjbDMGYmYfWUnitPcv2qZWxdc3Zo9JIfkpqVkSmt+WpSidxFpPGrtIIsdxCHgStV9ZXAcuA9IvJKYA3wDVV9GfAN7znAG4GXeX+XA5/MUDbDMNpAO+3lkG6FHVQmUC2h4U/lxYEChXzz02NOhEvPPDH1+27ato/BD9497X6l8au0gszyIFT1SeBJ7/EvROT7QBG4AHidd9hngW8Cf+mNf05VFdgmIgMicpx3HsMweoxOtOBMs8IOUyZKdVKPi2qa359n4mCZPpFE30VFlWuHTwOY6h3hyv6DZUZu28W6zXt4tlTm+IECZ526gNt3jNfIHeVXaQVt8UGIyCJgELgPeElg0v8x8BLvcREI1tB9whszDKMHabe9HNKtsKOUQNIkvvbNS3hk/XlOVVtFYPCDd3Pztn28uDAn9a6kXFEmSuUph/TtO8a56IwixYFCol+lFWSeSS0iLwRuB65Q1Z9LwB6nqioiqcIHRORyqiYoFi5c2EpRDcNoIe22l0M1cqk+IzrfJxw8dJiT19xVU25DcKvMWo/f7MeflK/YOBZ5rOqR2kzBGk2NUipXuPfBp9m65uymz+VCpjsIEclTVQ43q+od3vBPROQ47/XjgKe88XEgaKw7wRurQVVvUNUhVR1asGBBdsIbhtEU7bSXByOX+uSIL0GASaqTs78Kv2LjGK/4X19tSDngncP3DQwPFtseATU+UWLZNdP9E1mQZRSTAJ8Gvq+qHwm8tBl4p/f4ncC/BMbf4UUzLQeeNf+DYfQuafIQmqE+cunAocrU5K8Q2gK0VJ5s6prB3IOwz5k1E6UyI7fuylxJZLmDWAH8T+BsERnz/s4F1gP/Q0R+ALzeew7wFeCHwEPAPwLvzlA2wzAyZniwONU4J629PE30UyfqKAV9KcHP2U7Kk9mX5xDtkgzCRhgaGtLt27d3WgzDmPWkLasdd55r7twzzV5fyOcilcvJa+5q2FzUDAI8sv68mrHRneOxPokghXyuJYrto6uWpb7XIrJDVYeSjrNMasMwmqJV5R/884Q5c+Oin1rl0yjkc7x9uXvgi1+iO7jbcV3R9+f7Yncdvh+lOFBgfn98KXArtWEYRtcSFc76/i/tTpUkl2Qqiop+apUP4KIzilw7fBoDjr0ZRMKVYxw5Ed6+fCF/c+HpbNiyNzSDu5DP8ZFLlvHI+vPYuuZszjv9uNhzZhk6bA2DDMNwJsyUFDVxHzhU4cCh6msuSXJJ4a9ROwX/fGGmqTTc++DTQHXid2H/wTLX3LnH2Uz09uULuXb4tGkJhH4Gt1LdMQTNc6M7x7l9R/LuoBdLbRiGMYMIWy2v3jhG/1y31XvSSjfOVOQS/fRcQmSSQOzuwJ9kJ1IoGReF5O8a/IzqqAzu4kBhWr0oVwd8z5XaMAxjZhE1saXp4Ry30g1LcoPqpL7u/CWxjliXiVSBNy09jpu2hVdL9SfZ4wcKLWsgJMDD151bMxZ1D/z8BpGqknKVI8tSG6YgDMNwIm5yd85KFli05q6a99SbVYK2eb8i6xUbx1i3eU+konAxsRQHClNmpDD5/Ul2ZOVirrx1V2j+RFrCVvZxE/9E6ciOJE455ESYVG0qYswFUxCGYTgRN7EpkM8J5Ur8pBqMqvcf+qaq7Y89M2WGGbl1F+W6CdpPDoNaP8bozvHEwnn5nHDg+cM1E3C9/MFz9gH1+5FCvi9Vgl1Q6QR9NwP9efqoZng3ggCXnnni1L3KEvNBGIbhxMjKxfE9E5Sa5jxvX77QOXlMgZu37ZuaSOuVg099cpjvFwlTDr6s8/vzoEQqB/8YP+Lqyk3TlRPAUflcqp4RvtK5enQ3qzeOTflu9h8sp1YOwesqcPuO8baU2rAdhGEYTgwPFtn+2DPcvG1fqDmpPKn0z53Dzg+cUzPuN+RJQjliXoojaE6K8z345quDhw5HKhyo7i5++dzhKYdz1E5k/8Eyc3PCoYRdkk9xoMDozvHI+5WG+vf7Dv+sTEs+toMwDMOZa4dP4/pVyyJfD5vc04RgunR26xOZWj0nKZPxiVJspFFxoMDRc+fEKpAgrsrBdxxv2LI3syzvLKvi+piCMIxZTJqaR8GKqVFtNMU7LkiaEMycSOKEWlGdyh5O284ziO8jeDbG9JSGsJpTrZjEozKpswptDWIKwjBmKWlKZNQfG2WG8c1EQVwznfO55A5tPr6JJU2HtnoUWL1prCUrfD+Hwc9+9k0/zU7iA4U8a9+8pC1VccMwBWEYs5Q0Hd/SVEytXzWHVXV9+/KF5PvqVv8avVoOY3yi1NQOAmqjqsJwPf8zB54PVayJjv0AYUrAD+tttCpus5iT2jBmKUkd34KhmWlW2b6PYHiwOK00x/Ve5dEV6++ZZvcvTyqq6aqcNrODSEKAly7o5wdPHUg8tlSe5IqNY9y6fR+P/qxUU4rksuULEx3V8/urO4WoirjBDnbtxBSEYcxSovIajveib8Kyml2oqLLamyy/s+/ZqXME6zFFKadnS2WuX7WMKzftSjX5N9o+NA4FHnJQDkG2PvzM1GP/81534WkMnXRMbITWL587XH1/m1qJumImJsOYpcR1fGu2CY9SnSyjTFhx7UiHB4tMptwZZLWPaPa8pXKFa+7cw/Bgka1rzubR9eeF1oNqR/OfRjAFYRizlDDb9kVnFJ1yEZphfKKU2I60HRE67WL/wTJXj+6eeh4VNTU+UXIqi95OrKOcYRgADZmVBgr52AzlKIRqw51gYbr6MteNmri6Fb/oYJICjuue1yqso5xhdBlpcg46QRqzklDtbzC29hw+umpZqhIUcKQI38TBMpctXzitzLW/u0kT1VTP/P789EipDjJRKnPVHbs569QFsWG/WTYASospCMNoA61qy5klcUldA4V8TZ2l61ctmyoWNzxY5LLlC6cpCZepOViDqZ7hwSI7P1BVQGkVhQDnnX4cG966tOlQ2FZSKle498GnY9uNQu130cmFhSkIw2gDaXIOOkWU3b84UGBs7Tns/MA5PLL+vCkndnDCGjrpGAYCk7gfVeSqJOLuw/Bgkf656QIu/YJ22x97JtNQ2Eb40URpymkdpST87yKqSdOiNikLC3M1jDaQlHPQbsJah4Y17BGOOE99B3LwmPGJUrUEt1BT6lsD/7qEoCbdh0buU6lc4eaI5kCdJKiIw+550Fkf1aQJ3Nq4NovtIAyjDcSFdbabKHMXUGP6CE7s/jFhPZjLkxrbB8KvqhrHvJhWoND4fWrX3qGP6aXOjw5pxZrvk5oSGUlZ0kmKMetdqEUxGUZGBFfp8wp5Dhw6XDORtiNaJWynEBVF49cTAvcS3a4Iya08g32b6xndOc7IbbsSGxJ1inxO2HDx0poorLCmRwBHz82Rz/XxbGl69FY9Lt+DAI+sPy+VvBbFZBgdpH6VPlEqT2uo0w7lELZTiJpwgqvVVpu+FHjy2fhz3hThrIbqSnvDxUtDV+XdQLlSm+gW1/TowKEKE6WyU7CCS6HDLHeh5oMwjAwIsx1HNdRppwylciXSJxA080St9gcKeZ4/PNlQfoJLy4V1m/fE1iPasGUvBw5l67c5em6Og4cqqc1TjSrYuOY/UX26fbKu6mo7CMPIgG5wSkddK2riC0aDhq1chepO6Kh8HwOFIzuhVjJRKoeGAvuhns2YvURqP2MY8/vzfOgtjfV6Dq7k067q434XwTId169a1taqrraDMIwMiCuE12kZopgIdF6LW7nuP1imkM/VVGbNqjRHqVzhr+54gFJ5simHs+8jACJ9A1D9bFfdsZt5KTPE61fyIysXx16nHtffRburutoOwjAyIKnWUCdlcO1QFozVj+qJ3Mh10nKwSeUwvz8/pRx830DcTqJUriBSVSru1Eo4PFhkw1uXhhbmq6fdv4s02A7CMDIguAIPs6d3UgYgNva+niRzWdx1Vm9069g24NVliusfnZaPejscmF7bKanvRHWX1OccNVUqT07LSQiu9oPRZAP9eVRximLqNJmFuYrIjcCbgKdU9VXe2DLgU8BRwGHg3ap6v4gI8DHgXOAg8Luq+p2ka1iYq2Gkw5+o/G5sFVWKEZNU8NgwgmGxUVw9uju2WY5fwM5vLtTKAn3Bz5WlGaz+mq3q6RAWotwqReIa5prlDuIzwCeAzwXG/g64RlW/KiLnes9fB7wReJn3dybwSe9fw5hRZPmf3uXawQm4ojq1cwhTDnGTtatZ5NrhI81yfjRRopDv42B5cup1v4AdhO9EDjx/uKFqseDWoKjVNHOduLyZdmRNh5GZglDVb4vIovph4MXe43nAj7zHFwCf0+p2ZpuIDIjIcar6ZFbyGUa7qZ902/2fPq4eVP314yq7Ru04goQpQqianOoJylDvhG12VxFsUNSOHUSjQQj1nzNMKcaFw2ZFu30QVwBbROTvqTrIf8MbLwKPB457whubpiBE5HLgcoCFCxdmKqxhtJI0E3QWpAm9jTpWiG6LGTRJhZXpOCrfF2lqirpe/a6ifgfiQjsUAzTnbHYttd7u2l3tjmL6Y2C1qp4IrAY+nfYEqnqDqg6p6tCCBQtaLqBhZEWncyPS1INKWzsqmLUN03MtSuVKvANaiCxn7UdTXb9qGaWUyiGKVreJyIk0lZPg+htod+2udiuIdwJ3eI9vBV7jPR4HTgwcd4I3Zhgzhk4X7EsTeps2TLfpHtZKYjnrDVv2tqz43ouPyqduchTHpGpTu0CX30AnwmHbrSB+BPx37/HZwA+8x5uBd0iV5cCz5n8wZhppJt0smsQkVQ5t9FhwWwEPFPKJdYVgumnK/+xx10g72fu1kKJI22SoWSUf9tvI90lba3eFkZkPQkRuoRqhdKyIPAGsBf4A+JiIzAGew/MlAF+hGuL6ENUw13dlJZdhdArX3IgsndlpMnHTHJvkBC7kc6w7fwkAV4Q4qqMI+mjirqFUzUaOictAdJ8KAZa/dD5bH37G+VzNruy7IW8mDCv3bRhdRlTMfitj7F1xDcsNizbyJ+D6qKe0OQl+OWuXiKZcn1BJoSWilISfI+LCQCHP2Nr2FGBsFd2QB2EYRgN02pnt47qT8ZVIqVxJTL6D8C5qcfSJMLpzfOpcV27aFTp5FxvIm1CqpTjqHeiuyiG4M5qJWC0mw+gyOu3M9nHpo10fvRSXfOdT799IoqLKyK27ppTEhy9ZOq1OUj5X7dT2bANJdWnLe/j+iU75BdqJ7SAMo8tI6lPcqmzspPO47GQaze0I+jdcTE7lSWXd5j1Hzlm/wPeeR/kp5vfnea7cWB8Ln3Z0AOw2TEEYRpdR77Cc5xWyW71xjHWb9ziVYEia/MPMR6s3jrH9sWem2n66lCxvhTnM1eQ0USoz+MG7UWVaGe3yZLWjW9i5hOouodrsqJLKkR1ktikHMBOTYXQlweSw5w9Psv9geap1aX2F0TizT1Rby7CVvwI3B9p+uoTlNmsOC/ovXMxN+w+WI30MvpJ7wZy+qVLjQSf0RKncsHLok/bWQOoWTEEYRhfTSAkGF99BXLc5/ziXXIionskTBw+x7Jq7Y/M4krKvG8FXohNeue5WxWg2qlh6HTMxGUYX00gJhjizj79ij5vvgu9PyoXwX1u3eU/Nyv7AoQoQH/3UbPZ1HAotK8sxm7EdhGF0MY2UYIh6z0B/vmbFHsU8hy5oQYYHixz9gvi1Zv0OBtofttsMLp3hZiKmIAyjiznr1PCClEfPzaUy+xTyOVRxWrGnrDIBuE32/g7GLyHS18iFOsSblh7XaRE6gpmYDKOLuffBp0PHB/rnsueD4VnVUWUbwnoxhDERkRcQFxnl0m9BqW1B6pqM1g1EfQ8zHVMQhtHFNBpGGuY7iGsfGiTMRBUXFjt00jEceP5w4nkh3hEtAoU5fZTKk8wr5BvuJJcF4xMlTrnqK1RUmd9DPaWbxUxMhtHFtDKrOiriKIgQXnguKiz2pm37GLlt17TJ/Oi5udR2+2rJb+H6VctYd/4SCvnump78HY8fahsVQjyT6K5vwDCMGtL2ZYgjGLYaxW+cckzoajhux1KflwFVE9jY2nNSl+EulStcc+ceRm7dNS0KqdVNfnwGCvmpUN6BQn5aGY8kwhzwMwUzMRlGF9PqMtC+6SmqvMXOfROh70vb09lXKI30go6qjZRVLsKzpXJNNdZg61RXeikiKw22gzCMLsfPqn5k/XlsXXN2S+zdURPawfIkV4/unjaedsfim8DCdkCdil2KMlnVm+v8+x2300o6x0whVkGIyJ0isjnqr11CGobRWuImtFvue3za2PBgkaPnJneDg2onNF+hhGVjX7Z8oVNnuTCilEvQTBTVDe6ofI58iJ3qwPOHQ30ILj4b6Ewr0HaRZGL6e+/fC4FfBW7ynl8K/CQroQzDyJaRlYsjO7tFhZ9+6C2nMXLbrlCfQxC/cB4cMWkFdz0r1t/TUAZ1WN8Gn6CZ6OQ1d4Ues/9gOdS/MFEqh2Z6B8174xOlqV4XsymKKVZBqOq3AETkw3Xdh+4UEWvlZhg9yvBgkdWbxgjTBVEr8PoJM464NqmN2OvzOWHtm5dEXlupKp6RlYsj/R45kUjlFlWiPE3b1ZmIqw/iaBF5qf9ERE4Gjs5GJMPoPYIZwlHF6TpJmHyXnbkw9NhLzzwx8hy+szxKiQQJi+4Z3TneUAb1houXMjxYZGTl4lAzERxRSmeduiA08ispMW+mOpqbwVVBrAa+KSLfFJFvAfcCV2QnlmH0Di7ltYPHtluRRMk3dNIxvH35wqnJPifC25cvnOoHEXcO1yzo4EreP0fYe/M5ISrtIehcHh4s8sKjog0fpXKFex98OrQKbZLTud4v0+1Kvx04hbmq6tdE5GXAqd7Qg6r6fHZiGUbv4NpVzbXHczvl27rm7CmFcPXobm6573Fu2raPnAiXnnni1GuNVl4N7jTizlGpKFG1V0vlSUZu3QVU71NSi9DxidJU86D6+xrVmKje0dyp76rbcNpBiEg/MAL8iaruAhaKyJsylcwwegTXchgufRqywEW+q0d3c9O2fVOr+4oqN23bNxXy2qj5JbhbiDtHUmFuv+Xo6M5xpzBZvxRIMGS3PlEwrrd0p76rbsM1Ue6fgR3Ar3vPx4FbgS9nIZRh9BJRTtH6stmtaM8ZxLU3tUvr0LDQVn/82uHTGkp4A2rMOo2ew2eiVE7sZRHE75A3dNKR7HBXp3Orv6texdUHcYqq/h1QBlDVg3Qu38Uwuooox+mBQ7Xx9a2sq5TG7+FSriPKp+CPu+YExF2jkXPUk3aCDnbIS0Mrv6texlVBHBKRAl4xRhE5BTAfhGEQ7TgtV7RmcmplXaU0JhCX1qFRUUk5kZq+0UGHdhw5ES46ozgtr+C6C09rqN+Ez0B/+sY9jaz6W/ld9TKuJqZ1wNeAE0XkZmAF8K6shDKMXiOqh0J9+05oTV2ltCaQJNPKpWeeyE3b9k0bX/7S+TXO2ooqhXyOVy+cx388/Eykuaeiyu07xmvMOz5zRCg32AtCtTpRp3GYp+2QB62vgdWruEYx3S0iO4DlVE1L71XVn2YqmWH0EC52fmgu8Sroc+jzsnqTrhd3jnmFPCJV5Xb8QIEVpxzDth/up6I6FcV074NPh+5U4pRD8LgrNx2JPoLqhFtuoures6VytRx4XQ/sOBrdscz2JDlwj2L6hqr+TFXvUtUvq+pPReQbWQtnGL1C1iYJlzyEpOvVn2OiVGb/wSN9DbY+/AwvLszho6uW8fB153Lt8GmRTmXXKb6iWuMbadbJe/xAgeHBImNrz+Gjq5bVmM2iiNrdGcnE7iBE5CigHzhWROZzxDH9YmB2q1bDCOBqknCNPKonKocgJ8KkqtO5XHIZ9h+srUuUi9ippKFUrrBu8x6GB4tNRzIdPHSYk9fcNfV5t66ptl0d3Tle0840yGxzLLeSJBPTH1LNmD6eapirryB+Dnwi7o0iciPwJuApVX1VYPxPgfcAFeAuVX2fN34V8Pve+J+p6pbUn8Yw2kjYZO9PWFHHN5p8FbXynlTlkfXnTcmyeuNYpLJwXb0Hk/xa1Td6olRmdOc4IysXRyarueAnydXfu6jw16gOeYYbsSYmVf2Yqp4M/IWqvlRVT/b+lqpqrIIAPgO8ITggImcBFwBLVXUJXrVYEXkl8DZgifeefxCR5uLhDCND0oSZ+sdfuWlXquSrYKmHqPpFxw8UnGVJs5L2lUmanghJ+ErHpeyFC8F7F6X8lNmV+dxqXMNcJ0VkwH8iIvNF5N1xb1DVbwPP1A3/MbDeL9Ohqk954xcAX1TV51X1EeAh4DWOshlG20kTZhpXgwhqJzdfKSxacxerN445+RxcZUmzkvYVz4HnDzu/Jwn/czbSkCfpnFHKr5UKbjbiqiD+QFWnehGq6n7gDxq43suB3xSR+0TkWyLya954EQimcj6B+TgMurdgWpQdPWwlm2T79ye34E4Awh3BOZFpuQxRsjRq6y/kc5x16gKuumO3c6SQC/WTeCsS5+I6183GvIVW45oHkRMRUVU/US4HzG3wesdQDZf9NWBTsIy4CyJyOXA5wMKF4eWKjZlBtxZM8+sBuTpE42z/QvVzrVh/DwcPHU60zVdU+eiqZdOS3MJ2GEK1R4LvI3HZDRQ9/0WjxfmiyOeEs05dUCPPyMrFXHfhaZFhtweePxyroIIKoD5IYMBr6rN641hk4T4jGVcF8TVgo4j8X+/5H3pjaXkCuMNTNPeLyCRwLNXaTsEi9Cd4Y9NQ1RuAGwCGhoYyamNudAOuVVLbTVqHaFzkjn+eNKv9oJIc3Tkeabry/REu568v8x3Vba5RKhXllvsfp+LlQIxPlBi5bRcbLl4a6divXyAAU4q5GOKI9/MWunVh0Yu4Koi/pKoU/th7/nXgnxq43ihwFnCviLyc6i7kp8Bm4Asi8hGqEVMvA+5v4PzGDKKbCqYFI5aiViVRDtGwyJ2oHYgLQf+CP/E1w/z+fI1yiNshNcokQF2CXLmiXHPnnshJu9Fs5m5dWPQirpnUk8AnvT8nROQW4HVUcyieANYCNwI3ish3gUPAO73dxB4R2QR8DzgMvEdVW7e/NXoS1+zkrAlbyYYR5RANm+iayQWAqpJshRmokM+x9s1LasbSVEz1KQ4UGvpsSb0d0mQz+0o8jX/IiCcpUW6Tql4iIrsJWVCo6ulR71XVSyNeenvE8R8CPhQnjzG7CFt5d8Lx6DIRJ8lVP9GtWH9P6EQ2UMhz9Avm1PgNwuzwx3sTchTz+/Ohk+/8/jz9c+dMW5G77JCiCPot/HakrcqfcMVFiVvCXHqSdhDv9f615kBG2+mWgmlJTuZG5IpSfuvOXxLbhc4/zp+Qw5SMP2GHvW/tm5dMk9N1hyTAUXWF8vyIp5HbdlGuHGk25MpAA4X0wkhS4hbR1BixCkJVn/T+faw94hhGLd1QMC3KbFIcKMRmTsfhqvySjotSHmmUq6upSqEm6sg/5zV37plSDkHEc2QcP1DgrFMXsPH+x2sK9eX7hHXnL5n2vkaIU+JhDm3DjSQT0y+I8VWp6otbLpFhdBlZmbrClF9UraawyS1JCTTbPa2eolcor/6cURFPqvDo+vOmng+ddExmu8EslLiRvIN4EYCI/DXwJPB5qjvNy4DjMpfOMLqArE1dQedqMHrIJTyzFTssF8dyowpxdOd4aoXVCHFKvNECiQaIOtgLRWSXqi5NGms3Q0NDun379k5/BANhAAAZUklEQVSKYBhN4WL/z3IVPLpzPLS3Qq5PmJzUKWXlKy6/T0QwLHbZNXdHJrQV8rlp3euyIkwRQLgZrl0ydSsiskNVh5KOc82DOCAilwFfpPo7uRQ40IR8htGTtHo16mL/rzcBJcngKmOUcjp6bo5DhycJjvqKoqI61XnOVxLrzl/CyK27QhsBtTP/IGyHsmL9PZYT0QSutZh+B7gE+In391ZvzDBmDWkruLrgYv8PhmcmyZBGxijl9Fx5MrHr2y33HSmdNjxYZMNbo40JzeYfNFOPq5uSLXsRJwWhqo+q6gWqeqyqLlDVYVV9NGPZDKOrSFPB1ZWk2Px623+SDGlkjJokXcJU648ZHixGJgo2k3/QrFKOurblRLjh2nL05SLyDS8DGhE5XUSuzlY0w+gusliNhlUh9Ts/BCu2usqQRsZmJslcSH+KLCqqNquUrcprc7iamP4RuAooA6jqA1Qb/BjGrCGL1WiwgY5fxvv6Vct4dP15bF1z9jQ7eZIMaWQcWbmYfC68EVESl5554rSxsM/SrDO4WaWchUyzCVcndb+q3i+1q4bWdRIxjB6gnfkQjcqQRsbhwWJoBFMQAQr5PkqHJ1ENj2Jq9LPUE+Zcb0U9rm5ItuxVXBXET0XkFLxgBhG5mGpehGHMGlqRDxEViukyFpzo4pLjAK65c89ULaZSucK6zXtqXvfliFMO7UwyiyrRfdEZRW7fMd7xelyzFdc8iJdS7cHwG8B+4BHgsk6X4LA8CKOXCAsrzfcJCDWlKvI5AaUmkihN7P7ozvGa2kjBa/nRRlFhqT4CXL9qGZBeITYSChxVvLC+EKAlurWGluVBiEgfMKSqrxeRo4E+Vf1FK4Q0jNlEmMM1bJIOq2tUH6kUN1lu2LI39BzlSWXDlr0ceP5wonK4bHm1W2PaxjuNNuuJ8zWYiahzJCoIVZ0UkfcBm1TVkuOMnqCZhLasSjM0G3vvT7ZRk29SPwRfhjibQbCwXSNJZo026+mW3h9GLa4+iH8Vkb8ANhLIoFbVZzKRyjCaoJmWk1m2q2y2UZBAbMinS8nuJBm2rjl7KjGtkcY7aaKOgop4XiFPPic1Ox/zNXQe1zDXVcC7gW8B2wN/htFVjO4c58pNuxqOnc8iGc4nLCY/irD/mFEr/zTd5Rb9SoH5/eE9GOb352sS06KIW9W7htnWJ8BNlMqgVRksHLV7cN1BvJKqgngt1d/pvwGfykoowwjDpQbRVXfsjswEdjHxRB0zPlHi5DV3JdY2ipOvPgKpL6Lz2kAhj0hyO06fpO5yQbY+/AwrTjmG+x/dP80xvvbNS5puvOMaZhvlj+mfO4edHzjH6bMY2eO6g/gs8Arg48D/pqowPpuVUIZRj0vJhaTJzcWeHXdMXKmHq0d3s3rjWGJJiOHBIlvXnM0j68/jw5csDc3yXXf+EiYclUO+T6byBVzZ9sP9bLh4aU3y2IaLlzI8WIxVNAOFfOKq3jUxzWok9QauO4hXqeorA8/vFZHvZSGQYYTh4vyMm1xc7dkjKxeHhojWX/fKTbtYvXGM4wcKLPqVAlsfnu6OS3LOxuU0JDmbfV541JzI7nJRVFQjI4PifBQHnnfLjXWJOjKndG/guoP4jogs95+IyJmYD8JoIy4rzqjJJSeSzp7t0FK5ojq1UwhTDmHyhRHcUQRLa7j6K/ydhr9yD6uRVE/cMXHXLU8qV2wcS11RNUjQAV4vhTmluw/XHcQZwH+IyD7v+UJgr4jsBlRVT89EOsPwcFlxRtm/0yiHDVv2Jpa6TsNAhEM4CVd/hf/5ff+HSyXWsDpKQf/JQH8+difSaGRXfYSYcqQRkfWN7k5cFcQbMpXCMBJwcX7GmWxccxtabQP/5XOHa9pupsmxCJpqwrKw/c9/9ehubt62z2Xjw4pTjplWR6n+3C7O8Uaa7oSZCX3lYH2juxMnBdHpkhqG4VoHKcz+nSa3IWqnkhNhUjVyJR+Fn73sKykXOeKUSHCVrwpXbBxzkmN+f561b14SOqG7hsjWk1aZmmO693CqxdStWC0mw4W4Oj/1K9eolbpvphrdOe48KQcRiFQuQTmSrh91TBTFgQJnnbqAex98OlKxnrzmLqfdR5LcSco7zfdgZItrLSZXJ7VhtIVm2ktGkWblmhSmOTxYZKCQ3q+gRHdqC8rhkqi3bvMeZ+UwsnIxt+8Yjw2/jXLuDxTyU13i4hzKrl3frHlP7+HqgzCMzMmqzEXakMqkMM115y8JXeU3YqaplyNJmY3uHI8t0e0jMFUFNalFqR9RFFRffj5GmO/EN2+t3jjGhi17OXjosFP9pVaUSzfaiykIo2totNBbEs00+nHxB/jjrrkLQerlSFJmriU/Llu+kOHBIqsjzGH1hf+SIop8pRmmxKOI2qGZQugdTEEYXUNWTsxGVq6jO8drmu7A9B1N2PtdfAO+wztMjrBEvXxOOOvUBbEF9OrxI5XinO6NRBSlcWhb0lvvYwrC6BqyzK5Ns3KNcwLH7WiCiihuIv/wJUvjZalzVVQqysb/fDw2uztIcSA5NyRqkk9Sxq7K2nwLM4PMnNQicqOIPCUi3w157UoRURE51nsuIvJxEXlIRB4QkVdnJZfRvbTKidmsoztplRxVunrF+numTDpRjuyBQj5WOYQl6k0S3kQojHxOpuWGXHRGcSp7OifCRWcUa5RIjXz9+ch7N7pznL6ILOyBQr6mSuwL5lj8y0wgy2/xM4Qk2InIicA5wL7A8BuBl3l/lwOfzFAuo0txLfQWh2tETRxJq+Sk0tXjEyUOHDpcbScawHf8NnPteoLz9fz+/FTRvaBst+8Yn4qgqqhy+45xzjp1QWhJjf0Hy6H3Lq5SbiGf401Lj+O58uTU2ESpnPq+G91HZiYmVf22iCwKeel64H3AvwTGLgA+p9WkjG0iMiAix6nqk1nJZ3QnzToxW+HojitY51y6uqLM78/TP3dOqoidNE2FXPIHou7HvQ8+zXUXnhYZxRQ81neMh+2q/DpXWQUYGJ2lrT4IEbkAGFfVXVK7VS0CjweeP+GNmYIwUtEKR3eY3R6qZpRg6CdUV+hRE/rEwXLq3gZh1873CQiJ3dbCIq5cej0nOb/DCuv5THqVYeOipYzepW0KQkT6gb+ial5q5jyXUzVDsXDhwhZIZswkXB3dceGrrlFPvtklTpb644ORUWEKJ+raSfJE5ZAM9OdDayu55F745ET41XlHxd7XqPsunmy2i+hN2rmDOAU4GfB3DydQLSP+GmAcCJaYPMEbm4aq3gDcANVSG1kKbPQOV4/u5pb7Ho+0kQdX2y4JeS6mrjhndtg168NXJ0plRm7dVXPduGsnObfDTDwvmNM3LWrJNffCp6KamEsysnIxqzeOTTNTqSebKYjepG2hBqq6W1X/m6ouUtVFVM1Ir1bVHwObgXd40UzLgWfN/2C4cvXobm7ati+yzlG9o7tVfafjVt4XnVGcds2wSCS/mF+zRMnybKkc6/gf3Tme2AioOFBwKkES1zPb6E0y20GIyC3A64BjReQJYK2qfjri8K8A5wIPAQeBd2UllzHzuOW+x0PHcyKhTty4vtMr1t/j7FSOW3nfvmOcoZOOcep2V/9ampLgSbIc703uUT20kxL7gruEpF1V0brEzTgy20Go6qWqepyq5lX1hHrl4O0kfuo9VlV9j6qeoqqnqaqVaDWciSqCFzUeN2EFQzxXbxxjUUwuxcjKxdNCWX3qdyRx1wy+1miYbiM5JFEmspxIQ2HGVoxv5mHZLEbPE9VCM2rctZ2nr16iJunhwSIvPCp6Ex7cGYysXEw+N12efF9tYluj5q9GckiidjWTqtNaoLrQijwWo7uwUhtGz3PpmSdy07Z9oeNhuJbECBIV0z8REiHkE9wZ+O9LimJqJkw3bQ5JFqVNrBjfzMIUhNHz+IXp/CimnAiXnnnitNaaQfyJLE2znLBJOi68s9604jJ5ZlmPqp5mqtwaswNTEMaM4Nrh02IVQhRpMpfDJumwSVY4Um47Le2ctK0/g5GEKQijYRqJtuk2orKm64mapIcHi2x/7JlUu5c4WjFpp/lezCRkxGEKwmiIrLq/tZtGM5d9oorhBUNc0yrSZibtmfK9GN2BaEQoYC8wNDSk27dbRGwnsAb0VZLuQ1iuQSGfyyy6x74XwwUR2aGqQ0nH2Q7CaIisur/1Gkn3od1VTlvxvcwE06HRGkxBGA2RRbRNL05M8wp5JkrRxfDarUib/V7MRGUEsUQ5oyFanTXbikY/wXM101EuzXUOHAqvY3TWqQuA6Ik5q/ITzX4vrapTZcwMbAdhNESz0Tb1u4UDzx9uiSkmixVw1M4mqgAfHKnF1O5cg2a/FzMdGkFMQRgN02i0TdgkHkXaianVNv84hRMnm39N3zHcTtNZM1FQ7UzUM7ofUxBG24nro1BP2omp1SvgOIWTlGTnX7OXcg0su9oIYgrCaDuuk3XUxBTnzI6atPtEGupsFlcaHIjs5ezL0mtYdrURxBSE0XbiVt7+hFuMmJiuHt3Nzdv2Tau0CtXJLSozuqLakC8iaZcQpRx6edXdSzseI1ssisloO3Hltn3lEFZqenTneI1y8AlG2fglp8NKfddH47hEO7mWBh8o5K3MtTHjsB2E0Xb8ifOKjWOhr0et2Dds2evU1nJ4sMjqiHP7x7lGO9WbXKKu/2ypzNjacyJeNYzexHYQRkcYHixSjLDRC4Su5uN8F/X2/qT8gzTx/sODRbauOZtH1p8XKXMv+hsMIwlTEEbHGFm5mLCebwqhE3XUJBzWeyEpYazRaCdrq2nMJkxBGB1jeLDoZDLyCZuco3ovJLW/bDTD2dpqGrMJ80EYHaWYIjGrlSGYzcT710f5+M5uCws1ZhqmIIyOknaidg3BTHJCt0rZWHE7YyZjCsLoKFklZrmU3GhFvH+7y3kbRjsxBWF0nCwSs9pVdM6K2xkzGXNSGzOSdpXZbnc5b8NoJ6YgjBlJu8JRLezVmMmYicmYkbSr6JwVtzNmMqIaFYne/QwNDen27ds7LUbL6MWWm4Zh9B4iskNVh5KOsx1El2DhkoZhdBuzTkE0ukrPenXfbeGStptpHruHRq+TmZNaRG4UkadE5LuBsQ0i8qCIPCAiXxKRgcBrV4nIQyKyV0RWZiGTv0of96py+qv0pKb2jb4vDd0ULtmOzzvTsXtozASyjGL6DPCGurGvA69S1dOB/wKuAhCRVwJvA5Z47/kHEUkuwp+SNBU8W/G+NHRTuGQ7Pu9Mx+6hMRPITEGo6reBZ+rG7lbVw97TbcAJ3uMLgC+q6vOq+gjwEPCaVsvU6Cq9Hav7bgqX7KbdTCO4NALKml6/h4YBnc2D+D3gq97jIvB44LUnvLGW0ugqvR2r+26qEtpNu5m0dItpp5fvoWH4dERBiMj7gcPAzQ2893IR2S4i259++ulU7210ld6u1X2wMU1Yy8120U27mbR0i2mnl++hYfi0PYpJRH4XeBPw23okCWMcODFw2Ane2DRU9QbgBqjmQaS5dqNJTbMtGaqXP2+3mHZ6+R4ahk+miXIisgj4sqq+ynv+BuAjwH9X1acDxy0BvkDV73A88A3gZapaqT9nkJmWKGc0z4r194T2lygOFNi65uwOSGQY3YdrolyWYa63AP8PWCwiT4jI7wOfAF4EfF1ExkTkUwCqugfYBHwP+BrwniTlYBhhmGnHMFqHldowZhyWoGYY8VipDWPWkkV/CcOYjZiCMFJjK3TDmB2YgjBSYUUFDWP2YA2DjFR0S56BYRjZYwrCSEW35BkYhpE9piCMVFgJCcOYPZiCMFLRDXkG3VCMzzBmA+akNlLR6RIS5iQ3jPZhCsJITSfzDLqt855hzGRMQRhA7+Q2mJPcMNqHKYg20O2Tby+ZbY4fKIQW4zMnuWG0HnNSZ0y3NLCJo5dyG7rBSW4YswVTEBnTC5NvL5ltuqnznmHMdMzElDG9MPn2mtnGivEZRnuwHUTG9EJimZltDMMIw3YQGTOycnGNAxgan3xdnN2NOMQ7ndtgGEZ3YgqiSZIm5FZNvi6RRs1EI5nZxjCMeqyjXBPUT8hQ3R1k4TR16bVs/ZgNw3Ch4z2pZwPtjFBycXb3gkPcMIzewRREE7RzQnZxdveCQ9wwjN7BFEQTtHNCdok0smgkwzBaiSmIJmjnhOySIGZJZIZhtBJzUjdJt9dZMgzDqMfVSW1hrk1i4aGGYcxUzMRkGIZhhGIKwjAMwwjFFIRhGIYRiikIwzAMIxRTEIZhGEYoPR3mKiJPA49lfJljgZ9mfI1W0Suy9oqcYLJmRa/I2ityQjpZT1LVBUkH9bSCaAcist0lXrgb6BVZe0VOMFmzoldk7RU5IRtZzcRkGIZhhGIKwjAMwwjFFEQyN3RagBT0iqy9IieYrFnRK7L2ipyQgazmgzAMwzBCsR2EYRiGEYopiAAi8qiI7BaRMRHZ7o0dIyJfF5EfeP/O7wI5F3sy+n8/F5ErRGSdiIwHxs/tkHw3ishTIvLdwFjofZQqHxeRh0TkARF5dRfIukFEHvTk+ZKIDHjji0SkFLi/n+oCWSO/cxG5yruve0VkZYfl3BiQ8VERGfPGO31PTxSRe0XkeyKyR0Te64131e81Rs5sf6uqan/eH/AocGzd2N8Ba7zHa4C/7bScdfLlgB8DJwHrgL/oApl+C3g18N2k+wicC3wVEGA5cF8XyHoOMMd7/LcBWRcFj+uS+xr6nQOvBHYBLwBOBh4Gcp2Ss+71DwMf6JJ7ehzwau/xi4D/8u5dV/1eY+TM9LdqO4hkLgA+6z3+LDDcQVnC+G3gYVXNOmHQGVX9NvBM3XDUfbwA+JxW2QYMiMhx7ZE0XFZVvVtVD3tPtwEntEueOCLuaxQXAF9U1edV9RHgIeA1mQkXIE5OERHgEuCWdsiShKo+qarf8R7/Avg+UKTLfq9Rcmb9WzUFUYsCd4vIDhG53Bt7iao+6T3+MfCSzogWyduo/c/2J95288ZuMIcFiLqPReDxwHFPeGPdwu9RXTH6nCwiO0XkWyLym50Sqo6w77xb7+tvAj9R1R8ExrrinorIImAQuI8u/r3WyRmk5b9VUxC1vFZVXw28EXiPiPxW8EWt7t26JuxLROYC5wO3ekOfBE4BlgFPUt3Kdx3ddh+jEJH3A4eBm72hJ4GFqjoI/DnwBRF5cafk8+iJ7zzApdQuaLrinorIC4HbgStU9efB17rp9xolZ1a/VVMQAVR13Pv3KeBLVLfkP/G3kN6/T3VOwmm8EfiOqv4EQFV/oqoVVZ0E/pE2mRQcibqP48CJgeNO8MY6ioj8LvAm4DJvgsAz1/zMe7yDql3/5R0TktjvvOvuq4jMAS4ENvpj3XBPRSRPddK9WVXv8Ia77vcaIWemv1VTEB4icrSIvMh/TNX5811gM/BO77B3Av/SGQlDqVmN1dlC30JV/m4h6j5uBt7hRYcsB54NbO07goi8AXgfcL6qHgyMLxCRnPf4pcDLgB92RsopmaK+883A20TkBSJyMlVZ72+3fHW8HnhQVZ/wBzp9Tz2fyKeB76vqRwIvddXvNUrOzH+r7fDA98If8FKqUR+7gD3A+73xXwG+AfwA+FfgmE7L6sl1NPAzYF5g7PPAbuABqj/k4zok2y1Ut7hlqjba34+6j1SjQf4P1RXObmCoC2R9iKqdecz7+5R37EXeb2MM+A7w5i6QNfI7B97v3de9wBs7Kac3/hngj+qO7fQ9fS1V89EDge/73G77vcbImelv1TKpDcMwjFDMxGQYhmGEYgrCMAzDCMUUhGEYhhGKKQjDMAwjFFMQhmEYRiimIAyjDq8SZjflkAAgIt8UkZ7oj2zMDExBGEYb8LKIDaOnMAVhzHpE5M9F5Lve3xXe8BwRuVlEvi8it4lIv3fseq8m/wMi8vfe2AIRuV1E/tP7W+GNrxORz4vIVuDzIrJNRJYErvtNERnysvhvFJH7veJqF3ivF0Tki54MXwIKbb0xxqzHVjXGrEZEzgDeBZxJNUv2PuBbwGKqGcBbReRG4N0i8s9Uy1mcqqoqXnMW4GPA9ar67yKyENgCvMJ77ZVUi0CWRGQ11VLXa70SGcep6nYR+RvgHlX9Pe+c94vIvwJ/CBxU1VeIyOlUM2INo23YDsKY7bwW+JKqHlDVXwJ3UC1J/biqbvWOuck77lngOeDTInIh4Ne+eT3wCal2SdsMvNirugmwWVVL3uNNwMXe40uA27zH5wBrvPd/EzgKWEi18c5NAKr6ANUyC4bRNmwHYRjh1NegUVU9LCKvodqk6WLgT4CzqS60lqvqc8E3VOurcSBwgnER+Zm3G1gF/JF/KHCRqu4Neb9hdAzbQRiznX8DhkWk36vi+xZvbKGI/Lp3zO8A/+7tCuap6leA1cBS7/W7gT/1Tygiy2Kut5Fq9c153q4AqiapP/UqdiIig974t71rIyKvAk5v6pMaRkpMQRizGq22cfwM1VLY9wH/BOynWgH1PSLyfWA+1cY8LwK+LCIPAP9OtRELwJ8BQ57j+nsc2RmEcRvVLoCbAmN/DeSBB0Rkj/cc75ov9GT4ILCjuU9rGOmwaq6GYRhGKLaDMAzDMEIxBWEYhmGEYgrCMAzDCMUUhGEYhhGKKQjDMAwjFFMQhmEYRiimIAzDMIxQTEEYhmEYofx/ZeuWCUJoYzAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model = output[\"model\"]\n",
    "model.set_weights(bestweights)\n",
    "\n",
    "pred = model.predict(tedat)\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.scatter(y.loc[te_idx], pred)\n",
    "plt.xlabel(\"observed\")\n",
    "plt.ylabel(\"predicted\")\n",
    "\n",
    "print(mse(y.loc[te_idx], pred)**.5)"
   ]
  },
  {
   "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}