From 7a91bfa05c57fc8a7c94b506da7658f10b21780f Mon Sep 17 00:00:00 2001 From: Clare Saunders Date: Wed, 4 Sep 2024 18:30:15 -0700 Subject: [PATCH] Add option to build camera in gbdes task --- python/lsst/drp/tasks/build_camera.py | 383 ++++++++++--------- python/lsst/drp/tasks/gbdesAstrometricFit.py | 102 ++++- tests/test_build_camera.py | 86 +++-- 3 files changed, 348 insertions(+), 223 deletions(-) diff --git a/python/lsst/drp/tasks/build_camera.py b/python/lsst/drp/tasks/build_camera.py index 3edb2ce1..c2122a66 100644 --- a/python/lsst/drp/tasks/build_camera.py +++ b/python/lsst/drp/tasks/build_camera.py @@ -20,38 +20,25 @@ # see . # import re -import time from collections import defaultdict import astshim as ast -import lsst.pex.config import numpy as np import scipy.optimize -from lsst.afw.cameraGeom import ( - ACTUAL_PIXELS, - FIELD_ANGLE, - FOCAL_PLANE, - PIXELS, - TAN_PIXELS, - Camera, - CameraSys, - DetectorConfig, -) -from lsst.afw.geom import TransformPoint2ToPoint2, transformConfig + +import lsst.pex.config +from lsst.afw.cameraGeom import FIELD_ANGLE, FOCAL_PLANE, Camera +from lsst.afw.geom import TransformPoint2ToPoint2 from lsst.geom import degrees from lsst.pipe.base import Task -from .gbdesAstrometricFit import _degreeFromNCoeffs, _nCoeffsFromDegree - __all__ = [ - "BuildCameraConfig", - "BuildCameraTask", + "BuildCameraFromAstrometryConfig", + "BuildCameraFromAstrometryTask", ] -cameraSysList = [FIELD_ANGLE, FOCAL_PLANE, PIXELS, TAN_PIXELS, ACTUAL_PIXELS] -cameraSysMap = dict((sys.getSysName(), sys) for sys in cameraSysList) - -Nfeval = 0 +# Set up global variable to use in minimization callback. +nFunctionEval = 0 def _z_function(params, x, y, order=4): @@ -108,11 +95,11 @@ def _z_function_dx(params, x, y, order=4): for j in range(i + 1): coeffs[j, i - j] = params[c] c += 1 - deriv_coeffs = np.zeros(coeffs.shape) - deriv_coeffs[:, :-1] = coeffs[:, 1:] - deriv_coeffs[:, :-1] *= np.arange(1, order + 1) + derivCoeffs = np.zeros(coeffs.shape) + derivCoeffs[:, :-1] = coeffs[:, 1:] + derivCoeffs[:, :-1] *= np.arange(1, order + 1) - z = np.polynomial.polynomial.polyval2d(x, y, deriv_coeffs.T) + z = np.polynomial.polynomial.polyval2d(x, y, derivCoeffs.T) return z @@ -141,15 +128,15 @@ def _z_function_dy(params, x, y, order=4): for j in range(i + 1): coeffs[j, i - j] = params[c] c += 1 - deriv_coeffs = np.zeros(coeffs.shape) - deriv_coeffs[:-1] = coeffs[1:] - deriv_coeffs[:-1] *= np.arange(1, order + 1).reshape(-1, 1) + derivCoeffs = np.zeros(coeffs.shape) + derivCoeffs[:-1] = coeffs[1:] + derivCoeffs[:-1] *= np.arange(1, order + 1).reshape(-1, 1) - z = np.polynomial.polynomial.polyval2d(x, y, deriv_coeffs.T) + z = np.polynomial.polynomial.polyval2d(x, y, derivCoeffs.T) return z -class BuildCameraConfig(lsst.pex.config.Config): +class BuildCameraFromAstrometryConfig(lsst.pex.config.Config): """Configuration for BuildCameraTask.""" tangentPlaneDegree = lsst.pex.config.Field( @@ -184,7 +171,7 @@ class BuildCameraConfig(lsst.pex.config.Config): plateScale = lsst.pex.config.Field( dtype=float, doc=("Scaling between camera coordinates in mm and angle on the sky in" " arcsec."), - default=1.0, + default=20.005867576692737, ) astInversionTolerance = lsst.pex.config.Field( dtype=float, @@ -201,29 +188,20 @@ class BuildCameraConfig(lsst.pex.config.Config): ) -class BuildCameraTask(Task): +class BuildCameraFromAstrometryTask(Task): """Build an `lsst.afw.cameraGeom.Camera` object out of the `gbdes` polynomials mapping from pixels to the tangent plane. Parameters ---------- - camera : `lsst.afw.cameraGeom.Camera` - Camera object from which to take pupil function, name, and other - properties. - detectorList : `list` [`int`] - List of detector ids. - visitList : `list` [`int`] - List of ids for visits that were used to train the input model. + """ - ConfigClass = BuildCameraConfig - _DefaultName = "buildCamera" + ConfigClass = BuildCameraFromAstrometryConfig + _DefaultName = "buildCameraFromAstrometry" - def __init__(self, camera, detectorList, visitList, **kwargs): + def __init__(self, **kwargs): super().__init__(**kwargs) - self.camera = camera - self.detectorList = detectorList - self.visitList = visitList # The gbdes model normalizes the pixel positions to the range -1 - +1. X = np.arange(-1, 1, 0.1) @@ -232,7 +210,7 @@ def __init__(self, camera, detectorList, visitList, **kwargs): self.x = x.ravel() self.y = y.ravel() - def run(self, mapParams, mapTemplate): + def run(self, mapParams, mapTemplate, detectorList, visitList, inputCamera, rotationAngle): """Convert the model parameters into a Camera object. Parameters @@ -243,6 +221,15 @@ def run(self, mapParams, mapTemplate): mapTemplate : `dict` Dictionary describing the format of the astrometric model, following the convention in `gbdes`. + detectorList : `list` [`int`] + List of detector ids. + visitList : `list` [`int`] + List of ids for visits that were used to train the input model. + camera : `lsst.afw.cameraGeom.Camera` + Camera object from which to take pupil function, name, and other + properties. + rotationAngle : `float` + Value in radians of the average rotation angle of the input visits. Returns ------- @@ -251,27 +238,29 @@ def run(self, mapParams, mapTemplate): """ # Normalize the model. - newParams, newIntX, newIntY = self._normalizeModel(mapParams, mapTemplate) + newParams, newTPx, newTPy = self._normalize_model(mapParams, mapTemplate, detectorList, visitList) if self.config.modelSplitting == "basic": # Put all of the camera distortion into the pixels->focal plane # part of the distortion model, with the focal plane->tangent plane # part only used for scaling between the focal plane and sky. - pixToFocalPlane, focalPlaneToTangentPlane = self._basicModel(newParams) + pixToFocalPlane, focalPlaneToTangentPlane = self._basic_model(newParams, detectorList) else: # Fit two polynomials, such that the first describes the pixel-> # focal plane part of the model, and the second describes the focal # plane->tangent plane part of the model, with the goal of # generating a more physically-motivated distortion model. - pixToFocalPlane, focalPlaneToTangentPlane = self._splitModel(newIntX, newIntY) + pixToFocalPlane, focalPlaneToTangentPlane = self._split_model(newTPx, newTPy, detectorList) # Turn the mappings into a Camera object. - camera = self._translateToAfw(pixToFocalPlane, focalPlaneToTangentPlane) + camera = self._translate_to_afw( + pixToFocalPlane, focalPlaneToTangentPlane, detectorList, inputCamera, rotationAngle + ) return camera - def _normalizeModel(self, mapParams, mapTemplate): + def _normalize_model(self, mapParams, mapTemplate, detectorList, visitList): """Normalize the camera mappings, such that they correspond to the average visit. @@ -291,13 +280,17 @@ def _normalizeModel(self, mapParams, mapTemplate): mapTemplate : `dict` Dictionary describing the format of the astrometric model, following the convention in `gbdes`. + detectorList : `list` [`int`] + List of detector ids. + visitList : `list` [`int`] + List of ids for visits that were used to train the input model. Returns ------- newDeviceArray : `np.ndarray` Array of NxM, where N is the number of detectors, and M is the number of coefficients for each per-detector mapping. - newIntX, newIntY : `np.ndarray` + newTPx, newTPy : `np.ndarray` Projection of `self.x` and `self.y` onto the tangent plane, given the normalized mapping. """ @@ -306,7 +299,7 @@ def _normalizeModel(self, mapParams, mapTemplate): deviceParams = [] visitParams = [] for element in mapTemplate["BAND/DEVICE"]["Elements"]: - for detector in self.detectorList: + for detector in detectorList: detectorTemplate = element.replace("DEVICE", str(detector)) detectorTemplate = detectorTemplate.replace("BAND", ".+") for k, params in mapParams.items(): @@ -315,7 +308,7 @@ def _normalizeModel(self, mapParams, mapTemplate): deviceArray = np.vstack(deviceParams) for element in mapTemplate["EXPOSURE"]["Elements"]: - for visit in self.visitList: + for visit in visitList: visitTemplate = element.replace("EXPOSURE", str(visit)) for k, params in mapParams.items(): if re.fullmatch(visitTemplate, k): @@ -326,6 +319,8 @@ def _normalizeModel(self, mapParams, mapTemplate): expoMean = expoArray.mean(axis=0) + # Shift the per-device part of the model to correspond with the mean + # per-visit behavior. newDeviceArray = np.zeros(deviceArray.shape) nCoeffsDev = deviceArray.shape[1] // 2 newDeviceArray[:, :nCoeffsDev] = ( @@ -337,24 +332,24 @@ def _normalizeModel(self, mapParams, mapTemplate): newDeviceArray[:, 0] += expoMean[0] newDeviceArray[:, nCoeffsDev] += expoMean[3] - # Then get the interim positions from the new device model: - newIntX = [] - newIntY = [] + # Then get the tangent plane positions from the new device model: + newTPx = [] + newTPy = [] for deviceMap in newDeviceArray: nCoeffsDev = len(deviceMap) // 2 - deviceDegree = _degreeFromNCoeffs(nCoeffsDev) + deviceDegree = int(-1.5 + 0.5 * (1 + 8 * nCoeffsDev) ** 0.5) intX = _z_function(deviceMap[:nCoeffsDev], self.x, self.y, order=deviceDegree) intY = _z_function(deviceMap[nCoeffsDev:], self.x, self.y, order=deviceDegree) - newIntX.append(intX) - newIntY.append(intY) + newTPx.append(intX) + newTPy.append(intY) - newIntX = np.array(newIntX).ravel() - newIntY = np.array(newIntY).ravel() + newTPx = np.array(newTPx).ravel() + newTPy = np.array(newTPy).ravel() - return newDeviceArray, newIntX, newIntY + return newDeviceArray, newTPx, newTPy - def _basicModel(self, modelParameters): + def _basic_model(self, modelParameters, detectorList): """This will just convert the pix->fp parameters into the right format, and return an identity mapping for the fp->tp part. @@ -363,6 +358,8 @@ def _basicModel(self, modelParameters): modelParameters : `np.ndarray` Array of NxM, where N is the number of detectors, and M is the number of coefficients for each per-detector mapping. + detectorList : `list` [`int`] + List of detector ids. Returns ------- @@ -374,17 +371,17 @@ def _basicModel(self, modelParameters): Dictionary giving the focal plane to tangent plane mapping. """ - nCoeffsFP = modelParameters.shape[1] // 2 + nCoeffsFp = modelParameters.shape[1] // 2 pixelsToFocalPlane = defaultdict(dict) - for d, det in enumerate(self.detectorList): - pixelsToFocalPlane[det]["x"] = modelParameters[d][:nCoeffsFP] - pixelsToFocalPlane[det]["y"] = modelParameters[d][nCoeffsFP:] + for d, det in enumerate(detectorList): + pixelsToFocalPlane[det]["x"] = modelParameters[d][:nCoeffsFp] + pixelsToFocalPlane[det]["y"] = modelParameters[d][nCoeffsFp:] focalPlaneToTangentPlane = {"x": [0, 1, 0], "y": [0, 0, 1]} return pixelsToFocalPlane, focalPlaneToTangentPlane - def _splitModel(self, targetX, targetY, pixToFPGuess=None, fpToTpGuess=None): + def _split_model(self, targetX, targetY, detectorList, pixToFpGuess=None, fpToTpGuess=None): """Fit a two-step model, with one polynomial per detector fitting the pixels to focal plane part, followed by a polynomial fitting the focal plane to tangent plane part. @@ -396,7 +393,9 @@ def _splitModel(self, targetX, targetY, pixToFPGuess=None, fpToTpGuess=None): ---------- targetX, targetY : `np.ndarray` Target x and y values in the tangent plane. - pixToFPGuess : `dict` [`int`, [`string`, `float`]] + detectorList : `list` [`int`] + List of detector ids. + pixToFpGuess : `dict` [`int`, [`string`, `float`]] Initial guess for the pixels to focal plane mapping to use in the model fit, in the form of a dictionary giving the per-detector pixel to focal plane mapping, with keys for each detector giving @@ -419,17 +418,17 @@ def _splitModel(self, targetX, targetY, pixToFPGuess=None, fpToTpGuess=None): tpDegree = self.config.tangentPlaneDegree fpDegree = self.config.focalPlaneDegree - nCoeffsFP = _nCoeffsFromDegree(fpDegree) - nCoeffsFP_tot = len(self.detectorList) * nCoeffsFP * 2 + nCoeffsFP = int((fpDegree + 2) * (fpDegree + 1) / 2) + nCoeffsFPTot = len(detectorList) * nCoeffsFP * 2 - nCoeffsTP = _nCoeffsFromDegree(tpDegree) + nCoeffsTP = int((tpDegree + 2) * (tpDegree + 1) / 2) # The constant and linear terms will be fixed to remove degeneracy with # the pix->fp part of the model nCoeffsFixed = 3 nCoeffsFree = nCoeffsTP - nCoeffsFixed - fx_params = np.zeros(nCoeffsTP * 2) - fx_params[1] = 1 - fx_params[nCoeffsTP + 2] = 1 + fixedParams = np.zeros(nCoeffsTP * 2) + fixedParams[1] = 1 + fixedParams[nCoeffsTP + 2] = 1 nX = len(self.x) # We need an array of the form [1, x, y, x^2, xy, y^2, ...] to use when @@ -445,26 +444,26 @@ def _splitModel(self, targetX, targetY, pixToFPGuess=None, fpToTpGuess=None): def two_part_function(params): # The function giving the split model. - intX_tot = [] - intY_tot = [] - for i in range(len(self.detectorList)): - intX = _z_function( + fpXAll = [] + fpYAll = [] + for i in range(len(detectorList)): + fpX = _z_function( params[2 * nCoeffsFP * i : (2 * i + 1) * nCoeffsFP], self.x, self.y, order=fpDegree ) - intY = _z_function( + fpY = _z_function( params[(2 * i + 1) * nCoeffsFP : 2 * nCoeffsFP * (1 + i)], self.x, self.y, order=fpDegree ) - intX_tot.append(intX) - intY_tot.append(intY) - intX_tot = np.array(intX_tot).ravel() - intY_tot = np.array(intY_tot).ravel() + fpXAll.append(fpX) + fpYAll.append(fpY) + fpXAll = np.array(fpXAll).ravel() + fpYAll = np.array(fpYAll).ravel() - tpParams = fx_params.copy() - tpParams[nCoeffsFixed:nCoeffsTP] = params[nCoeffsFP_tot : nCoeffsFP_tot + nCoeffsFree] - tpParams[nCoeffsFixed + nCoeffsTP :] = params[nCoeffsFP_tot + nCoeffsFree :] + tpParams = fixedParams.copy() + tpParams[nCoeffsFixed:nCoeffsTP] = params[nCoeffsFPTot : nCoeffsFPTot + nCoeffsFree] + tpParams[nCoeffsFixed + nCoeffsTP :] = params[nCoeffsFPTot + nCoeffsFree :] - tpX = _z_function(tpParams[:nCoeffsTP], intX_tot, intY_tot, order=tpDegree) - tpY = _z_function(tpParams[nCoeffsTP:], intX_tot, intY_tot, order=tpDegree) + tpX = _z_function(tpParams[:nCoeffsTP], fpXAll, fpYAll, order=tpDegree) + tpY = _z_function(tpParams[nCoeffsTP:], fpXAll, fpYAll, order=tpDegree) return tpX, tpY def min_function(params): @@ -480,141 +479,141 @@ def jac(params): dX = -2 * (targetX - tpX) dY = -2 * (targetY - tpY) jacobian = np.zeros(len(params)) - fp_params = params[:nCoeffsFP_tot] - tp_params = fx_params.copy() - tp_params[nCoeffsFixed:nCoeffsTP] = params[nCoeffsFP_tot : nCoeffsFP_tot + nCoeffsFree] - tp_params[nCoeffsFixed + nCoeffsTP :] = params[nCoeffsFP_tot + nCoeffsFree :] - intX_tot = [] - intY_tot = [] + fpParams = params[:nCoeffsFPTot] + tpParams = fixedParams.copy() + tpParams[nCoeffsFixed:nCoeffsTP] = params[nCoeffsFPTot : nCoeffsFPTot + nCoeffsFree] + tpParams[nCoeffsFixed + nCoeffsTP :] = params[nCoeffsFPTot + nCoeffsFree :] + fpXAll = [] + fpYAll = [] # Fill in the derivatives for the pix->fp terms. - for i in range(len(self.detectorList)): + for i in range(len(detectorList)): dX_i = dX[nX * i : nX * (i + 1)] dY_i = dY[nX * i : nX * (i + 1)] - intX = _z_function( - fp_params[2 * nCoeffsFP * i : (2 * i + 1) * nCoeffsFP], self.x, self.y, order=fpDegree + fpX = _z_function( + fpParams[2 * nCoeffsFP * i : (2 * i + 1) * nCoeffsFP], self.x, self.y, order=fpDegree ) - intY = _z_function( - fp_params[(2 * i + 1) * nCoeffsFP : 2 * nCoeffsFP * (1 + i)], + fpY = _z_function( + fpParams[(2 * i + 1) * nCoeffsFP : 2 * nCoeffsFP * (1 + i)], self.x, self.y, order=fpDegree, ) - intX_tot.append(intX) - intY_tot.append(intY) + fpXAll.append(fpX) + fpYAll.append(fpY) - dTpX_dfpX = _z_function_dx(tp_params[:nCoeffsTP], intX, intY, order=tpDegree) - dTpX_dfpY = _z_function_dy(tp_params[:nCoeffsTP], intX, intY, order=tpDegree) - dTpY_dfpX = _z_function_dx(tp_params[nCoeffsTP:], intX, intY, order=tpDegree) - dTpY_dfpY = _z_function_dy(tp_params[nCoeffsTP:], intX, intY, order=tpDegree) + dTpX_dFpX = _z_function_dx(tpParams[:nCoeffsTP], fpX, fpY, order=tpDegree) + dTpX_dFpY = _z_function_dy(tpParams[:nCoeffsTP], fpX, fpY, order=tpDegree) + dTpY_dFpX = _z_function_dx(tpParams[nCoeffsTP:], fpX, fpY, order=tpDegree) + dTpY_dFpY = _z_function_dy(tpParams[nCoeffsTP:], fpX, fpY, order=tpDegree) - dTpX_part = np.concatenate([xyArray * dTpX_dfpX, xyArray * dTpX_dfpY]) - dTpY_part = np.concatenate([xyArray * dTpY_dfpX, xyArray * dTpY_dfpY]) + dTpX_part = np.concatenate([xyArray * dTpX_dFpX, xyArray * dTpX_dFpY]) + dTpY_part = np.concatenate([xyArray * dTpY_dFpX, xyArray * dTpY_dFpY]) jacobian_i = (dX_i * dTpX_part + dY_i * dTpY_part).sum(axis=1) jacobian[2 * nCoeffsFP * i : 2 * nCoeffsFP * (i + 1)] = jacobian_i - intX_tot = np.array(intX_tot).ravel() - intY_tot = np.array(intY_tot).ravel() + fpXAll = np.array(fpXAll).ravel() + fpYAll = np.array(fpYAll).ravel() # Fill in the derivatives for the fp->tp terms. for j in range(nCoeffsFree): tParams = np.zeros(nCoeffsTP) tParams[nCoeffsFixed + j] = 1 - tmpZ = _z_function(tParams, intX_tot, intY_tot, order=tpDegree) - jacobian[nCoeffsFP_tot + j] = (dX * tmpZ).sum() - jacobian[nCoeffsFP_tot + nCoeffsFree + j] = (dY * tmpZ).sum() + tmpZ = _z_function(tParams, fpXAll, fpYAll, order=tpDegree) + jacobian[nCoeffsFPTot + j] = (dX * tmpZ).sum() + jacobian[nCoeffsFPTot + nCoeffsFree + j] = (dY * tmpZ).sum() return jacobian def hessian(params): # The Hessian of min_function, which is needed for the minimizer. hessian = np.zeros((len(params), len(params))) - fp_params = params[:nCoeffsFP_tot] - tp_params = fx_params.copy() - tp_params[nCoeffsFixed:nCoeffsTP] = params[nCoeffsFP_tot : nCoeffsFP_tot + nCoeffsFree] - tp_params[nCoeffsFixed + nCoeffsTP :] = params[nCoeffsFP_tot + nCoeffsFree :] - intX_tot = [] - intY_tot = [] + fpParams = params[:nCoeffsFPTot] + tpParams = fixedParams.copy() + tpParams[nCoeffsFixed:nCoeffsTP] = params[nCoeffsFPTot : nCoeffsFPTot + nCoeffsFree] + tpParams[nCoeffsFixed + nCoeffsTP :] = params[nCoeffsFPTot + nCoeffsFree :] + fpXAll = [] + fpYAll = [] # Loop over the detectors to fill in the d(pix->fp)**2 and # d(pix->fp)d(fp->tp) cross terms. - for i in range(len(self.detectorList)): - intX = _z_function( - fp_params[2 * nCoeffsFP * i : (2 * i + 1) * nCoeffsFP], self.x, self.y, order=fpDegree + for i in range(len(detectorList)): + fpX = _z_function( + fpParams[2 * nCoeffsFP * i : (2 * i + 1) * nCoeffsFP], self.x, self.y, order=fpDegree ) - intY = _z_function( - fp_params[(2 * i + 1) * nCoeffsFP : 2 * nCoeffsFP * (1 + i)], + fpY = _z_function( + fpParams[(2 * i + 1) * nCoeffsFP : 2 * nCoeffsFP * (1 + i)], self.x, self.y, order=fpDegree, ) - intX_tot.append(intX) - intY_tot.append(intY) - dTpX_dfpX = _z_function_dx(tp_params[:nCoeffsTP], intX, intY, order=tpDegree) - dTpX_dfpY = _z_function_dy(tp_params[:nCoeffsTP], intX, intY, order=tpDegree) - dTpY_dfpX = _z_function_dx(tp_params[nCoeffsTP:], intX, intY, order=tpDegree) - dTpY_dfpY = _z_function_dy(tp_params[nCoeffsTP:], intX, intY, order=tpDegree) + fpXAll.append(fpX) + fpYAll.append(fpY) + dTpX_dFpX = _z_function_dx(tpParams[:nCoeffsTP], fpX, fpY, order=tpDegree) + dTpX_dFpY = _z_function_dy(tpParams[:nCoeffsTP], fpX, fpY, order=tpDegree) + dTpY_dFpX = _z_function_dx(tpParams[nCoeffsTP:], fpX, fpY, order=tpDegree) + dTpY_dFpY = _z_function_dy(tpParams[nCoeffsTP:], fpX, fpY, order=tpDegree) - dTpX_part = np.concatenate([xyArray * dTpX_dfpX, xyArray * dTpX_dfpY]) - dTpY_part = np.concatenate([xyArray * dTpY_dfpX, xyArray * dTpY_dfpY]) + dTpX_part = np.concatenate([xyArray * dTpX_dFpX, xyArray * dTpX_dFpY]) + dTpY_part = np.concatenate([xyArray * dTpY_dFpX, xyArray * dTpY_dFpY]) - for l in range(2 * nCoeffsFP): + for k in range(2 * nCoeffsFP): for m in range(2 * nCoeffsFP): - hessian[2 * nCoeffsFP * i + l, 2 * nCoeffsFP * i + m] = ( - 2 * (dTpX_part[l] * dTpX_part[m] + dTpY_part[l] * dTpY_part[m]).sum() + hessian[2 * nCoeffsFP * i + k, 2 * nCoeffsFP * i + m] = ( + 2 * (dTpX_part[k] * dTpX_part[m] + dTpY_part[k] * dTpY_part[m]).sum() ) for j in range(nCoeffsFree): - tParams = np.zeros(nCoeffsTP) - tParams[nCoeffsFixed + j] = 1 - tmpZ = _z_function(tParams, intX, intY, order=tpDegree) + tmpParams = np.zeros(nCoeffsTP) + tmpParams[nCoeffsFixed + j] = 1 + tmpZ = _z_function(tmpParams, fpX, fpY, order=tpDegree) - hessX = 2 * (dTpX_part[l] * tmpZ).sum() - hessY = 2 * (dTpY_part[l] * tmpZ).sum() + hessX = 2 * (dTpX_part[k] * tmpZ).sum() + hessY = 2 * (dTpY_part[k] * tmpZ).sum() # dTP_x part - hessian[2 * nCoeffsFP * i + l, nCoeffsFP_tot + j] = hessX - hessian[nCoeffsFP_tot + j, 2 * nCoeffsFP * i + l] = hessX + hessian[2 * nCoeffsFP * i + k, nCoeffsFPTot + j] = hessX + hessian[nCoeffsFPTot + j, 2 * nCoeffsFP * i + k] = hessX # dTP_y part - hessian[2 * nCoeffsFP * i + l, nCoeffsFP_tot + nCoeffsFree + j] = hessY - hessian[nCoeffsFP_tot + nCoeffsFree + j, 2 * nCoeffsFP * i + l] = hessY + hessian[2 * nCoeffsFP * i + k, nCoeffsFPTot + nCoeffsFree + j] = hessY + hessian[nCoeffsFPTot + nCoeffsFree + j, 2 * nCoeffsFP * i + k] = hessY - intX_tot = np.array(intX_tot).ravel() - intY_tot = np.array(intY_tot).ravel() - tmpZ_array = np.zeros((nCoeffsFree, nX * len(self.detectorList))) + fpXAll = np.array(fpXAll).ravel() + fpYAll = np.array(fpYAll).ravel() + tmpZArray = np.zeros((nCoeffsFree, nX * len(detectorList))) # Finally, get the d(fp->tp)**2 terms for j in range(nCoeffsFree): tParams = np.zeros(nCoeffsTP) tParams[nCoeffsFixed + j] = 1 - tmpZ_array[j] = _z_function(tParams, intX_tot, intY_tot, order=tpDegree) + tmpZArray[j] = _z_function(tParams, fpXAll, fpYAll, order=tpDegree) for j in range(nCoeffsFree): for m in range(nCoeffsFree): # X-Y cross terms are zero - hess = 2 * (tmpZ_array[j] * tmpZ_array[m]).sum() - hessian[nCoeffsFP_tot + j, nCoeffsFP_tot + m] = hess - hessian[nCoeffsFP_tot + nCoeffsFree + j, nCoeffsFP_tot + nCoeffsFree + m] = hess + hess = 2 * (tmpZArray[j] * tmpZArray[m]).sum() + hessian[nCoeffsFPTot + j, nCoeffsFPTot + m] = hess + hessian[nCoeffsFPTot + nCoeffsFree + j, nCoeffsFPTot + nCoeffsFree + m] = hess return hessian - global Nfeval - Nfeval = 0 + global nFunctionEval + nFunctionEval = 0 def callbackMF(params): - global Nfeval - self.log.info(f"Iteration {Nfeval}, function value {min_function(params)}") - Nfeval += 1 - - initGuess = np.zeros(nCoeffsFP_tot + 2 * nCoeffsFree) - if pixToFPGuess: - useVar = min(nCoeffsFP, len(pixToFPGuess[0]["x"])) - for i, det in enumerate(self.detectorList): - initGuess[2 * nCoeffsFP * i : 2 * nCoeffsFP * i + useVar] = pixToFPGuess[det]["x"][:useVar] - initGuess[(2 * i + 1) * nCoeffsFP : (2 * i + 1) * nCoeffsFP + useVar] = pixToFPGuess[det][ + global nFunctionEval + self.log.info(f"Iteration {nFunctionEval}, function value {min_function(params)}") + nFunctionEval += 1 + + initGuess = np.zeros(nCoeffsFPTot + 2 * nCoeffsFree) + if pixToFpGuess: + useVar = min(nCoeffsFP, len(pixToFpGuess[0]["x"])) + for i, det in enumerate(detectorList): + initGuess[2 * nCoeffsFP * i : 2 * nCoeffsFP * i + useVar] = pixToFpGuess[det]["x"][:useVar] + initGuess[(2 * i + 1) * nCoeffsFP : (2 * i + 1) * nCoeffsFP + useVar] = pixToFpGuess[det][ "y" ][:useVar] if fpToTpGuess: useVar = min(nCoeffsTP, len(fpToTpGuess["x"])) - initGuess[nCoeffsFP_tot : nCoeffsFP_tot + useVar - nCoeffsFixed] = fpToTpGuess["x"][ + initGuess[nCoeffsFPTot : nCoeffsFPTot + useVar - nCoeffsFixed] = fpToTpGuess["x"][ nCoeffsFixed:useVar ] - initGuess[nCoeffsFP_tot + nCoeffsFree : nCoeffsFP_tot + nCoeffsFree + useVar - nCoeffsFixed] = ( + initGuess[nCoeffsFPTot + nCoeffsFree : nCoeffsFPTot + nCoeffsFree + useVar - nCoeffsFixed] = ( fpToTpGuess["y"][nCoeffsFixed:useVar] ) @@ -631,21 +630,21 @@ def callbackMF(params): # Convert parameters to a dictionary. pixToFP = {} - for i, det in enumerate(self.detectorList): + for i, det in enumerate(detectorList): pixToFP[det] = { "x": res.x[2 * nCoeffsFP * i : (2 * i + 1) * nCoeffsFP], "y": res.x[(2 * i + 1) * nCoeffsFP : 2 * nCoeffsFP * (1 + i)], } - fpToTpAll = fx_params.copy() - fpToTpAll[nCoeffsFixed:nCoeffsTP] = res.x[nCoeffsFP_tot : nCoeffsFP_tot + nCoeffsFree] - fpToTpAll[nCoeffsFixed + nCoeffsTP :] = res.x[nCoeffsFP_tot + nCoeffsFree :] + fpToTpAll = fixedParams.copy() + fpToTpAll[nCoeffsFixed:nCoeffsTP] = res.x[nCoeffsFPTot : nCoeffsFPTot + nCoeffsFree] + fpToTpAll[nCoeffsFixed + nCoeffsTP :] = res.x[nCoeffsFPTot + nCoeffsFree :] fpToTp = {"x": fpToTpAll[:nCoeffsTP], "y": fpToTpAll[nCoeffsTP:]} return pixToFP, fpToTp - def makeAstPolyMapCoeffs(self, order, xCoeffs, yCoeffs): - """Convert polynomial coefficients in gbdes format to AST PolyMap + def make_ast_polymap_coeffs(self, order, xCoeffs, yCoeffs): + """Convert polynomial coefficients in gbdes format to AST PolyMap input format. Paramaters @@ -659,8 +658,8 @@ def makeAstPolyMapCoeffs(self, order, xCoeffs, yCoeffs): where N is the polynomial order. Returns ------- - focalPlaneToPupil: `lsst.afw.geom.TransformPoint2ToPoint2` - Transform from focal plane to field angle coordinates + polyArray : `np.ndarray` + Array formatted for AST PolyMap input. """ N = len(xCoeffs) @@ -676,7 +675,9 @@ def makeAstPolyMapCoeffs(self, order, xCoeffs, yCoeffs): return polyArray - def _translateToAfw(self, pixToFocalPlane, focalPlaneToTangentPlane): + def _translate_to_afw( + self, pixToFocalPlane, focalPlaneToTangentPlane, detectorList, inputCamera, rotationAngle + ): """Convert the model parameters to a Camera object. Parameters @@ -687,6 +688,8 @@ def _translateToAfw(self, pixToFocalPlane, focalPlaneToTangentPlane): polynomials. focalPlaneToTangentPlane : `dict` [`string`, `float`] Dictionary giving the focal plane to tangent plane mapping. + rotationAngle : `float` + Value in radians of the average rotation angle of the input visits. Returns ------- @@ -695,21 +698,25 @@ def _translateToAfw(self, pixToFocalPlane, focalPlaneToTangentPlane): """ if self.config.modelSplitting == "basic": tpDegree = 1 - fpDegree = _degreeFromNCoeffs(len(pixToFocalPlane[self.detectorList[0]]["x"])) + nCoeffsFP = len(pixToFocalPlane[detectorList[0]]["x"]) + fpDegree = int(-1.5 + 0.5 * (1 + 8 * nCoeffsFP) ** 0.5) else: tpDegree = self.config.tangentPlaneDegree fpDegree = self.config.focalPlaneDegree scaleConvert = (1 * degrees).asArcseconds() / self.config.plateScale - cameraBuilder = Camera.Builder(self.camera.getName()) - cameraBuilder.setPupilFactoryName(self.camera.getPupilFactoryName()) + cameraBuilder = Camera.Builder(inputCamera.getName()) + cameraBuilder.setPupilFactoryName(inputCamera.getPupilFactoryName()) # Convert fp->tp to AST format: - forwardCoeffs = self.makeAstPolyMapCoeffs( + forwardCoeffs = self.make_ast_polymap_coeffs( tpDegree, focalPlaneToTangentPlane["x"], focalPlaneToTangentPlane["y"] ) - rotateAndFlipCoeffs = self.makeAstPolyMapCoeffs(1, [0, 0, -1], [0, -1, 0]) + # Reverse rotation from input visits and flip x-axis. + cosRot = np.cos(-rotationAngle) + sinRot = np.sin(-rotationAngle) + rotateAndFlipCoeffs = self.make_ast_polymap_coeffs(1, [0, cosRot, -sinRot], [0, -sinRot, cosRot]) ccdZoom = ast.ZoomMap(2, 1 / scaleConvert) ccdToSky = ast.PolyMap( @@ -731,18 +738,20 @@ def _translateToAfw(self, pixToFocalPlane, focalPlaneToTangentPlane): cameraBuilder.setTransformFromFocalPlaneTo(FIELD_ANGLE, focalPlaneToTPMapping) # Convert pix->fp to AST format: - for detector in self.camera: + for detector in inputCamera: d = detector.getId() if d not in pixToFocalPlane: - # TODO: just grab detector map from the obs_ package. + # This camera will not include detectors that were not used in + # astrometric fit. continue detectorBuilder = cameraBuilder.add(detector.getName(), detector.getId()) - # TODO: Add all the detector details like pixel size, etc. + detectorBuilder.setBBox(detector.getBBox()) + detectorBuilder.setPixelSize(detector.getPixelSize()) for amp in detector.getAmplifiers(): detectorBuilder.append(amp.rebuild()) - normCoeffs = self.makeAstPolyMapCoeffs( + normCoeffs = self.make_ast_polymap_coeffs( 1, [-1.0, 2 / detector.getBBox().getWidth(), 0], [-1.0, 0, 2 / detector.getBBox().getHeight()] ) normMap = ast.PolyMap( @@ -751,7 +760,7 @@ def _translateToAfw(self, pixToFocalPlane, focalPlaneToTangentPlane): "IterInverse=1, TolInverse=%s, NIterInverse=%s" % (self.config.astInversionTolerance / 2.0, self.config.astInversionMaxIter), ) - forwardDetCoeffs = self.makeAstPolyMapCoeffs( + forwardDetCoeffs = self.make_ast_polymap_coeffs( fpDegree, pixToFocalPlane[d]["x"], pixToFocalPlane[d]["y"] ) ccdToFP = ast.PolyMap( @@ -765,6 +774,6 @@ def _translateToAfw(self, pixToFocalPlane, focalPlaneToTangentPlane): ccdToFPMapping = TransformPoint2ToPoint2(fullDetMap) detectorBuilder.setTransformFromPixelsTo(FOCAL_PLANE, ccdToFPMapping) - output_camera = cameraBuilder.finish() + outputCamera = cameraBuilder.finish() - return output_camera + return outputCamera diff --git a/python/lsst/drp/tasks/gbdesAstrometricFit.py b/python/lsst/drp/tasks/gbdesAstrometricFit.py index 2a0fad78..7c380716 100644 --- a/python/lsst/drp/tasks/gbdesAstrometricFit.py +++ b/python/lsst/drp/tasks/gbdesAstrometricFit.py @@ -44,6 +44,8 @@ ) from lsst.meas.algorithms.sourceSelector import sourceSelectorRegistry +from .build_camera import BuildCameraFromAstrometryTask + __all__ = [ "calculate_apparent_motion", "GbdesAstrometricFitConnections", @@ -349,6 +351,13 @@ class GbdesAstrometricFitConnections( storageClass="ArrowNumpyDict", dimensions=("instrument", "physical_filter"), ) + inputCamera = pipeBase.connectionTypes.PrerequisiteInput( + doc="Input camera to use when constructing camera from astrometric model.", + name="camera", + storageClass="Camera", + dimensions=("instrument",), + isCalibration=True, + ) outputWcs = pipeBase.connectionTypes.Output( doc=( "Per-tract, per-visit world coordinate systems derived from the fitted model." @@ -388,7 +397,13 @@ class GbdesAstrometricFitConnections( doc="Camera parameters to use for 'device' part of model", name="gbdesAstrometricFit_cameraModel", storageClass="ArrowNumpyDict", - dimensions=("instrument", "physical_filter"), + dimensions=("instrument", "skymap", "tract", "physical_filter"), + ) + camera = pipeBase.connectionTypes.Output( + doc="Camera object constructed using the per-detector part of the astrometric model", + name="gbdesAstrometricFitCamera", + storageClass="Camera", + dimensions=("instrument", "skymap", "tract", "physical_filter"), ) dcrCoefficients = pipeBase.connectionTypes.Output( doc="Per-visit coefficients for DCR correction.", @@ -412,6 +427,9 @@ def __init__(self, *, config=None): self.prerequisiteInputs.remove("inputCameraModel") if not self.config.saveCameraModel: self.outputs.remove("outputCameraModel") + if not self.config.saveCameraObject: + self.prerequisiteInputs.remove("inputCamera") + self.outputs.remove("camera") class GbdesAstrometricFitConfig( @@ -541,6 +559,14 @@ class GbdesAstrometricFitConfig( doc="Save the 'device' part of the model to be used as input in future runs.", default=False, ) + buildCamera = pexConfig.ConfigurableField( + target=BuildCameraFromAstrometryTask, doc="Subtask to build camera from astrometric model." + ) + saveCameraObject = pexConfig.Field( + dtype=bool, + doc="Build and output an lsst.afw.cameraGeom.Camera object using the fit per-detector model.", + default=False, + ) def setDefaults(self): # Use only stars because aperture fluxes of galaxies are biased and @@ -600,6 +626,13 @@ def validate(self): "saveCameraModel and useInputCameraModel cannot both be true.", ) + if self.saveCameraObject and (self.exposurePolyOrder != 1): + raise pexConfig.FieldValidationError( + GbdesAstrometricFitConfig.saveCameraObject, + self, + "If saveCameraObject is True, exposurePolyOrder must be set to 1.", + ) + class GbdesAstrometricFitTask(pipeBase.PipelineTask): """Calibrate the WCS across multiple visits of the same field using the @@ -613,6 +646,8 @@ def __init__(self, **kwargs): super().__init__(**kwargs) self.makeSubtask("sourceSelector") self.makeSubtask("referenceSelector") + if self.config.saveCameraObject: + self.makeSubtask("buildCamera") def runQuantum(self, butlerQC, inputRefs, outputRefs): # We override runQuantum to set up the refObjLoaders @@ -661,6 +696,8 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): butlerQC.put(output.modelParams, outputRefs.modelParams) if self.config.saveCameraModel: butlerQC.put(output.cameraModelParams, outputRefs.outputCameraModel) + if self.config.saveCameraObject: + butlerQC.put(output.camera, outputRefs.camera) if self.config.useColor: butlerQC.put(output.colorParams, outputRefs.dcrCoefficients) @@ -673,6 +710,7 @@ def run( refObjectLoader=None, inputCameraModel=None, colorCatalog=None, + inputCamera=None, ): """Run the WCS fit for a given set of visits @@ -694,6 +732,8 @@ def run( Parameters to use for the device part of the model. colorCatalog : `lsst.afw.table.SimpleCatalog` Catalog containing object coordinates and magnitudes. + inputCamera : `lsst.afw.cameraGeom.Camera`, optional + Camera to be used as template when constructing new camera. Returns ------- @@ -714,6 +754,8 @@ def run( needed as input for future runs. ``colorParams`` : `dict` [`int`, `np.ndarray`] DCR parameters fit in RA and Dec directions for each visit. + ``camera`` : `lsst.afw.cameraGeom.Camera` + Camera object constructed from the per-detector model. """ self.log.info("Gather instrument, exposure, and field info") @@ -811,12 +853,13 @@ def run( ) self.log.info("WCS fitting done") - outputWcss, cameraParams, colorParams = self._make_outputs( + outputWcss, cameraParams, colorParams, camera = self._make_outputs( wcsf, inputVisitSummaries, exposureInfo, mapTemplate, inputCameraModel=(inputCameraModel if self.config.useInputCameraModel else None), + inputCamera=(inputCamera if self.config.buildCamera else None), ) outputCatalog = wcsf.getOutputCatalog() starCatalog = wcsf.getStarCatalog() @@ -830,6 +873,7 @@ def run( modelParams=modelParams, cameraModelParams=cameraParams, colorParams=colorParams, + camera=camera, ) def _prep_sky(self, inputVisitSummaries, epoch, fieldName="Field"): @@ -1753,7 +1797,9 @@ def _make_afw_wcs(self, mapDict, centerRA, centerDec, doNormalizePixels=False, x outWCS = afwgeom.SkyWcs(frameDict) return outWCS - def _make_outputs(self, wcsf, visitSummaryTables, exposureInfo, mapTemplate, inputCameraModel=None): + def _make_outputs( + self, wcsf, visitSummaryTables, exposureInfo, mapTemplate, inputCameraModel=None, inputCamera=None + ): """Make a WCS object out of the WCS models. Parameters @@ -1793,6 +1839,8 @@ def _make_outputs(self, wcsf, visitSummaryTables, exposureInfo, mapTemplate, inp when used as input for future runs. colorFits : `dict` [`int`, `np.ndarray`], optional DCR parameters fit in RA and Dec directions for each visit. + camera : `lsst.afw.cameraGeom.Camera`, optional + Camera object constructed from the per-detector model. """ # Get the parameters of the fit models mapParams = wcsf.mapCollection.getParamDict() @@ -1805,6 +1853,26 @@ def _make_outputs(self, wcsf, visitSummaryTables, exposureInfo, mapTemplate, inp for k, params in mapParams.items(): if re.fullmatch(detectorTemplate, k): cameraParams[k] = params + if self.config.saveCameraObject: + # Get the average rotation angle of the input visits. + rotations = [ + visTable[0].visitInfo.boresightRotAngle.asRadians() for visTable in visitSummaryTables + ] + rotationAngle = np.mean(rotations) + if inputCamera is None: + raise RuntimeError( + "inputCamera must be provided to _make_outputs in order to build output camera." + ) + camera = self.buildCamera.run( + mapParams, + mapTemplate, + exposureInfo.detectors, + exposureInfo.visits, + inputCamera, + rotationAngle, + ) + else: + camera = None if self.config.useInputCameraModel: if inputCameraModel is None: raise RuntimeError( @@ -1906,7 +1974,7 @@ def _make_outputs(self, wcsf, visitSummaryTables, exposureInfo, mapTemplate, inp colorDec = np.array([colorFits[vis][1] for vis in colorVisits]) colorFits = {"visit": colorVisits, "raCoefficient": colorRA, "decCoefficient": colorDec} - return catalogs, cameraParams, colorFits + return catalogs, cameraParams, colorFits, camera def _compute_model_params(self, wcsf): """Get the WCS model parameters and covariance and convert to a @@ -2049,6 +2117,13 @@ class GbdesGlobalAstrometricFitConnections( storageClass="ArrowNumpyDict", dimensions=("instrument", "physical_filter"), ) + inputCamera = pipeBase.connectionTypes.PrerequisiteInput( + doc="Input camera to use when constructing camera from astrometric model.", + name="camera", + storageClass="Camera", + dimensions=("instrument",), + isCalibration=True, + ) outputWcs = pipeBase.connectionTypes.Output( doc=( "Per-visit world coordinate systems derived from the fitted model. These catalogs only contain " @@ -2094,6 +2169,11 @@ class GbdesGlobalAstrometricFitConnections( doc="Per-visit coefficients for DCR correction.", name="gbdesGlobalAstrometricFit_dcrCoefficients", storageClass="ArrowNumpyDict", + ) + camera = pipeBase.connectionTypes.Output( + doc="Camera object constructed using the per-detector part of the astrometric model", + name="gbdesGlobalAstrometricFitCamera", + storageClass="Camera", dimensions=("instrument", "physical_filter"), ) @@ -2112,6 +2192,9 @@ def __init__(self, *, config=None): self.prerequisiteInputs.remove("inputCameraModel") if not self.config.saveCameraModel: self.outputs.remove("outputCameraModel") + if not self.config.saveCameraObject: + self.prerequisiteInputs.remove("inputCamera") + self.outputs.remove("camera") class GbdesGlobalAstrometricFitConfig( @@ -2201,6 +2284,8 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): butlerQC.put(output.modelParams, outputRefs.modelParams) if self.config.saveCameraModel: butlerQC.put(output.cameraModelParams, outputRefs.outputCameraModel) + if self.config.saveCameraObject: + butlerQC.put(output.camera, outputRefs.camera) if self.config.useColor: butlerQC.put(output.colorParams, outputRefs.dcrCoefficients) @@ -2214,6 +2299,7 @@ def run( refObjectLoader=None, inputCameraModel=None, colorCatalog=None, + inputCamera=None, ): """Run the WCS fit for a given set of visits @@ -2237,6 +2323,8 @@ def run( Parameters to use for the device part of the model. colorCatalog : `lsst.afw.table.SimpleCatalog` Catalog containing object coordinates and magnitudes. + inputCamera : `lsst.afw.cameraGeom.Camera`, optional + Camera to be used as template when constructing new camera. Returns ------- @@ -2257,6 +2345,8 @@ def run( needed as input for future runs. ``colorParams`` : `dict` [`int`, `np.ndarray`] DCR parameters fit in RA and Dec directions for each visit. + ``camera`` : `lsst.afw.cameraGeom.Camera` + Camera object constructed from the per-detector model. """ self.log.info("Gather instrument, exposure, and field info") @@ -2335,12 +2425,13 @@ def run( ) self.log.info("WCS fitting done") - outputWcss, cameraParams, colorParams = self._make_outputs( + outputWcss, cameraParams, colorParams, camera = self._make_outputs( wcsf, inputVisitSummaries, exposureInfo, mapTemplate, inputCameraModel=(inputCameraModel if self.config.useInputCameraModel else None), + inputCamera=(inputCamera if self.config.buildCamera else None), ) outputCatalog = wcsf.getOutputCatalog() starCatalog = wcsf.getStarCatalog() @@ -2354,6 +2445,7 @@ def run( modelParams=modelParams, cameraModelParams=cameraParams, colorParams=colorParams, + camera=camera, ) def _prep_sky(self, inputVisitSummaries): diff --git a/tests/test_build_camera.py b/tests/test_build_camera.py index ac1e77fd..92a76121 100644 --- a/tests/test_build_camera.py +++ b/tests/test_build_camera.py @@ -30,13 +30,17 @@ import yaml from lsst.afw.cameraGeom.testUtils import FIELD_ANGLE, PIXELS, CameraWrapper from lsst.afw.table import ExposureCatalog -from lsst.drp.tasks.build_camera import BuildCameraConfig, BuildCameraTask, _z_function +from lsst.drp.tasks.build_camera import ( + BuildCameraFromAstrometryConfig, + BuildCameraFromAstrometryTask, + _z_function, +) from scipy.optimize import minimize TESTDIR = os.path.abspath(os.path.dirname(__file__)) -class TestBuildCamera(lsst.utils.tests.TestCase): +class TestBuildCameraFromAstrometry(lsst.utils.tests.TestCase): def setUp(self): @@ -44,7 +48,8 @@ def setUp(self): self.camera = CameraWrapper(isLsstLike=False).camera self.detectorList = [10, 20, 30, 40, 50, 60, 70, 80, 90, 100] self.visitList = np.arange(10) - self.task = BuildCameraTask(self.camera, self.detectorList, self.visitList) + self.rotationAngle = (270 * u.degree).to(u.radian).value + self.task = BuildCameraFromAstrometryTask() datadir = os.path.join(TESTDIR, "data") with open(os.path.join(datadir, "sample_wcs.yaml"), "r") as f: @@ -59,11 +64,13 @@ def setUp(self): ) self.mapParams[key] = polyCoefficients - visitSummary = ExposureCatalog.readFits(os.path.join(datadir, f"visitSummary_1176.fits")) + visitSummary = ExposureCatalog.readFits(os.path.join(datadir, "visitSummary_1176.fits")) gbdesTask = lsst.drp.tasks.gbdesAstrometricFit.GbdesAstrometricFitTask() _, self.mapTemplate = gbdesTask.make_yaml(visitSummary) - _, tangentPlaneX, tangentPlaneY = self.task._normalizeModel(self.mapParams, self.mapTemplate) + _, tangentPlaneX, tangentPlaneY = self.task._normalize_model( + self.mapParams, self.mapTemplate, self.detectorList, self.visitList + ) # There is a rotation and flip between the gbdes camera model and the # afw camera model, which is applied by-hand here. self.originalFieldAngleX = (-tangentPlaneY * u.degree).to(u.radian).value @@ -73,7 +80,7 @@ def setUp(self): self.xPix = (self.task.x + 1) * bbox.endX / 2 self.yPix = (self.task.y + 1) * bbox.endY / 2 - def test_normalizeModel(self): + def test_normalize_model(self): deviceParams = [] for element in self.mapTemplate["BAND/DEVICE"]["Elements"]: @@ -96,43 +103,45 @@ def test_normalizeModel(self): visitParams.append(identityVisitParams) expoArray = np.vstack(visitParams) - # Get the interim positions from the original device maps: - origIntX = [] - origIntY = [] + # Get the tangent plane positions from the original device maps: + origTpX = [] + origTpY = [] for deviceMap in deviceArray: nCoeffsDev = len(deviceMap) // 2 deviceDegree = lsst.drp.tasks.gbdesAstrometricFit._degreeFromNCoeffs(nCoeffsDev) intX = _z_function(deviceMap[:nCoeffsDev], self.task.x, self.task.y, order=deviceDegree) intY = _z_function(deviceMap[nCoeffsDev:], self.task.x, self.task.y, order=deviceDegree) - origIntX.append(intX) - origIntY.append(intY) + origTpX.append(intX) + origTpY.append(intY) - origIntX = np.array(origIntX).ravel() - origIntY = np.array(origIntY).ravel() + origTpX = np.array(origTpX).ravel() + origTpY = np.array(origTpY).ravel() # Get the interim positions with the new device maps: - newDeviceParams, newIntX, newIntY = self.task._normalizeModel(self.mapParams, self.mapTemplate) + _, newIntX, newIntY = self.task._normalize_model( + self.mapParams, self.mapTemplate, self.detectorList, self.visitList + ) # Now fit the per-visit parameters with the constraint that the new # tangent plane positions match the old tangent plane positions, and # then check they are sufficiently close to the old values. newExpoArrayEmp = np.zeros(expoArray.shape) for e, expo in enumerate(expoArray): - origMapX = expo[0] + origIntX * expo[1] + origIntY * expo[2] - origMapY = expo[3] + origIntX * expo[4] + origIntY * expo[5] + origMapX = expo[0] + origTpX * expo[1] + origTpY * expo[2] + origMapY = expo[3] + origTpX * expo[4] + origTpY * expo[5] def min_function(params): - tp_x = params[0] + params[1] * newIntX + params[2] * newIntY - tp_y = params[3] + params[4] * newIntX + params[5] * newIntY + tpX = params[0] + params[1] * newIntX + params[2] * newIntY + tpY = params[3] + params[4] * newIntX + params[5] * newIntY - diff = ((origMapX - tp_x)) ** 2 + ((origMapY - tp_y)) ** 2 + diff = ((origMapX - tpX)) ** 2 + ((origMapY - tpY)) ** 2 return diff.sum() def jac(params): - tp_x = params[0] + params[1] * newIntX + params[2] * newIntY - tp_y = params[3] + params[4] * newIntX + params[5] * newIntY - dX = 2 * (origMapX - tp_x) - dY = 2 * (origMapY - tp_y) + tpX = params[0] + params[1] * newIntX + params[2] * newIntY + tpY = params[3] + params[4] * newIntX + params[5] * newIntY + dX = 2 * (origMapX - tpX) + dY = 2 * (origMapY - tpY) jacobian = [ -dX.sum(), -(dX * newIntX).sum(), @@ -158,11 +167,18 @@ def jac(params): def test_run_with_basic_model(self): - config = BuildCameraConfig() + config = BuildCameraFromAstrometryConfig() config.modelSplitting = "basic" - task = BuildCameraTask(self.camera, self.detectorList, self.visitList, config=config) - - camera = task.run(self.mapParams, self.mapTemplate) + task = BuildCameraFromAstrometryTask(config=config) + + camera = task.run( + self.mapParams, + self.mapTemplate, + self.detectorList, + self.visitList, + self.camera, + self.rotationAngle, + ) testX, testY = [], [] for dev in camera: faX, faY = ( @@ -180,11 +196,18 @@ def test_run_with_basic_model(self): def test_run_with_splitModel(self): - config = BuildCameraConfig() + config = BuildCameraFromAstrometryConfig() config.modelSplitting = "physical" config.modelSplittingTolerance = 1e-6 - task = BuildCameraTask(self.camera, self.detectorList, self.visitList, config=config) - camera = task.run(self.mapParams, self.mapTemplate) + task = BuildCameraFromAstrometryTask(config=config) + camera = task.run( + self.mapParams, + self.mapTemplate, + self.detectorList, + self.visitList, + self.camera, + self.rotationAngle, + ) testX, testY = [], [] for dev in camera: @@ -198,7 +221,8 @@ def test_run_with_splitModel(self): testX = np.concatenate(testX) testY = np.concatenate(testY) - # This does not reconstruct the input field angles perfectly. + # The "physical" model splitting is not expected to + # reconstruct the input field angles perfectly. self.assertFloatsAlmostEqual(self.originalFieldAngleX, testX, atol=1e-4) self.assertFloatsAlmostEqual(self.originalFieldAngleY, testY, atol=1e-4)