Going to look at faster ways to make wedges.
% matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
Set up the model parameters.
length = 100.# x range
width = 100. # y range
depth = 240. # z range
Define the plotting function.
def plot_slices(data, interpolation='nearest', cmap='Greys', vmin=0., vmax=1.):
plt.figure(figsize=(16,12))
for j in range(3):
plt.subplot(3,4,j+1)
plt.title('y = ' + str(int(j*(length/3))))
plt.imshow(data[:,:,int(j*(length/3))], interpolation=interpolation, cmap=cmap, vmax=vmax, vmin=vmin)
plt.subplot(3,4,4)
plt.title('y = ' + str(int(length-1)))
plt.imshow(data[:,:,-1], interpolation=interpolation, cmap=cmap, vmax=vmax, vmin=vmin)
for j in range(3):
plt.subplot(3,4,5+j)
plt.title('x = ' + str(int(j*(width/3))))
plt.imshow(data[:,int(j*(width/3)),:], interpolation=interpolation, cmap=cmap, vmax=vmax, vmin=vmin)
plt.subplot(3,4,8)
plt.title('x = ' + str(int(width-1)))
plt.imshow(data[:,-1,:], interpolation=interpolation, cmap=cmap, vmax=vmax, vmin=vmin)
for j in range(3):
plt.subplot(3,4,9+j)
plt.title('z = ' + str(1 + int(1 + depth/3 + j*depth/(3*3))))
plt.imshow(data[1 + int(depth/3 + j*depth/(3*3)),:,:], interpolation=interpolation, cmap=cmap, vmax=vmax, vmin=vmin)
plt.subplot(3,4,12)
plt.title('z = ' + str(int(2*depth/3 - 1)))
plt.imshow(data[int(2*depth/3 - 1),:,:], interpolation=interpolation, cmap=cmap, vmax=vmax, vmin=vmin)
plt.show()
Let's start with the wedge code from version 1 of Spectral_wedge.ipynb
. The main difference is that we will need to use the rock properties directly, rather than building an array of int
s first.
rocks = [(2800, 2850), (2800, 2850), (2400, 2450), (2400,2450)]
def make_ai(rocks):
rx = np.array(rocks)
return rx[:,0] * rx[:,1] / 10e6
ai = make_ai(rocks)
print(ai)
[0.798 0.798 0.588 0.588]
def make_earth(depth, length, width, ai):
earth = np.ones((depth,length,width))
for xslice in range(int(length)):
thickness = (depth/3)*(1 + xslice/length)
value = (ai[1]*xslice/length) + (ai[2]*(length-xslice)/length)
earth[:,xslice,:] *= value
# Draw the bottom of the wedge
value = ai[0]
for (z,y),v in np.ndenumerate(earth[:,xslice,:]):
if z > thickness:
earth[z,xslice,y] = value
# Put the top on
value = ai[3]
for (z,y),v in np.ndenumerate(earth[:,xslice,:]):
if z < depth/3:
earth[z,xslice,y] = value
return earth
%timeit earth = make_earth(depth, length, width, ai)
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-7-e75eaee6f09e> in <module>() ----> 1 get_ipython().run_line_magic('timeit', 'earth = make_earth(depth, length, width, ai)') ~/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py in run_line_magic(self, magic_name, line, _stack_depth) 2093 kwargs['local_ns'] = sys._getframe(stack_depth).f_locals 2094 with self.builtin_trap: -> 2095 result = fn(*args,**kwargs) 2096 return result 2097 <decorator-gen-61> in timeit(self, line, cell, local_ns) ~/anaconda3/lib/python3.6/site-packages/IPython/core/magic.py in <lambda>(f, *a, **k) 185 # but it's overkill for just that one bit of state. 186 def magic_deco(arg): --> 187 call = lambda f, *a, **k: f(*a, **k) 188 189 if callable(arg): ~/anaconda3/lib/python3.6/site-packages/IPython/core/magics/execution.py in timeit(self, line, cell, local_ns) 1096 for index in range(0, 10): 1097 number = 10 ** index -> 1098 time_number = timer.timeit(number) 1099 if time_number >= 0.2: 1100 break ~/anaconda3/lib/python3.6/site-packages/IPython/core/magics/execution.py in timeit(self, number) 158 gc.disable() 159 try: --> 160 timing = self.inner(it, self.timer) 161 finally: 162 if gcold: <magic-timeit> in inner(_it, _timer) <ipython-input-6-ae58368e8d06> in make_earth(depth, length, width, ai) 1 def make_earth(depth, length, width, ai): 2 ----> 3 earth = np.ones((depth,length,width)) 4 5 for xslice in range(int(length)): ~/anaconda3/lib/python3.6/site-packages/numpy/core/numeric.py in ones(shape, dtype, order) 186 187 """ --> 188 a = empty(shape, dtype, order) 189 multiarray.copyto(a, 1, casting='unsafe') 190 return a TypeError: 'float' object cannot be interpreted as an integer
earth = make_earth(depth, length, width, ai)
earth.shape
(240, 100, 100)
You can index an array with an array-like object — a list or another array.
a = np.array([1,2,3,4,5,6,7,8]) # A 1D array
b = [1,3,5,7]
print a[b]
d = np.array([[1,2,3,4],[5,6,7,8]]) # A 2D array...
e = [[0,1],[1,3]] # Needs 2 indices
print d[e]
[2 4 6 8] [2 8]
Now let's try using index lists, building the indices with nested list comprehensions.
def use_index_lists(depth, length, width, ai):
earth = np.ones((depth,length,width))
# Bottom indices
indices = [[[x for x in range(int(length))] for y in range(int(width))] for z in range(int(depth/2))]
earth[indices] = 2.
return earth
# Only getting started, let's test it...
%timeit earth_index_lists = use_index_lists(depth, length, width, ai)
earth_index_lists = use_index_lists(depth, length, width, ai)
earth_index_lists.shape
That was unbelievably slow. Forget it.
It occurs to me I have too many loops. I think I need one, because the values I want vary with index.
I just realized that instead of enumerating and deciding if the z-index is greater than thickness, say, I can just slice with thickness as the index. Doh! NumPy is awesome...
def one_loop(depth, length, width, ai):
earth = np.ones((depth,length,width))
d = depth/3
for xslice in range(int(length)):
xl = xslice/length
thickness = d*(1 + xl)
mix = (ai[1]*xl) + (ai[2]*(length-xslice)/length)
# Draw wedge
earth[:,xslice,:] *= mix
# Put the top on
earth[:d,:,:] = ai[3]
return earth
%timeit one_loop(depth, length, width, ai)
100 loops, best of 3: 13.8 ms per loop
Holy crap, that's a bit better.
earth = one_loop(depth, length, width, ai)
plot_slices(earth, interpolation='nearest', cmap='jet')
The last tricky bit is the outer loop. It's only there to loop over the x-slices. The tricky part is that the wedge thickness calculation depends on the index of the slice. I think numpy.indices
can provide this aspect and avoid the loop...
zs, xs, ys = np.indices((4,3,2))
print "x values:\n", xs
print "\ny values:\n", ys
print "\nz values:\n", zs
x values: [[[0 0] [1 1] [2 2]] [[0 0] [1 1] [2 2]] [[0 0] [1 1] [2 2]] [[0 0] [1 1] [2 2]]] y values: [[[0 1] [0 1] [0 1]] [[0 1] [0 1] [0 1]] [[0 1] [0 1] [0 1]] [[0 1] [0 1] [0 1]]] z values: [[[0 0] [0 0] [0 0]] [[1 1] [1 1] [1 1]] [[2 2] [2 2] [2 2]] [[3 3] [3 3] [3 3]]]
def no_loop(depth, length, width, ai):
# Make an earth full of ones
earth = np.ones((depth,length,width))
# Set up the gradient, which we'll then crop, so to speak
mix = np.array([(ai[1]*x/length) + (ai[2]*(length-x)/length) for x in range(int(length))])
earth *= mix[:, np.newaxis]
# Now make the wedge
d = depth/3
thickness = [d*(1 + x/length) for x in range(int(length))]
# Make a set of index arrays
zs = np.indices(earth.shape)[0]
# Make a boolean array map of stuff below the wedge
below = zs > thickness
# I don't know why I have to do this but I do
below = np.swapaxes(below, axis1=1, axis2=2)
# Set those values to the lowermost rock AI
earth[below] = ai[0]
# Put the top layer on
earth[:d,:,:] = ai[3]
return earth
%timeit earth = no_loop(depth, length, width, ai)
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-9-56ebd4328593> in <module>() ----> 1 get_ipython().run_line_magic('timeit', 'earth = no_loop(depth, length, width, ai)') ~/anaconda3/lib/python3.6/site-packages/IPython/core/interactiveshell.py in run_line_magic(self, magic_name, line, _stack_depth) 2093 kwargs['local_ns'] = sys._getframe(stack_depth).f_locals 2094 with self.builtin_trap: -> 2095 result = fn(*args,**kwargs) 2096 return result 2097 <decorator-gen-61> in timeit(self, line, cell, local_ns) ~/anaconda3/lib/python3.6/site-packages/IPython/core/magic.py in <lambda>(f, *a, **k) 185 # but it's overkill for just that one bit of state. 186 def magic_deco(arg): --> 187 call = lambda f, *a, **k: f(*a, **k) 188 189 if callable(arg): ~/anaconda3/lib/python3.6/site-packages/IPython/core/magics/execution.py in timeit(self, line, cell, local_ns) 1096 for index in range(0, 10): 1097 number = 10 ** index -> 1098 time_number = timer.timeit(number) 1099 if time_number >= 0.2: 1100 break ~/anaconda3/lib/python3.6/site-packages/IPython/core/magics/execution.py in timeit(self, number) 158 gc.disable() 159 try: --> 160 timing = self.inner(it, self.timer) 161 finally: 162 if gcold: <magic-timeit> in inner(_it, _timer) <ipython-input-8-28ca7de9b565> in no_loop(depth, length, width, ai) 2 3 # Make an earth full of ones ----> 4 earth = np.ones((depth,length,width)) 5 6 # Set up the gradient, which we'll then crop, so to speak ~/anaconda3/lib/python3.6/site-packages/numpy/core/numeric.py in ones(shape, dtype, order) 186 187 """ --> 188 a = empty(shape, dtype, order) 189 multiarray.copyto(a, 1, casting='unsafe') 190 return a TypeError: 'float' object cannot be interpreted as an integer
def no_loop(depth, length, width, ai):
# Make an earth full of ones
earth = np.ones((depth,length,width))
# Set up the gradient, which we'll then crop, so to speak
mix = np.array([(ai[1]*x/length) + (ai[2]*(length-x)/length) for x in range(int(length))])
earth *= mix[:, np.newaxis]
# Now make the wedge
d = depth/3
thickness = [d*(1 + x/length) for x in range(int(length))]
# Make a set of index arrays
zs = np.indices(earth.shape)[0]
# Make a boolean array map of stuff below the wedge
below = zs > thickness
# I don't know why I have to do this but I do
below = np.swapaxes(below, axis1=1, axis2=2)
# Set those values to the lowermost rock AI
earth[below] = ai[0]
# Put the top layer on
earth[:d,:,:] = ai[3]
return earth
depth, length, width, ai
(240.0, 100.0, 100.0, array([0.798, 0.798, 0.588, 0.588]))
earth = no_loop(depth, length, width, ai)
plot_slices(earth, interpolation='nearest', cmap='jet')
Amazing — the for
loop is 3 times faster! And I think it's the indices
part that is slowing it down, because without that section (the line with np.indices
and the three lines after it), it only takes 9.5 ms to execute.
Assuming I'm not doing something stupid in this code, it seems it's not always faster with indexing.