LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
EllipseSolver.cc
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010 LSST Corporation.
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 
25 #include "lsst/meas/algorithms/shapelet/EllipseSolver.h"
26 #include "lsst/meas/algorithms/shapelet/dbg.h"
27 #include "lsst/meas/algorithms/shapelet/PsiHelper.h"
28 
29 //#define JTEST
30 #define MAX_X_DEV 4.
31 
32 #ifdef JTEST
33 #include "TestHelper.h"
34 bool shouldShowTests = true;
35 bool shouldThrow = false;
36 std::string lastSuccess = "";
37 std::ostream* testout = &std::cout;
38 #endif
39 
40 namespace lsst {
41 namespace meas {
42 namespace algorithms {
43 namespace shapelet {
44 
45  struct EllipseSolver3::ESImpl3
46  {
47 
48  ESImpl3(
49  const BVec& b0, int order,
50  bool isFixedCen, bool isFixedGamma, bool isFixedMu);
51 
52  void calculateF(const DVector& x, DVector& f) const;
53  void calculateJ(const DVector& x, const DVector& f, DMatrix& J) const;
54 
55  void doF3(const DVector& x, DVector& f) const;
56  void doJ3(const DVector& x, const DVector& f, DMatrix& J) const;
57 
58  void doF2(const DVector& x, DVector& f) const;
59  void doJ2(const DVector& x, const DVector& f, DMatrix& J) const;
60 
61  void doF1(const DVector& x, DVector& f) const;
62  void doJ1(const DVector& x, const DVector& f, DMatrix& J) const;
63 
64  int b0order, bxorder, bxorderp2;
65  int b0size, bxsize, bxsizep2;
66  mutable BVec b0;
67  mutable DVector bx;
68  mutable DVector bxsave;
69  mutable DMatrix Daug;
70  mutable DMatrix Saug;
71  mutable DMatrix Taug;
72  mutable DMatrixView D;
73  mutable DMatrixView S;
74  mutable DMatrixView T;
75  mutable DVector dbdE;
76  mutable DVector Db0;
77  mutable DVector SDb0;
78  mutable DVector GDb0;
79  mutable DVector GthDb0;
80  bool fixcen, fixgam, fixmu, numeric_j, zerob11;
81  DMatrix U;
82  mutable DVector xinit;
83  mutable DVector xx;
84  mutable DVector ff;
85  mutable DMatrix jj;
86  mutable DVector x_short;
87  mutable DVector f_short;
88  mutable double fixuc,fixvc,fixg1,fixg2,fixm;
89  DMatrix Gx;
90  DMatrix Gy;
91  DMatrix Gg1;
92  DMatrix Gg2;
93  DMatrix Gth;
94  DMatrix Gmu;
95  };
96 
97  EllipseSolver3::EllipseSolver3(
98  const BVec& b0, int order,
99  bool fixcen, bool fixgam, bool fixmu) :
100  _pimpl(new ESImpl3(b0,order,fixcen,fixgam,fixmu))
101  {}
102 
103  EllipseSolver3::~EllipseSolver3()
104  { delete _pimpl; }
105 
106  void EllipseSolver3::calculateF(const DVector& x, DVector& f) const
107  { _pimpl->calculateF(x,f); }
108 
109  void EllipseSolver3::calculateJ(
110  const DVector& x, const DVector& f, DMatrix& J) const
111  {
112 #ifdef ALWAYS_NUMERIC_J
113  NLSolver::calculateJ(x,f,J);
114 #else
115  if (_pimpl->numeric_j) NLSolver::calculateJ(x,f,J);
116  else _pimpl->calculateJ(x,f,J);
117 #endif
118  }
119 
120  EllipseSolver3::ESImpl3::ESImpl3(
121  const BVec& _b0, int _order,
122  bool _fixcen, bool _fixgam, bool _fixmu) :
123  b0order(_b0.getOrder()), bxorder(_order), bxorderp2(bxorder+2),
124  b0size(_b0.size()), bxsize((bxorder+1)*(bxorder+2)/2),
125  bxsizep2((bxorderp2+1)*(bxorderp2+2)/2),
126  b0(_b0), bx(bxsize), bxsave(bxsize),
127  Daug(bxsize,bxsizep2), Saug(bxsize,bxsizep2),
128  Taug(bxsize,bxsizep2),
129  D(TMV_colRange(Daug,0,bxsize)), S(TMV_colRange(Saug,0,bxsize)),
130  T(TMV_colRange(Taug,0,bxsize)),
131  dbdE(6), Db0(bxsize), SDb0(bxsize), GDb0(bxsizep2), GthDb0(bxsizep2),
132  fixcen(_fixcen), fixgam(_fixgam), fixmu(_fixmu),
133  numeric_j(false), zerob11(true),
134  U((fixcen?0:2)+(fixgam?0:2)+(fixmu?0:1),5),
135  xinit(5), xx(5), ff(5), jj(5,5),
136  x_short(U.TMV_colsize()), f_short(U.TMV_colsize()),
137  Gx(bxsizep2,bxsize), Gy(bxsizep2,bxsize),
138  Gg1(bxsizep2,bxsize), Gg2(bxsizep2,bxsize),
139  Gth(bxsizep2,bxsize), Gmu(bxsizep2,bxsize)
140  {
141  //xdbg<<"EllipseSolver3: \n";
142  //xdbg<<"b0 = "<<b0.vec()<<std::endl;
143  //xdbg<<"b0.order, sigma = "<<b0.getOrder()<<" "<<b0.getSigma()<<std::endl;
144  //xdbg<<"order = "<<_order<<std::endl;
145  //xdbg<<"fixcen,gam,mu = "<<fixcen<<" "<<fixgam<<" "<<fixmu<<std::endl;
146  //xdbg<<"bxorder = "<<bxorder<<" "<<bxsize<<std::endl;
147  //xdbg<<"bxorderp2 = "<<bxorderp2<<" "<<bxsizep2<<std::endl;
148  U.setZero();
149  xinit.setZero();
150 
151  int k=0;
152  if (!fixcen) { U(k++,0) = 1.; U(k++,1) = 1.; }
153  if (!fixgam) { U(k++,2) = 1.; U(k++,3) = 1.; }
154  if (!fixmu) { U(k,4) = 1.; }
155 
156  Daug.setZero();
157  Saug.setZero();
158  Taug.setZero();
159 
160  setupGx(Gx,bxorderp2,bxorder);
161  setupGy(Gy,bxorderp2,bxorder);
162  setupGg1(Gg1,bxorderp2,bxorder);
163  setupGg2(Gg2,bxorderp2,bxorder);
164  setupGth(Gth,bxorderp2,bxorder);
165  setupGmu(Gmu,bxorderp2,bxorder);
166 
167  if (fixcen) {
168  T.TMV_diag().TMV_setAllTo(1.);
169  }
170  if (fixgam) {
171  S.TMV_diag().TMV_setAllTo(1.);
172  }
173  if (fixmu) {
174  D.TMV_diag().TMV_setAllTo(1.);
175  }
176  }
177 
178  void EllipseSolver3::ESImpl3::doF3(const DVector& x, DVector& f) const
179  {
180  Assert(x.size() == 5);
181  Assert(f.size() == 5);
182 
183  //xdbg<<"Start doF3\n";
184  //xdbg<<"x = "<<x<<std::endl;
185 
186 #ifdef MAX_X_DEV
187  if ((x-xinit).TMV_normInf() > MAX_X_DEV) {
188  dbg<<"bad x-xinit: "<<(x-xinit)<<std::endl;
189  xdbg<<"x = "<<x<<std::endl;
190  f = 2.* bx.TMV_subVector(1,6) / bx(0);
191  return;
192  }
193 #endif
194 
195  std::complex<double> zc(x[0],x[1]);
196  std::complex<double> g(x[2],x[3]);
197  double mu = x[4];
198 
199  if (norm(g) > 0.99) {
200  dbg<<"bad gsq = "<<norm(g)<<std::endl;
201  xdbg<<"x = "<<x<<std::endl;
202  f = 2.* bx.TMV_subVector(1,6) / bx(0);
203  return;
204  }
205 
206  bxsave = bx;
207  bool bxset = false;
208 
209  //xdbg<<"Start with b0 = "<<b0.vec()<<std::endl;
210  //xdbg<<"bx = "<<bx<<std::endl;
211  //TODO: Switch to doing calculate??Transform calls for only the portion
212  // of the matrix that I acutally use.
213  // This has implications for doJ though...
214  if (!fixmu) {
215  //xdbg<<"calc Mu:\n";
216  calculateMuTransform(mu,bxorder,Daug);
217  bx = TMV_colRange(D,0,b0size) * b0.vec();
218  bxset = true;
219  //xdbg<<"mu = "<<mu<<" bx => "<<bx<<std::endl;
220  } else if (fixm != 0.) {
221  bx = TMV_colRange(D,0,b0size) * b0.vec();
222  bxset = true;
223  }
224  if (!fixgam) {
225  //xdbg<<"calc G:\n";
226  calculateGTransform(g,bxorder,Saug);
227  if (bxset) {
228  bx = S * bx;
229  } else {
230  bx = TMV_colRange(S,0,b0size) * b0.vec();
231  bxset = true;
232  }
233  //xdbg<<"g = "<<g<<" bx => "<<bx<<std::endl;
234  } else if (fixg1 != 0. || fixg2 != 0.) {
235  if (bxset) {
236  bx = S * bx;
237  } else {
238  bx = TMV_colRange(S,0,b0size) * b0.vec();
239  bxset = true;
240  }
241  }
242  if (!fixcen) {
243  //xdbg<<"calc Z:\n";
244  calculateZTransform(zc,bxorder,Taug);
245  //xdbg<<"after calc Z\n";
246  if (bxset) {
247  //xdbg<<"T = "<<T<<std::endl;
248  bx = T * bx;
249  } else {
250  //xdbg<<"T = "<<TMV_colRange(T,0,b0size)<<std::endl;
251  bx = TMV_colRange(T,0,b0size) * b0.vec();
252  bxset = true;
253  }
254  //xdbg<<"zc = "<<zc<<" bx => "<<bx<<std::endl;
255  } else if (fixuc != 0. || fixvc != 0.) {
256  if (bxset) {
257  bx = T * bx;
258  } else {
259  bx = TMV_colRange(T,0,b0size) * b0.vec();
260  bxset = true;
261  }
262  }
263 
264  if (bx(0) <= 0.) {
265  dbg<<"bad bx(0): "<<bx(0)<<std::endl;
266  xdbg<<"x = "<<x<<std::endl;
267  xdbg<<"bx = "<<bx<<std::endl;
268  f = 2.* bxsave.TMV_subVector(1,6) / bxsave(0);
269  bx = bxsave;
270  return;
271  }
272 
273  f = bx.TMV_subVector(1,6) / bx(0);
274 
275  if (!zerob11) f(4) = 0.;
276  }
277 
278  void EllipseSolver3::ESImpl3::doF2(const DVector& x, DVector& f) const
279  {
280  Assert(x.size() == 5);
281  Assert(f.size() == 5);
282 
283  //xdbg<<"Start doF3\n";
284  //xdbg<<"x = "<<x<<std::endl;
285 
286 #ifdef MAX_X_DEV
287  if ((x-xinit).TMV_normInf() > MAX_X_DEV) {
288  dbg<<"bad x-xinit: "<<(x-xinit)<<std::endl;
289  xdbg<<"x = "<<x<<std::endl;
290  f = 2.* bx.TMV_subVector(1,6) / bx(0);
291  return;
292  }
293 #endif
294 
295  std::complex<double> zc(x[0],x[1]);
296  double mu = x[4];
297 
298  bxsave = bx;
299  bool bxset = false;
300 
301  //xdbg<<"Start with b0 = "<<b0.vec()<<std::endl;
302  //xdbg<<"bx = "<<bx<<std::endl;
303  if (!fixmu) {
304  //xdbg<<"calc Mu:\n";
305  calculateMuTransform(mu,bxorder,Daug);
306  bx = TMV_colRange(D,0,b0size) * b0.vec();
307  bxset = true;
308  //xdbg<<"mu = "<<mu<<" bx => "<<bx<<std::endl;
309  } else if (fixm != 0.) {
310  bx = TMV_colRange(D,0,b0size) * b0.vec();
311  bxset = true;
312  }
313  if (!fixcen) {
314  //xdbg<<"calc Z:\n";
315  calculateZTransform(zc,bxorder,Taug);
316  //xdbg<<"after calc Z\n";
317  if (bxset) {
318  //xdbg<<"T = "<<T<<std::endl;
319  bx = T * bx;
320  } else {
321  //xdbg<<"T = "<<TMV_colRange(T,0,b0size)<<std::endl;
322  bx = TMV_colRange(T,0,b0size) * b0.vec();
323  bxset = true;
324  }
325  //xdbg<<"zc = "<<zc<<" bx => "<<bx<<std::endl;
326  } else if (fixuc != 0. || fixvc != 0.) {
327  if (bxset) {
328  bx = T * bx;
329  } else {
330  bx = TMV_colRange(T,0,b0size) * b0.vec();
331  bxset = true;
332  }
333  }
334 
335  if (bx(0) <= 0.) {
336  dbg<<"bad bx(0): "<<bx(0)<<std::endl;
337  xdbg<<"x = "<<x<<std::endl;
338  xdbg<<"bx = "<<x<<std::endl;
339  f = 2.* bxsave.TMV_subVector(1,6) / bxsave(0);
340  bx = bxsave;
341  return;
342  }
343 
344  //xdbg<<"Done: bx = "<<bx<<std::endl;
345  f = bx.TMV_subVector(1,6) / bx(0);
346 
347  if (!zerob11) f(4) = 0.;
348  }
349 
350  void EllipseSolver3::ESImpl3::doF1(const DVector& x, DVector& f) const
351  {
352  Assert(x.size() == 5);
353  Assert(f.size() == 5);
354 
355  //xdbg<<"Start doF3\n";
356  //xdbg<<"x = "<<x<<std::endl;
357 
358 #ifdef MAX_X_DEV
359  if ((x-xinit).TMV_normInf() > MAX_X_DEV) {
360  dbg<<"bad x-xinit: "<<(x-xinit)<<std::endl;
361  xdbg<<"x = "<<x<<std::endl;
362  f = 2.* bx.TMV_subVector(1,6) / bx(0);
363  return;
364  }
365 #endif
366 
367  std::complex<double> zc(x[0],x[1]);
368 
369  bxsave = bx;
370 
371  //xdbg<<"Start with b0 = "<<b0.vec()<<std::endl;
372  //xdbg<<"bx = "<<bx<<std::endl;
373  if (!fixcen) {
374  //xdbg<<"calc Z:\n";
375  calculateZTransform(zc,bxorder,Taug);
376  //xdbg<<"after calc Z\n";
377  bx = TMV_colRange(T,0,b0size) * b0.vec();
378  //xdbg<<"zc = "<<zc<<" bx => "<<bx<<std::endl;
379  } else if (fixuc != 0. || fixvc != 0.) {
380  bx = TMV_colRange(T,0,b0size) * b0.vec();
381  }
382 
383  if (bx(0) <= 0.) {
384  dbg<<"bad bx(0): "<<bx(0)<<std::endl;
385  xdbg<<"x = "<<x<<std::endl;
386  xdbg<<"bx = "<<x<<std::endl;
387  f = 2.* bxsave.TMV_subVector(1,6) / bxsave(0);
388  bx = bxsave;
389  return;
390  }
391 
392  f = bx.TMV_subVector(1,6) / bx(0);
393 
394  if (!zerob11) f(4) = 0.;
395  }
396 
397  void EllipseSolver3::ESImpl3::doJ3(
398  const DVector& x, const DVector& f, DMatrix& J) const
399  {
400  // bx = T S D b0
401  //
402  // dbx/dzc = (dT/dzc) S D b0
403  // dbx/dg = T (dS/dg) D b0
404  // dbx/dmu = T S (dD/dmu) b0
405  //
406  // f = bx(1:6) / bx(0)
407  // df/dE = dbx/dE(1:6) / bx(0) - (1/bx(0)^2) (dbx(0)/dE) bx(1:6)
408  // = ( dbx/dE(1:6) - dbx/dE(0) f ) / bx(0)
409 
410  //xdbg<<"Start doJ3\n";
411  //xdbg<<"x = "<<x<<std::endl;
412  //xdbg<<"f = "<<f<<std::endl;
413 
414  Assert(x.size() == 5);
415  Assert(J.TMV_rowsize() == 5);
416  Assert(J.TMV_colsize() == 5);
417 
418  std::complex<double> zc(x[0],x[1]);
419  std::complex<double> g(x[2],x[3]);
420  double mu = x[4];
421  double fact = 1./(1.-norm(g));
422 
423  if (!fixcen || !fixgam) {
424  // Both cen and gam use this:
425  Db0 = TMV_colRange(D,0,b0size) * b0.vec();
426  }
427 
428  // Leave off one factor of bx(0) for now.
429  // We apply it at the end to the whole matrix.
430  if (!fixcen) {
431  // dT/dx = T Gx;
432  // dT/dy = T Gy;
433  // db/dx = T Gx S D b0
434  // db/dy = T Gy S D b0
435  augmentZTransformCols(zc,bxorder,Taug);
436  SDb0 = S * Db0;
437  GDb0 = Gx * SDb0;
438  dbdE = TMV_rowRange(Taug,0,6) * GDb0;
439  //dbg<<"dbdE = "<<dbdE<<std::endl;
440  J.col(0) = dbdE.TMV_subVector(1,6) - dbdE(0) * f;
441  //dbg<<"J(0) = "<<J.col(0)<<std::endl;
442  GDb0 = Gy * SDb0;
443  dbdE = TMV_rowRange(Taug,0,6) * GDb0;
444  //dbg<<"dbdE = "<<dbdE<<std::endl;
445  J.col(1) = dbdE.TMV_subVector(1,6) - dbdE(0) * f;
446  //dbg<<"J(1) = "<<J.col(1)<<std::endl;
447  } else {
448  J.col(0).setZero();
449  //dbg<<"J(0) = "<<J.col(0)<<std::endl;
450  J.col(1).setZero();
451  //dbg<<"J(1) = "<<J.col(1)<<std::endl;
452  }
453 
454  if (!fixgam) {
455  // dS/dg1 = fact * S * (Gg1 + g2 * Gth);
456  // dS/dg2 = fact * S * (Gg2 - g1 * Gth);
457  augmentGTransformCols(g,bxorder,Saug);
458  double g1 = real(g);
459  double g2 = imag(g);
460  GDb0 = Gg1 * Db0;
461  GthDb0 = Gth * Db0;
462  GDb0 += g2*GthDb0;
463  SDb0 = fact * Saug * GDb0;
464  dbdE = TMV_rowRange(T,0,6) * SDb0;
465  //dbg<<"dbdE = "<<dbdE<<std::endl;
466  J.col(2) = dbdE.TMV_subVector(1,6) - dbdE(0) * f;
467  //dbg<<"J(2) = "<<J.col(2)<<std::endl;
468  GDb0 = Gg2 * Db0;
469  GDb0 -= g1*GthDb0;
470  SDb0 = fact * Saug * GDb0;
471  dbdE = TMV_rowRange(T,0,6) * SDb0;
472  //dbg<<"dbdE = "<<dbdE<<std::endl;
473  J.col(3) = dbdE.TMV_subVector(1,6) - dbdE(0) * f;
474  //dbg<<"J(3) = "<<J.col(3)<<std::endl;
475  } else {
476  J.col(2).setZero();
477  //dbg<<"J(2) = "<<J.col(2)<<std::endl;
478  J.col(3).setZero();
479  //dbg<<"J(3) = "<<J.col(3)<<std::endl;
480  }
481 
482  if (!fixmu) {
483  // dD/dmu = D Gmu
484  augmentMuTransformCols(mu,bxorder,Daug);
485  GDb0 = TMV_colRange(Gmu,0,b0size) * b0.vec();
486  Db0 = Daug * GDb0;
487  SDb0 = S * Db0;
488  dbdE = TMV_rowRange(T,0,6) * SDb0;
489  //dbg<<"dbdE = "<<dbdE<<std::endl;
490  J.col(4) = dbdE.TMV_subVector(1,6) - dbdE(0) * f;
491  //dbg<<"J(4) = "<<J.col(4)<<std::endl;
492  } else {
493  J.col(4).setZero();
494  //dbg<<"J(4) = "<<J.col(4)<<std::endl;
495  }
496  J /= bx(0);
497  if (!zerob11) J.row(4).setZero();
498  //xdbg<<"J = "<<J<<std::endl;
499 
500 #ifdef JTEST
501  double dE = 1.e-6;
502  dbg<<"dE = "<<dE<<std::endl;
503  testout = dbgout;
504 
505  // This section was only used for debugging purposes, but
506  // I'm leaving it in, since the derivations of Gx[i] are
507  // helpful in understanding some of the code following this section.
508 
509  //
510  // dD/dmu:
511  //
512  DMatrix Dbig(bxsizep2,bxsizep2,0.);
513  DMatrix D0(bxsize,bxsize,0.);
514  DMatrix D1(bxsize,bxsize,0.);
515  DMatrix D2(bxsize,bxsize,0.);
516  calculateMuTransform(mu,bxorderp2,Dbig);
517  calculateMuTransform(mu,bxorder,D0);
518  calculateMuTransform(mu-dE,bxorder,D1);
519  calculateMuTransform(mu+dE,bxorder,D2);
520  DMatrix Gmubig(bxsizep2,bxsize,0.);
521  setupGmu(Gmubig,bxorderp2,bxorder);
522 
523  DMatrix dDdmu_num = (D2-D1)/(2.*dE);
524  DMatrix d2D_num = (D2+D1-2.*D0)/(dE*dE);
525  DMatrix dDdmu = TMV_rowRange(Dbig,0,bxsize) * Gmubig;
526  dbg<<"Gmu = "<<Gmubig.TMV_subMatrix(0,6,0,6)<<std::endl;
527  dbg<<"D = "<<Dbig.TMV_subMatrix(0,6,0,6)<<std::endl;
528  dbg<<"dDdmu = "<<dDdmu.TMV_subMatrix(0,6,0,6)<<std::endl;
529  dbg<<"dDdmu_num = "<<dDdmu_num.TMV_subMatrix(0,6,0,6)<<std::endl;
530  dbg<<"Norm(dD/dmu - numeric dD/dmu) = "<<Norm(dDdmu-dDdmu_num)<<std::endl;
531  dbg<<"dE*Norm(d2D_num) = "<<dE*Norm(d2D_num)<<std::endl;
532  test(Norm(dDdmu-dDdmu_num) < 30.*dE*Norm(d2D_num),"dDdmu");
533 
534  //
535  // dS/dg1
536  //
537  DMatrix Sbig(bxsizep2,bxsizep2,0.);
538  DMatrix S0(bxsize,bxsize,0.);
539  DMatrix S1(bxsize,bxsize,0.);
540  DMatrix S2(bxsize,bxsize,0.);
541  calculateGTransform(g,bxorderp2,Sbig);
542  calculateGTransform(g,bxorder,S0);
543  calculateGTransform(g-std::complex<double>(dE,0),bxorder,S1);
544  calculateGTransform(g+std::complex<double>(dE,0),bxorder,S2);
545  DMatrix Gg1big(bxsizep2,bxsize,0.);
546  DMatrix Gthbig(bxsizep2,bxsize,0.);
547  setupGg1(Gg1big,bxorderp2,bxorder);
548  setupGth(Gthbig,bxorderp2,bxorder);
549 
550  DMatrix dSdg1_num = (S2-S1)/(2.*dE);
551  DMatrix d2S_num = (S2+S1-2.*S0)/(dE*dE);
552  DMatrix dSdg1 = fact * TMV_rowRange(Sbig,0,bxsize) * (Gg1big + imag(g) * Gthbig);
553  dbg<<"Gg1 = "<<Gg1big.TMV_subMatrix(0,6,0,6)<<std::endl;
554  dbg<<"Gth = "<<Gthbig.TMV_subMatrix(0,6,0,6)<<std::endl;
555  dbg<<"S = "<<Sbig.TMV_subMatrix(0,6,0,6)<<std::endl;
556  dbg<<"dSdg1 = "<<dSdg1.TMV_subMatrix(0,6,0,6)<<std::endl;
557  dbg<<"dSdg1_num = "<<dSdg1_num.TMV_subMatrix(0,6,0,6)<<std::endl;
558  dbg<<"Norm(dS/dg1 - numeric dS/dg1) = "<<Norm(dSdg1-dSdg1_num)<<std::endl;
559  dbg<<"dE*Norm(d2S_num) = "<<dE*Norm(d2S_num)<<std::endl;
560  test(Norm(dSdg1-dSdg1_num) < 30.*dE*Norm(d2S_num),"dSdg1");
561 
562  //
563  // dS/dg2
564  //
565  calculateGTransform(g-std::complex<double>(0,dE),bxorder,S1);
566  calculateGTransform(g+std::complex<double>(0,dE),bxorder,S2);
567  DMatrix Gg2big(bxsizep2,bxsize,0.);
568  setupGg2(Gg2big,bxorderp2,bxorder);
569 
570  DMatrix dSdg2_num = (S2-S1)/(2.*dE);
571  d2S_num = (S2+S1-2.*S0)/(dE*dE);
572  DMatrix dSdg2 = fact * TMV_rowRange(Sbig,0,bxsize) * (Gg2big - real(g) * Gthbig);
573  dbg<<"Gg2 = "<<Gg2big.TMV_subMatrix(0,6,0,6)<<std::endl;
574  dbg<<"Gth = "<<Gthbig.TMV_subMatrix(0,6,0,6)<<std::endl;
575  dbg<<"S = "<<Sbig.TMV_subMatrix(0,6,0,6)<<std::endl;
576  dbg<<"dSdg2 = "<<dSdg2.TMV_subMatrix(0,6,0,6)<<std::endl;
577  dbg<<"dSdg2_num = "<<dSdg2_num.TMV_subMatrix(0,6,0,6)<<std::endl;
578  dbg<<"Norm(dS/dg2 - numeric dS/dg2) = "<<Norm(dSdg2-dSdg2_num)<<std::endl;
579  dbg<<"dE*Norm(d2S_num) = "<<dE*Norm(d2S_num)<<std::endl;
580  test(Norm(dSdg2-dSdg2_num) < 30.*dE*Norm(d2S_num),"dSdg2");
581 
582  //
583  // dT/dx
584  //
585  DMatrix Tbig(bxsizep2,bxsizep2,0.);
586  DMatrix T0(bxsize,bxsize,0.);
587  DMatrix T1(bxsize,bxsize,0.);
588  DMatrix T2(bxsize,bxsize,0.);
589  calculateZTransform(zc,bxorderp2,Tbig);
590  calculateZTransform(zc,bxorder,T0);
591  calculateZTransform(zc-std::complex<double>(dE,0),bxorder,T1);
592  calculateZTransform(zc+std::complex<double>(dE,0),bxorder,T2);
593  DMatrix Gxbig(bxsizep2,bxsize,0.);
594  setupGx(Gxbig,bxorderp2,bxorder);
595 
596  DMatrix dTdx_num = (T2-T1)/(2.*dE);
597  DMatrix d2T_num = (T2+T1-2.*T0)/(dE*dE);
598  DMatrix dTdx = TMV_rowRange(Tbig,0,bxsize) * Gxbig;
599  dbg<<"Gx = "<<Gxbig.TMV_subMatrix(0,6,0,6)<<std::endl;
600  dbg<<"T = "<<Tbig.TMV_subMatrix(0,6,0,6)<<std::endl;
601  dbg<<"dTdx = "<<dTdx.TMV_subMatrix(0,6,0,6)<<std::endl;
602  dbg<<"dTdx_num = "<<dTdx_num.TMV_subMatrix(0,6,0,6)<<std::endl;
603  dbg<<"Norm(dT/dx - numeric dT/dx) = "<<Norm(dTdx-dTdx_num)<<std::endl;
604  dbg<<"dE*Norm(d2T_num) = "<<dE*Norm(d2T_num)<<std::endl;
605  test(Norm(dTdx-dTdx_num) < 30.*dE*Norm(d2T_num),"dTdx");
606 
607  //
608  // dT/dy
609  //
610  calculateZTransform(zc-std::complex<double>(0,dE),bxorder,T1);
611  calculateZTransform(zc+std::complex<double>(0,dE),bxorder,T2);
612  DMatrix Gybig(bxsizep2,bxsize,0.);
613  setupGy(Gybig,bxorderp2,bxorder,0.);
614 
615  DMatrix dTdy_num = (T2-T1)/(2.*dE);
616  d2T_num = (T2+T1-2.*T0)/(dE*dE);
617  DMatrix dTdy = TMV_rowRange(Tbig,0,bxsize) * Gybig;
618  dbg<<"Gy = "<<Gybig.TMV_subMatrix(0,6,0,6)<<std::endl;
619  dbg<<"T = "<<Tbig.TMV_subMatrix(0,6,0,6)<<std::endl;
620  dbg<<"dTdy = "<<dTdy.TMV_subMatrix(0,6,0,6)<<std::endl;
621  dbg<<"dTdy_num = "<<dTdy_num.TMV_subMatrix(0,6,0,6)<<std::endl;
622  dbg<<"Norm(dT/dy - numeric dT/dy) = "<<Norm(dTdy-dTdy_num)<<std::endl;
623  dbg<<"dE*Norm(d2T_num) = "<<dE*Norm(d2T_num)<<std::endl;
624  test(Norm(dTdy-dTdy_num) < 30.*dE*Norm(d2T_num),"dTdy");
625 
626  //
627  // J
628  //
629 
630  DMatrix J_num(5,5,0.);
631  DVector x0 = x;
632  DVector f0(5);
633  doF3(x0,f0);
634  DVector bx0 = bx.TMV_subVector(0,6);
635  dbg<<"x0 = "<<x0<<std::endl;
636  dbg<<"b0 = "<<bx0<<std::endl;
637  dbg<<"f0 = "<<f0<<std::endl;
638 
639  for(int k=0;k<5;k++) {
640  dbg<<"k = "<<k<<std::endl;
641 
642  DVector f1(5);
643  DVector x1 = x; x1(k) -= dE;
644  doF3(x1,f1);
645  DVector b1 = bx.TMV_subVector(0,6);
646  dbg<<"x1 = "<<x1<<std::endl;
647  dbg<<"b1 = "<<b1<<std::endl;
648  dbg<<"f1 = "<<f1<<std::endl;
649 
650  DVector f2(5);
651  DVector x2 = x; x2(k) += dE;
652  doF3(x2,f2);
653  DVector b2 = bx.TMV_subVector(0,6);
654  dbg<<"x2 = "<<x2<<std::endl;
655  dbg<<"b2 = "<<b2<<std::endl;
656  dbg<<"f2 = "<<f2<<std::endl;
657 
658  DVector dbdE_num = (b2-b1)/(2.*dE);
659  dbg<<"dbdE_num = "<<dbdE_num<<std::endl;
660  DVector d2bdE = (b2+b1-2.*bx0)/(dE*dE);
661  dbg<<"d2bdE = "<<d2bdE<<std::endl;
662 
663  DVector dbdE(6);
664  double g1 = real(g);
665  double g2 = imag(g);
666  Db0 = TMV_colRange(D,0,b0size) * b0.vec();
667  dbg<<"Db0 = "<<Db0<<std::endl;
668  switch (k) {
669  case 0: SDb0 = S * Db0;
670  dbg<<"SDb0 = "<<SDb0<<std::endl;
671  GDb0 = Gx * SDb0;
672  dbg<<"GDb0 = "<<GDb0<<std::endl;
673  dbdE = TMV_rowRange(Tbig,0,6) * GDb0;
674  dbg<<"dbdx = "<<dbdE<<std::endl;
675  dbg<<"dbdx = dTdx S D b0 = "<<
676  TMV_rowRange(dTdx,0,6) * S * Db0<<std::endl;
677  break;
678  case 1: SDb0 = S * Db0;
679  dbg<<"SDb0 = "<<SDb0<<std::endl;
680  GDb0 = Gy * SDb0;
681  dbg<<"GDb0 = "<<GDb0<<std::endl;
682  dbdE = TMV_rowRange(Tbig,0,6) * GDb0;
683  dbg<<"dbdy = "<<dbdE<<std::endl;
684  dbg<<"dbdy = dTdy S D b0 = "<<
685  TMV_rowRange(dTdy,0,6) * S * Db0<<std::endl;
686  break;
687  case 2: GDb0 = Gg1 * Db0;
688  dbg<<"GDb0 = "<<GDb0<<std::endl;
689  GthDb0 = Gth * Db0;
690  dbg<<"GthDb0 = "<<GthDb0<<std::endl;
691  GDb0 += g2*GthDb0;
692  dbg<<"GDb0 = "<<GDb0<<std::endl;
693  SDb0 = fact * TMV_rowRange(Sbig,0,bxsize) * GDb0;
694  dbg<<"SDb0 = "<<SDb0<<std::endl;
695  dbdE = TMV_rowRange(T,0,6) * SDb0;
696  dbg<<"dbdg1 = "<<dbdE<<std::endl;
697  dbg<<"dbdg1 = T dSdg1 D b0 = "<<
698  TMV_rowRange(T,0,6) * dSdg1 * Db0<<std::endl;
699  break;
700  case 3: GDb0 = Gg2 * Db0;
701  dbg<<"GDb0 = "<<GDb0<<std::endl;
702  GthDb0 = Gth * Db0;
703  dbg<<"GthDb0 = "<<GthDb0<<std::endl;
704  GDb0 -= g1*GthDb0;
705  dbg<<"GDb0 = "<<GDb0<<std::endl;
706  SDb0 = fact * TMV_rowRange(Sbig,0,bxsize) * GDb0;
707  dbg<<"SDb0 = "<<SDb0<<std::endl;
708  dbdE = TMV_rowRange(T,0,6) * SDb0;
709  dbg<<"dbdg2 = "<<dbdE<<std::endl;
710  dbg<<"dbdg2 = T dSdg2 D b0 = "<<
711  TMV_rowRange(T,0,6) * dSdg2 * Db0<<std::endl;
712  break;
713  case 4: GDb0 = TMV_colRange(Gmu,0,b0size) * b0.vec();
714  dbg<<"GDb0 = "<<GDb0<<std::endl;
715  Db0 = TMV_rowRange(Dbig,0,bxsize) * GDb0;
716  dbg<<"Db0 = "<<Db0<<std::endl;
717  SDb0 = S * Db0;
718  dbg<<"SDb0 = "<<SDb0<<std::endl;
719  dbdE = TMV_rowRange(T,0,6) * SDb0;
720  dbg<<"dbdmu = "<<dbdE<<std::endl;
721  dbg<<"dbdmu = T S dDdmu b0 = "<<
722  TMV_rowRange(T,0,6) * S * TMV_colRange(dDdmu,0,b0size) * b0.vec()
723  <<std::endl;
724  break;
725  }
726  dbg<<"dbdE = "<<dbdE<<std::endl;
727  dbg<<"Norm(dbdE-dbdE_num) = "<<Norm(dbdE-dbdE_num)<<std::endl;
728  dbg<<"dE*Norm(d2bdE) = "<<dE*Norm(d2bdE)<<std::endl;
729  test(Norm(dbdE-dbdE_num) < 30.*dE*Norm(d2bdE),"dbdE");
730 
731  DVector dfdE = dbdE.TMV_subVector(1,6) - dbdE(0) * f0;
732  dfdE /= bx(0);
733  dbg<<"dfdE from dbdE = "<<dfdE<<std::endl;
734  dbg<<"dfdE from dbdE_num = "<<
735  (dbdE_num.TMV_subVector(1,6) - dbdE_num(0) * f0)/bx(0)<<std::endl;;
736  dbg<<"dfdE in J = "<<J.col(k)<<std::endl;
737  J_num.col(k) = (f2-f1)/(2.*dE);
738  dbg<<"dfdE_num = "<<J_num.col(k)<<std::endl;
739  DVector d2f_num = (f2+f1-2.*f0)/(dE*dE);
740  dbg<<"d2f_num = "<<d2f_num<<std::endl;
741  dbg<<"Norm(dfdE_num-dfdE) ="<<Norm(J_num.col(k)-dfdE)<<std::endl;
742  dbg<<"dE*Norm(d2f_num) ="<<dE*Norm(d2f_num)<<std::endl;
743  test(Norm(J_num.col(k)-dfdE) < 30.*dE*Norm(d2f_num),"dfdE");
744  }
745 
746  dbg<<"J_num = "<<J_num<<std::endl;
747  dbg<<"analytic: "<<J<<std::endl;
748  if (fixcen) TMV_colRange(J_num,0,2).setZero();
749  if (fixgam) TMV_colRange(J_num,2,4).setZero();
750  if (fixmu) J_num.col(4).setZero();
751  dbg<<"J_num => "<<J_num<<std::endl;
752  dbg<<"Norm(diff) = "<<Norm(J-J_num)<<std::endl;
753  dbg<<"dE*Norm(J) = "<<dE*Norm(J_num)<<std::endl;
754  test(Norm(J_num-J) < dE*Norm(J_num),"dfdE");
755  doF3(x,f0); // return everything to original values
756 #endif
757  }
758 
759  void EllipseSolver3::ESImpl3::doJ2(
760  const DVector& x, const DVector& f, DMatrix& J) const
761  {
762  //xdbg<<"Start doJ2\n";
763  //xdbg<<"x = "<<x<<std::endl;
764  //xdbg<<"f = "<<f<<std::endl;
765 
766  Assert(x.size() == 5);
767  Assert(J.TMV_rowsize() == 5);
768  Assert(J.TMV_colsize() == 5);
769 
770  std::complex<double> zc(x[0],x[1]);
771  double mu = x[4];
772 
773  J.setZero();
774  if (!fixcen) {
775  augmentZTransformCols(zc,bxorder,Taug);
776  Db0 = TMV_colRange(D,0,b0size) * b0.vec();
777  //dbg<<"Db0 = "<<Db0<<std::endl;
778  GDb0 = Gx * Db0;
779  dbdE = TMV_rowRange(Taug,0,6) * GDb0;
780  //dbg<<"dbdE = "<<dbdE<<std::endl;
781  J.col(0) = dbdE.TMV_subVector(1,6) - dbdE(0) * f;
782  //dbg<<"J(0) = "<<J.col(0)<<std::endl;
783  GDb0 = Gy * Db0;
784  dbdE = TMV_rowRange(Taug,0,6) * GDb0;
785  //dbg<<"dbdE = "<<dbdE<<std::endl;
786  J.col(1) = dbdE.TMV_subVector(1,6) - dbdE(0) * f;
787  //dbg<<"J(1) = "<<J.col(1)<<std::endl;
788  } else {
789  //dbg<<"J(0) = "<<J.col(0)<<std::endl;
790  //dbg<<"J(1) = "<<J.col(1)<<std::endl;
791  }
792 
793  if (!fixmu) {
794  // dD/dmu = D Gmu
795  augmentMuTransformCols(mu,bxorder,Daug);
796  GDb0 = TMV_colRange(Gmu,0,b0size) * b0.vec();
797  //dbg<<"GDb0 = "<<GDb0<<std::endl;
798  Db0 = Daug * GDb0;
799  //dbg<<"Db0 = "<<Db0<<std::endl;
800  SDb0 = S * Db0;
801  //dbg<<"SDb0 = "<<SDb0<<std::endl;
802  dbdE = TMV_rowRange(T,0,6) * SDb0;
803  //dbg<<"dbdE = "<<dbdE<<std::endl;
804  J.col(4) = dbdE.TMV_subVector(1,6) - dbdE(0) * f;
805  //dbg<<"J(4) = "<<J.col(4)<<std::endl;
806  } else {
807  //dbg<<"J(4) = "<<J.col(4)<<std::endl;
808  }
809  //xdbg<<"bx(0) = "<<bx(0)<<std::endl;
810  J /= bx(0);
811  if (!zerob11) J.row(4).setZero();
812  //xdbg<<"J = "<<J<<std::endl;
813  }
814 
815  void EllipseSolver3::ESImpl3::doJ1(
816  const DVector& x, const DVector& f, DMatrix& J) const
817  {
818  //xdbg<<"Start doJ1\n";
819  //xdbg<<"x = "<<x<<std::endl;
820  //xdbg<<"f = "<<f<<std::endl;
821 
822  Assert(x.size() == 5);
823  Assert(J.TMV_rowsize() == 5);
824  Assert(J.TMV_colsize() == 5);
825 
826  std::complex<double> zc(x[0],x[1]);
827 
828  J.setZero();
829 
830  augmentZTransformCols(zc,bxorder,Taug);
831  GDb0 = TMV_colRange(Gx,0,b0size) * b0.vec();
832  dbdE = TMV_rowRange(Taug,0,6) * GDb0;
833  //dbg<<"dbdE = "<<dbdE<<std::endl;
834  J.col(0) = dbdE.TMV_subVector(1,6) - dbdE(0) * f;
835  //dbg<<"J(0) = "<<J.col(0)<<std::endl;
836  GDb0 = TMV_colRange(Gy,0,b0size) * b0.vec();
837  dbdE = TMV_rowRange(Taug,0,6) * GDb0;
838  //dbg<<"dbdE = "<<dbdE<<std::endl;
839  J.col(1) = dbdE.TMV_subVector(1,6) - dbdE(0) * f;
840  //dbg<<"J(1) = "<<J.col(1)<<std::endl;
841 
842  J /= bx(0);
843  if (!zerob11) J.row(4).setZero();
844  //xdbg<<"J = "<<J<<std::endl;
845  }
846 
847  void EllipseSolver3::ESImpl3::calculateF(const DVector& x, DVector& f) const
848  {
849  Assert(x.size() == U.TMV_colsize());
850  Assert(f.size() == U.TMV_colsize());
851 
852  Assert(xx.size() == U.TMV_rowsize());
853  EIGEN_Transpose(xx) = EIGEN_Transpose(x)*U;
854  if (fixcen) { xx[0] = fixuc; xx[1] = fixvc; }
855  if (fixgam) { xx[2] = fixg1; xx[3] = fixg2; }
856  if (fixmu) { xx[4] = fixm; }
857 
858  if (!fixgam || fixg1 != 0. || fixg2 != 0.)
859  doF3(xx,ff);
860  else if (!fixmu || fixm != 0.)
861  doF2(xx,ff);
862  else
863  doF1(xx,ff);
864 
865  f = U*ff;
866  }
867 
868  void EllipseSolver3::ESImpl3::calculateJ(
869  const DVector& x, const DVector& f, DMatrix& j) const
870  {
871  Assert(x.size() == U.TMV_colsize());
872  Assert(f.size() == U.TMV_colsize());
873  Assert(j.TMV_rowsize() == U.TMV_colsize());
874  Assert(j.TMV_colsize() == U.TMV_colsize());
875 
876  EIGEN_Transpose(xx) = EIGEN_Transpose(x)*U;
877  if (fixcen) { xx[0] = fixuc; xx[1] = fixvc; }
878  if (fixgam) { xx[2] = fixg1; xx[3] = fixg2; }
879  if (fixmu) { xx[4] = fixm; }
880 
881  if (!fixgam || fixg1 != 0. || fixg2 != 0.)
882  doJ3(xx,ff,jj);
883  else if (!fixmu || fixm != 0.)
884  doJ2(xx,ff,jj);
885  else
886  doJ1(xx,ff,jj);
887 
888  j = U*jj*U.transpose();
889  }
890 
891  void EllipseSolver3::useNumericJ() { _pimpl->numeric_j = true; }
892 
893  void EllipseSolver3::dontZeroB11() { _pimpl->zerob11 = false; this->useSVD(); }
894 
895  void EllipseSolver3::getCovariance(DMatrix& cov) const
896  {
897  DMatrix cov1(_pimpl->x_short.size(),_pimpl->x_short.size());
898  NLSolver::getCovariance(cov1);
899  cov = _pimpl->U.transpose()*cov1*_pimpl->U;
900  dbg<<"getCovariance:\n";
901  dbg<<"cov1 = "<<cov1<<std::endl;
902  dbg<<"full cov = "<<cov<<std::endl;
903  }
904 
905  void EllipseSolver3::getInverseCovariance(DMatrix& invcov) const
906  {
907  DMatrix invcov1(_pimpl->x_short.size(),_pimpl->x_short.size());
908  NLSolver::getInverseCovariance(invcov1);
909  invcov = _pimpl->U.transpose()*invcov1*_pimpl->U;
910  dbg<<"getInverseCovariance:\n";
911  dbg<<"invcov1 = "<<invcov1<<std::endl;
912  dbg<<"full invcov = "<<invcov<<std::endl;
913  }
914 
915  void EllipseSolver3::callF(const DVector& x, DVector& f) const
916  {
917  Assert(x.size() == 5);
918  Assert(f.size() == 5);
919  _pimpl->xinit = x;
920  if (_pimpl->fixcen) { _pimpl->fixuc = x[0]; _pimpl->fixvc = x[1]; }
921  if (_pimpl->fixgam) { _pimpl->fixg1 = x[2]; _pimpl->fixg2 = x[3]; }
922  if (_pimpl->fixmu) { _pimpl->fixm = x[4]; }
923 
924  _pimpl->x_short = _pimpl->U * x;
925 
926  calculateF(_pimpl->x_short,_pimpl->f_short);
927 
928  EIGEN_Transpose(f) = EIGEN_Transpose(_pimpl->f_short)*_pimpl->U;
929  }
930 
931  bool EllipseSolver3::solve(DVector& x, DVector& f) const
932  {
933  Assert(x.size() == 5);
934  Assert(f.size() == 5);
935  _pimpl->xinit = x;
936  _pimpl->bx.TMV_setAllTo(1.); // In case we bail out with f = 2*bx/b(0)
937  if (_pimpl->fixcen) { _pimpl->fixuc = x[0]; _pimpl->fixvc = x[1]; }
938  if (_pimpl->fixgam) { _pimpl->fixg1 = x[2]; _pimpl->fixg2 = x[3]; }
939  if (_pimpl->fixmu) { _pimpl->fixm = x[4]; }
940  if (_pimpl->fixcen && (x[0] != 0. || x[1] != 0.)) {
941  calculateZTransform(
942  std::complex<double>(x[0],x[1]),_pimpl->bxorder,_pimpl->Taug);
943  }
944  if (_pimpl->fixgam && (x[2] != 0. || x[3] != 0.)) {
945  calculateGTransform(
946  std::complex<double>(x[2],x[3]),_pimpl->bxorder,_pimpl->Saug);
947  }
948  if (_pimpl->fixmu && x[4] != 0.) {
949  calculateMuTransform(x[4],_pimpl->bxorder,_pimpl->Daug);
950  }
951 
952  _pimpl->x_short = _pimpl->U * x;
953 
954  bool ret;
955  try {
956  ret = NLSolver::solve(_pimpl->x_short,_pimpl->f_short);
957  } catch (...) {
958  xdbg<<"Caught exception during NLSolver::solve"<<std::endl;
959  ret = false;
960  }
961 
962  EIGEN_Transpose(x) = EIGEN_Transpose(_pimpl->x_short)*_pimpl->U;
963  if (_pimpl->fixcen) { x[0] = _pimpl->fixuc; x[1] = _pimpl->fixvc; }
964  if (_pimpl->fixgam) { x[2] = _pimpl->fixg1; x[3] = _pimpl->fixg2; }
965  if (_pimpl->fixmu) { x[4] = _pimpl->fixm; }
966  EIGEN_Transpose(f) = EIGEN_Transpose(_pimpl->f_short)*_pimpl->U;
967 
968  return ret;
969  }
970 
971  bool EllipseSolver3::testJ(
972  const DVector& x, DVector& f,
973  std::ostream* os, double relerr) const
974  {
975  Assert(x.size() == 5);
976  Assert(f.size() == 5);
977  _pimpl->xinit = x;
978  if (_pimpl->fixcen) { _pimpl->fixuc = x[0]; _pimpl->fixvc = x[1]; }
979  if (_pimpl->fixgam) { _pimpl->fixg1 = x[2]; _pimpl->fixg2 = x[3]; }
980  if (_pimpl->fixmu) { _pimpl->fixm = x[4]; }
981 
982  _pimpl->x_short = _pimpl->U * x;
983 
984  return NLSolver::testJ(_pimpl->x_short,_pimpl->f_short,os,relerr);
985  }
986 
987 
988 }}}}
lsst::meas::algorithms::shapelet::NLSolver::getCovariance
virtual void getCovariance(DMatrix &cov) const
Definition: NLSolver.cc:1638
MAX_X_DEV
#define MAX_X_DEV
Definition: EllipseSolver.cc:30
lsst::meas::algorithms::shapelet::EllipseSolver3::ESImpl3::calculateJ
void calculateJ(const DVector &x, const DVector &f, DMatrix &J) const
Definition: EllipseSolver.cc:868
lsst::meas::algorithms::shapelet::augmentMuTransformCols
void augmentMuTransformCols(double mu, int order, DMatrix &D)
Definition: BVec.cc:316
lsst::meas::algorithms::shapelet::setupGg1
void setupGg1(DMatrix &Gg1, int order1, int order2)
Definition: PsiHelper.cc:449
lsst::meas::algorithms::shapelet::EllipseSolver3::ESImpl3::ESImpl3
ESImpl3(const BVec &b0, int order, bool isFixedCen, bool isFixedGamma, bool isFixedMu)
Definition: EllipseSolver.cc:120
lsst::meas::algorithms::shapelet::setupGth
void setupGth(DMatrix &Gth, int order1, int order2)
Definition: PsiHelper.cc:579
lsst::meas::algorithms::shapelet::NLSolver::solve
virtual bool solve(DVector &x, DVector &f) const
Definition: NLSolver.cc:1620
lsst.afw.math::details::norm
T norm(const T &x)
Definition: Integrate.h:191
lsst::meas::algorithms::shapelet::EllipseSolver3::getInverseCovariance
void getInverseCovariance(DMatrix &invcov) const
Definition: EllipseSolver.cc:905
lsst::meas::algorithms::shapelet::calculateGTransform
void calculateGTransform(std::complex< double > g, int order1, int order2, DMatrix &S)
Definition: BVec.cc:468
lsst::meas::algorithms::shapelet::EllipseSolver3::solve
bool solve(DVector &x, DVector &f) const
Definition: EllipseSolver.cc:931
lsst.afw.display::x0
int const x0
Definition: saturated.cc:45
TMV_colsize
#define TMV_colsize()
Definition: MyMatrix.h:175
lsst::meas::algorithms::shapelet::augmentZTransformCols
void augmentZTransformCols(std::complex< double > z, int order, DMatrix &T)
Definition: BVec.cc:156
lsst::meas::algorithms::shapelet::calculateZTransform
void calculateZTransform(std::complex< double > z, int order1, int order2, DMatrix &T)
Definition: BVec.cc:84
lsst::meas::algorithms::shapelet::NLSolver::calculateJ
virtual void calculateJ(const DVector &x, const DVector &f, DMatrix &j) const
Definition: NLSolver.cc:1201
TMV_colRange
#define TMV_colRange(m, j1, j2)
Definition: MyMatrix.h:203
lsst::meas::algorithms::shapelet::EllipseSolver3::ESImpl3::doF1
void doF1(const DVector &x, DVector &f) const
Definition: EllipseSolver.cc:350
PsiHelper.h
lsst::meas::algorithms::shapelet::EllipseSolver3::calculateJ
void calculateJ(const DVector &x, const DVector &f, DMatrix &df) const
Definition: EllipseSolver.cc:109
lsst::meas::algorithms::shapelet::EllipseSolver3::ESImpl3::doF2
void doF2(const DVector &x, DVector &f) const
Definition: EllipseSolver.cc:278
lsst::meas::algorithms::shapelet::setupGx
void setupGx(DMatrix &Gx, int order1, int order2)
Definition: PsiHelper.cc:370
lsst::meas::algorithms::shapelet::DMatrixView
Eigen::Block< DMatrix > DMatrixView
Definition: MyMatrix.h:133
EllipseSolver.h
lsst::meas::algorithms::shapelet::setupGy
void setupGy(DMatrix &Gy, int order1, int order2)
Definition: PsiHelper.cc:409
lsst::meas::algorithms::shapelet::augmentGTransformCols
void augmentGTransformCols(std::complex< double > g, int order, DMatrix &S)
Definition: BVec.cc:577
x
int x
Definition: GaussianCentroid.cc:36
lsst::meas::algorithms::shapelet::EllipseSolver3::ESImpl3::calculateF
void calculateF(const DVector &x, DVector &f) const
Definition: EllipseSolver.cc:847
lsst::meas::algorithms::shapelet::setupGg2
void setupGg2(DMatrix &Gg2, int order1, int order2)
Definition: PsiHelper.cc:500
lsst::meas::algorithms::shapelet::EllipseSolver3::EllipseSolver3
EllipseSolver3(const BVec &b0, int order, bool fixcen=false, bool fixgam=false, bool fixmu=false)
Definition: EllipseSolver.cc:97
lsst::meas::algorithms::shapelet::NLSolver::testJ
virtual bool testJ(const DVector &, DVector &, std::ostream *os=0, double relerr=0.) const
Definition: NLSolver.cc:1224
lsst::meas::algorithms::shapelet::NLSolver::getInverseCovariance
virtual void getInverseCovariance(DMatrix &invcov) const
Definition: NLSolver.cc:1675
lsst::meas::algorithms::shapelet::EllipseSolver3::ESImpl3::doJ1
void doJ1(const DVector &x, const DVector &f, DMatrix &J) const
Definition: EllipseSolver.cc:815
TMV_rowRange
#define TMV_rowRange(m, i1, i2)
Definition: MyMatrix.h:204
lsst::meas::algorithms::shapelet::EllipseSolver3::testJ
bool testJ(const DVector &x, DVector &f, std::ostream *os=0, double relerr=0.) const
Definition: EllipseSolver.cc:971
xdbg
#define xdbg
Definition: dbg.h:70
lsst::meas::algorithms::shapelet::EllipseSolver3::calculateF
void calculateF(const DVector &x, DVector &f) const
Definition: EllipseSolver.cc:106
dbg
#define dbg
Definition: dbg.h:69
lsst.afw.math::details::real
T real(const T &x)
Definition: Integrate.h:193
lsst::meas::algorithms::shapelet::EllipseSolver3::ESImpl3::doJ2
void doJ2(const DVector &x, const DVector &f, DMatrix &J) const
Definition: EllipseSolver.cc:759
dbgout
std::ostream *const dbgout
Definition: dbg.h:64
EIGEN_Transpose
#define EIGEN_Transpose(m)
Definition: MyMatrix.h:211
lsst::meas::algorithms::shapelet::EllipseSolver3::callF
void callF(const DVector &x, DVector &f) const
Definition: EllipseSolver.cc:915
Assert
#define Assert(x)
Definition: dbg.h:73
lsst::meas::algorithms::shapelet::EllipseSolver3::ESImpl3::doJ3
void doJ3(const DVector &x, const DVector &f, DMatrix &J) const
Definition: EllipseSolver.cc:397
lsst::meas::algorithms::shapelet::EllipseSolver3::ESImpl3::doF3
void doF3(const DVector &x, DVector &f) const
Definition: EllipseSolver.cc:178
lsst::meas::algorithms::shapelet::calculateMuTransform
void calculateMuTransform(double mu, int order1, int order2, DMatrix &D)
Definition: BVec.cc:227
lsst::meas::algorithms::shapelet::setupGmu
void setupGmu(DMatrix &Gmu, int order1, int order2)
Definition: PsiHelper.cc:551
Generated on Wed Sep 16 2015 13:35:31 for LSSTApplications by   doxygen 1.8.5