Merge branch '#670' into #628

CHANGELOG.md

@@ -108,7 +108,10 @@ Bug Fixes and Improvements
 - Fixed a bug in rendering profiles that involve two separate shifts. (#645)
 - Fixed a bug if drawImage was given odd nx, ny parameters, the drawn profile
   was not correctly centered in the image. (#645)
-
+- Added intermediate results cache to `ChromaticObject.drawImage` and
+  `ChromaticConvolution.drawImage`, which, for some applications, can
+  significantly speed up (anywhere from 10% to 2000%) the rendering of groups
+  of similar (same SED, same Bandpass, same PSF) chromatic profiles.
 
 Updates to config options
 -------------------------

galsim/chromatic.py

@@ -84,18 +84,21 @@ class ChromaticObject(object):
     #    profile at a particular wavelength.
     # 2) Define a `withScaledFlux` method, which scales the flux at all wavelengths by a fixed
     #    multiplier.
-    # 3) Initialize a `separable` attribute.  This marks whether (`separable = True`) or not
+    # 3) Potentially define their own `__repr__` and `__str__` methods.  Note that the default
+    #    assumes that `.obj` is the only attribute of significance, but this isn't always
+    #    appropriate, (e.g. ChromaticSum, ChromaticConvolution).
+    # 4) Initialize a `separable` attribute.  This marks whether (`separable = True`) or not
     #    (`separable = False`) the given chromatic profile can be factored into a spatial profile
     #    and a spectral profile.  Separable profiles can be drawn quickly by evaluating at a single
     #    wavelength and adjusting the flux via a (fast) 1D integral over the spectral component.
     #    Inseparable profiles, on the other hand, need to be evaluated at multiple wavelengths
     #    in order to draw (slow).
-    # 4) Separable objects must initialize an `SED` attribute, which is a callable object (often a
+    # 5) Separable objects must initialize an `SED` attribute, which is a callable object (often a
     #    `galsim.SED` instance) that returns the _relative_ flux of the profile at a given
     #    wavelength. (The _absolute_ flux is controlled by both the `SED` and the `.flux` attribute
     #    of the underlying chromaticized GSObject(s).  See `galsim.Chromatic` docstring for details
     #    concerning normalization.)
-    # 5) Initialize a `wave_list` attribute, which specifies wavelengths at which the profile (or
+    # 6) Initialize a `wave_list` attribute, which specifies wavelengths at which the profile (or
     #    the SED in the case of separable profiles) will be evaluated when drawing a
     #    ChromaticObject.  The type of `wave_list` should be a numpy array, and may be empty, in
     #    which case either the Bandpass object being drawn against, or the integrator being used
@@ -105,6 +108,7 @@ class ChromaticObject(object):
     # attribute representing a manipulated `GSObject` or `ChromaticObject`, or an `objlist`
     # attribute in the case of compound classes like `ChromaticSum` and `ChromaticConvolution`.
 
+
     def __init__(self, obj):
         self.obj = obj
         if isinstance(obj, galsim.GSObject):
@@ -119,6 +123,25 @@ def __init__(self, obj):
             raise TypeError("Can only directly instantiate ChromaticObject with a GSObject "+
                             "or ChromaticObject argument.")
 
+    @staticmethod
+    def _get_multiplier(sed, bandpass, wave_list):
+        wave_list = np.array(wave_list)
+        if len(wave_list) > 0:
+            multiplier = np.trapz(sed(wave_list) * bandpass(wave_list), wave_list)
+        else:
+            multiplier = galsim.integ.int1d(lambda w: sed(w) * bandpass(w),
+                                            bandpass.blue_limit, bandpass.red_limit)
+        return multiplier
+
+    @staticmethod
+    def resize_multiplier_cache(maxsize):
+        """ Resize the cache (default size=10) containing the integral over the product of an SED
+        and a Bandpass, which is used by ChromaticObject.drawImage().
+
+        @param maxsize  The new number of products to cache.
+        """
+        ChromaticObject._multiplier_cache.resize(maxsize)
+
     def __repr__(self):
         return 'galsim.ChromaticObject(%r)'%self.obj
 
@@ -225,12 +248,19 @@ def drawImage(self, bandpass, image=None, integrator='trapezoidal', **kwargs):
 
         By default, the above two integrators will use the trapezoidal rule for integration.  The
         midpoint rule for integration can be specified instead by passing an integrator that has
-        been initialized with the `rule=galsim.integ.midpt` argument.  Finally, when creating a
+        been initialized with the `rule=galsim.integ.midpt` argument.  When creating a
         ContinuousIntegrator, the number of samples `N` is also an argument.  For example:
 
             >>> integrator = galsim.ContinuousIntegrator(rule=galsim.integ.midpt, N=100)
             >>> image = chromatic_obj.drawImage(bandpass, integrator=integrator)
 
+        Finally, this method uses a cache to avoid recomputing the integral over the product of
+        the bandpass and object SED when possible (i.e., for separable profiles).  Because the
+        cache size is finite, users may find that it is more efficient when drawing many images
+        to group images using the same SEDs and bandpasses together in order to hit the cache more
+        often.  The default cache size is 10, but may be resized using the
+        `ChromaticObject.resize_multiplier_cache()` method.
+
         @param bandpass         A Bandpass object representing the filter against which to
                                 integrate.
         @param image            Optionally, the Image to draw onto.  (See GSObject.drawImage()
@@ -265,11 +295,7 @@ def drawImage(self, bandpass, image=None, integrator='trapezoidal', **kwargs):
         wave_list = self._getCombinedWaveList(bandpass)
 
         if self.separable:
-            if len(wave_list) > 0:
-                multiplier = np.trapz(self.SED(wave_list) * bandpass(wave_list), wave_list)
-            else:
-                multiplier = galsim.integ.int1d(lambda w: self.SED(w) * bandpass(w),
-                                                bandpass.blue_limit, bandpass.red_limit)
+            multiplier = ChromaticObject._multiplier_cache(self.SED, bandpass, tuple(wave_list))
             prof0 *= multiplier/self.SED(bandpass.effective_wavelength)
             image = prof0.drawImage(image=image, **kwargs)
             return image
@@ -712,6 +738,9 @@ def offset_func(w):
 
         return galsim.Transform(self, offset=offset)
 
+ChromaticObject._multiplier_cache = galsim.utilities.LRU_Cache(
+    ChromaticObject._get_multiplier, maxsize=10)
+
 
 class InterpolatedChromaticObject(ChromaticObject):
     """A ChromaticObject that uses interpolation of predrawn images to speed up subsequent
@@ -788,8 +817,8 @@ def __init__(self, obj, waves, oversample_fac=1.0):
 
         # Finally, now that we have an image scale and size, draw all the images.  Note that
         # `no_pixel` is used (we want the object on its own, without a pixel response).
-        self.ims = [ obj.drawImage(scale=scale, nx=im_size, ny=im_size, method='no_pixel') \
-                         for obj in objs ]
+        self.ims = [ obj.drawImage(scale=scale, nx=im_size, ny=im_size, method='no_pixel')
+                     for obj in objs ]
 
     def __repr__(self):
         s = 'galsim.InterpolatedChromaticObject(%r,%r'%(self.original, self.waves)
@@ -1648,6 +1677,48 @@ def __init__(self, *args, **kwargs):
         for obj in self.objlist:
             self.wave_list = np.union1d(self.wave_list, obj.wave_list)
 
+    @staticmethod
+    def _get_effective_prof(sep_SED, insep_profs, bandpass,
+                            iimult, wave_list, wmult, integrator,
+                            gsparams):
+            # Collapse inseparable profiles into one effective profile
+            SED = lambda w: reduce(lambda x,y:x*y, [s(w) for s in sep_SED], 1)
+            insep_obj = galsim.Convolve(insep_profs, gsparams=gsparams)
+            # Find scale at which to draw effective profile
+            iiscale = insep_obj.evaluateAtWavelength(bandpass.effective_wavelength).nyquistScale()
+            if iimult is not None:
+                iiscale /= iimult
+            # Create the effective bandpass.
+            wave_list = np.union1d(wave_list, bandpass.wave_list)
+            wave_list = wave_list[wave_list >= bandpass.blue_limit]
+            wave_list = wave_list[wave_list <= bandpass.red_limit]
+            effective_bandpass = galsim.Bandpass(
+                galsim.LookupTable(wave_list, bandpass(wave_list) * SED(wave_list),
+                                   interpolant='linear'))
+            # If there's only one inseparable profile, let it draw itself.
+            if len(insep_profs) == 1:
+                effective_prof_image = insep_profs[0].drawImage(
+                    effective_bandpass, wmult=wmult, scale=iiscale, integrator=integrator,
+                    method='no_pixel')
+            # Otherwise, use superclass ChromaticObject to draw convolution of inseparable profiles.
+            else:
+                effective_prof_image = ChromaticObject.drawImage(
+                        insep_obj, effective_bandpass, wmult=wmult, scale=iiscale,
+                        integrator=integrator, method='no_pixel')
+
+            effective_prof = galsim.InterpolatedImage(effective_prof_image, gsparams=gsparams)
+            return effective_prof
+
+    @staticmethod
+    def resize_effective_prof_cache(maxsize):
+        """ Resize the cache containing effective profiles, (i.e., wavelength-integrated products
+        of separable profile SEDs, inseparable profiles, and Bandpasses), which are used by
+        ChromaticConvolution.drawImage().
+
+        @param maxsize  The new number of effective profiles to cache.
+        """
+        ChromaticConvolution._effective_prof_cache.resize(maxsize)
+
     def _findSED(self):
         # pull out the non-trivial seds
         sedlist = [ obj.SED for obj in self.objlist if obj.SED != galsim.SED('1') ]
@@ -1692,6 +1763,14 @@ def drawImage(self, bandpass, image=None, integrator='trapezoidal', iimult=None,
         Works by finding sums of profiles which include separable portions, which can then be
         integrated before doing any convolutions, which are pushed to the end.
 
+        This method uses a cache to avoid recomputing 'effective' profiles, which are the
+        wavelength-integrated products of inseparable profiles, the spectral components of
+        separable profiles, and the bandpass.  Because the cache size is finite, users may find
+        that it is more efficient when drawing many images to group images using the same
+        SEDs, bandpasses, and inseparable profiles (generally PSFs) together in order to hit the
+        cache more often.  The default cache size is 10, but may be resized using the
+        `ChromaticConvolution.resize_effective_prof_cache()` method.
+
         @param bandpass         A Bandpass object representing the filter against which to
                                 integrate.
         @param image            Optionally, the Image to draw onto.  (See GSObject.drawImage()
@@ -1810,42 +1889,20 @@ def drawImage(self, bandpass, image=None, integrator='trapezoidal', iimult=None,
         # insep_profs should never be empty, since separable cases were farmed out to
         # ChromaticObject.drawImage() above.
 
-        # Collapse inseparable profiles into one effective profile
-        SED = lambda w: reduce(lambda x,y:x*y, [s(w) for s in sep_SED], 1)
-        insep_obj = galsim.Convolve(insep_profs, gsparams=self.gsparams)
-        # Find scale at which to draw effective profile
-        iiscale = insep_obj.evaluateAtWavelength(bandpass.effective_wavelength).nyquistScale()
-        if iimult is not None:
-            iiscale /= iimult
-        # Create the effective bandpass.
-        wave_list = np.union1d(wave_list, bandpass.wave_list)
-        wave_list = wave_list[wave_list >= bandpass.blue_limit]
-        wave_list = wave_list[wave_list <= bandpass.red_limit]
-        effective_bandpass = galsim.Bandpass(
-            galsim.LookupTable(wave_list, bandpass(wave_list) * SED(wave_list),
-                               interpolant='linear'))
-        # If there's only one inseparable profile, let it draw itself.
         wmult = kwargs.get('wmult', 1)
-        if len(insep_profs) == 1:
-            effective_prof_image = insep_profs[0].drawImage(
-                effective_bandpass, wmult=wmult, scale=iiscale, integrator=integrator,
-                method='no_pixel')
-        # Otherwise, use superclass ChromaticObject to draw convolution of inseparable profiles.
-        else:
-            effective_prof_image = ChromaticObject.drawImage(
-                    insep_obj, effective_bandpass, wmult=wmult, scale=iiscale,
-                    integrator=integrator, method='no_pixel')
-
-        # Image -> InterpolatedImage
-        # It could be useful to cache this result if drawing more than one object with the same
-        # PSF+SED combination.  This naturally happens in a ring test or when fitting the
-        # parameters of a galaxy profile to an image when the PSF is constant.
-        effective_prof = galsim.InterpolatedImage(effective_prof_image, gsparams=self.gsparams)
+        # Collapse inseparable profiles into one effective profile
+        effective_prof = ChromaticConvolution._effective_prof_cache(
+            tuple(sep_SED), tuple(insep_profs), bandpass, iimult, tuple(wave_list),
+            wmult, integrator, self.gsparams)
+
         # append effective profile to separable profiles (which should all be GSObjects)
         sep_profs.append(effective_prof)
         # finally, convolve and draw.
         final_prof = galsim.Convolve(sep_profs, gsparams=self.gsparams)
-        return final_prof.drawImage(image=image,**kwargs)
+        return final_prof.drawImage(image=image, **kwargs)
+
+ChromaticConvolution._effective_prof_cache = galsim.utilities.LRU_Cache(
+    ChromaticConvolution._get_effective_prof, maxsize=10)
 
 
 class ChromaticDeconvolution(ChromaticObject):

galsim/utilities.py

@@ -417,3 +417,100 @@ def _gammafn(x):
                0.00000000000000122678, -0.00000000000000011813, 0.00000000000000000119,
                0.00000000000000000141, -0.00000000000000000023, 0.00000000000000000002
              )
+
+
+class LRU_Cache:
+    """ Simplified Least Recently Used Cache.
+    Mostly stolen from http://code.activestate.com/recipes/577970-simplified-lru-cache/,
+    but added a method for dynamic resizing.  The least recently used cached item is
+    overwritten on a cache miss.
+
+    @param user_function   A python function to cache.
+    @param maxsize         Maximum number of inputs to cache.  [Default: 1024]
+
+    Usage
+    -----
+    >>> def slow_function(*args) # A slow-to-evaluate python function
+    >>>    ...
+    >>>
+    >>> v1 = slow_function(*k1)  # Calling function is slow
+    >>> v1 = slow_function(*k1)  # Calling again with same args is still slow
+    >>> cache = galsim.utilities.LRU_Cache(slow_function)
+    >>> v1 = cache(*k1)  # Returns slow_function(*k1), slowly the first time
+    >>> v1 = cache(*k1)  # Returns slow_function(*k1) again, but fast this time.
+
+    Methods
+    -------
+    >>> cache.resize(maxsize) # Resize the cache, either upwards or downwards.  Upwards resizing
+                              # is non-destructive.  Downwards resizing will remove the least
+                              # recently used items first.
+    """
+    def __init__(self, user_function, maxsize=1024):
+        # Link layout:     [PREV, NEXT, KEY, RESULT]
+        self.root = root = [None, None, None, None]
+        self.user_function = user_function
+        self.cache = cache = {}
+
+        last = root
+        for i in range(maxsize):
+            key = object()
+            cache[key] = last[1] = last = [last, root, key, None]
+        root[0] = last
+
+    def __call__(self, *key):
+        cache = self.cache
+        root = self.root
+        link = cache.get(key)
+        if link is not None:
+            # Cache hit: move link to last position
+            link_prev, link_next, _, result = link
+            link_prev[1] = link_next
+            link_next[0] = link_prev
+            last = root[0]
+            last[1] = root[0] = link
+            link[0] = last
+            link[1] = root
+            return result
+        # Cache miss: evaluate and insert new key/value at root, then increment root
+        #             so that just-evaluated value is in last position.
+        result = self.user_function(*key)
+        root[2] = key
+        root[3] = result
+        oldroot = root
+        root = self.root = root[1]
+        root[2], oldkey = None, root[2]
+        root[3], oldvalue = None, root[3]
+        del cache[oldkey]
+        cache[key] = oldroot
+        return result
+
+    def resize(self, maxsize):
+        """ Resize the cache.  Increasing the size of the cache is non-destructive, i.e.,
+        previously cached inputs remain in the cache.  Decreasing the size of the cache will
+        necessarily remove items from the cache if the cache is already filled.  Items are removed
+        in least recently used order.
+
+        @param maxsize  The new maximum number of inputs to cache.
+        """
+        oldsize = len(self.cache)
+        if maxsize == oldsize:
+            return
+        else:
+            root = self.root
+            cache = self.cache
+            if maxsize < oldsize:
+                for i in range(oldsize - maxsize):
+                    # Delete root.next
+                    current_next_link = root[1]
+                    new_next_link = root[1] = root[1][1]
+                    new_next_link[0] = root
+                    del cache[current_next_link[2]]
+            elif maxsize > oldsize:
+                for i in range(maxsize - oldsize):
+                    # Insert between root and root.next
+                    key = object()
+                    cache[key] = link = [root, root[1], key, None]
+                    root[1][0] = link
+                    root[1] = link
+            else:
+                raise ValueError("Invalid maxsize: {0:}".format(maxsize))

tests/test_chromatic.py

@@ -1313,6 +1313,9 @@ def evaluateAtWavelength(self, wave):
             ret = ret.shear(g1 = this_shear)
             return ret
 
+        def __repr__(self):
+            return 'galsim.ChromaticGaussian(%r)'%self.sigma
+
     # For this test, we're going to use the ChromaticGaussian defined above.  This class is simple
     # enough that evaluation of both the interpolated and exact versions is very fast, so it won't
     # slow down the tests too much to do both ways.  Note that for initial tests (fair comparison