basin-inflow-lstm.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "9faa064e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.simplefilter(\"ignore\")\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import tensorflow as tf\n",
    "from tensorflow import keras\n",
    "import keras_tuner as kt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "93ea6530",
   "metadata": {},
   "outputs": [],
   "source": [
    "DATA_PATH = 'Reservoir_Project/Data'\n",
    "HP_TUNING_PATH = 'Reservoir_Project/Hyperparameter_Tuning'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "dd49ef72",
   "metadata": {},
   "outputs": [],
   "source": [
    "basin_inflow = pd.read_excel(f'{DATA_PATH}/Custom/basin_inflow_no_dates.xlsx', index_col=0)\n",
    "basin_inflow_train = pd.read_excel(f'{DATA_PATH}/Custom/basin_inflow_train.xlsx', index_col=0)\n",
    "basin_inflow_validation = pd.read_excel(f'{DATA_PATH}/Custom/basin_inflow_validation.xlsx', index_col=0)\n",
    "basin_inflow_test = pd.read_excel(f'{DATA_PATH}/Custom/basin_inflow_test.xlsx', index_col=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "748d978f",
   "metadata": {},
   "source": [
    "### Create data window"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d373e853",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" TensorFlow utility class for producing data windows from time series data \"\"\"\n",
    "\n",
    "class WindowGenerator():\n",
    "  def __init__(self, input_width, label_width, shift,\n",
    "               train_df, val_df, test_df,\n",
    "               label_columns=None):\n",
    "    # Store the raw data.\n",
    "    self.train_df = train_df\n",
    "    self.val_df = val_df\n",
    "    self.test_df = test_df\n",
    "    \n",
    "    self.column_indices = {name: i for i, name in\n",
    "                           enumerate(train_df.columns)}\n",
    "    \n",
    "    # Work out the label column indices.\n",
    "    self.label_columns = label_columns\n",
    "    if label_columns is not None:\n",
    "      self.label_columns_indices = {name: i for i, name in\n",
    "                                    enumerate(label_columns)}\n",
    "\n",
    "    # Work out the window parameters.\n",
    "    self.input_width = input_width\n",
    "    self.label_width = label_width\n",
    "    self.shift = shift\n",
    "\n",
    "    self.total_window_size = input_width + shift\n",
    "\n",
    "    self.input_slice = slice(0, input_width)\n",
    "    self.input_indices = np.arange(self.total_window_size)[self.input_slice]\n",
    "\n",
    "    self.label_start = self.total_window_size - self.label_width\n",
    "    self.labels_slice = slice(self.label_start, None)\n",
    "    self.label_indices = np.arange(self.total_window_size)[self.labels_slice]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "aeb825cc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<__main__.WindowGenerator at 0x299736610>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\"\"\"\n",
    "Window\n",
    "- Given 60 days of history predict 30 days ahead. A season is about 90 days in a CA Water Year.\n",
    "- Window size: 90\n",
    "\"\"\"\n",
    "\n",
    "window = WindowGenerator(input_width=60, label_width=1, shift=30,\n",
    "                     train_df=basin_inflow_train, val_df=basin_inflow_validation, \n",
    "                     test_df=basin_inflow_test, label_columns=['INFLOW'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "379f5c38",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" Create window of inputs and labels \"\"\"\n",
    "\n",
    "def split_window(self, features):\n",
    "  inputs = features[:, self.input_slice, :]\n",
    "  labels = features[:, self.labels_slice, :]\n",
    "  if self.label_columns is not None:\n",
    "    labels = tf.stack(\n",
    "        [labels[:, :, self.column_indices[name]] for name in self.label_columns],\n",
    "        axis=-1)\n",
    "\n",
    "  # set shapes after slicing\n",
    "  inputs.set_shape([None, self.input_width, None])\n",
    "  labels.set_shape([None, self.label_width, None])\n",
    "\n",
    "  return inputs, labels\n",
    "\n",
    "WindowGenerator.split_window = split_window"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3e37dcdb",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" Create dataset of sliding windows over a time series dataframe \"\"\"\n",
    "\n",
    "def make_dataset(self, data):\n",
    "  data = np.array(data, dtype=np.float32)\n",
    "  ds = tf.keras.utils.timeseries_dataset_from_array(\n",
    "      data=data,\n",
    "      targets=None,\n",
    "      sequence_length=self.total_window_size,\n",
    "      sequence_stride=1,\n",
    "      shuffle=True,\n",
    "      batch_size=32,)\n",
    "\n",
    "  # (input_window, label_window) pairs \n",
    "  ds = ds.map(self.split_window)\n",
    "\n",
    "  return ds\n",
    "\n",
    "WindowGenerator.make_dataset = make_dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "5f27daf2",
   "metadata": {},
   "outputs": [],
   "source": [
    "@property\n",
    "def train(self):\n",
    "  return self.make_dataset(self.train_df)\n",
    "\n",
    "@property\n",
    "def val(self):\n",
    "  return self.make_dataset(self.val_df)\n",
    "\n",
    "@property\n",
    "def test(self):\n",
    "  return self.make_dataset(self.test_df)\n",
    "\n",
    "WindowGenerator.train = train\n",
    "WindowGenerator.val = val\n",
    "WindowGenerator.test = test"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fccd065a",
   "metadata": {},
   "source": [
    "### Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "6b1cd773",
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_model(hp):\n",
    "  lstm = keras.Sequential()\n",
    "\n",
    "  ### Tuning hyperparameters ### \n",
    "\n",
    "  # Number of lstm layer units\n",
    "  hp_units = hp.Int('units', min_value=32, max_value=128, step=32)\n",
    "\n",
    "  # Dropout rate applied to input values\n",
    "  hp_dropout = hp.Choice(\"dropout\", [0.2, 0.3, 0.4, 0.5])\n",
    "\n",
    "  # Recurrent dropout rate applied to hidden cell states between time steps\n",
    "  hp_recurrent_dropout = hp.Choice(\"recurrent_dropout\", [0.2, 0.3, 0.4, 0.5])\n",
    "  \n",
    "  # L2 regularization \n",
    "  hp_l2_reg = hp.Choice(\"l2\", [0.001, 0.01, 0.02, 0.05])\n",
    "\n",
    "  # Optimizer learning rate\n",
    "  hp_learning_rate = hp.Choice('learning_rate', values=[0.01, 0.001, 0.0001])\n",
    "\n",
    "  ### Layers ### \n",
    "  \n",
    "  lstm.add(tf.keras.layers.LSTM(units=hp_units, dropout=hp_dropout, recurrent_dropout=hp_recurrent_dropout, \n",
    "                  kernel_regularizer=tf.keras.regularizers.l2(l2=hp_l2_reg),\n",
    "                  return_sequences=True, input_shape=[None, 39]))\n",
    "    \n",
    "  lstm.add(tf.keras.layers.LSTM(units=hp_units, dropout=hp_dropout, recurrent_dropout=hp_recurrent_dropout, \n",
    "                  kernel_regularizer=tf.keras.regularizers.l2(l2=hp_l2_reg)))\n",
    "    \n",
    "  lstm.add(keras.layers.Dense(1))    \n",
    "    \n",
    "  ### Compile ### \n",
    "           \n",
    "  lstm.compile(loss=tf.keras.losses.MeanAbsoluteError(),\n",
    "               optimizer=tf.keras.optimizers.legacy.Adam(learning_rate=hp_learning_rate),\n",
    "               metrics=['mean_absolute_error'])\n",
    "\n",
    "  return lstm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "7a1aa7cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "tuner = kt.Hyperband(build_model,\n",
    "                     objective='val_mean_absolute_error',\n",
    "                     max_epochs=10,\n",
    "                     factor=3,\n",
    "                     directory=HP_TUNING_PATH,\n",
    "                     project_name='reservoir_model_hp_tuning')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "5cb0f6f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "MAX_EPOCHS = 10\n",
    "early_stopping = tf.keras.callbacks.EarlyStopping(monitor='val_loss', patience=5, mode='min')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70e57c45",
   "metadata": {},
   "outputs": [],
   "source": [
    "tuner.search(window.train, epochs=MAX_EPOCHS, validation_data=window.val, callbacks=[early_stopping])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "3562e549",
   "metadata": {},
   "outputs": [],
   "source": [
    "best_hps=tuner.get_best_hyperparameters(num_trials=1)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a42a5d8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# find the optimal number of training epochs\n",
    "\n",
    "model = tuner.hypermodel.build(best_hps)\n",
    "history = model.fit(window.train, epochs=MAX_EPOCHS, validation_data=window.val, callbacks=[early_stopping])\n",
    "\n",
    "val_mae_per_epoch = history.history['val_mean_absolute_error']\n",
    "best_epoch = val_mae_per_epoch.index(min(val_mae_per_epoch)) + 1 # epoch of lowest validation MAE "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcd1da1b",
   "metadata": {},
   "outputs": [],
   "source": [
    "hypermodel = tuner.hypermodel.build(best_hps)\n",
    "history = hypermodel.fit(window.train, epochs=best_epoch, validation_data=window.val, callbacks=[early_stopping])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "790b60f1",
   "metadata": {},
   "source": [
    "### Predict next 30 days of inflow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "055cf12c",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" Inverse scaling for cfs prediction \"\"\"\n",
    "\n",
    "df_size = len(basin_inflow)\n",
    "basin_inflow_test = basin_inflow[int(df_size * 0.9):] # last 10%\n",
    "basin_inflow_test.reset_index(drop=True, inplace=True)\n",
    "# obtain absolute maximum inflow value for inverse scaling\n",
    "INFLOW_ABS_MAX = basin_inflow_test['INFLOW'].abs().max()\n",
    "\n",
    "def inverse_scaling(predictions):\n",
    "    next_30_days = predictions[0:30]\n",
    "    \n",
    "    print('Next 30 Days of Inflow\\n')\n",
    "    for idx, target in enumerate(next_30_days):\n",
    "        day = idx + 1\n",
    "        print(\"{}. {:.3f} cfs\".format(day, abs(target[0] * INFLOW_ABS_MAX)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "69b6b751",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\" 60-day windows of basin weather and river data \"\"\"\n",
    "\n",
    "windowTestBatch1 = window.test.take(1)\n",
    "windowTestBatch2 = window.test.take(2)\n",
    "windowTestBatch3 = window.test.take(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "393d47aa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1/1 [==============================] - 0s 302ms/step\n",
      "2/2 [==============================] - 0s 23ms/step\n",
      "3/3 [==============================] - 0s 23ms/step\n"
     ]
    }
   ],
   "source": [
    "\"\"\" Predicted Inflows \"\"\"\n",
    "\n",
    "predictedInflow1 = hypermodel.predict(windowTestBatch1)\n",
    "predictedInflow2 = hypermodel.predict(windowTestBatch2)\n",
    "predictedInflow3 = hypermodel.predict(windowTestBatch3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "84c8e0aa",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Next 30 Days of Inflow\n",
      "\n",
      "1. 6096.670 cfs\n",
      "2. 7168.289 cfs\n",
      "3. 7115.677 cfs\n",
      "4. 6249.618 cfs\n",
      "5. 6506.479 cfs\n",
      "6. 9250.617 cfs\n",
      "7. 8197.594 cfs\n",
      "8. 13122.252 cfs\n",
      "9. 5557.840 cfs\n",
      "10. 8920.194 cfs\n",
      "11. 9449.152 cfs\n",
      "12. 13048.285 cfs\n",
      "13. 19556.594 cfs\n",
      "14. 19501.059 cfs\n",
      "15. 9529.258 cfs\n",
      "16. 9934.721 cfs\n",
      "17. 7019.022 cfs\n",
      "18. 3893.511 cfs\n",
      "19. 12713.444 cfs\n",
      "20. 6187.136 cfs\n",
      "21. 7742.878 cfs\n",
      "22. 10527.569 cfs\n",
      "23. 11802.981 cfs\n",
      "24. 7219.112 cfs\n",
      "25. 7104.785 cfs\n",
      "26. 6460.319 cfs\n",
      "27. 9677.672 cfs\n",
      "28. 7686.907 cfs\n",
      "29. 18663.553 cfs\n",
      "30. 19070.869 cfs\n"
     ]
    }
   ],
   "source": [
    "inverse_scaling(predictedInflow1)"
   ]
  }
 ],
 "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.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}