pyraeus Package

pyraeus Package

astro Module

pyraeus.astro.apply_extinction_law(flux, k, ebv=0)

Apply inter-galactic extinction laws given E(B-V) for an SED

Examples

>>> apply_extinction_law(np.array([0.0064, 0.0063, 0.0065]), np.array([ 4.45, 4.43, 4.41]), ebv=1)
array([ 0.00010621,  0.0001065 ,  0.00011192])

pyraeus.astro.compute_int_mag_AB(wave_res, magc_res, filt_res)

pyraeus.astro.compute_pivot_wave(wave_res, filt_res)

pyraeus.astro.do_redshift(z, wave, flux)

Redshift the SED spectra

Examples

>>> do_redshift(0.1, np.array([5000, 5010, 5020]), np.array([0.0064, 0.0063, 0.0065]))
(array([ 5500.,  5511.,  5522.]), array([ 0.00581818,  0.00572727,  0.00590909]))

pyraeus.astro.flux_to_mag_AB(flux, wave, system='cgs')

Compute the AB magnitude from the flux

Examples

>>> print flux_to_mag_AB(5.29695457906e-14)
-15.41

pyraeus.astro.get_cosmo(key=None)

Get astropy cosmology object from global variable. See set_cosmo for more details.

pyraeus.astro.mag_AB_abs(wave, flux, trans)

Compute the absolute magnitude in AB system

See Eq. 3.1 of O. Ilbert thesis

\(M_{AB} = -2.5 * \log_{10}( \frac{\int(F(\lambda,z)T(\lambda)d\lambda)} {\int(T(\lambda)(c/\lambda^2)d\lambda)} ) - M_{AB,0}\)

\(M_{AB} = -2.5 * \log_{10}(F(\nu) / \text{area}) - M_{AB,0}\)

Examples

>>> import numpy as np
>>> print mag_AB_abs(np.array([5000, 5010, 5020]), np.array([0.0064, 0.0063, 0.0065]), np.array([1, 1, 1]))
-15.418265234

pyraeus.astro.mag_AB_to_flux(mab)

Compute the flux from the AB magnitude

Examples

>>> print mag_AB_to_flux(-15.41)
5.29695457906e-14

pyraeus.astro.set_cosmo(key)

Set global cosmology variable. Possible values are ‘default’ (H0=70, Om0=0.3) or ‘WMAP7’.

catalog Module

class pyraeus.catalog.Catalog(mag, err_mag, zref=None)

Bases: object

A catalog of observed objects.

Build a catalog of observed objects, from an array of magnitudes and associated errors.

Example

>>> import numpy as np
>>> C = Catalog(np.array([(17.696,  17.261,  17.265), (20.524,  18.928,  17.865)]), np.array([(0.01,  0.01,  0.01), (0.01,  0.01,  0.01)]) )
>>> print C.mag
[[ 17.696  17.261  17.265]
 [ 20.524  18.928  17.865]]

fit_redshift(lib, mode='broadcast', z_interp=False)

Fit redshift using an input template library.

Best fit redshift and associated chi2 are stored in the zphot and chi2 attributes. PDF(z) is stored in the pdz attribute.

static load(filename, mag_col, err_col)

Load catalog from ASCII file. Magnitudes and associated error are read from list of column index.

save(output_fn, lst_filt, print_mags=False, output_fmt='commented_header')

Saving outputs of the catalog. ascii fmt. Output formats from astropy.io (default fixed_width) Output variables : ‘obj ID’,’phot z’, ‘best chi2’, ‘SED id’

Optionnal: mags mag_err output

save_stat(filename, data)

Save specific data to disk

extlaw Module

class pyraeus.extlaw.ExtLaw(filename, wrange=None)

Bases: object

An ExtLaw contains an extinction law function k(wavelength)

wavelength in angstrom

Load extinction law from file (ASCII 2 columns, [w,h]).

If wrange is set, SED is resampled to match the input wavelength range.

Example

>>> print ExtLaw ("ancil_data/ext/calzetti.dat").k
[   0.       190.18576  160.93358 ...,   -0.19588   -0.19622   -0.19657]

filt Module

class pyraeus.filt.Filter(filename, wrange=None, name=None)

Bases: object

An filter contains a window function (between 0 and 1), denoted \(T \left( \lambda \right)\)

Load filter from file (ASCII 2 columns, [w,t]).

If wrange is set, filter is resampled to match the input wavelength range.

Example

>>> filt = Filter ("ancil_data/filter/Yp.pb")
>>> print filt.w[15:20]
[ 9500.  9600.  9700.  9800.  9900.]
>>> print filt.t[15:20]
[ 0.37830019  0.41279066  0.44283501  0.47776948  0.50711225]

pyraeus.filt.load_filter(lst_file, wrange=None, numref=0)

Load a list of filter files.

The filter are resampled to match a common wavelength range (numref filter taken as default).

Example

>>> lst_filt = load_filter (["ancil_data/filter/Jp.pb", "ancil_data/filter/Yp.pb" ], numref=1)
>>> print lst_filt[0].w[15:20]
[ 9500.  9600.  9700.  9800.  9900.]
>>> print lst_filt[1].w[15:20]
[ 9500.  9600.  9700.  9800.  9900.]

fit Module

pyraeus.fit.chi2(obs_flux, obs_eflux, temp_flux)

Compute the chi2 of template flux wrt the observed flux and error.

\(\chi^2 = \sum_i( \frac{F_{obs,i}-sF_{temp,i}} {\sigma_i})^2\)

with

\(s = \sum_j( \frac{F_{obs,j}F_{temp,j}} {\sigma_j^2}) / \sum_j( \frac{F_{temp,j}^2} {\sigma_j^2})\)

Parameters obs_flux, obs_eflux, temp_flux: 1d array of same size

Returns scalar

Example

>>> import numpy as np
>>> print chi2(np.arange(3), np.ones((3)), np.arange(0,6,2))
0.0

pyraeus.fit.chi2_broadcast(obs_flux, obs_eflux, temp_flux)

pyraeus.fit.chi2_inline(obs_flux, obs_eflux, temp_flux)

Compute the chi2 for a list of templates, in C.

Parameters obs_flux, obs_eflux: 1d array of same size temp_flux: 2d array, with number of columns equal to obs_flux size.

Returns 1d array

pyraeus.fit.cross_prod(*args)

Return the cross product of 2 lists as a list of tuples.

pyraeus.fit.run_minuit(obj, fname, guess)

Run minuit to minimize obj.fname (theta) with theta=guess as a starting point.

Parameters obj: an object for storing ancillary data fname: the object method name to minimize guess: a vector of parameter value corresponding to the starting point, used as the arguments of the fname criterion to minimize

Return m: a Minuit object parfit: the best fit parameter values parerr: errors on best fit parameter values

library Module

class pyraeus.library.TemplateLibrary(lst_sed, lst_filt, lst_madau=None, zmin=0, zmax=6, zstep=0.1, set_mag=False)

Bases: object

A magnitude library is a table with columns corresponding to magnitudes of each filter and rows corresponding to various sed model and redshift values.

Build a magnitude library.

By default, the magnitude table is only set to nan. If set_mag is set to True, magnitudes are computed and recorded.

Example

>>> from pyraeus import load_sed,load_filter
>>> lst_sed  = load_sed (["ancil_data/sed/Ell1_A_0.sed", "ancil_data/sed/Ell2_A_0.sed"])
>>> lst_filt = load_filter (["ancil_data/filter/Jp.pb", "ancil_data/filter/Yp.pb"],wrange=lst_sed[0].w)
>>> lib = TemplateLibrary (lst_sed, lst_filt)
>>> print lib.rec[0:5,1]
[ nan  nan  nan  nan  nan]

idx(zval, sedid)

Return the array index of a given z and sed.

Example

>>> from pyraeus import load_sed,load_filter
>>> lst_sed  = load_sed (["ancil_data/sed/Ell1_A_0.sed", "ancil_data/sed/Ell2_A_0.sed"])
>>> lst_filt = load_filter (["ancil_data/filter/Jp.pb", "ancil_data/filter/Yp.pb"],wrange=lst_sed[0].w)
>>> lib = TemplateLibrary (lst_sed, lst_filt, set_mag=True) 
>>> z   = 4
>>> sedid = lst_sed[1].id
>>> ii = lib.idx (z, sedid)
>>> print lib.rec[ii,:]
[ 34.36351145  35.3524804 ]

static load(filename)

Load object from disk.

Parameters

filaneme: pickle filename

save(filename, fmt='pkl')

Save object to disk.

Parameters

filaneme: string , including extension fmt: format option in [‘pkl’, ‘bin’, ‘ascii’]

Examples

>>> from pyraeus import load_sed,load_filter
>>> lst_sed  = load_sed (["ancil_data/sed/Ell1_A_0.sed", "ancil_data/sed/Ell2_A_0.sed"])
>>> lst_filt = load_filter (["ancil_data/filter/Jp.pb", "ancil_data/filter/Yp.pb"],wrange=lst_sed[0].w)
>>> lib = TemplateLibrary (lst_sed, lst_filt)
>>> lib.save ("lib.pkl")
>>> lib = TemplateLibrary.load ("lib.pkl")
>>> lib.save ("lib.txt", fmt='ascii')
>>> lib.save ("lib.bin", fmt='bin')

set_mag()

Compute and record magnitude.

Example

>>> from pyraeus import load_sed,load_filter
>>> lst_sed  = load_sed (["ancil_data/sed/Ell1_A_0.sed", "ancil_data/sed/Ell2_A_0.sed"])
>>> lst_filt = load_filter (["ancil_data/filter/Jp.pb", "ancil_data/filter/Yp.pb"],wrange=lst_sed[0].w)
>>> lib = TemplateLibrary (lst_sed, lst_filt,zmin=0, zmax=2, zstep=1)
>>> lib.set_mag() 
>>> print lib.rec[0:5,1]
[ 27.79665888  29.7357644   27.63744741  29.70090341]

madaulaw Module

class pyraeus.madaulaw.LstMadauLaw(lst_files, lst_z, wrange=None, name=None)

Bases: object

An LstMadauLaw contains the Madau law function exp(-tau(lambda_rest)) for several z

Load Madau laws from list of files and redhsift.

get_opacity(z)

Return nearest z opacity law.

class pyraeus.madaulaw.MadauLaw(filename, z, wrange=None, name=None)

Bases: object

An MadauLaw contains the Madau law function exp(-tau(lambda_rest)) for a given z

rest wavelength in angstrom

Load Madau law from file (ASCII 2 columns, [w,exp(-tau)]).

If wrange is set, opacity is resampled to match the input wavelength range.

Example

>>> print MadauLaw ("ancil_data/opa/tau02.out").etau
[ 1.        1.        1.       ...,  0.993293  0.993275  0.993258]

perfphotoz Module

pyraeus.perfphotoz.acuracy(delta_z_wo)

Comput the redshift accuracy on the data without outliers (?)

1. Root-mean-square (RMS) scatter: sigma_rms = standard_deviation(delta_z_wo)

2. Normalized median absolute deviation (NMAD): sigma_nmad = 1.48 * median(abs(delta_z_norm))

pyraeus.perfphotoz.bias(delta_z_wo)

Comput the redshift bias on the data without outliers

e_z = abs(median(delta_z_wo))

pyraeus.perfphotoz.error(z_cat, z_reco)

Comput the redshift error

1. Formula from Hildebrandt et al., AA, 523, A31 (2010):

   delta_z = z_cat - z_reco

2. Formula from Ilbert et al., AA, 457, 841 (2006):

   delta_z = delta_z / (1. + z_cat)

Examples

>>> import numpy as np
>>> print error(np.array([1.0,1.4]),np.array([1.1,1.5]))
(array([-0.1, -0.1]), array([-0.05      , -0.04166667]))

pyraeus.perfphotoz.outliers(delta_z, rate_outliers)

Comput the rate of outliers in %

outliers = abs(delta_z) > 0.10 - 0.15 eta = outliers * 100 / nb_outliers

pyraeus.perfphotoz.reject_outliers(z_cat, delta_z, rate_outliers)

Reject the outlier objects from the data

See the function outliers(arg1,arg2) for formula.

sed Module

class pyraeus.sed.Sed(filename, wrange=None)

Bases: object

An SED contains a list of fluxes (erg/s/cm^2/angstrom) associated to a wavelength (angstrom), denoted \(F \left( \lambda \right)\)

Load SED from file (ASCII 2 columns, [w,h]).

If wrange is set, SED is resampled to match the input wavelenght range. The interpolation is performed on a log/log scale.

Example

>>> sed0 = Sed ("ancil_data/sed/Ell1_A_0.sed")
>>> print sed0.w[500:505]
[ 1000.    1002.31  1004.62  1006.93  1009.25]
>>> print sed0.h[500:505]
[  2.13015000e-05   2.16627000e-05   2.20212000e-05   2.21186000e-05
   2.21636000e-05]
>>> sed1 = Sed ("ancil_data/sed/Ell1_A_0.sed", wrange=range(1000,1008,2))
>>> print sed1.h[0:5]
[  2.13015000e-05   2.16139224e-05   2.19244821e-05   2.20793622e-05]

add_extlaw(filename, ebv)

Add extension law attribut

pyraeus.sed.load_sed(lst_file, wrange=None, refnum=0)

Load a list of SED from files.

The SED are resampled to match a common wavelength range (first sed taken as default). The interpolation is performed on a log/log scale.

Example

>>> lst_sed = load_sed (["ancil_data/sed/Ell1_A_0.sed", "ancil_data/sed/SunLCB.sed"])
>>> print lst_sed[0].w[500:505]
[ 1000.    1002.31  1004.62  1006.93  1009.25]
>>> print lst_sed[1].w[500:505]
[ 1000.    1002.31  1004.62  1006.93  1009.25]
>>> lst_sed = load_sed (["ancil_data/sed/Ell1_A_0.sed", "ancil_data/sed/SunLCB.sed"], wrange=range(1000,1008,2))
>>> print lst_sed[0].w[0:4]
[ 1000.  1002.  1004.  1006.]
>>> lst_sed = load_sed (["ancil_data/sed/Ell1_A_0.sed", "ancil_data/sed/SunLCB.sed"], refnum=1)
>>> print lst_sed[0].w[500:505]
[ 6370.  6390.  6410.  6430.  6450.]

tools Module

pyraeus.tools.extrap1d(interpolator)

Return an extrapolator

from http://stackoverflow.com/questions/2745329/how-to-make-scipy-interpolate-give-an-extrapolated-result-beyond-the-input-range

Example

>>> x = np.arange(0,10)
>>> y = np.exp(-x/3.0)
>>> f_i = interp1d(x, y)
>>> f_x = extrap1d(f_i)
>>> print f_x([9,10])
[ 0.04978707  0.03009069]

pyraeus.tools.load_param_file(filename)

Return parameters keys and values as a dictionnary

Example

>>> print load_param_file ("test/example.conf")
{'h0': '70. # Hubble constante at present time in km.s-1.Mpc-1', 'om0': '0.3 # Omega matter', 'sed': ['ancil_data/sed/Ell1_A_0.sed', 'ancil_data/sed/Ell2_A_0.sed'], 'l0': '0.7 # Omega lambda'}

pyraeus.tools.resample(w, h, wrange, log=False, union=False, extrapol=False)

Resample a function of wavelength \(h \left( \lambda \right)\) to a given wavelength range

If union is True, the wavelength range is the union of the two ranges. If log option is set to True, the interpolation is performed on the log/log scale.

Example

>>> w,h = resample ( np.arange(0,10,2), np.arange(0,10,2), np.arange(0,10,1))
>>> print h
[  0.   1.   2.   3.   4.   5.   6.   7.   8.  nan]
>>> w,h = resample ( np.arange(0,10,2), np.arange(0,10,2), np.arange(4,10,1), union=True)
>>> print h
[  0.   2.   4.   5.   6.   7.   8.  nan]
>>> w,h = resample ( np.array([1,5,10]), np.array([1,10,100]), np.arange(1,10,2), log=True)
>>> print h
[  1.           4.8151098   10.          30.57924984  70.46880495]