Skip to content

Commit

Permalink
Merge pull request #37 from simonsobs/v0.4.0b1-fixes
Browse files Browse the repository at this point in the history 
V0.4.0b1 fixes
  • Loading branch information
@mattyowl
mattyowl authored Oct 1, 2020
2 parents ac7b264 + fd92446 commit 032c942
Show file tree
Hide file tree
Showing 9 changed files with 132 additions and 116 deletions.
45 changes: 40 additions & 5 deletions bin/nemoCosmo
View file
Original file line number Diff line number Diff line change
Expand Up @@ -131,7 +131,7 @@ def calcMass(tab, massOptions, tckQFitDict, fRelWeightsDict, mockSurvey, verbose
# Cuts on z, fixed_y_c for forced photometry mode (invalid objects will be listed but without a mass)
if row['fixed_y_c'] > 0 and np.isnan(row['redshift']) == False:
# Corrected for mass function steepness
massDict=signals.calcM500Fromy0(row['fixed_y_c']*1e-4, row['fixed_err_y_c']*1e-4,
massDict=signals.calcMass(row['fixed_y_c']*1e-4, row['fixed_err_y_c']*1e-4,
row['redshift'], row['redshiftErr'],
tenToA0 = massOptions['tenToA0'],
B0 = massOptions['B0'],
Expand Down Expand Up @@ -185,6 +185,11 @@ if __name__ == '__main__':
parser.add_argument("-f", "--footprint", dest = "footprint", help="""Footprint to use, e.g., DES,
HSC, KiDS (default: full). Note that the catalog will not be trimmed to match
the footprint.""", default = None)
parser.add_argument("-R", "--relativistic-correction", dest = "relativisticCorrection",
action = "store_true", default = False, help = """Apply relativistic correction
to the signal modelling and completeness calculation.""")
parser.add_argument("-D", "--mass-def", dest = "massDef", default = "M500c", help = """Mass
definition to use (e.g., M500c or M200m).""")
parser.add_argument("-m", "--max-samples", dest="maxSamples", help="""Maximum number of samples. If
given, overrides the value given in the Cobaya configuration file.""", type = int,
default = defaultMaxSamples)
Expand All @@ -209,6 +214,15 @@ if __name__ == '__main__':
footprintLabel=args.footprint
cosmoOutDir=args.cosmoOutDir
numToBurn=args.numToBurn
massDef=args.massDef
if massDef == "M500c":
rhoType="critical"
delta=500
elif massDef == "M200m":
rhoType="matter"
delta=200
else:
raise Exception("massDef should be either M500c or M200m (use -D M500c or -D M200m)")

if args.likelihoodType == "mass":
lnprob=lnprob_M500
Expand Down Expand Up @@ -266,16 +280,23 @@ if __name__ == '__main__':
info['sampler']['mcmc']['max_samples']=maxSamples

print("Setting up SNR > %.2f selection function (footprint: %s)" % (SNRCut, footprintLabel))
if args.relativisticCorrection == True:
print("Relativistic correction will be applied")
else:
print("Relativistic correction neglected")
selFn=completeness.SelFn(selFnDir, SNRCut, footprintLabel = footprintLabel, zStep = 0.1,
enableDrawSample = True, mockOversampleFactor = 1)
enableDrawSample = True, mockOversampleFactor = 1, downsampleRMS = True,
applyRelativisticCorrection = args.relativisticCorrection,
delta = delta, rhoType = rhoType)

# Cut catalog according to given SNR cut and survey footprint
tab=tab[np.where(tab['fixed_SNR'] > SNRCut)]
inMask=selFn.checkCoordsInAreaMask(tab['RADeg'], tab['decDeg'])
tab=tab[inMask]
print("%d clusters with fixed_SNR > %.1f in footprint %s" % (len(tab), SNRCut, str(footprintLabel)))

# NOTE: Since we use such coarse z binning, we completely ignore photo-z errors here
# Otherwise, signals.calcPM500 throws exception when z grid is not fine enough
# Otherwise, signals.calcPMass throws exception when z grid is not fine enough
tab['redshiftErr'][:]=0.0

# Define binning on (log10(fixed_y_c), redshift) grid
Expand All @@ -300,11 +321,25 @@ if __name__ == '__main__':

H0, Om0, Ob0, sigma8, ns = 68.0, 0.31, 0.049, 0.81, 0.965 # WebSky cosmology
#tenToA0, B0, Mpivot, sigma_int = 2.65e-05, 0.0, 3.0e+14, 0.2 # WebSky scaling relation - assumed scatter
tenToA0, B0, Mpivot, sigma_int = 3.02e-05, 0.0, 3.0e+14, 0.0 # WebSky scaling relation - no scatter
if massDef == 'M500c':
tenToA0, B0, Mpivot, sigma_int = 3.02e-05, 0.0, 3.0e+14, 0.0 # WebSky scaling relation - no scatter, M500c
elif massDef == 'M200m':
tenToA0, B0, Mpivot, sigma_int = 1.7e-05, 0.0, 3.0e+14, 0.0 # WebSky scaling relation - no scatter, M200m
#lnprob(H0, Om0, sigma8, Ob0, ns, p, B0, Mpivot, sigma_int)
testsToRun=['H0', 'Om0', 'sigma8', 'tenToA0', 'B0', 'sigma_int']
#testsToRun=['tenToA0', 'B0', 'sigma_int']
print("Running likelihood tests")
# For checking if recovered parameters beat the theoretical max likelihood ones
# If they do, we have bugs to fix...
theoreticalMaxLogLike=lnprob(H0, Om0, sigma8, Ob0, ns, tenToA0, B0, Mpivot, sigma_int)
print("Theoretical max log likelihood = %.3f" % (theoreticalMaxLogLike))
print("Num clusters expected = %1f" % ((selFn.compMz*selFn.mockSurvey.clusterCount).sum()))
print("Num clusters in catalog = %d" % (len(tab)))
mode=input("Drop into interactive mode [y] or run likelihood tests for each parameter [n] ? ")
if mode == "y":
import IPython
IPython.embed()
sys.exit()
testLogLike=lnprob(70.5, 0.304, 0.825, Ob0, ns, 2.78e-05, 0.306, Mpivot, 0.0) # Any test parameters
if 'H0' in testsToRun:
label="H0"; var=H0
pRange=np.linspace(60, 80, 21)
Expand Down
122 changes: 40 additions & 82 deletions bin/nemoMass
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,9 +97,19 @@ def addForcedPhotometry(pathToCatalog, config, zColumnName = None, zErrColumnNam
return tab

#------------------------------------------------------------------------------------------------------------
def calcMass(tab, massOptions, tckQFitDict, fRelWeightsDict):
def calcMass(tab, massOptions, tckQFitDict, fRelWeightsDict, mockSurvey):
"""Calculates masses for cluster data in table.
"""

label=mockSurvey.mdefLabel
colNames=['%s' % (label), '%sUncorr' % (label)]
if 'rescaleFactor' in massOptions.keys():
colNames.append('%sCal' % (label))
for c in colNames:
tab['%s' % (c)]=np.zeros(len(tab))
tab['%s_errPlus' % (c)]=np.zeros(len(tab))
tab['%s_errMinus' % (c)]=np.zeros(len(tab))

count=0
for row in tab:
Expand All @@ -112,93 +122,50 @@ def calcMass(tab, massOptions, tckQFitDict, fRelWeightsDict):
# Cuts on z, fixed_y_c for forced photometry mode (invalid objects will be listed but without a mass)
if row['fixed_y_c'] > 0 and np.isnan(row['redshift']) == False:
# Corrected for mass function steepness
massDict=signals.calcM500Fromy0(row['fixed_y_c']*1e-4, row['fixed_err_y_c']*1e-4,
massDict=signals.calcMass(row['fixed_y_c']*1e-4, row['fixed_err_y_c']*1e-4,
row['redshift'], row['redshiftErr'],
tenToA0 = massOptions['tenToA0'],
B0 = massOptions['B0'],
Mpivot = massOptions['Mpivot'],
sigma_int = massOptions['sigma_int'],
tckQFit = tckQFitDict[tileName], mockSurvey = mockSurvey,
applyMFDebiasCorrection = True,
applyRelativisticCorrection = True,
applyRelativisticCorrection = massOptions['relativisticCorrection'],
fRelWeightsDict = fRelWeightsDict[tileName])
row['M500']=massDict['M500']
row['M500_errPlus']=massDict['M500_errPlus']
row['M500_errMinus']=massDict['M500_errMinus']
# Without relativistic correction (may be useful)
massDict_noRel=signals.calcM500Fromy0(row['fixed_y_c']*1e-4, row['fixed_err_y_c']*1e-4,
row['redshift'], row['redshiftErr'],
tenToA0 = massOptions['tenToA0'],
B0 = massOptions['B0'],
Mpivot = massOptions['Mpivot'],
sigma_int = massOptions['sigma_int'],
tckQFit = tckQFitDict[tileName], mockSurvey = mockSurvey,
applyMFDebiasCorrection = True,
applyRelativisticCorrection = False,
fRelWeightsDict = fRelWeightsDict[tileName])
row['M500NoRel']=massDict_noRel['M500']
row['M500NoRel_errPlus']=massDict_noRel['M500_errPlus']
row['M500NoRel_errMinus']=massDict_noRel['M500_errMinus']
row['%s' % (label)]=massDict['%s' % (label)]
row['%s_errPlus' % (label)]=massDict['%s_errPlus' % (label)]
row['%s_errMinus' % (label)]=massDict['%s_errMinus' % (label)]
# Uncorrected for mass function steepness
unCorrMassDict=signals.calcM500Fromy0(row['fixed_y_c']*1e-4, row['fixed_err_y_c']*1e-4,
unCorrMassDict=signals.calcMass(row['fixed_y_c']*1e-4, row['fixed_err_y_c']*1e-4,
row['redshift'], row['redshiftErr'],
tenToA0 = massOptions['tenToA0'],
B0 = massOptions['B0'],
Mpivot = massOptions['Mpivot'],
sigma_int = massOptions['sigma_int'],
tckQFit = tckQFitDict[tileName], mockSurvey = mockSurvey,
applyMFDebiasCorrection = False,
applyRelativisticCorrection = True,
applyRelativisticCorrection = massOptions['relativisticCorrection'],
fRelWeightsDict = fRelWeightsDict)
row['M500Uncorr']=unCorrMassDict['M500']
row['M500Uncorr_errPlus']=unCorrMassDict['M500_errPlus']
row['M500Uncorr_errMinus']=unCorrMassDict['M500_errMinus']
# Uncorrected for mass function steepness - without relativistic correction
unCorrMassDict_noRel=signals.calcM500Fromy0(row['fixed_y_c']*1e-4, row['fixed_err_y_c']*1e-4,
row['redshift'], row['redshiftErr'],
tenToA0 = massOptions['tenToA0'],
B0 = massOptions['B0'],
Mpivot = massOptions['Mpivot'],
sigma_int = massOptions['sigma_int'],
tckQFit = tckQFitDict[tileName], mockSurvey = mockSurvey,
applyMFDebiasCorrection = False,
applyRelativisticCorrection = False,
fRelWeightsDict = fRelWeightsDict)
row['M500UncorrNoRel']=unCorrMassDict_noRel['M500']
row['M500UncorrNoRel_errPlus']=unCorrMassDict_noRel['M500_errPlus']
row['M500UncorrNoRel_errMinus']=unCorrMassDict_noRel['M500_errMinus']
row['%sUncorr' % (label)]=unCorrMassDict['%s' % (label)]
row['%sUncorr_errPlus' % (label)]=unCorrMassDict['%s_errPlus' % (label)]
row['%sUncorr_errMinus' % (label)]=unCorrMassDict['%s_errMinus' % (label)]
# Mass conversion
row['M200m']=signals.convertM500cToM200m(massDict['M500']*1e14, row['redshift'])/1e14
row['M200m_errPlus']=(row['M500_errPlus']/row['M500'])*row['M200m']
row['M200m_errMinus']=(row['M500_errMinus']/row['M500'])*row['M200m']
row['M200mUncorr']=signals.convertM500cToM200m(unCorrMassDict['M500']*1e14, row['redshift'])/1e14
row['M200mUncorr_errPlus']=(row['M500Uncorr_errPlus']/row['M500Uncorr'])*row['M200mUncorr']
row['M200mUncorr_errMinus']=(row['M500Uncorr_errMinus']/row['M500Uncorr'])*row['M200mUncorr']
# Mass conversion - no relativistic correction
row['M200mNoRel']=signals.convertM500cToM200m(massDict_noRel['M500']*1e14, row['redshift'])/1e14
row['M200mNoRel_errPlus']=(row['M500NoRel_errPlus']/row['M500NoRel'])*row['M200mNoRel']
row['M200mNoRel_errMinus']=(row['M500NoRel_errMinus']/row['M500NoRel'])*row['M200mNoRel']
row['M200mUncorrNoRel']=signals.convertM500cToM200m(unCorrMassDict_noRel['M500']*1e14, row['redshift'])/1e14
row['M200mUncorrNoRel_errPlus']=(row['M500UncorrNoRel_errPlus']/row['M500UncorrNoRel'])*row['M200mUncorrNoRel']
row['M200mUncorrNoRel_errMinus']=(row['M500UncorrNoRel_errMinus']/row['M500UncorrNoRel'])*row['M200mUncorrNoRel']
#row['M200m']=signals.convertM500cToM200m(massDict['M500']*1e14, row['redshift'])/1e14
#row['M200m_errPlus']=(row['M500_errPlus']/row['M500'])*row['M200m']
#row['M200m_errMinus']=(row['M500_errMinus']/row['M500'])*row['M200m']
#row['M200mUncorr']=signals.convertM500cToM200m(unCorrMassDict['M500']*1e14, row['redshift'])/1e14
#row['M200mUncorr_errPlus']=(row['M500Uncorr_errPlus']/row['M500Uncorr'])*row['M200mUncorr']
#row['M200mUncorr_errMinus']=(row['M500Uncorr_errMinus']/row['M500Uncorr'])*row['M200mUncorr']
# Re-scaling (e.g., using richness-based weak-lensing mass calibration)
row['M500Cal']=massDict['M500']/massOptions['rescaleFactor']
row['M500Cal_errPlus']=np.sqrt(np.power(row['M500_errPlus']/row['M500'], 2) + \
np.power(massOptions['rescaleFactorErr']/massOptions['rescaleFactor'], 2))*row['M500Cal']
row['M500Cal_errMinus']=np.sqrt(np.power(row['M500_errMinus']/row['M500'], 2) + \
np.power(massOptions['rescaleFactorErr']/massOptions['rescaleFactor'], 2))*row['M500Cal']
row['M200mCal']=signals.convertM500cToM200m(row['M500Cal']*1e14, row['redshift'])/1e14
row['M200mCal_errPlus']=(row['M500Cal_errPlus']/row['M500Cal'])*row['M200mCal']
row['M200mCal_errMinus']=(row['M500Cal_errMinus']/row['M500Cal'])*row['M200mCal']
# Re-scaling (e.g., using richness-based weak-lensing mass calibration) - no relativistic correction
row['M500CalNoRel']=massDict_noRel['M500']/massOptions['rescaleFactor']
row['M500CalNoRel_errPlus']=np.sqrt(np.power(row['M500NoRel_errPlus']/row['M500NoRel'], 2) + \
np.power(massOptions['rescaleFactorErr']/massOptions['rescaleFactor'], 2))*row['M500CalNoRel']
row['M500CalNoRel_errMinus']=np.sqrt(np.power(row['M500NoRel_errMinus']/row['M500NoRel'], 2) + \
np.power(massOptions['rescaleFactorErr']/massOptions['rescaleFactor'], 2))*row['M500CalNoRel']
row['M200mCalNoRel']=signals.convertM500cToM200m(row['M500CalNoRel']*1e14, row['redshift'])/1e14
row['M200mCalNoRel_errPlus']=(row['M500CalNoRel_errPlus']/row['M500CalNoRel'])*row['M200mCalNoRel']
row['M200mCalNoRel_errMinus']=(row['M500CalNoRel_errMinus']/row['M500CalNoRel'])*row['M200mCalNoRel']
if 'rescaleFactor' in massOptions.keys():
row['%sCal' % (label)]=massDict['%s' % (label)]/massOptions['rescaleFactor']
row['%sCal_errPlus' % (label)]=np.sqrt(np.power(row['%s_errPlus' % (label)]/row['%s' % (label)], 2) + \
np.power(massOptions['rescaleFactorErr']/massOptions['rescaleFactor'], 2))*row['%sCal' % (label)]
row['%sCal_errMinus' % (label)]=np.sqrt(np.power(row['%s_errMinus' % (label)]/row['%s' % (label)], 2) + \
np.power(massOptions['rescaleFactorErr']/massOptions['rescaleFactor'], 2))*row['%sCal' % (label)]
#row['M200mCal']=signals.convertM500cToM200m(row['M500Cal']*1e14, row['redshift'])/1e14
#row['M200mCal_errPlus']=(row['M500Cal_errPlus']/row['M500Cal'])*row['M200mCal']
#row['M200mCal_errMinus']=(row['M500Cal_errMinus']/row['M500Cal'])*row['M200mCal']

return tab

Expand Down Expand Up @@ -260,15 +227,15 @@ if __name__ == '__main__':
nemoTab['redshiftErr']=zTab['redshiftErr']
tab=nemoTab
if outFileName is None:
outFileName=optimalCatalogFileName.replace("_optimalCatalog.fits", "_M500.fits")
outFileName=optimalCatalogFileName.replace("_optimalCatalog.fits", "_mass.fits")

else:

# Load another catalog (e.g., a mock, for testing)
optimalCatalogFileName=catFileName
tab=atpy.Table().read(optimalCatalogFileName)
if outFileName is None:
outFileName=catFileName.replace(".fits", "_M500.fits")
outFileName=catFileName.replace(".fits", "_mass.fits")

# Enter forced photometry mode if we can't find the columns we need
# If we're doing forced photometry, we're going to need to set-up so we can find the filtered maps
Expand Down Expand Up @@ -324,16 +291,7 @@ if __name__ == '__main__':
else:
ns=0.95
mockSurvey=MockSurvey.MockSurvey(minMass, areaDeg2, zMin, zMax, H0, Om0, Ob0, sigma8, ns)

colsToAdd=['M500', 'M500Uncorr', 'M200m', 'M200mUncorr',
'M500NoRel', 'M500UncorrNoRel', 'M200mNoRel', 'M200mUncorrNoRel']
if 'rescaleFactor' in list(massOptions.keys()):
colsToAdd=colsToAdd+['M500Cal', 'M200mCal', 'M500CalNoRel', 'M200mCalNoRel']
for c in colsToAdd:
tab.add_column(atpy.Column(np.zeros(len(tab)), c))
tab.add_column(atpy.Column(np.zeros(len(tab)), c+"_errPlus"))
tab.add_column(atpy.Column(np.zeros(len(tab)), c+"_errMinus"))


# Seems like multiprocessing doesn't work under slurm...
# So use MPI instead: just divide up catalog among processes
tab.add_column(atpy.Column(np.arange(len(tab)), "sortIndex"))
Expand All @@ -345,7 +303,7 @@ if __name__ == '__main__':
endIndex=len(tab)
tab=tab[startIndex:endIndex]

tab=calcMass(tab, massOptions, tckQFitDict, fRelWeightsDict)
tab=calcMass(tab, massOptions, tckQFitDict, fRelWeightsDict, mockSurvey)

if config.MPIEnabled == True:
tabList=config.comm.gather(tab, root = 0)
Expand Down
Loading

0 comments on commit 032c942

Please sign in to comment.