I have been asking related questions on this stack for the past week to try to isolate things I didn't understand about using the @jit decorator with Numba in Python. However, I'm hitting the wall, so I'll just write the whole problem.
The issue at hand is to calculate the minimal distance between pairs of a large number of segments. The segments are represented by their beginning and endpoints in 3D. Mathematically, every segment is parameterized as [AB] = A + (B-A)*s, where s in [0,1], and A and B are the beginning and endpoints of the segment. For two such segments, the minimal distance can be computed and the formula is given here.
I have already exposed this problem on another thread, and the answer given dealt with replacing double loops of my code by vectorizing the problem, which however would hit memory issues for large sets of segments. Therefore, I've decided to stick with the loops, and use numba's jit instead.
Since the solution to the problem requires a lot of dot products, and numpy's dot product is not supported by numba, I started by implementing my own 3D dot product.
import numpy as np
from numba import jit, autojit, double, float64, float32, void, int32
def my_dot(a,b):
res = a[0]*b[0]+a[1]*b[1]+a[2]*b[2]
return res
dot_jit = jit(double(double[:], double[:]))(my_dot) #I know, it's not of much use here.
The function computing the minimal distance of all pairs in N segments takes as input an Nx6 array (6 coordinates)
def compute_stuff(array_to_compute):
N = len(array_to_compute)
con_mat = np.zeros((N,N))
for i in range(N):
for j in range(i+1,N):
p0 = array_to_compute[i,0:3]
p1 = array_to_compute[i,3:6]
q0 = array_to_compute[j,0:3]
q1 = array_to_compute[j,3:6]
s = ( dot_jit((p1-p0),(q1-q0))*dot_jit((q1-q0),(p0-q0)) - dot_jit((q1-q0),(q1-q0))*dot_jit((p1-p0),(p0-q0)))/( dot_jit((p1-p0),(p1-p0))*dot_jit((q1-q0),(q1-q0)) - dot_jit((p1-p0),(q1-q0))**2 )
t = ( dot_jit((p1-p0),(p1-p0))*dot_jit((q1-q0),(p0-q0)) -dot_jit((p1-p0),(q1-q0))*dot_jit((p1-p0),(p0-q0)))/( dot_jit((p1-p0),(p1-p0))*dot_jit((q1-q0),(q1-q0)) - dot_jit((p1-p0),(q1-q0))**2 )
con_mat[i,j] = np.sum( (p0+(p1-p0)*s-(q0+(q1-q0)*t))**2 )
return con_mat
fast_compute_stuff = jit(double[:,:](double[:,:]))(compute_stuff)
So, compute_stuff(arg) takes as argument an 2D np.array ( double[:,:]) , performs a bunch of numba-supported (?) operations, and returns another 2D np.array (double[:,:]). However,
v = np.random.random( (100,6) )
%timeit compute_stuff(v)
%timeit fast_compute_stuff(v)
I get 134 and 123 ms per loop. Can you shed some light on why I fail to speed up my function ? Any feedback would be much appreciated.
np.dot
using numba's JIT compiler.np.dot
is just a thin wrapper that calls BLAS*gemm/*gemv
functions, which are heavily optimized and often multithreaded. Your best bet is probably to make sure that numpy is linked against the fastest multithreaded BLAS library you can get your hands on (probably either Intel's MKL or OpenBLAS).np.dot
anyway, so shouldn't expect much performance benefit from JITing away the overhead from the nestedfor
loops.