Skip to content

Instantly share code, notes, and snippets.

@epifanio
Created September 24, 2018 03:16
Show Gist options
  • Save epifanio/47340242f5be2de2c50577bf82c37143 to your computer and use it in GitHub Desktop.
Save epifanio/47340242f5be2de2c50577bf82c37143 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"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
}
@epifanio
Copy link
Author

epifanio commented Sep 24, 2018

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment