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

Advantages

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

Disadvantages

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

Python

Advantages

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

Disadvantages

  • No Simulink
  • 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 cell capabilities

  • Markdown
  • LaTeX
  • Basic HTML

Code cells

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

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
#addition, subtraction multiplication, devision, addition, exponentiation, modulo
#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

print('The tea is ready!')
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