In [1]:
%load_ext cythonmagic

import numpy as np
X = np.random.random((1000, 3))
D = np.empty((1000, 1000))
In [2]:
# Pure python version

def pairwise_python(X, D):
    M = X.shape[0]
    N = X.shape[1]
    for i in range(M):
        for j in range(M):
            d = 0.0
            for k in range(N):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = np.sqrt(d)
In [3]:
# numba version

import numpy as np
from numba import double
from numba.decorators import jit

@jit(arg_types=[double[:,:], double[:,:]])
def pairwise_numba(X, D):
    M = X.shape[0]
    N = X.shape[1]
    for i in range(M):
        for j in range(M):
            d = 0.0
            for k in range(N):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = np.sqrt(d)
CONVERTING double[:, :]
CONVERTING double
['d']
CONVERTING double[:, :]
CONVERTING double
['d']
str_to_llvmtype(): str = 'f64'
{'blocks': {0: <llvm.core.BasicBlock object at 0x2e6efd0>,
            29: <llvm.core.BasicBlock object at 0x2e7e050>,
            39: <llvm.core.BasicBlock object at 0x2e7e090>,
            42: <llvm.core.BasicBlock object at 0x2e7e0d0>,
            45: <llvm.core.BasicBlock object at 0x2e7e150>,
            48: <llvm.core.BasicBlock object at 0x2e7e190>,
            58: <llvm.core.BasicBlock object at 0x2e7e1d0>,
            61: <llvm.core.BasicBlock object at 0x2e7e210>,
            64: <llvm.core.BasicBlock object at 0x2e7e290>,
            73: <llvm.core.BasicBlock object at 0x2e7e2d0>,
            83: <llvm.core.BasicBlock object at 0x2e7e310>,
            86: <llvm.core.BasicBlock object at 0x2e7e350>,
            89: <llvm.core.BasicBlock object at 0x2e7e3d0>,
            136: <llvm.core.BasicBlock object at 0x2e7e390>,
            165: <llvm.core.BasicBlock object at 0x2e7e250>,
            169: <llvm.core.BasicBlock object at 0x2e7e110>},
 'blocks_dom': {0: set([0]),
                29: set([0, 29]),
                39: set([0, 29, 39, 42, 45, 48, 61, 165]),
                42: set([0, 29, 42]),
                45: set([0, 29, 42, 45]),
                48: set([0, 29, 42, 45, 48]),
                58: set([0, 29, 42, 45, 48, 58, 61, 64, 73, 86, 136]),
                61: set([0, 29, 42, 45, 48, 61]),
                64: set([0, 29, 42, 45, 48, 61, 64]),
                73: set([0, 29, 42, 45, 48, 61, 64, 73]),
                83: set([0, 29, 42, 45, 48, 61, 64, 73, 83, 86, 89]),
                86: set([0, 29, 42, 45, 48, 61, 64, 73, 86]),
                89: set([0, 29, 42, 45, 48, 61, 64, 73, 86, 89]),
                136: set([0, 29, 42, 45, 48, 61, 64, 73, 86, 136]),
                165: set([0, 29, 42, 45, 48, 61, 165]),
                169: set([0, 29, 42, 169])},
 'blocks_in': {0: set(),
               29: set([0]),
               39: set([165]),
               42: set([29, 39]),
               45: set([42]),
               48: set([45]),
               58: set([136]),
               61: set([48, 58]),
               64: set([61]),
               73: set([64]),
               83: set([89]),
               86: set([73, 83]),
               89: set([86]),
               136: set([86]),
               165: set([61]),
               169: set([42])},
 'blocks_out': {0: set([29]),
                29: set([42]),
                39: set([42]),
                42: set([45, 169]),
                45: set([48]),
                48: set([61]),
                58: set([61]),
                61: set([64, 165]),
                64: set([73]),
                73: set([86]),
                83: set([86]),
                86: set([89, 136]),
                89: set([83]),
                136: set([58]),
                165: set([39]),
                169: set()},
 'blocks_reaching': {0: set([0]),
                     29: set([0, 29]),
                     39: set([0,
                              29,
                              39,
                              42,
                              45,
                              48,
                              58,
                              61,
                              64,
                              73,
                              83,
                              86,
                              89,
                              136,
                              165]),
                     42: set([0,
                              29,
                              39,
                              42,
                              45,
                              48,
                              58,
                              61,
                              64,
                              73,
                              83,
                              86,
                              89,
                              136,
                              165]),
                     45: set([0,
                              29,
                              39,
                              42,
                              45,
                              48,
                              58,
                              61,
                              64,
                              73,
                              83,
                              86,
                              89,
                              136,
                              165]),
                     48: set([0,
                              29,
                              39,
                              42,
                              45,
                              48,
                              58,
                              61,
                              64,
                              73,
                              83,
                              86,
                              89,
                              136,
                              165]),
                     58: set([0,
                              29,
                              39,
                              42,
                              45,
                              48,
                              58,
                              61,
                              64,
                              73,
                              83,
                              86,
                              89,
                              136,
                              165]),
                     61: set([0,
                              29,
                              39,
                              42,
                              45,
                              48,
                              58,
                              61,
                              64,
                              73,
                              83,
                              86,
                              89,
                              136,
                              165]),
                     64: set([0,
                              29,
                              39,
                              42,
                              45,
                              48,
                              58,
                              61,
                              64,
                              73,
                              83,
                              86,
                              89,
                              136,
                              165]),
                     73: set([0,
                              29,
                              39,
                              42,
                              45,
                              48,
                              58,
                              61,
                              64,
                              73,
                              83,
                              86,
                              89,
                              136,
                              165]),
                     83: set([0,
                              29,
                              39,
                              42,
                              45,
                              48,
                              58,
                              61,
                              64,
                              73,
                              83,
                              86,
                              89,
                              136,
                              165]),
                     86: set([0,
                              29,
                              39,
                              42,
                              45,
                              48,
                              58,
                              61,
                              64,
                              73,
                              83,
                              86,
                              89,
                              136,
                              165]),
                     89: set([0,
                              29,
                              39,
                              42,
                              45,
                              48,
                              58,
                              61,
                              64,
                              73,
                              83,
                              86,
                              89,
                              136,
                              165]),
                     136: set([0,
                               29,
                               39,
                               42,
                               45,
                               48,
                               58,
                               61,
                               64,
                               73,
                               83,
                               86,
                               89,
                               136,
                               165]),
                     165: set([0,
                               29,
                               39,
                               42,
                               45,
                               48,
                               58,
                               61,
                               64,
                               73,
                               83,
                               86,
                               89,
                               136,
                               165]),
                     169: set([0,
                               29,
                               39,
                               42,
                               45,
                               48,
                               58,
                               61,
                               64,
                               73,
                               83,
                               86,
                               89,
                               136,
                               165,
                               169])},
 'blocks_reads': {0: set([0]),
                  29: set([2]),
                  39: set(),
                  42: set([4]),
                  45: set(),
                  48: set([2]),
                  58: set(),
                  61: set([5]),
                  64: set(),
                  73: set([3]),
                  83: set(),
                  86: set([7]),
                  89: set([0, 4, 5, 6, 7, 8]),
                  136: set([1, 4, 5, 6]),
                  165: set(),
                  169: set()},
 'blocks_writer': {0: {2: 10, 3: 23},
                   29: {4: 38},
                   39: {4: 39},
                   42: {4: 42},
                   45: {},
                   48: {5: 57},
                   58: {5: 58},
                   61: {5: 61},
                   64: {6: 67},
                   73: {7: 82},
                   83: {7: 83},
                   86: {6: 86, 7: 86},
                   89: {6: 130, 8: 116},
                   136: {},
                   165: {},
                   169: {}},
 'blocks_writes': {0: set([0, 1, 2, 3]),
                   29: set([4]),
                   39: set([4]),
                   42: set([4]),
                   45: set(),
                   48: set([5]),
                   58: set([5]),
                   61: set([5]),
                   64: set([6]),
                   73: set([7]),
                   83: set([7]),
                   86: set([6, 7]),
                   89: set([6, 8]),
                   136: set(),
                   165: set(),
                   169: set()},
 'translator': <numba.translate.Translate object at 0x2e6ed90>}
op_LOAD_ATTR(): 3 106 shape <Variable(val=<llvm.core.Argument object at 0x2e6eed0>, _llvm=<llvm.core.Argument object at 0x2e6eed0>, typ='arr[f64]')> arr[f64]
op_LOAD_ATTR(): { i64, i32*, i8*, i32, i64*, i64*, i8*, i8*, i32, i8*, i8*, i8*, i64* }*
op_LOAD_ATTR(): 16 106 shape <Variable(val=<llvm.core.Argument object at 0x2e6eed0>, _llvm=<llvm.core.Argument object at 0x2e6eed0>, typ='arr[f64]')> arr[f64]
op_LOAD_ATTR(): { i64, i32*, i8*, i32, i64*, i64*, i8*, i8*, i32, i8*, i8*, i8*, i64* }*
('op_CALL_FUNCTION():', <Variable(val=<built-in function range>, _llvm=None, typ=['func'])>)
str_to_llvmtype(): str = 'i64'
add_phi_incomming(): reaching_defs = {29: {0: 0, 1: 0, 2: 0, 3: 0, 4: 29},
 39: {0: 0, 1: 0, 2: 0, 3: 0, 4: 39, 5: 61}}
    crnt_block=42, pred=39, local=4
op_BINARY_ADD(): <Variable(val=<llvm.core.PHINode object at 0x2e7ea90>, _llvm=<llvm.core.PHINode object at 0x2e7ea90>, typ='i64')> + <Variable(val=<llvm.core.ConstantInt object at 0x2e7e510>, _llvm=<llvm.core.ConstantInt object at 0x2e7e510>, typ='i64')>
resolve_type(): arg1 = <Variable(val=<llvm.core.PHINode object at 0x2e7ea90>, _llvm=<llvm.core.PHINode object at 0x2e7ea90>, typ='i64')>, arg2 = <Variable(val=<llvm.core.ConstantInt object at 0x2e7e510>, _llvm=<llvm.core.ConstantInt object at 0x2e7e510>, typ='i64')>
resolve_type() ==> 'i64'
resolve_type(): arg1 = <Variable(val=<llvm.core.PHINode object at 0x2e7ea90>, _llvm=<llvm.core.PHINode object at 0x2e7ea90>, typ='i64')>, arg2 = <Variable(val=<llvm.core.Instruction object at 0x2e7e490>, _llvm=<llvm.core.Instruction object at 0x2e7e490>, typ='i64')>
resolve_type() ==> 'i64'
('op_CALL_FUNCTION():', <Variable(val=<built-in function range>, _llvm=None, typ=['func'])>)
str_to_llvmtype(): str = 'i64'
add_phi_incomming(): reaching_defs = {48: {0: 0, 1: 0, 2: 0, 3: 0, 4: 42, 5: 48},
 58: {0: 0, 1: 0, 2: 0, 3: 0, 4: 42, 5: 58, 6: 86, 7: 86}}
    crnt_block=61, pred=58, local=5
op_BINARY_ADD(): <Variable(val=<llvm.core.PHINode object at 0x2e7e4d0>, _llvm=<llvm.core.PHINode object at 0x2e7e4d0>, typ='i64')> + <Variable(val=<llvm.core.ConstantInt object at 0x2e7ea50>, _llvm=<llvm.core.ConstantInt object at 0x2e7ea50>, typ='i64')>
resolve_type(): arg1 = <Variable(val=<llvm.core.PHINode object at 0x2e7e4d0>, _llvm=<llvm.core.PHINode object at 0x2e7e4d0>, typ='i64')>, arg2 = <Variable(val=<llvm.core.ConstantInt object at 0x2e7ea50>, _llvm=<llvm.core.ConstantInt object at 0x2e7ea50>, typ='i64')>
resolve_type() ==> 'i64'
resolve_type(): arg1 = <Variable(val=<llvm.core.PHINode object at 0x2e7e4d0>, _llvm=<llvm.core.PHINode object at 0x2e7e4d0>, typ='i64')>, arg2 = <Variable(val=<llvm.core.Instruction object at 0x2e7e490>, _llvm=<llvm.core.Instruction object at 0x2e7e490>, typ='i64')>
resolve_type() ==> 'i64'
('op_CALL_FUNCTION():', <Variable(val=<built-in function range>, _llvm=None, typ=['func'])>)
str_to_llvmtype(): str = 'f64'
add_phi_incomming(): reaching_defs = {73: {0: 0, 1: 0, 2: 0, 3: 0, 4: 42, 5: 61, 6: 64, 7: 73},
 83: {0: 0, 1: 0, 2: 0, 3: 0, 4: 42, 5: 61, 6: 89, 7: 83, 8: 89}}
    crnt_block=86, pred=83, local=6
str_to_llvmtype(): str = 'i64'
add_phi_incomming(): reaching_defs = {73: {0: 0, 1: 0, 2: 0, 3: 0, 4: 42, 5: 61, 6: 64, 7: 73},
 83: {0: 0, 1: 0, 2: 0, 3: 0, 4: 42, 5: 61, 6: 89, 7: 83, 8: 89}}
    crnt_block=86, pred=83, local=7
op_BINARY_ADD(): <Variable(val=<llvm.core.PHINode object at 0x2e7e590>, _llvm=<llvm.core.PHINode object at 0x2e7e590>, typ='i64')> + <Variable(val=<llvm.core.ConstantInt object at 0x2e7eb90>, _llvm=<llvm.core.ConstantInt object at 0x2e7eb90>, typ='i64')>
resolve_type(): arg1 = <Variable(val=<llvm.core.PHINode object at 0x2e7e590>, _llvm=<llvm.core.PHINode object at 0x2e7e590>, typ='i64')>, arg2 = <Variable(val=<llvm.core.ConstantInt object at 0x2e7eb90>, _llvm=<llvm.core.ConstantInt object at 0x2e7eb90>, typ='i64')>
resolve_type() ==> 'i64'
resolve_type(): arg1 = <Variable(val=<llvm.core.PHINode object at 0x2e7e590>, _llvm=<llvm.core.PHINode object at 0x2e7e590>, typ='i64')>, arg2 = <Variable(val=<llvm.core.Instruction object at 0x2e7e410>, _llvm=<llvm.core.Instruction object at 0x2e7e410>, typ='i64')>
resolve_type() ==> 'i64'
op_BINARY_SUBSCR(): arr_var.typ = arr[f64]
str_to_llvmtype(): str = 'f64'
  %28 = load double* %27
op_BINARY_SUBSCR(): arr_var.typ = arr[f64]
str_to_llvmtype(): str = 'f64'
  %39 = load double* %38
resolve_type(): arg1 = <Variable(val=<llvm.core.Instruction object at 0x2e7ef90>, _llvm=<llvm.core.Instruction object at 0x2e7ef90>, typ='f64')>, arg2 = <Variable(val=<llvm.core.Instruction object at 0x2e7edd0>, _llvm=<llvm.core.Instruction object at 0x2e7edd0>, typ='f64')>
resolve_type() ==> 'f64'
resolve_type(): arg1 = <Variable(val=<llvm.core.Instruction object at 0x2e7ed10>, _llvm=<llvm.core.Instruction object at 0x2e7ed10>, typ='f64')>, arg2 = <Variable(val=<llvm.core.Instruction object at 0x2e7ed10>, _llvm=<llvm.core.Instruction object at 0x2e7ed10>, typ='f64')>
resolve_type() ==> 'f64'
op_BINARY_ADD(): <Variable(val=<llvm.core.PHINode object at 0x2e7eb10>, _llvm=<llvm.core.PHINode object at 0x2e7eb10>, typ='f64')> + <Variable(val=<llvm.core.Instruction object at 0x2e7ef10>, _llvm=<llvm.core.Instruction object at 0x2e7ef10>, typ='f64')>
resolve_type(): arg1 = <Variable(val=<llvm.core.PHINode object at 0x2e7eb10>, _llvm=<llvm.core.PHINode object at 0x2e7eb10>, typ='f64')>, arg2 = <Variable(val=<llvm.core.Instruction object at 0x2e7ef10>, _llvm=<llvm.core.Instruction object at 0x2e7ef10>, typ='f64')>
resolve_type() ==> 'f64'
op_LOAD_ATTR(): 140 106 sqrt <Variable(val=<module 'numpy' from '/home/vanderplas/PythonEnv/pydev/local/lib/python2.7/site-packages/numpy/__init__.pyc'>, _llvm=None, typ=None)> None
('op_CALL_FUNCTION():', <Variable(val=<ufunc 'sqrt'>, _llvm=None, typ=None)>)
str_to_llvmtype(): str = 'f64'
op_STORE_SUBSCR(): 161 60 None
op_STORE_SUBSCR(): <Variable(val=<llvm.core.Argument object at 0x2e6ef90>, _llvm=<llvm.core.Argument object at 0x2e6ef90>, typ='arr[f64]')>[<Variable(val=(<Variable(val=<llvm.core.PHINode object at 0x2e7ea90>, _llvm=<llvm.core.PHINode object at 0x2e7ea90>, typ='i64')>, <Variable(val=<llvm.core.PHINode object at 0x2e7e4d0>, _llvm=<llvm.core.PHINode object at 0x2e7e4d0>, typ='i64')>), _llvm=None, typ='tuple')>] = <Variable(val=<llvm.core.CallOrInvokeInstruction object at 0x2e7ea50>, _llvm=<llvm.core.CallOrInvokeInstruction object at 0x2e7ea50>, typ='f64')>
op_STORE_SUBSCR(): arr_lval = '{ i64, i32*, i8*, i32, i64*, i64*, i8*, i8*, i32, i8*, i8*, i8*, i64* }* %D', arr_ltype = '{ i64, i32*, i8*, i32, i64*, i64*, i8*, i8*, i32, i8*, i8*, i8*, i64* }*'
str_to_llvmtype(): str = 'f64'
; ModuleID = 'pairwise_numba_mod_2e6ed90'

define double @pairwise_numba({ i64, i32*, i8*, i32, i64*, i64*, i8*, i8*, i32, i8*, i8*, i8*, i64* }* %X, { i64, i32*, i8*, i32, i64*, i64*, i8*, i8*, i32, i8*, i8*, i8*, i64* }* %D) {
Entry:
  %0 = getelementptr { i64, i32*, i8*, i32, i64*, i64*, i8*, i8*, i32, i8*, i8*, i8*, i64* }* %X, i32 0, i32 4
  %1 = load i64** %0
  %2 = getelementptr i64* %1, i32 0
  %3 = load i64* %2
  %4 = getelementptr { i64, i32*, i8*, i32, i64*, i64*, i8*, i8*, i32, i8*, i8*, i8*, i64* }* %X, i32 0, i32 4
  %5 = load i64** %4
  %6 = getelementptr i64* %5, i32 1
  %7 = load i64* %6
  br label %BLOCK_29

BLOCK_29:                                         ; preds = %Entry
  br label %BLOCK_42

BLOCK_39:                                         ; preds = %BLOCK_165
  %8 = add i64 %9, 1
  br label %BLOCK_42

BLOCK_42:                                         ; preds = %BLOCK_39, %BLOCK_29
  %9 = phi i64 [ 0, %BLOCK_29 ], [ %8, %BLOCK_39 ]
  %10 = icmp slt i64 %9, %3
  br i1 %10, label %BLOCK_45, label %BLOCK_169

BLOCK_169:                                        ; preds = %BLOCK_42
  ret double 0.000000e+00

BLOCK_45:                                         ; preds = %BLOCK_42
  br label %BLOCK_48

BLOCK_48:                                         ; preds = %BLOCK_45
  br label %BLOCK_61

BLOCK_58:                                         ; preds = %BLOCK_136
  %11 = add i64 %12, 1
  br label %BLOCK_61

BLOCK_61:                                         ; preds = %BLOCK_58, %BLOCK_48
  %12 = phi i64 [ 0, %BLOCK_48 ], [ %11, %BLOCK_58 ]
  %13 = icmp slt i64 %12, %3
  br i1 %13, label %BLOCK_64, label %BLOCK_165

BLOCK_165:                                        ; preds = %BLOCK_61
  br label %BLOCK_39

BLOCK_64:                                         ; preds = %BLOCK_61
  br label %BLOCK_73

BLOCK_73:                                         ; preds = %BLOCK_64
  br label %BLOCK_86

BLOCK_83:                                         ; preds = %BLOCK_89
  %14 = add i64 %16, 1
  br label %BLOCK_86

BLOCK_86:                                         ; preds = %BLOCK_83, %BLOCK_73
  %15 = phi double [ 0.000000e+00, %BLOCK_73 ], [ %53, %BLOCK_83 ]
  %16 = phi i64 [ 0, %BLOCK_73 ], [ %14, %BLOCK_83 ]
  %17 = icmp slt i64 %16, %7
  br i1 %17, label %BLOCK_89, label %BLOCK_136

BLOCK_136:                                        ; preds = %BLOCK_86
  %18 = call double @llvm.sqrt.f64(double %15)
  %19 = getelementptr { i64, i32*, i8*, i32, i64*, i64*, i8*, i8*, i32, i8*, i8*, i8*, i64* }* %D, i32 0, i32 5
  %20 = load i64** %19
  %21 = getelementptr i64* %20, i32 0
  %22 = load i64* %21
  %23 = mul i64 %9, %22
  %24 = getelementptr { i64, i32*, i8*, i32, i64*, i64*, i8*, i8*, i32, i8*, i8*, i8*, i64* }* %D, i32 0, i32 2
  %25 = load i8** %24
  %26 = getelementptr i8* %25, i64 %23
  %27 = bitcast i8* %26 to double*
  %28 = getelementptr double* %27, i64 %12
  store double %18, double* %28
  br label %BLOCK_58

BLOCK_89:                                         ; preds = %BLOCK_86
  %29 = getelementptr { i64, i32*, i8*, i32, i64*, i64*, i8*, i8*, i32, i8*, i8*, i8*, i64* }* %X, i32 0, i32 5
  %30 = load i64** %29
  %31 = getelementptr i64* %30, i32 0
  %32 = load i64* %31
  %33 = mul i64 %9, %32
  %34 = getelementptr { i64, i32*, i8*, i32, i64*, i64*, i8*, i8*, i32, i8*, i8*, i8*, i64* }* %X, i32 0, i32 2
  %35 = load i8** %34
  %36 = getelementptr i8* %35, i64 %33
  %37 = bitcast i8* %36 to double*
  %38 = getelementptr double* %37, i64 %16
  %39 = load double* %38
  %40 = getelementptr { i64, i32*, i8*, i32, i64*, i64*, i8*, i8*, i32, i8*, i8*, i8*, i64* }* %X, i32 0, i32 5
  %41 = load i64** %40
  %42 = getelementptr i64* %41, i32 0
  %43 = load i64* %42
  %44 = mul i64 %12, %43
  %45 = getelementptr { i64, i32*, i8*, i32, i64*, i64*, i8*, i8*, i32, i8*, i8*, i8*, i64* }* %X, i32 0, i32 2
  %46 = load i8** %45
  %47 = getelementptr i8* %46, i64 %44
  %48 = bitcast i8* %47 to double*
  %49 = getelementptr double* %48, i64 %16
  %50 = load double* %49
  %51 = fsub double %39, %50
  %52 = fmul double %51, %51
  %53 = fadd double %15, %52
  br label %BLOCK_83
}

declare double @llvm.sqrt.f64(double) nounwind readonly

In [4]:
%%cython

cimport cython
from libc.math cimport sqrt

@cython.boundscheck(False)
@cython.wraparound(False)
def pairwise_cython(double[:, ::1] X, double[:, ::1] D):
    cdef int M = X.shape[0]
    cdef int N = X.shape[1]
    cdef double tmp, d
    for i in range(M):
        for j in range(M):
            d = 0.0
            for k in range(N):
                tmp = X[i, k] - X[j, k]
                d += tmp * tmp
            D[i, j] = sqrt(d)
In [5]:
%timeit pairwise_python(X, D)
%timeit pairwise_numba(X, D)
%timeit pairwise_cython(X, D)
1 loops, best of 3: 12.1 s per loop
100 loops, best of 3: 15.5 ms per loop
100 loops, best of 3: 9.86 ms per loop