numpy/3. Nearest Points Exercise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Nearest Points Exercise\n",
    "\n",
    "Tamás Gál (tamas.gal@fau.de)\n",
    "\n",
    "The latest version of this notebook is available at [https://github.com/Asterics2020-Obelics](https://github.com/Asterics2020-Obelics/School2017/tree/master/numpy)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Python version: 3.6.5 (default, Jun  1 2018, 14:48:24) \n",
      "[GCC 4.2.1 Compatible Apple LLVM 9.1.0 (clang-902.0.39.1)]\n",
      "NumPy version: 1.14.3\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import sys\n",
    "\n",
    "print(\"Python version: {0}\\n\"\n",
    "      \"NumPy version: {1}\"\n",
    "      .format(sys.version, np.__version__))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Given an array of points (in 3D), find the nearest point for each one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.52682391, 0.84736992, 0.68428599],\n",
       "       [0.00342933, 0.55981111, 0.1826532 ],\n",
       "       [0.60919571, 0.78728847, 0.12705444],\n",
       "       ...,\n",
       "       [0.05618889, 0.1495999 , 0.14557865],\n",
       "       [0.64582629, 0.58587899, 0.52161598],\n",
       "       [0.14804983, 0.49087426, 0.35089677]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 500\n",
    "n_dims = 3\n",
    "points = np.random.random((N, n_dims))\n",
    "points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.07678683, 0.87280216, 0.24566237],\n",
       "       [0.21943355, 0.63271636, 0.77344128],\n",
       "       [0.1216955 , 0.44489452, 0.08678564],\n",
       "       [0.59164763, 0.09284068, 0.53764511]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 4\n",
    "n_dims = 3\n",
    "points = np.random.random((N, n_dims))\n",
    "points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 0, 0, 1])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distances = np.sum((points.reshape(N, 1, n_dims) - points)**2, axis=2)\n",
    "distances[np.diag_indices_from(distances)] = np.inf\n",
    "np.argmin(distances, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Discussion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.07678683, 0.87280216, 0.24566237],\n",
       "       [0.21943355, 0.63271636, 0.77344128],\n",
       "       [0.1216955 , 0.44489452, 0.08678564],\n",
       "       [0.59164763, 0.09284068, 0.53764511]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4, 3)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "points.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[0.07678683, 0.87280216, 0.24566237]],\n",
       "\n",
       "       [[0.21943355, 0.63271636, 0.77344128]],\n",
       "\n",
       "       [[0.1216955 , 0.44489452, 0.08678564]],\n",
       "\n",
       "       [[0.59164763, 0.09284068, 0.53764511]]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reshaped_points = points.reshape(N, 1, n_dims)\n",
    "reshaped_points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4, 1, 3)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reshaped_points.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[ 0.        ,  0.        ,  0.        ],\n",
       "        [-0.14264672,  0.2400858 , -0.52777891],\n",
       "        [-0.04490867,  0.42790764,  0.15887673],\n",
       "        [-0.5148608 ,  0.77996148, -0.29198274]],\n",
       "\n",
       "       [[ 0.14264672, -0.2400858 ,  0.52777891],\n",
       "        [ 0.        ,  0.        ,  0.        ],\n",
       "        [ 0.09773805,  0.18782183,  0.68665564],\n",
       "        [-0.37221408,  0.53987568,  0.23579617]],\n",
       "\n",
       "       [[ 0.04490867, -0.42790764, -0.15887673],\n",
       "        [-0.09773805, -0.18782183, -0.68665564],\n",
       "        [ 0.        ,  0.        ,  0.        ],\n",
       "        [-0.46995213,  0.35205385, -0.45085947]],\n",
       "\n",
       "       [[ 0.5148608 , -0.77996148,  0.29198274],\n",
       "        [ 0.37221408, -0.53987568, -0.23579617],\n",
       "        [ 0.46995213, -0.35205385,  0.45085947],\n",
       "        [ 0.        ,  0.        ,  0.        ]]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "differences = reshaped_points - points\n",
    "differences"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4, 1, 3)  -  (4, 3)  =>  (4, 4, 3)\n"
     ]
    }
   ],
   "source": [
    "print(reshaped_points.shape, \" - \", points.shape, \" => \", differences.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.        , 0.35653986, 0.21036355, 0.95867548],\n",
       "       [0.35653986, 0.        , 0.51632574, 0.48560891],\n",
       "       [0.21036355, 0.51632574, 0.        , 0.54807118],\n",
       "       [0.95867548, 0.48560891, 0.54807118, 0.        ]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distances = np.sum(differences**2, axis=2)\n",
    "distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4, 4)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distances.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.        , 0.35653986, 0.21036355, 0.95867548],\n",
       "       [0.35653986, 0.        , 0.51632574, 0.48560891],\n",
       "       [0.21036355, 0.51632574, 0.        , 0.54807118],\n",
       "       [0.95867548, 0.48560891, 0.54807118, 0.        ]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# get rid of self-distances first ;)\n",
    "distances[np.diag_indices_from(distances)] = np.inf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[       inf, 0.35653986, 0.21036355, 0.95867548],\n",
       "       [0.35653986,        inf, 0.51632574, 0.48560891],\n",
       "       [0.21036355, 0.51632574,        inf, 0.54807118],\n",
       "       [0.95867548, 0.48560891, 0.54807118,        inf]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 0, 0, 1])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nearest_points = np.argmin(distances, axis=1)\n",
    "nearest_points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Comparing to `scikit-learn`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 0, 0, 1])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nearest_points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2, 0, 0, 1])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.neighbors import NearestNeighbors\n",
    "_, indices = NearestNeighbors().fit(points).kneighbors(points, 2)\n",
    "indices[:, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "152 µs ± 4.25 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit _, indices = NearestNeighbors().fit(points).kneighbors(points, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "35.9 µs ± 1.36 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "distances = np.sum((points.reshape(N, 1, n_dims) - points)**2, axis=2)\n",
    "distances[np.diag_indices_from(distances)] = np.inf\n",
    "np.argmin(distances, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Our implementation is faster just because of the function call overhead of scikit-learn. Try bigger `N`s and see how they actually compare to each other."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}