Python programozási segédlet

Ajánlott programcsomag: Anaconda, ebből a Spyder környezetet érdemes használni, mely egy interaktív python shell. A Python programozási nyelv és környezet több implementációja is létezik, ebből mi a "hivatalos" kiadást használjuk, mert ehhez adottak a szükséges numerikus csomagok (numpy, scipy, matplotlib). Ezen felül a Python két verziója is elterjedt (jelenleg 2.7, illetve 3.4), melyek szintaktikailag különbözőek. Aki most kezdi a Python programozást, érdemes rögtön a 3-as verzióval kezdeni, de előfordulhat, hogy mások által írt scriptek futtatásához 2-es verzióra van szükség. Ezeket adott esetben vagy át kell írni 3-as verzióra, vagy a 2-es verziót is installálni kell. Szerencsére a két verzió vígan megfér egymás mellett.

Alapműveletek

Csomagok inicializálása

import pylab
import numpy
import matplotlib.pyplot as plt

Adatfájl betöltése

gal = pylab.loadtxt("galaxies.csv",skiprows=1,delimiter=',') 

Adatfeldolgozáshoz nem a python beépített tömbjeit érdemes használni, hanem a numpy array osztályát. Az array osztály támogatja a pythonból ismert szeletelő műveleteket is.

z = numpy.array(gal[:,1])
m = numpy.array(gal[:,2]) 

Adatok ábrázolása

pylab.plot(z,m,'.',markersize=1) 

Adatok szűrése:

z = z[m<limit_m]
m = m[m<limit_m] 

Hisztogram készítése

hist,bins = numpy.histogram(M, weights=w, bins=50)
plt.plot(bins[0:-1],hist)

Függvényillesztés

Függvények definiálása

def schechter(x, phi, M, alpha):
    a = 0.4 * (M - x)
    b = math.pow(10, a)
    return 0.4 * math.log(10) * phi * math.pow(b, alpha + 1) * math.exp(-b)

Nem-lineáris függvény illesztése

from scipy.optimize import curve_fit
popt,pcov = curve_fit(log_schechter, bins, hist, p, sigma=error) 

Numerikus integrálás

Függvény numerikus integrálásához (kvadratúra) először szükséges a függvény definitálása:

def ez(z):
    O_M = 0.3
    O_K = 0.0
    O_L = 0.7
    return 1 / math.sqrt((1 + z)**3 * O_M + (1 + z)**2 * O_K + (1 + z) * O_L)
    

Majd erre meghívjuk a szükséges numerikus függvényt. Ez a függvény eredményül nem egyetlen számot, hanem egy kételemű vektort ad vissza, mely első (nulladik) eleme az integrál értéke, a másik elem pedig a becsült numerikus hiba.

import math
from scipy.integrate import quad

def E_z(z):
    O_M = 0.3
    O_K = 0.0
    O_L = 0.7
    return 1 / math.sqrt((1 + z)**3 * O_M + (1 + z)**2 * O_K + (1 + z) * O_L)
    
def D_C(z):
    D_H = 3000 / 0.7
    return D_H * quad(E_z, 0, z)[0]
    
print(D_C(0.1))

Műveletek numpy.array objektumokkal

A numpy csomag erőssége, hogy az elemi műveleteket nem csak skalár értékeken lehet elvégezni, hanem teljes vektorokon, mátrixokon stb. is.

import numpy

a = numpy.array([0.1, 0.2, 0.3])
b = 2 * a
print(b)

Természetesen a fenti megvalósítható egy ciklussal is, de a numpy csomag nagy része nem Python nyelven, hanem C-ben íródott, így nagyságrendekkel gyorsabb, mintha Python ciklusokat használnánk. Az alábbi program működik, de ilyet ne használjunk.

a = [0.1, 0.2, 0.3]
b = [0,0,0]
for i in range(0, len(a)):
    b[i] = 2 * a[i]
print(b)

Egy fontos probléma, hogy miként lehet egy skalár értékű, skalár paraméterű függvényt egy teljes adatvektor minden elemére meghívni úgy, hogy kerüljük a for ciklusok használatát. A célra használhatjuk a numpy.vectorize függvényét, melynek paraméterül magát a skalár függvényt kell meadni (zárójelek nélkül!).

import math
import numpy
from scipy.integrate import quad

def E_z(z):
    O_M = 0.3
    O_K = 0.0
    O_L = 0.7
    return 1 / math.sqrt((1 + z)**3 * O_M + (1 + z)**2 * O_K + (1 + z) * O_L)
    
def D_C(z):
    D_H = 3000 / 0.7
    return D_H * quad(E_z, 0, z)[0]

zz = numpy.array([0.1, 0.2, 0.3])
rr = numpy.vectorize(D_C)(zz)
print(rr)

Ha a vektorizált függvényt gyakran használjuk, akkor érdemes egy függvényértékű változóban eltárolni:

D_C_vec = vectorize(D_C)

zz = numpy.array([0.1, 0.2, 0.3])
rr = D_C_vec(zz)
print(rr)

Az alábbi összetett példa bemutatja, hogy hogyan kell többparaméteres függvényeket csak az egyik paraméter szerint integrálni, illetve többparaméteres függvényeket vektorizálni.

import math
import numpy
from scipy.integrate import quad

def E_z(z, O_M = 0.3, O_K = 0.0, O_L = 0.7):
    return 1 / math.sqrt((1 + z)**3 * O_M + (1 + z)**2 * O_K + (1 + z) * O_L)
    
def D_C(z, O_M = 0.3, O_K = 0.0, O_L = 0.7, H_0 = 70):
    D_H = 3000 / H_0 * 100
    g = lambda x: E_z(x, O_M, O_K, O_L)
    return D_H * quad(g, 0, z)[0]  

zz = numpy.array([0.1, 0.2, 0.3])
D_C_vec = numpy.vectorize(D_C)
print(numpy.vectorize(D_C)(zz, 0.28, 0, 0.72, 69))

Kozmológiai függvények

Az astropy.cosmology csomag számos standard kozmológiát tartalmaz a legfontosabb mérésekből (WMAP, Planck) származó paraméterekkel.

from astropy.cosmology import WMAP9 as cosmo
DM = cosmo.distmod(z).value
M = numpy.subtract(m, DM.value)