numba_parallel.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "!rm -rf nbspatial_src2.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "!rm -rf nbspatial.cpython*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Writing nbspatial_src2.py\n"
     ]
    }
   ],
   "source": [
    "%%file nbspatial_src2.py\n",
    "from numba import jit, njit\n",
    "from numba.pycc import CC\n",
    "import numba\n",
    "cc = CC('nbspatial')\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "@cc.export('array_tracing', 'b1[:](f8[:,:], f8[:,:])')\n",
    "@njit(parallel=True)\n",
    "def array_tracing(xy, poly):\n",
    "    D = np.empty(len(xy), dtype=numba.boolean)\n",
    "    n = len(poly)\n",
    "    for i in numba.prange(1, len(D) - 1):\n",
    "        inside = False\n",
    "        p2x = 0.0\n",
    "        p2y = 0.0\n",
    "        xints = 0.0\n",
    "        p1x,p1y = poly[0]\n",
    "        x = xy[i][0]\n",
    "        y = xy[i][1]\n",
    "        for i in range(n+1):\n",
    "            p2x,p2y = poly[i % n]\n",
    "            if y > min(p1y,p2y):\n",
    "                if y <= max(p1y,p2y):\n",
    "                    if x <= max(p1x,p2x):\n",
    "                        if p1y != p2y:\n",
    "                            xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x\n",
    "                        if p1x == p2x or x <= xints:\n",
    "                            inside = not inside\n",
    "            p1x,p1y = p2x,p2y\n",
    "        D[i] = inside\n",
    "    return D\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    cc.compile()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "In file included from /usr/local/lib/python3.7/dist-packages/numba/pycc/modulemixin.c:17:0:\n",
      "/usr/local/lib/python3.7/dist-packages/numba/pycc/../_dynfunc.c: In function ‘dup_string’:\n",
      "/usr/local/lib/python3.7/dist-packages/numba/pycc/../_dynfunc.c:238:9: warning: assignment discards ‘const’ qualifier from pointer target type [-Wdiscarded-qualifiers]\n",
      "     tmp = PyString_AsString(strobj);\n",
      "         ^\n",
      "\u001b[0m"
     ]
    }
   ],
   "source": [
    "!python3.7 nbspatial_src2.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "# regular polygon for testing\n",
    "lenpoly = 10000\n",
    "polygon = np.array([[np.sin(x)+0.5,np.cos(x)+0.5] for x in np.linspace(0,2*np.pi,lenpoly)[:-1]])\n",
    "\n",
    "# random points set of points to test \n",
    "N = 1000000\n",
    "points = np.array([np.random.random(N), np.random.random(N)]).reshape(N,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nbspatial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 31s, sys: 512 ms, total: 1min 32s\n",
      "Wall time: 1min 31s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ True,  True,  True, ..., False, False, False])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "nbspatial.array_tracing(points, polygon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numba\n",
    "from numba import njit, jit\n",
    "\n",
    "@njit(parallel=True)\n",
    "def array_tracing2(xy, poly):\n",
    "    D = np.empty(len(xy), dtype=numba.boolean)\n",
    "    n = len(poly)\n",
    "    for i in numba.prange(1, len(D) - 1):\n",
    "        inside = False\n",
    "        p2x = 0.0\n",
    "        p2y = 0.0\n",
    "        xints = 0.0\n",
    "        p1x,p1y = poly[0]\n",
    "        x = xy[i][0]\n",
    "        y = xy[i][1]\n",
    "        for i in range(n+1):\n",
    "            p2x,p2y = poly[i % n]\n",
    "            if y > min(p1y,p2y):\n",
    "                if y <= max(p1y,p2y):\n",
    "                    if x <= max(p1x,p2x):\n",
    "                        if p1y != p2y:\n",
    "                            xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x\n",
    "                        if p1x == p2x or x <= xints:\n",
    "                            inside = not inside\n",
    "            p1x,p1y = p2x,p2y\n",
    "        D[i] = inside\n",
    "    return D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1min 45s, sys: 693 ms, total: 1min 46s\n",
      "Wall time: 9.48 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([ True,  True,  True, ..., False, False, False])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "array_tracing2(points, polygon) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

by adding a new module with the following function, which operates on a single point instead of an array:

%%file numba_parallel.py
import numba
from numba import jit
import numpy as np

@cc.export('ray_tracing1',  'b1(f8, f8, f8[:,:])')
@jit(nopython=True, nogil=True)
def ray_tracing1(x,y,poly):
    n = len(poly)
    inside = False
    p2x = 0.0
    p2y = 0.0
    xints = 0.0
    p1x,p1y = poly[0]
    for i in range(n+1):
        p2x,p2y = poly[i % n]
        if y > min(p1y,p2y):
            if y <= max(p1y,p2y):
                if x <= max(p1x,p2x):
                    if p1y != p2y:
                        xints = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
                    if p1x == p2x or x <= xints:
                        inside = not inside
        p1x,p1y = p2x,p2y

    return inside

I can define then a new function that make use of numba.prange to loop over all the points in parallel:

import numba

@numba.njit(parallel=True)
def parallel_tracing(pp, poly):
    D = np.empty(len(pp), dtype=numba.boolean)
    for i in numba.prange(1, len(D) - 1):
        D[i] = nbspatial.ray_tracing1(pp[i][0], pp[i][1], polygon)
    return D