Merge ff32c95 into 4aecf75

bin.src/imsim.py

@@ -2,11 +2,13 @@
 """
 This is the imSim program, used to drive GalSim to simulate the LSST.  Written
 for the DESC collaboration and LSST project.  This version of the program can
-read phoSim instance files as is. It leverages the LSST Sims GalSim interface
+read phoSim instance files as is.  It leverages the LSST Sims GalSim interface
 code found in sims_GalSimInterface.
 """
 import os
 import argparse
+import warnings
+from astropy._erfa import ErfaWarning
 import desc.imsim
 
 parser = argparse.ArgumentParser()
@@ -72,16 +74,18 @@
 
 apply_sensor_model = not args.disable_sensor_model
 
-image_simulator \
-    = desc.imsim.ImageSimulator(args.instcat, psf,
-                                numRows=args.numrows,
-                                config=args.config_file,
-                                seed=args.seed,
-                                outdir=args.outdir,
-                                sensor_list=sensor_list,
-                                apply_sensor_model=apply_sensor_model,
-                                create_centroid_file=args.create_centroid_file,
-                                file_id=args.file_id,
-                                log_level=args.log_level)
+with warnings.catch_warnings():
+    warnings.filterwarnings('ignore', 'ERFA', ErfaWarning)
+    image_simulator \
+        = desc.imsim.ImageSimulator(args.instcat, psf,
+                                    numRows=args.numrows,
+                                    config=args.config_file,
+                                    seed=args.seed,
+                                    outdir=args.outdir,
+                                    sensor_list=sensor_list,
+                                    apply_sensor_model=apply_sensor_model,
+                                    create_centroid_file=args.create_centroid_file,
+                                    file_id=args.file_id,
+                                    log_level=args.log_level)
 
-image_simulator.run(processes=args.processes)
+    image_simulator.run(processes=args.processes)

bin.src/make_flats.py

@@ -54,4 +54,6 @@
     my_flat = desc.imsim.make_flat(gs_det, counts_per_iter, niter, rng,
                                    logger=logger, wcs=wcs)
     ccd_id = "R{}_S{}".format(det_name[2:5:2], det_name[8:11:2])
-    my_flat.write('flat_{}_{}_{}.fits'.format(visit, ccd_id, obs_md.bandpass))
+    prefix = config['persistence']['eimage_prefix']
+    my_flat.write('{}_{}_{}_{}.fits'.format(prefix, visit, ccd_id,
+                                            obs_md.bandpass))

python/desc/imsim/ImageSimulator.py

@@ -10,6 +10,7 @@
 import shutil
 import sqlite3
 import numpy as np
+from astropy._erfa import ErfaWarning
 from lsst.afw.cameraGeom import WAVEFRONT, GUIDER
 from lsst.sims.photUtils import BandpassDict
 from lsst.sims.GalSimInterface import make_galsim_detector
@@ -88,15 +89,33 @@ def __init__(self, instcat, psf, numRows=None, config=None, seed=267,
         self.camera_wrapper = LSSTCameraWrapper()
         if sensor_list is None:
             sensor_list = self._get_all_sensors()
-        self.logger.debug("parsing instance catalog for %d sensors",
+        self.logger.debug("parsing instance catalog for %d sensor(s)",
                           len(sensor_list))
+        checkpoint_files = self._gather_checkpoint_files(sensor_list, file_id)
         self.obs_md, self.phot_params, sources \
-            = parsePhoSimInstanceFile(instcat, sensor_list, numRows=numRows)
+            = parsePhoSimInstanceFile(instcat, sensor_list, numRows=numRows,
+                                      checkpoint_files=checkpoint_files,
+                                      log_level=log_level)
         self.gs_obj_dict = sources[1]
         self.apply_sensor_model = apply_sensor_model
         self.file_id = file_id
         self._make_gs_interpreters(seed, sensor_list, file_id)
 
+    def _gather_checkpoint_files(self, sensor_list, file_id=None):
+        """
+        Gather any checkpoint files that have been created for the
+        desired sensors and return a dictionary with the discovered
+        filenames.
+        """
+        if file_id is None:
+            return None
+        checkpoint_files = dict()
+        for det_name in sensor_list:
+            filename = self.checkpoint_file(file_id, det_name)
+            if os.path.isfile(filename):
+                checkpoint_files[det_name] = filename
+        return checkpoint_files
+
     def _make_gs_interpreters(self, seed, sensor_list, file_id):
         """
         Create a separate GalSimInterpreter for each sensor so that
@@ -106,9 +125,6 @@ def _make_gs_interpreters(self, seed, sensor_list, file_id):
         Also extract GsObjectLists from gs_obj_dict for only the
         sensors in sensor_list so that the memory in the underlying
         InstCatTrimmer object in gs_obj_dict can be recovered.
-
-        TODO: Find a good way to pass a different seed to each
-        gs_interpreter or have them share the random number generator.
         """
         bp_dict = BandpassDict.loadTotalBandpassesFromFiles(bandpassNames=self.obs_md.bandpass)
         noise_and_background \
@@ -216,7 +232,8 @@ def run(self, processes=1, wait_time=None, node_id=0):
         # Set the CHECKPOINT_SUMMARY variable so that the summary info
         # from the subprocesses can be persisted.
         global CHECKPOINT_SUMMARY
-        if self.file_id is not None and CHECKPOINT_SUMMARY is None:
+        if (self.file_id is not None and processes > 1 and
+            CHECKPOINT_SUMMARY is None):
             visit = self.obs_md.OpsimMetaData['obshistID']
             CHECKPOINT_SUMMARY \
                 = CheckpointSummary(db_file='ckpt_{}_{}.sqlite3'.format(visit, node_id))
@@ -329,14 +346,17 @@ def __call__(self, gs_objects):
         with warnings.catch_warnings():
             warnings.filterwarnings('ignore', 'Automatic n_photons',
                                     UserWarning)
+            warnings.filterwarnings('ignore', 'ERFA function', ErfaWarning)
             nan_fluxes = 0
+            starting_for_loop = True
             for gs_obj in gs_objects:
+                if starting_for_loop:
+                    logger.info("drawing %d objects", len(gs_objects))
+                    starting_for_loop = False
                 if gs_obj.uniqueId in gs_interpreter.drawn_objects:
                     continue
                 flux = gs_obj.flux(IMAGE_SIMULATOR.obs_md.bandpass)
                 if not np.isnan(flux):
-                    if not gs_interpreter.drawn_objects:
-                        logger.info("drawing %d objects", len(gs_objects))
                     logger.debug("%s  %s  %s", gs_obj.uniqueId, flux,
                                  gs_obj.galSimType)
                     gs_interpreter.drawObject(gs_obj)

python/desc/imsim/flats.py

@@ -41,7 +41,7 @@ def make_flat(gs_det, counts_per_iter, niter, rng, buf=2, logger=None,
     -------
     galsim.ImageF
     """
-    config = desc.imsim.read_config()
+    config = desc.imsim.get_config()
 
     if logger is None:
         logger = desc.imsim.get_logger('DEBUG', name='make_flat')

python/desc/imsim/imSim.py

@@ -152,12 +152,16 @@ def metadata_from_file(file_name):
     return commands
 
 
-def sources_from_list(object_lines, obs_md, phot_params, file_name):
+def sources_from_list(object_lines, obs_md, phot_params, file_name,
+                      target_chip=None, log_level='INFO'):
     camera = get_obs_lsstSim_camera()
     config = get_config()
+    logger = get_logger(log_level, name=(target_chip if target_chip is
+                                         not None else 'sources_from_list'))
 
     num_objects = len(object_lines)
 
+    logger.debug('allocating object arrays')
     # RA, Dec in the coordinate system expected by PhoSim
     ra_phosim = np.zeros(num_objects, dtype=float)
     dec_phosim = np.zeros(num_objects, dtype=float)
@@ -186,6 +190,7 @@ def sources_from_list(object_lines, obs_md, phot_params, file_name):
     unique_id = [None]*num_objects
     object_type = np.zeros(num_objects, dtype=int)
 
+    logger.debug('looping over %s objects', num_objects)
     i_obj = -1
     for line in object_lines:
         params = line.strip().split()
@@ -308,6 +313,8 @@ def sources_from_list(object_lines, obs_md, phot_params, file_name):
 
     my_sed_dirs = sed_dirs(file_name)
 
+    logger.debug('constructing GalSimCelestialObjects for %s objects',
+                 num_objects)
     gs_object_arr = []
     for i_obj in range(num_objects):
         if not object_is_valid[i_obj]:
@@ -375,9 +382,13 @@ def sources_from_list(object_lines, obs_md, phot_params, file_name):
     assert len(x_pupil) == len(gs_object_arr)
     assert len(y_pupil) == len(gs_object_arr)
 
+    logger.debug('down-selecting by chip')
     out_obj_dict = {}
     for det in lsst_camera():
         chip_name = det.getName()
+        if target_chip is not None and chip_name != target_chip:
+            # Only consider the target_chip.
+            continue
         pixel_corners = getCornerPixels(chip_name, lsst_camera())
         x_min = pixel_corners[0][0]
         x_max = pixel_corners[2][0]
@@ -395,6 +406,8 @@ def sources_from_list(object_lines, obs_md, phot_params, file_name):
 
         out_obj_dict[chip_name] = gs_object_arr[on_chip]
 
+    if target_chip is not None:
+        logger.debug('objects remaining %s', len(out_obj_dict[target_chip]))
     return gs_object_arr, out_obj_dict
 
 def get_image_dirs():
@@ -492,28 +505,35 @@ def phosim_obs_metadata(phosim_commands):
     # Set the OpsimMetaData attribute with the obshistID info.
     obs_md.OpsimMetaData = {'obshistID': phosim_commands['obshistid']}
     obs_md.OpsimMetaData['FWHMgeom'] = fwhm_geom
-    obs_md.OpsimMetaData['FWHMeff'] =  fwhm_eff
+    obs_md.OpsimMetaData['FWHMeff'] = fwhm_eff
     obs_md.OpsimMetaData['rawSeeing'] = phosim_commands['seeing']
     obs_md.OpsimMetaData['altitude'] = phosim_commands['altitude']
     obs_md.OpsimMetaData['seed'] = phosim_commands['seed']
     return obs_md
 
 
-def parsePhoSimInstanceFile(fileName, sensor_list, numRows=None):
+def parsePhoSimInstanceFile(fileName, sensor_list, numRows=None,
+                            checkpoint_files=None, log_level='INFO'):
     """
     Read a PhoSim instance catalog into a Pandas dataFrame. Then use
     the information that was read-in to build and return a command
     dictionary and object dataFrame.
 
     Parameters
     ----------
-    fileName : str
+    fileName: str
         The instance catalog filename.
     sensor_list: list
         List of sensors for which to extract object lists.
-    numRows : int, optional
+    numRows: int [None]
         The number of rows to read from the instance catalog.
-        If None (the default), then all of the rows will be read in.
+        If None, then all of the rows will be read in.
+    checkpoint_files: dict [None]
+        Checkpoint files keyed by sensor name, e.g., "R:2,2 S:1,1".
+        The instance catalog lines corresponding to drawn_objects in
+        the checkpoint files will be skipped on ingest.
+    log_level: str ['INFO']
+        Logging level.
 
     Returns
     -------
@@ -522,13 +542,18 @@ def parsePhoSimInstanceFile(fileName, sensor_list, numRows=None):
         original DataFrames containing the header lines and object
         lines
     """
-
+    logger = get_logger(log_level, 'parsePhoSimInstanceFile')
     commands = metadata_from_file(fileName)
     obs_metadata = phosim_obs_metadata(commands)
     phot_params = photometricParameters(commands)
-    instcats = InstCatTrimmer(fileName, sensor_list, numRows=numRows)
+    logger.debug('creating InstCatTrimmer object')
+    instcats = InstCatTrimmer(fileName, sensor_list, numRows=numRows,
+                              checkpoint_files=checkpoint_files,
+                              log_level=log_level)
     gs_object_dict = {detname: GsObjectList(instcats[detname], instcats.obs_md,
-                                            phot_params, instcats.instcat_file)
+                                            phot_params, instcats.instcat_file,
+                                            chip_name=detname,
+                                            log_level=log_level)
                       for detname in sensor_list}
 
     return PhoSimInstanceCatalogContents(obs_metadata,
@@ -543,20 +568,23 @@ class GsObjectList:
     until items in the list are accessed.
     """
     def __init__(self, object_lines, obs_md, phot_params, file_name,
-                 chip_name=None):
+                 chip_name=None, log_level='INFO'):
         self.object_lines = object_lines
         self.obs_md = obs_md
         self.phot_params = phot_params
         self.file_name = file_name
         self.chip_name = chip_name
+        self.log_level = log_level
         self._gs_objects = None
 
     @property
     def gs_objects(self):
         if self._gs_objects is None:
             obj_arr, obj_dict \
                 = sources_from_list(self.object_lines, self.obs_md,
-                                    self.phot_params, self.file_name)
+                                    self.phot_params, self.file_name,
+                                    target_chip=self.chip_name,
+                                    log_level=self.log_level)
             if self.chip_name is not None:
                 try:
                     self._gs_objects = obj_dict[self.chip_name]
@@ -754,8 +782,6 @@ def add_cosmic_rays(gs_interpreter, phot_params):
         gs_interpreter.detectorImages[name] = \
             galsim.Image(crs.paint(imarr, exptime=exptime), wcs=image.wcs)
 
-    return None
-
 
 def add_treering_info(detectors, tr_filename=None):
     """