# Numpy - multidimensional data arrays¶

## Author: Steven Christe¶

Numpy is the fundamental library for scientific computing in Python. It contains list-like objects that work like arrays, matrices, and data tables. It also provides linear algebra, Fourier transforms, random number generation, and tools for integrating C/C++ and Fortran code.

Python was never designed originally for scientific computing, and contains many high-level abstractions necessary to enable its enormously flexible object-oriented interface. In Python, storing most integers requires more than just 4-8 bytes. It also requires at least a couple pointers per-integer. Performing a calculation on two numbers requires one or two bytecode operations, each of which can take dozens of CPU instructions for each pass through the Python eval loop. And when it comes to looping and index operations of Python lists the situation is even worst. NumPy provides efficient objects for dealing with arrays of numbers and overloads operators so you don't have to ever (hopefully) use loops!

This ipython notebook file can be found at

or can be viewed (if you don't have ipython installed) here

This is heavily cribbed from J.R. Johansson ([email protected]) http://dml.riken.jp/~rob/

The full version of the original talk is available at http://github.com/jrjohansson/scientific-python-lectures.

In [158]:
%pylab inline

Populating the interactive namespace from numpy and matplotlib

WARNING: pylab import has clobbered these variables: ['indices']
%pylab --no-import-all prevents importing * from pylab and numpy

In [159]:
import numpy as np

In [2]:
import sys # only need to get python version
print("Python version: " + sys.version)
print("Numpy version: " + np.__version__)

Python version: 2.7.6 |Anaconda 1.8.0 (x86_64)| (default, Jan 10 2014, 11:23:15)
[GCC 4.0.1 (Apple Inc. build 5493)]
Numpy version: 1.8.1


In the numpy package the terminology used for vectors, matrices and higher-dimensional data sets is array.

## Creating numpy arrays¶

There are a number of ways to initialize new numpy arrays, for example from

• a Python list or tuples
• using functions that are dedicated to generating numpy arrays, such as arange, linspace, etc.

### From lists¶

For example, to create new vector and matrix arrays from Python lists we can use the numpy.array function.

In [3]:
# a vector: the argument to the array function is a Python list
v = np.array([1,2,3,4])

v

Out[3]:
array([1, 2, 3, 4])
In [148]:
# a matrix: the argument to the array function is a nested Python list
M = np.array([[1, 6], [3, 4]])

M

Out[148]:
array([[1, 6],
[3, 4]])

The v and M objects are both of the type ndarray that the numpy module provides.

In [149]:
type(v), type(M)

Out[149]:
(numpy.ndarray, numpy.ndarray)

The difference between the v and M arrays is only their shapes. We can get information about the shape of an array by using the ndarray.shape property.

In [150]:
v.shape

Out[150]:
(4,)
In [7]:
M.shape

Out[7]:
(2, 2)

The number of elements in the array is available through the ndarray.size property:

In [8]:
M.size

Out[8]:
4

Equivalently, we could use the function numpy.shape and numpy.size

In [9]:
np.shape(M)

Out[9]:
(2, 2)
In [10]:
np.size(M)

Out[10]:
4

Features of a numpy.ndarray

• Numpy arrays are statically typed and homogeneous. The type of the elements is determined when array is created.
• Numpy arrays are memory efficient.
• Because of the static typing, fast implementation of mathematical functions such as multiplication and addition of numpy arrays can be implemented in a compiled language (C and Fortran is used).

Using the dtype (data type) property of an ndarray, we can see what type the data of an array has:

In [11]:
M.dtype

Out[11]:
dtype('int64')

We get an error if we try to assign a value of the wrong type to an element in a numpy array:

If we want, we can explicitly define the type of the array data when we create it, using the dtype keyword argument:

In [12]:
M = np.array([[1, 2], [3, 4]], dtype=np.float16)

M

Out[12]:
array([[ 1.,  2.],
[ 3.,  4.]], dtype=float16)

Common type that can be used with dtype are: int, float, complex, bool, object, etc.

### Array-generating functions¶

For larger arrays it is inpractical to initialize the data manually, using explicit python lists. Instead we can use one of the many functions in numpy that generates arrays of different forms. Some of the more common are:

#### arange¶

In [13]:
# create a range

x = np.arange(0, 10, 1) # arguments: start, stop, step

x

Out[13]:
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
In [14]:
x = np.arange(-1, 1, 0.1)

x

Out[14]:
array([ -1.00000000e+00,  -9.00000000e-01,  -8.00000000e-01,
-7.00000000e-01,  -6.00000000e-01,  -5.00000000e-01,
-4.00000000e-01,  -3.00000000e-01,  -2.00000000e-01,
-1.00000000e-01,  -2.22044605e-16,   1.00000000e-01,
2.00000000e-01,   3.00000000e-01,   4.00000000e-01,
5.00000000e-01,   6.00000000e-01,   7.00000000e-01,
8.00000000e-01,   9.00000000e-01])

#### linspace and logspace¶

In [152]:
# using linspace, both end points ARE included
np.linspace(0, 10, 25)

Out[152]:
array([  0.        ,   0.41666667,   0.83333333,   1.25      ,
1.66666667,   2.08333333,   2.5       ,   2.91666667,
3.33333333,   3.75      ,   4.16666667,   4.58333333,
5.        ,   5.41666667,   5.83333333,   6.25      ,
6.66666667,   7.08333333,   7.5       ,   7.91666667,
8.33333333,   8.75      ,   9.16666667,   9.58333333,  10.        ])
In [16]:
np.logspace(0, 10, 10, base=e)

Out[16]:
array([  1.00000000e+00,   3.03773178e+00,   9.22781435e+00,
2.80316249e+01,   8.51525577e+01,   2.58670631e+02,
7.85771994e+02,   2.38696456e+03,   7.25095809e+03,
2.20264658e+04])

#### mgrid¶

In [17]:
x, y = np.mgrid[0:5, 0:5] # similar to meshgrid in MATLAB

In [18]:
x

Out[18]:
array([[0, 0, 0, 0, 0],
[1, 1, 1, 1, 1],
[2, 2, 2, 2, 2],
[3, 3, 3, 3, 3],
[4, 4, 4, 4, 4]])
In [19]:
y

Out[19]:
array([[0, 1, 2, 3, 4],
[0, 1, 2, 3, 4],
[0, 1, 2, 3, 4],
[0, 1, 2, 3, 4],
[0, 1, 2, 3, 4]])

#### random data¶

In [20]:
from numpy import random

In [21]:
# uniform random numbers in [0,1]
random.rand(5,5)

Out[21]:
array([[ 0.32846187,  0.42213273,  0.01968231,  0.21261055,  0.51482732],
[ 0.08183575,  0.00824206,  0.62834203,  0.77968281,  0.126293  ],
[ 0.38743163,  0.86545303,  0.40312758,  0.8053914 ,  0.48047428],
[ 0.61395831,  0.89044911,  0.14621425,  0.48195724,  0.38045612],
[ 0.20086805,  0.27765159,  0.08753675,  0.18159165,  0.64932571]])
In [22]:
# standard normal distributed random numbers
random.randn(5,5)

Out[22]:
array([[-0.78404982, -0.57134026,  1.49107161, -1.26857542,  0.44456653],
[-0.02629558,  0.08901156, -1.08758579,  0.22911653, -1.98229464],
[ 0.33490693,  0.65502718,  0.23791072, -0.22042591, -0.01485002],
[ 1.49117271, -0.76377386,  0.54285327, -0.78657823,  1.72420386],
[ 1.77250365, -1.02862661,  0.07318933, -0.06165841, -1.93811428]])

#### zeros and ones¶

In [23]:
np.zeros((3,3))

Out[23]:
array([[ 0.,  0.,  0.],
[ 0.,  0.,  0.],
[ 0.,  0.,  0.]])
In [24]:
np.ones((3,3))

Out[24]:
array([[ 1.,  1.,  1.],
[ 1.,  1.,  1.],
[ 1.,  1.,  1.]])

## More properties of the numpy arrays¶

In [25]:
M.itemsize # bytes per element

Out[25]:
2
In [26]:
M.nbytes # number of bytes

Out[26]:
8
In [27]:
M.ndim # number of dimensions

Out[27]:
2

## Manipulating arrays¶

### Indexing¶

We can index elements in an array using the square bracket and indices:

In [154]:
# v is a vector, and has only one dimension, taking one index
v[0]

Out[154]:
1
In [29]:
# M is a matrix, or a 2 dimensional array, taking two indices
M[1,1]

Out[29]:
4.0

The same thing can be achieved with using : instead of an index:

In [30]:
M[1,:] # row 1

Out[30]:
array([ 3.,  4.], dtype=float16)
In [31]:
M[:,1] # column 1

Out[31]:
array([ 2.,  4.], dtype=float16)

We can assign new values to elements in an array using indexing:

In [32]:
M[0,0] = 1

In [33]:
M

Out[33]:
array([[ 1.,  2.],
[ 3.,  4.]], dtype=float16)
In [34]:
# also works for rows and columns
M[1,:] = 0
M[:,1] = -1

In [35]:
M

Out[35]:
array([[ 1., -1.],
[ 0., -1.]], dtype=float16)

### Index slicing¶

Index slicing is the technical name for the syntax M[lower:upper:step] to extract part of an array:

In [36]:
A = np.array([1,2,3,4,5])
A

Out[36]:
array([1, 2, 3, 4, 5])
In [37]:
A[1:3]

Out[37]:
array([2, 3])

Array slices are mutable: if they are assigned a new value the original array from which the slice was extracted is modified:

In [38]:
A[1:3] = [-2,-3]

A

Out[38]:
array([ 1, -2, -3,  4,  5])

We can omit any of the three parameters in M[lower:upper:step]:

In [39]:
A[::] # lower, upper, step all take the default values

Out[39]:
array([ 1, -2, -3,  4,  5])
In [40]:
A[::2] # step is 2, lower and upper defaults to the beginning and end of the array

Out[40]:
array([ 1, -3,  5])
In [41]:
A[:3] # first three elements

Out[41]:
array([ 1, -2, -3])
In [42]:
A[3:] # elements from index 3

Out[42]:
array([4, 5])

Negative indices counts from the end of the array (positive index from the begining):

In [43]:
A = np.array([1,2,3,4,5])

In [44]:
A[-1] # the last element in the array

Out[44]:
5
In [45]:
A[-3:] # the last three elements

Out[45]:
array([3, 4, 5])

Index slicing works exactly the same way for multidimensional arrays:

In [46]:
A = np.array([[n+m*10 for n in range(5)] for m in range(5)])

A

Out[46]:
array([[ 0,  1,  2,  3,  4],
[10, 11, 12, 13, 14],
[20, 21, 22, 23, 24],
[30, 31, 32, 33, 34],
[40, 41, 42, 43, 44]])
In [47]:
# a block from the original array
A[1:4, 1:4]

Out[47]:
array([[11, 12, 13],
[21, 22, 23],
[31, 32, 33]])
In [48]:
# strides
A[::2, ::2]

Out[48]:
array([[ 0,  2,  4],
[20, 22, 24],
[40, 42, 44]])
In [157]:
A[1]

Out[157]:
array([3, 4])

### Fancy indexing¶

Fancy indexing is the name for when an array or list is used in-place of an index:

In [49]:
row_indices = [1, 2, 3]
A[row_indices]

Out[49]:
array([[10, 11, 12, 13, 14],
[20, 21, 22, 23, 24],
[30, 31, 32, 33, 34]])
In [50]:
col_indices = [1, 2, -1] # remember, index -1 means the last element
A[row_indices, col_indices]

Out[50]:
array([11, 22, 34])

We can also index masks: If the index mask is an Numpy array of with data type bool, then an element is selected (True) or not (False) depending on the value of the index mask at the position each element:

In [51]:
B = array([n for n in range(5)])
B

Out[51]:
array([0, 1, 2, 3, 4])
In [52]:
row_mask = array([True, False, True, False, False])

Out[52]:
array([0, 2])
In [53]:
# same thing

Out[53]:
array([0, 2])

This feature is very useful to conditionally select elements from an array, using for example comparison operators:

In [54]:
x = np.arange(0, 10, 0.5)
x

Out[54]:
array([ 0. ,  0.5,  1. ,  1.5,  2. ,  2.5,  3. ,  3.5,  4. ,  4.5,  5. ,
5.5,  6. ,  6.5,  7. ,  7.5,  8. ,  8.5,  9. ,  9.5])
In [55]:
mask = (5 < x) * (x < 7.5)


Out[55]:
array([False, False, False, False, False, False, False, False, False,
False, False,  True,  True,  True,  True, False, False, False,
False, False], dtype=bool)
In [56]:
x[mask]

Out[56]:
array([ 5.5,  6. ,  6.5,  7. ])

## Functions for extracting data from arrays and creating arrays¶

### where¶

The index mask can be converted to position index using the where function

In [57]:
indices = np.where(mask)

indices

Out[57]:
(array([11, 12, 13, 14]),)
In [58]:
x[indices] # this indexing is equivalent to the fancy indexing x[mask]

Out[58]:
array([ 5.5,  6. ,  6.5,  7. ])

### diag¶

With the diag function we can also extract the diagonal and subdiagonals of an array:

In [59]:
np.diag(A)

Out[59]:
array([ 0, 11, 22, 33, 44])
In [60]:
np.diag(A, -1)

Out[60]:
array([10, 21, 32, 43])

### take¶

The take function is similar to fancy indexing described above:

In [61]:
v2 = np.arange(-3,3)
v2

Out[61]:
array([-3, -2, -1,  0,  1,  2])
In [62]:
row_indices = [1, 3, 5]
v2[row_indices] # fancy indexing

Out[62]:
array([-2,  0,  2])
In [63]:
v2.take(row_indices)

Out[63]:
array([-2,  0,  2])

But take also works on lists and other objects:

In [64]:
np.take([-3, -2, -1,  0,  1,  2], row_indices)

Out[64]:
array([-2,  0,  2])

### choose¶

Constructs and array by picking elements form several arrays:

In [65]:
which = [1, 0, 1, 0]
choices = [[-2,-2,-2,-2], [5,5,5,5]]

np.choose(which, choices)

Out[65]:
array([ 5, -2,  5, -2])

## Linear algebra¶

Vectorizing code is the key to writing efficient numerical calculation with Python/Numpy. That means that as much as possible of a program should be formulated in terms of matrix and vector operations, like matrix-matrix multiplication.

### Scalar-array operations¶

We can use the usual arithmetic operators to multiply, add, subtract, and divide arrays with scalar numbers.

In [66]:
v1 = np.arange(0, 5)

In [67]:
v1 * 2

Out[67]:
array([0, 2, 4, 6, 8])
In [68]:
v1 + 2

Out[68]:
array([2, 3, 4, 5, 6])
In [69]:
A * 2, A + 2

Out[69]:
(array([[ 0,  2,  4,  6,  8],
[20, 22, 24, 26, 28],
[40, 42, 44, 46, 48],
[60, 62, 64, 66, 68],
[80, 82, 84, 86, 88]]),
array([[ 2,  3,  4,  5,  6],
[12, 13, 14, 15, 16],
[22, 23, 24, 25, 26],
[32, 33, 34, 35, 36],
[42, 43, 44, 45, 46]]))

### Element-wise array-array operations¶

When we add, subtract, multiply and divide arrays with each other, the default behaviour is element-wise operations:

In [70]:
A * A # element-wise multiplication

Out[70]:
array([[   0,    1,    4,    9,   16],
[ 100,  121,  144,  169,  196],
[ 400,  441,  484,  529,  576],
[ 900,  961, 1024, 1089, 1156],
[1600, 1681, 1764, 1849, 1936]])
In [71]:
v1 * v1

Out[71]:
array([ 0,  1,  4,  9, 16])

If we multiply arrays with compatible shapes, we get an element-wise multiplication of each row:

In [72]:
A.shape, v1.shape

Out[72]:
((5, 5), (5,))
In [73]:
A * v1

Out[73]:
array([[  0,   1,   4,   9,  16],
[  0,  11,  24,  39,  56],
[  0,  21,  44,  69,  96],
[  0,  31,  64,  99, 136],
[  0,  41,  84, 129, 176]])

### Matrix algebra¶

What about matrix mutiplication? There are two ways. We can either use the dot function, which applies a matrix-matrix, matrix-vector, or inner vector multiplication to its two arguments:

In [74]:
np.dot(A, A)

Out[74]:
array([[ 300,  310,  320,  330,  340],
[1300, 1360, 1420, 1480, 1540],
[2300, 2410, 2520, 2630, 2740],
[3300, 3460, 3620, 3780, 3940],
[4300, 4510, 4720, 4930, 5140]])
In [75]:
np.dot(A, v1)

Out[75]:
array([ 30, 130, 230, 330, 430])
In [76]:
np.dot(v1, v1)

Out[76]:
30

Alternatively, we can cast the array objects to the type matrix. This changes the behavior of the standard arithmetic operators +, -, * to use matrix algebra.

In [77]:
M = np.matrix(A)
v = np.matrix(v1).T # make it a column vector

In [78]:
v

Out[78]:
matrix([[0],
[1],
[2],
[3],
[4]])
In [79]:
M * M

Out[79]:
matrix([[ 300,  310,  320,  330,  340],
[1300, 1360, 1420, 1480, 1540],
[2300, 2410, 2520, 2630, 2740],
[3300, 3460, 3620, 3780, 3940],
[4300, 4510, 4720, 4930, 5140]])
In [80]:
M * v

Out[80]:
matrix([[ 30],
[130],
[230],
[330],
[430]])
In [81]:
# inner product
v.T * v

Out[81]:
matrix([[30]])
In [82]:
# with matrix objects, standard matrix algebra applies
v + M*v

Out[82]:
matrix([[ 30],
[131],
[232],
[333],
[434]])

If we try to add, subtract or multiply objects with incomplatible shapes we get an error:

In [83]:
v = np.matrix([1,2,3,4,5,6]).T

In [84]:
shape(M), shape(v)

Out[84]:
((5, 5), (6, 1))
In [85]:
M * v

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
----> 1 M * v

/Users/schriste/anaconda/lib/python2.7/site-packages/numpy/matrixlib/defmatrix.pyc in __mul__(self, other)
328         if isinstance(other,(N.ndarray, list, tuple)) :
329             # This promotes 1-D vectors to row vectors
--> 330             return N.dot(self, asmatrix(other))
331         if isscalar(other) or not hasattr(other, '__rmul__') :
332             return N.dot(self, other)

ValueError: objects are not aligned

See also the related functions: inner, outer, cross, kron, tensordot. Try for example help(kron).

### Array/Matrix transformations¶

Above we have used the .T to transpose the matrix object v. We could also have used the transpose function to accomplish the same thing.

Other mathematical functions that transforms matrix objects are:

In [ ]:
C = np.matrix([[1j, 2j], [3j, 4j]])
C

In [86]:
np.conjugate(C)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-86-934f2dd5e2ef> in <module>()
----> 1 np.conjugate(C)

NameError: name 'C' is not defined

Hermitian conjugate: transpose + conjugate

In [87]:
C.H

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-87-bdd18b98227d> in <module>()
----> 1 C.H

NameError: name 'C' is not defined

We can extract the real and imaginary parts of complex-valued arrays using real and imag:

In [88]:
np.real(C) # same as: C.real

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-88-b672483d4bc4> in <module>()
----> 1 np.real(C) # same as: C.real

NameError: name 'C' is not defined
In [89]:
np.imag(C) # same as: C.imag

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-89-1ed1d19d3c2c> in <module>()
----> 1 np.imag(C) # same as: C.imag

NameError: name 'C' is not defined

Or the complex argument and absolute value

In [90]:
np.angle(C+1) # heads up MATLAB Users, angle is used instead of arg

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-90-3775d8ef0c87> in <module>()
----> 1 np.angle(C+1) # heads up MATLAB Users, angle is used instead of arg

NameError: name 'C' is not defined
In [91]:
np.abs(C)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-91-56f6751545e7> in <module>()
----> 1 np.abs(C)

NameError: name 'C' is not defined

### Matrix computations¶

#### Inverse¶

In [92]:
inv(C) # equivalent to C.I

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-92-0b30bd464f9b> in <module>()
----> 1 inv(C) # equivalent to C.I

NameError: name 'C' is not defined
In [93]:
C.I * C

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-93-33ab9fba3975> in <module>()
----> 1 C.I * C

NameError: name 'C' is not defined

#### Determinant¶

In [ ]:
det(C)

In [ ]:
det(C.I)


#### sum, prod, and trace¶

In [ ]:
d = arange(0, 10)
d

In [ ]:
# sum up all elements
sum(d)

In [94]:
# product of all elements
prod(d+1)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-94-e8cd400f9a6b> in <module>()
1 # product of all elements
----> 2 prod(d+1)

NameError: name 'd' is not defined
In [95]:
# cummulative sum
cumsum(d)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-95-242a18f7884c> in <module>()
1 # cummulative sum
----> 2 cumsum(d)

NameError: name 'd' is not defined
In [96]:
# cummulative product
cumprod(d+1)

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-96-bcb7dc412ffa> in <module>()
1 # cummulative product
----> 2 cumprod(d+1)

NameError: name 'd' is not defined
In [97]:
# same as: diag(A).sum()
trace(A)

Out[97]:
110

### Calculations with higher-dimensional data¶

When functions such as min, max, etc., is applied to a multidimensional arrays, it is sometimes useful to apply the calculation to the entire array, and sometimes only on a row or column basis. Using the axis argument we can specify how these functions should behave:

In [98]:
m = rand(3,3)
m

Out[98]:
array([[ 0.45974672,  0.90724345,  0.23833268],
[ 0.90236963,  0.51838115,  0.64488138],
[ 0.49125775,  0.50524105,  0.3813743 ]])
In [99]:
# global max
m.max()

Out[99]:
0.90724344539482737
In [100]:
# max in each column
m.max(axis=0)

Out[100]:
array([ 0.90236963,  0.90724345,  0.64488138])
In [101]:
# max in each row
m.max(axis=1)

Out[101]:
array([ 0.90724345,  0.90236963,  0.50524105])

Many other functions and methods in the array and matrix classes accept the same (optional) axis keyword argument.

## Reshaping, resizing and stacking arrays¶

The shape of an Numpy array can be modified without copying the underlaying data, which makes it a fast operation even for large arrays.

In [102]:
A

Out[102]:
array([[ 0,  1,  2,  3,  4],
[10, 11, 12, 13, 14],
[20, 21, 22, 23, 24],
[30, 31, 32, 33, 34],
[40, 41, 42, 43, 44]])
In [103]:
n, m = A.shape

In [104]:
B = A.reshape((1,n*m))
B

Out[104]:
array([[ 0,  1,  2,  3,  4, 10, 11, 12, 13, 14, 20, 21, 22, 23, 24, 30, 31,
32, 33, 34, 40, 41, 42, 43, 44]])
In [105]:
B[0,0:5] = 5 # modify the array

B

Out[105]:
array([[ 5,  5,  5,  5,  5, 10, 11, 12, 13, 14, 20, 21, 22, 23, 24, 30, 31,
32, 33, 34, 40, 41, 42, 43, 44]])
In [106]:
A # and the original variable is also changed. B is only a different view of the same data

Out[106]:
array([[ 5,  5,  5,  5,  5],
[10, 11, 12, 13, 14],
[20, 21, 22, 23, 24],
[30, 31, 32, 33, 34],
[40, 41, 42, 43, 44]])

We can also use the function flatten to make a higher-dimensional array into a vector. But this function create a copy of the data.

In [107]:
B = A.flatten()

B

Out[107]:
array([ 5,  5,  5,  5,  5, 10, 11, 12, 13, 14, 20, 21, 22, 23, 24, 30, 31,
32, 33, 34, 40, 41, 42, 43, 44])
In [108]:
B[0:5] = 10

B

Out[108]:
array([10, 10, 10, 10, 10, 10, 11, 12, 13, 14, 20, 21, 22, 23, 24, 30, 31,
32, 33, 34, 40, 41, 42, 43, 44])
In [109]:
A # now A has not changed, because B's data is a copy of A's, not refering to the same data

Out[109]:
array([[ 5,  5,  5,  5,  5],
[10, 11, 12, 13, 14],
[20, 21, 22, 23, 24],
[30, 31, 32, 33, 34],
[40, 41, 42, 43, 44]])

## Adding a new dimension: newaxis¶

With newaxis, we can insert new dimensions in an array, for example converting a vector to a column or row matrix:

In [110]:
v = np.array([1,2,3])

In [111]:
np.shape(v)

Out[111]:
(3,)
In [112]:
# make a column matrix of the vector v
v[:, newaxis]

Out[112]:
array([[1],
[2],
[3]])
In [113]:
# column matrix
v[:,newaxis].shape

Out[113]:
(3, 1)
In [114]:
# row matrix
v[newaxis,:].shape

Out[114]:
(1, 3)

## Stacking and repeating arrays¶

Using function repeat, tile, vstack, hstack, and concatenate we can create larger vectors and matrices from smaller ones:

### tile and repeat¶

In [115]:
a = np.array([[1, 2], [3, 4]])

In [116]:
# repeat each element 3 times
np.repeat(a, 3)

Out[116]:
array([1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 4])
In [117]:
# tile the matrix 3 times
np.tile(a, 3)

Out[117]:
array([[1, 2, 1, 2, 1, 2],
[3, 4, 3, 4, 3, 4]])

### concatenate¶

In [118]:
b = np.array([[5, 6]])

In [119]:
np.concatenate((a, b), axis=0)

Out[119]:
array([[1, 2],
[3, 4],
[5, 6]])
In [120]:
np.concatenate((a, b.T), axis=1)

Out[120]:
array([[1, 2, 5],
[3, 4, 6]])

### hstack and vstack¶

In [121]:
np.vstack((a,b))

Out[121]:
array([[1, 2],
[3, 4],
[5, 6]])
In [122]:
np.hstack((a,b.T))

Out[122]:
array([[1, 2, 5],
[3, 4, 6]])

## Copy and "deep copy"¶

To achieve high performance, assignments in Python usually do not copy the underlaying objects. This is important for example when objects are passed between functions, to avoid an excessive amount of memory copying when it is not necessary (techincal term: pass by reference).

In [123]:
A = np.array([[1, 2], [3, 4]])

A

Out[123]:
array([[1, 2],
[3, 4]])
In [124]:
# now B is referring to the same array data as A
B = A

In [125]:
# changing B affects A
B[0,0] = 10

B

Out[125]:
array([[10,  2],
[ 3,  4]])
In [126]:
A

Out[126]:
array([[10,  2],
[ 3,  4]])

If we want to avoid this behavior, so that when we get a new completely independent object B copied from A, then we need to do a so-called "deep copy" using the function copy:

In [127]:
B = copy(A)

In [128]:
# now, if we modify B, A is not affected
B[0,0] = -5

B

Out[128]:
array([[-5,  2],
[ 3,  4]])
In [129]:
A

Out[129]:
array([[10,  2],
[ 3,  4]])

## Iterating over array elements¶

Generally, we want to avoid iterating over the elements of arrays whenever we can (at all costs). The reason is that in a interpreted language like Python (or MATLAB), iterations are really slow compared to vectorized operations.

However, sometimes iterations are unavoidable. For such cases, the Python for loop is the most convenient way to iterate over an array:

In [130]:
v = np.array([1,2,3,4])

for element in v:
print(element)

1
2
3
4

In [131]:
M = np.array([[1,2], [3,4]])

for row in M:
print("row", row)

for element in row:
print(element)

('row', array([1, 2]))
1
2
('row', array([3, 4]))
3
4


When we need to iterate over each element of an array and modify its elements, it is convenient to use the enumerate function to obtain both the element and its index in the for loop:

In [132]:
for row_idx, row in enumerate(M):
print("row_idx", row_idx, "row", row)

for col_idx, element in enumerate(row):
print("col_idx", col_idx, "element", element)

# update the matrix M: square each element
M[row_idx, col_idx] = element ** 2

('row_idx', 0, 'row', array([1, 2]))
('col_idx', 0, 'element', 1)
('col_idx', 1, 'element', 2)
('row_idx', 1, 'row', array([3, 4]))
('col_idx', 0, 'element', 3)
('col_idx', 1, 'element', 4)

In [133]:
# each element in M is now squared
M

Out[133]:
array([[ 1,  4],
[ 9, 16]])

## Vectorizing functions¶

As mentioned several times by now, to get good performance we should try to avoid looping over elements in our vectors and matrices, and instead use vectorized algorithms. The first step in converting a scalar algorithm to a vectorized algorithm is to make sure that the functions we write work with vector inputs.

In [134]:
def Theta(x):
"""
Scalar implemenation of the Heaviside step function.
"""
if x >= 0:
return 1
else:
return 0

In [135]:
Theta(np.array([-3,-2,-1,0,1,2,3]))

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-135-d55419725688> in <module>()
----> 1 Theta(np.array([-3,-2,-1,0,1,2,3]))

<ipython-input-134-9a0cb13d93d4> in Theta(x)
3     Scalar implemenation of the Heaviside step function.
4     """
----> 5     if x >= 0:
6         return 1
7     else:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

OK, that didn't work because we didn't write the Theta function so that it can handle with vector input...

To get a vectorized version of Theta we can use the Numpy function vectorize. In many cases it can automatically vectorize a function:

In [136]:
Theta_vec = np.vectorize(Theta)

In [137]:
Theta_vec(np.array([-3,-2,-1,0,1,2,3]))

Out[137]:
array([0, 0, 0, 1, 1, 1, 1])

We can also implement the function to accept vector input from the beginning (requires more effort but might give better performance):

In [138]:
def Theta(x):
"""
Vector-aware implemenation of the Heaviside step function.
"""
return 1 * (x >= 0)

In [139]:
Theta(np.array([-3,-2,-1,0,1,2,3]))

Out[139]:
array([0, 0, 0, 1, 1, 1, 1])
In [140]:
# still works for scalars as well
Theta(-1.2), Theta(2.6)

Out[140]:
(0, 1)

## Using arrays in conditions¶

When using arrays in conditions in for example if statements and other boolean expressions, one need to use one of any or all, which requires that any or all elements in the array evalutes to True:

In [141]:
M

Out[141]:
array([[ 1,  4],
[ 9, 16]])
In [142]:
if (M > 5).any():
print("at least one element in M is larger than 5")
else:
print("no element in M is larger than 5")

at least one element in M is larger than 5

In [143]:
if (M > 5).all():
print("all elements in M are larger than 5")
else:
print("all elements in M are not larger than 5")

all elements in M are not larger than 5


## Type casting¶

Since Numpy arrays are statically typed, the type of an array does not change once created. But we can explicitly cast an array of some type to another using the astype functions (see also the similar asarray function). This always create a new array of new type:

In [144]:
M.dtype

Out[144]:
dtype('int64')
In [145]:
M2 = M.astype(float)

M2

Out[145]:
array([[  1.,   4.],
[  9.,  16.]])
In [146]:
M2.dtype

Out[146]:
dtype('float64')
In [147]:
M3 = M.astype(bool)

M3

Out[147]:
array([[ True,  True],
[ True,  True]], dtype=bool)

# SciPy¶

In [147]: