Efficient Numerical Code With Numpy Vectorization

C and fortran paradigm: For loops

• Python loops have python overhead
• Numpy operations are lower level

Ex: adding two arrays

n = 100
x = np.random.rand(n)
y = np.random.rand(n)
f = []
for i in range(100):
    f.append(x[i] + y[i])

Numpy paradigm: Array operations

• Useful when performing the same operation to many numerical pieces of data.
• Ex: Adding all the elements from two lists of numbers

x = np.rand(100)
y = np.rand(100)
f = x + y

Simple operations

• Addition two lists of numbers
• numpy array functions
  • np.sum(), np.cumsum(), np.exp()

Some Numpy version have additional optimizations:

Under the hood, complex operations will automatically be run on multiple processors.

• SIMD
• Intel MKL
• To get these benefits in anaconda:

conda install mkl

Some Types of Numerical Data Structures

• Grid data - possibly regularly spaced
• Point data

Diffusion Example

• Finite Difference Examples: diffusion.py

  • Grid data -http://en.wikipedia.org/wiki/Finite_difference_method#Example:_The_heat_equation
  • Equation: $$ u_j^{n+1} = (1-2r) u_j^n + r u_{n-1} + r u_{j+1}^n $$ Where $j$ is the spatial index, and $n$ is the time index

• 1-D

  • For loop
  • array version
  • Show 2-d plot of all results to compare to gaussian motion
• 2-D
  • For loop: diff2d_loop()
  • array version: diff2d_vec()
• Draw overlapping boxes to show how slice indexing should work
• %timeit and %lprun useful for measuring timing
• lineprofiler: https://github.com/rkern/line_profiler

Monte Carlo Example

• Gaussian motion
• List of Particles

Data reshaping

Often data is not in correct shape for easily performing array manipulations

• Useful array member functions:
  • .T
  • .ravel()
  • .reshape(newShape)
• Array views
• Rotating 2-D grid of points
• Rotating 2-D list of points from gaussian motion