Merge pull request #611 from desihub/ADM-mock-randoms

Updates to faciliate the mocks run on pixel-level imaging data.

doc/changes.rst

@@ -5,8 +5,14 @@ desitarget Change Log
 0.38.1 (unreleased)
 -------------------
 
-* No changes yet.
-
+* Help the mocks run on pixel-level imaging data [`PR #611`_]. Includes:
+    * New :func:`geomask.get_brick_info()` function to look up the
+      brick names associated with each HEALPixel.
+    * In :func:`randoms.quantities_at_positions_in_a_brick()`, add a
+      `justlist` option to list the (maximal) required input files.
+    * Minor bug fixes and documentation updates.
+
+.. _`PR #611`: https://github.com/desihub/desitarget/pull/611
 
 0.38.0 (2020-04-23)
 -------------------

py/desitarget/QA.py

@@ -1998,8 +1998,11 @@ def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.'
         # ADM Write out the densest pixels.
         if makeplots:
             html.write('<pre>')
-            html.write('Densest HEALPixels (NSIDE={}; links are to LS viewer):\n'.format(
-                densdict[objtype]["NSIDE"]))
+            resol = 0
+            if isinstance(densdict[objtype]["NSIDE"], int):
+                resol = hp.nside2resol(densdict[objtype]["NSIDE"], arcmin=True)
+            html.write('Densest HEALPixels (NSIDE={}; ~{:.0f} arcmin across; links are to LS viewer):\n'.format(
+                densdict[objtype]["NSIDE"], resol))
             ras, decs, dens = densdict[objtype]["RA"], densdict[objtype]["DEC"], densdict[objtype]["DENSITY"]
             for i in range(len(ras)):
                 ender = "   "
@@ -2226,7 +2229,8 @@ def make_qa_page(targs, mocks=False, makeplots=True, max_bin_area=1.0, qadir='.'
                             down = 0.9
                         qasystematics_scatterplot(pixmap, sysname, objtype, qadir=qadir,
                                                   downclip=down, upclip=up, nbins=10, xlabel=plotlabel)
-
+                log.info('Made systematics plots for {}...t ={:.1f}s'.format(
+                    sysname, time()-start))
         log.info('Done with systematics...t = {:.1f}s'.format(time()-start))
 
     # ADM html postamble for main page.

py/desitarget/geomask.py

@@ -863,6 +863,66 @@ def add_hp_neighbors(nside, pixnum):
     return pixnum
 
 
+def brick_names_touch_hp(nside, numproc=1):
+    """Determine which of a set of brick names touch a set of HEALPixels.
+
+    Parameters
+    ----------
+    nside : :class:`int`
+        (NESTED) HEALPixel nside.
+    numproc : :class:`int`, optional, defaults to 1
+        The number of parallel processes to use.
+
+    Returns
+    -------
+    :class:`list`
+        A list of lists of input brick names that touch each HEALPixel
+        at `nside`. So, e.g. for `nside=2` the returned list will have
+        48 entries, and, for example, output[0] will be a list of names
+        of bricks that touch HEALPixel 0.
+
+    Notes
+    -----
+        - Runs in about 65 (10) seconds for numproc=1 (32) at nside=2.
+        - Runs in about 325 (20) seconds for numproc=1 (32) at nside=64.
+    """
+    t0 = time()
+    # ADM grab the standard table of bricks.
+    bricktable = brick.Bricks(bricksize=0.25).to_table()
+
+    def _make_lookupdict(indexes):
+        """for a set of indexes that correspond to bricktable rows, make
+        a look-up dictionary of which pixels touch each brick"""
+
+        lookupdict = {bt["BRICKNAME"]: hp_in_box(
+            nside, [bt["RA1"], bt["RA2"], bt["DEC1"], bt["DEC2"]]
+        ) for bt in bricktable[indexes]}
+
+        return lookupdict
+
+    # ADM split the length of the bricktable into arrays of indexes.
+    indexes = np.array_split(np.arange(len(bricktable)), numproc)
+
+    if numproc > 1:
+        pool = sharedmem.MapReduce(np=numproc)
+        with pool:
+            lookupdict = pool.map(_make_lookupdict, indexes)
+        lookupdict = {key: val for lud in lookupdict for key, val in lud.items()}
+    else:
+        lookupdict = _make_lookupdict(indexes[0])
+
+    # ADM change the pixels-in-brick look-up table to a
+    # ADM bricks-in-pixel look-up table.
+    bricksperpixel = [[] for pix in range(hp.nside2npix(nside))]
+    for brickname in lookupdict:
+        for pixel in lookupdict[brickname]:
+            bricksperpixel[pixel].append(brickname)
+
+    log.info("Done...t = {:.1f}s".format(time()-t0))
+
+    return bricksperpixel
+
+
 def sweep_files_touch_hp(nside, pixlist, infiles):
     """Determine which of a set of sweep files touch a set of HEALPixels.
 

py/desitarget/gfa.py

@@ -297,7 +297,7 @@ def all_gaia_in_tiles(maglim=18, numproc=4, allsky=False,
     else:
         # ADM this is critical for, e.g., unit tests for which the
         # ADM Gaia "00000" pixel file might not exist.
-        dummyfile = find_gaia_files_hp(_get_gaia_nside(), pixlist[0],
+        dummyfile = find_gaia_files_hp(nside, pixlist[0],
                                        neighbors=False)[0]
     dummygfas = np.array([], gaia_in_file(dummyfile).dtype)