DynamicFinalModels.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5be20428",
   "metadata": {},
   "source": [
    "# Dynamic Models - Final Model Building and Training"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d4a9f7f",
   "metadata": {},
   "source": [
    "Import Julia machine learning packages, plotting package, and file loader"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6ce70021",
   "metadata": {},
   "outputs": [],
   "source": [
    "using Flux,Statistics,Plots,MLDataUtils,DelimitedFiles\n",
    "using Flux.Data: DataLoader\n",
    "include(\"utilityfunc.jl\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bb554f0",
   "metadata": {},
   "source": [
    "Define the upper and lower bounds on the features and define the scaling and unscaling functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3f3390b0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "inv_oscaler (generic function with 1 method)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Define scaler functions to preprocess data\n",
    "conc_ub=Float32(.6)\n",
    "#set lower and upper bounds on Lpt, Kt, Rs, ti,ci (Lpt is first element, Kt is second element, Rs is the third element, ti is the fourth element, ci is the fifth element\n",
    "i_lb=Float32[5e-7,7e-7,5.0,0.0,0.0] \n",
    "i_ub=Float32[5e-6,5e-6,30,(300/3600),conc_ub]\n",
    "iscaler(x)=(x-i_lb)./(i_ub-i_lb) #Min-Max Normalization on Features\n",
    "inv_iscaler(x)=x.*(i_ub.-i_lb).+i_lb\n",
    "oscaler(x)= x ./ conc_ub #Min-Max Normalization on Target\n",
    "inv_oscaler(x)=x.*conc_ub"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "dc038a7a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "iscalerbatch (generic function with 1 method)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function iscalerbatch(x) #Min-Max Normalization on Matrices\n",
    "    nfeat,sl = size(x)\n",
    "    out=Array{Float32,2}(undef,nfeat,sl)\n",
    "    for i=1:sl\n",
    "        out[:,i]=iscaler(x[:,i])\n",
    "    end\n",
    "    return out\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "62018489",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "inv_iscalerbatch (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function inv_iscalerbatch(x)\n",
    "    nfeat,sl = size(x)\n",
    "    out=Array{Float32,2}(undef,nfeat,sl)\n",
    "    for i=1:sl\n",
    "        out[:,i]=inv_iscaler(x[:,i])\n",
    "    end\n",
    "    return out\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ae89d5e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_samp=Int64(1e4); # number of samples\n",
    "nbatch = n_samp; # number of batches\n",
    "nfeat=5; #5 features\n",
    "ntarg=1; #1 target values\n",
    "sl=20; #sequence length of 20"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cefd60b",
   "metadata": {},
   "source": [
    "Load in the data and scale the features and output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e3d07b6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "d = readdlm(\"MLFinalData_v3.csv\",',',Float32) #entire dataset\n",
    "X= iscalerbatch(d[1:5,:]); #Divide into input and scale\n",
    "Y= oscaler.(d[6,:]); #Divide into outputput and scale"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91aa1512",
   "metadata": {},
   "source": [
    "Split the dataset into training, validation, and testsets using an 80/10/10 split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0b116486",
   "metadata": {},
   "outputs": [],
   "source": [
    "(x_train, y_train), (x_val, y_val),(x_test, y_test) = splitobs((X, Y), at = (0.8,.1) );"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b64c9be5",
   "metadata": {},
   "outputs": [],
   "source": [
    "YVAL = y_val[1:end]\n",
    "YTEST = y_test[1:end];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "f26501ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = DataLoader((x_train,y_train)); #Package data in Flux.jl's dataloader"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecd0e61d",
   "metadata": {},
   "source": [
    "Set up learning rate modes and number of epochs. Here we use a exponentially decaying learning rate for the MLP model and a constant learning rate for the recurrent models. We train each model for a 100 epochs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c94382a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "decay(epoch)=1e-3*exp(-.023*epoch)\n",
    "lr=5e-5;\n",
    "n_epochs = 100;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d6565d6",
   "metadata": {},
   "source": [
    "# MLP Code"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "423c7bf8",
   "metadata": {},
   "source": [
    "Set up MLP model - Here we use 1 hidden layer with 8 nodes, the ReLU activation function, and the MSE loss function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "7d9d47cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_hidden_mlp=8;\n",
    "activation_func=Flux.relu;\n",
    "m_MLP = Chain(Dense(nfeat,n_hidden_mlp),Dense(n_hidden_mlp,n_hidden_mlp,activation_func),Dense(n_hidden_mlp,ntarg));\n",
    "function loss_MLP(x, y)\n",
    "     Flux.mse(m_MLP(x), y)\n",
    "end\n",
    "ps_MLP = params(m_MLP);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2f55bb4",
   "metadata": {},
   "source": [
    "Train the model using the ADAM Optimizer for a 100 epochs with early stopping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca333375",
   "metadata": {},
   "outputs": [],
   "source": [
    "e_log_m = zeros(n_epochs);\n",
    "for i=1:n_epochs\n",
    "    Flux.train!(loss_MLP, ps_MLP, data, ADAM(decay(i)))\n",
    "    e_log_m[i]=cust_mse(x_val,YVAL,m_MLP)\n",
    "    if (i%1 == 0 || i ==1)\n",
    "        @show(i,e_log_m[i])\n",
    "    end\n",
    "    if (e_log_m[i]<1e-7) #early stopping criteria of 1e-7\n",
    "        break\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03de02c9",
   "metadata": {},
   "source": [
    "Evaluate performance metrics on the test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1165b1b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "@show cust_mse(x_test,YTEST,m_MLP)\n",
    "@show mean_error_func(x_test,YTEST,m_MLP)\n",
    "@show RMSE(x_test,YTEST,m_MLP)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b9e0417",
   "metadata": {},
   "source": [
    "# RNN Code"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc5531f5",
   "metadata": {},
   "source": [
    "Set up the vanilla RNN model - Here we use 1 hidden layer with 32 nodes and the MSE loss function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "37189689",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_hidden_rnn=32;\n",
    "#activation_func=Flux.tanh\n",
    "m_RNN = Chain(Dense(nfeat,n_hidden_rnn),RNN(n_hidden_rnn,n_hidden_rnn),Dense(n_hidden_rnn,ntarg))\n",
    "function loss_RNN(x, y)\n",
    "    Flux.reset!(m_RNN)\n",
    "     Flux.mse(m_RNN(x), y)\n",
    "end\n",
    "ps_RNN = params(m_RNN);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "505c6c4d",
   "metadata": {},
   "source": [
    "Train the model using the ADAM Optimizer for a 100 epochs with early stopping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b84fd29",
   "metadata": {},
   "outputs": [],
   "source": [
    "e_log_r = zeros(n_epochs)\n",
    "for i=1:n_epochs\n",
    "    Flux.train!(loss_RNN, ps_RNN, data, ADAM(5e-5))\n",
    "    if (i%1 == 0 || i ==1)\n",
    "        @show(i,cust_mse(x_val,YVAL,m_RNN))\n",
    "    end\n",
    "    if (e_log_r[i]<5e-7)\n",
    "        break\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eafb5146",
   "metadata": {},
   "source": [
    "Evaluate performance metrics on the test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f01b831",
   "metadata": {},
   "outputs": [],
   "source": [
    "@show cust_mse(x_test,YTEST,m_RNN)\n",
    "@show mean_error_func(x_test,YTEST,m_RNN)\n",
    "@show RMSE(x_test,YTEST,m_RNN)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53572559",
   "metadata": {},
   "source": [
    "# LSTM Code"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "89c8a20b",
   "metadata": {},
   "source": [
    "Set up the LSTM model - Here we use 1 hidden layer with 8 nodes and the MSE loss function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "5f4c0f5a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "loss_LSTM (generic function with 1 method)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_hidden_lstm=8;\n",
    "m_LSTM = Chain(Dense(nfeat,n_hidden_lstm),LSTM(n_hidden_lstm,n_hidden_lstm),Dense(n_hidden_lstm,ntarg))\n",
    "function loss_LSTM(x, y)\n",
    "    Flux.reset!(m_LSTM)\n",
    "     Flux.mse(m_LSTM(x), y)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "9b01c4d8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dr (generic function with 1 method)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dr(e)=5e-4*exp(-.016*e)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06050c8a",
   "metadata": {},
   "source": [
    "Train the model using the ADAM Optimizer for a 100 epochs with early stopping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7588c5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ps_LSTM = params(m_LSTM)\n",
    "e_log_l = zeros(n_epochs)\n",
    "for i=1:n_epochs\n",
    "    Flux.train!(loss_LSTM, ps_LSTM, data, ADAM(dr(i)))\n",
    "    e_log_l[i]=cust_mse(x_val,YVAL,m_LSTM)\n",
    "    if (i%1 == 0 || i ==1)\n",
    "        @show(i,e_log_l[i])\n",
    "    end\n",
    "    if (e_log_l[i]<1e-7)\n",
    "        break\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "2fd4cb01",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cust_mse(x_test, YTEST, m_LSTM) = 2.212111f-7\n",
      "mean_error_func(x_test, YTEST, m_LSTM) = 10.310354f0\n",
      "RMSE(x_test, YTEST, m_LSTM) = 0.00047033082f0\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.00047033082f0"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@show cust_mse(x_test,YTEST,m_LSTM)\n",
    "@show mean_error_func(x_test,YTEST,m_LSTM)\n",
    "@show RMSE(x_test,YTEST,m_LSTM)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73e33dde",
   "metadata": {},
   "source": [
    "Export the models (Uncomment lines to save models)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e2352aee",
   "metadata": {},
   "outputs": [],
   "source": [
    "using BSON: @save\n",
    "#@save \"m_MLP.bson\" m_MLP\n",
    "#@save \"m_RNN_f.bson\" m_RNN\n",
    "#@save \"m_LSTM.bson\" m_LSTM"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.6.1",
   "language": "julia",
   "name": "julia-1.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}