Coverage for python/lsst/drp/tasks/build_camera.py: 97%

1# This file is part of drp_tasks. 

2# 

3# LSST Data Management System 

4# This product includes software developed by the 

5# LSST Project (http://www.lsst.org/). 

6# See COPYRIGHT file at the top of the source tree. 

7# 

8# This program is free software: you can redistribute it and/or modify 

9# it under the terms of the GNU General Public License as published by 

10# the Free Software Foundation, either version 3 of the License, or 

11# (at your option) any later version. 

12# 

13# This program is distributed in the hope that it will be useful, 

14# but WITHOUT ANY WARRANTY; without even the implied warranty of 

15# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 

16# GNU General Public License for more details. 

17# 

18# You should have received a copy of the LSST License Statement and 

19# the GNU General Public License along with this program. If not, 

20# see <https://www.lsstcorp.org/LegalNotices/>. 

21# 

22import re 

23from collections import defaultdict 

24 

25import astshim as ast 

26import numpy as np 

27import scipy.optimize 

28 

29import lsst.pex.config 

30from lsst.afw.cameraGeom import FIELD_ANGLE, FOCAL_PLANE, PIXELS, Camera 

31from lsst.afw.geom import TransformPoint2ToPoint2 

32from lsst.geom import degrees 

33from lsst.pipe.base import Struct, Task 

34 

35__all__ = [ 

36 "BuildCameraFromAstrometryConfig", 

37 "BuildCameraFromAstrometryTask", 

38] 

39 

40# Set up global variable to use in minimization callback. 

41nFunctionEval = 0 

42 

43 

44def _z_function(params, x, y, order=4): 

45 """Convenience function for calculating the value of a 2-d polynomial with 

46 coefficients given by an array of parameters. 

47 

48 Parameters 

49 ---------- 

50 params : `np.ndarray` 

51 Coefficientss for polynomial with terms ordered 

52 [1, x, y, x^2, xy, y^2, ...]. 

53 x, y : `np.ndarray` 

54 x and y values at which to evaluate the 2-d polynomial. 

55 order : `int`, optional 

56 Order of the given polynomial. 

57 

58 Returns 

59 ------- 

60 z : `np.ndarray` 

61 Value of polynomial evaluated at (x, y). 

62 """ 

63 coeffs = np.zeros((order + 1, order + 1)) 

64 c = 0 

65 for i in range(order + 1): 

66 for j in range(i + 1): 

67 coeffs[j, i - j] = params[c] 

68 c += 1 

69 z = np.polynomial.polynomial.polyval2d(x, y, coeffs.T) 

70 return z 

71 

72 

73def _z_function_dx(params, x, y, order=4): 

74 """Convenience function for calculating the derivative with respect to x of 

75 a 2-d polynomial with coefficients given by an array of parameters. 

76 

77 Parameters 

78 ---------- 

79 params : `np.ndarray` 

80 Coefficientss for polynomial with terms ordered 

81 [1, x, y, x^2, xy, y^2, ...]. 

82 x, y : `np.ndarray` 

83 x and y values at which to evaluate the 2-d polynomial. 

84 order : `int`, optional 

85 Order of the given polynomial. 

86 

87 Returns 

88 ------- 

89 z : `np.ndarray` 

90 Derivative of polynomial evaluated at (x, y). 

91 """ 

92 coeffs = np.zeros((order + 1, order + 1)) 

93 c = 0 

94 for i in range(order + 1): 

95 for j in range(i + 1): 

96 coeffs[j, i - j] = params[c] 

97 c += 1 

98 derivCoeffs = np.zeros(coeffs.shape) 

99 derivCoeffs[:, :-1] = coeffs[:, 1:] 

100 derivCoeffs[:, :-1] *= np.arange(1, order + 1) 

101 

102 z = np.polynomial.polynomial.polyval2d(x, y, derivCoeffs.T) 

103 return z 

104 

105 

106def _z_function_dy(params, x, y, order=4): 

107 """Convenience function for calculating the derivative with respect to y of 

108 a 2-d polynomial with coefficients given by an array of parameters. 

109 

110 Parameters 

111 ---------- 

112 params : `np.ndarray` 

113 Coefficientss for polynomial with terms ordered 

114 [1, x, y, x^2, xy, y^2, ...]. 

115 x, y : `np.ndarray` 

116 x and y values at which to evaluate the 2-d polynomial. 

117 order : `int`, optional 

118 Order of the given polynomial. 

119 

120 Returns 

121 ------- 

122 z : `np.ndarray` 

123 Derivative of polynomial evaluated at (x, y). 

124 """ 

125 coeffs = np.zeros((order + 1, order + 1)) 

126 c = 0 

127 for i in range(order + 1): 

128 for j in range(i + 1): 

129 coeffs[j, i - j] = params[c] 

130 c += 1 

131 derivCoeffs = np.zeros(coeffs.shape) 

132 derivCoeffs[:-1] = coeffs[1:] 

133 derivCoeffs[:-1] *= np.arange(1, order + 1).reshape(-1, 1) 

134 

135 z = np.polynomial.polynomial.polyval2d(x, y, derivCoeffs.T) 

136 return z 

137 

138 

139class BuildCameraFromAstrometryConfig(lsst.pex.config.Config): 

140 """Configuration for BuildCameraTask.""" 

141 

142 tangentPlaneDegree = lsst.pex.config.Field( 

143 dtype=int, 

144 doc=( 

145 "Order of polynomial mapping between the focal plane and the tangent plane. Only used if " 

146 "modelSplitting='physical'." 

147 ), 

148 default=9, 

149 ) 

150 focalPlaneDegree = lsst.pex.config.Field( 

151 dtype=int, 

152 doc=( 

153 "Order of polynomial mapping between the pixels and the focal plane. Only used if " 

154 "modelSplitting='physical'." 

155 ), 

156 default=3, 

157 ) 

158 modelSplitting = lsst.pex.config.ChoiceField( 

159 dtype=str, 

160 doc="How to split the camera model into pixel to focal plane and focal plane to tangent plane parts.", 

161 default="basic", 

162 allowed={ 

163 "basic": "Put all the mapping except scaling into the pixel to focal plane part", 

164 "physical": ( 

165 "Fit a more physical mapping where as much as possible of the model variability goes" 

166 "into the focal plane to tangent plane part of the model, to imitate the effects of " 

167 "the telescope optics" 

168 ), 

169 }, 

170 ) 

171 plateScale = lsst.pex.config.Field( 

172 dtype=float, 

173 doc=("Scaling between camera coordinates in mm and angle on the sky in" " arcsec."), 

174 default=20.005867576692737, 

175 ) 

176 astInversionTolerance = lsst.pex.config.Field( 

177 dtype=float, 

178 doc="Tolerance for AST map inversion.", 

179 default=0.005, 

180 ) 

181 astInversionMaxIter = lsst.pex.config.Field( 

182 dtype=int, doc="Maximum iterations fro AST map inversion.", default=10 

183 ) 

184 modelSplittingTolerance = lsst.pex.config.Field( 

185 dtype=float, 

186 doc="Average error in model splitting minimization acceptable for convergence.", 

187 default=1e-8, 

188 ) 

189 

190 

191class BuildCameraFromAstrometryTask(Task): 

192 """Build an `lsst.afw.cameraGeom.Camera` object out of the `gbdes` 

193 polynomials mapping from pixels to the tangent plane. 

194 

195 Parameters 

196 ---------- 

197 

198 """ 

199 

200 ConfigClass = BuildCameraFromAstrometryConfig 

201 _DefaultName = "buildCameraFromAstrometry" 

202 

203 def __init__(self, **kwargs): 

204 super().__init__(**kwargs) 

205 

206 # The gbdes model normalizes the pixel positions to the range -1 - +1. 

207 self.X = np.linspace(-1, 1, 10) 

208 self.Y = np.linspace(-1, 1, 10) 

209 x, y = np.meshgrid(self.X, self.Y) 

210 self.x = x.ravel() 

211 self.y = y.ravel() 

212 

213 def run(self, mapParams, mapTemplate, detectorList, visitList, inputCamera): 

214 """Convert the model parameters into a Camera object. 

215 

216 Parameters 

217 ---------- 

218 mapParams : `dict` [`float`] 

219 Parameters that describe the per-detector and per-visit parts of 

220 the astrometric model. 

221 mapTemplate : `dict` 

222 Dictionary describing the format of the astrometric model, 

223 following the convention in `gbdes`. 

224 detectorList : `list` [`int`] 

225 List of detector ids. 

226 visitList : `list` [`int`] 

227 List of ids for visits that were used to train the input model. 

228 camera : `lsst.afw.cameraGeom.Camera` 

229 Camera object from which to take pupil function, name, and other 

230 properties. 

231 rotationAngle : `float` 

232 Value in radians of the average rotation angle of the input visits. 

233 

234 Returns 

235 ------- 

236 camera : `lsst.afw.cameraGeom.Camera` 

237 Camera object with transformations set by the input mapParams. 

238 """ 

239 

240 # Align the model with the obs_* model camera. 

241 recenteredParams = self._recenter_model(mapParams, mapTemplate, detectorList, inputCamera) 

242 

243 if self.config.modelSplitting == "basic": 

244 # Put all of the camera distortion into the pixels->focal plane 

245 # part of the distortion model, with the focal plane->tangent plane 

246 # part only used for scaling between the focal plane and sky. 

247 pixToFocalPlane, focalPlaneToTangentPlane = self._basic_model( 

248 recenteredParams.newDeviceArray, detectorList 

249 ) 

250 

251 else: 

252 # Fit two polynomials, such that the first describes the pixel-> 

253 # focal plane part of the model, and the second describes the focal 

254 # plane->tangent plane part of the model, with the goal of 

255 # generating a more physically-motivated distortion model. 

256 pixToFocalPlane, focalPlaneToTangentPlane = self._split_model( 

257 recenteredParams.newTPx, recenteredParams.newTPy, detectorList 

258 ) 

259 

260 xScale = int(mapTemplate["BAND/DEVICE/poly"]["XMax"]) - int(mapTemplate["BAND/DEVICE/poly"]["XMin"]) 

261 yScale = int(mapTemplate["BAND/DEVICE/poly"]["YMax"]) - int(mapTemplate["BAND/DEVICE/poly"]["YMin"]) 

262 

263 # Turn the mappings into a Camera object. 

264 camera = self._translate_to_afw( 

265 pixToFocalPlane, 

266 focalPlaneToTangentPlane, 

267 detectorList, 

268 inputCamera, 

269 xScale, 

270 yScale, 

271 ) 

272 

273 return camera 

274 

275 def _recenter_model(self, mapParams, mapTemplate, detectorList, camera, transform="affine"): 

276 """Normalize and recenter the camera mappings, such that they 

277 correspond to the average visit. 

278 

279 In order to break degeneracies, `gbdes` sets one visit to have the 

280 identity for the per-visit part of its model. This visit is thus 

281 imprinted in the per-detector part of the model. This function 

282 normalizes and recenters the model such that the parameters correspond 

283 instead to the average per-visit behavior. This step assumes that the 

284 per-visit part of the model is a first-order polynomial, in which case 

285 this conversion can be done without loss of information. If `normalize` 

286 is false, the model is just recentered so that the average offset is 

287 zero. 

288 

289 Parameters 

290 ---------- 

291 mapParams : `dict` [`float`] 

292 Parameters that describe the per-detector and per-visit parts of 

293 the astrometric model. 

294 mapTemplate : `dict` 

295 Dictionary describing the format of the astrometric model, 

296 following the convention in `gbdes`. 

297 detectorList : `list` [`int`] 

298 List of detector ids. 

299 camera: `lsst.afw.cameraGeom.Camera` 

300 Existing camera model ontoo which the new model will be normalized. 

301 transform: `str` 

302 Method for recentering the camera. 

303 - ``"affine"``: The camera parameters are shifted by a linear 

304 function to minimize the difference between this camera model 

305 and a preexisting model. 

306 - ``"rigid"`: The camera parameters are shifted by a 

307 translation and rotation to minimize the difference between 

308 this camera model and a preexisting model. TODO: DM-56786 This 

309 needs to be combined with a model for other causes of shifts, 

310 like refraction or camera flexure, in order to be 

311 scientifically useful. 

312 

313 Returns 

314 ------- 

315 result : `lsst.pipe.base.Struct` 

316 A struct object with the following attributes: 

317 ``newDeviceArray`` : `np.ndarray` 

318 Array of NxM, where N is the number of detectors, and M is the 

319 number of coefficients for each per-detector mapping. 

320 ``newTPx``, ``newTPy`` : `np.ndarray` 

321 Projection of `self.x` and `self.y` onto the tangent plane for 

322 each detector, given the normalized mapping. 

323 ``inputCoords`` : `np.ndarray` 

324 Projection of `self.x` and `self.y` onto the tangent plane for 

325 each detector using the input camera. 

326 ``targetCoords`` : `np.ndarray` 

327 Projection of `self.x` and `self.y` onto the tangent plane for 

328 each detector using the target preexisting camera. 

329 """ 

330 # Get the per-device part of the model. 

331 deviceParams = [] 

332 for element in mapTemplate["BAND/DEVICE"]["Elements"]: 

333 for detector in detectorList: 

334 detectorTemplate = element.replace("DEVICE", str(detector)) 

335 detectorTemplate = detectorTemplate.replace("BAND", ".+") 

336 for k, params in mapParams.items(): 

337 if re.fullmatch(detectorTemplate, k): 

338 deviceParams.append(params) 

339 deviceArray = np.vstack(deviceParams) 

340 

341 inputTPx = [] 

342 inputTPy = [] 

343 targetTPx = [] 

344 targetTPy = [] 

345 xMax = max([det.getBBox().endX for det in camera]) - 1 

346 yMax = max([det.getBBox().endY for det in camera]) - 1 

347 

348 for d, deviceMap in enumerate(deviceArray): 

349 detector = detectorList[d] 

350 obsDet = camera[detector] 

351 detBBox = obsDet.getBBox() 

352 

353 nCoeffsDev = len(deviceMap) // 2 

354 deviceDegree = int(-1.5 + 0.5 * (1 + 8 * nCoeffsDev) ** 0.5) 

355 

356 x = ((detBBox.endX - 1) / xMax) * (self.x + 1) - 1 

357 y = ((detBBox.endY - 1) / yMax) * (self.y + 1) - 1 

358 

359 intX = _z_function(deviceMap[:nCoeffsDev], x, y, order=deviceDegree) 

360 intY = _z_function(deviceMap[nCoeffsDev:], x, y, order=deviceDegree) 

361 inputTPx.append(intX) 

362 inputTPy.append(intY) 

363 

364 # Get the tangent plane mapping for the same points from the target 

365 # camera model. 

366 dX = np.linspace(detBBox.beginX, detBBox.endX - 1, len(self.X)) 

367 dY = np.linspace(detBBox.beginY, detBBox.endY - 1, len(self.Y)) 

368 detX, detY = np.meshgrid(dX, dY) 

369 detX = detX.ravel() 

370 detY = detY.ravel() 

371 mapping = obsDet.getTransform(PIXELS, FIELD_ANGLE).getMapping() 

372 faX, faY = mapping.applyForward(np.array([detX, detY])) 

373 targetTPx.append(np.degrees(faX)) 

374 targetTPy.append(np.degrees(faY)) 

375 

376 inputCoords = np.array([np.concatenate(inputTPx), np.concatenate(inputTPy)]) 

377 targetCoords = np.array([np.concatenate(targetTPx), np.concatenate(targetTPy)]) 

378 

379 if transform == "rigid": 

380 # Find the best rigid transformation matching the new model to the 

381 # input camera model using the method described in 

382 # https://igl.ethz.ch/projects/ARAP/svd_rot.pdf. 

383 inputAve = inputCoords.mean(axis=1) 

384 targetAve = targetCoords.mean(axis=1) 

385 

386 xVec = inputCoords - inputAve.reshape(-1, 1) 

387 yVec = targetCoords - targetAve.reshape(-1, 1) 

388 

389 # Now do SVD: 

390 SMatrix = xVec.dot(yVec.T) 

391 U, S, VTrans = np.linalg.svd(SMatrix) 

392 

393 rotationMatrix = (VTrans.T).dot(U.T) 

394 translation = targetAve.reshape(-1, 1) - rotationMatrix.dot(inputAve.reshape(-1, 1)) 

395 

396 # Convert the tangent plane coordinates already calculated. 

397 matchedTPCoords = rotationMatrix.dot(inputCoords) + translation 

398 

399 # Convert the device map coefficients. 

400 nCoeffsDev = deviceArray.shape[1] // 2 

401 newDeviceArray = np.zeros(deviceArray.shape) 

402 newDeviceArray[:, :nCoeffsDev] = ( 

403 deviceArray[:, :nCoeffsDev] * rotationMatrix[0, 0] 

404 + deviceArray[:, nCoeffsDev:] * rotationMatrix[0, 1] 

405 ) 

406 newDeviceArray[:, 0] += translation[0] 

407 newDeviceArray[:, nCoeffsDev:] = ( 

408 deviceArray[:, :nCoeffsDev] * rotationMatrix[1, 0] 

409 + deviceArray[:, nCoeffsDev:] * rotationMatrix[1, 1] 

410 ) 

411 newDeviceArray[:, nCoeffsDev] += translation[1] 

412 

413 elif transform == "affine": 413 ↛ 439line 413 didn't jump to line 439 because the condition on line 413 was always true

414 # Find the transformation parameters by minimizing 

415 # (target coords - transforms X input coords)**2. 

416 x = np.vstack([np.ones_like(inputCoords[0]), inputCoords]) 

417 A = x.dot(x.T) 

418 AInv = np.linalg.inv(A) 

419 bx = (x * targetCoords[0]).sum(axis=1) 

420 by = (x * targetCoords[1]).sum(axis=1) 

421 transformX = AInv.dot(bx) 

422 transformY = AInv.dot(by) 

423 

424 matchedTPCoords = np.array([transformX, transformY]).dot(x) 

425 

426 # Convert the device map coefficients. 

427 nCoeffsDev = deviceArray.shape[1] // 2 

428 newDeviceArray = np.zeros(deviceArray.shape) 

429 newDeviceArray[:, :nCoeffsDev] = ( 

430 deviceArray[:, :nCoeffsDev] * transformX[1] + deviceArray[:, nCoeffsDev:] * transformX[2] 

431 ) 

432 newDeviceArray[:, 0] += transformX[0] 

433 newDeviceArray[:, nCoeffsDev:] = ( 

434 deviceArray[:, :nCoeffsDev] * transformY[1] + deviceArray[:, nCoeffsDev:] * transformY[2] 

435 ) 

436 newDeviceArray[:, nCoeffsDev] += transformY[0] 

437 

438 else: 

439 raise RuntimeError("Transform '%s' is not supported", transform) 

440 

441 newTPx, newTPy = matchedTPCoords 

442 

443 return Struct( 

444 newDeviceArray=newDeviceArray, 

445 newTPx=newTPx, 

446 newTPy=newTPy, 

447 inputCoords=inputCoords, 

448 targetCoords=targetCoords, 

449 ) 

450 

451 def _basic_model(self, modelParameters, detectorList): 

452 """This will just convert the pix->fp parameters into the right format, 

453 and return an identity mapping for the fp->tp part. 

454 

455 Parameters 

456 ---------- 

457 modelParameters : `np.ndarray` 

458 Array of NxM, where N is the number of detectors, and M is the 

459 number of coefficients for each per-detector mapping. 

460 detectorList : `list` [`int`] 

461 List of detector ids. 

462 

463 Returns 

464 ------- 

465 pixelsToFocalPlane : `dict` [`int`, [`string`, `float`]] 

466 Dictionary giving the per-detector pixel to focal plane mapping, 

467 with keys for each detector giving the x and y transformation 

468 polynomials. 

469 focalPlaneToTangentPlane : `dict` [`string`, `float`] 

470 Dictionary giving the focal plane to tangent plane mapping. 

471 """ 

472 

473 nCoeffsFp = modelParameters.shape[1] // 2 

474 pixelsToFocalPlane = defaultdict(dict) 

475 for d, det in enumerate(detectorList): 

476 pixelsToFocalPlane[det]["x"] = modelParameters[d][:nCoeffsFp] 

477 pixelsToFocalPlane[det]["y"] = modelParameters[d][nCoeffsFp:] 

478 

479 focalPlaneToTangentPlane = {"x": [0, 1, 0], "y": [0, 0, 1]} 

480 

481 return pixelsToFocalPlane, focalPlaneToTangentPlane 

482 

483 def _split_model(self, targetX, targetY, detectorList, pixToFpGuess=None, fpToTpGuess=None): 

484 """Fit a two-step model, with one polynomial per detector fitting the 

485 pixels to focal plane part, followed by a polynomial fitting the focal 

486 plane to tangent plane part. 

487 

488 The linear parameters of the focal plane to tangent plane model are 

489 fixed to be an identity map to remove degeneracies. 

490 

491 Parameters 

492 ---------- 

493 targetX, targetY : `np.ndarray` 

494 Target x and y values in the tangent plane. 

495 detectorList : `list` [`int`] 

496 List of detector ids. 

497 pixToFpGuess : `dict` [`int`, [`string`, `float`]] 

498 Initial guess for the pixels to focal plane mapping to use in the 

499 model fit, in the form of a dictionary giving the per-detector 

500 pixel to focal plane mapping, with keys for each detector giving 

501 the x and y transformation polynomials. 

502 focalPlaneToTangentPlane : `dict` [`string`, `float`] 

503 Initial guess for the focal plane to tangent plane mapping to use 

504 in the model fit, in the form of a dictionary giving the focal 

505 plane to tangent plane mapping. 

506 

507 Returns 

508 ------- 

509 pixelsToFocalPlane : `dict` [`int`, [`string`, `float`]] 

510 Dictionary giving the per-detector pixel to focal plane mapping, 

511 with keys for each detector giving the x and y transformation 

512 polynomials. 

513 focalPlaneToTangentPlane : `dict` [`string`, `float`] 

514 Dictionary giving the focal plane to tangent plane mapping. 

515 """ 

516 

517 tpDegree = self.config.tangentPlaneDegree 

518 fpDegree = self.config.focalPlaneDegree 

519 

520 nCoeffsFP = int((fpDegree + 2) * (fpDegree + 1) / 2) 

521 nCoeffsFPTot = len(detectorList) * nCoeffsFP * 2 

522 

523 nCoeffsTP = int((tpDegree + 2) * (tpDegree + 1) / 2) 

524 # The constant and linear terms will be fixed to remove degeneracy with 

525 # the pix->fp part of the model 

526 nCoeffsFixed = 3 

527 nCoeffsFree = nCoeffsTP - nCoeffsFixed 

528 fixedParams = np.zeros(nCoeffsTP * 2) 

529 fixedParams[1] = 1 

530 fixedParams[nCoeffsTP + 2] = 1 

531 

532 nX = len(self.x) 

533 # We need an array of the form [1, x, y, x^2, xy, y^2, ...] to use when 

534 # calculating derivatives. 

535 xyArray = np.zeros((nCoeffsFP, nX)) 

536 c = 0 

537 for ii in range(fpDegree + 1): 

538 for jj in range(ii + 1): 

539 coeffs = np.zeros((fpDegree + 1, fpDegree + 1)) 

540 coeffs[jj, ii - jj] = 1 

541 xyArray[c] = np.polynomial.polynomial.polyval2d(self.x, self.y, coeffs.T) 

542 c += 1 

543 

544 def two_part_function(params): 

545 # The function giving the split model. 

546 fpXAll = [] 

547 fpYAll = [] 

548 for i in range(len(detectorList)): 

549 fpX = _z_function( 

550 params[2 * nCoeffsFP * i : (2 * i + 1) * nCoeffsFP], self.x, self.y, order=fpDegree 

551 ) 

552 fpY = _z_function( 

553 params[(2 * i + 1) * nCoeffsFP : 2 * nCoeffsFP * (1 + i)], self.x, self.y, order=fpDegree 

554 ) 

555 fpXAll.append(fpX) 

556 fpYAll.append(fpY) 

557 fpXAll = np.array(fpXAll).ravel() 

558 fpYAll = np.array(fpYAll).ravel() 

559 

560 tpParams = fixedParams.copy() 

561 tpParams[nCoeffsFixed:nCoeffsTP] = params[nCoeffsFPTot : nCoeffsFPTot + nCoeffsFree] 

562 tpParams[nCoeffsFixed + nCoeffsTP :] = params[nCoeffsFPTot + nCoeffsFree :] 

563 

564 tpX = _z_function(tpParams[:nCoeffsTP], fpXAll, fpYAll, order=tpDegree) 

565 tpY = _z_function(tpParams[nCoeffsTP:], fpXAll, fpYAll, order=tpDegree) 

566 return tpX, tpY 

567 

568 def min_function(params): 

569 # The least-squares function to be minimized. 

570 tpX, tpY = two_part_function(params) 

571 diff = ((targetX - tpX)) ** 2 + ((targetY - tpY)) ** 2 

572 result = diff.sum() 

573 return result 

574 

575 def jac(params): 

576 # The Jacobian of min_function, which is needed for the minimizer. 

577 tpX, tpY = two_part_function(params) 

578 dX = -2 * (targetX - tpX) 

579 dY = -2 * (targetY - tpY) 

580 jacobian = np.zeros(len(params)) 

581 fpParams = params[:nCoeffsFPTot] 

582 tpParams = fixedParams.copy() 

583 tpParams[nCoeffsFixed:nCoeffsTP] = params[nCoeffsFPTot : nCoeffsFPTot + nCoeffsFree] 

584 tpParams[nCoeffsFixed + nCoeffsTP :] = params[nCoeffsFPTot + nCoeffsFree :] 

585 fpXAll = [] 

586 fpYAll = [] 

587 # Fill in the derivatives for the pix->fp terms. 

588 for i in range(len(detectorList)): 

589 dX_i = dX[nX * i : nX * (i + 1)] 

590 dY_i = dY[nX * i : nX * (i + 1)] 

591 fpX = _z_function( 

592 fpParams[2 * nCoeffsFP * i : (2 * i + 1) * nCoeffsFP], self.x, self.y, order=fpDegree 

593 ) 

594 fpY = _z_function( 

595 fpParams[(2 * i + 1) * nCoeffsFP : 2 * nCoeffsFP * (1 + i)], 

596 self.x, 

597 self.y, 

598 order=fpDegree, 

599 ) 

600 fpXAll.append(fpX) 

601 fpYAll.append(fpY) 

602 

603 dTpX_dFpX = _z_function_dx(tpParams[:nCoeffsTP], fpX, fpY, order=tpDegree) 

604 dTpX_dFpY = _z_function_dy(tpParams[:nCoeffsTP], fpX, fpY, order=tpDegree) 

605 dTpY_dFpX = _z_function_dx(tpParams[nCoeffsTP:], fpX, fpY, order=tpDegree) 

606 dTpY_dFpY = _z_function_dy(tpParams[nCoeffsTP:], fpX, fpY, order=tpDegree) 

607 

608 dTpX_part = np.concatenate([xyArray * dTpX_dFpX, xyArray * dTpX_dFpY]) 

609 dTpY_part = np.concatenate([xyArray * dTpY_dFpX, xyArray * dTpY_dFpY]) 

610 jacobian_i = (dX_i * dTpX_part + dY_i * dTpY_part).sum(axis=1) 

611 jacobian[2 * nCoeffsFP * i : 2 * nCoeffsFP * (i + 1)] = jacobian_i 

612 

613 fpXAll = np.array(fpXAll).ravel() 

614 fpYAll = np.array(fpYAll).ravel() 

615 

616 # Fill in the derivatives for the fp->tp terms. 

617 for j in range(nCoeffsFree): 

618 tParams = np.zeros(nCoeffsTP) 

619 tParams[nCoeffsFixed + j] = 1 

620 tmpZ = _z_function(tParams, fpXAll, fpYAll, order=tpDegree) 

621 jacobian[nCoeffsFPTot + j] = (dX * tmpZ).sum() 

622 jacobian[nCoeffsFPTot + nCoeffsFree + j] = (dY * tmpZ).sum() 

623 return jacobian 

624 

625 def hessian(params): 

626 # The Hessian of min_function, which is needed for the minimizer. 

627 hessian = np.zeros((len(params), len(params))) 

628 

629 fpParams = params[:nCoeffsFPTot] 

630 tpParams = fixedParams.copy() 

631 tpParams[nCoeffsFixed:nCoeffsTP] = params[nCoeffsFPTot : nCoeffsFPTot + nCoeffsFree] 

632 tpParams[nCoeffsFixed + nCoeffsTP :] = params[nCoeffsFPTot + nCoeffsFree :] 

633 fpXAll = [] 

634 fpYAll = [] 

635 

636 # Loop over the detectors to fill in the d(pix->fp)**2 and 

637 # d(pix->fp)d(fp->tp) cross terms. 

638 for i in range(len(detectorList)): 

639 fpX = _z_function( 

640 fpParams[2 * nCoeffsFP * i : (2 * i + 1) * nCoeffsFP], self.x, self.y, order=fpDegree 

641 ) 

642 fpY = _z_function( 

643 fpParams[(2 * i + 1) * nCoeffsFP : 2 * nCoeffsFP * (1 + i)], 

644 self.x, 

645 self.y, 

646 order=fpDegree, 

647 ) 

648 fpXAll.append(fpX) 

649 fpYAll.append(fpY) 

650 dTpX_dFpX = _z_function_dx(tpParams[:nCoeffsTP], fpX, fpY, order=tpDegree) 

651 dTpX_dFpY = _z_function_dy(tpParams[:nCoeffsTP], fpX, fpY, order=tpDegree) 

652 dTpY_dFpX = _z_function_dx(tpParams[nCoeffsTP:], fpX, fpY, order=tpDegree) 

653 dTpY_dFpY = _z_function_dy(tpParams[nCoeffsTP:], fpX, fpY, order=tpDegree) 

654 

655 dTpX_part = np.concatenate([xyArray * dTpX_dFpX, xyArray * dTpX_dFpY]) 

656 dTpY_part = np.concatenate([xyArray * dTpY_dFpX, xyArray * dTpY_dFpY]) 

657 

658 for k in range(2 * nCoeffsFP): 

659 for m in range(2 * nCoeffsFP): 

660 hessian[2 * nCoeffsFP * i + k, 2 * nCoeffsFP * i + m] = ( 

661 2 * (dTpX_part[k] * dTpX_part[m] + dTpY_part[k] * dTpY_part[m]).sum() 

662 ) 

663 

664 for j in range(nCoeffsFree): 

665 tmpParams = np.zeros(nCoeffsTP) 

666 tmpParams[nCoeffsFixed + j] = 1 

667 tmpZ = _z_function(tmpParams, fpX, fpY, order=tpDegree) 

668 

669 hessX = 2 * (dTpX_part[k] * tmpZ).sum() 

670 hessY = 2 * (dTpY_part[k] * tmpZ).sum() 

671 # dTP_x part 

672 hessian[2 * nCoeffsFP * i + k, nCoeffsFPTot + j] = hessX 

673 hessian[nCoeffsFPTot + j, 2 * nCoeffsFP * i + k] = hessX 

674 # dTP_y part 

675 hessian[2 * nCoeffsFP * i + k, nCoeffsFPTot + nCoeffsFree + j] = hessY 

676 hessian[nCoeffsFPTot + nCoeffsFree + j, 2 * nCoeffsFP * i + k] = hessY 

677 

678 fpXAll = np.array(fpXAll).ravel() 

679 fpYAll = np.array(fpYAll).ravel() 

680 tmpZArray = np.zeros((nCoeffsFree, nX * len(detectorList))) 

681 # Finally, get the d(fp->tp)**2 terms 

682 for j in range(nCoeffsFree): 

683 tParams = np.zeros(nCoeffsTP) 

684 tParams[nCoeffsFixed + j] = 1 

685 tmpZArray[j] = _z_function(tParams, fpXAll, fpYAll, order=tpDegree) 

686 for j in range(nCoeffsFree): 

687 for m in range(nCoeffsFree): 

688 # X-Y cross terms are zero 

689 hess = 2 * (tmpZArray[j] * tmpZArray[m]).sum() 

690 hessian[nCoeffsFPTot + j, nCoeffsFPTot + m] = hess 

691 hessian[nCoeffsFPTot + nCoeffsFree + j, nCoeffsFPTot + nCoeffsFree + m] = hess 

692 return hessian 

693 

694 global nFunctionEval 

695 nFunctionEval = 0 

696 

697 def callbackMF(params): 

698 global nFunctionEval 

699 self.log.info(f"Iteration {nFunctionEval}, function value {min_function(params)}") 

700 nFunctionEval += 1 

701 

702 initGuess = np.zeros(nCoeffsFPTot + 2 * nCoeffsFree) 

703 if pixToFpGuess: 703 ↛ 704line 703 didn't jump to line 704 because the condition on line 703 was never true

704 useVar = min(nCoeffsFP, len(pixToFpGuess[0]["x"])) 

705 for i, det in enumerate(detectorList): 

706 initGuess[2 * nCoeffsFP * i : 2 * nCoeffsFP * i + useVar] = pixToFpGuess[det]["x"][:useVar] 

707 initGuess[(2 * i + 1) * nCoeffsFP : (2 * i + 1) * nCoeffsFP + useVar] = pixToFpGuess[det][ 

708 "y" 

709 ][:useVar] 

710 if fpToTpGuess: 710 ↛ 711line 710 didn't jump to line 711 because the condition on line 710 was never true

711 useVar = min(nCoeffsTP, len(fpToTpGuess["x"])) 

712 initGuess[nCoeffsFPTot : nCoeffsFPTot + useVar - nCoeffsFixed] = fpToTpGuess["x"][ 

713 nCoeffsFixed:useVar 

714 ] 

715 initGuess[nCoeffsFPTot + nCoeffsFree : nCoeffsFPTot + nCoeffsFree + useVar - nCoeffsFixed] = ( 

716 fpToTpGuess["y"][nCoeffsFixed:useVar] 

717 ) 

718 

719 self.log.info(f"Initial value of least squares function: {min_function(initGuess)}") 

720 res = scipy.optimize.minimize( 

721 min_function, 

722 initGuess, 

723 callback=callbackMF, 

724 method="Newton-CG", 

725 jac=jac, 

726 hess=hessian, 

727 options={"xtol": self.config.modelSplittingTolerance}, 

728 ) 

729 

730 # Convert parameters to a dictionary. 

731 pixToFP = {} 

732 for i, det in enumerate(detectorList): 

733 pixToFP[det] = { 

734 "x": res.x[2 * nCoeffsFP * i : (2 * i + 1) * nCoeffsFP], 

735 "y": res.x[(2 * i + 1) * nCoeffsFP : 2 * nCoeffsFP * (1 + i)], 

736 } 

737 

738 fpToTpAll = fixedParams.copy() 

739 fpToTpAll[nCoeffsFixed:nCoeffsTP] = res.x[nCoeffsFPTot : nCoeffsFPTot + nCoeffsFree] 

740 fpToTpAll[nCoeffsFixed + nCoeffsTP :] = res.x[nCoeffsFPTot + nCoeffsFree :] 

741 fpToTp = {"x": fpToTpAll[:nCoeffsTP], "y": fpToTpAll[nCoeffsTP:]} 

742 

743 return pixToFP, fpToTp 

744 

745 def make_ast_polymap_coeffs(self, order, xCoeffs, yCoeffs): 

746 """Convert polynomial coefficients in gbdes format to AST PolyMap input 

747 format. 

748 

749 Paramaters 

750 ---------- 

751 order: `int` 

752 Polynomial order 

753 xCoeffs, yCoeffs: `list` of `float` 

754 Forward or inverse polynomial coefficients for the x and y axes 

755 of output, in this order: 

756 x0y0, x1y0, x0y1, x2y0, x1y1, x0y2, ...xNy0, xN-1y1, ...x0yN 

757 where N is the polynomial order. 

758 Returns 

759 ------- 

760 polyArray : `np.ndarray` 

761 Array formatted for AST PolyMap input. 

762 """ 

763 N = len(xCoeffs) 

764 

765 polyArray = np.zeros((2 * N, 4))