Code/31C-MO_GRU-Publish-04-25-2016.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Yield Oracles:  Predicting Soybean Yields in Brazil Using Satellite and Weather Data\n",
    "\n",
    "This code was written as part of a project to satisfy the [Capstone Course](https://datascience.berkeley.edu/academics/curriculum/synthetic-capstone-course/) requirement of the [UC Berkeley Master of Information and Data Science](https://datascience.berkeley.edu) program.  Our project was completed in April 2016.  For more details and to see our results, please visit our [website](http://amitavadas.github.io/capstone/index.html).  We look forward to your feedback!\n",
    "\n",
    "For questions and comments about this code, please contact the author, Marguerite Oneto, by email at marguerite.oneto@ischool.berkeley.edu."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run A Gated Recurrent Units (GRU) Model Using Theano\n",
    "This code is based on an excellent tutorial by Denny Britz giving an introduction to Gated Recurrent Units.  Please see the links below for more information.  \n",
    "Reference:  http://www.wildml.com/2015/10/recurrent-neural-network-tutorial-part-4-implementing-a-grulstm-rnn-with-python-and-theano/   \n",
    "Reference:  https://github.com/dennybritz/rnn-tutorial-gru-lstm/blob/master/gru_theano.py   \n",
    "Data for Unit Test:  https://github.com/maoneto/W210/blob/master/Code/data/reddit-comments-2015-trunc.csv   \n",
    "Data for Yield Prediction Test:  https://github.com/maoneto/W210/blob/master/Code/data/train_trajectories_11_images_max_mv_trunc10000.csv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GRU Unit Test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define GRUTheano Class for Unit Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting gru_theano_ut.py\n"
     ]
    }
   ],
   "source": [
    "%%writefile gru_theano_ut.py\n",
    "\n",
    "import numpy as np\n",
    "import theano as theano\n",
    "import theano.tensor as T\n",
    "from theano.gradient import grad_clip\n",
    "import time\n",
    "import operator\n",
    "\n",
    "class GRUTheanoUT:\n",
    "    \n",
    "    def __init__(self, word_dim, hidden_dim=128, bptt_truncate=-1):\n",
    "        # Assign instance variables\n",
    "        self.word_dim = word_dim\n",
    "        self.hidden_dim = hidden_dim\n",
    "        self.bptt_truncate = bptt_truncate\n",
    "        # Initialize the network parameters\n",
    "        E = np.random.uniform(-np.sqrt(1./word_dim), np.sqrt(1./word_dim), (hidden_dim, word_dim))\n",
    "        U = np.random.uniform(-np.sqrt(1./hidden_dim), np.sqrt(1./hidden_dim), (6, hidden_dim, hidden_dim))\n",
    "        W = np.random.uniform(-np.sqrt(1./hidden_dim), np.sqrt(1./hidden_dim), (6, hidden_dim, hidden_dim))\n",
    "        V = np.random.uniform(-np.sqrt(1./hidden_dim), np.sqrt(1./hidden_dim), (word_dim, hidden_dim))\n",
    "        b = np.zeros((6, hidden_dim))\n",
    "        c = np.zeros(word_dim)\n",
    "        # Theano: Created shared variables\n",
    "        self.E = theano.shared(name='E', value=E.astype(theano.config.floatX))\n",
    "        self.U = theano.shared(name='U', value=U.astype(theano.config.floatX))\n",
    "        self.W = theano.shared(name='W', value=W.astype(theano.config.floatX))\n",
    "        self.V = theano.shared(name='V', value=V.astype(theano.config.floatX))\n",
    "        self.b = theano.shared(name='b', value=b.astype(theano.config.floatX))\n",
    "        self.c = theano.shared(name='c', value=c.astype(theano.config.floatX))\n",
    "        # SGD / rmsprop: Initialize parameters\n",
    "        self.mE = theano.shared(name='mE', value=np.zeros(E.shape).astype(theano.config.floatX))\n",
    "        self.mU = theano.shared(name='mU', value=np.zeros(U.shape).astype(theano.config.floatX))\n",
    "        self.mV = theano.shared(name='mV', value=np.zeros(V.shape).astype(theano.config.floatX))\n",
    "        self.mW = theano.shared(name='mW', value=np.zeros(W.shape).astype(theano.config.floatX))\n",
    "        self.mb = theano.shared(name='mb', value=np.zeros(b.shape).astype(theano.config.floatX))\n",
    "        self.mc = theano.shared(name='mc', value=np.zeros(c.shape).astype(theano.config.floatX))\n",
    "        # We store the Theano graph here\n",
    "        self.theano = {}\n",
    "        self.__theano_build__()\n",
    "    \n",
    "    def __theano_build__(self):\n",
    "        E, V, U, W, b, c = self.E, self.V, self.U, self.W, self.b, self.c\n",
    "        \n",
    "        x = T.ivector('x')\n",
    "        y = T.ivector('y')\n",
    "        \n",
    "        def forward_prop_step(x_t, s_t1_prev, s_t2_prev):\n",
    "            # This is how we calculated the hidden state in a simple RNN. No longer!\n",
    "            # s_t = T.tanh(U[:,x_t] + W.dot(s_t1_prev))\n",
    "            \n",
    "            # Word embedding layer\n",
    "            x_e = E[:,x_t]\n",
    "            \n",
    "            # GRU Layer 1\n",
    "            z_t1 = T.nnet.hard_sigmoid(U[0].dot(x_e) + W[0].dot(s_t1_prev) + b[0])\n",
    "            r_t1 = T.nnet.hard_sigmoid(U[1].dot(x_e) + W[1].dot(s_t1_prev) + b[1])\n",
    "            c_t1 = T.tanh(U[2].dot(x_e) + W[2].dot(s_t1_prev * r_t1) + b[2])\n",
    "            s_t1 = (T.ones_like(z_t1) - z_t1) * c_t1 + z_t1 * s_t1_prev\n",
    "            \n",
    "            # GRU Layer 2\n",
    "            z_t2 = T.nnet.hard_sigmoid(U[3].dot(s_t1) + W[3].dot(s_t2_prev) + b[3])\n",
    "            r_t2 = T.nnet.hard_sigmoid(U[4].dot(s_t1) + W[4].dot(s_t2_prev) + b[4])\n",
    "            c_t2 = T.tanh(U[5].dot(s_t1) + W[5].dot(s_t2_prev * r_t2) + b[5])\n",
    "            s_t2 = (T.ones_like(z_t2) - z_t2) * c_t2 + z_t2 * s_t2_prev\n",
    "            \n",
    "            # Final output calculation\n",
    "            # Theano's softmax returns a matrix with one row, we only need the row\n",
    "            o_t = T.nnet.softmax(V.dot(s_t2) + c)[0]\n",
    "\n",
    "            return [o_t, s_t1, s_t2]\n",
    "        \n",
    "        [o, s, s2], updates = theano.scan(\n",
    "            forward_prop_step,\n",
    "            sequences=x,\n",
    "            truncate_gradient=self.bptt_truncate,\n",
    "            outputs_info=[None, \n",
    "                          dict(initial=T.zeros(self.hidden_dim)),\n",
    "                          dict(initial=T.zeros(self.hidden_dim))])\n",
    "        \n",
    "        prediction = T.argmax(o, axis=1)\n",
    "        o_error = T.sum(T.nnet.categorical_crossentropy(o, y))\n",
    "        \n",
    "        # Total cost (could add regularization here)\n",
    "        cost = o_error\n",
    "        \n",
    "        # Gradients\n",
    "        dE = T.grad(cost, E)\n",
    "        dU = T.grad(cost, U)\n",
    "        dW = T.grad(cost, W)\n",
    "        db = T.grad(cost, b)\n",
    "        dV = T.grad(cost, V)\n",
    "        dc = T.grad(cost, c)\n",
    "        \n",
    "        # Assign functions\n",
    "        self.predict = theano.function([x], o)\n",
    "        self.predict_class = theano.function([x], prediction)\n",
    "        self.ce_error = theano.function([x, y], cost)\n",
    "        self.bptt = theano.function([x, y], [dE, dU, dW, db, dV, dc])\n",
    "        \n",
    "        # SGD parameters\n",
    "        learning_rate = T.scalar('learning_rate')\n",
    "        decay = T.scalar('decay')\n",
    "        \n",
    "        # rmsprop cache updates\n",
    "        mE = decay * self.mE + (1 - decay) * dE ** 2\n",
    "        mU = decay * self.mU + (1 - decay) * dU ** 2\n",
    "        mW = decay * self.mW + (1 - decay) * dW ** 2\n",
    "        mV = decay * self.mV + (1 - decay) * dV ** 2\n",
    "        mb = decay * self.mb + (1 - decay) * db ** 2\n",
    "        mc = decay * self.mc + (1 - decay) * dc ** 2\n",
    "        \n",
    "        self.sgd_step = theano.function(\n",
    "            [x, y, learning_rate, theano.Param(decay, default=0.9)],\n",
    "            [], \n",
    "            updates=[(E, E - learning_rate * dE / T.sqrt(mE + 1e-6)),\n",
    "                     (U, U - learning_rate * dU / T.sqrt(mU + 1e-6)),\n",
    "                     (W, W - learning_rate * dW / T.sqrt(mW + 1e-6)),\n",
    "                     (V, V - learning_rate * dV / T.sqrt(mV + 1e-6)),\n",
    "                     (b, b - learning_rate * db / T.sqrt(mb + 1e-6)),\n",
    "                     (c, c - learning_rate * dc / T.sqrt(mc + 1e-6)),\n",
    "                     (self.mE, mE),\n",
    "                     (self.mU, mU),\n",
    "                     (self.mW, mW),\n",
    "                     (self.mV, mV),\n",
    "                     (self.mb, mb),\n",
    "                     (self.mc, mc)\n",
    "                    ])\n",
    "        \n",
    "        \n",
    "    def calculate_total_loss(self, X, Y):\n",
    "        return np.sum([self.ce_error(x,y) for x,y in zip(X,Y)])\n",
    "    \n",
    "    def calculate_loss(self, X, Y):\n",
    "        # Divide calculate_loss by the number of words\n",
    "        num_words = np.sum([len(y) for y in Y])\n",
    "        return self.calculate_total_loss(X,Y)/float(num_words)  \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define Theano Utility Functions from utils.py for Unit Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting utils_ut.py\n"
     ]
    }
   ],
   "source": [
    "%%writefile utils_ut.py\n",
    "\n",
    "import csv\n",
    "import itertools\n",
    "import numpy as np\n",
    "import nltk\n",
    "import time\n",
    "import sys\n",
    "import operator\n",
    "import io\n",
    "import array\n",
    "from datetime import datetime\n",
    "from gru_theano_ut import GRUTheanoUT\n",
    "\n",
    "SENTENCE_START_TOKEN = \"SENTENCE_START\"\n",
    "SENTENCE_END_TOKEN = \"SENTENCE_END\"\n",
    "UNKNOWN_TOKEN = \"UNKNOWN_TOKEN\"\n",
    "\n",
    "def load_data(filename, vocabulary_size=2000, min_sent_characters=0):\n",
    "\n",
    "    word_to_index = []\n",
    "    index_to_word = []\n",
    "\n",
    "    # Read the data and append SENTENCE_START and SENTENCE_END tokens\n",
    "    print(\"Reading CSV file...\")\n",
    "    with open(filename, 'rt') as f:\n",
    "        reader = csv.reader(f, skipinitialspace=True)\n",
    "        reader.next()\n",
    "        # Split full comments into sentences\n",
    "        sentences = itertools.chain(*[nltk.sent_tokenize(x[0].decode(\"utf-8\").lower()) for x in reader])\n",
    "        # Filter sentences\n",
    "        sentences = [s for s in sentences if len(s) >= min_sent_characters]\n",
    "        sentences = [s for s in sentences if \"http\" not in s]\n",
    "        # Append SENTENCE_START and SENTENCE_END\n",
    "        sentences = [\"%s %s %s\" % (SENTENCE_START_TOKEN, x, SENTENCE_END_TOKEN) for x in sentences]\n",
    "    print(\"Parsed %d sentences.\" % (len(sentences)))\n",
    "\n",
    "    # Tokenize the sentences into words\n",
    "    tokenized_sentences = [nltk.word_tokenize(sent) for sent in sentences]\n",
    "\n",
    "    # Count the word frequencies\n",
    "    word_freq = nltk.FreqDist(itertools.chain(*tokenized_sentences))\n",
    "    print(\"Found %d unique words tokens.\" % len(word_freq.items()))\n",
    "\n",
    "    # Get the most common words and build index_to_word and word_to_index vectors\n",
    "    vocab = sorted(word_freq.items(), key=lambda x: (x[1], x[0]), reverse=True)[:vocabulary_size-2]\n",
    "    print(\"Using vocabulary size %d.\" % vocabulary_size)\n",
    "    print(\"The least frequent word in our vocabulary is '%s' and appeared %d times.\" % (vocab[-1][0], vocab[-1][1]))\n",
    "\n",
    "    sorted_vocab = sorted(vocab, key=operator.itemgetter(1))\n",
    "    index_to_word = [\"<MASK/>\", UNKNOWN_TOKEN] + [x[0] for x in sorted_vocab]\n",
    "    word_to_index = dict([(w, i) for i, w in enumerate(index_to_word)])\n",
    "\n",
    "    # Replace all words not in our vocabulary with the unknown token\n",
    "    for i, sent in enumerate(tokenized_sentences):\n",
    "        tokenized_sentences[i] = [w if w in word_to_index else UNKNOWN_TOKEN for w in sent]\n",
    "\n",
    "    # Create the training data\n",
    "    X_train = np.asarray([[word_to_index[w] for w in sent[:-1]] for sent in tokenized_sentences])\n",
    "    y_train = np.asarray([[word_to_index[w] for w in sent[1:]] for sent in tokenized_sentences])\n",
    "\n",
    "    return X_train, y_train, word_to_index, index_to_word\n",
    "\n",
    "\n",
    "def train_with_sgd(model, X_train, y_train, learning_rate=0.001, nepoch=20, decay=0.9,\n",
    "    callback_every=10000, callback=None):\n",
    "    num_examples_seen = 0\n",
    "    for epoch in range(nepoch):\n",
    "        # For each training example...\n",
    "        for i in np.random.permutation(len(y_train)):\n",
    "            # One SGD step\n",
    "            model.sgd_step(X_train[i], y_train[i], learning_rate, decay)\n",
    "            num_examples_seen += 1\n",
    "            # Optionally do callback\n",
    "            if (callback and callback_every and num_examples_seen % callback_every == 0):\n",
    "                callback(model, num_examples_seen)            \n",
    "    return model\n",
    "\n",
    "def save_model_parameters_theano(model, outfile):\n",
    "    np.savez(outfile,\n",
    "        E=model.E.get_value(),\n",
    "        U=model.U.get_value(),\n",
    "        W=model.W.get_value(),\n",
    "        V=model.V.get_value(),\n",
    "        b=model.b.get_value(),\n",
    "        c=model.c.get_value())\n",
    "    print \"Saved model parameters to %s.\" % outfile\n",
    "\n",
    "def load_model_parameters_theano(path, modelClass=GRUTheanoUT):\n",
    "    npzfile = np.load(path)\n",
    "    E, U, W, V, b, c = npzfile[\"E\"], npzfile[\"U\"], npzfile[\"W\"], npzfile[\"V\"], npzfile[\"b\"], npzfile[\"c\"]\n",
    "    hidden_dim, word_dim = E.shape[0], E.shape[1]\n",
    "    print \"Building model model from %s with hidden_dim=%d word_dim=%d\" % (path, hidden_dim, word_dim)\n",
    "    sys.stdout.flush()\n",
    "    model = modelClass(word_dim, hidden_dim=hidden_dim)\n",
    "    model.E.set_value(E)\n",
    "    model.U.set_value(U)\n",
    "    model.W.set_value(W)\n",
    "    model.V.set_value(V)\n",
    "    model.b.set_value(b)\n",
    "    model.c.set_value(c)\n",
    "    return model \n",
    "\n",
    "def gradient_check_theano(model, x, y, h=0.001, error_threshold=0.01):\n",
    "    # Overwrite the bptt attribute. We need to backpropagate all the way to get the correct gradient\n",
    "    model.bptt_truncate = 1000\n",
    "    # Calculate the gradients using backprop\n",
    "    bptt_gradients = model.bptt(x, y)\n",
    "    # List of all parameters we want to chec.\n",
    "    model_parameters = ['E', 'U', 'W', 'b', 'V', 'c']\n",
    "    # Gradient check for each parameter\n",
    "    for pidx, pname in enumerate(model_parameters):\n",
    "        # Get the actual parameter value from the mode, e.g. model.W\n",
    "        parameter_T = operator.attrgetter(pname)(model)\n",
    "        parameter = parameter_T.get_value()\n",
    "        print \"Performing gradient check for parameter %s with size %d.\" % (pname, np.prod(parameter.shape))\n",
    "        # Iterate over each element of the parameter matrix, e.g. (0,0), (0,1), ...\n",
    "        it = np.nditer(parameter, flags=['multi_index'], op_flags=['readwrite'])\n",
    "        while not it.finished:\n",
    "            ix = it.multi_index\n",
    "            # Save the original value so we can reset it later\n",
    "            original_value = parameter[ix]\n",
    "            # Estimate the gradient using (f(x+h) - f(x-h))/(2*h)\n",
    "            parameter[ix] = original_value + h\n",
    "            parameter_T.set_value(parameter)\n",
    "            gradplus = model.calculate_total_loss([x],[y])\n",
    "            parameter[ix] = original_value - h\n",
    "            parameter_T.set_value(parameter)\n",
    "            gradminus = model.calculate_total_loss([x],[y])\n",
    "            estimated_gradient = (gradplus - gradminus)/(2*h)\n",
    "            parameter[ix] = original_value\n",
    "            parameter_T.set_value(parameter)\n",
    "            # The gradient for this parameter calculated using backpropagation\n",
    "            backprop_gradient = bptt_gradients[pidx][ix]\n",
    "            # calculate The relative error: (|x - y|/(|x| + |y|))\n",
    "            relative_error = np.abs(backprop_gradient - estimated_gradient)/(np.abs(backprop_gradient) + np.abs(estimated_gradient))\n",
    "            # If the error is to large fail the gradient check\n",
    "            if relative_error > error_threshold:\n",
    "                print \"Gradient Check ERROR: parameter=%s ix=%s\" % (pname, ix)\n",
    "                print \"+h Loss: %f\" % gradplus\n",
    "                print \"-h Loss: %f\" % gradminus\n",
    "                print \"Estimated_gradient: %f\" % estimated_gradient\n",
    "                print \"Backpropagation gradient: %f\" % backprop_gradient\n",
    "                print \"Relative Error: %f\" % relative_error\n",
    "                return \n",
    "            it.iternext()\n",
    "        print \"Gradient check for parameter %s passed.\" % (pname)\n",
    "\n",
    "\n",
    "def print_sentence(s, index_to_word):\n",
    "    sentence_str = [index_to_word[x] for x in s[1:-1]]\n",
    "    print(\" \".join(sentence_str))\n",
    "    sys.stdout.flush()\n",
    "\n",
    "def generate_sentence(model, index_to_word, word_to_index, min_length=5):\n",
    "    # We start the sentence with the start token\n",
    "    new_sentence = [word_to_index[SENTENCE_START_TOKEN]]\n",
    "    # Repeat until we get an end token\n",
    "    while not new_sentence[-1] == word_to_index[SENTENCE_END_TOKEN]:\n",
    "        next_word_probs = model.predict(new_sentence)[-1]\n",
    "        samples = np.random.multinomial(1, next_word_probs)\n",
    "        sampled_word = np.argmax(samples)\n",
    "        new_sentence.append(sampled_word)\n",
    "        # Seomtimes we get stuck if the sentence becomes too long, e.g. \"........\" :(\n",
    "        # And: We don't want sentences with UNKNOWN_TOKEN's\n",
    "        if len(new_sentence) > 100 or sampled_word == word_to_index[UNKNOWN_TOKEN]:\n",
    "            return None\n",
    "    if len(new_sentence) < min_length:\n",
    "        return None\n",
    "    return new_sentence\n",
    "\n",
    "def generate_sentences(model, n, index_to_word, word_to_index):\n",
    "    for i in range(n):\n",
    "        sent = None\n",
    "        while not sent:\n",
    "            sent = generate_sentence(model, index_to_word, word_to_index)\n",
    "        print_sentence(sent, index_to_word)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train GRU for Unit Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reading CSV file...\n",
      "Parsed 94 sentences.\n",
      "Found 632 unique words tokens.\n",
      "Using vocabulary size 200.\n",
      "The least frequent word in our vocabulary is 'we' and appeared 1 times.\n",
      "SGD Step time: 58.140039 milliseconds\n",
      "\n",
      "2016-04-25T18:34:47.266062 (25)\n",
      "--------------------------------------------------\n",
      "Loss: 4.135619\n",
      "myself SENTENCE_START not while\n",
      "much in fact now that there similar % second\n",
      "in too all 're anyone however connections think someone tone overlap to\n",
      "when SENTENCE_START cost que more connections in the with with\n",
      "que é from who .\n",
      "s : from his right wrong n't anyone can did n't same\n",
      "no lot be its so in\n",
      ". it back % character . work .\n",
      "then of each much welcome down see ? by , leia season\n",
      "n't like last $ . 's a $ the\n",
      "Saved model parameters to GRU-2016-04-25-18-33-200-48-128.dat.\n",
      "\n",
      "\n",
      "\n",
      "2016-04-25T18:34:48.804555 (50)\n",
      "--------------------------------------------------\n",
      "Loss: 3.976999\n",
      "for have n't é )\n",
      "from yeah leia a how where with get\n",
      "from : lidstrom multiply they\n",
      "right with award\n",
      ". so though v.\n",
      "still is yt good on\n",
      "writing but up on stream\n",
      "was season tone 2 guys on be fact . result but\n",
      "while riven that the similar with have or\n",
      "leia guys other , with end someone\n",
      "Saved model parameters to GRU-2016-04-25-18-33-200-48-128.dat.\n",
      "\n",
      "\n",
      "\n",
      "2016-04-25T18:34:50.175742 (75)\n",
      "--------------------------------------------------\n",
      "Loss: 3.921327\n",
      "more what just\n",
      "the similar happens they\n",
      "result to get\n",
      "`` have but define but\n",
      "seeing connections yellow\n",
      "paul s in the\n",
      "would for : .\n",
      "other happens think it the he if\n",
      "SENTENCE_START good leia\n",
      "hit have said ?\n",
      "Saved model parameters to GRU-2016-04-25-18-33-200-48-128.dat.\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#! /usr/bin/env python\n",
    "\n",
    "import sys\n",
    "import os\n",
    "import time\n",
    "import numpy as np\n",
    "from utils_ut import *\n",
    "from datetime import datetime\n",
    "from gru_theano_ut import GRUTheanoUT\n",
    "\n",
    "LEARNING_RATE = float(os.environ.get(\"LEARNING_RATE\", \"0.001\"))\n",
    "VOCABULARY_SIZE = int(os.environ.get(\"VOCABULARY_SIZE\", \"200\"))\n",
    "EMBEDDING_DIM = int(os.environ.get(\"EMBEDDING_DIM\", \"48\"))\n",
    "HIDDEN_DIM = int(os.environ.get(\"HIDDEN_DIM\", \"128\"))\n",
    "NEPOCH = int(os.environ.get(\"NEPOCH\", \"1\"))\n",
    "MODEL_OUTPUT_FILE = os.environ.get(\"MODEL_OUTPUT_FILE\")\n",
    "INPUT_DATA_FILE = os.environ.get(\"INPUT_DATA_FILE\", \"./data/reddit-comments-2015-trunc.csv\")\n",
    "PRINT_EVERY = int(os.environ.get(\"PRINT_EVERY\", \"25\"))\n",
    "\n",
    "if not MODEL_OUTPUT_FILE:\n",
    "  ts = datetime.now().strftime(\"%Y-%m-%d-%H-%M\")\n",
    "  MODEL_OUTPUT_FILE = \"GRU-%s-%s-%s-%s.dat\" % (ts, VOCABULARY_SIZE, EMBEDDING_DIM, HIDDEN_DIM)\n",
    "\n",
    "# Load data\n",
    "x_train, y_train, word_to_index, index_to_word = load_data(INPUT_DATA_FILE, VOCABULARY_SIZE)\n",
    "\n",
    "# Build model\n",
    "model = GRUTheanoUT(VOCABULARY_SIZE, hidden_dim=HIDDEN_DIM, bptt_truncate=-1)\n",
    "\n",
    "# Print SGD step time\n",
    "t1 = time.time()\n",
    "model.sgd_step(x_train[10], y_train[10], LEARNING_RATE)\n",
    "t2 = time.time()\n",
    "print \"SGD Step time: %f milliseconds\" % ((t2 - t1) * 1000.)\n",
    "sys.stdout.flush()\n",
    "\n",
    "# We do this every few examples to understand what's going on\n",
    "def sgd_callback(model, num_examples_seen):\n",
    "  dt = datetime.now().isoformat()\n",
    "  loss = model.calculate_loss(x_train[:10000], y_train[:10000])\n",
    "  print(\"\\n%s (%d)\" % (dt, num_examples_seen))\n",
    "  print(\"--------------------------------------------------\")\n",
    "  print(\"Loss: %f\" % loss)\n",
    "  generate_sentences(model, 10, index_to_word, word_to_index)\n",
    "  save_model_parameters_theano(model, MODEL_OUTPUT_FILE)\n",
    "  print(\"\\n\")\n",
    "  sys.stdout.flush()\n",
    "\n",
    "for epoch in range(NEPOCH):\n",
    "  train_with_sgd(model, x_train, y_train, learning_rate=LEARNING_RATE, nepoch=1, decay=0.9, \n",
    "    callback_every=PRINT_EVERY, callback=sgd_callback)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GRU for Yield Prediction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create Training, Validation, and Test Datasets for Yield Prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X_train and X_test Shape:\n",
      "(8000, 7) (2000, 7)\n",
      "Y_train, Y_test, log_Y_train, log_Y_test, log_Y_bar Shape:\n",
      "(8000, 7) (2000, 7) (8000, 7) (2000, 7)\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "import csv\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from sklearn.cross_validation import train_test_split\n",
    "\n",
    "X = []\n",
    "Y = []\n",
    "\n",
    "with open('train_trajectories_11_images_max_mv_trunc10000.csv', 'r') as csvfile:\n",
    "    datareader = csv.reader(csvfile, delimiter=',')\n",
    "    for row in datareader:\n",
    "        label = row.pop()  # pop the last element in the list which is the label (yield)\n",
    "        if float(label) != 0.0:\n",
    "            X.append(row)\n",
    "            Y.append(len(X[0])*[label]) # output at each t, o_t, is the yield\n",
    "X = np.array(X).astype(np.float)\n",
    "Y = np.array(Y).astype(np.float)\n",
    "\n",
    "# Break labeled examples into train and test sets\n",
    "X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.20, random_state=22)\n",
    "\n",
    "log_Y_train = np.log(Y_train)\n",
    "log_Y_test = np.log(Y_test)\n",
    "\n",
    "print \"X_train and X_test Shape:\"\n",
    "print X_train.shape, X_test.shape\n",
    "\n",
    "print \"Y_train, Y_test, log_Y_train, log_Y_test, log_Y_bar Shape:\"\n",
    "print Y_train.shape, Y_test.shape, log_Y_train.shape, log_Y_test.shape\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define GRUTheano Class for Yield Prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting gru_theano_yp.py\n"
     ]
    }
   ],
   "source": [
    "%%writefile gru_theano_yp.py\n",
    "#!/usr/bin/env python\n",
    "\n",
    "import numpy as np\n",
    "import theano as theano\n",
    "import theano.tensor as T\n",
    "from theano.gradient import grad_clip\n",
    "import time\n",
    "import operator\n",
    "\n",
    "theano.exception_verbosity='high'\n",
    "theano.mode='FAST_COMPILE'\n",
    "theano.allow_gc=False\n",
    "theano.optimizer='fast_compile'\n",
    "theano.config.compute_test_value = 'off'\n",
    "\n",
    "class GRUTheanoYP:\n",
    "    \n",
    "    def __init__(self, x_dim, hidden_dim=128, bptt_truncate=-1):\n",
    "        # Assign instance variables\n",
    "        self.x_dim = x_dim\n",
    "        self.hidden_dim = hidden_dim\n",
    "        self.bptt_truncate = bptt_truncate\n",
    "        # Initialize the network parameters\n",
    "        E = np.random.uniform(-np.sqrt(1./x_dim), np.sqrt(1./x_dim), (hidden_dim, x_dim))\n",
    "        U = np.random.uniform(-np.sqrt(1./hidden_dim), np.sqrt(1./hidden_dim), (6, hidden_dim, hidden_dim))\n",
    "        W = np.random.uniform(-np.sqrt(1./hidden_dim), np.sqrt(1./hidden_dim), (6, hidden_dim, hidden_dim))\n",
    "        V = np.random.uniform(-np.sqrt(1./hidden_dim), np.sqrt(1./hidden_dim), (x_dim, hidden_dim))\n",
    "        b = np.zeros((6, hidden_dim, 1))\n",
    "        c = np.zeros((x_dim, 1))\n",
    "        # Theano: Created shared variables\n",
    "        self.E = theano.shared(name='E', value=E.astype(theano.config.floatX))\n",
    "        self.U = theano.shared(name='U', value=U.astype(theano.config.floatX))\n",
    "        self.W = theano.shared(name='W', value=W.astype(theano.config.floatX))\n",
    "        self.V = theano.shared(name='V', value=V.astype(theano.config.floatX))\n",
    "        self.b = theano.shared(name='b', value=b.astype(theano.config.floatX))\n",
    "        self.c = theano.shared(name='c', value=c.astype(theano.config.floatX))\n",
    "        # SGD / rmsprop: Initialize parameters\n",
    "        self.mE = theano.shared(name='mE', value=np.zeros(E.shape).astype(theano.config.floatX))\n",
    "        self.mU = theano.shared(name='mU', value=np.zeros(U.shape).astype(theano.config.floatX))\n",
    "        self.mV = theano.shared(name='mV', value=np.zeros(V.shape).astype(theano.config.floatX))\n",
    "        self.mW = theano.shared(name='mW', value=np.zeros(W.shape).astype(theano.config.floatX))\n",
    "        self.mb = theano.shared(name='mb', value=np.zeros(b.shape).astype(theano.config.floatX))\n",
    "        self.mc = theano.shared(name='mc', value=np.zeros(c.shape).astype(theano.config.floatX))\n",
    "        # We store the Theano graph here\n",
    "        self.theano = {}\n",
    "        self.__theano_build__()\n",
    "    \n",
    "    def __theano_build__(self):\n",
    "        E, V, U, W, b, c = self.E, self.V, self.U, self.W, self.b, self.c\n",
    "        \n",
    "        x = T.vector('x')\n",
    "        y = T.vector('y')\n",
    "        \n",
    "        def forward_prop_step(x_t, s_t1_prev, s_t2_prev):\n",
    "            # This is how we calculated the hidden state in a simple RNN. No longer!\n",
    "            # s_t = T.tanh(U[:,x_t] + W.dot(s_t1_prev))\n",
    "            \n",
    "            # Embedding layer: x_t is the NDVI value of the pixel\n",
    "            x_e = T.mul(E, x_t)\n",
    "            \n",
    "            # GRU Layer 1\n",
    "            z_t1 = T.nnet.hard_sigmoid(T.dot(U[0], x_e) + T.dot(W[0], s_t1_prev) + b[0])\n",
    "            r_t1 = T.nnet.hard_sigmoid(T.dot(U[1], x_e) + T.dot(W[1], s_t1_prev) + b[1])\n",
    "            c_t1 = T.tanh(T.dot(U[2], x_e) + W[2].dot(s_t1_prev * r_t1) + b[2])\n",
    "            s_t1 = (T.ones_like(z_t1) - z_t1) * c_t1 + z_t1 * s_t1_prev\n",
    "            # For debugging ...\n",
    "            # print T.shape(s_t1).eval({x_t: 0.8888, s_t1_prev: np.asarray([[0], [0], [0], [0], [0]])})     \n",
    "            # GRU Layer 2\n",
    "            z_t2 = T.nnet.hard_sigmoid(U[3].dot(s_t1) + W[3].dot(s_t2_prev) + b[3])\n",
    "            r_t2 = T.nnet.hard_sigmoid(U[4].dot(s_t1) + W[4].dot(s_t2_prev) + b[4])\n",
    "            c_t2 = T.tanh(U[5].dot(s_t1) + W[5].dot(s_t2_prev * r_t2) + b[5])\n",
    "            s_t2 = (T.ones_like(z_t2) - z_t2) * c_t2 + z_t2 * s_t2_prev\n",
    "            \n",
    "            # Final output calculation\n",
    "            # Theano's softmax returns a matrix with one row, we only need the row\n",
    "            # o_t = T.nnet.softmax(V.dot(s_t2) + c)[0]\n",
    "            # We do not use softmax because we are predicting a continuous variable, not doing classification:\n",
    "            o_t = (V.dot(s_t2) + c)[0]\n",
    "\n",
    "            return [o_t, s_t1, s_t2]\n",
    "        \n",
    "        [o, s, s2], updates = theano.scan(\n",
    "            forward_prop_step,\n",
    "            sequences=[x],\n",
    "            truncate_gradient=self.bptt_truncate,\n",
    "            outputs_info=[None, \n",
    "                          dict(initial=T.zeros_like(E)),\n",
    "                          dict(initial=T.zeros_like(E))])\n",
    "        \n",
    "        # prediction = T.argmax(o, axis=1)\n",
    "        # o_error = T.sum(T.nnet.categorical_crossentropy(o, y))\n",
    "        # Again, we are predicting a continuous variable, not doing classification.\n",
    "        # We use the sum of squared errors (SSE) as our cost function:\n",
    "        prediction = o\n",
    "        o_error = T.sum(T.sqr(o - T.reshape(y, [T.shape(y)[0],1,1], ndim=3)))\n",
    "        \n",
    "        # Total cost (could add regularization here)\n",
    "        cost = o_error\n",
    "        \n",
    "        # Gradients\n",
    "        dE = T.grad(cost, E)\n",
    "        dU = T.grad(cost, U)\n",
    "        dW = T.grad(cost, W)\n",
    "        db = T.grad(cost, b)\n",
    "        dV = T.grad(cost, V)\n",
    "        dc = T.grad(cost, c)\n",
    "        \n",
    "        # Assign functions\n",
    "        self.predict = theano.function([x], o)\n",
    "        # We do not need this because we are not doing classification\n",
    "        # self.predict_class = theano.function([x], prediction)\n",
    "        self.sse_error = theano.function([x, y], cost)\n",
    "        self.bptt = theano.function([x, y], [dU, dW, db, dV, dc])\n",
    "        \n",
    "        # SGD parameters\n",
    "        learning_rate = T.scalar('learning_rate')\n",
    "        decay = T.scalar('decay')\n",
    "        \n",
    "        # rmsprop cache updates\n",
    "        mE = decay * self.mE + (1 - decay) * dE ** 2\n",
    "        mU = decay * self.mU + (1 - decay) * dU ** 2\n",
    "        mW = decay * self.mW + (1 - decay) * dW ** 2\n",
    "        mV = decay * self.mV + (1 - decay) * dV ** 2\n",
    "        mb = decay * self.mb + (1 - decay) * db ** 2\n",
    "        mc = decay * self.mc + (1 - decay) * dc ** 2\n",
    "        \n",
    "        self.sgd_step = theano.function(\n",
    "            [x, y, learning_rate, theano.Param(decay, default=0.9)],\n",
    "            [], \n",
    "            updates=[(U, U - learning_rate * dU / T.sqrt(mU + 1e-6)),\n",
    "                     (W, W - learning_rate * dW / T.sqrt(mW + 1e-6)),\n",
    "                     (V, V - learning_rate * dV / T.sqrt(mV + 1e-6)),\n",
    "                     (b, b - learning_rate * db / T.sqrt(mb + 1e-6)),\n",
    "                     (c, c - learning_rate * dc / T.sqrt(mc + 1e-6)),\n",
    "                     (self.mU, mU),\n",
    "                     (self.mW, mW),\n",
    "                     (self.mV, mV),\n",
    "                     (self.mb, mb),\n",
    "                     (self.mc, mc)\n",
    "                    ])\n",
    "        \n",
    "        \n",
    "    def calculate_total_loss(self, X, Y):\n",
    "        return np.sum([self.sse_error(x,y) for x,y in zip(X,Y)])\n",
    "    \n",
    "    def calculate_loss(self, X, Y):\n",
    "        # Divide calculate_loss by the number of examples\n",
    "        num_examples = np.sum([len(y) for y in Y])\n",
    "        return self.calculate_total_loss(X,Y)/float(num_examples)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define Theano Utility Functions from utils.py for Yield Prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Overwriting utils_yp.py\n"
     ]
    }
   ],
   "source": [
    "%%writefile utils_yp.py\n",
    "\n",
    "import csv\n",
    "import itertools\n",
    "import numpy as np\n",
    "import nltk\n",
    "import time\n",
    "import sys\n",
    "import operator\n",
    "import io\n",
    "import array\n",
    "from datetime import datetime\n",
    "from gru_theano_yp import GRUTheanoYP\n",
    "\n",
    "def shuffle_data(p, X, y):\n",
    "    # shuffle it\n",
    "    shuffle = np.random.permutation(np.arange(X.shape[0]))\n",
    "    X, y = X[shuffle], y[shuffle]\n",
    "    # divide \n",
    "    n_train = np.round(X.shape[0]*p)\n",
    "    return X[:n_train], y[:n_train], X[n_train:], y[n_train:]\n",
    "\n",
    "# Used to shrink the size of the training and test datasets for debugging\n",
    "def shuffle_data2(X, y, n_train, n_test):\n",
    "    # shuffle it\n",
    "    shuffle = np.random.permutation(np.arange(X.shape[0]))\n",
    "    X, y = X[shuffle], y[shuffle]\n",
    "    # divide \n",
    "    return X[:n_train], y[:n_train], X[n_train:n_train + n_test], y[n_train:n_train + n_test]\n",
    "\n",
    "def save_model_parameters_theano(model, outfile):\n",
    "    np.savez(outfile,\n",
    "        E=model.E.get_value(),\n",
    "        U=model.U.get_value(),\n",
    "        W=model.W.get_value(),\n",
    "        V=model.V.get_value(),\n",
    "        b=model.b.get_value(),\n",
    "        c=model.c.get_value())\n",
    "#     print \"Saved model parameters to %s.\" % outfile\n",
    "\n",
    "def load_model_parameters_theano(path, modelClass=GRUTheanoYP):\n",
    "    npzfile = np.load(path)\n",
    "    E, U, W, V, b, c = npzfile[\"E\"], npzfile[\"U\"], npzfile[\"W\"], npzfile[\"V\"], npzfile[\"b\"], npzfile[\"c\"]\n",
    "    hidden_dim, x_dim = E.shape[0], E.shape[1]\n",
    "    print \"Building model from %s with hidden_dim=%d x_dim=%d\" % (path, hidden_dim, x_dim)\n",
    "    sys.stdout.flush()\n",
    "    model = modelClass(x_dim, hidden_dim=hidden_dim)\n",
    "    model.E.set_value(E)\n",
    "    model.U.set_value(U)\n",
    "    model.W.set_value(W)\n",
    "    model.V.set_value(V)\n",
    "    model.b.set_value(b)\n",
    "    model.c.set_value(c)\n",
    "    return model \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Train GRU for Yield Prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Start Time: 2016-04-25 18:40:06\n",
      "SGD Step time: 25.578022 milliseconds\n",
      "Epochs: 100, TrainingSamples: 8000, ModelTime: 20160425184114\n",
      "1) ExamplesSeen=7200, TrainLoss=0.9803, ValLoss=0.9588, ValR2=-0.6773, ValRMSE=0.3701, TrainTime=2.59 min, EndTime=06:43:57 (saved)\n",
      "2) ExamplesSeen=14400, TrainLoss=0.6720, ValLoss=0.6408, ValR2=-0.0658, ValRMSE=0.3026, TrainTime=2.99 min, EndTime=06:47:04 (saved)\n",
      "3) ExamplesSeen=21600, TrainLoss=0.6137, ValLoss=0.5375, ValR2=-0.0164, ValRMSE=0.2771, TrainTime=3.05 min, EndTime=06:50:14 (saved)\n",
      "4) ExamplesSeen=28800, TrainLoss=0.6190, ValLoss=0.6038, ValR2=-0.0367, ValRMSE=0.2937, TrainTime=2.97 min, EndTime=06:53:19 \n",
      "5) ExamplesSeen=36000, TrainLoss=0.6061, ValLoss=0.6145, ValR2=-0.0263, ValRMSE=0.2963, TrainTime=3.05 min, EndTime=06:56:28 \n",
      "6) ExamplesSeen=43200, TrainLoss=0.6232, ValLoss=0.6342, ValR2=-0.0574, ValRMSE=0.3010, TrainTime=2.99 min, EndTime=06:59:33 \n",
      "7) ExamplesSeen=50400, TrainLoss=0.6122, ValLoss=0.5663, ValR2=-0.0163, ValRMSE=0.2844, TrainTime=3.16 min, EndTime=07:02:49 \n",
      "8) ExamplesSeen=57600, TrainLoss=1.3032, ValLoss=1.4270, ValR2=-1.0802, ValRMSE=0.4515, TrainTime=3.20 min, EndTime=07:06:07 \n",
      "9) ExamplesSeen=64800, TrainLoss=0.7038, ValLoss=0.7538, ValR2=-0.2333, ValRMSE=0.3282, TrainTime=3.17 min, EndTime=07:09:23 \n",
      "10) ExamplesSeen=72000, TrainLoss=2.8250, ValLoss=2.8096, ValR2=-3.7463, ValRMSE=0.6335, TrainTime=3.17 min, EndTime=07:12:39 \n",
      "11) ExamplesSeen=79200, TrainLoss=0.6055, ValLoss=0.6124, ValR2=-0.0283, ValRMSE=0.2958, TrainTime=3.16 min, EndTime=07:15:55 \n",
      "12) ExamplesSeen=86400, TrainLoss=0.6043, ValLoss=0.6098, ValR2=-0.0069, ValRMSE=0.2951, TrainTime=3.16 min, EndTime=07:19:11 \n",
      "13) ExamplesSeen=93600, TrainLoss=0.6268, ValLoss=0.6937, ValR2=-0.0511, ValRMSE=0.3148, TrainTime=3.17 min, EndTime=07:22:27 \n",
      "14) ExamplesSeen=100800, TrainLoss=0.6992, ValLoss=0.6789, ValR2=-0.2196, ValRMSE=0.3114, TrainTime=3.16 min, EndTime=07:25:43 \n",
      "15) ExamplesSeen=108000, TrainLoss=2.2451, ValLoss=2.2353, ValR2=-2.9077, ValRMSE=0.5651, TrainTime=3.16 min, EndTime=07:28:59 \n",
      "16) ExamplesSeen=115200, TrainLoss=1.1983, ValLoss=1.2674, ValR2=-1.0014, ValRMSE=0.4255, TrainTime=3.29 min, EndTime=07:32:23 \n",
      "17) ExamplesSeen=122400, TrainLoss=0.6704, ValLoss=0.6520, ValR2=-0.1228, ValRMSE=0.3052, TrainTime=3.40 min, EndTime=07:35:54 \n",
      "18) ExamplesSeen=129600, TrainLoss=1.0103, ValLoss=0.9660, ValR2=-0.6862, ValRMSE=0.3715, TrainTime=3.44 min, EndTime=07:39:27 \n",
      "19) ExamplesSeen=136800, TrainLoss=1.3847, ValLoss=1.4458, ValR2=-1.4236, ValRMSE=0.4545, TrainTime=3.43 min, EndTime=07:43:00 \n",
      "20) ExamplesSeen=144000, TrainLoss=0.7284, ValLoss=0.6904, ValR2=-0.1877, ValRMSE=0.3140, TrainTime=3.47 min, EndTime=07:46:34 \n",
      "21) ExamplesSeen=151200, TrainLoss=0.6332, ValLoss=0.6367, ValR2=-0.0754, ValRMSE=0.3016, TrainTime=3.47 min, EndTime=07:50:09 \n",
      "22) ExamplesSeen=158400, TrainLoss=0.6326, ValLoss=0.6213, ValR2=-0.0619, ValRMSE=0.2979, TrainTime=3.44 min, EndTime=07:53:43 \n",
      "23) ExamplesSeen=165600, TrainLoss=0.6200, ValLoss=0.6160, ValR2=-0.0414, ValRMSE=0.2966, TrainTime=3.38 min, EndTime=07:57:12 \n",
      "24) ExamplesSeen=172800, TrainLoss=0.7124, ValLoss=0.7279, ValR2=-0.2801, ValRMSE=0.3225, TrainTime=3.35 min, EndTime=08:00:40 \n",
      "25) ExamplesSeen=180000, TrainLoss=0.6362, ValLoss=0.6619, ValR2=-0.0590, ValRMSE=0.3075, TrainTime=3.32 min, EndTime=08:04:05 \n",
      "26) ExamplesSeen=187200, TrainLoss=0.7171, ValLoss=0.6968, ValR2=-0.1671, ValRMSE=0.3155, TrainTime=3.37 min, EndTime=08:07:34 \n",
      "27) ExamplesSeen=194400, TrainLoss=0.8331, ValLoss=0.8048, ValR2=-0.3041, ValRMSE=0.3391, TrainTime=3.39 min, EndTime=08:11:04 \n",
      "28) ExamplesSeen=201600, TrainLoss=0.6768, ValLoss=0.6649, ValR2=-0.1226, ValRMSE=0.3082, TrainTime=3.40 min, EndTime=08:14:35 \n",
      "29) ExamplesSeen=208800, TrainLoss=0.7954, ValLoss=0.8067, ValR2=-0.3481, ValRMSE=0.3395, TrainTime=3.27 min, EndTime=08:17:57 \n",
      "30) ExamplesSeen=216000, TrainLoss=0.6135, ValLoss=0.6470, ValR2=-0.0390, ValRMSE=0.3040, TrainTime=3.25 min, EndTime=08:21:18 \n",
      "31) ExamplesSeen=223200, TrainLoss=0.6334, ValLoss=0.5841, ValR2=-0.0458, ValRMSE=0.2889, TrainTime=3.28 min, EndTime=08:24:41 \n",
      "32) ExamplesSeen=230400, TrainLoss=0.6509, ValLoss=0.5885, ValR2=-0.0696, ValRMSE=0.2899, TrainTime=3.28 min, EndTime=08:28:04 \n",
      "33) ExamplesSeen=237600, TrainLoss=1.0427, ValLoss=1.0250, ValR2=-0.7392, ValRMSE=0.3827, TrainTime=3.23 min, EndTime=08:31:23 \n",
      "34) ExamplesSeen=244800, TrainLoss=0.6964, ValLoss=0.7240, ValR2=-0.1709, ValRMSE=0.3216, TrainTime=3.18 min, EndTime=08:34:40 \n",
      "35) ExamplesSeen=252000, TrainLoss=0.6814, ValLoss=0.6839, ValR2=-0.1253, ValRMSE=0.3126, TrainTime=3.18 min, EndTime=08:37:57 \n",
      "36) ExamplesSeen=259200, TrainLoss=0.8927, ValLoss=0.8359, ValR2=-0.4416, ValRMSE=0.3456, TrainTime=3.17 min, EndTime=08:41:14 \n",
      "37) ExamplesSeen=266400, TrainLoss=0.6187, ValLoss=0.6161, ValR2=-0.0509, ValRMSE=0.2967, TrainTime=3.16 min, EndTime=08:44:29 \n",
      "38) ExamplesSeen=273600, TrainLoss=0.6687, ValLoss=0.6684, ValR2=-0.0980, ValRMSE=0.3090, TrainTime=3.19 min, EndTime=08:47:46 \n",
      "39) ExamplesSeen=280800, TrainLoss=1.2828, ValLoss=1.2485, ValR2=-1.1082, ValRMSE=0.4223, TrainTime=3.26 min, EndTime=08:51:08 \n",
      "40) ExamplesSeen=288000, TrainLoss=0.5997, ValLoss=0.6388, ValR2=-0.0104, ValRMSE=0.3021, TrainTime=3.25 min, EndTime=08:54:29 \n",
      "41) ExamplesSeen=295200, TrainLoss=0.6340, ValLoss=0.6536, ValR2=-0.0714, ValRMSE=0.3056, TrainTime=3.26 min, EndTime=08:57:50 \n",
      "42) ExamplesSeen=302400, TrainLoss=0.6020, ValLoss=0.6100, ValR2=-0.0068, ValRMSE=0.2952, TrainTime=3.24 min, EndTime=09:01:11 \n",
      "43) ExamplesSeen=309600, TrainLoss=0.6860, ValLoss=0.6976, ValR2=-0.2293, ValRMSE=0.3157, TrainTime=3.27 min, EndTime=09:04:33 \n",
      "44) ExamplesSeen=316800, TrainLoss=1.0441, ValLoss=1.0488, ValR2=-0.7068, ValRMSE=0.3871, TrainTime=3.21 min, EndTime=09:07:52 \n",
      "45) ExamplesSeen=324000, TrainLoss=0.7382, ValLoss=0.7020, ValR2=-0.2234, ValRMSE=0.3167, TrainTime=3.21 min, EndTime=09:11:11 \n",
      "46) ExamplesSeen=331200, TrainLoss=0.7713, ValLoss=0.7931, ValR2=-0.3101, ValRMSE=0.3366, TrainTime=3.31 min, EndTime=09:14:35 \n",
      "47) ExamplesSeen=338400, TrainLoss=0.6280, ValLoss=0.6366, ValR2=-0.0459, ValRMSE=0.3016, TrainTime=3.26 min, EndTime=09:17:57 \n",
      "48) ExamplesSeen=345600, TrainLoss=0.6144, ValLoss=0.6033, ValR2=-0.0279, ValRMSE=0.2936, TrainTime=3.26 min, EndTime=09:21:18 \n",
      "49) ExamplesSeen=352800, TrainLoss=0.6724, ValLoss=0.6476, ValR2=-0.1247, ValRMSE=0.3042, TrainTime=3.27 min, EndTime=09:24:40 \n",
      "50) ExamplesSeen=360000, TrainLoss=0.6977, ValLoss=0.6236, ValR2=-0.1735, ValRMSE=0.2985, TrainTime=3.27 min, EndTime=09:28:03 \n",
      "51) ExamplesSeen=367200, TrainLoss=0.9264, ValLoss=0.9168, ValR2=-0.5670, ValRMSE=0.3619, TrainTime=3.26 min, EndTime=09:31:24 \n",
      "52) ExamplesSeen=374400, TrainLoss=0.7697, ValLoss=0.7399, ValR2=-0.2403, ValRMSE=0.3251, TrainTime=3.27 min, EndTime=09:34:46 \n",
      "53) ExamplesSeen=381600, TrainLoss=0.6423, ValLoss=0.6454, ValR2=-0.0717, ValRMSE=0.3037, TrainTime=3.24 min, EndTime=09:38:06 \n",
      "54) ExamplesSeen=388800, TrainLoss=0.6948, ValLoss=0.6729, ValR2=-0.1693, ValRMSE=0.3100, TrainTime=3.24 min, EndTime=09:41:27 \n",
      "55) ExamplesSeen=396000, TrainLoss=0.6071, ValLoss=0.5962, ValR2=-0.0138, ValRMSE=0.2919, TrainTime=3.21 min, EndTime=09:44:45 \n",
      "56) ExamplesSeen=403200, TrainLoss=1.0773, ValLoss=1.1789, ValR2=-0.8432, ValRMSE=0.4104, TrainTime=3.24 min, EndTime=09:48:06 \n",
      "57) ExamplesSeen=410400, TrainLoss=0.6226, ValLoss=0.6389, ValR2=-0.0411, ValRMSE=0.3021, TrainTime=3.24 min, EndTime=09:51:26 \n",
      "58) ExamplesSeen=417600, TrainLoss=0.6177, ValLoss=0.6190, ValR2=-0.0320, ValRMSE=0.2974, TrainTime=3.23 min, EndTime=09:54:46 \n",
      "59) ExamplesSeen=424800, TrainLoss=0.6126, ValLoss=0.5951, ValR2=-0.0216, ValRMSE=0.2916, TrainTime=3.27 min, EndTime=09:58:08 \n",
      "60) ExamplesSeen=432000, TrainLoss=0.8821, ValLoss=0.8633, ValR2=-0.4670, ValRMSE=0.3512, TrainTime=3.28 min, EndTime=10:01:31 \n",
      "61) ExamplesSeen=439200, TrainLoss=0.7365, ValLoss=0.7636, ValR2=-0.2666, ValRMSE=0.3303, TrainTime=3.29 min, EndTime=10:04:54 \n",
      "62) ExamplesSeen=446400, TrainLoss=0.6013, ValLoss=0.6393, ValR2=-0.0102, ValRMSE=0.3022, TrainTime=3.29 min, EndTime=10:08:18 \n",
      "63) ExamplesSeen=453600, TrainLoss=0.6398, ValLoss=0.6732, ValR2=-0.0812, ValRMSE=0.3101, TrainTime=3.27 min, EndTime=10:11:40 \n",
      "64) ExamplesSeen=460800, TrainLoss=0.6538, ValLoss=0.6511, ValR2=-0.0647, ValRMSE=0.3050, TrainTime=3.29 min, EndTime=10:15:03 \n",
      "65) ExamplesSeen=468000, TrainLoss=0.6681, ValLoss=0.6425, ValR2=-0.1089, ValRMSE=0.3030, TrainTime=3.30 min, EndTime=10:18:27 \n",
      "66) ExamplesSeen=475200, TrainLoss=0.7962, ValLoss=0.8180, ValR2=-0.3240, ValRMSE=0.3418, TrainTime=3.13 min, EndTime=10:21:40 \n",
      "67) ExamplesSeen=482400, TrainLoss=1.0521, ValLoss=1.1469, ValR2=-0.6872, ValRMSE=0.4048, TrainTime=3.13 min, EndTime=10:24:54 \n",
      "68) ExamplesSeen=489600, TrainLoss=0.6255, ValLoss=0.6121, ValR2=-0.0292, ValRMSE=0.2957, TrainTime=3.24 min, EndTime=10:28:15 \n",
      "69) ExamplesSeen=496800, TrainLoss=0.6116, ValLoss=0.6102, ValR2=-0.0250, ValRMSE=0.2953, TrainTime=3.28 min, EndTime=10:31:38 \n",
      "70) ExamplesSeen=504000, TrainLoss=0.6812, ValLoss=0.6912, ValR2=-0.1381, ValRMSE=0.3142, TrainTime=3.28 min, EndTime=10:35:00 \n",
      "71) ExamplesSeen=511200, TrainLoss=0.6169, ValLoss=0.6332, ValR2=-0.0259, ValRMSE=0.3008, TrainTime=3.23 min, EndTime=10:38:20 \n",
      "72) ExamplesSeen=518400, TrainLoss=0.6727, ValLoss=0.6522, ValR2=-0.1497, ValRMSE=0.3052, TrainTime=3.28 min, EndTime=10:41:43 \n",
      "73) ExamplesSeen=525600, TrainLoss=0.6210, ValLoss=0.5830, ValR2=-0.0539, ValRMSE=0.2886, TrainTime=3.26 min, EndTime=10:45:04 \n",
      "74) ExamplesSeen=532800, TrainLoss=0.7657, ValLoss=0.7392, ValR2=-0.1842, ValRMSE=0.3250, TrainTime=3.27 min, EndTime=10:48:27 \n",
      "75) ExamplesSeen=540000, TrainLoss=0.6036, ValLoss=0.6139, ValR2=-0.0132, ValRMSE=0.2961, TrainTime=3.25 min, EndTime=10:51:48 \n",
      "76) ExamplesSeen=547200, TrainLoss=0.6061, ValLoss=0.5781, ValR2=-0.0124, ValRMSE=0.2874, TrainTime=3.24 min, EndTime=10:55:08 \n",
      "77) ExamplesSeen=554400, TrainLoss=0.6695, ValLoss=0.6484, ValR2=-0.0948, ValRMSE=0.3043, TrainTime=3.23 min, EndTime=10:58:28 \n",
      "78) ExamplesSeen=561600, TrainLoss=0.7782, ValLoss=0.7547, ValR2=-0.2326, ValRMSE=0.3283, TrainTime=3.31 min, EndTime=11:01:53 \n",
      "79) ExamplesSeen=568800, TrainLoss=0.6107, ValLoss=0.5860, ValR2=-0.0168, ValRMSE=0.2893, TrainTime=3.26 min, EndTime=11:05:15 \n",
      "80) ExamplesSeen=576000, TrainLoss=0.7498, ValLoss=0.7782, ValR2=-0.2770, ValRMSE=0.3334, TrainTime=3.27 min, EndTime=11:08:37 \n",
      "81) ExamplesSeen=583200, TrainLoss=0.6121, ValLoss=0.6348, ValR2=-0.0235, ValRMSE=0.3011, TrainTime=3.29 min, EndTime=11:12:00 \n",
      "82) ExamplesSeen=590400, TrainLoss=0.8285, ValLoss=0.8519, ValR2=-0.4161, ValRMSE=0.3488, TrainTime=3.27 min, EndTime=11:15:23 \n",
      "83) ExamplesSeen=597600, TrainLoss=0.6011, ValLoss=0.6011, ValR2=-0.0111, ValRMSE=0.2930, TrainTime=3.27 min, EndTime=11:18:45 \n",
      "84) ExamplesSeen=604800, TrainLoss=0.6761, ValLoss=0.6469, ValR2=-0.1276, ValRMSE=0.3040, TrainTime=3.28 min, EndTime=11:22:08 \n",
      "85) ExamplesSeen=612000, TrainLoss=0.6434, ValLoss=0.6398, ValR2=-0.0678, ValRMSE=0.3023, TrainTime=3.27 min, EndTime=11:25:30 \n",
      "86) ExamplesSeen=619200, TrainLoss=0.6087, ValLoss=0.5794, ValR2=-0.0126, ValRMSE=0.2877, TrainTime=3.29 min, EndTime=11:28:54 \n",
      "87) ExamplesSeen=626400, TrainLoss=0.6064, ValLoss=0.6047, ValR2=-0.0084, ValRMSE=0.2939, TrainTime=3.26 min, EndTime=11:32:16 \n",
      "88) ExamplesSeen=633600, TrainLoss=0.8032, ValLoss=0.7420, ValR2=-0.2693, ValRMSE=0.3256, TrainTime=3.23 min, EndTime=11:35:36 \n",
      "89) ExamplesSeen=640800, TrainLoss=0.6477, ValLoss=0.6773, ValR2=-0.1558, ValRMSE=0.3110, TrainTime=3.27 min, EndTime=11:38:58 \n",
      "90) ExamplesSeen=648000, TrainLoss=0.6927, ValLoss=0.7621, ValR2=-0.2232, ValRMSE=0.3300, TrainTime=3.26 min, EndTime=11:42:19 \n",
      "91) ExamplesSeen=655200, TrainLoss=0.6946, ValLoss=0.6593, ValR2=-0.1533, ValRMSE=0.3069, TrainTime=3.20 min, EndTime=11:45:37 \n",
      "92) ExamplesSeen=662400, TrainLoss=0.6913, ValLoss=0.6071, ValR2=-0.1556, ValRMSE=0.2945, TrainTime=3.27 min, EndTime=11:49:00 \n",
      "93) ExamplesSeen=669600, TrainLoss=0.6099, ValLoss=0.6568, ValR2=-0.0314, ValRMSE=0.3063, TrainTime=3.30 min, EndTime=11:52:23 \n",
      "94) ExamplesSeen=676800, TrainLoss=0.6166, ValLoss=0.6084, ValR2=-0.0296, ValRMSE=0.2948, TrainTime=3.26 min, EndTime=11:55:45 \n",
      "95) ExamplesSeen=684000, TrainLoss=0.6232, ValLoss=0.6107, ValR2=-0.0323, ValRMSE=0.2954, TrainTime=3.26 min, EndTime=11:59:07 \n",
      "96) ExamplesSeen=691200, TrainLoss=0.6043, ValLoss=0.5982, ValR2=-0.0079, ValRMSE=0.2923, TrainTime=3.29 min, EndTime=12:02:30 \n",
      "97) ExamplesSeen=698400, TrainLoss=0.6155, ValLoss=0.5462, ValR2=-0.0102, ValRMSE=0.2793, TrainTime=3.28 min, EndTime=12:05:53 \n",
      "98) ExamplesSeen=705600, TrainLoss=0.6149, ValLoss=0.6126, ValR2=-0.0152, ValRMSE=0.2958, TrainTime=3.28 min, EndTime=12:09:16 \n",
      "99) ExamplesSeen=712800, TrainLoss=0.7542, ValLoss=0.7707, ValR2=-0.2946, ValRMSE=0.3318, TrainTime=3.29 min, EndTime=12:12:40 \n",
      "100) ExamplesSeen=720000, TrainLoss=0.6116, ValLoss=0.5889, ValR2=-0.0234, ValRMSE=0.2900, TrainTime=3.33 min, EndTime=12:16:06 \n",
      "Finished.  Total train time = 5.58 hours\n",
      "End Time: 2016-04-26 00:16:06\n"
     ]
    }
   ],
   "source": [
    "#! /usr/bin/env python\n",
    "# Make sure the latest code updates are loaded\n",
    "%reload_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import sys\n",
    "import os\n",
    "import time\n",
    "import numpy as np\n",
    "from utils_yp import *\n",
    "from datetime import datetime\n",
    "from gru_theano_yp import GRUTheanoYP\n",
    "\n",
    "LEARNING_RATE = float(os.environ.get(\"LEARNING_RATE\", \"0.01\"))\n",
    "DECAY_RATE = float(os.environ.get(\"DECAY_RATE\", \"0.9\"))\n",
    "X_DIM = int(os.environ.get(\"X_DIM\", \"1\"))  # number of features in the x_t vector\n",
    "HIDDEN_DIM = int(os.environ.get(\"HIDDEN_DIM\", \"128\"))\n",
    "NEPOCH = int(os.environ.get(\"NEPOCH\", \"100\"))\n",
    "\n",
    "def train_with_sgd(model, X_train, y_train, X_test, y_test, nepoch=20, learning_rate=0.001, decay=0.9):\n",
    "    # Set start of training time\n",
    "    start_time = time.time()\n",
    "    # Set model time\n",
    "    modeltime = datetime.now().strftime(\"%Y%m%d%H%M%S\")\n",
    "    print 'Epochs: %d, TrainingSamples: %s, ModelTime: %s' %(nepoch, y_train.shape[0], modeltime)      \n",
    "    # Set epoch variables\n",
    "    min_loss = 100000\n",
    "    train_losses = []\n",
    "    val_losses = []\n",
    "    num_examples_seen = 0\n",
    "    for epoch in range(nepoch):\n",
    "        epoch_start = time.time()\n",
    "        # Divide the data\n",
    "        train_data, train_labels, val_data, val_labels = shuffle_data(0.9, X_train, y_train) # Use this for full sample training  \n",
    "        # train_data, train_labels, val_data, val_labels = shuffle_data2(X_train, y_train, 10000, 1000) # Use this for testing and debugging  \n",
    "        # For each training example...\n",
    "        for i in range(len(train_data)):\n",
    "            # One SGD step\n",
    "            model.sgd_step(train_data[i], train_labels[i], learning_rate, decay)\n",
    "            num_examples_seen += 1\n",
    "        epoch_time = time.time() - epoch_start\n",
    "        # Calculate training loss\n",
    "        train_loss = model.calculate_loss(train_data, train_labels)\n",
    "        train_losses.append((num_examples_seen, train_loss))\n",
    "        # Calculate validation loss, R-squared, and RMSE\n",
    "        val_loss = model.calculate_loss(val_data, val_labels)\n",
    "        val_losses.append((num_examples_seen, val_loss))\n",
    "        val_predictions = []\n",
    "        for i in range(len(val_data)):\n",
    "            val_predictions.append(model.predict(val_data[i]))\n",
    "        val_predictions = np.reshape(val_predictions, (len(val_predictions), len(val_predictions[0])))\n",
    "        val_R_squared = 1 - np.sum(np.square(val_predictions - val_labels))/np.sum(np.square(val_labels - np.mean(val_labels)))\n",
    "        val_rmse = np.sqrt(np.mean(np.square(val_predictions - val_labels)))\n",
    "        # If validation loss is a new minimum, save predictions and model   \n",
    "        is_saved = ''\n",
    "        if val_loss < min_loss:\n",
    "            min_loss = val_loss\n",
    "            # Make and save predictions\n",
    "            X_test_predictions = []\n",
    "            for i in range(len(X_test)):\n",
    "                X_test_predictions.append(np.append(model.predict(X_test[i]).reshape(len(X_test[i])), float(y_test[i][0])))\n",
    "            predictions_and_labels = np.asarray(X_test_predictions)\n",
    "            filename = \"./predictions/GRUs/pred-%s.txt\" %(modeltime)\n",
    "            np.savetxt(filename, predictions_and_labels, fmt='%.18f', delimiter=',',)\n",
    "            # Save model parameters\n",
    "            filename = \"./models/GRUs/GRU-%s.npz\" % (modeltime)\n",
    "            save_model_parameters_theano(model, filename)\n",
    "            is_saved = '(saved)'\n",
    "        # Print epoch stats\n",
    "        print '%d) ExamplesSeen=%d, TrainLoss=%.4f, ValLoss=%.4f, ValR2=%.4f, ValRMSE=%.4f, TrainTime=%.2f min, EndTime=%s %s' %(epoch+1, num_examples_seen, train_loss, val_loss, val_R_squared, val_rmse, epoch_time/60, time.strftime(\"%I:%M:%S\"), is_saved)\n",
    "    print 'Finished.  Total train time = %.2f hours' %((time.time() - start_time)/3600)   \n",
    "    \n",
    "    return train_losses, val_losses\n",
    "\n",
    "# Start Time\n",
    "print 'Start Time: %s' %(datetime.now().strftime(\"%Y-%m-%d %H:%M:%S\"))\n",
    "\n",
    "# Build model\n",
    "model = GRUTheanoYP(x_dim=X_DIM, hidden_dim=HIDDEN_DIM, bptt_truncate=-1)\n",
    "\n",
    "# Measure and Print SGD step time\n",
    "t1 = time.time()\n",
    "model.sgd_step(X_train[10], Y_train[10], LEARNING_RATE, DECAY_RATE)\n",
    "t2 = time.time()\n",
    "print \"SGD Step time: %f milliseconds\" % ((t2 - t1) * 1000.)\n",
    "\n",
    "# Train the model\n",
    "train_losses, val_losses = train_with_sgd(model, X_train, log_Y_train, X_test, log_Y_test, nepoch=NEPOCH, learning_rate=LEARNING_RATE, decay=DECAY_RATE)\n",
    "\n",
    "#End Time\n",
    "print 'End Time: %s' %(datetime.now().strftime(\"%Y-%m-%d %H:%M:%S\")) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Baseline Model:  The Mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R-Squared: -0.0002, MAE: 0.7594, RMSE: 0.9684, MAPE: 26.20%\n"
     ]
    }
   ],
   "source": [
    "Y_train_1d = Y_train[:,0]\n",
    "y_hat = len(Y_test)*[np.mean(Y_train_1d)]\n",
    "y_test = Y_test[:,0]\n",
    "\n",
    "y_hat = np.array(y_hat).astype(np.float)\n",
    "y_test = np.array(y_test).astype(np.float)\n",
    "y_bar = np.mean(y_test)\n",
    "\n",
    "R_squared = 1 - np.sum(np.square(y_hat - y_test))/np.sum(np.square(y_test - y_bar))\n",
    "mae = np.mean(np.abs(y_hat - y_test))\n",
    "rmse = np.sqrt(np.mean(np.square(y_hat - y_test)))\n",
    "mape = np.mean(abs(np.divide((y_hat - y_test), y_test)))\n",
    "\n",
    "print 'R-Squared: %.4f, MAE: %.4f, RMSE: %.4f, MAPE: %.2f%s' %(R_squared, mae, rmse, mape*100, '%')  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate Performance Stats on Holdout Sample\n",
    "These are the predictions on a holdout sample (X_test above) using the model with the lowest validation loss.    \n",
    "- t-0 = prediction at harvest time   \n",
    "- t-6 = prediction three months out from harvest   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Log-Transformed Yield Results:\n",
      "t-6) R-Squared: 0.0000, MAE: 0.2453, RMSE: 0.2978, MAPE: 26.09%\n",
      "t-5) R-Squared: -0.0001, MAE: 0.2454, RMSE: 0.2978, MAPE: 26.00%\n",
      "t-4) R-Squared: 0.0000, MAE: 0.2454, RMSE: 0.2978, MAPE: 26.06%\n",
      "t-3) R-Squared: -0.0002, MAE: 0.2454, RMSE: 0.2978, MAPE: 26.00%\n",
      "t-2) R-Squared: -0.0000, MAE: 0.2454, RMSE: 0.2978, MAPE: 26.06%\n",
      "t-1) R-Squared: -0.0002, MAE: 0.2454, RMSE: 0.2978, MAPE: 26.00%\n",
      "t-0) R-Squared: -0.0000, MAE: 0.2454, RMSE: 0.2978, MAPE: 26.06%\n",
      "\n",
      "Values Transformed Back to Original Units Results:\n",
      "t-6) R-Squared: -0.0196, MAE: 0.7597, RMSE: 0.9778, MAPE: 25.18%\n",
      "t-5) R-Squared: -0.0235, MAE: 0.7600, RMSE: 0.9797, MAPE: 25.08%\n",
      "t-4) R-Squared: -0.0209, MAE: 0.7597, RMSE: 0.9784, MAPE: 25.14%\n",
      "t-3) R-Squared: -0.0235, MAE: 0.7600, RMSE: 0.9797, MAPE: 25.08%\n",
      "t-2) R-Squared: -0.0209, MAE: 0.7597, RMSE: 0.9785, MAPE: 25.14%\n",
      "t-1) R-Squared: -0.0235, MAE: 0.7600, RMSE: 0.9797, MAPE: 25.08%\n",
      "t-0) R-Squared: -0.0209, MAE: 0.7597, RMSE: 0.9785, MAPE: 25.14%\n"
     ]
    }
   ],
   "source": [
    "y_hat = []\n",
    "y_test = []\n",
    "with open('./predictions/GRUs/pred-20160425154620.txt', 'r') as csvfile:\n",
    "    datareader = csv.reader(csvfile, delimiter=',')\n",
    "    for row in datareader:\n",
    "        label = row.pop()\n",
    "        y_hat.append(row)\n",
    "        y_test.append(label)\n",
    "\n",
    "# Results keeping predictions as logs of yield\n",
    "y_hat = np.array(y_hat).astype(np.float)\n",
    "y_test = np.array(y_test).astype(np.float)\n",
    "y_bar = np.mean(y_test)\n",
    "\n",
    "R_squared_list_tr = []\n",
    "mae_list_tr = []\n",
    "rmse_list_tr = []\n",
    "mape_list_tr = []\n",
    "\n",
    "print 'Log-Transformed Yield Results:'\n",
    "for j in range(len(y_hat[0])):\n",
    "    y_hat_last = y_hat[:, j]\n",
    "    R_squared = 1 - np.sum(np.square(y_hat_last - y_test))/np.sum(np.square(y_test - y_bar))\n",
    "    mae = np.mean(np.abs(y_hat_last - y_test))\n",
    "    rmse = np.sqrt(np.mean(np.square(y_hat_last - y_test)))\n",
    "    count = 0\n",
    "    sums = 0\n",
    "    for i in range(len(y_test)):\n",
    "        if y_test[i] != 0: \n",
    "            error = np.divide(abs(y_hat_last[i] - y_test[i]), abs(y_test[i]))\n",
    "            count += 1\n",
    "            sums += error\n",
    "    mean_abs_pct_error = sums/count\n",
    "    R_squared_list_tr.append(R_squared)\n",
    "    mae_list_tr.append(mae)\n",
    "    rmse_list_tr.append(rmse)\n",
    "    mape_list_tr.append(mean_abs_pct_error)\n",
    "    print 't-%d) R-Squared: %.4f, MAE: %.4f, RMSE: %.4f, MAPE: %.2f%s' %(6-j, R_squared, mae, rmse, mean_abs_pct_error*100, '%')  \n",
    "    \n",
    "# Results transforming predictions and yields back to original units\n",
    "y_hat = np.array(y_hat).astype(np.float)\n",
    "y_test = np.array(y_test).astype(np.float)\n",
    "\n",
    "y_hat_exp = np.exp(y_hat)\n",
    "y_test_exp = np.exp(y_test)\n",
    "y_bar_exp = np.mean(y_test_exp)\n",
    "\n",
    "R_squared_list_exp = []\n",
    "mae_list_exp = []\n",
    "rmse_list_exp = []\n",
    "mape_list_exp = []\n",
    "\n",
    "print '\\nValues Transformed Back to Original Units Results:'\n",
    "for j in range(len(y_hat_exp[0])):\n",
    "    y_hat_last = y_hat_exp[:, j]\n",
    "    R_squared = 1 - np.sum(np.square(y_hat_last - y_test_exp))/np.sum(np.square(y_test_exp - y_bar_exp))\n",
    "    mae = np.mean(np.abs(y_hat_last - y_test_exp))\n",
    "    rmse = np.sqrt(np.mean(np.square(y_hat_last - y_test_exp)))\n",
    "    count = 0\n",
    "    sums = 0\n",
    "    for i in range(len(y_test_exp)):\n",
    "        if y_test_exp[i] != 0: \n",
    "            error = np.divide(abs(y_hat_last[i] - y_test_exp[i]), abs(y_test_exp[i]))\n",
    "            count += 1\n",
    "            sums += error\n",
    "    mean_abs_pct_error = sums/count\n",
    "    R_squared_list_exp.append(R_squared)\n",
    "    mae_list_exp.append(mae)\n",
    "    rmse_list_exp.append(rmse)\n",
    "    mape_list_exp.append(mean_abs_pct_error)\n",
    "    print 't-%d) R-Squared: %.4f, MAE: %.4f, RMSE: %.4f, MAPE: %.2f%s' %(len(y_hat_exp[0])-1-j, R_squared, mae, rmse, mean_abs_pct_error*100, '%')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Graph Losses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAscAAAF/CAYAAAClh2SrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmUXGd95//3c2uvrt43dau1WJa8Y7CxWWzAwhnAgEPI\nOXMcAxlC4jkk+cUzY2Z+hJmEATuBw4HYOVl+wEBChhgDCSeZELANJGDEQOIlGNy25UWWJVlSd0ut\n7uqtuqrurar7/P6oVkvdXV1d3aqlS/68ztFRLbfvfXS7Wv2pb33v8xhrLSIiIiIiAk6jByAiIiIi\nslkoHIuIiIiILFA4FhERERFZoHAsIiIiIrJA4VhEREREZIHCsYiIiIjIgrLh2BgTNcY8aox5whjz\njDHmU6ts92fGmBeMMcPGmKtqM1QRERERkdoKlnvSWps1xrzZWps2xgSBnxhj3mCt/cnpbYwx7wB2\nW2v3GGNeC3weeF1thy0iIiIiUn1rtlVYa9MLN8NAAEgu2+RdwF8vbPso0GGM6a/mIEVERERE6mHN\ncGyMcYwxTwAngR9aa59ZtslW4NhZ948DQ9UbooiIiIhIfVRSOfatta+iGHjfZIzZW2Izs/zLqjA2\nEREREZG6KttzfDZr7Ywx5gHgGmDfWU+NANvOuj+08NgSxhgFZhERERGpOWvt8sJtxcqGY2NMD5C3\n1k4bY2LAW4C7lm32LeB24G+MMa8Dpq21J1cZ6EbHKU3uzjvv5M4772z0MGST0utDVqPXhqxGrw1Z\njTEbzsXA2pXjAeCvjTEOxRaMr1hrf2CM+U0Aa+0XrLUPGmPeYYw5CMwDv35OIxIRERERaZC1pnJ7\nCri6xONfWHb/9iqPS0RERESk7rRCntTF3r17Gz0E2cT0+pDV6LUhq9FrQ2rF1KsP2Bhj1XMsIiIi\nIrVkjKndBXkiIiIim8W5Xmgl559aFF4VjkVERKRp6FNoOa1Wb5bUcywiIiIiskDhWERERERkgcKx\niIiIiMgChWMRERGRTeId73gHX/nKVxo9jJc1hWMRERGRc5BIJGhtbaW1tRXHcYjH44v3v/71r69r\nXw8++CD/4T/8hw2NY+fOnfzgBz/Y0NfKGZqtQkREROQcpFKpxdsXXHABX/rSl7jxxhtXbJfP5wkG\naxe9jDGa7q4KVDkWERERqYF9+/YxNDTEZz7zGQYGBrjtttuYnp7m5ptvpq+vj66uLn7xF3+RkZGR\nxa/Zu3cvX/rSlwD48pe/zBve8AY+/OEP09XVxa5du/jud7+77nG4rssdd9zB1q1b2bp1Kx/60Ifw\nPA+AiYkJbr75Zjo7O+nu7uZNb3rT4td9+tOfZmhoiLa2Ni655BIeeuihczwjzUHhWERERKRGTp48\nydTUFEePHuULX/gCvu9z2223cfToUY4ePUosFuP2229f3H559fexxx7jkksuYXJykt/93d/ltttu\nW/cYPvnJT/LYY48xPDzM8PAwjz32GJ/4xCcAuOeee9i2bRsTExOMj4/zqU99CoDnn3+ez372s/z0\npz9ldnaWf/qnf2Lnzp3ndjKahMKxiIiISI04jsNdd91FKBQiGo3S1dXFL//yLxONRkkkEvze7/0e\nP/rRj1b9+h07dnDbbbdhjOH9738/Y2NjjI+Pr2sMX/va1/jYxz5GT08PPT09fPzjH1+86C8cDjM2\nNsaRI0cIBAJcf/31AAQCAVzXZf/+/eRyObZv386uXbs2fiKaiHqORURE5Lywz+yryn722r1V2Q9A\nb28v4XB48X46neZDH/oQ3/ve95iamgKKPcvW2pL9wlu2bFm8HY/HF7fv6+ureAyjo6Ps2LFj8f72\n7dsZHR0F4MMf/jB33nknb33rWwH44Ac/yEc+8hF2797Nn/zJn3DnnXeyf/9+3va2t/HHf/zHDAwM\nrONf35wUjkVEROS8UM1QWy3LA+8999zDgQMHeOyxx+jr6+OJJ57g6quvXjUcV8Pg4CBHjhzh0ksv\nBeDo0aMMDg4CxZk27r77bu6++27279/PjTfeyLXXXsuNN97Ie97zHt7znvcwNzfHb/7mb/KRj3yE\ne++9tyZj3EzUVlEjT//7p8kcyjR6GCIiIrKJpFIpYrEY7e3tJJNJ7rrrrqru3/M8stns4p98Ps97\n3vMePvGJTzAxMcHExAR/8Ad/sDhd3P3338/Bgwex1tLW1kYgECAQCHDgwAEeeughXNclEokQjUYJ\nBAJVHetmpXBcI4efzOOOeY0ehoiIiDTQ8mrwHXfcQSaToaenh+uuu463v/3tq1aMS03NtlZ1+R3v\neAfxeHzxzx/8wR/w0Y9+lGuuuYYrr7ySK6+8kmuuuYaPfvSjABw8eJC3vOUttLa2ct111/E7v/M7\n3HDDDbiuy//4H/+D3t5eBgYGmJiYWLxY73xnrLX1OZAxtl7H2gwujab4X3/mc8MH2xo9FBERkfOC\nMYaXU5aQ8lZ7PSw8vuEeFVWOayRTcEjN6gdYREREpJkoHNeI5zu46UaPQkRERETWQ+G4RjxryGZU\nORYRERFpJgrHNeJZh2xa4VhERESkmSgc14D1LTkcXFWORURERJqKwnENFFwfDwc3U5vJvEVERESk\nNhSOa8BL+fgYsllVjkVERESaicJxDWTniqHYzTZ4ICIiIiKyLgrHNZCe9QFwXVWORUREpDzHcTh0\n6BAAv/3bv80nPvGJirZdr69+9au87W1v29DXvpwoHNfAYuXYVc+xiIjI+e6mm27i4x//+IrH//Ef\n/5GBgQF83694X5///OcXl3Y+F0eOHMFxnCXHft/73sf3vve9c973cvv27WPbtm1V32+jKBzXQHqu\n+ELMug0eiIiIiNTcBz7wAe67774Vj3/lK1/hV3/1V3GcxsUtLbe9fgrHNZBNna4cN3ggIiIiUnO/\n9Eu/xOTkJD/+8Y8XH5uamuKBBx7g/e9/P4899hivf/3r6ezsZHBwkP/0n/4TuVyu5L4+8IEP8D//\n5/9cvP9Hf/RHDA4OMjQ0xF/91V8t2faBBx7gqquuor29ne3bt3PXXXctPvemN70JgI6ODtra2njk\nkUf48pe/zBvf+MbFbf71X/+Va6+9lo6ODl7zmtfw8MMPLz63d+9ePvaxj/GGN7yBtrY23va2tzE5\nObnuc/Pss8+yd+9eOjs7ueKKK/j2t7+9+NyDDz7I5ZdfTltbG0NDQ9xzzz0ATExMcPPNN9PZ2Ul3\ndzdvetOb6hryFY5rILMQjj2vwQMRERGRmovFYtxyyy3ce++9i4994xvf4NJLL+UVr3gFwWCQP/3T\nP2VycpKHH36YH/zgB3zuc58ruS9jDMYU2zK/+93vcs899/D973+fAwcO8P3vf3/JtolEgvvuu4+Z\nmRkeeOABPv/5z/OP//iPAItBfWZmhtnZWV73utct+dpkMsk73/lO7rjjDpLJJP/1v/5X3vnOdzI1\nNbW4zde//nW+/OUvMz4+jud53H333es6L7lcjl/8xV/kpptu4tSpU/z5n/8573vf+3jhhRcAuO22\n2/jiF7/I7Ows+/fv58YbbwTgnnvuYdu2bUxMTDA+Ps6nPvWpxXNSDwrHNZCdX6gc59RzLCIi8nLw\na7/2a/zd3/0d3kJl7N577+XXfu3XALj66qt5zWteg+M47Nixgw9+8IP86Ec/WnOf3/jGN/iN3/gN\nLrvsMuLx+JLKMMANN9zA5ZdfDsArXvEKbr311sX9rlVpfeCBB7j44ot53/veh+M43HrrrVxyySV8\n61vfAooh/dd//dfZvXs30WiUW265hSeeeGJd5+SRRx5hfn6e//7f/zvBYJA3v/nN3HzzzXzta18D\nIBwOs3//fmZnZ2lvb+eqq65afHxsbIwjR44QCAS4/vrr13Xcc6VwXAOLlePSn5iIiIhIDRhTnT8b\ncf3119PT08M//MM/8OKLL/Jv//ZvvPe97wXgwIED3HzzzQwMDNDe3s7v//7vV9SiMDY2tuRCt+3b\nty95/tFHH+XNb34zfX19dHR08IUvfKHi1ofR0dEV+9uxYwejo6OL97ds2bJ4OxaLkUqlKtr32cdY\nfqHejh07GBkZAeDv//7vefDBB9m5cyd79+7lkUceAeDDH/4wu3fv5q1vfSsXXnghn/70p9d13HOl\ncFwD2XTxb1WORURE6sfa6vzZqPe///3ce++93Hfffdx000309vYCxenZLrvsMg4ePMjMzAyf/OQn\nK5rBYmBggKNHjy7eP/s2wHvf+17e/e53c/z4caanp/mt3/qtxf2u1YawdetWXnrppSWPvfTSS2zd\nurWif2sltm7dyrFjx5ZUsV966SWGhoYAuOaaa/jmN7/JqVOnePe7380tt9wCFNtF7r77bl588UW+\n9a1v8cd//Mc89NBDVRvXWhSOa+B0W4UqxyIiIi8f73//+/nnf/5n/vIv/3KxpQIglUrR2tpKPB7n\nueee4/Of//yq+7DWLobJW265hS9/+cs8++yzpNPpFW0VqVSKzs5OwuEwjz32GF/72tcWQ3Fvby+O\n4/Diiy+WPM7b3/52Dhw4wNe//nXy+Tx/+7d/y3PPPcfNN9+8ZCzr4bou2Wx28c+1115LPB7nM5/5\nDLlcjn379nH//fdz6623ksvl+OpXv8rMzAyBQIDW1lYCgQAA999/PwcPHsRaS1tbG4FAYPG5elA4\nroFMxhIxBby8KsciIiIvFzt27OD6668nnU7zrne9a/Hxu+++m6997Wu0tbXxwQ9+kFtvvXVJZXf5\n7dP3b7rpJu644w5uvPFGLrroIn7hF35hybaf+9zn+NjHPkZbWxt/+Id/yK/8yq8sPhePx/n93/99\nrr/+erq6unj00UeX7Lu7u5v777+fe+65h56eHu6++27uv/9+urq61hzXcsYYRkZGiMVixONx4vE4\nLS0tjIyM8O1vf5vvfOc79Pb2cvvtt/OVr3yFiy66CID77ruPCy64gPb2dr74xS/y1a9+FYCDBw/y\nlre8hdbWVq677jp+53d+hxtuuGH935ANMvWaGsMYY18uc+39yQem+MRXW9jdkuWR6bZGD0dEROS8\nYIzRvL2yaLXXw8LjG65QqnJcA9k0JII+XkGVYxEREZFmonBcA65rSYQVjkVERESajcJxDWQy0BpV\nOBYRERFpNgrHNeC60Bq3eL7CsYiIiEgzUTiuAdeF1haL5+v0ioiIiDQTpbcacF1DR6slt/ELJUVE\nRESkARSOayDrQVsb5KxOr4iIiEgzCTZ6AOcj1zO0dxg8ivPvrbWEo4iIiFRGv1Ol1hSOa8DNQXuH\nxcPB5i0mpB9kERGRc6UFQKQe9Ll/Dbg5Q2uboYBDPuM3ejgiIiIiUqGy4dgYs80Y80NjzH5jzNPG\nmP9cYpu9xpgZY8zPF/58tHbDbQ5eDqJxCOGTndW7XBEREZFmsVZbRQ74kLX2CWNMAnjcGPPP1tpn\nl233I2vtu2ozxObj5g2xFkPI+GTmfNoaPSARERERqUjZyrG19oS19omF2yngWWCwxKZqqj2LlzdE\nWwxhY8nMqa1CREREpFlU3HNsjNkJXAU8uuwpC1xnjBk2xjxojLmsesNrTm7BEEsYQo7FnW/0aERE\nRESkUhXNVrHQUvF3wH9ZqCCf7WfANmtt2hjzduCbwEXVHWZz8QqnK8c+mZQqxyIiIiLNYs1wbIwJ\nAX8P3Get/eby5621c2fd/o4x5nPGmC5rbXL5tnfeeefi7b1797J3794NDntzc32HWCuEA5bsvC7I\nExEREamVffv2sW/fvqrtz5SbM9AUZ9r+a2DSWvuhVbbpB8attdYY8xrgG9banSW2sy+X+QkHgln2\n/V/DL721wP/60wJ7b2tt9JBEREREXhaMMVhrN3w93FqV4+uBXwWeNMb8fOGx3wO2A1hrvwD8e+C3\njTF5IA3cutHBnC9yvkNUlWMRERGRplM2HFtrf8LaM1p8FvhsNQfV7FxriLcZhWMRERGRJqPlo2sg\nh0OszRAOWty0wrGIiIhIs1A4rjK/YPFYaKsIQjbT6BGJiIiISKUqnudYKuOlfRwswaAhErK4Csci\nIiIiTUPhuMrS0z5hinMbR0KQzaitQkRERKRZqK2iyrJzlvDC5CHhsMVVOBYRERFpGgrHVZae8wk7\nxXQcCUMm2+ABiYiIiEjFFI6rLJs6UzmOhMFTOBYRERFpGgrHVZaZtYQDxVaKSBiyboMHJCIiIiIV\n0wV5VZZO+UScYjgORyyuwrGIiIhI01A4rjJ3HsLBYjiORlE4FhEREWkiaquoskzKEgkUb0ciBs8z\njR2QiIiIiFRM4bjKsmlLeOGsRmLg5ho7HhERERGpnNoqqiwzb4mEircjUfC8xo5HRERERCqnynGV\nuWdVjqMxg5tTW4WIiIhIs1DluMoyGYiGF6Zyi4GXVzgWERERaRYKx1XmZiC80FYRjRs89RyLiIiI\nNA2F4yrLZiAaWagcx8EtqHIsIiIi0izUc1xl2awlEi7ejsYNXr6x4xERERGRyikcV5mbhWikeDue\nMOQKjR2PiIiIiFRObRVVlnWLU7gBRFscPLVViIiIiDQNheMqc12ILFSOIwnwfIVjERERkWahcFxl\nWc8QjRVvxxKOwrGIiIhIE1E4rjLPKy7+ARBNGDyrUywiIiLSLJTcqmxJ5bjVaXg4zh7NNvT4IiIi\nIs1E4bjKvFxxZTyAaKshZxvXVmGt5Yd7fkohqykzRERERCqhcFxl2Zwh1lIMxLFWhxwO1tqGjKWQ\n8bnFey3pcYVjERERkUooHFeZlzdE4ws9xzFDDgffbUw4zkz5pAgxcVwrkYiIiIhUQuG4ytz8mcqx\n40AAn+yc35CxzCWLFePkaGOOLyIiItJsFI6rzMsboi1n+ozDWDINCsfzyWLFOnlC4VhERESkElo+\nusq8giGWOHM/ZPyGheO5qYXK8UmFYxEREZFKqHJcZa5viJ4VjsOOJTvXmJ7j+enicacmGnN8ERER\nkWajynGVeb5DrPXM/ZBjyaYaFI6nFsLxpMKxiIiISCVUOa4yzzfEWs/qOXZ8svONCaep2WI7xcxU\nQw4vIiIi0nRUOa4yzy6tHIcdSybVmJ7f9Gzx7+mZxi1EIiIiItJMVDmuMs8a4u1nTmskaHHnGzOW\n+YULAafnGnN8ERERkWajynEVWQs5HGJtZ9oowgFLNt2gnuMUhPCZndd7IBEREZFKKDVVkedaAlgC\n0bN6joONC8fpFPSEc8ym1VYhIiIiUgmF4ypKz/qE8DFmk4TjNPTH88xm9W0WERERqYRSUxVlZn3C\nZmkQDgfBzTRmPOk09LcXmM0FGjMAERERkSajcFxFmVlL2CydmSISsmTTjRlPOgNbei2pvMKxiIiI\nSCUUjqsom7KEnWWV4xBkM41pq8hkYHAA5nyFYxEREZFKKBxXUXrOJ+IsqxyHwc02aDyuoWeLwceQ\nnmvMXMsiIiIizUThuIqyKUtoWeU4ErYNC8cZ15Bog4TJkzxeaMwgRERERJqIwnEVZVKWSGBZOI6A\n6zZoPJ6hpc2hNVBgckThWERERGQtWgSkirLzlvCycBwOg+s2puc4mzO0tBtaQwUmxzTXsYiIiMha\nVDmuomzKEgkuqxxHG1g5zju0tBvaIj5TJxoT0EVERESaicJxFWXmi4t+nC0aAddrTNU2mze0dBra\nYz7JcV2QJyIiIrKWsuHYGLPNGPNDY8x+Y8zTxpj/vMp2f2aMecEYM2yMuao2Q938shlLNLS8cmxw\nvQaNp+DQ2hmgrcUyNaHKsYiIiMha1uo5zgEfstY+YYxJAI8bY/7ZWvvs6Q2MMe8Adltr9xhjXgt8\nHnhd7Ya8eWXmLeEV4Ri8RoVj3ylWjlst08nGjEFERESkmZStHFtrT1hrn1i4nQKeBQaXbfYu4K8X\ntnkU6DDG9NdgrJuemy3Oa3y2SAzcXGPaKlzr0NoToKMNpqcbMgQRERGRplJxz7ExZidwFfDosqe2\nAsfOun8cGDrXgTWjbGZlOI7FDV6+/uHYz/lkKV6Q195hmZmt+xBEREREmk5F4XihpeLvgP+yUEFe\nscmy+y/LBlc3A9HwsgvyYgYvV/+x5Od9cjjEYobOLsP0nKZyExEREVnLmvMcG2NCwN8D91lrv1li\nkxFg21n3hxYeW+HOO+9cvL1371727t27jqFuftmsJRJZ+lgkDm4DKsepZIEIBmMCdPYYZtOamERE\nRETOP/v27WPfvn1V21/ZcGyMMcCXgGestX+yymbfAm4H/sYY8zpg2lp7stSGZ4fj85HrQmxZOI62\nGLwGLE6XmvSJOj4QoKvPYTajyrGIiIicf5YXXO+6665z2t9alePrgV8FnjTG/Hzhsd8DtgNYa79g\nrX3QGPMOY8xBYB749XMaURPLutDZtfSxaNzg5es/ltSUTzRQnNu4c4th1g3UfxAiIiIiTaZsOLbW\n/oQK+pKttbdXbURNzHWLU7edrVGV4/kpuxiOuwcCzOXUViEiIiKyFiWmKnI9Qyy29LFoAjy/AT3H\n0z7RhdX6uocCzBXWbC8XERERedlTOK6irGeIRJcG4WjCaUg4np+2xBYWJGnpC2AtZLN1H4aIiIhI\nU1E4riIvB7H40seiCUPOr/9pnp+1xELFtopQW5AEeaanXpYz7ImIiIhUTOG4irI5QyS+rHLcYvBs\nAyrHc3Zx5gwTMCRMnsmRBjQ/i4iIiDQRheMq8vIQXVY5jrU65GwDKsdzlljkTKW4NVhQOBYRERFZ\ng8JxFbk5Q3RZ5TjW5uCtWECw9tIpiEfPCsdhn+QJv+7jEBEREWkmCsdV5BUMscTyC/IMORz8fH37\nfdPzdsnMGW1Rn+QJ9RyLiIiIlKNwXEVuwVkRjkMhgwG8+fpWbdNpQ/zscBzzSZ5S5VhERESkHIXj\nKvIKhljLyhaKED6Z2XqHY0u85cz99oRlelKVYxEREZFyFI6ryPMN0dYS4dhYMnN1DscZszQct1qm\nknUdgoiIiEjTUTiuIneVcBwxPtlUfau2mSy0nBWOO9phZqauQxARERFpOgrHVZTzHeKtK09pyLFk\n5+ocjl1D/Kyg3tEJ07P1nzVDREREpJkoHFeRZx1iJcJx2PHJ1PuCPM/Q0nYmDHd2G2bmFY5FRERE\nylE4rhJrIYch1l4qHNu6t1VkvaWV485ew2xa324RERGRcpSWqsTzIIAlGCsRjgMWd76+48nkDC3t\nZ8JxV5/DbEaVYxEREZFyFI6rJJuxhPEx4ZUBNBywZOfrXDnOOyQ6znx7OwcMc16grmMQERERaTYK\nx1WSnvUJ4WNMiXActGTqHI4zBYeWjjNj6R4IMJdXOBYREREpR+G4SjKzPmFTOgCHg5BN17lyXDAk\nus58exN9QayFbLauwxARERFpKgrHVZKZs0Sc0jNShIMWN1Pf8WR9h9buM5XiYEeAhMkzPV3fcYiI\niIg0E4XjKsmmLKFVKseRkCWbqV/l2FqLS2BJ5TjYHiRh80xNaQlpERERkdUoHFdJZs4SCaxSOQ6B\nW8dw7Gd9sji0JM70HDshh4TJkzxR3/mWRURERJqJwnGVpOd8wk6ZynEde33dGR8whEJLH28NFUiO\nFeo3EBEREZEmo3BcJdl5Szi4SjgOg1vHcJxKFoiYAssnzmgN+0yOqXIsIiIishqF4yrJzkMksMps\nFXUPxz7REhcHtsV8psbVcywiIiKyGoXjKsnOWyKrVY4j4Lr1G0tqyidaIqi3xy1TEwrHIiIiIqtR\nOK6STNoSDpZ+rt7heH7aEitxcWB7q8/UpMKxiIiIyGoUjqvETVsi4dLBMxoF16vfWFLTPtHQyrF0\ntMH0VP3GISIiItJsFI6rJJOGcIlAChCJglfHcJyegVhoZeW4ox1mZlcuby0iIiIiRQrHVeJmLdFw\n6eci9a4cz/rESlSxO7pgZq5+4xARERFpNqt0ycp6ZTOG6KptFQY3V7+xzM9ZYpGVY+nsNszMq3Is\nIiIishpVjqskmylO2VZKJAZern6hND0HscjKxzv7DLNpfctFREREVqOkVCWua4lGSz8XiYNXz8rx\nvCUeXVk57up3mHX1LRcRERFZjZJSlWTd4pRtpURjhmy+jpXjeYjHVj7eNeAwlwvUbRwiIiIizUbh\nuEpcl1Urx9G4watnOE5DLL7y8URvECxk67han4iIiEgzUTiukqxriKwWjlsMXqF+4TiThpaWlY8H\n24MknDzT03UbioiIiEhTUTiuEi8HsRKtDFD/cJzOGuIlwnGgPUAChWMRERGR1SgcV0nWM0RipQNw\ntMWQq2flOAvxxMrjBduDtPgKxyIiIiKrUTiuEi8Hsfgq4Thh8Pw6Vo5dQ0vryuM5UYcEeZKnVq6e\nJyIiIiIKx1Xj5gzREhfBAcQSTl3DcdYztLStPJ4xhtZQgeSYwrGIiIhIKQrHVeLmDdGW0gE4koCc\nrd+pzuQMLe2lx9IW8UmeUDgWERERKUXLR1dJMRyXfi7W6uCVXlm6JjI5Z/VwHPOZOlXHwYiIiIg0\nEYXjKvEKhlhL6epwrNUhV8c8mi0YEp2lx9LeYpmaUDgWERERKUVtFVXiFQyxEjNEAMTbHTwcrK1P\nKM0WnNXDcatlKqlwLCIiIlKKwnGVuAVDdJVwHIoYLJB36xeOW7pKj6WjDU3lJiIiIrIKheMq8Xxn\n1cqxMRDCJzNbnwvhstahtStQ8rmOTpierd/MGSIiIiLNROG4SjxriJaYPu20sLFk52pfOfZzPlkC\nq16Q19ENM3MKxyIiIiKlKBxXSc53iLeufjpDxiczV/vKcSHt4+IQX2VBks4ew2xa4VhERESklDXD\nsTHmr4wxJ40xT63y/F5jzIwx5ucLfz5a/WFubtaCh0O0TDgOG0s2VfvKcXqqQAhLoHRXBV19DrMZ\nvScSERERKaWSqdz+N/DnwL1ltvmRtfZd1RlS8/G8Yk9xMF6mcuzYulSOU0mfiOOz2vueri0Os56D\ntcVeaBERERE5Y80SorX2x8DUGpu9rGNWNlsMx06kTOXY8cmmaj+WVNIn6qwewlt6AhiKYxYRERGR\nparx+boFrjPGDBtjHjTGXFaFfTaVTNoSxseEy1yQF7Bk5mtfOZ6ftkQDqx8n2B6kNVDQdG4iIiIi\nJVRjhbzCppmrAAAgAElEQVSfAdustWljzNuBbwIXldrwzjvvXLy9d+9e9u7dW4XDN15mzi+G4zJ9\nCuGAJTtf+57j1LRPNLj6cYLtQRImz/R0iIGBmg9HREREpKb27dvHvn37qra/cw7H1tq5s25/xxjz\nOWNMl7U2uXzbs8Px+SQ96xNyygffeoXj9IwlFlq9chxoD9Bic6oci4iIyHlhecH1rrvuOqf9nXNb\nhTGm3yyUTI0xrwFMqWB8PnNTEDHlWybCQXDTtR9LatYnFlqjcuwrHIuIiIiUsmbl2BjzdeAGoMcY\ncwz4OBACsNZ+Afj3wG8bY/JAGri1dsPdnNJzPuG1KsdBSzZTh8rxLMTCqx8n0BKgpZAnObn6jBYi\nIiIiL1drhmNr7XvWeP6zwGerNqImlE1ZwoHywTcSsrjp2ofj+VlLLLL6cYxjaA0XmDpZ+7GIiIiI\nNBuVDqsgk7KEy8wQARAOQTZT+7GkUxCLlt+mLeqTPFn7mTNEREREmo3CcRVk5y2RSirHdZhbeH7e\nEo+VH0t7zGfqVHUrx9bChReiXmYRERFpagrHVZCZt4TXaFCJhMHN1qHnOA3xWPlt2hOWqcnqjmVq\nCg4dghMnqrpbERERkbpSOF7Dvn3woQ+V38adt0TKzC0MEAnbuqxKl05DLF5+m/Y2y/Raax6u05Ej\nxb9PnarufkVERETqSeF4DQcOwOOPl98mk4ZwmenTACIR6tJWkclAS0v51bw72mFmprrHPR2OJyaq\nu18RERGRelI4XkMyCSdPlt8mk7ZE1wrHYfDcKg5sFemsIZ4ov01np2F6tnyAXq/Dh4t/q3IsIiIi\nzUzheA2VhGM3Uwy/5USiBter3rhWk8kaWhJrVI67YXa+uuH4yBGIx1U5FhERkeamcLyGZLLYglCu\nXzibsUTWbKuw9QnHHsTbym/T2WuYSVc/HF/1Cl+VYxEREWlqCsdrSC4shD0+vvo2btYQiZTfT7Fy\nXN1AWkrGMyTayh+ns99hNutgqzhhxeHDlq2PjzJ+vFC9nYqIiIjUmcLxGk6H43KtFdmsXTscx6hP\n5Tjn0NJe/tva0h3AMeWr4ethLRw5DBfnZxkf1cp7IiIi0rwUjteQTMLWreXDsetCdI1wHI0a3Fzt\nK8fZvKGlo/xxgu1BWoOFqi3YMTkJIccySEZtFSIiItLUFI7XMDkJl166RuXYhcgaSzZH4uDlqju2\nUjJ5h0RH+W9rsD1IwqleOD58GLa25ekgx0Sy9m8ARERERGpljXXdJJlcOxy7riG6RjiOxgxuvrpj\nKyXrO7R0rV05TpCvWjg+cgQGglm6On2SVZ4iTkRERKSeVDkuI5Mp9tNecMEalWNv7cpxNG7w8nVo\nq/AdWrsCZbcJtAdI2FxVw3G/l2HghjZ8v7hKn4iIiEgzUjguI5mEri7o71+jcuwZYvHywbce4dha\ni2sdEl1rt1W0FKoXjg8dsvRMpei8sYPOSEFzHYuIiEjTUjguo9Jw7OWKs1GUE2kBr1DbcOxnfbI4\nay4CEmgN0JLLMTVVnZklDj/vMxjziF8ap8PJ6aI8ERERaVoKx2UUw7Glp71Qvq0iV0HluMXUPBzn\nUz45nDX7n52gQ2uowNS4X5XjHn4Rdu0xhPvDtFmFYxEREWleCsdlJJPQ4nrMffIFTpxYfTsvD7GW\nNcJxwpDzaxuOU8kCEeNjKjhMW9SSHD/3yrG1cHTMsOuVAcJbwrTnXLVViIiISNNSOC4jmYQ28kRG\n5pmbA2+VRTzcnCESL7+vWIuDV/Nw7BN1KqsGt7X4TE2cezgeH4dowKf/VXFC3SHa8h7jJ6pTkRYR\nERGpN4XjMpJJaCVHfsyjt3f1JaS9glm7ctzGpgrHHQmYnjz3cHzkCAwGXeKXxjGOoSvuc/KowrEs\ndeQTR5j5l5lGD0NERGRNCsdlJJPQWsjhjrn099tV+47dgkN0jYvgYgmHnK3t6Z6ftkSDlQXe9jZb\nldkqDh+GPi9D/LJi6by70zJ+XOFYlpr6/hSpJ1KNHoaIiMiaFI7LSCYh4XlQgN7O1cOxVzDE1grH\nrQ6erXHleNonFqwsmHZ2WpLT5z6eF/fn2WKyhPvDAPR2W06dqM4sGHL+8MY8vJOr9CWJiIhsIgrH\nZSST0JJ2wUBPy+ozVni+IZYofypjbQ65Gp/u+RlLNFRZMN3SCyeT5z6eg08W2DHoYxauAuzZYjg1\noVXyZCmFYxERaRYKx2UkkxBPZYntidEdzpcPx23lA2EkYcjjUMjXrqo6P2uJhSvbf98Wy+ScQ6Fw\nbsc8/ILlgj1n/u39Qw6TVahIy/kjn8rz6bkL+f6Ta8wxKCIisgkoHJeRTFri01kSr0rQ6XhlwrFD\ntLV8IHQcQwifbKp2/bjpucrDcbwrSGvEP+dp146OOey+Mrh4v3+Hw+ScXlZyhjfmcYQWXhwrv6y5\niIjIZqAUU0ZysjhbRWxPjE6/dDj2fcjhEGtd+1SG8MnO1rByPGeJRyvbf7A9SG8sz9jYxo/n+zA6\nG+Si68KLj/XtCjLnnntFWs4f3phH0olwYlr/3YiIyOan31ZlJJPQ0+8Q2Rqh3c2WDMeeVwy9geja\npzJsLOnZGlaO5yFW4SfXgbYAvZHcOYXjEycgTp6eq85M8hwbDNMaKJBMbny/cn5xRz0mbJhTqeDa\nG4uIiDSYwvEqPA+yWWjvd4gMRmhPlQ7H2SyE8XEqCMch45OpZVvFPMRjlW0bbA/SHfAYHd348Q7u\nz7OFLNHtZxJ5uD9Mu8lplTxZNPZCjpx1mLBhCvP6SEFERDY3heNVTE1BR4slsiVMeDBMYmq+ZDh2\n3YVwHKmgcuxYsqnatVWk5y3xeOVtFd3m3CrHLzzsMdSexzhn+q3DW8K0+x6nTm18v3J+OXbQJxL0\nSToRvHHNWCEiIpubwvEqkknoiBcI9YWIDEaIj6eZnoZ8ful2mbQlhI8Jrz1DQ6jW4ThjiMcrmyki\n1Bui08ueWzgeLrB9YGklPNgRpM3PcXJEFUIpGjkKV+zIM2nDms5NREQ2PYXjVSST0B4qEO4PE+oP\n4U/m6OqyKyqimZRP2JyZ57eccMCSna9lOIZYS2XbRrZFaJ+eZ2xs4+M5fNCya8/Sf7cxhq54gROH\nFY6laPSE5ZLdloI1TB/NNXo4co5GRoqfmImInK8UjleRTEJbIE+4L4wTdAh1h+jrXrlKXnrGJ2Qq\nC5hhx69pOM64hpY1Vuo7rdhW4TFybOPjOTpquPDKldNzdbf5jL+kJaSlaGwiwNBOQ29LnuMv6E1T\ns/ut34L/838aPQoRkdpROF5FMgmtNk+oPwRAeCBMT6u/Ihy7KYg4lQXBcLC2leOMC/HWyrY1xjC4\nFcZGNn68Y9MhLn59ZMXjPZ2W8VEtIS1FJ2YDbN8ToK/d5/hhvWlqdseOwUsvNXoUIiK1o3C8imQS\nEnmPcF9xDt/wYJie2MpV8tIpn7BTYeU4AG662iM9I+MaWtZYjORsgzsdTp4y2A3kWC9VYDwf5qI3\nlgjHvYZT4wrHAoVMgVO5EEN7HAZ6fMaON3pEcq5GRooBWUTkfKVwvIpkEhKuR7i/GI4jgxG6QyvD\ncSZlCQcqDMdBSzZdw8pxztDSXnk4br8gQixkNzQn8aF/ydAezBNPrHwJ9Q0aJia1hLSAd8IjGYwy\nNGTYsgXGTuh10cxcFyYm4Lje5Kzp0CH43OcaPQoR2QiF41Ukk9CSdgn1LbRVDIbpZOUqedl1hONI\nqNbh2CHRXvm3NLItQl/LxlbJe/5fPba25Us+1ztkmJxRCBLwRj1O2TBbt8LWbXBCb5qa2un/K1Q5\nXtt3vwt/8ReNHoWIbITC8SomJyyJrEuoqxiOI4MROnLuysrxvCUSrDAcBy3ZTO3CcTZvaOmoPHxE\ntkXoDuY2tBDIi8MFdgyW7h/dckGA5LxeWgKzR11ShQC9vbB1Z4Dx2ZUXcErzGB2FCy5QOK7E8LDO\nk0izUoJZxeRJS0ebxQSKYTM8GKYtnVl5QV6aisNxOAxuptojPSNbcEh0Vv4tjW6P0uW7G6ocH3rB\n54LdpZ/bsjvIlBvcUC+znF+OPV+gt6WA48DQxQ7jaS0h3cxGRuBVr4LZWcjU8P+y88HwMExOQrqG\n15mISG0oHK9i8pSlu/tMFTYyGKFtZmU4XlflOGRxs9Uc5VJZ36Gla32V4053YwuBvDTqcOErSged\nzp1hsFa/FITjh3y2dBanb9t2UZCJfAg/pxkrmtXICAwNwdatxdtSWqEATz8Nvb3qzxZpRgrHq5ia\ngu7+s5ZFHgyTmJznxIml27lpSyRcYTgOg5utYVuF79DaVfnH1pGhCB2pNKPrnHbNz/mMzAS56LWh\nks+Ht4TpIKclpIXjRy2DfcXX19YhQ9JEyI1rIZBmNToKg4OwbZtaBso5eBD6+uCyy3SeRJqRwvEq\npmYM3YNnhePeMImZDMmkpXDWOgaZDIRLZ8QVwmHI1mhlKT/nkyWwrtkqArEAffECo0fWV8nLvJjh\nRCDGhReXDuKBRIB2cpw4WvqCPXn5GD1h2Lq1eLujA3I4TB/REtLNamSkWDVWOC5veBhe+UqdJ5Fm\npXBcQqEAqayhe+hM24AJGGJ9ITrailMZneZm1lE5joCXrc3V+n7Gx8WhpWV9+9+yxTLy0voqx7NP\npZnww2zbVvp5Ywyd0QJjLygcv9ydmAywdWfxvxljoDuS49jzel00q9PheGhIoa8chWOR5qZwXML0\nNLSGfSL9S0vCkcEIvR1LV8nLZIrtEpWIRMD1atNWkZ0pYDGEKqxin7Z1h8PYyfUF6kOPunS3+ERW\nrv+xqKvFZ3ydFWk5/5ycddh+0ZlPGPpbC4y8qNdFsxodhdy9L9GdTauXtozh4eKFiwrHIs1J4biE\nZBLaQoXFBUBOCw+G6WkpLAnHngvRMiHxbJEIZN3aVI7nkz4RU8Csc/dbL3Q4mVzfKnkvPJFn+5by\nAae73efkMYWglzPf8znlhdl+yZlw3NfhM3JE05g0I2uLlePozyfpmplX6CtDlWOR5qZwXEIyCW1O\nfnHp6NMigxG6IktXycuuo3IcjYJbo3bL1JRPNLD+MNp1YZiAKU7NVKnDL1guuLD8Nj1dML7OC/3k\n/OKd8JgIRBjaduYd20CfZXRMr4tmNDMDgQAEjqbpTqUV+lYxOVn8/3THDsvWfl/nSaQJKRyXkExC\nq80RWtZWER4M0+XkloZjtxh6KxGJGrwahuNYhSv1nS26PUpvuPKFQKxvOTrqsPvK8vPV9vYv7c2W\n2rL+5guc7qjHhB9hcPDMYwMDcGKdbTzV9P/enufUqc13rprB6CgMbLHkp/O0T84p9K1ieBiuvBKS\nD04y95FndJ5EmtCa4dgY81fGmJPGmKfKbPNnxpgXjDHDxpirqjvE+ksmIZHzSlaOOwtLV8lzPYhU\nGo5j4OZqEwxS0xurHEe2RegyXsVzHWdfynIyHGPXxeVfOn1bHSaS6x6ObMDMwzMM/8Jwo4exwvgL\nHqGAJZE489jWHYaTU415T57Pw59/1uGHX9HqFRsxMgIDHQWcuEPs2DyZDMzPN3pUm8/plorph6YJ\nHpgjn1/fJ3Mi0niV/Jb638BNqz1pjHkHsNtauwf4IPD5Ko2tYSYnLQnXI9S3snLc5maXVY5N5ZXj\nGHg1muI1PQ3R0PorYpHtEbpylS8Ekn42zXg4zs6d5bfr3+6QnNMHE/WQ+nmK5GNzm656fPT5Av2J\npTNTbN3lcHKuMUtIP/9EHg+Hn/9I8yxvxOgo9MXztL+hHfdolqEhq4vySjgdjmd+PIM34rJtyKp6\nLNJk1kwv1tofA1NlNnkX8NcL2z4KdBhj+qszvMaYGPVpD+YJRJf+Eo8MRmify66oHEdjle03GjO4\n+RpVjmd9YhsJxwMROj2X0Qovnpt/Zp6xQoQLLii/Xf+uAJPpxoSgl5uf/NDnjvQryL5Uw+UXN+D4\nIcuWrqWvq+2XBDmVXeeUKlXy+PdcHCxPPd2Qwze9kRHocTzil8YJ9YTY2qt+2lKGh+GKi/LMPzNP\ndGeUwS6dJ5FmU43S3lbg7B/948BQFfbbMBOjPp1tK4NmeCBM6/T80nCcM0TjlQXeaNzg1aitIj1r\niVU43/LZTMDQ35rn+MHC2hsDM/vTnEoHGVrjOzywJ8i0p3BcD08/A8/TxtTPNtdn3MePWwb6lr4m\nt10WZDIfakiV+8lHC7y2ZYbnRhsTzpvdyAh057PELogRuzDGlkReoW8Zz4Pnn4ftqVkSVyVouaKF\nLXFP50mkyVTrc+/liW9zfb67TpMnLJ2dKx8P9YRom89y8uSZf56Xq7xyHImDV6hR5XjOEots7LT3\n91mOH67saw/9LEd/j11zPuX+i8Kk/CC5XFO/FJrC4WNOsV3goc218tzoCcPWZW+ievodsgSYPV7/\n1ob9zxre/ZY8x7IRMrOVvRmUM0ZHoTOVJnpBlOiuKP3hnNoqlnnuOdixA7zHZuh4Ywex3TH6HIVj\nkWZTfsqByowAZ6+VNrTw2Ap33nnn4u29e/eyd+/eKhy++iYnLN09Kx83jqF/i2HiJPg+OM76Ksex\nFoNXo8XB0imIRzcWRAeHDA9WMFuF9S2HD1h2vXrtf2+4NUCryTF+BLbuUaWuVgrpAsfmw7TFfH76\nmOWdjR7QWU4mHV6/a+n7b2OgJ5Tj+LMF2rdXOAdilTw/GuL3PlRg4Dseww8WeN2tibW/SBaNjED7\n1DzRC3qI7YrRezLD0WM6h2db7Df+yQzb/ts2skey9Pwsw4vH2hs9NJHz2r59+9i3b1/V9leNcPwt\n4Hbgb4wxrwOmrbUnS214djjezJJThu5XlS6qt2wNkZi1JJOGnh5w84ZYhUs2R1tMzSrH6RTEKrww\ncLmtuxxOPrX2uDKHMozHW7hgd2X/hs5QntHnrMJxDaUPpBkNtXDzvysw/JPNdZ5PzAXZfvHK1pqe\nWJ6jzxe4/G31G0s2CyPpEK98G1y8xeNn3y/wulvrd/zzwciIpX0ytVg57vp+hn9RRXSJ4WG48gqf\nuU/P0XZdGyZo6JqZYp/Ok0hNLS+43nXXXee0v0qmcvs68K/AxcaYY8aY3zDG/KYx5jcBrLUPAoeM\nMQeBLwD/zzmNaBOYmjX0DJY+NZHBCL1tZ5aQdvOGSEtl+40lHDy/NjM4pOchHt/Y1267OMDJmbXH\nNf/kPJPdrWtejHdaZ7TA+CF9fF1L88+mGclHeO9/DPDMdAzf2xyrEvp5n1NeiG2XrAzH/W0FRg7V\nd5xPPZpjq5OldVeEKy6zDD++Oc5TsygU4NQp6E0UCCaCxHbF6JxKqa1imeFhuKglS3RXlFBHiNju\nGB3jmhNapNmsWTm21r6ngm1ur85wNofptEPv9tIXk4UHw3RH85w8GeTyy8FbR+U40gI5v0aV47Ql\n3rqxtorePREKBUMqxZI5aZebf2qe8Ug31+2sbL/drT4ntFRwTR173CMWhjf/O4cjxJl+Kk3Xqxv/\nUXfuZI4JJ8K2HSvfdPV3+owcre/r4mff99jT6WJMC696Y5C/+LRel+tx8iR0tloSuyIARHdF6Rid\n45je+y6ythiOd7x+ivY3FtsoItsjdE6mODZtsdZgGrf+jYisgyaiXcb3Yc516NlZOhxHBiN0hc4s\nIe35hliiwp7jhINXq3CcMcQr7H1eLrojQndw7YVAUk+meDEV5eKLK9tvd4dlfEQVulp6/ok8OwcL\nxOOwvTXH499xGz0kAFLHXGb8IP0lJnUc7LecOFHf8Tz5mM+lO4tJ7tp3RnlhNrppquzNYHQUtrQV\niF5Q7N0K94eJZ1xyOcvcXIMHt0mMjRUDcmQ4SfsbiuHYCTp07QwRDhYXlxKR5lDXcLzZFikoZXYW\noo5PfLD0xULhwTCd1jvTVlEwxForO43RVlOzcJzJQnyDBcPItghdBXfNcDw9nObgiSBXXlnZfnt7\nYPzE5v+eN7ODL8DuPcXbr7wwx2M/2RylvOPP5OmKFAiW+GxqYAjGTtX3ffkzzxuueGXxZ2/P5QGm\nCXHi0c019d1mNjICvZEc0Z3FcGyMIX5hjMFeLXBxWvFiPMvsv8wshmOA2J4Yg10FnSeRJlLX31CZ\nQ5t/2dZkEtpMbsXS0adFBiN05M4sIZ3zHaKVVo5bHXK2Nqc84xpaKhzHcqHuEF3W4/jh1YNVYb7A\nc8cDXLCr8t7m3i0wObGhIUkFrG85POpw0ZXFBPrqa+Dnz2yOuaWPHijQ31p6apatOxxOTtc3HD93\nIsyr3lg8T4EAXNjl8dMHNteiKZvZyAj0GHexcgzF1oqB9rz6jhcMD8NlQ3mCrUGiQ2fOU2x3jC2x\nnMKxSBOp62+o1M9S9TzchiST0OrnViwdfVp4IEx7OnOmcuwb4m2VV45zNTrlGRda2jYWjo0x9LcV\nOP7c6vPMze+f52hfF1dfXfkx+oYcJqbVZFcr7nGXsWCci64ovqZe/9YQT5+INHhURccO+Qx0lW5b\nGNodYHy+GhPlVGZ2FqbdAJe/5cyE5Jdd6PPkw83VVmF927BP30ZHodMtLgByWmxXjC1hhb7Thofh\nQpNa7Dc+LbY7Rp9xdZ5Emkhdw/GJhzd/OD51okCrzRPsKP3LOzIYoXW2GI59HwoYIpVWjtuKizXY\nGvx+y3jOhsMxQH+3z/EXVw8LqSdTvBhv56qrKt9n3w6HydTmqGSej9LPpxkLxbnwwuL9a98e4Ugu\nxvx4jSbTXofRERjcUvqFPnRpgAm3fuH4iR/n2BlIEz2rVeqV1zo8/VxzvXF76Q9f4sXffbEhxx4Z\ngc659IrKcY/NKvQtGB6G7aemVoTj+J44PW5G50mkidQ1HP9sk/RDlnPqSIH2aAGzymXFwa4gHZ7L\nyTGL60IISyBa2WkMhA0BfLLp6lesMjlDS/vGf9kPDsDIsdVT+/xT8zzvtqwrHA/sDjKV0TWftZJ+\nLs1xL8ru3cX78RbD9pjLv3278e0CY+OGrdtKvx4HLg6T8QNks/Wpgv7sIY893d6Sn+mr3xLm+Ykw\nfq55qsf3fDHIJ/6/MIVM/f8fPX7c0jE9T3THWe0Cu2J0Z9IKfUAmA4cPQ8/Tp5b0G0Oxctw5M6/z\nJNJE6ppcnnzGwdaibFpFp44V6GxZ/RemMYb+/uKVya4LYePjVBiOjTGEsWRnq38OsjlDS8c5hOMd\nDidOrv71s8PzPHcyvL5wvCfEVD7UFBdiNqOx4Sx5DL29Zx57xZDHIz+o/9LMy40lAwxdUPrnIpQI\n0Gk8jh+oT8h76nGfyy5c+jP9qtcEOGxamN+frssYzpU34fHYyTjfYwvjf1//Rv7Ro5YtXT5O5Mz3\nNLorStf0vHqOgaefhot2+TjzeeKXLL0oI7IjQtfsPMde0v+DIs2iruH4BZsg+1Ljq1rlTI5aOtvK\n/yfWP2Q4NVGsFoRY+gtjLSHjk5mrQeW44JDo2Pi3c2hPgJNTpb/eWstzwwV6e6Gjo/J99g06zBLC\nm2x8WDsfHXiywAVb/SVzp159peVnw+f+Y50aPrcWqJOpINtLLAByWm84x7EyPe7V9MwLDlcse1O3\nZQuYoOHQD5tjxoqph6Z53mmDeIDv/clM3Y8/Mgrblr3Zie6MLixwodA3PAwXd7u0v6F9xaeOTtBh\naNByVHO+izSNuobjF8Ntm/6ivMlxn86u8tu0DYWJhYtztYbtesOxrUk4zhYcEp0brxxvvzzIqVX6\ng70xjxcKCV597fr2GY+DYyzJQ96GxyWre/Eg7L546ff8tW8O8ORL57aMtDfh8dOrfrrh2WVswXLK\nC7HjitX7intbChyrU+X4+VNhrrph6TkxBi4dyvOzHzbHG7en/iFFLA6/dbvhgf2Jus78k05DNmvo\nX7YMfCAWYLDH59hRanIdBUAhW9j0nzbCwsV4hbkVLRWn7bg4wOgJg988XTwiL2t1DceH5yMk/21z\nzxg/OQHdPeVDZngwTE9LgSOHLSF8TLjyUBo2PtlUDdoqfIdE18YvfttyWZhMwSFT4ndu6skURzo7\nuOqq9YfvzkiBky9s/l7zZpOfy3NsLsyeK5b+CL/m5ihH5iNkMht/jU390xRYmH14dkNfn5vIMWEi\nbNu5+n8v/e0FRo/UPimcPGnJ52H33pXzD15xpeGpJ2s+hKp45IcFXn013Po+hx85fYx8qX6rqIyO\nQl8iT2xXdMVzvbvD4MNMjYrZz//H5zn2mc3frDs8DENjyRUX453WeXGURMRnfLzOAxORDalrON62\nxefJ/9v4K+nLmZo29GwpHwIjgxG6I3leOmwJG3/Vi/dKCTu26uHYWotrHRJdG/92RrdH6LIeY2Mr\nxzb/1Dwv2NZ19Ruf1hUvMHZoc3/Pm1HmQIYTrQn2XLT0tde2I8z2YIbHz6HvePLBSWIXx5h5eGOJ\n59QBFwu0ta2+zZYen5E69Kr+/Ec5dgXTRErMW37VDUGePR7Ez2/ucl7mcIb9s3Fed2OAiy+Gge2G\n734xjS3Up6I6MgK9odySmSpOi10YZaCjNnMdW2uZemiakc+ObOrvkbXw5LBl+/gUiatKr8QU26O5\njkWaSV3D8VWvNgw/xab+mGxqztC9tfxpCQ+E6XRyC+F4ff+WUA3CsZ/1yeJseBEQgGAiSHfA49iz\nK4NsanieZ5NRrr56/fvtbvUZP7p5f7E1q/TzacYCZ6ZxO9vlvS4Pf3djrSzWt0x9b4qdH9/J7CMb\nqxwf3V+gL5an3HvGgS0wdnJDu1+Xn+/LcVFv6TcKr3pNgJeCCdLPbu6L8qZ+MMWLrR285jXFE/re\nDwT4oekj+U/1WY94dBS6fXfJHMenxXbVLvS5x1z+aGoXXzXbmbx/svoHqJIjRyAWtGx7bQwnVPp3\nh+Y6FmkudQ3Hr77O4UA+gTe6eXtQZ9IOfdvLz8EaHgzTWXA5+pIh7Kwv+EUCPtn56objwryPS4DY\nyt9d69KXKHDsmZXh+ODjHuFo8SKm9erpgvHRzftmqFmln0tzLBtZnMbtbFddUuCnj21sv3M/nSPU\nF6RIS1cAACAASURBVKLnl3tIP5umkF5/S8yxgwW2tJX/tGBwm+HkRO3/+3nq55bL9pT+Gb38cjiS\njzO9yVu9Jv5pimdnY1xzTfH+r/wK7Mt0cfQv1ljvvUpGRqArmylZOY7uitJraxP6Zv51hp86Xfz1\n+ACP3r15l9ocHoZLOrOrtlRAsXLcndVcx5vBf/tvxdmmRMqpazi++mrDoWgbc49v3l9G0/8/e/cd\nHkW5/QH8O1vTe29A6L0JUkRQlKYidrBgwd4Qe7l67e1nF9tVsCJ4EVQUBClSEjoJJIEE0nvb9N3N\ntpnz+2MSIGST7Ca7JFzP53l4SDazM+/uvjtz5i3nNSsR2qf9sbvaKC38TSbkFwIapXOBn0ZJMBlc\n25LaWCtCCYKyi+tthAdIKDpjfLBklZCarcKY8zrXKh0SBlRWcHDsarq0RtSblYiObv23cZMUOJLV\nuUU2qjZUIXhOMJQeSngP80bDQee/q4W5EiKD26/jMX0UKK93/wIxx7KVGD7Wft318wOC/SSkbze7\nvRydRRIhZbMZYeFAUNNE4d69gf6DFdi6CbBUur+hoTBfQpCpEdro1qsvesZ7IthNC1xk/mWAEUo8\n+xzw2oFwGI73zBb+I0eAPub6NifjAYBHLw8EG42csaKbVVcDH3wAbNrU3SVhPd1ZDY5HjQKOG71Q\n10ODYyKg3qZEWN/2Z/trojTw15tQWCxA62xwrAJMLj7H66sleDjZgm1PRDih6IxJUsbjRuT4+Xc6\nOA6LBHRV59ZKZGcqeLsAloqe1duRlSahV4wEhZ1v8HkzNcirVdudXNmR6j+roZ0WhNGjAdvIgE4N\nrSgplodNtCdmgBKVRvcGx0RAVrUGYy5uPd642dCBEo7s77nDfvQpepzQ+mP8xJYf9PwbBSRGxqD8\ne/ePTSnKkhARLEFQtv4ee8R7ILjePbmOd+8kjBsp4bEnFdD5+uDbx89+CjtHHE4ixFXWwG9C24Ps\nFWoFooJF5B/nycndacsWQKEg7NzZ3SVhPd1ZDY5DQgA/H8KJhJ7ZUqOvJ6hA8I1tPzhW+asQSGZU\n6ARoVM4GxwST0bWtB/oaCR7Krl/go2IElJ3R3WRINSBb49yy0acLi1Wiqu7cDY5rd9Ui56kclC7r\nOf1wJBFy8oVWadyaBY72QhyMSD7kXJ2wVFpgzDDizyJ/HD4MHPYI7lTGipIKBWLi2v/MIweoobcp\nYXHjPUd+PsFTsiFmYutMFc1GTVLhWI7irE1uc1bt1lrkhAZj3BlpFK+7Dthe4Yf8L8vcPoejuIAQ\n08Zqh5pwDYKtJhTkuvYGQzSKOJSvxZSZSmg0wNKPgZc2BKK+oucFl4cPShjeX4TKp/3eml69BRTk\n9sx69k+xcSPhco9ybP+r59Uj1rOc9bV9R40CkpLP9lEdU55jhZ/C1uakimaCICCsaVWyzgTHZhcH\nx4YacklwHB2vQJmu5UXQkGLAcYNnpybjAUBEHyWqDefmEtJEhNxncxF1fxRKl5X2mJX+TAUmlHp4\no/8g+++rykeFIT5G7N3kXMaK6k3VCLw4EF8tV2DWLCCx1Ad1e+qcDr7KaxWI7dtBxpcINQIFK8rc\nmJEs+W8r4tVGqIPavtkdeZ4C+RpfGHtol33Nlhqkm71bBcfR0cDIMQL21Pmjfl/nJk46qrRcQGw/\n+3VNEATExgIFea49ZsPBBqRrA3DBVPm4MxdoMC6yEc/e2rM+p4oKuWds2PTW47HP1GuQEkVl525D\nwbmOCNi4nnCVPh+11fJYesbactajlrGTFDjR6AVzWc9rPa7ItsFP7dgdZUSUfJLzcDI41qoIJhfn\n79fXSvB0shz2xA1WoeKMcaD5BxrRKCrQu3fn9hneR4laS89Pl2VP9aZqWKusSJnWH2aNGrU7ax1+\nblqaPLbNHYwZRpT7+tidjNdsZLwVBxKce8+r/6xGybAwlJQA778PbN2tBCA4vapluUGFuEHtt6Kp\n/FUIhhmFue5rwTm8y4qBEe3fIAwbBuQqfHrkPAjJLEGXUI+MYpXdnpsbbhCQEBaNMjfmPCYCymoV\niBvS9ufZa4ACxaWCSxcCKd9RhyyzV4ubgv/7P2DFZg8cPdozblIBYM0aYEpQHYKmtj3euFnvkWro\n6pWwcWbLbnH0KKCySRgwABjh0YBdu7q7RPZJNgnGEz3rJvCf6KwHx2PGCMjx9oc+ueetlFeZLyLA\n08HguKnbWOPkvCeNGjC7ODg21BE81F2/YMQOV6HS1LKV7fARAaOGU7tpudoTHqlArVIDq+7cWIms\nGREh97lc1N/WH9ddL2BL/75OBSGvvCDhmafILYsjGDOMKBE87aZxazZmLJB8zPExvSQSqjdVY21B\nEO68Exg4ENBoBFQMC3VqaAURocKiQVw7q+MBcotjqIcNRRnuC45TDgNDBrb/vRg4EChpVKNyX9vn\nI1GUJ/KcbfV761ESF4w+fQT42Emfe801wI5cLxSsroJN756Iq6oK8FBICBjQejJes5ABWqgVhJoa\n1x13/182DIgT4e196rFBNwTgjsBi3LtQdNuKfM5a+SPhwrqSdifjNfMb5IkAtY0zJXSTjRuBCV61\n6PNSbwwzVGP7lp7ZYFPxYwWSpyRDsvbM8v1TdMuwihMm7x7ZUlNZJCLQx7EK6R+ngbdShFbjZMux\nhmAyuXhYRT3B0wXBccwILfSSEiaj/B5Ya6zIqPfEmImdryYhIUAd1LCU9awJbR3RrdWBALy6IQB3\n3QV8f9gf5et0sNZ2HOSXlMjdd8OFOqxa6fqreOPxRhQ12k/j1mzMRRrkV6lgdLABov5APaRwT/z3\ndxXuuENeXnnmTCDJM8SpSXmmCitqoUZMBxlfADl1YFGW+4Lj9DwlRozrIGe5BoiPIaTtbju4fOIJ\n4JJLXF26jtVsqUFebEirIRXNwsKA8ycIONIvBpWrK91ShpISIERptZvjuJlHvAfCvVyX65iIsDdJ\ngUlTW96RCwoBDz6phC7Xih9/dM2xuqKoCEg7TLhorBWa0LYnfTbz7OeJUHCu4+6ycQNhlK4cwXOC\nMWmUDTt7aHCs+00HW52tR+f2/ic468Fxr16ABQrk7XZx86kLVJUQAttZ1et0mkgNAgULtO3P3Wv9\nPDVgdq6XukOGeoKntutBmEqrQIDCisJUOZA1pBqQ4xuAMWM6P04uMBAwSkoYis6d4JhEQu6/cpFx\n+QDodAI+/RSIjBJwZFgvVPzY8fqvH79pxXShHNcHlGHZUtcHf3XHGlFWr0SvXm1vEzTGG71VRhw5\n4tg+qzdUY0+vWEyeDMTGyo/NmAHsLvdxquW4KM0KX6UITcexAsIDRBQXuKcJ0GYD8mo1GDW944KM\nGCsgLV2wO6Y8IQFYtQrQ6YCkJHeUtG01W2qQIfi1GRwDTTmP1eFumzBaXAwE20x2cxw384z3RJjC\ndUFfY1Yjjkp+mHJp696HmDsj8Ij1BJ54zD29Ms5YtYowRaVD/+fiHNreo7cHQqyNyHfx5MX/BWXf\nlUFsdO5cWVgIbN7s2LYGA7B3LzBljA0qPxUmXOmBwlIBVT0s/hRNImq21CD+9XiUfsVdDN3prAfH\nggCMGE5IPtRD+sVOU1UhISjYsW21UVoE2CydaDkGTK4OjvWuCY4BIMzThsIUuRXNkGpAptW705Px\nAEChAPw1Iiqyz53ZweU/lAPBGryy0hfvvguoVMDixcDP+sgOgxCzGfjqKwH33S7hyru0yM8D0tNd\nW76cdBHRkQR1Ozdmnv09McBWj/2Jjl2Iq/+sxs/5Qbj77lOPXXwxsD9dheq0RocvXAXpNoR7OTaE\nJiKMUOKmSTGZJwghMCNsrHeH244Yo0Ce1rfVOD+jEbjjdsK/L6rANX1rsGyZe8pqj63OBkOaAanF\nmnaD46uuAnYd00B3wgJDusHl5SjIEhFCZqhD265sHvEeCLGYXBYc1+2uRxr8MGlS67+pA9WYer0W\nU2P1eP551xyvs1Z8YcOM4BoEzgh0aHuFRoFIHxtyU86NQcfOBqudVbOtBhm3ZkD3i3MLvXz4IbBw\nIRzKeLNjBzA0qBFxlwcAAEJnBmKYugEJCZ0psfvUbq2Fz0gfRN0bhfo99TAVujhY+B/z6qtApXs6\nzc5+cAwAYycqkFHrCWtVzxqHWlUFBIc61kqqidIgEBZo2x6KZ5dW6/qWY6Me8PJwUXDsJ6LwuHzy\nLj1oRIVJhYEDu7bPIB8J5XnnRnAsWSTkvZiHv8cOQN++AmbOlB+/5hogR6dCeokaDYfbHhL0/Udm\n9LbpMe2NcETND8OlinIsX+66G0FbvQ359Wr0G9B+PVWoFRgWacb+7R1fiC3lFqRmAKV1Ssyeferx\ngABgxAgBx2PCoE9ybI5AwQkR4X6OfdZRUUBZhXtm7yf/bUG8RyNU/h1PChg+HMj38IX+0KnXSCLh\n0av06F2sw2RzBS5IysKqldSp3NGdUbujFupx/sjKFjBiRNvbBQYCU6cKSJvcF4Xvur6/vuCYDRGB\nIoR2Jh149PZAkMGIQhf1AhzdZIRKKyCujQbZ6AeicWtpBn76iXD0qEsO6bQTJwiF+cBVrwa2+96c\nKSZCQr4bx9m7irXKit2Ru91yw3U6EgnZj2Uj5JoQp3J2EwG//AJ4eMiTIjuyaRMw1lyFoFnySjo+\no30wXKrFtvU960ZF95sOIVeGQOmlRNgNYSj7xrF5LkTA888DBw+6uYA9yPHj8mtetco9+++W4HjM\nGAG5/gFoSOpZ445r6gQERzp2otNGaREGM3w9nWw51sqti65k0FOXl45uFh5CKMqRWxsPHyAM6yd1\neeW9YH8J5UXnRldi6VelsPb1xbsrvfDOO6ceV6uB++8X8EdkfJsT84iAD96UcM/1ZqgD1PAe6o25\noVX4/mty2Qx143EjKkP80a9/x/V0zHDCocMdb1e9qRp/hfXGokUCVGfEkjNnAsneoajb41gfdnE+\nITLEsc86KlZAWY17TkGHE0QMinTs5nv4cCCr0RMNSQ0gIlRtqMKy/hlYs12LL37xwLCfh2HYDX4Y\nGmzC2rVuKW4rNVtqUDwwDEOGoMMb8Pnzga2GYOjW6mAqcu2dd2G2hKgOFnRReioR6WtD/gnXfMd3\nJwATxrY9Cdh3jC9CYxS4a6YRr7/ukkM67du3TLjIuwqR14U69bzYPoLL0965Q+nyUkgmCWXfujHX\nIuReOotGiYcqh+J4otnhLFZpaYBkI7y6xISlSzve/s/fJYwTq+AzUp7ZKigETJlA2LGl5/Rgk0So\n+r0KwVfK3deRd0aibHmZQylEf/gB+OQT4J575AnE/wTLlxMGBjTip+/c84K7JTgeNQrIsnk73Bp1\nttQ2CAiJduwt0URpcCdycP1I53KMarWA2cXDb41GwMtFwXFkFFBaRCCJcCRTiTETut6yFxIMVJb2\n/OBYNIrIfzUfP0X0x9VXA0OHtvz73XcDf2X7IOOHKrtdjttXNaK6VsAtH58amzPuZl9EaizYuNE1\nZTRmGFHm7d1upopmIyerkF+uhKGDxp/i32vwZ3kgFi1q/bcZM4DdOsfHHReXCIiKdGhTxPRVorLB\nPavkpaYShgx27MLXqxdgsCmRu74ORy45gqNLcvCGeQA+/1GFPjN9AcitlZdWFWDZV2fnYlqzpQaZ\nXgHtDqlodsUVQMI+BTzmR6HoXdcuVVdSDMT06vgcEBsDFOR0/Ttuq7chudQDF85qv15EPxCN2SU5\n+OsvIDOzy4d12k8/ATffq7K7amB7ep8DuY5JJJR8VoKI9wai/Ptyty2QIxpF5P4rF3+OGozEPQIS\n+vZCxaqO53QAcqvxtOA69P7oCAoLqd35ALm5QI2OMP4yLQTFqff+gmu0OFGoREMPaaNrONAAVZAK\n5mAvrFkD+IzxgSpAhZpt7aeByckBHn0U2LYN8PbGWR3+1V1sNuDbrwiP16YhLVWeOOxq3RIcDxoE\nlBlUKN/n3i4bZ9UalQiLcyw3m8pXBW8fAR7ezr2FWg/Hxkg5o9EIeHU8tNIh0b0VKK0QYMozIVvl\ni7ETuh68hIa5b1yQKxUvLUbdyFCs2qjBSy+1/ntoKHD1tQI2hfSyOz7u3afNWDS7EdrAU3Uo7Pow\nzDAV42sXDa0wZhhRQp7tZqpoFjjaG329TDh8uO1tSCSs/VOJCeeT3W7sceOACr0K2QmNDi0GUlYp\ndLg6XrPIfirUWxSwumF0VUaBCiMdrLuCAAwdQsgVvRF2Qxh+u+w8jL9QiWuuOfU6fEb6YMaQRqQc\nImRnu768pzMXm+WhLqVah4JjX1/g6quBlRSHsm/LYKl03QmmtFJAzICOz3Fx/QQUl3Q96KvfX490\nbQAmTWn/mKHXhQIp9bhnvuWstx7v/ckAvUmBy/4d5PRz40erUVrr3mXTu6p6YzUM/h4Y80QETvgH\noXqze/IYFr5bCNPoIHz6swe++QbYUBWMsu8cG1qx9mfC2KwCqFTAwmlGfPJJ29tu2gRM8K9HyOyW\nn1fE7EAMFBqQmNgzWo+bh1S89568AuYnnwiIvDOy3Yl5Vitw003As8/KjY4ffSQPNeiO1JNn08aN\nQAQ1Yvq/QzCBqvDfFa5vPe6W4FilAoYMICTt6znt/0SEOosCoX0dP3FpojRQaJ17Cz08ALPFtS0H\nRpPrguPYAUpU1ChgSDUgS+nbpcl4zcJ7KVFRTG7LxXq6jAxg/HigoMC559nqbCh8pxCfi33w2GNA\neLj97RYvBtZUhaHgy5YnrKytBuwo9MEjX7RMd+I91BszQ2uxZbOc8aCrjMeNKDRoHGo59h7ujf7W\nemzZ0vY29fvq8QdF4t6H7dd7pRKYfimwvzEA5sKOuzxLa5WI6evYd8IjSo1Apc3lq+Q1NgKlejWG\nX+z4hIARoxWwPjIIecOi8ONKAR9/3HqbPoujMNNPh6+/dmFh7ajZWoPAiwNx4KDgUHAMAG+9Bfyw\nVoWKi2JQ/JHrZjmW1SnRe1jHKXl6D1WhpFrR5fzDJX83oMDi0eFy9UoPJaLuicKVhgKsWye3Dp4t\ny/9twLwLzVB5OH/57D3OA3UWpcuH1rlS8dJi7B7SB15ewAplL5R/6/hYYEeZS80o+qAIyzV9ceed\nwIIFALQKpBSqYTjWfqNZXh5QlCvhwplK9H6hNy7OycHatWgz88TGPwmjqsoReGnLiZOe8Z4Y5d2A\nrWt7xtwn3W86eM0IxhdfyMHf228DfynDUb2xGhad/RveV14B/Pzk6xIgB8jXXAO88MJZLHg3+PIj\nG2ZZSxH3VByuOK8RK79wfWzRbev6jpmoQLrOw6G8sWeDaBDRQGqERjseHGujtBC0zgW6Gg/XD6to\nNAnw9nFNwB07RIUKgwq6QwYUNGoxbFjX9xk1UI3GGB/kvZDX9Z114OmnAa1SwowZzrVW57+ej9zz\nYnD4uAqPPNL2diNGAIOGK7D+kCcac07NznrvAQOunGBCsJ3602d+EKbF6LFihTOvxL6GdCMKdErE\nx3e8rTZGi7nKEnzzNeHWW2F3kYY93zagXOGJyy5rez8zZghI9nNsMZAKowpxQxz7DmnCNAgms8sX\nRUg/SohCIwJGeDn8nOHDgf37gdtvB5YulfNznylkXghmW0rw9ZeSW8f11WytgWJSEEpK5F42R4SF\nySvIvZwRh/xPi2Gr7/rFwmoF6i1KxI7q+CYjaJAWWgV1+QZw9182DO8rOjTROeqBKJh/KcNdC0W8\n+WbXjusofboB6zN9seiVzrVGePf1QDAsKMzrmcPMjFlG1B9swOoUP3z3HXBM54E9f5hcfp3O+3ce\nqub0wpZEFZ57Tu69uekmAQlxvTqcmPfLWsJEoQq9lsQg9LpQeOQ1YPYFVrvDCaxW4O+twNQBZru5\nqKdM6hnjjo1ZRtiqbVh33A/jx8vD2TZtAp59SYUjI3vJGZTOsGsX8OWXwDffyFmhmr3yCrB6NRxO\n43muqagA/t4hYMHdSig9lbj2JX9k5CpRXOTaz7H7guMxAvICAnrMSnnWcivqBRWCgx0PMjvTcuzp\nJcDi4vsBo0mAl69r9hU3VA2dqEbSNgv6RTt2kepISAhgG+KP8h/LUX/AuTHazkhMBA4minj+UAKu\nmGjGnDlwaDxZ9ZZqlH5fjg9KYvHGG+hwcuPiJQJ+9e6F0uVyVFe1twH/zQzAkx/bf2LY9WG4pKoA\nX3/dtS8viYSibEJICODlQNwnCAJGjgR2floLPz95qeR161pu880vatxyta3VRLzTzZgB7Kv1RW0H\nwTERodKqQZ8RDiQ5BqAOViNIMqOo0LWBQtI2C/p5NULl4/jylcOHA999B4weLbe82KNQK3DBQwEI\nITM2bXJRYc9ARKjZUoOc4CCMHo12P5cz3XILEBatwO+x/VHyedcH4ZWUEAIFK7z7tZ3juJlnvCfC\nlV3LdUwSYX+aChdc7Ng5VRuhRcjVIbjeoxirV+OsLK6x4clKaANVOG9S54ZGKLQKhGstyD7YM/O+\nl3xWgsrZvWAwypl6Hn9SwH/9+6LyJ9eNi9On6VH5qw7vZEfj5Zfllk9AHh6wId8PxT+UtzsJbfVy\nG6aH18Nvoh8UagWiH4zGPCrCp5+2noy2Zw8Q52tB3yvsL2Bw0XwPpBWoXJ5e1VlVv1Uh6PJgfPiR\ngCVL5McGDwb++AN4KSUG6z40tBjWVlsrf9+//BKIPGOOR3Aw8NJLwEMPocesJOlK3/5HxGToMHCJ\n/MLDpvtjsl8tvn/DtUtud1twPGoUkEU+PWZSXl2hBQTBqawPwZcHw3ecc1GphydgsTnfyitJsJtG\niiSCUU/w8XNNy3FEpIA6qHFgv/wZuUJoKFBVr0S/d/vh+J3H3bIsJhHwxBIJt1myMeC5OFy1OQmj\nBouYN6/97CDmUjPSb8nAnmtHQalVYP78jo91xRVArVKDrf/Rg0TCl/fWYFA/CSPG2r9geg/1xrhg\nA6rLCcnJnXyBAEx5JpT5+6BvP8c/a5/hPjBvrsJ7L8qrij36qHwRqqoCanPN+FMXhPv+1X4wGxcn\nT6rct6X9u7q6YhssUCA01sGxvkp5CeniTNc2wx5MEDEwxrmW05Ej5eE49oZTnC7yrkjM1BfhP5+5\noQ6LhIqVFVBoFDhS1H5+Y3sEAfjiC+Cb/DDs+T9dh3lqJZvUbh7V/DQbQhQWqPw6jtA94j0QajOh\nqAvzAY3pRhwV/DHlUscDz9hHY2H6pgi33yrh7bc7f2xHmApNWLNZjRvvULSZScMRUQEicg73jB7T\n04lGEWXflmG9LQJ33CG3Rt59N5Bs8MXuT2tddpzsx7Nx9PJBqDcocMcdpx7v3x/o00+BZHUQanfa\nP15lJZCaocBVz/qeTKEXdXcUohKLER4sYf36lttv3AiMp1Mp3M4UMycAvcmAvbu7tyVf95sOqbGR\nUKnk/PLNxo4Ffv5VwPOFffH3cjn4IwLuuw+47DLg8stPbUtEkMwSbHU23DrXgoZqCd+8Z0LD4QbU\n76+HPk0PU4HcC+CuSZbuRgR8uVTCjRc0wiNGvmkXBAE33CxgtYtTunVbcDxiBJBdp0X1wZ4xVVSX\nK8JfIzp10gtfEI7AaY4lgG+m9RRgtjp3ZpUkOSALDAQmT5YH32/aJLeK5r6QC7NCieBRjncht0el\nAvzVIg7YAjFuqmsmjoSFASdOAJo5YdBEatySj/W33wDdCQsW3qVE73/3RuSiCNyTk4LAAMJNN9lP\nb0Mi4Y95uXjaazSWb/bCl1+27J5qi1IJPLREgTViDPJeysMP6YF47LX2A8zw60Mxr3dtl8arGjOM\nqAzxc2gyXrPIOyNhzDBib5+98HroANZcmgXfOiOGDSE8cJuIEWEmxPfv+EXPulzAzhMeEE1tB1z5\nqVaEqC1OfYfCfEUUZrvuwpSZCfx3qxZzpzjXMhcYCOzbJ9/ItUcbocW1l4v4ewtQ7qKhmJZKC/Lf\nzMfevntR9GER+n/WH4cOOT7e+HR9+gD/ekHAuzQApcvbHsxt09uQNi8NB4YcQGOu/eTNeSkWhPs6\ndpOhCdcghEzI68KNTk1iHY6Jvpg40fHneA/1hs8oHyyMq8SKFXD5EJ3T5f1fIXaow3HT7V07L0ZH\nEvIznK/z5eXArFlyGjN3KP+xHOpxAVi7UYXbbpMf8/YGFj8GfHUiBMbjXW+Zq95UjbosM97cFoQP\nPkCrNKE33wz8HRDT5tCK1Z+ZME5dg9ibTn1R1UFqhM0Pw/w4Xau0bhv/kDDGoIPv+fYbsdTBaowN\nNWLLj903CNyis0B/RI+vEnzxyCNodf6cOlXAuzfX4PoHtUhNBb7/HkhJQYtUow2HGrB/4H7s8tuF\nPXF7kDT6AO7VHcNTTyuQdMsJZD6QiWM3HEPy5GTs7b0XOzQ7sNN7J3ZH7kbKnJR2z+s9yb49hMYq\nG654tWXcdd3LAcis0SB7jwsT0RPRWfknH6qlgfEifd8rpdXj3WHrC+XUL8Dk9uOsfb2WxgXWO/Wc\n994jmjCBqLaWaPNmouefJ7rwQiIvrUiDNQ0UGS7RX3+5royDAhvJQ2GjXbtcsz9JInrwQaJhw4hO\nJBppV/AuMmQaXLNzIrJaifpHW+m9qHSyGW3yMUWJUi5PobT7Mmn6dKK775bL0UynI7plXB0FqK30\nwfsSWSzOHbOmhsjfS6RXkEKxwVay2drfviG1gX6OPEQhIRKZOlnNCt4toLvGVNFrrzn/XNEiUt3e\nOsp/K5+OzDlCn3gdpnihgVY8UuXQ89evJxrtU0+1u2vb3GbdO/U02r/BqXI9P7iAbp7V6NRz2mK1\nEo0/T6IlwblUtcmx19UZtXtqaY5vBb39ttTxxm2QJInq9tbRsVuO0a6AXZR+WzrVHag7+ffYWKLM\nzM7t22YjGjPYSs8EZ5FoEVv93VRsogOjD1D6onTKey2Pki5MIsnW+rW8urCOFvTVOXzce8MKacnt\n5s4Vmoh+vyqbeoVYnX5e1V9VtG/oPnroIYkefbTTh2+XudJMH/kcoRFDWr+fznrpqmq6cXiNPwQM\n+AAAIABJREFUU8+x2YimTxVpamwDhYdKdPBgl4vRgiRJdGDUAVr6aANddlnLv9XWEgV62Gjbvfld\nO4ZNov3D99O/bm6gefPsb1NRQeTvJ9GmgMST5/LTTYttoA+urWj1uD5dT9tCEyksTKL0dPmxsjIi\nPy+Rkq5KbbdcX1xTQhfEG9svuyiRpcbJi4SDSr8ppV8vOUHh4USNbZwKTSUm+rdXOkVFShQSQnTk\nSFO5JIkKPyykhNAEKv+pvNXzbr6Z6OmnW+9PEiWyNljJVGyi1KtT6cSDJ1z4itznlhlGeiCmyO7f\n5g2pp+cvOlU3mmLOTses3dZyDABjxgs4WqpBY/ZZWnaqHZXFIgJ93d+14uEtwCI63rSWnAy8/jqw\nYgXg7w9ccgnw8svAujfr8LvvXrz3hQJLHhNcklWiWXiQBLOkwMiRrtmfIMgpZhYuBKbP90Tjbf1w\n4p4TDqUGc8RXS23wqTTgtlURUHrKTRGCQsCg7wdBv6UKn19bjuRkOcWN1Sp3nQ/qJ8Fw1IC0ZBGL\nHxHaXYrZnoAAYMGNwBuqoXjwcUWHC6V4D/VGbIANg+NE/P57516nMcOIYsmxNG5nUqgV8DvfD3FP\nxmHE+hG4t244jiQD898KcOj5U6cCx03eKPm77Z6egiwREf7OtUBEhBFKXJACDADefBNQFBlx1zwz\ngmY4n2bLUX7n+2FeZBW+/Fjs1Jg+U74JSROTcOzGY/Ae4Y3zs87HoK8Hwe88eVxkeTmg18OhjCT2\nKJXA8pUqfF4fh9TPW86Q06fqkTQxCao54VjdayCeT4+DVQQK32/dm1OYS4iKcvy4MdGE/KzOn0P3\n7BEw4Xzn39DASwIhKAXcNV7umWlrIm5dnTxhd+ZMYPlyoN7B6Q+SWULBGwVIjInDglu6fsnsPViJ\nojLn9vPayxJqDxvwZmgGHvPKwuxZhN27HXtu0UdFyP13brvjeOv31EPUi/jpgHerfOf+/sDdt4p4\n/zttl7rjS5eVotbLA5/96Y3/+z/724SGAhdMEXAwOhpVv7dMP1GVa8H+Qg/c+KZ/q+d5D/JG8Hk+\nWDBOfzKt2+bNwPggPcLntH8umH6jB5LzNe2mlMxakoU9UXuQ80yOyycn6n7TYbU1GvfdJ2ezskcb\nqcU1l1jxzGV1eO89uefdWm3F0auPouy7MozZMwZh14e1et7bb8vjks/MBS4oBKh8VNBGaTFw2UBU\nra9C5ZqenW/VaAR+26bGnU/b76W9eYkG63ZpIBpc1ArelcjamX+w03L89ttEiy5qoN2xu8lw3HUt\niZ3xyZximjGs/btHV9j2VQMN9tY7tK1eTzRwINGKFS0fb8xrpMTIRNL94XirjjMWXmulfr263kJi\nz4oVRKGhEn02IINKlpd0eX96PVGYl4V+vt5+q4Y+TU8JIQmUu6WBBg4kiokhuvhCkb4PSabqLdVd\nOnZGBlFYGFGVg42UOS/k0P/NLKM5czp3vKQLk2hYX6vLW40cdcEQE308yf77vG9lA0UqG+ndRW23\nLNvzy00FNCSq862NzQ4eJAr2E+m3+CSy6TtoxneB4uUl1MvbRImJzj3PcMJAu3vtpoJ3C0gS7bc8\n//470aWXdr2MSxYY6SLfqpPHqfqrilYFHaTbLjVQQADRHXcQzZhB9MBtVkoISaCGlJat/nP71dL7\nCx3/jnx3TRGdH9+5XgCLzkJzVKW09OPOnXdKvymlw5cepnvuIXrmmTP2bSFaupQoPFx+zatWEc2b\nR+TvTzR/PtGGDXKvw8ntayyk26Cj7GezKWlKEu3w3kF7xxykkCCJcnM7VbwWdq02Un+NY9cAIqJt\nWyUK8bDQ1jnpJIkSZS7JpI8HnaDQEIm2bm37eZIkUdZTWbRv8D5KmpJEqVenks1g/7tx9MajtP3p\nEgoPJ7u9aDodkb/SSik/ONfi3axqcxUlhCbQwqss9Pjj7W+7ciXRtGEmSrm8Za/yp/MraHJ02z1T\nVX9V0boBSRQYKFFdHdHNN0n0uHcmNRa0XydtRhvFK/SUuM1+r0XJ8hLa228vNaQ0UPod6ZQQmkAF\n7xaQrbHr5xmb0UbrfHZTgL9E5a0bfluoXFdJhyYeIiKi2t21tLvXbsp8JJNEU/vfmbffJpo9u2XP\n6Znq9tVRQmgCGXNaxkDtPeds+8+rjTRRW223N4yIyGQi8lNb6dBbpUTU9Zbjbg2ON2+WhweULCuh\nxKhEakh1rkvWlV4+r4jmT3V/cLx7pZ76ejh2I3DnnUQLF7Z8zFpvpf3D91PBewVuKJ3sueeIbrjB\nbbunbduIQgJFesE3g8xlXQuMnl1ooIu9dGRtaLs7tuLnCtrdazflHDHTxg0SJU1Lppx/53TpuM2s\nTvQCN6Q20LbovRQQIFGR/Z6hNpkrzLQrOIF8vCWq6dz1qcteecJM13iXtnr852erKEAw08cPOzdc\niIho/7OFFOLV8k20WokOHCB65x2iJ54gqu4gPjMaiQb1E+kFn+NUn+R8GTrDZrTRvd65tPBaxyuA\nPk1PidGJVPxlcbvbvfAC0bPPdrWEREajRHHaRvr6yVpa90wlTdVWUrC/SM89R1Ta9DFWVxP17Uv0\n8aIa2j9yP4nmUxeecYH1tOb1ujb23trOf5VSjF/nvs+6P3TU26uRDh/u1NNJNIuUGJlIaRv0FBQk\n37BKknyjMXAg0fTp1GrfOh3R0qUSjRslUliAjW4bWU3v9z5OK7z20/6pyZTzfA5Vbaoia72V/vxT\nHtrmCqUFIvnBQqK14xuBsjKiCF8rLR184tSQMUmijDsz6MvRJyg0VKL161s/T7SKlL4onbaPTaIV\nX1lpX6JIx245RgfPO0imkpbjusxlZtoVsIsee9hGTzzRdlnunV5PNw1wvkGmZlcNJYQk0I6v6ik8\nXB6m0R6DgSggQKJffPaQuUKuT6JZpEs9K+iD59oOdCVJon1D99HcKWb66CP5GvNLv2SHynhDTCW9\neFvrG5a6vXWUEJJA+qOn/qZP01PKFSm0O243lX5bandIkqMqf6+k+3sV0x13dLytaJXr+ImHT1BC\nWAJV/lrp0DHMZqLx44luu03+uS0F7xXQwfEHSTSL1NBAdO21RBMnktPXKncZH6Gnzxa0HlJzuusv\nMdFjkXkkSdK5HRzrdER+fkSiSFS2oowSwhPO2sXtTI/EF9ED17lm7GN7Dq3TU6ym4+OsXi1ftOpP\nezskm0Qpc1Mo484Mktx4S5eWRrR3r9t2T0REKSlEkb4WenRkWafvTktzrOSvsNCBbzqOFrOfzqbk\ni5Ip+zn5/66c0DpLkiTaN3gf3Xalie66q/27ckmUqG5fHeW+mEsHxx+knX47aef1xyk4uPtu5ZOT\nJYpRGKmxsPFkGd++rJwCFRZa/1XnbiwLl5WQSpBo506i118nmjVLPicMHUp0331Ed91FFBdHtHNn\n2/t4ZLFElwRVUf677rthtGfPA7nkp7W1+I62pf5QPSVGJFLp92W0fLn8mhYtItq3r3U9mD2b6Jdf\nXFPG316pITVEilI10rv/MpPeToNlSgpRSIhE/73wBGU/k33y8Ti1kQ7+4fjnWrhGR2pBJLETjb8H\nl+SSj1bscOx+e/Jey6Njtx6j228nuv12oosvJho8mOiPP069x5IoUf2heir8qJDSbkijxOhESghP\noHWXHqfFl9bSxedbKT5eIq2WKD5eblm//36iyZOJPvig82U7nSQRaSCSLrX991YUiS4c3Ei3Bhad\nDBJP7sMmUdoNafT9BZkUFibRzz+f+lthjo2eH1VEk4PryNdXopkziUJDiZKSJMp9JZd2x+2mhiOn\nGqJyX8ml1DsyKCKCTo7XtafgqIV8BQsVHnf8hrDugNwimfdrFY0cSfSf/zj2vFtvJXpqdCkVflwo\nH/vbMvJVWamkg87G4i+Ladn5WRQYSBQfYqHMRx0buP/JjZV0cZ+WjXOmEhMlRidS2ZpKevllov79\niaZOJbrlFrkB6f3HDfTxoExaNSCFyjd0bo5Dyu3HKczPRikOTr3KfTmXki5IosY852IVvZ5o7lz5\nO9FW44okSZRyRQr9vSiXRoyQg+lXXyWKiiKXzT/qrGP7zBQgmKmhuP2b73XrJBrlWUc122vO7eCY\nSL4I/vGH/HPFzxWUEJZAtXuc65p1RHY20fLl1OaFbGFwMb3wsPsn5B3bZqQIVfvHyc+XT2b797d8\nPOvJLEqeltyidedclnfCRvFqA00bYKTXFtbR9rcrqWxVOZWvLqeKNRVU8UsFNRxuaLOFZeHwarpx\niGPdvpJNosMzDlNCWEKrlpOzKeeFHDp0XzaNGUP01FMtAyNrg5XKV5XTsVuOUUJoAu0bso+yHs+i\n6m3VcstYItH553db0UkUiYK1Ftq3tJIsdVa6Z0AZRXuY6HBC51v/dRt01NfHSKNHEy1eTLR2rTwp\nh4jIVGqishVl9NsaG0VEyBNRz2yp37aNKNzHSjsuTmtzmIK7NOY10gVqHd2zSGw3QK7dXUsJYQmU\n/UMlLVggn/O2byd64w05ABsxQu72r6mR60NICFFhoWvKKIkSbXqwiIyl7X9GK1cS9Y6T6M+wfVSb\nWEuSTSJPWKmqzPFoVX9MT3FqI02fTvT559RhN/HpPhqeTVNHda0XyVJloV0Bu+hogokGDiT69NOW\n9UV/VE+HJh2ivf33UsZdGVT6TSkZMg12GxrMZqITJ+SJqB9+SPTYY3JjjqvEeppo75ft39Q/u9BA\nI9W1VHfMfk+jaBbpyJwj9PPsbIqIkOihh4jOHyeRr8pKl8XV0KoVp+rlf/8rDykrKCAqW1lGCaEJ\npFuvI9Eq0u6Y3bTyPSNNntxxuRfE6+i+mY718DakNFBCeAKVrq2kOXPknlBHG0I2byYa2c9KB8cf\nJEmSaOmA43TewI4nxNmMNtoVkkBD+ou0ILLc4Ym5JzY2kJ/CevLGTjSJdGjCITrwRD5deqncw33o\nENHWrURff0300kvyze2MGRL1CrfRWI9a2jMzhQwZjg8PlUSJnvM7ThdPdv8wMCJ5UufixfINY04b\nHafb/rBSsMJMryxqOPlZbdggxyOffNJ9wyweuLCGbh3W8WdpMskT5bfOST/3g+MNG+Q7suYZ/Lo/\ndJQQkkA1O7red3zsGNErrxCNGkUUFibR+BE2unSqSGZz6094rraEPn7b+ZnSzsrZ30iBirYvAlYr\n0QUXEL355qnHDJkGynwsk/b03UMWnXtmzHaXom119NbkIro6vppivM0UrLXQrJgaen5kEa2blkF7\nBu6jHd476NDkQ5S5JJPKVpaRMctIST/Ukp/CQkUnHP/MrA3WVmOqzraGVHmMfWWFRMOHE73wtI3K\nfiij1HmptNN3Jx2eeZiKPi0iY27rcn77LdGNN3ZDoU9z9agGenpUCc32r6QRoUYqKejajVr9wXo6\nMOpAi8dsjTbKez2PdgXvoqQLkigxOpGS3i6lSy+VaOLEUyf22lqimHCR3gk82uXhOZ21dWY6XRZR\nRaG+Nnr3cSOZGlq+H9V/V1NCaAL99UENxccT3XuvPAykmSgSbdlCdP318hjYa68liojonovQo48S\nXTTKTAnxe6hkr548BZtT5bAZbbRRvZNWfm6mG64Tyd+faNo0OfBvr8VPtIp0izqfnn2860HC8QeO\nU/az2S0eE00i5byQQwkhCVT0SdFZv4myZ3yUgZbNyqeG1Aa75dnwpYGCBDMd+6X9YS02o42SpibR\nhhty6dEHbPRx/+OUcmeG3Z6xd96Rb8xqauQbtsSIREq7IY0OTTpEV1whB30dSV5WRX4qa4c3Cobj\nBkqMSqTSH8vogQfkMfTOZASy2YgiIyX6IfgQlX5TSlf5ldGbbzj2ueX8K4fWzsqiNd57HR4XLNnk\nXrGDm80kSRKlL0qnb6ZkU0yMRM880/7wOZuN6LaFEk3o20h/BSVS5iOZZKnu+MXWJNbSAA+93WEx\n7vTRR0SRka17h//zH3kOzdr39JQQltBirHZmplx3br+97YwaokkkY7aR6vbV2c000lkWo0ihChPt\nW+PYjcdN80Va7JnV5eBYIDo7yaAFQaC2jnXFFcCUKcCTT8q/12ytwbH5x9Dvw35Q+atgqbDAUm6B\ntdx68meFRgHfcb7wO98PvuN9oQnRoKEBSE0F/vwTWLMGqK8hzBxixBSxAn0OF0MdqcETOf3gJ1jx\n+sRS+Iz0hs9IH3gP9cbcCWbc92MIbpjvmpnzbSnLsKD/EAXqbCq7OXVffhnYuRPYsNoG3c8VKPu2\nDI2ZjQhbEIbYJbHw6NXxalXnsrw84O+/5X/btwNlZUBIMCHER0SAYEWA2QSf6kYc03th5o1qvLWi\nc8u4dhciwoGhBxB6XSjy95hw69ZeuHpAHZ56WkDw3GCoA9tOm/HCC3Lmj5deOosFPsMXzxlw/+te\nmDHChJ93e8Dbu2vfF1ORCUnnJ2FS8SQQESr/W4nsp7LhO8YX8W/Hw6ufF+r31SPn2Rw05pvx18Qh\n+GSTDz78UMCG3yQY/qjAlz+rETwr2EWv0DnWWiuqN1Rjz28mvPFnIEr1KizuX4q5cwgecVrkvV6A\nv68ZjaVrvfDZZ22vvgfIWRa++07++bHHzk75T2ezAZdeCgys0uESWxkey+2H/Ebnzjepc1PRkNQA\na6UVFqUSST4h2CGGIqHeHwODzbh+tB5zJ1vgF6OCJlQDdZgaVp0VM65V47U1fpg1q2uvwZhlRPLE\nZEzImwCltxK1u2px/K7j8Brkhf5L+59cOKC73XSlFUN0lbikohBWnRV+E/zgN9kP/pP8UeehxfgL\nlfjoqUbMf63jbDK2ehuOXHIEjZmNiH4wGr1f7n1ygYzTEQEPPwwcOyZfI8WiRqRdnQaPB+Mx9clg\nFBbKeY3bI1klzPOvwJDbg/HmJ/bPVY15jTh84WH0frE3fqqPxLJlQEKCnPXCGY8/Dhh3VuO6I6mY\n73EBEg4q0b9/x88zl5qxt9deBF4SiBEbRjh8vGv61OC8KUrcPK4Bb70q4WfE4OtvBMye3fFzJQm4\n804g+7iEDwZmwbi+Us63f3ckFCr7mUl+vLEEz28KQWalxqH8+q70++/AHXfIiwbNnQssWSJn91i3\nDhgwACh4qwC6dTqM2j4KCrUColFE1XEz7npEhYJC4LNrKxDUYIS5yCz/KzbDVmuDJlIDVaAKplwT\nAqcHImReCIIvC4Y62H5dISJYyizQH9bDUmIBFHImDSjl/wWlgPVrRby70Q9pdY5d6//4A3jpnkYc\nLPECEXX6AtUjguOsLGDCBDmxdXPqoNpdtch5OgcqPxXU4WpowjTQhMsnU1WoBukngKRtVqQeJqQX\nq5ELb9SSBv0ibTg/0oiJ1SXoV12N4BmBCJ4djMAZgdBGaqHXE6ZPJUzqb8aSsZUwpBigP6LH/bkD\n8cYvfrjkEve+D/pyK+IiCLWCGp5awMuT4O3Z9L8XUFAk4KcLs6DcXo7A6YGIuDUCQbODoFB3a9a9\nbmOxyGupV1TIKa6af9bXSXjiaQV8fLq7hM4r+aIE1RurEXptKExjgjH9ChUefli+cLXnppvkNFQL\nF56dctpTX0/4bqkV9z2l6TB9nSMki4RdPrswetdoZD2aBalRQt/3+tpdXKd6SzVyn81FRp0HXjQO\ngqXGht9vL8KIjzuZ88zFiIA/14l44lGCt2jDorhyrLZEoVGpxo8/Ar16dXcJO1ZRAZw3lnB+bRkK\nvP2wr6xzN59EBLFBhFVnhbXSioYSCzZuVWLFVg8k52kwu1cd5gXrEG+uR2OFDbPLxqOoQoEAxzIL\ntivt6jT4jvWFqcCEqvVV6P9Rf4Re3cHqLmfZ//0f8NRTcuourYagUUjQSBJUFhH6RgHXX2zGJ1vt\nL3lsj7XairqEOoTMDWl3O1GUb9D8/IBvv5Vvtt94A8jNBf7zH8eO9feiPNzwUwwunwv8+2ELgrQi\nRIMIyShBbBCR/WQ2Yh6JwcHYGDzwQNMSznEOv5STDh8G5l0m4Sl9KpZGjcDRdMfjnMyHM+E71hcR\nt0Y4/JwPb6/BqjUKKMwixKH+WP2bErGxjpf3ZICcDfz0ph5lz2XBUm6B/yR/KLwUUHgqoPRSyj97\nKLHgKU9c/aAnlrzlxLK8LpSUJAfGQUFAdDSwciVOfv9IIqRelgpjhhG2OhtEowhtjBaaaC1WNETj\nh8wgXDjEArNKCYtCATOUMIsCjEYBVivQv4+EwR4G9CnTITqlDNHneSLkyhD4TfRDY2Yj9If10B/R\nQ39YD5IIPqN84BHnAZC8OBdJBIiA3izg9r974/4HgQdecex9MpuByHBCTZ3i3A+OAXnVt8JCefWX\n9kgScOutcuvqeecBw4cDw4YS+nk3IqSsDsaDDVCHqhE0Owh+4/wgKFu/NzqdvNLcQw8BDz4oPzZ6\nNLBsGVyaL9hu+W0Sjl1/DOYqKxrNAozN/ywCGq0KRIZKGHl7EMLmh7V5t8X+t+Tny3mEn31WXq61\nLeefD7z/PjBp0tkr29mQEJQAhYcCfV7rg4iFEXa/s82ICLpfdEh/Nh9mtQrTD46AQtuzbhxFUT6P\nvfUWcPXVcku/quMVmHuM/fuBKRcQ5l1sxU8b21/5sTMKC+U8w8uWAeHhcmv1unWuW/mtLrEOyVOS\nEXl3JOLfjIc6oGeeR0VRvpCbzYDJdOqf1UoYNkxwW2ui0QhcdJF8o/3SS/KyzT/+KC+f7ghDugEJ\n09OwvD4Wf5lCcWdECa6Lq4LWWwmFtwJBs4JQMjYac+bILdTnnde5chIBw4YBgQGEaRcJePXVzu3H\nURk7TRgyVYv7rzXh/R89nc59D8jxyV13yXmF168nWPfVwJRrkm8cGkXoawm/JnlhRbIfbCQgOU8D\nbx/39la3p7BQ/u7de2/r1QrFRhHGDCO0MVqoQ9QteiP27AGOHwe8vABPT/lf889Kpdw7cfAgcOAA\nkJxMCPaWMNjTgNGWKiyYZITvKB/4jPKBz0gfaKI0dns6zGZ5eew+feQWbmdWXl24EPj+e6FLwXG3\njzlu1tBAFB1NlJDQ9jaSJA8onzxZTvfSFTk58izM5lm+cXHkkhyWjHVGZqY8Yea77+TfzWY53dbR\no/JM4d9+IwoIcG6S07miNrG23VR89kg2ySV5Rpl9P/0k/3Mnm02ec3LVVUQvv+zafXfXGPRzRXm5\nPBl00SJ55dLOjnFPTSW66CKi4cOJduyQH8vLk6+tv/7a9XK+/ro8M+ps5XbPTel6vRFF+X2dMkWO\na4jknPiLFxMFBclZIzZupE5ldTkX2Wzy/K9vv5UnlC9Y0Pa45dOfc8MNRFdfTZ3KYPP77/S/Mea4\n2Y8/Au++K7dc2OuyfeMNeZudO4HA1r2uTktOlu+e16wB5swBiovl7ibGukNGBnDxxUBDg9yCFBgI\nBAfL3V5BQfId9IcfOncHzRhj9pw4AUycKM9lWLy48/shAlavlscIX3CBPDzyzjuBRx7pehkLC4H7\n75dbN8+l854kAffcI/eG+PjI78miRfJj58LwKndpbJR7/ktKgF9/BULsjAIiknv1jx6Vex7aWjWw\nPWYz4OHRtZbjHhUcEwEXXgjcckvr7uWvvgJeew1ITIRTS5p2ZPNm4OabgepqeXzrufQFZP97jEZ5\nYpSvL9dFxph7lZXJAYorhv0YDPIS7oA8sfyffv6SJOC994CICOC66wCttrtL1DNIEvCvf8k3VOvX\nyxMAT/fKK8DatfKEfGcncZ5OEP6HgmNAHoQ/c6bcitbcOvzrr8B99wE7drR+I11hxQrgxRdbrz/O\nGGOMMcZca9kyeZ7N6tVyoygAfP458M47cnaTCMfnUtr1PxccA3IgrFYDH30kB8TXXgts3AiMHeu+\n8omi/aEcjDHGGGPMtbZsAW68UR5O6+kpD+/ZuRPo64IERP+TwXFVFTB4sNwl8eijcoqR6dPdXEDG\nGGOMMXbWHD0qZ6VoaJCD5VGjXLNftwfHgiDMAvABACWAr4jorTP+Pg3AbwBymh5aQ0Stkq44ExwD\nwKefymnWfvpJHq/DGGOMMcb+t1RWyil2Bw923T7dGhwLgqAEcBzAJQCKARwAsICI0k/bZhqAR4lo\nbgcFdSo4FkX5jmKE4wvcMMYYY4yxf7iuBscdpRkfDyCLiPKIyApgFYAr7ZWjswVoi1LJgTFjjDHG\nGDu7OgqOowEUnvZ7UdNjpyMAkwRBOCIIwgZBEIa4soCMMcYYY4ydLR1lN3RkHEQSgFgiMgqCMBvA\nrwDckHCNMcYYY4wx9+ooOC4GEHva77GQW49PIqKG037+UxCETwVBCCKi6jN39uKLL578edq0aZg2\nbVoniswYY4wxxphs+/bt2L59u8v219GEPBXkCXnTAZQA2I/WE/LCAVQQEQmCMB7Af4mot519OTUh\njzHGGGOMMWd1dUJeuy3HRGQTBOFBAJsgp3JbRkTpgiDc0/T3LwBcC+A+QRBsAIwA5ne2MIwxxhhj\njHWnHrkICGOMMcYYY53h7lRujDHGGGOM/WNwcMwYY4wxxlgTDo4ZY4wxxhhrwsExY4wxxhhjTTg4\nZowxxhhjrAkHx4wxxhhjjDXh4JgxxhhjjLEmHBwzxhhjjDHWhINjxhhjjDHGmnBwzBhjjDHGWBMO\njhljjDHGGGvCwTFjjDHGGGNNODhmjDHGGGOsCQfHjDHGGGOMNeHgmDHGGGOMsSYcHDPGGGOMMdaE\ng2PGGGOMMcaacHDMGGOMMcZYEw6OGWOMMcYYa8LBMWOMMcYYY004OGaMMcYYY6wJB8eMMcYYY4w1\n4eCYMcYYY4yxJhwcM8YYY4wx1oSDY8YYY4wxxppwcMwYY4wxxlgTDo4ZY4wxxhhrwsExY4wxxhhj\nTTg4ZowxxhhjrAkHx4wxxhhjjDXh4JgxxhhjjLEmHBwzxhhjjDHWhINjxhhjjDHGmnBwzBhjjDHG\nWBMOjhljjDHGGGvCwTFjjDHGGGNNODhmjDHGGGOsCQfHjDHGGGOMNeHgmDHGGGOMsSbnRKeXAAAN\nZklEQVQcHDPGGGOMMdaEg2PGGGOMMcaacHDMGGOMMcZYEw6OGWOMMcYYa8LBMWOMMcYYY004OGaM\nMcYYY6wJB8eMMcYYY4w14eCYMcYYY4yxJh0Gx4IgzBIEIUMQhExBEJ5qY5uPmv5+RBCE0a4vJmOM\nMcYYY+7XbnAsCIISwFIAswAMAbBAEITBZ2wzB0A/IuoP4G4An7mprOwctn379u4uAuvBuH6wtnDd\nYG3husHcpaOW4/EAsogoj4isAFYBuPKMbeYC+BYAiGgfgABBEMJdXlJ2TuOTGGsP1w/WFq4brC1c\nN5i7dBQcRwMoPO33oqbHOtomputFY4wxxhhj7OzqKDgmB/cjdPJ5jDHGGGOM9RgCUdtxrCAIEwC8\nSESzmn5/BoBERG+dts3nALYT0aqm3zMATCWi8jP2xQEzY4wxxhhzOyI6s+HWYaoO/n4QQH9BEHoD\nKAFwA4AFZ2yzDsCDAFY1BdO1ZwbGXS0kY4wxxhhjZ0O7wTER2QRBeBDAJgBKAMuIKF0QhHua/v4F\nEW0QBGGOIAhZAAwAbnd7qRljjDHGGHODdodVMMYYY4wx9k/CK+QxlxMEIVYQhL8FQTgqCEKaIAgP\nNz0eJAjCZkEQTgiC8JcgCAHdXVbWPQRBUAqCkCwIwu9Nv3PdYBAEIUAQhJ8FQUgXBOGYIAjnc91g\ngDznqemakioIwo+CIGi5bvxzCYKwXBCEckEQUk97rM360FR/MpsWtZvR0f45OGbuYAWwhIiGApgA\n4IGmxWOeBrCZiAYA2Nr0O/tnWgzgGE5ltuG6wQDgQwAbiGgwgBEAMsB14x+vad7TXQDGENFwyMM8\n54Prxj/Z15AXqDud3fogCMIQyHPmhjQ951NBENqNfzk4Zi5HRGVEdLjpZz2AdMj5sE8uGNP0/7zu\nKSHrToIgxACYA+ArnEoDyXXjH04QBH8AU4hoOSDPeSGiOnDdYEA95EYXL0EQVAC8ICcJ4LrxD0VE\nuwDUnPFwW/XhSgArichKRHkAsiAvctcmDo6ZWzXd8Y8GsA9A+GmZTMoB8EqK/0zvA3gCgHTaY1w3\nWB8AlYIgfC0IQpIgCF8KguANrhv/eERUDeBdAAWQg+JaItoMrhuspbbqQxTkBeqa2VvQrgUOjpnb\nCILgA2ANgMVE1HD630ieCcqzQf9hBEG4HEAFESWj9eJBALhu/IOpAIwB8CkRjYGc/ahFNznXjX8m\nQRD6AngEQG/IgY6PIAg3n74N1w12OgfqQ7t1hYNj5haCIKghB8bfE9GvTQ+XC4IQ0fT3SAAV3VU+\n1m0mAZgrCEIugJUALhYE4Xtw3WBya04RER1o+v1nyMFyGdeNf7zzAOwmoioisgFYC2AiuG6wltq6\njhQDiP3/9u492MqqjOP49weCouIFNbGEIFInGnTKMtRMU3OaUnLygoIMWNnNwhrzgrdpxmbKy6jF\njDoaChOKoglCM14AQywUlNtRCG+BmYGWRV4wRc7TH+vZnZfT3ocjFw/C7zPzDmuvd79rrb32q+fZ\nz17vfivP2zfrGnJwbJucJAFjgCURcV1l1xRgeJaHA5NbH2tbt4i4KCJ6RURfygU1D0XEMHxubPMi\nYiXwoqT9s+pYYDEwFZ8b27qlwEBJ3fLvy7GUC3p9blhVo78jU4DTJHWV1BfYD5jbVkP+nWPb5CR9\nHpgFNNHy1cUoysk4EegNLAdOjYhVHTFG63iSjgTOjYhBknrgc2ObJ+kgyoWaXYHnKTeV6ozPjW2e\npPMpAU8zMB/4FtAdnxvbJEkTgCOBPSnriy8D7qXB+SDpIuAbwLuUpZ4PtNm+g2MzMzMzs8LLKszM\nzMzMkoNjMzMzM7Pk4NjMzMzMLDk4NjMzMzNLDo7NzMzMzJKDYzMzMzOz5ODYzP5H0lpJCyrb+Zuw\n7T6SntyI4wdIuiXL51XG+KSkdyXtJqmXpN9LWizpKUkjK8f3kDRN0jOSHpS0W2XfKEnPSloq6bgG\n/c/M/bV+J27oa9kYksZKOmkj29hR0m2SmnL+HpG006YaY53+jpL075y3JZIu20Rt1bajN+FYR0ga\n3WDfDEndN1VfZrZl2q6jB2BmW5TVEfGpjh5EA+cBowEi4irgKgBJxwM/iohVknYAfhwRCyXtDMyT\n9GBELAUuBKZFxJWSLsjHF0rqDwwG+gMfAaZL2j8imlv1H8CQiJj/PrzWtgQtN9fZUOcAKyJiKICk\n/YA1Gzuw9ZgVESdI2hFYKGlqRCxY30GSOkfE2lbVD0fEoM0zzDbn9g7gLOCazdS3mW0BnDk2s/WS\ntFzSFZlpnCOpX9b3kfSQpEWSpkvqlfV7S5okaWFuA7OpzpJuyqzuAxnMImlkZnsX5Z2PWve/PTAw\nIh6vM7whwAQotyCOiIVZfgP4EyXgBRgEjMvyOODELH8NmBARayJiOfAccEijqagztsmShmX5O5LG\nZ/ksSXPz9d8tqVvWj5V0vaRHJT2fmdBxmVG9tdLuG5KuybmaLmnP1uOQdHBmtJ+QdL+knu2ZT6An\n8Lfag4h4NiLeyWPPyPd4gaQbJXXK+uMkzZY0T9LEWqY5z42fZn2TpAMazF2tr9XAPKCfpH6S7svx\nz6odm3N0o6THgCva+T70ycz++JzLuypzfoyk+Tm+MZK6Zv1nJf0x36PH8gMVwIdzXM9IqvY/hXLb\nczPbmkWEN2/evBERUG6tuaCynZL1y4BRWR4GTM3yVGBYls8EJmX5TmBkljsBuwB9KNnJAyvPGZrl\nl4AuWd6lzrgG1vpsVb8j8CqwW519fYAXgJ3z8b8q+1R7TMlGD63s+zVwUp32ZgJLK3NzRdZ/CHgW\nOAJ4ujYWoEfl2MuBH2T5VuD2LA8CXgM+mWN6ojI/zcDpWb4UGF05/utAF2A2sEfWDwbGtHM+D6Lc\ncnV2ju3jWf8JSgDYOR9fn+/3nsDDQLesvwC4tHJunJ3l7wE31+nvqMo5s0ce0x+YUen7c8CMLI/N\ncahBW6tY9zztm+93M3BoPm8McC6wA/CXSj/jKJnz2i2qD876nSm3qh6R9d2B7Sm3od230v+fgZ06\n+r9Vb968bb7NyyrMrOqtaLysopaBvAO4NssDacnAjgeuzPIXgTMAoixPeE1SD2BZRDTlc+ZRAhqA\nJuB2SZOByXX6/iiwok79CcAfImJVtTIzgHcD50TJIK8jIkJSW1+f19tXd1lFRLySa2gfAk6sjGWA\npJ8Bu1ICr/srh03Nf58CVkbE4hz3YsqcNFECvTvzeeOBe6ovETiAElRPlwQlsKtlg9ucz4hYJOlj\nwHHAscDjkg4FjgEOBp7INncAVlIC1/7A7KzvSgmsa2pjm08J3Os5QtL8fF0/p3xwOQy4K9ustQtl\nru+KiEbv0SMRcUK1QlIf4MWIeDSrxgMjgWmU8+65rB8HnE0JzFdExLyckzeynaAE6a/n4yWU8++v\nefzLQC/KByUz2wo5ODazDVENWv7vK+426t+ulNcC3bL8VeALlGD3YkkDYt11ptGgvdNoCdpLp1IX\n4LfA+IioBoYvS+oZESsl7QO8kvUvUYKdmn2z7r04EPgHLUs4oGQ/B0XEk5KGUzKeNe/kv82sOyfN\n1P//sqgfsC+OiMPq1K9vPomIN4FJwCRJzcBXclzjIuKidTov67qnRcSQOn1ReQ1rG4wfWgW0knah\nZO8bfRhb3aC+La3Py3pz1uh8rWp9nnZuR7tmtpXwmmMza6/BlX9rWcPZtKzBHArMyvIMylfsSOqc\ngVBdKmnD3hExk3KR3K5A619OeIGyTrZ63K6UAPDeVm2NAZZExHWt2pgCDM/ycFoyqlOA0yR1ldQX\n2A+Y22i4dcZ/CPBl4NPATzKDCSVbvDKD9TN47wFVJ+CULA8BHqnsC8oSjr1q67kldZHUvz3zKekw\nSbtnuSslK7yc8r6dLGmv3NdDUm/gMeBwtaw130nlIr4NFhGvAcsknZxtStKBG9Mm0Fst69trc/Y0\n0Kc2dsoykZlZv4+kz2T/3SV1pn7wXK3bm5YssplthZw5NrOqbpKqvyBwXyWLuLukRcB/gNOz7ofA\nrZLOo2Riz8z6c4CbJH2Tknn7LuXr6NYBYlCycr/JYFfALzNwqlpEWUZQdSLwQES8Vak7nBKINlVe\nx6iIuB/4BTAxx7QcOBUgIpao/CzbEsqa6++38XX+bZJq/f0dOB64CRgRESsknQvcAhxNWSc8J583\nhxIsV193vXLVm8Ahki6hzN3g6s6IWJOB5a9y7rajLHd5hvXPZz/ghgykOwG/i4h7ALK/B/NCvDU5\nH3MljQAmqFwcCXAxZa31OsNq8Hoa1Q/NcVxCWUM9gbIkpK15CcoSjep5ejllmc7TwNkqP/m3GLgh\nIt6WdCZl+cZ2lA8+N+b8DQZG54V7q4EvNRhr5Nz0BF7NrLuZbaXU+G+AmVkhaRnlwqV/duAYxlKC\nnTkdNYb3k6TXI8K/qdtOmbGfGhEDNmMf36ZcjHftep9sZh9YXlZhZu2xJXyKvpqSgd5WbAlz/kGz\nuedsMHDzZu7DzDqYM8dmZmZmZsmZYzMzMzOz5ODYzMzMzCw5ODYzMzMzSw6OzczMzMySg2MzMzMz\ns+Tg2MzMzMws/RelflGQ91Z04QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x11b337b90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def plot_losses(train_losses, val_losses, epoch_examples, epochs, min_y, max_y):\n",
    "    x1 = [a/epoch_examples for (a,b)  in train_losses]\n",
    "    y1 = [b for (a,b)  in train_losses]\n",
    "    x2 = [a/epoch_examples for (a,b)  in val_losses]\n",
    "    y2 = [b for (a,b)  in val_losses]\n",
    "    plt.figure(figsize=(12, 6))\n",
    "    plt.axis((1, epochs, min_y, max_y))\n",
    "    plt.plot(x1, y1, \"-m\", x2, y2, \"-b\")\n",
    "    plt.legend(('Train Loss', 'Validation Loss'))\n",
    "    xlabel = 'Epochs (%d Examples Seen Per Epoch)' %(epoch_examples)\n",
    "    plt.xlabel(xlabel)\n",
    "    plt.show()\n",
    "\n",
    "plot_losses(train_losses, val_losses, 7200, 100, 0, 3.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Code for Reloading the Saved Best Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Building model from ./models/GRUs/GRU-20160425154620.npz with hidden_dim=128 x_dim=1\n"
     ]
    }
   ],
   "source": [
    "model = load_model_parameters_theano(\"./models/GRUs/GRU-20160425154620.npz\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}