Scientific Computing with Python

_intro.md

Scientific Computing with Python

Why Python?

• Programming is fun again!
• Free (as in beer and as in speech)
• Steep learning curve
• Highly readable, easy to code
• Batteries included
• Package management
• Scales pretty well (ie. namespaces, proper OO)

Difference to Matlab

• Indexing starts at 0, not 1. Look:

When to use

Use Matlab when...

• you only do throw-away scripts
• you need to build GUIs to interact (but don't need to deploy them)
• you want print-ready graphs

Use R when...

• you do data analysis, and data analysis only

Use python when...

• you want code that scales well with your project
• you collaborate with others on the code
• you want parallel computing, Sockets, interface C, C++ or Fortran,
• you want to parse a string

Setting up your Environment:

• We'll need python, ipython, matplotlib, numpy, scipy and PIL (the python image library)
• For nerds:
  • install python2.7
  • install pip (curl https://raw.github.com/pypa/pip/master/contrib/get-pip.py | python)
  • using pip, install the other stuff (ie. pip install matplotlib)
• Alternative distributions (haven't used any of them):
  • Spyeder
  • Enthought Python
  • Sage (Supports symbolic computation, not 100% python compatible though)
  • Python(x,y), Windows only
• At any rate, we can the interactive ipython console with ipython qtconsole --pylab inline
• I we feel fancy, we'll run the console in the browser: ipython notebook --pylab inline Why?
  • Everything you do is saved automatically.
  • Start it from a folder in your dropbox and you have the same environment on all your machines
  • You can even tell your code cells to be just Markdown-Formatted text-cells. This way you can annotate your code, log results, etc. A proper lab book!

Further Reading:

• Numpy for Matlab Users
• Cheat Sheet for Matlab and R users
• The Excellent SciPy Lectures
• Some interesting matplotlib recipes

_why_python?

# Scientific Computing with Python

## Why Python?

* Programming is fun again!
* Free (as in beer and as in speech)
* Steep learning curve
* Highly readable, easy to code
* Batteries included
* Package management
* Scales pretty well (ie. namespaces, proper OO)

### Difference to Matlab

* Indexing starts at `0`, not `1`. Look:

## When to use

### Use Matlab when...

* you only do throw-away scripts
* you need to build GUIs to interact (but don't need to deploy them)
* you want print-ready graphs

### Use R when...

* you do data analysis, and data analysis _only_

### Use python when...

* you want code that scales well with your project
* you collaborate with others on the code
* you want parallel computing, Sockets, interface C, C++ or Fortran,
* you want to [parse a string](http://www.youtube.com/watch?v=1lBeungEnx4&t=5m9s)

## Setting up your Environment:

* We'll need python, ipython, matplotlib, numpy, scipy and PIL (the python image library)
* For nerds:
  * install python2.7
  * install pip (`curl https://raw.github.com/pypa/pip/master/contrib/get-pip.py | python`)
  * using pip, install the other stuff (ie. `pip install matplotlib`)
* Alternative distributions (haven't used any of them):
  * [Spyeder](http://code.google.com/p/spyderlib/)
  * [Enthought Python](http://www.enthought.com/products/edudownload.php)
  * [Sage](http://sagemath.org/) (Supports symbolic computation, not 100% python compatible though)
  * [Python(x,y)](http://code.google.com/p/pythonxy/wiki/Welcome), Windows only
* At any rate, we can the interactive ipython console with `ipython qtconsole --pylab inline`
* I we feel fancy, we'll run the console in the browser: `ipython notebook --pylab inline` Why?
  * Everything you do is saved automatically.
  * Start it from a folder in your dropbox and you have the same environment on all your machines
  * You can even tell your code cells to be just Markdown-Formatted text-cells. This way you can annotate your code, log results, etc. A proper lab book!

## Further Reading:

* [Numpy for Matlab Users](http://www.scipy.org/NumPy_for_Matlab_Users)
* Cheat [Sheet for Matlab](http://mathesaurus.sourceforge.net/matlab-numpy.html) and [R users](http://mathesaurus.sourceforge.net/r-numpy.html)
* The Excellent [SciPy Lectures](http://scipy-lectures.github.com)
* Some interesting [matplotlib recipes](http://www.scipy.org/Cookbook/Matplotlib)

a_case_for_python.py

import this
love = this
this is love   # True
love is True   # False
love is False  # True
love is not True or False, love is love # (True, True)
import antigravity
import inspect
print inspect.getsource(antigravity)
# import webbrowser
# webbrowser.open("http://xkcd.com/353/")

basic_numpy.py

# If you haven't started ipython with --pylab, then do this:
# from pylab import *

# Constructing Arrays
x = array([0,2,3,4]) # Matlab: a = [0,2,3,4];

# Basic plotting
x = arange(11)
a = 2 * pi/10
y = sin(a*x)
plot(x, y)
# If you're not on --pylab inline, do this to show the plot in a new window:
# show()

# Matlab: X = [0, 1, 2, 3 ; 10, 11, 12, 13 ; 20, 21, 22, 23]
X = array([[0, 1, 2, 3], [10,11,12,13], [20,21,22,23]])
print X

# Array Slicing
X[1, 2] # Row, Column

# Caveat: Pass-by-Reference!
b = X[:,2]
print b
b[1] = 99
print X

# Work-Around: copy first.
Y = X.copy()
Y[1,2] = 42
print Y
print X

# Other ways of crating arrays:
print zeros(5)
print ones((2,4))
print identity(3)

# Gotchas:
a = array(1,2,3,4)

# Logical indexing:
x = rand(100)
x[x < .5] = 0
scatter(range(100), x)

# Polar plots
r = rand(200)
phi = rand(200) * 2 * pi
polar(phi, r, 'ro')

# Hexbins FTW!
num_data = 5000.0
x = arange(num_data) / num_data + .2
y = x * random.normal(x, 10)
scatter(x, y)
hexbin(x, y, gridsize=30)

basic_python.py

# Functions: indent blocks with 4 spaces
def hello_world(name, times=3):
    print "Say", name*times
 
hello_world("ni")

# Data Types: 
print [1, 2, "dog", True] # Lists
print {1,1,2,3,5}         # Sets
print {"Barca": 1, "Real Madrid": 2}

# List indexing and slicing
a = range(6)
print a
print a[2]        
print a[1:2]      
print a[:2]       
print a[3:]       
print a[-1]       
print a[2:-2]     
print a[-2:1:-1]  

# List comprehension:
[x**3 for x in a]

[x**3 for x in a if x % 2 == 1]

# Can also construct sets this way:
word = "supercalifragilisticexpialidocious"
{letter for letter in word} # Strings are sequences -> we can iterate over them

# Or even dicts
{letter: word.count(letter) for letter in word}

# Loops are always loops over sequences:
for item in ["eggs", "sugar", "milk"]:
    print "buy", item

for index, item in enumerate(["eggs", "sugar", "milk"]):
    print "({number:02}) {stuff}".format(number=index, stuff=item)

# Being clever:
for v in enumerate(["eggs", "sugar", "milk"]):
    print "({}) {}".format(*v)
In [68]:

# More functions:
def flex(*args, **kwargs):
    print ", ".join([str(arg) for arg in args]) # List comprehension: convert items to string
    for key in kwargs:
        print key, "=", kwargs[key]
    return 23, 42
 
x, y = flex("hello", 99 ,dog="Bark")
print x + 7

image_processing.py

# Read and display image
upf = imread("upf.png")
imshow(upf)

# Convert to grayscayle and display upright
upf_grey = mean(upf, cmap=cm.gray)
imshow(upf_grey)

# Make a cheap line filter
sobel_h = array([[1,2,1], [0,0,0], [-1,-2,-1]])
sobel_v = sobel_h.transpose()
from scipy.signal import convolve2d
edges_h = convolve2d(upf_grey, sobel_h)
edges_v = convolve2d(upf_grey, sobel_v)
edges = hypot(edges_h, edges_v)
imshow(edges)

# But actually, scipy does that:
from scipy import ndimage
edges_h = ndimage.sobel(upf_grey, axis=0)
edges_v = ndimage.sobel(upf_grey, axis=1)
edges = hypot(edges_h, edges_v)
imshow(edges)

# make a feature-mappy thing
edge_feature = ndimage.gaussian_filter(edges, 10)
imshow(edge_feature)

import Image
luminance = Image.fromarray(upf_grey).convert('L')
luminance_feature = ndimage.gaussian_filter(luminance, 10)
imshow(luminance_feature)

attention = 0.3 * luminance_feature + 0.7 * edge_feature
imshow(attention)

from scipy.ndimage import filters 
maxima = filters.maximum_filter(attention, 10) 
minima = filters.minimum_filter(attention, 10) 
diff = (maxima - minima > 25)  * 1.0 # arbitrary threshold 
maxima_smooth = ndimage.gaussian_filter(diff, 15) 
imshow(maxima_smooth)