1 Documentation and Resources

[Pages:11]1 Documentation and Resources

* Download: o Requirements: Python, IPython, Numpy, Scipy, Matplotlib

o Windows: google "windows download (Python,IPython,Numpy,Scipy,Matplotlib" o Debian based: sudo apt-get install python ipython numpy scipy matplotlib

* Screencast for IPython o * Documentation for Plotting o \#-subplot o * Documentation for Arrays/Matrices/Linear Algebra o o

2 Getting Started with Numerical Computations in Python

The following listings contain about 90% of what you are going to need to do serious numerical programming in Python. The standard libraries are Numpy and Scipy. The are pretty mature but often the case for Open Source, never really finished. Especially the sparse matrix representations included in SciPy apparently has still some limitations, as can be seen in Listing 5.

1 #!/ usr/bin/env python from numpy import from numpy . l i n a l g import

Listing 1: Numpy examples

# defining matrices 6 A = array ([[1 ,2 ,3] ,[4 ,5 ,6] , [7 ,8 ,9]] , dtype = float )

print '3x3 identity matrix' print eye (3) print '3x3 matrix filled with 1' print ones ((3 ,3)) 11 print ' repeat matrix along the first axis ' print rep eat (A, 2 , a x i s =1) print 'array with 3 elements equally spaced from 0 to 5' print linspace (0 ,5 ,3) print 'populate matrix' 16 G = a r r a y ( [ [ s i n ( x)+ c o s ( y ) f o r x in l i n s p a c e ( 0 , 2 pi , 4 ) ]

for y in linspace (2 pi ,4 pi , 4 ) ] ) print G print 'functions operate elementwise on arrays' G = s i n (G)

21

# modifying matrices print 'transpose of a 2dim array' G transposed = G.T 26 print ' second column of A ' print G[ : , 1 ] print 'second column and all rows but the last ' print G[: -1 ,1] x = linspace (0 ,5 ,10) 31 print ' whole vector x ' print x print 'every second element' print x [ : : 2 ]

1

2 GETTING STARTED WITH NUMERICAL COMPUTATIONS IN PYTHON

print 'every second element in reverse ' 36 print x [ : : - 2 ]

# linalg procedures print 'rank(G)=' , 41 print ran k (G) print 'matrix multiplication if 2D arrays A \cdot B' print dot ( ones ((3 ,3)) , diag ( [ 1 , 2 , 3] ) ) print 'inner product' print dot ( [ 1 , 2 , 3 ] , [ 4 , 5 , 6 ] ) # works a l s o on l i s t s 46 print ' dyadic product v v^T ' print outer ( [1 ,2 ,3 ] ,[ 4 , 5 ] )

print 'elementwise multiplication A[i,j]* D[i,j]' print G G transposed

51

print 'inverse of G'

try :

i n v (G) # t h i s wont work since rank(G)=2 and G a 4x4 matrix

except :

56

print 'matrix inversion failed'

print 'solve A x=b' b = array ([1 ,2 ,3]) A = diag ( ones ( 3 ) ) + diag ( ones ( 2 ) , k=1) 61 x = s o l v e (A, b ) print x

+ d i a g ( o n e s ( 2 ) , k=-1)

#!/ usr/bin/env python 2 from pylab import

import numpy a s npy

Listing 2: Pylab examples

# f i g u r e with two overlayed plots x = linspace (0 ,10 ,100) 7 y = sin (x) z = cos (x) figure () p1=p l o t ( x , y , 'b - ' ) p2=p l o t ( x , z , 'r. ' ) 12 x l a b e l ( 'lala ' ) ylabel ('dumdedum') legend ( ( p1 , p2 ) , ( 'foo ' , 'bar' ) ) text (2 ,0.5 , r 'some text with Latex code $ i arrow ( 0 ,0 ,2 ,1) 17 g r i d ( 'on ' ) t i t l e ( 'foo') sav efig ( 'lala.eps')

\frac{d}{dt} \Psi =

H \Psi $')

# figure with subplots 22 f i g u r e ( )

subplot (2 ,2 ,1) p1=p l o t ( x , y , 'b - ' ) subplot (2 ,2 ,2) p2=p l o t ( x , y , 'r - ' ) 27 s u b p l o t ( 2 , 2 , 3 ) p3=p l o t ( x , y , 'g - ' ) subplot (2 ,2 ,4) p4=p l o t ( x , y , 'k - ' ) sav efig ( 'subplots.eps')

32

# f i g u r e with contour plots of z=f (x , y)= x^T A x # U orthonormal matrix u = array ([5. ,3])

2

u = u/norm ( u ) 37 U = ( ey e ( 2 ) - 2 o u t e r ( u , u ) )

# A indefinite matrix A = dot ( dot (U. T, diag ( [ 3 , - 2 ] ) ) ,U) f = lambda x : dot ( x , dot (A, x ) ) x = npy . r [ -10:10:100 j ] 42 N = sh ap e ( x ) [ 0 ] z = zeros ((N,N)) for n in range (N) :

for m in range (N) : z [m, n]= f ( [ x [ n ] , x [m] ] )

47

figure () imshow ( z , o r i g i n=' lower ' , e x t e n t =[min ( x ) , max( x ) , min ( x ) , max( x ) ] ) colorbar () 52 s a v e f i g ( ' imshow . eps ' ) ## Contour plot figure () levels = linspace ( -300,300 ,20) c o n t o u r ( z , l e v e l s = l e v e l s , o r i g i n=' lower ' , e x t e n t =[min ( x ) , max( x ) , min ( x ) , max( x ) ] ) 57 c o l o r b a r ( format=FormatStrFormatter ( r '$10 ^{% d} \ frac {\ pi }{2} $ ' ) ) sav efig ( 'contour.eps')

## combined imshow with contour figure () 62 imshow ( z , o r i g i n=' lower ' , e x t e n t =[min ( x ) , max ( x ) , min ( x ) , max ( x ) ] ) c o n t o u r ( z , o r i g i n=' lower ' , e x t e n t =[min ( x ) , max( x ) , min ( x ) , max( x ) ] ) hot ( ) #change coloring colorbar () sav efig ( 'imshow_and_contour .eps')

67

figure () ### ContourF plot c o n t o u r f ( z , o r i g i n=' lower ' , e x t e n t =[min ( x ) , max( x ) , min ( x ) , max( x ) ] ) 72 autumn ( ) colorbar () sav efig ( 'contourf.eps')

77 ### 3D PLOTS from numpy import import pylab as p import m a t p l o t l i b . axes3d as p3

82 x = y = l i n s p a c e ( 0 , 1 0 , 1 0 0 ) X,Y = meshgrid (x , y) Z = exp ( -((X-7.)2 + (Y- 5 . ) 2 ) ) + 2 exp ( -2((X-5.)2 + (Y- 5 . ) 2 ) )

f i g=p . f i g u r e ( ) 87 ax = p3 . Axes3D ( f i g )

ax . plot wireframe (X,Y, Z) ax . s e t x l a b e l ( 'X') ax . s e t y l a b e l ( 'Y') ax . s e t z l a b e l ( 'Z')

92

a = b = c = array ( linspace (0 ,3 ,20)) ax . s c a t t e r 3 D ( a , b , c , c o l o r='r ' ) sav efig ( 'plot3d.eps') #p . show ( )

97

show ( ) #pops up some windows

Listing 3: Pitfalls examples

3

2 GETTING STARTED WITH NUMERICAL COMPUTATIONS IN PYTHON

1 #!/ usr/bin/env python from numpy import

# a l l objects are passed as references # f o r example d i c t i o n a r i e s ( hash maps, c a l l e d map in c++) 6 a ={ 'foo ' : 'bar ' , 2 : 1 } print a c=a c [ 'foo' ] = 'lala ' print c 11 print a

# or arrays A = ones ((2 ,2)) B=A 16 B[ 0 , 0 ] = 1 2 . 3 2 3 print A

# a l l ? No, int , f l o a t , etc are immutable in Python

def f (x ) :

21

x = 23.

def g(x ) :

x [0] = 123.

a=2 26 b = a r r a y ( [ 1 , 2 , 3 . ] )

f (a) print a #output 2 g(b) print b #output [ 1 2 3 . ,

2. ,3.]

31

# mutable datatypes as default arguments

def f (x =[1 ,2 ,3]):

print x

36

x[0]=2423.

f ( ) #output [ 1 , 2 , 3 ]

f ( ) #output [ 2 4 2 3 . , 2 , 3 ]

41 # sol uti on : copy ! There are se ve ral ways # 1 : by f a r the best way: use copy-method of the c l a s s C = A. copy () C[0 ,1]=123456 print A

46 # 2 : use copy constructor x = [1 ,2 ,3] # this is a list y = l i s t (x ) # copies the l i s t

51 # at the moment: i nte ge r d i v i s i o n by standard ( changes in Python 3000) print 2/3 # output : 0 print 2 . / 3 #output 0.666666667

1 #!/ usr/bin/env python

Listing 4: general stuff

from numpy import

print 'output formatting ' 6 print 'fixed number of integer digits'

print ' an int with at least 5 digits length : %5d '%-23 print 'fill whitespaces with 0' print ' %05 d '%43 print 'floating point'

4

11 print ' %2.7 f '% -1.23254234 print '%2.7f'%-1232323.2 print 'exponential representation ' print '%e'%(6.231023) print ' insert a string %s '%' INPUTSTRING '

16

# elementary datatypes print 'dictionary (=hash map)' mydict = {'1' :0 , 2: 'lala ' , 'myarray ' : array ([ 1 , 2 , 3 ] ) } 21 print ' access with ' print mydict [ '1' ] , mydict [ 2 ] , mydict [ 'myarray' ]

print 'list ' mylist = [1 , 'two' ,3]

26

class myclass : a=1 b=2

myobject = myclass () 31 my ob ject . a = 2

myobject .b = 3 myobject . c = 4 # this

also

works

# looping and i f then e l s e 36 f o r i in r a n g e ( 1 0 ) : # very slow ! ! !

print i

for i , x in enumerate ( mylist ) : print 'The % d\' th element is %s! '%( i , s t r ( x ) )

41

i =0

while True : # slow ! ! !

print i

i = i +1

46

i f i ==9:

break

e l i f i == 8 :

print '8 yay!'

else :

51

print '.'

print 'normal functions '

def f (A, x ) :

56

return dot (A, x )

x = linspace (0 ,9 ,10)

A = ones ((10 ,10))

print f (A, x)

61 print ' lambda functions ( anonymous functions ) ' f = lambda x : x 2

print 'pass functions as argument of functions '

def f (g ) :

66

return g ( 1 . 5 )

print f (lambda x : x 2 )

print 'variable number of arguments , variable number of keywords ' 71 def f o o ( a r g s , keywords ) :

print "Number of arguments :" , len ( args ) print "Arguments are: " , args print "Keyword passed arguments are: " , keywords

76 f o o ( 3 , 2 , [ 2 , 3 ] , l a l a = a r r a y ( [ 1 , 2 , 3 ] ) )

5

3 OBJECT ORIENTED PROGRAMMING

#!/ usr/bin/env python

Listing 5: Sparse Matrices

from numpy import 4 from numpy . l i n a l g import

from s c i p y import from s c i p y . s p a r s e import # for sparse matrix containers from s c i p y . l i n s o l v e import # for spsolve from pylab import

9

Asp = l i l m a t r i x ( ( 1 0 , 1 0 ) )

#the case k=-1 i s buggy in SciPy 0 . 6 . 0

# either download Bleeding Edge with

# svn co http ://svn . scipy . org/svn/scipy/trunk

14 #or change sparse . py !

#def setdiag ( s e l f , values , k=0):

#""" F i l l s the diagonal elements { a i i } with the values from the

#given sequence . I f k != 0 , f i l l s the o f f -diagonal elements

#{ a { i , i+k}} instead .

19

#"""

#M, N = s e l f . shape

#i f len ( values ) > min(M, N+k ) :

#r a i s e ValueError , "sequence of target values i s too long"

#i f k>=0:

24

#f o r i , v in enumerate( values ) :

#s e l f [ i , i+k ] = v

#e l s e :

#f o r i , v in enumerate( values ) :

#s e l f [ i-k , i ] = v

29

#return

Asp . s e t d i a g ( o n e s ( 9 ) , k=-1) Asp . s e t d i a g ( ones ( 1 0 ) , k=0) Asp . s e t d i a g ( ones ( 9 ) , k=1)

34

b = ones (10) x = s p s o l v e (Asp , b ) print x

3 Object Oriented Programming

1 #!/ usr/bin/env python import numpy

Listing 6: Algorithmic differentiation

class adouble( object ):

""" derive the class from object """

6

def i n i t ( s e l f , x , dx = 0 ) :

""" this function is automatically called when >>> a = adouble (0),

i.e. runs some initialization . Usually called CONSTRUCTOR

"""

self .x = x

11

s e l f . dx = dx

def s t r ( s e l f ) :

""" string representation of adouble:

i.e. the result you see when you do: >>> print adouble(0) """

16

return '[%f ,%f] '%( s e l f . x , s e l f . dx )

def mul ( s e lf , rhs ) : """multiplication between doubles and/or adboules a*b is equivalent to a.__mul__(b)

6

21

26

31

if

36

41

""" i f numpy . i s s c a l a r ( r h s ) :

return adouble ( s e l f . x rhs , s e l f . dx rhs ) else :

return adouble ( s e l f . x rhs . x , s e l f . x rhs . dx + s e l f . dx rhs . x)

def rmul ( s e lf , lhs ) : """ 2 * adouble(0) is equivalent to 2.__mul__(adouble(0)), but integers dont have this functionality . therefore this must be a right multipilcation , i.e. adouble (0). __mul(2) """ return s e l f lh s

n a m e == " __main__ " :

### t h i s i s only executed i f t h i s s c r i p t i s c a l l e d

###by $python a l g o r i t h m i c d i f f e r e n t i a t i o n . py

a = adouble (2 ,1)

b = adouble (3)

print a b print a 2. print 2 a

1.0

foo fboaor

0.5

some

text

with

Latex

code

i

d dt

=H

dumdedum

0.0

-0.5

-1.00

2

4

6

8

10

lala

Fig. 3.1: Output produced by Listing 2.

4 Global, Local and Immutable Data Types

Python distinguishes between so-called mutable (dicts, lists, all user defined classes) and immutable data types (int, double, tuple, string). That means that once an immutable data type has been defined, it IS this value and not a reference to a memory location. This in particular is the difference between a list and a tuple:

In [ 1 ] : mytuple = (1 ,2 ,3)

In [ 2 ] : mytuple [ 0 ] = 23

3 ---------------------------------------------------------------------------

Traceback ( most recent c a l l l a s t )

/ u s r /math/ s 3 s u n u 2 / w a l t e r/ in ()

: 'tuple ' o b j e c t does not support item assignment

7

4 GLOBAL, LOCAL AND IMMUTABLE DATA TYPES

1.0

1.0

0.5

0.5

0.0

0.0

-0.5

-0.5

-1.0

-1.0

0

2

4

6

8

10

0

2

4

6

8

10

1.0

1.0

0.5

0.5

0.0

0.0

-0.5

-0.5

-1.0

-1.0

0

2

4

6

8

10

0

2

4

6

8

10

Fig. 3.2: Output produced by Listing 2.

This distinction results in subtle and unintuitive behavior as show in Listing 7. The key difference is that all mutable objects "are" a "reference" to a memory address. So while the "reference" cannot change the memory it is pointing to can.

Listing 7: Local, Global and Immutable Data Types

#!/ usr/bin/env python

a = 12 #immutable data type

b = [ 1 , 2 , 3 ] #mutable data type

4 def f ( ) :

a = 234 #l o c a l l y defined

f ()

print a #output 12

def f ( ) :

9

global a

a = 237 #now uses the global variable

f ()

print a

def f ( ) :

14

b[0]= -23 #not l o c a l l y defined , uses the global variable

f ()

print b #output [ -23 ,2 ,3]

def f ( ) :

b = { 1 : 3 } #does not change the global variable b

19 f ( )

print b #output [ -23 ,2 ,3]

def f ( ) :

global b

b = { 1 : 3 } #now overwrites global variable

24 f ( )

print b #output {1: 3}

b = [1 ,2 ,3]

def f ( ) :

29

print b

f ( ) #output [ 1 , 2 , 3 ]

b = [4 ,5 ,6]

8

 ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download