03-End_to_End_RAPIDS.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "12b84d8e-0f54-48c3-8a5f-c48161743f79",
   "metadata": {},
   "source": [
    "# Introduction to End-to-End RAPIDS Workflows"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f12fb5e1-f301-484f-abaf-acfcb2c7bfc7",
   "metadata": {},
   "source": [
    "This tutorial will teach developers how to build an end-to-end workflow with cuDF, cuML, and accelerated XGBoost. You will have the chance to ingest data, conduct ETL, perform EDA, train an XGBoost model, and use SHAP to gain insights into the predictions made by the model. \n",
    "\n",
    "\n",
    "We're going to be working with data from the [CitiBike data set](https://console.cloud.google.com/marketplace/product/city-of-new-york/nyc-citi-bike?pli=1&project=nv-ai-infra). CitiBike is a bike rental company which operates in NYC. Bikes are 'stored' at docking stations around the city, and users can rent a bike and return it to any docking station. We will use the historical information to attempt to predict the duration of a user's ride, given their starting station, as well as some other information. \n",
    "\n",
    "\n",
    "Before we begin, we're going to check what kind of GPU we have using [nvidia-smi](https://developer.nvidia.com/nvidia-system-management-interface). `nvidia-smi` has a whole range of functions described at the link. We are just going to use it to see general information about our GPU."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ebf90750-1bd8-4a09-b382-4bb276dd0e3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "!nvidia-smi"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f888a4c0-5535-4298-8c04-f43b2cce333f",
   "metadata": {},
   "source": [
    "Here we see that we have a [16 GB card](https://www.google.com/search?q=mib+to+gb&ei=3s07YsmeELHt9AP4z5zgDg&ved=0ahUKEwjJhbPQ0N32AhWxNn0KHfgnB-wQ4dUDCA4&uact=5&oq=mib+to+gb&gs_lcp=Cgdnd3Mtd2l6EAMyBwgAEEcQsAMyBwgAEEcQsAMyBwgAEEcQsAMyBwgAEEcQsAMyBwgAEEcQsAMyBwgAEEcQsAMyBwgAEEcQsAMyBwgAEEcQsAMyBwgAELADEEMyBwgAELADEEMyBwgAELADEEMyBwgAELADEEMyCggAEOQCELADGAEyCggAEOQCELADGAEyCggAEOQCELADGAEyDwguENQCEMgDELADEEMYAjIPCC4Q1AIQyAMQsAMQQxgCSgQIQRgASgQIRhgBUABYAGCABmgBcAF4AIABAIgBAJIBAJgBAMgBEcABAdoBBggBEAEYCdoBBggCEAEYCA&sclient=gws-wiz). If we had multiple cards, we would use `dask_cudf`. This will be covered in another notebook. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e378859f-24b0-4a09-a1c3-207d4b2a388d",
   "metadata": {},
   "source": [
    "## Importing the data\n",
    "\n",
    "Before we begin, we need to install a couple of packages.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a629b283-6316-4f5c-8e4e-458e3b1496b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install google-cloud-bigquery"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9e8a2cd-3648-47de-b139-372f2d5c7b3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install db_dtypes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d269a758-f92a-4983-8ba3-ed834e3cd0aa",
   "metadata": {},
   "source": [
    "The CitiBike data is available for download directly from an BigQuery. In the following cell, we import the data from 2014 only. \n",
    "\n",
    "\n",
    "_You can change the years and the number of years in the cell below by altering the `WHERE` statement._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6c87a6e-246d-43bc-a302-0a3b62870dcf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import time\n",
    "import cupy as cp\n",
    "import cudf \n",
    "from google.cloud import bigquery\n",
    "\n",
    "#os.environ.setdefault(\"GCLOUD_PROJECT\", \"hotornot-1078\")\n",
    "\n",
    "query = \"\"\"\n",
    "SELECT * \n",
    "FROM `bigquery-public-data.new_york_citibike.citibike_trips` \n",
    "WHERE EXTRACT(YEAR from starttime) = 2014\n",
    "\"\"\"\n",
    "client = bigquery.Client()\n",
    "job = client.query(query)\n",
    "pd_df = job.to_dataframe()\n",
    "df = cudf.from_pandas(pd_df)\n",
    "del(pd_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2fe90e95-bc1d-44ec-a0ba-ff4a1a5c5fd3",
   "metadata": {},
   "source": [
    "Let's look at the data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1113a504-2bc6-49a8-98a2-e4ac97381e9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a72167fa-7e9b-4393-83a2-d0148f01f0d5",
   "metadata": {},
   "source": [
    "Let's see the data types we have from the import."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ba3c8e1-1703-4f5f-a899-94b8dab073e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.dtypes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "834a78ad-b15f-4a85-91f3-033cc9d25364",
   "metadata": {},
   "source": [
    "Let's take a quick look at the data. We see that when `df.describe()` is given mixed types, we should tell it to include all the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc5dd7fc-0632-4626-b9d9-ee90358b81cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.describe(include='all')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73f2d6bd-5f2c-45c6-a082-94e12d43764d",
   "metadata": {},
   "source": [
    "## Data Cleaning and Feature Engineering"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46b157af-574d-42af-89b7-d048653f0238",
   "metadata": {},
   "source": [
    "The data contains some redundant information - `start_station_id` and `end_station_id` are both captured by the station names and latitude/longitude data. We drop this redundant information. \n",
    "\n",
    "We also remove all information about the end station. We wish to predict the duration of the user's ride at the point of pick up, and their bike drop-off destination would not be known to us at that time. \n",
    "\n",
    "We don't expect the `bike_id` to give us insight into ride duration so we remove that from the data set."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a655e7f0-0bc2-4343-aa06-bfe0e51ee0cf",
   "metadata": {},
   "source": [
    "We drop infromation based on `tripduration`, starting with observations where `tripduration` is negative - Bikes can do a lot of things, but they can't travel back in time!  \n",
    "\n",
    "We remove any trips lasting less than five minutes, as these are likely to indicate a malfunctioning bike which is quickly returned, rather than a real journey. \n",
    "\n",
    "We also drop all rides that lasted longer than 10 hours from our data  -  The citi bikes are supposed to be used for relatively short trips round the city, and are not suitable for long journeys. We don't want this data to skew our model.\n",
    "\n",
    "Finally, we drop all recorded rides that contain missing data for any of the remaining columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f036645b-93c8-4ccf-a9e7-70a28680e4b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = df.drop(['start_station_id', 'end_station_id', 'end_station_name', 'bikeid', 'stoptime', 'end_station_latitude', 'end_station_longitude'], axis=1) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e99ac38-bf0f-491f-b330-8d85fda30a79",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['tripduration'] = df['tripduration'].where(df['tripduration']>300)\n",
    "df['tripduration'] = df['tripduration'].where(df['tripduration']<=36000)\n",
    "df = df.dropna()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "36ae7a39-c3a4-4ec2-9c44-7e922de24d72",
   "metadata": {},
   "source": [
    "Next, we grab some things from the time fields that will be useful as features for our model. We're doing to create a variable grouping the time of day into one of six periods, the day of week, the month, and then we're going to drop those time variables. The exact second a bike was rented or returned likely has limited explanatory value. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f30c9bb-7787-4272-874d-e464492b7938",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['start_hour_of_the_day'] = df['starttime'].dt.hour\n",
    "df['dow'] = df['starttime'].dt.dayofweek\n",
    "df['month'] = df['starttime'].dt.month\n",
    "df = df.drop(['starttime'], axis=1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5c2b01c-71e8-4da1-854f-0c2bfa1d0d73",
   "metadata": {},
   "source": [
    "We're going to use cuML for the next bit of ETL to encode labels into numbers for our analysis. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ce325f5-3aef-4e66-a5dd-01ff3c3979ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "import cuml\n",
    "\n",
    "le = cuml.LabelEncoder()\n",
    "df['start_station_name'] = le.fit_transform(df['start_station_name'])\n",
    "df['usertype'] = le.fit_transform(df['usertype'])\n",
    "df['gender'] = le.fit_transform(df['gender'])\n",
    "df['customer_plan'] = le.fit_transform(df['customer_plan'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71317b31-dcd3-48cb-a15e-dbeac1d8bd65",
   "metadata": {},
   "source": [
    "Given that we are aiming to predict the length of the ride in seconds, it seems unfair to include both the hour at which the journey starts and the hour at which the journey stops in our feature vectors - let's remove this now, and see how well we can predict trip duration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cfe29b1-1f6f-40dc-8e44-9831a0c3e1f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c938040e-a0f8-4f37-8005-eabb97f5c3ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edd61d9e-7631-46c8-b792-c78840ab1451",
   "metadata": {},
   "source": [
    "Now that our data is cleaned up, it's time to see how well we can predict trip duration. We'll start by making a simple XGBoost model and then we will move onto an ensemble with some other methods with cuML."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48d34673-1e08-4dd4-973e-947a8f311c38",
   "metadata": {},
   "source": [
    "## XGBoost Prediction Model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26ed63a0-6efb-4316-9bcc-6bc758d45c2f",
   "metadata": {},
   "source": [
    "First, we want to split our data into train and test sets. We do this with cuML. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e726248f-b027-482a-9368-d0e3c4127987",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = cuml.train_test_split(df, 'tripduration', train_size=0.8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72428e6d-36fe-4c8c-abc0-d87737e188a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import xgboost as xgb\n",
    "\n",
    "dtrain = xgb.DMatrix(X_train, y_train)\n",
    "dtest = xgb.DMatrix(X_test, y_test) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36a1ea3d-9a23-41ee-9017-8864f053b956",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = {\n",
    "    'learning_rate': 0.01,\n",
    "    'max_depth': 5,\n",
    "    'objective': 'reg:squarederror',\n",
    "    'subsample': 0.8,\n",
    "    'disable_default_eval_metric':True, \n",
    "    'tree_method':'gpu_hist' \n",
    "}\n",
    "\n",
    "trained_model = xgb.train(\n",
    "    params,\n",
    "    dtrain,\n",
    "    num_boost_round=1000,\n",
    "    evals=[(dtrain, 'train')]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf0fdf0d-d20f-42c7-9d0f-2369202336db",
   "metadata": {},
   "source": [
    "We'll save the trained model, so that we can re-load it later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0654d51-4be7-4d6d-add6-b9f835646854",
   "metadata": {},
   "outputs": [],
   "source": [
    "trained_model.save_model(\"xgb.model\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7321699-beb4-4327-8475-bb929b591970",
   "metadata": {},
   "source": [
    "Now let's see how well our fitted model looks out in the wild."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5cb7c566-0ff6-4a26-ac16-5c567eefbb2e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "prediction = trained_model.predict(dtest).astype('int64')\n",
    "print(\"RMSE: {}\".format(cp.sqrt(cuml.metrics.mean_squared_error(y_test.values, prediction))))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87d1e350-00f1-489c-91a2-db4e6682472a",
   "metadata": {},
   "source": [
    "Looks like our model's predictions are off 13 minutes with our quick model - why not see if you can change the parameter values and improve the model's performance. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2b7bd16-9fcc-41b6-a7d7-34783a5fc4d1",
   "metadata": {},
   "source": [
    "## Model Explainability with SHAP\n",
    "\n",
    "When using complex models, such as XGBoost, it's not always straightforward to understand the predictions made by the model. In this section we use Shapley Additive Explanation (SHAP) values to gain insight into the Machine Learning model.\n",
    "\n",
    "Computing SHAP values is a computationally expensive procedure, but we accelerate the procedure by running on NVIDIA GPUs. To save more time, we compute SHAP values on a subset of our data.\n",
    "\n",
    "Much of the code in this section is taken from this great [blog](https://medium.com/rapids-ai/gpu-accelerated-shap-values-with-xgboost-1-3-and-rapids-587fad6822) on GPU-Accelerated SHAP Values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11622a9e-8da0-4152-b05a-7208b9ebdaf5",
   "metadata": {},
   "outputs": [],
   "source": [
    "shap_sample = xgb.DMatrix(X_test.sample(frac=0.01))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08291f46-95a6-4740-94b8-d459a37f304f",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "trained_model.set_param({\"predictor\": \"gpu_predictor\"})\n",
    "shap_values = trained_model.predict(shap_sample, pred_contribs=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed6a3333-f838-4c32-ad98-3c63a422b451",
   "metadata": {},
   "source": [
    "We can aggregate and visualse these SHAP values to see which of the features in our data had the most impact on the predictions made by our model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c78d7688-7e7c-4771-95e1-65c28f3d5c49",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def plot_feature_importance(feature_names, shap_values):\n",
    " # Get the mean absolute contribution for each feature\n",
    " aggregate = np.mean(np.abs(shap_values[:, 0:-1]), axis=0)\n",
    " # sort by magnitude\n",
    " z = [(x, y) for y, x in sorted(zip(aggregate, feature_names), reverse=True)]\n",
    " z = list(zip(*z))\n",
    " plt.bar(z[0], z[1])\n",
    " plt.xticks(rotation=90)\n",
    " plt.tight_layout()\n",
    " plt.show()\n",
    "\n",
    "\n",
    "plot_feature_importance(X_test.columns, shap_values)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22f578dc-c462-4989-8872-be17ed270199",
   "metadata": {},
   "source": [
    "This shows us that the most important features in predicting ride duration are the location of the pick up point. \n",
    "\n",
    "\n",
    "We can also use SHAP to consider the importance of interactions between features. This is more computationally expensive again, but can bring valuable insights. The following cell will take around 100 seconds to run. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef2fb4f0-89ba-4dd1-ac95-bff88e6d9e69",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "shap_interactions = trained_model.predict(shap_sample, pred_interactions=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cdd4e2a-d667-40c3-b8e5-0d2f08105e7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_top_k_interactions(feature_names, shap_interactions, k):\n",
    " # Get the mean absolute contribution for each feature interaction\n",
    " aggregate_interactions = np.mean(np.abs(shap_interactions[:, :-1, :-1]), axis=0)\n",
    " interactions = []\n",
    " for i in range(aggregate_interactions.shape[0]):\n",
    "     for j in range(aggregate_interactions.shape[1]):\n",
    "         if j < i:\n",
    "             interactions.append(\n",
    "             (feature_names[i] + \"-\" + feature_names[j], aggregate_interactions[i][j] * 2))\n",
    " # sort by magnitude\n",
    " interactions.sort(key=lambda x: x[1], reverse=True)\n",
    " interaction_features, interaction_values = map(tuple, zip(*interactions))\n",
    " plt.bar(interaction_features[:k], interaction_values[:k])\n",
    " plt.xticks(rotation=90)\n",
    " plt.tight_layout()\n",
    " plt.show()\n",
    "\n",
    "\n",
    "plot_top_k_interactions(X_test.columns, shap_interactions, 10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70363ba3-213c-431a-ba43-b50d3c760190",
   "metadata": {},
   "source": [
    "Here we see (unsurprisingly) that the interactions between the starting longitude and latitude greatly influence the predictions, followed by a location and starting time of ride. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64e0f44a-a5da-4f1b-b932-619df627d818",
   "metadata": {},
   "source": [
    "## Accelerating Inference \n",
    "\n",
    "Throughout this notebook we have run most of our computation on the GPU. In this Section, we compare the speed it takes to make predictions on a CPU vs the GPU. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a39eceeb-ad1d-435d-b42f-a6a629d3f665",
   "metadata": {},
   "outputs": [],
   "source": [
    "xgb_features = xgb.DMatrix(X_test.astype(\"float32\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6551d3bf-d3bf-4373-a3cd-4e2c2270af9a",
   "metadata": {},
   "source": [
    "### CPU\n",
    "\n",
    "We first re-load the model from file, as XGBoost caches the results of previous predictions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c920c434-5a91-4b7d-baa1-8e042abe147d",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "model = xgb.Booster(model_file=\"xgb.model\")\n",
    "model.set_param({\"predictor\": \"cpu_predictor\"})\n",
    "predictions = model.predict(xgb_features)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ffcdce86-9436-46cf-a5e8-e5f680c007a1",
   "metadata": {},
   "source": [
    "### GPU\n",
    "\n",
    "Now we can again reload the model, and this time run the same predictions on the GPU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cdd7cbc4-6e27-4678-99c8-a014d2a07267",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "model = xgb.Booster(model_file=\"xgb.model\")\n",
    "model.set_param({\"predictor\": \"gpu_predictor\"})\n",
    "predictions = model.predict(xgb_features)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca2fc191-57c6-4c13-88ca-09539fcca314",
   "metadata": {},
   "source": [
    "So you can see that the GPU allows us to make predictions in a fraction of the time taken on CPU. This is ideal for situations requiring real-time inference. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56b014a6-5a26-4e1a-92c1-40b84fe4bf65",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "In this notebook you've seen how we can use cuML, cuDF and XGBoost to explore and clean data, compute feature vectors and train a machine learning model to predict ride duration on the CitiBike Data Set. \n",
    "\n",
    "To find out more, check out [RAPIDS.ai](http://rapids.ai)."
   ]
  }
 ],
 "metadata": {
  "environment": {
   "kernel": "conda-env-rapids-22.04-py",
   "name": "common-cu110.m91",
   "type": "gcloud",
   "uri": "gcr.io/deeplearning-platform-release/base-cu110:m91"
  },
  "kernelspec": {
   "display_name": "Python [conda env:rapids-22.04]",
   "language": "python",
   "name": "conda-env-rapids-22.04-py"
  },
  "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}