diff --git a/python/lsst/drp/tasks/build_camera.py b/python/lsst/drp/tasks/build_camera.py index 2371bc44..e4cfeb9c 100644 --- a/python/lsst/drp/tasks/build_camera.py +++ b/python/lsst/drp/tasks/build_camera.py @@ -37,7 +37,7 @@ ] # Set up global variable to use in minimization callback. -Nfeval = 0 +nFunctionEval = 0 def _z_function(params, x, y, order=4): @@ -94,11 +94,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 @@ -127,11 +127,11 @@ 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 @@ -197,7 +197,7 @@ class BuildCameraFromAstrometryTask(Task): """ ConfigClass = BuildCameraFromAstrometryConfig - _DefaultName = "buildCamera" + _DefaultName = "buildCameraFromAstrometry" def __init__(self, **kwargs): super().__init__(**kwargs) @@ -237,29 +237,29 @@ def run(self, mapParams, mapTemplate, detectorList, visitList, inputCamera, rota """ # Normalize the model. - newParams, newIntX, newIntY = self._normalizeModel(mapParams, mapTemplate, detectorList, visitList) + 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, detectorList) + 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, detectorList) + pixToFocalPlane, focalPlaneToTangentPlane = self._split_model(newTPx, newTPy, detectorList) # Turn the mappings into a Camera object. - camera = self._translateToAfw( + camera = self._translate_to_afw( pixToFocalPlane, focalPlaneToTangentPlane, detectorList, inputCamera, rotationAngle ) return camera - def _normalizeModel(self, mapParams, mapTemplate, detectorList, visitList): + def _normalize_model(self, mapParams, mapTemplate, detectorList, visitList): """Normalize the camera mappings, such that they correspond to the average visit. @@ -279,13 +279,17 @@ def _normalizeModel(self, mapParams, mapTemplate, detectorList, visitList): 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. """ @@ -314,6 +318,8 @@ def _normalizeModel(self, mapParams, mapTemplate, detectorList, visitList): 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] = ( @@ -325,24 +331,24 @@ def _normalizeModel(self, mapParams, mapTemplate, detectorList, visitList): 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 = 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, detectorList): + 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. @@ -351,6 +357,8 @@ def _basicModel(self, modelParameters, detectorList): 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 ------- @@ -362,17 +370,17 @@ def _basicModel(self, modelParameters, detectorList): 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(detectorList): - pixelsToFocalPlane[det]["x"] = modelParameters[d][:nCoeffsFP] - pixelsToFocalPlane[det]["y"] = modelParameters[d][nCoeffsFP:] + 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, detectorList, 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. @@ -384,7 +392,9 @@ def _splitModel(self, targetX, targetY, detectorList, pixToFPGuess=None, fpToTpG ---------- 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 @@ -408,16 +418,16 @@ def _splitModel(self, targetX, targetY, detectorList, pixToFPGuess=None, fpToTpG fpDegree = self.config.focalPlaneDegree nCoeffsFP = int((fpDegree + 2) * (fpDegree + 1) / 2) - nCoeffsFP_tot = len(detectorList) * nCoeffsFP * 2 + nCoeffsFPTot = len(detectorList) * nCoeffsFP * 2 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 @@ -433,26 +443,26 @@ def _splitModel(self, targetX, targetY, detectorList, pixToFPGuess=None, fpToTpG def two_part_function(params): # The function giving the split model. - intX_tot = [] - intY_tot = [] + fpXAll = [] + fpYAll = [] for i in range(len(detectorList)): - intX = _z_function( + 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): @@ -468,141 +478,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(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(detectorList)): - 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) - 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(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 + global nFunctionEval + self.log.info(f"Iteration {nFunctionEval}, function value {min_function(params)}") + nFunctionEval += 1 - initGuess = np.zeros(nCoeffsFP_tot + 2 * nCoeffsFree) - if pixToFPGuess: - useVar = min(nCoeffsFP, len(pixToFPGuess[0]["x"])) + 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][ + 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] ) @@ -625,15 +635,15 @@ def callbackMF(params): "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 @@ -647,8 +657,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) @@ -664,7 +674,7 @@ def makeAstPolyMapCoeffs(self, order, xCoeffs, yCoeffs): return polyArray - def _translateToAfw( + def _translate_to_afw( self, pixToFocalPlane, focalPlaneToTangentPlane, detectorList, inputCamera, rotationAngle ): """Convert the model parameters to a Camera object. @@ -699,13 +709,13 @@ def _translateToAfw( cameraBuilder.setPupilFactoryName(inputCamera.getPupilFactoryName()) # Convert fp->tp to AST format: - forwardCoeffs = self.makeAstPolyMapCoeffs( + forwardCoeffs = self.make_ast_polymap_coeffs( tpDegree, focalPlaneToTangentPlane["x"], focalPlaneToTangentPlane["y"] ) # Reverse rotation from input visits and flip x-axis. - cosRot = np.cos(rotationAngle) - sinRot = np.sin(rotationAngle) - rotateAndFlipCoeffs = self.makeAstPolyMapCoeffs(1, [0, cosRot, -sinRot], [0, -sinRot, cosRot]) + 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( @@ -730,7 +740,8 @@ def _translateToAfw( for detector in inputCamera: d = detector.getId() if d not in pixToFocalPlane: - # This camera will not include detectors that were not used in astrometric fit. + # This camera will not include detectors that were not used in + # astrometric fit. continue detectorBuilder = cameraBuilder.add(detector.getName(), detector.getId()) @@ -739,7 +750,7 @@ def _translateToAfw( 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( @@ -748,7 +759,7 @@ def _translateToAfw( "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( @@ -762,6 +773,6 @@ def _translateToAfw( 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 fedeac9c..b07dc20d 100644 --- a/python/lsst/drp/tasks/gbdesAstrometricFit.py +++ b/python/lsst/drp/tasks/gbdesAstrometricFit.py @@ -2102,6 +2102,7 @@ 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="gbdesAstrometricFitCamera", diff --git a/tests/test_build_camera.py b/tests/test_build_camera.py index 41819914..92a76121 100644 --- a/tests/test_build_camera.py +++ b/tests/test_build_camera.py @@ -48,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 = BuildCameraFromAstrometryTask(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: @@ -63,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 @@ -77,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"]: @@ -100,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(), @@ -164,9 +169,16 @@ def test_run_with_basic_model(self): config = BuildCameraFromAstrometryConfig() config.modelSplitting = "basic" - task = BuildCameraFromAstrometryTask(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 = ( @@ -187,8 +199,15 @@ def test_run_with_splitModel(self): config = BuildCameraFromAstrometryConfig() config.modelSplitting = "physical" config.modelSplittingTolerance = 1e-6 - task = BuildCameraFromAstrometryTask(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: