diff --git a/simqso/lumfun.py b/simqso/lumfun.py index baf7702..75099bc 100644 --- a/simqso/lumfun.py +++ b/simqso/lumfun.py @@ -217,10 +217,11 @@ def logPhi(self,M,z,par=None): np.log10(10**(0.4*(alpha+1)*(M-Mstar)) + \ 10**(0.4*( beta+1)*(M-Mstar))) def _sample(self,Mrange,zrange,p,**kwargs): + zin = kwargs.pop('zin',None) + verbose = kwargs.pop('verbose',0) eps_M,eps_z = 5e-2,2e-2 nM = int(-np.diff(Mrange(zrange)) / eps_M) nz = int(np.diff(zrange) / eps_z) - zin = kwargs.get('zin') if zin is None: # integrate across redshift to get the dN/dz distribution skyfrac = kwargs.get('skyArea',skyDeg2) / skyDeg2 @@ -234,7 +235,8 @@ def _sample(self,Mrange,zrange,p,**kwargs): Ntot = np.sum(zsamp) zfun = interp1d(np.cumsum(zsamp)/Ntot,zbins) Ntot = np.int(np.round(Ntot * skyfrac * 4*np.pi)) - print('integration returned ',Ntot,' objects') + if verbose > 0: + print('integration returned ',Ntot,' objects') else: # redshifts supplied by user zfun = lambda x: zin @@ -252,11 +254,13 @@ def _sample(self,Mrange,zrange,p,**kwargs): N_M = np.sum(Msamp) Mfun = interp1d(np.cumsum(Msamp)/N_M,Mbins) M[i] = Mfun(y[i]) - if Ntot > 1e4 and ((i+1)%(Ntot//10))==0: + if verbose > 1 and Ntot > 1e4 and ((i+1)%(Ntot//10))==0: print(i+1,' out of ',Ntot) return M,z def _fast_sample(self,Mrange,zrange,p,**kwargs): - print('using fast sample') + verbose = kwargs.pop('verbose',0) + if verbose > 1: + print('using fast sample for QLF') skyfrac = kwargs.get('skyArea',skyDeg2) / skyDeg2 eps_M,eps_z = 0.05,0.10 magLimPad = 0.2 @@ -284,10 +288,11 @@ def _fast_sample(self,Mrange,zrange,p,**kwargs): M,z = np.hstack(Mz) M += dM * (np.random.rand(len(M)) - 0.5) z += dz * (np.random.rand(len(M)) - 0.5) - print('to generate {} quasars'.format(len(M))) + if verbose > 1: + print('to generate {} quasars'.format(len(M))) return M,z def sample_from_fluxrange(self,mrange,zrange,kcorr,p=None,**kwargs): - fast = kwargs.get('fast_sample',False) + fast = kwargs.pop('fast_sample',False) m2M = AppToAbsMag(self.cosmo,kcorr) _mrange = np.array(mrange[::-1]) _Mrange = lambda z: _mrange - m2M(_mrange,z) diff --git a/simqso/sqgrids.py b/simqso/sqgrids.py index 4d11a06..7dd4ad7 100644 --- a/simqso/sqgrids.py +++ b/simqso/sqgrids.py @@ -1251,7 +1251,7 @@ def generateBEffEmissionLines(M1450,**kwargs): lineCatalog['name'][useLines]) return lines -def generateVdBCompositeEmLines(minEW=1.0,noFe=False): +def generateVdBCompositeEmLines(minEW=1.0,noFe=False,verbose=0): tmplfits = os.path.join(datadir,'simqso_templates.fits') all_lines = Table(fits.getdata(tmplfits,'VdB01CompEmLines')) # blended lines are repeated in the table @@ -1263,8 +1263,9 @@ def generateVdBCompositeEmLines(minEW=1.0,noFe=False): if noFe: isFe = lines['ID'].find('Fe') == 0 lines = lines[~isFe] - print('using the following lines from VdB template: ', end=' ') - print(','.join(list(lines['ID']))) + if verbose > 0: + print('using the following lines from VdB template: ', end=' ') + print(','.join(list(lines['ID']))) c = ConstSampler lineList = [ [c(l['OWave']),c(l['EqWid']),c(l['Width'])] for l in lines ] lines = GaussianEmissionLinesTemplateVar(lineList) diff --git a/simqso/sqrun.py b/simqso/sqrun.py index fed6351..258d7fd 100644 --- a/simqso/sqrun.py +++ b/simqso/sqrun.py @@ -141,17 +141,15 @@ def buildForest(wave,z,simParams,outputDir): return tgrid -def buildContinuumModels(qsoGrid,simParams): +def buildContinuumModels(qsoGrid,simParams,verbose=0): continuumParams = simParams['QuasarModelParams']['ContinuumParams'] reseed(continuumParams) slopes = continuumParams['PowerLawSlopes'][::2] breakpts = continuumParams['PowerLawSlopes'][1::2] - print('... building continuum grid') + if verbose > 0: + print('... building continuum grid') cmodel = continuumParams['ContinuumModel'] - if cmodel in ['GaussianPLawDistribution','FixedPLawDistribution', - 'BrokenPowerLaw']: - if cmodel in ['GaussianPLawDistribution','FixedPLawDistribution']: - print('WARNING: %s continuum is deprecated' % cmodel) + if cmodel == 'BrokenPowerLaw': slopeVars = [ grids.GaussianSampler(*s) for s in slopes ] continuumVars = [ grids.BrokenPowerLawContinuumVar(slopeVars, breakpts) ] @@ -179,8 +177,9 @@ def buildEmissionLineGrid(qsoGrid,simParams): emLineParams['EmissionLineModel']) qsoGrid.addVar(emLineGrid) -def buildDustGrid(qsoGrid,simParams): - print('... building dust extinction grid') +def buildDustGrid(qsoGrid,simParams,verbose=0): + if verbose > 0: + print('... building dust extinction grid') dustParams = simParams['QuasarModelParams']['DustExtinctionParams'] reseed(dustParams) if dustParams['DustExtinctionModel'] == 'Fixed E(B-V)': @@ -202,8 +201,8 @@ def buildDustGrid(qsoGrid,simParams): qsoGrid.addVar(dustVar) -def buildFeatures(qsoGrid,wave,simParams,forest=None): - buildContinuumModels(qsoGrid,simParams) +def buildFeatures(qsoGrid,wave,simParams,forest=None,verbose=0): + buildContinuumModels(qsoGrid,simParams,verbose=verbose) qsoParams = simParams['QuasarModelParams'] if 'EmissionLineParams' in qsoParams: buildEmissionLineGrid(qsoGrid,simParams) @@ -213,7 +212,7 @@ def buildFeatures(qsoGrid,wave,simParams,forest=None): feGrid = grids.VW01FeTemplateGrid(qsoGrid.z,wave,scales=scalings) qsoGrid.addVar(grids.FeTemplateVar(feGrid)) if 'DustExtinctionParams' in qsoParams: - buildDustGrid(qsoGrid,simParams) + buildDustGrid(qsoGrid,simParams,verbose=verbose) if forest is not None: if isinstance(forest,hiforest.CachedIGMTransmissionGrid): losMap = forest.losMap @@ -342,9 +341,10 @@ def buildSpectraBySightLine(wave,qsoGrid,procMap=map, Input wavelength grid. ''' photoCache = qsoGrid.getPhotoCache(wave) - print('simulating ',qsoGrid.nObj,' quasar spectra') - print('units are ',qsoGrid.units) - print('max number iterations: ',maxIter) + if verbose > 0: + print('simulating ',qsoGrid.nObj,' quasar spectra') + print('units are ',qsoGrid.units) + print('max number iterations: ',maxIter) verby = 0 if not verbose else qsoGrid.nObj//(5*verbose) if qsoGrid.units == 'luminosity' or photoCache is None: nIter = 1 @@ -404,8 +404,9 @@ def buildSpectraBulk(wave,qsoGrid,procMap=map, Input wavelength grid. ''' photoCache = qsoGrid.getPhotoCache(wave) - print('simulating ',qsoGrid.nObj,' quasar spectra') - print('units are ',qsoGrid.units) + if verbose > 0: + print('simulating ',qsoGrid.nObj,' quasar spectra') + print('units are ',qsoGrid.units) if qsoGrid.units == 'luminosity' or photoCache is None: nIter = 1 fluxBand = None @@ -424,7 +425,8 @@ def buildSpectraBulk(wave,qsoGrid,procMap=map, build_one_spec = partial(buildSpecWithPhot,wave,qsoGrid.cosmo, specFeatures,photoCache,iterNum=iterNum, saveSpectra=saveSpectra) - print('buildSpectra iteration ',iterNum,' out of ',nIter) + if verbose > 1: + print('buildSpectra iteration ',iterNum,' out of ',nIter) specOut = list(procMap(build_one_spec,qsoGrid)) specOut = _regroup(specOut) synMag,synFlux,spectra = specOut @@ -439,7 +441,8 @@ def buildSpectraBulk(wave,qsoGrid,procMap=map, if nIter > 1: # find the largest mag offset dm = synMag[:,fluxBand] - qsoGrid.appMag - print('--> delta mag mean = %.7f, rms = %.7f, |max| = %.7f' % \ + if verbose > 1: + print('--> delta mag mean = %.7f, rms = %.7f, |max| = %.7f' % \ (dm.mean(),dm.std(),np.abs(dm).max())) qsoGrid.absMag[:] -= dm dmagMax = np.abs(dm).max() @@ -526,19 +529,23 @@ def qsoSimulation(simParams,**kwargs): outputDir,retParams=True, clean=True) except IOError: - print(simParams['FileName']+' output not found') + if verbose > 0: + print(simParams['FileName']+' output not found') if 'GridFileName' in simParams: - print('restoring grid from ',simParams['GridFileName']) + if verbose > 0: + print('restoring grid from ',simParams['GridFileName']) try: qsoGrid = readSimulationData(simParams['GridFileName'], outputDir) except IOError: - print(simParams['GridFileName'],' not found, generating') + if verbose > 0: + print(simParams['GridFileName'],' not found, generating') qsoGrid = buildQsoGrid(simParams) qsoGrid.write(simParams,outputDir, simParams['GridFileName']+'.fits') else: - print('generating QSO grid') + if verbose > 0: + print('generating QSO grid') qsoGrid = buildQsoGrid(simParams) if not forestOnly: if not noWriteOutput and 'GridFileName' in simParams: @@ -585,7 +592,7 @@ def qsoSimulation(simParams,**kwargs): # # add the quasar model variables to the grid (does the random sampling) # - buildFeatures(qsoGrid,wave,simParams,forest) + buildFeatures(qsoGrid,wave,simParams,forest,verbose=verbose) timerLog('Generate Features') # # Use continuum and emission line distributions to build the components @@ -599,7 +606,8 @@ def qsoSimulation(simParams,**kwargs): # map the simulated photometry to observed values with uncertainties # if not noPhotoMap: - print('mapping photometry') + if verbose > 0: + print('mapping photometry') reseed(simParams['PhotoMapParams']) photoData = sqphoto.calcObsPhot(qsoGrid.synFlux,qsoGrid.photoMap) qsoGrid.addData(photoData)