# PME Python Workshop¶

## Contents:¶

Note: This document is written for python 3.5

# Introduction¶

## What is Python?¶

Python is a general-purpose, high-level programming language. Its design philosophy emphasizes code readability. The scientific Python ecosystem is maturing fast and Python is an appealing alternative, because it is free, open-source, and becoming ever more powerful.

Why Python?

## Difference between Python and Matlab¶

### Matlab¶

• All in one package (sort of)
• Well establised in companies and universities alike

• Closed source
• Expensive
• Other institutions with which you share code should pay for a license

### Python¶

• Free
• General-purpose
• Opensource
• Large and still growning community

• Slow migration to Python 3
• Less mature libraries(toolboxes)

## Python for scientific research¶

Python uses scientific libraries for specific tasks like linear algebra, optimization, plotting, symbolic math, etc. In Matlab all these libraries (toolboxes) are always imported automatically. In case of Python you only import what you need.

### Scientific Python¶

 NumPy Base N-dimensional array package SciPy Fundamental library for scientific computing Matplotlib Comprehensive 2D Plotting IPython Enhanced Interactive Console SymPy Symbolic mathematics pandas Data structures & analysis

Great overview of the current and past scientific community

## Jupyter notebook, A.K.A. IPython notebook¶

Python is not restricted to the use of one graphical user interface e.g. Notepad, Notepad++, IDLE, Pycharm, Spyder, Jupyter notebook etc. Spyder and Jupyter are used most often in scientific research. We will be using Jupyter in this workshop.

The notebook is combination of scripting (code cells) and documentation (Markdown cells).
Jupyter Notebook Users Manual

• Markdown
• LaTeX
• Basic HTML

### Code cells¶

• Python code
• Magic functions (build-in shortcut commands)
In [ ]:
#Explanation of the magic commands
%magic

In [ ]:
#Cheatsheet
%quickref


## PME and Python¶

Hybrida - Alejandro Aragón (SOM)

# Python basics¶

## Hello World¶

In [ ]:
print("Hello", "World!") # built-in function for output


## Data types¶

#### Basic data types¶

• Booleans
• Integers
• Floats
• Compex numbers
• Stings
In [ ]:
# Run a cell by <Shift><Enter>, <Ctrl><Enter> or <Alt><Enter>.
# Find out the difference!

# Booleans
a = True
b = False

# Integers
A = -10
B = 0
C = 10982734018926734981726512087541234924092874

# Floating-point values (Floats)
α = 3. # \alpha<Tab>
β = .5 # \beta<Tab>
γ = 1e5 # \gamma<Tab>

# Complex numbers
Σ = 1+2j # \Sigma<Tab>
Ψ = 2.+0.j # \Psi<Tab>
Ω = complex(1.) # \Omega<Tab>

# Strings

descriptive_name = 'This is a string'
canBeUsefull = "This is 'also' \"a\" \n string" #start with a lowercase to avoid potential namespace issues.
# Python recommends UpperCamelCase for class names, CAPITALIZED_WITH_UNDERSCORES for constants,
# and lowercase_separated_by_underscores for other names.

In [ ]:
print(canBeUsefull)


#### Collection data types¶

• Tuples
• Lists
In [ ]:
# Tuples (immutable)
#immutable: its state cannot be modified after it is created. foo[0] = 36 will result in an error

foo = 1,2
fos = ("hello",'World',2,5j,3.14)

# Lists (mutable)
# mutable: its state be modified. qux[1] = 18 will delete [A, B, C] and place the interger 18 at its position

qux = [[a, b],
[A, B, C],
[α, β, γ],
[Σ, Ψ, Ω],
[descriptive_name, canBeUsefull]] # Autocompletion, try: descr<Tab>

%whos

In [ ]:
# tuples are great for swapping!
a, b = 4, 6
a,b = b,a
print("swapping tuple: ", a,b)
# in other programming languages this would require a more involved notation:
# tmp = a
# a = b
# b = tmp


#### Python does not work with variables but with object references instead!¶

In [ ]:
x = 4
print(id(x)) # id() returns the "identity" of an object. This is the address of the object in memory.

y = x
print(id(y))
x = 10
print(id(x))

y = 11
print(id(y))

y

In [ ]:
a = [3, 4]
b = a

print("id of a before mutation of b", id(a[0]))
print("id of b before mutation of b", id(b[0]))

b[0] = -1

print("id of a after mutation of b", id(a[0]))
print("id of b after mutation of b", id(b[0]))

print("After changing b[0] to -1, a[0] =", a[0])

print("Before reassigning a, a =", a, "and b =", b)
a = [5, 6]
print("After reassigning a, a =", a, "and b =", b)


## Logical operators¶

Logic operators add relations between data types and their instances.

• The identity operator is, is not
• The comparison operators >, >=, <, <=, ==, !=
• The membership operator in
• The logical operators and, or, not

#### The identity operator:¶

In [ ]:
string1 = "Hello  World!"
string2 = string1
string1 is string2

In [ ]:
string1 is not None # None is a built-in null object, and is often used as a place-marking value to signify “unknown” or “nonexistent”


#### The comparison operators:¶

In [ ]:
integer1 = integer2 = 15
integer1<integer2, integer1<=integer2, integer1==integer2, integer1!=integer2, integer1>=integer2, integer1>integer2
# less, equal or less, equal, not equal, equal or more, more

In [ ]:
11 < integer1 == integer2 < 20 # operators can be chained!


#### The membership operator:¶

In [ ]:
'anu' in "Peanut butter"


#### The logic operators:¶

In [ ]:
T = True
F = False
T and T, T and F, F and T, F and F

In [ ]:
T or T, T or F, F or T, F or F

In [ ]:
not True , not False , not 0, not 1, not 2, not "hello"


## Arithmetic operators¶

• Basic operators +, −, ∗, **, % (modulo operator) and /
• Compound assignment operators +=, −=, ∗=, and /=

#### Basic operators +, −, ∗, ** and / :¶

In [ ]:
1 + 1, 1 - 1, 2*2, 10/2, 'First' + 'Last', 3**2, 10%4
#modulo opperator % : divides left hand operand by right hand operand and returns remainder


#### Compound assignment operators +=, −=, ∗=, **=, %= and /= :¶

In [ ]:
n = 0
n += 2
n *= 2
n -= 1
n **= 2
n %= 5
n /= 2

n

In [ ]:
colors = [" red ", " yellow "]
colors += [" green "]
print(colors)
colors += "blue"
colors

In [ ]:
n = 4
print(id(n))  # The built-in id() function returns the object’s address in memory
n += 3
print(id(n))
print(n)

In [ ]:
n = None
del n
%whos


## Functions¶

Create functions to perform set of operations repeatedly

def functionName(arguments) :
Code block

In [ ]:
def mid(num1,num2):
average = 0.5*(num1+num2)
return average

mid(4,19)

In [ ]:
def getstring(var1, var2, var3):
var1 = 'Hello'
var2 = 'World'
var3 = '!'
return var1,var2,var3 #defines a tuple on the fly

k, l, m = 0, 0, 0
k, l, m = getstring(k,l,m)
print(k, l, m)

In [ ]:
#internal variables do not show up as saved variables
%whos


#### Included libraries¶

Python has some basic functions included like:

• len(a) : length of a list or string
• int(), float(), str(), boolean() : type conversion
• range() : Produces a sequence of integers from start to stop by a certain step
• input() : asks the user for an input and outputs it as a string

For specialized tasks, libraries of pre-defined functions are available. These libraries are called modules in Python and can be imported using the following syntax:

from module import function as functionNameYouLike
from module import * # not recommended! Loads the module into the local namespace, no prefix will be nessecary to use  the functions of the module


or

import module as moduleNameYouLike
import module.function as functionNameYouLike


## Control flow statements¶

Control flow statements determine execution order/paths in a program

• The if ... else statement
• The while loop statement
• The for ... in loop statement
• The try ... except... else statement

#### The if statement:¶

if boolean_expression1:
# Execute this block (called a suite in Python)
suite1
elif boolean_expression2:
# Execute this block
suite2
...
elif boolean_expressionN:
suiteN
else:
# Execute this block when the other expressions do not hold
else_suite

In [ ]:
#input() asks the user for an input and outputs it as a string
#int() converts the input to an interger
#print() displays the input

temp = int(input('What is the water temperature in degrees Celcius? '))
if temp < 0:
print("It's frozen!") #Put your cursur between the brackets of the function print() and try <Shift><Tab>. Also try pressing it multiple times.
elif temp > 100:
print("It's a gas!")
else :
print('At',temp ,'degrees Celcius water is liquid. Go with the flow...' )


#### The while loop statement:¶

while boolean_expression:
suite


Two important built-in statements for loops:

break # breaks out of the (innermost) loop
continue # switches control to the start of the loop

In [ ]:
temperature = 100

while temperature >= 80: # first while loop code
print(temperature)
temperature = temperature - 1


In [ ]:
while True:
n = input("Please enter 'hello': ")
if n.strip() == 'hello':
break

In [ ]:
var = 10
while var > 0:
var -= 1
if var == 5:
continue
print('Current variable value :', var)
print("Good bye!")


#### The for ... in loop statement:¶

for variable in iterable:
Execute this block


As before,

break breaks out of the (innermost) loop
continue switches control to the start of the loop

In [ ]:
for letter in 'lalelilolu':
if letter in 'aeiou':
print(letter, 'is a vowel')


#### The try ... except... else statement:¶

"It’s easier to ask forgiveness than permission” is the Pythonic way to go! Other programming languages “look before you leap”, where checks are made in advance.

try:
# you do your operations here
try_suite
except Exception1 as Variable1:
# if there is Exception1, then execute this block
exception_suite1
...
except ExceptionN as VariableN:
# if there is ExceptionN, then execute this block
exception_suiteN
else:
# executed when the try block’s suite has finished normally
else_suite
finally:
# always executed at the end
finally_suite

In [ ]:
age = input('Enter your age : ')
try:
i = int(age)
except ValueError as err:
print(err)
else:
print('Your age is : ', i)


# Exercise¶

If we list all the natural numbers below 10 that are multiples of 3 or 5, we get 3, 5, 6 and 9. The sum of these multiples is 23.

Find the sum of all the multiples of 3 or 5 below 1000.

Tip: use the modulo % AND use the help function, OR the question mark range?, to find out more about range()


Done already? Check your answer here. Do you have time left? Don't worry, there is enough practice material over at project Euler.

# Frequently used scientific libraries¶

## NumPy¶

NumPy is the fundamental package for scientific (numerical) computing with Python. It includes discrete Fourier transforms, more complex linear algebra operations, random number capabilities, size/shape/type testing of arrays, splitting and joining arrays, histograms, creating arrays of numbers spaced in various ways, creating and evaluating functions on grid arrays, treating arrays with special (NaN, Inf) values, set operations, creating various kinds of special matrices, and evaluating special mathematical functions (e.g., Bessel functions).

Numpy is written in C to increase the speed of calculations.

Frequently used Matlab functions and their numpy equivalents

In [ ]:
import numpy as np
a_np = np.array([3, 4, 8]) #Put cursor between brackets of np.array(), hold <Shift> and press <Tab> twice or three times.
b_np = np.array([7, 2, 39], dtype=float)
c_np = np.linspace(1,10,4) # start, end, num-points
d_np = np.arange(1, 9, 2) # start, end (exclusive), step

print(a_np, b_np, c_np, d_np)

In [ ]:
A_np = np.array([[ 1,  2,  3,  4],
[ 5,  6,  7,  8],
[ 9, 10, 11, 12],
[13, 14, 15, 16]])
B_np = np.random.rand(4,4)
print("A_np =\n", A_np,"\n B_np =\n", B_np)


#### Dimensions¶

In [ ]:
print('Number of dimensions =',np.ndim(A_np),
', Number of elements =', np.size(A_np),
', Shape =', np.shape(A_np))


#### Indexing and slicing¶

The items of an array can be accessed and assigned to the same way as other Python sequences (e.g. lists):

In [ ]:
a_np[0], a_np[2], a_np[-1], d_np[0:3:2] # [start:end:step]


Indices begin at 0, like other Python sequences (and C/C++). In contrast, in Fortran or Matlab, indices begin at 1. The usual python idiom for reversing a sequence is supported:

In [ ]:
print(a_np)
print(a_np[::-1])


For multidimensional arrays, indexes are tuples of integers:

In [ ]:
A_np[1, 1] # second line, second column

In [ ]:
A_np[2, 1] = 99 # third line, second column
A_np

In [ ]:
A_np[1:2] = c_np #second line
A_np

In [ ]:
A_np[:,2] = d_np # third column
A_np


#### Elementwise multiplication¶

In [ ]:
print(A_np*B_np,
'\n\n',
a_np*b_np)


#### Multiplication (or dot product depending on the shape)¶

In [ ]:
print("A_npB_np =\n", A_np.dot(B_np),
'\n\na_np*b_np =\n', a_np.dot(b_np))

In [ ]:
print("A_npB_np =\n", A_np@B_np,
'\n\na_np*b_np =\n', a_np@b_np) #Optional short hand notation since Python 3.5


#### Cross product¶

In [ ]:
np.cross(a_np,b_np)


## SciPy¶

SciPy greatly extends the functionality of the NumPy routines.A large community of developers continually builds new functionality into SciPy. A good rule of thumb is: if you are thinking about implementing a numerical routine into your code, check the SciPy documentation website first. Chances are, if it's a common task, someone will have added it to SciPy. We will not cover what SciPy has to offer in detail, but in the table below we mention a subset of its capabilities:

import scipy as sc

Module Code for
sc.constants Many mathematical and physical constants.
sc.special Special functions for mathematical physics, such as iry, elliptic, bessel, gamma, beta, hypergeometric, parabolic cylinder, mathieu, spheroidal wave, struve, and kelvin functions.
sc.integrate Functions for performing numerical integration using trapezoidal, Simpson's, Romberg, and other methods. Also provides methods for integration of ordinary differential equations.
sc.optimize Standard minimization/maximization routines that operate on generic user-defined objective functions. Algorithms include: Nelder-Mead Simplex, Powell's, conjugate gradient, BFGS, least-squares, constrained optimizers, simulated annealing, brute force, Brent's method, Newton's method, bisection method, Broyden, Anderson, and line search.
sc.linalg Much broader base of linear algebra routines than NumPy. Offers more control for using special, faster routines for specific cases (e.g. tridiagonal matrices). Methods include: inverse, determinant, solving a linear system of equations, computing norms and pseudo/generalized inverses, eigenvalue/eigenvector decomposition, singular value decomposition, LU decomposition, Cholesky decomposition, QR decomposition, Schur decomposition, and various other mathematical operations on matrices.
sc.sparse Routines for working with large, sparse matrices.
sc.interpolate Routines and classes for interpolation objects that can be used with discrete numeric data. Linear and spline interpolation available for one -and two-dimensional data sets.
sc.fftpack Fast Fourier transform routines and processing.
sc.signal Signal processing routines, such as convolution, correlation, finite fourier transforms, B-spline smoothing, filtering, etc.
sc.stats Huge library of various statistical distributions and statistical functions for operating on sets of data.

## Matplotlib¶

Matplotlib is a plotting library for the Python programming language and its numerical mathematics extension NumPy.

In [ ]:
%pylab inline


#### Cosine¶

In [ ]:
X = np.linspace(0, 6*pi, 100)
Y = np.cos(X)
plt.plot(X, Y, 'o')       # line plot


#### Random image¶

In [ ]:
image = np.random.rand(30, 30)
plt.imshow(image, cmap=plt.cm.hot)
plt.colorbar()


## SymPy¶

SymPy is a symbolic computation library. SymPy is written entirely in Python and does not depend on any additional libraries.

#### Example: calculate the momentum of a rotating rigid body¶

Consider the Euler equations for the rotational dynamics of a rigid body about its center of mass given by:

$$\dot{\bf p}= J \frac{\partial H}{\partial {\bf p}} = \left[\begin{matrix}0 & - p_{3} & p_{2}\\p_{3} & 0 & - p_{1}\\- p_{2} & p_{1} & 0\end{matrix}\right] \left[\begin{matrix} \frac{\partial H}{\partial p_1}\\\frac{\partial H}{\partial p_2}\\\frac{\partial H}{\partial p_3}\end{matrix}\right]$$

and Hamiltonian, a function comparable to the Lagrangian, the total kinetic energy:

$$H = \frac{1}{2} \left(\frac{p_{1}^{2}}{I_{1}} + \frac{p_{2}^{2}}{I_{2}} + \frac{p_{3}^{2}}{I_{3}} \right)$$

where $p_k$, for $k = 1 , 2 , 3$, are the angular momenta relative to the axis of a frame fixed to the body and $I_1 > I_2 > I_3$ are the principle moments of inertia.

#### Show the system is lossless, meaning $\dot{ H} = 0$¶

In [ ]:
import sympy as sp
sp.init_printing()

In [ ]:
t, g_2, g_3, I_1, I_2, I_3 = sp.symbols('t g_2 g_3 I_1 I_2 I_3') # Create variables

p_1 = sp.Function("p_1")(t) #create functions
p_2 = sp.Function("p_2")(t) #time dependency is only nessacery because of differentiation later on
p_3 = sp.Function("p_3")(t)

p = sp.Matrix([
[p_1],
[p_2],
[p_3]
]) #set up momentum vector

H = sp.Rational(1, 2)*( p_1**2/I_1 + p_2**2/I_2 + p_3**2/I_3) #set up the hamiltonian. Similar to a Lagrangian

J = sp.Matrix([
[0, -p_3, p_2],
[p_3, 0, -p_1],
[-p_2, p_1, 0]
]) #set up the J matrix

delH = sp.Matrix([
[sp.diff(H, p_1)],
[sp.diff(H, p_2)],
[sp.diff(H, p_3)]
]) # set up delH/delp

pd = J*delH ; pd #calculate the derivative of p

In [ ]:
Hd = sp.diff(H,t); Hd #calculate the derivative of the Hamiltonian

In [ ]:
Sd = sp.simplify(sp.Matrix([H]).jacobian(p) * pd)[0] #the [0] is to convert the matrix datatype to a single dimension datatype.
print("Q: Is the system lossless? A: This is",Sd==0)


## pandas¶

pandas is an library providing high-performance, easy-to-use data structures and data analysis tools.

Introduction to pandas

Extensive pandas demo by the Scipy comunity

Quick pandas guide