Mathematical and Scientific Computing in Python

a talk by Sean Vig at the Py-Camp Lightning Talk presented 10/21/2013

Transcript

1. Mathematical and Scientific Computing in Python Sean Vig Py-Camp Lightning Talk - 10/21/2013
2. My Background Research - Condensed matter physics Electron scattering on highly correlated electron systems Data processing and visualization
3. Science and Computation Computers have become necessary Data collection Data processing Theoretical modeling What do scientists need?
4. What do scientists need? Abstracted tasks Symbolic manipulation Exploratory discovery e.g. Mathematica, Matlab, Maple, IDL Sharing Scalable
5. 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
6. Scientific Python Stack http://scipy.org/topical-software.html
7. 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 ] ] )
8. 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
9. 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
10. 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 ] )
11. 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
12. 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 ( )
13. 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)
14. 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
15. 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
16. "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
17. Working with Other Languages Hooks to call other languages Cython (C) F2PY (FORTRAN) RPy (R)
18. 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
19. 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
20. 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
21. 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 ) )
22. 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
23. 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
24. 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
25. 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/)