From 7acd236ae2c38a29b02c327a7ee1514c7dd89f4e Mon Sep 17 00:00:00 2001 From: James Chiang Date: Fri, 23 Aug 2019 14:59:30 -0700 Subject: [PATCH 01/10] use generator for gs_object_list to save memory --- python/desc/imsim/imSim.py | 207 +++++++++++++++++++++++++++++-------- 1 file changed, 162 insertions(+), 45 deletions(-) diff --git a/python/desc/imsim/imSim.py b/python/desc/imsim/imSim.py index b74e5c0b..38ebe92d 100644 --- a/python/desc/imsim/imSim.py +++ b/python/desc/imsim/imSim.py @@ -611,67 +611,184 @@ def parsePhoSimInstanceFile(fileName, sensor_list, numRows=None, log_level=log_level, sort_magnorm=sort_magnorm) gs_object_dict = {detname: GsObjectList(instcats[detname], instcats.obs_md, phot_params, instcats.instcat_file, - chip_name=detname, - log_level=log_level) + detname, log_level=log_level) for detname in sensor_list} return PhoSimInstanceCatalogContents(obs_metadata, phot_params, ([], gs_object_dict)) - class GsObjectList: - """ - List-like class to provide access to lists of objects from an - instance catalog, deferring creation of GalSimCelestialObjects - until items in the list are accessed. - """ - def __init__(self, object_lines, obs_md, phot_params, file_name, - chip_name=None, log_level='INFO'): - self.object_lines = object_lines + def __init__(self, object_lines, obs_md, phot_params, instcat_file, + chip_name, log_level='INFO'): self.obs_md = obs_md + self.logger = get_logger(log_level, name=chip_name) + self._find_objects_on_chip(object_lines, chip_name) 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, - 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] - except KeyError: - self._gs_objects = [] - else: - self._gs_objects = obj_arr - return self._gs_objects + self.sed_dirs = sed_dirs(instcat_file) + config = get_config() + self.gamma2_sign = config['wl_params']['gamma2_sign'] + + def _find_objects_on_chip(self, object_lines, chip_name): + num_lines = len(object_lines) + ra_phosim = np.zeros(num_lines, dtype=float) + dec_phosim = np.zeros(num_lines, dtype=float) + mag_norm = 55.*np.ones(num_lines, dtype=float) + for i, line in enumerate(object_lines): + if not line.startswith('object'): + raise RuntimeError('Trying to process non-object entry from ' + 'the instance catalog.') + tokens = line.split() + ra_phosim[i] = float(tokens[2]) + dec_phosim[i] = float(tokens[3]) + mag_norm[i] = float(tokens[4]) + ra_appGeo, dec_appGeo \ + = PhoSimAstrometryBase._appGeoFromPhoSim(np.radians(ra_phosim), + np.radians(dec_phosim), + self.obs_md) + ra_obs_rad, dec_obs_rad \ + = _observedFromAppGeo(ra_appGeo, dec_appGeo, + obs_metadata=self.obs_md, + includeRefraction=True) + x_pupil, y_pupil = _pupilCoordsFromObserved(ra_obs_rad, dec_obs_rad, + self.obs_md) + on_chip_dict = _chip_downselect(mag_norm, x_pupil, y_pupil, + self.logger, [chip_name]) + index = on_chip_dict[chip_name] + self.object_lines = [] + for i in index[0]: + self.object_lines.append(object_lines[i]) + self.x_pupil = list(x_pupil[index]) + self.y_pupil = list(y_pupil[index]) + self.bp_dict = BandpassDict.loadTotalBandpassesFromFiles() def reset(self): - """ - Reset the ._gs_objects attribute to None in order to recover - memory devoted to the GalSimCelestialObject instances. - """ - self._gs_objects = None + pass def __len__(self): - try: - return len(self._gs_objects) - except TypeError: - return len(self.object_lines) + return len(self.object_lines) def __iter__(self): - for gs_obj in self.gs_objects: - yield gs_obj + return self - def __getitem__(self, index): - return self.gs_objects[index] + def __next__(self): + try: + while True: + gs_obj = self._make_gs_object(self.object_lines.pop(0), + self.x_pupil.pop(0), + self.y_pupil.pop(0)) + if gs_obj is not None: + return gs_obj + except IndexError: + raise StopIteration() + + def _make_gs_object(self, object_line, x_pupil, y_pupil): + params = object_line.strip().split() + unique_id = params[1] + ra_phosim = float(params[2]) + dec_phosim = float(params[3]) + mag_norm = float(params[4]) + sed_name = params[5] + redshift = float(params[6]) + gamma1 = float(params[7]) + gamma2 = self.gamma2_sign*float(params[8]) + kappa = float(params[9]) + internal_av = 0 + internal_rv = 0 + galactic_av = 0 + galactic_rv = 0 + fits_file = None + pixel_scale = 0 + rotation_angle = 0 + npoints = 0 + semi_major_arcsec = 0 + semi_minor_arcsec = 0 + position_angle_degrees = 0 + sersic_index = 0 + if params[12].lower() == 'point': + gs_type = 'pointSource' + i_gal_dust_model = 14 + if params[13].lower() != 'none': + i_gal_dust_model = 16 + internal_av = float(params[14]) + internal_rv =float(params[15]) + if params[i_gal_dust_model].lower() != 'none': + galactic_av = float(params[i_gal_dust_model+1]) + galactic_rv = float(params[i_gal_dust_model+2]) + elif params[12].lower() == 'sersic2d': + gs_type = 'sersic' + semi_major_arcsec = float(params[13]) + semi_minor_arcsec = float(params[14]) + position_angle_degrees = float(params[15]) + sersic_index = float(params[16]) + i_gal_dust_model = 18 + if params[17].lower() != 'none': + i_gal_dust_model = 20 + internal_av = float(params[18]) + internal_rv = float(params[19]) + if params[i_gal_dust_model].lower() != 'none': + galactic_av = float(params[i_gal_dust_model+1]) + galactic_rv = float(params[i_gal_dust_model+2]) + elif params[12].lower() == 'knots': + gs_type = 'RandomWalk' + semi_major_arcsec = float(params[13]) + semi_minor_arcsec = float(params[14]) + position_angle_degrees = float(params[15]) + npoints = int(params[16]) + i_gal_dust_model = 18 + if params[17].lower() != 'none': + i_gal_dust_model = 20 + internal_av = float(params[18]) + internal_rv = float(params[19]) + if params[i_gal_dust_model].lower() != 'none': + galactic_av = float(params[i_gal_dust_model+1]) + galactic_rv = float(params[i_gal_dust_model+2]) + elif (params[12].endswith('.fits') or params[12].endswith('.fits.gz')): + gs_type = 'FitsImage' + fits_file = find_file_path(params[12], get_image_dirs()) + pixel_scale = float(params[13]) + rotation_angle = float(params[14]) + i_gal_dust_model = 16 + if params[15].lower() != 'none': + i_gal_dust_model = 18 + internal_av = float(params[16]) + internal_rv = float(params[17]) + if params[i_gal_dust_model].lower() != 'none': + galactic_av = float(params[i_gal_dust_model+1]) + galactic_rv = float(params[i_gal_dust_model+2]) + else: + raise RuntimeError("Do not know how to handle " + "object type: %s" % params[12]) + + object_is_valid = (mag_norm < 50.0 and + (galactic_av != 0 or galactic_rv != 0) and + not (gs_type == 'sersic' and + semi_major_arcsec < semi_minor_arcsec) and + not (gs_type == 'RandomWalk' and npoints <=0)) + if not object_is_valid: + return None + sed_obj = SedWrapper(find_file_path(sed_name, self.sed_dirs), + mag_norm, redshift, internal_av, internal_rv, + galactic_av, galactic_rv, self.bp_dict) + position_angle_radians = np.radians(360 - position_angle_degrees) + gs_object = GalSimCelestialObject(gs_type, x_pupil, y_pupil, + radiansFromArcsec(semi_major_arcsec), + radiansFromArcsec(semi_minor_arcsec), + radiansFromArcsec(semi_major_arcsec), + position_angle_radians, + sersic_index, + sed_obj, + self.bp_dict, + self.phot_params, + npoints, + fits_file, + pixel_scale, + rotation_angle, + gamma1=gamma1, + gamma2=gamma2, + kappa=kappa, + uniqueId=unique_id) + return gs_object class ImSimConfiguration(object): From 8fea54fba533038bcd5b8d58a5228c83b542d9f6 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Fri, 23 Aug 2019 22:07:20 +0000 Subject: [PATCH 02/10] remove check for warning messages for instcat parsing --- tests/test_instcat_parser.py | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/tests/test_instcat_parser.py b/tests/test_instcat_parser.py index f4ccffb9..e58d5c6a 100644 --- a/tests/test_instcat_parser.py +++ b/tests/test_instcat_parser.py @@ -400,10 +400,6 @@ def test_validate_phosim_object_list(self): desired_warning_dex = i_ww break - if desired_warning_dex<0: - raise RuntimeError("Expected warning about bad sources " - "not issued") - message = wa[desired_warning_dex].message.args[0] # these are the objects that should be omitted @@ -412,12 +408,6 @@ def test_validate_phosim_object_list(self): 811883374596, 956090392580, 34307989098523, 34304522113056]]) - self.assertIn('Omitted 6 suspicious objects', message) - self.assertIn('4 had galactic_Av', message) - self.assertIn('1 had mag_norm', message) - self.assertIn('1 had semi_major_axis', message) - self.assertIn('1 had n_points', message) - self.assertEqual(len(my_objs), 18) for obj in instcat_contents.sources[0]: self.assertNotIn(obj.uniqueId, bad_unique_ids) From 0cd4a446b7269d2706d6c3caf4cbd74f744cced0 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Fri, 23 Aug 2019 15:29:54 -0700 Subject: [PATCH 03/10] remove the rest of warnings handling code from the instcat parser tests --- tests/test_instcat_parser.py | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/tests/test_instcat_parser.py b/tests/test_instcat_parser.py index e58d5c6a..2467169c 100644 --- a/tests/test_instcat_parser.py +++ b/tests/test_instcat_parser.py @@ -389,18 +389,6 @@ def test_validate_phosim_object_list(self): my_objs = set() for sensor in sensors: [my_objs.add(x) for x in instcat_contents.sources[1][sensor]] - self.assertGreater(len(wa), 0) - - # we must detect which warning is the warning we are actually - # testing, because PALPY keeps raising ERFAWarnings over our - # request for dates in the future - desired_warning_dex = -1 - for i_ww, ww in enumerate(wa): - if 'Omitted 6 suspicious objects' in ww.message.args[0]: - desired_warning_dex = i_ww - break - - message = wa[desired_warning_dex].message.args[0] # these are the objects that should be omitted bad_unique_ids = set([str(x) for x in From 0090dfde5049ca6a33ec912e8510797f36da92fb Mon Sep 17 00:00:00 2001 From: James Chiang Date: Tue, 27 Aug 2019 14:35:07 -0700 Subject: [PATCH 04/10] reset chunk_size so that smaller numRows is used --- python/desc/imsim/trim.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/python/desc/imsim/trim.py b/python/desc/imsim/trim.py index 5753d629..b3b0e9aa 100644 --- a/python/desc/imsim/trim.py +++ b/python/desc/imsim/trim.py @@ -207,6 +207,8 @@ def _process_objects(self, sensor_list, chunk_size, radius=0.18, """ num_gals = defaultdict(lambda: 0) self.update({sensor: [] for sensor in sensor_list}) + if numRows is not None: + chunk_size = min(chunk_size, numRows) with desc.imsim.fopen(self.instcat_file, mode='rt') as fd: nread = 0 while numRows is None or nread < numRows: From c3772c527e41ba882f9ff6c106907f051f885226 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Tue, 17 Sep 2019 15:20:59 -0700 Subject: [PATCH 05/10] silence afw.geom deprecation warnings --- python/desc/imsim/camera_info.py | 8 ++-- python/desc/imsim/camera_readout.py | 64 ++++++++++++++--------------- 2 files changed, 36 insertions(+), 36 deletions(-) diff --git a/python/desc/imsim/camera_info.py b/python/desc/imsim/camera_info.py index 42cf5fae..0b05cf03 100644 --- a/python/desc/imsim/camera_info.py +++ b/python/desc/imsim/camera_info.py @@ -2,7 +2,7 @@ Class to encapsulate info from lsst.obs.lsst.imsim.ImsimMapper().camera. """ import astropy.time -import lsst.afw.geom as afw_geom +import lsst.geom as lsst_geom from lsst.obs.lsst.imsim import ImsimMapper __all__ = ['CameraInfo', 'getHourAngle'] @@ -65,15 +65,15 @@ def mosaic_section(amp_info): Returns ------- - lsst.afw.geom.Box2I + lsst.geom.Box2I """ yseg, xseg = (int(x) for x in amp_info.getName()[-2:]) width = amp_info.getBBox().getWidth() height = amp_info.getBBox().getHeight() xmin = xseg*width ymin = 0 if yseg == 1 else height - return afw_geom.Box2I(afw_geom.Point2I(xmin, ymin), - afw_geom.Extent2I(width, height)) + return lsst_geom.Box2I(lsst_geom.Point2I(xmin, ymin), + lsst_geom.Extent2I(width, height)) def getHourAngle(observatory, mjd, ra): diff --git a/python/desc/imsim/camera_readout.py b/python/desc/imsim/camera_readout.py index 245d62d2..90d32161 100644 --- a/python/desc/imsim/camera_readout.py +++ b/python/desc/imsim/camera_readout.py @@ -23,7 +23,7 @@ from astropy.io import fits import astropy.time import galsim -import lsst.afw.geom as afwGeom +import lsst.geom as lsst_geom import lsst.afw.image as afwImage with warnings.catch_warnings(): warnings.simplefilter('ignore') @@ -567,7 +567,7 @@ def _noao_section_keyword(bbox, flipx=False, flipy=False): Parameters ---------- - bbox : lsst.afw.geom.Box2I + bbox : lsst.geom.Box2I Bounding box. flipx : bool Flag to indicate that data should be flipped in the x-direction. @@ -598,16 +598,16 @@ def set_itl_bboxes(amp): lsst.afw.table.tableLib.AmpInfoRecord The updated AmpInfoRecord. """ - amp.setRawBBox(afwGeom.Box2I(afwGeom.Point2I(0, 0), - afwGeom.Extent2I(544, 2048))) - amp.setRawDataBBox(afwGeom.Box2I(afwGeom.Point2I(3, 0), - afwGeom.Extent2I(509, 2000))) - amp.setRawHorizontalOverscanBBox(afwGeom.Box2I(afwGeom.Point2I(512, 0), - afwGeom.Extent2I(48, 2000))) - amp.setRawVerticalOverscanBBox(afwGeom.Box2I(afwGeom.Point2I(0, 2000), - afwGeom.Extent2I(544, 48))) - amp.setRawPrescanBBox(afwGeom.Box2I(afwGeom.Point2I(0, 0), - afwGeom.Extent2I(3, 2000))) + amp.setRawBBox(lsst_geom.Box2I(lsst_geom.Point2I(0, 0), + lsst_geom.Extent2I(544, 2048))) + amp.setRawDataBBox(lsst_geom.Box2I(lsst_geom.Point2I(3, 0), + lsst_geom.Extent2I(509, 2000))) + amp.setRawHorizontalOverscanBBox(lsst_geom.Box2I(lsst_geom.Point2I(512, 0), + lsst_geom.Extent2I(48, 2000))) + amp.setRawVerticalOverscanBBox(lsst_geom.Box2I(lsst_geom.Point2I(0, 2000), + lsst_geom.Extent2I(544, 48))) + amp.setRawPrescanBBox(lsst_geom.Box2I(lsst_geom.Point2I(0, 0), + lsst_geom.Extent2I(3, 2000))) return amp @@ -626,16 +626,16 @@ def set_e2v_bboxes(amp): lsst.afw.table.tableLib.AmpInfoRecord The updated AmpInfoRecord. """ - amp.setRawBBox(afwGeom.Box2I(afwGeom.Point2I(0, 0), - afwGeom.Extent2I(542, 2022))) - amp.setRawDataBBox(afwGeom.Box2I(afwGeom.Point2I(10, 0), - afwGeom.Extent2I(522, 2002))) - amp.setRawHorizontalOverscanBBox(afwGeom.Box2I(afwGeom.Point2I(522, 0), - afwGeom.Extent2I(20, 2002))) - amp.setRawVerticalOverscanBBox(afwGeom.Box2I(afwGeom.Point2I(0, 2002), - afwGeom.Extent2I(542, 20))) - amp.setRawPrescanBBox(afwGeom.Box2I(afwGeom.Point2I(0, 0), - afwGeom.Extent2I(10, 2002))) + amp.setRawBBox(lsst_geom.Box2I(lsst_geom.Point2I(0, 0), + lsst_geom.Extent2I(542, 2022))) + amp.setRawDataBBox(lsst_geom.Box2I(lsst_geom.Point2I(10, 0), + lsst_geom.Extent2I(522, 2002))) + amp.setRawHorizontalOverscanBBox(lsst_geom.Box2I(lsst_geom.Point2I(522, 0), + lsst_geom.Extent2I(20, 2002))) + amp.setRawVerticalOverscanBBox(lsst_geom.Box2I(lsst_geom.Point2I(0, 2002), + lsst_geom.Extent2I(542, 20))) + amp.setRawPrescanBBox(lsst_geom.Box2I(lsst_geom.Point2I(0, 0), + lsst_geom.Extent2I(10, 2002))) return amp @@ -654,16 +654,16 @@ def set_phosim_bboxes(amp): lsst.afw.table.tableLib.AmpInfoRecord The updated AmpInfoRecord. """ - amp.setRawBBox(afwGeom.Box2I(afwGeom.Point2I(0, 0), - afwGeom.Extent2I(519, 2001))) - amp.setRawDataBBox(afwGeom.Box2I(afwGeom.Point2I(4, 1), - afwGeom.Extent2I(509, 2000))) - amp.setRawHorizontalOverscanBBox(afwGeom.Box2I(afwGeom.Point2I(513, 1), - afwGeom.Extent2I(6, 2000))) - amp.setRawVerticalOverscanBBox(afwGeom.Box2I(afwGeom.Point2I(0, 2001), - afwGeom.Extent2I(519, 0))) - amp.setRawPrescanBBox(afwGeom.Box2I(afwGeom.Point2I(0, 1), - afwGeom.Extent2I(4, 2000))) + amp.setRawBBox(lsst_geom.Box2I(lsst_geom.Point2I(0, 0), + lsst_geom.Extent2I(519, 2001))) + amp.setRawDataBBox(lsst_geom.Box2I(lsst_geom.Point2I(4, 1), + lsst_geom.Extent2I(509, 2000))) + amp.setRawHorizontalOverscanBBox(lsst_geom.Box2I(lsst_geom.Point2I(513, 1), + lsst_geom.Extent2I(6, 2000))) + amp.setRawVerticalOverscanBBox(lsst_geom.Box2I(lsst_geom.Point2I(0, 2001), + lsst_geom.Extent2I(519, 0))) + amp.setRawPrescanBBox(lsst_geom.Box2I(lsst_geom.Point2I(0, 1), + lsst_geom.Extent2I(4, 2000))) return amp PixelParameters = namedtuple('PixelParameters', From 8bd6893e0eea82604fe1fa2f4f8254ceee0e71da Mon Sep 17 00:00:00 2001 From: James Chiang Date: Tue, 17 Sep 2019 15:50:15 -0700 Subject: [PATCH 06/10] use lsstdesc/stack-sims:w_2019_37-sims_w_2019_37 docker image in travis-ci build --- .travis.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.travis.yml b/.travis.yml index 6621e7de..e1438232 100644 --- a/.travis.yml +++ b/.travis.yml @@ -7,7 +7,7 @@ dist: trusty compiler: gcc env: - - OS_TYPE=centos OS_VERSION=7 DOCKER_IMAGE=lsstdesc/stack-sims:w_2019_23-sims_w_2019_23 + - OS_TYPE=centos OS_VERSION=7 DOCKER_IMAGE=lsstdesc/stack-sims:w_2019_37-sims_w_2019_37 services: - docker From dd2a632395a70e0fac47f30fd667a87f7864e2b0 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Tue, 17 Sep 2019 18:43:30 -0700 Subject: [PATCH 07/10] do setup -t current lsst_sims to use correct version of throughputs; clone master of sims_GalSimInterface to suppress afw.geom deprecation warnings --- travis_test.sh | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/travis_test.sh b/travis_test.sh index 29d37653..7501fce7 100755 --- a/travis_test.sh +++ b/travis_test.sh @@ -1,7 +1,7 @@ #!/bin/bash source scl_source enable devtoolset-6 source loadLSST.bash -setup lsst_sims +setup -t current lsst_sims pip install nose pip install coveralls pip install pylint @@ -16,6 +16,14 @@ setup -r . -j scons lib python shebang examples doc policy python/lsst/obs/lsst/version.py cd .. + +# Clone sims_GalSimInterface and set up master to suppress +# lsst.afw.geom deprecation warnings. +git clone https://github.com/lsst/sims_GalSimInterface.git +cd sims_GalSimInterface +setup -r . -j +cd .. + eups declare imsim -r ${TRAVIS_BUILD_DIR} -t current setup imsim cd ${TRAVIS_BUILD_DIR} From b89dc77ed959564831cd5ba2ed7018a1f65d488b Mon Sep 17 00:00:00 2001 From: James Chiang Date: Thu, 19 Sep 2019 13:08:09 -0700 Subject: [PATCH 08/10] remove sources_from_list --- python/desc/imsim/imSim.py | 266 +---------------------------------- tests/test_instcat_parser.py | 165 ++-------------------- 2 files changed, 17 insertions(+), 414 deletions(-) diff --git a/python/desc/imsim/imSim.py b/python/desc/imsim/imSim.py index 38ebe92d..228aa19c 100644 --- a/python/desc/imsim/imSim.py +++ b/python/desc/imsim/imSim.py @@ -54,7 +54,6 @@ __all__ = ['PhosimInstanceCatalogParseError', 'photometricParameters', 'phosim_obs_metadata', - 'sources_from_list', 'metadata_from_file', 'read_config', 'get_config', 'get_logger', 'get_image_dirs', 'get_obs_lsstSim_camera', @@ -62,7 +61,7 @@ '_POINT_SOURCE', '_SERSIC_2D', '_RANDOM_WALK', '_FITS_IMAGE', 'parsePhoSimInstanceFile', 'add_treering_info', 'airmass', 'FWHMeff', 'FWHMgeom', 'make_psf', - 'save_psf', 'load_psf', 'TracebackDecorator'] + 'save_psf', 'load_psf', 'TracebackDecorator', 'GsObjectList'] class PhosimInstanceCatalogParseError(RuntimeError): @@ -159,268 +158,6 @@ def uss_mem(): return my_process.memory_full_info().uss/1024.**3 -def sources_from_list(object_lines, obs_md, phot_params, file_name, - target_chip=None, log_level='INFO'): - """. - - Parse the object lines from a phosim-style instance catalog and - repackage as GalSimCelestialObjects after applying on-chip - selections and consistency cuts on mag_norm, extinction - parameters, and galaxy shape and knot parameters. - - Parameters - ---------- - object_lines: list - List of object line entries from the instance catalog. - obs_md: ObservationMetaData - Visit-specific metadata from the instance catalog. - phot_params: PhotometricParameters - Visit-specific photometric data needed for computing object - fluxes. - file_name: str - The filename path of the instance catalog, used for inferring - the location of the Dynamic SEDs used by the sprinkled objects. - target_chip: str [None] - For multiprocessing mode, this is the name of the chip - (e.g., "R:2,2 S:1,1") being simulated. If None, then - of the object lines are partitioned among all 189 science - sensors in the LSST focalplane. Since this is a costly - calculation, doing it only for the target chip saves a lot - of compute time. - log_level: str ['INFO'] - Logging level. - - Returns - ------- - list, dict: A list of all the GalSimCelestialObjects that pass - the focalplane-level selections and a dict of those objects keyed - by chip name. - """ - 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, %s GB', uss_mem()) - # 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) - - sed_name = [None]*num_objects - mag_norm = 55.0*np.ones(num_objects, dtype=float) - gamma1 = np.zeros(num_objects, dtype=float) - gamma2 = np.zeros(num_objects, dtype=float) - gamma2_sign = config['wl_params']['gamma2_sign'] - kappa = np.zeros(num_objects, dtype=float) - - internal_av = np.zeros(num_objects, dtype=float) - internal_rv = np.zeros(num_objects, dtype=float) - galactic_av = np.zeros(num_objects, dtype=float) - galactic_rv = np.zeros(num_objects, dtype=float) - semi_major_arcsec = np.zeros(num_objects, dtype=float) - semi_minor_arcsec = np.zeros(num_objects, dtype=float) - position_angle_degrees = np.zeros(num_objects, dtype=float) - sersic_index = np.zeros(num_objects, dtype=float) - npoints = np.zeros(num_objects, dtype=int) - redshift = np.zeros(num_objects, dtype=float) - pixel_scale = np.zeros(num_objects, dtype=float) - rotation_angle = np.zeros(num_objects, dtype=float) - fits_image_file = dict() - - unique_id = [None]*num_objects - object_type = np.zeros(num_objects, dtype=int) - - logger.debug('looping over %s objects; %s GB', num_objects, uss_mem()) - i_obj = -1 - for line in object_lines: - params = line.strip().split() - if params[0] != 'object': - continue - i_obj += 1 - unique_id[i_obj] = params[1] - ra_phosim[i_obj] = float(params[2]) - dec_phosim[i_obj] = float(params[3]) - mag_norm[i_obj] = float(params[4]) - sed_name[i_obj] = params[5] - redshift[i_obj] = float(params[6]) - gamma1[i_obj] = float(params[7]) - gamma2[i_obj] = gamma2_sign*float(params[8]) - kappa[i_obj] = float(params[9]) - if params[12].lower() == 'point': - object_type[i_obj] = _POINT_SOURCE - i_gal_dust_model = 14 - if params[13].lower() != 'none': - i_gal_dust_model = 16 - internal_av[i_obj] = float(params[14]) - internal_rv[i_obj] =float(params[15]) - if params[i_gal_dust_model].lower() != 'none': - galactic_av[i_obj] = float(params[i_gal_dust_model+1]) - galactic_rv[i_obj] = float(params[i_gal_dust_model+2]) - elif params[12].lower() == 'sersic2d': - object_type[i_obj] = _SERSIC_2D - semi_major_arcsec[i_obj] = float(params[13]) - semi_minor_arcsec[i_obj] = float(params[14]) - position_angle_degrees[i_obj] = float(params[15]) - sersic_index[i_obj] = float(params[16]) - i_gal_dust_model = 18 - if params[17].lower() != 'none': - i_gal_dust_model = 20 - internal_av[i_obj] = float(params[18]) - internal_rv[i_obj] = float(params[19]) - if params[i_gal_dust_model].lower() != 'none': - galactic_av[i_obj] = float(params[i_gal_dust_model+1]) - galactic_rv[i_obj] =float(params[i_gal_dust_model+2]) - elif params[12].lower() == 'knots': - object_type[i_obj] = _RANDOM_WALK - semi_major_arcsec[i_obj] = float(params[13]) - semi_minor_arcsec[i_obj] = float(params[14]) - position_angle_degrees[i_obj] = float(params[15]) - npoints[i_obj] = int(params[16]) - i_gal_dust_model = 18 - if params[17].lower() != 'none': - i_gal_dust_model = 20 - internal_av[i_obj] = float(params[18]) - internal_rv[i_obj] = float(params[19]) - if params[i_gal_dust_model].lower() != 'none': - galactic_av[i_obj] = float(params[i_gal_dust_model+1]) - galactic_rv[i_obj] = float(params[i_gal_dust_model+2]) - elif (params[12].endswith('.fits') or params[12].endswith('.fits.gz')): - object_type[i_obj] = _FITS_IMAGE - fits_image_file[i_obj] = find_file_path(params[12], get_image_dirs()) - pixel_scale[i_obj] = float(params[13]) - rotation_angle[i_obj] = float(params[14]) - i_gal_dust_model = 16 - if params[15].lower() != 'none': - i_gal_dust_model = 18 - internal_av[i_obj] = float(params[16]) - internal_rv[i_obj] = float(params[17]) - if params[i_gal_dust_model].lower() != 'none': - galactic_av[i_obj] = float(params[i_gal_dust_model+1]) - galactic_rv[i_obj] = float(params[i_gal_dust_model+2]) - else: - raise RuntimeError("Do not know how to handle " - "object type: %s" % params[12]) - - logger.debug("computing pupil coords, %s GB", uss_mem()) - ra_appGeo, dec_appGeo \ - = PhoSimAstrometryBase._appGeoFromPhoSim(np.radians(ra_phosim), - np.radians(dec_phosim), - obs_md) - - ra_obs_rad, dec_obs_rad \ - = _observedFromAppGeo(ra_appGeo, dec_appGeo, - obs_metadata=obs_md, - includeRefraction=True) - - semi_major_radians = radiansFromArcsec(semi_major_arcsec) - semi_minor_radians = radiansFromArcsec(semi_minor_arcsec) - # Account for PA sign difference wrt phosim convention. - position_angle_radians = np.radians(360. - position_angle_degrees) - - x_pupil, y_pupil = _pupilCoordsFromObserved(ra_obs_rad, - dec_obs_rad, - obs_md) - - bp_dict = BandpassDict.loadTotalBandpassesFromFiles() - - object_is_valid = np.array([True]*num_objects) - - invalid_objects = np.where(np.logical_or(np.logical_or( - mag_norm>50.0, - np.logical_and(galactic_av==0.0, galactic_rv==0.0)), - np.logical_or( - np.logical_and(object_type==_SERSIC_2D, - semi_major_arcsec 0: - message = "\nOmitted %d suspicious objects from " % len(invalid_objects[0]) - message += "the instance catalog:\n" - n_bad_mag_norm = len(np.where(mag_norm>50.0)[0]) - message += " %d had mag_norm > 50.0\n" % n_bad_mag_norm - n_bad_av = len(np.where(np.logical_and(galactic_av==0.0, galactic_rv==0.0))[0]) - message += " %d had galactic_Av == galactic_Rv == 0\n" % n_bad_av - n_bad_axes = len(np.where(np.logical_and(object_type==_SERSIC_2D, - semi_major_arcsec Date: Thu, 19 Sep 2019 13:20:48 -0700 Subject: [PATCH 09/10] remove GsObjectsList.reset method --- python/desc/imsim/ImageSimulator.py | 3 --- python/desc/imsim/imSim.py | 3 --- 2 files changed, 6 deletions(-) diff --git a/python/desc/imsim/ImageSimulator.py b/python/desc/imsim/ImageSimulator.py index 90bef76f..a26f918c 100644 --- a/python/desc/imsim/ImageSimulator.py +++ b/python/desc/imsim/ImageSimulator.py @@ -405,9 +405,6 @@ def __call__(self, gs_objects): if nan_fluxes > 0: logger.info("%s objects had nan fluxes", nan_fluxes) - # Recover the memory devoted to the GalSimCelestialObject instances. - gs_objects.reset() - add_cosmic_rays(gs_interpreter, IMAGE_SIMULATOR.phot_params) full_well = int(IMAGE_SIMULATOR.config['ccd']['full_well']) apply_channel_bleeding(gs_interpreter, full_well) diff --git a/python/desc/imsim/imSim.py b/python/desc/imsim/imSim.py index 228aa19c..df6510a0 100644 --- a/python/desc/imsim/imSim.py +++ b/python/desc/imsim/imSim.py @@ -400,9 +400,6 @@ def _find_objects_on_chip(self, object_lines, chip_name): self.y_pupil = list(y_pupil[index]) self.bp_dict = BandpassDict.loadTotalBandpassesFromFiles() - def reset(self): - pass - def __len__(self): return len(self.object_lines) From 79100d0a848841d6c26d452f9dececa839651cae Mon Sep 17 00:00:00 2001 From: James Chiang Date: Fri, 20 Sep 2019 13:43:18 -0700 Subject: [PATCH 10/10] restore flux, coordinate, and parameter tests for instance catalog reading --- tests/test_instcat_parser.py | 147 +++++++++++++++++++++++++++++++++-- 1 file changed, 140 insertions(+), 7 deletions(-) diff --git a/tests/test_instcat_parser.py b/tests/test_instcat_parser.py index 03bce20e..e8e075de 100644 --- a/tests/test_instcat_parser.py +++ b/tests/test_instcat_parser.py @@ -22,13 +22,20 @@ from lsst.sims.GalSimInterface import LSSTCameraWrapper def sources_from_list(lines, obs_md, phot_params, file_name): + """Return a two-item tuple containing + * a list of GalSimCelestialObjects for each object entry in `lines` + * a dictionary of these objects disaggregated by chip_name. + """ target_chips = [det.getName() for det in lsst_camera()] + gs_object_dict = dict() out_obj_dict = dict() for chip_name in target_chips: out_obj_dict[chip_name] \ - = desc.imsim.GsObjectList(lines, obs_md, phot_params, file_name, - chip_name) - return out_obj_dict + = [_ for _ in desc.imsim.GsObjectList(lines, obs_md, phot_params, + file_name, chip_name)] + for gsobj in out_obj_dict[chip_name]: + gs_object_dict[gsobj.uniqueId] = gsobj + return list(gs_object_dict.values()), out_obj_dict class InstanceCatalogParserTestCase(unittest.TestCase): """ @@ -143,11 +150,61 @@ def test_object_extraction_stars(self): dtype=truth_dtype, delimiter=';') truth_data.sort() + + gs_object_arr, gs_object_dict \ + = sources_from_list(lines, obs_md, phot_params, self.phosim_file) + + id_arr = [None]*len(gs_object_arr) + for i_obj in range(len(gs_object_arr)): + id_arr[i_obj] = gs_object_arr[i_obj].uniqueId + id_arr = sorted(id_arr) + np.testing.assert_array_equal(truth_data['uniqueId'], id_arr) + + ######## test that pupil coordinates are correct to within + ######## half a milliarcsecond + + x_pup_test, y_pup_test = _pupilCoordsFromRaDec(truth_data['raJ2000'], + truth_data['decJ2000'], + pm_ra=truth_data['pmRA'], + pm_dec=truth_data['pmDec'], + v_rad=truth_data['v_rad'], + parallax=truth_data['parallax'], + obs_metadata=obs_md) + + for gs_obj in gs_object_arr: + i_obj = np.where(truth_data['uniqueId'] == gs_obj.uniqueId)[0][0] + dd = np.sqrt((x_pup_test[i_obj]-gs_obj.xPupilRadians)**2 + + (y_pup_test[i_obj]-gs_obj.yPupilRadians)**2) + dd = arcsecFromRadians(dd) + self.assertLess(dd, 0.0005) + + ######## test that fluxes are correctly calculated + + bp_dict = BandpassDict.loadTotalBandpassesFromFiles() + imsim_bp = Bandpass() + imsim_bp.imsimBandpass() + phot_params = PhotometricParameters(nexp=1, exptime=30.0) + + for gs_obj in gs_object_arr: + i_obj = np.where(truth_data['uniqueId'] == gs_obj.uniqueId)[0][0] + sed = Sed() + full_sed_name = os.path.join(os.environ['SIMS_SED_LIBRARY_DIR'], + truth_data['sedFilename'][i_obj]) + sed.readSED_flambda(full_sed_name) + fnorm = sed.calcFluxNorm(truth_data['magNorm'][i_obj], imsim_bp) + sed.multiplyFluxNorm(fnorm) + sed.resampleSED(wavelen_match=bp_dict.wavelenMatch) + a_x, b_x = sed.setupCCM_ab() + sed.addDust(a_x, b_x, A_v=truth_data['Av'][i_obj], + R_v=truth_data['Rv'][i_obj]) + + for bp in ('u', 'g', 'r', 'i', 'z', 'y'): + flux = sed.calcADU(bp_dict[bp], phot_params)*phot_params.gain + self.assertAlmostEqual(flux/gs_obj.flux(bp), 1.0, 10) + ######## test that objects are assigned to the right chip in ######## gs_object_dict - gs_object_dict = sources_from_list(lines, obs_md, phot_params, - self.phosim_file) unique_id_dict = {} for chip_name in gs_object_dict: local_unique_id_list = [] @@ -200,11 +257,87 @@ def test_object_extraction_galaxies(self): truth_data.sort() + gs_object_arr, gs_object_dict \ + = sources_from_list(lines, obs_md, phot_params, galaxy_phosim_file) + + id_arr = [None]*len(gs_object_arr) + for i_obj in range(len(gs_object_arr)): + id_arr[i_obj] = gs_object_arr[i_obj].uniqueId + id_arr = sorted(id_arr) + np.testing.assert_array_equal(truth_data['uniqueId'], id_arr) + + ######## test that galaxy parameters are correctly read in + + g1 = truth_data['gamma1']/(1.0-truth_data['kappa']) + g2 = truth_data['gamma2']/(1.0-truth_data['kappa']) + mu = 1.0/((1.0-truth_data['kappa'])**2 - (truth_data['gamma1']**2 + truth_data['gamma2']**2)) + for gs_obj in gs_object_arr: + i_obj = np.where(truth_data['uniqueId'] == gs_obj.uniqueId)[0][0] + self.assertAlmostEqual(gs_obj.mu/mu[i_obj], 1.0, 6) + self.assertAlmostEqual(gs_obj.g1/g1[i_obj], 1.0, 6) + self.assertAlmostEqual(gs_obj.g2/g2[i_obj], 1.0, 6) + self.assertGreater(np.abs(gs_obj.mu), 0.0) + self.assertGreater(np.abs(gs_obj.g1), 0.0) + self.assertGreater(np.abs(gs_obj.g2), 0.0) + + self.assertAlmostEqual(gs_obj.halfLightRadiusRadians, + truth_data['majorAxis'][i_obj], 13) + self.assertAlmostEqual(gs_obj.minorAxisRadians, + truth_data['minorAxis'][i_obj], 13) + self.assertAlmostEqual(gs_obj.majorAxisRadians, + truth_data['majorAxis'][i_obj], 13) + self.assertAlmostEqual(2.*np.pi - gs_obj.positionAngleRadians, + truth_data['positionAngle'][i_obj], 7) + self.assertAlmostEqual(gs_obj.sindex, + truth_data['sindex'][i_obj], 10) + + ######## test that pupil coordinates are correct to within + ######## half a milliarcsecond + + x_pup_test, y_pup_test = _pupilCoordsFromRaDec(truth_data['raJ2000'], + truth_data['decJ2000'], + obs_metadata=obs_md) + + for gs_obj in gs_object_arr: + i_obj = np.where(truth_data['uniqueId'] == gs_obj.uniqueId)[0][0] + dd = np.sqrt((x_pup_test[i_obj]-gs_obj.xPupilRadians)**2 + + (y_pup_test[i_obj]-gs_obj.yPupilRadians)**2) + dd = arcsecFromRadians(dd) + self.assertLess(dd, 0.0005) + + ######## test that fluxes are correctly calculated + + bp_dict = BandpassDict.loadTotalBandpassesFromFiles() + imsim_bp = Bandpass() + imsim_bp.imsimBandpass() + phot_params = PhotometricParameters(nexp=1, exptime=30.0) + + for gs_obj in gs_object_arr: + i_obj = np.where(truth_data['uniqueId'] == gs_obj.uniqueId)[0][0] + sed = Sed() + full_sed_name = os.path.join(os.environ['SIMS_SED_LIBRARY_DIR'], + truth_data['sedFilename'][i_obj]) + sed.readSED_flambda(full_sed_name) + fnorm = sed.calcFluxNorm(truth_data['magNorm'][i_obj], imsim_bp) + sed.multiplyFluxNorm(fnorm) + + a_x, b_x = sed.setupCCM_ab() + sed.addDust(a_x, b_x, A_v=truth_data['internalAv'][i_obj], + R_v=truth_data['internalRv'][i_obj]) + + sed.redshiftSED(truth_data['redshift'][i_obj], dimming=True) + sed.resampleSED(wavelen_match=bp_dict.wavelenMatch) + a_x, b_x = sed.setupCCM_ab() + sed.addDust(a_x, b_x, A_v=truth_data['galacticAv'][i_obj], + R_v=truth_data['galacticRv'][i_obj]) + + for bp in ('u', 'g', 'r', 'i', 'z', 'y'): + flux = sed.calcADU(bp_dict[bp], phot_params)*phot_params.gain + self.assertAlmostEqual(flux/gs_obj.flux(bp), 1.0, 6) + ######## test that objects are assigned to the right chip in ######## gs_object_dict - gs_object_dict = sources_from_list(lines, obs_md, phot_params, - galaxy_phosim_file) unique_id_dict = {} for chip_name in gs_object_dict: local_unique_id_list = []