Merge pull request #64 from simonsobs/dev

Dev

bin/nemoMass

@@ -89,7 +89,7 @@ def addForcedPhotometry(pathToCatalog, config, zColumnName = None, zErrColumnNam
         if mapFilter['label'] == config.parDict['photFilter']:
             trimmedList.append(mapFilter)
     config.parDict['mapFilters']=trimmedList
-    forcedTab=pipelines.filterMapsAndMakeCatalogs(config)
+    forcedTab=pipelines.filterMapsAndMakeCatalogs(config, useCachedFilteredMaps = True)
 
     # Need to graft the redshifts back on
     zMatched, forcedMatched, rDeg=catalogs.crossMatch(zTab, forcedTab)
@@ -322,7 +322,7 @@ if __name__ == '__main__':
             if key not in tab.keys():
                 forcedPhotometry=True
         config=startUp.NemoConfig(parDictFileName, MPIEnabled = args.MPIEnabled, divideTilesByProcesses = False,
-                                  setUpMaps = forcedPhotometry, writeTileDir = False, verbose = False,
+                                  setUpMaps = forcedPhotometry, writeTileInfo = False, verbose = False,
                                   strictMPIExceptions = strictMPIExceptions)
 
     # Remaining set-up

nemo/filters.py

@@ -77,7 +77,7 @@ def filterMaps(unfilteredMapsDictList, filterParams, tileName, diagnosticsDir =
 
     f=filterParams
     label=f['label']+"#"+tileName
-        
+
     print("... making filtered map %s" % (label))
     filterClass=eval('%s' % (f['class']))
     filterObj=filterClass(f['label'], unfilteredMapsDictList, f['params'], tileName = tileName,

nemo/maps.py

@@ -26,7 +26,6 @@
 import yaml
 import pickle
 from pixell import enmap, curvedsky, utils, powspec
-import sharp
 import nemo
 try:
     import reproject
@@ -1323,7 +1322,7 @@ def simNoiseMap(shape, noiseLevel, wcs = None, lKnee = None, alpha = -3, noiseMo
         def _mod_noise_map(ivar, Nl):
             map1 = enmap.rand_gauss(ivar.shape, ivar.wcs)
             lmax = len(Nl)-1
-            ainfo = sharp.alm_info(lmax)
+            ainfo = curvedsky.alm_info(lmax)
             alm = curvedsky.map2alm(map1, ainfo=ainfo)
             map2 = curvedsky.alm2map(alm, np.zeros_like(map1))
             map1 -= map2

nemo/photometry.py

@@ -351,7 +351,7 @@ def measureFluxes(catalog, filteredMapDict, diagnosticsDir, photFilteredMapDict
                     obj[prefix+"err_fluxJy"]=deltaTToJyPerSr(obj[prefix+'err_deltaT_c'], obsFreqGHz)*beamSolidAngle_nsr*1.e-9
 
 #------------------------------------------------------------------------------------------------------------
-def makeForcedPhotometryCatalog(filteredMapDict, inputCatalog, useInterpolator = True,
+def makeForcedPhotometryCatalog(filteredMapDict, inputCatalog, useInterpolator = True,\
                                 DS9RegionsPath = None):
     """Make a catalog on which we can do forced photometry at the locations of objects in it.
     

nemo/pipelines.py

@@ -30,8 +30,8 @@
 
 #------------------------------------------------------------------------------------------------------------
 def filterMapsAndMakeCatalogs(config, rootOutDir = None, useCachedFilters = False, useCachedRMSMap = False,\
-                              measureFluxes = True, invertMap = False, verbose = True, writeAreaMask = False,\
-                              writeFlagMask = False):
+                              useCachedFilteredMaps = False, measureFluxes = True, invertMap = False, \
+                              verbose = True, writeAreaMask = False, writeFlagMask = False):
     """Runs the map filtering and catalog construction steps according to the given configuration.
     
     Args:
@@ -42,6 +42,9 @@ def filterMapsAndMakeCatalogs(config, rootOutDir = None, useCachedFilters = Fals
             read from disk, rather than re-calculated (used by source injection simulations).
         useCachedRMSMap (:obj:`bool`, optional): If True, use the previously estimated noise map, which has
             been saved to disk. This is only useful for source injection simulations.
+        useCachedFilteredMaps (:obj:`bool`, optional): If True, and previously written maps are found, these
+            will be read from disk, rather than re-made. This is useful for performing forced photometry using
+            external catalogs without having to run nemo again. Otherwise, this should be set to False.
         measureFluxes (bool, optional): If True, measure fluxes. If False, just extract S/N values for 
             detected objects.
         invertMap (bool, optional): If True, multiply all maps by -1; needed by 
@@ -69,6 +72,7 @@ def filterMapsAndMakeCatalogs(config, rootOutDir = None, useCachedFilters = Fals
                 writeFlagMask=True
             config.filterSetOptions[setNum]['catalog']=_filterMapsAndMakeCatalogs(config, verbose = True,
                                                                                   useCachedFilters = False,
+                                                                                  useCachedFilteredMaps = False,
                                                                                   writeAreaMask = writeAreaMask,
                                                                                   writeFlagMask = writeFlagMask)
 
@@ -93,6 +97,7 @@ def filterMapsAndMakeCatalogs(config, rootOutDir = None, useCachedFilters = Fals
     else:
         # Default single pass behaviour (also used by source injection tests)
         catalog=_filterMapsAndMakeCatalogs(config, rootOutDir = rootOutDir, useCachedFilters = useCachedFilters,
+                                           useCachedFilteredMaps = useCachedFilteredMaps,
                                            useCachedRMSMap = useCachedRMSMap, measureFluxes = measureFluxes,
                                            invertMap = invertMap, verbose = verbose,
                                            writeAreaMask = writeAreaMask, writeFlagMask = writeFlagMask)
@@ -104,8 +109,8 @@ def filterMapsAndMakeCatalogs(config, rootOutDir = None, useCachedFilters = Fals
 
 #------------------------------------------------------------------------------------------------------------
 def _filterMapsAndMakeCatalogs(config, rootOutDir = None, useCachedFilters = False, useCachedRMSMap = False,\
-                               measureFluxes = True, invertMap = False, verbose = True, writeAreaMask = False,\
-                               writeFlagMask = False):
+                               useCachedFilteredMaps = False, measureFluxes = True, invertMap = False, \
+                               verbose = True, writeAreaMask = False, writeFlagMask = False):
     """Runs the map filtering and catalog construction steps according to the given configuration.
     
     Args:
@@ -114,6 +119,9 @@ def _filterMapsAndMakeCatalogs(config, rootOutDir = None, useCachedFilters = Fal
             output filtered maps and catalogs are written.
         useCachedFilters (:obj:`bool`, optional): If True, and previously made filters are found, they will be
             read from disk, rather than re-calculated (used by source injection simulations).
+        useCachedFilteredMaps (:obj:`bool`, optional): If True, and previously written maps are found, these
+            will be read from disk, rather than re-made. This is useful for performing forced photometry using
+            external catalogs without having to run nemo again. Otherwise, this should be set to False.
         measureFluxes (bool, optional): If True, measure fluxes. If False, just extract S/N values for 
             detected objects.
         invertMap (bool, optional): If True, multiply all maps by -1; needed by 
@@ -181,10 +189,28 @@ def _filterMapsAndMakeCatalogs(config, rootOutDir = None, useCachedFilters = Fal
             else:
                 DS9RegionsPath=None
 
-            filteredMapDict=filters.filterMaps(config.unfilteredMapsDictList, f, tileName, 
-                                               diagnosticsDir = config.diagnosticsDir, selFnDir = config.selFnDir,
-                                               verbose = True, undoPixelWindow = undoPixelWindow,
-                                               useCachedFilter = useCachedFilters)
+            filteredMapFileName=filteredMapsDir+os.path.sep+tileName+os.path.sep+"%s_filteredMap.fits"  % (label)
+            SNMapFileName=filteredMapsDir+os.path.sep+tileName+os.path.sep+"%s_SNMap.fits" % (label)
+            if useCachedFilteredMaps == True and os.path.exists(filteredMapFileName):
+                print("... loading cached filtered map %s ..." % (filteredMapFileName))
+                filteredMapDict={}
+                with pyfits.open(filteredMapFileName) as img:
+                    filteredMapDict['data']=img[0].data
+                    filteredMapDict['wcs']=astWCS.WCS(img[0].header, mode = 'pyfits')
+                    filteredMapDict['mapUnits']=filteredMapDict['wcs'].header['BUNIT']
+                    if 'BEAMNSR' in filteredMapDict['wcs'].header.keys():
+                        filteredMapDict['beamSolidAngle_nsr']=filteredMapDict['wcs'].header['BEAMNSR']
+                        filteredMapDict['obsFreqGHz']=filteredMapDict['wcs'].header['FREQGHZ']
+                with pyfits.open(SNMapFileName) as img:
+                    filteredMapDict['SNMap']=img[0].data
+                filteredMapDict['surveyMask'], wcs=completeness.loadAreaMask(tileName, config.selFnDir)
+                filteredMapDict['label']=f['label']
+                filteredMapDict['tileName']=tileName
+            else:
+                filteredMapDict=filters.filterMaps(config.unfilteredMapsDictList, f, tileName,
+                                                   diagnosticsDir = config.diagnosticsDir, selFnDir = config.selFnDir,
+                                                   verbose = True, undoPixelWindow = undoPixelWindow,
+                                                   useCachedFilter = useCachedFilters)
 
             if useCachedRMSMap == True and photFilter is not None: # i.e., only an option for cluster insertion sims
                 # This is messy:
@@ -205,8 +231,6 @@ def _filterMapsAndMakeCatalogs(config, rootOutDir = None, useCachedFilters = Fal
                 filteredMapDict['data'][mask]=0 # just in case we rely elsewhere on zero == no data
 
             if 'saveFilteredMaps' in f['params'] and f['params']['saveFilteredMaps'] == True:
-                filteredMapFileName=filteredMapsDir+os.path.sep+tileName+os.path.sep+"%s_filteredMap.fits"  % (label)
-                SNMapFileName=filteredMapsDir+os.path.sep+tileName+os.path.sep+"%s_SNMap.fits" % (label)
                 maps.saveFITS(filteredMapFileName, filteredMapDict['data'], filteredMapDict['wcs'])
                 maps.saveFITS(SNMapFileName, filteredMapDict['SNMap'], filteredMapDict['wcs'])
 

tests/configs/quickstart-clusters-Q.yml

@@ -0,0 +1,121 @@
+# Nemo config file
+# YAML format 
+# - use null to return None in Python
+# - note that YAML is fussy about large numbers: use e.g. 1.0e+14 for M500MSun (not 1e14)
+
+unfilteredMaps:
+    - {mapFileName: "maps/f150_1_10_8_map.fits",
+       weightsFileName: "maps/f150_1_10_8_ivar.fits",
+       obsFreqGHz: 149.6, units: 'uK',
+       beamFileName: "maps/s16_pa2_f150_nohwp_night_beam_profile_jitter.txt"}
+    - {mapFileName: "maps/f090_1_10_8_map.fits",
+       weightsFileName: "maps/f090_1_10_8_ivar.fits",
+       obsFreqGHz: 97.8, units: 'uK',
+       beamFileName: "maps/s16_pa3_f090_nohwp_night_beam_profile_jitter.txt"}
+
+# Masks
+#surveyMask: "maps/Jun2020/AdvACTSurveyMask_v7_galLatCut_S18-dust-artifacts-extended-post10mJy.fits"
+
+# Instead of giving point source mask, can give catalog instead
+maskPointSourcesFromCatalog: 
+    - "PS_S18_f150_auto_rArcmin.fits"
+
+# Detection/catalog options
+thresholdSigma: 4.0
+minObjPix: 1
+findCenterOfMass: True
+useInterpolator: True
+rejectBorder: 0
+objIdent: 'ACT-CL'
+longNames: False
+removeRings: False
+
+# Photometry options
+photFilter: 'Arnaud_M2e14_z0p4'
+fitQ: True
+
+# Optionally override the GNFW parameters - if not present, Arnaud et al. (2010) parameters are used
+# The example below is for the Planck Pressure Profile (PPP)
+#GNFWParams: {P0: 6.41, c500: 1.81, gamma: 0.31, alpha: 1.33, beta: 4.13, tol: 1e-7, npts: 100}
+
+# Mass measurement options - used by nemoMass and nemoSelFn scripts
+# Writes out .fits file to nemoOutDir/nemoOutDir_M500.fits
+# redshiftCatalog:  A .fits table containing name, RADeg, decDeg, redshift, redshiftErr columns
+# tenToA0, B0, Mpivot, sigma_int: Fixed scaling relation options (see H13 or ACTPol paper)
+# rescaleFactor, rescaleFactorErr: For MCal masses, as in the ACTPol paper (i.e., just rescales M500 results by 1/rescaleFactor)
+massOptions: {tenToA0: 4.95e-5, 
+              B0: 0.08, 
+              Mpivot: 3.0e+14, 
+              sigma_int: 0.2,
+              relativisticCorrection: True,
+              rescaleFactor: 0.69, 
+              rescaleFactorErr: 0.07,
+              redshiftCatalog: "DR5_cluster-catalog_v1.1.fits"}
+
+# Selection function options - this calculation can also be enabled with the nemo -S switch
+# NOTE: could eventually add 'completenessFraction' to 'massLimitMaps', which is why that's a dictionary list
+# Use selFnFootprints to calculate average completeness in given sky areas - e.g., overlap with optical surveys
+calcSelFn: False
+selFnOptions: {fixedSNRCut: 5.0, 
+               massLimitMaps: [{z: 0.5}]}
+
+# Filter definitions:
+# allFilters is a dictionary of parameters that will be copied into all mapFilters
+# (these can be overridden by keys with the same name in mapFilters)
+allFilters: {class: "ArnaudModelMatchedFilter",
+             params: {noiseParams: {method: "dataMap",
+                                    noiseGridArcmin: 40.},
+                      saveFilteredMaps: True,
+                      saveRMSMap: True,
+                      savePlots: False,
+                      saveDS9Regions: False,
+                      outputUnits: 'yc',
+                      edgeTrimArcmin: 0.0}
+            }
+
+# mapFilters is a list of all the different filters to apply
+# (keys in mapFilters with the same name as those in allFilters take priority)
+mapFilters:
+    #- {label: "Arnaud_M1e14_z0p2",
+       #params: {M500MSun: 1.0e+14, z: 0.2}}
+    #- {label: "Arnaud_M2e14_z0p2",
+       #params: {M500MSun: 2.0e+14, z: 0.2}}
+    #- {label: "Arnaud_M4e14_z0p2",
+       #params: {M500MSun: 4.0e+14, z: 0.2}}
+    #- {label: "Arnaud_M8e14_z0p2",
+       #params: {M500MSun: 8.0e+14, z: 0.2}}
+    #- {label: "Arnaud_M1e14_z0p4",
+       #params: {M500MSun: 1.0e+14, z: 0.4}}
+    - {label: "Arnaud_M2e14_z0p4",
+       params: {M500MSun: 2.0e+14, z: 0.4,
+                saveFilteredMaps: True,
+                savePlots: True}}
+    #- {label: "Arnaud_M4e14_z0p4",
+       #params: {M500MSun: 4.0e+14, z: 0.4}}
+    #- {label: "Arnaud_M8e14_z0p4",
+       #params: {M500MSun: 8.0e+14, z: 0.4}}              
+    #- {label: "Arnaud_M1e14_z0p8",
+       #params: {M500MSun: 1.0e+14, z: 0.8}}
+    #- {label: "Arnaud_M2e14_z0p8",
+       #params: {M500MSun: 2.0e+14, z: 0.8}}
+    #- {label: "Arnaud_M4e14_z0p8",
+       #params: {M500MSun: 4.0e+14, z: 0.8}}
+    #- {label: "Arnaud_M8e14_z0p8",
+       #params: {M500MSun: 8.0e+14, z: 0.8}}
+    #- {label: "Arnaud_M1e14_z1p2",
+       #params: {M500MSun: 1.0e+14, z: 1.2}}
+    #- {label: "Arnaud_M2e14_z1p2",
+       #params: {M500MSun: 2.0e+14, z: 1.2}}
+    #- {label: "Arnaud_M4e14_z1p2",
+       #params: {M500MSun: 4.0e+14, z: 1.2}}
+    #- {label: "Arnaud_M8e14_z1p2",
+       #params: {M500MSun: 8.0e+14, z: 1.2}}
+
+# Source injection test - this can also be enabled with the nemo -I switch
+sourceInjectionIterations: 10
+sourcesPerTile: 50
+sourceInjectionModels:
+    - {redshift: 0.8, M500: 2.0e+14}
+    - {redshift: 0.4, M500: 2.0e+14}
+    - {redshift: 0.2, M500: 2.0e+14}
+    - {redshift: 0.1, M500: 2.0e+14}

tests/lib/NemoTests.py

@@ -120,7 +120,7 @@ def setup_quickstart(self):
         if os.path.exists(self.inMapFileName) == False:
             print(">>> Downloading quickstart files ...")
             os.chdir(self.cacheDir)
-            os.system("wget https://astro.ukzn.ac.za/~mjh/nemo-quickstart-maps.tar.gz")
+            os.system("wget 'https://www.dropbox.com/scl/fi/5sjapfshk8g3kxgy4zc56/nemo-quickstart-maps.tar.gz?rlkey=swpdttpgvkffbgr9pcmebslw2&dl=0' -O nemo-quickstart-maps.tar.gz")
             os.system("tar -zxvf nemo-quickstart-maps.tar.gz")
             os.remove("nemo-quickstart-maps.tar.gz")
             os.chdir(thisDir)
@@ -163,8 +163,11 @@ def set_config(self, configFileName):
         self.mocksDir=configFileName.replace(".yml", "")+os.path.sep+"mocks"
 
 
-    def run_nemo(self):
-        self._run_command(["nemo", self.configFileName])
+    def run_nemo(self, forcedPhotometryCatalog = None):
+        args=['nemo', self.configFileName]
+        if forcedPhotometryCatalog is not None:
+            args=args+['-f', forcedPhotometryCatalog]
+        self._run_command(args)
 
 
     def run_nemo_injection_test(self):
@@ -175,10 +178,12 @@ def run_parallel_nemo(self):
         self._run_command(["mpiexec", "-np", "4", "nemo", self.configFileName, "-M"])
 
 
-    def run_nemo_mass(self, catalogFileName = None):
+    def run_nemo_mass(self, catalogFileName = None, forcedPhotometry = False):
         args=['nemoMass', self.configFileName]
-        if catalogFileName != None:
+        if catalogFileName is not None:
             args=args+['-c', catalogFileName]
+        if forcedPhotometry == True:
+            args=args+['-F']
         self._run_command(args)
 
 

tests/quick.robot

@@ -9,6 +9,18 @@ Quickstart clusters tutorial runs
 Quickstart sources tutorial runs
     Run quickstart sources
 
+Mass estimation works
+    Run quickstart clusters with Q
+    Run nemo mass
+
+Forced photometry using nemo works
+    Set config              configs/quickstart-clusters-Q.yml
+    Run nemo    DR5_cluster-catalog_v1.1.fits
+
+Forced photometry using nemoMass works
+    Run quickstart clusters with Q
+    Run nemo mass   DR5_cluster-catalog_v1.1.fits   True
+
 Source injection test runs
     Generate simulated source maps
     Set config      configs/sim_ptsrc_f090.yml
@@ -57,6 +69,10 @@ Run quickstart clusters
     Set config              ../examples/quickstart/quickstart-clusters.yml
     Run nemo
 
+Run quickstart clusters with Q
+    Set config              configs/quickstart-clusters-Q.yml
+    Run nemo
+
 Run quickstart sources
     Set config              ../examples/quickstart/quickstart-sources.yml
     Run nemo
@@ -115,7 +131,9 @@ Clean up
     Remove directory        testsCache/quickstart-clusters              True
     Remove directory        testsCache/quickstart-sources               True
     Remove directory        testsCache/quickstart-multipass             True
-
+    Remove directory        testsCache/quickstart-multipass-Q           True
+    Remove file             testsCache/*forcedCatalog*
+    Remove file             testsCache/*_mass.fits
 
 *** Settings ***
 #Documentation              To be added here