Defining ufuncs using vectorize
You have been able to define your own NumPy ufuncs for quite some time, but it's a little involved.
You can read through the documentation, the example they post there is a ufunc to perform
It looks like this:
static void double_logit(char **args, npy_intp *dimensions,
npy_intp* steps, void* data)
{
npy_intp i;
npy_intp n = dimensions[0];
char *in = args[0], *out = args[1];
npy_intp in_step = steps[0], out_step = steps[1];
double tmp;
for (i = 0; i < n; i++) {
/*BEGIN main ufunc computation*/
tmp = *(double *)in;
tmp /= 1-tmp;
*((double *)out) = log(tmp);
/*END main ufunc computation*/
in += in_step;
out += out_step;
}
}
And note, that's just for a double. If you want floats, long doubles, etc... you have to write all of those, too. And then create a setup.py file to install it. And I left out a bunch of boilerplate stuff to set up the import hooks, etc...
Say "thank you" to the NumPy devs
We can use Numba to define ufuncs without all of the pain.
import numpy
import math
Let's define a function that operates on two inputs
def trig(a, b):
return math.sin(a**2) * math.exp(b)
trig(1, 1)
2.2873552871788423
Seems reasonable. However, the math library only works on scalars. If we try to pass in arrays, we'll get an error.
a = numpy.ones((5,5))
b = numpy.ones((5,5))
trig(a, b)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
<ipython-input-5-6bc62cd8d328> in <module>()
----> 1 trig(a, b)
<ipython-input-2-27083e35e9e9> in trig(a, b)
1 def trig(a, b):
----> 2 return math.sin(a**2) * math.exp(b)
TypeError: only length-1 arrays can be converted to Python scalars
from numba import vectorize
vec_trig = vectorize()(trig)
vec_trig(a, b)
array([[ 2.28735529, 2.28735529, 2.28735529, 2.28735529, 2.28735529],
[ 2.28735529, 2.28735529, 2.28735529, 2.28735529, 2.28735529],
[ 2.28735529, 2.28735529, 2.28735529, 2.28735529, 2.28735529],
[ 2.28735529, 2.28735529, 2.28735529, 2.28735529, 2.28735529],
[ 2.28735529, 2.28735529, 2.28735529, 2.28735529, 2.28735529]])
And just like that, the scalar function trig is now a NumPy ufunc called vec_trig
Note that this is a "Dynamic UFunc" with no signature given.
How does it compare to just using NumPy? Let's check
def numpy_trig(a, b):
return numpy.sin(a**2) * numpy.exp(b)
a = numpy.random.random((1000, 1000))
b = numpy.random.random((1000, 1000))
%timeit vec_trig(a, b)
10 loops, best of 3: 29.4 ms per loop
%timeit numpy_trig(a, b)
10 loops, best of 3: 32.5 ms per loop
What happens if we do specify a signature? Is there a speed boost?
vec_trig = vectorize('float64(float64, float64)')(trig)
%timeit vec_trig(a, b)
10 loops, best of 3: 29.5 ms per loop
No, not really. But(!), if we have a signature, then we can add the target kwarg.
vec_trig = vectorize('float64(float64, float64)', target='parallel')(trig)
%timeit vec_trig(a, b)
100 loops, best of 3: 6.24 ms per loop
Automatic multicore operations!
Note: target='parallel' is not always the best option. There is overhead in setting up the threading, so if the individual scalar operations that make up a ufunc are simple you'll probably get better performance in serial. If the individual operations are more expensive (like trig!) then parallel is (usually) a good option.
Passing multiple signatures
If you use multiple signatures, they have to be listed in order of most specific -> least specific
@vectorize(['int32(int32, int32)',
'int64(int64, int64)',
'float32(float32, float32)',
'float64(float64, float64)'])
def trig(a, b):
return math.sin(a**2) * math.exp(b)
trig(1, 1)
2
trig(1., 1.)
2.2873552871788423
trig.ntypes
4
Exercise: Clipping an array
Yes, NumPy has a clip ufunc already, but let's pretend it doesn't.
Create a Numba vectorized ufunc that takes a vector a, a lower limit amin and an upper limit amax. It should return the vector a with all values clipped such that :
# %load snippets/clip.py
a = numpy.random.random((5000))
amin = .2
amax = .6
%timeit vec_truncate_serial(a, amin, amax)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-24-1221b70b3424> in <module>()
----> 1 get_ipython().magic('timeit vec_truncate_serial(a, amin, amax)')
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/interactiveshell.py in magic(self, arg_s)
2156 magic_name, _, magic_arg_s = arg_s.partition(' ')
2157 magic_name = magic_name.lstrip(prefilter.ESC_MAGIC)
-> 2158 return self.run_line_magic(magic_name, magic_arg_s)
2159
2160 #-------------------------------------------------------------------------
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/interactiveshell.py in run_line_magic(self, magic_name, line)
2077 kwargs['local_ns'] = sys._getframe(stack_depth).f_locals
2078 with self.builtin_trap:
-> 2079 result = fn(*args,**kwargs)
2080 return result
2081
<decorator-gen-58> in timeit(self, line, cell)
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/magic.py in <lambda>(f, *a, **k)
186 # but it's overkill for just that one bit of state.
187 def magic_deco(arg):
--> 188 call = lambda f, *a, **k: f(*a, **k)
189
190 if callable(arg):
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/magics/execution.py in timeit(self, line, cell)
1042 number = 1
1043 for _ in range(1, 10):
-> 1044 time_number = timer.timeit(number)
1045 worst_tuning = max(worst_tuning, time_number / number)
1046 if time_number >= 0.2:
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/magics/execution.py in timeit(self, number)
137 gc.disable()
138 try:
--> 139 timing = self.inner(it, self.timer)
140 finally:
141 if gcold:
<magic-timeit> in inner(_it, _timer)
NameError: name 'vec_truncate_serial' is not defined
%timeit vec_truncate_par(a, amin, amax)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-25-a1fc14977f98> in <module>()
----> 1 get_ipython().magic('timeit vec_truncate_par(a, amin, amax)')
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/interactiveshell.py in magic(self, arg_s)
2156 magic_name, _, magic_arg_s = arg_s.partition(' ')
2157 magic_name = magic_name.lstrip(prefilter.ESC_MAGIC)
-> 2158 return self.run_line_magic(magic_name, magic_arg_s)
2159
2160 #-------------------------------------------------------------------------
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/interactiveshell.py in run_line_magic(self, magic_name, line)
2077 kwargs['local_ns'] = sys._getframe(stack_depth).f_locals
2078 with self.builtin_trap:
-> 2079 result = fn(*args,**kwargs)
2080 return result
2081
<decorator-gen-58> in timeit(self, line, cell)
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/magic.py in <lambda>(f, *a, **k)
186 # but it's overkill for just that one bit of state.
187 def magic_deco(arg):
--> 188 call = lambda f, *a, **k: f(*a, **k)
189
190 if callable(arg):
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/magics/execution.py in timeit(self, line, cell)
1042 number = 1
1043 for _ in range(1, 10):
-> 1044 time_number = timer.timeit(number)
1045 worst_tuning = max(worst_tuning, time_number / number)
1046 if time_number >= 0.2:
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/magics/execution.py in timeit(self, number)
137 gc.disable()
138 try:
--> 139 timing = self.inner(it, self.timer)
140 finally:
141 if gcold:
<magic-timeit> in inner(_it, _timer)
NameError: name 'vec_truncate_par' is not defined
%timeit numpy.clip(a, amin, amax)
100000 loops, best of 3: 4.78 µs per loop
a = numpy.random.random((100000))
%timeit vec_truncate_serial(a, amin, amax)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-28-1221b70b3424> in <module>()
----> 1 get_ipython().magic('timeit vec_truncate_serial(a, amin, amax)')
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/interactiveshell.py in magic(self, arg_s)
2156 magic_name, _, magic_arg_s = arg_s.partition(' ')
2157 magic_name = magic_name.lstrip(prefilter.ESC_MAGIC)
-> 2158 return self.run_line_magic(magic_name, magic_arg_s)
2159
2160 #-------------------------------------------------------------------------
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/interactiveshell.py in run_line_magic(self, magic_name, line)
2077 kwargs['local_ns'] = sys._getframe(stack_depth).f_locals
2078 with self.builtin_trap:
-> 2079 result = fn(*args,**kwargs)
2080 return result
2081
<decorator-gen-58> in timeit(self, line, cell)
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/magic.py in <lambda>(f, *a, **k)
186 # but it's overkill for just that one bit of state.
187 def magic_deco(arg):
--> 188 call = lambda f, *a, **k: f(*a, **k)
189
190 if callable(arg):
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/magics/execution.py in timeit(self, line, cell)
1042 number = 1
1043 for _ in range(1, 10):
-> 1044 time_number = timer.timeit(number)
1045 worst_tuning = max(worst_tuning, time_number / number)
1046 if time_number >= 0.2:
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/magics/execution.py in timeit(self, number)
137 gc.disable()
138 try:
--> 139 timing = self.inner(it, self.timer)
140 finally:
141 if gcold:
<magic-timeit> in inner(_it, _timer)
NameError: name 'vec_truncate_serial' is not defined
%timeit vec_truncate_par(a, amin, amax)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-29-a1fc14977f98> in <module>()
----> 1 get_ipython().magic('timeit vec_truncate_par(a, amin, amax)')
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/interactiveshell.py in magic(self, arg_s)
2156 magic_name, _, magic_arg_s = arg_s.partition(' ')
2157 magic_name = magic_name.lstrip(prefilter.ESC_MAGIC)
-> 2158 return self.run_line_magic(magic_name, magic_arg_s)
2159
2160 #-------------------------------------------------------------------------
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/interactiveshell.py in run_line_magic(self, magic_name, line)
2077 kwargs['local_ns'] = sys._getframe(stack_depth).f_locals
2078 with self.builtin_trap:
-> 2079 result = fn(*args,**kwargs)
2080 return result
2081
<decorator-gen-58> in timeit(self, line, cell)
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/magic.py in <lambda>(f, *a, **k)
186 # but it's overkill for just that one bit of state.
187 def magic_deco(arg):
--> 188 call = lambda f, *a, **k: f(*a, **k)
189
190 if callable(arg):
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/magics/execution.py in timeit(self, line, cell)
1042 number = 1
1043 for _ in range(1, 10):
-> 1044 time_number = timer.timeit(number)
1045 worst_tuning = max(worst_tuning, time_number / number)
1046 if time_number >= 0.2:
/home/gil/anaconda/lib/python3.5/site-packages/IPython/core/magics/execution.py in timeit(self, number)
137 gc.disable()
138 try:
--> 139 timing = self.inner(it, self.timer)
140 finally:
141 if gcold:
<magic-timeit> in inner(_it, _timer)
NameError: name 'vec_truncate_par' is not defined
%timeit numpy.clip(a, amin, amax)
1000 loops, best of 3: 212 µs per loop
Exercise: Create logit ufunc
Recall from above that this is a ufunc which performs this operation:
# %load snippets/logit.py
logit(a)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-32-43208ff16126> in <module>()
----> 1 logit(a)
NameError: name 'logit' is not defined
Performance of vectorize vs. regular array-wide operations
@vectorize
def discriminant(a, b, c):
return b**2 - 4 * a * c
a = numpy.arange(10000)
b = numpy.arange(10000)
c = numpy.arange(10000)
%timeit discriminant(a, b, c)
The slowest run took 3182.90 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 11.7 µs per loop
%timeit b**2 - 4 * a * c
The slowest run took 5.46 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 39.4 µs per loop
What's going on?
- Each array operation creates a temporary copy
- Each of these arrays are loaded into and out of cache a whole bunch
del a, b, c