Benchmarking of PR15049 "Faster manhattan_distances() for sparse matrices"

sparse_manhattan.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "sklearn version: 0.21.2\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "import numpy as np\n",
    "import scipy.sparse as ss\n",
    "\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "np.random.seed(0)\n",
    "\n",
    "N = 10\n",
    "\n",
    "X = ss.random(2000,3000)\n",
    "Y = ss.random(2000,3000)\n",
    "\n",
    "import sklearn\n",
    "print(\"sklearn version:\",sklearn.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython -f\n",
    "# distutils: extra_compile_args=-fopenmp\n",
    "# distutils: extra_link_args=-fopenmp\n",
    "#cython: boundscheck=False\n",
    "#cython: cdivision=True\n",
    "#cython: wraparound=False\n",
    "#\n",
    "# Author: Andreas Mueller <amueller@ais.uni-bonn.de>\n",
    "#         Lars Buitinck\n",
    "#         Paolo Toccaceli\n",
    "#\n",
    "# License: BSD 3 clause\n",
    "\n",
    "import numpy as np\n",
    "cimport numpy as np\n",
    "from cython cimport floating\n",
    "from cython.parallel cimport prange\n",
    "from libc.math cimport fabs\n",
    "\n",
    "np.import_array()\n",
    "\n",
    "def sparse_manhattan(floating[::1] X_data, int[:] X_indices, int[:] X_indptr,\n",
    "                      floating[::1] Y_data, int[:] Y_indices, int[:] Y_indptr,\n",
    "                      double[:, ::1] D, int num_threads):\n",
    "    \"\"\"Pairwise L1 distances for CSR matrices.\n",
    "\n",
    "    Usage:\n",
    "    >>> D = np.zeros(X.shape[0], Y.shape[0])\n",
    "    >>> _sparse_manhattan(X.data, X.indices, X.indptr,\n",
    "    ...                   Y.data, Y.indices, Y.indptr,\n",
    "    ...                   D)\n",
    "    \"\"\"\n",
    "    cdef np.npy_intp px, py, i, j, ix, iy\n",
    "    cdef double d = 0.0\n",
    "\n",
    "    cdef int m = D.shape[0]\n",
    "    cdef int n = D.shape[1]\n",
    "    \n",
    "    cdef int last_idx = 0\n",
    "\n",
    "    # We scan the matrices row by row.\n",
    "    # Given row px in X and row py in Y, we find the positions (i and j\n",
    "    # respectively), in .indices where the indices for the two rows start.\n",
    "    # If the indices (ix and iy) are the same, the corresponding data values\n",
    "    # are processed and the cursors i and j are advanced.\n",
    "    # If not, the lowest index is considered. Its associated data value is\n",
    "    # processed and its cursor is advanced.\n",
    "    # We proceed like this until one of the cursors hits the end for its row.\n",
    "    # Then we process all remaining data values in the other row.\n",
    "\n",
    "    # Below the avoidance of inplace operators is intentional.\n",
    "    # When prange is used, the inplace operator has a special meaning, i.e. it\n",
    "    # signals a \"reduction\"\n",
    "\n",
    "    for px in prange(m, nogil=True, num_threads=num_threads):\n",
    "        for py in range(n):\n",
    "            i = X_indptr[px]\n",
    "            j = Y_indptr[py]\n",
    "            d = 0.0\n",
    "            while i < X_indptr[px + 1] and j < Y_indptr[py + 1]:\n",
    "                ix = X_indices[i]\n",
    "                iy = Y_indices[j]\n",
    "\n",
    "                if ix == iy:\n",
    "                    d = d + fabs(X_data[i] - Y_data[j])\n",
    "                    i = i + 1\n",
    "                    j = j + 1\n",
    "                elif ix < iy:\n",
    "                    d = d + fabs(X_data[i])\n",
    "                    i = i + 1\n",
    "                else:\n",
    "                    d = d + fabs(Y_data[j])\n",
    "                    j = j + 1\n",
    "\n",
    "            if i == X_indptr[px + 1]:\n",
    "                last_idx = Y_indptr[py + 1]\n",
    "                while j < last_idx:\n",
    "                    d = d + fabs(Y_data[j])\n",
    "                    j = j + 1\n",
    "            else:\n",
    "                last_idx = X_indptr[px + 1]\n",
    "                while i < last_idx:\n",
    "                    d = d + fabs(X_data[i])\n",
    "                    i = i + 1\n",
    "\n",
    "            D[px, py] = d\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.sparse import csr_matrix\n",
    "\n",
    "def cython_manhattan_distances(X,Y,num_threads):\n",
    "    X = csr_matrix(X, copy=False)\n",
    "    Y = csr_matrix(Y, copy=False)\n",
    "    X.sum_duplicates()   # this also sorts indices in-place\n",
    "    Y.sum_duplicates()\n",
    "    D = np.zeros((X.shape[0], Y.shape[0]))\n",
    "    sparse_manhattan(X.data, X.indices, X.indptr,\n",
    "                      Y.data, Y.indices, Y.indptr,\n",
    "                      D,num_threads)\n",
    "    return D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics.pairwise import manhattan_distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 34.4 s, sys: 69.8 ms, total: 34.4 s\n",
      "Wall time: 34.5 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "for i in range(N):\n",
    "    D = manhattan_distances(X, Y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 11.8 s, sys: 12.1 ms, total: 11.8 s\n",
      "Wall time: 11.9 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "for i in range(N):\n",
    "    cython_manhattan_distances(X,Y,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 12 s, sys: 4.36 ms, total: 12 s\n",
      "Wall time: 6.04 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "for i in range(N):\n",
    "    cython_manhattan_distances(X,Y,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 14.4 s, sys: 3.64 ms, total: 14.4 s\n",
      "Wall time: 5.16 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "for i in range(N):\n",
    "    cython_manhattan_distances(X,Y,3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 16.3 s, sys: 40 ms, total: 16.4 s\n",
      "Wall time: 4.25 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "for i in range(N):\n",
    "    Dc = cython_manhattan_distances(X,Y,4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.all(np.isclose(D,Dc))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "model name\t: Intel(R) Core(TM) i5-7300U CPU @ 2.60GHz\n"
     ]
    }
   ],
   "source": [
    "! grep -m 1 \"model name\" /proc/cpuinfo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4\n"
     ]
    }
   ],
   "source": [
    "! grep -c \"core id\" /proc/cpuinfo"
   ]
  },
  {
   "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}