Add makePhot method that is like drawPhot, but without the image

galsim/gsobject.py

@@ -2067,6 +2067,101 @@ def _calculate_nphotons(self, n_photons, poisson_flux, max_extra_noise, rng):
         return iN, g
 
 
+    def makePhot(self, n_photons=0, rng=None, max_extra_noise=0., poisson_flux=None,
+                 sensor=None, surface_ops=(), maxN=None, orig_center=PositionI(0,0),
+                 local_wcs=None):
+        """
+        Make photons for a profile.
+
+        This is equivalent to drawPhot, except that the photons are not placed onto
+        an image.  Instead, it just returns the PhotonArray.
+
+        .. note::
+
+            The (x,y) positions returned are in the same units as the distance units
+            of the GSObject being rendered.  If you want (x,y) in pixel coordinates, you
+            should call this function for the profile in image coordinates::
+
+                >>> photons = image.wcs.toImage(obj).makePhot()
+
+            Or if you just want a simple pixel scale conversion from sky coordinates to image
+            coordinates, you can instead do
+
+                >>> photons = obj.makePhot()
+                >>> photons.scaleXY(1./pixel_scale)
+
+        Parameters:
+            n_photons:      If provided, the number of photons to use for photon shooting.
+                            If not provided (i.e. ``n_photons = 0``), use as many photons as
+                            necessary to result in an image with the correct Poisson shot
+                            noise for the object's flux.  For positive definite profiles, this
+                            is equivalent to ``n_photons = flux``.  However, some profiles need
+                            more than this because some of the shot photons are negative
+                            (usually due to interpolants).  [default: 0]
+            rng:            If provided, a random number generator to use for photon shooting,
+                            which may be any kind of `BaseDeviate` object.  If ``rng`` is None, one
+                            will be automatically created, using the time as a seed.
+                            [default: None]
+            max_extra_noise: If provided, the allowed extra noise in each pixel when photon
+                            shooting.  This is only relevant if ``n_photons=0``, so the number of
+                            photons is being automatically calculated.  In that case, if the image
+                            noise is dominated by the sky background, then you can get away with
+                            using fewer shot photons than the full ``n_photons = flux``.
+                            Essentially each shot photon can have a ``flux > 1``, which increases
+                            the noise in each pixel.  The ``max_extra_noise`` parameter specifies
+                            how much extra noise per pixel is allowed because of this approximation.
+                            A typical value for this might be ``max_extra_noise = sky_level / 100``
+                            where ``sky_level`` is the flux per pixel due to the sky.  Note that
+                            this uses a "variance" definition of noise, not a "sigma" definition.
+                            [default: 0.]
+            poisson_flux:   Whether to allow total object flux scaling to vary according to
+                            Poisson statistics for ``n_photons`` samples when photon shooting.
+                            [default: True, unless ``n_photons`` is given, in which case the default
+                            is False]
+            surface_ops:    A list of operators that can modify the photon array that will be
+                            applied in order before accumulating the photons on the sensor.
+                            [default: ()]
+            orig_center:    The position of the image center in the original image coordinates.
+                            [default: (0,0)]
+            local_wcs:      The local wcs in the original image. [default: None]
+
+        Returns:
+            - a `PhotonArray` with the data about the photons.
+        """
+        from .sensor import Sensor
+        # Make sure the type of n_photons is correct and has a valid value:
+        if n_photons < 0.:
+            raise GalSimRangeError("Invalid n_photons < 0.", n_photons, 0., None)
+
+        if poisson_flux is None:
+            # If n_photons is given, poisson_flux = False
+            poisson_flux = (n_photons == 0.)
+
+        # Check that either n_photons is set to something or flux is set to something
+        if n_photons == 0. and self.flux == 1.:
+            galsim_warn(
+                    "Warning: drawImage for object with flux == 1, area == 1, and "
+                    "exptime == 1, but n_photons == 0.  This will only shoot a single photon.")
+
+        Ntot, g = self._calculate_nphotons(n_photons, poisson_flux, max_extra_noise, rng)
+
+        try:
+            photons = self.shoot(Ntot, rng)
+        except (GalSimError, NotImplementedError) as e:
+            raise GalSimNotImplementedError(
+                    "Unable to draw this GSObject with photon shooting.  Perhaps it "
+                    "is a Deconvolve or is a compound including one or more "
+                    "Deconvolve objects.\nOriginal error: %r"%(e))
+
+        if g != 1.:
+            photons.scaleFlux(g)
+
+        for op in surface_ops:
+            op.applyTo(photons, local_wcs)
+
+        return photons
+
+
     def drawPhot(self, image, gain=1., add_to_image=False,
                  n_photons=0, rng=None, max_extra_noise=0., poisson_flux=None,
                  sensor=None, surface_ops=(), maxN=None, orig_center=PositionI(0,0),
@@ -2145,8 +2240,8 @@ def drawPhot(self, image, gain=1., add_to_image=False,
             raise GalSimRangeError("Invalid n_photons < 0.", n_photons, 0., None)
 
         if poisson_flux is None:
-            if n_photons == 0.: poisson_flux = True
-            else: poisson_flux = False
+            # If n_photons is given, poisson_flux = False
+            poisson_flux = (n_photons == 0.)
 
         # Check that either n_photons is set to something or flux is set to something
         if n_photons == 0. and self.flux == 1.:

tests/test_airy.py

@@ -233,22 +233,26 @@ def test_airy_shoot():
     obj = galsim.Airy(lam_over_diam=0.01, flux=1.e4)
     im = galsim.Image(1000,1000, scale=1)
     im.setCenter(0,0)
-    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng)
+    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng.duplicate())
     print('obj.flux = ',obj.flux)
     print('added_flux = ',added_flux)
     print('photon fluxes = ',photons.flux.min(),'..',photons.flux.max())
     print('image flux = ',im.array.sum())
     assert np.isclose(added_flux, obj.flux)
     assert np.isclose(im.array.sum(), obj.flux)
+    photons2 = obj.makePhot(poisson_flux=False, rng=rng)
+    assert photons2 == photons, "Airy makePhot not equivalent to drawPhot"
 
     obj = galsim.Airy(lam_over_diam=0.01, flux=1.e4, obscuration=0.4)
-    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng)
+    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng.duplicate())
     print('obj.flux = ',obj.flux)
     print('added_flux = ',added_flux)
     print('photon fluxes = ',photons.flux.min(),'..',photons.flux.max())
     print('image flux = ',im.array.sum())
     assert np.isclose(added_flux, obj.flux)
     assert np.isclose(im.array.sum(), obj.flux)
+    photons2 = obj.makePhot(poisson_flux=False, rng=rng)
+    assert photons2 == photons, "Airy makePhot not equivalent to drawPhot"
 
 
 @timer

tests/test_box.py

@@ -250,31 +250,37 @@ def test_box_shoot():
     obj = galsim.Box(width=1.3, height=2.4, flux=1.e4)
     im = galsim.Image(100,100, scale=1)
     im.setCenter(0,0)
-    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng)
+    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng.duplicate())
     print('obj.flux = ',obj.flux)
     print('added_flux = ',added_flux)
     print('photon fluxes = ',photons.flux.min(),'..',photons.flux.max())
     print('image flux = ',im.array.sum())
     assert np.isclose(added_flux, obj.flux)
     assert np.isclose(im.array.sum(), obj.flux)
+    photons2 = obj.makePhot(poisson_flux=False, rng=rng)
+    assert photons2 == photons, "Box makePhot not equivalent to drawPhot"
 
     obj = galsim.Pixel(scale=9.3, flux=1.e4)
-    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng)
+    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng.duplicate())
     print('obj.flux = ',obj.flux)
     print('added_flux = ',added_flux)
     print('photon fluxes = ',photons.flux.min(),'..',photons.flux.max())
     print('image flux = ',im.array.sum())
     assert np.isclose(added_flux, obj.flux)
     assert np.isclose(im.array.sum(), obj.flux)
+    photons2 = obj.makePhot(poisson_flux=False, rng=rng)
+    assert photons2 == photons, "Pixel makePhot not equivalent to drawPhot"
 
     obj = galsim.TopHat(radius=4.7, flux=1.e4)
-    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng)
+    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng.duplicate())
     print('obj.flux = ',obj.flux)
     print('added_flux = ',added_flux)
     print('photon fluxes = ',photons.flux.min(),'..',photons.flux.max())
     print('image flux = ',im.array.sum())
     assert np.isclose(added_flux, obj.flux)
     assert np.isclose(im.array.sum(), obj.flux)
+    photons2 = obj.makePhot(poisson_flux=False, rng=rng)
+    assert photons2 == photons, "TopHat makePhot not equivalent to drawPhot"
 
 
 @timer

tests/test_convolve.py

@@ -152,6 +152,7 @@ def test_convolve():
     assert_raises(galsim.GalSimError, deconv.xValue, galsim.PositionD(0,0))
     assert_raises(galsim.GalSimError, deconv.drawReal, myImg)
     assert_raises(galsim.GalSimError, deconv.drawPhot, myImg, n_photons=10)
+    assert_raises(galsim.GalSimError, deconv.makePhot, n_photons=10)
 
 
 @timer

tests/test_draw.py

@@ -1125,9 +1125,17 @@ def test_shoot():
     np.testing.assert_equal(image4.array, 0)
 
     # Warns if flux is 1 and n_photons not given.
+    psf = galsim.Gaussian(sigma=3)
     with assert_warns(galsim.GalSimWarning):
-        psf = galsim.Gaussian(sigma=3)
         psf.drawImage(method='phot')
+    with assert_warns(galsim.GalSimWarning):
+        psf.drawPhot(image4)
+    with assert_warns(galsim.GalSimWarning):
+        psf.makePhot()
+    # With n_photons=1, it's fine.
+    psf.drawImage(method='phot', n_photons=1)
+    psf.drawPhot(image4, n_photons=1)
+    psf.makePhot(n_photons=1)
 
     # Check negative flux shooting with poisson_flux=True
     # The do_shoot test in galsim_test_helpers checks negative flux with a fixed number of photons.
@@ -1187,6 +1195,16 @@ def test_drawImage_area_exptime():
     # But not if we explicitly tell it to shoot 1 photon
     with assert_raises(AssertionError):
         assert_warns(galsim.GalSimWarning, obj1.drawImage, method='phot', n_photons=1)
+    # Likewise for makePhot
+    with assert_warns(galsim.GalSimWarning):
+        obj1.makePhot()
+    with assert_raises(AssertionError):
+        assert_warns(galsim.GalSimWarning, obj1.makePhot, n_photons=1)
+    # And drawPhot
+    with assert_warns(galsim.GalSimWarning):
+        obj1.drawPhot(im1)
+    with assert_raises(AssertionError):
+        assert_warns(galsim.GalSimWarning, obj1.drawPhot, im1, n_photons=1)
 
 
 @timer
@@ -1601,13 +1619,25 @@ def test_direct_scale():
                                    "drawReal made different image when off-center")
 
     obj.drawImage(im1, method='phot', rng=rng.duplicate())
-    obj.drawPhot(im2, rng=rng.duplicate())
-    obj.drawPhot(im3, rng=rng.duplicate())
+    _, phot1 = obj.drawPhot(im2, rng=rng.duplicate())
+    _, phot2 = obj.drawPhot(im3, rng=rng.duplicate())
+    phot3 = obj.makePhot(rng=rng.duplicate())
+    phot3.scaleXY(1./scale)
+    phot4 = im3.wcs.toImage(obj).makePhot(rng=rng.duplicate())
     print('phot: max diff = ',np.max(np.abs(im1.array - im2.array)))
     np.testing.assert_almost_equal(im1.array, im2.array, 15,
                                    "drawPhot made different image than method='phot'")
     np.testing.assert_almost_equal(im3.array, im2[im3.bounds].array, 15,
                                    "drawPhot made different image when off-center")
+    assert phot2 == phot1, "drawPhot made different photons than method='phot'"
+    assert phot3 == phot1, "makePhot made different photons than method='phot'"
+    # phot4 has a different order of operations for the math, so it doesn't come out exact.
+    np.testing.assert_almost_equal(phot4.x, phot3.x, 15,
+                                   "two ways to have makePhot apply scale have different x")
+    np.testing.assert_almost_equal(phot4.y, phot3.y, 15,
+                                   "two ways to have makePhot apply scale have different y")
+    np.testing.assert_almost_equal(phot4.flux, phot3.flux, 15,
+                                   "two ways to have makePhot apply scale have different flux")
 
     # Check images with invalid wcs raise ValueError
     im4 = galsim.ImageD(65, 65)
@@ -1621,6 +1651,7 @@ def test_direct_scale():
     # Also some other errors from drawPhot
     assert_raises(ValueError, obj.drawPhot, im2, n_photons=-20)
     assert_raises(TypeError, obj.drawPhot, im2, sensor=5)
+    assert_raises(ValueError, obj.makePhot, n_photons=-20)
 
 if __name__ == "__main__":
     test_drawImage()

tests/test_exponential.py

@@ -219,13 +219,15 @@ def test_exponential_shoot():
     obj = galsim.Exponential(half_light_radius=3.5, flux=1.e4)
     im = galsim.Image(100,100, scale=1)
     im.setCenter(0,0)
-    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng)
+    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng.duplicate())
     print('obj.flux = ',obj.flux)
     print('added_flux = ',added_flux)
     print('photon fluxes = ',photons.flux.min(),'..',photons.flux.max())
     print('image flux = ',im.array.sum())
     assert np.isclose(added_flux, obj.flux)
     assert np.isclose(im.array.sum(), obj.flux)
+    photons2 = obj.makePhot(poisson_flux=False, rng=rng)
+    assert photons2 == photons, "Exponential makePhot not equivalent to drawPhot"
 
 
 @timer

tests/test_gaussian.py

@@ -339,13 +339,15 @@ def test_gaussian_shoot():
     obj = galsim.Gaussian(fwhm=3.5, flux=1.e4)
     im = galsim.Image(100,100, scale=1)
     im.setCenter(0,0)
-    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng)
+    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng.duplicate())
     print('obj.flux = ',obj.flux)
     print('added_flux = ',added_flux)
     print('photon fluxes = ',photons.flux.min(),'..',photons.flux.max())
     print('image flux = ',im.array.sum())
     assert np.isclose(added_flux, obj.flux)
     assert np.isclose(im.array.sum(), obj.flux)
+    photons2 = obj.makePhot(poisson_flux=False, rng=rng)
+    assert photons2 == photons, "Gaussian makePhot not equivalent to drawPhot"
 
 
 @timer

tests/test_interpolatedimage.py

@@ -1515,7 +1515,7 @@ def test_ii_shoot():
         flux = 1.e4
     for interp in interp_list:
         obj = galsim.InterpolatedImage(image_in, x_interpolant=interp, scale=3.3, flux=flux)
-        added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng)
+        added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng.duplicate())
         print('interp = ',interp)
         print('obj.flux = ',obj.flux)
         print('added_flux = ',added_flux)
@@ -1533,6 +1533,8 @@ def test_ii_shoot():
             rtol = 1.e-7
         assert np.isclose(added_flux, obj.flux, rtol=rtol)
         assert np.isclose(im.array.sum(), obj.flux, rtol=rtol)
+        photons2 = obj.makePhot(poisson_flux=False, rng=rng)
+        assert photons2 == photons, "InterpolatedImage makePhot not equivalent to drawPhot"
 
 
 @timer

tests/test_kolmogorov.py

@@ -331,13 +331,15 @@ def test_kolmogorov_shoot():
     obj = galsim.Kolmogorov(fwhm=1.5, flux=1.e4)
     im = galsim.Image(500,500, scale=1)
     im.setCenter(0,0)
-    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng)
+    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng.duplicate())
     print('obj.flux = ',obj.flux)
     print('added_flux = ',added_flux)
     print('photon fluxes = ',photons.flux.min(),'..',photons.flux.max())
     print('image flux = ',im.array.sum())
     assert np.isclose(added_flux, obj.flux)
     assert np.isclose(im.array.sum(), obj.flux)
+    photons2 = obj.makePhot(poisson_flux=False, rng=rng)
+    assert photons2 == photons, "Kolmogorov makePhot not equivalent to drawPhot"
 
 
 @timer

tests/test_moffat.py

@@ -439,23 +439,27 @@ def test_moffat_shoot():
     obj = galsim.Moffat(fwhm=3.5, beta=4.7, flux=1.e4)
     im = galsim.Image(500,500, scale=1)
     im.setCenter(0,0)
-    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng)
+    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng.duplicate())
     print('obj.flux = ',obj.flux)
     print('added_flux = ',added_flux)
     print('photon fluxes = ',photons.flux.min(),'..',photons.flux.max())
     print('image flux = ',im.array.sum())
     assert np.isclose(added_flux, obj.flux)
     assert np.isclose(im.array.sum(), obj.flux)
+    photons2 = obj.makePhot(poisson_flux=False, rng=rng)
+    assert photons2 == photons, "Moffat makePhot not equivalent to drawPhot"
 
     # Note: low beta has large wings, so don't go too low.  Also, reduce fwhm a bit.
     obj = galsim.Moffat(fwhm=1.5, beta=1.9, flux=1.e4)
-    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng)
+    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng.duplicate())
     print('obj.flux = ',obj.flux)
     print('added_flux = ',added_flux)
     print('photon fluxes = ',photons.flux.min(),'..',photons.flux.max())
     print('image flux = ',im.array.sum())
     assert np.isclose(added_flux, obj.flux)
     assert np.isclose(im.array.sum(), obj.flux)
+    photons2 = obj.makePhot(poisson_flux=False, rng=rng)
+    assert photons2 == photons, "Moffat makePhot not equivalent to drawPhot"
 
 
 @timer

tests/test_photon_array.py

@@ -358,10 +358,18 @@ def applyTo(self, photon_array, local_wcs=None):
     rng2 = galsim.BaseDeviate(1234)
     sampler2 = galsim.WavelengthSampler(sed, bandpass, rng2)
     obj.drawImage(im2, method='phot', n_photons=nphotons, use_true_center=False,
-                  surface_ops=[sampler2,clip600], rng=rng2)
+                  surface_ops=[sampler2,clip600], rng=rng2, save_photons=True)
     print('sum = ',im1.array.sum(),im2.array.sum())
     np.testing.assert_array_equal(im1.array, im2.array)
 
+    # Equivalent version just getting photons back
+    rng2.seed(1234)
+    photons = obj.makePhot(n_photons=nphotons, surface_ops=[sampler2,clip600], rng=rng2)
+    print('phot.x = ',photons.x)
+    print('im2.photons.x = ',im2.photons.x)
+    assert photons == im2.photons
+
+
 @timer
 def test_photon_angles():
     """Test the photon_array function

tests/test_second_kick.py

@@ -265,13 +265,15 @@ def test_sk_shoot():
     obj = galsim.SecondKick(lam=500, r0=0.2, diam=4, flux=1.e4)
     im = galsim.Image(500,500, scale=1)
     im.setCenter(0,0)
-    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng)
+    added_flux, photons = obj.drawPhot(im, poisson_flux=False, rng=rng.duplicate())
     print('obj.flux = ',obj.flux)
     print('added_flux = ',added_flux)
     print('photon fluxes = ',photons.flux.min(),'..',photons.flux.max())
     print('image flux = ',im.array.sum())
     assert np.isclose(added_flux, obj.flux)
     assert np.isclose(im.array.sum(), obj.flux)
+    photons2 = obj.makePhot(poisson_flux=False, rng=rng)
+    assert photons2 == photons, "SecondKick makePhot not equivalent to drawPhot"
 
 
 @timer