cudf/notebooks_Apply_Operations_in_cuDF.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Background\n",
    "\n",
    "The RAPIDS cuDF library is a GPU DataFrame manipulation library based on Apache Arrow that accelerates loading, filtering, and manipulation of data for model training by providing a familiar Pandas-like API. While the Pandas library consists of a variety of methods to cover most data science use cases, it does not always cover all of them. Because accelerated, customized data transformation have been found to be very valuable, cuDF provides two special methods that serve this particular purpose:  `apply_rows` and `apply_chunk` functions, which utilize the Numba library to accelerate the data transformation via GPU in parallel.\n",
    "\n",
    "In this tutorial, I am going to show a few examples of how to use them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Difference between `apply_rows` and `apply_chunks`\n",
    "\n",
    "`apply_rows` is a special case of `apply_chunks`, which processes each of the rows of the Dataframe independently in parallel. Under the hood, the `apply_rows` method will optimally divide the long columns into chunks, and assign chunks into different GPU blocks to compute. Here is one example using `apply_rows` to double the input array and also print out the GPU block/grid allocation information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "999000.0\n",
      "999000.0\n"
     ]
    }
   ],
   "source": [
    "import cudf\n",
    "import numpy as np\n",
    "from numba import cuda\n",
    " \n",
    "df = cudf.DataFrame()\n",
    "df['in1'] = np.arange(1000, dtype=np.float64)\n",
    " \n",
    "def kernel(in1, out):\n",
    "    for i, x in enumerate(in1):\n",
    "        print('tid:', cuda.threadIdx.x, 'bid:', cuda.blockIdx.x,\n",
    "              'array size:', in1.size, 'block threads:', cuda.blockDim.x)\n",
    "        out[i] = x * 2.0\n",
    " \n",
    "outdf = df.apply_rows(kernel,\n",
    "                      incols=['in1'],\n",
    "                      outcols=dict(out=np.float64),\n",
    "                      kwargs=dict())\n",
    " \n",
    "print(outdf['in1'].sum()*2.0)\n",
    "print(outdf['out'].sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the output.txt, we can see that the for-loop in the kernel function is unrolled by the compiler automatically. It uses 14 CUDA blocks. Each CUDA block uses 64 threads to do the computation. Each of the thread in the block most of the time deals one element in the input array and sometimes deals with two elements. The order of processing row element is not defined. \n",
    "\n",
    "We implement the same array double logic with the `apply_chunks` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "9900.0\n",
      "9900.0\n"
     ]
    }
   ],
   "source": [
    "import cudf\n",
    "import numpy as np\n",
    "from numba import cuda\n",
    " \n",
    " \n",
    "df = cudf.DataFrame()\n",
    "df['in1'] = np.arange(100, dtype=np.float64)\n",
    " \n",
    " \n",
    "def kernel(in1, out):\n",
    "    print('tid:', cuda.threadIdx.x, 'bid:', cuda.blockIdx.x,\n",
    "          'array size:', in1.size, 'block threads:', cuda.blockDim.x)\n",
    "    for i in range(cuda.threadIdx.x, in1.size, cuda.blockDim.x):\n",
    "        out[i] = in1[i] * 2.0\n",
    " \n",
    "outdf = df.apply_chunks(kernel,\n",
    "                        incols=['in1'],\n",
    "                        outcols=dict(out=np.float64),\n",
    "                        kwargs=dict(),\n",
    "                        chunks=16,\n",
    "                        tpb=8)\n",
    " \n",
    "print(outdf['in1'].sum()*2.0)\n",
    "print(outdf['out'].sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the output.txt, we can see `apply_chunks` has more control than the `apply_rows` method. It can specify how to divide the long array into chunks, map each of the array chunks to different GPU blocks to process (chunks argument), and assign the number of thread in the block (tpb argument). The for-loop is no longer automatically unrolled in the kernel function as `apply_rows` method, but instead stays as the for-loop for that GPU thread. Each kernel corresponds to each thread in one block, and each has full access to all the elements in that chunk of the array. In this example, the chunk size is 16, and it uniformly cuts the 100 elements into 7 chunks (except the last one), and assigns them to 7 blocks. Each block has 8 threads to process this length 16 subarray (or length 4 for the last block). \n",
    "\n",
    "## Performance benchmark compare\n",
    "\n",
    "Here we compare the benefits of using cuDF `apply_rows` vs pandas `apply` method with the following python code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cudf\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import time\n",
    " \n",
    " \n",
    "data_length = 1e6\n",
    "df = cudf.DataFrame()\n",
    "df['in1'] = np.arange(data_length, dtype=np.float64)\n",
    " \n",
    " \n",
    "def kernel(in1, out):\n",
    "    for i, x in enumerate(in1):\n",
    "        out[i] = x * 2.0\n",
    " \n",
    "start = time.time()\n",
    "df = df.apply_rows(kernel,\n",
    "                   incols=['in1'],\n",
    "                   outcols=dict(out=np.float64),\n",
    "                   kwargs=dict())\n",
    "end = time.time()\n",
    "print('cuDF time', end-start)\n",
    "assert(np.isclose(df['in1'].sum()*2.0, df['out'].sum()))\n",
    " \n",
    " \n",
    "df = pd.DataFrame()\n",
    "df['in1'] = np.arange(data_length, dtype=np.float64)\n",
    "start = time.time()\n",
    "df['out'] = df.in1.apply(lambda x: x*2)\n",
    "end = time.time()\n",
    "print('pandas time', end-start)\n",
    "assert(np.isclose(df['in1'].sum()*2.0, df['out'].sum()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We change the data_length from 1e4 to 1e7, here is the computation time spent in cuDF and pandas.\n",
    "\n",
    "| data length | 1e3         | 1e4         | 1e5         | 1e6         | 1e7         | 1e8         |\n",
    "| ----------- | ----------- | ----------- | ----------- | ----------- | ----------- | ----------- |\n",
    "| cuDF Time(s)| 0.1750      | 0.1840      | 0.1750      | 0.1720      | 0.1770      | 0.2490      |\n",
    "| Pandas Time(s)| 0.0006    | 0.0022      | 0.0180      | 0.2500      | 2.1300      | 21.400      |\n",
    "| Speed Up    | 0.003x      | 0.011x      | 0.103x      | 1.453x      | **12.034x**     | **85.944x**      |\n",
    "\n",
    "As we can see, cuDF has an overhead of launching GPU kernels (mostly the kernel compilation time), but the computation time remains relatively constant in this test due to the massive number of cores in the P100 card. Notice that the CPU computation scales linearly with the length of the array due to the series computation nature of the `apply` function. cuDF has the computation advantage once the array size exceeds one million. \n",
    "\n",
    "## Realistic application\n",
    "\n",
    "In the financial service industry, data scientists usually need to compute features from time series data. The most popular method to process the time series data is to compute moving average. The following example shows how to utilize `apply_chunks` to speed up moving average computation for a long array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "cuDF time 0.7422566413879395\n",
      "pandas time 13.539819955825806\n"
     ]
    }
   ],
   "source": [
    "import cudf\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from numba import cuda\n",
    "import time\n",
    " \n",
    "data_length = int(5e8)\n",
    "average_window = 4\n",
    "df = cudf.DataFrame()\n",
    "threads_per_block = 128\n",
    "trunk_size = 10240\n",
    "df['in1'] = np.arange(data_length, dtype=np.float64)\n",
    " \n",
    " \n",
    "def kernel1(in1, out, average_length):\n",
    "    for i in range(cuda.threadIdx.x,\n",
    "               average_length-1, cuda.blockDim.x):\n",
    "        out[i] = np.inf\n",
    "    for i in range(cuda.threadIdx.x + average_length - 1,\n",
    "                   in1.size, cuda.blockDim.x):\n",
    "        summ = 0.0\n",
    "        for j in range(i - average_length + 1,\n",
    "                       i + 1):\n",
    "            summ += in1[j]\n",
    "        out[i] = summ / np.float64(average_length)\n",
    " \n",
    "def kernel2(in1, out, average_length):\n",
    "    if in1.size - average_length + cuda.threadIdx.x - average_length + 1 < 0 :\n",
    "        return\n",
    "    for i in range(in1.size - average_length + cuda.threadIdx.x,\n",
    "                   in1.size, cuda.blockDim.x):\n",
    "        summ = 0.0\n",
    "        for j in range(i - average_length + 1,\n",
    "                       i + 1):\n",
    "            #print(i,j, in1.size)\n",
    "            summ += in1[j]\n",
    "        out[i] = summ / np.float64(average_length)\n",
    " \n",
    " \n",
    "start = time.time()\n",
    "df = df.apply_chunks(kernel1,\n",
    "                     incols=['in1'],\n",
    "                     outcols=dict(out=np.float64),\n",
    "                     kwargs=dict(average_length=average_window),\n",
    "                     chunks=list(range(0, data_length,\n",
    "                                       trunk_size))+ [data_length],\n",
    "                     tpb=threads_per_block)\n",
    " \n",
    "df = df.apply_chunks(kernel2,\n",
    "                     incols=['in1', 'out'],\n",
    "                     outcols=dict(),\n",
    "                     kwargs=dict(average_length=average_window),\n",
    "                     chunks=[0]+list(range(average_window, data_length,\n",
    "                                           trunk_size))+ [data_length],\n",
    "                     tpb=threads_per_block)\n",
    "end = time.time()\n",
    "print('cuDF time', end-start)\n",
    " \n",
    "pdf = pd.DataFrame()\n",
    "pdf['in1'] = np.arange(data_length, dtype=np.float64)\n",
    "start = time.time()\n",
    "pdf['out'] = pdf.rolling(average_window).mean()\n",
    "end = time.time()\n",
    "print('pandas time', end-start)\n",
    " \n",
    "assert(np.isclose(pdf.out.as_matrix()[average_window:].mean(),\n",
    "       df.out.to_array()[average_window:].mean()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above code, we divide the array into subarrays of size \"trunk_size\", and send those subarrays to GPU blocks to compute moving average. However, there is no history for the elements at the beginning of the subarray. To fix this, we shift the chunk division by an offset of \"average_window\". Then we call the kernel2 to compute the moving average of those missing records only.. Note, in kernel2, we didn't define outcols as it will create a new GPU memory buffer and overwrite the old \"out\" values. Instead, we reuse the out array as input. For an array of 5e8 length, cuDF uses 0.74s to do the computation while pandas use 13.54s. \n",
    "\n",
    "This code is not optimized for performance. There are a few things we can do to make it faster. First, we can use shared memory to load the array and reduce the IO when doing the summation. Secondly, there is a lot of redundant summation done by the threads. We can maintain an accumulation summation array to help reduce the redundancy. This is outside the scope of this tutorial."
   ]
  },
  {
   "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": 2
}