From 92d71518d3d5cce0d909059c0d482298edeb3ee1 Mon Sep 17 00:00:00 2001 From: yuedong0607 Date: Fri, 3 Apr 2026 03:33:57 -0400 Subject: [PATCH 01/11] remove galsim.roman calls; calling from romanisim.models instead --- roman_imsim/bandpass.py | 4 +-- roman_imsim/noise.py | 41 +++++++++++++++++------------- roman_imsim/psf.py | 56 +++++++++++++++++++++++++++++++++++------ roman_imsim/sca.py | 6 ++--- roman_imsim/scafile.py | 4 +-- roman_imsim/skycat.py | 8 +++--- roman_imsim/stamp.py | 4 +-- roman_imsim/utils.py | 8 +++--- roman_imsim/wcs.py | 8 +++--- 9 files changed, 93 insertions(+), 46 deletions(-) diff --git a/roman_imsim/bandpass.py b/roman_imsim/bandpass.py index 96472722..dd735f1a 100644 --- a/roman_imsim/bandpass.py +++ b/roman_imsim/bandpass.py @@ -1,4 +1,4 @@ -import galsim.roman as roman +import romanisim.models as models from galsim.config import BandpassBuilder, GetAllParams, RegisterBandpassType @@ -25,7 +25,7 @@ def buildBandpass(self, config, base, logger): kwargs, safe = GetAllParams(config, base, req=req) name = kwargs["name"] - bandpass = roman.getBandpasses(red_limit=2000)[name] + bandpass = models.bandpass.getBandpasses(red_limit=2000)[name] return bandpass, safe diff --git a/roman_imsim/noise.py b/roman_imsim/noise.py index 1d7bd711..f59d0926 100644 --- a/roman_imsim/noise.py +++ b/roman_imsim/noise.py @@ -1,5 +1,5 @@ import galsim -import galsim.roman as roman +import romanisim.models as models from galsim.config import NoiseBuilder, RegisterNoiseType from astropy.time import Time @@ -43,6 +43,7 @@ def addNoise(self, config, base, image, rng, current_var, draw_method, logger): "ipc": bool, "read_noise": bool, "sky_subtract": bool, + "use_crds": bool, } params, safe = galsim.config.GetAllParams(config, base, req={}, opt=opt, ignore=[]) @@ -56,6 +57,8 @@ def addNoise(self, config, base, image, rng, current_var, draw_method, logger): ipc = params.get("ipc", False) read_noise = params.get("read_noise", False) sky_subtract = params.get("sky_subtract", True) + + use_crds = params.get("use_crds", False) base["current_noise_image"] = base["current_image"] wcs = base["wcs"] @@ -69,17 +72,17 @@ def addNoise(self, config, base, image, rng, current_var, draw_method, logger): # value added to sky_image. So technically, this includes things that aren't just sky. # E.g. includes dark_current and thermal backgrounds. sky_image = image.copy() - sky_level = roman.getSkyLevel(bp, world_pos=wcs.toWorld(image.true_center), date=date) + sky_level = models.backgrounds.getSkyLevel(bp, world_pos=wcs.toWorld(image.true_center), date=date) logger.debug("Adding sky_level = %s", sky_level) if stray_light: - logger.debug("Stray light fraction = %s", roman.stray_light_fraction) - sky_level *= 1.0 + roman.stray_light_fraction + logger.debug("Stray light fraction = %s", models.parameters.stray_light_fraction) + sky_level *= 1.0 + models.parameters.stray_light_fraction wcs.makeSkyImage(sky_image, sky_level) # The other background is the expected thermal backgrounds in this band. # These are provided in e-/pix/s, so we have to multiply by the exposure time. if thermal_background: - tb = roman.thermal_backgrounds[filter_name] * exptime + tb = models.backgrounds.thermal_backgrounds[filter_name] * exptime logger.debug("Adding thermal background: %s", tb) sky_image += tb @@ -112,39 +115,41 @@ def addNoise(self, config, base, image, rng, current_var, draw_method, logger): if reciprocity_failure: logger.debug("Applying reciprocity failure") - roman.addReciprocityFailure(image) + models.nonlinearity.addReciprocityFailure(img=image) if dark_current: - dc = roman.dark_current * exptime - logger.debug("Adding dark current: %s", dc) - sky_image += dc - dark_noise = galsim.noise.DeviateNoise(galsim.random.PoissonDeviate(rng, dc)) - image.addNoise(dark_noise) + logger.debug("Adding dark current: %s") + dc = models.DarkCurrent(usecrds=use_crds) + dc.apply(img=image, exptime=exptime) if nonlinearity: logger.debug("Applying classical nonlinearity") - roman.applyNonlinearity(image) + non_linear = models.Nonlinearity(usecrds=use_crds) + non_linear.apply(img=image, electrons=False) # Mosby and Rauscher say there are two read noises. One happens before IPC, the other # one after. # TODO: Add read_noise1 if ipc: logger.debug("Applying IPC") - roman.applyIPC(image) + IPC = models.IPC(usecrds=use_crds) + IPC.apply(img=image) if read_noise: - logger.debug("Adding read noise %s", roman.read_noise) - image.addNoise(galsim.GaussianNoise(rng, sigma=roman.read_noise)) + logger.debug("Adding read noise") + rn = models.ReadNoise(usecrds=use_crds) + rn.apply(img=image) - logger.debug("Applying gain %s", roman.gain) - image /= roman.gain + logger.debug("Applying gain") + gain = models.Gain(usecrds=use_crds) + gain.apply(img=image) # Make integer ADU now. image.quantize() if sky_subtract: logger.debug("Subtracting sky image") - sky_image /= roman.gain + gain.apply(img=sky_image) sky_image.quantize() image -= sky_image diff --git a/roman_imsim/psf.py b/roman_imsim/psf.py index 033dbd10..e859a44c 100644 --- a/roman_imsim/psf.py +++ b/roman_imsim/psf.py @@ -57,7 +57,6 @@ """ import galsim -import galsim.roman as roman from galsim.config import ( InputLoader, RegisterInputType, @@ -66,6 +65,7 @@ RegisterObjectType, ) from galsim.errors import GalSimConfigValueError +import romanisim.models as models ########################## # PSF Interpolator Input # @@ -77,6 +77,48 @@ class PSFInterpolator: """Base class for PSF interpolator""" + +class RomanPSF(object): + """Class building needed Roman PSFs.""" + + def __init__( + self, + SCA=None, + WCS=None, + n_waves=None, + bpass=None, + extra_aberrations=None, + logger=None, + ): + + logger = galsim.config.LoggerWrapper(logger) + + if n_waves == -1: + if bpass.name == "W146": + n_waves = 10 + else: + n_waves = 5 + + corners = [ + galsim.PositionD(1, 1), + galsim.PositionD(1, models.parameters.n_pix), + galsim.PositionD(models.parameters.n_pix, 1), + galsim.PositionD(models.parameters.n_pix, models.parameters.n_pix), + ] + cc = galsim.PositionD(models.parameters.n_pix / 2, models.parameters.n_pix / 2) + tags = ["ll", "lu", "ul", "uu"] + self.PSF = {} + pupil_bin = 8 + self.PSF[pupil_bin] = {} + for tag, SCA_pos in tuple(zip(tags, corners)): + self.PSF[pupil_bin][tag] = self._psf_call( + SCA, bpass, SCA_pos, WCS, pupil_bin, n_waves, logger, extra_aberrations + ) + for pupil_bin in [4, 2, "achromatic"]: + self.PSF[pupil_bin] = self._psf_call( + SCA, bpass, cc, WCS, pupil_bin, n_waves, logger, extra_aberrations + ) + def _parse_pupil_bin(self, pupil_bin): if pupil_bin == "achromatic": return 8 @@ -86,7 +128,7 @@ def _parse_pupil_bin(self, pupil_bin): def _psf_call(self, SCA, bpass, SCA_pos, WCS, pupil_bin, n_waves, logger, extra_aberrations): if pupil_bin == 8: - psf = roman.getPSF( + psf = models.psf_utils.getPSF( SCA, bpass.name, SCA_pos=SCA_pos, @@ -100,7 +142,7 @@ def _psf_call(self, SCA, bpass, SCA_pos, WCS, pupil_bin, n_waves, logger, extra_ extra_aberrations=extra_aberrations, ) else: - psf = roman.getPSF( + psf = models.psf_utils.getPSF( SCA, bpass.name, SCA_pos=SCA_pos, @@ -260,12 +302,12 @@ def getPSF(self, pupil_bin, pos): if pupil_bin != 8: return psf - wll = (self._image_xsize - pos.x) * (self._image_ysize - pos.y) - wlu = (self._image_xsize - pos.x) * (pos.y - 1) - wul = (pos.x - 1) * (self._image_ysize - pos.y) + wll = (models.parameters.n_pix - pos.x) * (models.parameters.n_pix - pos.y) + wlu = (models.parameters.n_pix - pos.x) * (pos.y - 1) + wul = (pos.x - 1) * (models.parameters.n_pix - pos.y) wuu = (pos.x - 1) * (pos.y - 1) return (wll * psf["ll"] + wlu * psf["lu"] + wul * psf["ul"] + wuu * psf["uu"]) / ( - (self._image_xsize - 1) * (self._image_ysize - 1) + (models.parameters.n_pix - 1) * (models.parameters.n_pix - 1) ) diff --git a/roman_imsim/sca.py b/roman_imsim/sca.py index 7d362440..7524a21a 100644 --- a/roman_imsim/sca.py +++ b/roman_imsim/sca.py @@ -2,7 +2,7 @@ import galsim import galsim.config -import galsim.roman as roman +import romanisim.models as models import numpy as np from astropy.time import Time from galsim.config import RegisterImageType @@ -67,7 +67,7 @@ def setup(self, config, base, image_num, obj_num, ignore, logger): if "bandpass" not in config: base["bandpass"] = galsim.config.BuildBandpass(base["image"], "bandpass", base, logger=logger) - return roman.n_pix, roman.n_pix + return models.parameters.n_pix, models.parameters.n_pix def buildImage(self, config, base, image_num, obj_num, logger): """Build an Image containing multiple objects placed at arbitrary locations. @@ -100,7 +100,7 @@ def buildImage(self, config, base, image_num, obj_num, logger): full_image.header["MJD-OBS"] = self.mjd full_image.header["DATE-OBS"] = Time(self.mjd, format="mjd").datetime.isoformat() full_image.header["FILTER"] = self.filter - full_image.header["ZPTMAG"] = 2.5 * np.log10(self.exptime * roman.collecting_area) + full_image.header["ZPTMAG"] = 2.5 * np.log10(self.exptime * models.parameters.collecting_area) base["current_image"] = full_image diff --git a/roman_imsim/scafile.py b/roman_imsim/scafile.py index de1bea8e..75574a64 100644 --- a/roman_imsim/scafile.py +++ b/roman_imsim/scafile.py @@ -4,7 +4,7 @@ import astropy.time import galsim -import galsim.roman as roman +import romanisim.models as models from galsim.config import OutputBuilder, RegisterOutputType @@ -29,7 +29,7 @@ def setup(self, config, base, file_num, logger): if "exptime" in config: base["exptime"] = galsim.config.ParseValue(config, "exptime", base, float)[0] else: - base["exptime"] = roman.exptime + base["exptime"] = models.parameters.exptime # Save the detector size, so the input catalogs can use it to figure out which # objects will be visible. diff --git a/roman_imsim/skycat.py b/roman_imsim/skycat.py index 0b2b8c9b..bf31eca1 100644 --- a/roman_imsim/skycat.py +++ b/roman_imsim/skycat.py @@ -3,7 +3,7 @@ """ import galsim -import galsim.roman as roman +import romanisim.models as models import numpy as np from galsim.config import ( InputLoader, @@ -68,11 +68,11 @@ def __init__( if xsize is not None: self.xsize = xsize else: - self.xsize = roman.n_pix + self.xsize = models.parameters.n_pix if ysize is not None: self.ysize = ysize else: - self.ysize = roman.n_pix + self.ysize = models.parameters.n_pix self.obj_types = obj_types self.edge_pix = edge_pix self.logger = galsim.config.LoggerWrapper(logger) @@ -184,7 +184,7 @@ def getObj(self, index, gsparams=None, rng=None, exptime=30): # Compute the flux or get the cached value. flux = ( - skycat_obj.get_roman_flux(self.bandpass.name, mjd=self.mjd) * self.exptime * roman.collecting_area + skycat_obj.get_roman_flux(self.bandpass.name, mjd=self.mjd) * self.exptime * models.parameters.collecting_area ) if np.isnan(flux): return None diff --git a/roman_imsim/stamp.py b/roman_imsim/stamp.py index 2554d180..f491fe06 100644 --- a/roman_imsim/stamp.py +++ b/roman_imsim/stamp.py @@ -1,6 +1,6 @@ import galsim import galsim.config -import galsim.roman as roman +import romanisim.models as models import numpy as np from galsim.config import RegisterStampType, StampBuilder @@ -100,7 +100,7 @@ def setup(self, config, base, xsize, ysize, ignore, logger): # psf = galsim.config.BuildGSObject(base, 'psf', logger=logger)[0]['achromatic'] # obj = galsim.Convolve(gal_achrom, psf).withFlux(self.flux) obj = gal_achrom.withGSParams(galsim.GSParams(stepk_minimum_hlr=20)) - image_size = obj.getGoodImageSize(roman.pixel_scale) + image_size = obj.getGoodImageSize(models.parameters.pixel_scale) # print('stamp setup3',process.memory_info().rss) base["pupil_bin"] = self.pupil_bin diff --git a/roman_imsim/utils.py b/roman_imsim/utils.py index 408f8c97..8b6d0219 100644 --- a/roman_imsim/utils.py +++ b/roman_imsim/utils.py @@ -2,7 +2,7 @@ import galsim import galsim.config -import galsim.roman as roman +import romanisim.models as models import numpy as np @@ -72,9 +72,9 @@ def getPSF(self, x=None, y=None, pupil_bin=8): raise ValueError( "x,y position for pupil_bin values other than 8 not supported. Using SCA center." ) - return self.PSF.getPSF(pupil_bin, galsim.PositionD(roman.n_pix / 2, roman.n_pix / 2)) + return self.PSF.getPSF(pupil_bin, galsim.PositionD(models.parameters.n_pix / 2, models.parameters.n_pix / 2)) if (x is None) | (y is None): - return self.PSF.getPSF(8, galsim.PositionD(roman.n_pix / 2, roman.n_pix / 2)) + return self.PSF.getPSF(8, galsim.PositionD(models.parameters.n_pix / 2, models.parameters.n_pix / 2)) return self.PSF.getPSF(8, galsim.PositionD(x, y)) def getWCS(self): @@ -148,7 +148,7 @@ def getPSF_Image( return psf.drawImage(self.bpass, image=stamp, wcs=wcs, method=method) photon_ops = [self.getPSF(x, y, pupil_bin)] + self.photon_ops if include_pixel and oversampling_factor > 1: - scale = galsim.roman.pixel_scale + scale = models.parameters.pixel_scale # An array (comb) of Dirac Delta functions. comb = galsim.Add( [ diff --git a/roman_imsim/wcs.py b/roman_imsim/wcs.py index 94700318..9e7951b0 100644 --- a/roman_imsim/wcs.py +++ b/roman_imsim/wcs.py @@ -1,5 +1,5 @@ import galsim -import galsim.roman as roman +import romanisim.models as models from astropy.time import Time from galsim.angle import Angle from galsim.celestial import CelestialCoord @@ -21,10 +21,10 @@ def buildWCS(self, config, base, logger): kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt) if "max_sun_angle" in kwargs: - roman.max_sun_angle = kwargs["max_sun_angle"] - roman.roman_wcs.max_sun_angle = kwargs["max_sun_angle"] + models.parameters.max_sun_angle = kwargs["max_sun_angle"] + models.wcs_utils.max_sun_angle = kwargs["max_sun_angle"] pointing = CelestialCoord(ra=kwargs["ra"], dec=kwargs["dec"]) - wcs = roman.getWCS( + wcs = models.wcs_utils.getWCS( world_pos=pointing, PA=kwargs["pa"], date=Time(kwargs["mjd"], format="mjd").datetime, From f74a5dc1e92fe816ed86ef303d50a0e6cfe75388 Mon Sep 17 00:00:00 2001 From: yuedong0607 Date: Mon, 20 Apr 2026 01:39:49 -0400 Subject: [PATCH 02/11] enabling the option for using per-SCA throughput --- roman_imsim/bandpass.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/roman_imsim/bandpass.py b/roman_imsim/bandpass.py index dd735f1a..19b657b5 100644 --- a/roman_imsim/bandpass.py +++ b/roman_imsim/bandpass.py @@ -22,10 +22,14 @@ def buildBandpass(self, config, base, logger): the constructed Bandpass object. """ req = {"name": str} - kwargs, safe = GetAllParams(config, base, req=req) + opt = {"SCA": int} + kwargs, safe = GetAllParams(config, base, req=req, opt=opt) name = kwargs["name"] - bandpass = models.bandpass.getBandpasses(red_limit=2000)[name] + if "SCA" in kwargs: + bandpass = models.bandpass.getBandpasses(red_limit=2000, sca=kwargs["SCA"])[name] + else: + bandpass = models.bandpass.getBandpasses(red_limit=2000)[name] return bandpass, safe From 1f670fefad5ae4b0c2b7007274606db565ddc62f Mon Sep 17 00:00:00 2001 From: yuedong0607 Date: Tue, 21 Apr 2026 05:11:09 -0400 Subject: [PATCH 03/11] add "per-effect" options to effects in noise class --- roman_imsim/config/hack.yaml | 175 +++++++++++++++++++++++++++++++++++ roman_imsim/noise.py | 48 +++++----- 2 files changed, 202 insertions(+), 21 deletions(-) create mode 100644 roman_imsim/config/hack.yaml diff --git a/roman_imsim/config/hack.yaml b/roman_imsim/config/hack.yaml new file mode 100644 index 00000000..e9556791 --- /dev/null +++ b/roman_imsim/config/hack.yaml @@ -0,0 +1,175 @@ +# Default settings for roman simulation +# Includes creation of noiseless oversampled images (including PSF) +# - processing of other detector and instrument effects are still handled in the +# python postprocessing layer to enable things not currently in galsim.roman + +modules: + + # Including galsim.roman in the list of modules to import will add a number of Roman-specific + # functions and classes that we will use here. + - roman_imsim + - galsim.roman + + # We need this for one of our Eval items. GalSim does not by default import datetime into + # the globals dict it uses when evaluating Eval items, so we can tell it to import it here. + - datetime + +# Define some other information about the images +image: + + # A special Image type that knows all the Roman SCA geometry, WCS, gain, etc. + # It also by default applies a number of detector effects, but these can be turned + # off if desired by setting some parameters (given below) to False. + type: roman_sca + + wcs: + type: RomanWCS + SCA: '@image.SCA' + ra: { type: ObSeqData, field: ra } + dec: { type: ObSeqData, field: dec } + pa: { type: ObSeqData, field: pa } + mjd: { type: ObSeqData, field: mjd } + + bandpass: + type: RomanBandpassTrimmed + SCA: '@image.SCA' + name: { type: ObSeqData, field: filter } + + # When you want to have multiple images generate the same random galaxies, then + # you can set up multiple random number generators with different update cadences + # by making random_seed a list. + # The default behavior is just to have the random seeds for each object go in order by + # object number across all images, but this shows how to set it up so we use two separate + # cadences. + # The first one behaves normally, which will be used for things like noise on the image. + # The second one sets the initial seed for each object to repeat to the same starting value + # at the start of each filter. If we were doing more than 3 total files, it would then + # move on to another sequence for the next 3 and so on. + random_seed: + # Used for noise and nobjects. + - { type: ObSeqData, field: visit } + + # Used for objects. Repeats sequence for each filter + # Note: Don't use $ shorthand here, since that will implicitly be evaluated once and then + # treated the same way as an integer (i.e. making a regular sequence starting from that + # value). Using an explicit dict with an Eval type means GalSim will leave it alone and + # evaluate it as is for each object. + + + # We're just doing one SCA here. + # If you wanted to do all of them in each of three filters (given below), you could use: + # + # SCA: + # type: Sequence + # first: 1 + # last: 18 + # repeat: 3 # repeat each SCA num 3 times before moving on, for the 3 filters. + # + SCA: 4 + mjd: { type: ObSeqData, field: mjd } + filter: { type: ObSeqData, field: filter } + exptime: { type: ObSeqData, field: exptime } + + # Photon shooting is way faster for chromatic objects than fft, especially when most of them + # are fairly faint. The cross-over point for achromatic objects is generally of order + # flux=1.e6 or so (depending on the profile). Most of our objects here are much fainter than + # that. The fft rendering for chromatic is a factor of 10 or so slower still, whereas + # chromatic photon shooting is only slightly slower than achromatic, so the difference + # is even more pronounced in this case. + draw_method: 'phot' + use_fft_bright: True + + noise: + # These are all by default turned on, but you can turn any of them off if desired: + type: RomanNoise + mjd: { type: ObSeqData, field: mjd } + stray_light: True + thermal_background: True + reciprocity_failure: True + dark_current: + turn_on: True + use_crds: True + nonlinearity: + turn_on: True + use_crds: True + ipc: + turn_on: True + use_crds: True + read_noise: + turn_on: True + use_crds: True + gain: + turn_on: True + use_crds: True + sky_subtract: False + # noise: + # type: NoNoise + +stamp: + type: Roman_stamp + world_pos: + type: SkyCatWorldPos + exptime: { type: ObSeqData, field: exptime } + skip_failures: True + photon_ops: + - { type: ChargeDiff } + +# psf: +# type: roman_psf +# # If omitted, it would figure this out automatically, because we are using the RomanSCA image +# # type. But if we weren't, you'd have to tell it which SCA to build the PSF for. +# SCA: '@image.SCA' +# # n_waves defines how finely to sample the PSF profile over the bandpass. +# # Using 10 wavelengths usually gives decent accuracy. +# n_waves: 10 + +# Define the galaxy type and positions to use +gal: + type: SkyCatObj + +input: + obseq_data: + file_name: "Roman_WAS_obseq_11_1_23.fits" + visit: 12909 + SCA: '@image.SCA' + roman_psf: + SCA: '@image.SCA' + n_waves: 5 + sky_catalog: + file_name: "skyCatalog/skyCatalog.yaml" + edge_pix: 512 + mjd: { type: ObSeqData, field: mjd } + exptime: { type: ObSeqData, field: exptime } + obj_types: ['diffsky_galaxy','star','snana','gaia_star'] + +output: + + nfiles: 1 + dir: output/RomanWAS_new/images/truth + file_name: + type: FormattedStr + format: "Roman_WAS_truth_%s_%i_%i.fits.gz" + items: + - { type: ObSeqData, field: filter } + - { type: ObSeqData, field: visit } + - '@image.SCA' + + truth: + dir: output/RomanWAS_new/truth + file_name: + type: FormattedStr + format: "Roman_WAS_index_%s_%i_%i.txt" + items: + - { type: ObSeqData, field: filter } + - { type: ObSeqData, field: visit } + - '@image.SCA' + columns: + object_id: "@object_id" + ra: "$sky_pos.ra.deg" + dec: "$sky_pos.dec.deg" + x: "$image_pos.x" + y: "$image_pos.y" + realized_flux: "@realized_flux" + flux: "@flux" + mag: "@mag" + obj_type: "@object_type" diff --git a/roman_imsim/noise.py b/roman_imsim/noise.py index f59d0926..7bc0df20 100644 --- a/roman_imsim/noise.py +++ b/roman_imsim/noise.py @@ -38,12 +38,12 @@ def addNoise(self, config, base, image, rng, current_var, draw_method, logger): "stray_light": bool, "thermal_background": bool, "reciprocity_failure": bool, - "dark_current": bool, - "nonlinearity": bool, - "ipc": bool, - "read_noise": bool, + "dark_current": dict, + "nonlinearity": dict, + "ipc": dict, + "read_noise": dict, "sky_subtract": bool, - "use_crds": bool, + "gain": dict, } params, safe = galsim.config.GetAllParams(config, base, req={}, opt=opt, ignore=[]) @@ -52,11 +52,13 @@ def addNoise(self, config, base, image, rng, current_var, draw_method, logger): stray_light = params.get("stray_light", False) thermal_background = params.get("thermal_background", False) reciprocity_failure = params.get("reciprocity_failure", False) - dark_current = params.get("dark_current", False) - nonlinearity = params.get("nonlinearity", False) - ipc = params.get("ipc", False) - read_noise = params.get("read_noise", False) + # dark_current = params.get("dark_current", False) + dark_current = params.get("dark_current", {'turn_on': False, 'use_crds': False}) + nonlinearity = params.get("nonlinearity", {'turn_on': False, 'use_crds': False}) + ipc = params.get("ipc", {'turn_on': False, 'use_crds': False}) + read_noise = params.get("read_noise", {'turn_on': False, 'use_crds': False}) sky_subtract = params.get("sky_subtract", True) + gain = params.get('gain', {'turn_on': False, 'use_crds': False}) use_crds = params.get("use_crds", False) @@ -117,39 +119,43 @@ def addNoise(self, config, base, image, rng, current_var, draw_method, logger): logger.debug("Applying reciprocity failure") models.nonlinearity.addReciprocityFailure(img=image) - if dark_current: + if dark_current['turn_on']: logger.debug("Adding dark current: %s") - dc = models.DarkCurrent(usecrds=use_crds) + # print(f'for drak_current, use_crds = {dark_current["use_crds"]}') + dc = models.DarkCurrent(usecrds=dark_current['use_crds']) dc.apply(img=image, exptime=exptime) - if nonlinearity: + if nonlinearity['turn_on']: logger.debug("Applying classical nonlinearity") - non_linear = models.Nonlinearity(usecrds=use_crds) + non_linear = models.Nonlinearity(usecrds=nonlinearity['use_crds']) non_linear.apply(img=image, electrons=False) # Mosby and Rauscher say there are two read noises. One happens before IPC, the other # one after. # TODO: Add read_noise1 - if ipc: + if ipc['turn_on']: logger.debug("Applying IPC") - IPC = models.IPC(usecrds=use_crds) + IPC = models.IPC(usecrds=ipc['use_crds']) IPC.apply(img=image) - if read_noise: + if read_noise['turn_on']: logger.debug("Adding read noise") - rn = models.ReadNoise(usecrds=use_crds) + rn = models.ReadNoise(usecrds=read_noise['use_crds']) rn.apply(img=image) - logger.debug("Applying gain") - gain = models.Gain(usecrds=use_crds) - gain.apply(img=image) + if gain['turn_on']: + logger.debug("Applying gain") + gain_model = models.Gain(usecrds=gain['use_crds']) + gain_model.apply(img=image) # Make integer ADU now. image.quantize() if sky_subtract: logger.debug("Subtracting sky image") - gain.apply(img=sky_image) + if gain['turn_on']: + logger.debug("Applying gain") + gain_model.apply(img=sky_image) sky_image.quantize() image -= sky_image From d86f6a90826864f188e0606e23a925ca44531ebd Mon Sep 17 00:00:00 2001 From: Yuedong Fang Date: Tue, 28 Apr 2026 03:01:09 -0400 Subject: [PATCH 04/11] include romanisim in project.toml --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index 3e616be4..7205bf76 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,6 +23,7 @@ dependencies = [ "skycatalogs", "packaging", "setuptools", + "roman_imsim", ] [project.urls] From 7fd70e4ef0f7ece38bcc055595e607ff8c0e07b0 Mon Sep 17 00:00:00 2001 From: Yuedong Fang Date: Tue, 28 Apr 2026 03:01:35 -0400 Subject: [PATCH 05/11] fix a typo --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 7205bf76..b505fadd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,7 +23,7 @@ dependencies = [ "skycatalogs", "packaging", "setuptools", - "roman_imsim", + "romanisim", ] [project.urls] From 04dc9bb9a4f76ace0b5dbd1c9b07bacaa1b7c830 Mon Sep 17 00:00:00 2001 From: Yuedong Fang Date: Thu, 30 Apr 2026 16:33:01 -0400 Subject: [PATCH 06/11] reformatted --- roman_imsim/bandpass.py | 4 +- roman_imsim/detector_physics.py | 293 +++++++++++++++++++++++--------- roman_imsim/noise.py | 50 +++--- roman_imsim/obseq.py | 20 ++- roman_imsim/psf.py | 21 ++- roman_imsim/sca.py | 24 ++- roman_imsim/scafile.py | 8 +- roman_imsim/skycat.py | 40 +++-- roman_imsim/stamp.py | 70 ++++++-- roman_imsim/utils.py | 46 +++-- roman_imsim/wcs.py | 1 - tests/test_end_to_end.py | 12 +- tests/test_imports.py | 8 +- 13 files changed, 441 insertions(+), 156 deletions(-) diff --git a/roman_imsim/bandpass.py b/roman_imsim/bandpass.py index 19b657b5..aec401c7 100644 --- a/roman_imsim/bandpass.py +++ b/roman_imsim/bandpass.py @@ -27,7 +27,9 @@ def buildBandpass(self, config, base, logger): name = kwargs["name"] if "SCA" in kwargs: - bandpass = models.bandpass.getBandpasses(red_limit=2000, sca=kwargs["SCA"])[name] + bandpass = models.bandpass.getBandpasses(red_limit=2000, sca=kwargs["SCA"])[ + name + ] else: bandpass = models.bandpass.getBandpasses(red_limit=2000)[name] diff --git a/roman_imsim/detector_physics.py b/roman_imsim/detector_physics.py index 7e58f8d4..3e29c366 100644 --- a/roman_imsim/detector_physics.py +++ b/roman_imsim/detector_physics.py @@ -75,22 +75,30 @@ def __init__(self, params, visit, SCA): self.exptime = obseq_data.ob["exptime"] self.bpass = roman.getBandpasses()[self.filter] self.WCS = roman.getWCS( - world_pos=galsim.CelestialCoord(ra=obseq_data.ob["ra"], dec=obseq_data.ob["dec"]), + world_pos=galsim.CelestialCoord( + ra=obseq_data.ob["ra"], dec=obseq_data.ob["dec"] + ), PA=obseq_data.ob["pa"], date=self.date, SCAs=self.sca, PA_is_FPA=True, )[self.sca] - self.radec = self.WCS.toWorld(galsim.PositionI(int(roman.n_pix / 2), int(roman.n_pix / 2))) + self.radec = self.WCS.toWorld( + galsim.PositionI(int(roman.n_pix / 2), int(roman.n_pix / 2)) + ) self.WCS = roman.getWCS( - world_pos=galsim.CelestialCoord(ra=obseq_data.ob["ra"], dec=obseq_data.ob["dec"]), + world_pos=galsim.CelestialCoord( + ra=obseq_data.ob["ra"], dec=obseq_data.ob["dec"] + ), PA=obseq_data.ob["pa"], date=self.date, SCAs=self.sca, PA_is_FPA=True, )[self.sca] - self.radec = self.WCS.toWorld(galsim.PositionI(int(roman.n_pix / 2), int(roman.n_pix / 2))) + self.radec = self.WCS.toWorld( + galsim.PositionI(int(roman.n_pix / 2), int(roman.n_pix / 2)) + ) class modify_image(object): @@ -98,7 +106,9 @@ class modify_image(object): Class to simulate non-idealities and noise of roman detector images. """ - def __init__(self, params, visit, sca, dither_from_file, sca_filepath=None, use_galsim=False): + def __init__( + self, params, visit, sca, dither_from_file, sca_filepath=None, use_galsim=False + ): """ Set up noise properties of image @@ -119,9 +129,13 @@ def __init__(self, params, visit, sca, dither_from_file, sca_filepath=None, use_ # Load sca file if applicable if sca_filepath is not None: - self.df = fio.FITS(sca_filepath + "/" + sca_number_to_file[self.pointing.sca]) + self.df = fio.FITS( + sca_filepath + "/" + sca_number_to_file[self.pointing.sca] + ) print("------- Using SCA files --------") - self.df = fio.FITS(sca_filepath + "/" + sca_number_to_file[self.pointing.sca]) + self.df = fio.FITS( + sca_filepath + "/" + sca_number_to_file[self.pointing.sca] + ) print("------- Using SCA files --------") else: self.df = None @@ -132,9 +146,15 @@ def __init__(self, params, visit, sca, dither_from_file, sca_filepath=None, use_ old_filename = os.path.join(self.params["output"]["dir"], imfilename) if not os.path.exists( - self.params["output"]["dir"].replace("truth", self.get_path_name(use_galsim=use_galsim)) + self.params["output"]["dir"].replace( + "truth", self.get_path_name(use_galsim=use_galsim) + ) ): - os.mkdir(self.params["output"]["dir"].replace("truth", self.get_path_name(use_galsim=use_galsim))) + os.mkdir( + self.params["output"]["dir"].replace( + "truth", self.get_path_name(use_galsim=use_galsim) + ) + ) new_filename = os.path.join(self.params["output"]["dir"], imfilename).replace( "truth", self.get_path_name(use_galsim=use_galsim) ) @@ -142,9 +162,15 @@ def __init__(self, params, visit, sca, dither_from_file, sca_filepath=None, use_ b = galsim.BoundsI(xmin=1, ymin=1, xmax=roman.n_pix, ymax=roman.n_pix) old_filename = os.path.join(self.params["output"]["dir"], imfilename) if not os.path.exists( - self.params["output"]["dir"].replace("truth", self.get_path_name(use_galsim=use_galsim)) + self.params["output"]["dir"].replace( + "truth", self.get_path_name(use_galsim=use_galsim) + ) ): - os.mkdir(self.params["output"]["dir"].replace("truth", self.get_path_name(use_galsim=use_galsim))) + os.mkdir( + self.params["output"]["dir"].replace( + "truth", self.get_path_name(use_galsim=use_galsim) + ) + ) new_filename = os.path.join(self.params["output"]["dir"], imfilename).replace( "truth", self.get_path_name(use_galsim=use_galsim) ) @@ -160,8 +186,12 @@ def __init__(self, params, visit, sca, dither_from_file, sca_filepath=None, use_ force_cvz = True self.setup_sky(im, self.pointing, rng, visit * sca, force_cvz=force_cvz) - img, err, dq, sky_mean, sky_noise = self.add_effects(im, None, self.pointing, use_galsim=use_galsim) - img, err, dq, sky_mean, sky_noise = self.add_effects(im, None, self.pointing, use_galsim=use_galsim) + img, err, dq, sky_mean, sky_noise = self.add_effects( + im, None, self.pointing, use_galsim=use_galsim + ) + img, err, dq, sky_mean, sky_noise = self.add_effects( + im, None, self.pointing, use_galsim=use_galsim + ) write_fits( old_filename, @@ -172,7 +202,15 @@ def __init__(self, params, visit, sca, dither_from_file, sca_filepath=None, use_ self.pointing.sca, sky_mean=sky_mean, ) - write_fits(old_filename, new_filename, img, sky_noise, dq, self.pointing.sca, sky_mean=sky_mean) + write_fits( + old_filename, + new_filename, + img, + sky_noise, + dq, + self.pointing.sca, + sky_mean=sky_mean, + ) def get_path_name(self, use_galsim=False): @@ -309,7 +347,13 @@ def add_effects_scafile(self, im, wt, pointing): self.sky_mean, sky_noise, ) - return im, self.sky[self.sky.bounds & im.bounds] - self.sky_mean, dq, self.sky_mean, sky_noise + return ( + im, + self.sky[self.sky.bounds & im.bounds] - self.sky_mean, + dq, + self.sky_mean, + sky_noise, + ) def add_effects_galsim(self, im, wt, pointing): """ @@ -439,7 +483,13 @@ def add_effects_simple(self, im, wt, pointing): self.sky_mean, sky_noise, ) - return im, self.sky[self.sky.bounds & im.bounds] - self.sky_mean, dq, self.sky_mean, sky_noise + return ( + im, + self.sky[self.sky.bounds & im.bounds] - self.sky_mean, + dq, + self.sky_mean, + sky_noise, + ) def set_diff(self, im=None): if self.params["save_diff"]: @@ -498,7 +548,9 @@ def bfe(self, im): # solve boundary shfit kernel aX components # ======================================================================= a_area = self.df["BFE"][:, :, :, :] # 5x5x32x32 - a_components = np.zeros((4, 2 * nbfe + 1, 2 * nbfe + 1, n_max, m_max)) # 4x5x5x32x32 + a_components = np.zeros( + (4, 2 * nbfe + 1, 2 * nbfe + 1, n_max, m_max) + ) # 4x5x5x32x32 # solve aR aT aL aB for each a for n in range(n_max): # m_max and n_max = 32 (binned in 128x128) @@ -509,7 +561,16 @@ def bfe(self, im): a = (a + np.fliplr(a) + np.flipud(a) + np.flip(a)) / 4.0 r = 0.5 * (3.25 / 4.25) ** (1.5) / 1.5 # source-boundary projection - B = (a[2, 2], a[3, 2], a[2, 3], a[3, 3], a[4, 2], a[2, 4], a[3, 4], a[4, 4]) + B = ( + a[2, 2], + a[3, 2], + a[2, 3], + a[3, 3], + a[4, 2], + a[2, 4], + a[3, 4], + a[4, 4], + ) A = np.array( [ @@ -562,11 +623,15 @@ def bfe(self, im): # The img is clipped by the saturation level here to cap the brighter fatter effect # and avoid unphysical behavior - array_pad = self.saturate(im.copy()).array[4:-4, 4:-4] # img of interest 4088x4088 + array_pad = self.saturate(im.copy()).array[ + 4:-4, 4:-4 + ] # img of interest 4088x4088 array_pad = np.pad( array_pad, [(4 + nbfe, 4 + nbfe), (4 + nbfe, 4 + nbfe)], mode="symmetric" ) # 4100x4100 array - array_pad = self.saturate(im.copy()).array[4:-4, 4:-4] # img of interest 4088x4088 + array_pad = self.saturate(im.copy()).array[ + 4:-4, 4:-4 + ] # img of interest 4088x4088 array_pad = np.pad( array_pad, [(4 + nbfe, 4 + nbfe), (4 + nbfe, 4 + nbfe)], mode="symmetric" ) # 4100x4100 array @@ -584,7 +649,6 @@ def bfe(self, im): for gj in range(num_grids): for gi in range(num_grids): - a_components_pad = np.zeros( ( 4, @@ -595,7 +659,13 @@ def bfe(self, im): ) ) # (4,5,5,sub_grid,sub_grid) a_components_pad = np.zeros( - (4, 2 * nbfe + 1, 2 * nbfe + 1, bin_size * n_sub + 2 * nbfe, bin_size * m_sub + 2 * nbfe) + ( + 4, + 2 * nbfe + 1, + 2 * nbfe + 1, + bin_size * n_sub + 2 * nbfe, + bin_size * m_sub + 2 * nbfe, + ) ) # (4,5,5,sub_grid,sub_grid) for comp in range(4): @@ -611,9 +681,7 @@ def bfe(self, im): gi * m_sub : (gi + 1) * m_sub, ] ) - ).dot( - t.T - ) # sub_grid*sub_grid + ).dot(t.T) # sub_grid*sub_grid a_components_pad[comp, j, i, :, :] = np.pad( tmp, [(nbfe, nbfe), (nbfe, nbfe)], mode="symmetric" ) @@ -635,14 +703,10 @@ def bfe(self, im): nbfe - dx : nbfe - dx + bin_size * m_sub, ] * array_pad[ - -dy - + nbfe - + gj * bin_size * n_sub : -dy + -dy + nbfe + gj * bin_size * n_sub : -dy + nbfe + (gj + 1) * bin_size * n_sub, - -dx - + nbfe - + gi * bin_size * m_sub : -dx + -dx + nbfe + gi * bin_size * m_sub : -dx + nbfe + (gi + 1) * bin_size * m_sub, ] @@ -657,12 +721,18 @@ def bfe(self, im): gi * bin_size * m_sub : (gi + 1) * bin_size * m_sub, ] *= 0.5 * ( array_pad[ - nbfe + gj * bin_size * n_sub : nbfe + (gj + 1) * bin_size * n_sub, - nbfe + gi * bin_size * m_sub : nbfe + (gi + 1) * bin_size * m_sub, + nbfe + gj * bin_size * n_sub : nbfe + + (gj + 1) * bin_size * n_sub, + nbfe + gi * bin_size * m_sub : nbfe + + (gi + 1) * bin_size * m_sub, ] + array_pad[ - dj + nbfe + gj * bin_size * n_sub : dj + nbfe + (gj + 1) * bin_size * n_sub, - di + nbfe + gi * bin_size * m_sub : di + nbfe + (gi + 1) * bin_size * m_sub, + dj + nbfe + gj * bin_size * n_sub : dj + + nbfe + + (gj + 1) * bin_size * n_sub, + di + nbfe + gi * bin_size * m_sub : di + + nbfe + + (gi + 1) * bin_size * m_sub, ] ) dQ_components[ @@ -671,12 +741,18 @@ def bfe(self, im): gi * bin_size * m_sub : (gi + 1) * bin_size * m_sub, ] *= 0.5 * ( array_pad[ - nbfe + gj * bin_size * n_sub : nbfe + (gj + 1) * bin_size * n_sub, - nbfe + gi * bin_size * m_sub : nbfe + (gi + 1) * bin_size * m_sub, + nbfe + gj * bin_size * n_sub : nbfe + + (gj + 1) * bin_size * n_sub, + nbfe + gi * bin_size * m_sub : nbfe + + (gi + 1) * bin_size * m_sub, ] + array_pad[ - dj + nbfe + gj * bin_size * n_sub : dj + nbfe + (gj + 1) * bin_size * n_sub, - di + nbfe + gi * bin_size * m_sub : di + nbfe + (gi + 1) * bin_size * m_sub, + dj + nbfe + gj * bin_size * n_sub : dj + + nbfe + + (gj + 1) * bin_size * n_sub, + di + nbfe + gi * bin_size * m_sub : di + + nbfe + + (gi + 1) * bin_size * m_sub, ] ) @@ -697,14 +773,20 @@ def get_eff_sky_bg(self, pointing, radec): radec : World coordinate position of image """ - sky_level = roman.getSkyLevel(pointing.bpass, world_pos=radec, date=pointing.date) + sky_level = roman.getSkyLevel( + pointing.bpass, world_pos=radec, date=pointing.date + ) sky_level *= (1.0 + roman.stray_light_fraction) * roman.pixel_scale**2 - sky_level = roman.getSkyLevel(pointing.bpass, world_pos=radec, date=pointing.date) + sky_level = roman.getSkyLevel( + pointing.bpass, world_pos=radec, date=pointing.date + ) sky_level *= (1.0 + roman.stray_light_fraction) * roman.pixel_scale**2 return sky_level - def translate_cvz(self, orig_radec, field_ra=9.5, field_dec=-44, cvz_ra=61.24, cvz_dec=-48.42): + def translate_cvz( + self, orig_radec, field_ra=9.5, field_dec=-44, cvz_ra=61.24, cvz_dec=-48.42 + ): ra = orig_radec.ra / galsim.degrees - field_ra dec = orig_radec.dec / galsim.degrees - field_dec @@ -737,10 +819,12 @@ def setup_sky(self, im, pointing, rng, rng_iter, force_cvz=False): self.dark_current_ = roman.dark_current * roman.exptime else: self.dark_current_ = ( - roman.dark_current * roman.exptime + self.df["DARK"][:, :].flatten() * roman.exptime + roman.dark_current * roman.exptime + + self.df["DARK"][:, :].flatten() * roman.exptime ) self.dark_current_ = ( - roman.dark_current * roman.exptime + self.df["DARK"][:, :].flatten() * roman.exptime + roman.dark_current * roman.exptime + + self.df["DARK"][:, :].flatten() * roman.exptime ) if self.df is None: self.gain = roman.gain @@ -753,9 +837,13 @@ def setup_sky(self, im, pointing, rng, rng_iter, force_cvz=False): radec = self.translate_cvz(pointing.radec) else: radec = pointing.radec - sky_level = roman.getSkyLevel(pointing.bpass, world_pos=radec, date=pointing.date) + sky_level = roman.getSkyLevel( + pointing.bpass, world_pos=radec, date=pointing.date + ) sky_level *= 1.0 + roman.stray_light_fraction - sky_level = roman.getSkyLevel(pointing.bpass, world_pos=radec, date=pointing.date) + sky_level = roman.getSkyLevel( + pointing.bpass, world_pos=radec, date=pointing.date + ) sky_level *= 1.0 + roman.stray_light_fraction # Make a image of the sky that takes into account the spatially variable pixel scale. Note # that makeSkyImage() takes a bit of time. If you do not care about the variable pixel @@ -775,7 +863,10 @@ def setup_sky(self, im, pointing, rng, rng_iter, force_cvz=False): # ~ (35, 3050, 0.008, 1.2E6) (e-/p). Hot pixels could be removed in further analysis using the # dq array. self.sky_mean = np.mean( - np.round((np.round(self.sky.array) + round(np.median(self.dark_current_))) / np.mean(self.gain)) + np.round( + (np.round(self.sky.array) + round(np.median(self.dark_current_))) + / np.mean(self.gain) + ) ) self.sky.addNoise(self.noise) @@ -811,7 +902,9 @@ def add_background(self, im): return im - def recip_failure(self, im, exptime=roman.exptime, alpha=roman.reciprocity_alpha, base_flux=1.0): + def recip_failure( + self, im, exptime=roman.exptime, alpha=roman.reciprocity_alpha, base_flux=1.0 + ): """ Introduce reciprocity failure to image. @@ -853,18 +946,27 @@ def dark_current(self, im): if self.df is None: self.im_dark = im.copy() dark_current_ = self.dark_current_ - dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(self.rng, dark_current_)) + dark_noise = galsim.DeviateNoise( + galsim.PoissonDeviate(self.rng, dark_current_) + ) im.addNoise(dark_noise) self.im_dark = im - self.im_dark else: - dark_current_ = self.dark_current_.clip(0) # opt for numpy random geneator instead for speed - self.im_dark = self.rng_np.poisson(dark_current_).reshape(im.array.shape).astype(im.dtype) + self.im_dark = ( + self.rng_np.poisson(dark_current_) + .reshape(im.array.shape) + .astype(im.dtype) + ) im.array[:, :] += self.im_dark - self.im_dark = self.rng_np.poisson(dark_current_).reshape(im.array.shape).astype(im.dtype) + self.im_dark = ( + self.rng_np.poisson(dark_current_) + .reshape(im.array.shape) + .astype(im.dtype) + ) im.array[:, :] += self.im_dark # NOTE: Sky level and dark current might appear like a constant background that can be @@ -955,14 +1057,26 @@ def add_persistence(self, im, pointing): dither_sca_array = np.loadtxt(self.params["dither_from_file"]).astype(int) # select adjacent exposures for the same sca (within 10*roman.exptime) - dither_list_selected = dither_sca_array[dither_sca_array[:, 1] == pointing.sca, 0] - dither_list_selected = dither_list_selected[np.abs(dither_list_selected - pointing.visit) < 10] - p_list = np.array([get_pointing(self.params, i, pointing.sca) for i in dither_list_selected]) + dither_list_selected = dither_sca_array[ + dither_sca_array[:, 1] == pointing.sca, 0 + ] + dither_list_selected = dither_list_selected[ + np.abs(dither_list_selected - pointing.visit) < 10 + ] + p_list = np.array( + [get_pointing(self.params, i, pointing.sca) for i in dither_list_selected] + ) dt_list = np.array([(pointing.date - p.date).total_seconds() for p in p_list]) p_pers = p_list[np.where((dt_list > 0) & (dt_list < roman.exptime * 10))] - dither_list_selected = dither_sca_array[dither_sca_array[:, 1] == pointing.sca, 0] - dither_list_selected = dither_list_selected[np.abs(dither_list_selected - pointing.visit) < 10] - p_list = np.array([get_pointing(self.params, i, pointing.sca) for i in dither_list_selected]) + dither_list_selected = dither_sca_array[ + dither_sca_array[:, 1] == pointing.sca, 0 + ] + dither_list_selected = dither_list_selected[ + np.abs(dither_list_selected - pointing.visit) < 10 + ] + p_list = np.array( + [get_pointing(self.params, i, pointing.sca) for i in dither_list_selected] + ) dt_list = np.array([(pointing.date - p.date).total_seconds() for p in p_list]) p_pers = p_list[np.where((dt_list > 0) & (dt_list < roman.exptime * 10))] @@ -973,7 +1087,9 @@ def add_persistence(self, im, pointing): pointing.date - p.date ).total_seconds() - roman.exptime / 2 # avg time since end of exposures self.params["output"]["file_name"]["items"] = [p.filter, p.visit, p.sca] - imfilename = ParseValue(self.params["output"], "file_name", self.params, str)[0] + imfilename = ParseValue( + self.params["output"], "file_name", self.params, str + )[0] fn = os.path.join(self.params["output"]["dir"], imfilename) # apply all the effects that occured before persistence on the previouse exposures @@ -985,11 +1101,16 @@ def add_persistence(self, im, pointing): x = x.clip(0) # remove negative stimulus - im.array[:, :] += galsim.roman.roman_detectors.fermi_linear(x.array, dt) * roman.exptime - im.array[:, :] += galsim.roman.roman_detectors.fermi_linear(x.array, dt) * roman.exptime + im.array[:, :] += ( + galsim.roman.roman_detectors.fermi_linear(x.array, dt) + * roman.exptime + ) + im.array[:, :] += ( + galsim.roman.roman_detectors.fermi_linear(x.array, dt) + * roman.exptime + ) else: - # setup parameters for persistence Q01 = self.df["PERSIST"].read_header()["Q01"] Q02 = self.df["PERSIST"].read_header()["Q02"] @@ -1005,10 +1126,12 @@ def add_persistence(self, im, pointing): pointing.date - p.date ).total_seconds() - roman.exptime / 2 # avg time since end of exposures fac_dt = ( - roman.exptime / 2.0 - ) / dt # linear time dependence (approximate until we get t1 and Delat t of the data) + (roman.exptime / 2.0) / dt + ) # linear time dependence (approximate until we get t1 and Delat t of the data) self.params["output"]["file_name"]["items"] = [p.filter, p.visit, p.sca] - imfilename = ParseValue(self.params["output"], "file_name", self.params, str)[0] + imfilename = ParseValue( + self.params["output"], "file_name", self.params, str + )[0] fn = os.path.join(self.params["output"]["dir"], imfilename) # apply all the effects that occured before persistence on the previouse exposures @@ -1023,7 +1146,7 @@ def add_persistence(self, im, pointing): # Do linear interpolation a = np.zeros(x.shape) - a += ((x < Q01)) * x / Q01 + a += (x < Q01) * x / Q01 a += ((x >= Q01) & (x < Q02)) * (Q02 - x) / (Q02 - Q01) im.array[:, :] += a * self.df["PERSIST"][0, :, :][0] * fac_dt @@ -1049,7 +1172,9 @@ def add_persistence(self, im, pointing): a = np.zeros(x.shape) a += ((x >= Q05) & (x < Q06)) * (x - Q05) / (Q06 - Q05) - a += ((x >= Q06)) * (x / Q06) ** alpha # avoid fractional power of negative values + a += (x >= Q06) * ( + x / Q06 + ) ** alpha # avoid fractional power of negative values im.array[:, :] += a * self.df["PERSIST"][5, :, :][0] * fac_dt return im @@ -1107,9 +1232,13 @@ def interpix_cap(self, im, kernel=roman.ipc_kernel): grid_size = 4096 // num_grids array_pad = im.array[4:-4, 4:-4] # it's an array instead of img - array_pad = np.pad(array_pad, [(5, 5), (5, 5)], mode="symmetric") # 4098x4098 array + array_pad = np.pad( + array_pad, [(5, 5), (5, 5)], mode="symmetric" + ) # 4098x4098 array array_pad = im.array[4:-4, 4:-4] # it's an array instead of img - array_pad = np.pad(array_pad, [(5, 5), (5, 5)], mode="symmetric") # 4098x4098 array + array_pad = np.pad( + array_pad, [(5, 5), (5, 5)], mode="symmetric" + ) # 4098x4098 array K = self.df["IPC"][:, :, :, :] # 3,3,512, 512 @@ -1126,13 +1255,16 @@ def interpix_cap(self, im, kernel=roman.ipc_kernel): for j in range(3): for i in range(3): tmp = (t.dot(K[j, i, :, :])).dot(t.T) # grid_sizexgrid_size - K_pad[j, i, :, :] = np.pad(tmp, [(1, 1), (1, 1)], mode="symmetric") + K_pad[j, i, :, :] = np.pad( + tmp, [(1, 1), (1, 1)], mode="symmetric" + ) tmp = (t.dot(K[j, i, :, :])).dot(t.T) # grid_sizexgrid_size - K_pad[j, i, :, :] = np.pad(tmp, [(1, 1), (1, 1)], mode="symmetric") + K_pad[j, i, :, :] = np.pad( + tmp, [(1, 1), (1, 1)], mode="symmetric" + ) for dy in range(-1, 2): for dx in range(-1, 2): - array_out[ gj * grid_size : (gj + 1) * grid_size, gi * grid_size : (gi + 1) * grid_size, @@ -1144,8 +1276,12 @@ def interpix_cap(self, im, kernel=roman.ipc_kernel): 1 - dx : 1 - dx + grid_size, ] * array_pad[ - 1 - dy + gj * grid_size : 1 - dy + (gj + 1) * grid_size, - 1 - dx + gi * grid_size : 1 - dx + (gi + 1) * grid_size, + 1 - dy + gj * grid_size : 1 + - dy + + (gj + 1) * grid_size, + 1 - dx + gi * grid_size : 1 + - dx + + (gi + 1) * grid_size, ] ) @@ -1174,12 +1310,16 @@ def add_read_noise(self, im): # use numpy random generator to draw 2-d noise map read_noise = self.df["READ"][2, :, :].flatten() # flattened 4096x4096 array self.im_read = ( - self.rng_np.normal(loc=0.0, scale=read_noise).reshape(im.array.shape).astype(im.dtype) + self.rng_np.normal(loc=0.0, scale=read_noise) + .reshape(im.array.shape) + .astype(im.dtype) ) im.array[:, :] += self.im_read read_noise = self.df["READ"][2, :, :].flatten() # flattened 4096x4096 array self.im_read = ( - self.rng_np.normal(loc=0.0, scale=read_noise).reshape(im.array.shape).astype(im.dtype) + self.rng_np.normal(loc=0.0, scale=read_noise) + .reshape(im.array.shape) + .astype(im.dtype) ) im.array[:, :] += self.im_read @@ -1258,7 +1398,6 @@ def finalize_sky_im(self, im, pointing): im = self.e_to_ADU(im) im.quantize() else: - bound_pad = galsim.BoundsI(xmin=1, ymin=1, xmax=4096, ymax=4096) im_pad = galsim.Image(bound_pad) im_pad.array[4:-4, 4:-4] = im.array[:, :] diff --git a/roman_imsim/noise.py b/roman_imsim/noise.py index 7bc0df20..5ff071bc 100644 --- a/roman_imsim/noise.py +++ b/roman_imsim/noise.py @@ -46,21 +46,21 @@ def addNoise(self, config, base, image, rng, current_var, draw_method, logger): "gain": dict, } - params, safe = galsim.config.GetAllParams(config, base, req={}, opt=opt, ignore=[]) + params, safe = galsim.config.GetAllParams( + config, base, req={}, opt=opt, ignore=[] + ) mjd = params.get("mjd", None) stray_light = params.get("stray_light", False) thermal_background = params.get("thermal_background", False) reciprocity_failure = params.get("reciprocity_failure", False) # dark_current = params.get("dark_current", False) - dark_current = params.get("dark_current", {'turn_on': False, 'use_crds': False}) - nonlinearity = params.get("nonlinearity", {'turn_on': False, 'use_crds': False}) - ipc = params.get("ipc", {'turn_on': False, 'use_crds': False}) - read_noise = params.get("read_noise", {'turn_on': False, 'use_crds': False}) + dark_current = params.get("dark_current", {"turn_on": False, "use_crds": False}) + nonlinearity = params.get("nonlinearity", {"turn_on": False, "use_crds": False}) + ipc = params.get("ipc", {"turn_on": False, "use_crds": False}) + read_noise = params.get("read_noise", {"turn_on": False, "use_crds": False}) sky_subtract = params.get("sky_subtract", True) - gain = params.get('gain', {'turn_on': False, 'use_crds': False}) - - use_crds = params.get("use_crds", False) + gain = params.get("gain", {"turn_on": False, "use_crds": False}) base["current_noise_image"] = base["current_image"] wcs = base["wcs"] @@ -68,16 +68,22 @@ def addNoise(self, config, base, image, rng, current_var, draw_method, logger): filter_name = bp.name exptime, _ = galsim.config.ParseValue(base["image"], "exptime", base, float) date = Time(mjd, format="mjd").to_datetime() if mjd is not None else None - logger.info("image %d: Start RomanSCA detector effects", base.get("image_num", 0)) + logger.info( + "image %d: Start RomanSCA detector effects", base.get("image_num", 0) + ) # Things that will eventually be subtracted (if sky_subtract) will have their expectation # value added to sky_image. So technically, this includes things that aren't just sky. # E.g. includes dark_current and thermal backgrounds. sky_image = image.copy() - sky_level = models.backgrounds.getSkyLevel(bp, world_pos=wcs.toWorld(image.true_center), date=date) + sky_level = models.backgrounds.getSkyLevel( + bp, world_pos=wcs.toWorld(image.true_center), date=date + ) logger.debug("Adding sky_level = %s", sky_level) if stray_light: - logger.debug("Stray light fraction = %s", models.parameters.stray_light_fraction) + logger.debug( + "Stray light fraction = %s", models.parameters.stray_light_fraction + ) sky_level *= 1.0 + models.parameters.stray_light_fraction wcs.makeSkyImage(sky_image, sky_level) @@ -119,33 +125,33 @@ def addNoise(self, config, base, image, rng, current_var, draw_method, logger): logger.debug("Applying reciprocity failure") models.nonlinearity.addReciprocityFailure(img=image) - if dark_current['turn_on']: + if dark_current["turn_on"]: logger.debug("Adding dark current: %s") # print(f'for drak_current, use_crds = {dark_current["use_crds"]}') - dc = models.DarkCurrent(usecrds=dark_current['use_crds']) + dc = models.DarkCurrent(usecrds=dark_current["use_crds"]) dc.apply(img=image, exptime=exptime) - if nonlinearity['turn_on']: + if nonlinearity["turn_on"]: logger.debug("Applying classical nonlinearity") - non_linear = models.Nonlinearity(usecrds=nonlinearity['use_crds']) + non_linear = models.Nonlinearity(usecrds=nonlinearity["use_crds"]) non_linear.apply(img=image, electrons=False) # Mosby and Rauscher say there are two read noises. One happens before IPC, the other # one after. # TODO: Add read_noise1 - if ipc['turn_on']: + if ipc["turn_on"]: logger.debug("Applying IPC") - IPC = models.IPC(usecrds=ipc['use_crds']) + IPC = models.IPC(usecrds=ipc["use_crds"]) IPC.apply(img=image) - if read_noise['turn_on']: + if read_noise["turn_on"]: logger.debug("Adding read noise") - rn = models.ReadNoise(usecrds=read_noise['use_crds']) + rn = models.ReadNoise(usecrds=read_noise["use_crds"]) rn.apply(img=image) - if gain['turn_on']: + if gain["turn_on"]: logger.debug("Applying gain") - gain_model = models.Gain(usecrds=gain['use_crds']) + gain_model = models.Gain(usecrds=gain["use_crds"]) gain_model.apply(img=image) # Make integer ADU now. @@ -153,7 +159,7 @@ def addNoise(self, config, base, image, rng, current_var, draw_method, logger): if sky_subtract: logger.debug("Subtracting sky image") - if gain['turn_on']: + if gain["turn_on"]: logger.debug("Applying gain") gain_model.apply(img=sky_image) sky_image.quantize() diff --git a/roman_imsim/obseq.py b/roman_imsim/obseq.py index 79c05263..9bf4b3c2 100644 --- a/roman_imsim/obseq.py +++ b/roman_imsim/obseq.py @@ -26,11 +26,17 @@ def __init__(self, file_name, visit, SCA, logger=None): def read_obseq(self): """Read visit info from the obseq file.""" if self.file_name is None: - raise ValueError("No obseq filename provided, trying to build from config information.") + raise ValueError( + "No obseq filename provided, trying to build from config information." + ) if self.visit is None: - raise ValueError("The visit must be set when reading visit info from an obseq file.") + raise ValueError( + "The visit must be set when reading visit info from an obseq file." + ) - self.logger.warning("Reading info from obseq file %s for visit %s", self.file_name, self.visit) + self.logger.warning( + "Reading info from obseq file %s for visit %s", self.file_name, self.visit + ) ob = fio.FITS(self.file_name)[-1][self.visit] @@ -62,5 +68,9 @@ def ObSeqData(config, base, value_type): return val, safe -RegisterInputType("obseq_data", InputLoader(ObSeqDataLoader, file_scope=True, takes_logger=True)) -RegisterValueType("ObSeqData", ObSeqData, [float, int, str, Angle], input_type="obseq_data") +RegisterInputType( + "obseq_data", InputLoader(ObSeqDataLoader, file_scope=True, takes_logger=True) +) +RegisterValueType( + "ObSeqData", ObSeqData, [float, int, str, Angle], input_type="obseq_data" +) diff --git a/roman_imsim/psf.py b/roman_imsim/psf.py index e859a44c..f4ba979b 100644 --- a/roman_imsim/psf.py +++ b/roman_imsim/psf.py @@ -125,7 +125,9 @@ def _parse_pupil_bin(self, pupil_bin): else: return pupil_bin - def _psf_call(self, SCA, bpass, SCA_pos, WCS, pupil_bin, n_waves, logger, extra_aberrations): + def _psf_call( + self, SCA, bpass, SCA_pos, WCS, pupil_bin, n_waves, logger, extra_aberrations + ): if pupil_bin == 8: psf = models.psf_utils.getPSF( @@ -306,9 +308,9 @@ def getPSF(self, pupil_bin, pos): wlu = (models.parameters.n_pix - pos.x) * (pos.y - 1) wul = (pos.x - 1) * (models.parameters.n_pix - pos.y) wuu = (pos.x - 1) * (pos.y - 1) - return (wll * psf["ll"] + wlu * psf["lu"] + wul * psf["ul"] + wuu * psf["uu"]) / ( - (models.parameters.n_pix - 1) * (models.parameters.n_pix - 1) - ) + return ( + wll * psf["ll"] + wlu * psf["lu"] + wul * psf["ul"] + wuu * psf["uu"] + ) / ((models.parameters.n_pix - 1) * (models.parameters.n_pix - 1)) def RegisterPSFInterpolatorType(interp_type, builder, input_type=None): @@ -346,7 +348,16 @@ def getKwargs(self, config, base, logger): } ignore = ["extra_aberrations"] - kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt, ignore=ignore) + # If SCA is in base, then don't require it in the config file. + # (Presumably because using Roman image type, which sets it there for convenience.) + if "SCA" in base: + opt["SCA"] = int + else: + req["SCA"] = int + + kwargs, safe = galsim.config.GetAllParams( + config, base, req=req, opt=opt, ignore=ignore + ) kwargs["extra_aberrations"] = galsim.config.ParseAberrations( "extra_aberrations", config, base, "RomanPSF" diff --git a/roman_imsim/sca.py b/roman_imsim/sca.py index 7524a21a..46559cc2 100644 --- a/roman_imsim/sca.py +++ b/roman_imsim/sca.py @@ -52,7 +52,9 @@ def setup(self, config, base, image_num, obj_num, ignore, logger): "draw_method": str, "use_fft_bright": bool, } - params = galsim.config.GetAllParams(config, base, req=req, opt=opt, ignore=ignore + extra_ignore)[0] + params = galsim.config.GetAllParams( + config, base, req=req, opt=opt, ignore=ignore + extra_ignore + )[0] self.sca = params["SCA"] base["SCA"] = self.sca @@ -61,11 +63,15 @@ def setup(self, config, base, image_num, obj_num, ignore, logger): self.exptime = params["exptime"] # If draw_method isn't in image field, it may be in stamp. Check. - self.draw_method = params.get("draw_method", base.get("stamp", {}).get("draw_method", "phot")) + self.draw_method = params.get( + "draw_method", base.get("stamp", {}).get("draw_method", "phot") + ) # If user hasn't overridden the bandpass to use, get the standard one. if "bandpass" not in config: - base["bandpass"] = galsim.config.BuildBandpass(base["image"], "bandpass", base, logger=logger) + base["bandpass"] = galsim.config.BuildBandpass( + base["image"], "bandpass", base, logger=logger + ) return models.parameters.n_pix, models.parameters.n_pix @@ -98,9 +104,13 @@ def buildImage(self, config, base, image_num, obj_num, logger): full_image.header["VERSION"] = "unknown" full_image.header["EXPTIME"] = self.exptime full_image.header["MJD-OBS"] = self.mjd - full_image.header["DATE-OBS"] = Time(self.mjd, format="mjd").datetime.isoformat() + full_image.header["DATE-OBS"] = Time( + self.mjd, format="mjd" + ).datetime.isoformat() full_image.header["FILTER"] = self.filter - full_image.header["ZPTMAG"] = 2.5 * np.log10(self.exptime * models.parameters.collecting_area) + full_image.header["ZPTMAG"] = 2.5 * np.log10( + self.exptime * models.parameters.collecting_area + ) base["current_image"] = full_image @@ -151,7 +161,9 @@ def buildImage(self, config, base, image_num, obj_num, logger): # avoid an error. But this isn't covered in the imsim test suite. continue - logger.debug("image %d: full bounds = %s", image_num, str(full_image.bounds)) + logger.debug( + "image %d: full bounds = %s", image_num, str(full_image.bounds) + ) logger.debug( "image %d: stamp %d bounds = %s", image_num, diff --git a/roman_imsim/scafile.py b/roman_imsim/scafile.py index 75574a64..a9315037 100644 --- a/roman_imsim/scafile.py +++ b/roman_imsim/scafile.py @@ -27,7 +27,9 @@ def setup(self, config, base, file_num, logger): logger.debug("file %d: seed = %d", file_num, seed) if "exptime" in config: - base["exptime"] = galsim.config.ParseValue(config, "exptime", base, float)[0] + base["exptime"] = galsim.config.ParseValue(config, "exptime", base, float)[ + 0 + ] else: base["exptime"] = models.parameters.exptime @@ -92,7 +94,9 @@ def buildImages(self, config, base, file_num, image_num, obj_num, ignore, logger if cosmic_ray_rate > 0: cosmic_ray_catalog = params.get("cosmic_ray_catalog", None) if cosmic_ray_catalog is None: - cosmic_ray_catalog = os.path.join(data_dir, "cosmic_rays_itl_2017.fits.gz") + cosmic_ray_catalog = os.path.join( + data_dir, "cosmic_rays_itl_2017.fits.gz" + ) if not os.path.isfile(cosmic_ray_catalog): raise FileNotFoundError(f"{cosmic_ray_catalog} not found") diff --git a/roman_imsim/skycat.py b/roman_imsim/skycat.py index bf31eca1..b2de68e0 100644 --- a/roman_imsim/skycat.py +++ b/roman_imsim/skycat.py @@ -79,7 +79,9 @@ def __init__( if obj_types is not None: self.logger.warning(f"Object types restricted to {obj_types}") - self.sca_center = wcs.toWorld(galsim.PositionD(self.xsize / 2.0, self.ysize / 2.0)) + self.sca_center = wcs.toWorld( + galsim.PositionD(self.xsize / 2.0, self.ysize / 2.0) + ) self._objects = None # import os, psutil @@ -111,10 +113,14 @@ def objects(self): vertices = [] for x, y in corners: sky_coord = self.wcs.toWorld(galsim.PositionD(x, y)) - vertices.append((sky_coord.ra / galsim.degrees, sky_coord.dec / galsim.degrees)) + vertices.append( + (sky_coord.ra / galsim.degrees, sky_coord.dec / galsim.degrees) + ) region = PolygonalRegion(vertices) sky_cat = skyCatalogs.open_catalog(self.file_name) - self._objects = sky_cat.get_objects_by_region(region, obj_type_set=self.obj_types, mjd=self.mjd) + self._objects = sky_cat.get_objects_by_region( + region, obj_type_set=self.obj_types, mjd=self.mjd + ) if not self._objects: self.logger.warning("No objects found on image.") # import os, psutil @@ -184,7 +190,9 @@ def getObj(self, index, gsparams=None, rng=None, exptime=30): # Compute the flux or get the cached value. flux = ( - skycat_obj.get_roman_flux(self.bandpass.name, mjd=self.mjd) * self.exptime * models.parameters.collecting_area + skycat_obj.get_roman_flux(self.bandpass.name, mjd=self.mjd) + * self.exptime + * models.parameters.collecting_area ) if np.isnan(flux): return None @@ -207,7 +215,9 @@ def getObj(self, index, gsparams=None, rng=None, exptime=30): gs_obj_list = [] for component in gsobjs: if faint: - gsobjs[component] = gsobjs[component].evaluateAtWavelength(self.bandpass) + gsobjs[component] = gsobjs[component].evaluateAtWavelength( + self.bandpass + ) gs_obj_list.append(gsobjs[component] * self._trivial_sed) else: if component in seds: @@ -234,7 +244,9 @@ def getObj(self, index, gsparams=None, rng=None, exptime=30): gs_object.flux = flux # Get the object type - if (skycat_obj.object_type == "diffsky_galaxy") | (skycat_obj.object_type == "galaxy"): + if (skycat_obj.object_type == "diffsky_galaxy") | ( + skycat_obj.object_type == "galaxy" + ): gs_object.object_type = "galaxy" if skycat_obj.object_type in {"star", "gaia_star"}: gs_object.object_type = "star" @@ -264,7 +276,9 @@ def getKwargs(self, config, base, logger): kwargs["logger"] = logger if "bandpass" not in config: - base["bandpass"] = galsim.config.BuildBandpass(base["image"], "bandpass", base, logger=logger)[0] + base["bandpass"] = galsim.config.BuildBandpass( + base["image"], "bandpass", base, logger=logger + )[0] kwargs["bandpass"] = base["bandpass"] # Sky catalog object lists are created per CCD, so they are @@ -289,11 +303,11 @@ def SkyCatObj(config, base, ignore, gsparams, logger): message = ( "skyCatalogs selection and SCA center do not agree: \n" "skycat.sca_center: " - f"{sca_center.ra/galsim.degrees:.5f}, " - f"{sca_center.dec/galsim.degrees:.5f}\n" + f"{sca_center.ra / galsim.degrees:.5f}, " + f"{sca_center.dec / galsim.degrees:.5f}\n" "world_center: " - f"{world_center.ra/galsim.degrees:.5f}, " - f"{world_center.dec/galsim.degrees:.5f} \n" + f"{world_center.ra / galsim.degrees:.5f}, " + f"{world_center.dec / galsim.degrees:.5f} \n" f"Separation: {sep:.2e} arcsec" ) raise RuntimeError(message) @@ -338,7 +352,9 @@ def SkyCatWorldPos(config, base, value_type): RegisterInputType("sky_catalog", SkyCatalogLoader(SkyCatalogInterface, has_nobj=True)) RegisterObjectType("SkyCatObj", SkyCatObj, input_type="sky_catalog") -RegisterValueType("SkyCatWorldPos", SkyCatWorldPos, [galsim.CelestialCoord], input_type="sky_catalog") +RegisterValueType( + "SkyCatWorldPos", SkyCatWorldPos, [galsim.CelestialCoord], input_type="sky_catalog" +) # This class was modified from https://github.com/LSSTDESC/imSim/. License info follows: diff --git a/roman_imsim/stamp.py b/roman_imsim/stamp.py index f491fe06..9e57c67b 100644 --- a/roman_imsim/stamp.py +++ b/roman_imsim/stamp.py @@ -83,7 +83,9 @@ def setup(self, config, base, xsize, ysize, ignore, logger): else: gal_achrom = gal.evaluateAtWavelength(bandpass.effective_wavelength) - if hasattr(gal_achrom, "original") and isinstance(gal_achrom.original, galsim.DeltaFunction): + if hasattr(gal_achrom, "original") and isinstance( + gal_achrom.original, galsim.DeltaFunction + ): # For bright stars, set the following stamp size limits if self.flux < 1e6: image_size = 500 @@ -105,12 +107,16 @@ def setup(self, config, base, xsize, ysize, ignore, logger): # print('stamp setup3',process.memory_info().rss) base["pupil_bin"] = self.pupil_bin logger.info("Object flux is %d", self.flux) - logger.info("Object %d will use stamp size = %s", base.get("obj_num", 0), image_size) + logger.info( + "Object %d will use stamp size = %s", base.get("obj_num", 0), image_size + ) # Determine where this object is going to go: # This is the same as what the base StampBuilder does: if "image_pos" in config: - image_pos = galsim.config.ParseValue(config, "image_pos", base, galsim.PositionD)[0] + image_pos = galsim.config.ParseValue( + config, "image_pos", base, galsim.PositionD + )[0] else: image_pos = None @@ -121,6 +127,31 @@ def setup(self, config, base, xsize, ysize, ignore, logger): return image_size, image_size, image_pos, world_pos + def buildPSF(self, config, base, gsparams, logger): + """Build the PSF object. + + For the Basic stamp type, this builds a PSF from the base['psf'] dict, if present, + else returns None. + + Parameters: + config: The configuration dict for the stamp field. + base: The base configuration dict. + gsparams: A dict of kwargs to use for a GSParams. More may be added to this + list by the galaxy object. + logger: A logger object to log progress. + + Returns: + the PSF + """ + if base.get("psf", {}).get("type", "roman_psf") != "roman_psf": + return galsim.config.BuildGSObject( + base, "psf", gsparams=gsparams, logger=logger + )[0] + + roman_psf = galsim.config.GetInputObj("roman_psf", config, base, "buildPSF") + psf = roman_psf.getPSF(self.pupil_bin, base["image_pos"]) + return psf + def getDrawMethod(self, config, base, logger): """Determine the draw method to use. @@ -133,7 +164,9 @@ def getDrawMethod(self, config, base, logger): method = galsim.config.ParseValue(config, "draw_method", base, str)[0] self.use_fft_bright = False if "use_fft_bright" in config: - self.use_fft_bright = galsim.config.ParseValue(config, "use_fft_bright", base, bool)[0] + self.use_fft_bright = galsim.config.ParseValue( + config, "use_fft_bright", base, bool + )[0] if method not in galsim.config.valid_draw_methods: raise galsim.GalSimConfigValueError( @@ -145,12 +178,16 @@ def getDrawMethod(self, config, base, logger): logger.info("Auto -> Use FFT drawing for object %d.", base["obj_num"]) return "fft" else: - logger.info("Auto -> Use photon shooting for object %d.", base["obj_num"]) + logger.info( + "Auto -> Use photon shooting for object %d.", base["obj_num"] + ) return "phot" else: # If user sets something specific for the method, rather than auto, # then respect their wishes. - logger.info("Use specified method=%s for object %d.", method, base["obj_num"]) + logger.info( + "Use specified method=%s for object %d.", method, base["obj_num"] + ) return method @classmethod @@ -163,11 +200,16 @@ def _fix_seds_24(cls, prof, bandpass): sed = prof.sed # TODO: This bit should probably be ported back to Galsim. # Something like sed.make_tabulated() - if not isinstance(sed._spec, galsim.LookupTable) or sed._spec.interpolant != "linear": + if ( + not isinstance(sed._spec, galsim.LookupTable) + or sed._spec.interpolant != "linear" + ): # Workaround for https://github.com/GalSim-developers/GalSim/issues/1228 f = np.broadcast_to(sed(wave_list), wave_list.shape) new_spec = galsim.LookupTable(wave_list, f, interpolant="linear") - new_sed = galsim.SED(new_spec, "nm", "fphotons" if sed.spectral else "1") + new_sed = galsim.SED( + new_spec, "nm", "fphotons" if sed.spectral else "1" + ) prof.sed = new_sed # Also recurse onto any components. @@ -252,7 +294,9 @@ def draw(self, prof, image, method, offset, config, base, logger): # print('stamp draw3b ',process.memory_info().rss) if not faint and "photon_ops" in config: - photon_ops = galsim.config.BuildPhotonOps(config, "photon_ops", base, logger) + photon_ops = galsim.config.BuildPhotonOps( + config, "photon_ops", base, logger + ) else: photon_ops = [] @@ -290,7 +334,9 @@ def draw(self, prof, image, method, offset, config, base, logger): if not faint and config.get("fft_photon_ops"): kwargs.update( { - "photon_ops": galsim.config.BuildPhotonOps(config, "fft_photon_ops", base, logger), + "photon_ops": galsim.config.BuildPhotonOps( + config, "fft_photon_ops", base, logger + ), "maxN": maxN, "rng": self.rng, "n_subsample": 1, @@ -308,7 +354,9 @@ def draw(self, prof, image, method, offset, config, base, logger): # and raise the error. logger.error("Caught error trying to draw using FFT:") logger.error("%s", e) - logger.error("You may need to add a gsparams field with maximum_fft_size to") + logger.error( + "You may need to add a gsparams field with maximum_fft_size to" + ) logger.error("either the psf or gal field to allow larger FFTs.") logger.info("prof = %r", prof) logger.info("fft_image = %s", fft_image) diff --git a/roman_imsim/utils.py b/roman_imsim/utils.py index 8b6d0219..788d4bbe 100644 --- a/roman_imsim/utils.py +++ b/roman_imsim/utils.py @@ -11,7 +11,9 @@ class roman_utils(object): Class to contain a variety of helper routines to work with the simulation data. """ - def __init__(self, config_file, visit=None, sca=None, image_name=None, setup_skycat=False): + def __init__( + self, config_file, visit=None, sca=None, image_name=None, setup_skycat=False + ): """ Setup information about a simulated Roman image. Parameters: @@ -36,10 +38,16 @@ def __init__(self, config_file, visit=None, sca=None, image_name=None, setup_sky self.skycat = galsim.config.GetInputObj( "sky_catalog", config["input"]["sky_catalog"], config, "sky_catalog" ) - self.PSF = galsim.config.GetInputObj("roman_psf", config["input"]["roman_psf"], config, "roman_psf") + self.PSF = galsim.config.GetInputObj( + "roman_psf", config["input"]["roman_psf"], config, "roman_psf" + ) self.wcs = galsim.config.BuildWCS(config["image"], "wcs", config) - self.bpass = galsim.config.BuildBandpass(config["image"], "bandpass", config, None)[0] - self.photon_ops = galsim.config.BuildPhotonOps(config["stamp"], "photon_ops", config, None) + self.bpass = galsim.config.BuildBandpass( + config["image"], "bandpass", config, None + )[0] + self.photon_ops = galsim.config.BuildPhotonOps( + config["stamp"], "photon_ops", config, None + ) self.rng = galsim.config.GetRNG(config, config["image"], None, "roman_sca") def check_input(self, visit, sca, image_name): @@ -54,7 +62,9 @@ def check_input(self, visit, sca, image_name): tmp = np.array(image_name[start:end].split("_")).astype(int) return tmp[0], tmp[1] if (visit is None) | (sca is None): - raise ValueError("Insufficient information to construct visit info - all inputs are None.") + raise ValueError( + "Insufficient information to construct visit info - all inputs are None." + ) return visit, sca def getPSF(self, x=None, y=None, pupil_bin=8): @@ -72,9 +82,19 @@ def getPSF(self, x=None, y=None, pupil_bin=8): raise ValueError( "x,y position for pupil_bin values other than 8 not supported. Using SCA center." ) - return self.PSF.getPSF(pupil_bin, galsim.PositionD(models.parameters.n_pix / 2, models.parameters.n_pix / 2)) + return self.PSF.getPSF( + pupil_bin, + galsim.PositionD( + models.parameters.n_pix / 2, models.parameters.n_pix / 2 + ), + ) if (x is None) | (y is None): - return self.PSF.getPSF(8, galsim.PositionD(models.parameters.n_pix / 2, models.parameters.n_pix / 2)) + return self.PSF.getPSF( + 8, + galsim.PositionD( + models.parameters.n_pix / 2, models.parameters.n_pix / 2 + ), + ) return self.PSF.getPSF(8, galsim.PositionD(x, y)) def getWCS(self): @@ -142,7 +162,9 @@ def getPSF_Image( dvdx=local_wcs.dvdx / oversampling_factor, dvdy=local_wcs.dvdy / oversampling_factor, ) - stamp = galsim.Image(stamp_size * oversampling_factor, stamp_size * oversampling_factor, wcs=wcs) + stamp = galsim.Image( + stamp_size * oversampling_factor, stamp_size * oversampling_factor, wcs=wcs + ) if not include_photonOps: psf = galsim.Convolve(point, self.getPSF(x, y, pupil_bin)) return psf.drawImage(self.bpass, image=stamp, wcs=wcs, method=method) @@ -157,8 +179,12 @@ def getPSF_Image( dy=j * scale / oversampling_factor, ) for i, j in product( - np.arange((1 - oversampling_factor) / 2, (1 + oversampling_factor) / 2), - np.arange((1 - oversampling_factor) / 2, (1 + oversampling_factor) / 2), + np.arange( + (1 - oversampling_factor) / 2, (1 + oversampling_factor) / 2 + ), + np.arange( + (1 - oversampling_factor) / 2, (1 + oversampling_factor) / 2 + ), ) ] ) diff --git a/roman_imsim/wcs.py b/roman_imsim/wcs.py index 9e7951b0..964e4d69 100644 --- a/roman_imsim/wcs.py +++ b/roman_imsim/wcs.py @@ -7,7 +7,6 @@ class RomanWCS(WCSBuilder): - def buildWCS(self, config, base, logger): req = { diff --git a/tests/test_end_to_end.py b/tests/test_end_to_end.py index 1d2fd66f..e75c0314 100644 --- a/tests/test_end_to_end.py +++ b/tests/test_end_to_end.py @@ -22,10 +22,18 @@ def test_compare_truth(): reference_truth = ( - Path(TEST_REPO) / REFERENCE_DIR / "RomanWAS_new" / "truth" / "Roman_WAS_index_J129_12909_4.txt" + Path(TEST_REPO) + / REFERENCE_DIR + / "RomanWAS_new" + / "truth" + / "Roman_WAS_index_J129_12909_4.txt" ) output_truth = ( - Path(TEST_REPO) / OUTPUT_DIR / "RomanWAS_new" / "truth" / "Roman_WAS_index_J129_12909_4.txt" + Path(TEST_REPO) + / OUTPUT_DIR + / "RomanWAS_new" + / "truth" + / "Roman_WAS_index_J129_12909_4.txt" ) output_table = ascii.read(output_truth) diff --git a/tests/test_imports.py b/tests/test_imports.py index 0d1bc016..2511b8e9 100644 --- a/tests/test_imports.py +++ b/tests/test_imports.py @@ -9,10 +9,14 @@ def test_all_modules_import(): package = roman_imsim errors = [] - for loader, module_name, is_pkg in pkgutil.walk_packages(package.__path__, package.__name__ + "."): + for loader, module_name, is_pkg in pkgutil.walk_packages( + package.__path__, package.__name__ + "." + ): try: importlib.import_module(module_name) except Exception as e: errors.append((module_name, repr(e))) - assert not errors, "Failed to import modules:\n" + "\n".join(f"{m}: {e}" for m, e in errors) + assert not errors, "Failed to import modules:\n" + "\n".join( + f"{m}: {e}" for m, e in errors + ) From 84935e889cb7394e6f78dc057c19cbc2ec7631d5 Mon Sep 17 00:00:00 2001 From: Yuedong Fang Date: Thu, 30 Apr 2026 17:01:45 -0400 Subject: [PATCH 07/11] formatted --- roman_imsim/_templates.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/roman_imsim/_templates.py b/roman_imsim/_templates.py index 7a2931be..9ee09d0b 100644 --- a/roman_imsim/_templates.py +++ b/roman_imsim/_templates.py @@ -4,4 +4,6 @@ from ._meta_data import config_dir -galsim.config.RegisterTemplate("roman_imsim_default", os.path.join(config_dir, "default.yaml")) +galsim.config.RegisterTemplate( + "roman_imsim_default", os.path.join(config_dir, "default.yaml") +) From 9cbf3e34a94746a762d56349922d19cd176c1bcc Mon Sep 17 00:00:00 2001 From: Yuedong Fang Date: Thu, 30 Apr 2026 17:57:23 -0400 Subject: [PATCH 08/11] code formatted with black --- roman_imsim/_templates.py | 4 +- roman_imsim/bandpass.py | 4 +- roman_imsim/detector_physics.py | 241 +++++++++----------------------- roman_imsim/noise.py | 16 +-- roman_imsim/obseq.py | 20 +-- roman_imsim/psf.py | 16 +-- roman_imsim/sca.py | 24 +--- roman_imsim/scafile.py | 8 +- roman_imsim/skycat.py | 28 +--- roman_imsim/stamp.py | 49 ++----- roman_imsim/utils.py | 40 ++---- tests/test_end_to_end.py | 12 +- tests/test_imports.py | 8 +- 13 files changed, 125 insertions(+), 345 deletions(-) diff --git a/roman_imsim/_templates.py b/roman_imsim/_templates.py index 9ee09d0b..7a2931be 100644 --- a/roman_imsim/_templates.py +++ b/roman_imsim/_templates.py @@ -4,6 +4,4 @@ from ._meta_data import config_dir -galsim.config.RegisterTemplate( - "roman_imsim_default", os.path.join(config_dir, "default.yaml") -) +galsim.config.RegisterTemplate("roman_imsim_default", os.path.join(config_dir, "default.yaml")) diff --git a/roman_imsim/bandpass.py b/roman_imsim/bandpass.py index aec401c7..19b657b5 100644 --- a/roman_imsim/bandpass.py +++ b/roman_imsim/bandpass.py @@ -27,9 +27,7 @@ def buildBandpass(self, config, base, logger): name = kwargs["name"] if "SCA" in kwargs: - bandpass = models.bandpass.getBandpasses(red_limit=2000, sca=kwargs["SCA"])[ - name - ] + bandpass = models.bandpass.getBandpasses(red_limit=2000, sca=kwargs["SCA"])[name] else: bandpass = models.bandpass.getBandpasses(red_limit=2000)[name] diff --git a/roman_imsim/detector_physics.py b/roman_imsim/detector_physics.py index 3e29c366..ccfe4323 100644 --- a/roman_imsim/detector_physics.py +++ b/roman_imsim/detector_physics.py @@ -75,30 +75,22 @@ def __init__(self, params, visit, SCA): self.exptime = obseq_data.ob["exptime"] self.bpass = roman.getBandpasses()[self.filter] self.WCS = roman.getWCS( - world_pos=galsim.CelestialCoord( - ra=obseq_data.ob["ra"], dec=obseq_data.ob["dec"] - ), + world_pos=galsim.CelestialCoord(ra=obseq_data.ob["ra"], dec=obseq_data.ob["dec"]), PA=obseq_data.ob["pa"], date=self.date, SCAs=self.sca, PA_is_FPA=True, )[self.sca] - self.radec = self.WCS.toWorld( - galsim.PositionI(int(roman.n_pix / 2), int(roman.n_pix / 2)) - ) + self.radec = self.WCS.toWorld(galsim.PositionI(int(roman.n_pix / 2), int(roman.n_pix / 2))) self.WCS = roman.getWCS( - world_pos=galsim.CelestialCoord( - ra=obseq_data.ob["ra"], dec=obseq_data.ob["dec"] - ), + world_pos=galsim.CelestialCoord(ra=obseq_data.ob["ra"], dec=obseq_data.ob["dec"]), PA=obseq_data.ob["pa"], date=self.date, SCAs=self.sca, PA_is_FPA=True, )[self.sca] - self.radec = self.WCS.toWorld( - galsim.PositionI(int(roman.n_pix / 2), int(roman.n_pix / 2)) - ) + self.radec = self.WCS.toWorld(galsim.PositionI(int(roman.n_pix / 2), int(roman.n_pix / 2))) class modify_image(object): @@ -106,9 +98,7 @@ class modify_image(object): Class to simulate non-idealities and noise of roman detector images. """ - def __init__( - self, params, visit, sca, dither_from_file, sca_filepath=None, use_galsim=False - ): + def __init__(self, params, visit, sca, dither_from_file, sca_filepath=None, use_galsim=False): """ Set up noise properties of image @@ -129,13 +119,9 @@ def __init__( # Load sca file if applicable if sca_filepath is not None: - self.df = fio.FITS( - sca_filepath + "/" + sca_number_to_file[self.pointing.sca] - ) + self.df = fio.FITS(sca_filepath + "/" + sca_number_to_file[self.pointing.sca]) print("------- Using SCA files --------") - self.df = fio.FITS( - sca_filepath + "/" + sca_number_to_file[self.pointing.sca] - ) + self.df = fio.FITS(sca_filepath + "/" + sca_number_to_file[self.pointing.sca]) print("------- Using SCA files --------") else: self.df = None @@ -146,15 +132,9 @@ def __init__( old_filename = os.path.join(self.params["output"]["dir"], imfilename) if not os.path.exists( - self.params["output"]["dir"].replace( - "truth", self.get_path_name(use_galsim=use_galsim) - ) + self.params["output"]["dir"].replace("truth", self.get_path_name(use_galsim=use_galsim)) ): - os.mkdir( - self.params["output"]["dir"].replace( - "truth", self.get_path_name(use_galsim=use_galsim) - ) - ) + os.mkdir(self.params["output"]["dir"].replace("truth", self.get_path_name(use_galsim=use_galsim))) new_filename = os.path.join(self.params["output"]["dir"], imfilename).replace( "truth", self.get_path_name(use_galsim=use_galsim) ) @@ -162,15 +142,9 @@ def __init__( b = galsim.BoundsI(xmin=1, ymin=1, xmax=roman.n_pix, ymax=roman.n_pix) old_filename = os.path.join(self.params["output"]["dir"], imfilename) if not os.path.exists( - self.params["output"]["dir"].replace( - "truth", self.get_path_name(use_galsim=use_galsim) - ) + self.params["output"]["dir"].replace("truth", self.get_path_name(use_galsim=use_galsim)) ): - os.mkdir( - self.params["output"]["dir"].replace( - "truth", self.get_path_name(use_galsim=use_galsim) - ) - ) + os.mkdir(self.params["output"]["dir"].replace("truth", self.get_path_name(use_galsim=use_galsim))) new_filename = os.path.join(self.params["output"]["dir"], imfilename).replace( "truth", self.get_path_name(use_galsim=use_galsim) ) @@ -186,12 +160,8 @@ def __init__( force_cvz = True self.setup_sky(im, self.pointing, rng, visit * sca, force_cvz=force_cvz) - img, err, dq, sky_mean, sky_noise = self.add_effects( - im, None, self.pointing, use_galsim=use_galsim - ) - img, err, dq, sky_mean, sky_noise = self.add_effects( - im, None, self.pointing, use_galsim=use_galsim - ) + img, err, dq, sky_mean, sky_noise = self.add_effects(im, None, self.pointing, use_galsim=use_galsim) + img, err, dq, sky_mean, sky_noise = self.add_effects(im, None, self.pointing, use_galsim=use_galsim) write_fits( old_filename, @@ -548,9 +518,7 @@ def bfe(self, im): # solve boundary shfit kernel aX components # ======================================================================= a_area = self.df["BFE"][:, :, :, :] # 5x5x32x32 - a_components = np.zeros( - (4, 2 * nbfe + 1, 2 * nbfe + 1, n_max, m_max) - ) # 4x5x5x32x32 + a_components = np.zeros((4, 2 * nbfe + 1, 2 * nbfe + 1, n_max, m_max)) # 4x5x5x32x32 # solve aR aT aL aB for each a for n in range(n_max): # m_max and n_max = 32 (binned in 128x128) @@ -623,15 +591,11 @@ def bfe(self, im): # The img is clipped by the saturation level here to cap the brighter fatter effect # and avoid unphysical behavior - array_pad = self.saturate(im.copy()).array[ - 4:-4, 4:-4 - ] # img of interest 4088x4088 + array_pad = self.saturate(im.copy()).array[4:-4, 4:-4] # img of interest 4088x4088 array_pad = np.pad( array_pad, [(4 + nbfe, 4 + nbfe), (4 + nbfe, 4 + nbfe)], mode="symmetric" ) # 4100x4100 array - array_pad = self.saturate(im.copy()).array[ - 4:-4, 4:-4 - ] # img of interest 4088x4088 + array_pad = self.saturate(im.copy()).array[4:-4, 4:-4] # img of interest 4088x4088 array_pad = np.pad( array_pad, [(4 + nbfe, 4 + nbfe), (4 + nbfe, 4 + nbfe)], mode="symmetric" ) # 4100x4100 array @@ -681,7 +645,9 @@ def bfe(self, im): gi * m_sub : (gi + 1) * m_sub, ] ) - ).dot(t.T) # sub_grid*sub_grid + ).dot( + t.T + ) # sub_grid*sub_grid a_components_pad[comp, j, i, :, :] = np.pad( tmp, [(nbfe, nbfe), (nbfe, nbfe)], mode="symmetric" ) @@ -703,10 +669,14 @@ def bfe(self, im): nbfe - dx : nbfe - dx + bin_size * m_sub, ] * array_pad[ - -dy + nbfe + gj * bin_size * n_sub : -dy + -dy + + nbfe + + gj * bin_size * n_sub : -dy + nbfe + (gj + 1) * bin_size * n_sub, - -dx + nbfe + gi * bin_size * m_sub : -dx + -dx + + nbfe + + gi * bin_size * m_sub : -dx + nbfe + (gi + 1) * bin_size * m_sub, ] @@ -721,18 +691,12 @@ def bfe(self, im): gi * bin_size * m_sub : (gi + 1) * bin_size * m_sub, ] *= 0.5 * ( array_pad[ - nbfe + gj * bin_size * n_sub : nbfe - + (gj + 1) * bin_size * n_sub, - nbfe + gi * bin_size * m_sub : nbfe - + (gi + 1) * bin_size * m_sub, + nbfe + gj * bin_size * n_sub : nbfe + (gj + 1) * bin_size * n_sub, + nbfe + gi * bin_size * m_sub : nbfe + (gi + 1) * bin_size * m_sub, ] + array_pad[ - dj + nbfe + gj * bin_size * n_sub : dj - + nbfe - + (gj + 1) * bin_size * n_sub, - di + nbfe + gi * bin_size * m_sub : di - + nbfe - + (gi + 1) * bin_size * m_sub, + dj + nbfe + gj * bin_size * n_sub : dj + nbfe + (gj + 1) * bin_size * n_sub, + di + nbfe + gi * bin_size * m_sub : di + nbfe + (gi + 1) * bin_size * m_sub, ] ) dQ_components[ @@ -741,18 +705,12 @@ def bfe(self, im): gi * bin_size * m_sub : (gi + 1) * bin_size * m_sub, ] *= 0.5 * ( array_pad[ - nbfe + gj * bin_size * n_sub : nbfe - + (gj + 1) * bin_size * n_sub, - nbfe + gi * bin_size * m_sub : nbfe - + (gi + 1) * bin_size * m_sub, + nbfe + gj * bin_size * n_sub : nbfe + (gj + 1) * bin_size * n_sub, + nbfe + gi * bin_size * m_sub : nbfe + (gi + 1) * bin_size * m_sub, ] + array_pad[ - dj + nbfe + gj * bin_size * n_sub : dj - + nbfe - + (gj + 1) * bin_size * n_sub, - di + nbfe + gi * bin_size * m_sub : di - + nbfe - + (gi + 1) * bin_size * m_sub, + dj + nbfe + gj * bin_size * n_sub : dj + nbfe + (gj + 1) * bin_size * n_sub, + di + nbfe + gi * bin_size * m_sub : di + nbfe + (gi + 1) * bin_size * m_sub, ] ) @@ -773,20 +731,14 @@ def get_eff_sky_bg(self, pointing, radec): radec : World coordinate position of image """ - sky_level = roman.getSkyLevel( - pointing.bpass, world_pos=radec, date=pointing.date - ) + sky_level = roman.getSkyLevel(pointing.bpass, world_pos=radec, date=pointing.date) sky_level *= (1.0 + roman.stray_light_fraction) * roman.pixel_scale**2 - sky_level = roman.getSkyLevel( - pointing.bpass, world_pos=radec, date=pointing.date - ) + sky_level = roman.getSkyLevel(pointing.bpass, world_pos=radec, date=pointing.date) sky_level *= (1.0 + roman.stray_light_fraction) * roman.pixel_scale**2 return sky_level - def translate_cvz( - self, orig_radec, field_ra=9.5, field_dec=-44, cvz_ra=61.24, cvz_dec=-48.42 - ): + def translate_cvz(self, orig_radec, field_ra=9.5, field_dec=-44, cvz_ra=61.24, cvz_dec=-48.42): ra = orig_radec.ra / galsim.degrees - field_ra dec = orig_radec.dec / galsim.degrees - field_dec @@ -819,12 +771,10 @@ def setup_sky(self, im, pointing, rng, rng_iter, force_cvz=False): self.dark_current_ = roman.dark_current * roman.exptime else: self.dark_current_ = ( - roman.dark_current * roman.exptime - + self.df["DARK"][:, :].flatten() * roman.exptime + roman.dark_current * roman.exptime + self.df["DARK"][:, :].flatten() * roman.exptime ) self.dark_current_ = ( - roman.dark_current * roman.exptime - + self.df["DARK"][:, :].flatten() * roman.exptime + roman.dark_current * roman.exptime + self.df["DARK"][:, :].flatten() * roman.exptime ) if self.df is None: self.gain = roman.gain @@ -837,13 +787,9 @@ def setup_sky(self, im, pointing, rng, rng_iter, force_cvz=False): radec = self.translate_cvz(pointing.radec) else: radec = pointing.radec - sky_level = roman.getSkyLevel( - pointing.bpass, world_pos=radec, date=pointing.date - ) + sky_level = roman.getSkyLevel(pointing.bpass, world_pos=radec, date=pointing.date) sky_level *= 1.0 + roman.stray_light_fraction - sky_level = roman.getSkyLevel( - pointing.bpass, world_pos=radec, date=pointing.date - ) + sky_level = roman.getSkyLevel(pointing.bpass, world_pos=radec, date=pointing.date) sky_level *= 1.0 + roman.stray_light_fraction # Make a image of the sky that takes into account the spatially variable pixel scale. Note # that makeSkyImage() takes a bit of time. If you do not care about the variable pixel @@ -863,10 +809,7 @@ def setup_sky(self, im, pointing, rng, rng_iter, force_cvz=False): # ~ (35, 3050, 0.008, 1.2E6) (e-/p). Hot pixels could be removed in further analysis using the # dq array. self.sky_mean = np.mean( - np.round( - (np.round(self.sky.array) + round(np.median(self.dark_current_))) - / np.mean(self.gain) - ) + np.round((np.round(self.sky.array) + round(np.median(self.dark_current_))) / np.mean(self.gain)) ) self.sky.addNoise(self.noise) @@ -902,9 +845,7 @@ def add_background(self, im): return im - def recip_failure( - self, im, exptime=roman.exptime, alpha=roman.reciprocity_alpha, base_flux=1.0 - ): + def recip_failure(self, im, exptime=roman.exptime, alpha=roman.reciprocity_alpha, base_flux=1.0): """ Introduce reciprocity failure to image. @@ -946,9 +887,7 @@ def dark_current(self, im): if self.df is None: self.im_dark = im.copy() dark_current_ = self.dark_current_ - dark_noise = galsim.DeviateNoise( - galsim.PoissonDeviate(self.rng, dark_current_) - ) + dark_noise = galsim.DeviateNoise(galsim.PoissonDeviate(self.rng, dark_current_)) im.addNoise(dark_noise) self.im_dark = im - self.im_dark @@ -956,17 +895,9 @@ def dark_current(self, im): dark_current_ = self.dark_current_.clip(0) # opt for numpy random geneator instead for speed - self.im_dark = ( - self.rng_np.poisson(dark_current_) - .reshape(im.array.shape) - .astype(im.dtype) - ) + self.im_dark = self.rng_np.poisson(dark_current_).reshape(im.array.shape).astype(im.dtype) im.array[:, :] += self.im_dark - self.im_dark = ( - self.rng_np.poisson(dark_current_) - .reshape(im.array.shape) - .astype(im.dtype) - ) + self.im_dark = self.rng_np.poisson(dark_current_).reshape(im.array.shape).astype(im.dtype) im.array[:, :] += self.im_dark # NOTE: Sky level and dark current might appear like a constant background that can be @@ -1057,26 +988,14 @@ def add_persistence(self, im, pointing): dither_sca_array = np.loadtxt(self.params["dither_from_file"]).astype(int) # select adjacent exposures for the same sca (within 10*roman.exptime) - dither_list_selected = dither_sca_array[ - dither_sca_array[:, 1] == pointing.sca, 0 - ] - dither_list_selected = dither_list_selected[ - np.abs(dither_list_selected - pointing.visit) < 10 - ] - p_list = np.array( - [get_pointing(self.params, i, pointing.sca) for i in dither_list_selected] - ) + dither_list_selected = dither_sca_array[dither_sca_array[:, 1] == pointing.sca, 0] + dither_list_selected = dither_list_selected[np.abs(dither_list_selected - pointing.visit) < 10] + p_list = np.array([get_pointing(self.params, i, pointing.sca) for i in dither_list_selected]) dt_list = np.array([(pointing.date - p.date).total_seconds() for p in p_list]) p_pers = p_list[np.where((dt_list > 0) & (dt_list < roman.exptime * 10))] - dither_list_selected = dither_sca_array[ - dither_sca_array[:, 1] == pointing.sca, 0 - ] - dither_list_selected = dither_list_selected[ - np.abs(dither_list_selected - pointing.visit) < 10 - ] - p_list = np.array( - [get_pointing(self.params, i, pointing.sca) for i in dither_list_selected] - ) + dither_list_selected = dither_sca_array[dither_sca_array[:, 1] == pointing.sca, 0] + dither_list_selected = dither_list_selected[np.abs(dither_list_selected - pointing.visit) < 10] + p_list = np.array([get_pointing(self.params, i, pointing.sca) for i in dither_list_selected]) dt_list = np.array([(pointing.date - p.date).total_seconds() for p in p_list]) p_pers = p_list[np.where((dt_list > 0) & (dt_list < roman.exptime * 10))] @@ -1087,9 +1006,7 @@ def add_persistence(self, im, pointing): pointing.date - p.date ).total_seconds() - roman.exptime / 2 # avg time since end of exposures self.params["output"]["file_name"]["items"] = [p.filter, p.visit, p.sca] - imfilename = ParseValue( - self.params["output"], "file_name", self.params, str - )[0] + imfilename = ParseValue(self.params["output"], "file_name", self.params, str)[0] fn = os.path.join(self.params["output"]["dir"], imfilename) # apply all the effects that occured before persistence on the previouse exposures @@ -1101,14 +1018,8 @@ def add_persistence(self, im, pointing): x = x.clip(0) # remove negative stimulus - im.array[:, :] += ( - galsim.roman.roman_detectors.fermi_linear(x.array, dt) - * roman.exptime - ) - im.array[:, :] += ( - galsim.roman.roman_detectors.fermi_linear(x.array, dt) - * roman.exptime - ) + im.array[:, :] += galsim.roman.roman_detectors.fermi_linear(x.array, dt) * roman.exptime + im.array[:, :] += galsim.roman.roman_detectors.fermi_linear(x.array, dt) * roman.exptime else: # setup parameters for persistence @@ -1126,12 +1037,10 @@ def add_persistence(self, im, pointing): pointing.date - p.date ).total_seconds() - roman.exptime / 2 # avg time since end of exposures fac_dt = ( - (roman.exptime / 2.0) / dt - ) # linear time dependence (approximate until we get t1 and Delat t of the data) + roman.exptime / 2.0 + ) / dt # linear time dependence (approximate until we get t1 and Delat t of the data) self.params["output"]["file_name"]["items"] = [p.filter, p.visit, p.sca] - imfilename = ParseValue( - self.params["output"], "file_name", self.params, str - )[0] + imfilename = ParseValue(self.params["output"], "file_name", self.params, str)[0] fn = os.path.join(self.params["output"]["dir"], imfilename) # apply all the effects that occured before persistence on the previouse exposures @@ -1172,9 +1081,7 @@ def add_persistence(self, im, pointing): a = np.zeros(x.shape) a += ((x >= Q05) & (x < Q06)) * (x - Q05) / (Q06 - Q05) - a += (x >= Q06) * ( - x / Q06 - ) ** alpha # avoid fractional power of negative values + a += (x >= Q06) * (x / Q06) ** alpha # avoid fractional power of negative values im.array[:, :] += a * self.df["PERSIST"][5, :, :][0] * fac_dt return im @@ -1232,13 +1139,9 @@ def interpix_cap(self, im, kernel=roman.ipc_kernel): grid_size = 4096 // num_grids array_pad = im.array[4:-4, 4:-4] # it's an array instead of img - array_pad = np.pad( - array_pad, [(5, 5), (5, 5)], mode="symmetric" - ) # 4098x4098 array + array_pad = np.pad(array_pad, [(5, 5), (5, 5)], mode="symmetric") # 4098x4098 array array_pad = im.array[4:-4, 4:-4] # it's an array instead of img - array_pad = np.pad( - array_pad, [(5, 5), (5, 5)], mode="symmetric" - ) # 4098x4098 array + array_pad = np.pad(array_pad, [(5, 5), (5, 5)], mode="symmetric") # 4098x4098 array K = self.df["IPC"][:, :, :, :] # 3,3,512, 512 @@ -1255,13 +1158,9 @@ def interpix_cap(self, im, kernel=roman.ipc_kernel): for j in range(3): for i in range(3): tmp = (t.dot(K[j, i, :, :])).dot(t.T) # grid_sizexgrid_size - K_pad[j, i, :, :] = np.pad( - tmp, [(1, 1), (1, 1)], mode="symmetric" - ) + K_pad[j, i, :, :] = np.pad(tmp, [(1, 1), (1, 1)], mode="symmetric") tmp = (t.dot(K[j, i, :, :])).dot(t.T) # grid_sizexgrid_size - K_pad[j, i, :, :] = np.pad( - tmp, [(1, 1), (1, 1)], mode="symmetric" - ) + K_pad[j, i, :, :] = np.pad(tmp, [(1, 1), (1, 1)], mode="symmetric") for dy in range(-1, 2): for dx in range(-1, 2): @@ -1276,12 +1175,8 @@ def interpix_cap(self, im, kernel=roman.ipc_kernel): 1 - dx : 1 - dx + grid_size, ] * array_pad[ - 1 - dy + gj * grid_size : 1 - - dy - + (gj + 1) * grid_size, - 1 - dx + gi * grid_size : 1 - - dx - + (gi + 1) * grid_size, + 1 - dy + gj * grid_size : 1 - dy + (gj + 1) * grid_size, + 1 - dx + gi * grid_size : 1 - dx + (gi + 1) * grid_size, ] ) @@ -1310,16 +1205,12 @@ def add_read_noise(self, im): # use numpy random generator to draw 2-d noise map read_noise = self.df["READ"][2, :, :].flatten() # flattened 4096x4096 array self.im_read = ( - self.rng_np.normal(loc=0.0, scale=read_noise) - .reshape(im.array.shape) - .astype(im.dtype) + self.rng_np.normal(loc=0.0, scale=read_noise).reshape(im.array.shape).astype(im.dtype) ) im.array[:, :] += self.im_read read_noise = self.df["READ"][2, :, :].flatten() # flattened 4096x4096 array self.im_read = ( - self.rng_np.normal(loc=0.0, scale=read_noise) - .reshape(im.array.shape) - .astype(im.dtype) + self.rng_np.normal(loc=0.0, scale=read_noise).reshape(im.array.shape).astype(im.dtype) ) im.array[:, :] += self.im_read diff --git a/roman_imsim/noise.py b/roman_imsim/noise.py index 5ff071bc..b2fcd4ae 100644 --- a/roman_imsim/noise.py +++ b/roman_imsim/noise.py @@ -46,9 +46,7 @@ def addNoise(self, config, base, image, rng, current_var, draw_method, logger): "gain": dict, } - params, safe = galsim.config.GetAllParams( - config, base, req={}, opt=opt, ignore=[] - ) + params, safe = galsim.config.GetAllParams(config, base, req={}, opt=opt, ignore=[]) mjd = params.get("mjd", None) stray_light = params.get("stray_light", False) @@ -68,22 +66,16 @@ def addNoise(self, config, base, image, rng, current_var, draw_method, logger): filter_name = bp.name exptime, _ = galsim.config.ParseValue(base["image"], "exptime", base, float) date = Time(mjd, format="mjd").to_datetime() if mjd is not None else None - logger.info( - "image %d: Start RomanSCA detector effects", base.get("image_num", 0) - ) + logger.info("image %d: Start RomanSCA detector effects", base.get("image_num", 0)) # Things that will eventually be subtracted (if sky_subtract) will have their expectation # value added to sky_image. So technically, this includes things that aren't just sky. # E.g. includes dark_current and thermal backgrounds. sky_image = image.copy() - sky_level = models.backgrounds.getSkyLevel( - bp, world_pos=wcs.toWorld(image.true_center), date=date - ) + sky_level = models.backgrounds.getSkyLevel(bp, world_pos=wcs.toWorld(image.true_center), date=date) logger.debug("Adding sky_level = %s", sky_level) if stray_light: - logger.debug( - "Stray light fraction = %s", models.parameters.stray_light_fraction - ) + logger.debug("Stray light fraction = %s", models.parameters.stray_light_fraction) sky_level *= 1.0 + models.parameters.stray_light_fraction wcs.makeSkyImage(sky_image, sky_level) diff --git a/roman_imsim/obseq.py b/roman_imsim/obseq.py index 9bf4b3c2..79c05263 100644 --- a/roman_imsim/obseq.py +++ b/roman_imsim/obseq.py @@ -26,17 +26,11 @@ def __init__(self, file_name, visit, SCA, logger=None): def read_obseq(self): """Read visit info from the obseq file.""" if self.file_name is None: - raise ValueError( - "No obseq filename provided, trying to build from config information." - ) + raise ValueError("No obseq filename provided, trying to build from config information.") if self.visit is None: - raise ValueError( - "The visit must be set when reading visit info from an obseq file." - ) + raise ValueError("The visit must be set when reading visit info from an obseq file.") - self.logger.warning( - "Reading info from obseq file %s for visit %s", self.file_name, self.visit - ) + self.logger.warning("Reading info from obseq file %s for visit %s", self.file_name, self.visit) ob = fio.FITS(self.file_name)[-1][self.visit] @@ -68,9 +62,5 @@ def ObSeqData(config, base, value_type): return val, safe -RegisterInputType( - "obseq_data", InputLoader(ObSeqDataLoader, file_scope=True, takes_logger=True) -) -RegisterValueType( - "ObSeqData", ObSeqData, [float, int, str, Angle], input_type="obseq_data" -) +RegisterInputType("obseq_data", InputLoader(ObSeqDataLoader, file_scope=True, takes_logger=True)) +RegisterValueType("ObSeqData", ObSeqData, [float, int, str, Angle], input_type="obseq_data") diff --git a/roman_imsim/psf.py b/roman_imsim/psf.py index f4ba979b..c0e7d5c4 100644 --- a/roman_imsim/psf.py +++ b/roman_imsim/psf.py @@ -125,9 +125,7 @@ def _parse_pupil_bin(self, pupil_bin): else: return pupil_bin - def _psf_call( - self, SCA, bpass, SCA_pos, WCS, pupil_bin, n_waves, logger, extra_aberrations - ): + def _psf_call(self, SCA, bpass, SCA_pos, WCS, pupil_bin, n_waves, logger, extra_aberrations): if pupil_bin == 8: psf = models.psf_utils.getPSF( @@ -246,10 +244,10 @@ def initPSF( self._image_xsize = image_xsize if self._image_xsize is None: - self._image_xsize = roman.n_pix + self._image_xsize = models.parameters.n_pix self._image_ysize = image_ysize if self._image_ysize is None: - self._image_ysize = roman.n_pix + self._image_ysize = models.parameters.n_pix corners = [ galsim.PositionD(1, 1), @@ -308,9 +306,9 @@ def getPSF(self, pupil_bin, pos): wlu = (models.parameters.n_pix - pos.x) * (pos.y - 1) wul = (pos.x - 1) * (models.parameters.n_pix - pos.y) wuu = (pos.x - 1) * (pos.y - 1) - return ( - wll * psf["ll"] + wlu * psf["lu"] + wul * psf["ul"] + wuu * psf["uu"] - ) / ((models.parameters.n_pix - 1) * (models.parameters.n_pix - 1)) + return (wll * psf["ll"] + wlu * psf["lu"] + wul * psf["ul"] + wuu * psf["uu"]) / ( + (models.parameters.n_pix - 1) * (models.parameters.n_pix - 1) + ) def RegisterPSFInterpolatorType(interp_type, builder, input_type=None): @@ -493,7 +491,7 @@ def BuildRomanPSF(config, base, ignore, gsparams, logger): base["image_pos"], ) else: - psf = roman.getPSF( + psf = models.psf_utils.getPSF( SCA=base["SCA"], bandpass=base["bandpass"].name, SCA_pos=base["image_pos"], diff --git a/roman_imsim/sca.py b/roman_imsim/sca.py index 46559cc2..7524a21a 100644 --- a/roman_imsim/sca.py +++ b/roman_imsim/sca.py @@ -52,9 +52,7 @@ def setup(self, config, base, image_num, obj_num, ignore, logger): "draw_method": str, "use_fft_bright": bool, } - params = galsim.config.GetAllParams( - config, base, req=req, opt=opt, ignore=ignore + extra_ignore - )[0] + params = galsim.config.GetAllParams(config, base, req=req, opt=opt, ignore=ignore + extra_ignore)[0] self.sca = params["SCA"] base["SCA"] = self.sca @@ -63,15 +61,11 @@ def setup(self, config, base, image_num, obj_num, ignore, logger): self.exptime = params["exptime"] # If draw_method isn't in image field, it may be in stamp. Check. - self.draw_method = params.get( - "draw_method", base.get("stamp", {}).get("draw_method", "phot") - ) + self.draw_method = params.get("draw_method", base.get("stamp", {}).get("draw_method", "phot")) # If user hasn't overridden the bandpass to use, get the standard one. if "bandpass" not in config: - base["bandpass"] = galsim.config.BuildBandpass( - base["image"], "bandpass", base, logger=logger - ) + base["bandpass"] = galsim.config.BuildBandpass(base["image"], "bandpass", base, logger=logger) return models.parameters.n_pix, models.parameters.n_pix @@ -104,13 +98,9 @@ def buildImage(self, config, base, image_num, obj_num, logger): full_image.header["VERSION"] = "unknown" full_image.header["EXPTIME"] = self.exptime full_image.header["MJD-OBS"] = self.mjd - full_image.header["DATE-OBS"] = Time( - self.mjd, format="mjd" - ).datetime.isoformat() + full_image.header["DATE-OBS"] = Time(self.mjd, format="mjd").datetime.isoformat() full_image.header["FILTER"] = self.filter - full_image.header["ZPTMAG"] = 2.5 * np.log10( - self.exptime * models.parameters.collecting_area - ) + full_image.header["ZPTMAG"] = 2.5 * np.log10(self.exptime * models.parameters.collecting_area) base["current_image"] = full_image @@ -161,9 +151,7 @@ def buildImage(self, config, base, image_num, obj_num, logger): # avoid an error. But this isn't covered in the imsim test suite. continue - logger.debug( - "image %d: full bounds = %s", image_num, str(full_image.bounds) - ) + logger.debug("image %d: full bounds = %s", image_num, str(full_image.bounds)) logger.debug( "image %d: stamp %d bounds = %s", image_num, diff --git a/roman_imsim/scafile.py b/roman_imsim/scafile.py index a9315037..75574a64 100644 --- a/roman_imsim/scafile.py +++ b/roman_imsim/scafile.py @@ -27,9 +27,7 @@ def setup(self, config, base, file_num, logger): logger.debug("file %d: seed = %d", file_num, seed) if "exptime" in config: - base["exptime"] = galsim.config.ParseValue(config, "exptime", base, float)[ - 0 - ] + base["exptime"] = galsim.config.ParseValue(config, "exptime", base, float)[0] else: base["exptime"] = models.parameters.exptime @@ -94,9 +92,7 @@ def buildImages(self, config, base, file_num, image_num, obj_num, ignore, logger if cosmic_ray_rate > 0: cosmic_ray_catalog = params.get("cosmic_ray_catalog", None) if cosmic_ray_catalog is None: - cosmic_ray_catalog = os.path.join( - data_dir, "cosmic_rays_itl_2017.fits.gz" - ) + cosmic_ray_catalog = os.path.join(data_dir, "cosmic_rays_itl_2017.fits.gz") if not os.path.isfile(cosmic_ray_catalog): raise FileNotFoundError(f"{cosmic_ray_catalog} not found") diff --git a/roman_imsim/skycat.py b/roman_imsim/skycat.py index b2de68e0..19274a4e 100644 --- a/roman_imsim/skycat.py +++ b/roman_imsim/skycat.py @@ -79,9 +79,7 @@ def __init__( if obj_types is not None: self.logger.warning(f"Object types restricted to {obj_types}") - self.sca_center = wcs.toWorld( - galsim.PositionD(self.xsize / 2.0, self.ysize / 2.0) - ) + self.sca_center = wcs.toWorld(galsim.PositionD(self.xsize / 2.0, self.ysize / 2.0)) self._objects = None # import os, psutil @@ -113,14 +111,10 @@ def objects(self): vertices = [] for x, y in corners: sky_coord = self.wcs.toWorld(galsim.PositionD(x, y)) - vertices.append( - (sky_coord.ra / galsim.degrees, sky_coord.dec / galsim.degrees) - ) + vertices.append((sky_coord.ra / galsim.degrees, sky_coord.dec / galsim.degrees)) region = PolygonalRegion(vertices) sky_cat = skyCatalogs.open_catalog(self.file_name) - self._objects = sky_cat.get_objects_by_region( - region, obj_type_set=self.obj_types, mjd=self.mjd - ) + self._objects = sky_cat.get_objects_by_region(region, obj_type_set=self.obj_types, mjd=self.mjd) if not self._objects: self.logger.warning("No objects found on image.") # import os, psutil @@ -215,9 +209,7 @@ def getObj(self, index, gsparams=None, rng=None, exptime=30): gs_obj_list = [] for component in gsobjs: if faint: - gsobjs[component] = gsobjs[component].evaluateAtWavelength( - self.bandpass - ) + gsobjs[component] = gsobjs[component].evaluateAtWavelength(self.bandpass) gs_obj_list.append(gsobjs[component] * self._trivial_sed) else: if component in seds: @@ -244,9 +236,7 @@ def getObj(self, index, gsparams=None, rng=None, exptime=30): gs_object.flux = flux # Get the object type - if (skycat_obj.object_type == "diffsky_galaxy") | ( - skycat_obj.object_type == "galaxy" - ): + if (skycat_obj.object_type == "diffsky_galaxy") | (skycat_obj.object_type == "galaxy"): gs_object.object_type = "galaxy" if skycat_obj.object_type in {"star", "gaia_star"}: gs_object.object_type = "star" @@ -276,9 +266,7 @@ def getKwargs(self, config, base, logger): kwargs["logger"] = logger if "bandpass" not in config: - base["bandpass"] = galsim.config.BuildBandpass( - base["image"], "bandpass", base, logger=logger - )[0] + base["bandpass"] = galsim.config.BuildBandpass(base["image"], "bandpass", base, logger=logger)[0] kwargs["bandpass"] = base["bandpass"] # Sky catalog object lists are created per CCD, so they are @@ -352,9 +340,7 @@ def SkyCatWorldPos(config, base, value_type): RegisterInputType("sky_catalog", SkyCatalogLoader(SkyCatalogInterface, has_nobj=True)) RegisterObjectType("SkyCatObj", SkyCatObj, input_type="sky_catalog") -RegisterValueType( - "SkyCatWorldPos", SkyCatWorldPos, [galsim.CelestialCoord], input_type="sky_catalog" -) +RegisterValueType("SkyCatWorldPos", SkyCatWorldPos, [galsim.CelestialCoord], input_type="sky_catalog") # This class was modified from https://github.com/LSSTDESC/imSim/. License info follows: diff --git a/roman_imsim/stamp.py b/roman_imsim/stamp.py index 9e57c67b..47fa8106 100644 --- a/roman_imsim/stamp.py +++ b/roman_imsim/stamp.py @@ -83,9 +83,7 @@ def setup(self, config, base, xsize, ysize, ignore, logger): else: gal_achrom = gal.evaluateAtWavelength(bandpass.effective_wavelength) - if hasattr(gal_achrom, "original") and isinstance( - gal_achrom.original, galsim.DeltaFunction - ): + if hasattr(gal_achrom, "original") and isinstance(gal_achrom.original, galsim.DeltaFunction): # For bright stars, set the following stamp size limits if self.flux < 1e6: image_size = 500 @@ -107,16 +105,12 @@ def setup(self, config, base, xsize, ysize, ignore, logger): # print('stamp setup3',process.memory_info().rss) base["pupil_bin"] = self.pupil_bin logger.info("Object flux is %d", self.flux) - logger.info( - "Object %d will use stamp size = %s", base.get("obj_num", 0), image_size - ) + logger.info("Object %d will use stamp size = %s", base.get("obj_num", 0), image_size) # Determine where this object is going to go: # This is the same as what the base StampBuilder does: if "image_pos" in config: - image_pos = galsim.config.ParseValue( - config, "image_pos", base, galsim.PositionD - )[0] + image_pos = galsim.config.ParseValue(config, "image_pos", base, galsim.PositionD)[0] else: image_pos = None @@ -144,9 +138,7 @@ def buildPSF(self, config, base, gsparams, logger): the PSF """ if base.get("psf", {}).get("type", "roman_psf") != "roman_psf": - return galsim.config.BuildGSObject( - base, "psf", gsparams=gsparams, logger=logger - )[0] + return galsim.config.BuildGSObject(base, "psf", gsparams=gsparams, logger=logger)[0] roman_psf = galsim.config.GetInputObj("roman_psf", config, base, "buildPSF") psf = roman_psf.getPSF(self.pupil_bin, base["image_pos"]) @@ -164,9 +156,7 @@ def getDrawMethod(self, config, base, logger): method = galsim.config.ParseValue(config, "draw_method", base, str)[0] self.use_fft_bright = False if "use_fft_bright" in config: - self.use_fft_bright = galsim.config.ParseValue( - config, "use_fft_bright", base, bool - )[0] + self.use_fft_bright = galsim.config.ParseValue(config, "use_fft_bright", base, bool)[0] if method not in galsim.config.valid_draw_methods: raise galsim.GalSimConfigValueError( @@ -178,16 +168,12 @@ def getDrawMethod(self, config, base, logger): logger.info("Auto -> Use FFT drawing for object %d.", base["obj_num"]) return "fft" else: - logger.info( - "Auto -> Use photon shooting for object %d.", base["obj_num"] - ) + logger.info("Auto -> Use photon shooting for object %d.", base["obj_num"]) return "phot" else: # If user sets something specific for the method, rather than auto, # then respect their wishes. - logger.info( - "Use specified method=%s for object %d.", method, base["obj_num"] - ) + logger.info("Use specified method=%s for object %d.", method, base["obj_num"]) return method @classmethod @@ -200,16 +186,11 @@ def _fix_seds_24(cls, prof, bandpass): sed = prof.sed # TODO: This bit should probably be ported back to Galsim. # Something like sed.make_tabulated() - if ( - not isinstance(sed._spec, galsim.LookupTable) - or sed._spec.interpolant != "linear" - ): + if not isinstance(sed._spec, galsim.LookupTable) or sed._spec.interpolant != "linear": # Workaround for https://github.com/GalSim-developers/GalSim/issues/1228 f = np.broadcast_to(sed(wave_list), wave_list.shape) new_spec = galsim.LookupTable(wave_list, f, interpolant="linear") - new_sed = galsim.SED( - new_spec, "nm", "fphotons" if sed.spectral else "1" - ) + new_sed = galsim.SED(new_spec, "nm", "fphotons" if sed.spectral else "1") prof.sed = new_sed # Also recurse onto any components. @@ -294,9 +275,7 @@ def draw(self, prof, image, method, offset, config, base, logger): # print('stamp draw3b ',process.memory_info().rss) if not faint and "photon_ops" in config: - photon_ops = galsim.config.BuildPhotonOps( - config, "photon_ops", base, logger - ) + photon_ops = galsim.config.BuildPhotonOps(config, "photon_ops", base, logger) else: photon_ops = [] @@ -334,9 +313,7 @@ def draw(self, prof, image, method, offset, config, base, logger): if not faint and config.get("fft_photon_ops"): kwargs.update( { - "photon_ops": galsim.config.BuildPhotonOps( - config, "fft_photon_ops", base, logger - ), + "photon_ops": galsim.config.BuildPhotonOps(config, "fft_photon_ops", base, logger), "maxN": maxN, "rng": self.rng, "n_subsample": 1, @@ -354,9 +331,7 @@ def draw(self, prof, image, method, offset, config, base, logger): # and raise the error. logger.error("Caught error trying to draw using FFT:") logger.error("%s", e) - logger.error( - "You may need to add a gsparams field with maximum_fft_size to" - ) + logger.error("You may need to add a gsparams field with maximum_fft_size to") logger.error("either the psf or gal field to allow larger FFTs.") logger.info("prof = %r", prof) logger.info("fft_image = %s", fft_image) diff --git a/roman_imsim/utils.py b/roman_imsim/utils.py index 788d4bbe..19e39170 100644 --- a/roman_imsim/utils.py +++ b/roman_imsim/utils.py @@ -11,9 +11,7 @@ class roman_utils(object): Class to contain a variety of helper routines to work with the simulation data. """ - def __init__( - self, config_file, visit=None, sca=None, image_name=None, setup_skycat=False - ): + def __init__(self, config_file, visit=None, sca=None, image_name=None, setup_skycat=False): """ Setup information about a simulated Roman image. Parameters: @@ -38,16 +36,10 @@ def __init__( self.skycat = galsim.config.GetInputObj( "sky_catalog", config["input"]["sky_catalog"], config, "sky_catalog" ) - self.PSF = galsim.config.GetInputObj( - "roman_psf", config["input"]["roman_psf"], config, "roman_psf" - ) + self.PSF = galsim.config.GetInputObj("roman_psf", config["input"]["roman_psf"], config, "roman_psf") self.wcs = galsim.config.BuildWCS(config["image"], "wcs", config) - self.bpass = galsim.config.BuildBandpass( - config["image"], "bandpass", config, None - )[0] - self.photon_ops = galsim.config.BuildPhotonOps( - config["stamp"], "photon_ops", config, None - ) + self.bpass = galsim.config.BuildBandpass(config["image"], "bandpass", config, None)[0] + self.photon_ops = galsim.config.BuildPhotonOps(config["stamp"], "photon_ops", config, None) self.rng = galsim.config.GetRNG(config, config["image"], None, "roman_sca") def check_input(self, visit, sca, image_name): @@ -62,9 +54,7 @@ def check_input(self, visit, sca, image_name): tmp = np.array(image_name[start:end].split("_")).astype(int) return tmp[0], tmp[1] if (visit is None) | (sca is None): - raise ValueError( - "Insufficient information to construct visit info - all inputs are None." - ) + raise ValueError("Insufficient information to construct visit info - all inputs are None.") return visit, sca def getPSF(self, x=None, y=None, pupil_bin=8): @@ -84,16 +74,12 @@ def getPSF(self, x=None, y=None, pupil_bin=8): ) return self.PSF.getPSF( pupil_bin, - galsim.PositionD( - models.parameters.n_pix / 2, models.parameters.n_pix / 2 - ), + galsim.PositionD(models.parameters.n_pix / 2, models.parameters.n_pix / 2), ) if (x is None) | (y is None): return self.PSF.getPSF( 8, - galsim.PositionD( - models.parameters.n_pix / 2, models.parameters.n_pix / 2 - ), + galsim.PositionD(models.parameters.n_pix / 2, models.parameters.n_pix / 2), ) return self.PSF.getPSF(8, galsim.PositionD(x, y)) @@ -162,9 +148,7 @@ def getPSF_Image( dvdx=local_wcs.dvdx / oversampling_factor, dvdy=local_wcs.dvdy / oversampling_factor, ) - stamp = galsim.Image( - stamp_size * oversampling_factor, stamp_size * oversampling_factor, wcs=wcs - ) + stamp = galsim.Image(stamp_size * oversampling_factor, stamp_size * oversampling_factor, wcs=wcs) if not include_photonOps: psf = galsim.Convolve(point, self.getPSF(x, y, pupil_bin)) return psf.drawImage(self.bpass, image=stamp, wcs=wcs, method=method) @@ -179,12 +163,8 @@ def getPSF_Image( dy=j * scale / oversampling_factor, ) for i, j in product( - np.arange( - (1 - oversampling_factor) / 2, (1 + oversampling_factor) / 2 - ), - np.arange( - (1 - oversampling_factor) / 2, (1 + oversampling_factor) / 2 - ), + np.arange((1 - oversampling_factor) / 2, (1 + oversampling_factor) / 2), + np.arange((1 - oversampling_factor) / 2, (1 + oversampling_factor) / 2), ) ] ) diff --git a/tests/test_end_to_end.py b/tests/test_end_to_end.py index e75c0314..1d2fd66f 100644 --- a/tests/test_end_to_end.py +++ b/tests/test_end_to_end.py @@ -22,18 +22,10 @@ def test_compare_truth(): reference_truth = ( - Path(TEST_REPO) - / REFERENCE_DIR - / "RomanWAS_new" - / "truth" - / "Roman_WAS_index_J129_12909_4.txt" + Path(TEST_REPO) / REFERENCE_DIR / "RomanWAS_new" / "truth" / "Roman_WAS_index_J129_12909_4.txt" ) output_truth = ( - Path(TEST_REPO) - / OUTPUT_DIR - / "RomanWAS_new" - / "truth" - / "Roman_WAS_index_J129_12909_4.txt" + Path(TEST_REPO) / OUTPUT_DIR / "RomanWAS_new" / "truth" / "Roman_WAS_index_J129_12909_4.txt" ) output_table = ascii.read(output_truth) diff --git a/tests/test_imports.py b/tests/test_imports.py index 2511b8e9..0d1bc016 100644 --- a/tests/test_imports.py +++ b/tests/test_imports.py @@ -9,14 +9,10 @@ def test_all_modules_import(): package = roman_imsim errors = [] - for loader, module_name, is_pkg in pkgutil.walk_packages( - package.__path__, package.__name__ + "." - ): + for loader, module_name, is_pkg in pkgutil.walk_packages(package.__path__, package.__name__ + "."): try: importlib.import_module(module_name) except Exception as e: errors.append((module_name, repr(e))) - assert not errors, "Failed to import modules:\n" + "\n".join( - f"{m}: {e}" for m, e in errors - ) + assert not errors, "Failed to import modules:\n" + "\n".join(f"{m}: {e}" for m, e in errors) From f8c02af935cbc0d263505e654c9e6735ac10f5b4 Mon Sep 17 00:00:00 2001 From: Yuedong Fang Date: Wed, 13 May 2026 15:32:08 -0400 Subject: [PATCH 09/11] temporary change: using the un-released main branch of romanisim --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index b505fadd..73a5d12b 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,7 +23,7 @@ dependencies = [ "skycatalogs", "packaging", "setuptools", - "romanisim", + "romanisim @ git+https://github.com/spacetelescope/romanisim.git@main", ] [project.urls] From 414da08645c8a37f806eb61e1f979e52455dcf88 Mon Sep 17 00:00:00 2001 From: Yuedong Fang Date: Thu, 14 May 2026 14:43:14 -0400 Subject: [PATCH 10/11] switch to use the newly released version of romanisim --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index 73a5d12b..b505fadd 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -23,7 +23,7 @@ dependencies = [ "skycatalogs", "packaging", "setuptools", - "romanisim @ git+https://github.com/spacetelescope/romanisim.git@main", + "romanisim", ] [project.urls] From a7dc540402bb655e1a370dbc09b6c0d315e35f2d Mon Sep 17 00:00:00 2001 From: Yuedong Fang Date: Thu, 28 May 2026 16:43:56 -0400 Subject: [PATCH 11/11] black reformatted. --- roman_imsim/psf.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/roman_imsim/psf.py b/roman_imsim/psf.py index c0e7d5c4..8f6449ad 100644 --- a/roman_imsim/psf.py +++ b/roman_imsim/psf.py @@ -353,9 +353,7 @@ def getKwargs(self, config, base, logger): else: req["SCA"] = int - kwargs, safe = galsim.config.GetAllParams( - config, base, req=req, opt=opt, ignore=ignore - ) + kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt, ignore=ignore) kwargs["extra_aberrations"] = galsim.config.ParseAberrations( "extra_aberrations", config, base, "RomanPSF"