blog_notebooks/plasticc/notebooks/rapids_lsst_gpu_only_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kaggle Comptetion: PLasTiCC LSST GPU Only Demo\n",
    "\n",
    "This is a 2 minutes end-to-end gpu-only demos the 8th place solution (8/1094) of Rapids.ai for the __[PLAsTiCC Astronomical Classification](https://www.kaggle.com/c/PLAsTiCC-2018/leaderboard)__.** More details can be found at our __[blog](https://medium.com/rapids-ai/make-sense-of-the-universe-with-rapids-ai-d105b0e5ec95)__\n",
    "\n",
    "**Note: this notebook is here for archival purposes and is not intended to illustrate best practices.  [Please use this updated PLAsTiCC notebook, shown at KDD 2019 instead](conference_notebooks/KDD_2019/plasticc)**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "GPU_id = 0\n",
    "os.environ['CUDA_VISIBLE_DEVICES'] = str(GPU_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cudf as gd\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import math\n",
    "import xgboost as xgb\n",
    "import seaborn as sns\n",
    "from functools import partial\n",
    "from sklearn.preprocessing import LabelEncoder\n",
    "from sklearn.model_selection import train_test_split\n",
    "from termcolor import colored\n",
    "from cudf_workaround import cudf_groupby_aggs\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "sns.set()\n",
    "print(gd.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Table of contents\n",
    "[1. Global variables](#global)<br>\n",
    "[2. Functions](#func)<br>\n",
    "[3. ETL & Visualizations](#etl)<br>\n",
    "[4. Model training](#train)<br>\n",
    "[5. Conclusions](#conclusions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"global\"></a>\n",
    "## 1. Global variables "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Original data download and description __[link](https://www.kaggle.com/c/PLAsTiCC-2018/data)__**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PATH = '../data'\n",
    "#PATH = '/raid/data/ml/lsst/input'\n",
    "#PATH = '../lsst/input'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Tested on V100 with 32 GB GPU memory. Please reset this variable if memory capacity is smaller, and the input data will be down sampled accordingly.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GPU_MEMORY = 32 # GB. \n",
    "#GPU_MEMORY = 16 # GB. Both 32 and 16 GB have been tested"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TEST_ROWS = 453653104 # number of rows in test data\n",
    "# no skip if your gpu has 32 GB memory\n",
    "# otherwise, skip rows porportionally\n",
    "OVERHEAD = 1.2 # cudf 0.7 introduces 20% memory overhead comparing to cudf 0.4\n",
    "SKIP_ROWS = int((1 - GPU_MEMORY/(32.0*OVERHEAD))*TEST_ROWS) \n",
    "GPU_RUN_TIME = {}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"func\"></a>\n",
    "## 2. Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def scatter(x,y,values,xlabel='x',ylabel='y',title=None):\n",
    "    colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k']\n",
    "    colors = np.array([colors[i] for i in values])\n",
    "    ps = []\n",
    "    bs = []\n",
    "    bands = ['passband_%s'%i for i in ['u', 'g', 'r', 'i', 'z','y']]\n",
    "    for i in sorted(np.unique(values)):\n",
    "        mask = values==i\n",
    "        if len(x[mask]):\n",
    "            p = plt.scatter(x[mask],y[mask],c=colors[mask])\n",
    "            ps.append(p)\n",
    "            bs.append(bands[i])\n",
    "    plt.legend(ps,bs,scatterpoints=1)\n",
    "    if title is not None:\n",
    "        plt.title(title)\n",
    "    \n",
    "    plt.xlim([np.min(x)-10,np.min(x)+1500])\n",
    "    plt.ylabel('y: %s'%ylabel)\n",
    "    plt.xlabel('x: %s'%xlabel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def multi_weighted_logloss(y_true, y_preds, classes, class_weights):\n",
    "    \"\"\"\n",
    "    refactor from\n",
    "    @author olivier https://www.kaggle.com/ogrellier\n",
    "    multi logloss for PLAsTiCC challenge\n",
    "    \"\"\"\n",
    "    y_p = y_preds.reshape(y_true.shape[0], len(classes), order='F')\n",
    "    y_ohe = pd.get_dummies(y_true)\n",
    "    y_p = np.clip(a=y_p, a_min=1e-15, a_max=1 - 1e-15)\n",
    "    y_p_log = np.log(y_p)\n",
    "    y_log_ones = np.sum(y_ohe.values * y_p_log, axis=0)\n",
    "    nb_pos = y_ohe.sum(axis=0).values.astype(float)\n",
    "    class_arr = np.array([class_weights[k] for k in sorted(class_weights.keys())])\n",
    "    y_w = y_log_ones * class_arr / nb_pos\n",
    "\n",
    "    loss = - np.sum(y_w) / np.sum(class_arr)\n",
    "    return loss\n",
    "\n",
    "def xgb_multi_weighted_logloss(y_predicted, y_true, classes, class_weights):\n",
    "    loss = multi_weighted_logloss(y_true.get_label(), y_predicted, \n",
    "                                  classes, class_weights)\n",
    "    return 'wloss', loss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### GPU ETL functions "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To save GPU memory, we drop the column as soon as it is done with groupby\n",
    "# this hits performance a little but avoids GPU OOM.\n",
    "def groupby_aggs(df,aggs,col):\n",
    "    res = None\n",
    "    for i,j in aggs.items():\n",
    "        for k in j:\n",
    "            #print(i,k)\n",
    "            tmp = df.groupby(col,as_index=False).agg({i:[k]})\n",
    "            if res is None:\n",
    "                res = tmp\n",
    "            else:\n",
    "                res = res.merge(tmp,on=[col],how='left')\n",
    "        df.drop_column(i)\n",
    "    return res\n",
    "\n",
    "def etl_gpu(df,df_meta):\n",
    "    aggs = {\n",
    "        'passband': ['mean'], \n",
    "        'detected': ['mean'],\n",
    "        'mjd':['max','min'],\n",
    "    }\n",
    "    agg_df = groupby_aggs(df,aggs,'object_id')\n",
    "    # at this step, columns ['passband','detected','mjd'] are deleted \n",
    "    \n",
    "    df['flux_ratio_sq'] = df['flux'] / df['flux_err']\n",
    "    df['flux_ratio_sq'] = df['flux_ratio_sq'].applymap(lambda x: math.pow(x,2))\n",
    "    df['flux_by_flux_ratio_sq'] = df['flux'] * df['flux_ratio_sq']\n",
    "    \n",
    "    aggs2 = {\n",
    "        'flux_ratio_sq':['sum'],\n",
    "        'flux_by_flux_ratio_sq':['sum'],\n",
    "        'flux': ['min', 'max', 'mean'],\n",
    "        'flux_err': ['min', 'max', 'mean'],\n",
    "    }\n",
    "    agg_df2 = groupby_aggs(df,aggs2,'object_id')\n",
    "    agg_df = agg_df.merge(agg_df2,on=['object_id'],how='left')\n",
    "    del agg_df2\n",
    "\n",
    "    agg_df['flux_diff'] = agg_df['max_flux'] - agg_df['min_flux']\n",
    "    agg_df['flux_dif2'] = (agg_df['max_flux'] - agg_df['min_flux']) / agg_df['mean_flux']\n",
    "    agg_df['flux_w_mean'] = agg_df['sum_flux_by_flux_ratio_sq'] / agg_df['sum_flux_ratio_sq']\n",
    "    agg_df['flux_dif3'] = (agg_df['max_flux'] - agg_df['min_flux']) / agg_df['flux_w_mean']\n",
    "    \n",
    "    agg_df['mjd_diff'] = agg_df['max_mjd'] - agg_df['min_mjd']\n",
    "    agg_df.drop_column('max_mjd')\n",
    "    agg_df.drop_column('min_mjd')\n",
    "    \n",
    "    for col in ['ra','decl','gal_l','gal_b']:\n",
    "        df_meta.drop_column(col)\n",
    "    \n",
    "    df_meta = df_meta.merge(agg_df,on=['object_id'],how='left')\n",
    "    return df_meta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"etl\"></a>\n",
    "## 3. ETL & Visualizations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data for ETL part 1\n",
    "**GPU load data**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "start = time.time()\n",
    "step = 'load data part1'\n",
    "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,\n",
    "            names=ts_cols,dtype=ts_dtypes,skiprows=1)\n",
    "test_gd = gd.read_csv('%s/test_set.csv'%PATH,\n",
    "            names=ts_cols,dtype=ts_dtypes,skiprows=1+SKIP_ROWS) # skip the header\n",
    "GPU_RUN_TIME[step] = time.time() - start"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "oid = 615\n",
    "train = train_gd.to_pandas()\n",
    "mask = train.object_id== oid\n",
    "scatter(train.loc[mask,'mjd'].values,\n",
    "                train.loc[mask,'flux'].values,\n",
    "                values=train.loc[mask,'passband'].values,\n",
    "                xlabel='time',ylabel='flux',title='object %d class 42'%oid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ETL part 1 with 100x  speedup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# to save memory, we need to move dataframe to cpu and only keep the columns we need\n",
    "test_gd = test_gd[['object_id','flux']]\n",
    "train_gd = train_gd[['object_id','flux']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# GPU\n",
    "step = 'ETL part1'\n",
    "start = time.time()\n",
    "aggs = {'flux':['skew']}\n",
    "test_gd = cudf_groupby_aggs(test_gd,group_id_col='object_id',aggs=aggs)\n",
    "train_gd = cudf_groupby_aggs(train_gd,group_id_col='object_id',aggs=aggs)\n",
    "GPU_RUN_TIME[step] = time.time() - start"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "test_flux_skew_gd = test_gd.sort_values(by='object_id')\n",
    "train_flux_skew_gd = train_gd.sort_values(by='object_id')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load data for the ETL part 2 with 11x speedup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# read data on gpu\n",
    "step = 'load data part2'\n",
    "start = time.time()\n",
    "ts_cols = ['object_id', 'mjd', 'passband', 'flux', 'flux_err', 'detected']\n",
    "ts_dtypes = ['int32', 'float32', 'int32', 'float32','float32','int32']\n",
    "\n",
    "test_gd = gd.read_csv('%s/test_set.csv'%PATH,\n",
    "            names=ts_cols,dtype=ts_dtypes,skiprows=1+SKIP_ROWS) # skip the header\n",
    "train_gd = gd.read_csv('%s/training_set.csv'%PATH,\n",
    "            names=ts_cols,dtype=ts_dtypes,skiprows=1)\n",
    "\n",
    "cols = ['object_id', 'ra', 'decl', 'gal_l', 'gal_b', 'ddf',\n",
    "       'hostgal_specz', 'hostgal_photoz', 'hostgal_photoz_err', \n",
    "       'distmod','mwebv', 'target']\n",
    "dtypes = ['int32']+['float32']*4+['int32']+['float32']*5+['int32']\n",
    "\n",
    "train_meta_gd = gd.read_csv('%s/training_set_metadata.csv'%PATH,\n",
    "            names=cols,dtype=dtypes,skiprows=1)\n",
    "del cols[-1],dtypes[-1]\n",
    "test_meta_gd = gd.read_csv('%s/test_set_metadata.csv'%PATH,\n",
    "            names=cols,dtype=dtypes,skiprows=1)\n",
    "GPU_RUN_TIME[step] = time.time() - start"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ETL part2 with 9x ~ 12x speedup "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# GPU\n",
    "start = time.time()\n",
    "step = 'ETL part2'\n",
    "train_final_gd = etl_gpu(train_gd,train_meta_gd)\n",
    "train_final_gd = train_final_gd.merge(train_flux_skew_gd,on=['object_id'],how='left')\n",
    "test_final_gd = etl_gpu(test_gd,test_meta_gd)\n",
    "del test_gd,test_meta_gd\n",
    "test_final_gd = test_final_gd.merge(test_flux_skew_gd,on=['object_id'],how='left')\n",
    "GPU_RUN_TIME[step] = time.time() - start"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"train\"></a>\n",
    "## 4. Model training"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### train and validation with 5x speedup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = train_final_gd['target'].to_array()\n",
    "classes = sorted(np.unique(y))    \n",
    "# Taken from Giba's topic : https://www.kaggle.com/titericz\n",
    "# https://www.kaggle.com/c/PLAsTiCC-2018/discussion/67194\n",
    "# with Kyle Boone's post https://www.kaggle.com/kyleboone\n",
    "class_weights = {c: 1 for c in classes}\n",
    "class_weights.update({c:2 for c in [64, 15]})\n",
    "\n",
    "lbl = LabelEncoder()\n",
    "y = lbl.fit_transform(y)\n",
    "\n",
    "cols = [i for i in test_final_gd.columns if i not in ['object_id','target']]\n",
    "for col in cols:\n",
    "    train_final_gd[col] = train_final_gd[col].fillna(0).astype('float32')\n",
    "\n",
    "for col in cols:\n",
    "    test_final_gd[col] = test_final_gd[col].fillna(0).astype('float32')\n",
    "    \n",
    "X = train_final_gd[cols].as_matrix()\n",
    "Xt = test_final_gd[cols].as_matrix()\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(X,y,test_size=0.1,stratify=y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cpu_params = {\n",
    "            'objective': 'multi:softprob', \n",
    "            'tree_method': 'hist', \n",
    "            'nthread': 16, \n",
    "            'num_class':14,\n",
    "            'max_depth': 7, \n",
    "            'silent':1,\n",
    "            'subsample':0.7,\n",
    "            'colsample_bytree': 0.7,}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "func_loss = partial(xgb_multi_weighted_logloss, \n",
    "                        classes=classes, \n",
    "                        class_weights=class_weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# GPU\n",
    "gpu_params = cpu_params.copy()\n",
    "gpu_params.update({'objective': 'multi:softprob',\n",
    "                   'tree_method': 'gpu_hist', \n",
    "                  })"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "step = 'training'\n",
    "start = time.time()\n",
    "dtrain = xgb.DMatrix(data=X_train, label=y_train)\n",
    "dvalid = xgb.DMatrix(data=X_test, label=y_test)\n",
    "dtest = xgb.DMatrix(data=Xt)\n",
    "watchlist = [(dvalid, 'eval'), (dtrain, 'train')]\n",
    "clf = xgb.train(gpu_params, dtrain=dtrain,\n",
    "                num_boost_round=60,evals=watchlist,\n",
    "                feval=func_loss,early_stopping_rounds=10,\n",
    "                verbose_eval=1000)\n",
    "yp = clf.predict(dvalid)\n",
    "gpu_loss = multi_weighted_logloss(y_test, yp, classes, class_weights)\n",
    "ysub = clf.predict(dtest)\n",
    "line = 'validation loss %.4f'%gpu_loss\n",
    "print(colored(line,'green'))\n",
    "GPU_RUN_TIME[step] = time.time() - start"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"conclusions\"></a>\n",
    "## 5. Conclustions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GPU_RUN_TIME"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "steps = ['load data part1','ETL part1','load data part2','ETL part2','training']\n",
    "GPU_RUN_TIME['Overall'] = sum([GPU_RUN_TIME[i] for i in steps])\n",
    "steps.append('Overall')\n",
    "gpu_time = [GPU_RUN_TIME[i] for i in steps]\n",
    "df = pd.DataFrame({'GPU': gpu_time}, index=steps)\n",
    "df.plot.bar(rot=0,figsize=(20,5), fontsize=15, title='Running time: seconds')"
   ]
  }
 ],
 "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
}