Scientific computing * Python 2.x

Module/Pakete math, pylab, numpy, scipy, cmath, python-tk, ipython/matplotlib installieren und obectivec probieren.

Run python with alpha numeric output

Welche python-Version und Module nutze ich? /usr/lib64/python2.5 Beachte 3/2.-3/2-1/2-0.5=0.0 !

# -*- coding: utf-8 -*-


print ('%.40f'%(1./3))
print ('{0:.40f}'.format(1./3))
python -c "x='python -c ';print(2**333)"
from sys import *
setrecursionlimit(2000)
Selbstreproduktion:
python -c "x='python -c %sx=%s; print x%%(chr(34),repr(x),chr(34))%s'; print x%(chr(34),repr(x),chr(34))"

Arbeiten in der Konsole:

if i==1: print(i);print(i*i)

for i in [1,3,9]:
	print(i)

for i in 2,4,7,8,99:
	print(i)

for i in [-2,4,77,'Nein ',(3,8),pi,2j,i,e,i*3j,i]: print('Aus ',i,' wird ',2*i)

for i in range(1,111):
	if i==9:print(i);print(i*i)

for i in range(1,7):
...  print(i)	# Einruecken!

g='15.9', G=float(g), 2*g,g+g, 2*G-G-G, int(G)

>>> import math
>>> help(math)
>>> 99-70*math.sqrt(2)-1/(99+70*math.sqrt(2))
>>> math.sqrt(-2), math.log(math.e)
>>> cmath.sqrt(-2)
>>> cmath.sqrt(1+1j)-cmath.exp(1j*cmath.pi/8)*cmath.sqrt(cmath.sqrt(2))

pylab - objektorientiertes Plot-modul

from pylab import *
help(pylab) matlab/scilab/octave-analog

• Plot: acorr, cohere, contour, csd, errorbar, hist, hold, pcolor, plotfile, psd, scatter, show, specgram, xcorr
• Matrix: range(min,max-,step), arange(min,max-,step), linspace, mgrid, ogrid, cumprod, cumsum, detrend, diag, diff, eig, eye, linspace, ones, rand, svd, zeros
• Wahrscheinlichkeit/Statistik: random, randn(N), randint(min,max-,N), levypdf, normpdf, corrcoef, mean, median, norm, std
• TSA: bartlett, blackman, cohere, csd, fft, hamming, hanning, hist, kaiser, psd, sinc
• Weiteres: append, date2num, angle, griddata, load, polyfit, polyval, roots, save, trapz

• copy: append
• Ma
• W

Definition der Funktion

def my_count(N):
	""" my_count(N) computes the partial sum of N integers """
	i=0
	s=0
	while(i 1:
        return x * fakultaet(x - 1)
        print(x)
    else:
        return 1

x=12; y=fakultaet(x); print("Fakultaet von",x,'ist', y)

# logmap
def map(p,x,n):
   i = 0
   orbit=[]
#   y = x
   while(i 0 """
# generiert gleichverteilte Pseudozufallszahlen auf dem Intervall [0,1)   
    import time
#    from time import *
    a=65539;b=0;m=2**31
#   print(a,b,m)
    x=time.time()
#    x=time()
# aktuelle Zeit time.asctime(time.localtime(time.time()))
    r=[]
    for i in range(1,N+1):
        x=(a*x+b)%m
        z=x/m
        r.append(z)
#       print(z)
    return r


#!/bin/python
# Iteration einer Map
# /usr/cnld/shw/TeX/introductions/Praktikum+Exercise/2010WS_ScientificComputing/5-example-codes/2_2_1DMap.py

""" Hier die Idee des Programms beschreiben
      Referenzen angeben
      Autor und Datum
      Wo ist die zu finden
                                              """

# Module laden

from pylab import *

# Parameter festlegen

TotaleSchrittzahl       = 10
Parameter               =  3.678
Startwert               = rand()

# Eigene Funktionen kreieren

def  map(zustand,parameter):
    """ logistische Abbildung 
        zustand muss im Intervall [0,1] liegen
        parameter muss aus [0,4] sein 
    """
#    zustand=zustand*parameter*(1.-zustand)
    zustand*=parameter*(1.-zustand)
    return zustand

def ErzeugeTrajektorie(Startwert, Parameter, TotaleSchrittzahl):
    """ Inizialisierung eines Arrays, in das die Trajektorie geschrieben wird """
    Trajektorie                 = zeros(TotaleSchrittzahl, dtype='float')
    x=Startwert
    for schritt in range(TotaleSchrittzahl):
        x=map(x,Parameter)
        Trajektorie[schritt]=x
    return Trajektorie

def ErzeugeTrajektoriePerListe(Startwert, Parameter, TotaleSchrittzahl):
    """ Inizialisierung einer Liste, in das die Trajektorie geschrieben wird """
    Trajektorie=[]
    x=Startwert
    for schritt in range(TotaleSchrittzahl):
        x=map(x,Parameter)
        Trajektorie.append(x)
    return Trajektorie

# Hauptprogramm

#Trajektorie = ErzeugeTrajektorie(Startwert, Parameter, TotaleSchrittzahl)
Trajektorie = ErzeugeTrajektoriePerListe(Startwert, Parameter, TotaleSchrittzahl)
plot(Trajektorie)
show()

Anonyme Funktionen: lambda

f=lambda x,y,z:x+y+z;l=sin(pi*f(1,1,1)/3.)
x=(lambda a="fee", b="fie",c="foe":a+b+c);x("wwe"); x("adc","hh"); x("X","","");
L=[(lambda x:x**2), (lambda x:x**3), (lambda x:x**4)]
for f in L: print(f(2))
print(L[0](3))
Eigenwert-Einzeiler: K=[1,2,3]; w1=1;w2=2
plot( K, imag(array( map( lambda k: eig(matrix([ [0,1,0,0],[-w1**2,0,k,0],[0,0,0,1],[k,0,-w2**2,0] ]))[0] , K )) ) ); show()

Zufallsvariable

from pylab import *
help(hist)
yf=y[y < 10]
subplot(3,2,x)
hist(randint(1,7,N),bins=6,cumulative=True)
x=linspace(0,50,1000)+2*pi*rand(1)
N=len(x)
T=13
s=sin(2*pi*x/T)
r=randn(N)
y=s+r
pylab.plot(x,s,'r.')
pylab.plot(x,y)
pylab.show()

Run python scripts

Shell command 		> python script01.py

Python command 		>>> execfile('LogisticMap.py') raus mit Ctrl-d oder quit()

gnuplot command 	gnuplot> plot "< python  LogisticMap.py" using 2 w lp raus mit quit

__________________________________________________________
# file: script01.py	# leitet Kommentazeile ein
print('Hallo Welt')	# print gibt auf Bildschirm aus
print(2 ** 99)
__________________________________________________________
__________________________________________________________
# file: Orbit.py
for i in range(7):
    print "  "*i,"Dieses Programm wurde soeben ins Leben gerufen"
__________________________________________________________

oder
__________________________________________________________
# file: Fibonacci-Reihe.py
a, b = 0, 1
print "Fibonacci-Reihe"
print "Start: ",a,b
while b < 10200003040:
        print b,
        a, b = b, a+b
__________________________________________________________
oder
__________________________________________________________
# file: Fakultaet.py
# Definition der Funktion x!

def fakultaet(x):
    if x > 1:
        return x * fakultaet(x - 1)
        print x
    else:
        return 1

        x=12
        print x
        y=fakultaet(x)
        print y
__________________________________________________________

oder
__________________________________________________________
# file: Count1.py 
# http://www-alt.physik.uni-muenchen.de/kurs/Computing/python_05/node22.html
sum = 0.
sum_tot = 1.0
for i in range(10):
  sum += 0.1
  print('i=',i,"    Summe=",sum);
if ( sum != sum_tot ):
  print " Python kann nicht rechnen: 1 - 1 = %g" %  (sum - sum_tot)
__________________________________________________________

import numpy
from numpy import mean

Plotten

from pylab import randn
r=randn(500)
pylab.plot(r)
pylab.show()
__________________________________________________________
# file: LogisticMap.py 
# gnuplot > plot " < python  LogisticMap.py " using 2 w l
# Python:>>> execfile('LogisticMap.py')
# http://board.gulli.com/thread/1269828-python-python-script-als-exe/
# 
# mathplotlib #
import pylab
i, N, r, x = 0, 100, 4, 0.123456789
print "# Logistic map"
print "# ",N," points for control parameter ",r
orbit=[]
y = x
while i < N:
        i, x = i+1, x * r * (1-x)
        print y,x
        y = x
        orbit.append(x)
# http://www.daniweb.com/code/snippet691.html
# pylab.xlabel("\alpha")
# pylab.ylabel("x_i")
# pylab.plot(orbit,'bd--')
# pylab.show()
__________________________________________________________
http://matplotlib.sourceforge.net/

from pylab import randn, hist
import matplotlib.pyplot as plt
x=randn(1000)
pylab.plot(x)
pylab.show()
plt.hist(x, 100)
pylab.show()
__________________________________________________________
Image plotten:
from pylab import imshow,randn,show  # function pylab.imshow
rows=100
columns=300
Gaussian_random_field=randn(100,300)
imshow(Gaussian_random_field)
show()

Subplots

subplot(nrows,ncols,plotn)
Rows=2;Cols=3
subplot(Rows,Cols,1);plot(X);subplot(Rows,Cols,2);plot(rabbits,foxes);subplot(Rows,Cols,3);
hist(rabbits,20);subplot(Rows,Cols,4);psd(rabbits);subplot(Rows,Cols,5);specgram(foxes);
subplot(Rows,Cols,6);xcorr(rabbits/1000.,foxes);show()

Kraut und Rüben-Daten lesen, schreiben und plotten - IO

plotfile(fname, cols=(0,), plotfuncs=None, comments='#', skiprows=0, checkrows=5, delimiter=',', **kwargs)
lo5.dat
8.331 13.291 18.063 150.483
8.218 13.630 17.879 146.93
8.767 14.960 19.103 167.476
9.331 15.232 20.690 193.058
10.076 15.583 21.442 216.05
11.268 15.997 21.670 244.178
......

from pylab import *
plotfile('lo5.dat',cols=(0,1,2,3),delimiter=' '); show()

from pylab import load
a=load('Bier.dat')
x=a[:,0]
y=a[:,1]
pylab.plot(x,y)
pylab.show()

# execfile('Histogram_from_data_file.py')
from pylab import *

datenList = [] # leere Liste
f = open("daten.txt")
for line in f:
    line = line.rstrip()
    datenList.append(float(line))

# datenList[0:4]

hist(datenList,10)
show()

# Plot_2_columns_from_data_file.py
from pylab import *

xList = [] # leere Liste
yList = []

f = open("Bier_without_header.dat")

for line in f:
    line = line.rstrip()
    parts = line.split()
    xList.append(float(parts[0]))
    yList.append(float(parts[1]))

plot(xList,yList, '-ro')        # help(plot)
# axis([0,10,0,1])
show()

# Plot_of_2_2_colum_data_files.py
from pylab import *

def readXYDataFromFile(filename):
    xList = []
    yList = []

    f = open(filename)

    for line in f:
        line = line.rstrip()
        parts = line.split()
        xList.append(float(parts[0]))
        yList.append(float(parts[1]))

    return (xList,yList)

(xListA, yListA) = readXYDataFromFile("linie1.txt")
(xListB, yListB) = readXYDataFromFile("linie2.txt")

plot(xListA, yListA,'-r' , label="red")
plot(xListB, yListB, 'b-.', linewidth=3., label="blue")
legend()

show()

__________________________________________________________
# file: IO.py
# IO operation
# Standard input
# num = int( raw_input("Zahl eingeben:") )
# Python:>>> execfile('IO1.py')
# Redirection to statndard output python IO1.py > OrbitOutput.dat
import pylab

num = 19
f = file('numbers.txt','w')
liste = ["%d\n" % x for x in range(num)]
f.writelines(liste)
f.close()

from pylab import *
num = 19
f = file('numbers.txt','w')
s = ['alt','neu','wo?','Ich','gehe','jetzt']
r = linspace(-1.1,pi,len(s))
liste = ["Nummer %d\n" % x for x in range(num)]
f.writelines(liste)
sliste= ["Wort %s %0.3f\n" % (s[i], r[i]) for i in range(len(s))]
f.writelines(sliste)
f.close()


i, N, r, x = 0, 100, 4, 0.123456789
print "# Logistic map"
print "# ",N," points for control parameter ",r
orbit=[]
y = x
#f = file('Orbit.dat','w')
while i < N:
        i, x = i+1, x * r * (1-x)
        print y,x
#       liste=[y,x]
#       f.writelines(liste)
        y = x
        orbit.append(x)

#liste=["%f\n", orbit]
# f.write(orbit)
# f.writelines(orbit)
#f.close()
# http://www.scipy.org/Cookbook/InputOutput
from numpy import *
from scipy.io import write_array
#from scipy.io import read_array
write_array("OrbitLogMap0.dat", orbit)

Input and Output: open, format, pickle

http://docs.python.org/tutorial/inputoutput.html

>>> f = open('/tmp/workfile', 'w')
>>> print f

r, w, r+, a

print 'The value of PI is approximately %5.3f.' % math.pi

pickle.dump(x, f)
x = pickle.load(f)


# Matplotlib (pylab)
from numpy import *
from pylab import load           # warning, the load() function of numpy will be shadowed
from pylab import save
save('OrbitLogMap1.dat', orbit)
#read_data = load("myfile.txt")

# numPy
savetxt('OrbitLogMap2.dat', orbit, fmt="%12.6G")    # save to file
# from numpy import *
# data = loadtxt('table.dat', unpack=True)

# Hash im Datenfile bedeutet standardmaessig Kommentar
 x,b=load("bier.d",skiprows=3,usecols=[0,1],unpack=True)
zeit, data = loadtxt("Sonne.dat",usecols=(0,1), unpack=True)


# http://www.daniweb.com/code/snippet691.html
# pylab.xlabel("\alpha")
# pylab.ylabel("x_i")
# pylab.plot(orbit,'bd--')
# pylab.show()
# 
# http://www.wspiegel.de/pykurs/python09.htm
f = file('PythonProgram.py','w')
programm = 'for i in range(7):\n'
programm = programm + '    print "  "*i,"Dieses Programm wurde soeben ins Leben gerufen"\n'
f.write(programm)
f.close()

__________________________________________________________
# Daten von URL lesen

from pylab import *
import urllib, urllib2

# URL der Datenquelle

URLdata="http://www.agnld.uni-potsdam.de/~shw/Lehre/lehrangebot/2004WS-NLD/Exercises/x.dat"
# alte Wetterdaten 

# Zieldateiname

fileName="YyY.dat"

# URL abfragen, oeffnen und runterladen

req = urllib2.Request(URLdata)
urllib2.urlopen(req)
urllib.urlretrieve(URLdata,fileName)

# File laden und plotten

Data=load(fileName)
plot(Data)
show()

Listen und Arrays

D=Array[Liste] oder L=list(D)

from pylab import arange

Shuffle: r=rand(N); r[argsort(rand(N))]

x=range(-3,12,4)
y=arange(-3,12,4)

a = np.array([[1,2,3], [4,5,6]])
np.reshape(a, 6)
    array([1, 2, 3, 4, 5, 6])
np.reshape(a, 6, order='F')
    array([1, 4, 2, 5, 3, 6])
np.reshape(a, (3,-1))       # the unspecified value is inferred to be 2
    array([[1, 2],
           [3, 4],
           [5, 6]])
 reshape(a,(3,2))
array([[1, 2],
       [3, 4],
       [5, 6]])

Listen:

a=range(n); a=a[-m:]+a[:-m]
Listen: x=[-1 3 7]; x=range(10); y=range(1,10); z =range(-2,10,2); r=2*x; a=[2,1, 4]; len(a); a.sort; a.append(-11); a.insert(2,77); x=[0]*3

Listenmethoden: s.append(x),s.extend(t),s.count,s.insert(i,x),s.pop([i]),s.remove(x),s.reverse(),s.sort([key[,revers]])

Tensorprodukt fuer Listen l1 = ['a', 'b', 'c'];l2 = [1, 2, 3];crossprod = [(x,y) for x in l1 for y in l2]

Arrays:

A=array(a) macht aus Liste ein Array.

a.shape;
Arrays, Matrizen: x=3*[range(4)]; x[2][3]; y=4*[[1]*2]; y=4*[[1]*2]; y[3][0]
a=[[1,2,3],[2,1,4],[0,1,0]];det(a)
 b=array(((1,2),(3,4))); b=array([[1,2],[3,4]])
Ausschneiden: x[:3] oder x[-2:]
Interne Funktionen: sum(x) 

Matrix Loesung von Gleichungssystemen http://www.scipy.org/Tentative_NumPy_Tutorial
from numpy import matrix
from numpy import linalg
 c=matrix[-7,5,6]
A=matrix([[1,2,2],[2,-2,1],[2,1,-2]])
C=matrix([[2,-1,-1],[1,0,2],[1,2,0]])
 C.I*A*C
a=matrix([[1,2,3],[2,1,4],[0,1,0]])
det(a)
c.T
c.I
Transponiere:
linalg.solve(a, c.T)

A=matrix([1,2,3])
Erzeuge symmetrische Matix:
A.T+A
Erzeuge schiefsymmetrische Matix:
A.T-A


Vektorprodukt in 3d
from numpy import * oder from pylab import *
cross(u,v)
a=array([1,-2,-1]), b=array([-3,-2,-1]), cross(a,cross(a,b))

Skalarprodukt
dot(u,v)
sum(A[i]*B[i] for i in range(len(A)))

Laenge eines Vektors sqrt(sum(u*u))

Kurvenfit

o=arange(3)+1	% Polynom 2. Ordnung
p=poly1d(o)	%
p(x)		% Polynom mit Argument x
Wurzeln p.r  Koeffizienten p.c Ordnung p.order Koeffizient k-ter Ordnung p[k]
Polynom als Ring p*p, (p**3 + 4) / p, ... help(poly1d)
Methoden: 
N=200; r=randn(N); x=arange(-N/2.,N/2.); plot(x,p(x)+500*r); show()

p_deg=polyfit(x,y,deg)
c_2=polyfit(x,y,2); p_2=poly1d(c_2);plot(x,y,'D',x,p_2(x)); show()

Kovarianz, Hauptkomponenteanalyse

Numerische Integration von odes

from numpy import *
from pylab import *
from scipy import integrate
a = 1.; b = 0.1; c = 1.5; d = 0.75	# Wahl der Parameter 
def dX_dt(X, t=0):			# Definition der rechten Seite
 return array([ a*X[0] -   b*X[0]*X[1] , -c*X[1] + d*b*X[0]*X[1] ])
t = linspace(0, 15,  1000)              # Zeitachse
X0 = array([10, 5])                     # Anfangsbedingungen: 10 rabbits and 5 foxes  

X, infodict = integrate.odeint(dX_dt, X0, t, full_output=True)
infodict['message']                     # >>> 'Integration successful.'
rabbits, foxes = X.T

plot(t, rabbits, 'r-', label='Rabbits')
plot(t, foxes  , 'b-', label='Foxes')
grid()
legend(loc='best')
xlabel('time')
ylabel('population')
title('Evolution of fox and rabbit populations'+str(round(pi,4))+' '+str(pi)[:4])
show()

Pfad fuer Python in der Bash-Shell setzen

Die Datei .bashrc durch das Kommando

__________________________________________________________
#Python
export PYTHONPATH=$PYTHONPATH:$HOME/bin/python
export PYTHONSTARTUP=$HOME/.pythonrc
__________________________________________________________

source ~/.bashrc

.pythonrc 
__________________________________________________________
# startup script for python to enable saving of interpreter history and
# enabling name completion

# import needed modules
import atexit
import os
import readline
import rlcompleter
import sys


# where is history saved
historyPath = os.path.expanduser("~/.pyhistory")

# handler for saving history
def save_history(historyPath=historyPath):
   import readline
   readline.write_history_file(historyPath)

# read history, if it exists
if os.path.exists(historyPath):
   readline.read_history_file(historyPath)

# register saving handler
atexit.register(save_history)

# enable completion
readline.parse_and_bind('tab: complete')

# fancier prompt
sys.ps1 = '>>>'
# cleanup
del os, atexit, readline, rlcompleter, save_history, historyPath

# personal module imports
from numpy import *
from numpy import matlib
import pylab
__________________________________________________________

Weiteres

CA

Pythonmodule kompilieren unter Windows

1. Laden Sie Swig herunter:
http://prdownloads.sourceforge.net/swig/swigwin-1.3.40.zip
und entpacken Sie das Archiv - zBsp. nach
C:\Programme\swig-1.3.40

2. Setzen Sie folgende Umgebungsvariablen:
PYTHON_INCLUDE = C:\Programme\Python25\include
Das Verzeichnis mit Python.h

PYTHON_LIB = C:\Programme\Python25\libs\python25.lib

NUMPY_INCLUDE = C:\Programme\Python25\Lib\site-packages\numpy\core\include
Das Verzeichnis mit numpy\arrayobject.h

SWIG = C:\Programme\swigwin-1.3.40\swig.exe

3. Laden Sie Visual Studio Express C++ 2008 herunter (kostenlos):
http://www.microsoft.com/germany/Express/download/webdownload.aspx

4. Oeffnen Sie die Projektdatei CompPhys.vcproj in Visual Studio

5. Stellen Sie den Compiler auf "Release"

6. Rufen Sie "Projekt - CompPhys Neu erstellen" auf

7. Testen Sie das Modul mittels test.py

Home page