Script 3.1. of the online course "An introduction to agent-based modeling with Python" by Claudius Gräbner

For the course homepage and links to the acompanying videos (in German) see: http://claudius-graebner.com/introabmdt.html

For the Enlish version (currently without videos) see: http://claudius-graebner.com/introabmen.html

Last updated July 18 2019

One of the most commonly used packages in Python is `numpy`

.
It stands for *Numeric Python* since it originates from the desire to use Python for numerical analysis.

Basically, the main features `numpy`

provides are the following:

- A new data type called
`array`

. It is well suited to represent vectors, and matrices. - A submodule that allows the creation and handling of random numbers
- All kinds of mathematical routines in the area of numerical and stochastic modeling

All data types and functions provided by `numpy`

are very *efficient*.
This meant they require relatively few computational ressources and they are fast.
Why is that?
Remember that we mentioned Python to be a relatively slow language.
That is true, but most demanding functions of `numpy`

are not written in Python, but in C or Fortran.
These are very demanding languages, which require much effort and time from programmers.
Yet, they are extremely fast.
Thereby, `numpy`

manages to exploit their speed, but to maintain the ease of the Python syntax.

For us, `numpy`

becomes important because many mathematical operations in the context of ABM involve:

- The use of matrix algebra
- The use of random processes

For both tasks, `numpy`

provides us with everything we need, so we will use it extensively when building ABMs.

Before we get started we need to import `numpy`

.
By convention, this is done using the alias `np`

:

In [2]:

```
import numpy as np
```

Arrays are the basic data type of `numpy`

(this means: the basic class; many other class objects, such as matrices, *inherit* from the class `ndarray`

):

In [2]:

```
a_1 = np.array([1,2])
type(a_1)
```

Out[2]:

In many ways, arrays behave similar to lists, yet there are also some important differences:

- Arrays can only contain elements that are of the same type
- Arrays are the only data type supported by more advanced mathematical functions
- Most operations on arrays are much faster than on lists

There are some more specific differences underlying this (such that arrays have constant memory size), and that explain the superiority of arrays over lists, yet we do not want to bother with them right now.

For you, point 1 above is probably the practically most relevant distinction for the moment.

Lets see how we can create arrays. There are several options.

First, we can take an existing container, e.g. a `list`

or a `tuple`

and transform it into an array:

In [9]:

```
l_1 = [1,2,3,4]
print(type(l_1))
t_1 = (5,6,7,8)
print(type(t_1))
a_1 = np.array(l_1)
print(type(a_1))
a_2 = np.array(t_1)
print(type(a_2))
```

Note that building arrays from `sets`

is a bit more problematic and requires more steps, which we do not cover here.

Indexing works for arrays just as for tuples and lists:

In [15]:

```
print(a_1)
a_1[0]
```

Out[15]:

It is also straightforward to create multidimensional arrays, such as matrices.

For example, to create the following matrix:

$\boldsymbol{A}=\begin{pmatrix} 1 & 3 \\ 6 & 8 \end{pmatrix}$

Just type:

In [16]:

```
a_3 = np.array([[1,3], [6,8]])
print(a_3)
```

At this point it is important to look at the very precise syntax of `np.array`

:
it **only takes one argument**: a container.
In case you want to specify a multidimensional array, this needs to be encolosed in brackets:

In [19]:

```
a_4 = np.array(([1,2], [3,4], [5,6]) )
a_5 = np.array([[1,2], [3,4], [5,6]] ) # This is equivalent
print(a_4)
print(a_5)
```

The following, however, will not work:

In [20]:

```
a_6 = np.array([1,2], [3,4], [5,6])
```

The reasons why the brackets are needed is the following:
`np.array`

creates an array from the first object it receives.

If you want to create a matrix you must provide more than one containes (one for each row).
So the two (or more containers) must be summarized into one arbitrary container before getting passsed to `np.array`

as its first argument.

If there were no brackts, `np.array`

would receive more than one argument, and it does not know how to handle this!

Aside from creating arrays from tuples, there are a number of convenience functions to create typical arrays.

If you want an array full of zeros you can use `np.zeros`

with the desired shape:

In [5]:

```
print( np.zeros( (3,4) ) ) # for a 3 * 4 matrix
print("\n") # print new line
print( np.zeros(4) ) # for a 4-element vector
```

Exactly the same logic applies to ones:

In [7]:

```
print( np.ones( (3,4) ) ) # for a 3 * 4 matrix
print("\n") # print new line
print( np.ones( (4) ) ) # for a 4-element vector
```

If you want to get a sequence of numbers from $x_0$ to $x_n$, you can use the `np.arrange`

function:
it takes $x_0$ and $x_n$ as the first two arguments, and a third argument that specifies the step length.
Note that $x_n$ will not be included by default:

In [67]:

```
np.arange(10, 30, 5 ) # returns a sequence from 10 to 30 in 5 steps, excluding 30
```

Out[67]:

If you want to specify the number of elements in the array, and let Python determine the step length, use `np.linspace`

instead:

In [68]:

```
np.linspace(5, 12, 11 ) # returns an evenly-spaced array with 11 elements from 5 to 12
```

Out[68]:

These functions are useful when you want to evaluate functions at many steps:

In [69]:

```
x = np.linspace(-5, 5, 10)
y = np.sin(x)
y
```

Out[69]:

Note that `np.linspace`

is recommended when using `floats`

.

We learned about attributes in the previous sessions. We know that there are two main types of attributes:

- Attribute properties
- Attribute functions

Typical properties of an array are, among others:

In [70]:

```
a = np.array(([1,2], [3,4]))
a.shape
```

Out[70]:

Thus, the `shape`

gives you rows and columns of the array. Since `a`

in the example corresponds to the following matrix:

$\boldsymbol{A}=\begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix}$

it is of shape $2\times2$.

Similar, the `ndim`

attribute gives us the dimension of the array.
For the case of `a`

this is obviously 2:

In [22]:

```
a.ndim
```

Out[22]:

To get the number of elements in the array we can check its `size`

:

In [23]:

```
a.size
```

Out[23]:

Typical attribute functions, i.e. methods of the array instantiations, are among others, the sum of the columns:

In [24]:

```
a.sum(axis=0)
```

Out[24]:

or the sum of the rows:

In [25]:

```
a.sum(axis=1)
```

Out[25]:

Similarly, we can get the minimum of the rows and colums:

In [27]:

```
a.min(axis=1) # for the rows
```

Out[27]:

In [28]:

```
a.min(axis=0) # for the columns
```

Out[28]:

Finally, the cumulative sum of an array is given by the following:

In [29]:

```
a.cumsum(axis=0) # for the columns
```

Out[29]:

In [30]:

```
a.cumsum(axis=1) # for the rows
```

Out[30]:

I leave it to you to figure out what the cumulative sum is;)

We will learn more attributes in the next section.

We now introduce those mathematical operations that we are concerned with most of the time. Many of them will be familiar to you from the 'Einführung in Wirtschaftsmathematik', but here we learn about some very simple tricks to do the calculations in no time!

Consider a matrix $\boldsymbol{A}=\begin{pmatrix} a_{11} a_{12} a_{13}\\ a_{21} a_{22} a_{23}\end{pmatrix}$.
We can the directly get the `shape`

of the matrix:

In [21]:

```
A = np.array([[1,2,3], [4,5,6]])
A.shape
```

Out[21]:

Often we would like to *transpose* the matrix:
if $\boldsymbol{A}$ is a $2\times 3$ matrix, then its transpose $\boldsymbol{A}^T$ is a $3\times 2$ matrix.

In other words, if

$$\boldsymbol{A}=\begin{pmatrix} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23}\end{pmatrix},$$

then (abusing matrix notation a bit)

$$\boldsymbol{A}^T=\begin{pmatrix} a_{11} & a_{21} \\ a_{12} & a_{22} \\ a_{13} & a_{23}\end{pmatrix},$$

In [22]:

```
A = np.array([[1,2,3], [4,5,6]])
print(A.shape)
A_t = A.T
print(A_t.shape)
A_t
```

Out[22]:

Scalar addition/substraction functions like this:

$$\boldsymbol{A} - 5 = \begin{pmatrix} a_{11}-5 & a_{12}-5 & a_{13}-5\\ a_{21}-5 & a_{22}-5 & a_{23}-5\end{pmatrix}$$

It works for all kind of matrices and is really simple to compute:

In [23]:

```
A - 5
```

Out[23]:

If you want to add or substract matrices you need to make sure they are of the same dimension. Then for $\boldsymbol{A}=\begin{pmatrix} a_{11} a_{12}\\ a_{21} a_{22}\end{pmatrix}$ and $\boldsymbol{B}=\begin{pmatrix} b_{11} b_{12}\\ b_{21} b_{22}\end{pmatrix}$ we get the sum (or difference) by:

$$\boldsymbol{A} + \boldsymbol{B} = \begin{pmatrix} a_{11}+b_{11} & a_{12}+b_{12}\\ a_{21}+b_{21} & a_{22}+b_{22}\end{pmatrix}$$

In [24]:

```
A = np.array([[1,2], [3,4]])
B = np.array([[5,6], [7,8]])
A + B
```

Out[24]:

Suppose you want to *multiply* a matrix with some scalar $b$.
If $\boldsymbol{A}=\begin{pmatrix} a_{11} a_{12}\\ a_{21} a_{22}\end{pmatrix}$ then:

$$b\boldsymbol{A}=\begin{pmatrix} ba_{11} & ba_{12}\\ ba_{21} & ba_{22}\end{pmatrix}$$

In [25]:

```
A = np.array([[1,2], [3,4]])
b = 2
b*A
```

Out[25]:

For matrix *multiplication* to be defined, the number of columns in the first operand must equal the number of rows of the second operand.

This is because matrix multiplication is defined as follows: if If $\boldsymbol{A}=\begin{pmatrix} a_{11} a_{12}\\ a_{21} a_{22}\end{pmatrix}$ (i.e. a $2\times 2$ matrix) and $\boldsymbol{b}=\begin{pmatrix} b_{1} b_{2}\end{pmatrix}$ (i.e. a $1\times 2$ vector):

$$\begin{pmatrix} a_{11} & a_{12}\\ a_{21} & a_{22}\end{pmatrix} \begin{pmatrix} b_{1} & b_{2}\end{pmatrix} = \begin{pmatrix} a_{11} b_{1} + a_{12} b_2\\ a_{21} b_{1} + a_{22} b_2 \end{pmatrix} = \begin{pmatrix} c_{1} & c_{2}\end{pmatrix}$$

or, in other words, if $\boldsymbol{A}\boldsymbol{b}=\boldsymbol{C}$, and if $\boldsymbol{A}$ is $n\times m$ and $\boldsymbol{b}$ is $m\times p$, then $\boldsymbol{C}$ will be $n\times p$ with:

$$c_{ij}=a_{i1}b_{1j} + ... + a_{im}b_{mj} = \sum_{k=1}^m a_{ik}b_{kj}.$$

There are several ways to do this computationally:

In [26]:

```
A = np.array([[1,2], [3,4]])
B = np.array([4,5])
print(A.dot(B))
print(np.dot(A, B))
print(A@B)
```

Many ABM somehow relate to game theoretic concepts, where matrix multiplication is convenient to compute payoffs (but we will learn about this using examples later on).

The *inverse* of a square matrix $\boldsymbol{A}$ is $\boldsymbol{A}^{-1}$ such that

$$\boldsymbol{A}\boldsymbol{A}^{-1}=\boldsymbol{I}$$

where $\boldsymbol{I}$ is the identity matrix, i.e. the matrix consisting of ones on its main diagonal, and zeros outside the diagonal:

In [27]:

```
import scipy.linalg
A = np.array([[1,2], [3,4]])
print(scipy.linalg.inv(A))
print("\n")
print(scipy.linalg.inv(A)@A)
```

A standard problem in linear algebra is to solve the equation:

$$\boldsymbol{A}x = b$$

for which the key challenge is to find the inverse of $\boldsymbol{A}$.

Aside from computing $\boldsymbol{A}^{-1}$ as described above, this problem can be solved directly:

In [28]:

```
A = np.array([[1,2], [3,4]])
b = np.array([5,6])
x = scipy.linalg.solve(A,b)
```

You can verify the solution:

In [29]:

```
print(b)
print(A@x)
```

There are many more things to be done with matrices and once you understood the basics of Python, all these are straightforward to implement.
The help page of `numpy`

provides a nice overview: https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.linalg.html

For those of you interested I also upload an *optional* script that explains you how to do the following:

- Finding roots for equations
- Compute fixed points
- Solve constrained and unconstrained optimization problems

Here we skip this because it is not immediately needed for building ABMs - but it might come in handy in some more mathematical courses during your studies.

The second major application of numpy is the creation and manipulation of random numbers.
There is much functionality provided by the numpy submodule `numpy.random`

.

Here we introduce the most important concepts frequently used when using ABM.
For a complete documentation of all objects, classes and functions provided by `numpy.random`

see here.

If simply want to obtain a few random integers, where each integer appears with equal probability we can use the `randint`

function:

In [59]:

```
np.random.randint(10, 20, 5) # Creates five random integers in the range of 10 and 20
```

Out[59]:

To get random floats between zero and one, we can use the function `random`

:

In [55]:

```
size_of_sample = 10
np.random.random(size_of_sample)
```

Out[55]:

Since we can modify the resulting floats, this is a very commonly used method. It is easy to combine it with if/else statements, e.g. to simulate the toss of a fair coin we could do:

In [56]:

```
if np.random.random(1) > 0.5:
print("Head!")
else:
print("Tails!")
```

We might also choose some elements of a list randomly. If all elemens shall be chosen with the same probability:

In [63]:

```
x = np.arange(1, 5, 1)
print(x)
np.random.choice(x, 2) # choose two elements out of x
```

Out[63]:

The probabilities can also be specified using a vector of probabilities:

In [66]:

```
x_probs = [0.1, 0.25, 0.35, 0.3]
np.random.choice(x, 2, p=x_probs)
```

Out[66]:

There are several variants, such as choosing with or without replacement. Its best to have a look at the documentation of the function using `help(np.random.choice)`

.

Another common task is to shuffle the elements of a list of an array.
The function `shuffle`

does this (and thereby modifies the container):

In [44]:

```
x = np.arange(1, 10)
print(x)
np.random.shuffle(x)
print(x)
```

Out[44]:

If you do not want to change the original container, but return a new object, you can use the function `permutation`

(which, since it does not modifies the container, also works with tuples):

In [49]:

```
x = np.arange(1, 10)
print(x)
np.random.permutation(x)
```

Out[49]:

In [47]:

```
print(x) # x has not been changed
```

Often we want to draw random numbers from a particular distribution. There are many distributions pre-defined in numpy, and they are all used following the same syntax:

To draw 100 random numbers from the Normal distribution with $\mu=1$ and $\sigma=2$:

In [38]:

```
np.random.normal(1, 2, 20)
```

Out[38]:

Similarly, we can draw 10 random numbers from the Poisson distribution with $\lambda=2$:

In [39]:

```
np.random.poisson(2, 10)
```

Out[39]:

The basic principle is always the same. See the official documentation for a list of all available distributions.

Numpy uses the Mersenne Twister pseudo-random number generator to generate random numbers. The results are not actual random numbers since the process is deterministic, but for basically all purposes they simulate true random processes sufficiently closely. For example, you must generate $2^{19937}-1$ random numbers before the sequence starts repeating itself...

While the process as such is super complicated, think of it as a linear function that starts from some initial conditions and generates numbers that are hard to predict. The process is deterministic in the sense that for the same initial conditions, the same sequence of pseudo random numbers gets created.

This actually comes with a huge advantage: by setting the initial conditions for the Mersenne Twister, we can make sure that two random processes yield the same result:

In [31]:

```
np.random.sample(5)
```

Out[31]:

In [32]:

```
np.random.sample(5)
```

Out[32]:

The two results are obviously different...

But if we set the initial conditions (the 'seed') explicitly...

In [33]:

```
np.random.seed(123)
np.random.sample(5)
```

Out[33]:

In [34]:

```
np.random.seed(123)
np.random.sample(5)
```

Out[34]:

...we get exactly the same result!

This is useful when we want to verify ABM: if we run the model twice and the results differ we do not know whether the difference is due to different realizations of random processes, or due to a design change we made. By using the same seed for both runs we are able to isolate the effect of design changes.

By the way, 'real' random numbers is hard.
An attempt is made by the `quantumrandom`

project:
https://pypi.org/project/quantumrandom/

It interacts with the ANU Quantum Random Numbers Server which generates random numbers in real-time by measuring the quantum fluctuations of the vacuum.