Switch to use romanisim

pyproject.toml

@@ -23,6 +23,7 @@ dependencies = [
   "skycatalogs",
   "packaging",
   "setuptools",
+  "romanisim",
 ]
 
 [project.urls]

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
 
 
@@ -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 = roman.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
 

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"

roman_imsim/detector_physics.py

@@ -172,7 +172,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 +317,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 +453,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"]:
@@ -509,7 +529,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(
                     [
@@ -584,7 +613,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 +623,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):
@@ -858,7 +892,6 @@ def dark_current(self, im):
             self.im_dark = im - self.im_dark
 
         else:
-
             dark_current_ = self.dark_current_.clip(0)
 
             # opt for numpy random geneator instead for speed
@@ -989,7 +1022,6 @@ def add_persistence(self, im, pointing):
                 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"]
@@ -1023,7 +1055,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 +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
@@ -1132,7 +1164,6 @@ def interpix_cap(self, im, kernel=roman.ipc_kernel):
 
                     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,
@@ -1258,7 +1289,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[:, :]