From e21299218d4089140a18b9fbb8873a56ec18eb7b Mon Sep 17 00:00:00 2001 From: "J. Xavier Prochaska" Date: Mon, 1 Dec 2014 11:09:31 -0800 Subject: [PATCH 01/18] first steps by JXP --- py/qso_template/__init__.py | 1 + py/qso_template/qso_pca.py | 125 ++++++++++++++++++++++++++++++++++++ 2 files changed, 126 insertions(+) create mode 100644 py/qso_template/__init__.py create mode 100644 py/qso_template/qso_pca.py diff --git a/py/qso_template/__init__.py b/py/qso_template/__init__.py new file mode 100644 index 000000000..7f366e162 --- /dev/null +++ b/py/qso_template/__init__.py @@ -0,0 +1 @@ +import qso_pca diff --git a/py/qso_template/qso_pca.py b/py/qso_template/qso_pca.py new file mode 100644 index 000000000..b8258cc3b --- /dev/null +++ b/py/qso_template/qso_pca.py @@ -0,0 +1,125 @@ +""" +#;+ +#; NAME: +#; qso_pca +#; Version 1.0 +#; +#; PURPOSE: +#; Module for generate QSO PCA templates +#; 24-Nov-2014 by JXP +#;- +#;------------------------------------------------------------------------------ +""" +from __future__ import print_function, absolute_import, division, unicode_literals + +import numpy as np +import os + +from astropy.io import fits +from xastropy.xutils import xdebug as xdb + +""" +## Tinkering about +# JXP on 01 Dec 2014 + +hdu = fits.open('spEigenQSO-55732.fits') +eigen = hdu[0].data # There are 4 eigenvectors +wav = 10.**(2.6534 + np.arange(13637)*0.0001) # Rest-frame, Ang +xdb.xplot(wav,eigen[0,:]) + # Looks sensible enough + +# QSO sample used to generate the PCA eigenvectors +qsos = hdu[1].data +xdb.xhist(qsos['REDSHIFT']) # Skewed to high z +""" + + +""" +Let's test in IDL on a BOSS spectrum + +cd /Users/xavier/BOSS/DR10/BOSSLyaDR10_spectra_v2.1/5870 + +IDL> boss_objinf, [5870,1000], /plot + Name RA(2000) DEC(2000) Mag(R) zobj +5870-56065-100 14:07:37.28 +21:48:38.5 21.06 2.180 + +Looks like a fine quasar + +Read the eigenspectra + +eigen = xmrdfits(getenv('IDLSPEC2D_DIR')+'/templates/spEigenQSO-55732.fits',0,head) +eigen_wave = 10.^(sxpar(head,'COEFF0') + dindgen(sxpar(head,'NAXIS1'))*sxpar(head,'COEFF1')) + + +BOSS spectrum + +boss_readspec, 'speclya-5870-56065-1000.fits', flux, wave, ivar=ivar +zqso = 2.180 +wrest = wave / (1+zqso) +npix = n_elements(wrest) + +Are they 'registered'? +mn = min(abs(wrest[0]-eigen_wave),imn) +x_splot, findgen(npix), flux, ytwo=eigen[imn:*,0]/5292*4.4 + Yes, they are.. + +Build the inputs to computechi2 + +objflux = flux +sqivar = ivar +imx = npix+imn-1 +starflux = eigen[imn:imx,*] + +chi2 = computechi2(objflux, sqivar, starflux, acoeff=acoeff, yfit=yfit) + ;; The fit looks ok + +IDL> print, acoeff + 0.00045629701 -0.00034376233 0.00047730298 0.00045008259 +""" + + +# Read Eigenspectrum + +hdu = fits.open(os.environ.get('IDLSPEC2D_DIR')+'/templates/spEigenQSO-55732.fits') +eigen = hdu[0].data +head = hdu[0].header +eigen_wave = 10.**(head['COEFF0'] + np.arange(head['NAXIS1'])*head['COEFF1']) + +#xdb.xplot(eigen_wave, eigen[0,:]) +#xdb.set_trace() + +#BOSS spectrum + +boss_hdu = fits.open('/Users/xavier/BOSS/DR10/BOSSLyaDR10_spectra_v2.1/5870/speclya-5870-56065-1000.fits.gz') +t = boss_hdu[1].data +flux = t['flux'] +wave = 10.**t['loglam'] +ivar = t['ivar'] + +zqso = 2.180 +wrest = wave / (1+zqso) +npix = len(wrest) +imx = npix+imn + +imn = np.argmin(np.abs(wrest[0]-eigen_wave)) +eigen_flux = eigen[:,imn:imx] + + +# Generate the matrices +""" +Now Python with simple approach following Dan F-M http://dan.iel.fm/emcee/current/user/line/ +""" + +C = np.diag(1./ivar) +A = eigen_flux.T +y = flux + +cov = np.linalg.inv(np.dot(A.T, np.linalg.solve(C, A))) +acoeff = np.dot(cov, np.dot(A.T, np.linalg.solve(C, y))) + +#In [48]: acoeff +#Out[48]: array([ 0.00045461, -0.00012764, 0.00010843, 0.0003281 ], dtype=float32) + +model = np.dot(A,acoeff) +xdb.xplot(wrest, flux, model) + # Looks good to me From d6b16d8bf5ca0df87e6e56f56ae73eada2c3b587 Mon Sep 17 00:00:00 2001 From: "J. Xavier Prochaska" Date: Mon, 1 Dec 2014 11:10:20 -0800 Subject: [PATCH 02/18] putting into tests --- py/qso_template/tests.py | 125 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 125 insertions(+) create mode 100644 py/qso_template/tests.py diff --git a/py/qso_template/tests.py b/py/qso_template/tests.py new file mode 100644 index 000000000..b8258cc3b --- /dev/null +++ b/py/qso_template/tests.py @@ -0,0 +1,125 @@ +""" +#;+ +#; NAME: +#; qso_pca +#; Version 1.0 +#; +#; PURPOSE: +#; Module for generate QSO PCA templates +#; 24-Nov-2014 by JXP +#;- +#;------------------------------------------------------------------------------ +""" +from __future__ import print_function, absolute_import, division, unicode_literals + +import numpy as np +import os + +from astropy.io import fits +from xastropy.xutils import xdebug as xdb + +""" +## Tinkering about +# JXP on 01 Dec 2014 + +hdu = fits.open('spEigenQSO-55732.fits') +eigen = hdu[0].data # There are 4 eigenvectors +wav = 10.**(2.6534 + np.arange(13637)*0.0001) # Rest-frame, Ang +xdb.xplot(wav,eigen[0,:]) + # Looks sensible enough + +# QSO sample used to generate the PCA eigenvectors +qsos = hdu[1].data +xdb.xhist(qsos['REDSHIFT']) # Skewed to high z +""" + + +""" +Let's test in IDL on a BOSS spectrum + +cd /Users/xavier/BOSS/DR10/BOSSLyaDR10_spectra_v2.1/5870 + +IDL> boss_objinf, [5870,1000], /plot + Name RA(2000) DEC(2000) Mag(R) zobj +5870-56065-100 14:07:37.28 +21:48:38.5 21.06 2.180 + +Looks like a fine quasar + +Read the eigenspectra + +eigen = xmrdfits(getenv('IDLSPEC2D_DIR')+'/templates/spEigenQSO-55732.fits',0,head) +eigen_wave = 10.^(sxpar(head,'COEFF0') + dindgen(sxpar(head,'NAXIS1'))*sxpar(head,'COEFF1')) + + +BOSS spectrum + +boss_readspec, 'speclya-5870-56065-1000.fits', flux, wave, ivar=ivar +zqso = 2.180 +wrest = wave / (1+zqso) +npix = n_elements(wrest) + +Are they 'registered'? +mn = min(abs(wrest[0]-eigen_wave),imn) +x_splot, findgen(npix), flux, ytwo=eigen[imn:*,0]/5292*4.4 + Yes, they are.. + +Build the inputs to computechi2 + +objflux = flux +sqivar = ivar +imx = npix+imn-1 +starflux = eigen[imn:imx,*] + +chi2 = computechi2(objflux, sqivar, starflux, acoeff=acoeff, yfit=yfit) + ;; The fit looks ok + +IDL> print, acoeff + 0.00045629701 -0.00034376233 0.00047730298 0.00045008259 +""" + + +# Read Eigenspectrum + +hdu = fits.open(os.environ.get('IDLSPEC2D_DIR')+'/templates/spEigenQSO-55732.fits') +eigen = hdu[0].data +head = hdu[0].header +eigen_wave = 10.**(head['COEFF0'] + np.arange(head['NAXIS1'])*head['COEFF1']) + +#xdb.xplot(eigen_wave, eigen[0,:]) +#xdb.set_trace() + +#BOSS spectrum + +boss_hdu = fits.open('/Users/xavier/BOSS/DR10/BOSSLyaDR10_spectra_v2.1/5870/speclya-5870-56065-1000.fits.gz') +t = boss_hdu[1].data +flux = t['flux'] +wave = 10.**t['loglam'] +ivar = t['ivar'] + +zqso = 2.180 +wrest = wave / (1+zqso) +npix = len(wrest) +imx = npix+imn + +imn = np.argmin(np.abs(wrest[0]-eigen_wave)) +eigen_flux = eigen[:,imn:imx] + + +# Generate the matrices +""" +Now Python with simple approach following Dan F-M http://dan.iel.fm/emcee/current/user/line/ +""" + +C = np.diag(1./ivar) +A = eigen_flux.T +y = flux + +cov = np.linalg.inv(np.dot(A.T, np.linalg.solve(C, A))) +acoeff = np.dot(cov, np.dot(A.T, np.linalg.solve(C, y))) + +#In [48]: acoeff +#Out[48]: array([ 0.00045461, -0.00012764, 0.00010843, 0.0003281 ], dtype=float32) + +model = np.dot(A,acoeff) +xdb.xplot(wrest, flux, model) + # Looks good to me From 29712e065c70923c5934d77637e3231d645afc66 Mon Sep 17 00:00:00 2001 From: "J. Xavier Prochaska" Date: Mon, 1 Dec 2014 11:32:52 -0800 Subject: [PATCH 03/18] a bit more --- py/qso_template/tests.py | 28 +++++++++++++++++++++------- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/py/qso_template/tests.py b/py/qso_template/tests.py index b8258cc3b..e2e9a08ad 100644 --- a/py/qso_template/tests.py +++ b/py/qso_template/tests.py @@ -16,7 +16,12 @@ import os from astropy.io import fits -from xastropy.xutils import xdebug as xdb + +flg_xdb = True +try: + from xastropy.xutils import xdebug as xdb +except ImportError: + flg_xdb = False """ ## Tinkering about @@ -80,7 +85,9 @@ # Read Eigenspectrum -hdu = fits.open(os.environ.get('IDLSPEC2D_DIR')+'/templates/spEigenQSO-55732.fits') +eigen_fil = os.environ.get('IDLSPEC2D_DIR')+'/templates/spEigenQSO-55732.fits' +print('Using these eigen spectra: {:s}'.format(eigen_fil)) +hdu = fits.open(eigen_fil) eigen = hdu[0].data head = hdu[0].header eigen_wave = 10.**(head['COEFF0'] + np.arange(head['NAXIS1'])*head['COEFF1']) @@ -90,7 +97,9 @@ #BOSS spectrum -boss_hdu = fits.open('/Users/xavier/BOSS/DR10/BOSSLyaDR10_spectra_v2.1/5870/speclya-5870-56065-1000.fits.gz') +boss_fil = '/Users/xavier/BOSS/DR10/BOSSLyaDR10_spectra_v2.1/5870/speclya-5870-56065-1000.fits.gz' +print('Fitting this spectrum: {:s}'.format(boss_fil)) +boss_hdu = fits.open(boss_fil) t = boss_hdu[1].data flux = t['flux'] wave = 10.**t['loglam'] @@ -99,9 +108,9 @@ zqso = 2.180 wrest = wave / (1+zqso) npix = len(wrest) -imx = npix+imn imn = np.argmin(np.abs(wrest[0]-eigen_wave)) +imx = npix+imn eigen_flux = eigen[:,imn:imx] @@ -114,12 +123,17 @@ A = eigen_flux.T y = flux -cov = np.linalg.inv(np.dot(A.T, np.linalg.solve(C, A))) -acoeff = np.dot(cov, np.dot(A.T, np.linalg.solve(C, y))) +alpha = np.dot(A.T, np.linalg.solve(C, A)) # Numerical Recipe notation +cov = np.linalg.inv(alpha) +beta = np.dot(A.T, np.linalg.solve(C, y)) +acoeff = np.dot(cov, beta) + +print('acoeff = {:g}, {:g}, {:g}, {:g}'.format(*acoeff)) #In [48]: acoeff #Out[48]: array([ 0.00045461, -0.00012764, 0.00010843, 0.0003281 ], dtype=float32) model = np.dot(A,acoeff) -xdb.xplot(wrest, flux, model) +if flg_xdb is True: + xdb.xplot(wrest, flux, model) # Looks good to me From 0b1357306a682fa3963e50ed3e24092103bf2c8b Mon Sep 17 00:00:00 2001 From: "J. Xavier Prochaska" Date: Mon, 1 Dec 2014 19:58:15 -0800 Subject: [PATCH 04/18] jxp: first pca code --- py/qso_template/qso_pca.py | 214 +++++++++++++++++++------------------ 1 file changed, 108 insertions(+), 106 deletions(-) diff --git a/py/qso_template/qso_pca.py b/py/qso_template/qso_pca.py index b8258cc3b..22b7b9100 100644 --- a/py/qso_template/qso_pca.py +++ b/py/qso_template/qso_pca.py @@ -16,110 +16,112 @@ import os from astropy.io import fits -from xastropy.xutils import xdebug as xdb +flg_xdb = True +try: + from xastropy.xutils import xdebug as xdb +except ImportError: + flg_xdb = False + +# +def read_qso_eigen(eigen_fil=None): + ''' + Input the QSO Eigenspectra + ''' + # File + if eigen_fil is None: + eigen_fil = os.environ.get('IDLSPEC2D_DIR')+'/templates/spEigenQSO-55732.fits' + print('Using these eigen spectra: {:s}'.format(eigen_fil)) + hdu = fits.open(eigen_fil) + eigen = hdu[0].data + head = hdu[0].header + # Rest wavelength + eigen_wave = 10.**(head['COEFF0'] + np.arange(head['NAXIS1'])*head['COEFF1']) + + # Return + return eigen, eigen_wave + +## +def fit_eigen(flux,ivar,eigen_flux): + ''' + Fit the spectrum with the eigenvectors. + Pass back the coefficients + ''' + #C = np.diag(1./ivar) + Cinv = np.diag(ivar) + A = eigen_flux.T + + #alpha = np.dot(A.T, np.linalg.solve(C, A)) # Numerical Recipe notation + alpha = np.dot(A.T, np.dot(Cinv,A)) + cov = np.linalg.inv(alpha) + #beta = np.dot(A.T, np.linalg.solve(C, y)) + beta = np.dot(A.T, np.dot(Cinv, flux)) + acoeff = np.dot(cov, beta) + + # Return + return acoeff + +## +def do_boss_lya(debug=False): + ''' + Generate PCA coeff for the BOSS Lya DR10 dataset, v2.1 + ''' + # Eigen + eigen, eigen_wave = read_qso_eigen() + + # Open the BOSS catalog file + boss_cat_fil = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_cat_v2.1.fits.gz' + bcat_hdu = fits.open(boss_cat_fil) + t_boss = bcat_hdu[1].data + nqso = len(t_boss) + + pca_val = np.zeros((nqso, 4)) + + # Loop us -- Should spawn on multiple CPU + #for ii in range(nqso): + datdir = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_spectra_v2.1/' + for ii in range(1000): + if (ii % 100) == 0: + print('ii = {:d}'.format(ii)) + # Spectrum file + pnm = str(t_boss['PLATE'][ii]) + fnm = str(t_boss['FIBERID'][ii]).rjust(4,str('0')) + mjd = str(t_boss['MJD'][ii]) + sfil = datdir+pnm+'/speclya-' + sfil = sfil+pnm+'-'+mjd+'-'+fnm+'.fits.gz' + # Read spectrum + spec_hdu = fits.open(sfil) + t = spec_hdu[1].data + flux = t['flux'] + wave = 10.**t['loglam'] + ivar = t['ivar'] + zqso = t_boss['z_pipe'][ii] + + wrest = wave / (1+zqso) + npix = len(wrest) + + # Pack + imn = np.argmin(np.abs(wrest[0]-eigen_wave)) + imx = npix+imn + eigen_flux = eigen[:,imn:imx] + + # FIT + acoeff = fit_eigen(flux, ivar, eigen_flux) + pca_val[ii,:] = acoeff + + # Check + if debug is True: + model = np.dot(eigen_flux.T,acoeff) + if flg_xdb is True: + xdb.xplot(wrest, flux, model) + xdb.set_trace() + + xdb.set_trace() + +## ################################# +## ################################# +## TESTING +## ################################# +if __name__ == '__main__': + # Run + do_boss_lya()#debug=True) -""" -## Tinkering about -# JXP on 01 Dec 2014 - -hdu = fits.open('spEigenQSO-55732.fits') -eigen = hdu[0].data # There are 4 eigenvectors -wav = 10.**(2.6534 + np.arange(13637)*0.0001) # Rest-frame, Ang -xdb.xplot(wav,eigen[0,:]) - # Looks sensible enough - -# QSO sample used to generate the PCA eigenvectors -qsos = hdu[1].data -xdb.xhist(qsos['REDSHIFT']) # Skewed to high z -""" - - -""" -Let's test in IDL on a BOSS spectrum - -cd /Users/xavier/BOSS/DR10/BOSSLyaDR10_spectra_v2.1/5870 - -IDL> boss_objinf, [5870,1000], /plot - Name RA(2000) DEC(2000) Mag(R) zobj -5870-56065-100 14:07:37.28 +21:48:38.5 21.06 2.180 - -Looks like a fine quasar - -Read the eigenspectra - -eigen = xmrdfits(getenv('IDLSPEC2D_DIR')+'/templates/spEigenQSO-55732.fits',0,head) -eigen_wave = 10.^(sxpar(head,'COEFF0') + dindgen(sxpar(head,'NAXIS1'))*sxpar(head,'COEFF1')) - - -BOSS spectrum - -boss_readspec, 'speclya-5870-56065-1000.fits', flux, wave, ivar=ivar -zqso = 2.180 -wrest = wave / (1+zqso) -npix = n_elements(wrest) - -Are they 'registered'? -mn = min(abs(wrest[0]-eigen_wave),imn) -x_splot, findgen(npix), flux, ytwo=eigen[imn:*,0]/5292*4.4 - Yes, they are.. - -Build the inputs to computechi2 - -objflux = flux -sqivar = ivar -imx = npix+imn-1 -starflux = eigen[imn:imx,*] - -chi2 = computechi2(objflux, sqivar, starflux, acoeff=acoeff, yfit=yfit) - ;; The fit looks ok - -IDL> print, acoeff - 0.00045629701 -0.00034376233 0.00047730298 0.00045008259 -""" - - -# Read Eigenspectrum - -hdu = fits.open(os.environ.get('IDLSPEC2D_DIR')+'/templates/spEigenQSO-55732.fits') -eigen = hdu[0].data -head = hdu[0].header -eigen_wave = 10.**(head['COEFF0'] + np.arange(head['NAXIS1'])*head['COEFF1']) - -#xdb.xplot(eigen_wave, eigen[0,:]) -#xdb.set_trace() - -#BOSS spectrum - -boss_hdu = fits.open('/Users/xavier/BOSS/DR10/BOSSLyaDR10_spectra_v2.1/5870/speclya-5870-56065-1000.fits.gz') -t = boss_hdu[1].data -flux = t['flux'] -wave = 10.**t['loglam'] -ivar = t['ivar'] - -zqso = 2.180 -wrest = wave / (1+zqso) -npix = len(wrest) -imx = npix+imn - -imn = np.argmin(np.abs(wrest[0]-eigen_wave)) -eigen_flux = eigen[:,imn:imx] - - -# Generate the matrices -""" -Now Python with simple approach following Dan F-M http://dan.iel.fm/emcee/current/user/line/ -""" - -C = np.diag(1./ivar) -A = eigen_flux.T -y = flux - -cov = np.linalg.inv(np.dot(A.T, np.linalg.solve(C, A))) -acoeff = np.dot(cov, np.dot(A.T, np.linalg.solve(C, y))) - -#In [48]: acoeff -#Out[48]: array([ 0.00045461, -0.00012764, 0.00010843, 0.0003281 ], dtype=float32) - -model = np.dot(A,acoeff) -xdb.xplot(wrest, flux, model) - # Looks good to me From 47d0a62260dcd4fd6e61ef7021e49bea1bd45a2e Mon Sep 17 00:00:00 2001 From: "J. Xavier Prochaska" Date: Tue, 2 Dec 2014 06:45:27 -0800 Subject: [PATCH 05/18] parallelized --- py/qso_template/qso_pca.py | 54 +++++++++++++++++++++++++++++++++----- 1 file changed, 48 insertions(+), 6 deletions(-) diff --git a/py/qso_template/qso_pca.py b/py/qso_template/qso_pca.py index 22b7b9100..b1cc1d2a3 100644 --- a/py/qso_template/qso_pca.py +++ b/py/qso_template/qso_pca.py @@ -14,6 +14,7 @@ import numpy as np import os +import multiprocessing as mp from astropy.io import fits flg_xdb = True @@ -61,7 +62,7 @@ def fit_eigen(flux,ivar,eigen_flux): return acoeff ## -def do_boss_lya(debug=False): +def do_boss_lya_parallel(istart, iend, output, debug=False): ''' Generate PCA coeff for the BOSS Lya DR10 dataset, v2.1 ''' @@ -74,12 +75,13 @@ def do_boss_lya(debug=False): t_boss = bcat_hdu[1].data nqso = len(t_boss) - pca_val = np.zeros((nqso, 4)) + pca_val = np.zeros((iend-istart, 4)) # Loop us -- Should spawn on multiple CPU #for ii in range(nqso): datdir = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_spectra_v2.1/' - for ii in range(1000): + jj = 0 + for ii in range(istart,iend): if (ii % 100) == 0: print('ii = {:d}'.format(ii)) # Spectrum file @@ -106,7 +108,8 @@ def do_boss_lya(debug=False): # FIT acoeff = fit_eigen(flux, ivar, eigen_flux) - pca_val[ii,:] = acoeff + pca_val[jj,:] = acoeff + jj += 1 # Check if debug is True: @@ -115,13 +118,52 @@ def do_boss_lya(debug=False): xdb.xplot(wrest, flux, model) xdb.set_trace() - xdb.set_trace() + #xdb.set_trace() + output.put((istart,iend,pca_val)) ## ################################# ## ################################# ## TESTING ## ################################# if __name__ == '__main__': + # Run - do_boss_lya()#debug=True) + #do_boss_lya()#debug=True) + # Parallel + boss_cat_fil = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_cat_v2.1.fits.gz' + bcat_hdu = fits.open(boss_cat_fil) + t_boss = bcat_hdu[1].data + nqso = len(t_boss) + nqso = 800 # Testing + + output = mp.Queue() + processes = [] + nproc = 4 + nsub = nqso // nproc + # Setup the Processes + for ii in range(nproc): + # Generate + istrt = ii * nsub + if ii == (nproc-1): + iend = nqso + else: + iend = (ii+1)*nsub + #xdb.set_trace() + process = mp.Process(target=do_boss_lya_parallel, + args=(istrt,iend,output)) + processes.append(process) + + # Run processes + for p in processes: + p.start() + + # Exit the completed processes + for p in processes: + p.join() + + # Get process results from the output queue + print('Got here') + results = [output.get() for p in processes] + # Combine + xdb.set_trace() From 62397d3a1ff500c3b6bb2cc5480d4c0266eaf0f8 Mon Sep 17 00:00:00 2001 From: "J. Xavier Prochaska" Date: Tue, 2 Dec 2014 08:18:56 -0800 Subject: [PATCH 06/18] jxp: moved boss stuff to a new module --- py/qso_template/fit_boss_qsos.py | 189 +++++++++++++++++++++++++++++++ 1 file changed, 189 insertions(+) create mode 100644 py/qso_template/fit_boss_qsos.py diff --git a/py/qso_template/fit_boss_qsos.py b/py/qso_template/fit_boss_qsos.py new file mode 100644 index 000000000..42f671763 --- /dev/null +++ b/py/qso_template/fit_boss_qsos.py @@ -0,0 +1,189 @@ +""" +#;+ +#; NAME: +#; qso_pca +#; Version 1.0 +#; +#; PURPOSE: +#; Module for generate QSO PCA templates +#; 24-Nov-2014 by JXP +#;- +#;------------------------------------------------------------------------------ +""" +from __future__ import print_function, absolute_import, division, unicode_literals + +import numpy as np +import os +import multiprocessing as mp + +from astropy.io import fits +flg_xdb = True +try: + from xastropy.xutils import xdebug as xdb +except ImportError: + flg_xdb = False + +# +def read_qso_eigen(eigen_fil=None): + ''' + Input the QSO Eigenspectra + ''' + # File + if eigen_fil is None: + eigen_fil = os.environ.get('IDLSPEC2D_DIR')+'/templates/spEigenQSO-55732.fits' + print('Using these eigen spectra: {:s}'.format(eigen_fil)) + hdu = fits.open(eigen_fil) + eigen = hdu[0].data + head = hdu[0].header + # Rest wavelength + eigen_wave = 10.**(head['COEFF0'] + np.arange(head['NAXIS1'])*head['COEFF1']) + + # Return + return eigen, eigen_wave + +## +def fit_eigen(flux,ivar,eigen_flux): + ''' + Fit the spectrum with the eigenvectors. + Pass back the coefficients + ''' + #C = np.diag(1./ivar) + Cinv = np.diag(ivar) + A = eigen_flux.T + + #alpha = np.dot(A.T, np.linalg.solve(C, A)) # Numerical Recipe notation + alpha = np.dot(A.T, np.dot(Cinv,A)) + cov = np.linalg.inv(alpha) + #beta = np.dot(A.T, np.linalg.solve(C, y)) + beta = np.dot(A.T, np.dot(Cinv, flux)) + acoeff = np.dot(cov, beta) + + # Return + return acoeff + +## +def do_boss_lya_parallel(istart, iend, output, debug=False, cut_Lya=True): + ''' + Generate PCA coeff for the BOSS Lya DR10 dataset, v2.1 + + Parameters + ---------- + cut_Lya: boolean (True) + Avoid using the Lya forest in the analysis + ''' + # Eigen + eigen, eigen_wave = read_qso_eigen() + + # Open the BOSS catalog file + boss_cat_fil = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_cat_v2.1.fits.gz' + bcat_hdu = fits.open(boss_cat_fil) + t_boss = bcat_hdu[1].data + nqso = len(t_boss) + + pca_val = np.zeros((iend-istart, 4)) + + # Loop us -- Should spawn on multiple CPU + #for ii in range(nqso): + datdir = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_spectra_v2.1/' + jj = 0 + for ii in range(istart,iend): + if (ii % 100) == 0: + print('ii = {:d}'.format(ii)) + # Spectrum file + pnm = str(t_boss['PLATE'][ii]) + fnm = str(t_boss['FIBERID'][ii]).rjust(4,str('0')) + mjd = str(t_boss['MJD'][ii]) + sfil = datdir+pnm+'/speclya-' + sfil = sfil+pnm+'-'+mjd+'-'+fnm+'.fits.gz' + # Read spectrum + spec_hdu = fits.open(sfil) + t = spec_hdu[1].data + flux = t['flux'] + wave = 10.**t['loglam'] + ivar = t['ivar'] + zqso = t_boss['z_pipe'][ii] + + wrest = wave / (1+zqso) + wlya = 1215. + + # Cut Lya forest? + if cut_Lya is True: + Ly_imn = np.argmin(np.abs(wrest-wlya)) + else: + Ly_imn = 0 + + # Pack + imn = np.argmin(np.abs(wrest[Ly_imn]-eigen_wave)) + npix = len(wrest[Ly_imn:]) + imx = npix+imn + eigen_flux = eigen[:,imn:imx] + + # FIT + acoeff = fit_eigen(flux[Ly_imn:], ivar[Ly_imn:], eigen_flux) + pca_val[jj,:] = acoeff + jj += 1 + + # Check + if debug is True: + model = np.dot(eigen.T,acoeff) + if flg_xdb is True: + xdb.xplot(wrest, flux, xtwo=eigen_wave, ytwo=model) + xdb.set_trace() + + #xdb.set_trace() + if output is not None: + output.put((istart,iend,pca_val)) + +## ################################# +## ################################# +## TESTING +## ################################# +if __name__ == '__main__': + + # Run + #do_boss_lya_parallel(0,10,None,debug=True,cut_Lya=True) + #xdb.set_trace() + + ## ############################ + # Parallel + boss_cat_fil = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_cat_v2.1.fits.gz' + bcat_hdu = fits.open(boss_cat_fil) + t_boss = bcat_hdu[1].data + nqso = len(t_boss) + nqso = 45 # Testing + + output = mp.Queue() + processes = [] + nproc = 4 + nsub = nqso // nproc + # Setup the Processes + for ii in range(nproc): + # Generate + istrt = ii * nsub + if ii == (nproc-1): + iend = nqso + else: + iend = (ii+1)*nsub + #xdb.set_trace() + process = mp.Process(target=do_boss_lya_parallel, + args=(istrt,iend,output)) + processes.append(process) + + # Run processes + for p in processes: + p.start() + + # Exit the completed processes + for p in processes: + p.join() + + # Get process results from the output queue + results = [output.get() for p in processes] + + # Bring together + #sorted(results, key=lambda result: result[0]) + #all_is = [ir[0] for ir in results] + pca_val = np.zeros((nqso, 4)) + for ir in results: + pca_val[ir[0]:ir[1],:] = ir[2] + xdb.set_trace() From 6a2add5daedf3ac96292278090f8788fce1eed14 Mon Sep 17 00:00:00 2001 From: "J. Xavier Prochaska" Date: Tue, 2 Dec 2014 15:44:00 -0800 Subject: [PATCH 07/18] jxp: boss qsos --- py/qso_template/fit_boss_qsos.py | 40 ++++++++++++++++++++++++++++---- 1 file changed, 35 insertions(+), 5 deletions(-) diff --git a/py/qso_template/fit_boss_qsos.py b/py/qso_template/fit_boss_qsos.py index 42f671763..b42394ed8 100644 --- a/py/qso_template/fit_boss_qsos.py +++ b/py/qso_template/fit_boss_qsos.py @@ -15,8 +15,12 @@ import numpy as np import os import multiprocessing as mp +import Queue + +#Queue.Queue(30000000) from astropy.io import fits + flg_xdb = True try: from xastropy.xutils import xdebug as xdb @@ -87,7 +91,7 @@ def do_boss_lya_parallel(istart, iend, output, debug=False, cut_Lya=True): datdir = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_spectra_v2.1/' jj = 0 for ii in range(istart,iend): - if (ii % 100) == 0: + if (ii % 1000) == 0: print('ii = {:d}'.format(ii)) # Spectrum file pnm = str(t_boss['PLATE'][ii]) @@ -131,8 +135,10 @@ def do_boss_lya_parallel(istart, iend, output, debug=False, cut_Lya=True): xdb.set_trace() #xdb.set_trace() + print('Done with my subset {:d}, {:d}'.format(istart,iend)) if output is not None: output.put((istart,iend,pca_val)) + #output.put(None) ## ################################# ## ################################# @@ -150,12 +156,13 @@ def do_boss_lya_parallel(istart, iend, output, debug=False, cut_Lya=True): bcat_hdu = fits.open(boss_cat_fil) t_boss = bcat_hdu[1].data nqso = len(t_boss) - nqso = 45 # Testing + #nqso = 4500 # Testing output = mp.Queue() processes = [] - nproc = 4 + nproc = 8 nsub = nqso // nproc + # Setup the Processes for ii in range(nproc): # Generate @@ -173,12 +180,16 @@ def do_boss_lya_parallel(istart, iend, output, debug=False, cut_Lya=True): for p in processes: p.start() + # Get process results from the output queue + print('Grabbing Output') + results = [output.get() for p in processes] + #xdb.set_trace() + # Exit the completed processes + print('Joining') for p in processes: p.join() - # Get process results from the output queue - results = [output.get() for p in processes] # Bring together #sorted(results, key=lambda result: result[0]) @@ -186,4 +197,23 @@ def do_boss_lya_parallel(istart, iend, output, debug=False, cut_Lya=True): pca_val = np.zeros((nqso, 4)) for ir in results: pca_val[ir[0]:ir[1],:] = ir[2] + + # Write to disk as a binary FITS table + col0 = fits.Column(name='PCA0',format='E',array=pca_val[:,0]) + col1 = fits.Column(name='PCA1',format='E',array=pca_val[:,1]) + col2 = fits.Column(name='PCA2',format='E',array=pca_val[:,2]) + col3 = fits.Column(name='PCA3',format='E',array=pca_val[:,3]) + cols = fits.ColDefs([col0, col1, col2, col3]) + tbhdu = fits.BinTableHDU.from_columns(cols) + + prihdr = fits.Header() + prihdr['OBSERVER'] = 'Edwin Hubble' + prihdr['COMMENT'] = "Here's some commentary about this FITS file." + prihdu = fits.PrimaryHDU(header=prihdr) + + thdulist = fits.HDUList([prihdu, tbhdu]) + thdulist.writeto('BOSS_DR10Lya_PCA_values.fits',clobber=True) + + # Done xdb.set_trace() + print('All done') From 49dd4a32246df00fd25ac1c9e72715adabc8c1dd Mon Sep 17 00:00:00 2001 From: "J. Xavier Prochaska" Date: Tue, 2 Dec 2014 17:08:34 -0800 Subject: [PATCH 08/18] more boss --- py/qso_template/boss_qsos_figs.py | 0 py/qso_template/fit_boss_qsos.py | 6 +++--- 2 files changed, 3 insertions(+), 3 deletions(-) create mode 100644 py/qso_template/boss_qsos_figs.py diff --git a/py/qso_template/boss_qsos_figs.py b/py/qso_template/boss_qsos_figs.py new file mode 100644 index 000000000..e69de29bb diff --git a/py/qso_template/fit_boss_qsos.py b/py/qso_template/fit_boss_qsos.py index b42394ed8..1ab8b724d 100644 --- a/py/qso_template/fit_boss_qsos.py +++ b/py/qso_template/fit_boss_qsos.py @@ -1,12 +1,12 @@ """ #;+ #; NAME: -#; qso_pca +#; fit_boss_qsos #; Version 1.0 #; #; PURPOSE: -#; Module for generate QSO PCA templates -#; 24-Nov-2014 by JXP +#; Module for Fitting PCA to the BOSS QSOs +#; 01-Dec-2014 by JXP #;- #;------------------------------------------------------------------------------ """ From 9bcc2464120ea54901d02f96ba100c45baa7985a Mon Sep 17 00:00:00 2001 From: "J. Xavier Prochaska" Date: Mon, 8 Dec 2014 05:16:56 -0800 Subject: [PATCH 09/18] boss qso figs by JXP --- py/qso_template/boss_qsos_figs.py | 234 ++++++++++++++++++++++++++++++ 1 file changed, 234 insertions(+) diff --git a/py/qso_template/boss_qsos_figs.py b/py/qso_template/boss_qsos_figs.py index e69de29bb..9bf5efa3e 100644 --- a/py/qso_template/boss_qsos_figs.py +++ b/py/qso_template/boss_qsos_figs.py @@ -0,0 +1,234 @@ +""" +#;+ +#; NAME: +#; boss_qso_figs +#; Version 1.0 +#; +#; PURPOSE: +#; Module for making figures with the BOSS QSO fit outpus +#; 02-Dec-2014 by JXP +#;- +#;------------------------------------------------------------------------------ +""" +from __future__ import print_function, absolute_import, division, unicode_literals + +import numpy as np +import os + +import matplotlib as mpl +mpl.rcParams['font.family'] = 'stixgeneral' +from matplotlib.backends.backend_pdf import PdfPages +from matplotlib import pyplot as plt +import matplotlib.gridspec as gridspec +from matplotlib.colors import LogNorm + +from astropy.io import fits + +from xastropy.plotting import utils as xputils + +try: + from xastropy.xutils import xdebug as xdb +except ImportError: + import pdb as xdb + +# ##################### ##################### +# ##################### ##################### +# Plots the PCA coeff against one another from BOSS DR10 QSOs +def fig_boss_pca_coeff(outfil=None, boss_fil=None): + + # Todo + # Include NHI on the label + # Imports + + # Read FITS table + if boss_fil is None: + boss_fil = 'BOSS_DR10Lya_PCA_values.fits.gz' + hdu = fits.open(boss_fil) + pca_coeff = hdu[1].data + #xdb.set_trace() + + # Initialize + #if 'xmnx' not in locals(): + # xmnx = (17.0, 20.4) + #ymnx = ((-1.9, 2.3), + # (-1.9, 1.4), + # (-1.1, 1.4), + # (-2.1, 0.9)) + + ms = 1. # point size + scl = 100. + + allxi = [0,0,0,1,1,2] + allyi = [1,2,3,2,3,3] + + rngi = ([-0.02, 0.3], + [-0.2,0.4], + [-0.2, 0.5], + [-0.2, 0.4]) + + # Start the plot + if outfil != None: + pp = PdfPages(outfil) + + plt.figure(figsize=(4, 5.5)) + plt.clf() + gs = gridspec.GridSpec(3, 2) + + + # Looping + for ii in range(6): + + # Axis + ax = plt.subplot(gs[ii]) + #ax = plt.subplot(gs[ii//2,ii%2]) + + ''' + ax.xaxis.set_minor_locator(plt.MultipleLocator(0.5)) + ax.xaxis.set_major_locator(plt.MultipleLocator(1.)) + #ax.yaxis.set_minor_locator(plt.MultipleLocator(0.5)) + ax.yaxis.set_major_locator(plt.MultipleLocator(0.5)) + ax.set_xlim(xmnx) + ax.set_ylim((-0.5, 0.5)) + ''' + + xi = allxi[ii] + xlbl=str('PCA'+str(xi)) + yi = allyi[ii] + ylbl=str('PCA'+str(yi)) + + + # Labels + ax.set_xlabel(xlbl) + ax.set_ylabel(ylbl) + + # Data + + # Values + #ax.scatter(pca_coeff[xlbl], pca_coeff[ylbl], color='black', s=ms)#, marker=mark) + ax.hist2d(scl*pca_coeff[xlbl], scl*pca_coeff[ylbl], bins=100, norm=LogNorm(), + range=[rngi[xi],rngi[yi]]) + + # Font size + xputils.set_fontsize(ax,8.) + + # Layout and save + plt.tight_layout(pad=0.2,h_pad=0.1,w_pad=0.25) + if outfil != None: + pp.savefig(bbox_inches='tight') + pp.close() + else: + plt.show() + +# ##################### ##################### +# ##################### ##################### +# Plots X vs PCA coeff from BOSS DR10 QSOs +def fig_boss_x_vs_pca(outfil=None, boss_fil=None, flg=0): + ''' + flg = 0: Redshift + flg = 1: imag + ''' + + # Todo + # Include NHI on the label + # Imports + + # Read PCA FITS table + if boss_fil is None: + boss_fil = 'BOSS_DR10Lya_PCA_values.fits.gz' + hdu = fits.open(boss_fil) + pca_coeff = hdu[1].data + + # Read BOSS catalog table + boss_cat_fil = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_cat_v2.1.fits.gz' + bcat_hdu = fits.open(boss_cat_fil) + t_boss = bcat_hdu[1].data + #xdb.set_trace() + if flg == 0: + xparm = t_boss['z_pipe'] + xmnx = [2., 4.] + xlbl=str(r'$z_{\rm QSO}$') + elif flg == 1: + tmp = t_boss['PSFMAG'] + xparm = tmp[:,3] # i-band mag + xlbl=str('i mag') + xmnx = [17.,22.5] + else: + raise ValueError('fig_boss_x_vs_pca: flg={:d} not allowed'.format(flg)) + + # Initialize + #if 'xmnx' not in locals(): + # xmnx = (17.0, 20.4) + ymnx = ((-0.02, 0.3), + (-0.2,0.4), + (-0.2, 0.5), + (-0.2, 0.4)) + + ms = 1. # point size + scl = 100. + + + # Start the plot + if outfil != None: + pp = PdfPages(outfil) + + plt.figure(figsize=(4, 4)) + plt.clf() + gs = gridspec.GridSpec(2, 2) + + + # Looping + for ii in range(4): + + # Axis + ax = plt.subplot(gs[ii]) + #ax = plt.subplot(gs[ii//2,ii%2]) + + + ylbl=str('PCA'+str(ii)) + + # Labels + ax.set_xlabel(xlbl) + ax.set_ylabel(ylbl) + + # Data + + # Values + #xdb.set_trace() + ax.hist2d(xparm, scl*pca_coeff[ylbl], bins=100, norm=LogNorm(), + range=[xmnx, [ymnx[ii][0], ymnx[ii][1]]]) + + # Font size + xputils.set_fontsize(ax,8.) + + # Layout and save + plt.tight_layout(pad=0.2,h_pad=0.1,w_pad=0.25) + if outfil != None: + pp.savefig(bbox_inches='tight') + pp.close() + else: + plt.show() + + +#### ########################## ######################### +#### ########################## ######################### +#### ########################## ######################### + +# Command line execution +if __name__ == '__main__': + + flg_fig = 0 + flg_fig += 1 # PCA vs PCA + flg_fig += 2**1 # PCA vs z + flg_fig += 2**2 # PCA vs i + + # PCA vs PCA + if (flg_fig % 2) == 1: + fig_boss_pca_coeff(outfil='fig_boss_pca.pdf') + + # PCA vs z + if (flg_fig % 2**2) >= 2**1: + fig_boss_x_vs_pca(outfil='fig_boss_z_vs_pca.pdf',flg=0) + + # PCA vs imag + if (flg_fig % 2**3) >= 2**2: + fig_boss_x_vs_pca(outfil='fig_boss_imag_vs_pca.pdf',flg=1) From 545937bae43c905c92a0ac8bc430d30e27e43f6e Mon Sep 17 00:00:00 2001 From: "J. Xavier Prochaska" Date: Mon, 8 Dec 2014 05:17:40 -0800 Subject: [PATCH 10/18] all files --- py/qso_template/__init__.py | 1 + py/qso_template/fit_boss_qsos.py | 17 +++++++++++++---- py/qso_template/tests.py | 7 +++++-- 3 files changed, 19 insertions(+), 6 deletions(-) diff --git a/py/qso_template/__init__.py b/py/qso_template/__init__.py index 7f366e162..9ff0a1d66 100644 --- a/py/qso_template/__init__.py +++ b/py/qso_template/__init__.py @@ -1 +1,2 @@ import qso_pca +import tests diff --git a/py/qso_template/fit_boss_qsos.py b/py/qso_template/fit_boss_qsos.py index 1ab8b724d..edea366bf 100644 --- a/py/qso_template/fit_boss_qsos.py +++ b/py/qso_template/fit_boss_qsos.py @@ -66,7 +66,7 @@ def fit_eigen(flux,ivar,eigen_flux): return acoeff ## -def do_boss_lya_parallel(istart, iend, output, debug=False, cut_Lya=True): +def do_boss_lya_parallel(istart, iend, cut_Lya, output, debug=False): ''' Generate PCA coeff for the BOSS Lya DR10 dataset, v2.1 @@ -86,6 +86,9 @@ def do_boss_lya_parallel(istart, iend, output, debug=False, cut_Lya=True): pca_val = np.zeros((iend-istart, 4)) + if cut_Lya is False: + print('do_boss: Not cutting the Lya Forest in the fit') + # Loop us -- Should spawn on multiple CPU #for ii in range(nqso): datdir = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_spectra_v2.1/' @@ -163,6 +166,8 @@ def do_boss_lya_parallel(istart, iend, output, debug=False, cut_Lya=True): nproc = 8 nsub = nqso // nproc + cut_Lya = False + # Setup the Processes for ii in range(nproc): # Generate @@ -173,7 +178,7 @@ def do_boss_lya_parallel(istart, iend, output, debug=False, cut_Lya=True): iend = (ii+1)*nsub #xdb.set_trace() process = mp.Process(target=do_boss_lya_parallel, - args=(istrt,iend,output)) + args=(istrt,iend,cut_Lya, output)) processes.append(process) # Run processes @@ -212,8 +217,12 @@ def do_boss_lya_parallel(istart, iend, output, debug=False, cut_Lya=True): prihdu = fits.PrimaryHDU(header=prihdr) thdulist = fits.HDUList([prihdu, tbhdu]) - thdulist.writeto('BOSS_DR10Lya_PCA_values.fits',clobber=True) + if cut_Lya is False: + outfil = 'BOSS_DR10Lya_PCA_values_nocut.fits' + else: + outfil = 'BOSS_DR10Lya_PCA_values.fits' + thdulist.writeto(outfil, clobber=True) # Done - xdb.set_trace() + #xdb.set_trace() print('All done') diff --git a/py/qso_template/tests.py b/py/qso_template/tests.py index e2e9a08ad..f2115856d 100644 --- a/py/qso_template/tests.py +++ b/py/qso_template/tests.py @@ -120,12 +120,15 @@ """ C = np.diag(1./ivar) +Cinv = np.diag(ivar) A = eigen_flux.T y = flux -alpha = np.dot(A.T, np.linalg.solve(C, A)) # Numerical Recipe notation +#alpha = np.dot(A.T, np.linalg.solve(C, A)) # Numerical Recipe notation +alpha = np.dot(A.T, np.dot(Cinv,A)) cov = np.linalg.inv(alpha) -beta = np.dot(A.T, np.linalg.solve(C, y)) +#beta = np.dot(A.T, np.linalg.solve(C, y)) +beta = np.dot(A.T, np.dot(Cinv, y)) acoeff = np.dot(cov, beta) print('acoeff = {:g}, {:g}, {:g}, {:g}'.format(*acoeff)) From f9f0d44c0cac5e8e28cb61888151a2aff62200cf Mon Sep 17 00:00:00 2001 From: "J. Xavier Prochaska" Date: Mon, 8 Dec 2014 06:46:24 -0800 Subject: [PATCH 11/18] some templates now; JXP --- py/qso_template/boss_qsos_figs.py | 209 ++++++++++++++++++++++++++++-- py/qso_template/desi_qso_templ.py | 195 ++++++++++++++++++++++++++++ py/qso_template/qso_pca.py | 42 ++++-- 3 files changed, 427 insertions(+), 19 deletions(-) create mode 100644 py/qso_template/desi_qso_templ.py diff --git a/py/qso_template/boss_qsos_figs.py b/py/qso_template/boss_qsos_figs.py index 9bf5efa3e..bb59b032c 100644 --- a/py/qso_template/boss_qsos_figs.py +++ b/py/qso_template/boss_qsos_figs.py @@ -26,6 +26,8 @@ from xastropy.plotting import utils as xputils +import fit_boss_qsos as fbq + try: from xastropy.xutils import xdebug as xdb except ImportError: @@ -42,7 +44,7 @@ def fig_boss_pca_coeff(outfil=None, boss_fil=None): # Read FITS table if boss_fil is None: - boss_fil = 'BOSS_DR10Lya_PCA_values.fits.gz' + boss_fil = 'BOSS_DR10Lya_PCA_values_nocut.fits.gz' hdu = fits.open(boss_fil) pca_coeff = hdu[1].data #xdb.set_trace() @@ -70,9 +72,9 @@ def fig_boss_pca_coeff(outfil=None, boss_fil=None): if outfil != None: pp = PdfPages(outfil) - plt.figure(figsize=(4, 5.5)) + plt.figure(figsize=(5.5, 4)) plt.clf() - gs = gridspec.GridSpec(3, 2) + gs = gridspec.GridSpec(2, 3) # Looping @@ -134,7 +136,7 @@ def fig_boss_x_vs_pca(outfil=None, boss_fil=None, flg=0): # Read PCA FITS table if boss_fil is None: - boss_fil = 'BOSS_DR10Lya_PCA_values.fits.gz' + boss_fil = 'BOSS_DR10Lya_PCA_values_nocut.fits.gz' hdu = fits.open(boss_fil) pca_coeff = hdu[1].data @@ -175,7 +177,6 @@ def fig_boss_x_vs_pca(outfil=None, boss_fil=None, flg=0): plt.clf() gs = gridspec.GridSpec(2, 2) - # Looping for ii in range(4): @@ -209,6 +210,188 @@ def fig_boss_x_vs_pca(outfil=None, boss_fil=None, flg=0): plt.show() +# ##################### ##################### +# ##################### ##################### +# Plots BOSS PCZ Eigenvectors +def fig_boss_eigen(outfil=None, boss_fil=None, flg=0): + ''' + flg = 0: Redshift + flg = 1: imag + ''' + + # Todo + # Include NHI on the label + # Imports + + # Read + eigen, eigen_wave = fbq.read_qso_eigen() + + # Initialize + ymnx = ((0.0, 10), + (-1, 1), + (-1, 1), + (-1, 1)) + + # Start the plot + if outfil != None: + pp = PdfPages(outfil) + + plt.figure(figsize=(8, 5)) + plt.clf() + gs = gridspec.GridSpec(4, 1) + #xdb.set_trace() + + # Looping + for ii in range(4): + + # Axis + ax = plt.subplot(gs[ii]) + #ax = plt.subplot(gs[ii//2,ii%2]) + + ylbl=str('Eigen'+str(ii)) + + # Labels + if ii == 3: + ax.set_xlabel('Rest Wavelength') + else: + ax.get_xaxis().set_ticks([]) + ax.set_ylabel(ylbl) + + ax.set_xlim((400., 8000)) + + # Data + + # Values + ax.plot( eigen_wave, eigen[ii,:], 'k-',drawstyle='steps-mid', linewidth=0.5) + + # Layout and save + plt.tight_layout(pad=0.2,h_pad=0.0,w_pad=0.25) + if outfil != None: + pp.savefig(bbox_inches='tight') + pp.close() + else: + plt.show() + +# ##################### ##################### +# ##################### ##################### +# Plots i vs zQSO +def fig_boss_i_vs_z(outfil=None, boss_fil=None, flg=0): + ''' + ''' + + # Todo + # Include NHI on the label + # Imports + + # Read BOSS catalog table + boss_cat_fil = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_cat_v2.1.fits.gz' + bcat_hdu = fits.open(boss_cat_fil) + t_boss = bcat_hdu[1].data + #xdb.set_trace() + zQSO = t_boss['z_pipe'] + xmnx = [2., 4.] + xlbl=str(r'$z_{\rm QSO}$') + tmp = t_boss['PSFMAG'] + imag = tmp[:,3] # i-band mag + ylbl=str('i mag') + ymnx = [17.,22.5] + + # Initialize + + ms = 1. # point size + scl = 100. + + + # Start the plot + if outfil != None: + pp = PdfPages(outfil) + + plt.figure(figsize=(4, 4)) + plt.clf() + gs = gridspec.GridSpec(1, 1) + + # Looping + ax = plt.subplot(gs[0]) + + # Labels + ax.set_xlabel(xlbl) + ax.set_ylabel(ylbl) + + #xdb.set_trace() + ax.hist2d(zQSO, imag, bins=100, norm=LogNorm(), range=[xmnx, ymnx]) + + # Font size + #xputils.set_fontsize(ax,8.) + + # Layout and save + plt.tight_layout(pad=0.2,h_pad=0.1,w_pad=0.25) + if outfil != None: + pp.savefig(bbox_inches='tight') + pp.close() + else: + plt.show() + +# ##################### ##################### +# ##################### ##################### +# Plots DESI templates at a range of z and imag +def fig_desi_templ_z_i(outfil=None, boss_fil=None, flg=0): + ''' + flg = 0: Redshift + flg = 1: imag + ''' + + # Todo + # Include NHI on the label + # Imports + + + # Initialize + ymnx = ((0.0, 10), + (-1, 1), + (-1, 1), + (-1, 1)) + + # Start the plot + if outfil != None: + pp = PdfPages(outfil) + + plt.figure(figsize=(8, 5)) + plt.clf() + gs = gridspec.GridSpec(4, 1) + #xdb.set_trace() + + # Looping + for ii in range(4): + + # Axis + ax = plt.subplot(gs[ii]) + #ax = plt.subplot(gs[ii//2,ii%2]) + + ylbl=str('Eigen'+str(ii)) + + # Labels + if ii == 3: + ax.set_xlabel('Rest Wavelength') + else: + ax.get_xaxis().set_ticks([]) + ax.set_ylabel(ylbl) + + ax.set_xlim((400., 8000)) + + # Data + + # Values + ax.plot( eigen_wave, eigen[ii,:], 'k-',drawstyle='steps-mid', linewidth=0.5) + + # Layout and save + plt.tight_layout(pad=0.2,h_pad=0.0,w_pad=0.25) + if outfil != None: + pp.savefig(bbox_inches='tight') + pp.close() + else: + plt.show() + + #### ########################## ######################### #### ########################## ######################### #### ########################## ######################### @@ -217,9 +400,11 @@ def fig_boss_x_vs_pca(outfil=None, boss_fil=None, flg=0): if __name__ == '__main__': flg_fig = 0 - flg_fig += 1 # PCA vs PCA - flg_fig += 2**1 # PCA vs z - flg_fig += 2**2 # PCA vs i + #flg_fig += 1 # PCA vs PCA + #flg_fig += 2**1 # PCA vs z + #flg_fig += 2**2 # PCA vs i + #flg_fig += 2**3 # Eigenvectors + flg_fig += 2**4 # z vs i for BOSS # PCA vs PCA if (flg_fig % 2) == 1: @@ -232,3 +417,11 @@ def fig_boss_x_vs_pca(outfil=None, boss_fil=None, flg=0): # PCA vs imag if (flg_fig % 2**3) >= 2**2: fig_boss_x_vs_pca(outfil='fig_boss_imag_vs_pca.pdf',flg=1) + + # Eigenvectors + if (flg_fig % 2**4) >= 2**3: + fig_boss_eigen(outfil='fig_boss_eigen.pdf') + + # z vs i + if (flg_fig % 2**5) >= 2**4: + fig_boss_i_vs_z(outfil='fig_boss_i_vs_z.pdf') diff --git a/py/qso_template/desi_qso_templ.py b/py/qso_template/desi_qso_templ.py new file mode 100644 index 000000000..49ae11903 --- /dev/null +++ b/py/qso_template/desi_qso_templ.py @@ -0,0 +1,195 @@ +""" +#;+ +#; NAME: +#; fit_boss_qsos +#; Version 1.0 +#; +#; PURPOSE: +#; Module for Fitting PCA to the BOSS QSOs +#; 01-Dec-2014 by JXP +#;- +#;------------------------------------------------------------------------------ +""" +from __future__ import print_function, absolute_import, division, unicode_literals + +import numpy as np +import os + +from astropy.io import fits + +import fit_boss_qsos as fbq + +flg_xdb = True +try: + from xastropy.xutils import xdebug as xdb +except ImportError: + flg_xdb = False + + +## +def mean_templ_zi(zimag, debug=False, i_wind=0.1, z_wind=0.05, + boss_pca_fil=None): + ''' + Generate 'mean' templates at given z,i + + Parameters + ---------- + zimag: list of tuples + Redshift, imag pairs for the templates + i_wind: float (0.1 mag) + Window for smoothing imag + z_wind: float (0.05 mag) + Window for smoothing redshift + ''' + # PCA values + if boss_pca_fil is None: + boss_pca_fil = 'BOSS_DR10Lya_PCA_values_nocut.fits.gz' + hdu = fits.open(boss_pca_fil) + pca_coeff = hdu[1].data + + # BOSS Eigenvectors + eigen, eigen_wave = fbq.read_qso_eigen() + npix = len(eigen_wave) + + # Open the BOSS catalog file + boss_cat_fil = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_cat_v2.1.fits.gz' + bcat_hdu = fits.open(boss_cat_fil) + t_boss = bcat_hdu[1].data + nqso = len(t_boss) + zQSO = t_boss['z_pipe'] + tmp = t_boss['PSFMAG'] + imag = tmp[:,3] # i-band mag + + # Output array + ntempl = len(zimag) + out_spec = np.zeros( (ntempl, npix) ) + + # Iterate on z,imag + for izi in zimag: + tt = zimag.index(izi) + # Find matches + idx = np.where( (np.fabs(imag-izi[1]) < i_wind) & + (np.fabs(zQSO-izi[0]) < z_wind))[0] + if len(idx) < 50: + raise ValueError('mean_templ_zi: Not enough QSOs! {:d}'.format(len(idx))) + + # Calculate median PCA values + PCA0 = np.median(pca_coeff['PCA0'][idx]) + PCA1 = np.median(pca_coeff['PCA1'][idx]) + PCA2 = np.median(pca_coeff['PCA2'][idx]) + PCA3 = np.median(pca_coeff['PCA3'][idx]) + acoeff = np.array( [PCA0, PCA1, PCA2, PCA3] ) + + # Make the template + out_spec[tt,:] = np.dot(eigen.T,acoeff) + if debug is True: + xdb.xplot(eigen_wave*(1.+izi[0]), out_spec[tt,:]) + xdb.set_trace() + + # Return + return out_spec + +# ##################### ##################### +# ##################### ##################### +# Plots DESI templates at a range of z and imag +def fig_desi_templ_z_i(outfil=None, boss_fil=None, flg=0): + ''' + flg = 0: Redshift + flg = 1: imag + ''' + + # Todo + # Include NHI on the label + # Imports + import matplotlib as mpl + mpl.rcParams['font.family'] = 'stixgeneral' + from matplotlib.backends.backend_pdf import PdfPages + from matplotlib import pyplot as plt + import matplotlib.gridspec as gridspec + from matplotlib.colors import LogNorm + + # Eigen (for wavelengths) + eigen, eigen_wave = fbq.read_qso_eigen() + + all_zi = [ [ (2.3, 18.5), (2.3, 19.5), (2.3, 20.5), (2.3, 21.5) ], + [ (2.5, 18.5), (2.5, 19.5), (2.5, 20.5), (2.5, 21.5) ], + [ (2.7, 19.5), (2.7, 20.5), (2.7, 21.5) ], + [ (3.2, 19.5), (3.2, 20.5), (3.2, 21.5) ] ] + + xmnx = (3600., 9000.) + + # Start the plot + if outfil != None: + pp = PdfPages(outfil) + + plt.figure(figsize=(8, 5)) + plt.clf() + gs = gridspec.GridSpec(4, 1) + #xdb.set_trace() + + # Looping + for ii in range(4): + + # Get the templates + ztempl = all_zi[ii][0][0] + spec = mean_templ_zi(all_zi[ii]) + + # Axis + ax = plt.subplot(gs[ii]) + #ax = plt.subplot(gs[ii//2,ii%2]) + + # Labels + if ii == 3: + ax.set_xlabel('Wavelength') + else: + ax.get_xaxis().set_ticks([]) + ax.set_ylabel('Flux') + + ax.set_xlim(xmnx) + + # Data + + # Values + for jj in range(len(all_zi[ii])): + ax.plot( eigen_wave*(1.+ ztempl), + spec[jj,:], '-',drawstyle='steps-mid', linewidth=0.5) + if jj == 0: + ymx = 1.05*np.max(spec[jj,:]) + ax.set_ylim((0., ymx)) + # Label + zlbl = 'z={:g}'.format(ztempl) + ax.text(7000., ymx*0.7, zlbl) + + # Layout and save + plt.tight_layout(pad=0.2,h_pad=0.0,w_pad=0.25) + if outfil != None: + pp.savefig(bbox_inches='tight') + pp.close() + else: + plt.show() + + +## ################################# +## ################################# +## TESTING +## ################################# +if __name__ == '__main__': + + # Run + flg_test = 0 + #flg_test += 1 # Mean templates with z,imag + flg_test += 2 # Mean templates with z,imag + + # Make Mean templates + if (flg_test % 2) == 1: + zimag = [ (2.3, 19.) ] + mean_templ_zi(zimag) + + # Mean template fig + if (flg_test % 2**2) >= 2**1: + fig_desi_templ_z_i(outfil='fig_desi_templ_z_i.pdf') + + + # Done + #xdb.set_trace() + print('All done') diff --git a/py/qso_template/qso_pca.py b/py/qso_template/qso_pca.py index b1cc1d2a3..42f671763 100644 --- a/py/qso_template/qso_pca.py +++ b/py/qso_template/qso_pca.py @@ -62,9 +62,14 @@ def fit_eigen(flux,ivar,eigen_flux): return acoeff ## -def do_boss_lya_parallel(istart, iend, output, debug=False): +def do_boss_lya_parallel(istart, iend, output, debug=False, cut_Lya=True): ''' Generate PCA coeff for the BOSS Lya DR10 dataset, v2.1 + + Parameters + ---------- + cut_Lya: boolean (True) + Avoid using the Lya forest in the analysis ''' # Eigen eigen, eigen_wave = read_qso_eigen() @@ -99,27 +104,35 @@ def do_boss_lya_parallel(istart, iend, output, debug=False): zqso = t_boss['z_pipe'][ii] wrest = wave / (1+zqso) - npix = len(wrest) + wlya = 1215. + # Cut Lya forest? + if cut_Lya is True: + Ly_imn = np.argmin(np.abs(wrest-wlya)) + else: + Ly_imn = 0 + # Pack - imn = np.argmin(np.abs(wrest[0]-eigen_wave)) + imn = np.argmin(np.abs(wrest[Ly_imn]-eigen_wave)) + npix = len(wrest[Ly_imn:]) imx = npix+imn eigen_flux = eigen[:,imn:imx] # FIT - acoeff = fit_eigen(flux, ivar, eigen_flux) + acoeff = fit_eigen(flux[Ly_imn:], ivar[Ly_imn:], eigen_flux) pca_val[jj,:] = acoeff jj += 1 # Check if debug is True: - model = np.dot(eigen_flux.T,acoeff) + model = np.dot(eigen.T,acoeff) if flg_xdb is True: - xdb.xplot(wrest, flux, model) + xdb.xplot(wrest, flux, xtwo=eigen_wave, ytwo=model) xdb.set_trace() #xdb.set_trace() - output.put((istart,iend,pca_val)) + if output is not None: + output.put((istart,iend,pca_val)) ## ################################# ## ################################# @@ -128,14 +141,16 @@ def do_boss_lya_parallel(istart, iend, output, debug=False): if __name__ == '__main__': # Run - #do_boss_lya()#debug=True) + #do_boss_lya_parallel(0,10,None,debug=True,cut_Lya=True) + #xdb.set_trace() + ## ############################ # Parallel boss_cat_fil = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_cat_v2.1.fits.gz' bcat_hdu = fits.open(boss_cat_fil) t_boss = bcat_hdu[1].data nqso = len(t_boss) - nqso = 800 # Testing + nqso = 45 # Testing output = mp.Queue() processes = [] @@ -163,7 +178,12 @@ def do_boss_lya_parallel(istart, iend, output, debug=False): p.join() # Get process results from the output queue - print('Got here') results = [output.get() for p in processes] - # Combine + + # Bring together + #sorted(results, key=lambda result: result[0]) + #all_is = [ir[0] for ir in results] + pca_val = np.zeros((nqso, 4)) + for ir in results: + pca_val[ir[0]:ir[1],:] = ir[2] xdb.set_trace() From 73ebef1e63ee602f1d15ad3352f8e99378ce87c8 Mon Sep 17 00:00:00 2001 From: "J. Xavier Prochaska" Date: Thu, 18 Dec 2014 13:51:38 -0800 Subject: [PATCH 12/18] jxp --- py/qso_template/fit_boss_qsos.py | 130 ++++++++++++++++++++++++++++--- 1 file changed, 118 insertions(+), 12 deletions(-) diff --git a/py/qso_template/fit_boss_qsos.py b/py/qso_template/fit_boss_qsos.py index edea366bf..70c31bc78 100644 --- a/py/qso_template/fit_boss_qsos.py +++ b/py/qso_template/fit_boss_qsos.py @@ -142,6 +142,93 @@ def do_boss_lya_parallel(istart, iend, cut_Lya, output, debug=False): if output is not None: output.put((istart,iend,pca_val)) #output.put(None) + +## +def do_sdss_lya_parallel(istart, iend, cut_Lya, output, debug=False): + ''' + Generate PCA coeff for the SDSS DR7 dataset, 0.50.)[0] + ivar[gd] = 1./sig[gd]**2 + zqso = t_sdss['z'][ii] + + wrest = wave / (1+zqso) + wlya = 1215. + + # Cut Lya forest? + if cut_Lya is True: + Ly_imn = np.argmin(np.abs(wrest-wlya)) + else: + Ly_imn = 0 + + # Pack + imn = np.argmin(np.abs(wrest[Ly_imn]-eigen_wave)) + npix = len(wrest[Ly_imn:]) + imx = npix+imn + eigen_flux = eigen[:,imn:imx] + + # FIT + acoeff = fit_eigen(flux[Ly_imn:], ivar[Ly_imn:], eigen_flux) + pca_val[jj,:] = acoeff + jj += 1 + + # Check + if debug is True: + model = np.dot(eigen.T,acoeff) + xdb.set_trace() + if flg_xdb is True: + xdb.xplot(wrest, flux, xtwo=eigen_wave, ytwo=model) + xdb.set_trace() + + #xdb.set_trace() + print('Done with my subset {:d}, {:d}'.format(istart,iend)) + if output is not None: + output.put((istart,iend,pca_val)) + #output.put(None) ## ################################# ## ################################# @@ -153,13 +240,26 @@ def do_boss_lya_parallel(istart, iend, cut_Lya, output, debug=False): #do_boss_lya_parallel(0,10,None,debug=True,cut_Lya=True) #xdb.set_trace() + flg = 0 # 0=BOSS, 1=SDSS + ## ############################ # Parallel - boss_cat_fil = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_cat_v2.1.fits.gz' - bcat_hdu = fits.open(boss_cat_fil) - t_boss = bcat_hdu[1].data - nqso = len(t_boss) - #nqso = 4500 # Testing + if flg == 0: + boss_cat_fil = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_cat_v2.1.fits.gz' + bcat_hdu = fits.open(boss_cat_fil) + t_boss = bcat_hdu[1].data + nqso = len(t_boss) + elif flg == 1: + sdss_cat_fil = os.environ.get('SDSSPATH')+'/DR7_QSO/dr7_qso.fits.gz' + scat_hdu = fits.open(sdss_cat_fil) + t_sdss = scat_hdu[1].data + nqso = len(t_sdss) + outfil = 'SDSS_DR7Lya_PCA_values_nocut.fits' + + nqso = 450 # Testing + + #do_sdss_lya_parallel(0, nqso, False, None) + #xdb.set_trace() output = mp.Queue() processes = [] @@ -177,18 +277,23 @@ def do_boss_lya_parallel(istart, iend, cut_Lya, output, debug=False): else: iend = (ii+1)*nsub #xdb.set_trace() - process = mp.Process(target=do_boss_lya_parallel, - args=(istrt,iend,cut_Lya, output)) + if flg == 0: + process = mp.Process(target=do_boss_lya_parallel, + args=(istrt,iend,cut_Lya, output)) + elif flg == 1: + process = mp.Process(target=do_sdss_lya_parallel, + args=(istrt,iend,cut_Lya, output)) processes.append(process) # Run processes for p in processes: p.start() + #xdb.set_trace() # Get process results from the output queue print('Grabbing Output') results = [output.get() for p in processes] - #xdb.set_trace() + xdb.set_trace() # Exit the completed processes print('Joining') @@ -217,10 +322,11 @@ def do_boss_lya_parallel(istart, iend, cut_Lya, output, debug=False): prihdu = fits.PrimaryHDU(header=prihdr) thdulist = fits.HDUList([prihdu, tbhdu]) - if cut_Lya is False: - outfil = 'BOSS_DR10Lya_PCA_values_nocut.fits' - else: - outfil = 'BOSS_DR10Lya_PCA_values.fits' + if not (outfil in locals()): + if cut_Lya is False: + outfil = 'BOSS_DR10Lya_PCA_values_nocut.fits' + else: + outfil = 'BOSS_DR10Lya_PCA_values.fits' thdulist.writeto(outfil, clobber=True) # Done From bec9ef2817021eb7463a20fef8a13051d13cc7d9 Mon Sep 17 00:00:00 2001 From: "J. Xavier Prochaska" Date: Thu, 18 Dec 2014 13:52:12 -0800 Subject: [PATCH 13/18] jxp --- py/qso_template/desi_qso_templ.py | 195 +++++++++++++++++++++++++++++- 1 file changed, 194 insertions(+), 1 deletion(-) diff --git a/py/qso_template/desi_qso_templ.py b/py/qso_template/desi_qso_templ.py index 49ae11903..90346fe98 100644 --- a/py/qso_template/desi_qso_templ.py +++ b/py/qso_template/desi_qso_templ.py @@ -19,6 +19,8 @@ import fit_boss_qsos as fbq +from xastropy.stats import basic as xstat_b + flg_xdb = True try: from xastropy.xutils import xdebug as xdb @@ -168,6 +170,187 @@ def fig_desi_templ_z_i(outfil=None, boss_fil=None, flg=0): else: plt.show() +## +def desi_qso_templates(z_wind=0.2, zmnx=(2.,4.), outfil=None, Ntempl=100, + boss_pca_fil=None): + ''' + Generate 'mean' templates at given z,i + + Parameters + ---------- + z_wind: float (0.2) + Window for sampling + zmnx: tuple ( (2,4) ) + Min/max for generation + Ntempl: int (100) + Number of draws per redshift window + ''' + # Cosmology + from astropy import cosmology + cosmo = cosmology.core.FlatLambdaCDM(70., 0.3) + # PCA values + if boss_pca_fil is None: + boss_pca_fil = 'BOSS_DR10Lya_PCA_values_nocut.fits.gz' + hdu = fits.open(boss_pca_fil) + pca_coeff = hdu[1].data + + # BOSS Eigenvectors + eigen, eigen_wave = fbq.read_qso_eigen() + npix = len(eigen_wave) + chkpix = np.where((eigen_wave > 900.) & (eigen_wave < 5000.) )[0] + lambda_912 = 911.76 + pix912 = np.argmin( np.abs(eigen_wave-lambda_912) ) + + # Open the BOSS catalog file + boss_cat_fil = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_cat_v2.1.fits.gz' + bcat_hdu = fits.open(boss_cat_fil) + t_boss = bcat_hdu[1].data + nqso = len(t_boss) + zQSO = t_boss['z_pipe'] + + + # Outfil + if outfil is None: + outfil = 'DESI_QSO_Templates_v1.1.fits' + + # Loop on redshift + z0 = np.arange(zmnx[0],zmnx[1],z_wind) + z1 = z0 + z_wind + + pca_list = ['PCA0', 'PCA1', 'PCA2', 'PCA3'] + PCA_mean = np.zeros(4) + PCA_sig = np.zeros(4) + PCA_rand = np.zeros( (4,Ntempl*2) ) + + final_spec = np.zeros( (npix, Ntempl * len(z0)) ) + final_wave = np.zeros( (npix, Ntempl * len(z0)) ) + final_z = np.zeros( Ntempl * len(z0) ) + + seed = -1422 + for ii in range(len(z0)): + # Random z values and wavelengths + zrand = np.random.uniform( z0[ii], z1[ii], Ntempl*2) + wave = np.outer(eigen_wave, 1+zrand) + + # MFP (Worseck+14) + mfp = 37. * ( (1+zrand)/5. )**(-5.4) # Physical Mpc + + # Grab PCA mean + sigma + idx = np.where( (zQSO >= z0[ii]) & (zQSO < z1[ii]) )[0] + + # Get PCA stats and random values + for ipca in pca_list: + jj = pca_list.index(ipca) + if jj == 0: # Use bounds for PCA0 [avoids negative values] + xmnx = xstat_b.perc( pca_coeff[ipca][idx], per=0.95 ) + PCA_rand[jj,:] = np.random.uniform( xmnx[0], xmnx[1], Ntempl*2) + else: + PCA_mean[jj] = np.mean(pca_coeff[ipca][idx]) + PCA_sig[jj] = np.std(pca_coeff[ipca][idx]) + # Draws + PCA_rand[jj,:] = np.random.uniform( PCA_mean[jj] - 2*PCA_sig[jj], + PCA_mean[jj] + 2*PCA_sig[jj], Ntempl*2) + + # Generate the templates (2*Ntempl) + spec = np.dot(eigen.T,PCA_rand) + + # Take first good Ntempl + + # Truncate, MFP, Fill + ngd = 0 + for kk in range(2*Ntempl): + # Any zero values? + mn = np.min(spec[chkpix,kk]) + if mn < 0.: + continue + # MFP + z912 = wave[0:pix912,kk]/lambda_912 - 1. + phys_dist = np.fabs( cosmo.lookback_distance(z912) - + cosmo.lookback_distance(zrand[kk]) ) # Mpc + #xdb.set_trace() + spec[0:pix912,kk] = spec[0:pix912,kk] * np.exp(-phys_dist.value/mfp[kk]) + # Write + final_spec[:, ii*Ntempl+ngd] = spec[:,kk] + final_wave[:, ii*Ntempl+ngd] = wave[:,kk] + final_z[ii*Ntempl+ngd] = zrand[kk] + ngd += 1 + if ngd == Ntempl: + break + + #xdb.set_trace() + return final_wave, final_spec, final_z + +# ##################### ##################### +# ##################### ##################### +# Plots DESI templates at a range of z and imag +def chk_desi_qso_templates(infil=None, outfil=None, Ntempl=100): + ''' + ''' + # Get the templates + if infil is None: + final_wave, final_spec, final_z = desi_qso_templates(Ntempl=Ntempl) + sz = final_spec.shape + npage = sz[1] // Ntempl + + # Imports + import matplotlib as mpl + mpl.rcParams['font.family'] = 'stixgeneral' + from matplotlib.backends.backend_pdf import PdfPages + from matplotlib import pyplot as plt + import matplotlib.gridspec as gridspec + from matplotlib.colors import LogNorm + + # Eigen (for wavelengths) + xmnx = (3600., 10000.) + + # Start the plot + if outfil != None: + pp = PdfPages(outfil) + + #xdb.set_trace() + + # Looping + for ii in range(npage): + #for ii in range(1): + i0 = ii * Ntempl + i1 = i0 + Ntempl + ymx = 0. + + plt.figure(figsize=(8, 5)) + plt.clf() + gs = gridspec.GridSpec(1, 1) + + # Axis + ax = plt.subplot(gs[0]) + #ax = plt.subplot(gs[ii//2,ii%2]) + + # Labels + ax.set_xlabel('Wavelength') + ax.set_ylabel('Flux') + ax.set_xlim(xmnx) + + # Data + #for jj in range(i0,i1): + for jj in range(i0,i0+20): + ax.plot( final_wave[:,jj], final_spec[:,jj], + '-',drawstyle='steps-mid', linewidth=0.5) + ymx = max( ymx, np.max(final_spec[:,jj]) ) + ax.set_ylim( (0., ymx*1.05) ) + # Label + zmin = np.min(final_z[i0:i1]) + zmax = np.max(final_z[i0:i1]) + zlbl = 'z=[{:g},{:g}]'.format(zmin,zmax) + ax.text(7000., ymx*0.7, zlbl) + + # Layout and save + plt.tight_layout(pad=0.2,h_pad=0.0,w_pad=0.25) + if outfil != None: + pp.savefig()#bbox_inches='tight') + plt.close() + else: + plt.show() + + pp.close() ## ################################# ## ################################# @@ -178,7 +361,9 @@ def fig_desi_templ_z_i(outfil=None, boss_fil=None, flg=0): # Run flg_test = 0 #flg_test += 1 # Mean templates with z,imag - flg_test += 2 # Mean templates with z,imag + #flg_test += 2 # Mean template fig + #flg_test += 2**2 # v1.1 templates for z=2-4 + flg_test += 2**3 # Check v1.1 templates for z=2-4 # Make Mean templates if (flg_test % 2) == 1: @@ -189,6 +374,14 @@ def fig_desi_templ_z_i(outfil=None, boss_fil=None, flg=0): if (flg_test % 2**2) >= 2**1: fig_desi_templ_z_i(outfil='fig_desi_templ_z_i.pdf') + # Make z=2-4 templates; v1.1 + if (flg_test % 2**3) >= 2**2: + aa,bb,cc = desi_qso_templates() + + # Check z=2-4 templates; v1.1 + if (flg_test % 2**4) >= 2**3: + chk_desi_qso_templates(outfil='chk_desi_qso_templates.pdf') + # Done #xdb.set_trace() From 0d231368b403f7ac0783872cadfcb58f713a6ac1 Mon Sep 17 00:00:00 2001 From: "J. Xavier Prochaska" Date: Fri, 19 Dec 2014 06:53:37 -0800 Subject: [PATCH 14/18] jxp --- py/qso_template/boss_qsos_figs.py | 16 +++- py/qso_template/desi_qso_templ.py | 118 +++++++++++++++++++++----- py/qso_template/fit_boss_qsos.py | 134 +++++++++++++++++++++++++----- py/qso_template/run_qso_fits.py | 42 ++++++++++ 4 files changed, 264 insertions(+), 46 deletions(-) create mode 100644 py/qso_template/run_qso_fits.py diff --git a/py/qso_template/boss_qsos_figs.py b/py/qso_template/boss_qsos_figs.py index bb59b032c..f8d08984f 100644 --- a/py/qso_template/boss_qsos_figs.py +++ b/py/qso_template/boss_qsos_figs.py @@ -36,7 +36,7 @@ # ##################### ##################### # ##################### ##################### # Plots the PCA coeff against one another from BOSS DR10 QSOs -def fig_boss_pca_coeff(outfil=None, boss_fil=None): +def fig_boss_pca_coeff(outfil=None, boss_fil=None,scl=100.): # Todo # Include NHI on the label @@ -58,7 +58,6 @@ def fig_boss_pca_coeff(outfil=None, boss_fil=None): # (-2.1, 0.9)) ms = 1. # point size - scl = 100. allxi = [0,0,0,1,1,2] allyi = [1,2,3,2,3,3] @@ -404,7 +403,9 @@ def fig_desi_templ_z_i(outfil=None, boss_fil=None, flg=0): #flg_fig += 2**1 # PCA vs z #flg_fig += 2**2 # PCA vs i #flg_fig += 2**3 # Eigenvectors - flg_fig += 2**4 # z vs i for BOSS + #flg_fig += 2**4 # z vs i for BOSS + + flg_fig += 2**10 # SDSS: PCA vs PCA # PCA vs PCA if (flg_fig % 2) == 1: @@ -425,3 +426,12 @@ def fig_desi_templ_z_i(outfil=None, boss_fil=None, flg=0): # z vs i if (flg_fig % 2**5) >= 2**4: fig_boss_i_vs_z(outfil='fig_boss_i_vs_z.pdf') + + # ####### + # SDSS FIGURES + + # ####### + # PCA vs PCA + sdss_fil = 'SDSS_DR7Lya_PCA_values_nocut.fits' + if (flg_fig % 2**11) >= 2**10: + fig_boss_pca_coeff(outfil='fig_sdss_pca.pdf', boss_fil=sdss_fil, scl=1.) diff --git a/py/qso_template/desi_qso_templ.py b/py/qso_template/desi_qso_templ.py index 90346fe98..edfb32200 100644 --- a/py/qso_template/desi_qso_templ.py +++ b/py/qso_template/desi_qso_templ.py @@ -57,7 +57,6 @@ def mean_templ_zi(zimag, debug=False, i_wind=0.1, z_wind=0.05, boss_cat_fil = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_cat_v2.1.fits.gz' bcat_hdu = fits.open(boss_cat_fil) t_boss = bcat_hdu[1].data - nqso = len(t_boss) zQSO = t_boss['z_pipe'] tmp = t_boss['PSFMAG'] imag = tmp[:,3] # i-band mag @@ -171,30 +170,38 @@ def fig_desi_templ_z_i(outfil=None, boss_fil=None, flg=0): plt.show() ## -def desi_qso_templates(z_wind=0.2, zmnx=(2.,4.), outfil=None, Ntempl=100, - boss_pca_fil=None): +def desi_qso_templates(z_wind=0.2, zmnx=(0.4,4.), outfil=None, Ntempl=500, + boss_pca_fil=None, wvmnx=(3500., 10000.), + sdss_pca_fil=None, no_write=False): ''' - Generate 'mean' templates at given z,i + Generate 10000 templates for DESI from z=0.4 to 4 Parameters ---------- z_wind: float (0.2) Window for sampling - zmnx: tuple ( (2,4) ) + zmnx: tuple ( (0.5,4) ) Min/max for generation - Ntempl: int (100) + Ntempl: int (500) Number of draws per redshift window ''' # Cosmology from astropy import cosmology cosmo = cosmology.core.FlatLambdaCDM(70., 0.3) + # PCA values if boss_pca_fil is None: boss_pca_fil = 'BOSS_DR10Lya_PCA_values_nocut.fits.gz' hdu = fits.open(boss_pca_fil) - pca_coeff = hdu[1].data + boss_pca_coeff = hdu[1].data - # BOSS Eigenvectors + if sdss_pca_fil is None: + sdss_pca_fil = 'SDSS_DR7Lya_PCA_values_nocut.fits' + hdu2 = fits.open(sdss_pca_fil) + sdss_pca_coeff = hdu2[1].data + + + # Eigenvectors eigen, eigen_wave = fbq.read_qso_eigen() npix = len(eigen_wave) chkpix = np.where((eigen_wave > 900.) & (eigen_wave < 5000.) )[0] @@ -205,9 +212,16 @@ def desi_qso_templates(z_wind=0.2, zmnx=(2.,4.), outfil=None, Ntempl=100, boss_cat_fil = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_cat_v2.1.fits.gz' bcat_hdu = fits.open(boss_cat_fil) t_boss = bcat_hdu[1].data - nqso = len(t_boss) - zQSO = t_boss['z_pipe'] + boss_zQSO = t_boss['z_pipe'] + # Open the SDSS catalog file + sdss_cat_fil = os.environ.get('SDSSPATH')+'/DR7_QSO/dr7_qso.fits.gz' + scat_hdu = fits.open(sdss_cat_fil) + t_sdss = scat_hdu[1].data + sdss_zQSO = t_sdss['z'] + if len(sdss_pca_coeff) != len(sdss_zQSO): + print('Need to run all SDSS models!') + sdss_zQSO = sdss_zQSO[0:len(sdss_pca_coeff)] # Outfil if outfil is None: @@ -228,6 +242,15 @@ def desi_qso_templates(z_wind=0.2, zmnx=(2.,4.), outfil=None, Ntempl=100, seed = -1422 for ii in range(len(z0)): + + # BOSS or SDSS? + if z0[ii] > 1.99: + zQSO = boss_zQSO + pca_coeff = boss_pca_coeff + else: + zQSO = sdss_zQSO + pca_coeff = sdss_pca_coeff + # Random z values and wavelengths zrand = np.random.uniform( z0[ii], z1[ii], Ntempl*2) wave = np.outer(eigen_wave, 1+zrand) @@ -237,6 +260,7 @@ def desi_qso_templates(z_wind=0.2, zmnx=(2.,4.), outfil=None, Ntempl=100, # Grab PCA mean + sigma idx = np.where( (zQSO >= z0[ii]) & (zQSO < z1[ii]) )[0] + print('Making z=({:g},{:g}) with {:d} input quasars'.format(z0[ii],z1[ii],len(idx))) # Get PCA stats and random values for ipca in pca_list: @@ -263,12 +287,14 @@ def desi_qso_templates(z_wind=0.2, zmnx=(2.,4.), outfil=None, Ntempl=100, mn = np.min(spec[chkpix,kk]) if mn < 0.: continue + # MFP - z912 = wave[0:pix912,kk]/lambda_912 - 1. - phys_dist = np.fabs( cosmo.lookback_distance(z912) - - cosmo.lookback_distance(zrand[kk]) ) # Mpc - #xdb.set_trace() - spec[0:pix912,kk] = spec[0:pix912,kk] * np.exp(-phys_dist.value/mfp[kk]) + if z0[ii] > 2.39: + z912 = wave[0:pix912,kk]/lambda_912 - 1. + phys_dist = np.fabs( cosmo.lookback_distance(z912) - + cosmo.lookback_distance(zrand[kk]) ) # Mpc + spec[0:pix912,kk] = spec[0:pix912,kk] * np.exp(-phys_dist.value/mfp[kk]) + # Write final_spec[:, ii*Ntempl+ngd] = spec[:,kk] final_wave[:, ii*Ntempl+ngd] = wave[:,kk] @@ -277,7 +303,52 @@ def desi_qso_templates(z_wind=0.2, zmnx=(2.,4.), outfil=None, Ntempl=100, if ngd == Ntempl: break - #xdb.set_trace() + if no_write is True: # Mainly for plotting + return final_wave, final_spec, final_z + + # Rebin + light = 2.99792458e5 # [km/s] + velpixsize = 10. # [km/s] + pixsize = velpixsize/light/np.log(10) # [pixel size in log-10 A] + minwave = np.log10(wvmnx[0]) # minimum wavelength [log10-A] + maxwave = np.log10(wvmnx[1]) # maximum wavelength [log10-A] + r_npix = np.round((maxwave-minwave)/pixsize+1) + + log_wave = minwave+np.arange(r_npix)*pixsize # constant log-10 spacing + + totN = Ntempl * len(z0) + rebin_spec = np.zeros((r_npix, totN)) + + from scipy.interpolate import interp1d + + for ii in range(totN): + # Interpolate (in log space) + f1d = interp1d(np.log10(final_wave[:,ii]), final_spec[:,ii]) + rebin_spec[:,ii] = f1d(log_wave) + #xdb.xplot(final_wave[:,ii], final_spec[:,ii], xtwo=10.**log_wave, ytwo=rebin_spec[:,ii]) + #xdb.set_trace() + + # Write + hdu = fits.PrimaryHDU(rebin_spec) + hdu.header.set('PROJECT', 'DESI QSO TEMPLATES') + hdu.header.set('VERSION', '1.1') + hdu.header.set('OBJTYPE', 'QSO') + hdu.header.set('DISPAXIS', 1, 'dispersion axis') + hdu.header.set('CRPIX1', 1, 'reference pixel number') + hdu.header.set('CRVAL1', minwave, 'reference log10(Ang)') + hdu.header.set('CDELT1', pixsize, 'delta log10(Ang)') + hdu.header.set('LOGLAM', 1, 'log10 spaced wavelengths?') + hdu.header.set('AIRORVAC', 'vac', ' wavelengths in vacuum (vac) or air') + hdu.header.set('VELSCALE', velpixsize, ' pixel size in km/s') + hdu.header.set('WAVEUNIT', 'Angstrom', ' wavelength units') + + prihdu = fits.PrimaryHDU(header=prihdr) + + table_hdu = fits.BinTableHDU.from_columns(np.array(full_tab.filled())) + thdulist = fits.HDUList([prihdu, table_hdu]) + print('Writing {:s} table, with {:d} rows'.format(outfil,len(full_tab))) + thdulist.writeto(outfil, clobber=True) + return final_wave, final_spec, final_z # ##################### ##################### @@ -288,7 +359,8 @@ def chk_desi_qso_templates(infil=None, outfil=None, Ntempl=100): ''' # Get the templates if infil is None: - final_wave, final_spec, final_z = desi_qso_templates(Ntempl=Ntempl) + final_wave, final_spec, final_z = desi_qso_templates(Ntempl=Ntempl, zmnx=(0.4,0.8), + no_write=True) sz = final_spec.shape npage = sz[1] // Ntempl @@ -331,7 +403,7 @@ def chk_desi_qso_templates(infil=None, outfil=None, Ntempl=100): # Data #for jj in range(i0,i1): - for jj in range(i0,i0+20): + for jj in range(i0,i0+15): ax.plot( final_wave[:,jj], final_spec[:,jj], '-',drawstyle='steps-mid', linewidth=0.5) ymx = max( ymx, np.max(final_spec[:,jj]) ) @@ -362,8 +434,8 @@ def chk_desi_qso_templates(infil=None, outfil=None, Ntempl=100): flg_test = 0 #flg_test += 1 # Mean templates with z,imag #flg_test += 2 # Mean template fig - #flg_test += 2**2 # v1.1 templates for z=2-4 - flg_test += 2**3 # Check v1.1 templates for z=2-4 + flg_test += 2**2 # v1.1 templates + #flg_test += 2**3 # Check v1.1 templates # Make Mean templates if (flg_test % 2) == 1: @@ -376,11 +448,11 @@ def chk_desi_qso_templates(infil=None, outfil=None, Ntempl=100): # Make z=2-4 templates; v1.1 if (flg_test % 2**3) >= 2**2: - aa,bb,cc = desi_qso_templates() + aa,bb,cc = desi_qso_templates(zmnx=(2.,2.4)) - # Check z=2-4 templates; v1.1 + # Check z=0.4-4 templates; v1.1 if (flg_test % 2**4) >= 2**3: - chk_desi_qso_templates(outfil='chk_desi_qso_templates.pdf') + chk_desi_qso_templates(outfil='chk_desi_qso_templates.pdf', Ntempl=20) # Done diff --git a/py/qso_template/fit_boss_qsos.py b/py/qso_template/fit_boss_qsos.py index 70c31bc78..8775c06d1 100644 --- a/py/qso_template/fit_boss_qsos.py +++ b/py/qso_template/fit_boss_qsos.py @@ -93,9 +93,11 @@ def do_boss_lya_parallel(istart, iend, cut_Lya, output, debug=False): #for ii in range(nqso): datdir = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_spectra_v2.1/' jj = 0 + print('istart = {:d}'.format(istart)) for ii in range(istart,iend): - if (ii % 1000) == 0: + if (ii % 100) == 0: print('ii = {:d}'.format(ii)) + #print('ii = {:d}'.format(ii)) # Spectrum file pnm = str(t_boss['PLATE'][ii]) fnm = str(t_boss['FIBERID'][ii]).rjust(4,str('0')) @@ -125,8 +127,11 @@ def do_boss_lya_parallel(istart, iend, cut_Lya, output, debug=False): imx = npix+imn eigen_flux = eigen[:,imn:imx] + # FIT - acoeff = fit_eigen(flux[Ly_imn:], ivar[Ly_imn:], eigen_flux) + tflux = flux[Ly_imn:] + tivar = ivar[Ly_imn:] + acoeff = fit_eigen(tflux, ivar, eigen_flux) pca_val[jj,:] = acoeff jj += 1 @@ -137,11 +142,14 @@ def do_boss_lya_parallel(istart, iend, cut_Lya, output, debug=False): xdb.xplot(wrest, flux, xtwo=eigen_wave, ytwo=model) xdb.set_trace() + #xdb.set_trace() print('Done with my subset {:d}, {:d}'.format(istart,iend)) if output is not None: output.put((istart,iend,pca_val)) #output.put(None) + else: + return pca_val ## def do_sdss_lya_parallel(istart, iend, cut_Lya, output, debug=False): @@ -219,7 +227,6 @@ def do_sdss_lya_parallel(istart, iend, cut_Lya, output, debug=False): # Check if debug is True: model = np.dot(eigen.T,acoeff) - xdb.set_trace() if flg_xdb is True: xdb.xplot(wrest, flux, xtwo=eigen_wave, ytwo=model) xdb.set_trace() @@ -229,17 +236,14 @@ def do_sdss_lya_parallel(istart, iend, cut_Lya, output, debug=False): if output is not None: output.put((istart,iend,pca_val)) #output.put(None) - -## ################################# -## ################################# -## TESTING -## ################################# -if __name__ == '__main__': - - # Run - #do_boss_lya_parallel(0,10,None,debug=True,cut_Lya=True) - #xdb.set_trace() + else: + return pca_val +def failed_parallel(): + ''' + Collision with np.dot + Might fix with OPENBLAS_NUM_THREADS=1 + ''' flg = 0 # 0=BOSS, 1=SDSS ## ############################ @@ -256,14 +260,13 @@ def do_sdss_lya_parallel(istart, iend, cut_Lya, output, debug=False): nqso = len(t_sdss) outfil = 'SDSS_DR7Lya_PCA_values_nocut.fits' - nqso = 450 # Testing + nqso = 40 # Testing - #do_sdss_lya_parallel(0, nqso, False, None) - #xdb.set_trace() + #do_boss_lya_parallel(0,nqso, False, None,debug=False) output = mp.Queue() processes = [] - nproc = 8 + nproc = 1 nsub = nqso // nproc cut_Lya = False @@ -289,18 +292,17 @@ def do_sdss_lya_parallel(istart, iend, cut_Lya, output, debug=False): for p in processes: p.start() - #xdb.set_trace() - # Get process results from the output queue print('Grabbing Output') results = [output.get() for p in processes] - xdb.set_trace() + # Get process results from the output queue # Exit the completed processes print('Joining') for p in processes: p.join() + xdb.set_trace() # Bring together #sorted(results, key=lambda result: result[0]) #all_is = [ir[0] for ir in results] @@ -332,3 +334,95 @@ def do_sdss_lya_parallel(istart, iend, cut_Lya, output, debug=False): # Done #xdb.set_trace() print('All done') + +def stack_fits(flg=0): + ''' + Splices together the various PCA fits for SDSS or BOSS + + flg: int (0) + 0=BOSS, 1=SDSS + ''' + import glob + from astropy.table.table import Table + + if flg == 0: + outroot = 'Output/BOSS_DR10Lya_PCA_values_nocut' + outfil = 'BOSS_DR10Lya_PCA_values_nocut.fits' + elif flg ==1: + outroot = 'Output/SDSS_DR7Lya_PCA_values_nocut' + outfil = 'SDSS_DR7Lya_PCA_values_nocut.fits' + + # Get all the files + files = glob.glob(outroot+'*') + + for ifil in files: + hdu = fits.open(ifil) + tab = Table(hdu[1].data) + # + if not 'full_tab' in locals(): + full_tab = tab + else: + #xdb.set_trace() + for row in tab: + full_tab.add_row(row) + # Write + prihdr = fits.Header() + if flg == 0: + prihdr['PROJECT'] = 'BOSS: z>2 quasars' + elif flg == 1: + prihdr['PROJECT'] = 'SDSS: Meant for z<2 quasars' + prihdr['COMMENT'] = 'PCA fits to the quasars' + prihdu = fits.PrimaryHDU(header=prihdr) + + table_hdu = fits.BinTableHDU.from_columns(np.array(full_tab.filled())) + thdulist = fits.HDUList([prihdu, table_hdu]) + print('Writing {:s} table, with {:d} rows'.format(outfil,len(full_tab))) + thdulist.writeto(outfil, clobber=True) + + +## ################ +if __name__ == '__main__': + ''' + flg: int (0) + 0=BOSS, 1=SDSS + ''' + import sys + + flg = int(sys.argv[1]) + istrt = int(sys.argv[2]) + iend = int(sys.argv[3]) + outfil = str(sys.argv[4]) + + cut_Lya = False + + # Run + #do_boss_lya_parallel(0,10, False, None,debug=True) + #xdb.set_trace() + #do_sdss_lya_parallel(0, 10, False, None, debug=True) + + ## ############################ + if flg == 0: + pca_val = do_boss_lya_parallel(istrt,iend,cut_Lya, None) + elif flg == 1: + pca_val = do_sdss_lya_parallel(istrt,iend,cut_Lya, None) + + # Write to disk as a binary FITS table + col0 = fits.Column(name='PCA0',format='E',array=pca_val[:,0]) + col1 = fits.Column(name='PCA1',format='E',array=pca_val[:,1]) + col2 = fits.Column(name='PCA2',format='E',array=pca_val[:,2]) + col3 = fits.Column(name='PCA3',format='E',array=pca_val[:,3]) + cols = fits.ColDefs([col0, col1, col2, col3]) + tbhdu = fits.BinTableHDU.from_columns(cols) + + prihdr = fits.Header() + prihdr['OBSERVER'] = 'Edwin Hubble' + prihdr['COMMENT'] = "Here's some commentary about this FITS file." + prihdu = fits.PrimaryHDU(header=prihdr) + + thdulist = fits.HDUList([prihdu, tbhdu]) + thdulist.writeto(outfil, clobber=True) + + # Done + #xdb.set_trace() + print('All done') + diff --git a/py/qso_template/run_qso_fits.py b/py/qso_template/run_qso_fits.py new file mode 100644 index 000000000..f4841a700 --- /dev/null +++ b/py/qso_template/run_qso_fits.py @@ -0,0 +1,42 @@ +"popen_ex.py" +from subprocess import Popen +import os +from astropy.io import fits + +flg = 1 # SDSS +nproc = 4 + +## ############################ +# "Parallel" +if flg == 0: + print('Running BOSS!') + boss_cat_fil = os.environ.get('BOSSPATH')+'/DR10/BOSSLyaDR10_cat_v2.1.fits.gz' + bcat_hdu = fits.open(boss_cat_fil) + t_boss = bcat_hdu[1].data + nqso = len(t_boss) + outroot = 'Output/BOSS_DR10Lya_PCA_values_nocut' +elif flg == 1: + print('Running SDSS!') + sdss_cat_fil = os.environ.get('SDSSPATH')+'/DR7_QSO/dr7_qso.fits.gz' + scat_hdu = fits.open(sdss_cat_fil) + t_sdss = scat_hdu[1].data + nqso = len(t_sdss) + outroot = 'Output/SDSS_DR7Lya_PCA_values_nocut' + +nqso = 20000 # Testing +nsub = nqso // nproc + +cut_Lya = False + +# Setup the Processes +for ii in range(nproc): + # Generate + istrt = ii * nsub + if ii == (nproc-1): + iend = nqso + else: + iend = (ii+1)*nsub + + outfil = outroot+str(ii)+'.fits' + Popen(['python', './fit_boss_qsos.py', str(flg), str(istrt), str(iend), outfil]) + From 9fc4c2a20c26dabf74a668864f455082ecb8d97c Mon Sep 17 00:00:00 2001 From: "J. Xavier Prochaska" Date: Fri, 19 Dec 2014 08:09:43 -0800 Subject: [PATCH 15/18] jxp --- py/qso_template/desi_qso_templ.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/py/qso_template/desi_qso_templ.py b/py/qso_template/desi_qso_templ.py index edfb32200..3181a17e3 100644 --- a/py/qso_template/desi_qso_templ.py +++ b/py/qso_template/desi_qso_templ.py @@ -342,12 +342,14 @@ def desi_qso_templates(z_wind=0.2, zmnx=(0.4,4.), outfil=None, Ntempl=500, hdu.header.set('VELSCALE', velpixsize, ' pixel size in km/s') hdu.header.set('WAVEUNIT', 'Angstrom', ' wavelength units') - prihdu = fits.PrimaryHDU(header=prihdr) - - table_hdu = fits.BinTableHDU.from_columns(np.array(full_tab.filled())) - thdulist = fits.HDUList([prihdu, table_hdu]) - print('Writing {:s} table, with {:d} rows'.format(outfil,len(full_tab))) - thdulist.writeto(outfil, clobber=True) + idval = range(totN) + col0 = fits.Column(name='idval',format='K', array=idval) + col1 = fits.Column(name='zQSO',format='E',array=final_z) + cols = fits.ColDefs([col0, col1]) + tbhdu = fits.BinTableHDU.from_columns(cols) + + hdulist = fits.HDUList([hdu, tbhdu]) + hdulist.writeto(outfil, clobber=True) return final_wave, final_spec, final_z From c32aedf805ac97780eac5e92cafe5f774763b7e7 Mon Sep 17 00:00:00 2001 From: "J. Xavier Prochaska" Date: Fri, 19 Dec 2014 08:52:10 -0800 Subject: [PATCH 16/18] jxp --- py/qso_template/desi_qso_templ.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/py/qso_template/desi_qso_templ.py b/py/qso_template/desi_qso_templ.py index 3181a17e3..235f474a2 100644 --- a/py/qso_template/desi_qso_templ.py +++ b/py/qso_template/desi_qso_templ.py @@ -328,8 +328,11 @@ def desi_qso_templates(z_wind=0.2, zmnx=(0.4,4.), outfil=None, Ntempl=500, #xdb.xplot(final_wave[:,ii], final_spec[:,ii], xtwo=10.**log_wave, ytwo=rebin_spec[:,ii]) #xdb.set_trace() + # Transpose for consistency + out_spec = np.array(rebin_spec.T, dtype='float32') + # Write - hdu = fits.PrimaryHDU(rebin_spec) + hdu = fits.PrimaryHDU(out_spec) hdu.header.set('PROJECT', 'DESI QSO TEMPLATES') hdu.header.set('VERSION', '1.1') hdu.header.set('OBJTYPE', 'QSO') @@ -343,10 +346,11 @@ def desi_qso_templates(z_wind=0.2, zmnx=(0.4,4.), outfil=None, Ntempl=500, hdu.header.set('WAVEUNIT', 'Angstrom', ' wavelength units') idval = range(totN) - col0 = fits.Column(name='idval',format='K', array=idval) - col1 = fits.Column(name='zQSO',format='E',array=final_z) + col0 = fits.Column(name='TEMPLATEID',format='K', array=idval) + col1 = fits.Column(name='Z',format='E',array=final_z) cols = fits.ColDefs([col0, col1]) tbhdu = fits.BinTableHDU.from_columns(cols) + tbhdu.header.set('EXTNAME','METADATA') hdulist = fits.HDUList([hdu, tbhdu]) hdulist.writeto(outfil, clobber=True) From 83e83ff9682c5021c0267212f9fa9b8e7c37dc3a Mon Sep 17 00:00:00 2001 From: "J. Xavier Prochaska" Date: Fri, 19 Dec 2014 11:45:42 -0800 Subject: [PATCH 17/18] jxp --- py/qso_template/desi_qso_templ.py | 13 ++++++++----- py/qso_template/fit_boss_qsos.py | 4 +++- py/qso_template/run_qso_fits.py | 4 ++-- 3 files changed, 13 insertions(+), 8 deletions(-) diff --git a/py/qso_template/desi_qso_templ.py b/py/qso_template/desi_qso_templ.py index 235f474a2..0ef1095dc 100644 --- a/py/qso_template/desi_qso_templ.py +++ b/py/qso_template/desi_qso_templ.py @@ -174,7 +174,7 @@ def desi_qso_templates(z_wind=0.2, zmnx=(0.4,4.), outfil=None, Ntempl=500, boss_pca_fil=None, wvmnx=(3500., 10000.), sdss_pca_fil=None, no_write=False): ''' - Generate 10000 templates for DESI from z=0.4 to 4 + Generate 9000 templates for DESI from z=0.4 to 4 Parameters ---------- @@ -196,7 +196,7 @@ def desi_qso_templates(z_wind=0.2, zmnx=(0.4,4.), outfil=None, Ntempl=500, boss_pca_coeff = hdu[1].data if sdss_pca_fil is None: - sdss_pca_fil = 'SDSS_DR7Lya_PCA_values_nocut.fits' + sdss_pca_fil = 'SDSS_DR7Lya_PCA_values_nocut.fits.gz' hdu2 = fits.open(sdss_pca_fil) sdss_pca_coeff = hdu2[1].data @@ -220,7 +220,7 @@ def desi_qso_templates(z_wind=0.2, zmnx=(0.4,4.), outfil=None, Ntempl=500, t_sdss = scat_hdu[1].data sdss_zQSO = t_sdss['z'] if len(sdss_pca_coeff) != len(sdss_zQSO): - print('Need to run all SDSS models!') + print('Need to finish running the SDSS models!') sdss_zQSO = sdss_zQSO[0:len(sdss_pca_coeff)] # Outfil @@ -302,6 +302,9 @@ def desi_qso_templates(z_wind=0.2, zmnx=(0.4,4.), outfil=None, Ntempl=500, ngd += 1 if ngd == Ntempl: break + if ngd != Ntempl: + print('Did not make enough!') + xdb.set_trace() if no_write is True: # Mainly for plotting return final_wave, final_spec, final_z @@ -365,7 +368,7 @@ def chk_desi_qso_templates(infil=None, outfil=None, Ntempl=100): ''' # Get the templates if infil is None: - final_wave, final_spec, final_z = desi_qso_templates(Ntempl=Ntempl, zmnx=(0.4,0.8), + final_wave, final_spec, final_z = desi_qso_templates(Ntempl=Ntempl, #zmnx=(0.4,0.8), no_write=True) sz = final_spec.shape npage = sz[1] // Ntempl @@ -454,7 +457,7 @@ def chk_desi_qso_templates(infil=None, outfil=None, Ntempl=100): # Make z=2-4 templates; v1.1 if (flg_test % 2**3) >= 2**2: - aa,bb,cc = desi_qso_templates(zmnx=(2.,2.4)) + aa,bb,cc = desi_qso_templates() #zmnx=(2.,2.4)) # Check z=0.4-4 templates; v1.1 if (flg_test % 2**4) >= 2**3: diff --git a/py/qso_template/fit_boss_qsos.py b/py/qso_template/fit_boss_qsos.py index 8775c06d1..501dfc3c4 100644 --- a/py/qso_template/fit_boss_qsos.py +++ b/py/qso_template/fit_boss_qsos.py @@ -335,7 +335,8 @@ def failed_parallel(): #xdb.set_trace() print('All done') -def stack_fits(flg=0): +# ######## +def splice_fits(flg=0): ''' Splices together the various PCA fits for SDSS or BOSS @@ -356,6 +357,7 @@ def stack_fits(flg=0): files = glob.glob(outroot+'*') for ifil in files: + print('Reading {:s}'.format(ifil)) hdu = fits.open(ifil) tab = Table(hdu[1].data) # diff --git a/py/qso_template/run_qso_fits.py b/py/qso_template/run_qso_fits.py index f4841a700..e501a9162 100644 --- a/py/qso_template/run_qso_fits.py +++ b/py/qso_template/run_qso_fits.py @@ -4,7 +4,7 @@ from astropy.io import fits flg = 1 # SDSS -nproc = 4 +nproc = 8 ## ############################ # "Parallel" @@ -23,7 +23,7 @@ nqso = len(t_sdss) outroot = 'Output/SDSS_DR7Lya_PCA_values_nocut' -nqso = 20000 # Testing +#nqso = 800 #20000 # Testing nsub = nqso // nproc cut_Lya = False From 489df193b2dbf53a8747da50de6d8c826b061e6c Mon Sep 17 00:00:00 2001 From: "J. Xavier Prochaska" Date: Mon, 5 Jan 2015 11:11:14 -0800 Subject: [PATCH 18/18] bunit --- py/qso_template/desi_qso_templ.py | 1 + 1 file changed, 1 insertion(+) diff --git a/py/qso_template/desi_qso_templ.py b/py/qso_template/desi_qso_templ.py index 0ef1095dc..7f4de30f2 100644 --- a/py/qso_template/desi_qso_templ.py +++ b/py/qso_template/desi_qso_templ.py @@ -347,6 +347,7 @@ def desi_qso_templates(z_wind=0.2, zmnx=(0.4,4.), outfil=None, Ntempl=500, hdu.header.set('AIRORVAC', 'vac', ' wavelengths in vacuum (vac) or air') hdu.header.set('VELSCALE', velpixsize, ' pixel size in km/s') hdu.header.set('WAVEUNIT', 'Angstrom', ' wavelength units') + hdu.header.set('BUNIT', '1e-17 erg/s/cm2/A', ' flux unit') idval = range(totN) col0 = fits.Column(name='TEMPLATEID',format='K', array=idval)