conference_notebooks/KDD_2019/plasticc/Part_1-1_RNN_Feature_Engineering.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 1.1 - Feature Engineering for RNN Model\n",
    "In this notebook, we will analyze and transform timeseries features so that we can feed them through a pre-trained RNN model in order to turn our variable-length timeseries sequences into fixed-length vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "import math\n",
    "import cudf as gd\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from numba import cuda\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from utils import scatter\n",
    "%matplotlib inline\n",
    "\n",
    "print(gd.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PATH = \"../../../../../data/plasticc_data\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read data\n",
    "\n",
    "To begin, we'll look at how to engineer features from the raw training data. We'll load it in from a CSV format, and represent it as a Dataframe.\n",
    "\n",
    "Thanks to the `read_csv` API, we can manually define all column names, types, etc., to optimize the reading/loading process. We'll have to skip the record header at the beginning of the file. This can be done by setting `skiprows=1`.\n",
    "\n",
    "After the read is complete, we can drop the columns we don't need to save memory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts_cols = ['object_id', 'mjd', 'passband', 'flux', 'flux_err', 'detected']\n",
    "ts_dtypes = ['int32', 'float32', 'int32', 'float32','float32','int32']\n",
    "\n",
    "train_gd = gd.read_csv('%s/training_set.csv'%PATH, names=ts_cols,dtype=ts_dtypes,skiprows=1)\n",
    "test_gd = gd.read_csv('%s/test_set_sample.csv'%PATH, names=ts_cols,dtype=ts_dtypes,skiprows=1)\n",
    "\n",
    "train_gd.drop_column('detected')\n",
    "test_gd.drop_column('detected')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_gd.head().to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normalization\n",
    "\n",
    "Let's histogram the `flux` and `flux_err` values to get a feel for their distributions. \n",
    "\n",
    "We can use `hist()` from Pandas to generate histograms. To access the Pandas API from cuDF, we'll have to call `to_pandas()` to copy the data to the host and convert it to a Pandas Dataframe.\n",
    "\n",
    "A quick look at the output suggests that most of the binned values occur near the zeroth bin. To exemplify the shape of the distribution, we can perform a logarithmic transformation.\n",
    "\n",
    "A couple more details worth noting here:\n",
    "- There are extreme negative values in `flux`, which cannot be log-transformed. We will need to clip them.\n",
    "- The `flux_err` is entirely non-negative; so, we can apply log-transform directly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1,2, figsize=(20, 5))\n",
    "ax = train_gd.flux.to_pandas().hist(ax=axs[0])\n",
    "\n",
    "ax.set_xlabel('flux', fontsize=18)\n",
    "ax.set_ylabel('Count', fontsize=18)\n",
    "ax.set_title('Original data')\n",
    "\n",
    "ax = train_gd.flux_err.to_pandas().hist(ax=axs[1])\n",
    "\n",
    "ax.set_xlabel('flux_err', fontsize=18)\n",
    "ax.set_ylabel('Count', fontsize=18)\n",
    "ax.set_title('Original data')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `applymap(func)` function can be used on cuDF Series objects. \n",
    "\n",
    "This enables us to apply the log transform and clipped log transform in parallel on the GPU for all the values in the Series. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def clip_and_log_transform(df):\n",
    "    df['flux'] = df['flux'].applymap(lambda x: math.log1p(x+10) if x>-10 else 0)\n",
    "    df['flux_err'] = df['flux_err'].applymap(lambda x: math.log1p(x))\n",
    "    return df\n",
    "\n",
    "train_gd = clip_and_log_transform(train_gd)\n",
    "test_gd = clip_and_log_transform(test_gd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's plot the histograms of the clipped and log-transformed values to get a better idea of the data distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1,2, figsize=(20, 5))\n",
    "ax = train_gd.flux.to_pandas().hist(ax=axs[0])\n",
    "ax.set_xlabel('flux', fontsize=18)\n",
    "ax.set_ylabel('Count', fontsize=18)\n",
    "ax.set_title('After clipping and log transformation')\n",
    "\n",
    "ax = train_gd.flux_err.to_pandas().hist(ax=axs[1])\n",
    "ax.set_xlabel('flux_err', fontsize=18)\n",
    "ax.set_ylabel('Count', fontsize=18)\n",
    "ax.set_title('After log transformation')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Differentiation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First let's take look at the scatter plot of the light curve for a single object.\n",
    "\n",
    "We want to learn how `flux` varies over timestamp `mjd`, namely $\\frac{\\nabla flux}{\\nabla mjd}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "oid = 615\n",
    "sample = train_gd[train_gd['object_id']==oid].to_pandas()\n",
    "scatter(sample['mjd'].values,\n",
    "        sample['flux'].values,\n",
    "        values=sample['passband'].values,\n",
    "        xlabel='time',ylabel='flux',title='object %d class 42'%oid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_gd = train_gd.sort_values(['object_id','mjd'])\n",
    "test_gd = test_gd.sort_values(['object_id','mjd'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use Panda's `groupby() - apply()` technique to group the dataframes by `object_id`. \n",
    "\n",
    "Rather than using `apply()` on the grouped Dataframe, we use cuDF's `apply_grouped()` so that we can run a function in parallel on each group. Our custom function will compute the deltas between timestamps.\n",
    "\n",
    "The following function, `compute_delta()` computes the differences for `flux` and `mjd` observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def delta_in_group(val,val_delta, default):\n",
    "    \"\"\"\n",
    "    A Numba CUDA function that diffs all the values in a single timeseries sequence\n",
    "    \"\"\"\n",
    "    for i in range(cuda.threadIdx.x, len(val), cuda.blockDim.x):\n",
    "        if i>0:\n",
    "            val_delta[i] = val[i]-val[i-1]\n",
    "        else:\n",
    "            val_delta[i] = default\n",
    "\n",
    "def compute_delta(df):\n",
    "    for col,d in zip(['flux','mjd'],[0,180]):\n",
    "        df = df.rename({col:'val'})\n",
    "        df = df.groupby('object_id',method=\"cudf\",\n",
    "                        as_index=False).apply_grouped(delta_in_group,\n",
    "                                  incols=['val'],\n",
    "                                  outcols={'val_delta': np.float32},\n",
    "                                  kwargs={'default':d},\n",
    "                                  tpb=32)\n",
    "        df = df.rename({'val': col})\n",
    "        df = df.rename({'val_delta':'%s_delta'%col})\n",
    "\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_gd = compute_delta(train_gd)\n",
    "test_gd = compute_delta(test_gd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_gd.drop_column('mjd')\n",
    "test_gd.drop_column('mjd')\n",
    "\n",
    "train_gd.drop_column('flux')\n",
    "test_gd.drop_column('flux')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_gd.head().to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the variance of `mjd_delta` and `flux_delta`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1,1)\n",
    "ax = train_gd.mjd_delta.to_pandas().hist()\n",
    "ax.set_xlabel('mjd_delta', fontsize=18)\n",
    "ax.set_ylabel('Count', fontsize=18)\n",
    "ax.set_title('Original')\n",
    "\n",
    "fig, axs = plt.subplots(1,1)\n",
    "ax = train_gd.flux_delta.to_pandas().hist()\n",
    "ax.set_xlabel('flux_delta', fontsize=18)\n",
    "ax.set_ylabel('Count', fontsize=18)\n",
    "ax.set_title('Original')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Independent Exercise\n",
    "\n",
    "The logarithmic transform lowered the variance of our `flux` values above, so `flux_delta` contains values with a fairly low variance. The variance of `mjd_delta` is pretty large. Lowering this will help to improve accuracy in our final predictions.\n",
    "\n",
    "Use the blank cell below to scalle all the values in the `mjd_delta` column by the value `90` to shrink the variance. This will give better results when training the RNN.\n",
    "\n",
    "Make sure both `train_gd` and `test_gd` are scaled"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = train_gd.mjd_delta.to_pandas().hist()\n",
    "ax.set_xlabel('mjd_delta', fontsize=18)\n",
    "ax.set_ylabel('Count', fontsize=18)\n",
    "ax.set_title('After scaling')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create step numbers for each light curve sequence\n",
    "\n",
    "For our final engineering step, let's mark each timeseries element with a step number, so that we can preserve the order of each sequence. \n",
    "\n",
    "Since the train and test Dataframes are already ordered by their object ids and the time steps, we can just number the ocurrences of object ids starting at 0. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_gd.head().to_pandas()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_pd = train_gd.to_pandas()\n",
    "test_pd = test_gd.to_pandas()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we have converted our reference dataframe to host memory and we only need the `object_id` column of our GPU Dataframes, we can loop through all the other columns and drop them. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cols = [col for col in train_gd.columns]\n",
    "for col in cols:\n",
    "    if col!='object_id':\n",
    "        train_gd.drop_column(col)\n",
    "        test_gd.drop_column(col)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can number the steps in the GPU Dataframe by using the `groupby() - apply_grouped()` strategy we used for computing the deltas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_order_in_group(row_id,step):\n",
    "    \"\"\"\n",
    "    A Numba CUDA function that assigns a step value to all the elements in a single timeseries sequence\n",
    "    \"\"\"\n",
    "    for i in range(cuda.threadIdx.x, len(row_id), cuda.blockDim.x):\n",
    "        step[i] = i\n",
    "\n",
    "def add_step(df):\n",
    "    df['row_id'] = np.arange(df.shape[0])\n",
    "    df = df.groupby('object_id',method=\"cudf\",\n",
    "                as_index=False).apply_grouped(get_order_in_group,\n",
    "                          incols=['row_id'],\n",
    "                          outcols={'step': np.float32},\n",
    "                          tpb=32)\n",
    "    df = df.sort_values('row_id')\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_gd = add_step(train_gd)\n",
    "test_gd = add_step(test_gd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will convert our final `Series` of numbered steps to to `Pandas.Series` and assign the steps to our Dataframes on host."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_pd['step'] = train_gd['step'].to_pandas().values\n",
    "test_pd['step'] = test_gd['step'].to_pandas().values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_pd.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_pd.to_pickle('train_rnn.pkl')\n",
    "test_pd.to_pickle('test_rnn.pkl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we process the light curve timseries data in order to train a RNN. \n",
    "- `object_id` can be used to sample sequences from the dataframe.\n",
    "- `mjd_delta` and `flux_delta` enable RNN to learn the temporal pattern of light curves.\n",
    "- `passband` is an original feature and will be embedded."
   ]
  },
  {
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}