In [1]:
import sys
print sys.version
import numpy; print 'numpy', numpy.version.full_version
import scipy; print 'scipy', scipy.version.full_version
import numexpr; print 'numexpr', numexpr.version.version
import Cython; print Cython.__version__
In [2]:
import numpy.random
a = numpy.random.rand(20000).astype(np.float32)
In [3]:
import scipy.misc
print(scipy.misc.logsumexp(a))
%timeit scipy.misc.logsumexp(a)
In [4]:
import scipy.weave
def lse_weave(u): # no timing difference taking it out of the function, ie leaving scipy.weave.inline() bare
code = """
int i;
float r = 0;
float largest_in_u = 0;
for (i=0; i<Nu[0]; i++) {
if (U1(i) > largest_in_u) {
largest_in_u = U1(i);
}
}
for (i=0; i<Nu[0]; i++) {
r += exp(U1(i) - largest_in_u);
}
return_val = largest_in_u + log(r);
"""
return scipy.weave.inline(code, ['u']) # no change using type_converters=converters.blitz
print(lse_weave(a))
%timeit lse_weave(a)
In [5]:
import numexpr
def lse_numexpr(a):
largest_in_a = np.max(a)
r = numexpr.evaluate("sum(exp(a - largest_in_a))")
return largest_in_a + np.log(r)
# np.log(b) faster than numexpr.evaluate("log(r)")
# numexpr.evaluate("log(sum(exp(a)))") will return an error
print(lse_numexpr(a))
%timeit lse_numexpr(a)
In [6]:
%load_ext cythonmagic
In [7]:
%%cython
cimport cython
import numpy as np
cimport numpy as np
from libc.math cimport exp, log # 40x speedup using this instead of np.exp, np.log which result in python numpy calls
DTYPE=np.float32
ctypedef np.float32_t DTYPE_t
@cython.boundscheck(False)
@cython.wraparound(False)
cpdef lse_cython(np.ndarray[DTYPE_t, ndim=1] a) :
cdef int i
cdef double result = 0.0
cdef double largest_in_a = a[0]
for i in range(1,a.shape[0]):
if (a[i] > largest_in_a):
largest_in_a = a[i]
for i in range(a.shape[0]):
result += exp(a[i] - largest_in_a)
return largest_in_a + log(result)
In [8]:
print(lse_cython(a))
%timeit lse_cython(a)
In [9]:
import sselogsumexp
sselogsumexp.logsumexp(a)
%timeit sselogsumexp.logsumexp(a)
In [10]:
print("%timeit scipy.misc.logsumexp(a)")
print(scipy.misc.logsumexp(a))
%timeit scipy.misc.logsumexp(a)
print(lse_weave(a))
print("%timeit lse_weave(a)")
%timeit lse_weave(a)
print(lse_numexpr(a))
print("%timeit lse_numexpr(a)")
%timeit lse_numexpr(a)
print(lse_cython(a))
print("%timeit lse_cython(a)")
%timeit lse_cython(a)
print sselogsumexp.logsumexp(a)
print("%timeit sselogsumexp.logsumexp(a)")
%timeit sselogsumexp.logsumexp(a)
In [10]: