Code to randomly insert DLAs into a sightline [v1.0]



In [50]:

    
# imports
from scipy import interpolate

from specdb.specdb import IgmSpec

from pyigm.surveys.dlasurvey import DLASurvey
from pyigm.fN import dla as pyi_fd
from pyigm.fN.fnmodel import FNModel

#from xastropy.xutils import xdebug as xdb

sys.path.append(os.path.abspath("../../src"))
import training_set as tset

Random numbers



In [43]:

    
rstate = np.random.RandomState(1234)

Spectra



In [3]:

    
igmsp = IgmSpec()









    



Using /u/xavier/local/Python/igmspec/DB/IGMspec_DB_v02.hdf5 for the catalog file
Using /u/xavier/local/Python/igmspec/DB/IGMspec_DB_v02.hdf5 for the DB file
Available surveys: [u'BOSS_DR12', u'HSTQSO', u'SDSS_DR7', u'KODIAQ_DR1', u'MUSoDLA', u'HD-LLS_DR1', u'2QZ', u'ESI_DLA', u'HDLA100', u'GGG', u'COS-Halos', u'HST_z2', u'COS-Dwarfs', u'XQ-100']
Database is igmspec
Created on 2016-Oct-25

Grab the sightlines



In [2]:

    
reload(tset)
slines, sdict = tset.grab_sightlines()









    



Using the DR5 sample for the sightlines
SDSS-DR5: Loading DLA file /Users/xavier/local/Python/pyigm/pyigm/data/DLA/SDSS_DR5/dr5_alldla.fits.gz
SDSS-DR5: Loading QSOs file /Users/xavier/local/Python/pyigm/pyigm/data/DLA/SDSS_DR5/dr5_dlagz_s2n4.fits
We have 5034 sightlines for analysis
Min z = 2.23244, Median z = 2.65984, Max z = 5.19597

Fiddling with random $z$

z~3



In [8]:

    
i3 = np.argmin(np.abs(slines['ZEM']-3.))
s3 = slines[i3]
spec3l, meta = igmsp.spec_from_coord((s3['RA'], s3['DEC']), isurvey=['SDSS_DR7'])
spec3 = spec3l[0]









    



Your search yielded 1 match[es]
Staged 1 spectra totalling 6.4e-05 Gb
Loaded spectra



In [11]:

    
spec3.wvmin.value/1215.67 - 1.









    Out[11]:





2.1209319596016654

zlya



In [48]:

    
zlya = spec3.wavelength.value/1215.67 - 1
dz = np.roll(zlya,-1)-zlya
dz[-1] = dz[-2]
dz









    Out[48]:





array([ 0.0007187 ,  0.00071887,  0.00071903, ...,  0.00174   ,
        0.0017404 ,  0.0017404 ])



In [49]:

    
### Cut on zem and 920A rest-frame
gdz = (zlya < s3['ZEM']) & (spec3.wavelength > 910.*u.AA*(1+s3['ZEM']))

$\ell(z)$



In [25]:

    
reload(pyi_fd)
lz = pyi_fd.lX(zlya[gdz], extrap=True, calc_lz=True)
lz









    Out[25]:





array([ 0.18331483,  0.18334042,  0.18336602, ...,  0.27237858,
        0.27241322,  0.27244785])

Cumul



In [37]:

    
cum_lz = np.cumsum(lz*dz[gdz])
tot_lz = cum_lz[-1]
tot_lz
#xdb.xplot(cum_lz)









    Out[37]:





0.18932704664780059

Random draw



In [39]:

    
ndla = np.random.poisson(tot_lz, 20)
ndla









    Out[39]:





array([0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1])

Random redshifts



In [41]:

    
fzdla = interpolate.interp1d(cum_lz/tot_lz, zlya[gdz],bounds_error=False,fill_value=np.min(zlya[gdz]))#



In [46]:

    
randz = rstate.random_sample(100)



In [47]:

    
zdla = fzdla(randz)
zdla









    Out[47]:





array([ 2.98564734,  2.91595431,  2.71686489,  2.95142175,  2.79594765,
        2.78945962,  2.16232933,  2.55063337,  2.4011178 ,  2.45195463,
        2.20491238,  2.82471283,  2.61489647,  2.74420171,  2.71902585,
        2.48622071,  2.55734776,  2.52803089,  2.58838005,  2.98862937,
        2.48859465,  2.13300936,  2.94525584,  2.91020805,  2.46795376,
        2.71887703,  2.47375047,  2.33323137,  2.34364447,  2.53217058,
        2.19429392,  2.73584985,  2.82060163,  2.95257784,  2.49142715,
        2.41543319,  2.48718822,  2.87657154,  2.93031168,  2.53935166,
        2.12247379,  2.16322234,  2.26281261,  2.69246489,  2.22777802,
        2.41565617,  2.17973902,  2.71164358,  2.14358104,  2.89692227,
        2.42448172,  2.82306824,  2.66294105,  2.15939305,  2.66487041,
        2.74242401,  2.30413911,  2.85346724,  2.70425838,  2.65633302,
        2.71835043,  2.76471922,  2.36052721,  2.84606566,  2.89796926,
        2.79835032,  2.76974259,  2.58098138,  2.89826008,  2.74722291,
        2.45040503,  2.62301031,  2.50877492,  2.15217171,  2.67452064,
        2.89337549,  2.96633594,  2.7621093 ,  2.197549  ,  2.523326  ,
        2.16315804,  2.29508922,  2.6647718 ,  2.29379411,  2.66580953,
        2.61922984,  2.90702007,  2.79239517,  2.93207794,  2.39681872,
        2.37655071,  2.77573314,  2.93062102,  2.98888596,  2.72771866,
        2.44718162,  2.69952   ,  2.8724076 ,  2.71709734,  2.23917044])

Fiddling with random NHI

Load $f(N)$ [just for the shape]



In [54]:

    
fN_model = FNModel.default_model()









    



Using P14 spline values to generate a default model
Loading: /Users/xavier/local/Python/pyigm/pyigm/data/fN/fN_spline_z24.fits.gz

Cumulative



In [57]:

    
lX, cum_lX, lX_NHI = fN_model.calculate_lox(fN_model.zpivot,
    20.3,NHI_max=22.5, cumul=True)
xdb.xplot(lX_NHI, cum_lX)

Interpolator



In [58]:

    
fNHI = interpolate.interp1d(cum_lX/cum_lX[-1], lX_NHI,bounds_error=False,fill_value=lX_NHI[0])#



In [59]:

    
randNHI = rstate.random_sample(100)



In [60]:

    
dla_NHI = fNHI(randNHI)
dla_NHI









    Out[60]:





array([ 20.48331731,  21.67512512,  20.59492666,  20.72369414,
        20.67841503,  20.4226547 ,  20.40168945,  20.654418  ,
        21.73114995,  20.56439854,  20.88243825,  20.40832829,
        20.34198151,  20.39645882,  20.39900071,  20.4817214 ,
        20.68728349,  20.96605553,  21.16749841,  20.79019633,
        21.10583448,  20.74783509,  20.8566711 ,  20.68772767,
        21.18949896,  20.33081183,  20.98149793,  20.38846215,
        20.33469602,  20.33487522,  20.33063393,  20.66801997,
        20.57315366,  22.24008286,  20.36533816,  20.51346323,
        21.98394868,  20.57005345,  21.17476473,  20.71862699,
        21.41103761,  21.17982791,  20.97538906,  20.31855475,
        21.46239687,  20.67137241,  20.42941226,  21.15508122,
        20.44042832,  20.38430537,  20.89732151,  21.53947189,
        20.6861296 ,  20.60572458,  20.31279553,  20.75622788,
        21.47579246,  21.49400243,  21.31255748,  20.41408149,
        20.77435049,  20.57955718,  20.54032425,  20.32141439,
        21.13921647,  20.83101798,  20.35225833,  21.18377088,
        20.98036536,  20.50603771,  21.04027285,  20.66033865,
        20.31225203,  20.59634708,  20.60699395,  20.45042549,
        21.06201559,  20.92422136,  20.73603896,  20.64716319,
        20.64228378,  20.46978483,  20.75520071,  20.93677031,
        20.47248507,  20.67383392,  20.35896704,  20.88194592,
        21.17649669,  20.71656616,  21.83211643,  20.41018322,
        20.30694166,  20.61601168,  21.31888774,  21.09247831,
        20.63002161,  21.22513789,  20.44167876,  20.49645255])



In [63]:

    
rstate.rand?

Debugging



In [1]:

    
from linetools.spectra import io as lsio
spec = lsio.readspec('/Users/xavier/Dropbox/MachineLearning/DLAs/training_20456_5000.hdf5')



In [2]:

    
spec.select = 2
spec.plot(xspec=True)









    



Out of bounds
button=1, x=379.000000, y=362.000000, xdata=5121.696774, ydata=13.2475






    



WARNING: UnitsWarning: The unit 'Angstrom' has been deprecated in the FITS standard. Suggested: 10**-1 nm. [astropy.units.format.utils]
/Users/xavier/local/Python/linetools/linetools/lists/linelist.py:374: RuntimeWarning: divide by zero encountered in log10
  self._data['log(w*f)'] = np.log10(qm_strength)






    



Loading abundances from Asplund2009
Abundances are relative by number on a logarithmic scale with H=12
button=3, x=630.000000, y=418.000000, xdata=5088.294212, ydata=16.1189
You chose: HI 1215 :: 1215.67 Angstrom :: 0.416
button=1, x=264.000000, y=161.000000, xdata=4556.757794, ydata=2.9413
button=3, x=264.000000, y=161.000000, xdata=4556.757794, ydata=2.9413
You chose: HI 1215 :: 1215.67 Angstrom :: 0.416
button=3, x=299.000000, y=244.000000, xdata=3932.275116, ydata=7.1971
You chose: HI 1215 :: 1215.67 Angstrom :: 0.416



In [4]:

    
hdrlist = spec.meta['headers']



In [13]:

    
meta2 = hdrlist[1]
meta2['zem']









    Out[13]:





3.744961977005005



In [8]:

    
import h5py
import json



In [7]:

    
hdf = h5py.File('/Users/xavier/Dropbox/MachineLearning/DLAs/training_20456_5000.hdf5','r')



In [9]:

    
meta = json.loads(hdf['meta'].value)



In [11]:

    
meta['headers'][2]









    Out[11]:





{u'AU': 0.086,
 u'BESTFLAG': 3,
 u'BESTID': u'588017949356327018',
 u'CAMCOL': 5,
 u'DATE-OBS': u'2004-05-13',
 u'DEC': 42.471892,
 u'DGMI': -0.056,
 u'EPOCH': 2000.0,
 u'FIBER': 504,
 u'FIELD': 240,
 u'FIRSTDEL': 0.0,
 u'FIRSTMAG': 0.0,
 u'FIRSTSN': 0.0,
 u'FTFLAG': 0,
 u'FTTFLAG': 0,
 u'GMAG': 19.062,
 u'GMAGERR': 0.014,
 u'GRATING': u'BOTH',
 u'GTMAG': 19.101,
 u'GTMAGERR': 0.019,
 u'GXFLAG': 0,
 u'GXTFLAG': 0,
 u'HMAG': 0.0,
 u'HMAGERR': 0.0,
 u'HZFLAG': 1,
 u'HZTFLAG': 1,
 u'IGM_ID': 185041,
 u'IMAG': 18.768,
 u'IMAGERR': 0.03,
 u'INSTR': u'SDSS',
 u'ITMAG': 18.728,
 u'ITMAGERR': 0.038,
 u'JMAG': 17.551,
 u'JMAGERR': 0.323,
 u'KMAG': 0.0,
 u'KMAGERR': 0.0,
 u'LZFLAG': 1,
 u'LZTFLAG': 1,
 u'MASSDEL': 0.0,
 u'MASSFLG': 18,
 u'MFLAG': 0,
 u'MIMAG': -27.67,
 u'MODEFLAG': 1,
 u'NPIX': 3848,
 u'OBJECT': 106,
 u'ONAME': u' SDSS J131252.22+422818.7 ',
 u'PLATE': 1460,
 u'R': 2000.0,
 u'RA': 198.21759,
 u'RASSCNT': -9.0,
 u'RASSDEL': 0.0,
 u'RASSSN': 0.0,
 u'RERUN': 41,
 u'RMAG': 18.806,
 u'RMAGERR': 0.021,
 u'RMJD': 52754,
 u'RTFLAG': 0,
 u'RTMAG': 18.81,
 u'RTMAGERR': 0.023,
 u'RTTFLAG': 0,
 u'RUN': 3893,
 u'SDSSJ': u'131252.22+422818.8',
 u'SMJD': 53138,
 u'SPECOID': u'411181694083661824',
 u'SPEC_FILE': u'spSpec-53138-1460-504.fit.gz',
 u'SPFLAG': 1,
 u'SRFLAG': 0,
 u'SRTFLAG': 0,
 u'STFLAG': 0,
 u'STTFLAG': 0,
 u'SURVEY_ID': 66057,
 u'TELESCOPE': u'SDSS 2.5-M',
 u'TFLAG': 3,
 u'UMAG': 20.633,
 u'UMAGERR': 0.072,
 u'USELFLAG': 1,
 u'UTMAG': 20.683,
 u'UTMAGERR': 0.067,
 u'WV_MAX': 9193.904748039808,
 u'WV_MIN': 3791.4034418330666,
 u'ZMAG': 18.789,
 u'ZMAGERR': 0.046,
 u'ZTMAG': 18.828,
 u'ZTMAGERR': 0.043,
 u'flag_zem': u'SDSS-HW',
 u'logNH': 20.164,
 u'sig_zem': 0.0,
 u'zem': 3.1747050285339355}



In [ ]: