In [1]:
######################################################################3
# Example: Timing
#
# Timing in Python is easily accomplished using the "time"
# function provided by the "time" module. A simple
# example is provided below.
import time
nsizes = 8
for i in range(nsizes):
list1 = []
list2 = []
nelem = 10**i
t0 = time.time()
for j in range(nelem):
list1.append(j)
list2.append(i+j)
t1 = time.time()
dt_create = t1-t0
t0 = time.time()
for j in range(nelem):
list2[j] = list2[j]*list1[j]
t1 = time.time()
dt_prod = t1-t0
nestr = '{:8}'.format(nelem)
dtcstr = '{:.4e}'.format(dt_create)
dtpstr = '{:.4e}'.format(dt_prod)
print('')
print('/////////////////////////////////////////////////////////////////')
print(' Time to create two lists of length '+nestr+' : '+dtcstr+' seconds')
print(' Time to perform element-wise product on two lists of length '+nestr+' : '+dtpstr+' seconds')
In [ ]:
#################################################################
# Example: Numpy array initialization
#
# There are various ways to initialize a Numpy array.
# Several examples are shown below.
import numpy as np
#Array dimensions for use later
nx = 10
ny = 20
#1-D array, 8-byte real, using np.empty means that
# values are not initialized to anything (in particular...)
earr = np.empty(nx,dtype='float64')
print('earr: ', earr)
print(' ')
#1-D array, 4-byte real, using np.zeros
# means that values are initialized to zero
zarr = np.zeros(nx,dtype='float32')
print('zarr: ', zarr)
print(' ')
#2-d array, 4-byte integer, values set to zero
# Row-major ordering (second index is fastest; like C/C++)
iarr2da = np.zeros((nx,ny),dtype='int32')
print('iarr2da: ')
print(iarr2da)
print(' ')
#2-d array, 4-byte integer, values set to zero
#Column-major ordering (first index is fastest; like Fortran)
iarr2db = np.zeros((nx,ny),dtype='int32', order ='F')
print('iarr2db: ')
print(iarr2db)
print(' ')
#1-d array, 4-byte real, values spaced with integer spacing
# in the range [istart,iend], with stepsize istep
# This is accomplished via np.arange
istart = 10
iend = 20
istep = 2
arrsp = np.arange(istart,iend,istep,dtype='float32')
print('arrsp: ', arrsp)
print(' ')
#1-d array, 8-byte real, values spaced with non-integer spacing
# in the range [istart,iend], with stepsize = (iend-istart)/nstep
istart = 100
iend = 200
nstep = 200
arrsp2 = np.linspace(istart,iend,nstep,dtype='float64')
print('arrsp2: ')
print(arrsp2)
print(' ')
# 1-d array, 8-byte real, initialized using values from an existing list
oned = [0,1.2, 5.6, 8.9]
arrinit = np.array(oned,dtype='float64')
print('arrinit: ', arrinit)
print(' ')
# 2-d array, 4-byte integer, initialized using values from existing 2-D list
twod = [ [ 1, 2], [3,4] ]
arrinit2d = np.array(twod,dtype='float64')
print('arrinit2d: ')
print(arrinit2d)
print(' ')
array_names = ['zarr', 'earr', 'iarr2da', 'iarr2db','arrsp', 'arrsp2', 'arrinit', 'arrinit2d']
arrays = [zarr,earr,iarr2da,iarr2db, arrsp, arrsp2, arrinit, arrinit2d]
for i,o in enumerate(arrays):
ndim = o.ndim
dtype = o.dtype
isize = o.itemsize
ssp = o.nbytes
ne = o.size
print('////////////////////////////////////////')
print(' Information for ndarray '+array_names[i])
print(' data type : ', dtype)
print(' number of elements : ', ne)
print(' element size (bytes) : ', isize)
print(' dimensions : ', ndim)
for j in range(ndim):
print(' dimension '+str(j)+' size :', o.shape[j])
print(' storage space (bytes): ', ssp)
if (ndim > 1):
print(' Element spacing along dimension 1 (bytes): ', o.strides[0])
print(' Element spacing along dimension 2 (bytes): ', o.strides[1])
print('')
In [ ]:
#//////////////////////////////////////////////////////////////
# Exercise 1:
# Rewrite the following program using numpy arrays and
# array operations where possible (instead of explicit loops).
n = 1000000
a = []
b = []
for i in range(1,n+1):
a.append(4*i)
b.append(i**2)
dsum = 0.0
for i in range(n):
dsum += a[i]*b[i]
print('dsum is: ', dsum)
In [ ]:
import numpy as np
n = 1000000
a = 4 * (np.arange(n) + 1)
b = np.square(np.arange(n) + 1)
#dsum = np.sum(a * b, dtype='float64')
dsum = np.sum(a.astype('float64') * b)
print("Our sum is : ", dsum)
In [ ]:
###########################################################################
# Example: In-place vs. out-of-place operations
#
# We can avoid unnecessary array copying, and save time,
# by using in-place operators where possible.
# Use a += 2 instead of a = a+2, a *=2 instead of a = a*2, etc.
import numpy as np
import time
npts = 10000000
ntrials = 4
a = np.zeros(npts)
b = np.zeros(npts)
print(' ')
print(' Timing in-place vs. out-of-place array operations')
print(' Number of elements: ', npts)
print(' Number of trials : ', ntrials)
print('')
# This appears to be in-place, but in fact a new array is made (a*2) and reassigned to a
t0 = time.time()
for i in range(ntrials):
a = a*2
t1 = time.time()
dt1 = t1-t0
# This is truly in-place
t0 = time.time()
for i in range(ntrials):
a *= 2
t1 = time.time()
dt2 = t1-t0
# And here, we have a clearly out-of-place operation (a*2 is calculated and then assigned to b)
t0 = time.time()
for i in range(ntrials):
b = a*2
t1 = time.time()
dt3 = t1-t0
tstr1 = '{:.4e}'.format(dt1)
tstr2 = '{:.4e}'.format(dt2)
tstr3 = '{:.4e}'.format(dt3)
print(' "In-place" ( a = a*2 ) : '+tstr1+' seconds')
print(' In-place ( a *= 2 ) : '+tstr2+' seconds')
print(' Out-of-place ( b = a*2 ) : '+tstr3+' seconds')
In [ ]:
##################################################################
# Example: Numpy array ordering
# In this example, we demonstrate the difference between
# row-major and column-major ordering of the array.
import numpy as np
dt = 'int32' # 4-byte integers
#Create a 1-D array with elements numbered 1 through 8
values = np.arange(1,9,1,dtype=dt)
#Reshape the 1-D array in two different ways
#Row-major ordering (default behavior; C-like).
#Values are loaded into row 0, then row 1, etc.
array2d_row_major = np.reshape(values,(4,2),order='C')
#Column-major ordering (Fortran-like)
#Values are loaded into column 0, then column 1, etc.
array2d_col_major = np.reshape(values,(4,2),order='F')
print('')
print('values: ')
print(values)
print('')
print('2-D reshape; row-major):')
print(array2d_row_major)
print('')
print('')
print('2-D reshape; column-major):')
print(array2d_col_major)
print('')
In [ ]:
####################################################################
# Example: Numpy array access patterns
#
# In general, try not to loop over array elements; use
# vectorized operations whenever possible. That said,
# sometimes loops cannot be avoided. When you HAVE to write a
# loop, the order in which you access the array can affect
# computation speed. If the array is row-major (default),
# it is most efficient to work row-by-row. If the array is
# column major, work column-by-column instead.
#
# In this example, we perform vector dot products with
# rows and columns of an arrays that are both row-major
# and column-major
from time import time
import numpy as np
nx = 64
nsizes = 6
ntests = 50
orders = ['C', 'F']
otitles = [' Row Major (C-style)', 'Column Major (Fortran-style)']
print('Vector-Matrix Multiplication Timings')
for i,o in enumerate(orders):
print(' ')
print('Array ordering: ', otitles[i])
for k in range(nsizes):
nxy = nx*(2**k)
nxys = 'nx: '+'{:4}'.format(nxy)
vec = np.zeros(nxy,dtype='float64', order=o)
mat = np.zeros((nxy,nxy),dtype='float64', order=o)
# first pass: vec = mat[:,j]
t0 = time()
for n in range(ntests):
for j in range(nxy):
vec2 = mat[:,j]
dsum = np.sum(vec2*vec)
t1 = time()
dt1 = (t1-t0)/ntests
# second pass: vec = mat[j,:]
t0 = time()
for n in range(ntests):
for j in range(nxy):
vec2 = mat[j,:]
dsum = np.sum(vec2*vec)
t1 = time()
dt2 = (t1-t0)/ntests
s2 = 1.0/dt2
s1 = 1.0/dt1
ds = (s2-s1)/s1
ds = ds*100
dratio = dt1/dt2
dss = '{:.4f}'.format(dratio)
dss = '(Column-wise Time)/(Row-wise time) = '+dss
print(nxys, ';', dss)
In [ ]:
##################################################################################
# Example: Writing & reading numpy arrays
#
# We can write array contents in unformatted binary using tofile.
# We can read from a file using fromfile.
#
# NOTE: Regardless of an array's ordering, tofile will
# ALWAYS write the array in row-major order.
# To see this, change the ordering of simple_array to 'F.'
# array2d[a,b] and array3d will remain unchanged.
import numpy as np
ofile = 'numpy_output.dat'
dt = 'int32' # 4-byte integers
simple_array = np.zeros((4,2),dtype=dt)
simple_array[0,0] = 1
simple_array[0,1] = 2
simple_array[1,0] = 3
simple_array[1,1] = 4
simple_array[2,0] = 5
simple_array[2,1] = 6
simple_array[3,0] = 7
simple_array[3,1] = 8
#Writing an array is easy - just specify a filename
simple_array.tofile(ofile)
#When reading unformatted binary, we specify
# (1) The filename
# (2) The datatype
# (3) The number of values to read
values = np.fromfile(ofile,dtype=dt,count=8)
#We can reshape the data as desired
array2da = np.reshape(values,(4,2))
array2db = np.reshape(values,(2,4))
array3d = np.reshape(values,(2,2,2))
print('')
print('values: ')
print(values)
print('')
print('2-D array; 4 rows, 2 columns):')
print(array2da)
print('')
print('')
print('2-D array; 2 rows, 4 columns):')
print(array2db)
print('')
print('')
print('3-D array:')
print(array3d)
print('')
In [ ]:
##########################################################################
# Example: Map & Reduction in Serial Python
#
# If we have a function designed to accept one argument,
# but want to evaluate it for multiple data values, we
# can use Python's "map" function as shorthand for a loop.
#
# If we have a function that we would like to call again
# and again, using results from one call in tandem with
# the subsequent call, we can use "reduce" function.
#
# The following program demonstrates the use of map & reduce.
from functools import reduce
def squarex(x):
return (x*x)
def add(x):
return (x+x)
def addtwo(x,y):
return x+y
def multtwo(x,y):
return x*y
arg_list = [1,3,9,27]
# Compute the square of 1,3,9, and 27
vals = map(squarex, [1,3,9,27])
for i,v in enumerate(vals):
str1 = '{:2d}'.format(arg_list[i])
str2 = '{:3d}'.format(v)
print(str1+' squared is '+str2+'.')
#Compute ((1+2)+3)
# -- first 1 and 2 are passed to addtwo
# -- next, (1+2 = 3) and 3 are passed to addtwo
# -- the final result is assigned to v2
res1 = reduce(addtwo,[1,2,3])
#Compute (((1*2) * 3) * 4)
# -- first 1 and 2 are passed to multtwo
# -- next, (1*2=2) and 3 are passed to addtwo
# -- next (2*3=6) and 4 are passed to addtwo
# -- the final result is assigned to res2
res2 = reduce(multtwo,[1,2,3,4])
print('')
print(' 1+2+3 = ', res1)
print('1*2*3*4 = ', res2)
In [ ]:
#///////////////////////////////////////////////////////////////
# Exercise 2:
# Rewrite the following code using the map and
# reduce functions.
from functools import reduce
n = 4
b = []
for i in range(1,n+1):
b.append(2*i)
print(b)
prod = 1
for i in range(n):
prod = prod*b[i]
print(prod)
In [ ]:
def times2(x):
return(x*2)
n = 4
idx = np.arange(1,n+1,1,dtype='int32')
b = map(times2, idx)
vals = [val for val in b] # N.B. !!!!!
print(vals)
# N.B. !!!!! b can only be iterated over once, after that it disappears!!!!
#for i, val in enumerate(b):
# print(i,val)
def product(x, y):
return(x * y)
# this is the same as: lambda x,y: x*y
prod = reduce(product, vals)
print("prod = ", prod)
In [ ]:
new = list(map(times2, idx))
print(new)
In [ ]: