Mathematical and Scientific Computing in Python Sean Vig Py-Camp Lightning Talk - 10/21/2013

My Background Research - Condensed matter physics Electron scattering on highly correlated electron systems Data processing and visualization

Science and Computation Computers have become necessary Data collection Data processing Theoretical modeling What do scientists need?

What do scientists need? Abstracted tasks Symbolic manipulation Exploratory discovery e.g. Mathematica, Matlab, Maple, IDL Sharing Scalable

What do scientists need? Abstracted tasks Symbolic manipulation Exploratory discovery e.g. Mathematica, Matlab, Maple, IDL Sharing Scalable Low level tasks Text processing Data storage Web interfaces SPEED

Scientific Python Stack http://scipy.org/topical-software.html

Multi-dimensional matrices, Workhorse of scientific computing I n [ 1 ] : I n [ 2 ] : Same syntax as Python, but more powerful manipulation and slicing I n [ 3 ] : I n [ 4 ] : i m p o r t n u m p y a s n p a = n p . a r a n g e ( 2 0 ) . r e s h a p e ( 4 , 5 ) a O u t [ 2 ] : a r r a y ( [ [ 0 , 1 , 2 , 3 , 4 ] , [ 5 , 6 , 7 , 8 , 9 ] , [ 1 0 , 1 1 , 1 2 , 1 3 , 1 4 ] , [ 1 5 , 1 6 , 1 7 , 1 8 , 1 9 ] ] ) a [ 2 , 3 ] O u t [ 3 ] : 13 a [ ( 1 , 3 , 0 ) , : : 2 ] O u t [ 4 ] : a r r a y ( [ [ 5 , 7 , 9 ] , [ 1 5 , 1 7 , 1 9 ] , [ 0 , 2 , 4 ] ] )

More intuitive syntax, faster I n [ 5 ] : I n [ 6 ] : I n [ 7 ] : A p y = r a n g e ( 1 0 0 0 0 0 0 ) B p y = r a n g e ( 1 0 0 0 0 0 0 ) C p y = r a n g e ( 1 0 0 0 0 0 0 ) A n p = n p . a r a n g e ( 1 0 0 0 0 0 0 ) B n p = n p . a r a n g e ( 1 0 0 0 0 0 0 ) C n p = n p . a r a n g e ( 1 0 0 0 0 0 0 ) % % t i m e i t # P u r e p y t h o n i m p l e m e n t a t i o n Z = [ ] f o r i d x i n x r a n g e ( l e n ( A p y ) ) : Z . a p p e n d ( A p y [ i d x ] + B p y [ i d x ] * C p y [ i d x ] ) 1 l o o p s , b e s t o f 3 : 1 8 4 m s p e r l o o p # N u m P y i m p l e m e n t a t i o n % t i m e i t Z = A n p + B n p * C n p 1 0 0 l o o p s , b e s t o f 3 : 3 . 3 3 m s p e r l o o p

Algorithms and tools BLAS, LAPACK, FFTPACK, QUADPACK, etc. Mathematical functions I n [ 8 ] : Integrate a bessel function j v ( 2 . 5 , x ) along the interval : I n [ 9 ] : i m p o r t s c i p y a s s p [0, 4.5] (x) dx ∫ 4.5 0 J 2.5 f r o m s c i p y i m p o r t i n t e g r a t e , s p e c i a l i n t e g r a t e . q u a d ( l a m b d a x : s p e c i a l . j v ( 2 . 5 , x ) , 0 , 4 . 5 ) O u t [ 9 ] : ( ) 1.11781793808, 7.86631718254e − 09

Drumhead example Bessel functions are a family of solutions to Bessel's differential equation: These show up when solving for vibrational modes. I n [ 1 0 ] : + x + ( − )y = 0 x 2 d 2 dx2 dy dx x 2 2 d e f d r u m h e a d _ h e i g h t ( n , k , d i s t a n c e , a n g l e , t ) : n t h _ z e r o = s p e c i a l . j n _ z e r o s ( n , k ) r e t u r n n p . c o s ( t ) * n p . c o s ( n * a n g l e ) * s p e c i a l . j n ( n , d i s t a n c e * n t h _ z e r o ) t h e t a = s p . r _ [ 0 : 2 * n p . p i : 5 0 j ] r a d i u s = s p . r _ [ 0 : 1 : 5 0 j ] x = n p . a r r a y ( [ r * n p . c o s ( t h e t a ) f o r r i n r a d i u s ] ) y = n p . a r r a y ( [ r * n p . s i n ( t h e t a ) f o r r i n r a d i u s ] ) # H e i g h t o f d r u m a s a f u n c t i o n o f x a n d y z = n p . a r r a y ( [ d r u m h e a d _ h e i g h t ( 1 , 1 , r , t h e t a , 0 . 5 ) f o r r i n r a d i u s ] )

Matlab-style plotting library Use gallery (http://matplotlib.org/gallery.html) I n [ 1 1 ] : i m p o r t p y l a b f r o m m p l _ t o o l k i t s . m p l o t 3 d i m p o r t A x e s 3 D f r o m m a t p l o t l i b i m p o r t c m

Plot drum head example from before I n [ 1 2 ] : f i g = p y l a b . f i g u r e ( ) a x = A x e s 3 D ( f i g ) a x . p l o t _ s u r f a c e ( x , y , z , r s t r i d e = 1 , c s t r i d e = 1 , c m a p = c m . j e t ) a x . s e t _ x l a b e l ( ' X ' ) a x . s e t _ y l a b e l ( ' Y ' ) a x . s e t _ z l a b e l ( ' Z ' ) p y l a b . s h o w ( )

Symbolic mathematics library Manipulate expressions symbolically I n [ 1 3 ] : I n [ 1 4 ] : I n [ 1 5 ] : I n [ 1 6 ] : f r o m s y m p y i m p o r t * x , y , z = s y m b o l s ( ' x y z ' ) i n i t _ p r i n t i n g ( ) x O u t [ 1 4 ] : x x * * 2 + 5 O u t [ 1 5 ] : + 5 x 2 c o s ( x ) O u t [ 1 6 ] : cos (x)

Calculus I n [ 1 7 ] : I n [ 1 8 ] : I n [ 1 9 ] : I n [ 2 0 ] : d i f f ( c o s ( x ) , x ) O u t [ 1 7 ] : − sin (x) e x p r = e x p ( x * y * z ) d i f f ( e x p r , x , y , y , z , z , z , z ) O u t [ 1 8 ] : ( + 14 + 52xyz + 48) x 3 y 2 x 3 y 3 z 3 x 2 y 2 z 2 e xyz i n t e g = I n t e g r a l ( x * * y * e x p ( - x ) , ( x , 0 , o o ) ) i n t e g O u t [ 1 9 ] : dx ∫ 0 x y e −x _ . d o i t ( ) O u t [ 2 0 ] : ⎧ ⎩ ⎨ ⎪ ⎪ (y + 1) dx ∫ 0 x y e −x for − y < 1 otherwise

Solvers I n [ 2 1 ] : I n [ 2 2 ] : x − y + 2 = 0 x + y − 3 = 0 s o l v e ( [ x - y + 2 , x + y - 3 ] , [ x , y ] ) O u t [ 2 1 ] : { } x : , 1 2 y : 5 2 + 2xf (x) − x = 0 df dx e −x 2 f = F u n c t i o n ( ' f ' ) d s o l v e ( f ( x ) . d i f f ( x ) + 2 * x * f ( x ) - x * e x p ( - x * * 2 ) , f ( x ) ) O u t [ 2 2 ] : f (x) = ( + ) C 1 x 2 2 e −x 2

"Interactive Python" 'Glue' of scientific Python More informative interactive environment "Magic" functions Improved introspection Tab completion Graphical notebook Rich media representations (Images, LaTeX, Javascript, Video) Inline graphics Markdown and LaTeX text formatting Parallelization tools

Working with Other Languages Hooks to call other languages Cython (C) F2PY (FORTRAN) RPy (R)

Working with Other Languages Hooks to call other languages Cython (C) F2PY (FORTRAN) RPy (R) IPython to the rescue Use 'magic' functions Functions that act on the code or Python shell https://github.com/ipython/ipython/wiki/Extensions-Index

Cython made easy I n [ 2 3 ] : I n [ 2 4 ] : I n [ 2 5 ] : % l o a d _ e x t c y t h o n m a g i c % % c y t h o n _ p y x i m p o r t f o o d e f f ( x ) : r e t u r n 4 . 0 * x f ( 1 0 ) O u t [ 2 5 ] : 40.0

Just how fast is it? Solving Poisson's equation (PDE) on a rectangular grid for fixed boundary conditions: I n [ 2 7 ] : d e f p y t h o n _ i t e r ( u ) : n x , n y = u . s h a p e f o r i i n x r a n g e ( 1 , n x - 1 ) : f o r j i n x r a n g e ( 1 , n y - 1 ) : u [ i , j ] = ( ( u [ i + 1 , j ] + u [ i - 1 , j ] ) * d y 2 + ( u [ i , j + 1 ] + u [ i , j - 1 ] ) * d x 2 ) / ( 2 * ( d x 2 + d y 2 ) ) % t i m e i t c a l c ( 1 0 0 , p y t h o n _ i t e r ) # 1 0 0 x 1 0 0 a r r a y , 1 0 0 i t e r a t i o n s 1 l o o p s , b e s t o f 3 : 3 . 8 5 s p e r l o o p

Just how fast is it? Solving Poisson's equation (PDE) on a rectangular grid for fixed boundary conditions: I n [ 2 7 ] : Python Cython, just add types I n [ 2 8 ] : d e f p y t h o n _ i t e r ( u ) : n x , n y = u . s h a p e f o r i i n x r a n g e ( 1 , n x - 1 ) : f o r j i n x r a n g e ( 1 , n y - 1 ) : u [ i , j ] = ( ( u [ i + 1 , j ] + u [ i - 1 , j ] ) * d y 2 + ( u [ i , j + 1 ] + u [ i , j - 1 ] ) * d x 2 ) / ( 2 * ( d x 2 + d y 2 ) ) % t i m e i t c a l c ( 1 0 0 , p y t h o n _ i t e r ) # 1 0 0 x 1 0 0 a r r a y , 1 0 0 i t e r a t i o n s 1 l o o p s , b e s t o f 3 : 3 . 8 5 s p e r l o o p % % c y t h o n c i m p o r t n u m p y a s n p d e f c y t h o n _ i t e r ( n p . n d a r r a y [ d o u b l e , n d i m = 2 ] u , d o u b l e d x 2 , d o u b l e d y 2 ) : c d e f u n s i g n e d i n t i , j f o r i i n x r a n g e ( 1 , u . s h a p e [ 0 ] - 1 ) : f o r j i n x r a n g e ( 1 , u . s h a p e [ 1 ] - 1 ) : u [ i , j ] = ( ( u [ i + 1 , j ] + u [ i - 1 , j ] ) * d y 2 + ( u [ i , j + 1 ] + u [ i , j - 1 ] ) * d x 2 ) / ( 2 * ( d x 2 + d y 2 ) )

Just how fast is it? Solving Poisson's equation (PDE) on a rectangular grid for fixed boundary conditions: I n [ 2 7 ] : Python Cython, just add types I n [ 2 8 ] : I n [ 2 9 ] : d e f p y t h o n _ i t e r ( u ) : n x , n y = u . s h a p e f o r i i n x r a n g e ( 1 , n x - 1 ) : f o r j i n x r a n g e ( 1 , n y - 1 ) : u [ i , j ] = ( ( u [ i + 1 , j ] + u [ i - 1 , j ] ) * d y 2 + ( u [ i , j + 1 ] + u [ i , j - 1 ] ) * d x 2 ) / ( 2 * ( d x 2 + d y 2 ) ) % t i m e i t c a l c ( 1 0 0 , p y t h o n _ i t e r ) # 1 0 0 x 1 0 0 a r r a y , 1 0 0 i t e r a t i o n s 1 l o o p s , b e s t o f 3 : 3 . 8 5 s p e r l o o p % % c y t h o n c i m p o r t n u m p y a s n p d e f c y t h o n _ i t e r ( n p . n d a r r a y [ d o u b l e , n d i m = 2 ] u , d o u b l e d x 2 , d o u b l e d y 2 ) : c d e f u n s i g n e d i n t i , j f o r i i n x r a n g e ( 1 , u . s h a p e [ 0 ] - 1 ) : f o r j i n x r a n g e ( 1 , u . s h a p e [ 1 ] - 1 ) : u [ i , j ] = ( ( u [ i + 1 , j ] + u [ i - 1 , j ] ) * d y 2 + ( u [ i , j + 1 ] + u [ i , j - 1 ] ) * d x 2 ) / ( 2 * ( d x 2 + d y 2 ) ) % t i m e i t c a l c ( 1 0 0 , c y t h o n _ i t e r , a r g s = ( d x 2 , d y 2 ) ) # 1 0 0 x 1 0 0 a r r a y , 1 0 0 i t e r a t i o n s 1 0 l o o p s , b e s t o f 3 : 2 1 . 6 m s p e r l o o p

NumPy can be just as fast... I n [ 2 8 ] : I n [ 2 9 ] : % % c y t h o n c i m p o r t n u m p y a s n p d e f c y t h o n _ i t e r ( n p . n d a r r a y [ d o u b l e , n d i m = 2 ] u , d o u b l e d x 2 , d o u b l e d y 2 ) : c d e f u n s i g n e d i n t i , j f o r i i n x r a n g e ( 1 , u . s h a p e [ 0 ] - 1 ) : f o r j i n x r a n g e ( 1 , u . s h a p e [ 1 ] - 1 ) : u [ i , j ] = ( ( u [ i + 1 , j ] + u [ i - 1 , j ] ) * d y 2 + ( u [ i , j + 1 ] + u [ i , j - 1 ] ) * d x 2 ) / ( 2 * ( d x 2 + d y 2 ) ) % t i m e i t c a l c ( 1 0 0 , c y t h o n _ i t e r , a r g s = ( d x 2 , d y 2 ) ) # 1 0 0 x 1 0 0 a r r a y , 1 0 0 i t e r a t i o n s 1 0 l o o p s , b e s t o f 3 : 2 1 . 6 m s p e r l o o p

NumPy can be just as fast... I n [ 2 8 ] : I n [ 2 9 ] : I n [ 3 0 ] : % % c y t h o n c i m p o r t n u m p y a s n p d e f c y t h o n _ i t e r ( n p . n d a r r a y [ d o u b l e , n d i m = 2 ] u , d o u b l e d x 2 , d o u b l e d y 2 ) : c d e f u n s i g n e d i n t i , j f o r i i n x r a n g e ( 1 , u . s h a p e [ 0 ] - 1 ) : f o r j i n x r a n g e ( 1 , u . s h a p e [ 1 ] - 1 ) : u [ i , j ] = ( ( u [ i + 1 , j ] + u [ i - 1 , j ] ) * d y 2 + ( u [ i , j + 1 ] + u [ i , j - 1 ] ) * d x 2 ) / ( 2 * ( d x 2 + d y 2 ) ) % t i m e i t c a l c ( 1 0 0 , c y t h o n _ i t e r , a r g s = ( d x 2 , d y 2 ) ) # 1 0 0 x 1 0 0 a r r a y , 1 0 0 i t e r a t i o n s 1 0 l o o p s , b e s t o f 3 : 2 1 . 6 m s p e r l o o p d e f n u m p y _ i t e r ( u ) : u [ 1 : - 1 , 1 : - 1 ] = ( ( u [ 2 : , 1 : - 1 ] + u [ : - 2 , 1 : - 1 ] ) * d y 2 + ( u [ 1 : - 1 , 2 : ] + u [ 1 : - 1 , : - 2 ] ) * d x 2 ) / ( 2 * ( d x 2 + d y 2 ) ) % t i m e i t c a l c ( 1 0 0 , n u m p y _ i t e r ) # 1 0 0 x 1 0 0 a r r a y , 1 0 0 i t e r a t i o n s 1 0 0 l o o p s , b e s t o f 3 : 1 2 . 4 m s p e r l o o p

Where to go from here... Check out tutorials/cookbooks NumPy (http://wiki.scipy.org/Tentative_NumPy_Tutorial) SciPy (http://wiki.scipy.org/Cookbook) matplotlib (http://matplotlib.org/gallery.html) SymPy (http://docs.sympy.org/latest/tutorial/index.html) Read the docs, modify examples, play around https://github.com/flacjacket/pycamp-lightning-talk (use http://nbviewer.ipython.org/)