Avoid gratuitous constructions of UniformDeviates in python

examples/fft_vs_geom_movie.py

@@ -228,7 +228,8 @@ def make_movie(args):
                     geom_psf = geom_psl.makePSF(
                             lam=args.lam, aper=geom_aper, t0=t0, exptime=args.time_step)
                     geom_img0 = geom_psf.drawImage(nx=args.nx, ny=args.nx, scale=scale,
-                                                   method='phot', n_photons=args.geom_nphot)
+                                                   method='phot', n_photons=args.geom_nphot,
+                                                   rng=rng)
 
                 t0 += args.time_step
 

galsim/convolve.py

@@ -400,18 +400,18 @@ def _drawReal(self, image):
             raise GalSimError("Cannot use real_space convolution for >2 profiles")
 
     @doc_inherit
-    def _shoot(self, photons, ud):
+    def _shoot(self, photons, rng):
         from .photon_array import PhotonArray
 
-        self.obj_list[0]._shoot(photons, ud)
+        self.obj_list[0]._shoot(photons, rng)
         # It may be necessary to shuffle when convolving because we do not have a
         # gaurantee that the convolvee's photons are uncorrelated, e.g., they might
         # both have their negative ones at the end.
         # However, this decision is now made by the convolve method.
         for obj in self.obj_list[1:]:
             p1 = PhotonArray(len(photons))
-            obj._shoot(p1, ud)
-            photons.convolve(p1, ud)
+            obj._shoot(p1, rng)
+            photons.convolve(p1, rng)
 
     @doc_inherit
     def _drawKImage(self, image):
@@ -774,12 +774,12 @@ def _prepareDraw(self):
         self.orig_obj._prepareDraw()
 
     @doc_inherit
-    def _shoot(self, photons, ud):
+    def _shoot(self, photons, rng):
         from .photon_array import PhotonArray
-        self.orig_obj._shoot(photons, ud)
+        self.orig_obj._shoot(photons, rng)
         photons2 = PhotonArray(len(photons))
-        self.orig_obj._shoot(photons2, ud)
-        photons.convolve(photons2, ud)
+        self.orig_obj._shoot(photons2, rng)
+        photons.convolve(photons2, rng)
 
 
 def AutoCorrelate(obj, real_space=None, gsparams=None, propagate_gsparams=True):
@@ -935,13 +935,13 @@ def _prepareDraw(self):
         self._orig_obj._prepareDraw()
 
     @doc_inherit
-    def _shoot(self, photons, ud):
+    def _shoot(self, photons, rng):
         from .photon_array import PhotonArray
-        self.orig_obj._shoot(photons, ud)
+        self.orig_obj._shoot(photons, rng)
         photons2 = PhotonArray(len(photons))
-        self.orig_obj._shoot(photons2, ud)
+        self.orig_obj._shoot(photons2, rng)
 
         # Flip sign of (x, y) in one of the results
         photons2.scaleXY(-1)
 
-        photons.convolve(photons2, ud)
+        photons.convolve(photons2, rng)

galsim/gsobject.py

@@ -270,7 +270,7 @@ def __init__(self):
     #     _xValue(self, pos)
     #     _kValue(self, kpos)
     #     _drawReal(self, image)
-    #     _shoot(self, photons, ud):
+    #     _shoot(self, photons, rng):
     #     _drawKImage(self, image)
     #
     # Required for real-space convolution
@@ -1518,7 +1518,6 @@ def drawImage(self, image=None, nx=None, ny=None, bounds=None, scale=None, wcs=N
         from .convolve import Convolve, Convolution, Deconvolve
         from .box import Pixel
         from .wcs import PixelScale
-        from .random import UniformDeviate
         from .photon_array import PhotonArray
 
         # Check that image is sane
@@ -1680,8 +1679,7 @@ def drawImage(self, image=None, nx=None, ny=None, bounds=None, scale=None, wcs=N
                 added_photons = prof.drawFFT(draw_image, add)
 
             if sensor is not None:
-                ud = UniformDeviate(rng)
-                photons = PhotonArray.makeFromImage(draw_image, rng=ud)
+                photons = PhotonArray.makeFromImage(draw_image, rng=rng)
                 for op in surface_ops:
                     op.applyTo(photons, local_wcs)
                 if imview.dtype in (np.float32, np.float64):
@@ -2079,7 +2077,6 @@ def drawPhot(self, image, gain=1., add_to_image=False,
             nphotons is the total flux of photons that landed inside the image bounds, and
             photons is the PhotonArray that was applied to the image.
         """
-        from .random import UniformDeviate
         from .sensor import Sensor
         from .image import ImageD
         # Make sure the type of n_photons is correct and has a valid value:
@@ -2096,8 +2093,6 @@ def drawPhot(self, image, gain=1., add_to_image=False,
                     "Warning: drawImage for object with flux == 1, area == 1, and "
                     "exptime == 1, but n_photons == 0.  This will only shoot a single photon.")
 
-        ud = UniformDeviate(rng)
-
         # Make sure the image is set up to have unit pixel scale and centered at 0,0.
         if image.wcs is None or not image.wcs.isPixelScale():
             raise GalSimValueError("drawPhot requires an image with a PixelScale wcs", image)
@@ -2107,7 +2102,7 @@ def drawPhot(self, image, gain=1., add_to_image=False,
         elif not isinstance(sensor, Sensor):
             raise TypeError("The sensor provided is not a Sensor instance")
 
-        Ntot, g = self._calculate_nphotons(n_photons, poisson_flux, max_extra_noise, ud)
+        Ntot, g = self._calculate_nphotons(n_photons, poisson_flux, max_extra_noise, rng)
 
         if gain != 1.:
             g /= gain
@@ -2129,7 +2124,7 @@ def drawPhot(self, image, gain=1., add_to_image=False,
             thisN = min(maxN, Nleft)
 
             try:
-                photons = self.shoot(thisN, ud)
+                photons = self.shoot(thisN, rng)
             except (GalSimError, NotImplementedError) as e:
                 raise GalSimNotImplementedError(
                         "Unable to draw this GSObject with photon shooting.  Perhaps it "
@@ -2170,23 +2165,24 @@ def shoot(self, n_photons, rng=None):
 
         @returns PhotonArray.
         """
-        from .random import UniformDeviate
+        from .random import BaseDeviate
         from .photon_array import PhotonArray
 
         photons = PhotonArray(n_photons)
         if n_photons == 0:
             # It's ok to shoot 0, but downstream can have problems with it, so just stop now.
             return photons
+        if rng is None:
+            rng = BaseDeviate()
 
-        ud = UniformDeviate(rng)
-        self._shoot(photons, ud)
+        self._shoot(photons, rng)
         return photons
 
-    def _shoot(self, photons, ud):
+    def _shoot(self, photons, rng):
         """Shoot photons into the given PhotonArray
 
         @param photons      A PhotonArray instance into which the photons should be placed.
-        @param ud           A UniformDeviate instance to use for the photon shooting,
+        @param rng          A BaseDeviate instance to use for the photon shooting,
         """
         raise NotImplementedError("%s does not implement shoot"%self.__class__.__name__)
 

galsim/interpolatedimage.py

@@ -685,9 +685,9 @@ def _kValue(self, kpos):
         return self._sbp.kValue(kpos._p)
 
     @doc_inherit
-    def _shoot(self, photons, ud):
+    def _shoot(self, photons, rng):
         with convert_cpp_errors():
-            self._sbp.shoot(photons._pa, ud._rng)
+            self._sbp.shoot(photons._pa, rng._rng)
 
     @doc_inherit
     def _drawReal(self, image):

galsim/phase_psf.py

@@ -1487,15 +1487,17 @@ def _drawReal(self, image):
         self._ii._drawReal(image)
 
     @doc_inherit
-    def _shoot(self, photons, ud):
+    def _shoot(self, photons, rng):
         from .photon_array import PhotonArray
+        from .random import UniformDeviate
 
         if not self._geometric_shooting:
             self._prepareDraw()
-            return self._ii._shoot(photons, ud)
+            return self._ii._shoot(photons, rng)
 
         n_photons = len(photons)
         t = np.empty((n_photons,), dtype=float)
+        ud = UniformDeviate(rng)
         ud.generate(t)
         t *= self.exptime
         t += self.t0
@@ -1519,8 +1521,8 @@ def _shoot(self, photons, ud):
 
         if self.second_kick:
             p2 = PhotonArray(len(photons))
-            self.second_kick._shoot(p2, ud)
-            photons.convolve(p2, ud)
+            self.second_kick._shoot(p2, rng)
+            photons.convolve(p2, rng)
 
     @doc_inherit
     def _drawKImage(self, image):
@@ -1956,8 +1958,8 @@ def _drawReal(self, image):
         self._psf._drawReal(image)
 
     @doc_inherit
-    def _shoot(self, photons, ud):
-        self._psf._shoot(photons, ud)
+    def _shoot(self, photons, rng):
+        self._psf._shoot(photons, rng)
 
     @doc_inherit
     def _drawKImage(self, image):

galsim/photon_array.py

@@ -23,7 +23,7 @@
 import numpy as np
 
 from . import _galsim
-from .random import UniformDeviate
+from .random import UniformDeviate, BaseDeviate
 from .utilities import lazy_property
 from .angle import radians, arcsec
 from .errors import GalSimError, GalSimRangeError, GalSimValueError, GalSimUndefinedBoundsError
@@ -217,8 +217,9 @@ def convolve(self, rhs, rng=None):
         if rhs.size() != self.size():
             raise GalSimIncompatibleValuesError("PhotonArray.convolve with unequal size arrays",
                                                 self_pa=self, rhs=rhs)
-        ud = UniformDeviate(rng)
-        self._pa.convolve(rhs._pa, ud._rng)
+        if rng is None:
+            rng = BaseDeviate()
+        self._pa.convolve(rhs._pa, rng._rng)
 
     def __repr__(self):
         s = "galsim.PhotonArray(%d, x=array(%r), y=array(%r), flux=array(%r)"%(
@@ -315,7 +316,6 @@ def makeFromImage(cls, image, max_flux=1., rng=None):
 
         @returns a PhotonArray
         """
-        ud = UniformDeviate(rng)
         max_flux = float(max_flux)
         if (max_flux <= 0):
             raise GalSimRangeError("max_flux must be positive", max_flux, 0.)
@@ -326,7 +326,9 @@ def makeFromImage(cls, image, max_flux=1., rng=None):
         N = int(np.prod(image.array.shape) + total_flux / max_flux)
         photons = cls(N)
 
-        N = photons._pa.setFrom(image._image, max_flux, ud._rng)
+        if rng is None:
+            rng = BaseDeviate()
+        N = photons._pa.setFrom(image._image, max_flux, rng._rng)
         photons._x = photons.x[:N]
         photons._y = photons.y[:N]
         photons._flux = photons.flux[:N]
@@ -415,6 +417,8 @@ class WavelengthSampler(object):
     def __init__(self, sed, bandpass, rng=None, npoints=None):
         self.sed = sed
         self.bandpass = bandpass
+        if rng is None:
+            rng = BaseDeviate()
         self.rng = rng
         self.npoints = npoints
 

galsim/randwalk.py

@@ -263,8 +263,7 @@ def _get_points(self):
 
         The most efficient way is to write into an image
         """
-        ud = UniformDeviate(self._rng)
-        photons = self._profile.shoot(self._npoints, ud)
+        photons = self._profile.shoot(self._npoints, self._rng)
         ar = np.column_stack([ photons.x, photons.y ])
 
         return ar
@@ -354,8 +353,8 @@ def _kValue(self, kpos):
         return self._sbp.kValue(kpos._p)
 
     @doc_inherit
-    def _shoot(self, photons, ud):
-        self._sbp.shoot(photons._pa, ud._rng)
+    def _shoot(self, photons, rng):
+        self._sbp.shoot(photons._pa, rng._rng)
 
     @doc_inherit
     def _drawKImage(self, image):

galsim/sum.py

@@ -303,7 +303,7 @@ def _drawReal(self, image):
                 image += im1
 
     @doc_inherit
-    def _shoot(self, photons, ud):
+    def _shoot(self, photons, rng):
         from .photon_array import PhotonArray
         from .random import BinomialDeviate
 
@@ -322,10 +322,10 @@ def _shoot(self, photons, ud):
             thisN = remainingN  # All of what's left, if this is the last summand...
             if i < len(self.obj_list)-1:
                 # otherwise, allocate a randomized fraction of the remaining photons to summand.
-                bd = BinomialDeviate(ud, remainingN, thisAbsoluteFlux/remainingAbsoluteFlux)
+                bd = BinomialDeviate(rng, remainingN, thisAbsoluteFlux/remainingAbsoluteFlux)
                 thisN = int(bd())
             if thisN > 0:
-                thisPA = obj.shoot(thisN, ud)
+                thisPA = obj.shoot(thisN, rng)
                 # Now rescale the photon fluxes so that they are each nominally fluxPerPhoton
                 # whereas the shoot() routine would have made them each nominally
                 # thisAbsoluteFlux/thisN

galsim/transform.py

@@ -490,8 +490,8 @@ def _drawReal(self, image):
             self._sbp.draw(image._image, image.scale)
 
     @doc_inherit
-    def _shoot(self, photons, ud):
-        self._original._shoot(photons, ud)
+    def _shoot(self, photons, rng):
+        self._original._shoot(photons, rng)
         photons.x, photons.y = self._fwd(photons.x, photons.y)
         photons.x += self.offset.x
         photons.y += self.offset.y

include/galsim/PhotonArray.h

@@ -207,9 +207,9 @@ namespace galsim {
          * totals, if the two photon streams are uncorrelated.
          *
          * @param[in] rhs PhotonArray to convolve with this one.  Must be same size.
-         * @param[in] ud  A UniformDeviate in case we need to shuffle.
+         * @param[in] rng  A BaseDeviate in case we need to shuffle.
          */
-        void convolve(const PhotonArray& rhs, UniformDeviate ud);
+        void convolve(const PhotonArray& rhs, BaseDeviate ud);
 
         /**
          * @brief Convolve this array with another, shuffling the order in which photons are
@@ -219,9 +219,9 @@ namespace galsim {
          * multiplied into the array is randomized to destroy any flux or position correlations.
          *
          * @param[in] rhs PhotonArray to convolve with this one.  Must be same size.
-         * @param[in] ud  A UniformDeviate used to shuffle the input photons.
+         * @param[in] rng  A BaseDeviate used to shuffle the input photons.
          */
-        void convolveShuffle(const PhotonArray& rhs, UniformDeviate ud);
+        void convolveShuffle(const PhotonArray& rhs, BaseDeviate rng);
 
         /**
          * @brief Add flux of photons to an image by binning into pixels.
@@ -250,12 +250,12 @@ namespace galsim {
          *
          * @param image     The image to use for the photon fluxes and positions.
          * @param maxFlux   The maximum flux that any photon should have.
-         * @param ud        A UniformDeviate in case we need to shuffle.
+         * @param rng       A BaseDeviate in case we need to shuffle.
          *
          * @returns the total number of photons set.
          */
         template <class T>
-        int setFrom(const BaseImage<T>& image, double maxFlux, UniformDeviate ud);
+        int setFrom(const BaseImage<T>& image, double maxFlux, BaseDeviate ud);
 
         /**
          * @brief Check if the current array has correlated photons.

include/galsim/SBProfile.h

@@ -295,9 +295,9 @@ namespace galsim {
          * alternative to rejection sampling.  See `OneDimensionalDeviate` documentation.
          *
          * @param[in] photons PhotonArray in which to write the photon information
-         * @param[in] ud UniformDeviate that will be used to draw photons from distribution.
+         * @param[in] rng BaseDeviate that will be used to draw photons from distribution.
          */
-        void shoot(PhotonArray& photons, UniformDeviate ud) const;
+        void shoot(PhotonArray& photons, BaseDeviate rng) const;
 
         /**
          * @brief Return expectation value of flux in positive photons when shoot() is called

include/galsim/Silicon.h

@@ -47,7 +47,7 @@ namespace galsim
         bool insidePixel(int ix, int iy, double x, double y, double zconv,
                          ImageView<T> target, bool* off_edge=0) const;
 
-        double calculateConversionDepth(const PhotonArray& photons, int i, UniformDeviate ud) const;
+        double calculateConversionDepth(const PhotonArray& photons, int i, BaseDeviate rng) const;
 
         template <typename T>
         void updatePixelDistortions(ImageView<T> target);
@@ -56,7 +56,7 @@ namespace galsim
         void addTreeRingDistortions(ImageView<T> target, Position<int> orig_center);
 
         template <typename T>
-        double accumulate(const PhotonArray& photons, UniformDeviate ud, ImageView<T> target,
+        double accumulate(const PhotonArray& photons, BaseDeviate rng, ImageView<T> target,
                           Position<int> orig_center, bool resume);
 
         template <typename T>

pysrc/PhotonArray.cpp

@@ -27,7 +27,7 @@ namespace galsim {
         wrapper
             .def("addTo", (double (PhotonArray::*)(ImageView<T>) const) &PhotonArray::addTo)
             .def("setFrom",
-                 (int (PhotonArray::*)(const BaseImage<T>&, double, UniformDeviate))
+                 (int (PhotonArray::*)(const BaseImage<T>&, double, BaseDeviate))
                  &PhotonArray::setFrom);
     }