Include intrinsic zernikes in danish output

Fixed adding intrinsics to Danish.

Streamlining zk4up and caches in instrument.

python/lsst/ts/wep/estimation/danish.py

@@ -99,12 +99,14 @@ def history(self) -> dict:
 
         The history is a dictionary that contains intermediate products
         from the Zernike fitting. The dict contains entries for "intra"
-        and/or "extra". Each of these is a dictionary that contains the
+        and/or "extra", plus the final zernike estimate under "zk".
+        The "intra" and "extra" entries are dictionaries that contain the
         following entries
             - "image" - the image that is being fit
             - "variance" - the background variance that was used for fitting
             - "lstsqResult" - dictionary of results returned by least_squares
-            - "zk" - the Zernike coefficients fit to the donut
+            - "zkResid" - the Zernike coefficients fit to the donut
+            - "zk" - zkResid + the intrinsic Zernikes
             - "model" - the final forward modeled donut image
         """
         return super().history
@@ -114,7 +116,7 @@ def _estimateSingleZk(
         image: Image,
         instrument: Instrument,
         factory: danish.DonutFactory,
-        x0: list,
+        jmax: int,
     ) -> Tuple[np.ndarray, dict]:
         """Estimate Zernikes (in meters) for a single donut stamp.
 
@@ -126,8 +128,8 @@ def _estimateSingleZk(
             The ts_wep Instrument
         factory : danish.DonutFactory
             The Danish donut factory
-        x0 : list
-            The initial guess for the model parameters
+        jmax : int
+            The maximum Zernike Noll index returned by Danish
 
         Returns
         -------
@@ -159,7 +161,7 @@ def _estimateSingleZk(
         model = danish.SingleDonutModel(
             factory,
             z_ref=zkRef,
-            z_terms=np.arange(4, len(x0) + 1),
+            z_terms=np.arange(4, jmax + 1),
             thx=np.deg2rad(image.fieldAngle[0]),
             thy=np.deg2rad(image.fieldAngle[1]),
             npix=img.shape[0],
@@ -169,6 +171,9 @@ def _estimateSingleZk(
         var = img[:4].std() ** 2
         img = img - img[:4].mean()  # Make a copy
 
+        # Create the initial guess for the model parameters
+        x0 = [0.0, 0.0, 1.0] + [0.0] * (jmax - 3)
+
         # Use scipy to optimize the parameters
         result = least_squares(
             model.chi,
@@ -179,7 +184,16 @@ def _estimateSingleZk(
         )
 
         # Unpack the parameters
-        dx, dy, fwhm, *zk = result.x
+        dx, dy, fwhm, *zkResid = result.x
+        zkResid = np.array(zkResid)
+
+        # Add the intrinsic telescope zernikes back into the result
+        zkIntrinsic = instrument.getIntrinsicZernikes(
+            *image.fieldAngle,
+            image.bandLabel,
+            jmax=jmax,
+        )
+        zk = zkResid + zkIntrinsic
 
         if self.saveHistory:
             # Save the image and variance
@@ -192,10 +206,11 @@ def _estimateSingleZk(
             hist["lstsqResult"] = dict(result)
 
             # Save the Zernike coefficients in meters
+            hist["zkResid"] = zkResid
             hist["zk"] = zk
 
             # Also add the forward modeled image
-            hist["model"] = model.model(dx, dy, fwhm, zk)
+            hist["model"] = model.model(dx, dy, fwhm, zkResid)
 
         else:
             hist = {}
@@ -253,9 +268,6 @@ def estimateZk(
             pixel_scale=instrument.pixelSize,
         )
 
-        # Create the initial guess for the model parameters
-        x0 = [0.0, 0.0, 1.0] + [0.0] * (jmax - 3)
-
         # Create an empty history
         hist = {}
 
@@ -264,7 +276,7 @@ def estimateZk(
             I1,
             instrument,
             factory,
-            x0,
+            jmax,
         )
 
         if I2 is not None:
@@ -273,15 +285,16 @@ def estimateZk(
                 I2,
                 instrument,
                 factory,
-                x0,
+                jmax,
             )
 
             # Average the Zernikes
             zk = np.mean([zk1, zk2], axis=0)
         else:
-            zk = np.array(zk1)
+            zk = zk1
 
         if self.saveHistory:
+            hist["zk"] = zk
             self._history = hist
 
         return zk

python/lsst/ts/wep/instrument.py

@@ -515,7 +515,6 @@ def _getIntrinsicZernikesCached(
         yAngle: float,
         band: Union[BandLabel, str],
         jmax: int,
-        return4Up: bool,
     ) -> np.ndarray:
         """Cached interior function for the getIntrinsicZernikes method.
 
@@ -533,8 +532,6 @@ def _getIntrinsicZernikesCached(
             "{band}".
         jmax : int, optional
             The maximum Noll index of the intrinsic Zernikes.
-        return4Up : bool
-            Whether to only return the coefficients for Noll indices >= 4.
 
         Returns
         -------
@@ -613,7 +610,6 @@ def _getOffAxisCoeffCached(
         defocalType: DefocalType,
         band: Union[BandLabel, str],
         jmax: int,
-        return4Up: bool,
     ) -> np.ndarray:
         """Cached interior function for the getOffAxisCoeff method.
 
@@ -634,8 +630,6 @@ def _getOffAxisCoeffCached(
             contains "{band}".
         jmax : int
             The maximum Noll index of the off-axis model Zernikes.
-        return4Up : bool
-            Whether to only return the coefficients for Noll indices >= 4.
 
         Returns
         -------