import numpy as np

X = np.random.random((1000, 3))
X_wide = np.random.random((1000, 100))

def pairwise_numpy(X):
    return np.sqrt(((X[:, None, :] - X) ** 2).sum(-1))
%timeit pairwise_numpy(X)

def pairwise_python(X):
    M = X.shape[0]
    N = X.shape[1]
    D = np.empty((M, M), dtype=np.float)
    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)
    return D

%timeit pairwise_python(X)

def pairwise_python2(X):
    n_samples = X.shape[0]
    result = np.zeros((n_samples, n_samples), dtype=X.dtype)
    for i in xrange(X.shape[0]):
        for j in xrange(X.shape[0]):
            result[i, j] = np.sqrt(np.sum((X[i, :] - X[j, :]) ** 2))
    return result

%timeit pairwise_python2(X)

#np.allclose(pairwise_python(X), pairwise_python2(X))

import numba

numba.__version__

from numba import double
from numba.decorators import jit, autojit

pairwise_numba = autojit(pairwise_python)

%timeit pairwise_numba(X)

%timeit pairwise_numba(X_wide)

pairwise_numba2 = autojit(pairwise_python2)

%timeit pairwise_numba2(X)

from parakeet import jit

pairwise_parakeet = jit(pairwise_python)

%timeit pairwise_parakeet(X)

%timeit pairwise_parakeet(X_wide)

pairwise_parakeet2 = jit(pairwise_python2)
%timeit pairwise_parakeet2(X)

%timeit pairwise_parakeet2(X_wide)

np.allclose(pairwise_parakeet(X), pairwise_parakeet2(X))

np.allclose(pairwise_parakeet(X_wide), pairwise_parakeet2(X_wide))

!cython --version

%load_ext cythonmagic

%%cython

import numpy as np
cimport cython
from libc.math cimport sqrt

@cython.boundscheck(False)
@cython.wraparound(False)
def pairwise_cython(double[:, ::1] X):
    cdef int M = X.shape[0]
    cdef int N = X.shape[1]
    cdef double tmp, d
    cdef double[:, ::1] D = np.empty((M, M), dtype=np.float64)
    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)
    return np.asarray(D)

%timeit pairwise_cython(X)

%timeit pairwise_cython(X_wide)

%%file pairwise_fortran.f

      subroutine pairwise_fortran(X,D,m,n)
          integer :: n,m
          double precision, intent(in) :: X(m,n)
          double precision, intent(out) :: D(m,m) 
          integer :: i,j,k
          double precision :: r 
          do i = 1,m 
              do j = 1,m 
                  r = 0
                  do k = 1,n 
                      r = r + (X(i,k) - X(j,k)) * (X(i,k) - X(j,k)) 
                  end do 
                  D(i,j) = sqrt(r) 
              end do 
          end do 
      end subroutine pairwise_fortran

# Compile the Fortran with f2py.
# We'll direct the output into /dev/null so it doesn't fill the screen
!f2py -c pairwise_fortran.f -m pairwise_fortran > /dev/null

from pairwise_fortran import pairwise_fortran
XF = np.asarray(X, order='F')
%timeit pairwise_fortran(XF)

XF_wide = np.asarray(X_wide, order='F')
%timeit pairwise_fortran(XF_wide)

from scipy.spatial.distance import cdist
%timeit cdist(X, X)

%timeit cdist(X_wide, X_wide)

from sklearn.metrics import euclidean_distances
%timeit euclidean_distances(X, X)

%timeit euclidean_distances(X_wide, X_wide)