performance_evaluation.ipynb

{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "c1186017-bf51-4733-905d-ae25c138cb1d",
      "metadata": {
        "id": "c1186017-bf51-4733-905d-ae25c138cb1d"
      },
      "outputs": [],
      "source": [
        "from models import models\n",
        "from utils import *\n",
        "from sklearn.metrics import mean_absolute_error\n",
        "from sklearn.metrics import mean_squared_error\n",
        "from sklearn.metrics import r2_score\n",
        "from models.keras_vision_transformer import transformer_layers\n",
        "from models.keras_vision_transformer import swin_layers\n",
        "from sklearn.multioutput import MultiOutputRegressor\n",
        "from sklearn.svm import SVR\n",
        "from sklearn import linear_model\n",
        "from sklearn.preprocessing import PolynomialFeatures, StandardScaler"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2e2307c6-9877-447b-b6b8-59ad8a8304c1",
      "metadata": {
        "id": "2e2307c6-9877-447b-b6b8-59ad8a8304c1"
      },
      "source": [
        "## DL Approaches"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "11d08c89-45b5-4e90-aa61-dd280d328d88",
      "metadata": {
        "id": "11d08c89-45b5-4e90-aa61-dd280d328d88"
      },
      "outputs": [],
      "source": [
        "\n",
        "MODEL_NAME = \"Attn-Unet-pretrained\"\n",
        "MODEL_NAME = \"SWIN-Unet\"\n",
        "MODEL_NAME = \"CASPIAN\"\n",
        "\n",
        "MODEL_NAME = \"CASPIAN_beta\"\n",
        "MODEL_NAME = \"CASPIAN_beta_v2\"\n",
        "#MODEL_NAME = \"CASPIAN_beta_v3\"\n",
        "#MODEL_NAME = \"CASPIAN_beta_v4\"\n",
        "\n",
        "\n",
        "output_1d = False\n",
        "the_mask = get_the_mask()\n",
        "\n",
        "amae = []\n",
        "armse = []\n",
        "ar2 = []\n",
        "delta_05 = []\n",
        "delta_01 = []\n",
        "acc = []\n",
        "artae = []\n",
        "\n",
        "amae_h = []\n",
        "armse_h = []\n",
        "ar2_h = []\n",
        "delta_05_h = []\n",
        "delta_01_h = []\n",
        "acc_h = []\n",
        "artae_h = []\n",
        "\n",
        "if \"_\" in MODEL_NAME:\n",
        "    splits = 2\n",
        "else:\n",
        "    splits = 4\n",
        "\n",
        "for split in range(1,splits):\n",
        "    #Loading the eval. datasets\n",
        "    ds_eval = {\n",
        "        'test': tf.data.Dataset.load(\"./data/test_ds_split_%d\" % split),\n",
        "        'holdout': tf.data.Dataset.load(\"./data/holdout_dataset\"),\n",
        "    }\n",
        "    if \"Attn\" in MODEL_NAME:\n",
        "        ds_eval['test'] = ds_eval['test'].map(lambda f,x,y, yf: expand_input_channels_full(f,x,y, yf))\n",
        "        ds_eval['holdout'] = ds_eval['holdout'].map(lambda f,x,y, yf: expand_input_channels_full(f,x,y, yf))\n",
        "\n",
        "    ds_eval['test'] = ds_eval['test'].batch(1)\n",
        "    ds_eval['holdout'] = ds_eval['holdout'].batch(1)\n",
        "\n",
        "    model = tf.keras.models.load_model(\"./models/trained_models/\"+MODEL_NAME+\"_split_{}\".format(str(split)),\n",
        "                                       compile=False,\n",
        "                                       custom_objects={\"patch_extract\": transformer_layers.patch_extract,\n",
        "                                                      \"patch_embedding\": transformer_layers.patch_embedding,\n",
        "                                                      \"SwinTransformerBlock\": swin_layers.SwinTransformerBlock,\n",
        "                                                      \"patch_merging\": transformer_layers.patch_merging,\n",
        "                                                      \"patch_expanding\": transformer_layers.patch_expanding})\n",
        "    model.summary()\n",
        "    ypred = []\n",
        "    ypred_h = []\n",
        "    ytest = []\n",
        "    yhold = []\n",
        "\n",
        "    for sample in ds_eval['test'].as_numpy_iterator():\n",
        "        scenario, input_grid, label, label_flat = sample\n",
        "        pred = model.predict(input_grid)[0, :, :, 0]\n",
        "        pred = pred[the_mask]\n",
        "        label = label[0, :, :][the_mask]\n",
        "        ypred.append(pred)\n",
        "        ytest.append(label)\n",
        "\n",
        "    for sample in ds_eval['holdout'].as_numpy_iterator():\n",
        "        scenario, input_grid, label, label_flat = sample\n",
        "        pred = model.predict(input_grid)[0, :, :, 0]\n",
        "        pred = pred[the_mask]\n",
        "        label = label[0, :, :][the_mask]\n",
        "        ypred_h.append(pred)\n",
        "        yhold.append(label)\n",
        "\n",
        "    ypred = np.array(ypred)\n",
        "    ypred_h = np.array(ypred_h)\n",
        "    ytest = np.array(ytest)\n",
        "    yhold = np.array(yhold)\n",
        "    print(ypred.shape, ypred_h.shape)\n",
        "    amae.extend(mean_absolute_error(ytest.T, ypred.T, multioutput='raw_values'))\n",
        "    amae_h.extend(mean_absolute_error(yhold.T, ypred_h.T, multioutput='raw_values'))\n",
        "\n",
        "    armse.extend(mean_squared_error(ytest.T, ypred.T, squared=False, multioutput='raw_values'))\n",
        "    armse_h.extend(mean_squared_error(yhold.T, ypred_h.T, squared=False, multioutput='raw_values'))\n",
        "\n",
        "    artae.extend(artae_score(ytest, ypred))\n",
        "    artae_h.extend(artae_score(yhold, ypred_h))\n",
        "\n",
        "    delta_05.extend(delta_score(ytest, ypred, 0.5))\n",
        "    delta_05_h.extend(delta_score(yhold, ypred_h, 0.5))\n",
        "\n",
        "    delta_01.extend(delta_score(ytest, ypred, 0.1))\n",
        "    delta_01_h.extend(delta_score(yhold, ypred_h, 0.1))\n",
        "\n",
        "    ar2.extend(r2_score(ytest.T, ypred.T, multioutput='raw_values'))\n",
        "    ar2_h.extend(r2_score(yhold.T, ypred_h.T, multioutput='raw_values'))\n",
        "\n",
        "    acc.extend(acc0(ytest, ypred))\n",
        "    acc_h.extend(acc0(yhold, ypred_h))\n",
        "\n",
        "\n",
        "\n",
        "print(\"\\n------- Results on the Test dataset -------------\")\n",
        "print(\"\\nAverage Mean Absolute Erorr: %.5f +- %.5f\" % (np.mean(amae), np.std(amae)))\n",
        "print(\"\\nAverage Root Mean Squared Erorr: %.5f +- %.5f\" % (np.mean(armse), np.std(armse)))\n",
        "print(\"\\nAverage Relative Total Absolute Error: %.5f +- %.5f\" % (np.mean(artae), np.std(artae)))\n",
        "print(\"\\nAverage delta_05_err: %.5f %% +- %.5f\" % (np.mean(delta_05), np.std(delta_05)))\n",
        "print(\"\\nAverage delta_01_err: %.5f %% +- %.5f\" % (np.mean(delta_01), np.std(delta_01)))\n",
        "print(\"\\nAverage R2 score: %.5f +- %.5f\" % (np.mean(ar2), np.std(ar2)))\n",
        "print(\"\\nAverage Accuracy (of 0s): %.5f %% +- %.5f\" % (np.mean(acc), np.std(acc)))\n",
        "\n",
        "print(\"\\n------- Results on the Holdout dataset -------------\")\n",
        "print(\"\\nAverage Mean Absolute Erorr: %.5f +- %.5f\" % (np.mean(amae_h), np.std(amae_h)))\n",
        "print(\"\\nAverage Root Mean Squared Erorr: %.5f +- %.5f\" % (np.mean(armse_h), np.std(armse_h)))\n",
        "print(\"\\nAverage Relative Total Absolute Error: %.5f +- %.5f\" % (np.mean(artae_h), np.std(artae_h)))\n",
        "print(\"\\nAverage delta_05_err: %.5f %% +- %.5f\" % (np.mean(delta_05_h), np.std(delta_05_h)))\n",
        "print(\"\\nAverage delta_01_err: %.5f %% +- %.5f\" % (np.mean(delta_01_h), np.std(delta_01_h)))\n",
        "print(\"\\nAverage R2 score: %.5f +- %.5f\" % (np.mean(ar2_h), np.std(ar2_h)))\n",
        "print(\"\\nAverage Accuracy (of 0s): %.5f %% +- %.5f\" % (np.mean(acc_h), np.std(acc_h)))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "62b5d46f-75a1-4429-bea9-f596fe36fc9b",
      "metadata": {
        "id": "62b5d46f-75a1-4429-bea9-f596fe36fc9b"
      },
      "source": [
        "## Traditional ML Approaches & Baseline"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "9a20bec4-e91c-4246-9a81-177c00892bad",
      "metadata": {
        "id": "9a20bec4-e91c-4246-9a81-177c00892bad"
      },
      "outputs": [],
      "source": [
        "mode = \"ML\"\n",
        "\n",
        "model_name = \"baseline\"\n",
        "model_name = \"linear\"\n",
        "model_name = \"krig\"\n",
        "model_name = \"svr\"\n",
        "model_name = \"poly_lasso\"\n",
        "\n",
        "amae = []\n",
        "armse = []\n",
        "ar2 = []\n",
        "delta_05 = []\n",
        "delta_01 = []\n",
        "acc = []\n",
        "artae = []\n",
        "\n",
        "amae_h = []\n",
        "armse_h = []\n",
        "ar2_h = []\n",
        "delta_05_h = []\n",
        "delta_01_h = []\n",
        "acc_h = []\n",
        "artae_h = []\n",
        "\n",
        "\n",
        "for split in range(1,4):\n",
        "    if \"basel\" in model_name:\n",
        "        xdata, ydata, xtest, ytest, xhold, yhold = tfDS_to_np (split=split, baseline=True)\n",
        "    else:\n",
        "        xdata, ydata, xtest, ytest, xhold, yhold = tfDS_to_np (split=split)\n",
        "\n",
        "    if model_name == \"krig\":\n",
        "        Z, pcs, means = PCA(ydata)\n",
        "        tt = np.dot(Z, pcs) + means\n",
        "        tt = tt - ydata\n",
        "        print(\"Loss due to PCA: %.15f \" % tt.sum())\n",
        "        model = Kriging(xdata, Z, corr='abs_exp', poly='linear')\n",
        "        ypred = np.dot(model.predict_values(xtest), pcs) + means\n",
        "        ypred_h = np.dot(model.predict_values(xhold), pcs) + means\n",
        "    elif model_name == \"svr\":\n",
        "        print (\"Fitting split %d .....\" % split)\n",
        "        base = SVR(C=5, epsilon=0.05, kernel=\"linear\")\n",
        "        model = MultiOutputRegressor(base, n_jobs=-1).fit(xdata, ydata)\n",
        "        ypred = model.predict(xtest)\n",
        "        ypred_h = model.predict(xhold)\n",
        "    elif model_name == 'linear':\n",
        "        print (\"Fitting split %d .....\" % split)\n",
        "        model = linear_model.LinearRegression().fit(xdata, ydata)\n",
        "        ypred = model.predict(xtest)\n",
        "        ypred_h = model.predict(xhold)\n",
        "    elif model_name == 'poly_lasso':\n",
        "        print (\"Fitting split %d .....\" % split)\n",
        "        poly = PolynomialFeatures(interaction_only=True)\n",
        "        poly_xdata = poly.fit_transform(xdata)\n",
        "        print(\"Shape of the input %s\" % str(poly_xdata.shape))\n",
        "        model = linear_model.Lasso(alpha=0.025, max_iter=2000).fit(poly_xdata, ydata)\n",
        "        poly_xtest = poly.fit_transform(xtest)\n",
        "        poly_xhold = poly.fit_transform(xhold)\n",
        "        ypred = model.predict(poly_xtest)\n",
        "        ypred_h = model.predict(poly_xhold)\n",
        "    else:\n",
        "        print (\"Fitting split %d .....\" % split)\n",
        "        avg_pwl = dataset_average_pwl(ydata)\n",
        "        ypred =  simple_predictor(xtest, avg_pwl)\n",
        "        ypred_h = simple_predictor(xhold, avg_pwl)\n",
        "    ypred [ypred < 0] = 0\n",
        "    ypred_h [ypred_h < 0] = 0\n",
        "\n",
        "    amae.extend(mean_absolute_error(ytest.T, ypred.T, multioutput='raw_values'))\n",
        "    amae_h.extend(mean_absolute_error(yhold.T, ypred_h.T, multioutput='raw_values'))\n",
        "\n",
        "    armse.extend(mean_squared_error(ytest.T, ypred.T, squared=False, multioutput='raw_values'))\n",
        "    armse_h.extend(mean_squared_error(yhold.T, ypred_h.T, squared=False, multioutput='raw_values'))\n",
        "\n",
        "    artae.extend(artae_score(ytest, ypred))\n",
        "    artae_h.extend(artae_score(yhold, ypred_h))\n",
        "\n",
        "    delta_05.extend(delta_score(ytest, ypred, 0.5))\n",
        "    delta_05_h.extend(delta_score(yhold, ypred_h, 0.5))\n",
        "\n",
        "    delta_01.extend(delta_score(ytest, ypred, 0.1))\n",
        "    delta_01_h.extend(delta_score(yhold, ypred_h, 0.1))\n",
        "\n",
        "    ar2.extend(r2_score(ytest.T, ypred.T, multioutput='raw_values'))\n",
        "    ar2_h.extend(r2_score(yhold.T, ypred_h.T, multioutput='raw_values'))\n",
        "\n",
        "    acc.extend(acc0(ytest, ypred))\n",
        "    acc_h.extend(acc0(yhold, ypred_h))\n",
        "\n",
        "\n",
        "\n",
        "print(\"\\n------- Results on the Test dataset -------------\")\n",
        "print(\"\\nAverage Mean Absolute Erorr: %.5f +- %.5f\" % (np.mean(amae), np.std(amae)))\n",
        "print(\"\\nAverage Root Mean Squared Erorr: %.5f +- %.5f\" % (np.mean(armse), np.std(armse)))\n",
        "print(\"\\nAverage Relative Total Absolute Error: %.5f +- %.5f\" % (np.mean(artae), np.std(artae)))\n",
        "print(\"\\nAverage delta_05_err: %.5f %% +- %.5f\" % (np.mean(delta_05), np.std(delta_05)))\n",
        "print(\"\\nAverage delta_01_err: %.5f %% +- %.5f\" % (np.mean(delta_01), np.std(delta_01)))\n",
        "print(\"\\nAverage R2 score: %.5f +- %.5f\" % (np.mean(ar2), np.std(ar2)))\n",
        "print(\"\\nAverage Accuracy (of 0s): %.5f %% +- %.5f\" % (np.mean(acc), np.std(acc)))\n",
        "\n",
        "print(\"\\n------- Results on the Holdout dataset -------------\")\n",
        "print(\"\\nAverage Mean Absolute Erorr: %.5f +- %.5f\" % (np.mean(amae_h), np.std(amae_h)))\n",
        "print(\"\\nAverage Root Mean Squared Erorr: %.5f +- %.5f\" % (np.mean(armse_h), np.std(armse_h)))\n",
        "print(\"\\nAverage Relative Total Absolute Error: %.5f +- %.5f\" % (np.mean(artae_h), np.std(artae_h)))\n",
        "print(\"\\nAverage delta_05_err: %.5f %% +- %.5f\" % (np.mean(delta_05_h), np.std(delta_05_h)))\n",
        "print(\"\\nAverage delta_01_err: %.5f %% +- %.5f\" % (np.mean(delta_01_h), np.std(delta_01_h)))\n",
        "print(\"\\nAverage R2 score: %.5f +- %.5f\" % (np.mean(ar2_h), np.std(ar2_h)))\n",
        "print(\"\\nAverage Accuracy (of 0s): %.5f %% +- %.5f\" % (np.mean(acc_h), np.std(acc_h)))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a29f1f9f-18b6-49c0-a5d3-2f0a2de4d6fc",
      "metadata": {
        "id": "a29f1f9f-18b6-49c0-a5d3-2f0a2de4d6fc"
      },
      "source": [
        "## Visualize Predictions"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "b98f5345-3df4-4d31-ac35-51034583c821",
      "metadata": {
        "id": "b98f5345-3df4-4d31-ac35-51034583c821"
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "from matplotlib.colors import LinearSegmentedColormap, BoundaryNorm, ListedColormap\n",
        "%matplotlib widget\n",
        "\n",
        "colors = [(1, 1, 1),\n",
        "          (0, 1, 0),\n",
        "          (1.0, 0, 0.0)]  # Color for 1 (red)\n",
        "\n",
        "# Create a ListedColormap\n",
        "incmap = ListedColormap(colors)\n",
        "\n",
        "colors = [(0.863, 0.824, 0.80),  # Very light gray for 0\n",
        "          (0.0, 0.0, 1.0),    # Blue for (0, 1]\n",
        "          (0.098, 0.098, 0.439),    # Dark blue for (1, 2]\n",
        "          (0.627, 0.125, 0.941),    # Violet for (2, 3]\n",
        "          (0.294, 0, 0.51),    # Dark violet for (3, 4]\n",
        "          (0.855, 0.439, 0.839)]    # Red for > 4\n",
        "outcmap = ListedColormap(colors)\n",
        "boundaries = [0, 4, 4.5, 5, 5.5, 7.0]\n",
        "outnorm = BoundaryNorm(boundaries, outcmap.N, clip=True)\n",
        "\n",
        "ds_eval = {\n",
        "    'test': tf.data.Dataset.load(\"./data/test_ds_split_%d\" % split),\n",
        "    'holdout': tf.data.Dataset.load(\"./data/holdout_dataset\"),\n",
        "}\n",
        "\n",
        "\n",
        "if mode != \"ML\":\n",
        "    scenario = \"00000001111110000.input.npy\"\n",
        "    plt.close()\n",
        "    fig, axs = plt.subplots(1,3, figsize=(12, 6), facecolor='w', edgecolor='k')\n",
        "    fig.tight_layout(pad=2, h_pad=1.5, w_pad=0)\n",
        "    axs = axs.ravel()\n",
        "\n",
        "    ds = ds_eval['holdout'].filter(find_scenario(scenario))\n",
        "    if (len(list(ds.as_numpy_iterator()))) != 0:\n",
        "\n",
        "        for sample in ds.as_numpy_iterator():\n",
        "            scenario, input_grid, label, label_flat = sample\n",
        "            pred = model.predict(input_grid)[0, :, :, 0]\n",
        "            pred[~the_mask] = 0\n",
        "            axs[0].imshow(input_grid[0, :, :], cmap=incmap, origin=\"lower\", interpolation='none', aspect='auto')\n",
        "            axs[0].set_title(scenario[0])\n",
        "            axs[1].imshow(label[0, :, :], cmap=outcmap, norm=outnorm, origin=\"lower\", interpolation='none', aspect='auto')\n",
        "            axs[1].set_title(\"Ground Truth\")\n",
        "            axs[2].imshow(pred, cmap=outcmap, norm=outnorm, origin=\"lower\", interpolation='none', aspect='auto')\n",
        "            axs[2].set_title(\"Prediction\")\n",
        "            print(\"Ground Truth MinMax: %f %f\" % (label.min(), label.max()))\n",
        "            print(\"Prediction MinMax: %f %f \" % (pred.min(), pred.max()))\n",
        "            fig = plt.figure(figsize=(10, 5))\n",
        "            ax = fig.gca()\n",
        "            ax.scatter(label[0, :, :][the_mask], pred[the_mask], s=1)\n",
        "            ax.plot([0, 7], [0,7], color=\"yellow\", alpha=0.5)\n",
        "            ax.set_xlabel(\"Actual Peak Water Levels\")\n",
        "            ax.set_ylabel(\"Predicted Peak Water Levels\")\n",
        "            plt.show()\n",
        "    else:\n",
        "        print(\"Scenario not found\")\n",
        "else:\n",
        "    model.predict(xhold[0].reshape(1, -1))\n",
        "    fig = plt.figure(figsize=(10, 5))\n",
        "    ax = fig.gca()\n",
        "    ax.scatter(yhold[0], ypred_h[0], s=1)\n",
        "    ax.plot([0, 7], [0,7], color=\"yellow\", alpha=0.5)\n",
        "    ax.set_xlabel(\"Actual Peak Water Levels\")\n",
        "    ax.set_ylabel(\"Predicted Peak Water Levels\")\n",
        "    plt.show()"
      ]
    }
  ],
  "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.10.11"
    },
    "colab": {
      "provenance": []
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}