Merge branch 'main' into photon_sampling

Author: Higgins00

.github/workflows/rebase_checker.yaml

@@ -0,0 +1,7 @@
+name: Check that 'main' is not merged into the development branch
+
+on: pull_request
+
+jobs:
+  call-workflow:
+    uses: lsst/rubin_workflows/.github/workflows/rebase_checker.yaml@main

README.md

@@ -1,6 +1,7 @@
-![Build](https://github.com/matroxel/roman_imsim/actions/workflows/build.yml/badge.svg)
+![Build](https://github.com/DukeCosmology/roman_imsim/actions/workflows/build.yml/badge.svg)
+
+# roman_imsim
 
-# wfirst_imsim
 
 Nonstandard dependencies are represented in setup.py.
 

config/default.yaml

@@ -70,13 +70,14 @@ image:
     filter: { type: ObSeqData, field: filter }
     exptime: { type: ObSeqData, field: exptime }
 
+    draw_method: 'phot'
     # 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 slighly slower than achromatic, so the difference
     # is even more pronounced in this case.
-    draw_method: 'auto'
+    use_fft_bright: True
 
     # These are all by default turned on, but you can turn any of them off if desired:
     ignore_noise: True

config/tds.yaml

@@ -71,13 +71,14 @@ image:
     filter: { type: ObSeqData, field: filter }
     exptime: { type: ObSeqData, field: exptime }
 
+    draw_method: 'phot'
     # 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 slighly slower than achromatic, so the difference
     # is even more pronounced in this case.
-    draw_method: 'auto'
+    use_fft_bright: True
 
     # These are all by default turned on, but you can turn any of them off if desired:
     ignore_noise: True

config/was.yaml

@@ -69,13 +69,14 @@ image:
     filter: { type: ObSeqData, field: filter }
     exptime: { type: ObSeqData, field: exptime }
 
+    draw_method: 'phot'
     # 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 slighly slower than achromatic, so the difference
     # is even more pronounced in this case.
-    draw_method: 'auto'
+    use_fft_bright: True
 
     # These are all by default turned on, but you can turn any of them off if desired:
     # ignore_noise: True

pyproject.toml

@@ -16,14 +16,14 @@ dependencies = [
   "galsim",
   "fitsio",
   "astropy",
-  "healpy",
   "numpy",
-  "matplotlib",
-  "scipy"
+  "skycatalogs",
+  "packaging",
+  "setuptools",
 ]
 
 [project.urls]
-Homepage = "https://github.com/matroxel/roman_imsim"
+Homepage = "https://github.com/DukeCosmology/roman_imsim"
 
 [tool.setuptools.packages.find]
 where = ["."]

roman_imsim/skycat.py

@@ -208,14 +208,10 @@ def getObj(self, index, gsparams=None, rng=None, exptime=30):
         for component in gsobjs:
             if faint:
                 gsobjs[component] = gsobjs[component].evaluateAtWavelength(self.bandpass)
-                gs_obj_list.append(
-                    gsobjs[component] * self._trivial_sed * self.exptime * roman.collecting_area
-                )
+                gs_obj_list.append(gsobjs[component] * self._trivial_sed)
             else:
                 if component in seds:
-                    gs_obj_list.append(
-                        gsobjs[component] * seds[component] * self.exptime * roman.collecting_area
-                    )
+                    gs_obj_list.append(gsobjs[component] * seds[component])
 
         if not gs_obj_list:
             return None
@@ -225,14 +221,22 @@ def getObj(self, index, gsparams=None, rng=None, exptime=30):
         else:
             gs_object = galsim.Add(gs_obj_list)
 
+        # This should catch both "star" and "gaia_star" objects
+        if "star" in skycat_obj.object_type:
+            # Cap (star) flux at 30M photons to avoid gross artifacts when trying
+            # to draw the Roman PSF in finite time and memory
+            flux_cap = 3e7
+            if flux > flux_cap:
+                flux = flux_cap
+
         # Give the object the right flux
+        gs_object = gs_object.withFlux(flux, self.bandpass)
         gs_object.flux = flux
-        gs_object.withFlux(gs_object.flux, self.bandpass)
 
         # Get the object type
         if (skycat_obj.object_type == "diffsky_galaxy") | (skycat_obj.object_type == "galaxy"):
             gs_object.object_type = "galaxy"
-        if skycat_obj.object_type == "star":
+        if skycat_obj.object_type in {"star", "gaia_star"}:
             gs_object.object_type = "star"
         if skycat_obj.object_type == "snana":
             gs_object.object_type = "transient"
@@ -251,6 +255,8 @@ def getKwargs(self, config, base, logger):
             "edge_pix": float,
             "obj_types": list,
             "mjd": float,
+            "xsize": int,
+            "ysize": int,
         }
         kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt)
         wcs = galsim.config.BuildWCS(base["image"], "wcs", base, logger=logger)

roman_imsim/stamp.py

@@ -58,18 +58,6 @@ def setup(self, config, base, xsize, ysize, ignore, logger):
             # or cached by the skyCatalogs code.
             gal.flux = gal.calculateFlux(bandpass)
         self.flux = gal.flux
-        # Cap (star) flux at 30M photons to avoid gross artifacts when trying
-        # to draw the Roman PSF in finite time and memory
-        flux_cap = 3e7
-        if self.flux > flux_cap:
-            if (
-                hasattr(gal, "original")
-                and hasattr(gal.original, "original")
-                and isinstance(gal.original.original, galsim.DeltaFunction)
-            ) or (isinstance(gal, galsim.DeltaFunction)):
-                gal = gal.withFlux(flux_cap, bandpass)
-                self.flux = flux_cap
-                gal.flux = flux_cap
         base["flux"] = gal.flux
         base["mag"] = -2.5 * np.log10(gal.flux) + bandpass.zeropoint
         # print('stamp setup2',process.memory_info().rss)
@@ -282,9 +270,12 @@ def draw(self, prof, image, method, offset, config, base, logger):
         # print('stamp draw2',process.memory_info().rss)
 
         if method == "phot":
+<<<<<<< photon_sampling
             # We already calculated realized_flux above.  Use that now and tell GalSim not
             # recalculate the Poisson realization of the flux.
             gal = gal.withFlux(self.realized_flux, bandpass)
+=======
+>>>>>>> main
             # print('stamp draw3b ',process.memory_info().rss)
 
             if not faint and "photon_ops" in config:
@@ -304,7 +295,11 @@ def draw(self, prof, image, method, offset, config, base, logger):
 
             # print('stamp draw3a',process.memory_info().rss)
             gal.drawImage(
+<<<<<<< photon_sampling
                 bandpass,
+=======
+                bandpass=bandpass,
+>>>>>>> main
                 method="phot",
                 offset=offset,
                 rng=self.rng,
@@ -314,7 +309,11 @@ def draw(self, prof, image, method, offset, config, base, logger):
                 photon_ops=photon_ops,
                 sensor=None,
                 add_to_image=True,
+<<<<<<< photon_sampling
                 poisson_flux=False,
+=======
+                poisson_flux=True,
+>>>>>>> main
             )
         else:
             fft_image = image.copy()
@@ -335,10 +334,9 @@ def draw(self, prof, image, method, offset, config, base, logger):
                 )
 
             # Go back to a combined convolution for fft drawing.
-            gal = gal.withFlux(self.flux, bandpass)
             prof = galsim.Convolve([gal] + psfs)
             try:
-                prof.drawImage(bandpass, **kwargs)
+                prof.drawImage(bandpass=bandpass, **kwargs)
             except galsim.errors.GalSimFFTSizeError as e:
                 # I think this shouldn't happen with the updates I made to how the image size
                 # is calculated, even for extremely bright things.  So it should be ok to
@@ -357,6 +355,7 @@ def draw(self, prof, image, method, offset, config, base, logger):
             # with FFT because we switched from phot to fft above.
             if self.use_fft_bright:
                 self.add_poisson_noise(fft_image)
+<<<<<<< photon_sampling
             # In case we had to make a bigger image, just copy the part we need.
             image += fft_image[image.bounds]
         # print('stamp draw3',process.memory_info().rss)
@@ -759,9 +758,17 @@ def draw(self, prof, image, method, offset, config, base, logger):
                 photons.flux=1
             else:
                 photons = galsim.PhotonArray(0)
+=======
+            # In case we had to make a bigger image, just copy the part we need.
+            image += fft_image[image.bounds]
+>>>>>>> main
         # print('stamp draw3',process.memory_info().rss)
         return photons
 
+    def add_poisson_noise(self, fft_image):
+        fft_image.array[fft_image.array < 0] = 0.0
+        fft_image.addNoise(galsim.PoissonNoise(rng=self.rng))
+
 
 # Pick the right function to be _fix_seds.
 if galsim.__version_info__ < (2, 5):

setup.py

@@ -7,7 +7,7 @@
     license="MIT",
     author="Michael Troxel",
     author_email="michael.troxel@duke.edu",
-    url="",
+    url="https://github.com/DukeCosmology/roman_imsim",
     packages=["roman_imsim"],
-    install_requires=["galsim", "fitsio", "astropy", "healpy", "numpy", "matplotlib", "scipy"],
+    install_requires=["galsim", "fitsio", "astropy", "numpy"],
 )