Merge 4146648 into 4a552df

bin.src/imsim.py

@@ -5,6 +5,7 @@
 read phoSim instance files as is. It leverages the LSST Sims GalSim interface
 code found in sims_GalSimInterface.
 """
+import os
 import argparse
 import desc.imsim
 
@@ -37,13 +38,24 @@
                     help='integer used to seed random number generator')
 parser.add_argument('--processes', type=int, default=1,
                     help='number of processes to use in multiprocessing mode')
+parser.add_argument('--psf_file', type=str, default=None,
+                    help="Pickle file containing for the persisted PSF. "
+                    "If the file exists, the psf will be loaded from that "
+                    "file, ignoring the --psf option; "
+                    "if not, a PSF will be created and saved to that filename.")
+
 args = parser.parse_args()
 
 commands = desc.imsim.metadata_from_file(args.instcat)
 
 obs_md = desc.imsim.phosim_obs_metadata(commands)
 
-psf = desc.imsim.make_psf(args.psf, obs_md, log_level=args.log_level)
+if args.psf_file is None or not os.path.isfile(args.psf_file):
+    psf = desc.imsim.make_psf(args.psf, obs_md, log_level=args.log_level)
+    if args.psf_file is not None:
+        desc.imsim.save_psf(psf, args.psf_file)
+else:
+    psf = desc.imsim.load_psf(args.psf_file)
 
 sensor_list = args.sensors.split('^') if args.sensors is not None \
     else args.sensors

python/desc/imsim/ImageSimulator.py

@@ -97,7 +97,7 @@ def _make_gs_interpreters(self, seed, sensor_list, file_id):
         """
         bp_dict = BandpassDict.loadTotalBandpassesFromFiles(bandpassNames=self.obs_md.bandpass)
         noise_and_background \
-            = make_sky_model(self.obs_md, self.phot_params,
+            = make_sky_model(self.obs_md, self.phot_params, seed=seed,
                              apply_sensor_model=self.apply_sensor_model)
         self.gs_interpreters = dict()
         for det in self.camera_wrapper.camera:

python/desc/imsim/atmPSF.py

@@ -46,11 +46,13 @@ class AtmosphericPSF(PSFbase):
     @param exptime      Exposure time in seconds.  default: 30.
     @param kcrit        Critical Fourier mode at which to split first and second kicks
                         in units of (1/r0).  default: 0.2
+    @param screen_size  Size of the phase screens in meters.  default: 819.2
+    @param screen_scale Size of phase screen "pixels" in meters.  default: 0.1
     @param doOpt        Add in optical phase screens?  default: True
     @param logger       Optional logger.  default: None
     """
     def __init__(self, airmass, rawSeeing, band, rng, t0=0.0, exptime=30.0, kcrit=0.2,
-                 doOpt=True, logger=None):
+                 screen_size=819.2, screen_scale=0.1, doOpt=True, logger=None):
         self.airmass = airmass
         self.rawSeeing = rawSeeing
 
@@ -62,6 +64,8 @@ def __init__(self, airmass, rawSeeing, band, rng, t0=0.0, exptime=30.0, kcrit=0.
         self.rng = rng
         self.t0 = t0
         self.exptime = exptime
+        self.screen_size = screen_size
+        self.screen_scale = screen_scale
         self.logger = logger
 
         self.atm = galsim.Atmosphere(**self._getAtmKwargs())
@@ -139,7 +143,7 @@ def _getAtmKwargs(self):
 
         return dict(r0_500=r0_500, L0=L0, speed=speeds, direction=directions,
                     altitude=altitudes, r0_weights=weights, rng=self.rng,
-                    screen_size=819.2, screen_scale=0.1)
+                    screen_size=self.screen_size, screen_scale=self.screen_scale)
 
     def _getPSF(self, xPupil=None, yPupil=None, gsparams=None):
         """

python/desc/imsim/imSim.py

@@ -6,6 +6,7 @@
 import sys
 import warnings
 from collections import namedtuple, defaultdict
+import pickle
 import logging
 import gc
 import copy
@@ -58,7 +59,8 @@
            'add_cosmic_rays',
            '_POINT_SOURCE', '_SERSIC_2D', '_RANDOM_WALK',
            'parsePhoSimInstanceFile',
-           'add_treering_info', 'airmass', 'FWHMeff', 'FWHMgeom', 'make_psf']
+           'add_treering_info', 'airmass', 'FWHMeff', 'FWHMgeom', 'make_psf',
+           'save_psf', 'load_psf']
 
 
 class PhosimInstanceCatalogParseError(RuntimeError):
@@ -459,6 +461,7 @@ def phosim_obs_metadata(phosim_commands):
     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
 
 
@@ -680,6 +683,8 @@ def read_config(config_file=None):
             my_config.set_from_config(section, key, value)
     return my_config
 
+read_config()
+
 
 def get_logger(log_level, name=None):
     """
@@ -854,7 +859,7 @@ def FWHMgeom(rawSeeing, band, altitude):
     return 0.822*FWHMeff(rawSeeing, band, altitude) + 0.052
 
 
-def make_psf(psf_name, obs_md, log_level='WARN', rng=None):
+def make_psf(psf_name, obs_md, log_level='WARN', rng=None, **kwds):
     """
     Make the requested PSF object.
 
@@ -870,10 +875,13 @@ def make_psf(psf_name, obs_md, log_level='WARN', rng=None):
         Logging level ('DEBUG', 'INFO', 'WARN', 'ERROR', 'CRITICAL').
     rng: galsim.BaseDeviate
         Instance of the galsim.baseDeviate random number generator.
+    **kwds: **dict
+        Additional keyword arguments to pass to the AtmosphericPSF,
+        i.e., screen_size(=819.2) and screen_scale(=0.1).
 
     Returns
     -------
-    lsst.sims.GalSimInterface.PSFbase:  Instance of a subclass of PSFbase.
+    lsst.sims.GalSimInterface.PSFbase: Instance of a subclass of PSFbase.
     """
     if psf_name.lower() == 'doublegaussian':
         return SNRdocumentPSF(obs_md.OpsimMetaData['FWHMgeom'])
@@ -888,11 +896,31 @@ def make_psf(psf_name, obs_md, log_level='WARN', rng=None):
                                           band=obs_md.bandpass)
     elif psf_name.lower() == 'atmospheric':
         if rng is None:
-            rng = galsim.UniformDeviate()
+            # Use the 'seed' value from the instance catalog for the rng
+            # used by the atmospheric PSF.
+            rng = galsim.UniformDeviate(obs_md.OpsimMetaData['seed'])
         logger = get_logger(log_level, 'psf')
         psf = AtmosphericPSF(airmass=my_airmass,
                              rawSeeing=rawSeeing,
                              band=obs_md.bandpass,
                              rng=rng,
-                             logger=logger)
+                             logger=logger, **kwds)
+    return psf
+
+def save_psf(psf, outfile):
+    """
+    Save the psf as a pickle file.
+    """
+    # Set any logger attribute to None since loggers cannot be persisted.
+    if hasattr(psf, 'logger'):
+        psf.logger = None
+    with open(outfile, 'wb') as output:
+        pickle.dump(psf, output)
+
+def load_psf(psf_file):
+    """
+    Load a psf from a pickle file.
+    """
+    with open(psf_file, 'rb') as fd:
+        psf = pickle.load(fd)
     return psf

tests/test_psf.py

@@ -0,0 +1,57 @@
+"""
+Unit tests for PSF-related code.
+"""
+import os
+import glob
+import unittest
+import desc.imsim
+
+
+def are_psfs_equal(psf1, psf2):
+    """
+    Test that two PSFs are equal.  For PSFs implemented in
+    sims_GalSimInterface, it is sufficient to compare the
+    ._cached_psf attributes.  For the AtmosphericPSF, only certain
+    attributes are meaningful to compare. See issue #117.
+    """
+    if type(psf1) != type(psf2):
+        return False
+    if not isinstance(psf1, desc.imsim.atmPSF.AtmosphericPSF):
+        # Compare cached galsim objects.
+        return psf1._cached_psf == psf2._cached_psf
+    # See issue #117 for an explanation of these comparisons:
+    return (psf1.atm[:-1] == psf2.atm[:-1]) and (psf1.aper == psf2.aper)
+
+
+class PsfTestCase(unittest.TestCase):
+    """
+    TestCase class for PSF-related functions.
+    """
+    def setUp(self):
+        self.test_dir = os.path.join(os.environ['IMSIM_DIR'], 'tests',
+                                     'psf_tests_dir')
+        os.makedirs(self.test_dir)
+
+    def tearDown(self):
+        for item in glob.glob(os.path.join(self.test_dir, '*')):
+            os.remove(item)
+        os.rmdir(self.test_dir)
+
+    def test_save_and_load_psf(self):
+        """
+        Test that the different imSim PSFs are saved and retrieved
+        correctly.
+        """
+        instcat = os.path.join(os.environ['IMSIM_DIR'], 'tests',
+                               'tiny_instcat.txt')
+        obs_md, _, _ = desc.imsim.parsePhoSimInstanceFile(instcat)
+        for psf_name in ("DoubleGaussian", "Kolmogorov", "Atmospheric"):
+            psf = desc.imsim.make_psf(psf_name, obs_md, screen_scale=0.4)
+            psf_file = os.path.join(self.test_dir, '{}.pkl'.format(psf_name))
+            desc.imsim.save_psf(psf, psf_file)
+            psf_retrieved = desc.imsim.load_psf(psf_file)
+            self.assertTrue(are_psfs_equal(psf, psf_retrieved))
+
+
+if __name__ == '__main__':
+    unittest.main()