notebooks/rnn/rnn.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "# Based on this paper - https://apps.dtic.mil/dtic/tr/fulltext/u2/a164453.pdf"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "import math\n",
    "\n",
    "data = pd.read_csv(\"../../data/clean_weather.csv\")\n",
    "data = data.ffill()"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [],
   "source": [
    "PREDICTORS = [\"tmax\", \"tmin\", \"rain\"]\n",
    "TARGET = \"tmax_tomorrow\"\n",
    "\n",
    "scaler = StandardScaler()\n",
    "data[PREDICTORS] = scaler.fit_transform(data[PREDICTORS])\n",
    "\n",
    "np.random.seed(0)\n",
    "split_data = np.split(data, [int(.7*len(data)), int(.85*len(data))])\n",
    "(train_x, train_y), (valid_x, valid_y), (test_x, test_y) = [[d[PREDICTORS].to_numpy(), d[[TARGET]].to_numpy()] for d in split_data]"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "outputs": [],
   "source": [
    "# Rnn\n",
    "# Input -> hidden\n",
    "# hidden -> hidden\n",
    "# hidden -> output"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "outputs": [],
   "source": [
    "def init_params(layer_conf):\n",
    "    layers = []\n",
    "    for i in range(1, len(layer_conf)):\n",
    "        np.random.seed(0)\n",
    "        k = 1/math.sqrt(layer_conf[i][\"hidden\"])\n",
    "        i_weight = np.random.rand(layer_conf[i-1][\"units\"], layer_conf[i][\"hidden\"]) * 2 * k - k\n",
    "\n",
    "        h_weight = np.random.rand(layer_conf[i][\"hidden\"], layer_conf[i][\"hidden\"]) * 2 * k - k\n",
    "        h_bias = np.random.rand(1, layer_conf[i][\"hidden\"]) * 2 * k - k\n",
    "\n",
    "        o_weight = np.random.rand(layer_conf[i][\"hidden\"], layer_conf[i][\"output\"]) * 2 * k - k\n",
    "        o_bias = np.random.rand(1, layer_conf[i][\"output\"]) * 2 * k - k\n",
    "\n",
    "        layers.append(\n",
    "            [i_weight, h_weight, h_bias, o_weight, o_bias]\n",
    "        )\n",
    "    return layers"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "outputs": [],
   "source": [
    "def forward(x, layers):\n",
    "    hiddens = []\n",
    "    outputs = []\n",
    "    for i in range(len(layers)):\n",
    "        i_weight, h_weight, h_bias, o_weight, o_bias = layers[i]\n",
    "        hidden = np.zeros((x.shape[0], i_weight.shape[1]))\n",
    "        output = np.zeros((x.shape[0], o_weight.shape[1]))\n",
    "        for j in range(x.shape[0]):\n",
    "            input_x = x[j,:][np.newaxis,:] @ i_weight\n",
    "            hidden_x = input_x + hidden[max(j-1,0),:][np.newaxis,:] @ h_weight + h_bias\n",
    "            # Activation.  tanh avoids outputs getting larger and larger.\n",
    "            hidden_x = np.tanh(hidden_x)\n",
    "            # Store hidden for use in backprop\n",
    "            hidden[j,:] = hidden_x\n",
    "\n",
    "            # Output layer\n",
    "            output_x = hidden_x @ o_weight + o_bias\n",
    "            output[j,:] = output_x\n",
    "        hiddens.append(hidden)\n",
    "        outputs.append(output)\n",
    "    return hiddens, outputs[-1]"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "outputs": [],
   "source": [
    "def mse(actual, predicted):\n",
    "    return np.mean((actual-predicted)**2)\n",
    "\n",
    "def mse_grad(actual, predicted):\n",
    "    return (predicted - actual)"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "outputs": [],
   "source": [
    "def backward(layers, x, lr, grad, hiddens):\n",
    "    for i in range(len(layers)):\n",
    "        i_weight, h_weight, h_bias, o_weight, o_bias = layers[i]\n",
    "        hidden = hiddens[i]\n",
    "        next_h_grad = None\n",
    "        i_weight_grad, h_weight_grad, h_bias_grad, o_weight_grad, o_bias_grad = [0] * 5\n",
    "\n",
    "        for j in range(x.shape[0] - 1, -1, -1):\n",
    "            # 1,1\n",
    "            # Add newaxis in the first dimension\n",
    "            out_grad = grad[j,:][np.newaxis, :]\n",
    "\n",
    "            # Output updates\n",
    "            # 3x1 @ 1x1 = 3,1\n",
    "            # np.newaxis creates a size 1 axis, in this case transposing matrix\n",
    "            o_weight_grad += hidden[j,:][:, np.newaxis] @ out_grad\n",
    "            # 1,1\n",
    "            o_bias_grad += out_grad\n",
    "\n",
    "            # Propagate gradient to hidden unit\n",
    "            # 1,1 @ 1,3 = 1,3\n",
    "            h_grad = out_grad @ o_weight.T\n",
    "\n",
    "            if j < x.shape[0] - 1:\n",
    "                # Then we multiply the gradient by the hidden weights to pull gradient from next hidden state to current hidden state\n",
    "                hh_grad = next_h_grad @ h_weight.T\n",
    "                # Add the gradients together to combine output contribution and hidden contribution\n",
    "                h_grad += hh_grad\n",
    "\n",
    "            # Pull the gradient across the current hidden nonlinearity\n",
    "            # derivative of tanh is 1 - tanh(x) ** 2\n",
    "            # So we take the output of tanh (next hidden state), and plug in\n",
    "            tanh_deriv = 1 - hidden[j][np.newaxis,:] ** 2\n",
    "\n",
    "            # next_h_grad @ np.diag(tanh_deriv_next) multiplies each element of next_h_grad by the deriv\n",
    "            # Effect is to pull value across nonlinearity\n",
    "            h_grad = np.multiply(h_grad, tanh_deriv)\n",
    "\n",
    "            # Store to compute h grad for previous sequence position\n",
    "            next_h_grad = h_grad.copy()\n",
    "\n",
    "            # If we're not at the very beginning\n",
    "            if j > 0:\n",
    "                # 3,1 @ 1,3\n",
    "                # Multiply input from previous layer by post-nonlinearity grad at current layer\n",
    "                h_weight_grad += hidden[j-1][:, np.newaxis] @ h_grad\n",
    "                h_bias_grad += h_grad\n",
    "\n",
    "            # 3,1 (with newaxis) @ h_grad\n",
    "            i_weight_grad += x[j,:][:,np.newaxis] @ h_grad\n",
    "\n",
    "        lr = lr / x.shape[0]\n",
    "        i_weight -= i_weight_grad * lr\n",
    "        h_weight -= h_weight_grad * lr\n",
    "        h_bias -= h_bias_grad * lr\n",
    "        o_weight -= o_weight_grad * lr\n",
    "        o_bias -= o_bias_grad * lr\n",
    "        layers[i] = [i_weight, h_weight, h_bias, o_weight, o_bias]\n",
    "    return layers"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Epoch: 0 train loss 3122.594400144508 valid loss 2171.3186862102025\n",
      "Epoch: 10 train loss 84.7851930495282 valid loss 81.34087884698094\n",
      "Epoch: 20 train loss 52.84166669240144 valid loss 52.30941167562456\n",
      "Epoch: 30 train loss 34.89142680576195 valid loss 34.42998362512595\n",
      "Epoch: 40 train loss 33.906747974192236 valid loss 33.789660221951614\n",
      "Epoch: 50 train loss 30.593193275313595 valid loss 30.568271740103427\n",
      "Epoch: 60 train loss 27.767753755456706 valid loss 27.466770094463552\n",
      "Epoch: 70 train loss 27.987707979297518 valid loss 27.581650000896172\n",
      "Epoch: 80 train loss 28.42712118570627 valid loss 27.791959394170433\n",
      "Epoch: 90 train loss 26.598198926417233 valid loss 25.806742323792438\n",
      "Epoch: 100 train loss 25.263986813543738 valid loss 24.435517510355645\n",
      "Epoch: 110 train loss 24.445111564692038 valid loss 23.612611329688526\n",
      "Epoch: 120 train loss 23.89688420942656 valid loss 23.07347371245211\n",
      "Epoch: 130 train loss 23.49986279360287 valid loss 22.690708423306877\n",
      "Epoch: 140 train loss 23.196481286413093 valid loss 22.402625853595673\n",
      "Epoch: 150 train loss 22.9567624295313 valid loss 22.177010971976852\n",
      "Epoch: 160 train loss 22.763861662989473 valid loss 21.99587590418936\n",
      "Epoch: 170 train loss 22.607700358322987 valid loss 21.848663372586287\n",
      "Epoch: 180 train loss 22.4819605672346 valid loss 21.72898008580287\n",
      "Epoch: 190 train loss 22.3825592047859 valid loss 21.632903358642164\n",
      "Epoch: 200 train loss 22.306774327704215 valid loss 21.557992202834164\n",
      "Epoch: 210 train loss 22.2525802166852 valid loss 21.502527457734892\n",
      "Epoch: 220 train loss 22.217889244207004 valid loss 21.464647837242342\n",
      "Epoch: 230 train loss 22.199595042067486 valid loss 21.441263527739522\n",
      "Epoch: 240 train loss 22.192934673016204 valid loss 21.427339933427824\n"
     ]
    }
   ],
   "source": [
    "epochs = 250\n",
    "lr = 1e-5\n",
    "\n",
    "layer_conf = [\n",
    "    {\"type\":\"input\", \"units\": 3},\n",
    "    {\"type\": \"rnn\", \"hidden\": 4, \"output\": 1}\n",
    "]\n",
    "layers = init_params(layer_conf)\n",
    "\n",
    "for epoch in range(epochs):\n",
    "    sequence_len = 7\n",
    "    epoch_loss = 0\n",
    "    for j in range(train_x.shape[0] - sequence_len):\n",
    "        seq_x = train_x[j:(j+sequence_len),]\n",
    "        seq_y = train_y[j:(j+sequence_len),]\n",
    "        hiddens, outputs = forward(seq_x, layers)\n",
    "        grad = mse_grad(seq_y, outputs)\n",
    "        params = backward(layers, seq_x, lr, grad, hiddens)\n",
    "        epoch_loss += mse(seq_y, outputs)\n",
    "\n",
    "    if epoch % 10 == 0:\n",
    "        sequence_len = 7\n",
    "        valid_loss = 0\n",
    "        for j in range(valid_x.shape[0] - sequence_len):\n",
    "            seq_x = valid_x[j:(j+sequence_len),]\n",
    "            seq_y = valid_y[j:(j+sequence_len),]\n",
    "            _, outputs = forward(seq_x, layers)\n",
    "            valid_loss += mse(seq_y, outputs)\n",
    "\n",
    "        print(f\"Epoch: {epoch} train loss {epoch_loss / len(train_x)} valid loss {valid_loss / len(valid_x)}\")"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "from sympy import symbols, exp, diff, simplify\n",
    "\n",
    "x = symbols(\"x\")\n",
    "tanh = (exp(2 * x) - 1) / (exp(2 * x) + 1)\n",
    "simplify(diff(tanh, x))"
   ],
   "metadata": {
    "collapsed": false
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}