Upgrade to Pro — share decks privately, control downloads, hide ads and more …

Mathematical and Scientific Computing in Python

Mathematical and Scientific Computing in Python

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

Justin Olmanson

October 30, 2005
Tweet

More Decks by Justin Olmanson

Other Decks in Programming

Transcript

  1. My Background Research - Condensed matter physics Electron scattering on

    highly correlated electron systems Data processing and visualization
  2. Science and Computation Computers have become necessary Data collection Data

    processing Theoretical modeling What do scientists need?
  3. 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
  4. 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 ] ] )
  5. 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
  6. 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
  7. 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 ] )
  8. 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
  9. 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 ( )
  10. 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)
  11. 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
  12. 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
  13. "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
  14. 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
  15. 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
  16. 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
  17. 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 ) )
  18. 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
  19. 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
  20. 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
  21. 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/)