Add bandpass option to SED.thin

galsim/sed.py

@@ -875,7 +875,8 @@ def calculateMagnitude(self, bandpass):
         flux = self.calculateFlux(bandpass)
         return -2.5 * np.log10(flux) + bandpass.zeropoint
 
-    def thin(self, rel_err=1.e-4, trim_zeros=True, preserve_range=True, fast_search=True):
+    def thin(self, rel_err=1.e-4, trim_zeros=True, preserve_range=True, fast_search=True,
+             bandpass=None):
         """Remove some tabulated values while keeping the integral over the set of tabulated values
         still accurate to ``rel_err``.
 
@@ -901,14 +902,58 @@ def thin(self, rel_err=1.e-4, trim_zeros=True, preserve_range=True, fast_search=
                               tends to yield a thinned representation that retains fewer samples
                               while still meeting the relative error requirement.
                               [default: True]
+            bandpass:         If set, then target the thinned SED to a particular bandpass.
+                              The integral being optimized will be the product of the SED
+                              and the bandpass throughput.  And if preserve_range is False,
+                              then the output SED will be truncated to the same limits as the
+                              bandpass limits. [default: None]
 
         Returns:
             the thinned `SED`.
         """
-        if len(self.wave_list) > 0:
-            rest_wave_native = self._get_rest_native_waves(self.wave_list)
+        def _bandpass_native_waves(waves):
+            # The thinning algorithm runs on the SED's native tabulation grid, but when a
+            # bandpass target is provided we need to evaluate the throughput at the same physical
+            # observed wavelengths.  This helper converts SED-native wavelengths into the native
+            # wavelength units used by bandpass._tp.
+            observed_nm = waves / self.wave_factor
+            observed_nm *= (1.0 + self.redshift)
+            return observed_nm * bandpass.wave_factor
+
+        if bandpass is not None:
+            if self.blue_limit > bandpass.red_limit or self.red_limit < bandpass.blue_limit:
+                raise GalSimError("Bandpass does not overlap the SED wavelength range.")
+            wave_list, _, _ = utilities.combine_wave_list(self, bandpass)
+            if preserve_range and not trim_zeros:
+                # If we want to preserve the range, add back the limits to each end.
+                front = [self.blue_limit] if self.blue_limit < wave_list[0] else []
+                back = [self.red_limit] if self.red_limit > wave_list[-1] else []
+                if front or back:
+                    wave_list = np.concatenate((front, wave_list, back))
+        else:
+            wave_list = self.wave_list
+
+        if len(wave_list) > 0:
+            rest_wave_native = self._get_rest_native_waves(wave_list)
             spec_native = self._spec(rest_wave_native)
 
+            if bandpass is not None:
+                # Identify the overlapping region in the SED's native wavelength units to
+                # determine which portion of the wave_list is within the bandpass limits.
+                band_native_limits = self._get_rest_native_waves(
+                        np.array([bandpass.blue_limit, bandpass.red_limit]))
+                native_blue_limit = np.min(band_native_limits)
+                native_red_limit = np.max(band_native_limits)
+                in_band = np.logical_and(rest_wave_native >= native_blue_limit,
+                                         rest_wave_native <= native_red_limit)
+
+                # Compute the product SED * bandpass, since that is the quantity whose integral we
+                # want to preserve for this observation.
+                tp_native = np.zeros_like(rest_wave_native, dtype=float)
+                bp_wave_native = _bandpass_native_waves(rest_wave_native[in_band])
+                tp_native[in_band] = bandpass._tp(bp_wave_native) / bandpass.wave_factor
+                spec_native *= tp_native
+
             # Note that this is thinning in native units, not nm and photons/nm.
             interpolant = (self.interpolant if not isinstance(self._spec, LookupTable)
                            else self._spec.interpolant)
@@ -917,6 +962,17 @@ def thin(self, rel_err=1.e-4, trim_zeros=True, preserve_range=True, fast_search=
                     trim_zeros=trim_zeros, preserve_range=preserve_range,
                     fast_search=fast_search, interpolant=interpolant)
 
+            if bandpass is not None:
+                # Convert the thinned product back into an SED by dividing out the bandpass
+                # wherever the throughput is non-zero.
+                in_band = np.logical_and(newx >= native_blue_limit, newx <= native_red_limit)
+                tp_native = np.zeros_like(newx, dtype=float)
+                bp_wave_native = _bandpass_native_waves(newx[in_band])
+                tp_native[in_band] = bandpass._tp(bp_wave_native) / bandpass.wave_factor
+                nz = tp_native != 0.  # Don't divide by 0.
+                assert np.all(newf[~nz] == 0.)
+                newf[nz] /= tp_native[nz]
+
             newspec = _LookupTable(newx, newf, interpolant=interpolant)
             return SED(newspec, self.wave_type, self.flux_type, redshift=self.redshift,
                        fast=self.fast)

tests/test_sed.py

@@ -1115,6 +1115,129 @@ def test_thin():
     print('s = ',s)
     np.testing.assert_equal(s.wave_list, [0,1])
 
+@timer
+def test_thin_bandpass():
+    # Check bandpass-targeted thinning on a synthetic case where we can make strong assertions
+    # about the intended semantics without depending on quirks of any particular throughput file.
+    sed_waves_nm = np.array(
+            [400, 450, 500, 510, 520, 540, 560, 580, 600, 620, 640, 660, 680, 700, 800, 900],
+            dtype=float)
+    sed_vals = np.array(
+            [0.2, 1.0, 3.0, 2.85, 2.7, 2.4, 2.2, 2.0, 1.8, 1.6, 1.4, 1.2, 1.0, 0.8, 0.4, 0.1],
+            dtype=float)
+    bp_waves_nm = np.array([500, 520, 540, 560, 580, 600, 620, 640, 660, 680, 700], dtype=float)
+    bp_vals = np.array([0.0, 2.e-4, 0.03, 0.25, 0.68, 1.0, 0.68, 0.25, 0.03, 2.e-4, 0.0],
+                       dtype=float)
+
+    for sed_wave_type, bp_wave_type in [('ang', 'nm'), ('nm', 'ang')]:
+        # Exercise both directions of the unit conversion logic in SED.thin:
+        # SED-native Angstroms with an nm bandpass, and vice versa.
+        sed_scale = 10. if sed_wave_type == 'ang' else 1.
+        bp_scale = 10. if bp_wave_type == 'ang' else 1.
+        sed_table = galsim.LookupTable(sed_waves_nm * sed_scale, sed_vals, interpolant='linear')
+        bp_table = galsim.LookupTable(bp_waves_nm * bp_scale, bp_vals, interpolant='linear')
+        s = galsim.SED(sed_table, sed_wave_type, 'fphotons')
+        bp = galsim.Bandpass(bp_table, bp_wave_type)
+        flux = s.calculateFlux(bp)
+        wave_list, _, _ = galsim.utilities.combine_wave_list(s, bp)
+        combined_in_band = wave_list[np.logical_and(wave_list >= bp.blue_limit,
+                                                    wave_list <= bp.red_limit)]
+        print('wave_types = ', sed_wave_type, bp_wave_type)
+        print('unthinned: ',s.wave_list)
+        print('in band: ',combined_in_band)
+
+        # Default preserve_range=True, trim_zeros=True keeps the full range of the
+        # bandpass, not the SED's original full range.
+        thin_s = s.thin(rel_err=1.e-3, bandpass=bp)
+        print('thinned: ',thin_s.wave_list)
+        thin_flux = thin_s.calculateFlux(bp)
+        thin_err = abs((thin_flux - flux) / flux)
+        assert thin_err < 1.e-3
+        np.testing.assert_allclose([thin_s.blue_limit, thin_s.red_limit], [500., 700.])
+        assert len(thin_s.wave_list) < len(s.wave_list)
+        assert len(thin_s.wave_list) < len(combined_in_band)
+
+        # preserve_range=False can trim the support further when the bandpass tails are
+        # insignificant for the requested error tolerance.
+        thin_s2 = s.thin(rel_err=1.e-3, preserve_range=False, bandpass=bp)
+        print('preserve_range=False: ',thin_s2.wave_list)
+        thin_flux2 = thin_s2.calculateFlux(bp.truncate(thin_s2.blue_limit, thin_s2.red_limit))
+        thin_err2 = abs((thin_flux2 - flux) / flux)
+        assert thin_err2 < 1.e-3
+        assert thin_s2.blue_limit > bp.blue_limit
+        assert thin_s2.red_limit < bp.red_limit
+
+        # preserve_range=True, trim_zeros=False should keep the original SED support even though
+        # the bandpass only matters over a narrower interval.
+        thin_s3 = s.thin(rel_err=1.e-3, preserve_range=True, trim_zeros=False, bandpass=bp)
+        print('trim_zeros=False: ',thin_s3.wave_list)
+        thin_flux3 = thin_s3.calculateFlux(bp)
+        thin_err3 = abs((thin_flux3 - flux) / flux)
+        assert thin_err3 < 1.e-3
+        np.testing.assert_allclose([thin_s3.blue_limit, thin_s3.red_limit],
+                                   [s.blue_limit, s.red_limit])
+        # This is not required in general, but is true here.
+        assert len(thin_s3.wave_list) > len(thin_s.wave_list)
+
+    # Non-overlapping wavelength ranges raises an exception.
+    bp = galsim.Bandpass(galsim.LookupTable([950, 980, 1000], [0.1, 1.0, 0.1]), 'nm')
+    assert_raises(galsim.GalSimError, s.thin, bandpass=bp)
+
+    # If the SED and bandpass already share the same support, then preserve_range=True and
+    # trim_zeros=False should not need to add any extra front/back points.
+    s = galsim.SED(galsim.LookupTable(bp_waves_nm, sed_vals[2:-3], interpolant='linear'),
+                   'nm', 'fphotons')
+    bp = galsim.Bandpass(galsim.LookupTable(bp_waves_nm, bp_vals, interpolant='linear'), 'nm')
+    thin_s = s.thin(rel_err=1.e-3, preserve_range=True, trim_zeros=False, bandpass=bp)
+    np.testing.assert_allclose([thin_s.blue_limit, thin_s.red_limit],
+                               [s.blue_limit, s.red_limit])
+
+    # Use a realistic tabulated SED and bandpass with mixed wavelength units plus redshift.
+    s = galsim.SED('CWW_E_ext.sed', wave_type='ang', flux_type='flambda',
+                   fast=False, interpolant='linear').atRedshift(0.2)
+    bp = galsim.Bandpass('LSST_r.dat', 'nm')
+    flux = s.calculateFlux(bp)
+    wave_list, _, _ = galsim.utilities.combine_wave_list(s, bp)
+    combined_in_band = wave_list[np.logical_and(wave_list >= bp.blue_limit,
+                                                wave_list <= bp.red_limit)]
+    print('CWW original: ',len(combined_in_band), s.wave_list)
+    print('in band: ',len(combined_in_band), combined_in_band)
+
+    thin_s = s.thin(rel_err=1.e-3, preserve_range=True, bandpass=bp)
+    np.testing.assert_allclose([thin_s.blue_limit, thin_s.red_limit],
+                               [bp.blue_limit, bp.red_limit])
+    assert len(thin_s.wave_list) < len(combined_in_band)
+    assert len(thin_s.wave_list) < len(s.wave_list)
+    thin_flux = thin_s.calculateFlux(bp)
+    thin_err = abs((thin_flux - flux) / flux)
+    print('thinned', len(thin_s.wave_list), thin_s.wave_list)
+    print('thin_err = ',thin_err)
+    assert thin_err < 1.e-3
+
+    # trim_zeros=False should retain the original redshifted SED support.
+    thin_s2 = s.thin(rel_err=1.e-3, preserve_range=True, trim_zeros=False, bandpass=bp)
+    print('trim_zeros=False', len(thin_s2.wave_list), thin_s2.wave_list)
+    np.testing.assert_allclose([thin_s2.blue_limit, thin_s2.red_limit],
+                               [s.blue_limit, s.red_limit])
+    thin_flux2 = thin_s2.calculateFlux(bp)
+    thin_err2 = abs((thin_flux2 - flux) / flux)
+    print('thin_err2 = ',thin_err2)
+    assert thin_err2 < 1.e-3
+
+    # Compare the thinned error when not using bp in the thin call.
+    thin_s3 = s.thin(rel_err=1.e-3, preserve_range=True, trim_zeros=False)
+    print('trim_zeros=False, no bandpass', len(thin_s3.wave_list), thin_s3.wave_list)
+    np.testing.assert_allclose([thin_s3.blue_limit, thin_s3.red_limit],
+                               [s.blue_limit, s.red_limit])
+    thin_flux3 = thin_s3.calculateFlux(bp)
+    thin_err3 = abs((thin_flux3 - flux) / flux)
+    print('thin_err3 = ',thin_err3)
+    # There are no hard guarantees about this result, but it should normally also be < rel_err.
+    assert thin_err3 < 1.e-3
+    # It would also be expected in practice that the bandpass-aware thinning does a better
+    # job, but there are no guarantees about that either.  It turns out to be true here.
+    assert thin_err2 < thin_err3
+
 @timer
 def test_broadcast():
     """ Check that constand SED broadcasts over waves.