# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# + [markdown] slideshow={"slide_type": "slide"}
# <div class="frontmatter text-center">
# <h1> Introduction to Data Science and Programming</h1>
# <h2>Class 11: Scientific programming with numpy</h2>
# <h3>IT University of Copenhagen, Fall 2019</h3>
# <h3>Instructor: <NAME></h3>
# </div>
# -

# # Source
# This notebook was adapted from:
# * Scientific Python course by <NAME>
# * a lecture of <NAME> (http://jrjohansson.github.io)
# * Python for Data Analysis by <NAME>

# + [markdown] slideshow={"slide_type": "slide"}
# ## Numpy -  multidimensional data arrays
# ### Introduction
# -

# The `numpy` package (module) is used in almost all numerical computation using Python. It is a package that provide high-performance vector, matrix and higher-dimensional data structures for Python. It is implemented in C and Fortran so when calculations are vectorized (formulated with vectors and matrices), performance is very good. 
#
# To use `numpy` you need to import the module. The standard way is:

import numpy as np

# In the `numpy` package the terminology used for vectors, matrices and higher-dimensional data sets is *array*. Numpy's array object is called `ndarray`, for N-dimensional array.
#
#

# + slideshow={"slide_type": "slide"}
# Ignore this for now, it is only for plotting
import matplotlib.pyplot as plt

# + [markdown] slideshow={"slide_type": "slide"}
# ## 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.
# * reading data from files

# + [markdown] slideshow={"slide_type": "slide"}
# ### From lists
# -

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

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

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

# + [markdown] slideshow={"slide_type": "slide"}
# The vector has 1 dimension, the matrix has 2. We learn this with `numpy.ndim`.
# -

np.ndim(v), np.ndim(M)

# + [markdown] slideshow={"slide_type": "-"}
# The `v` and `M` objects are both of the type `ndarray` that the `numpy` module provides.
# -

type(v), type(M)

# + [markdown] slideshow={"slide_type": "slide"}
# The difference between the `v` and `M` arrays is their shapes. We can get information about the shape of an array by using the `ndarray.shape` property.
# -

v.shape

M.shape

# + [markdown] slideshow={"slide_type": "slide"}
# The number of elements in the array is available through the `ndarray.size` property:
# -

v.size

M.size

# + [markdown] slideshow={"slide_type": "slide"}
# Equivalently, we could use the function `numpy.shape` and `numpy.size`
# -

np.shape(M)

np.size(M)

# + [markdown] slideshow={"slide_type": "slide"}
# So far the `numpy.ndarray` looks awfully much like a Python list (or nested list). Why not simply use Python lists for computations instead of creating a new array type? 
#
# There are several reasons:
#
# * Python lists are very general. They can contain any kind of object. They are dynamically typed. They do not support mathematical functions such as matrix and dot multiplications, etc. Implementing such functions for Python lists would not be very efficient because of the dynamic typing.
# * Numpy arrays are **statically typed** and **homogeneous**. The type of the elements is determined when the 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:
# -

v.dtype

# + [markdown] slideshow={"slide_type": "slide"}
# We get an error if we try to assign a value of the wrong type to an element in a numpy array:
# -

v[0] = "hello"

# + [markdown] slideshow={"slide_type": "slide"}
# If we want, we can explicitly define the type of the array data when we create it, using the `dtype` keyword argument: 
# -

M = np.array([[1, 2], [3, 4]], dtype=complex)
M

# Common data types that can be used with `dtype` are: `int`, `float`, `complex`, `bool`, `object`, etc.
#
# We can also explicitly define the bit size of the data types, for example: `int64`, `int16`, `float128`, `complex128`.

# + [markdown] slideshow={"slide_type": "slide"}
# ### Using 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 generate arrays of different forms. Some of the more common are:

# #### arange

# create a range
x = np.arange(0, 10, 1) # arguments: start, stop, step. Like the function range for lists!
x

# + slideshow={"slide_type": "slide"}
x = np.arange(-0.5, 0.5, 0.1) #note that here we can use floats and non-integer steps. You could not do this with lists
x

# + [markdown] slideshow={"slide_type": "slide"}
# #### linspace and logspace
# -

# using linspace, both end points ARE included
np.linspace(0, 10, 25)

# + slideshow={"slide_type": "slide"}
np.logspace(0, 1, 5, base=10)
# -

import math
np.logspace(0, 1, 5, base=math.exp(1))

# + [markdown] slideshow={"slide_type": "slide"}
# #### Random data
# -

# uniform random numbers in [0,1]
np.random.rand(5,5)

# + [markdown] slideshow={"slide_type": "slide"}
# standard normal distributed random numbers $\mu = 0$ and $\sigma^2=1$
# -

np.random.randn(5,5)

# + [markdown] slideshow={"slide_type": "slide"}
# There is a huge variety of functions that you can use: 
# <img src="files/random1.png" width="700px"/>
#
#

# + [markdown] slideshow={"slide_type": "slide"}
# and you can generate samples from all the major distributions.
# Have a look at the documentation at https://docs.scipy.org/doc/numpy-1.14.1/reference/routines.random.html:
# <img src="files/random2.png" width="700px"/>

# + [markdown] slideshow={"slide_type": "slide"}
# #### zeros and ones
# -

# You can also create arrays filled with the same element:

np.zeros((3,3))

np.ones((3,3))

# + slideshow={"slide_type": "slide"}
np.zeros((3,3))+2 # We will see more matrix and vector operations later
# -

np.ones((3,3))*4 # We will see more matrix and vector operations later

# + [markdown] slideshow={"slide_type": "slide"}
# The function `numpy.empty` initializes an "empty" array.

# + slideshow={"slide_type": "-"}
np.empty((3,3))
# -

#  <span style="color:red">It’s not safe to assume that np.empty will return an array of all zeros. In many cases it will return uninitialized garbage values!</span>

# + [markdown] slideshow={"slide_type": "slide"}
# ## Manipulating arrays
# ### Indexing
# We can index elements in an array using square brackets and indices:
# -

v = np.array([1,2,3,4])
M = np.array([[1, 2], [3, 4]])
# v is a vector, and has only one dimension, taking one index
v[0]

# + slideshow={"slide_type": "slide"}
# M is a matrix, or a 2 dimensional array, taking two indices 
M[0,1]
# -

# <img src="files/indexing.png" width="400px"/>

# + [markdown] slideshow={"slide_type": "slide"}
# If we omit an index of a multidimensional array it returns the whole row (or, in general, a N-1 dimensional array) 
# -

M

M[1]

# + [markdown] slideshow={"slide_type": "slide"}
# *** 
# ### Quick exercise
# If I have a list of lists, what is the syntax to access the first element of the first list? Pay attention to the difference with arrays!
# ***

# + [markdown] slideshow={"slide_type": "slide"}
# ### Index slicing
# -

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

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


# It works in the same way as for **lists**.

A[1:3]

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

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

# + slideshow={"slide_type": "slide"}
A[::2] # step is 2, lower and upper defaults to the beginning and end of the array
# -

A[:3] # first three elements

A[3:] # elements from index 3

# + [markdown] slideshow={"slide_type": "slide"}
# Negative indices counts from the end of the array (positive index from the beginning):
# -

A[-1:]

A[-2:]

# + [markdown] slideshow={"slide_type": "slide"}
# Index slicing works exactly the same way for multidimensional arrays:

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

A
# -

# a block from the original array
A[1:4, 1:4]

# + [markdown] slideshow={"slide_type": "slide"}
# <img src="files/slicing.png" width="400px"/>

# + [markdown] slideshow={"slide_type": "slide"}
# The devil is in the details!
# -

A[4:,:]


A[4,:]

# + [markdown] slideshow={"slide_type": "slide"}
# ### An index slice only creates a view!
# An important distinction from lists is that array slices are *views* on the original array: <span style="color:red">The data is not copied, and any modifications to the view will be reflected in the source array!</span>
# -

arr = np.arange(10)
arr

arr_slice = arr[5:8]
arr_slice

arr_slice[:] = 666
arr

arr[:] = 999
arr_slice

# <span style="color:red">If you want a copy instead, use:</span> `arr[5:8].copy()`

# + [markdown] slideshow={"slide_type": "slide"}
# ### Fancy indexing
# Fancy indexing is the name for when an array or list is used in-place of an index: 
# -

A

row_indices = [1, 2, 3]
A[row_indices,:] # this selects the second, third and fourth row of A, and all its columns

A[row_indices] #this is equivalent to the expression above

# <span style="color:red">Fancy indexing, unlike slicing, always copies the data into a new array.</span>

# + slideshow={"slide_type": "slide"}
row_indices = [1, 2, 3]
A

# + slideshow={"slide_type": "-"}
col_indices = [1, 2, -1] # remember, index -1 means the last element
A[row_indices, col_indices]
# -

# *** 
# ### Quiz
# If this is matrix `A`:  
# `[[ 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]]`
#
# Then what is the output of `A[[1,2,3], [1,-1,2]]`?
# ***

# +
# A[[1,2,3], [1,-1,2]] # Uncomment to see the solution.
# -

# #### Masks

# + [markdown] slideshow={"slide_type": "slide"}
# We can also use index masks: If the index mask is an Numpy array of data type `bool`, then an element is selected (True) or not (False) depending on the value of the index mask at the position of each element: 
# -

B = np.array([n for n in range(5)])
B

row_mask = np.array([True, False, True, False, False])
B[row_mask]

# same thing
row_mask = np.array([1,0,1,0,0], dtype=bool) #1 is true, 0 is false
B[row_mask]

# + [markdown] slideshow={"slide_type": "slide"}
# This feature is very useful to conditionally select elements from an array, using for example comparison operators:
# -

x = np.arange(-1, 5, 0.5)
x

mask = (2 < x) & (x < 4.5)
# Always use parentheses for mask conditions. Only then can you join them with & and | 
mask

x[mask]

#  <span style="color:red">The Python keywords `and` and `or` do not work with boolean arrays. Use instead `&` and `| `.</span>

# + [markdown] slideshow={"slide_type": "slide"}
# ## Functions for extracting data from arrays and creating arrays
# -

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

x = np.arange(-1, 5, 0.5)
mask = (2 < x) & (x < 4.5)
indices = np.where(mask)
print(type(indices))
indices

x[indices] # this indexing is equivalent to the fancy indexing x[mask]

# + slideshow={"slide_type": "slide"}
x
# -

# Setting values with boolean arrays works in a common-sense way. To set all of the negative values to 0 we need only do:

x[x < 0] = 0
x

# + [markdown] slideshow={"slide_type": "slide"}
# ## File I/O
# -

# ### Comma-separated values (CSV)

# A very common file format for data files is comma-separated values (CSV), or related formats such as TSV (tab-separated values). 

# !head stockholm_temperatures.dat # head is a shell command that displays the beginning of a file

# The file stockholm_temperatures.dat contains the temperature in Stockholm since 1800 until 2011. The first three columns are respectively year, month and day, and the last column is the temperature.

# + [markdown] slideshow={"slide_type": "slide"}
# To read data from such files into Numpy arrays we can use the `numpy.loadtxt` function. For example:
# -

data = np.loadtxt('files/stockholm_temperatures.dat') 

data.shape

# + slideshow={"slide_type": "slide"}
# Ignore this code for now - we will explain it later.
fig, ax = plt.subplots(figsize=(14,4))
ax.plot(data[:,0]+data[:,1]/12.0+data[:,2]/365, data[:,3])
ax.set_title('Temperatures in Stockholm')
ax.set_xlabel('Year')
ax.set_ylabel('Temperature (C)');

# + [markdown] slideshow={"slide_type": "slide"}
# Using `numpy.savetxt` we can store a Numpy array to a file in CSV format:
# -

M = np.random.rand(3,3)
M

np.savetxt("random-matrix.csv", M)

# !head random-matrix.csv

np.savetxt("random-matrix.csv", M, fmt='%.5f') # fmt specifies the format
# !head random-matrix.csv

# + [markdown] slideshow={"slide_type": "slide"}
# ### Numpy's native file format (uncompressed)
# -

# `np.save` and `np.load` are the two workhorse functions for efficiently saving and loading array data on disk. Arrays are saved by default in an uncompressed raw binary format with file extension `.npy`.

np.save("random-matrix.npy", M)
# !file random-matrix.npy # file is a shell command that displays the file type

np.load("random-matrix.npy")

# + [markdown] slideshow={"slide_type": "slide"}
# ### Numpy's native file format (compressed)
# -

# You save multiple arrays in a zip archive using `np.savez` and passing the arrays as keyword arguments:

np.savez('array_archive.npz', a=M, b=data)

# When loading an .npz file, you get back a dict-like object which loads the individual arrays:

arch = np.load('array_archive.npz')
arch['a']

# + [markdown] slideshow={"slide_type": "slide"}
# ### Special values
# -

a = np.array([-1, 0, 1, 2])
a = a/0 
a # [-1/0 0/0 1/0 2/0]

# + slideshow={"slide_type": "slide"}
np.nan == np.nan # nan is not equal to anything, not even nan

# + slideshow={"slide_type": "-"}
np.isnan(a) # nan is nan

# + slideshow={"slide_type": "-"}
np.isinf(a) # nan is not infinite

# + slideshow={"slide_type": "-"}
np.isfinite(a) # nan is not finite

# + [markdown] slideshow={"slide_type": "slide"}
# #### Selecting subsets of a real data set
# -

# We can select subsets of the data in an array using indexing, fancy indexing, and the other methods of extracting data from an array (described above), and run computations later.
#
# For example, if we want to calculate the average temperature in 1971 only, we can create a mask in the following way:

# + slideshow={"slide_type": "slide"}
# reminder, the temperature dataset is stored in the data variable:
np.shape(data)
# -

mask = (data[:,0] == 1971)
data[mask,0]

# + slideshow={"slide_type": "slide"}
data[mask,3]

# + slideshow={"slide_type": "slide"}
print("This is the mean temperature in Stockholm in 1971: "+str(np.mean(data[mask,3])))

# + [markdown] slideshow={"slide_type": "slide"}
# If we are interested in the average temperature only in a particular month, say February, then we can create a index mask and use it to select only the data for that month using:
# -

np.unique(data[:,1]) # the month column takes values from 1 to 12

mask_feb = (data[:,1] == 2)

# the temperature data is in column 3
np.mean(data[mask_feb,3])

# + [markdown] slideshow={"slide_type": "slide"}
# #### More functions: sum, prod, and trace
# -

d = np.arange(0, 10)
d

# sum up all elements
np.sum(d)

# product of all elements
np.prod(d+1)

# cummulative sum
np.cumsum(d)

# cummulative product
np.cumprod(d+1)

# + [markdown] slideshow={"slide_type": "slide"}
# When you have two dimensional objects, you can specificy along which dimension (axis) you want to perform the sum (or mean, or the maximum, etc.)
# <img src="files/sum_axis.png">
# -

x = np.array([[1, 1], [2, 2]])
print(x.sum(axis=0))   # columns (first dimension)
print(x[:, 0].sum(), x[:, 1].sum())
print(x.sum(axis=1))   # rows (second dimension)
print(x[0, :].sum(), x[1, :].sum())

# + [markdown] slideshow={"slide_type": "slide"}
# ## 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.

A

n, m = A.shape

B = A.reshape((1,n*m))
B

# + slideshow={"slide_type": "slide"}
B[0,0:5] = 5 # modify the array
B
# -

A # and the original variable is also changed. B is only a different view of the same data

# + [markdown] slideshow={"slide_type": "slide"}
# We can also use the function `flatten` to make a higher-dimensional array into a vector. But this function creates a copy of the data.
# -

B = A.flatten()
B

B[0:5] = 10
B

A # now A has not changed, because B's data is a copy of A's, not refering to the same data

# + [markdown] slideshow={"slide_type": "slide"}
# ## 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:

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

v.shape

# + slideshow={"slide_type": "slide"}
# make a column matrix of the vector v
v[:,np.newaxis]
# -

# column matrix
v[:,np.newaxis].shape

# row matrix
v[np.newaxis,:].shape

# + [markdown] slideshow={"slide_type": "slide"}
# ## 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

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

# repeat each element 3 times
np.repeat(a, 3)

# tile the matrix 3 times 
np.tile(a, 3)

# + [markdown] slideshow={"slide_type": "slide"}
# ### concatenate
# -

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

np.concatenate((a, b), axis=0)

np.concatenate((a, b.T), axis=1)

# + [markdown] slideshow={"slide_type": "slide"}
# ### hstack and vstack
# -

np.vstack((a,b))

np.hstack((a,b.T))

# + [markdown] slideshow={"slide_type": "slide"}
# ## Using arrays in conditions
# -

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

M

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")

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

# + [markdown] slideshow={"slide_type": "slide"}
# ## Using numpy on a photo
# -

# This code loads a photo.
from skimage import io
photo = io.imread('files/balloons.jpg')
type(photo)

# + slideshow={"slide_type": "-"}
photo.shape

# + slideshow={"slide_type": "slide"}
np.ndim(photo)
# -

photo

# + slideshow={"slide_type": "slide"}
# This code plots a photo. Ignore this for now - we will cover plotting in a later class.
plt.imshow(photo);

# + slideshow={"slide_type": "slide"}
plt.imshow(photo[20:120, 350:420]);

# + slideshow={"slide_type": "slide"}
plt.imshow(photo[::-1]);

# + slideshow={"slide_type": "slide"}
plt.imshow(photo[::2, ::2]);

# + slideshow={"slide_type": "slide"}
# when np.where is used with 3 arguments, it replaces all True elements with the second, and all False elements with the third
photo_masked = np.where(photo > 155, photo, 0) 
plt.imshow(photo_masked);

# + slideshow={"slide_type": "slide"}
photo_onlyblue = np.copy(photo)
photo_onlyblue[:, :, :-1] = 0
plt.imshow(photo_onlyblue);
# -

photo_bw = np.copy(photo)
mask_dark = photo.mean(axis=2) < 127
photo_bw[mask_dark, :] = 0
photo_bw[~mask_dark, :] = 255
plt.imshow(photo_bw);

# + [markdown] slideshow={"slide_type": "slide"}
# ## 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.

v1 = np.arange(0, 5)

v1 * 2

v1 + 2

# + slideshow={"slide_type": "slide"}
A * 2
# -

A + 2

# + [markdown] slideshow={"slide_type": "slide"}
# ### Element-wise array-array operations
# -

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

A * A # element-wise multiplication

v1 * v1

# + [markdown] slideshow={"slide_type": "slide"}
# If we multiply arrays with compatible shapes, we get an element-wise multiplication of each row:
# -

A.shape, v1.shape

A * v1

# + [markdown] slideshow={"slide_type": "slide"}
# ### 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: 

np.dot(A, A)

np.dot(A, v1)

np.dot(v1, v1)

# + [markdown] slideshow={"slide_type": "slide"}
# Alternatively, we can cast the array objects to the type `matrix`. This changes the behavior of the standard arithmetic operators `+, -, *` to use matrix algebra.
# -

M = np.matrix(A)
v = np.matrix(v1).T # make it a column vector

v

M * M

# + slideshow={"slide_type": "slide"}
M * v
# -

# inner product
v.T * v

# with matrix objects, standard matrix algebra applies
v + M*v

# + [markdown] slideshow={"slide_type": "slide"}
# If we try to add, subtract or multiply objects with incompatible shapes we get an error:
# -

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

M.shape, v.shape

M * v

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

StarcoderdataJupyterScriptsDedupFiltered