Pretrained.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/jlab-sensing/MFC_Modeling/blob/main/Pretrained.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "#  Pretrained Models\n",
        "###  In order to run the code in this notebook, you must download the `ucscMFCDataset` directory and `stanfordMFCDataset.zip`, which expands into the directory `rocket4`, from [Hugging Face](https://huggingface.co/datasets/adunlop621/Soil_MFC/tree/main), and store them in the same directory as this notebook. You can also find several pretrained models in the at this link, with the naming conventions described in the [README](https://github.com/jlab-sensing/MFC_Modeling#:~:text=Repository%20files%20navigation-,README,-MFC_Modeling)\n"
      ],
      "metadata": {
        "id": "zWfnXLVWrBcK"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "zcF-Dv23N_ON",
        "collapsed": true
      },
      "outputs": [],
      "source": [
        "!pip install --upgrade hepml\n",
        "!pip install arrow\n",
        "!pip install keras_lr_finder\n",
        "# reload modules before executing user code\n",
        "#%load_ext autoreload\n",
        "# reload all modules every time before executing Python code\n",
        "#%autoreload 2\n",
        "# render plots in notebook\n",
        "%matplotlib inline\n",
        "import datetime\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "import glob\n",
        "from sklearn.model_selection import train_test_split\n",
        "import matplotlib.pyplot as plt\n",
        "import seaborn as sns\n",
        "from hepml.core import plot_regression_tree\n",
        "sns.set(color_codes=True)\n",
        "sns.set_palette(sns.color_palette(\"muted\"))\n",
        "from sklearn.model_selection import train_test_split\n",
        "from sklearn.metrics import mean_squared_error\n",
        "from sklearn.ensemble import RandomForestRegressor\n",
        "from sklearn.metrics import mean_absolute_percentage_error as MAPE\n",
        "from sklearn.metrics import mean_absolute_error"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "##  1. Load and Format Data"
      ],
      "metadata": {
        "id": "uUFONhAaeCDl"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "###  1.1 Load and format dataset 1\n",
        "\n",
        "### Remember to download `stanfordMFCDataset.zip`, which expands into the directory `rocket4`, from [Hugging Face](https://huggingface.co/datasets/adunlop621/Soil_MFC/tree/main), and store it in the same directory as this notebook before executing the following code."
      ],
      "metadata": {
        "id": "qtRveSm6rzh3"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "B7PoS8fwOS-Q"
      },
      "outputs": [],
      "source": [
        "#Load teros data\n",
        "import glob\n",
        "teros_files = glob.glob(\"rocket4/TEROSoutput*.csv\")\n",
        "X = pd.DataFrame()\n",
        "for f in teros_files:\n",
        "  try:\n",
        "    csv = pd.read_csv(f, index_col=False).dropna()\n",
        "    X = pd.concat([X, csv])\n",
        "  except:\n",
        "    continue"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "zaaHwFRvXtN4"
      },
      "outputs": [],
      "source": [
        "#Load power data\n",
        "power_files = glob.glob(\"rocket4/soil*.csv\")\n",
        "y = pd.DataFrame()\n",
        "for f in sorted(power_files, key=lambda x: int(x.split('.')[0].split('_')[-1])):\n",
        "#in power_files:\n",
        "  try:\n",
        "    csv = pd.read_csv(f, on_bad_lines='skip', skiprows=10).dropna(how='all')\n",
        "    csv = csv.rename({'Unnamed: 0': 'timestamp'}, axis='columns')\n",
        "    y = pd.concat([y,csv])\n",
        "  except:\n",
        "    continue\n",
        "y[\"timestamp\"] = y[\"timestamp\"].round(decimals = 1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ViS0TjC4n6lk"
      },
      "outputs": [],
      "source": [
        "#Convert current to amps, voltage to volts\n",
        "y[\"I1L [10pA]\"] = np.abs(y[\"I1L [10pA]\"] * 1E-11)\n",
        "y[\"V1 [10nV]\"] = np.abs(y[\"V1 [10nV]\"] * 1E-8)\n",
        "y[\"I1H [nA]\"] = np.abs(y[\"I1H [nA]\"] * 1E-9)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "z0aqozsZpSOB",
        "outputId": "e21bd865-0dc3-4a90-84ab-ddb85cd6befb"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "<ipython-input-5-38f7a5b5fe7c>:14: SettingWithCopyWarning: \n",
            "A value is trying to be set on a copy of a slice from a DataFrame.\n",
            "Try using .loc[row_indexer,col_indexer] = value instead\n",
            "\n",
            "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
            "  df.timestamp = df.timestamp.dt.tz_localize('UTC').dt.tz_convert('US/Pacific')\n"
          ]
        }
      ],
      "source": [
        "#Sort data by timestamp, convert to datetime\n",
        "X = X.sort_values(['timestamp'])\n",
        "y = y.sort_values(['timestamp'])\n",
        "X['timestamp'] = pd.to_datetime(X['timestamp'], unit='s')\n",
        "y['timestamp'] = pd.to_datetime(y['timestamp'], unit='s')\n",
        "\n",
        "#Merge data by timestamp\n",
        "uncut_df = pd.merge_asof(left=X,right=y,direction='nearest',tolerance=pd.Timedelta('1 sec'), on = 'timestamp').dropna(how='all')\n",
        "\n",
        "#Isolate data from cell0\n",
        "df = uncut_df.loc[uncut_df['sensorID'] == 0]\n",
        "\n",
        "#Localize timestamp\n",
        "df.timestamp = df.timestamp.dt.tz_localize('UTC').dt.tz_convert('US/Pacific')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "CBAfOHB61Jwc"
      },
      "outputs": [],
      "source": [
        "#Use only data from before/after deployment date\n",
        "df = df.loc[(df['timestamp'] < '2021-06-04')] #Deployment date\n",
        "\n",
        "#Power drop\n",
        "df = df.drop(df[(df.timestamp >= '2021-05-16 14:46:00') & (df.timestamp <= '2021-05-16 14:46:15')].index)\n",
        "\n",
        "#Drop data outages\n",
        "#df = df.drop(df[(df.timestamp > '2021-11-11') & (df.timestamp < '2021-11-22 01:00:00')].index)\n",
        "#df = df.drop(df[(df.timestamp > '2022-01-27')].index)\n",
        "#df = df.set_index('timestamp')\n",
        "#df = df[:-1]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "p8fFsoLNfrFA"
      },
      "outputs": [],
      "source": [
        "df = df.set_index('timestamp')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "CDOfybleg4cn"
      },
      "outputs": [],
      "source": [
        "#Get time since deployement\n",
        "df['tsd'] = (df.index - df.index[0]).days\n",
        "df['hour'] = (df.index).hour"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "-Pfyv0fM3te1"
      },
      "outputs": [],
      "source": [
        "#Calculate power\n",
        "df[\"power\"] = np.abs(np.multiply(df.iloc[:, 7], df.iloc[:, 8]))\n",
        "#df[\"power\"] = np.abs(np.multiply(df[\"I1L [10pA]\"], df[\"V1 [10nV]\"]))\n",
        "\n",
        "#Convert to microWatts\n",
        "df['power'] = df['power']*1E6"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "jA-WVzVh2-lf"
      },
      "outputs": [],
      "source": [
        "#Convert to microamps, millivolts\n",
        "df[\"I1L [10pA]\"] = np.abs(df[\"I1L [10pA]\"] * 1E6)\n",
        "df[\"V1 [10nV]\"] = np.abs(df[\"V1 [10nV]\"] * 1E3)\n",
        "df[\"I1H [nA]\"] = np.abs(df[\"I1H [nA]\"] * 1E6)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "XxAlm4FXh57m"
      },
      "outputs": [],
      "source": [
        "df = df.reset_index()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "8sHtXZCVrsQJ"
      },
      "outputs": [],
      "source": [
        "#Add power time series\n",
        "df['power - 1h'] = df['power'].shift(1).dropna()\n",
        "df['power - 2h'] = df['power'].shift(2).dropna()\n",
        "df['power - 3h'] = df['power'].shift(3).dropna()\n",
        "#df['power - 2h'] = df['power'].shift(2).dropna()\n",
        "#df['previous_power - 3'] = df['power'].shift(3).dropna()\n",
        "#df['previous_power - 4'] = df['power'].shift(4).dropna()\n",
        "\n",
        "#Add teros time series\n",
        "df['EC - 1h'] = df['EC'].shift(1).dropna()\n",
        "df['EC - 2h'] = df['EC'].shift(2).dropna()\n",
        "df['EC - 3h'] = df['EC'].shift(3).dropna()\n",
        "\n",
        "df['temp - 1h'] = df['temp'].shift(1).dropna()\n",
        "df['temp - 2h'] = df['temp'].shift(2).dropna()\n",
        "df['temp - 3h'] = df['temp'].shift(3).dropna()\n",
        "\n",
        "df['raw_VWC - 1h'] = df['raw_VWC'].shift(1).dropna()\n",
        "df['raw_VWC - 2h'] = df['raw_VWC'].shift(2).dropna()\n",
        "df['raw_VWC - 3h'] = df['raw_VWC'].shift(3).dropna()\n",
        "\n",
        "#Add voltage and current time series\n",
        "df['V1 - 1h'] = df['V1 [10nV]'].shift(1).dropna()\n",
        "df['V1 - 2h'] = df['V1 [10nV]'].shift(2).dropna()\n",
        "df['V1 - 3h'] = df['V1 [10nV]'].shift(3).dropna()\n",
        "\n",
        "df['I1L - 1h'] = df['I1L [10pA]'].shift(1).dropna()\n",
        "df['I1L - 2h'] = df['I1L [10pA]'].shift(2).dropna()\n",
        "df['I1L - 3h'] = df['I1L [10pA]'].shift(3).dropna()\n",
        "\n",
        "df['I1H - 1h'] = df['I1H [nA]'].shift(1).dropna()\n",
        "df['I1H - 2h'] = df['I1H [nA]'].shift(2).dropna()\n",
        "df['I1H - 3h'] = df['I1H [nA]'].shift(3).dropna()\n",
        "df = df.dropna()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "u4SM1_EvGS6y"
      },
      "outputs": [],
      "source": [
        "#df = df.rename(columns={'power': 'power [μW]'})\n",
        "df = df.rename(columns={'I1L [10pA]': 'Current (uA)', 'V1 [10nV]' : 'Voltage (mV)', 'power' : 'Power (uW)'})\n",
        "df = df.set_index('timestamp')"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "###  1.2 Load And Format Dataset 2\n",
        "\n",
        "### Remember to download the `ucscMFCDataset` directory from [Hugging Face](https://huggingface.co/datasets/adunlop621/Soil_MFC/tree/main), and store it in the same directory as this notebook before executing the following code."
      ],
      "metadata": {
        "id": "wS8WW-PZei1g"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "KrhS93hE4RrM",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "7777ce6f-aab1-4fc3-fa5b-821e30057830"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Mounted at /content/drive\n"
          ]
        }
      ],
      "source": [
        "#Load testing data\n",
        "csv = pd.read_csv(\"ucscMFCDataset/small_scaffold_8cm_2.csv\", index_col=False)#.dropna()\n",
        "csv['timestamp'] = csv['timestamp'].str.replace('T', ' ')\n",
        "csv = csv.drop(csv[(csv.timestamp >= '2024-04-19 11:00:00') & (csv.timestamp <= '2024-04-24 08:00:00')].index)\n",
        "#csv = csv.drop(csv[(csv.timestamp >= '2024-04-18 11:00:00') & (csv.timestamp <= '2024-04-25 08:00:00')].index)\n",
        "#csv = csv.drop(csv[(csv.timestamp >= '2024-04-21 20:00:00') & (csv.timestamp <= '2024-04-23 19:00:00')].index)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "FP7Cp6c752TJ"
      },
      "outputs": [],
      "source": [
        "X_test = pd.DataFrame()\n",
        "X_test['timestamp'] = csv['timestamp']\n",
        "X_test['raw_VWC'] = csv['Raw VWC']\n",
        "X_test['temp'] = csv['Temperature (C)']\n",
        "X_test['EC'] = csv['EC (uS/cm)']\n",
        "X_test = X_test.dropna()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "PQWz5Fej8ZWK"
      },
      "outputs": [],
      "source": [
        "y_test = pd.DataFrame()\n",
        "y_test['timestamp'] = csv['timestamp']\n",
        "y_test['Voltage (mV)'] = csv['Voltage (mV)']\n",
        "y_test['Current (uA)'] = csv['Current (uA)']\n",
        "y_test['Power (uW)'] = csv['Power (uW)']\n",
        "y_test = y_test.dropna()\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Eu1jPYNjkX6k"
      },
      "outputs": [],
      "source": [
        "X_test['timestamp'] = pd.to_datetime(X_test['timestamp'])\n",
        "y_test['timestamp'] = pd.to_datetime(y_test['timestamp'])\n",
        "\n",
        "#Merge data by timestamp\n",
        "df_test = pd.merge_asof(left=X_test,right=y_test,direction='nearest',tolerance=pd.Timedelta('3 sec'), on = 'timestamp').dropna()\n",
        "\n",
        "#Localize timestamp\n",
        "df_test.timestamp = df_test.timestamp.dt.tz_convert('UTC').dt.tz_convert('US/Pacific')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ynFP8WB8mEvM"
      },
      "outputs": [],
      "source": [
        "df_test = df_test.set_index('timestamp')\n",
        "#Get time since deployement\n",
        "df_test['tsd'] = (df_test.index - df_test.index[0]).days\n",
        "df_test['hour'] = (df_test.index).hour"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "_adeOWEumdzD"
      },
      "outputs": [],
      "source": [
        "#Add power time series\n",
        "df_test['power - 1h'] = df_test['Power (uW)'].shift(1).dropna()\n",
        "df_test['power - 2h'] = df_test['Power (uW)'].shift(2).dropna()\n",
        "df_test['power - 3h'] = df_test['Power (uW)'].shift(3).dropna()\n",
        "#df['power - 2h'] = df['power'].shift(2).dropna()\n",
        "#df['previous_power - 3'] = df['power'].shift(3).dropna()\n",
        "#df['previous_power - 4'] = df['power'].shift(4).dropna()\n",
        "\n",
        "#Add teros time series\n",
        "df_test['EC - 1h'] = df_test['EC'].shift(1).dropna()\n",
        "df_test['EC - 2h'] = df_test['EC'].shift(2).dropna()\n",
        "df_test['EC - 3h'] = df_test['EC'].shift(3).dropna()\n",
        "\n",
        "df_test['temp - 1h'] = df_test['temp'].shift(1).dropna()\n",
        "df_test['temp - 2h'] = df_test['temp'].shift(2).dropna()\n",
        "df_test['temp - 3h'] = df_test['temp'].shift(3).dropna()\n",
        "\n",
        "df_test['raw_VWC - 1h'] = df_test['raw_VWC'].shift(1).dropna()\n",
        "df_test['raw_VWC - 2h'] = df_test['raw_VWC'].shift(2).dropna()\n",
        "df_test['raw_VWC - 3h'] = df_test['raw_VWC'].shift(3).dropna()\n",
        "\n",
        "#Add voltage and current time series\n",
        "df_test['V1 - 1h'] = df_test['Voltage (mV)'].shift(1).dropna()\n",
        "df_test['V1 - 2h'] = df_test['Voltage (mV)'].shift(2).dropna()\n",
        "df_test['V1 - 3h'] = df_test['Voltage (mV)'].shift(3).dropna()\n",
        "\n",
        "df_test['I1L - 1h'] = df_test['Current (uA)'].shift(1).dropna()\n",
        "df_test['I1L - 2h'] = df_test['Current (uA)'].shift(2).dropna()\n",
        "df_test['I1L - 3h'] = df_test['Current (uA)'].shift(3).dropna()\n",
        "\n",
        "df_test = df_test.dropna()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "GsVfeqKjz58s"
      },
      "outputs": [],
      "source": [
        "df = df.rename(columns={'I1L [10pA]': 'Current (uA)', 'V1 [10nV]' : 'Voltage (mV)', 'power' : 'Power (uW)'})"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "SlMWcsxzDC9H"
      },
      "outputs": [],
      "source": [
        "#All Data (For Type 1 and 2 Models)\n",
        "X = pd.concat([df[\"power - 1h\"], df[\"power - 2h\"], df[\"power - 3h\"], df[\"V1 - 1h\"], df[\"V1 - 2h\"], df[\"V1 - 3h\"], df[\"I1L - 1h\"], df[\"I1L - 2h\"], df[\"I1L - 3h\"],df[\"EC - 1h\"], df[\"EC - 2h\"], df[\"EC - 3h\"], df[\"raw_VWC - 1h\"], df[\"raw_VWC - 2h\"], df[\"raw_VWC - 3h\"], df[\"temp - 1h\"], df[\"temp - 2h\"], df[\"temp - 3h\"], df[\"tsd\"], df[\"hour\"]], axis = 1)\n",
        "\n",
        "#Electricity Data Omitted (For Type 1A and 2A Models)\n",
        "#X = pd.concat([df[\"EC - 1h\"], df[\"EC - 2h\"], df[\"EC - 3h\"], df[\"raw_VWC - 1h\"], df[\"raw_VWC - 2h\"], df[\"raw_VWC - 3h\"], df[\"temp - 1h\"], df[\"temp - 2h\"], df[\"temp - 3h\"], df[\"tsd\"], df[\"hour\"]], axis = 1)\n",
        "\n",
        "#Environmental Data Omitted (For Type 1B and 2B Models)\n",
        "#X = pd.concat([df[\"power - 1h\"], df[\"power - 2h\"], df[\"power - 3h\"], df[\"V1 - 1h\"], df[\"V1 - 2h\"], df[\"V1 - 3h\"], df[\"I1L - 1h\"], df[\"I1L - 2h\"], df[\"I1L - 3h\"]], axis = 1)\n",
        "\n",
        "y = pd.concat([df['Power (uW)'], df['Voltage (mV)'], df['Current (uA)']], axis = 1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "olQJQQB70EMm"
      },
      "outputs": [],
      "source": [
        "#All Data (For Type 1 and 2 Models)\n",
        "X_new = pd.concat([df_test[\"power - 1h\"], df_test[\"power - 2h\"], df_test[\"power - 3h\"], df_test[\"V1 - 1h\"], df_test[\"V1 - 2h\"], df_test[\"V1 - 3h\"], df_test[\"I1L - 1h\"], df_test[\"I1L - 2h\"], df_test[\"I1L - 3h\"],df_test[\"EC - 1h\"], df_test[\"EC - 2h\"], df_test[\"EC - 3h\"], df_test[\"raw_VWC - 1h\"], df_test[\"raw_VWC - 2h\"], df_test[\"raw_VWC - 3h\"], df_test[\"temp - 1h\"], df_test[\"temp - 2h\"], df_test[\"temp - 3h\"], df_test[\"tsd\"], df_test[\"hour\"]], axis = 1)\n",
        "\n",
        "#Electricity Data Omitted (For Type 1A and 2A Models)\n",
        "#X_new = pd.concat([df_test[\"EC - 1h\"], df_test[\"EC - 2h\"], df_test[\"EC - 3h\"], df_test[\"raw_VWC - 1h\"], df_test[\"raw_VWC - 2h\"], df_test[\"raw_VWC - 3h\"], df_test[\"temp - 1h\"], df_test[\"temp - 2h\"], df_test[\"temp - 3h\"], df_test[\"tsd\"], df_test[\"hour\"]], axis = 1)\n",
        "\n",
        "#Environmental Data Omitted (For Type 1B and 2B Models)\n",
        "#X_new = pd.concat([df_test[\"power - 1h\"], df_test[\"power - 2h\"], df_test[\"power - 3h\"], df_test[\"V1 - 1h\"], df_test[\"V1 - 2h\"], df_test[\"V1 - 3h\"], df_test[\"I1L - 1h\"], df_test[\"I1L - 2h\"], df_test[\"I1L - 3h\"]], axis = 1)\n",
        "\n",
        "y_new = pd.concat([df_test['Power (uW)'], df_test['Voltage (mV)'], df_test['Current (uA)']], axis = 1)"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "###  1.3 Runtime Simulation Code"
      ],
      "metadata": {
        "id": "peqVeLHjspvL"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#New runtime calculation\n",
        "import math\n",
        "from dateutil import parser\n",
        "from matplotlib import pyplot as plt\n",
        "from datetime import datetime, timedelta\n",
        "\n",
        "def internal_R_v3(R=2000): #return internal resistance of v3 cells in ohms\n",
        "    #https://www.jstage.jst.go.jp/article/jwet/20/1/20_21-087/_pdf\n",
        "    v0_oc = 48.5e-3 #48.5 mV\n",
        "    v0_cc = 4.8e-3\n",
        "    v0_r = R*((v0_oc/v0_cc)-1)\n",
        "\n",
        "    v1_oc = 43.8e-3\n",
        "    v1_cc = 20.9e-3\n",
        "    v1_r = R*((v1_oc/v1_cc)-1)\n",
        "\n",
        "    v2_oc = 45.2e-3\n",
        "    v2_cc = 23.5e-3\n",
        "    v2_r = R*((v2_oc/v2_cc)-1)\n",
        "\n",
        "    return (v0_r+v1_r+v2_r)/3\n",
        "\n",
        "def internal_R_v0(R=2000): #return internal resistance of v0 cells in ohms\n",
        "    v3_oc = 41.7e-3 #41.7mV\n",
        "    v3_cc = 5.1e-3\n",
        "    v3_r = R*((v3_oc/v3_cc)-1)\n",
        "\n",
        "    v4_oc = 48.7e-3\n",
        "    v4_cc = 16.8e-3\n",
        "    v4_r = R*((v4_oc/v4_cc)-1)\n",
        "\n",
        "    v5_oc = 39.1e-3\n",
        "    v5_cc = 16.9e-3\n",
        "    v5_r = R*((v5_oc/v5_cc)-1)\n",
        "\n",
        "    return (v3_r+v4_r+v5_r)/3\n",
        "\n",
        "def SMFC_current(v, R):\n",
        "    return v/R\n",
        "\n",
        "#MODEL\n",
        "def cap_leakage(E_cap_tn, timestep):\n",
        "    #Spec for KEMET T491\n",
        "    return 0.01e-6 * E_cap_tn * timestep\n",
        "\n",
        "def Matrix_Power(V, R):\n",
        "    #efficiency interpolated from https://www.analog.com/media/en/technical-documentation/data-sheets/ADP5091-5092.pdf\n",
        "    #given I_in = 100 uA and SYS = 3V\n",
        "    #V is the voltage (V) of the SMFC we captured\n",
        "    #R is the resistance (ohms) of the load we used to get that voltage trace\n",
        "    #Eta = -292.25665*V**4 + 784.30311*V**3 - 770.71691*V**2 + 342.00502*V + 15.83307\n",
        "    #Eta = Eta/100\n",
        "    Eta = 0.60\n",
        "    Pmax = (V**2)/R\n",
        "    Pout = Eta*Pmax\n",
        "    #assert((Eta > 0) & (Eta < 1))\n",
        "    #assert(Pout < 12000e-6)\n",
        "    return Pout\n",
        "\n",
        "def update_capEnergy(e0, V_applied, R, C, dt):\n",
        "    # e0: initial energy stored\n",
        "    # V_applied: voltage from SMFC\n",
        "    # R: internal resistance of SMFC\n",
        "    # C: capacitance of capacitor\n",
        "    # dt: time step since last data point\n",
        "    e_cap = e0 + Matrix_Power(V_applied, R)*dt - cap_leakage(e0, dt)\n",
        "    v_cap = math.sqrt(2*e_cap/C)\n",
        "    if e_cap < 0: #Not charging if leakage is greater than energy\n",
        "        e_cap = 0\n",
        "\n",
        "    return e_cap, v_cap #output final e and v\n",
        "\n",
        "def Advanced_energy():\n",
        "    #Now representing \"Advanced\"\n",
        "    #startup time of 2500 ms\n",
        "    t = 2500e-3\n",
        "    e = 2.4 * 128e-3 * t\n",
        "    e_startup = 2.4 * 128e-3 * 5e-3\n",
        "    return e+e_startup\n",
        "\n",
        "def Minimal_energy():\n",
        "    #Now representing \"Minimal\"\n",
        "    t = 0.888e-3 #tentative time\n",
        "    e = 0.9 * 4.8e-3 * t #this uses average current\n",
        "    e_startup = 0#assume negligible, no known startup time given\n",
        "    return  e + e_startup\n",
        "\n",
        "def Analog_energy():\n",
        "    #Now representing Analog\n",
        "    t = 1e-3 #estimated operating time\n",
        "    e = 0.11 * 2.15e-6 * t\n",
        "    e_startup = 0 #analog device, no startup needed :)\n",
        "    return e + e_startup\n",
        "\n",
        "#STEP 3:\n",
        "# For each day:\n",
        "#   on_Minimal, on_Advanced, on_Analog = 0\n",
        "#   For each time step (like every 60 s given our logging freq):\n",
        "#       - Update the energy in our capacitor (put fcn in models.py) given (1) input voltage, (2) time step, (3) capacitance (prob 10 uF), this will be an integral\n",
        "#       - Check if energy is enough to turn on (1) 1 uJ load, (2) 10 uJ load, and (3) 20 uJ load (will tweak later to reflect real energy cost of each system)\n",
        "#       - If so, add to on_Minimal, on_Advanced, and on_Analog and reset capacitor energy to 0 J (might tweak this value)\n",
        "#   Append on_Minimal, on_Advanced, on_Analog to on_Minimal_list, on_Advanced_list, on_Analog_list. This will be a list of how many sensor readings we are able to take with each of these systems every day given the energy we got\n",
        "#STEP 4: Visualize the daily # of readings with 3 bar graphs, y axis is # of readings and x axis is days.\n",
        "#   - Given 3 lists of integer values, plot them on bar graphs\n",
        "\n",
        "def group_util(test_date1, test_date2, N):\n",
        "    diff = (test_date2 - test_date1) / N\n",
        "    return [test_date1 + diff * idx for idx in range(N)] + [test_date2]\n",
        "\n",
        "def oracle_simulate(v_list, C_h, time_frame_seconds):\n",
        "    #Calculate maximum energy\n",
        "    total_E = 0\n",
        "    for i in range(len(v_list) - 1):\n",
        "        t = (v_list.index[i+1] - v_list.index[i]).total_seconds()\n",
        "        if t > 9000:\n",
        "          print(\"Discontinuity\")\n",
        "          print(v_list.index[i+1], v_list.index[i])\n",
        "          print(v_list['Voltage (mV)'][i+1], v_list['Voltage (mV)'][i])\n",
        "          #total_E, ignore = update_capEnergy(total_E, V_applied=(v_list['Voltage (mV)'][i+1] + v_list['Voltage (mV)'][i])/2, R=internal_R_v0(), C=C_h[0], dt = t)\n",
        "        else:\n",
        "          total_E, ignore = update_capEnergy(total_E, V_applied=max(v_list['Voltage (mV)'][i], v_list['Voltage (mV)'][i+1]), R=internal_R_v0(), C=C_h[0], dt = t)\n",
        "    print(\"Oracle activations:\", math.floor(total_E/Minimal_energy()))\n",
        "    return(math.floor(total_E/Minimal_energy()))\n",
        "\n",
        "def naive_simulate(t_list, v_list, v_list_naive, v_list_fine, C_h):\n",
        "    # t_list: list of decimal time stamps in unit of days (e.g. 71.85893518518519 day), same length as v_list\n",
        "    # v_list: list of voltage values from SFMC\n",
        "    # C_h: capacitance of the capacitor being filled up by harvester\n",
        "\n",
        "    #assume capacitor is completely discharged at start\n",
        "    e_minimal_stored = 0\n",
        "    e_minimal_stored_theo = 0\n",
        "\n",
        "    #Initialize evaluation metrics\n",
        "    false_act = 0\n",
        "    max_act = 0\n",
        "    pred_act = 0\n",
        "    succ_act = 0\n",
        "\n",
        "    total_E = 0\n",
        "    total_E_naive = 0\n",
        "\n",
        "    #Calculate maximum energy\n",
        "    #for i in range(len(v_list_fine) - 1):\n",
        "    #    t = (v_list_fine.index[i+1] - v_list_fine.index[i]).total_seconds()\n",
        "    #    total_E, ignore = update_capEnergy(total_E, V_applied=v_list_fine['V1 [10nV]'][i], R=internal_R_v0(), C=C_h[0], dt = t)\n",
        "    #print(total_E/Minimal_energy())\n",
        "    v = v_list_naive.mean()\n",
        "    #for each voltage data point\n",
        "    for jj in range(len(v_list) - 1): #last data point was at 71.85893518518519 day\n",
        "        t = (v_list.index[jj+1] - v_list.index[jj]).total_seconds()\n",
        "        if t <= time_frame_seconds:\n",
        "          #Total predicted vs. actual energy stored\n",
        "          #Predict energy stored during scheduled sub-interval\n",
        "          total_E, ignore = update_capEnergy(total_E, V_applied=v_list[jj], R=internal_R_v0(), C=C_h[0], dt = t)\n",
        "          total_E_naive, ignore = update_capEnergy(total_E_naive, V_applied=v, R=internal_R_v0(), C=C_h[0], dt = t)\n",
        "\n",
        "          E_Minimal_pred, v_minimal_pred = update_capEnergy(e_minimal_stored, V_applied=v, R=internal_R_v0(), C=C_h[0], dt = t) #set dt as length of prediction interval, in seconds\n",
        "          pred_act += math.floor(E_Minimal_pred/Minimal_energy()) #Update number of activations predicted\n",
        "          itn = 0\n",
        "          if math.floor(E_Minimal_pred/Minimal_energy()) > 0:\n",
        "              minimal_intervals = [date for date in group_util(v_list.index[jj], v_list.index[jj] + timedelta(seconds=t), math.floor(E_Minimal_pred/Minimal_energy()))]\n",
        "              #Calculate desired interval\n",
        "              int_len = time_frame_seconds /  math.floor(E_Minimal_pred/Minimal_energy())\n",
        "              for i in range(len(minimal_intervals) - 1):\n",
        "                  #Determine actual energy stored during scheduled sub-interval\n",
        "                  start = v_list_fine.index.searchsorted(minimal_intervals[i])\n",
        "                  end =  v_list_fine.index.searchsorted(minimal_intervals[i+1])\n",
        "\n",
        "                  E_Minimal, ignore = update_capEnergy(e_minimal_stored, V_applied=v_list_fine.iloc[start:end]['Voltage (mV)'].mean(), R=internal_R_v0(), C=C_h[0], dt = int_len)\n",
        "                  if not math.isnan(v_list_fine.iloc[start:end]['Voltage (mV)'].mean()):\n",
        "                    if E_Minimal < Minimal_energy():\n",
        "                        false_act += 1\n",
        "                        e_minimal_stored = max(0, E_Minimal - Minimal_energy())\n",
        "                        itn += 1\n",
        "\n",
        "                    elif E_Minimal >= Minimal_energy():\n",
        "                        succ_act += 1\n",
        "                        e_minimal_stored = max(0, E_Minimal - Minimal_energy())\n",
        "                        itn+= 1\n",
        "\n",
        "                    else:\n",
        "                      print('Error')\n",
        "                      print(e_minimal_stored, v)\n",
        "\n",
        "                  #Unit test\n",
        "                  #else:\n",
        "                  #  print(\"?\")\n",
        "                  #  print(v_list_fine.index[start])\n",
        "                  #  print(v_list_fine.index[end])\n",
        "                  #  print(minimal_intervals[i], minimal_intervals[i+1])\n",
        "\n",
        "              #Unit test\n",
        "              #if itn != math.floor(E_Minimal_pred/Minimal_energy()):\n",
        "              #    print(\"itn not matching\")\n",
        "              #    print(itn, math.floor(E_Minimal_pred/Minimal_energy()))\n",
        "              #    continue\n",
        "\n",
        "          else:\n",
        "              e_minimal_stored, ignore = update_capEnergy(e_minimal_stored, V_applied=v_list[jj], R=internal_R_v0(), C=C_h[0], dt = t)\n",
        "              #Added this\n",
        "              #start = v_list_fine.index.searchsorted(v_list.index[jj])\n",
        "              #end =  v_list_fine.index.searchsorted(v_list.index[jj+1])\n",
        "              #for h in range(start, end):\n",
        "              #    v = v_list_fine.iloc[h]['V1 [mV]']\n",
        "              #    interval_length = ((v_list_fine.index[h+1]) - (v_list_fine.index[h])).total_seconds()\n",
        "              #    E_Minimal, ignore = update_capEnergy(e_minimal_stored, V_applied=v, R=internal_R_v0(), C=C_h[0], dt = interval_length)\n",
        "              #    e_minimal_stored = E_Minimal\n",
        "\n",
        "\n",
        "        else:\n",
        "          print(\"It's over 9000!\", v_list.index[jj], v_list.index[jj+1])\n",
        "\n",
        "    print(\"Naive total_E activations:\", total_E/Minimal_energy())\n",
        "    print(\"Naive total_E_pred activations:\", total_E_naive/Minimal_energy())\n",
        "    return pred_act, false_act, succ_act, total_E_naive\n",
        "\n",
        "def getMax(c_list, input_list):\n",
        "    max_value = max(input_list)\n",
        "    i = [index for index, item in enumerate(input_list) if item == max_value][0]\n",
        "    return i, max_value, c_list[i]\n",
        "\n",
        "\n",
        "#SMFC\n",
        "import csv\n",
        "from collections import defaultdict\n",
        "from scipy.signal import butter, lfilter\n",
        "import matplotlib.pyplot as plt\n",
        "from datetime import datetime\n",
        "\n",
        "def butter_lowpass(cutoff, fs, order=5):\n",
        "        return butter(order, cutoff, fs=fs, btype='low', analog=False)\n",
        "\n",
        "def butter_lowpass_filter(data, cutoff, fs, order=5):\n",
        "    b, a = butter_lowpass(cutoff, fs, order=order)\n",
        "    y = lfilter(b, a, data)\n",
        "    return y\n",
        "\n",
        "def getMFC_data(y_test, test_pred):\n",
        "    unix_time = y_test.index\n",
        "    d0 = unix_time[0]\n",
        "    days = []\n",
        "    for d in unix_time:\n",
        "        day = d\n",
        "        day_from_start = day-d0\n",
        "        decimal_day = day_from_start.total_seconds()/(24 * 3600)\n",
        "        days.append(decimal_day)\n",
        "\n",
        "    return days\n",
        "\n",
        "def simulate(t_list, v_list, v_list_pred, v_list_fine, C_h):\n",
        "    # t_list: list of decimal time stamps in unit of days (e.g. 71.85893518518519 day), same length as v_list\n",
        "    # v_list: list of voltage values from SFMC\n",
        "    # C_h: capacitance of the capacitor being filled up by harvester\n",
        "\n",
        "    #assume capacitor is completely discharged at start\n",
        "    e_minimal_stored = 0\n",
        "    e_minimal_stored_theo = 0\n",
        "\n",
        "    #Initialize evaluation metrics\n",
        "    false_act = 0\n",
        "    max_act = 0\n",
        "    pred_act = 0\n",
        "    succ_act = 0\n",
        "\n",
        "    total_E = 0\n",
        "    total_E_pred = 0\n",
        "\n",
        "    #Calculate maximum energy\n",
        "    #for i in range(len(v_list_fine) - 1):\n",
        "    #    t = (v_list_fine.index[i+1] - v_list_fine.index[i]).total_seconds()\n",
        "    #    total_E, ignore = update_capEnergy(total_E, V_applied=v_list_fine['V1 [10nV]'][i], R=internal_R_v0(), C=C_h[0], dt = t)\n",
        "    #print(total_E/Minimal_energy())\n",
        "    #for each voltage data point\n",
        "    for jj in range(len(v_list) - 1): #last data point was at 71.85893518518519 day\n",
        "        t = (v_list.index[jj+1] - v_list.index[jj]).total_seconds()\n",
        "        total_E, ignore = update_capEnergy(total_E, V_applied=v_list[jj], R=internal_R_v0(), C=C_h[0], dt = t)\n",
        "        total_E_pred, ignore = update_capEnergy(total_E_pred, V_applied=v_list_pred[jj], R=internal_R_v0(), C=C_h[0], dt = t)\n",
        "        if t <= time_frame_seconds:\n",
        "          #Total predicted vs. actual energy stored\n",
        "          #Predict energy stored during scheduled sub-interval\n",
        "          E_Minimal_pred, v_minimal_pred = update_capEnergy(e_minimal_stored, V_applied=v_list_pred[jj], R=internal_R_v0(), C=C_h[0], dt = t) #set dt as length of prediction interval, in seconds\n",
        "          pred_act += math.floor(E_Minimal_pred/Minimal_energy()) #Update number of activations predicted\n",
        "          itn = 0\n",
        "          if math.floor(E_Minimal_pred/Minimal_energy()) > 0:\n",
        "              minimal_intervals = [date for date in group_util(v_list_pred.index[jj], v_list_pred.index[jj] + timedelta(seconds=t), math.floor(E_Minimal_pred/Minimal_energy()))]\n",
        "              #Calculate desired interval\n",
        "              int_len = time_frame_seconds /  math.floor(E_Minimal_pred/Minimal_energy())\n",
        "              for i in range(len(minimal_intervals) - 1):\n",
        "                  #Determine actual energy stored during scheduled sub-interval\n",
        "                  start = v_list_fine.index.searchsorted(minimal_intervals[i])\n",
        "                  end =  v_list_fine.index.searchsorted(minimal_intervals[i+1])\n",
        "                  v = v_list_fine.iloc[start:end]['Voltage (mV)'].mean()\n",
        "\n",
        "                  #interval_length = ((v_list_fine.index[end]) - (v_list_fine.index[start])).total_seconds()\n",
        "                  #if interval_length > int_len:\n",
        "                  #  print('interval_length > int_len')\n",
        "                  #  print('interval_length, int_len:', interval_length, int_len)\n",
        "                  #  print(v_list_fine.index[start], v_list_fine.index[end])\n",
        "                  #else:\n",
        "                  #  print('interval_length <= int_len')\n",
        "                  #  print('interval_length, int_len:', interval_length, int_len)\n",
        "                  #  print(v_list_fine.index[start], v_list_fine.index[end])\n",
        "\n",
        "                  E_Minimal, ignore = update_capEnergy(e_minimal_stored, V_applied=v, R=internal_R_v0(), C=C_h[0], dt = int_len)\n",
        "                  if not math.isnan(v_list_fine.iloc[start:end]['Voltage (mV)'].mean()):\n",
        "                    if E_Minimal < Minimal_energy():\n",
        "                        false_act += 1\n",
        "                        e_minimal_stored = max(0, E_Minimal - Minimal_energy())\n",
        "                        itn += 1\n",
        "\n",
        "                    elif E_Minimal >= Minimal_energy():\n",
        "                        succ_act += 1\n",
        "                        e_minimal_stored = max(0, E_Minimal - Minimal_energy())\n",
        "                        itn+= 1\n",
        "\n",
        "                    else:\n",
        "                      print('Error')\n",
        "                      print(e_minimal_stored, v)\n",
        "\n",
        "                  #Unit test\n",
        "                  #else:\n",
        "                  #  print(\"?\")\n",
        "                  #  print(v_list_fine.index[start])\n",
        "                  #  print(v_list_fine.index[end])\n",
        "                  #  print(minimal_intervals[i], minimal_intervals[i+1])\n",
        "\n",
        "              #Unit test\n",
        "              #if itn != math.floor(E_Minimal_pred/Minimal_energy()):\n",
        "              #    print(\"itn not matching\")\n",
        "              #    print(itn, math.floor(E_Minimal_pred/Minimal_energy()))\n",
        "              #    continue\n",
        "\n",
        "          else:\n",
        "              e_minimal_stored, ignore = update_capEnergy(e_minimal_stored, V_applied=v_list[jj], R=internal_R_v0(), C=C_h[0], dt = t)\n",
        "              #Added this\n",
        "              #start = v_list_fine.index.searchsorted(v_list.index[jj])\n",
        "              #end =  v_list_fine.index.searchsorted(v_list.index[jj+1])\n",
        "              #for h in range(start, end):\n",
        "              #    v = v_list_fine.iloc[h]['V1 [mV]']\n",
        "              #    interval_length = ((v_list_fine.index[h+1]) - (v_list_fine.index[h])).total_seconds()\n",
        "              #    E_Minimal, ignore = update_capEnergy(e_minimal_stored, V_applied=v, R=internal_R_v0(), C=C_h[0], dt = interval_length)\n",
        "              #    e_minimal_stored = E_Minimal\n",
        "\n",
        "\n",
        "        else:\n",
        "          print(\"It's over 9000!\", v_list.index[jj], v_list.index[jj+1])\n",
        "\n",
        "    print(\"Runtime total_E activations:\", total_E/Minimal_energy())\n",
        "    print(\"Runtime total_E_pred activations:\", total_E_pred/Minimal_energy())\n",
        "    return pred_act, false_act, succ_act, total_E, total_E_pred"
      ],
      "metadata": {
        "id": "2ejNP9wmsvw5"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "##  2. Load Models"
      ],
      "metadata": {
        "id": "gubEvcHAeOfn"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "###  2.1 Load pretrained Type 1 Models, validate on Dataset 1\n",
        "###  In order to use pretrained models it is neccesary to download the ```trained_models``` directory from [Hugging Face](https://huggingface.co/datasets/adunlop621/Soil_MFC/tree/main) and store it in the same directory as this notebook"
      ],
      "metadata": {
        "id": "PcqsZz8zz8Az"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from keras.models import load_model\n",
        "from sklearn.model_selection import TimeSeriesSplit\n",
        "from keras.models import Sequential\n",
        "from keras.layers import Dense\n",
        "from keras.layers import LSTM\n",
        "from keras import backend as K\n",
        "\n",
        "power_mape = []\n",
        "voltage_mape = []\n",
        "current_mape = []\n",
        "\n",
        "E_actual_list = []\n",
        "E_pred_list = []\n",
        "\n",
        "max_act_list = []\n",
        "pred_act_list = []\n",
        "succ_act_list = []\n",
        "\n",
        "pred_act_naive_list = []\n",
        "false_act_naive_list = []\n",
        "succ_act_naive_list = []\n",
        "\n",
        "#Set parameters\n",
        "batchsize_list = [300, 150, 50, 20, 8]\n",
        "time_frame_list = ['3min', '5min', '15min', '30min', '60min']\n",
        "time_frame_seconds_list = [180, 300, 900, 1800, 3600]\n",
        "n = 0\n",
        "\n",
        "for j in range(len(batchsize_list)):\n",
        "  if j >=0:\n",
        "\n",
        "    X_train, X_test = train_test_split(X, test_size=0.3, shuffle=False)\n",
        "    y_train, y_test = train_test_split(y, test_size=0.3, shuffle=False)\n",
        "\n",
        "    X_valid, X_test = train_test_split(X_test, test_size=0.5, shuffle=False)\n",
        "    y_valid, y_test = train_test_split(y_test, test_size=0.5, shuffle=False)\n",
        "\n",
        "    batchsize = batchsize_list[j]\n",
        "    time_frame = time_frame_list[j]\n",
        "    time_frame_seconds = time_frame_seconds_list[j]\n",
        "\n",
        "    E_actual = 0\n",
        "    for i in range(len(y_test) - 1):\n",
        "      t = (y_test.index[i+1] - y_test.index[i]).total_seconds()\n",
        "      if t < 180:\n",
        "        E_actual += y_test['Power (uW)'][i] * t\n",
        "\n",
        "    #Set dataset bounds\n",
        "    train_bound_lower = y_train.index[0]\n",
        "    train_bound_upper = y_train.index[-1]\n",
        "    valid_bound_lower = y_valid.index[0]\n",
        "    valid_bound_upper = y_valid.index[-1]\n",
        "    test_bound_lower = y_test.index[0]\n",
        "    test_bound_upper = y_test.index[-1]\n",
        "\n",
        "    v_test = y_test.loc[(((y_test.index >= test_bound_lower)) & (y_test.index <= test_bound_upper))]['Voltage (mV)']\n",
        "    #v_test = v_test.drop(v_test[(v_test.index > '2021-11-11') & (v_test.index < '2021-11-22 01:00:00')].index)\n",
        "    v_test = pd.DataFrame(v_test)/1E5\n",
        "    v_avg_true = v_test['Voltage (mV)'].resample(time_frame).mean().dropna()\n",
        "    C0 = [0.007000000000000006, 0.007000000000000006, 0.007000000000000006]\n",
        "\n",
        "\n",
        "    #Resample data\n",
        "    X_train = X_train.resample(time_frame).mean().dropna()\n",
        "    X_valid = X_valid.resample(time_frame).mean().dropna()\n",
        "    X_test = X_test.resample(time_frame).mean().dropna()\n",
        "\n",
        "    y_train = y_train.resample(time_frame).mean().dropna()\n",
        "    y_valid = y_valid.resample(time_frame).mean().dropna()\n",
        "    y_test = y_test.resample(time_frame).mean().dropna()\n",
        "\n",
        "    print(v_test)\n",
        "\n",
        "    #Reshape data\n",
        "    X_train = X_train.values.reshape((X_train.shape[0], 1, X_train.shape[1]))\n",
        "    X_valid = X_valid.values.reshape((X_valid.shape[0], 1, X_valid.shape[1]))\n",
        "    X_test = X_test.values.reshape((X_test.shape[0], 1, X_test.shape[1]))\n",
        "\n",
        "    def quantile_loss(y_true, y_pred, quantile = 0.05):\n",
        "      error = y_true - y_pred\n",
        "      return tf.reduce_mean(K.maximum(quantile * error, (quantile - 1) * error), axis=-1)\n",
        "\n",
        "    # design network\n",
        "    model = Sequential()\n",
        "    model.add(LSTM(200, input_shape=(X_train.shape[1], X_train.shape[2]), activation='relu'))\n",
        "    model.add(Dense(100, activation='relu'))\n",
        "    model.add(Dense(3))\n",
        "    model.compile(loss=quantile_loss, metrics=['mape'], optimizer='adam')\n",
        "\n",
        "    #Select model to load\n",
        "    file = 'trained_models/type1B_' + time_frame_list[j] + '_quant50'\n",
        "    print(file)\n",
        "    model = load_model(file, custom_objects={'quantile_loss': quantile_loss})\n",
        "\n",
        "\n",
        "    train_pred = model.predict(X_train, batch_size=batchsize)\n",
        "    test_pred = model.predict(X_test, batch_size=batchsize)\n",
        "\n",
        "    days  = getMFC_data(y_test, test_pred)\n",
        "\n",
        "    #Prepare data for runtime simulation\n",
        "    y_test['power pred'] = test_pred[:, 0]\n",
        "    y_test['Voltage (mV) pred'] = test_pred[:, 1]\n",
        "    y_test['Current (uA) pred'] = test_pred[:, 2]\n",
        "\n",
        "    v_avg_pred = y_test['Voltage (mV) pred']/1E5\n",
        "    days  = getMFC_data(y_test, test_pred)\n",
        "\n",
        "\n",
        "    #Remove first and last entries of averaged data to prevent overestimation of available energy\n",
        "    v_avg_true = v_avg_true[1:][:-1]\n",
        "    v_avg_pred = v_avg_pred[1:][:-1]\n",
        "\n",
        "\n",
        "    #Run oracle model\n",
        "    max_act = oracle_simulate(v_test, C0, time_frame_seconds)\n",
        "\n",
        "    #Call simulate function\n",
        "    pred_act, false_act, succ_act, total_E, total_E_pred = simulate(days, v_avg_true, v_avg_pred, v_test, C0)\n",
        "\n",
        "    #Run naive model\n",
        "    v_valid = y_valid.loc[(((y_valid.index >= valid_bound_lower)) & (y_valid.index < valid_bound_upper))]['Voltage (mV)']/1E5\n",
        "    pred_act_naive, false_act_naive, succ_act_naive, total_E = naive_simulate(days, v_avg_true, v_valid, v_test, C0)\n",
        "    print(\"Dataset, train set, and test set size:\", len(y_train) + len(y_valid) + len(y_test), len(y_train), len(y_test))\n",
        "    print('Timeframe:', time_frame)\n",
        "\n",
        "    print('Minimal Application')\n",
        "    if succ_act_naive > 0:\n",
        "      print(\"Naive vs. DL succesful activations:\", succ_act/succ_act_naive)\n",
        "    else:\n",
        "      print(\"No succesful naive activation\")\n",
        "\n",
        "    #print('Predicted vs. Actual percent difference: %.3f%%' % ((total_E * 100 / total_E_pred) - 100))\n",
        "    print('Maximum possible activations:', max_act)\n",
        "    print('Predicted activations:', pred_act)\n",
        "    print('Successful activations: %d, %.3f%%' % (succ_act, succ_act * 100/pred_act))\n",
        "    print('Failed activations: %d, %.3f%%' % (false_act, false_act * 100/pred_act))\n",
        "    print('Missed activations: %d, %.3f%%' % (max_act - succ_act, (max_act - succ_act) * 100/max_act))\n",
        "    print('Successful activations compared to max: %d, %.3f%%' % (succ_act, succ_act * 100/max_act))\n",
        "\n",
        "    #Naive model\n",
        "    print('Naive predicted activations (usual actual energy average):', pred_act_naive)\n",
        "    print('Naive successful activations (usual actual energy average): %d, %.3f%%' % (succ_act_naive, succ_act_naive * 100/pred_act_naive))\n",
        "    print('Naive failed activations (usual actual energy average): %d, %.3f%%' % (false_act_naive, false_act_naive * 100/pred_act_naive))\n",
        "    print('Naive missed activations (usual actual energy average): %d, %.3f%%' % (max_act - succ_act_naive, (max_act - succ_act_naive) * 100/max_act))\n",
        "\n",
        "    print('Naive successful activations (compared to max): %d, %.3f%%' % (succ_act_naive, succ_act_naive * 100/max_act))\n",
        "\n",
        "    print('Voltage overestimation rate: %.3f%%' % ((y_test['Voltage (mV)'].values <= y_test['Voltage (mV) pred']).mean() * 100))\n",
        "    print(\"Test MAPE power: %3f\" %  MAPE(y_test['Power (uW)'].values.ravel(), y_test['power pred']))\n",
        "    print(\"Test MAPE voltage: %3f\" % MAPE(y_test['Voltage (mV)'], y_test['Voltage (mV) pred']))\n",
        "    print(\"Test MAPE current: %3f\" % MAPE(y_test['Current (uA)'], y_test['Current (uA) pred']))\n",
        "\n",
        "    E_pred = 0\n",
        "    for i in range(len(y_test) - 1):\n",
        "      t = (y_test.index[i+1] - y_test.index[i]).total_seconds()\n",
        "      if t <= time_frame_seconds + 50:\n",
        "        E_pred += y_test['power pred'][i] * t\n",
        "\n",
        "    print('Predicted vs. Actual Total Energy Percent Difference: %.3f%%' % ((E_pred - E_actual) * 100 / E_actual))\n",
        "    print(\"Test MAPE power:\\n\", MAPE(y_test['Power (uW)'].values.ravel(), y_test['power pred'].values.ravel()))\n",
        "\n",
        "    V_actual = y_test['Voltage (mV)'].mean()\n",
        "    V_pred = y_test['Voltage (mV) pred'].mean()\n",
        "    print('Predicted vs. Actual Total Voltage Percent Difference: %.3f%%' % ((V_pred - V_actual) * 100 / V_actual))\n",
        "    print(\"Test MAPE voltage:\\n\", MAPE(y_test['Voltage (mV)'].values.ravel(), y_test['Voltage (mV) pred'].values.ravel()))\n",
        "\n",
        "\n",
        "\n",
        "    print('Voltage overestimation rate: %.3f%%' % ((y_test['Voltage (mV)'].values <= y_test['Voltage (mV) pred']).mean() * 100))\n",
        "    print(\"Test MAPE power: %3f\" %  MAPE(y_test['Power (uW)'].values.ravel(), y_test['power pred']))\n",
        "    print(\"Test MAPE voltage: %3f\" % MAPE(y_test['Voltage (mV)'], y_test['Voltage (mV) pred']))\n",
        "    print(\"Test MAPE current: %3f\" % MAPE(y_test['Current (uA)'], y_test['Current (uA) pred']))"
      ],
      "metadata": {
        "id": "LiB28Wr0z_MP"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Graph Selected Data"
      ],
      "metadata": {
        "id": "KFLthf114Z-F"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from keras.models import load_model\n",
        "\n",
        "#Select pretrained models to graph\n",
        "\n",
        "def quantile_loss(y_true, y_pred, quantile = 0.5):\n",
        "    error = y_true - y_pred\n",
        "    return tf.reduce_mean(K.maximum(quantile * error, (quantile - 1) * error), axis=-1)\n",
        "\n",
        "model = load_model(\"trained_models/type1A_60min_quant50\", custom_objects={'quantile_loss': quantile_loss})\n",
        "predictions = model.predict(np.concatenate((X_train, X_valid, X_test)))\n",
        "mv1 = df\n",
        "mv1[\"power_pred_med\"] = predictions[:, 0]\n",
        "mv1[\"voltage_pred_med\"] = predictions[:, 1]\n",
        "mv1[\"current_pred_med\"] = predictions[:, 2]\n",
        "\n",
        "def quantile_loss(y_true, y_pred, quantile = 0.05):\n",
        "    error = y_true - y_pred\n",
        "    return tf.reduce_mean(K.maximum(quantile * error, (quantile - 1) * error), axis=-1)\n",
        "\n",
        "model = load_model(\"trained_models/type1A_60min_quant5\", custom_objects={'quantile_loss': quantile_loss})\n",
        "predictions = model.predict(np.concatenate((X_train, X_valid, X_test)))\n",
        "mv1 = df\n",
        "mv1[\"power_pred_lower\"] = predictions[:, 0]\n",
        "mv1[\"voltage_pred_lower\"] = predictions[:, 1]\n",
        "mv1[\"current_pred_lower\"] = predictions[:, 2]\n",
        "\n",
        "def quantile_loss(y_true, y_pred, quantile = 0.75):\n",
        "    error = y_true - y_pred\n",
        "    return tf.reduce_mean(K.maximum(quantile * error, (quantile - 1) * error), axis=-1)\n",
        "\n",
        "model = load_model(\"trained_models/type1A_60min_quant75\", custom_objects={'quantile_loss': quantile_loss})\n",
        "predictions = model.predict(np.concatenate((X_train, X_valid, X_test)))\n",
        "mv1 = df\n",
        "mv1[\"power_pred_upper\"] = predictions[:, 0]\n",
        "mv1[\"voltage_pred_upper\"] = predictions[:, 1]\n",
        "mv1[\"current_pred_upper\"] = predictions[:, 2]\n",
        "\n",
        "\n",
        "mv1 = mv1.loc[(mv1.index > '2021-12-12') & (mv1.index < '2021-12-14')]\n",
        "mv2 = mv1"
      ],
      "metadata": {
        "id": "oW04Uicd4e05"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import matplotlib as mpl\n",
        "mpl.use('Agg')\n",
        "#mpl.rc('font', **font)\n",
        "import matplotlib.pyplot as plt\n",
        "import matplotlib.gridspec as gridspec\n",
        "import matplotlib.dates as md\n",
        "import datetime\n",
        "import numpy as np\n",
        "from pytz import timezone\n",
        "import pandas as pd\n",
        "import arrow\n",
        "import os\n",
        "from pandas.plotting import register_matplotlib_converters\n",
        "register_matplotlib_converters()\n",
        "%matplotlib inline\n",
        "\n",
        "# Limits for graphs\n",
        "VOLTAGE_LIM = 0.2\n",
        "CURRENT_LIM = 40\n",
        "POWER_LIM = 4\n",
        "\n",
        "line_width = 0.5\n",
        "plt.title(\"1 hour prediction interval\")\n",
        "\n",
        "plt.close()\n",
        "plt.xlabel(\"Time\")\n",
        "fig, (ax1, ax2, ax3) = plt.subplots(3,figsize=(4,3), sharex=False)\n",
        "fig.autofmt_xdate()\n",
        "\n",
        "volt_color= 'tab:blue'\n",
        "\n",
        "amp_color = 'tab:red'\n",
        "\n",
        "\n",
        "volt_color1= 'tab:blue'\n",
        "volt_style1 = 'dashed'\n",
        "volt_color2= 'tab:green'\n",
        "volt_style2 = 'dashed'\n",
        "volt_color3= 'tab:red'\n",
        "volt_style3 = 'dashed'\n",
        "volt_color4= 'tab:orange'\n",
        "volt_style4 = 'dashed'\n",
        "\n",
        "#amp_color1 = 'tab:red'\n",
        "#amp_style1='dashed'\n",
        "#amp_color2 = 'tab:orange'\n",
        "#amp_style2='dashdot'\n",
        "\n",
        "ax1.fmt_xdata = md.DateFormatter('%m-%d')\n",
        "ax1.xaxis.set_major_formatter(md.DateFormatter('%m-%d'))\n",
        "ax1.grid(True)\n",
        "\n",
        "ax2.fmt_xdata = md.DateFormatter('%m-%d')\n",
        "ax2.xaxis.set_major_formatter(md.DateFormatter('%m-%d'))\n",
        "ax2.grid(True)\n",
        "\n",
        "ax3.fmt_xdata = md.DateFormatter('%m-%d')\n",
        "ax3.xaxis.set_major_formatter(md.DateFormatter('%m-%d'))\n",
        "ax3.grid(True)\n",
        "\n",
        "ax1.set_ylabel(\"Power (μW)\", fontsize=7.5, labelpad = 2.5)\n",
        "ax2.set_ylabel(\"Voltage (mV)\", fontsize=7.5, labelpad = 2.5)\n",
        "ax3.set_ylabel(\"Current (μA)\", fontsize=7.5, labelpad = 2.5)\n",
        "#ax3.set_xlabel(\"Date\", fontsize=10)\n",
        "\n",
        "ax1.set_ylim(0, 1.5)\n",
        "ax2.set_ylim(0, 60)\n",
        "ax3.set_ylim(0, 40)\n",
        "\n",
        "#ax3.set_xlim([datetime.date(2021, 12, 12), datetime.date(2021, 12, 14)])\n",
        "ax1.set_xlim([mv1.index[0], mv1.index[-1]])\n",
        "ax2.set_xlim([mv1.index[0], mv1.index[-1]])\n",
        "ax3.set_xlim([mv1.index[0], mv1.index[-1]])\n",
        "\n",
        "ax1.plot(mv1.index, mv2[\"power_pred_lower\"] * 1E-3, color='tab:red', ls = volt_style3, linewidth=line_width)\n",
        "ax1.plot(mv1.index, mv2[\"power_pred_upper\"] * 1E-3, color='tab:orange', ls = volt_style4, linewidth=line_width)\n",
        "ax1.plot(mv1.index, mv2[\"power\"] * 1E-3, color='tab:blue', ls = 'solid', linewidth=line_width)\n",
        "ax1.plot(mv1.index, mv2[\"power_pred_med\"] * 1E-3, color='tab:green', ls = volt_style2, linewidth=line_width)\n",
        "ax1.fill_between(mv1.index, mv2[\"power_pred_lower\"] * 1E-3, mv2[\"power_pred_upper\"] * 1E-3, color='grey', alpha=0.5)\n",
        "#ax1.legend(['lower bound', 'upper bound', 'actual', 'median predictions'], loc='lower center', prop={'size': 5.5}, ncol=2)\n",
        "\n",
        "ax2.plot(mv1.index, mv2[\"voltage_pred_lower\"] * 1E-2, color='tab:red', ls = volt_style3, linewidth=line_width)\n",
        "ax2.plot(mv1.index, mv2[\"voltage_pred_upper\"] * 1E-2, color='tab:orange', ls = volt_style4, linewidth=line_width)\n",
        "ax2.plot(mv1.index, mv2['V1 [mV]'] * 1E-2, color='tab:blue', ls = 'solid', linewidth=line_width)\n",
        "ax2.plot(mv1.index, mv2[\"voltage_pred_med\"] * 1E-2, color='tab:green', ls = volt_style2, linewidth=line_width)\n",
        "ax2.fill_between(mv1.index, mv2[\"voltage_pred_lower\"] * 1E-2, mv2[\"voltage_pred_upper\"] * 1E-2, color='grey', alpha=0.5)\n",
        "ax2.legend(['lower bound', 'upper bound', 'ground truth', 'median prediction'], loc='lower left', prop={'size': 6.6}, ncol=2, columnspacing=0.5)\n",
        "\n",
        "ax3.plot(mv1.index, mv2[\"current_pred_lower\"] * 1E-2, color='tab:red', ls = volt_style3, linewidth=line_width)\n",
        "ax3.plot(mv1.index, mv2[\"current_pred_upper\"] * 1E-2, color='tab:orange', ls = volt_style4, linewidth=line_width)\n",
        "ax3.plot(mv1.index, mv2['I1L [μA]'] * 1E-2, color='tab:blue', ls = 'solid', linewidth=line_width)\n",
        "ax3.plot(mv1.index, mv2[\"current_pred_med\"] * 1E-2, color='tab:green', ls = volt_style2, linewidth=line_width)\n",
        "ax3.fill_between(mv1.index, mv2[\"current_pred_lower\"] * 1E-2, mv2[\"current_pred_upper\"] * 1E-2, color='grey', alpha=0.5)\n",
        "#ax3.legend(['lower bound', 'upper bound', 'actual', 'median predictions'], loc='upper right', prop={'size': 5.5}, ncol=1)\n",
        "\n",
        "\n",
        "#Plot error\n",
        "#ax3.plot(mv1['timestamp'], mv1['error']/mv1['power'], color=volt_color2, ls = volt_style2)\n",
        "#ax3.legend(['error'], loc='upper right', prop={'size': 6})\n",
        "\n",
        "ax3.tick_params(axis='x', labelsize=7.5, rotation=0, pad = 0.1)\n",
        "ax3.set_xticks(list(ax3.get_xticks()) + [ax3.get_xlim()[0], ax3.get_xlim()[1]])\n",
        "for label in ax3.get_xticklabels():\n",
        "    label.set_horizontalalignment('center')\n",
        "\n",
        "ax1.tick_params(axis='y', labelsize=7.5, rotation=0, pad = 0.1)\n",
        "ax2.tick_params(axis='y', labelsize=7.5, rotation=0, pad = 0.1)\n",
        "ax3.tick_params(axis='y', labelsize=7.5, rotation=0, pad = 0.1)\n",
        "\n",
        "plt.tight_layout(pad=0.3, w_pad=0.5, h_pad=0.1)\n",
        "#plt.subplots_adjust(hspace=0.15)\n",
        "plt.savefig('twobat.pdf')\n",
        "plt.close()"
      ],
      "metadata": {
        "id": "TQQu1Odj4iZ1"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "###  2.2 Load pretrained Type 2 Models, validate on Dataset 2\n",
        "\n",
        "###  In order to use pretrained models it is neccesary to download the ```trained_models``` directory from [Hugging Face](https://huggingface.co/datasets/adunlop621/Soil_MFC/tree/main) and store it in the same directory as this notebook"
      ],
      "metadata": {
        "id": "umLPmZ1-uH_b"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "\n",
        "from keras.models import load_model\n",
        "from sklearn.model_selection import TimeSeriesSplit\n",
        "from keras.models import Sequential\n",
        "from keras.layers import Dense\n",
        "from keras.layers import LSTM\n",
        "from keras import backend as K\n",
        "\n",
        "power_mape = []\n",
        "voltage_mape = []\n",
        "current_mape = []\n",
        "\n",
        "E_actual_list = []\n",
        "E_pred_list = []\n",
        "\n",
        "max_act_list = []\n",
        "pred_act_list = []\n",
        "succ_act_list = []\n",
        "\n",
        "pred_act_naive_list = []\n",
        "false_act_naive_list = []\n",
        "succ_act_naive_list = []\n",
        "\n",
        "#Set parameters\n",
        "batchsize_list = [300, 150, 50, 20, 8]\n",
        "time_frame_list = ['3min', '5min', '15min', '30min', '60min']\n",
        "time_frame_seconds_list = [180, 300, 900, 1800, 3600]\n",
        "n = 0\n",
        "\n",
        "for j in range(len(batchsize_list)):\n",
        "  if j>=0:\n",
        "    batchsize = batchsize_list[j]\n",
        "    time_frame = time_frame_list[j]\n",
        "    time_frame_seconds = time_frame_seconds_list[j]\n",
        "\n",
        "    #Split train and test sets\n",
        "    X_valid, X_test = train_test_split(X_new, test_size=0.5, shuffle=False)\n",
        "    y_valid, y_test = train_test_split(y_new, test_size=0.5, shuffle=False)\n",
        "\n",
        "    X_train = X\n",
        "    y_train = y\n",
        "\n",
        "    E_actual = 0\n",
        "    for i in range(len(y_test) - 1):\n",
        "      t = (y_test.index[i+1] - y_test.index[i]).total_seconds()\n",
        "      if t < 180:\n",
        "        E_actual += y_test['Power (uW)'][i] * t\n",
        "\n",
        "    #Set dataset bounds\n",
        "    train_bound_lower = y_train.index[0]\n",
        "    train_bound_upper = y_train.index[-1]\n",
        "    valid_bound_lower = y_valid.index[0]\n",
        "    valid_bound_upper = y_valid.index[-1]\n",
        "    test_bound_lower = y_test.index[0]\n",
        "    test_bound_upper = y_test.index[-1]\n",
        "\n",
        "    #Resample data\n",
        "    X_train = X_train.resample(time_frame).mean().dropna()\n",
        "    X_valid = X_valid.resample(time_frame).mean().dropna()\n",
        "    X_test = X_test.resample(time_frame).mean().dropna()\n",
        "\n",
        "    y_train = y_train.resample(time_frame).mean().dropna()\n",
        "    y_valid = y_valid.resample(time_frame).mean().dropna()\n",
        "    y_test = y_test.resample(time_frame).mean().dropna()\n",
        "\n",
        "    #Reshape data\n",
        "    X_train = X_train.values.reshape((X_train.shape[0], 1, X_train.shape[1]))\n",
        "    X_valid = X_valid.values.reshape((X_valid.shape[0], 1, X_valid.shape[1]))\n",
        "    X_test = X_test.values.reshape((X_test.shape[0], 1, X_test.shape[1]))\n",
        "\n",
        "\n",
        "\n",
        "    def quantile_loss(y_true, y_pred, quantile = 0.05):\n",
        "      error = y_true - y_pred\n",
        "      return tf.reduce_mean(K.maximum(quantile * error, (quantile - 1) * error), axis=-1)\n",
        "\n",
        "    file = 'trained_models/type2_' + time_frame_list[j] + '_quant5'\n",
        "    print(file)\n",
        "    model = load_model(file, custom_objects={'quantile_loss': quantile_loss})\n",
        "    train_pred = model.predict(X_train, batch_size=batchsize)\n",
        "    test_pred = model.predict(X_test, batch_size=batchsize)\n",
        "\n",
        "\n",
        "    #Prepare data for runtime simulation\n",
        "    y_test['power pred'] = test_pred[:, 0]\n",
        "    y_test['Voltage (mV) pred'] = test_pred[:, 1]\n",
        "    y_test['Current (uA) pred'] = test_pred[:, 2]\n",
        "\n",
        "    days  = getMFC_data(y_test, test_pred)\n",
        "\n",
        "    v_test = y_new.loc[(((y_new.index >= test_bound_lower)) & (y_new.index <= test_bound_upper))]['Voltage (mV)']\n",
        "    #v_test = v_test.drop(v_test[(v_test.index > '2021-11-11') & (v_test.index < '2021-11-22 01:00:00')].index)\n",
        "    v_test = pd.DataFrame(v_test)/1E4\n",
        "    v_avg_true = v_test['Voltage (mV)'].resample(time_frame).mean().dropna()\n",
        "    v_avg_pred = y_test['Voltage (mV) pred']/1E4\n",
        "    C0 = [0.007000000000000006, 0.007000000000000006, 0.007000000000000006]\n",
        "\n",
        "    #Remove first and last entries of averaged data to prevent overestimation of available energy\n",
        "    v_avg_true = v_avg_true[1:][:-1]\n",
        "    v_avg_pred = v_avg_pred[1:][:-1]\n",
        "\n",
        "    #Run oracle model\n",
        "    max_act = oracle_simulate(v_test, C0, time_frame_seconds)\n",
        "\n",
        "    #Call simulate function\n",
        "    pred_act, false_act, succ_act, total_E, total_E_pred = simulate(days, v_avg_true, v_avg_pred, v_test, C0)\n",
        "\n",
        "    #Run naive model\n",
        "    v_valid = y_new.loc[(((y_new.index >= valid_bound_lower)) & (y_new.index < valid_bound_upper))]['Voltage (mV)']/1E4\n",
        "    pred_act_naive, false_act_naive, succ_act_naive, total_E = naive_simulate(days, v_avg_true, v_valid, v_test, C0)\n",
        "    print(\"Dataset, train set, and test set size:\", len(y_train) + len(y_valid) + len(y_test), len(y_train), len(y_test))\n",
        "    print('Timeframe:', time_frame)\n",
        "\n",
        "    #print('Predicted vs. Actual percent difference: %.3f%%' % ((total_E * 100 / total_E_pred) - 100))\n",
        "    print('Maximum possible activations:', max_act)\n",
        "    print('Predicted activations:', pred_act)\n",
        "    print('Successful activations: %d, %.3f%%' % (succ_act, succ_act * 100/pred_act))\n",
        "    print('Failed activations: %d, %.3f%%' % (false_act, false_act * 100/pred_act))\n",
        "    print('Missed activations: %d, %.3f%%' % (max_act - succ_act, (max_act - succ_act) * 100/max_act))\n",
        "\n",
        "\n",
        "\n",
        "    print('Voltage overestimation rate: %.3f%%' % ((y_test['Voltage (mV)'].values <= y_test['Voltage (mV) pred']).mean() * 100))\n",
        "    print(\"Test MAPE power: %3f\" %  MAPE(y_test['Power (uW)'].values.ravel(), y_test['power pred']))\n",
        "    print(\"Test MAPE voltage: %3f\" % MAPE(y_test['Voltage (mV)'], y_test['Voltage (mV) pred']))\n",
        "    print(\"Test MAPE current: %3f\" % MAPE(y_test['Current (uA)'], y_test['Current (uA) pred']))\n",
        "    print('Successful activations compared to max: %d, %.3f%%' % (succ_act, succ_act * 100/max_act))\n",
        "\n",
        "    #Naive model\n",
        "    print('Naive predicted activations (usual actual energy average):', pred_act_naive)\n",
        "    print('Naive successful activations (usual actual energy average): %d, %.3f%%' % (succ_act_naive, succ_act_naive * 100/pred_act_naive))\n",
        "    print('Naive failed activations (usual actual energy average): %d, %.3f%%' % (false_act_naive, false_act_naive * 100/pred_act_naive))\n",
        "    print('Naive missed activations (usual actual energy average): %d, %.3f%%' % (max_act - succ_act_naive, (max_act - succ_act_naive) * 100/max_act))\n",
        "\n",
        "    print('Naive successful activations (compared to max): %d, %.3f%%' % (succ_act_naive, succ_act_naive * 100/max_act))\n",
        "\n",
        "    E_pred = 0\n",
        "    for i in range(len(y_test) - 1):\n",
        "      t = (y_test.index[i+1] - y_test.index[i]).total_seconds()\n",
        "      if t <= time_frame_seconds + 50:\n",
        "        E_pred += y_test['power pred'][i] * t\n",
        "\n",
        "    print('Predicted vs. Actual Total Energy Percent Difference: %.3f%%' % ((E_pred - E_actual) * 100 / E_actual))\n",
        "    print(\"Test MAPE power:\\n\", MAPE(y_test['Power (uW)'].values.ravel(), y_test['power pred'].values.ravel()))\n",
        "\n",
        "    V_actual = y_test['Voltage (mV)'].mean()\n",
        "    V_pred = y_test['Voltage (mV) pred'].mean()\n",
        "    print('Predicted vs. Actual Total Voltage Percent Difference: %.3f%%' % ((V_pred - V_actual) * 100 / V_actual))\n",
        "    print(\"Test MAPE voltage:\\n\", MAPE(y_test['Voltage (mV)'].values.ravel(), y_test['Voltage (mV) pred'].values.ravel()))\n",
        "\n",
        "\n",
        "\n",
        "    print('Voltage overestimation rate: %.3f%%' % ((y_test['Voltage (mV)'].values <= y_test['Voltage (mV) pred']).mean() * 100))\n",
        "    print(\"Test MAPE power: %3f\" %  MAPE(y_test['Power (uW)'].values.ravel(), y_test['power pred']))\n",
        "    print(\"Test MAPE voltage: %3f\" % MAPE(y_test['Voltage (mV)'], y_test['Voltage (mV) pred']))\n",
        "    print(\"Test MAPE current: %3f\" % MAPE(y_test['Current (uA)'], y_test['Current (uA) pred']))\n",
        "    print(v_test)"
      ],
      "metadata": {
        "id": "5k2ZgWSWAPp6",
        "collapsed": true
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Graph Selected Data"
      ],
      "metadata": {
        "id": "oTdyJAYj40w3"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#Prep Data for Graphing\n",
        "from keras.models import load_model\n",
        "\n",
        "time_frame = '30min'\n",
        "\n",
        "X_test = pd.concat([df_test[\"power - 1h\"], df_test[\"power - 2h\"], df_test[\"power - 3h\"], df_test[\"V1 - 1h\"], df_test[\"V1 - 2h\"], df_test[\"V1 - 3h\"], df_test[\"I1L - 1h\"], df_test[\"I1L - 2h\"], df_test[\"I1L - 3h\"],df_test[\"EC - 1h\"], df_test[\"EC - 2h\"], df_test[\"EC - 3h\"], df_test[\"raw_VWC - 1h\"], df_test[\"raw_VWC - 2h\"], df_test[\"raw_VWC - 3h\"], df_test[\"temp - 1h\"], df_test[\"temp - 2h\"], df_test[\"temp - 3h\"], df_test[\"tsd\"], df_test[\"hour\"]], axis = 1)\n",
        "#X_test = pd.concat([df_test[\"EC - 1h\"], df_test[\"EC - 2h\"], df_test[\"EC - 3h\"], df_test[\"raw_VWC - 1h\"], df_test[\"raw_VWC - 2h\"], df_test[\"raw_VWC - 3h\"], df_test[\"temp - 1h\"], df_test[\"temp - 2h\"], df_test[\"temp - 3h\"], df_test[\"tsd\"], df_test[\"hour\"]], axis = 1)\n",
        "#X_test = pd.concat([df_test[\"power - 1h\"], df_test[\"power - 2h\"], df_test[\"power - 3h\"], df_test[\"V1 - 1h\"], df_test[\"V1 - 2h\"], df_test[\"V1 - 3h\"], df_test[\"I1L - 1h\"], df_test[\"I1L - 2h\"], df_test[\"I1L - 3h\"]], axis = 1)\n",
        "#X_new = pd.concat([df_test[\"V1 - 1h\"], df_test[\"V1 - 2h\"]], axis = 1)\n",
        "y_test = pd.concat([df_test['Power (uW)'], df_test['Voltage (mV)'], df_test['Current (uA)']], axis = 1)\n",
        "\n",
        "#Resample data\n",
        "X_test = X_test.resample(time_frame).mean().dropna()\n",
        "\n",
        "y_test = y_test.resample(time_frame).mean().dropna()\n",
        "\n",
        "#Reshape data\n",
        "X_test = X_test.values.reshape((X_test.shape[0], 1, X_test.shape[1]))\n",
        "\n",
        "mv1 = y_test\n",
        "\n",
        "def quantile_loss(y_true, y_pred, quantile = 0.5):\n",
        "    error = y_true - y_pred\n",
        "    return tf.reduce_mean(K.maximum(quantile * error, (quantile - 1) * error), axis=-1)\n",
        "\n",
        "file = 'trained_models/lstm16_' + time_frame + '_quant50'\n",
        "model = load_model(file, custom_objects={'quantile_loss': quantile_loss})\n",
        "predictions = model.predict(X_test)\n",
        "\n",
        "mv1[\"power_pred_med\"] = predictions[:, 0]\n",
        "mv1[\"voltage_pred_med\"] = predictions[:, 1]\n",
        "mv1[\"current_pred_med\"] = predictions[:, 2]\n",
        "\n",
        "def quantile_loss(y_true, y_pred, quantile = 0.05):\n",
        "    error = y_true - y_pred\n",
        "    return tf.reduce_mean(K.maximum(quantile * error, (quantile - 1) * error), axis=-1)\n",
        "\n",
        "file = 'trained_models/lstm16_' + time_frame + '_quant5'\n",
        "model = load_model(file, custom_objects={'quantile_loss': quantile_loss})\n",
        "predictions = model.predict(X_test)\n",
        "\n",
        "mv1[\"power_pred_lower\"] = predictions[:, 0]\n",
        "mv1[\"voltage_pred_lower\"] = predictions[:, 1]\n",
        "mv1[\"current_pred_lower\"] = predictions[:, 2]\n",
        "\n",
        "def quantile_loss(y_true, y_pred, quantile = 0.95):\n",
        "    error = y_true - y_pred\n",
        "    return tf.reduce_mean(K.maximum(quantile * error, (quantile - 1) * error), axis=-1)\n",
        "\n",
        "file = 'trained_models/lstm16_' + time_frame + '_quant95'\n",
        "model = load_model(file, custom_objects={'quantile_loss': quantile_loss})\n",
        "predictions = model.predict(X_test)\n",
        "mv1[\"power_pred_upper\"] = predictions[:, 0]\n",
        "mv1[\"voltage_pred_upper\"] = predictions[:, 1]\n",
        "mv1[\"current_pred_upper\"] = predictions[:, 2]\n",
        "\n",
        "\n",
        "mv1 = mv1.loc[(mv1.index > '2024-04-24 12:00:00') & (mv1.index < '2024-05-02')]\n",
        "mv2 = mv1"
      ],
      "metadata": {
        "id": "FG1eBXH-43M7"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "import matplotlib as mpl\n",
        "from matplotlib.ticker import FixedLocator\n",
        "mpl.use('Agg')\n",
        "#mpl.rc('font', **font)\n",
        "import matplotlib.pyplot as plt\n",
        "import matplotlib.gridspec as gridspec\n",
        "import matplotlib.dates as md\n",
        "import datetime\n",
        "import numpy as np\n",
        "from pytz import timezone\n",
        "import pandas as pd\n",
        "import arrow\n",
        "import os\n",
        "from pandas.plotting import register_matplotlib_converters\n",
        "register_matplotlib_converters()\n",
        "%matplotlib inline\n",
        "\n",
        "# Limits for graphs\n",
        "VOLTAGE_LIM = 400\n",
        "CURRENT_LIM = 150\n",
        "POWER_LIM = 40\n",
        "\n",
        "line_width = 0.5\n",
        "plt.title(\"1 hour prediction interval\")\n",
        "\n",
        "plt.close()\n",
        "plt.xlabel(\"Time\")\n",
        "fig, (ax1, ax2, ax3) = plt.subplots(3,figsize=(4,3), sharex=False)\n",
        "fig.autofmt_xdate()\n",
        "\n",
        "volt_color= 'tab:blue'\n",
        "\n",
        "amp_color = 'tab:red'\n",
        "\n",
        "\n",
        "volt_color1= 'tab:blue'\n",
        "volt_style1 = 'dashed'\n",
        "volt_color2= 'tab:green'\n",
        "volt_style2 = 'dashed'\n",
        "volt_color3= 'tab:red'\n",
        "volt_style3 = 'dashed'\n",
        "volt_color4= 'tab:orange'\n",
        "volt_style4 = 'dashed'\n",
        "\n",
        "#amp_color1 = 'tab:red'\n",
        "#amp_style1='dashed'\n",
        "#amp_color2 = 'tab:orange'\n",
        "#amp_style2='dashdot'\n",
        "\n",
        "ax1.fmt_xdata = md.DateFormatter('%m-%d')\n",
        "ax1.xaxis.set_major_formatter(md.DateFormatter('%m-%d'))\n",
        "ax1.grid(True)\n",
        "\n",
        "ax2.fmt_xdata = md.DateFormatter('%m-%d')\n",
        "ax2.xaxis.set_major_formatter(md.DateFormatter('%m-%d'))\n",
        "ax2.grid(True)\n",
        "\n",
        "ax3.fmt_xdata = md.DateFormatter('%m-%d')\n",
        "ax3.xaxis.set_major_formatter(md.DateFormatter('%m-%d'))\n",
        "ax3.grid(True)\n",
        "\n",
        "ax1.set_ylabel(\"Power (μW)\", fontsize=7.5, labelpad = 10)\n",
        "ax2.set_ylabel(\"Voltage (mV)\", fontsize=7.5, labelpad = 7)\n",
        "ax3.set_ylabel(\"Current (μA)\", fontsize=7.5, labelpad = 7)\n",
        "#ax3.set_xlabel(\"Date\", fontsize=10)\n",
        "\n",
        "ax1.set_ylim(-10, 60)\n",
        "ax2.set_ylim(-25, 450)\n",
        "ax3.set_ylim(-10, 200)\n",
        "\n",
        "#ax3.set_xlim([datetime.date(2021, 12, 12), datetime.date(2021, 12, 14)])\n",
        "ax1.set_xlim([mv1.index[0], mv1.index[-1]])\n",
        "ax2.set_xlim([mv1.index[0], mv1.index[-1]])\n",
        "ax3.set_xlim([mv1.index[0], mv1.index[-1]])\n",
        "\n",
        "ax1.plot(mv1.index, mv2[\"power_pred_lower\"], color='tab:red', ls = volt_style3, linewidth=line_width)\n",
        "ax1.plot(mv1.index, mv2[\"power_pred_upper\"], color='tab:orange', ls = volt_style4, linewidth=line_width)\n",
        "ax1.plot(mv1.index, mv2[\"Power (uW)\"], color='tab:blue', ls = 'solid', linewidth=line_width)\n",
        "ax1.plot(mv1.index, mv2[\"power_pred_med\"], color='tab:green', ls = volt_style2, linewidth=line_width)\n",
        "ax1.fill_between(mv1.index, mv2[\"power_pred_lower\"], mv2[\"power_pred_upper\"], color='grey', alpha=0.5)\n",
        "#ax1.legend(['lower bound', 'upper bound', 'actual', 'median predictions'], loc='lower center', prop={'size': 5.5}, ncol=2)\n",
        "\n",
        "ax2.plot(mv1.index, mv2[\"voltage_pred_lower\"], color='tab:red', ls = volt_style3, linewidth=line_width)\n",
        "ax2.plot(mv1.index, mv2[\"voltage_pred_upper\"], color='tab:orange', ls = volt_style4, linewidth=line_width)\n",
        "ax2.plot(mv1.index, mv2['Voltage (mV)'], color='tab:blue', ls = 'solid', linewidth=line_width)\n",
        "ax2.plot(mv1.index, mv2[\"voltage_pred_med\"], color='tab:green', ls = volt_style2, linewidth=line_width)\n",
        "ax2.fill_between(mv1.index, mv2[\"voltage_pred_lower\"], mv2[\"voltage_pred_upper\"], color='grey', alpha=0.5)\n",
        "ax2.legend(['lower bound', 'upper bound', 'ground truth', 'median prediction'], loc='lower left', prop={'size': 6.6}, ncol=2, columnspacing=0.5)\n",
        "\n",
        "ax3.plot(mv1.index, mv2[\"current_pred_lower\"], color='tab:red', ls = volt_style3, linewidth=line_width)\n",
        "ax3.plot(mv1.index, mv2[\"current_pred_upper\"], color='tab:orange', ls = volt_style4, linewidth=line_width)\n",
        "ax3.plot(mv1.index, mv2['Current (uA)'], color='tab:blue', ls = 'solid', linewidth=line_width)\n",
        "ax3.plot(mv1.index, mv2[\"current_pred_med\"], color='tab:green', ls = volt_style2, linewidth=line_width)\n",
        "ax3.fill_between(mv1.index, mv2[\"current_pred_lower\"], mv2[\"current_pred_upper\"], color='grey', alpha=0.5)\n",
        "#ax3.legend(['lower bound', 'upper bound', 'actual', 'median predictions'], loc='upper right', prop={'size': 5.5}, ncol=1)\n",
        "\n",
        "\n",
        "#Plot error\n",
        "#ax3.plot(mv1['timestamp'], mv1['error']/mv1['power'], color=volt_color2, ls = volt_style2)\n",
        "#ax3.legend(['error'], loc='upper right', prop={'size': 6})\n",
        "\n",
        "ax3.tick_params(axis='x', labelsize=7.5, rotation=30, pad = 0.1,)\n",
        "ax3.set_xticks(list(ax3.get_xticks()))\n",
        "for label in ax3.get_xticklabels():\n",
        "    label.set_horizontalalignment('center')\n",
        "\n",
        "ax1.tick_params(axis='y', labelsize=7.5, rotation=0, pad = 0.1)\n",
        "ax2.tick_params(axis='y', labelsize=7.5, rotation=0, pad = 0.1)\n",
        "ax3.tick_params(axis='y', labelsize=7.5, rotation=0, pad = 0.1)\n",
        "\n",
        "plt.tight_layout(pad=0.3, w_pad=0.5, h_pad=0.18)\n",
        "#plt.subplots_adjust(hspace=0.15)\n",
        "plt.savefig('twobat.pdf')\n",
        "plt.close()"
      ],
      "metadata": {
        "id": "pfPA-38C44AE"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "###  2.3 Load and graph pretrained SNN models\n",
        "\n",
        "###  In order to use pretrained models it is neccesary to download the ```trained_models``` directory from [Hugging Face](https://huggingface.co/datasets/adunlop621/Soil_MFC/tree/main) and store it in the same directory as this notebook"
      ],
      "metadata": {
        "id": "V7juu2pL0QXm"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from keras.models import load_model\n",
        "\n",
        "time_frame = '60min'\n",
        "batchsize = 8\n",
        "\n",
        "X = pd.concat([df[\"power - 1h\"], df[\"power - 2h\"], df[\"power - 3h\"], df[\"V1 - 1h\"], df[\"V1 - 2h\"], df[\"V1 - 3h\"], df[\"I1L - 1h\"], df[\"I1L - 2h\"], df[\"I1L - 3h\"],df[\"EC - 1h\"], df[\"EC - 2h\"], df[\"EC - 3h\"], df[\"raw_VWC - 1h\"], df[\"raw_VWC - 2h\"], df[\"raw_VWC - 3h\"], df[\"temp - 1h\"], df[\"temp - 2h\"], df[\"temp - 3h\"], df[\"tsd\"], df[\"hour\"]], axis = 1)\n",
        "#X = pd.concat([df[\"power - 1h\"], df[\"power - 2h\"], df[\"power - 3h\"], df[\"V1 - 1h\"], df[\"V1 - 2h\"], df[\"V1 - 3h\"], df[\"I1L - 1h\"], df[\"I1L - 2h\"], df[\"I1L - 3h\"], df[\"tsd\"], df[\"hour\"]], axis = 1)\n",
        "y = pd.concat([df[\"Power (uW)\"], df['Voltage (mV)'], df['Current (uA)']], axis = 1)\n",
        "\n",
        "#Normalize Data\n",
        "X_normalized = ((X - X.min()) / (X.max() - X.min()))\n",
        "\n",
        "#Split train and test sets\n",
        "X_train, X_test = train_test_split(X_normalized, test_size=0.3, shuffle=False)\n",
        "y_train, y_test = train_test_split(y, test_size=0.3, shuffle=False)\n",
        "\n",
        "X_valid, X_test = train_test_split(X_test, test_size=0.5, shuffle=False)\n",
        "y_valid, y_test = train_test_split(y_test, test_size=0.5, shuffle=False)\n",
        "\n",
        "#Resample data\n",
        "\n",
        "X_valid = X_valid.resample(time_frame).mean().dropna()\n",
        "y_valid = y_valid.resample(time_frame).mean().dropna()\n",
        "\n",
        "X_test = X_test.resample(time_frame).mean().dropna()\n",
        "y_test = y_test.resample(time_frame).mean().dropna()\n",
        "\n",
        "#Define mv1\n",
        "mv1 = y_valid\n",
        "\n",
        "#Reshape data\n",
        "X_train = X_train.values.reshape((X_train.shape[0], 1, X_train.shape[1]))\n",
        "X_valid = X_valid.values.reshape((X_valid.shape[0], 1, X_valid.shape[1]))\n",
        "X_test = X_test.values.reshape((X_test.shape[0], 1, X_test.shape[1]))\n",
        "\n",
        "# convert to tensor\n",
        "X_train = torch.tensor(X_train)\n",
        "y_train = torch.tensor(y_train.values)\n",
        "X_valid = torch.tensor(X_valid)\n",
        "y_valid = torch.tensor(y_valid.values)\n",
        "X_test = torch.tensor(X_test)\n",
        "y_test = torch.tensor(y_test.values)\n",
        "\n",
        "# make datasets\n",
        "train_dataset = TensorDataset(X_train, y_train)\n",
        "valid_dataset = TensorDataset(X_valid, y_valid)\n",
        "test_dataset = TensorDataset(X_test, y_test)\n",
        "\n",
        "# Create DataLoaders\n",
        "train_loader = DataLoader(train_dataset, batch_size=batchsize, shuffle=False)\n",
        "valid_loader = DataLoader(valid_dataset, batch_size=batchsize, shuffle=False)\n",
        "test_loader = DataLoader(test_dataset, batch_size=batchsize, shuffle=False)\n",
        "\n",
        "num_steps = 50\n",
        "num_inputs = X_train.shape[2]\n",
        "\n",
        "# create new inctance of the SNN Class\n",
        "model = Net(num_inputs, num_steps).to(device)\n",
        "\n",
        "model.load_state_dict(torch.load(\"trained_models/snn_60min_quant50\", map_location=torch.device('cpu')))\n",
        "model.eval()\n",
        "actuals = []\n",
        "predictions = []\n",
        "\n",
        "with torch.no_grad():\n",
        "    for data, targets in valid_loader:\n",
        "\n",
        "        # prepare data\n",
        "        data = data.to(device)\n",
        "        targets = targets.to(device)\n",
        "\n",
        "        data = data.float()\n",
        "        targets = targets.float()\n",
        "\n",
        "        _, _, _, output = model(data)\n",
        "\n",
        "        output = output.cpu()\n",
        "        output = output.squeeze(1).detach()\n",
        "\n",
        "        prediction = output[-1]\n",
        "\n",
        "        actuals.append(targets)\n",
        "        predictions.append(prediction)\n",
        "\n",
        "# Convert lists to tensors\n",
        "actuals = torch.cat(actuals, dim=0)\n",
        "predictions = torch.cat(predictions, dim=0)\n",
        "\n",
        "mv1[\"power_pred_med\"] = predictions[:, 0]\n",
        "mv1[\"voltage_pred_med\"] = predictions[:, 1]\n",
        "mv1[\"current_pred_med\"] = predictions[:, 2]\n",
        "#mv1 = mv1.loc[(mv1.index > '2022-01-04') & (mv1.index < '2022-01-06')]\n",
        "mv1 = mv1.loc[(mv1.index > '2021-12-12') & (mv1.index < '2021-12-14')]\n",
        "mv2 = mv1"
      ],
      "metadata": {
        "id": "9_jKnA030V9a"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Other Performance Metrics for Pretrained SNNs"
      ],
      "metadata": {
        "id": "jEBBNiE_03PT"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from keras.models import load_model\n",
        "\n",
        "#Set parameters\n",
        "batchsize_list = [300, 150, 50, 20, 8]\n",
        "time_frame_list = ['3min', '5min', '15min', '30min', '60min']\n",
        "time_frame_seconds_list = [180, 300, 900, 1800, 3600]\n",
        "n = 0\n",
        "\n",
        "snn_power_mape_list = []\n",
        "snn_volt_mape_list = []\n",
        "snn_curr_mape_list = []\n",
        "\n",
        "for j in range(len(batchsize_list)):\n",
        "    batchsize = batchsize_list[j]\n",
        "    time_frame = time_frame_list[j]\n",
        "    time_frame_seconds = time_frame_seconds_list[j]\n",
        "\n",
        "    X = pd.concat([df[\"power - 1h\"], df[\"power - 2h\"], df[\"power - 3h\"], df[\"V1 - 1h\"], df[\"V1 - 2h\"], df[\"V1 - 3h\"], df[\"I1L - 1h\"], df[\"I1L - 2h\"], df[\"I1L - 3h\"],df[\"EC - 1h\"], df[\"EC - 2h\"], df[\"EC - 3h\"], df[\"raw_VWC - 1h\"], df[\"raw_VWC - 2h\"], df[\"raw_VWC - 3h\"], df[\"temp - 1h\"], df[\"temp - 2h\"], df[\"temp - 3h\"], df[\"tsd\"], df[\"hour\"]], axis = 1)\n",
        "    #X = pd.concat([df[\"power - 1h\"], df[\"power - 2h\"], df[\"power - 3h\"], df[\"V1 - 1h\"], df[\"V1 - 2h\"], df[\"V1 - 3h\"], df[\"I1L - 1h\"], df[\"I1L - 2h\"], df[\"I1L - 3h\"], df[\"tsd\"], df[\"hour\"]], axis = 1)\n",
        "    y = pd.concat([df[\"Power (uW)\"], df['Voltage (mV)'], df['Current (uA)']], axis = 1)\n",
        "\n",
        "    #Normalize Data\n",
        "    X_normalized = ((X - X.min()) / (X.max() - X.min()))\n",
        "\n",
        "    #Split train and test sets\n",
        "    X_train, X_test = train_test_split(X_normalized, test_size=0.3, shuffle=False)\n",
        "    y_train, y_test = train_test_split(y, test_size=0.3, shuffle=False)\n",
        "\n",
        "    X_valid, X_test = train_test_split(X_test, test_size=0.5, shuffle=False)\n",
        "    y_valid, y_test = train_test_split(y_test, test_size=0.5, shuffle=False)\n",
        "\n",
        "    #Calculate actual energy generated in test set\n",
        "    E_actual = 0\n",
        "    for i in range(len(y_test) - 1):\n",
        "      t = (y_test.index[i+1] - y_test.index[i]).total_seconds()\n",
        "      if t < 180:\n",
        "        E_actual += y_test['Power (uW)'][i] * t\n",
        "\n",
        "    #Resample data\n",
        "\n",
        "    X_valid = X_valid.resample(time_frame).mean().dropna()\n",
        "    y_valid = y_valid.resample(time_frame).mean().dropna()\n",
        "\n",
        "    X_test = X_test.resample(time_frame).mean().dropna()\n",
        "    y_test = y_test.resample(time_frame).mean().dropna()\n",
        "\n",
        "    #Define mv1\n",
        "    mv1 = y_test\n",
        "\n",
        "    #Reshape data\n",
        "    X_train = X_train.values.reshape((X_train.shape[0], 1, X_train.shape[1]))\n",
        "    X_valid = X_valid.values.reshape((X_valid.shape[0], 1, X_valid.shape[1]))\n",
        "    X_test = X_test.values.reshape((X_test.shape[0], 1, X_test.shape[1]))\n",
        "\n",
        "    # convert to tensor\n",
        "    X_train = torch.tensor(X_train)\n",
        "    y_train = torch.tensor(y_train.values)\n",
        "    X_valid = torch.tensor(X_valid)\n",
        "    y_valid = torch.tensor(y_valid.values)\n",
        "    X_test = torch.tensor(X_test)\n",
        "    y_test = torch.tensor(y_test.values)\n",
        "\n",
        "    # make datasets\n",
        "    train_dataset = TensorDataset(X_train, y_train)\n",
        "    valid_dataset = TensorDataset(X_valid, y_valid)\n",
        "    test_dataset = TensorDataset(X_test, y_test)\n",
        "\n",
        "    # Create DataLoaders\n",
        "    train_loader = DataLoader(train_dataset, batch_size=batchsize, shuffle=False)\n",
        "    valid_loader = DataLoader(valid_dataset, batch_size=batchsize, shuffle=False)\n",
        "    test_loader = DataLoader(test_dataset, batch_size=batchsize, shuffle=False)\n",
        "\n",
        "    num_steps = 50\n",
        "    num_inputs = X_train.shape[2]\n",
        "\n",
        "    # create new inctance of the SNN Class\n",
        "    model = Net(num_inputs, num_steps).to(device)\n",
        "\n",
        "    file = 'trained_models/snn_' + time_frame_list[j] + '_quant50.pth'\n",
        "    print(file)\n",
        "\n",
        "    checkpoint = torch.load(file, map_location=torch.device('cpu'))\n",
        "    model.load_state_dict(checkpoint['state_dict'])\n",
        "\n",
        "    #model.load_state_dict(torch.load(file, map_location=torch.device('cpu'))['model_state_dict'])\n",
        "    model.eval()\n",
        "    actuals = []\n",
        "    predictions = []\n",
        "\n",
        "    with torch.no_grad():\n",
        "        for data, targets in test_loader:\n",
        "\n",
        "            # prepare data\n",
        "            data = data.to(device)\n",
        "            targets = targets.to(device)\n",
        "\n",
        "            data = data.float()\n",
        "            targets = targets.float()\n",
        "\n",
        "            _, _, _, output = model(data)\n",
        "\n",
        "            output = output.cpu()\n",
        "            output = output.squeeze(1).detach()\n",
        "\n",
        "            prediction = output[-1]\n",
        "\n",
        "            actuals.append(targets)\n",
        "            predictions.append(prediction)\n",
        "\n",
        "    # Convert lists to tensors\n",
        "    actuals = torch.cat(actuals, dim=0)\n",
        "    predictions = torch.cat(predictions, dim=0)\n",
        "\n",
        "    mv1[\"power_pred_med_\" + time_frame_list[j]] = predictions[:, 0].numpy()\n",
        "    mv1[\"voltage_pred_med_\" + time_frame_list[j]] = predictions[:, 1].numpy()\n",
        "    mv1[\"current_pred_med_\" + time_frame_list[j]] = predictions[:, 2].numpy()\n",
        "\n",
        "    print('Voltage overestimation rate: %.3f%%' % ((mv1['Voltage (mV)'].values <= mv1[\"voltage_pred_med_\" + time_frame_list[j]]).mean() * 100))\n",
        "    print(\"Test MAPE power: %3f\" %  MAPE(mv1['Power (uW)'].values.ravel(), mv1[\"power_pred_med_\" + time_frame_list[j]]))\n",
        "    print(\"Test MAPE voltage: %3f\" % MAPE(mv1['Voltage (mV)'], mv1[\"voltage_pred_med_\" + time_frame_list[j]]))\n",
        "    print(\"Test MAPE current: %3f\" % MAPE(mv1['Current (uA)'], mv1[\"current_pred_med_\" + time_frame_list[j]]))\n",
        "\n",
        "    E_pred = 0\n",
        "    for i in range(len(mv1) - 1):\n",
        "      t = (mv1.index[i+1] - mv1.index[i]).total_seconds()\n",
        "      if t <= time_frame_seconds + 50:\n",
        "        E_pred += mv1[\"power_pred_med_\" + time_frame_list[j]][i] * t\n",
        "\n",
        "    print('Predicted vs. Actual Total Energy Percent Difference: %.3f%%' % ((E_pred - E_actual) * 100 / E_actual))\n",
        "\n",
        "    V_actual = mv1['Voltage (mV)'].mean()\n",
        "    V_pred = mv1[\"voltage_pred_med_\" + time_frame_list[j]].mean()\n",
        "    print('Predicted vs. Actual Total Voltage Percent Difference: %.3f%%' % ((V_pred - V_actual) * 100 / V_actual))\n",
        "    #print(mv1)\n",
        "    #mv1 = mv1.loc[(mv1.index > '2022-01-04') & (mv1.index < '2022-01-06')]\n",
        "    #mv1 = mv1.loc[(mv1.index > '2021-12-12') & (mv1.index < '2021-12-14')]\n",
        "    #mv2 = mv1"
      ],
      "metadata": {
        "id": "_xTlvH2q04S6"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### Get total spikes for pretrained SNN models"
      ],
      "metadata": {
        "id": "XJ4aJz3V1JDP"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from keras.models import load_model\n",
        "import time\n",
        "from sklearn.model_selection import TimeSeriesSplit\n",
        "from keras.models import Sequential\n",
        "from keras.layers import Dense\n",
        "from keras.layers import LSTM\n",
        "from keras import backend as K\n",
        "from torch.utils.data import DataLoader, TensorDataset\n",
        "\n",
        "#Set parameters\n",
        "batchsize_list = [300, 150, 50, 20, 8]\n",
        "time_frame_list = ['3min', '5min', '15min', '30min', '60min']\n",
        "time_frame_seconds_list = [180, 300, 900, 1800, 3600]\n",
        "n = 0\n",
        "\n",
        "power_mape_list = []\n",
        "volt_mape_list = []\n",
        "curr_mape_list = []\n",
        "total_spk_snn_list = []\n",
        "\n",
        "for j in range(len(batchsize_list)):\n",
        "  if j>=0:\n",
        "    batchsize = batchsize_list[j]\n",
        "    time_frame = time_frame_list[j]\n",
        "    time_frame_seconds = time_frame_seconds_list[j]\n",
        "\n",
        "    X = pd.concat([df[\"power - 1h\"], df[\"power - 2h\"], df[\"power - 3h\"], df[\"V1 - 1h\"], df[\"V1 - 2h\"], df[\"V1 - 3h\"], df[\"I1L - 1h\"], df[\"I1L - 2h\"], df[\"I1L - 3h\"],df[\"EC - 1h\"], df[\"EC - 2h\"], df[\"EC - 3h\"], df[\"raw_VWC - 1h\"], df[\"raw_VWC - 2h\"], df[\"raw_VWC - 3h\"], df[\"temp - 1h\"], df[\"temp - 2h\"], df[\"temp - 3h\"], df[\"tsd\"], df[\"hour\"]], axis = 1)\n",
        "    #X = pd.concat([df[\"power - 1h\"], df[\"power - 2h\"], df[\"power - 3h\"], df[\"V1 - 1h\"], df[\"V1 - 2h\"], df[\"V1 - 3h\"], df[\"I1L - 1h\"], df[\"I1L - 2h\"], df[\"I1L - 3h\"], df[\"tsd\"], df[\"hour\"]], axis = 1)\n",
        "    y = pd.concat([df[\"Power (uW)\"], df['Voltage (mV)'], df['Current (uA)']], axis = 1)\n",
        "\n",
        "    #Normalize Data\n",
        "    X_normalized = ((X - X.min()) / (X.max() - X.min()))\n",
        "\n",
        "    #Split train and test sets\n",
        "    X_train, X_test = train_test_split(X_normalized, test_size=0.3, shuffle=False)\n",
        "    y_train, y_test = train_test_split(y, test_size=0.3, shuffle=False)\n",
        "\n",
        "    X_valid, X_test = train_test_split(X_test, test_size=0.5, shuffle=False)\n",
        "    y_valid, y_test = train_test_split(y_test, test_size=0.5, shuffle=False)\n",
        "\n",
        "    #Resample data\n",
        "    X_valid = X_valid.resample(time_frame).mean().dropna()\n",
        "    y_valid = y_valid.resample(time_frame).mean().dropna()\n",
        "\n",
        "    X_test = X_test.resample(time_frame).mean().dropna()\n",
        "    y_test = y_test.resample(time_frame).mean().dropna()\n",
        "\n",
        "    #Reshape data\n",
        "    X_train = X_train.values.reshape((X_train.shape[0], 1, X_train.shape[1]))\n",
        "    X_valid = X_valid.values.reshape((X_valid.shape[0], 1, X_valid.shape[1]))\n",
        "    X_test = X_test.values.reshape((X_test.shape[0], 1, X_test.shape[1]))\n",
        "\n",
        "    # convert to tensor\n",
        "    X_train = torch.tensor(X_train)\n",
        "    y_train = torch.tensor(y_train.values)\n",
        "    X_valid = torch.tensor(X_valid)\n",
        "    y_valid = torch.tensor(y_valid.values)\n",
        "    X_test = torch.tensor(X_test)\n",
        "    y_test = torch.tensor(y_test.values)\n",
        "\n",
        "    # make datasets\n",
        "    train_dataset = TensorDataset(X_train, y_train)\n",
        "    valid_dataset = TensorDataset(X_valid, y_valid)\n",
        "    test_dataset = TensorDataset(X_test, y_test)\n",
        "\n",
        "    # Create DataLoaders\n",
        "    train_loader = DataLoader(train_dataset, batch_size=batchsize, shuffle=False)\n",
        "    valid_loader = DataLoader(valid_dataset, batch_size=batchsize, shuffle=False)\n",
        "    test_loader = DataLoader(test_dataset, batch_size=batchsize, shuffle=False)\n",
        "\n",
        "    num_steps = 50\n",
        "    num_inputs = X_train.shape[2]\n",
        "\n",
        "    # create new inctance of the SNN Class\n",
        "    model = Net(num_inputs, num_steps).to(device)\n",
        "\n",
        "    file = 'trained_models/snn_' + time_frame_list[j] + '_quant50.pth'\n",
        "\n",
        "    checkpoint = torch.load(file, map_location=torch.device('cpu'))\n",
        "    model.load_state_dict(checkpoint['state_dict'])\n",
        "\n",
        "\n",
        "    #model.load_state_dict(torch.load(file, map_location=torch.device('cpu'))['model_state_dict'])\n",
        "    model.eval()\n",
        "    actuals = []\n",
        "    predictions = []\n",
        "\n",
        "    total_spk_snn = 0\n",
        "    n = 0\n",
        "    with torch.no_grad():\n",
        "        for data, targets in test_loader:\n",
        "            n += 1\n",
        "            if n==1:\n",
        "              # prepare data\n",
        "              data = data.to(device)\n",
        "              targets = targets.to(device)\n",
        "\n",
        "              data = data.float()\n",
        "              targets = targets.float()\n",
        "\n",
        "              data = data[0]\n",
        "\n",
        "              spk1, spk2, spk3, output = model(data)\n",
        "\n",
        "              total_spk_snn = (spk1.sum() + spk2.sum() + spk3.sum() + output.sum()).item() # count up the total number of spikes\n",
        "\n",
        "    total_spk_snn_list.append(total_spk_snn)"
      ],
      "metadata": {
        "id": "wGu7e_w_1Mhh"
      },
      "execution_count": null,
      "outputs": []
    }
  ],
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}