demos/rnn/lstm_foundations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/Humboldt-WI/adams/blob/master/demos/rnn/lstm_foundations.ipynb) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Long Short-Term Memory (LSTM)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having looked into the computations within a RNN cell using numpy,  the previous demo equipped us with an understanding of recurrent computations. It turns out that RNNs have trouble learning long-term dependencies due to the vanishing gradient problem. Long short-term memory cells are designed to mitigate this problem and, together with GRUs, represent a state-of-the-art approach for sequential data. \n",
    "\n",
    "This demo splits into two parts. First, we revisit the LSTM cell architecture and the corresponding calculations. Using plain numpy for that purpose should give us a good idea of how things work. Thereafter, to prepare for building proper LSTMs in Keras, which we do in a [subsequent demo](https://colab.research.google.com/github/Humboldt-WI/adams/blob/master/demos/rnn/lstm_fin_forecasting.ipynb)), we discuss some peculiarities about LSTMs and their implementation in Keras. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LSTM Cell\n",
    "The architecture of LSTM nodes can be daunting at first, so we will look at the pieces one by one. Here is a well-known picture from our favorite textbook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"lstm_3.svg\" />\n",
    "\n",
    "Source: Dive into Deep Learning http://d2l.ai/chapter_recurrent-modern/lstm.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "See the line that passes from state to state at the top with a plus in the middle? That's the hidden state and we see that the input information $x_t$ gets added to it after being processed in a more complex way. Imagine the hidden state is our long-term memory that stores information that may not be immediatly useful. \n",
    "\n",
    "Before information is added to the hidden state, the hidden state gets multiplied by the output of a sigmoid activated layer. We call this layer the **forget gate**, because it reduces some values of the hidden state by multiplying them with a value between 0 and 1.\n",
    "\n",
    "The input that is added to the hidden state is processed in two steps, which make up the **input gate**. The tanh activated layer processes and transforms the information. Almost as in our RNN, although the input processing in our RNN did not use an activation after the input layer. The processed input is multiplied with another sigmoid activated layer, just like in the forget gate. And just like in the forget gate, the sigmoid layer here determines which parts of the processed input will be added to the hidden state and which will be forgotten.     \n",
    "\n",
    "To the right, we see that the hidden state is passed on in two ways. On the upper level, it gets passed on as the hidden state. On the lower level, it is put through a tanh activation and filtered by yet another sigmoid layer. The result of this **output gate** passed on as the output at time step $t$ and passed on to be combined with the input of the next time step $t+1$. Image the output state that is passed on as our short-term memory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step-by-step LSTM"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our LSTM is to process input data of the form $ X_t \\in \\mathbb{R} ^{n\\times d}$, where $n$ is the size of the mini-batch and $d$ the input dimensions. We use the same settings as in the RNN demo and set $n=1, d=1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "BATCH_SIZE = 1\n",
    "INPUT_DIM = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will further assume that the we are interested in predicting a scalar value such as the stock price one day ahead, and thus set the number of nodes in the output layer to one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "NO_OUTPUT_NODES = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the number of nodes in the hidden state, we also consider the same setting of $h=5$ as in the RNN demo. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "NO_STATE_NODES = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With these settings, we can start with defining the data structures for our LSTM. These are the hidden state matrix $H_t \\in \\mathbb{R} ^{n\\times h} $, the input gate $I_t \\in \\mathbb{R} ^{n\\times h} $, the forget gate $F_t \\in \\mathbb{R} ^{n\\times h} $, and the output gate $O_t \\in \\mathbb{R} ^{n\\times h} $.\n",
    "<br>\n",
    "We also need a couple of weight matrices $W$ and biases $b$ to compute the values of the gates as follows:\n",
    "$$ I_t = \\sigma \\left( X_t W_{xi} + H_{t-1}W_{hi} + b_i \\right) $$\n",
    "$$ F_t = \\sigma \\left( X_t W_{xf} + H_{t-1}W_{hf} + b_f \\right) $$\n",
    "$$ O_t = \\sigma \\left( X_t W_{xo} + H_{t-1}W_{ho} + b_o \\right) $$\n",
    "<br>\n",
    "where $W_{xi}, W_{xf}, W_{xo} \\in \\mathbb{R} ^{d\\times h}$ and $W_{hi}, W_{hf}, W_{ho} \\in \\mathbb{R} ^{h\\times h}$ denote the weight matrices and $b_i, b_f, b_o \\in \\mathbb{R}^{1\\times h}$ the biases. \n",
    "\n",
    "<br>In order to compute predictions from our LSTM, we need yet another set of parameters, $W_{hy} \\in \\mathbb{R}^{h\\times c}$ and $b_{y} \\in \\mathbb{R}^{2\\times c}$, where $c$ denotes the number of output units. The prediction equation, \n",
    "with $g(\\cdot)$ denoting the output layer activation function, is then: \n",
    "$$Y_t = g \\left( H_tW_{hy} + b_y\\right) $$ \n",
    "\n",
    "Many parameters to initialize.  Let's go."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Define two anonymous functions to simplify the code\n",
    "init_x = lambda: np.random.randn(INPUT_DIM, NO_STATE_NODES)\n",
    "init_h = lambda: np.random.randn(NO_STATE_NODES, NO_STATE_NODES)\n",
    "\n",
    "# Initialize weight matrices\n",
    "W_xi, W_xf, W_xo = init_x(), init_x(), init_x() \n",
    "W_hi, W_hf, W_ho = init_h(), init_h(), init_h()\n",
    "\n",
    "# Initialize biases\n",
    "b_i = np.zeros(NO_STATE_NODES)\n",
    "b_f = np.zeros(NO_STATE_NODES)\n",
    "b_o = np.zeros(NO_STATE_NODES)\n",
    "\n",
    "W_hy = np.random.randn(NO_STATE_NODES, NO_OUTPUT_NODES)\n",
    "b_y = np.zeros(NO_OUTPUT_NODES)\n",
    "\n",
    "# And finally we have the hidden state, which we initialize to zero\n",
    "H = np.zeros(NO_STATE_NODES) # Since our hidden state gets update, we drop the 'index' t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that the LSTM maintains a seperate cell memory, $C_t \\in \\mathbb{R}^{n\\times h}$, in addition to the hidden state. We will also compute a candidate value of the new cell state, denoted by $\\tilde{C_t}$. Finally, we need a new set of weights and biases to compute the new candidate memory as:\n",
    "$$ \\tilde{C}_t = tanH \\left( X_tW_{xc} + H_{t-1}W_{hc} + b_c\\right) $$\n",
    "\n",
    "where $W_{xc} \\in \\mathbb{R}^{d \\times h}, W_{hc} \\in \\mathbb{R}^{h \\times h}, b_c  \\in \\mathbb{R}^{1 \\times h}$. <br>\n",
    "Here is the code to initialize the cell memory and its associated parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (Candidate) Cell memory\n",
    "C_tilde = C =  np.zeros((BATCH_SIZE, NO_STATE_NODES))\n",
    "\n",
    "# Parameters related to the cell memory\n",
    "W_xc = np.random.randn(INPUT_DIM, NO_STATE_NODES)\n",
    "W_hc = np.random.randn(NO_STATE_NODES, NO_STATE_NODES)\n",
    "b_c = np.zeros(NO_STATE_NODES)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The actual processing is similar to the previous [RNN example](https://github.com/Humboldt-WI/adams/tree/master/demos/rnn/rnn_foundations.ipynb). To show this, the following cell re-creates the synthetic data set we were using for our RNN in the previous demo. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#-----------------------------------------------------------\n",
    "# Copy-pasted from the rnn_foundations notebook\n",
    "#-----------------------------------------------------------\n",
    "import pandas as pd\n",
    "# load data\n",
    "data = pd.read_csv('AMZN.csv')\n",
    "data.set_index('Date', inplace=True)\n",
    "# Compute returns\n",
    "returns_series = data[\"Adj Close\"].diff(periods=1).to_numpy().reshape([-1,1])\n",
    "# Partitioning\n",
    "train_size  = int(len(returns_series) * 0.80)\n",
    "train, test = returns_series[:train_size], returns_series[train_size:]\n",
    "# Scaling\n",
    "from sklearn.preprocessing import MinMaxScaler\n",
    "#Standardize the values to a range from zero to 1\n",
    "scaler = MinMaxScaler()\n",
    "train  = scaler.fit_transform(train)\n",
    "test   = scaler.transform(test)\n",
    "\n",
    "# Helper function to create AR structure (note that Panda's shift() or numpy's roll() function could be used to accelerate the calculations)\n",
    "def create_dataset(time_series, window_size):\n",
    "    \"\"\" Function to create data set with lagged response values for estimating an autoregressive model\"\"\"\n",
    "    dataX, dataY = [], []\n",
    "    \n",
    "    for i in range(0,len(time_series) - window_size - 1):\n",
    "        x = time_series[i:i + window_size] # Remember: x[a:b-1]\n",
    "        dataX.append(x)\n",
    "        y = time_series[i + window_size] # Remember: x[b]\n",
    "        dataY.append(y)\n",
    "           \n",
    "    return np.array(dataX), np.array(dataY)\n",
    "\n",
    "WINDOW_SIZE=22\n",
    "\n",
    "X_train, y_train  = create_dataset(train, WINDOW_SIZE)\n",
    "X_test, y_test   = create_dataset(test, WINDOW_SIZE)\n",
    "\n",
    "#Reshape input to be [samples, time steps, features]\n",
    "X_train  = np.reshape(X_train, (X_train.shape[0],X_train.shape[1], 1))\n",
    "X_test   = np.reshape(X_test, (X_test.shape[0],X_test.shape[1], 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We begin with the first point in our input sequence, which is $X_{t-22}$ in our example, put it into our LSTM cell, perform all the internal computations to update our hidden state and the cell memory, and then proceed with the next point. So we first extract the minibatch of data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[       nan, 0.27370381, 0.28942156, 0.28371118, 0.27234689,\n",
       "        0.259343  , 0.28167578, 0.28913885, 0.2686719 , 0.26363992,\n",
       "        0.25917335, 0.28772543, 0.2778877 , 0.27065071, 0.2799796 ,\n",
       "        0.27596536, 0.28133657, 0.29513198, 0.26578842, 0.26929379,\n",
       "        0.25860797, 0.27749189]])"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "minibatch = X_train[0:BATCH_SIZE,:,:].reshape([-1,22]) # reshaping as in the RNN case\n",
    "minibatch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then start with the first value and compute the gates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[nan]])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = minibatch[:,0:1]\n",
    "X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Input gate\n",
    "$ I_t = \\sigma \\left( X_t W_{xi} + H_{t-1}W_{hi} + b_i \\right) $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[nan, nan, nan, nan, nan]])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.special import expit as sigmoid\n",
    "\n",
    "I = sigmoid(np.dot(X, W_xi) + np.dot(H, W_hi) + b_i )\n",
    "I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Forget gate\n",
    "$ F_t = \\sigma \\left( X_t W_{xf} + H_{t-1}W_{hf} + b_f \\right) $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[nan, nan, nan, nan, nan]])"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F = sigmoid(np.dot(X, W_xf) + np.dot(H, W_hf) + b_f )\n",
    "F"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Output gate\n",
    "$ O_t = \\sigma \\left( X_t W_{xo} + H_{t-1}W_{ho} + b_0 \\right) $"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[nan, nan, nan, nan, nan]])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "O = sigmoid(np.dot(X, W_xo) + np.dot(H, W_ho) + b_o )\n",
    "O"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cell state\n",
    "Next, we can proceed with computing our new candidate memory: $ \\tilde{C}_t = tanH \\left( X_tW_{xc} + H_{t-1}W_{hc} + b_c\\right) $, "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[nan, nan, nan, nan, nan]])"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C_tilde = np.tanh(np.dot(X,W_xc) + np.dot(H,W_hc) + b_c)\n",
    "C_tilde"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which, in turn, allows us to update the cell memory according to: $C_t = F_t \\odot C_{t-1} + I_t \\odot \\tilde{C}_t$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[nan, nan, nan, nan, nan]])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C = F * C + I * C_tilde\n",
    "C"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and finally our hidden state: $H_t = O_t \\odot tanH \\left(C_t \\right)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = O * np.tanh(C)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We would proceed with repeating the process using the next elements of our input, that is time-lagged observations $x_{t-21}, x_{t-20}, ..., x_{t}$ and update the state of our LSTM cell, as represented by its memory and hidden state. At any point in time, we could make a prediction using the current hidden state.\n",
    "#### Make a prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[nan]])"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Yhat = np.dot(H, W_hy) + b_y\n",
    "Yhat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This completes our journey through coding the forward pass in a LSTM \"by hand\". For comparison purpose, we could look at the next value in our time series to assess the previous prediction. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.27370381]])"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "minibatch\n",
    "minibatch[:,1:2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Apparently, our prediction is pretty far off. That should not come as a surprise since we did not train our LSTM. \n",
    "Based on the prediction and the true value, we could apply our loss function, backpropagate the error and update the weights based on their negative gradient. \n",
    "Backpropagation for RNNs and LSTMs works exactly as backpropagation for other neural network architectures on the *unrolled* neural network graph (the one without the loop). It is not easy to work out on your own though. A step-by-step walkthrough is available on our blog if you are intersted https://humboldt-wi.github.io/blog/research/information_systems_1819/02lstmgruandbeyond/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A note on notation\n",
    "The notation we used for the LSTM follow the lecture slides and draws upon chapter 9.2 of [Dive into Deep Learning (D2L)](http://d2l.ai/chapter_recurrent-modern/lstm.html). When it comes to LSTM, [colah's blog post on Understanding LSTM Networks](http://colah.github.io/posts/2015-08-Understanding-LSTMs/) is maybe the most famous source. I highly recommend reading the post, and when you do, please note that it uses a slightly different notation. Consider for example our input gate, which we defined as $ I_t = \\sigma \\left( X_t W_{xi} + H_{t-1}W_{hi} + b_i \\right) $. Colah's blog defines the same gate as $ I_t = \\sigma \\left( W_i \\left[ X_t, H_{t-1} \\right] +b_i \\right) $; same for other gates. Notation can be confusing when it differs across sources. Since Colah's blog is so popular, it is useful to see the equivalence between the two versions. The blog post assumes we concatenate the current input with the hidden state matrix of the previous step and then apply **one weight matrix** $W_i$. We, on the other hand, following D2L, consider **two weight matrices** $W_{xi}$ and $W_{hi}$ for the gates. You can easly convince yourself that both versions are equivalent. While writing down the mathematical notation is a bit tedious a little bit of numpy might just do the trick. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[nan, nan, nan, nan, nan]])"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# This is what we consider\n",
    "W_xi\n",
    "X\n",
    "W_hi\n",
    "H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.91185717, -0.82777974, -1.53530328,  2.18516906, -1.07571311],\n",
       "       [-1.15318561,  0.39309418, -0.94894074, -0.05277698,  1.15676662],\n",
       "       [ 1.10273649, -0.91311492,  1.7057517 ,  0.50535371, -1.7855959 ],\n",
       "       [-1.20352411,  1.02903543,  1.43188882,  0.46212581,  1.34001901],\n",
       "       [-0.79696491,  1.02946739,  1.5338228 ,  1.54174453, -2.71035077],\n",
       "       [ 0.47497536, -0.28324661,  1.3181224 , -0.06929086, -1.54099183]])"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# This is what Colah's blog considers\n",
    "np.hstack((X, H))\n",
    "np.vstack((W_xi, W_hi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.38777878e-17,  1.38777878e-17,  0.00000000e+00,\n",
       "         0.00000000e+00,  0.00000000e+00]])"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# And yes, it is the same\n",
    "d2l = np.dot(X, W_xi) + np.dot(H, W_hi)\n",
    "colah = np.dot(np.hstack((X, H)), np.vstack((W_xi, W_hi)))\n",
    "d2l-colah"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Concluding remarks and notes\n",
    "\n",
    "### Notation\n",
    "The notation we used for the LSTM follow the lecture slides and draws upon chapter 9.2 of [Dive into Deep Learning (D2L)](http://d2l.ai/chapter_recurrent-modern/lstm.html). When it comes to LSTM, [colah's blog post on Understanding LSTM Networks](http://colah.github.io/posts/2015-08-Understanding-LSTMs/) is maybe the most famous source. I highly recommend reading the post, and when you do, please note that it uses a slightly different notation. Consider for example our input gate, which we defined as $ I_t = \\sigma \\left( X_t W_{xi} + H_{t-1}W_{hi} + b_i \\right) $. Colah's blog defines the same gate as $ I_t = \\sigma \\left( W_i \\left[ X_t, H_{t-1} \\right] +b_i \\right) $; same for other gates. Notation can be confusing when it differs across sources. Since Colah's blog is so popular, it is useful to see the equivalence between the two versions. The blog post assumes we concatenate the current input with the hidden state matrix of the previous step and then apply **one weight matrix** $W_i$. We, on the other hand, following D2L, consider **two weight matrices** $W_{xi}$ and $W_{hi}$ for the gates. You can easly convince yourself that both versions are equivalent. While writing down the mathematical notation is a bit tedious a little bit of numpy might just do the trick. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.09024775,  0.47879233,  0.16190862, -1.94446208,  0.47213937]])"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([[0.27370381]])"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([[-0.58085862, -0.0363312 ,  1.08331669,  1.19785893,  0.69495794],\n",
       "       [ 1.63013443,  1.37963973,  1.00650678, -1.41765707, -1.09447956],\n",
       "       [-0.93631225,  0.11963822, -0.0278419 , -0.36064368, -1.21106708],\n",
       "       [ 1.75247636, -1.48884166,  0.47270972, -1.16958587, -0.05062323],\n",
       "       [ 0.88035388, -1.58002961,  0.37723041, -0.76339754, -1.88443342]])"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([[ 0.07363884, -0.08464281, -0.08260762,  0.02893675,  0.02347467]])"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# This is what we consider\n",
    "W_xi\n",
    "X\n",
    "W_hi\n",
    "H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.27370381,  0.07363884, -0.08464281, -0.08260762,  0.02893675,\n",
       "         0.02347467]])"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "text/plain": [
       "array([[-0.09024775,  0.47879233,  0.16190862, -1.94446208,  0.47213937],\n",
       "       [-0.58085862, -0.0363312 ,  1.08331669,  1.19785893,  0.69495794],\n",
       "       [ 1.63013443,  1.37963973,  1.00650678, -1.41765707, -1.09447956],\n",
       "       [-0.93631225,  0.11963822, -0.0278419 , -0.36064368, -1.21106708],\n",
       "       [ 1.75247636, -1.48884166,  0.47270972, -1.16958587, -0.05062323],\n",
       "       [ 0.88035388, -1.58002961,  0.37723041, -0.76339754, -1.88443342]])"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# This is what Colah's blog considers\n",
    "np.hstack((X, H))\n",
    "np.vstack((W_xi, W_hi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1.38777878e-17,  1.38777878e-17,  0.00000000e+00,\n",
       "         0.00000000e+00,  0.00000000e+00]])"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# And yes, it is the same\n",
    "d2l = np.dot(X, W_xi) + np.dot(H, W_hi)\n",
    "colah = np.dot(np.hstack((X, H)), np.vstack((W_xi, W_hi)))\n",
    "d2l-colah"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Why Tanh activation?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The original RNN uses tanh rather than the - at that time typical - sigmoid activation, because the tanh gradient lies in [0;1] rather than [0; .25]. Increasing the values of the activation function derivatives helps mitigate the vanishing gradient problem for lower layers and helps control the size of values in the hidden state. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABCHUlEQVR4nO3dd3iTZffA8e8pFJAhskQ2iKgsQajgBhcqiqivCrhw8qr4qu/PAW4UBy7coqgoooIT4VVUcIIbiohskCFl7z1aen5/nMS00EKbNn2S9nyu67mSO3mS3OnIyb3OLaqKc845l19JQVfAOedcYvIA4pxzLioeQJxzzkXFA4hzzrmoeABxzjkXldJBV6AoVa9eXRs2bBh0NZxzLqGkpqauVtUau99eogJIw4YNmTRpUtDVcM65hCIii3K63buwnHPORcUDiHPOuah4AHHOOReVEjUGkpP09HTS0tLYvn170FVJGOXKlaNu3bokJycHXRXnXIBKfABJS0ujUqVKNGzYEBEJujpxT1VZs2YNaWlpNGrUKOjqOOcCFGgXlogMEZGVIjItl/tFRJ4TkXkiMlVE2mS57wwRmR26r2+0ddi+fTvVqlXz4JFHIkK1atW8xeacC3wM5E3gjL3cfybQJHT0AgYBiEgp4MXQ/c2AHiLSLNpKePDIH/95Oecg4C4sVR0vIg33ckpX4C21nPO/iMgBIlILaAjMU9X5ACIyInTujBhX2bnCs3Qp/O9/cNZZULcuLFsG48fDGWdA5cpB187FO1XIyICdO+1IT48cGRmRy/Bx6KGF/ncV72MgdYDFWcppodtyur19Tk8gIr2w1gv169ePTS0LYP369bz77rvccMMNUT3+iiuu4Oyzz+aCCy7457aKFSuyefPmwqqiK0zTpsFBB0H16jBjBlx3HfzyiwWQceOgZ087p3JlWLECSpWyc13iSU+H9ethwwY7Nm60Y9Mm2Lw5cmzZAlu3Ri63bbPL7dvt+o4ddn3HjsgRDhr52c/piy/g9NML9S3GewDJqa9E93L7njeqDgYGA6SkpMTd7lnr16/npZdeijqAuASyeDEccQT06wf33QfHHQd//QXhyQg9ekDz5nD44VZ+6CEYOhQWLIBq1QKrtgvZts1ajcuXW2txxQpYudKONWtg9WpYuzZybNmSt+ctX37Po1w52G8/qFLFLsuVg7Jl7ShTJnIZPpKTI5fho3TpyGXp0tCmzb7rkk/xHkDSgHpZynWBpUCZXG5POH379uWvv/6idevWnHTSSUydOpV169aRnp7OQw89RNeuXVm4cCFnnnkmxx9/PD/99BN16tRh1KhR7LfffkFX3+1Lejr8+iscfzzUqwfDh8Opp9p9++0HBx8cOTc5Gdq2jZSvuw5OPDESPObMsW4IFxsbN8LcuTB/vgXthQvh778t8KelWVDYnQhUrQo1atjvqUED+6CuUsWOAw6w1mTlylCpEuy/vx0VK1q5fHlICnooOnrxHkBGAzeGxjjaAxtUdZmIrAKaiEgjYAnQHbi4wK92yy0wZUqBnyab1q3hmWdyvXvAgAFMmzaNKVOmkJGRwdatW9l///1ZvXo1Rx99NOeccw4Ac+fOZfjw4bz66qtcdNFFfPTRR1x66aUA3H777Tz00EOFW29XOO6/H556yloadetCt255f2zz5nYA/P47HHUUDBtmLRUXHVVrOUybZseMGTBrlh2rVmU/t0oVCwgNG9oXgDp1oHZtqFXLuiFr1rTuxdLx9TG6eXP2xtHKlfbWevSINHYLS6DvXESGAx2B6iKSBtwPJAOo6svAGKAzMA/YClwZui9DRG4EvgRKAUNUdXqRv4FCpqrcddddjB8/nqSkJJYsWcKKFSsAaNSoEa1btwagbdu2LFy48J/HPfHEE3uMgbg4cdtt0KqVBY+CaNYM7rrLBtxd3i1bBr/9ZsfkyXasXBm5v0YN6zLs2hWaNIFDDoHGje2Tdv/9g6v3blRh3TpYssR60bIeWXvUli+34ZOcHHFEMQsgqrrXr1Kh2Ve9c7lvDBZgCs9eWgpF4Z133mHVqlWkpqaSnJxMw4YN/1lvUbZs2X/OK1WqFNu2bQuqmm5ffvoJXn8dXnnFujfy0+rITdmy8OCDdj093Vo2t91mz++MqnVBffst/PCDHeEvWqVLW2uuc2frFTjiCGjRwgJIHNixw3rLFi6ERYvsetbes7S0nANDlSqRBtExx1ijKHwceKAdNWrYEYse7/hqe5VAlSpVYtOmTQBs2LCBAw88kOTkZL799lsWLcoxg7KLd5MmwYQJ1md+4IGF//xTpsDAgfaBeMklhf/8iWT9epu99vnndpmWZrfXrAknnAA33QTt28ORR8bmEzSPVG2cfd48682cN8+GWsLDLUuXZp9QlZRkgaFePYt3Z59tPWhZj1q1bGw9SB5AAlatWjWOO+44WrRowVFHHcWsWbNISUmhdevWHB6ejeMSQ2am/effdBNcfTVUqBCb1znqKJg92/rnS6IlS2DkSDu+/x527bLB6lNPhVNOgZNPtu6oABa87thhcx1mzbJf0ezZVp4zx2JdmIj1ah58sFW7USMbamnY0H6tderYnIp4J5qfecQJLiUlRXffUGrmzJk0bdo0oBolLv+57WbmTOuqevdd6xopKn/+CY88Am+8EfzX0Vhatw7ef99+vhMm2Nf1ww+Hc8+1caGjjy7Swez0dAsKf/5pY/HTp9vx11/2PSKsXj047DCLZ+HjkEMsUCTSr0tEUlU1ZffbvQXiXGHIyLB5+LFqdeRm5kz7QE1Ls0+m4kQVvvkGXnvNWhs7dljQ6NcPLrwQiugLzLp11msYPv74wyZvpafb/aVK2ezqI46A7t2tWocfbrcV9Z9DUfMA4lxBqFp/RMuWMHFi0XebXHSRfQMvTp9UGzfaAsoXX7Q+oCpV4Npr4corbSwjhj/j9eshNdWGsSZOtElbCxZE7q9Vy8YkzjjDfuUtW1oLI8sclxLFA4hz0VK11CPNmkHfvoH0uQMWPFThySehfv3CmfUVhOXL4dlnYdAgS/3Rvj289Za1NmLQ35ORYd1PP/9s2WR+/dXiVdjBB0NKCvTqZXHryCNjMycikXkAcS5a4WR1u3YFXROrx+jR1o2VaAFk+XJ49FGb9rxzJ/zrX3D77dCuXaG+zObNFigmTIAff7SAEU4Zd+CBFq8uu8xetm1bnyGdFx5AnItWmTLwzjtB18IkJ8OYMZYiI1GsX2+B4/nnLXD07Al33lloYzmbNtlSkO++syM11WJ9UpKNV/TsCccea+snGjYMrgGZyDyAOJdfq1fDf/5jXUZ16gRdm4hKlexyxQqbldWnT3x+Kqanw8svwwMP2FqZSy6xhZEFDBw7dlh31Fdf2dj7b79ZwEhOtlZFnz6WWuyYY+JqkXlCS9wsXsVEOBtvtK644go+/PDDbLeFU5ksXLiQjh07FqR6LidTpsDYsXvmTooXw4fbh/OsWUHXZE8//GD9QzfdZCleUlMtv1cUwUPVptI+95zNI6haFU46CQYMsPv69LFgsn69vezDD1s2cw8ehcdbIAHzdO4J6NRTLd9EvHYX3XSTLV2Op2m969fDrbfCkCE20D9ypOWfymcLaft2y1QyZowd8+fb7U2a2CSt006Djh19P66i4gEkYLFM516qVCmq+khg4fnzT0tQdNZZ8Rs8wDr5w8Hjhx9s5XqQ80w/+8ymMq1YAXfcYXuh5GPa8cqV8OmnNkdg3DjLCVW+vC06v+02m1Jb2EkCXR6paok52rZtq7ubMWNG9hs6dFB94w27vnOnlYcNs/KWLVYeMcLK69db+aOPrLxqlZVHj7bysmV7vN7uFixYoM2bN1dV1fT0dN2wYUPoqVZp48aNNTMzUxcsWKClSpXS33//XVVVL7zwQh0WqlPPnj21YcOG2qpVq3+OChUq7PN1C2qPn1tJcPHFqjVrqm7eHHRN8mbGDFUR1QEDgnn9zZtVr71WFVRbtFCdODHPD50/X/Wpp1SPO87eAqjWr6/au7fq55+rbtsWw3q7PQCTNIfPVG+BxBH1dO7xbcgQy/aaKIv2mja11B+hPWWK1B9/2LLs2bOt1fHgg/tsBc2dCx9+aMfkyXZbq1Y2vt61q12PxzkBJZkHkN19913kenJy9nL58tnLlStnL1evnr180EH5emlP5x6n5s+32VZlyxZtnqvC0L27Xaan21G+fOxf8803bTfFqlWtz+mUU3I9deFCGDEC3nsvspfb0UfDE0/A+edn37DRxR8PIAHzdO5xLj3dOtmbNoVRo4KuTXTS0y21eZs2UIAZf/u0c6ft6jlokGXEHTEix/02Vq2K5EX86Se7rX17y1B/wQWWgNAlhqB3JDwDeBbbVfA1VR2w2/23A+END0oDTYEaqrpWRBYCm4BdQIbmkCkyEXg69ziXnGyfbOE1FokoOdkG/g87LHavsWYNnHeeLfO+4w6bM5slO+727TYIPmwYfPGFLZw/4ghbR9i9uy3kc4knsHTuIlIKmAOcBqQBE4Eeqjojl/O7AP9V1ZND5YVAiqquzutrejr3wlMifm7hRIlu7+bMsQC1eLF1X4W6zVQtXcgbb1gX1YYNtgfGJZfY0bJlsNV2eReP6dzbAfNUdT6AiIwAugI5BhCgBzC8iOrmSrrt26FDB7j5Zrj44qBrU3jefttGqx94oHCe7+efLXiULm0LNI45hpUrLQfikCGWbX6//Sy91RVX2BqNUqUK56Vd8IJciV4HWJylnBa6bQ8iUh44A/goy80KjBWRVBHplduLiEgvEZkkIpNWxevKYRd/NmywJcvVqgVdk8L166+W5yO8mUVBjB1riyqrVSPzp1/4cuMxXHCBzTe4/XbbJPDVVy1X4rBhNpbuwaN4CbIFklPfQG79aV2AH1V1bZbbjlPVpSJyIDBORGap6vg9nlB1MDAYrAuroJV2JUTNmvYBWdy6sJ54wpJAJhXwu+OHH8LFF7OiyfG83nUUr55WiYULbSLizTfbjr7FvYfTBdsCSQOyzreoCyzN5dzu7NZ9papLQ5crgZFYl5hzBZORYWsW1q0rfsEDbF+NpCTLY/7551E9hb7/ARO6vUCPymOoN/dr7n60Eo0a2aSrtDTLMenBo2QIMoBMBJqISCMRKYMFidG7nyQilYEOwKgst1UQkUrh60AnYFqR1NoVbz//bAHk66+Drkls3XOPzZoKLVTNiy1bYPC/J9G626GcmPkdX2ScQu/ewqxZ1ivWrVvJ3ZmvpAqsC0tVM0TkRuBLbBrvEFWdLiLXhe5/OXTqecBYVd2S5eE1gZFi3xBLA++q6hdFV3tXbJ1wgmWxbdw46JrE1j332E5/NWvu89RFi2x32dcG7WTd5hRaV5jDa49tp8eV5YpkXaKLX4Gmc1fVMap6qKo2VtWHQ7e9nCV4oKpvqmr33R43X1VbhY7m4ccWJ9dccw0zZuQ2Ia1wdO7cmfXr1+9xe79+/XjyySdj+tpxJ5wbHCwRYXHsvsqqenU47ji7vmXLHner2q59F1xgq8EHPpXJqVtHM6HZv5m8rDZX9/bg4Xw/kLj12muv0axZs5i+xpgxYzjggANi+hoJY/hw29s8vDS6pHjvPVvFt9gmRGZk2E3t28Pxx1vX1O2XLmNBuWa83+wBjp/wKFLJc6054wEkDmzZsoWzzjqLVq1a0aJFC9577z06duxIeNHj66+/zqGHHkrHjh259tprufHGGwHbTOr666/npJNO4uCDD+b777/nqquuomnTplxxxRX/PP/w4cNp2bIlLVq0oE+fPv/c3rBhQ1avtnWYDz/8MIcddhinnnoqs2fPLro3Hy/OPBP697dETCVJu3bQuTOb0svxzDPWc9e9u23f8eKLsHj8AgZ83op6NXfCl1/6RuEuG8+FlcUtt0QSuhWW1q3hmWf2fs4XX3xB7dq1+eyzzwDLiTVo0CAAli5dSv/+/Zk8eTKVKlXi5JNPplWrVv88dt26dXzzzTeMHj2aLl268OOPP/Laa69x1FFHMWXKFA488ED69OlDamoqVapUoVOnTnzyySece+65/zxHamoqI0aM4PfffycjI4M2bdrQtm3bwv1BxDNVqFLF9uMuYVaUb8SzdYYyqK0FjRNOsC3Kzz4bktavhWPPtH1hv/gCatcOurouzngLJA60bNmSr776ij59+jBhwgQqZ9lO7bfffqNDhw5UrVqV5ORkLrzwwmyP7dKlCyJCy5YtqVmzJi1btiQpKYnmzZuzcOFCJk6cSMeOHalRowalS5fmkksuYfz47MtlJkyYwHnnnUf58uXZf//9OSeI9N9BGTsWOnXK12yk4mDePEuY26CBbQF7yrHb+OXshxg/ah3nnANJ6TtsltaCBfDJJ3DooUFX2cUhb4Fksa+WQqwceuihpKamMmbMGO688046der0z337ylUWTvOelJSULeV7UlISGRkZlC6dt1+xFPdB49ysWWOrzkvIHqhTpljA+OADy7HYs6ft6tdk00zo+Dj8dpQF1BtugPHjLWXuCScEXW0Xp7wFEgeWLl1K+fLlufTSS7ntttuYHN5NB2jXrh3ff/8969atIyMjg48++mgvz7Sn9u3b8/3337N69Wp27drF8OHD6dChQ7ZzTjzxREaOHMm2bdvYtGkT//vf/wrlfSWEHj0svUe5ckHXJKZ++slSVh15pO0lfvvt1rh45RXbT5w2bWwg/fTTbfBjyBCb6tujR9BVd3HMWyBx4M8//+T2228nKSmJ5ORkBg0axG233QZAnTp1uOuuu2jfvj21a9emWbNm2bq49qVWrVo8+uijnHTSSagqnTt3pmvXrtnOadOmDd26daN169Y0aNCAE0rCN87Jk2HZMvtULaatL1WbRfXQQ7bPWfXqdr13b8tTtYfwBmk33wydOxdewkVXfOW0z21xPfK0J3oc2rRpk6ranulnn322fvzxxwHXKDF+bnvVo4dqrVq2z30xk5mp+tlnqkcfbXuJ166tOnBgHrZyT0tTrVLFHtSvX5HU1SUGfE/0xNWvXz+++uortm/fTqdOnbLNoHJRevNNG0kuRqvhVG3Tpv79ITXVBsgHDYIrr8xDipGMDJu/u3OnbaD1738XSZ1dYvMAkgBK3KrwWFqxwvpvypa1hYPFQGamTZR68EH44w9bOf7663DZZTZQnif33AM//ADvvBPZ/yQjwy7zOBHDlTw+iO5KDlXL/3TqqXY9wWVmwscf28D4v/4FW7fC0KEwezZcdVU+gseYMfDYY9bqCAePDRtsUeXAgTGrv0t8/tXClRwi0LevfdIm8MB5uMXxwAMwdaot0Xj7beuByveGTcuX21aBRxyRfR77/vtDSorlBXMuFx5AXMkQ3t+8c+egaxI1VRg1Cvr1s66qQw+1nf569Ihypz9VGyDZtMlmX2WdyiwCL7+c60OdA+/CciWBqvXxDBkSdE2iogr/+x+0bWuLw7dssT3Hp0+HSy8twDaxL7xgKUqeeir38SBVe7HXXou6/q748gASp+IlnfsVV1zBd999F9N6xNymTbYD386dQdckX1Tt8719ezjnHBuWePNNmDnTBsgLNLY9c6atJjzrLLj++r2f+/77tnS9GIwbucLlXVhx6rUi+MY3ZsyYmL9GXNh/f8skmyBUbUPE++6zDRIbNLAGwOWX52NgfG8yMmzco2JFm661t/EgEZuZVbFiQo8budgItAUiImeIyGwRmScifXO4v6OIbBCRKaHjvrw+NpHEczr3ypUrU6ZMmSL4KcRAZiY88khkf/ME+AAcPx46doTTTrPMIi+/bPtcXX11IQUPgCeegN9+g5deytOOhFSubP1kmzbBuHGFVAlXLOS0urAoDmwb27+Ag4EywB9As93O6Qh8Gs1jczryshK9QwfVN96w6zt3WnnYMCtv2WLlESOsvH69lT/6yMqrVll59GgrL1u2t7WdER9++KFec801/5TXr1+vHTp00IkTJ+qSJUu0QYMGumbNGt25c6cef/zx2rt3b1VV7dmzp3br1k0zMzP1k08+0UqVKunUqVN1165d2qZNG/399991yZIlWq9ePV25cqWmp6frSSedpCNHjlRV1QYNGuiqVat00qRJ2qJFC92yZYtu2LBBGzdurE888cRe65wQK9F/+UW1dGnVoUODrsk+/fST6imn2CLwWrVUn3tOddu2GLzQn3+qJierXnhh/h973XWq5currl5d+PVycY1cVqIH2QJpB8xT2552JzAC6LqPxxTGY+OOp3OPkfbtYdo0GzCIU6mpNjHs2GNtSu7AgfDXX/Cf/8Qgv+OuXbZA5IADLGFifvXrZ12B1aoVcsVcogpyDKQOsDhLOQ1on8N5x4jIH8BS4DZVnZ6PxyIivYBeAPXr199npbKOFycnZy+XL5+9HM49F1a9evbyQQft8+UAT+de6HbutMDRpg0cdljQtcnR1Kk2xjFqlG3yN2AA3HgjVKgQwxd9/nmYONG2761RI/+Pr1kz0uW1ciUceGDh1s8lnCBbIDl9Yu3+aTkZaKCqrYDngU/y8Vi7UXWwqqaoakqNaP5pioCncy9kTz5pW7XOnRt0TfYwfTpcdBG0amVfNh54wNKq9+kT4+CxcKGlK+ncGbp1K9hzffmljexPmFAoVXOJK8gWSBpQL0u5LtbK+IeqbsxyfYyIvCQi1fPy2ETi6dwLWe/eUKtWaKOL+DB7tuWqGj7cAsXdd8Ott9pOujGnGpmqO2hQwScTHHccXHttsckl5gogp4GRojiw4DUfaERkILz5buccBEjoejvgb6z1sc/H5nR4OvfCE5c/t02bVDMygq5FNnPnql5+uWpSko0/33GHTbYoUu+9Z6Pzzz5b+M+dmWmHK9aIt0F0Vc0AbgS+BGYC76vqdBG5TkSuC512ATAtNAbyHNA99H5yfGzRv4ui0a9fP1q3bk2LFi1o1KiRp3PPiaolAjznnLhY8DZ/vo1XH364rcH773+tq+qxx2ysrMhs2mQvfuSR1jIr7Oc++2xfpV6CBbqQUFXHAGN2u+3lLNdfAF7I62OLK0/nnkddu0J6eqDrPRYsgIcftqy4pUvbbKo+ffI+oaLQ9etnOy9+/HEBcp7kokKFwn9Ol1B8JTrWjVesZiHFmMbBN/xswokSr746sCrMn29rFocOtc/UG26wwFG7dmBVsqlezz5r4xXtc5ykWDBJSTaNzP93SqwSnwurXLlyrFmzJv4+FOOUqrJmzRrKFfoihSitXm1px7/9NpCXnz/f4tZhh1lK9euus3Uczz4bcPBQtXnBBxxgkS1WwsHjm2+sueX/RyVKiW+B1K1bl7S0NFatWhV0VRJGuXLlqFu3btDVMJs22Tfh/fcv0pedN8+6qoYNs66q66+3rUYCDRpZvf++TbN95ZWiWfj366/w1Vewdq0vNCxBpCR9805JSdFwfilXjIS7sIrAzJkWOIYPhzJlrMVxxx02azhubN1qo/fVq9vCwaIYp8jMhG3bYryYxQVFRFJVNWX320t8F5ZLUB99BLfcYplliyB4TJliCwCbN4eRI+H//s8GzJ9+Os6CB9hUr8WLrR+tqAa5k5IseGRkwP33Q1pa0byuC5QHEJeYUlMto+yuXTF9mZ9/hi5dbBbsl1/CnXfaou4nnghwZtXeLFoEjz9u+9sGsSD0778tqn7wQdG/tityJX4MxCWoRx6B7dshS/6vwqJqWcsffdTSjVSrZqvI//MfG5OOa3fdZZePPx7M6x98sOVrqVdv3+e6hOctEJc4tm2zXZXmz7dyIc8E27XLxp5TUuD0020fjoEDrcVx770JEDx++w3efddypAT5AR5+7QULbNtcV2x5AHGJY+5c+Pxz+PPPQn3abdts46bDDrM8g5s32+Lq+fNtEXfFioX6crGhagMzNWvaApR48MILlnJ45cqga+JixLuwXOI44ghbZFFIU3ZXr7ZN+V54AVatsgS+jz9uC9oTboH1xx/Djz/atN1KlYKujQnnqPe078WWt0Bc/HvsMdu7GwoleMyebes26te3CUPt2tlYxy+/wPnnJ2DwSE+3RSjNm1sCrniRnAyNGtn1IUNsrYgrVrwF4uLbrl3w/fc2kn3VVVFP2VW1xdLPPAOffmpj75ddZl1UCZ+V/NVXbWXjp5/aqsZ4s22bTXo4+ujYpFRxgfGFhC5+hRcI7txpl8nJ+X6KrVttXPm552zopEYNa33ccENkc72EtmkTHHIING1q6VziNS/VkiX2wy9TJuiauCj4QkKXWN57D846yyJAmTL5Dh4LFtgK8Xr1LJdgUpL1ovz9t+0CWCyCB9g0sZUrrZsvXoMHQJ069nvcsiX7TDqX0OKwvesc1u2xbVu+kvNlZsLYsTYw/umnFjS6doWbb7Y1dfH8+RqVlStt+95//StxuoYWL7YVmV262JoRl9AC7cISkTOAZ4FSwGuqOmC3+y8BwnMSNwPXq+ofofsWApuAXUBGTs2r3XkXVgJYsyaSjC8z06LAPqxcCW++aROQ5s+3ST+9esG//w3xkvMxJm6+GV58EWbMgEMPDbo2ebdxY2QyRHp6VF2TrmjFXReWiJQCXgTOBJoBPURk9+HMBUAHVT0C6A8M3u3+k1S1dV6Ch0sAb79t/fkzZ1p5L8EjM9OSv3brZkGiTx/rJRk+3L7k9u9fzIPHokW2eOXKKxMreEAkeEycaEkfp04Ntj4uakGOgbQD5qnqfFXdCYwAumY9QVV/UtV1oeIvQHH+SHAnnmg5nMJTP3OweLFlwz3kEDjtNEs50ru3fQkfP94eHk/jtDt2WLADWLcOJk+228DGlceOtZ46sPGZceMi969caUkcMzKsvHNn5DoPPmiX991XFG8jNg44wKJ8jRpB18RFKcgAUgdYnKWcFrotN1cDn2cpKzBWRFJFpFduDxKRXiIySUQm+Z4fcWjnTmt5qNrCjEGD9khRsm0bjBgBZ5wBDRrAPfdAw4Y2u2rpUsvd17Rp7KuqattdbN9u5YULbex66VIrjx9vQxFz5lj5vffsrcyebeVPP4W2bS0Igg0FnH66LWIEGDMGOnWC9eut/O67lsRx40YrP/ec9fZsSp0Db77J6ycO5fge9f4JON9/b9OUwwFr/XqbgxC3mjSxBTi1atkP98cfg66Ry6cgA0hOQ5o5DsiIyElYAMmao+E4VW2DdYH1FpETc3qsqg5W1RRVTanh33Tiz7BhtiDj55+z3ZyZaR/I115rWW979LCerXvvtcXo33xjtxVmOixV+zAPf2D//be9/uTJVv7hBxuemTDByosW2fq9cI/bfvvZl+rwB3irVrb8oWpVK3fsCJ98Esnie9ZZ9pkZnhF23nn2GuHzu3SxrPXhHp/jjrMZZBUeuw/224+y551JhQqRFteoUfbzCff89etnzx0e5hwyxOobtmKFpW0JVHhmw1tvwfHHR364LjGoaiAHcAzwZZbyncCdOZx3BPAXcOhenqsfcNu+XrNt27bq4sSuXZHLr7765+apU1XvvFO1QQNVUK1QQfXyy1W//jrykILKzLTLTZtU779fdfx4K//1l73m669becEC1Zo1VT/5xMqrVqkOHKg6f76Vd+5U3by5cOqUZ1OmWCXvvnuPu3btUl23LlL+/nvVl1+OlG++WfXooyPl889XPfzwSHnw4Mh7V438nIrEjh1WgSJ9UZdXwCTN6bM3pxuL4sCmEM8HGgFlgD+A5rudUx+YBxy72+0VgEpZrv8EnLGv1/QAEic++US1VSvVNWtUVXXOHNX+/VVbtLC/yFKlVE8/XfXttwv2AZ2ZqTp5sur06VbeuVO1cWN7LVX7zEpOVn38cSunp6s+/bTqzJnRv2bMde2qWrmy6tq1BX6qL79UHTEiUj7lFNUzz4yUTzxR9eqrI+Xx41UXLizwy+7b8uX2BzBjRhG8mMuL3AJIYOtAVDVDRG4EvsSm8Q5R1ekicl3o/peB+4BqwEtiTd3wdN2awMjQbaWBd1X1iwDehotGzZrMLdeSD54ozQdf2EAxWBfNCy/AhRdGn3/vmWcse+4111jvyFln2bjCm2/a+MGZZ0ZSl5QpY1044S6g0qVtk8O4NWmS9VM98ABUqVLgp+vUKXv5q69sSCrslFMi3Wuq1sV2/vkwODQX8tZb7TlOP73AVclu4UKYNSt7ZVxc8lQmrkjolD+YMnIBI/VcRo6EadPs9qOPtoBx4YV528IiPR2WLbPxdrA9ydeutX08ADp0sA+9cHn8eJvoUyzWrHXubAkJFywotIzEeaVqs24rVrQAvHWrTWS49VabQr1tm+Vy7N8fLrnEUpilpdnvKaoFnFnXh7z3nkWpuN+QpfjKbR2Ir0R3MbNzp80MGj0aRr9Rh7+3tCIpSTnhBOGZZ+zb7L6CxqxZ1kLp3t3KPXvaePuCBVZu0CD7l/Gvv86eT/DEHKdWJKCff7a9UB59tMiDB1gQaNcuUi5f3gbhwzsKb9xoXwbCLZbZsy2gvP22BZTVq23K8umnR9aJ7lU4eKSlWeqT//7X0sO7+JJTv1ZxPXwMJPaWLLGB2H/9S7VShQwF1XLlVLucsUNff2ajrlix98d//73qVVepZmRY+Z57bExk61Yrf/ON6rvvlsCx1k6dVGvUsJH/BLBypeqgQaqLFln5k09sfOvnn608aZLqTTepLl2ahyebNCnyvhcuLJTxH5c/xNsgehCHB5DCt327zZC64w4bF7fODtU6tXdprzJD9H8dn9QtW7I/Zv16G8BWVR07VrVpU9XFi608dKjNfPr7bysvWWKfGSUuYGT144/2Q33iiaBrErWdO20C2fbtVn77bdXy5VVXr7bykCGq7dtHZpFl/RvJpmNH1WbNCm9KnsuT3AKIZ+N1+bJrl62LePJJW9hXpYoNtg4cCFWqKI9eO58//oDFaUm88r86nDzsSn74wVZdg3VpHXBAZM1YtWq2niy84O3SS2H58kjXVu3a1k1V7BIh5ke/fjar4Prrg65J1JKTbV1M2bJWvuQS6/YKd2dVqGDXK1e2cv/+9pbDXWQTJ9rBwIHWlZWUZN9VZs0q8vfissgpqhTXw1sg+ZeRYVNhn37aZpBWqRJpZRx+uOp//qM6erTqxo2q+uqruokKeu+Vi/W77+zx8+fbueH1CGvWqD78sK25cHnwww/2A3zyyaBrUqS++Ub12Wcj5S5drKUa9uKLqq/1nqwqYk1gF1Pk0gLxWVgum61b7Zvejz/a8cMPkZXZjRvbLKfjj4eUFGjZEjK/G89pt7XizO6Vue2GraR/OIoqN3Tn3nuFPn0s1Hz3naXwCGDsN/GddpolG5w/376ml1CLF9ugfUpoHtBJJ8H++6UzqsNAuPVWzr2gNEdUTePB+3dBgwYsWWIZUvKQzNnlgc/CcntQtZ1Qf/vN9gP/+Wf4449Iwr5mzSxdSNWqcNhhNgMKLNfhscfCO0MzSLricmqkv0XFiidC+fIkX96DNd0iXRUi9s/uovDjj7Y448knS3TwAOvSzDpj79tvYceOZChr2Y2qVVUqfz4Cpr2P/vobzZtbd+gLL9j5r75qX3yKImdaSeItkBJC1XI3pabaMWmSHetCuY7Ll7dpmsccY+ssSpWC116z+04+2TaS+/VX4NNPGXbvHA4acAunnZ5kCzoaN7ZEUK5wnX46/P67zVku4QEkTxYvhlWryDiiDW+9ns6hHz7C8Y91YXX9NtSoYXH41ltt8ej558Ptt1sDLyPD/r7D4y9uT94CKUF27LAEf1OnWovi999tLUU4WJQqZYPTF1wARx0FX3xhrZBvv7X777gjy46jqgz5v+lUadMIqACbNnFZ0jvQsjtQG1q0KPo3WBL88ostnHj8cQ8eeRVqppQGrjp6Btz3Eqw7nmpHQtr0DZTN3AYcxKpVtvg0Pd0eNmOGDfB//LGttl+2zBb8n3tuJPGly0VOAyPF9Shug+g7d1repo8+Un3wQdVu3WyGY6lSkYHucuVUGzVSrVpV9YUXVH/7zfLwiahu22bP88UXdl+2qbLhaZKpqfZEr70Wub1Ez6ktImeeqVq9esKs+4hL27dH/o4fecT+MZYt2+O0tDS7O5zna+RI+5OfONHK48aptmunOneulZcutX+LnTtj/xbiBb4OJDEDSGam6ooVNhnnjTdU+/a12VCNG6uWLh0JFKBat67NVunc2cpjx1qCwFGj7DHhOferV/+Tx3BP6emqxx0XyfaamWmLM3zxVtH59Vf7BT76aNA1KT7mzLGVjWE33qh6xRU5npqZaeuQwutQvvpK9dRTI/8/zz1nv57ly608cqTqZZdFEn8uX25Bpjh9z8otgHgXVhxIT7fu2wULrOto3jz4809bO7FwYWQWFNh8+oMOsvMvvti6ybdts+nxH3wARxxh2R+mT7eB7tKl4Zxz7AjbI5XEI49YromBA+0BRx4ZSR4lYqkkXNHp399mLvTuHXRNio8mTewIq1Ile86b3r1tELBnT0SyD9ifcoodYeedZ9snhxN+Ll1qPY7hYcCBA+HZZ21Gowi88YZ1IT/7rN3/99/WjVxnb9vnJQgPIDGmoV3s0tLsCAeK1FTYsMGmJqalRTb9Afvj2rXLpsleeqkN7r36qq2f6tnTxjI+/dQW8tWqZY/5978jj69bdx/7gb/9Nnz2mW0gDrZ3anilH8Dzzxfa+3f59Pvv9st96CGoVCno2hRf4S2BwZK2TZwY2Vp31y77dnbttXDqqXs8dPf/rxtusCOsWzcbGgxPIZ4zxwJMWJ8+9nLz5ln5gQds98inn7byL79Yhug2bQr+NmPNA0iUdu2CNWssAKxYYQPWa9fat44lS2zG0o4d9ocR3gI1LBwgGja0Xepq1bLtS6+6yo5q1Wwnug4d7BywRkJYjRpw5ZX7qODGjTb4WqqUNU3697eR8nLlrOILF1ply5e3HOguPvTvb0v1b7wx6JqUHGXK2P9G+FvcsmUWyJcts/LSpTbjZMAAy86ZkWHnhhM+7qZNm+wf/o8+mv3+m2+272xha9faEXb77dY4Ck9qOftsa+0MGWLlF16wXogLLrDynDnWYK1ePcr3XxA59WsV1yPaMZChQ1WPPFK1eXPVli0tV1PWsYesR/nyqk2a2PjEYYep3nab7WLXsaN1aael2eru77+P9KEW2IoVqm++aRnsVFU//NAqE95J6YsvbHAkLa2QXtDFxNSp9nu7//6ga+JUI4MYU6eqnnBCZFT9q69Uy5SJZIZcsED1gw9C6RgKbtYs1T/+iJT79bNMEGGHH656ySWRcoMGtmtn2IknWraHsAcftM3DCgIfRI8+gNxzj/29VK+ueu65qr16qbZpYxln33/fgsHzz9tlWI6J4KK1ebONhC9YYOXZs1VTUiIpHH7+2X6Vo0dbef581Yce8oCRaC66SLVSJZ+wEO9mzrTsoatWWfmVV+z/LzyN6/337VM8PFNl9uzIjJZCkJmZ/fNl1CjVn36KlC+/3HYHDp9bsWKOOyDnS1wGEOAMYDa2bW3fHO4X4LnQ/VOBNnl9bE5HQWZhFcrvPus+4D/9FNlce9s2m14V3ht89WrVevVUX33Vymlp9qsKzyJZtszSe4cTTm3bZn+kJWleYXEzY4bNrb7rrqBr4vJryxZLGBf+//7gA2uxhD807rvPfrfh/89HHrFmRLiF88knkX2WVe1/efLkSLmA07kyMyPbI0Qr3wEEGAM0zO3+gh7YNrZ/AQcT2RO92W7ndAY+DwWSo4Ff8/rYnI4CTeOdMCH7L3XYMPtWEXbPParvvBMpn3226lNPRcoHHGDfWlTtN1qqVOTDIiNDtWxZ1cces3J6un2NCD9/RoYt4AjnunbFzyWXqFaoEPlW64qP5cstJX/YiBGq114bKd90k2qdOpHyFVfYnPywSy9VbdEiUr7nHtUrr4yUBw3Knmxz1CibWxz2yy8WlAogmgByETAHuBtIzu28aA/gGODLLOU7gTt3O+cVoEeW8mygVl4em9NRoADSrJnqBRdEyo0bq158caTcooX9IYSdd57qM89Eyv36qY4ZEymPGxdpgagWr0njLn/mzFFNSrIBM1cyZd3fZPp01W+/jZSHD8++F8zdd9uua2EXXqh68smRcseO1gIKO/ZY6+EogNwCyF5zYYlIBeC+UHfRMCAzy+D7wH2Mz++ViFwAnKGq14TKlwHtVfXGLOd8CgxQ1R9C5a+BPkDDfT02y3P0AnoB1K9fv+2iRYuiq/CUKTarKTyXfNUqm/hdsWJ0z+dc2FVX2ZTqBQs8d4YruPXrbaZYeFrW9On2OdWgQdRPGW0urHRgC1AWqESWAFIIctoiaPdolts5eXms3ag6GBgMlkwxPxXMpnXr7OXwnHHnCmLBAnjrLVvI5sHDFYYDDshebt48Zi+VawARkTOAgcBobPB6ayG/dhqQZb0ndYGleTynTB4e61z8GzDA1urccUfQNXEu3/a23crdwIWq2jcGwQNgItBERBqJSBmgOxasshoNXC7maGCDqi7L42Odi2+LF1uei6uvLh55LVyJk2sLRFVPiOULq2qGiNwIfInNqhqiqtNF5LrQ/S9jM8E6Y1N1twJX7u2xsayvc4Xuscfssm/fYOvhXJR8QynngrB0qSWsvOwyS3TmXBzLbRDddwx2LghPPGEzZe68M+iaOBc1DyDOFbUVK+Dlly3VcjhtvnMJyAOIc0XtiScshfhddwVdE+cKxAOIc0Vp5UoYNMj2mzj00KBr41yBeABxrig99ZRtIXn33UHXxLkC8wDiXFFZvRpefBG6d4fDDw+6Ns4VmAcQ54rKwIG2C+Q99wRdE+cKhQcQ54rCmjW21/yFF0KzZkHXxrlC4QHEuaLw1FOwZQvcd1/QNXGu0HgAcS7WVq+21sdFF8U0M6pzRc0DiHOxNnCgtz5cseQBxLlYCrc+unXzsQ9X7HgAcS6WnnzSWx+u2PIA4lysrFhhrY+LL4amTYOujXOFzgOIc7Hy2GOwYwfcf3/QNXEuJjyAOBcLS5dazqvLL4cmTYKujXMxEUgAEZGqIjJOROaGLqvkcE49EflWRGaKyHQRuTnLff1EZImITAkdnYv2HTi3D488Yvt93Htv0DVxLmaCaoH0Bb5W1SbA16Hy7jKAW1W1KXA00FtEsk5jeVpVW4eOMbGvsnN5tGgRDB4MV10FjRoFXRvnYiaoANIVGBq6PhQ4d/cTVHWZqk4OXd8EzATqFFUFnYtav36QlOStD1fsBRVAaqrqMrBAARy4t5NFpCFwJPBrlptvFJGpIjIkpy6wLI/tJSKTRGTSqlWrCqHqzu3FzJnw1lvQuzfUrRt0bZyLqZgFEBH5SkSm5XB0zefzVAQ+Am5R1Y2hmwcBjYHWwDLgqdwer6qDVTVFVVNq1KgR3ZtxLq/uuw/Kl4e+OfXKOle8lI7VE6vqqbndJyIrRKSWqi4TkVrAylzOS8aCxzuq+nGW516R5ZxXgU8Lr+bORSk1FT780IKIf1lxJUBQXVijgZ6h6z2BUbufICICvA7MVNWBu91XK0vxPGBajOrpXN7deSdUrQq33hp0TZwrEkEFkAHAaSIyFzgtVEZEaotIeEbVccBlwMk5TNd9XET+FJGpwEnAf4u4/s5lN26cHffcA/vvH3RtnCsSoqpB16HIpKSk6KRJk4KuhituMjMhJQXWroXZs6Fs2aBr5FyhEpFUVU3Z/faYjYE4V2K89x78/jsMG+bBw5UonsrEuYLYsQPuvhtatbKkic6VIN4Cca4gXnwRFiyAL76wxYPOlSD+F+9ctNasgf794fTT7XCuhPEA4ly0+veHjRtt0yjnSiAPIM5FY+5c6766+mpo0SLo2jgXCA8gzkXj9tttxtWDDwZdE+cC44PozuXXuHEwahQ8+igcdFDQtXEuMN4CcS4/0tPh5puhcWP4rydAcCWbt0Ccy4+XXrKU7aNG+aJBV+J5C8S5vFq5Eu6/Hzp1gi5dgq6Nc4HzAOJcXvXpA1u2wDPPgEjQtXEucB5AnMuLCRPgzTfhttugadOga+NcXPAA4ty+pKfDDTdA/fqWrt05B/ggunP79vzzMG0ajBwJFSoEXRvn4oa3QJzbm4UL4d574ayzoGvXoGvjXFwJJICISFURGScic0OXVXI5b2Fo58EpIjIpv493rkBU4brrLMvuSy/5wLlzuwmqBdIX+FpVmwBfh8q5OUlVW++2G1Z+Hu9cdN59F778Eh55xMY/nHPZBBVAugJDQ9eHAucW8eOd27vVq+GWW6B9extAd87tIagAUlNVlwGELg/M5TwFxopIqoj0iuLxiEgvEZkkIpNWrVpVSNV3xV7v3rBhA7z6KpQqFXRtnItLMZuFJSJfATllmrs7H09znKouFZEDgXEiMktVx+enHqo6GBgMkJKSovl5rCuh3n/fjocfhpYtg66Nc3ErZgFEVU/N7T4RWSEitVR1mYjUAlbm8hxLQ5crRWQk0A4YD+Tp8c7l24oV1mWVkgJ33BF0bZyLa0F1YY0Geoau9wRG7X6CiFQQkUrh60AnYFpeH+9cvoVnXW3eDEOHQmlfJuXc3gQVQAYAp4nIXOC0UBkRqS0iY0Ln1AR+EJE/gN+Az1T1i7093rkCef11+OQTeOghaNYs6No4F/dEteQMC6SkpOikSZP2faIreWbPhjZt4JhjYOxYW/vhnANARFJ3W0oB+Ep052DnTrjkEihXzrquPHg4lyfeyevcnXdCaip8/DHUqRN0bZxLGP5Vy5VsI0fCwIG27uO884KujXMJxQOIK7nmz4crr7Qpu089FXRtnEs4HkBcybRtG1x0kSVIfP9939/cuSj4GIgreVShVy+YPBlGjYJGjYKukXMJyQOIK3mefhrefhv694cuXYKujXMJy7uwXMkydizcfjucfz7cdVfQtXEuoXkAcSXHn3/ChRdC8+a+3sO5QuD/Qa5kWLrUtqWtWBE++8wunXMF4mMgrvjbuNHGOtauhQkToF69oGvkXLHgAcQVb9u3w7nnwh9/2IyrI48MukbOFRseQFzxlZEB3bvDt9/arKuzzgq6Rs4VKz4G4oqnXbvgiius1fHCC5Ys0TlXqDyAuOJn1y7o2RPeeQceecTyXDnnCl0gAUREqorIOBGZG7qsksM5h4nIlCzHRhG5JXRfPxFZkuW+zkX+Jlx8ysiIBI+HH7ZMu865mAiqBdIX+FpVmwBfh8rZqOpsVW2tqq2BtsBWYGSWU54O36+qY3Z/vCuBtm+3dR7hlocvFHQupoIKIF2BoaHrQ4Fz93H+KcBfqroolpVyCWzTJhsk/+QTePZZb3k4VwSCCiA1VXUZQOjywH2c3x0YvtttN4rIVBEZklMXmCtBliyBDh3g++9h2DC46aaga+RciRCzACIiX4nItByOrvl8njLAOcAHWW4eBDQGWgPLgFw3cxCRXiIySUQmrVq1Kv9vxMW3KVOgfXuYOxdGj4ZLLw26Rs6VGDFbB6Kqp+Z2n4isEJFaqrpMRGoBK/fyVGcCk1V1RZbn/ue6iLwKfLqXegwGBgOkpKRoPt6Ci3cffwyXXw5VqsAPP0CrVkHXyLkSJagurNFAz9D1nsCovZzbg926r0JBJ+w8YFqh1s7Ft4wM6NsX/vUvS4z4668ePJwLQFABZABwmojMBU4LlRGR2iLyz4wqESkfuv/j3R7/uIj8KSJTgZOA/xZNtV3gli2D00+Hxx6Df/8bxo+H2rWDrpVzJVIgqUxUdQ02s2r325cCnbOUtwLVcjjvsphW0MWnzz6z1eVbtsDrr8NVVwVdI+dKNF+J7uLfpk1w/fVw9tlQpw6kpnrwcC4OeABx8W3cOGjRAl55Bf7v/+CXX6Bp06Br5ZzDA4iLV8uWwcUXQ6dOUK6czbJ66im77pyLCx5AXHzZscMCxeGHw0cfwX332VqPY48NumbOud34fiAuPmRmwocfWgqS+fPhzDMtJUmTJkHXzDmXC2+BuGCp2p4dbdpAt25QoQJ8+SWMGePBw7k45wHEBSMjA4YPt8Bx7rk2NXfYMPj9dxv3cM7FPe/CckVr7Vp44w3bJXDhQhvreOMNy2FV2v8cnUsk/h/rYk8VfvrJFv+NGAHbtsEJJ8Azz0CXLpDkDWHnEpEHEBc78+ZZN9Xbb8OcOVCxIlx2mW0xe8QRQdfOOVdAHkBc4VGFWbMsS+7IkbZiHODEE2121QUXWBBxzhULHkBcwWzebAkNP//cZk7Nn2+3t2sHjz8O3btDvXrB1tE5FxMeQFz+bNhg4xk//ADffQe//WYzqvbbD04+GW69Fbp2tZxVzrlizQOIy93WrTB1KkyeDBMn2r4bs2ZZV1WpUtC2Ldx2G5xyChx/vKcZca6E8QDiLNvtvHkWHGbOhOnT4c8/7TYNbeJYvbptHdujBxx3nF2vUCHYejvnAuUBpCTYtAmWLIHFi+Hvv2HRIliwwI6//oLlyyPnJiXBIYfYLKlLLoHWrW2xX926IBLYW3DOxZ9AAoiIXAj0A5oC7VR1Ui7nnQE8C5QCXlPV8M6FVYH3gIbAQuAiVV0X84rHA1Vbtb1unS3KW7sWVq+GNWtg1SpYudKO5cvtWLrUBrqzSkqygNCoEXTubClDDjnE0qQfcgiULRvMe3POJZSgWiDTgPOBV3I7QURKAS9iW9qmARNFZLSqzgD6Al+r6gAR6Rsq94l9tfdh1y5IT7eMsjt32mXWY/v2yLFtmx1bt9qxZUv2Y/Nmazls3Bg5NmywY9eu3OtwwAFw4IFw0EFw5JGWlLBOHTvq1bOjTh1ITi6yH4tzrngKakvbmQCy9y6RdsA8VZ0fOncE0BWYEbrsGDpvKPAdsQwg/fvDu+/abKP0dDuyXg8fmZkFe52kJFsnUbGijS9UqmRH/fqw//5QubIFiPBRrRpUqWLjE9WrW7lMmUJ4w845t2/xPAZSB1icpZwGtA9dr6mqywBUdZmIHJjbk4hIL6AXQP369aOrSa1aNiZQurQdycmRy/BRpkzketmyVi5bNnKUK2dH2bI25bVcOShf3q6XL28Bo0wZH2dwziWMmAUQEfkKOCiHu+5W1VF5eYocbtP81kNVBwODAVJSUvL9eACuucYO55xz/4hZAFHVUwv4FGlA1iXMdYGloesrRKRWqPVRC1hZwNdyzjmXT/GcBnUi0EREGolIGaA7MDp032igZ+h6TyAvLRrnnHOFKJAAIiLniUgacAzwmYh8Gbq9toiMAVDVDOBG4EtgJvC+qk4PPcUA4DQRmYvN0hpQ1O/BOedKOlGNblggEaWkpOikSTkuOXHOOZcLEUlV1ZTdb4/nLiznnHNxzAOIc865qHgAcc45FxUPIM4556JSogbRRWQVsCgGT10dWB2D5y0qiV5/SPz3kOj1h8R/D4lef4jde2igqjV2v7FEBZBYEZFJOc1QSBSJXn9I/PeQ6PWHxH8PiV5/KPr34F1YzjnnouIBxDnnXFQ8gBSOwUFXoIASvf6Q+O8h0esPif8eEr3+UMTvwcdAnHPORcVbIM4556LiAcQ551xUPIAUEhHpLyJTRWSKiIwVkdpB1yk/ROQJEZkVeg8jReSAoOuUXyJyoYhMF5FMEUmY6ZgicoaIzBaReSLSN+j65JeIDBGRlSIyLei6RENE6onItyIyM/T3c3PQdcoPESknIr+JyB+h+j9QZK/tYyCFQ0T2V9WNoes3Ac1U9bqAq5VnItIJ+EZVM0TkMQBVjd0+8zEgIk2BTOAV4DZVjfvUyyJSCpiDbUuQhu2D00NVZwRasXwQkROBzcBbqtoi6PrkV2hTulqqOllEKgGpwLmJ8jsQEQEqqOpmEUkGfgBuVtVfYv3a3gIpJOHgEVKBKLbfDZKqjg3twQLwC7YDZEJR1ZmqOjvoeuRTO2Ceqs5X1Z3ACKBrwHXKF1UdD6wNuh7RUtVlqjo5dH0Ttv9QnWBrlXdqNoeKyaGjSD5/PIAUIhF5WEQWA5cA9wVdnwK4Cvg86EqUEHWAxVnKaSTQh1dxIyINgSOBXwOuSr6ISCkRmYJt7z1OVYuk/h5A8kFEvhKRaTkcXQFU9W5VrQe8g+2mGFf2Vf/QOXcDGdh7iDt5eQ8JRnK4LaFar8WFiFQEPgJu2a1HIe6p6i5VbY31HLQTkSLpSixdFC9SXKjqqXk89V3gM+D+GFYn3/ZVfxHpCZwNnKJxOjiWj99BokgD6mUp1wWWBlSXEis0dvAR8I6qfhx0faKlqutF5DvgDCDmkxq8BVJIRKRJluI5wKyg6hINETkD6AOco6pbg65PCTIRaCIijUSkDNAdGB1wnUqU0CD068BMVR0YdH3yS0RqhGdNish+wKkU0eePz8IqJCLyEXAYNgtoEXCdqi4JtlZ5JyLzgLLAmtBNvyTSLDIAETkPeB6oAawHpqjq6YFWKg9EpDPwDFAKGKKqDwdbo/wRkeFARyyV+ArgflV9PdBK5YOIHA9MAP7E/n8B7lLVMcHVKu9E5AhgKPb3kwS8r6oPFslrewBxzjkXDe/Ccs45FxUPIM4556LiAcQ551xUPIA455yLigcQ55xzUfEA4lxAQllgF4hI1VC5SqjcIOi6OZcXHkCcC4iqLgYGAQNCNw0ABqvqouBq5Vze+ToQ5wIUSqGRCgwBrgWODGXldS7ueS4s5wKkqukicjvwBdDJg4dLJN6F5VzwzgSWAQm3GZMr2TyAOBcgEWmN7UZ4NPDf0O54ziUEDyDOBSSUBXYQtv/E38ATwJPB1sq5vPMA4lxwrgX+VtVxofJLwOEi0iHAOjmXZz4LyznnXFS8BeKccy4qHkCcc85FxQOIc865qHgAcc45FxUPIM4556LiAcQ551xUPIA455yLyv8DBW445My/d8cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "in_array = np.linspace(-np.pi, np.pi, 120) \n",
    "\n",
    "tan_h = np.tanh(in_array) \n",
    "sig = sigmoid(in_array)  \n",
    " \n",
    "plt.plot(in_array, tan_h, color = 'red', label='tanH')\n",
    "plt.plot(in_array, 1-tan_h**2, linestyle=':', color = 'red', label=\"tanH'\")\n",
    "plt.plot(in_array, sig, color = 'blue', label='sigmoid')\n",
    "plt.plot(in_array, sig*(1-sig), linestyle=':', color = 'blue',  label=\"sigmoid'\")\n",
    "plt.xlabel(\"X\") \n",
    "plt.ylabel(\"Y\") \n",
    "plt.legend(loc='best')\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LSTM peculiarities\n",
    "\n",
    "Having looked into the computations within a LSTM cell by writing numpy code in the previous part, our next task is to move on to Keras and build proper LSTMs. It turns out that this task is not quite as easy as when defining and training a feedforward network. Sequence data involves some additional challenges. The implementation of the LSTM in Keras also exhibits some peculiarities that we have to be aware of. \n",
    "\n",
    "The purpose of this part is to review some of the new challenges mentioned above and to provide some introductory content to help get you started. We will only briefly touch on selected concepts and revisit these in the next tutorial on [financial time series forecasting using an LSTM](https://colab.research.google.com/github/Humboldt-WI/adams/blob/master/demos/rnn/lstm_fin_forecasting.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Peculiarities #1: input/output\n",
    "\n",
    "LSTMs are flexible in terms of the structure of the input and the output sequence. This picture from [A. Karpathy's blog](http://karpathy.github.io/2015/05/21/rnn-effectiveness/) reminds us of the possibilities. \n",
    "\n",
    "<img src=\"RNN_structures.jpg\" style=\"width: 800px;\"/>\n",
    "Source: http://karpathy.github.io/2015/05/21/rnn-effectiveness/\n",
    "\n",
    "Applications like:\n",
    "1. Image recognition (binary outcome like cat/dog) - strictly speaking its not a sequence\n",
    "2. Image captioning (outcome \"Tall man in red shirt holds beer) - not really for Keras\n",
    "3. Sentiment analysis, typing suggestions or our time series prediction (using a bunch of past values to predict one in the future)\n",
    "4. Translation or also can be prediction of several steps in time series\n",
    "5. May be used for example for ongoing video classification\n",
    "\n",
    "Additionally, we have an architecture called **Bidirectional LSTM** that not only preserves information from the past and passes it to the future but also passes information from the future to the past (used widely for the tasks where context matters a lot)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Peculiarities #2: train set structure\n",
    "In the manual implementation, we have touched on how the data can be structured before feeding it into the model. Ultimately, we always want to have a set of features and a target. Although LSTMs are designed as models for sequential data, it turns out that they do not work well if the length of the input sequence is too long. As a rule of thumb, many sources mention a range of 200 - 400 time steps, which an LSTM can still process. When saying \"can\", you want to note that the quality of a predictive model can always be evaluated by measuring forecast accuracy. In addition to decreased accuracy, issues with too long sequences also concern the speed of training and the vanishing gradient problem. It is common practice to split a long sequence into subsequences and using the resulting \"chunks of data\" as input for the LSTM. This is the approach we adopted in the last tutorial. Remember this picture:\n",
    "\n",
    "<img src=\"time_series_one_step.png\" style=\"width: 500px;\"/>\n",
    "\n",
    "An important questions arises: what should be the length of the subsequence or time window? It will depend a lot on field knowledge. For example, in the case of sentiment analysis, the question is how many words you need on average to convey a sentiment. With ordinary time series, you can examine auto-correlation or partial auto-correlation to derive a judgment.\n",
    "\n",
    "Next, relating the previous picture to the picture on possible LSTM setups, we note that we could also use a LSTM to forecast multiple steps ahead and not just one. This would look as follows:  \n",
    "\n",
    "<img src=\"time_series_sequence.png\" style=\"width: 500px;\"/>\n",
    "\n",
    "We decide on the sequence length and keep it constant. Our input and output sequences are now of the same length. Each new sequence starts one time step ahead of the sample sequence. The remaining values, i.e. the remainder of the full sequence length divided by the chosen length, are dropped. Keep in mind that then your model turns into seq2seq type. By making the network stateful, we use the hidden state at the end of the previous sequence as the starting point of the next sequence. Since the hidden state summarizes information from earlier observation, this increases the information available for each prediction beyond the window size. \n",
    "\n",
    "Last, note that a sequence learning task could also have an entirely different form of time series or text data. To illustrate this, let's say we are expecting the model to figure out that it has to count the number of ones in an input sequence and give it as prediction."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "     X       y\n",
    "[0,1,1,1,2] [3]\n",
    "[0,0,1,1,2] [2]\n",
    "[0,0,0,1,2] [1]\n",
    "[1,0,0,1,2] [2]\n",
    "[0,0,0,0,1] [0]\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the order of those observations does not matter, although it is a sequence problem. We could shuffle the observations as we often do when working with non-sequential data, hoping that shuffling gives better generalization.\n",
    "\n",
    "In any case, remember that Keras expects the input data for an LSTM to be a 3d tensor with dimensions:\n",
    "- Number of samples or batch size\n",
    "- Number of time steps\n",
    "- Number of features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Peculiarities #3: state\n",
    "Keras supports two types of LSTMs: **stateful** or **stateless**\n",
    "\n",
    "These terms are somewhat confusing. Isn't the very idea of RNNs and latent variable models to keep a hidden state in which information from processing past time steps accumulates? Of course it is. In Keras, the difference between a stateful and a stateless LSTM concerns the point when the hidden state is reset. The maybe more common case is a stateless LSTM. There, the hidden state is reset after processing one batch (see below for peculiarities wrt batches). \n",
    "\n",
    "\n",
    "Recall that longer time series are commonly split into chunks. This resetting approach enforces what we said above. The time window or length of a subsequence in which a longer time series is split is crucial. A *stateless* LSTM will not maintain state information across multiple batches. The information that you deem relevant to forecast must be included in the subsequence. \n",
    "\n",
    "*Stateful* LSTMs are Keras's solution toward maintaining state when working with longer sequences. Conceptually, *stateful* LSTM are conceptually more suitable for such settings. The splitting of the time series into chunks is somewhat arbitrary and done because not doing it would probably break LSTM training. You would not bother with splitting a time series into chunks when using an econometric model like ARIMA.   \n",
    "\n",
    "Before demonstrating *stateful* and *stateless* LSTMs in the [next demo](https://colab.research.google.com/github/Humboldt-WI/adams/blob/master/demos/rnn/lstm_fin_forecasting.ipynb), let's note that a stateful recurrent model is one for which the internal states (memories) obtained after processing a batch of samples are reused as initial states for the samples of the next batch. This allows processing longer sequences while keeping computational complexity manageable. Again, the hidden state is passed from batch to batch, not within the batch. Within a batch the sub-sequences are treated as independent.\n",
    "\n",
    "If the model is stateless, the cell states are reset after each batch. It is considered to be more efficient in implementations than stateful and the clear choice when every observation does not depend on the previous one (e.g., sentence classification). Additionally, stateless LSTMs can be implemented with and without shuffling of observations.\n",
    "\n",
    "\n",
    "This [blog post](http://philipperemy.github.io/keras-stateful-lstm/) offers a very good explanation of stateful versus stateless and why the two are distinguished at all. It also offers a nice demo. Same more demos are available in [this tutorial](https://machinelearningmastery.com/stateful-stateless-lstm-time-series-forecasting-python/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Peculiarities #4: batches \n",
    "\n",
    "When it comes to time series forecasting, there are many suggestions what should be the batch size for an LSTM. One of the rules one can observe is your batch size is the same as your output (according to Keras functionality). In case your output is next day prediction, you might want to consider batch size=1. More specifically, the Keras implementation enforces a constant batch size for training and prediction (see, e.g., [here for a discussion](https://stackoverflow.com/questions/43702481/why-does-keras-lstm-batch-size-used-for-prediction-have-to-be-the-same-as-fittin)). Thus, if you want to make prediction one day into the future, a common setting in stock prediction settings, then you should train your LSTM with batch size equal to one. Conversely, if you train your LSTM with batch size equal to, say, 10, you must provide exactly ten samples to the *.predict()* function. In practice, the latter implies that you would need to wait for ten days before you can make a prediction. I know this sounds very serious. Relax, there are ways to get round this issue. One is to use Tensorflow but you can [also *persuade* Keras to let you change the batch size from training to prediction].(https://machinelearningmastery.com/use-different-batch-sizes-training-predicting-python-keras/)  \n",
    "\n",
    "In general, the advice is to keep batch size a value that divides without remainder into the train and validation set sizes, so that no data gets discarded. Keep in mind that the LSTM will be chaining together line 1 of batch n with line 1 of batch n+1 (see below). So try to keep the batch size somewhat meaningful. \n",
    "\n",
    "However, when it comes to text processing, the situation changes. [Jeremy Howard](https://www.youtube.com/watch?v=H3g26EVADgY&feature=youtu.be&t=17m50sin) advocates that splitting a big string of text into chunks, stacking it and then creating batches from the first \"slice\", second \"slice\" and etc. offers nice parallelization capacities and does not really harm the training process. Keep in mind that in case you have short sequences (that are still longer than your window size) this method might bring in additional distortion.\n",
    "\n",
    "Example (mind that numbers are only used for simplicity, should rather be words: we have a \"long\" string [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17......1000]\n",
    "Let's say we want to have batch=10, then we split the string into 10 chunks (that would make 100 numbers per chunk) and stack them:\n",
    "\n",
    "[1.....100]\n",
    "\n",
    "[101...200]\n",
    "\n",
    "[201...300]\n",
    "\n",
    "...\n",
    "\n",
    "[900...1000]\n",
    "\n",
    "(10 rows, 100 columns)\n",
    "\n",
    "Now, if you remember in order to transform it into a supervised task we need to decide on the \"window\" size and the output type. Let's say we will go for a many-to-many architecture, so we choose the window=10 and same size output without overlap. \n",
    "Then our stack would look like this:\n",
    "\n",
    "[1..10][11...20][21...30]......[91..100] \n",
    "\n",
    "[101..110][111...120][121...130]......[191..200]\n",
    "\n",
    ".........\n",
    "\n",
    "[901..910][911...920][921...930]......[991..1000]\n",
    "\n",
    "That is [1..10] will be used to predict [2...11], [11...20] will be used to predict [12...21] and so forth.\n",
    "\n",
    "\n",
    "Then, the first batch would be the first slice/column:\n",
    "\n",
    "[1..10]\n",
    "\n",
    "[101..110]\n",
    "\n",
    "....\n",
    "\n",
    "[901..910]\n",
    "\n",
    "\n",
    "Second batch - second column and etc. where we use a stateful network to continue with the hidden state from the previous batch.\n",
    "\n",
    "Structuring the data in this way, the sequence is lost at the end points of each row (i.e. 100 to 101, 200 to 201). If we consider the sequence to be not hundreds but thousands of words, that may be something we are willing to sacrifice, because we gain a parallelization of the process: the 10 chunks/rows are training simultaneously with every batch. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Peculiarities #5: Resetting ###\n",
    "In case you are using a stateful LSTM, make sure you reset the state after every epoch, otherwise the NN will treat it as a continuation of the time series.\n"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "epochs=30\n",
    "\n",
    "for i in range(epoch):\n",
    "        model.fit(X, y, epochs=1, batch_size=n_batch, verbose=1, shuffle=False)\n",
    "        model.reset_states()\n",
    "        \n",
    "# Also you might want to consider the online forecast\n",
    "for i in range(len(X)):\n",
    "        testX, testy = X[i], y[i]\n",
    "        testX = testX.reshape(1, 1, 1)\n",
    "        yhat = model.predict(testX, batch_size=1)\n",
    "        print('>Expected=%.1f, Predicted=%.1f' % (testy, yhat))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.11"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}