Make the RomanOptics chromatic=True option more effective and more efficient.

devel/roman/diagnose_bilinear_vs_direct.py

@@ -96,20 +96,15 @@ def run_scan(model, scas, grid, stamp_size, scale, margin, aberrations):
 
                 prof = galsim.roman.getPSF(
                     int(sca),
-                    model.filter,
+                    model.filter_name,
                     SCA_pos=star.image_pos,
                     pupil_bin=piff.roman_psf.pupil_bin,
                     wcs=star.image.wcs,
                     extra_aberrations=model._make_extra_aberrations(aberrations),
                     wavelength=None if model.chromatic else model.bandpass.effective_wavelength,
                 )
                 direct_image = star.image.copy()
-                prof.drawImage(
-                    direct_image,
-                    method=model._method,
-                    center=star.image_pos,
-                    bandpass=model.bandpass if model.chromatic else None,
-                )
+                model._draw_profile_to_image(prof, direct_image, center=star.image_pos, star=star)
 
                 max_abs, rms, frac_peak, frac_l2 = point_metrics(
                     model_star.image.array, direct_image.array
@@ -301,12 +296,12 @@ def main():
 
     piff.roman_psf.pupil_bin = int(args.pupil_bin)
     model = piff.Roman(
-        filter="H158",
         chromatic=False,
         max_zernike=args.max_zernike,
         aberration_interp="constant",
         aberration_prior_sigma=1.0e6,
     )
+    model.set_bandpass(galsim.roman.getBandpasses()["H158"])
 
     scas = parse_scas(args.scas)
     aberrations = parse_aberrations(args.aberrations, model.param_len)

devel/roman/piff.yaml

@@ -1,3 +1,5 @@
+# Note: this run takes ~ 80 minutes to run.
+
 input:
 
     # This file is quite large, so not saved in the repo.
@@ -9,6 +11,10 @@ input:
 
     cat_file_name: 'stars_13906_test_v251105.parquet'
 
+    bandpass:
+        type: RomanBandpass
+        name: H158
+
     ra_col: ra
     dec_col: dec
     properties:
@@ -26,16 +32,15 @@ select:
     hsm_size_reject: 5    # number of inter-quartile ranges from median to reject size estimate.
     max_pixel_cut: 50000  # reject stars with a pixel value higher than this. (sort of)
     min_sep: 1.0          # reject stars with a neighbor closer than 1 arcsec away.
-    #nstars: 100
 
 output:
     dir: 'output'
-    file_name: 'ffov_13906_11.piff'
+    file_name: 'ffov_13906_17.piff'
 
     stats:
         - 
           type: StarImages
-          file_name: 'stars_11.png'
+          file_name: 'stars_17.png'
           nplot: 0  # all stars
 
 psf:
@@ -46,12 +51,11 @@ psf:
 
       -
         type: RomanOptics
-        filter: H158
         chromatic: False
         max_zernike: 12    # Fit Zernike coefficients 4-12 inclusive.
         aberration_prior_sigma: 3.e-3  # Gaussian prior on all Zernike coefficients in absolute units
         aberration_interp: constant    # linear is possible but doesn't seem to work better.
-        nominal_interp: bilinear       # five_points is possible but doesn't work better on
+        nominal_interp: bilinear       # five_point is possible but doesn't work better on
                                        # this sim. May be different on real data.
         nproc: 6
 

devel/roman/piff_chromatic.yaml

@@ -0,0 +1,88 @@
+# Note: this run takes ~ 380 minutes to run.
+
+input:
+
+    # This file is quite large, so not saved in the repo.
+    # It and the cat_file are available from here:
+    #   https://drive.google.com/drive/folders/1akvHjdKSTMppPTcNfXjEy9SVq_hFBDBc
+    image_file_name: 'ffov_13906_test_v251020.fits'
+    image_hdu: [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18]  # hdu = SCA
+    badpix_zeros: True
+
+    cat_file_name: 'stars_13906_test_v251105.parquet'
+
+    # For now just use a single SED for all stars.  This is the SED for a K type star.
+    sed_file_name: 'seds/XSL_DR3_release/xsl_spectrum_X0066_merged_scl.fits'
+    sed_wave_type: 'nm'
+    sed_flux_type: 'erg cm**(-2) s**(-1) angstrom**(-1)'
+    sed_max_samples: 2  # Just keep a linear approximation to the sed across the bandpass.
+    bandpass:
+        type: RomanBandpass
+        name: H158
+
+    ra_col: ra
+    dec_col: dec
+    properties:
+        sca: '$@image_num + 1'
+    noise: BKGNDVAR
+    gain: EQVGAIN
+    sky: median
+
+    stamp_size: 25
+
+select:
+
+    min_snr: 50           # reject very faint stars
+    max_snr_weight: 500   # don't over-weight very bright stars
+    hsm_size_reject: 5    # number of inter-quartile ranges from median to reject size estimate.
+    max_pixel_cut: 50000  # reject stars with a pixel value higher than this. (sort of)
+    min_sep: 1.0          # reject stars with a neighbor closer than 1 arcsec away.
+
+output:
+    dir: 'output'
+    file_name: 'ffov_13906_18.piff'
+
+    stats:
+        - 
+          type: StarImages
+          file_name: 'stars_18.png'
+          nplot: 0  # all stars
+
+psf:
+
+    type: Sum
+
+    components:
+
+      -
+        type: RomanOptics
+        chromatic: True
+        max_zernike: 12    # Fit Zernike coefficients 4-12 inclusive.
+        aberration_prior_sigma: 3.e-3  # Gaussian prior on all Zernike coefficients.
+        aberration_interp: constant    # linear is possible but doesn't seem to work better.
+        nominal_interp: bilinear       # five_point is possible but doesn't work better on
+                                       # this sim. May be different on real data.
+        nproc: 6
+
+      - 
+        type: Simple
+        model:
+            type: PixelGrid
+            scale: 0.11
+            size: 5
+        interp:
+            type: BasisPolynomial
+            order: 2
+
+    outliers:
+      -
+        type: Chisq
+        nsigma: 10
+        max_remove: 0.03
+      -
+        type: Centroid
+        max_offset: 0.2 # arcsec
+
+
+
+verbose: 2

piff/config.py

@@ -159,11 +159,12 @@ def process(config, logger=None):
             config['select'][key] = config['input'].pop(key)
 
     # read in the input images
-    objects, wcs, pointing = Input.process(config['input'], logger=logger)
+    objects, wcs, pointing, bandpass = Input.process(config['input'], logger=logger)
     stars = Select.process(config.get('select',{}), objects, logger=logger)
 
     psf = PSF.process(config['psf'], logger=logger)
-    psf.fit(stars, wcs, pointing, logger=logger)
+    psf.set_context(wcs, pointing, bandpass)
+    psf.fit(stars, logger=logger)
 
     # Attach these for reference
     psf.initial_objects = objects

piff/convolvepsf.py

@@ -106,6 +106,13 @@ def interp_property_names(self):
             names.update(c.interp_property_names)
         return names
 
+    def set_context(self, wcs, pointing, bandpass):
+        super().set_context(wcs, pointing, bandpass)
+        if isinstance(self.components, list):
+            # If it's a list, then building components has been completed.
+            for comp in self.components:
+                comp.set_context(wcs, pointing, bandpass)
+
     @classmethod
     def parseKwargs(cls, config_psf, logger):
         """Parse the psf field of a configuration dict and return the kwargs to use for