step2_anomaly_screening.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import os\n",
    "import pickle\n",
    "from datetime import datetime"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define Which Input Files to Use\n",
    "The default settings will use the input files recently produced in Step 1) using the notebook `get_eia_demand_data.ipynb`. For those interested in reproducing the exact results included in the repository, you will need to point to the files containing the original `raw` EIA demand data that we querried on 10 Sept 2019.\n",
    "\n",
    "Use the `clean_subregions` flag to clean and prepare the BAs and subregions for imputation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "check_reproducibility = os.getenv('CHECK_REPRODUCIBILITY', default=None)\n",
    "\n",
    "clean_subregions = False\n",
    "#clean_subregions = True\n",
    "\n",
    "\n",
    "\n",
    "if not os.path.exists('./data'):\n",
    "    os.mkdir('./data')\n",
    "\n",
    "if check_reproducibility == 'TRUE':\n",
    "    # download published files from Zenodo for comparison if not already present\n",
    "    # v1.1 was used for Scientific Data paper: https://zenodo.org/record/3690240\n",
    "    if not os.path.exists('./data/truggles-EIA_Cleaned_Hourly_Electricity_Demand_Data-1b31ad5'):\n",
    "        !curl -o ./data/eia_data.zip https://zenodo.org/record/3690240/files/truggles/EIA_Cleaned_Hourly_Electricity_Demand_Data-v1.1.zip?download=1\n",
    "        !unzip ./data/eia_data.zip -d ./data\n",
    "    input_path = './data/truggles-EIA_Cleaned_Hourly_Electricity_Demand_Data-1b31ad5/data/release_2019_Oct/original_eia_files'\n",
    "    assert(os.path.exists(input_path)), f\"Were the files accessed from Zenodo and unzipped correctly?\"\n",
    "\n",
    "# Running with step1 files\n",
    "else:\n",
    "    input_path = './data'\n",
    "    if clean_subregions:\n",
    "        input_path = './data_subregions'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define The Screening Algorithms\n",
    "\n",
    "The following cell defines all of the screening algorithms used. The specific selection parameters are defined in a following cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_rolling_dem(df, short_hour_window):\n",
    "    df[\"rollingDem\"] = df[\"demand (MW)\"].rolling(\n",
    "        short_hour_window * 2, min_periods=1, center=True\n",
    "    ).median()\n",
    "    return df\n",
    "\n",
    "\n",
    "def add_rolling_dem_long(df, nDays):\n",
    "    df[\"rollingDemLong\"] = df[\"demand (MW)\"].rolling(\n",
    "        nDays * 24 * 2, min_periods=1, center=True\n",
    "    ).median()\n",
    "    return df\n",
    "\n",
    "\n",
    "def add_demand_minus_rolling_dem(df):\n",
    "    diff = df['demand (MW)'] - df['rollingDem']\n",
    "    df = df.assign(dem_minus_rolling=diff)\n",
    "    return df\n",
    "\n",
    "\n",
    "def add_demand_rel_diff_wrt_hourly(df):\n",
    "    diff = df['demand (MW)'] / (df['rollingDem'] * df['hourly_median_dem_dev'])\n",
    "    df = df.assign(dem_rel_diff_wrt_hourly=diff)\n",
    "    diff2 = df['demand (MW)'] / (df['rollingDemLong'] * df['hourly_median_dem_dev'])\n",
    "    df = df.assign(dem_rel_diff_wrt_hourly_long=diff2)\n",
    "    return df\n",
    "\n",
    "\n",
    "def add_delta_demand_rel_diff_wrt_hourly(df):\n",
    "    diff3 = df['dem_rel_diff_wrt_hourly'].diff()\n",
    "    df = df.assign(dem_rel_diff_wrt_hourly_delta_pre=diff3)\n",
    "    diff4 = df['dem_rel_diff_wrt_hourly'].diff(periods=-1)\n",
    "    df = df.assign(dem_rel_diff_wrt_hourly_delta_post=diff4)\n",
    "    return df\n",
    "\n",
    "\n",
    "# This is a global value so does not need to be added to the df\n",
    "def calculate_relative_demand_difference_IQR(df):\n",
    "    iqr_relative_deltas = np.nanpercentile(df['dem_rel_diff_wrt_hourly_delta_pre'], 75) - \\\n",
    "            np.nanpercentile(df['dem_rel_diff_wrt_hourly_delta_pre'], 25)\n",
    "    return iqr_relative_deltas\n",
    "\n",
    "    \n",
    "def add_demand_minus_rolling_dem_iqr(df, irq_hours):\n",
    "    rolling_window = df[\"dem_minus_rolling\"].rolling(iqr_hours * 2, min_periods=1, center=True)\n",
    "    df[\"dem_minus_rolling_IQR\"] = rolling_window.quantile(0.75) - rolling_window.quantile(0.25)\n",
    "    return df\n",
    "\n",
    "\n",
    "def add_hourly_median_dem_deviations(df):\n",
    "    # Create a df to hold all values to take nanmedian later\n",
    "    vals_dem_minus_rolling = df['dem_minus_rolling']\n",
    "    # Loop over nDays days on each side\n",
    "    for i in range(-nDays, nDays+1):\n",
    "        # Already initialized with zero value\n",
    "        if i == 0:\n",
    "            continue\n",
    "        vals_dem_minus_rolling = pd.concat(\n",
    "            [vals_dem_minus_rolling, df.shift(periods=i*24)['dem_minus_rolling']], axis=1)\n",
    "\n",
    "    df['vals_dem_minus_rolling'] = vals_dem_minus_rolling.median(axis=1, skipna=True)\n",
    "    # 1+vals to make it a scale factor\n",
    "    return df.assign(hourly_median_dem_dev=1.+df['vals_dem_minus_rolling']/df['rollingDemLong'])\n",
    "\n",
    "\n",
    "def add_deltas(df):\n",
    "    diff = df['demand (MW)'].diff()\n",
    "    df = df.assign(delta_pre=diff)\n",
    "    diff = df['demand (MW)'].diff(periods=-1)\n",
    "    df = df.assign(delta_post=diff)\n",
    "    return df\n",
    "\n",
    "\n",
    "def add_rolling_delta_iqr(df, iqr_hours):\n",
    "    rolling_window = df[\"delta_pre\"].rolling(iqr_hours * 2, min_periods=1, center=True)\n",
    "    df[\"delta_rolling_IQR\"] = rolling_window.quantile(0.75) - rolling_window.quantile(0.25)\n",
    "    return df\n",
    "\n",
    "\n",
    "def add_categories(df):\n",
    "    df['category'] = np.where(df['demand (MW)'].isna(), 'MISSING', 'OKAY')\n",
    "    return df\n",
    "\n",
    "\n",
    "def filter_neg_and_zeros(df):\n",
    "    df['category'] = np.where(df['demand (MW)'] <= 0., 'NEG_OR_ZERO', df['category'])\n",
    "    filtered = np.where(df['demand (MW)'] <= 0., df['demand (MW)'], np.nan)\n",
    "    df['negAndZeroFiltered'] = filtered\n",
    "    df['demand (MW)'] = df['demand (MW)'].mask(df['demand (MW)'] <= 0.)\n",
    "    return df\n",
    "\n",
    "    \n",
    "def filter_extrem_demand(df, multiplier):\n",
    "    med = np.nanmedian(df['demand (MW)'])\n",
    "    filtered = df['demand (MW)'].where(df['demand (MW)'] < med * multiplier)\n",
    "    df['globalDemandFiltered'] = np.where(df['demand (MW)'] != filtered, df['demand (MW)'], np.nan)\n",
    "    df['category'] = df['category'].mask(((df['demand (MW)'] != filtered) & \\\n",
    "                    (df['demand (MW)'].notna())), other='GLOBAL_DEM')\n",
    "    df['demand (MW)'] = filtered\n",
    "    return df\n",
    "\n",
    "\n",
    "def filter_global_plus_minus_one(df):\n",
    "    globalDemPlusMinusFiltered = [np.nan for _ in df.index]\n",
    "    for idx in df.index:\n",
    "        if df.loc[idx, 'category'] == 'GLOBAL_DEM':\n",
    "            if df.loc[idx-1, 'category'] == 'OKAY':\n",
    "                df.loc[idx-1, 'category'] = 'GLOBAL_DEM_PLUS_MINUS'\n",
    "                globalDemPlusMinusFiltered[idx-1] = df.loc[idx-1, 'demand (MW)']\n",
    "                df.loc[idx-1, 'demand (MW)'] = np.nan\n",
    "            if df.loc[idx+1, 'category'] == 'OKAY':\n",
    "                df.loc[idx+1, 'category'] = 'GLOBAL_DEM_PLUS_MINUS'\n",
    "                globalDemPlusMinusFiltered[idx+1] = df.loc[idx+1, 'demand (MW)']\n",
    "                df.loc[idx+1, 'demand (MW)'] = np.nan\n",
    "    df['globalDemPlusMinusFiltered'] = globalDemPlusMinusFiltered\n",
    "    return df\n",
    "    \n",
    "\n",
    "def filter_local_demand(df, multiplier_up, multiplier_down):\n",
    "    \n",
    "    # Filter in two steps to provide different labels for the categories\n",
    "    filtered = df['demand (MW)'].where(\n",
    "            (df['demand (MW)'] < df['rollingDem'] * df['hourly_median_dem_dev'] + \\\n",
    "                     multiplier_up * df['dem_minus_rolling_IQR']))\n",
    "    df['localDemandFilteredUp'] = np.where(df['demand (MW)'] != filtered, df['demand (MW)'], np.nan)\n",
    "    df['category'] = df['category'].mask(((df['demand (MW)'] != filtered) & \\\n",
    "                    (df['demand (MW)'].notna())), other='LOCAL_DEM_UP')\n",
    "    df['demand (MW)'] = filtered\n",
    "    \n",
    "    filtered = df['demand (MW)'].where(\n",
    "            (df['demand (MW)'] > df['rollingDem'] * df['hourly_median_dem_dev'] - \\\n",
    "                     multiplier_down * df['dem_minus_rolling_IQR']))\n",
    "    df['localDemandFilteredDown'] = np.where(df['demand (MW)'] != filtered, df['demand (MW)'], np.nan)\n",
    "    df['category'] = df['category'].mask(((df['demand (MW)'] != filtered) & \\\n",
    "                    (df['demand (MW)'].notna())), other='LOCAL_DEM_DOWN')\n",
    "    df['demand (MW)'] = filtered\n",
    "    \n",
    "    return df\n",
    "\n",
    "\n",
    "# Filter on a multiplier of the IQR and set\n",
    "# the associated 'demand (MW)' value to NAN.\n",
    "# Only consider \"double deltas\", hours with\n",
    "# large deltas on both sides\n",
    "def filter_deltas(df, multiplier):\n",
    "    \n",
    "    filtered = df['demand (MW)'].mask(\n",
    "            ((df['delta_pre'] > df['delta_rolling_IQR'] * multiplier) & \\\n",
    "            (df['delta_post'] > df['delta_rolling_IQR'] * multiplier)) | \\\n",
    "            ((df['delta_pre'] < -1. * df['delta_rolling_IQR'] * multiplier) & \\\n",
    "            (df['delta_post'] < -1. * df['delta_rolling_IQR'] * multiplier)))\n",
    "\n",
    "    df['deltaFiltered'] = np.where(df['demand (MW)'] != filtered, df['demand (MW)'], np.nan)\n",
    "    df['category'] = df['category'].mask(((df['demand (MW)'] != filtered) & \\\n",
    "                    (df['demand (MW)'].notna())), other='DELTA')\n",
    "    df['demand (MW)'] = filtered\n",
    "    return df\n",
    "\n",
    "\n",
    "# March through all hours recording previous \"good\"\n",
    "# demand value and its index.  Calculate deltas between\n",
    "# this value and next \"good\" hour.  If delta is LARGE\n",
    "# mark NAN.\n",
    "# Go forwards once, then backwards once to get all options.\n",
    "def filter_single_sided_deltas(df, multiplier, rel_multiplier, iqr_relative_deltas):\n",
    "\n",
    "\n",
    "    # Go through forwards first, then reverse\n",
    "    prev_good_index = np.nan\n",
    "    \n",
    "    deltaSingleFiltered = []\n",
    "    for idx in df.index:\n",
    "        deltaSingleFiltered.append(np.nan)\n",
    "        if np.isnan(df.loc[idx, 'demand (MW)']):\n",
    "            continue\n",
    "        \n",
    "        \n",
    "        # Initialize first good entry, this will never be flagged\n",
    "        if np.isnan(prev_good_index):\n",
    "            prev_good_index = idx\n",
    "            \n",
    "        \n",
    "        # Check deltas demand and relative wrt hourly adjustment\n",
    "        prev_good_delta_dem = abs(df.loc[prev_good_index, 'demand (MW)'] - df.loc[idx, 'demand (MW)'])\n",
    "        prev_good_delta_dem_rel_diff_wrt_hourly = abs(df.loc[prev_good_index, \n",
    "                        'dem_rel_diff_wrt_hourly'] - df.loc[idx, 'dem_rel_diff_wrt_hourly'])\n",
    "\n",
    "        \n",
    "        # delta_rolling_IQR is over 5 days on each side so should be\n",
    "        # similar regardless of which hours' we use. If delta is\n",
    "        # large, mark this hour anomalous\n",
    "        if (prev_good_delta_dem > df.loc[idx, 'delta_rolling_IQR'] * multiplier) and \\\n",
    "                (prev_good_delta_dem_rel_diff_wrt_hourly > rel_multiplier * iqr_relative_deltas):\n",
    "            \n",
    "            \n",
    "            # If the previous \"good\" value was farther from expected values, then consider current hour good\n",
    "            # and the previous hour will be caught on the way back through the reverse direction.\n",
    "            # The max deviation from the rolling 4 day dem and the rolling 10 day dem is taken\n",
    "            # to help catch cases where a large deviation pulls the rolling 4 day dem to center\n",
    "            # on its values.  i.e. SCL 2016 Dec 15.\n",
    "            prev_max = max(abs(1. - df.loc[prev_good_index, 'dem_rel_diff_wrt_hourly']),\n",
    "                            abs(1. - df.loc[prev_good_index, 'dem_rel_diff_wrt_hourly_long']))\n",
    "            current_max = max(abs(1. - df.loc[idx, 'dem_rel_diff_wrt_hourly']),\n",
    "                            abs(1. - df.loc[idx, 'dem_rel_diff_wrt_hourly_long']))         \n",
    "            if abs(current_max) < abs(prev_max):\n",
    "                prev_good_index = idx\n",
    "            \n",
    "            # else, continue to filter this hour\n",
    "            else:\n",
    "                deltaSingleFiltered[-1] = df.loc[idx, 'demand (MW)']\n",
    "                df.loc[idx, 'demand (MW)'] = np.nan\n",
    "                df.loc[idx, 'category'] = 'SINGLE_DELTA'\n",
    "        else:\n",
    "            prev_good_index = idx\n",
    "\n",
    "    \n",
    "    df['deltaSingleFilteredFwd'] = deltaSingleFiltered\n",
    "    \n",
    "    \n",
    "    ### Go through reversed, ~ copy of above code ###\n",
    "    prev_good_index = np.nan\n",
    "    \n",
    "    deltaSingleFiltered = []\n",
    "    for idx in reversed(df.index):\n",
    "        deltaSingleFiltered.append(np.nan)\n",
    "        if np.isnan(df.loc[idx, 'demand (MW)']):\n",
    "            continue\n",
    "        \n",
    "        \n",
    "        # Initialize first good entry, this will never be flagged\n",
    "        if np.isnan(prev_good_index):\n",
    "            prev_good_index = idx\n",
    "            \n",
    "        \n",
    "        # Check deltas demand and relative wrt hourly adjustment\n",
    "        prev_good_delta_dem = abs(df.loc[prev_good_index, 'demand (MW)'] - df.loc[idx, 'demand (MW)'])\n",
    "        prev_good_delta_dem_rel_diff_wrt_hourly = abs(df.loc[prev_good_index, \n",
    "                        'dem_rel_diff_wrt_hourly'] - df.loc[idx, 'dem_rel_diff_wrt_hourly'])\n",
    "\n",
    "        \n",
    "        # delta_rolling_IQR is over 5 days on each side so should be\n",
    "        # similar regardless of which hours' we use. If delta is\n",
    "        # large, mark this hour anomalous\n",
    "        if (prev_good_delta_dem > df.loc[idx, 'delta_rolling_IQR'] * multiplier) and \\\n",
    "                (prev_good_delta_dem_rel_diff_wrt_hourly > rel_multiplier * iqr_relative_deltas):\n",
    "            \n",
    "            \n",
    "            deltaSingleFiltered[-1] = df.loc[idx, 'demand (MW)']\n",
    "            df.loc[idx, 'demand (MW)'] = np.nan\n",
    "            df.loc[idx, 'category'] = 'SINGLE_DELTA'\n",
    "        else:\n",
    "            prev_good_index = idx\n",
    "\n",
    "    to_app = [val for val in reversed(deltaSingleFiltered)]\n",
    "    df['deltaSingleFilteredBkw'] = to_app\n",
    "    \n",
    "    return df\n",
    "\n",
    "\n",
    "def filter_runs(df):\n",
    "    \n",
    "    d1 = df['demand (MW)'].diff(periods=1)\n",
    "    d2 = df['demand (MW)'].diff(periods=2)\n",
    "\n",
    "    # cannot compare a dtyped [float64] array with a scalar of type [bool]\n",
    "    filtered = df['demand (MW)'].mask((d1 == 0) & (d2 == 0))\n",
    "    df['runFiltered'] = np.where(df['demand (MW)'] != filtered, df['demand (MW)'], np.nan)\n",
    "    df['demand (MW)'] = filtered\n",
    "    df['category'] = np.where(df['runFiltered'].notna(), 'IDENTICAL_RUN', df['category'])\n",
    "    return df\n",
    "    \n",
    "\n",
    "def filter_anomalous_regions(df, width, anomalous_pct):\n",
    "    \n",
    "    percent_good_data_cnt = [0. for _ in df.index]\n",
    "    percent_good_data_pre = [0. for _ in df.index]\n",
    "    percent_good_data_post = [0. for _ in df.index]\n",
    "    df['len_good_data'] = [0 for _ in df.index]\n",
    "    data_quality_cnt = []\n",
    "    data_quality_short = []\n",
    "    start_good_data = np.nan\n",
    "    end_good_data = np.nan\n",
    "    for idx in df.index:\n",
    "            \n",
    "        # Remove the oldest item in the list\n",
    "        if len(data_quality_short) > width:\n",
    "            data_quality_short.pop(0)\n",
    "        if len(data_quality_cnt) > 2 * width:\n",
    "            data_quality_cnt.pop(0)\n",
    "        \n",
    "        # Add new item and don't count MISSING as 'bad' data\n",
    "        if df.loc[idx, 'category'] == 'OKAY' or df.loc[idx, 'category'] == 'MISSING':\n",
    "            data_quality_cnt.append(1)\n",
    "            data_quality_short.append(1)\n",
    "            # Track length of good data chunks\n",
    "            if np.isnan(start_good_data):\n",
    "                start_good_data = idx\n",
    "            end_good_data = idx\n",
    "        else:\n",
    "            data_quality_cnt.append(0)\n",
    "            data_quality_short.append(0)\n",
    "            # Fill in length of good data chunk\n",
    "            if not (np.isnan(start_good_data) or np.isnan(end_good_data)):\n",
    "                len_good = end_good_data - start_good_data + 1\n",
    "                df.loc[start_good_data:end_good_data, 'len_good_data'] = len_good\n",
    "            start_good_data = np.nan\n",
    "            end_good_data = np.nan\n",
    "\n",
    "        \n",
    "        # centered measurements have length 2 * width\n",
    "        if len(data_quality_cnt) > 2 * width:\n",
    "            percent_good_data_cnt[idx-width] = np.mean(data_quality_cnt)\n",
    "        # left and right / pre and post measurements have length = width + 1\n",
    "        if len(data_quality_short) > width:\n",
    "            percent_good_data_pre[idx] = np.mean(data_quality_short)\n",
    "            percent_good_data_post[idx-width] = np.mean(data_quality_short)\n",
    "\n",
    "\n",
    "\n",
    "    \n",
    "    anomalousRegionsFiltered = [np.nan for _ in df.index]\n",
    "    for idx in df.index:\n",
    "        if percent_good_data_cnt[idx] <= anomalous_pct:\n",
    "            for j in range(idx-width, idx+width):\n",
    "                if j < 1 or j >= len(df.index):\n",
    "                    continue\n",
    "                if df.loc[j, 'category'] == 'OKAY':\n",
    "                    # If this is the start or end of continuous good data, don't filter\n",
    "                    if percent_good_data_pre[j] == 1.0 or percent_good_data_post[j] == 1.0:\n",
    "                        continue\n",
    "                    if df.loc[j, 'len_good_data'] > width:\n",
    "                        continue\n",
    "                    df.loc[j, 'category'] = 'ANOMALOUS_REGION'\n",
    "                    anomalousRegionsFiltered[j] = df.loc[j, 'demand (MW)']\n",
    "                    df.loc[j, 'demand (MW)'] = np.nan\n",
    "    \n",
    "\n",
    "    df['anomalousRegionsFiltered'] = anomalousRegionsFiltered\n",
    "    return df\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def return_all_regions():\n",
    "    return ['AEC', 'AECI', 'CPLE', 'CPLW',\n",
    "    'DUK', 'FMPP', 'FPC',\n",
    "    'FPL', 'GVL', 'HST', 'ISNE',\n",
    "    'JEA', 'LGEE', 'MISO', 'NSB',\n",
    "    'NYIS', 'OVEC', 'PJM', 'SC',\n",
    "    'SCEG', 'SEC', 'SOCO',\n",
    "    'SPA', 'SWPP', 'TAL', 'TEC',\n",
    "    'TVA', 'ERCO',\n",
    "    'AVA', 'AZPS', 'BANC', 'BPAT',\n",
    "    'CHPD', 'CISO', 'DOPD',\n",
    "    'EPE', 'GCPD', 'IID',\n",
    "    'IPCO', 'LDWP', 'NEVP', 'NWMT',\n",
    "    'PACE', 'PACW', 'PGE', 'PNM',\n",
    "    'PSCO', 'PSEI', 'SCL', 'SRP',\n",
    "    'TEPC', 'TIDC', 'TPWR', 'WACM',\n",
    "    'WALC', 'WAUW']\n",
    "\n",
    "# append subregions to regions list\n",
    "def append_sub_regions(l):\n",
    "    subregions = [\n",
    "        'CISO-PGAE', 'CISO-SCE', 'CISO-SDGE', 'CISO-VEA',\n",
    "\n",
    "        'ERCO-COAS', 'ERCO-EAST', 'ERCO-FWES', 'ERCO-NCEN',\n",
    "        'ERCO-NRTH', 'ERCO-SCEN', 'ERCO-SOUT', 'ERCO-WEST',\n",
    "\n",
    "        'ISNE-4001', 'ISNE-4002', 'ISNE-4003', 'ISNE-4004',\n",
    "        'ISNE-4005', 'ISNE-4006', 'ISNE-4007', 'ISNE-4008',\n",
    "\n",
    "        'MISO-0001', 'MISO-0004', 'MISO-0006', 'MISO-0027', 'MISO-0035', 'MISO-8910',\n",
    "\n",
    "        'NYIS-ZONA', 'NYIS-ZONB', 'NYIS-ZONC', 'NYIS-ZOND', 'NYIS-ZONE',\n",
    "        'NYIS-ZONF', 'NYIS-ZONG', 'NYIS-ZONH', 'NYIS-ZONI', 'NYIS-ZONJ', 'NYIS-ZONK',\n",
    "\n",
    "        'PJM-AE', 'PJM-AEP', 'PJM-AP', 'PJM-ATSI', 'PJM-BC', 'PJM-CE', 'PJM-DAY',\n",
    "        'PJM-DEOK', 'PJM-DOM', 'PJM-DPL', 'PJM-DUQ', 'PJM-EKPC', 'PJM-JC',\n",
    "        'PJM-ME', 'PJM-PE', 'PJM-PEP', 'PJM-PL', 'PJM-PN', 'PJM-PS', 'PJM-RECO',\n",
    "\n",
    "        'PNM-KAFB', 'PNM-KCEC', 'PNM-LAC', 'PNM-NTUA', 'PNM-PNM', 'PNM-TSGT',\n",
    "\n",
    "        'SWPP-CSWS', 'SWPP-EDE', 'SWPP-GRDA', 'SWPP-INDN', 'SWPP-KACY',\n",
    "        'SWPP-KCPL', 'SWPP-LES', 'SWPP-MPS', 'SWPP-NPPD', 'SWPP-OKGE', 'SWPP-OPPD',\n",
    "        'SWPP-SECI', 'SWPP-SPRM', 'SWPP-SPS', 'SWPP-WAUE', 'SWPP-WFEC', 'SWPP-WR',    \n",
    "    ]\n",
    "    for subR in subregions:\n",
    "        l.append(subR)\n",
    "    return l\n",
    "\n",
    "# Remobes BAs with subregions to avoid double counting in the imputation step.\n",
    "# This excludes the ERCO subregions because of their 2019 data start date.\n",
    "def remove_regions(l):\n",
    "    \n",
    "    to_exclude = [\n",
    "        # because of 2019 subregion start dates\n",
    "        'ERCO-COAS', 'ERCO-EAST', 'ERCO-FWES', 'ERCO-NCEN',\n",
    "        'ERCO-NRTH', 'ERCO-SCEN', 'ERCO-SOUT', 'ERCO-WEST',\n",
    "        \n",
    "        # because subregions are so small there are too many\n",
    "        # identical run filters that result in many anomalous regions\n",
    "        'PNM-KAFB', 'PNM-KCEC', 'PNM-LAC', 'PNM-NTUA', 'PNM-PNM', 'PNM-TSGT',\n",
    "        \n",
    "        # This PJM region has last data entry April 2019\n",
    "        'PJM-RECO',\n",
    "        \n",
    "        # BAs with subregions\n",
    "        'CISO', 'ISNE', 'MISO', 'NYIS', 'PJM', 'SWPP',\n",
    "    ]\n",
    "    \n",
    "    for rm in to_exclude:\n",
    "        if rm in l:\n",
    "            l.remove(rm)\n",
    "    \n",
    "    return l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Screening selection parameters\n",
    "For a detailed description of each, please see the associated Data Descriptor linked in the repository's README."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "short_hour_window = 24 # 48 hour moving median (M_{t,48hr})\n",
    "iqr_hours = 24*5 # width in hours of IQR values of relative deviations from diurnal cycle template (IQR_{dem,t})\n",
    "nDays = 10 # Used for normalized hourly demand template (h_{t,diurnal}) and 480 hour moving median (M_{t,480hr})\n",
    "global_dem_cut = 10 # threshold selection for global demand filter\n",
    "local_dem_cut_up = 3.5 # upwards threshold for local demand filter\n",
    "local_dem_cut_down = 2.5 # downwards threshold for local demand filter\n",
    "delta_multiplier = 2 # selection threshold for double-sided delta filter\n",
    "delta_single_multiplier = 5 # selection threshold for single-sided delta filter\n",
    "rel_multiplier = 15 # other selection threshold for single-sided delta filter\n",
    "anomalous_regions_width = 24 # width in hours of anomalous region filter\n",
    "anomalous_pct = .85 # required pct of good data in anomalous region filter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run the screening code over each EIA balancing authority"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "regions = return_all_regions()\n",
    "\n",
    "if clean_subregions:\n",
    "    regions = append_sub_regions(regions)\n",
    "\n",
    "\n",
    "regions.sort()\n",
    "for region in regions:\n",
    "    print(f'Processing {region}')\n",
    "\n",
    "\n",
    "    df = pd.read_csv(f'{input_path}/{region}.csv',\n",
    "                    dtype={'demand (MW)':np.float64},\n",
    "                    na_values=['MISSING', 'EMPTY'])\n",
    "\n",
    "\n",
    "    # Drop unused columns\n",
    "    df = df.drop(['forecast demand (MW)'], axis=1)\n",
    "        \n",
    "    # Convert date/time\n",
    "    df['date_time'] = pd.to_datetime(df['date_time'])\n",
    "        \n",
    "    # Add category labels to track which algo screens an hourly value\n",
    "    df = add_categories(df)\n",
    "\n",
    "    # Mark missing and empty values\n",
    "    df = df.assign(missing=df['demand (MW)'].isna())\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "    #---------------------------------------------\n",
    "    # Screening Step 1\n",
    "    #---------------------------------------------\n",
    "    \n",
    "    # Set all negative and zero values to NAN\n",
    "    # (negative or zero filter)\n",
    "    df = filter_neg_and_zeros(df)\n",
    "\n",
    "    # Set last demand values in runs of 3+ to NAN\n",
    "    # (identical run filter)\n",
    "    df = filter_runs(df)\n",
    "\n",
    "    # Global demand filter on 10x the median value\n",
    "    # (global demand filter)\n",
    "    df = filter_extrem_demand(df, global_dem_cut)\n",
    "\n",
    "    # Filter +/- 1 hour from any global deman filtered hours\n",
    "    # (global demand plus/minus 1 hour filter)\n",
    "    df = filter_global_plus_minus_one(df)\n",
    "\n",
    "\n",
    "\n",
    "    \n",
    "    #---------------------------------------------\n",
    "    # Calculate demand characteristics for Step 2\n",
    "    #---------------------------------------------\n",
    "    \n",
    "    # 48 hour moving median (M_{t,48hr})\n",
    "    df = add_rolling_dem(df, short_hour_window)\n",
    "\n",
    "    # 480 hour moving median (M_{t,480hr})\n",
    "    df = add_rolling_dem_long(df, nDays)\n",
    "\n",
    "    # demand minus moving median (Delta(d_{t},M_{t,48hr}))\n",
    "    df = add_demand_minus_rolling_dem(df)\n",
    "\n",
    "    # IQR values of relative deviations from diurnal cycle template (IQR_{dem,t})\n",
    "    df = add_demand_minus_rolling_dem_iqr(df, iqr_hours)\n",
    "\n",
    "    # demand deltas (delta(d_{t-1},d_{t}))\n",
    "    df = add_deltas(df)\n",
    "\n",
    "    # IQR values of demand deltas (IQR_{delta,t})\n",
    "    df = add_rolling_delta_iqr(df, iqr_hours)\n",
    "\n",
    "    # normalized hourly demand template (h_{t,diurnal})\n",
    "    df = add_hourly_median_dem_deviations(df)\n",
    "\n",
    "    # Demand deviation from hourly diurnal template (r_{t})\n",
    "    # This adds both the short and long moving medians\n",
    "    df = add_demand_rel_diff_wrt_hourly(df)\n",
    "\n",
    "    # Hour-to-hour differences between hourly diurnal templates\n",
    "    # (delta(r_{t-1},r_{t}))\n",
    "    # This adds differences for both the short and long moving medians\n",
    "    df = add_delta_demand_rel_diff_wrt_hourly(df)\n",
    "\n",
    "    # Calculate the global IQR for the hour-to-hour differences \n",
    "    # between hourly diurnal templates (IQR_{r})\n",
    "    # This is a global value and is not added to the dataframe\n",
    "    iqr_relative_deltas = calculate_relative_demand_difference_IQR(df)\n",
    "\n",
    "\n",
    "    \n",
    "    #---------------------------------------------\n",
    "    # Screening Step 2\n",
    "    #---------------------------------------------\n",
    "    \n",
    "    # (local demand filter)\n",
    "    df = filter_local_demand(df, local_dem_cut_up, local_dem_cut_down)\n",
    "    \n",
    "    # (double-sided delta filter)\n",
    "    df = filter_deltas(df, delta_multiplier)\n",
    "    \n",
    "    # (single-sided delta filter)\n",
    "    df = filter_single_sided_deltas(df, delta_single_multiplier,\n",
    "                                rel_multiplier, iqr_relative_deltas)\n",
    "    \n",
    "    # (anomalous regions filter)\n",
    "    df = filter_anomalous_regions(df, anomalous_regions_width, anomalous_pct)\n",
    "    \n",
    "    # Save as csv for easy viewing and pickle for computatinal ease\n",
    "    print(f'Saving pickle {input_path}/pickle_{region}.pkl')\n",
    "    pickle_file = open(f'{input_path}/pickle_{region}.pkl', 'wb') \n",
    "    pickle.dump(df, pickle_file)\n",
    "    pickle_file.close()\n",
    "    df.to_csv(f'{input_path}/csv_{region}.csv')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prepare a csv File for the MICE Imputation Process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "prep_final_output = True\n",
    "print(f\"prep_final_output {prep_final_output}\")\n",
    "\n",
    "regions = return_all_regions()\n",
    "\n",
    "if clean_subregions:\n",
    "    regions = append_sub_regions(regions)\n",
    "    regions = remove_regions(regions)\n",
    "\n",
    "print(f\"Length of REGIONS: {len(regions)}\")\n",
    "\n",
    "regions.sort()\n",
    "print(regions)\n",
    "for i, region in enumerate(regions):\n",
    "    if not prep_final_output:\n",
    "        break\n",
    "    print(f'Loading from pickle pickle_{region}.pkl')\n",
    "    pickle_in = open(f'{input_path}/pickle_{region}.pkl','rb')\n",
    "    if i == 0: # Load first instance to master\n",
    "        master = pickle.load(pickle_in)\n",
    "        cols = master.columns.tolist()\n",
    "        cols.remove('date_time')\n",
    "        master['date_time'] = pd.to_datetime(master['date_time'])\n",
    "        master[region] = master['demand (MW)']\n",
    "        master[region+'_category'] = master['category']\n",
    "        master = master.drop(cols, axis=1)\n",
    "        continue\n",
    "        \n",
    "\n",
    "    df = pickle.load(pickle_in)\n",
    "    master[region] = df['demand (MW)']\n",
    "    master[region+'_category'] = df['category']\n",
    "\n",
    "if prep_final_output:\n",
    "    print(master.head(5))\n",
    "    master.to_csv(f'{input_path}/csv_MASTER.csv', index=False, na_rep='NA')\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}