imcgreer / simqso Goto Github PK

View Code? Open in Web Editor NEW

12.0 12.0 11.0 6.43 MB

Module for generating simulated quasar spectra

License: BSD 3-Clause "New" or "Revised" License

Python 99.32% IDL 0.42% Makefile 0.25%

simqso's People

Contributors

Stargazers

Watchers

Forkers

yalsayyad findaz jtschindler desihub swayamtrupta sbailey astroland93 muskaan-chaudhary

simqso's Issues

Lyman-alpha continuum different than data

With desi we observe a large difference between the Lyman-alpha forest continuum in SIMQSO and in BOSS+eBOSS data. The following plot gives the stack for both cases. Would it be possible for you @imcgreer to update the lines parameters to reproduce BOSS+eBOSS stack and its variation?

Edit: The differences can be observed in many places of the spectrum but are more important
there.
Edit2: The flux is normalized in this plot.

add support for pysynphot

add BASS + MzLS filters

For DESI we need simqso to know about the BASS + MzLS filter profiles, which can be found here--
https://github.com/dkirkby/speclite/tree/master/speclite/data/filters

mag-dependent k-corr within sims

in order to use emission line kcorr and reduce number of iterations.

notebook for color-z and k-correction

reorganize forest sightlines

make a container object (with special case for grid sightlines?) that can generate sightlines on-the-fly, so they aren't required to be precomputed

definition of speed of light

@moustakas and @imcgreer, Doing a simple grep '3e5' -r * --include \*.py, I see the following in this repository and its fork https://github.com/moustakas/simqso:

imcgreer/simqso/simqso/sqgrids.py:        c_kms = 3e5
moustakas/simqso/simqso/sqgrids.py:        c_kms = 3e5

Would it be possible to change to the more universal from scipy.constants import speed_of_light or c = 299792458. m/s?
Maybe other grep would show other definition of the speed of light, could you check for that also?

package trying to read file not included in data

Hi,

When I was trying to generate spectra like in the example I got the following error: [Error] No such file or directory: data/VW01_Fe/Fe_UVOPT_V01_T06_BR92.asc.

It looks like the folder VW01_Fe was removed in a previous commit but it is still used by the package.

Is there any workaround to this problem?

power law continuum model incorrectly handling breakpoints

FixedPLContinuumGrid does nothing to the input breakpoints.
GaussianPLContinuumGrid adds '0' to the beginning so num(breaks)=num(slopes)

But the way the iteration works, it should be asserted in setPowerLawContinuum that num(breaks) = num(slopes)-1, then breaks = [0,breaks,np.inf]

(what is currently happening with default models? last slope is being dropped?)

emission lines parameters

I have used SIMQSO spectra to assert desi's redshift fitter performance.
In PR desihub/redrock#104, I seem to see that the parameters you use for all emission lines blue of the Lya (including Lya) have been extracted from absorbed spectra. Is it true? How do you get the parameters giving your Lya emission line in SIMQSO? Would it be easy to get these parameters from unabsorbed data? i.e. data corrected for absorption.

request to add DECaLS photometric system

The filter curves are here --
https://github.com/dkirkby/speclite/tree/master/speclite/data/filters
and documented here --
http://speclite.readthedocs.io/en/stable/filters.html#decam-filters

This is needed so that DESI can use simqso.

use random number seed to ensure reproducibility

The random sampling that occurs in lumfun and many other places should take as an option a numpy.random.RandomState() object or an input seed, to ensure reproducibility.

Let me know if you want me to add this @imcgreer or whether you would rather. Also, this should probably happen after #15 has been merged.

sqgrids.QsoSimObjects.read should not ignore EXTNAME

sqgrids.QsoSimObjects.read should read the header from the specified EXTNAME on this line instead of reading the first HDU.

Also, I'm not sure why the header is read again here

add warm dust component

multiple instances of QsoSimGrid results in crashes

Something is being held in global memory. This patch results in a crash if multiple sims are executed:

diff --git a/sdss/ebosscore.py b/sdss/ebosscore.py
index 4576766..ffb3d26 100644
--- a/sdss/ebosscore.py
+++ b/sdss/ebosscore.py
@@ -257,10 +257,10 @@ def run_colorz_sim(model,nm=7,nz=500):
    zrange = (0.9,4.0)
    mbins = np.linspace(*tuple(mrange+(nm,)))
    zbins = np.linspace(*tuple(zrange+(nz,)))
-   M,z = np.meshgrid(mbins,zbins,indexing='ij')
-   M = grids.AbsMagVar(grids.FixedSampler(M.flatten()),restWave=1450)
-   z = grids.RedshiftVar(grids.FixedSampler(z.flatten()))
-   qsos = grids.QsoSimPoints([M,z],cosmo=dr9cosmo,units='luminosity')
+   M = grids.AbsMagVar(grids.FixedSampler(mbins),restWave=1450)
+   z = grids.RedshiftVar(grids.FixedSampler(zbins))
+   qsos = grids.QsoSimGrid([M,z],(None,None),1,fixed_vars=['absMag','z'],
+                           cosmo=dr9cosmo,units='luminosity')
    qsos = runsim(model,None,'sdss_forest_grid.fits',qsos,
                  medianforest=True,const=True,nophot=True)
    synmags = np.array(qsos.data['synMag'].reshape(nm,nz,-1))

Wrong datatype assumed for sqgrids.py - blackbody_lambda function

Hi, I'm new to this package and trying to replicate the examples. I'm running astropy v5.0.4.

The example I'm following is this: https://github.com/imcgreer/simqso/blob/master/examples/SimpleSpecExample.ipynb
I haven't tried other examples yet. Also, I didn't do the %pylab inline in the beginning as I usually don't use pylab at all, and hence I imported numpy and matplotlib by myself and edited some parts of the example to correctly refer to intended numpy/matplotlib functions, which shouldn't be causing the problem I experienced.

Upon calling the 9th cell in the example;

# ready to generate spectra
_,spectra = buildSpectraBulk(wave,qsos,saveSpectra=True)

I get the error

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-9-8d131b44e14a> in <module>
      1 # ready to generate spectra
----> 2 _,spectra = buildSpectraBulk(wave,qsos,saveSpectra=True)

/opt/miniconda3/lib/python3.8/site-packages/simqso/sqrun.py in buildSpectraBulk(wave, qsoGrid, procMap, maxIter, verbose, saveSpectra)
    431         if verbose > 1:
    432             print('buildSpectra iteration ',iterNum,' out of ',nIter)
--> 433         specOut = list(procMap(build_one_spec,qsoGrid))
    434         specOut = _regroup(specOut)
    435         synMag,synFlux,spectra = specOut

/opt/miniconda3/lib/python3.8/site-packages/simqso/sqrun.py in buildSpecWithPhot(wave, cosmo, specFeatures, photoCache, objData, iterNum, saveSpectra)
    386 def buildSpecWithPhot(wave,cosmo,specFeatures,photoCache,
    387                       objData,iterNum=None,saveSpectra=False):
--> 388     sp = buildQsoSpectrum(wave,cosmo,specFeatures,objData,
    389                           iterNum=iterNum)
    390     if photoCache is None:

/opt/miniconda3/lib/python3.8/site-packages/simqso/sqrun.py in buildQsoSpectrum(wave, cosmo, specFeatures, obj, iterNum, save_components)
    250         if isinstance(feature,grids.ContinuumVar):
    251             assocvals = _getpar(feature.get_associated_var(),obj)
--> 252             spec = feature.add_to_spec(spec,_getpar(feature,obj),
    253                                        assocvals=assocvals,
    254                                        fluxNorm=fluxNorm)

/opt/miniconda3/lib/python3.8/site-packages/simqso/sqgrids.py in add_to_spec(self, spec, par, **kwargs)
    472         par : sampled values of the variable that are passed to render()
    473         '''
--> 474         spec.f_lambda[:] += self.render(spec.wave,spec.z,par,**kwargs)
    475         return spec
    476 

/opt/miniconda3/lib/python3.8/site-packages/simqso/sqgrids.py in render(self, wave, z, par, fluxNorm, assocvals)
    755         fdisk = self.assocVar.total_flux(assocvals,fluxNorm,z)
    756         flux_bb = fracdust * fdisk
--> 757         Blam = self._calc_bb(Tdust,wave,z)
    758         return flux_bb * np.pi * Blam / L_bb
    759 

/opt/miniconda3/lib/python3.8/site-packages/simqso/sqgrids.py in _calc_bb(self, Tdust, wave, z)
    740             rfwave = np.concatenate([ (lampeak-dwv1)[::-1], rfwv2 ])
    741             self.rfwave[Tdust] = rfwave
--> 742             bvals = blackbody_lambda(self.rfwave[Tdust],Tdust).value
    743             self.Blam[Tdust] = interp1d(self.rfwave[Tdust],bvals,
    744                                         kind='cubic')

/opt/miniconda3/lib/python3.8/site-packages/simqso/sqgrids.py in blackbody_lambda(in_x, temperature)
     44             return flux / u.sr  # Add per steradian to output flux unit
     45         bb_nu = blackbody_nu(in_x, temperature) * u.sr  # Remove sr for conversion
---> 46         flux = bb_nu.to(FLAM, u.spectral_density(in_x))
     47         return flux / u.sr  # Add per steradian to output flux unit
     48 

/opt/miniconda3/lib/python3.8/site-packages/astropy/units/quantity.py in to(self, unit, equivalencies, copy)
    846             # Avoid using to_value to ensure that we make a copy. We also
    847             # don't want to slow down this method (esp. the scalar case).
--> 848             value = self._to_value(unit, equivalencies)
    849         else:
    850             # to_value only copies if necessary

/opt/miniconda3/lib/python3.8/site-packages/astropy/units/quantity.py in _to_value(self, unit, equivalencies)
    800         if not self.dtype.names or isinstance(self.unit, StructuredUnit):
    801             # Standard path, let unit to do work.
--> 802             return self.unit.to(unit, self.view(np.ndarray),
    803                                 equivalencies=equivalencies)
    804 

/opt/miniconda3/lib/python3.8/site-packages/astropy/units/core.py in to(self, other, value, equivalencies)
   1133             return UNITY
   1134         else:
-> 1135             return self._get_converter(Unit(other),
   1136                                        equivalencies=equivalencies)(value)
   1137 

/opt/miniconda3/lib/python3.8/site-packages/astropy/units/core.py in convert(v)
    988         def make_converter(scale1, func, scale2):
    989             def convert(v):
--> 990                 return func(_condition_arg(v) / scale1) * scale2
    991             return convert
    992 

/opt/miniconda3/lib/python3.8/site-packages/astropy/units/equivalencies.py in iconverter(x)
    205 
    206     def iconverter(x):
--> 207         return x / (wav.to_value(si.AA, spectral()) ** 2 / c_Aps)
    208 
    209     def converter_f_nu_to_nu_f_nu(x):

AttributeError: 'numpy.ndarray' object has no attribute 'to_value'

I think this concerns one of the functions that is added recently ( #33 ). Could you check this?

comparing bytes and strings

The following code does not work,

simqso/simqso/hiforest.py

Line 29 in f066628

if b'HI' in linelist.ION[i]])

it should be replaced by the following

if 'HI' in linelist.ION[i]])

Traceback (most recent call last):
  File "rrarchetypes", line 229, in <module>
    flux, subtype = get_quasars(args.nb,wave,args.seed)
  File "rrarchetypes", line 168, in get_quasars
    data['QSO']['FLUX'], data['QSO']['WAVE'], data['QSO']['META'], data['QSO']['OBJMETA'] = SIMQSO(basewave_max=10e4).make_templates(nmodel=data['QSO']['NB'],
  File "<HOME>/desisim/py/desisim/templates.py", line 2352, in __init__
    from simqso.sqbase import ContinuumKCorr, fixed_R_dispersion
  File "<HOME>/simqso/simqso/__init__.py", line 5, in <module>
    from .sqrun import qsoSimulation
  File "<HOME>/simqso/simqso/sqrun.py", line 12, in <module>
    from . import hiforest
  File "<HOME>/simqso/simqso/hiforest.py", line 37, in <module>
    transitionParams = _getlinelistdata()
  File "<HOME>/simqso/simqso/hiforest.py", line 28, in _getlinelistdata
    Hlines = np.array([i for i in range(linelist.size) 
  File "<HOME>/simqso/simqso/hiforest.py", line 29, in <listcomp>
    if b'HI' in linelist.ION[i]])
TypeError: 'in <string>' requires string as left operand, not bytes

ValueError: low >= high when adding Lyman-alpha forest to BOSS/DR9 QLF

from astropy.cosmology import Planck13
from simqso.sqgrids import *
from simqso import sqbase
from simqso.sqmodels import BOSS_DR9_PLEpivot,get_BossDr9_model_vars
wave = sqbase.fixed_R_dispersion(1000,20e4,1000)
nqso = 10
np.random.seed(12345)
zin = 2.0 + np.random.rand(nqso)
kcorr = sqbase.ContinuumKCorr('DECam-r',1450,effWaveBand='SDSS-r')
qsos = generateQlfPoints(BOSS_DR9_PLEpivot(cosmo=Planck13),
                         (17,22),(2.0,3.0),
                         kcorr=kcorr,zin=zin,
                         qlfseed=12345,gridseed=67890)
sedVars = get_BossDr9_model_vars(qsos,wave,0,noforest=False)
qsos.addVars(sedVars)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-13-e6c2e1a93eb5> in <module>()
     14                          qlfseed=12345,gridseed=67890)
     15 sedVars = get_BossDr9_model_vars(qsos,wave,0,noforest=False)
---> 16 qsos.addVars(sedVars)

~/repos/simqso/simqso/sqgrids.py in addVars(self, newVars, noVals)
    950         '''
    951         for var in newVars:
--> 952             self.addVar(var,noVals=noVals)
    953     def addData(self,data):
    954         self.data = hstack([self.data,data])

~/repos/simqso/simqso/sqgrids.py in addVar(self, var, noVals)
    942         self.varNames.append(var.name)
    943         if not noVals:
--> 944             vals = var(self.nObj)
    945             if vals is not None:
    946                 self.data[var.name] = vals

~/repos/simqso/simqso/sqgrids.py in __call__(self, n, **kwargs)
    366         if self.seed:
    367             np.random.seed(self.seed)
--> 368         vals = self.sampler(n,**kwargs)
    369         if vals is not None:
    370             vals = np.array(vals).astype(self.dtype)

~/repos/simqso/simqso/sqgrids.py in __call__(self, n, **kwargs)
     52         pass
     53     def __call__(self,n,**kwargs):
---> 54         return self.sample(n,**kwargs)
     55     def __str__(self):
     56         s = str((self.low,self.high))

~/repos/simqso/simqso/sqgrids.py in sample(self, n, **kwargs)
    109         super(RandomSubSampler,self).__init__(0,n)
    110     def sample(self,n,**kwargs):
--> 111         return np.random.randint(self.low,self.high,n)
    112 
    113 class ConstSampler(Sampler):

mtrand.pyx in mtrand.RandomState.randint()

ValueError: low >= high

create first release?

@imcgreer thank you for all the recent changes. Would you be willing to create a v1.0 tag and release of simqso so that my development in desihub/desisim#293 can be done against a static version of your code?

add vanden berk template

Include in data table, and also use it to generate an emission line k-correction, which can be used as the default to reduce the number of iterations for flux grids.

write/read metadata to/from existing FITS file and HDU name?

It would be very helpful if the sqgrids.QsoSimObject.write method could optionally write to an existing FITS file and to a specified HDU name. Similarly, it would be great if the corresponding read method could read a particular HDU / extension name.

In case it's useful, in DESI we use this simple wrapper on various astropy convenience methods --
https://github.com/desihub/desispec/blob/master/py/desispec/io/util.py#L147-L209

Thanks!

print only when verbose=True

There are many places in the code that print out helpful messages (e.g., sqrun.buildSpectraBulk). However, it would be great if these messages could be suppressed by setting verbose=False (or similar).

add balmer continuum

have some code in old version...

cleanup handling of forest sightlines

There are two cases for generating sightlines: 1) [the default] create a subsample of forest sightlines and randomly assign them to quasars, e.g., 10 sightlines for 50 quasars means each sightline is used by an average of 5 quasars. This mode is implemented by grouping the quasars by sightline, sorting by redshift, and then generating the spectra in order so that the absorbers get added in increasing redshift. 2) each quasar has its own sightline, in which case this grouping/sorting isn't necessary.

Implementing both of these modes in IGMTransmissionGrid is messy. The second mode could be implemented by using all_spec() from that class, and then passing the resulting Table to CachedIGMTransmissionGrid. But that's kind of messy too.

Unable to fix emission line profiles with VariedEmissionLineGrid

VariedEmissionLineGrid has a 'fixProfiles' option. It is meant to use the emission line grid to derive line profiles but without any scatter in the profile parameters. It is not currently implemented.

self.fixed = kwargs.get('fixLineProfiles',False)

hiforest.py:        self.zmin = wave.min() / 1215.7 - 1.01
sqgrids.py:                                       [1215.7])
sqgrids.py:    >>> v = GaussianEmissionLineVar([GaussianSampler(1215.7,0.1),GaussianSampler(100.,10.),GaussianSampler(10.,1.)])
sqgrids.py:    array([[ 1215.645,   113.125,     9.099],
sqgrids.py:           [ 1215.987,   109.654,     9.312],
sqgrids.py:           [ 1215.74 ,   101.765,    10.822]])

Does that mean that the rest-frame wavelength of Lyman-alpha emission line and the definition of zmin is not using lambda_lya = 1215.67 A? but is using 1215.7?

Thanks for the package, it works very well.