model3_NRM-G.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "GeneralizedModel___test.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "machine_shape": "hm"
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ROFaYRSBFRhs"
      },
      "source": [
        "This file is now set for all 62 upstream gauges, end-to-end and available to run directly.\n",
        "\n",
        "I run a test with several epochs, and the results are in the folder as well."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "L2eEDiiEsQ3T"
      },
      "source": [
        "The first several lines make the codes connected to my personal Google's drive.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rt3oqBGQsFHx"
      },
      "source": [
        "from google.colab import drive\n",
        "drive.mount('/content/drive')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UI7fJfFEsJKq"
      },
      "source": [
        "%cd '/content/drive/My Drive/Benchmark'"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "HpfzP3SUseuS",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        },
        "outputId": "3da15645-283b-449a-fd71-5d2caed4977d"
      },
      "source": [
        "# This example is developed using tensorflow 1.15.\n",
        "%tensorflow_version 1.x"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "TensorFlow 1.x selected.\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "rEDiMkviso6B"
      },
      "source": [
        "import numpy as np\n",
        "import pandas as pd\n",
        "    \n",
        "def sequence_data_genearter():\n",
        "\t\"\"\"\n",
        "  station_list = [519,521,522,523,525,526,527,532,534,535,536,537,538,539,541,542,\n",
        "\t\t\t\t\t543,546,549,551,552,553,554,556,557,562,564,565,568,569,572,574,\n",
        "\t\t\t\t\t590,595,599,601,604,605,607,616,617,621,624,626,631,640,641,642,\n",
        "\t\t\t\t\t644,648,653,655,656,659,661,662,663,665,668,669,670,671,540,544,\n",
        "\t\t\t\t\t563,566,570,571,573,575,579,585,586,587,588,589,591,592,593,594,\n",
        "\t\t\t\t\t596,597,598,600,602,603,606,608,609,610,611,612,1688,613,614,615,\n",
        "\t\t\t\t\t618,619,620,622,625,627,628,629,630,632,634,635,636,637,638,639,\n",
        "\t\t\t\t\t643,645,646,649,654,657,658,660,664,666,667,673,1609] # 125 staion IDs in IFIS\n",
        "  \"\"\"          \n",
        "  # IN THIS TEST, my model can work on the following 62 watersheds.\n",
        "\tstation_list = [553,542,522,536,569,543,538,568,574,535,671,539,546,617,534,653,554,552,551,616,527,607,525,565,644,661,562,656,665,642,641,526,556,599,659,670,621,668,648,662,601,624,541,532,557,669,523,537,521,590,549,572,640,631,519,595,605,655,663,564,604,626]\n",
        "\n",
        "\n",
        "\t# convert series to supervised learning sequence data\n",
        "\tdef series_to_supervised(data, n_in=1, n_out=1, dropnan=True):\n",
        "\t\tn_vars = 1 if type(data) is list else data.shape[1]\n",
        "\t\tdf = pd.DataFrame(data)\n",
        "\t\tcols, names = list(), list()\n",
        "\t\t# input sequence (t-n, ... t-1)\n",
        "\t\tfor i in range(n_in, 0, -1):\n",
        "\t\t\tcols.append(df.shift(i))\n",
        "\t\t\tnames += [('var%d(t-%d)' % (j+1, i)) for j in range(n_vars)]\n",
        "\t\t# forecast sequence (t, t+1, ... t+n)\n",
        "\t\tfor i in range(0, n_out):\n",
        "\t\t\tcols.append(df.shift(-i))\n",
        "\t\t\tif i == 0:\n",
        "\t\t\t\tnames += [('var%d(t)' % (j+1)) for j in range(n_vars)]\n",
        "\t\t\telse:\n",
        "\t\t\t\tnames += [('var%d(t+%d)' % (j+1, i)) for j in range(n_vars)]\n",
        "\t\t# put it all together\n",
        "\t\tagg = pd.concat(cols, axis=1)\n",
        "\t\tagg.columns = names\n",
        "\t\t# drop rows with NaN values\n",
        "\t\tif dropnan:\n",
        "\t\t\tagg.dropna(inplace=True)\n",
        "\t\treturn agg\n",
        "\t\n",
        "\tstations_train = pd.DataFrame()\n",
        "\tstations_valid = pd.DataFrame()\n",
        "\tstations_test = pd.DataFrame()\n",
        "\t\n",
        "\tphy = pd.read_csv('./input/Features_Watershed.csv') # slope, travel time, area, and soil type data\n",
        "\tET = pd.read_csv('./input/ET_Iowa.csv') # pretreated min-max scaled, hourly ET data\n",
        "\t\n",
        "\tfor station_id in station_list:\n",
        "\n",
        "\t\twq = pd.read_csv('./input/USGS/'+str(station_id)+'_wq.csv') # watershed discharge Q\n",
        "\t\tpcp = pd.read_csv('./input/PCP/'+str(station_id)+'_pcp_SMA6.csv') # hourly precipitation amount\n",
        "\t\t\n",
        "\t\tfor datatype in ['train','valid','test']:\n",
        "\t\t\tif datatype == 'train':\n",
        "\t\t\t\twq_station = wq[:35064] # the first 4 years as training\n",
        "\t\t\t\tpcp_station = pcp[:35064]\n",
        "\t\t\t\tET_station = ET[:35064]\n",
        "\t\t\telif datatype == 'valid':\n",
        "\t\t\t\twq_station = wq[35064:52608] # the middle 2 years as validation\n",
        "\t\t\t\tpcp_station = pcp[35064:52608]\n",
        "\t\t\t\tET_station = ET[35064:52608]\n",
        "\t\t\telif datatype == 'test':\n",
        "\t\t\t\twq_station = wq[52608:] # the last 1 year as test, the final evaluation\n",
        "\t\t\t\tpcp_station = pcp[52608:]\n",
        "\t\t\t\tET_station = ET[52608:]\n",
        "\t\t\t\n",
        "\t\t\tdataset = pd.concat([wq_station, pcp_station, ET_station], axis=1) # combine the feature of discharge, rainfall, and ET.\n",
        "\t\t\tdataset = dataset.iloc[:,[3,4,1]] # feature shape = [?, 3], '?' is used to represent the instances as in TensorFlow\n",
        "\t\t\tdataset = dataset.fillna(-9999).values # fill NAs using negative values\n",
        "\t\t\tdataset = dataset.astype('float32') # convert to float\n",
        "\t\t\t\n",
        "\t\t\t# reframe the data as a supervised learning task\n",
        "\t\t\treframed = series_to_supervised(dataset, hours_history, hours_forecast) # feature shape = [?, 576], (72+120)*3=576\n",
        "\t\t\treframed_nanlist = reframed[(reframed < -1 ).any(1)].index.tolist() # find rows with any negative value (allowance of 1)\n",
        "\t\t\treframed = reframed.drop(reframed_nanlist) # remove rows with negative value\n",
        "\t\t\t\n",
        "\t\t\t# combining physical data into reframed data\n",
        "\t\t\tphy_station = phy.loc[phy.ifc_id==station_id] # feature shape = [1, 16], the station name and 15 feature values\n",
        "\t\t\tphy_station2 = pd.DataFrame(pd.np.repeat(phy_station.values,reframed.shape[0],axis=0),columns=phy_station.columns) # feature shape = [?, 16]\n",
        "\t\t\treframed = pd.concat([reframed.reset_index(drop=True), phy_station2.reset_index(drop=True)], axis=1, ignore_index=True) # feature shape = [?, 592]\n",
        "\t\t\t\n",
        "\t\t\tif datatype == 'train':\n",
        "\t\t\t\tstations_train = pd.concat([stations_train, reframed], axis=0)\n",
        "\t\t\telif datatype == 'valid':\n",
        "\t\t\t\tstations_valid = pd.concat([stations_valid, reframed], axis=0)\n",
        "\t\t\telif datatype == 'test':\n",
        "\t\t\t\tstations_test = pd.concat([stations_test, reframed], axis=0)\n",
        "\n",
        "\treturn stations_train.values, stations_valid.values, stations_test.values\n",
        "\n",
        "def min_max_normalization(dataset):\n",
        "\n",
        "\t# min-max normalization\n",
        "\t# the max precipitation and discharge should be calculated from the training+validation input.\n",
        "  # Thus, in my generalized model on 62 watereshed. the maximum hourly rainfall in the training+validation period among all 62 watersheds is 899.9.\n",
        "  # And, the maximum hourly discharge is 48775 in the training+validation period among all 62 watersheds.\n",
        "  # If you train a different model, the values here should be different. And the same, the values for the de-normalization in the post-process should be different as well.\n",
        "  \n",
        "\tPCPmax = 899.9092407226562\n",
        "\tPCPmin = 0.0\n",
        "\tQmax = 48775.0\n",
        "\tQmin = 0.0\n",
        "\t\n",
        "\t# normalize the Q and PCP\n",
        "\tPCPlist = list(range(0, dataset.shape[1], 3))[:hours_history+hours_forecast] # from column 0, every 3 features are PCP.\n",
        "\tQlist = list(range(2, dataset.shape[1], 3))[:hours_history+hours_forecast] # from column 2, every 3 features are Q.\n",
        "\tdataset[:, PCPlist] = (dataset[:, PCPlist] - PCPmin) / (PCPmax-PCPmin)\n",
        "\tdataset[:, Qlist] = (dataset[:, Qlist] - Qmin) / (Qmax-Qmin)\n",
        "\t\n",
        "\treturn dataset\n",
        "\n",
        "def X_y_split(train, valid, test):\n",
        "\n",
        "\t# randomize the training data\n",
        "\tnp.random.shuffle(train)\n",
        "\t\n",
        "\t# for the history 72 hours, we know PCP, ET, and Q, which is X1. shape = [?, 72, 3]\n",
        "\t# for the future 120 hours, we know PCP and ET, which is X2. shape = [?, 120, 2]\n",
        "\t# for each watershed, we know physical features, which is X3. shape = [?, 1, 15]\n",
        "\t# for the future 120 hours, we want to get Q, which is y. shape = [?, 120]\n",
        "\t\n",
        "\t# My method is generating two GRU layers, one for history and one for future.\n",
        "\tPCP_history = list(range(0, train.shape[1], 3))[:hours_history]\n",
        "\tPCP_forecast = list(range(0, train.shape[1], 3))[hours_history:hours_history+hours_forecast]\n",
        "\tET_history = list(range(1, train.shape[1], 3))[:hours_history]\n",
        "\tET_forecast = list(range(1, train.shape[1], 3))[hours_history:hours_history+hours_forecast]\n",
        "\tQ_history = list(range(2, train.shape[1], 3))[:hours_history]\n",
        "\tQ_forecast = list(range(2, train.shape[1], 3))[hours_history:hours_history+hours_forecast]\n",
        "\t\n",
        "\tX1_list = PCP_history + ET_history + Q_history\n",
        "\tX1_list.sort()\n",
        "\tX2_list = PCP_forecast + ET_forecast\n",
        "\tX2_list.sort()\n",
        "\n",
        "\t# split the data into X1, X2, X3 and y.\n",
        "\ttrain_X1, train_X2, train_X3, train_y, train_id = train[:,X1_list], train[:,X2_list], train[:,-15:], train[:,Q_forecast], train[:,-16:-15]\n",
        "\tvalid_X1, valid_X2, valid_X3, valid_y, valid_id = valid[:,X1_list], valid[:,X2_list], valid[:,-15:], valid[:,Q_forecast], valid[:,-16:-15]\n",
        "\ttest_X1, test_X2, test_X3, test_y, test_id = test[:,X1_list], test[:,X2_list], test[:,-15:], test[:,Q_forecast], test[:,-16:-15]\n",
        "\n",
        "\t# reshape X1 and X2 into 3D [samples, timesteps, features]\n",
        "\ttrain_X1 = train_X1.reshape(train_X1.shape[0], hours_history, 3)\n",
        "\tvalid_X1 = valid_X1.reshape(valid_X1.shape[0], hours_history, 3)\n",
        "\ttest_X1 = test_X1.reshape(test_X1.shape[0], hours_history, 3)\n",
        "\ttrain_X2 = train_X2.reshape(train_X2.shape[0], hours_forecast, 2)\n",
        "\tvalid_X2 = valid_X2.reshape(valid_X2.shape[0], hours_forecast, 2)\n",
        "\ttest_X2 = test_X2.reshape(test_X2.shape[0], hours_forecast, 2)\n",
        "\t\n",
        "\t# expand X3 to 3D, shape from [?, 15] to [?, 1, 15].\n",
        "\ttrain_X3, valid_X3, test_X3 = np.expand_dims(train_X3, axis=1), np.expand_dims(valid_X3, axis=1), np.expand_dims(test_X3, axis=1)\n",
        "\treturn [train_X1, train_X2, train_X3, train_y, train_id, valid_X1, valid_X2, valid_X3, valid_y, valid_id, test_X1, test_X2, test_X3, test_y, test_id]\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ep2mLBGOtbBT",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 35
        },
        "outputId": "58d62595-5fda-453d-84ac-9d0c8582b46a"
      },
      "source": [
        "# parameters\n",
        "hours_forecast = 120 # 120 hours forecast is the current goal.\n",
        "hours_history = 72 # 72 hours of history is enough for predicting the future.\n",
        "\n",
        "train, valid, test = sequence_data_genearter()\n",
        "\n",
        "train = min_max_normalization(train)\n",
        "valid = min_max_normalization(valid)\n",
        "test = min_max_normalization(test)\n",
        "\n",
        "train_X1, train_X2, train_X3, train_y, train_id, valid_X1, valid_X2, valid_X3, valid_y, valid_id, test_X1, test_X2, test_X3, test_y, test_id = X_y_split(train, valid, test)\n",
        "del train, valid, test\n"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "/usr/local/lib/python3.6/dist-packages/ipykernel_launcher.py:80: FutureWarning: The pandas.np module is deprecated and will be removed from pandas in a future version. Import numpy directly instead\n"
          ],
          "name": "stderr"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8cWsIgxvUjAY"
      },
      "source": [
        "Design of Model"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kQVssbP9tiVJ"
      },
      "source": [
        "import pandas as pd\n",
        "from tensorflow import keras\n",
        "from tensorflow.keras.models import Model\n",
        "from tensorflow.keras.layers import Input, Dense, Dropout, CuDNNGRU, Flatten, TimeDistributed, Lambda, concatenate\n",
        "import tensorflow.keras.backend as K\n",
        "from tensorflow.keras.callbacks import ReduceLROnPlateau, EarlyStopping, ModelCheckpoint\n",
        "\n",
        "def NRM_generalized_basic():\n",
        "\n",
        "\t# design network\n",
        "\tdim_dense=[128, 64, 64, 32, 32]\n",
        "\tdrop=0.2\n",
        "\n",
        "\tphy_input = Input(shape=(train_X3.shape[1],train_X3.shape[2]))\n",
        "\tphy_input_zeros = Lambda(lambda x: x * 0)(phy_input)\n",
        "\t\n",
        "\tencoder_phy = Lambda(lambda x: K.repeat_elements(x, hours_history, axis=1))(phy_input_zeros) # simple encoder without physical data\n",
        "\tdecoder_phy = Lambda(lambda x: K.repeat_elements(x, hours_forecast, axis=1))(phy_input) # physical data is necessary in decoder phase\n",
        "\n",
        "\tencoder_input = Input(shape=(train_X1.shape[1],train_X1.shape[2]))\n",
        "\tencoder_input_phy = concatenate([encoder_input,encoder_phy],axis=-1)\n",
        "\t\n",
        "\tencoder_rnn1 = CuDNNGRU(32, return_state=True, return_sequences=True)\n",
        "\tencoder_output1, encoder_hc1 = encoder_rnn1(encoder_input_phy)\n",
        "\tencoder_rnn2 = CuDNNGRU(32, return_state=True, return_sequences=True)\n",
        "\tencoder_output2, encoder_hc2 = encoder_rnn2(encoder_output1)\n",
        "\tencoder_rnn3 = CuDNNGRU(32, return_state=True, return_sequences=True)\n",
        "\tencoder_output3, encoder_hc3 = encoder_rnn3(encoder_output2)\n",
        "\tencoder_rnn4 = CuDNNGRU(32, return_state=True, return_sequences=True)\n",
        "\tencoder_output4, encoder_hc4 = encoder_rnn4(encoder_output3)\n",
        "\tencoder_rnn5 = CuDNNGRU(32, return_state=True)\n",
        "\tencoder_output5, encoder_hc5 = encoder_rnn5(encoder_output4)\n",
        "\t\n",
        "\tdecoder_input = Input(shape=(train_X2.shape[1],train_X2.shape[2]))\n",
        "\tdecoder_input_phy = concatenate([decoder_input,decoder_phy],axis=-1)\n",
        "\tdecoder_rnn1 = CuDNNGRU(32, return_sequences=True)\n",
        "\tdecoder_rnn2 = CuDNNGRU(32, return_sequences=True)\n",
        "\tdecoder_rnn3 = CuDNNGRU(32, return_sequences=True)\n",
        "\tdecoder_rnn4 = CuDNNGRU(32, return_sequences=True)\n",
        "\tdecoder_rnn5 = CuDNNGRU(32, return_sequences=True)\n",
        "\tx = decoder_rnn1(decoder_input_phy, initial_state=encoder_hc1)\n",
        "\tx = decoder_rnn2(x, initial_state=encoder_hc2)\n",
        "\tx = decoder_rnn3(x, initial_state=encoder_hc3)\n",
        "\tx = decoder_rnn4(x, initial_state=encoder_hc4)\n",
        "\tx = decoder_rnn5(x, initial_state=encoder_hc5)\n",
        "\n",
        "\tfor dim in dim_dense:\n",
        "\t\tx = TimeDistributed(Dense(dim, activation='relu'))(x)\n",
        "\t\tx = TimeDistributed(Dropout(drop))(x)\n",
        "\tmain_out = TimeDistributed(Dense(1, activation='relu'))(x)\n",
        "\tmain_out = Flatten()(main_out)\n",
        "\t\n",
        "\tmodel = Model(inputs=[encoder_input, decoder_input, phy_input], outputs=main_out)\n",
        "\tmodel.summary()\n",
        "\treturn model"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YYTZGoWcv4KB"
      },
      "source": [
        "Other settings of the model.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "5INnZHszuesQ",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "outputId": "20d5972c-e95c-4f31-d159-e4e8f31b3017"
      },
      "source": [
        "model = NRM_generalized_basic()\n",
        "\n",
        "testname = './NRM_generalized_basic'\n",
        "\n",
        "# some technologies used in training\n",
        "reduce_lr = ReduceLROnPlateau(monitor='val_loss', factor=0.1,\n",
        "\t\t\t\t\t\t\t\t  patience=5, cooldown=200, min_lr=1e-8)\n",
        "earlystoping = EarlyStopping(monitor='val_loss', min_delta=0,\n",
        "\t\t\t\t\t\t\t\t\t\t\t patience=10, verbose=1, mode='auto')\n",
        "checkpoint = ModelCheckpoint(testname+'.h5', monitor='val_loss', verbose=1,\n",
        "\t\t\t\t\t\t\t save_best_only=True, save_weights_only=True, mode='min')\n",
        "\n",
        "# setup the loss function and optimizer\n",
        "# optimizer\n",
        "RMSprop=keras.optimizers.RMSprop(lr=0.00003)\n",
        "\n",
        "# loss function\n",
        "# I used this customized loss function. MSE would work as well as a default method.\n",
        "def nseloss(y_true, y_pred):\n",
        "  return K.mean(K.sum((y_pred-y_true)**2,axis=0)/K.sum((y_true-K.mean(y_true))**2,axis=0))\n",
        "\n",
        "model.compile(optimizer=RMSprop, loss=nseloss)\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "WARNING:tensorflow:From /tensorflow-1.15.2/python3.6/tensorflow_core/python/ops/resource_variable_ops.py:1630: calling BaseResourceVariable.__init__ (from tensorflow.python.ops.resource_variable_ops) with constraint is deprecated and will be removed in a future version.\n",
            "Instructions for updating:\n",
            "If using Keras pass *_constraint arguments to layers.\n",
            "Model: \"model\"\n",
            "__________________________________________________________________________________________________\n",
            "Layer (type)                    Output Shape         Param #     Connected to                     \n",
            "==================================================================================================\n",
            "input_1 (InputLayer)            [(None, 1, 15)]      0                                            \n",
            "__________________________________________________________________________________________________\n",
            "lambda (Lambda)                 (None, 1, 15)        0           input_1[0][0]                    \n",
            "__________________________________________________________________________________________________\n",
            "input_2 (InputLayer)            [(None, 72, 3)]      0                                            \n",
            "__________________________________________________________________________________________________\n",
            "lambda_1 (Lambda)               (None, 72, 15)       0           lambda[0][0]                     \n",
            "__________________________________________________________________________________________________\n",
            "input_3 (InputLayer)            [(None, 120, 2)]     0                                            \n",
            "__________________________________________________________________________________________________\n",
            "lambda_2 (Lambda)               (None, 120, 15)      0           input_1[0][0]                    \n",
            "__________________________________________________________________________________________________\n",
            "concatenate (Concatenate)       (None, 72, 18)       0           input_2[0][0]                    \n",
            "                                                                 lambda_1[0][0]                   \n",
            "__________________________________________________________________________________________________\n",
            "concatenate_1 (Concatenate)     (None, 120, 17)      0           input_3[0][0]                    \n",
            "                                                                 lambda_2[0][0]                   \n",
            "__________________________________________________________________________________________________\n",
            "cu_dnngru (CuDNNGRU)            [(None, 72, 32), (No 4992        concatenate[0][0]                \n",
            "__________________________________________________________________________________________________\n",
            "cu_dnngru_5 (CuDNNGRU)          (None, 120, 32)      4896        concatenate_1[0][0]              \n",
            "                                                                 cu_dnngru[0][1]                  \n",
            "__________________________________________________________________________________________________\n",
            "cu_dnngru_1 (CuDNNGRU)          [(None, 72, 32), (No 6336        cu_dnngru[0][0]                  \n",
            "__________________________________________________________________________________________________\n",
            "cu_dnngru_6 (CuDNNGRU)          (None, 120, 32)      6336        cu_dnngru_5[0][0]                \n",
            "                                                                 cu_dnngru_1[0][1]                \n",
            "__________________________________________________________________________________________________\n",
            "cu_dnngru_2 (CuDNNGRU)          [(None, 72, 32), (No 6336        cu_dnngru_1[0][0]                \n",
            "__________________________________________________________________________________________________\n",
            "cu_dnngru_7 (CuDNNGRU)          (None, 120, 32)      6336        cu_dnngru_6[0][0]                \n",
            "                                                                 cu_dnngru_2[0][1]                \n",
            "__________________________________________________________________________________________________\n",
            "cu_dnngru_3 (CuDNNGRU)          [(None, 72, 32), (No 6336        cu_dnngru_2[0][0]                \n",
            "__________________________________________________________________________________________________\n",
            "cu_dnngru_8 (CuDNNGRU)          (None, 120, 32)      6336        cu_dnngru_7[0][0]                \n",
            "                                                                 cu_dnngru_3[0][1]                \n",
            "__________________________________________________________________________________________________\n",
            "cu_dnngru_4 (CuDNNGRU)          [(None, 32), (None,  6336        cu_dnngru_3[0][0]                \n",
            "__________________________________________________________________________________________________\n",
            "cu_dnngru_9 (CuDNNGRU)          (None, 120, 32)      6336        cu_dnngru_8[0][0]                \n",
            "                                                                 cu_dnngru_4[0][1]                \n",
            "__________________________________________________________________________________________________\n",
            "time_distributed (TimeDistribut (None, 120, 128)     4224        cu_dnngru_9[0][0]                \n",
            "__________________________________________________________________________________________________\n",
            "time_distributed_1 (TimeDistrib (None, 120, 128)     0           time_distributed[0][0]           \n",
            "__________________________________________________________________________________________________\n",
            "time_distributed_2 (TimeDistrib (None, 120, 64)      8256        time_distributed_1[0][0]         \n",
            "__________________________________________________________________________________________________\n",
            "time_distributed_3 (TimeDistrib (None, 120, 64)      0           time_distributed_2[0][0]         \n",
            "__________________________________________________________________________________________________\n",
            "time_distributed_4 (TimeDistrib (None, 120, 64)      4160        time_distributed_3[0][0]         \n",
            "__________________________________________________________________________________________________\n",
            "time_distributed_5 (TimeDistrib (None, 120, 64)      0           time_distributed_4[0][0]         \n",
            "__________________________________________________________________________________________________\n",
            "time_distributed_6 (TimeDistrib (None, 120, 32)      2080        time_distributed_5[0][0]         \n",
            "__________________________________________________________________________________________________\n",
            "time_distributed_7 (TimeDistrib (None, 120, 32)      0           time_distributed_6[0][0]         \n",
            "__________________________________________________________________________________________________\n",
            "time_distributed_8 (TimeDistrib (None, 120, 32)      1056        time_distributed_7[0][0]         \n",
            "__________________________________________________________________________________________________\n",
            "time_distributed_9 (TimeDistrib (None, 120, 32)      0           time_distributed_8[0][0]         \n",
            "__________________________________________________________________________________________________\n",
            "time_distributed_10 (TimeDistri (None, 120, 1)       33          time_distributed_9[0][0]         \n",
            "__________________________________________________________________________________________________\n",
            "flatten (Flatten)               (None, 120)          0           time_distributed_10[0][0]        \n",
            "==================================================================================================\n",
            "Total params: 80,385\n",
            "Trainable params: 80,385\n",
            "Non-trainable params: 0\n",
            "__________________________________________________________________________________________________\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7iZ0wbeHIEuG",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 1000
        },
        "outputId": "80a62f63-a016-4296-e902-3ee1f5443d4d"
      },
      "source": [
        "# train the model with large batci size for 5 epochs to get a good initialization weights with dense layers freezed.\n",
        "history = model.fit([train_X1, train_X2, train_X3], train_y, epochs=100, batch_size=64,\n",
        "\t\t\t\t\tvalidation_data=([valid_X1, valid_X2, valid_X3], valid_y), callbacks=[reduce_lr, earlystoping, checkpoint], verbose=1)\n",
        "\n",
        "# save training loss into local file\n",
        "loss_train = history.history['loss']\n",
        "loss_valid = history.history['val_loss']\n",
        "loss_train = pd.DataFrame({'TrainLoss':loss_train})\n",
        "loss_valid = pd.DataFrame({'TestLoss':loss_valid})\n",
        "loss_epoches = pd.concat([loss_train, loss_valid], axis=1)\n",
        "loss_name = testname + '-loss.csv'\n",
        "loss_epoches.to_csv(loss_name, index = True)\n"
      ],
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "text": [
            "WARNING:tensorflow:From /tensorflow-1.15.2/python3.6/tensorflow_core/python/ops/math_grad.py:1375: where (from tensorflow.python.ops.array_ops) is deprecated and will be removed in a future version.\n",
            "Instructions for updating:\n",
            "Use tf.where in 2.0, which has the same broadcast rule as np.where\n",
            "Train on 1375835 samples, validate on 766268 samples\n",
            "Epoch 1/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.6070\n",
            "Epoch 00001: val_loss improved from inf to 0.36188, saving model to ./NRM_generalized_basic.h5\n",
            "1375835/1375835 [==============================] - 1112s 808us/sample - loss: 0.6070 - val_loss: 0.3619\n",
            "Epoch 2/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.3444\n",
            "Epoch 00002: val_loss improved from 0.36188 to 0.28124, saving model to ./NRM_generalized_basic.h5\n",
            "1375835/1375835 [==============================] - 1081s 786us/sample - loss: 0.3444 - val_loss: 0.2812\n",
            "Epoch 3/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.3086\n",
            "Epoch 00003: val_loss improved from 0.28124 to 0.26835, saving model to ./NRM_generalized_basic.h5\n",
            "1375835/1375835 [==============================] - 1071s 778us/sample - loss: 0.3086 - val_loss: 0.2684\n",
            "Epoch 4/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.2822\n",
            "Epoch 00004: val_loss did not improve from 0.26835\n",
            "1375835/1375835 [==============================] - 1068s 776us/sample - loss: 0.2822 - val_loss: 0.4239\n",
            "Epoch 5/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.2630\n",
            "Epoch 00005: val_loss improved from 0.26835 to 0.26587, saving model to ./NRM_generalized_basic.h5\n",
            "1375835/1375835 [==============================] - 1076s 782us/sample - loss: 0.2630 - val_loss: 0.2659\n",
            "Epoch 6/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.2502\n",
            "Epoch 00006: val_loss did not improve from 0.26587\n",
            "1375835/1375835 [==============================] - 1077s 783us/sample - loss: 0.2502 - val_loss: 0.2758\n",
            "Epoch 7/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.2421\n",
            "Epoch 00007: val_loss improved from 0.26587 to 0.23244, saving model to ./NRM_generalized_basic.h5\n",
            "1375835/1375835 [==============================] - 1064s 773us/sample - loss: 0.2421 - val_loss: 0.2324\n",
            "Epoch 8/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.2348\n",
            "Epoch 00008: val_loss did not improve from 0.23244\n",
            "1375835/1375835 [==============================] - 1063s 773us/sample - loss: 0.2348 - val_loss: 0.2603\n",
            "Epoch 9/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.2283\n",
            "Epoch 00009: val_loss did not improve from 0.23244\n",
            "1375835/1375835 [==============================] - 1063s 773us/sample - loss: 0.2283 - val_loss: 0.2417\n",
            "Epoch 10/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.2239\n",
            "Epoch 00010: val_loss improved from 0.23244 to 0.23171, saving model to ./NRM_generalized_basic.h5\n",
            "1375835/1375835 [==============================] - 1064s 773us/sample - loss: 0.2239 - val_loss: 0.2317\n",
            "Epoch 11/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.2189\n",
            "Epoch 00011: val_loss improved from 0.23171 to 0.23107, saving model to ./NRM_generalized_basic.h5\n",
            "1375835/1375835 [==============================] - 1058s 769us/sample - loss: 0.2189 - val_loss: 0.2311\n",
            "Epoch 12/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.2147\n",
            "Epoch 00012: val_loss did not improve from 0.23107\n",
            "1375835/1375835 [==============================] - 1061s 771us/sample - loss: 0.2147 - val_loss: 0.2385\n",
            "Epoch 13/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.2122\n",
            "Epoch 00013: val_loss did not improve from 0.23107\n",
            "1375835/1375835 [==============================] - 1053s 765us/sample - loss: 0.2122 - val_loss: 0.2324\n",
            "Epoch 14/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.2081\n",
            "Epoch 00014: val_loss did not improve from 0.23107\n",
            "1375835/1375835 [==============================] - 1050s 763us/sample - loss: 0.2081 - val_loss: 0.2339\n",
            "Epoch 15/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.2062\n",
            "Epoch 00015: val_loss improved from 0.23107 to 0.22793, saving model to ./NRM_generalized_basic.h5\n",
            "1375835/1375835 [==============================] - 1048s 762us/sample - loss: 0.2062 - val_loss: 0.2279\n",
            "Epoch 16/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.2042\n",
            "Epoch 00016: val_loss improved from 0.22793 to 0.21296, saving model to ./NRM_generalized_basic.h5\n",
            "1375835/1375835 [==============================] - 1054s 766us/sample - loss: 0.2042 - val_loss: 0.2130\n",
            "Epoch 17/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.2025\n",
            "Epoch 00017: val_loss did not improve from 0.21296\n",
            "1375835/1375835 [==============================] - 1043s 758us/sample - loss: 0.2025 - val_loss: 0.2551\n",
            "Epoch 18/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.2010\n",
            "Epoch 00018: val_loss improved from 0.21296 to 0.20443, saving model to ./NRM_generalized_basic.h5\n",
            "1375835/1375835 [==============================] - 1043s 758us/sample - loss: 0.2010 - val_loss: 0.2044\n",
            "Epoch 19/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.1983\n",
            "Epoch 00019: val_loss did not improve from 0.20443\n",
            "1375835/1375835 [==============================] - 1043s 758us/sample - loss: 0.1983 - val_loss: 0.2116\n",
            "Epoch 20/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.1961\n",
            "Epoch 00020: val_loss did not improve from 0.20443\n",
            "1375835/1375835 [==============================] - 1043s 758us/sample - loss: 0.1961 - val_loss: 0.2098\n",
            "Epoch 21/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1946\n",
            "Epoch 00021: val_loss did not improve from 0.20443\n",
            "1375835/1375835 [==============================] - 1043s 758us/sample - loss: 0.1946 - val_loss: 0.2139\n",
            "Epoch 22/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1920\n",
            "Epoch 00022: val_loss did not improve from 0.20443\n",
            "1375835/1375835 [==============================] - 1041s 757us/sample - loss: 0.1920 - val_loss: 0.2344\n",
            "Epoch 23/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.1896\n",
            "Epoch 00023: val_loss improved from 0.20443 to 0.20323, saving model to ./NRM_generalized_basic.h5\n",
            "1375835/1375835 [==============================] - 1048s 761us/sample - loss: 0.1896 - val_loss: 0.2032\n",
            "Epoch 24/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1886\n",
            "Epoch 00024: val_loss did not improve from 0.20323\n",
            "1375835/1375835 [==============================] - 1041s 757us/sample - loss: 0.1886 - val_loss: 0.2192\n",
            "Epoch 25/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.1867\n",
            "Epoch 00025: val_loss improved from 0.20323 to 0.19779, saving model to ./NRM_generalized_basic.h5\n",
            "1375835/1375835 [==============================] - 1043s 758us/sample - loss: 0.1867 - val_loss: 0.1978\n",
            "Epoch 26/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1848\n",
            "Epoch 00026: val_loss did not improve from 0.19779\n",
            "1375835/1375835 [==============================] - 1051s 764us/sample - loss: 0.1848 - val_loss: 0.2092\n",
            "Epoch 27/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1843\n",
            "Epoch 00027: val_loss did not improve from 0.19779\n",
            "1375835/1375835 [==============================] - 1044s 759us/sample - loss: 0.1843 - val_loss: 0.2128\n",
            "Epoch 28/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.1823\n",
            "Epoch 00028: val_loss did not improve from 0.19779\n",
            "1375835/1375835 [==============================] - 1056s 767us/sample - loss: 0.1823 - val_loss: 0.2044\n",
            "Epoch 29/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1810\n",
            "Epoch 00029: val_loss did not improve from 0.19779\n",
            "1375835/1375835 [==============================] - 1054s 766us/sample - loss: 0.1810 - val_loss: 0.2319\n",
            "Epoch 30/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.1788\n",
            "Epoch 00030: val_loss did not improve from 0.19779\n",
            "1375835/1375835 [==============================] - 1060s 770us/sample - loss: 0.1788 - val_loss: 0.2015\n",
            "Epoch 31/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.1672\n",
            "Epoch 00031: val_loss improved from 0.19779 to 0.19397, saving model to ./NRM_generalized_basic.h5\n",
            "1375835/1375835 [==============================] - 1055s 767us/sample - loss: 0.1672 - val_loss: 0.1940\n",
            "Epoch 32/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1666\n",
            "Epoch 00032: val_loss did not improve from 0.19397\n",
            "1375835/1375835 [==============================] - 1058s 769us/sample - loss: 0.1666 - val_loss: 0.1979\n",
            "Epoch 33/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1660\n",
            "Epoch 00033: val_loss improved from 0.19397 to 0.19200, saving model to ./NRM_generalized_basic.h5\n",
            "1375835/1375835 [==============================] - 1061s 771us/sample - loss: 0.1660 - val_loss: 0.1920\n",
            "Epoch 34/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1658\n",
            "Epoch 00034: val_loss did not improve from 0.19200\n",
            "1375835/1375835 [==============================] - 1057s 768us/sample - loss: 0.1658 - val_loss: 0.1956\n",
            "Epoch 35/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1658\n",
            "Epoch 00035: val_loss did not improve from 0.19200\n",
            "1375835/1375835 [==============================] - 1053s 765us/sample - loss: 0.1658 - val_loss: 0.1978\n",
            "Epoch 36/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1657\n",
            "Epoch 00036: val_loss did not improve from 0.19200\n",
            "1375835/1375835 [==============================] - 1052s 765us/sample - loss: 0.1657 - val_loss: 0.1959\n",
            "Epoch 37/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1653\n",
            "Epoch 00037: val_loss improved from 0.19200 to 0.19151, saving model to ./NRM_generalized_basic.h5\n",
            "1375835/1375835 [==============================] - 1056s 767us/sample - loss: 0.1653 - val_loss: 0.1915\n",
            "Epoch 38/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.1653\n",
            "Epoch 00038: val_loss did not improve from 0.19151\n",
            "1375835/1375835 [==============================] - 1055s 767us/sample - loss: 0.1653 - val_loss: 0.1989\n",
            "Epoch 39/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1646\n",
            "Epoch 00039: val_loss did not improve from 0.19151\n",
            "1375835/1375835 [==============================] - 1054s 766us/sample - loss: 0.1646 - val_loss: 0.1953\n",
            "Epoch 40/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1647\n",
            "Epoch 00040: val_loss did not improve from 0.19151\n",
            "1375835/1375835 [==============================] - 1057s 768us/sample - loss: 0.1647 - val_loss: 0.1936\n",
            "Epoch 41/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.1645\n",
            "Epoch 00041: val_loss improved from 0.19151 to 0.18889, saving model to ./NRM_generalized_basic.h5\n",
            "1375835/1375835 [==============================] - 1052s 765us/sample - loss: 0.1645 - val_loss: 0.1889\n",
            "Epoch 42/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1638\n",
            "Epoch 00042: val_loss did not improve from 0.18889\n",
            "1375835/1375835 [==============================] - 1055s 767us/sample - loss: 0.1638 - val_loss: 0.1980\n",
            "Epoch 43/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1638\n",
            "Epoch 00043: val_loss did not improve from 0.18889\n",
            "1375835/1375835 [==============================] - 1051s 764us/sample - loss: 0.1638 - val_loss: 0.1923\n",
            "Epoch 44/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1640\n",
            "Epoch 00044: val_loss did not improve from 0.18889\n",
            "1375835/1375835 [==============================] - 1052s 765us/sample - loss: 0.1640 - val_loss: 0.1911\n",
            "Epoch 45/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1635\n",
            "Epoch 00045: val_loss improved from 0.18889 to 0.18758, saving model to ./NRM_generalized_basic.h5\n",
            "1375835/1375835 [==============================] - 1049s 762us/sample - loss: 0.1635 - val_loss: 0.1876\n",
            "Epoch 46/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.1634\n",
            "Epoch 00046: val_loss did not improve from 0.18758\n",
            "1375835/1375835 [==============================] - 1055s 767us/sample - loss: 0.1634 - val_loss: 0.1892\n",
            "Epoch 47/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1634\n",
            "Epoch 00047: val_loss did not improve from 0.18758\n",
            "1375835/1375835 [==============================] - 1054s 766us/sample - loss: 0.1634 - val_loss: 0.1899\n",
            "Epoch 48/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1632\n",
            "Epoch 00048: val_loss did not improve from 0.18758\n",
            "1375835/1375835 [==============================] - 1054s 766us/sample - loss: 0.1632 - val_loss: 0.1980\n",
            "Epoch 49/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1628\n",
            "Epoch 00049: val_loss did not improve from 0.18758\n",
            "1375835/1375835 [==============================] - 1057s 768us/sample - loss: 0.1628 - val_loss: 0.2003\n",
            "Epoch 50/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1626\n",
            "Epoch 00050: val_loss did not improve from 0.18758\n",
            "1375835/1375835 [==============================] - 1049s 763us/sample - loss: 0.1626 - val_loss: 0.1915\n",
            "Epoch 51/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.1626\n",
            "Epoch 00051: val_loss did not improve from 0.18758\n",
            "1375835/1375835 [==============================] - 1056s 767us/sample - loss: 0.1626 - val_loss: 0.1906\n",
            "Epoch 52/100\n",
            "1375744/1375835 [============================>.] - ETA: 0s - loss: 0.1628\n",
            "Epoch 00052: val_loss did not improve from 0.18758\n",
            "1375835/1375835 [==============================] - 1049s 763us/sample - loss: 0.1629 - val_loss: 0.1933\n",
            "Epoch 53/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1619\n",
            "Epoch 00053: val_loss did not improve from 0.18758\n",
            "1375835/1375835 [==============================] - 1047s 761us/sample - loss: 0.1619 - val_loss: 0.1941\n",
            "Epoch 54/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1622\n",
            "Epoch 00054: val_loss did not improve from 0.18758\n",
            "1375835/1375835 [==============================] - 1056s 768us/sample - loss: 0.1622 - val_loss: 0.1936\n",
            "Epoch 55/100\n",
            "1375808/1375835 [============================>.] - ETA: 0s - loss: 0.1620\n",
            "Epoch 00055: val_loss did not improve from 0.18758\n",
            "1375835/1375835 [==============================] - 1049s 762us/sample - loss: 0.1620 - val_loss: 0.1946\n",
            "Epoch 00055: early stopping\n"
          ],
          "name": "stdout"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "a4eZJ4YWwkDK"
      },
      "source": [
        "Evaluation\n",
        "\n",
        "I am using four statistics now. NSE, KGE, bias, and r (np.corrcoef).\n",
        "The most popular one is NSE, a traditional one. The second popular metric is KGE, a new metric proposed in 2009.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "4E7iCOEQwiYc"
      },
      "source": [
        "def nse(y_true, y_pred):\n",
        "\treturn 1-np.sum((y_pred-y_true)**2)/np.sum((y_true-np.mean(y_true))**2)\n",
        "\n",
        "def kge(y_true, y_pred):\n",
        "\tkge_r = np.corrcoef(y_true,y_pred)[1][0]\n",
        "\tkge_a = np.std(y_pred)/np.std(y_true)\n",
        "\tkge_b = np.mean(y_pred)/np.mean(y_true)\n",
        "\treturn 1-np.sqrt((kge_r-1)**2+(kge_a-1)**2+(kge_b-1)**2)\n",
        "\n",
        "def bias(y_true, y_pred):\n",
        "\treturn np.sum(y_pred)/np.sum(y_true)-1\n",
        "\n",
        "model.load_weights('./NRM_generalized_basic.h5')\n",
        "\n",
        "Q_predict = model.predict([test_X1, test_X2, test_X3])\n",
        "\n",
        "\n",
        "# post-treatment\n",
        "# I used the min-max scalling in the beginning, so I scale back.\n",
        "# Thus, the output was 0-1 scaled, we need to rescaled back to the real streamflow rates in cfs.\n",
        "# In addition, current unit is cfs, can convert to cms here if needed.\n",
        "\n",
        "\n",
        "# for detail, see function min_max_normalization()\n",
        "Q_predict_cfs = Q_predict*48775.0 \n",
        "test_y_cfs = test_y*48775.0"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JEq1R6rwy2s7"
      },
      "source": [
        "station_list = [553,542,522,536,569,543,538,568,574,535,671,539,546,617,534,653,554,552,551,616,527,607,525,565,644,661,562,656,665,642,641,526,556,599,659,670,621,668,648,662,601,624,541,532,557,669,523,537,521,590,549,572,640,631,519,595,605,655,663,564,604,626]\n",
        "\n",
        "for station_id in station_list:\n",
        "  # locate the index of the station\n",
        "  station_idx = np.argwhere(test_id.flatten() == station_id).flatten()\n",
        "\n",
        "  Q_predict_station = Q_predict_cfs[station_idx]\n",
        "  Q_true_station = test_y_cfs[station_idx]\n",
        "\n",
        "  # Save files for later analysis, optional\n",
        "  # np.savetxt(str(station_id)+'_time_series_true.csv',Q_true_station, delimiter=',')\n",
        "  # np.savetxt(str(station_id)+'_time_series_pred.csv',Q_predict_station, delimiter=',')\n",
        "\n",
        "  NSE_test_eachHour = []\n",
        "  r_test_eachHour = []\n",
        "  bias_test_eachHour = []\n",
        "  KGE_test_eachHour = []\n",
        "  for x in range(hours_forecast):\n",
        "    valuePred_test=Q_predict_station[:,x]\n",
        "    valueTrue_test=Q_true_station[:,x]\n",
        "    NSE_test_eachHour.append(nse(valueTrue_test,valuePred_test))\n",
        "    r_test_eachHour.append(np.corrcoef(valueTrue_test,valuePred_test)[0][1])\n",
        "    bias_test_eachHour.append(bias(valueTrue_test,valuePred_test))\n",
        "    KGE_test_eachHour.append(kge(valueTrue_test,valuePred_test))\n",
        "\n",
        "  NSE_test_eachHour=pd.DataFrame(NSE_test_eachHour)\n",
        "  NSE_test_eachHour.columns = ['NSE']\n",
        "  r_test_eachHour=pd.DataFrame(r_test_eachHour)\n",
        "  r_test_eachHour.columns = ['r']\n",
        "  bias_test_eachHour=pd.DataFrame(bias_test_eachHour)\n",
        "  bias_test_eachHour.columns = ['bias']\n",
        "  KGE_test_eachHour=pd.DataFrame(KGE_test_eachHour)\n",
        "  KGE_test_eachHour.columns = ['KGE']\n",
        "\n",
        "  evaluation_result = pd.concat([NSE_test_eachHour, KGE_test_eachHour, r_test_eachHour, bias_test_eachHour], axis=1)\n",
        "  evaluation_name = testname+'-'+str(station_id)+'-evaluation.csv'\n",
        "  evaluation_result.to_csv(evaluation_name, index = True)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HiiTH1onfuiX"
      },
      "source": [
        "Post evaluation...\n",
        "\n",
        "These are basic evaluation for each watershed and each hour.\n",
        "We can then, calculate the median of these watersheds among 120 prediction hours. (1*120 values, each value is the median of the 62 watersheds)\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hrKdSdgYftSQ"
      },
      "source": [
        "# up\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "from pandas import DataFrame, concat\n",
        "\n",
        "eval_station_list = station_list\n",
        "\n",
        "NSE_median = []\n",
        "KGE_median = []\n",
        "r_median = []\n",
        "bias_median = []\n",
        "for i in range(120):\t\n",
        "    NSEs = []\n",
        "    KGEs = []\n",
        "    rs = []\n",
        "    biases = []\n",
        "    for station_id in eval_station_list:\n",
        "        a = pd.read_csv('./NRM_generalized_basic-'+str(station_id)+'-evaluation.csv')\n",
        "        result_nse = a['NSE'][i]\n",
        "        result_kge = a['KGE'][i]\n",
        "        result_r = a['r'][i]\n",
        "        result_bias = a['bias'][i]\n",
        "        NSEs.append(result_nse)\n",
        "        KGEs.append(result_kge)\n",
        "        rs.append(result_r)\n",
        "        biases.append(result_bias)\n",
        "    NSE_median.append(np.median(NSEs))\n",
        "    KGE_median.append(np.median(KGEs))\n",
        "    r_median.append(np.median(rs))\n",
        "    bias_median.append(np.median(biases))\n",
        "print(NSE_median)"
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}