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
PsiHelper.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 <vector>
26 
29 #include "lsst/afw/geom/Angle.h"
30 
31 namespace afwGeom = lsst::afw::geom;
32 
33 namespace lsst {
34 namespace meas {
35 namespace algorithms {
36 namespace shapelet {
37 
38  // Here is the order of p,q along the indices of psi:
39  //
40  // k 0 1,2 3,4 5 6,7 8,9 10,11 12,13 14 15,16 17,18 19,20 21,22 23,24 25,26 27
41  // p 0 1 2 1 3 2 4 3 2 5 4 3 6 5 4 3
42  // q 0 0 0 1 0 1 0 1 2 0 1 2 0 1 2 3
43  // n 0 1 2 2 3 3 4 4 4 5 5 5 6 6 6 6
44  // m 0 1 2 0 3 1 4 2 0 5 3 1 6 4 2 0
45  //
46 
47 #ifdef USE_TMV
48  // The Eigen version required pretty significant rewriting before it was
49  // remotely efficient, so I separate this out with an ifdef, rather than
50  // use the TMV and EIGEN macros.
51 
52  void makePsi(DMatrix& psi, CDVectorView z, int order, const DVector* coeff)
53  {
54  // For p>=q:
55  //
56  // psi_pq = (pi p! q!)^-1/2 z^m exp(-r^2/2) K_pq(r^2)
57  //
58  // where K_pq(r^2) satisfies the recursion relation:
59  // K_pq = (r^2 - (N-1)) K_p-1,q-1 - (p-1)(q-1) K_p-2,q-2
60  // with K_p0 = 1
61  //
62  // The recursion for psi_pq can then be derived as:
63  //
64  // psi_pq = (pq)^-1/2 (r^2-(N-1)) psi_p-1,q-1
65  // - sqrt( (p-1)(q-1)/(pq) ) psi_p-2,q-2
66  //
67  // psi_p0 = (z/sqrt(p)) psi_p-1,0
68  //
69  // psi_00 = 1/sqrt(pi) exp(-r^2/2)
70 
71  Assert(int(psi.rowsize()) >= (order+1)*(order+2)/2);
72  Assert(psi.colsize() == z.size());
73  if (coeff) Assert(psi.colsize() == coeff->size());
74  Assert(psi.iscm());
75  Assert(!psi.isconj());
76 
77  // Setup rsq, z vectors and set psi_00
78  DVector rsq(z.size());
79  double* rsqit = rsq.ptr();
80  double* psi00it = psi.ptr();
81  const std::complex<double>* zit = z.cptr();
82  const int zsize = z.size();
83  for(int i=0;i<zsize;++i) {
84  rsqit[i] = std::norm(zit[i]);
85  psi00it[i] = afwGeom::INVSQRTPI * exp(-(rsqit[i])/2.);
86  }
87  if (coeff) psi.col(0) *= DiagMatrixViewOf(*coeff);
88 
89  DVector zr = z.realPart();
90  DVector zi = z.imagPart();
91  if (order >= 1) {
92  // Set psi_10
93  // All m > 0 elements are intrinsically complex.
94  // However, we are fitting to a real intensity pattern
95  // with complex coefficients of the complex shapelets.
96  // Since psi_pq = psi_qp* (* = complex conjugate),
97  // we know that b_pq must be b_qp*
98  // b_pq psi_pq + b_pq* psi_pq* = 2 b_pqr psi_pqr - 2 b_pqi psi_pqi
99  // So the values we want for the real fitter are
100  // really 2 real(psi_pq) and -2 imag(psi_pq)
101  // Putting the 2's here carries through to the rest of the
102  // elements via the recursion.
103  psi.col(1) = 2. * DiagMatrixViewOf(zr) * psi.col(0);
104  psi.col(2) = -2. * DiagMatrixViewOf(zi) * psi.col(0);
105  }
106  for(int N=2,k=3;N<=order;++N) {
107  // Set psi_N0
108  // The signs of these are not what you naively think due to
109  // the +2, -2 discussed above. You just have to follow through
110  // what the complex psi_N0 is, and what value is stored in the
111  // psi_N-1,0 location, and what needs to get stored here.
112  psi.col(k) = sqrt(1./N) * DiagMatrixViewOf(zr) * psi.col(k-N);
113  psi.col(k) += sqrt(1./N) * DiagMatrixViewOf(zi) * psi.col(k-N+1);
114  psi.col(k+1) = -sqrt(1./N) * DiagMatrixViewOf(zi) * psi.col(k-N);
115  psi.col(k+1) += sqrt(1./N) * DiagMatrixViewOf(zr) * psi.col(k-N+1);
116  k+=2;
117 
118  // Set psi_pq with q>0
119  // The rsq part of this calculation can be done in batch, which
120  // speeds things up a bit.
121  psi.colRange(k,k+N-1) =
122  DiagMatrixViewOf(rsq) * psi.colRange(k-2*N-1,k-N-2);
123  psi.colRange(k,k+N-1) -= (N-1.) * psi.colRange(k-2*N-1,k-N-2);
124  // The other calculation steps are different for each component:
125  for(int m=N-2,p=N-1,q=1;m>=0;--p,++q,m-=2) {
126  double pq = p*q;
127  if (m==0) {
128  psi.col(k) /= sqrt(pq);
129  if (q > 1) psi.col(k) -= sqrt(1.-(N-1.)/pq)*psi.col(k+2-4*N);
130  ++k;
131  } else {
132  psi.colRange(k,k+2) /= sqrt(pq);
133  if (q > 1)
134  psi.colRange(k,k+2) -=
135  sqrt(1.-(N-1.)/pq)*psi.colRange(k+2-4*N,k+4-4*N);
136  k+=2;
137  }
138  }
139  }
140  }
141 
142  void makePsi(DVector& psi, std::complex<double> z, int order)
143  {
144  double rsq = std::norm(z);
145  psi(0) = afwGeom::INVSQRTPI * exp(-rsq/2.);
146 
147  double zr = std::real(z);
148  double zi = std::imag(z);
149  if (order >= 1) {
150  psi(1) = 2. * zr * psi(0);
151  psi(2) = -2. * zi * psi(0);
152  }
153  for(int N=2,k=3;N<=order;++N) {
154  psi(k) = sqrt(1./N) * zr * psi(k-N);
155  psi(k) += sqrt(1./N) * zi * psi(k-N+1);
156  psi(k+1) = -sqrt(1./N) * zi * psi(k-N);
157  psi(k+1) += sqrt(1./N) * zr * psi(k-N+1);
158  k+=2;
159 
160  psi.subVector(k,k+N-1) = rsq * psi.subVector(k-2*N-1,k-N-2);
161  psi.subVector(k,k+N-1) -= (N-1.) * psi.subVector(k-2*N-1,k-N-2);
162  for(int m=N-2,p=N-1,q=1;m>=0;--p,++q,m-=2) {
163  double pq = p*q;
164  if (m==0) {
165  psi(k) /= sqrt(pq);
166  if (q > 1) psi(k) -= sqrt(1.-(N-1.)/pq)*psi(k+2-4*N);
167  ++k;
168  } else {
169  psi.subVector(k,k+2) /= sqrt(pq);
170  if (q > 1)
171  psi.subVector(k,k+2) -=
172  sqrt(1.-(N-1.)/pq)*psi.subVector(k+2-4*N,k+4-4*N);
173  k+=2;
174  }
175  }
176  }
177  }
178 
179  void augmentPsi(DMatrix& psi, CDVectorView z, int order)
180  {
181  Assert(int(psi.rowsize()) >= (order+3)*(order+4)/2);
182  Assert(psi.colsize() == z.size());
183  Assert(order >= 1);
184  //Assert(psi.iscm());
185  //Assert(!psi.isconj());
186 
187  DVector rsq(z.size());
188  double* rsqit = rsq.ptr();
189  const std::complex<double>* zit = z.cptr();
190  const int zsize = z.size();
191  for(int i=0;i<zsize;++i) {
192  rsqit[i] = std::norm(zit[i]);
193  }
194 
195  DVector zr = z.realPart();
196  DVector zi = z.imagPart();
197  for(int N=order+1,k=N*(N+1)/2;N<=order+2;++N) {
198  psi.col(k) = sqrt(1./N) * DiagMatrixViewOf(zr) * psi.col(k-N);
199  psi.col(k) += sqrt(1./N) * DiagMatrixViewOf(zi) * psi.col(k-N+1);
200  psi.col(k+1) = -sqrt(1./N) * DiagMatrixViewOf(zi) * psi.col(k-N);
201  psi.col(k+1) += sqrt(1./N) * DiagMatrixViewOf(zr) * psi.col(k-N+1);
202  k+=2;
203 
204  psi.colRange(k,k+N-1) =
205  DiagMatrixViewOf(rsq) * psi.colRange(k-2*N-1,k-N-2);
206  psi.colRange(k,k+N-1) -= (N-1.) * psi.colRange(k-2*N-1,k-N-2);
207  for(int m=N-2,p=N-1,q=1;m>=0;--p,++q,m-=2) {
208  double pq = p*q;
209  if (m==0) {
210  psi.col(k) /= sqrt(pq);
211  if (q > 1) psi.col(k) -= sqrt(1.-(N-1.)/pq)*psi.col(k+2-4*N);
212  ++k;
213  } else {
214  psi.colRange(k,k+2) /= sqrt(pq);
215  if (q > 1)
216  psi.colRange(k,k+2) -=
217  sqrt(1.-(N-1.)/pq)*psi.colRange(k+2-4*N,k+4-4*N);
218  k+=2;
219  }
220  }
221  }
222  }
223 
224 #else
225 
226  void makePsi(DMatrix& psi, CDVectorView z, int order, const DVector* coeff)
227  {
228  // For p>=q:
229  //
230  // psi_pq = (pi p! q!)^-1/2 z^m exp(-r^2/2) K_pq(r^2)
231  //
232  // where K_pq(r^2) satisfies the recursion relation:
233  // K_pq = (r^2 - (N-1)) K_p-1,q-1 - (p-1)(q-1) K_p-2,q-2
234  // with K_p0 = 1
235  //
236  // The recursion for psi_pq can then be derived as:
237  //
238  // psi_pq = (pq)^-1/2 (r^2-(N-1)) psi_p-1,q-1
239  // - sqrt( (p-1)(q-1)/(pq) ) psi_p-2,q-2
240  //
241  // psi_p0 = (z/sqrt(p)) psi_p-1,0
242  //
243  // psi_00 = 1/sqrt(pi) exp(-r^2/2)
244 
245  Assert(int(psi.TMV_rowsize()) >= (order+1)*(order+2)/2);
246  Assert(psi.TMV_colsize() == z.size());
247  if (coeff) Assert(psi.TMV_colsize() == coeff->size());
248  //Assert(psi.iscm());
249  //Assert(!psi.isconj());
250 
251  // Setup rsq, z vectors and set psi_00
252  DVector rsq(z.size());
253  double* rsqit = TMV_ptr(rsq);
254  double* psi00it = TMV_ptr(psi);
255  const std::complex<double>* zit = TMV_cptr(z);
256  const int zsize = z.size();
257  for(int i=0;i<zsize;++i) {
258  rsqit[i] = std::norm(zit[i]);
259  psi00it[i] = afwGeom::INVSQRTPI * exp(-(rsqit[i])/2.);
260  }
261  if (coeff) psi.col(0).array() = coeff->array() * psi.col(0).array();
262 
263  DVector zr = z.TMV_realPart();
264  DVector zi = z.TMV_imagPart();
265  if (order >= 1) {
266  // Set psi_10
267  // All m > 0 elements are intrinsically complex.
268  // However, we are fitting to a real intensity pattern
269  // with complex coefficients of the complex shapelets.
270  // Since psi_pq = psi_qp* (* = complex conjugate),
271  // we know that b_pq must be b_qp*
272  // b_pq psi_pq + b_pq* psi_pq* = 2 b_pqr psi_pqr - 2 b_pqi psi_pqi
273  // So the values we want for the real fitter are
274  // really 2 real(psi_pq) and -2 imag(psi_pq)
275  // Putting the 2's here carries through to the rest of the
276  // elements via the recursion.
277  psi.col(1).array() = zr.array() * psi.col(0).array();
278  psi.col(2).array() = (-zi).array() * psi.col(0).array();
279  psi.col(1) *= 2.;
280  psi.col(2) *= 2.;
281  }
282  for(int N=2,k=3;N<=order;++N) {
283  // Set psi_N0
284  // The signs of these are not what you naively think due to
285  // the +2, -2 discussed above. You just have to follow through
286  // what the complex psi_N0 is, and what value is stored in the
287  // psi_N-1,0 location, and what needs to get stored here.
288  psi.col(k).array() = zr.array() * psi.col(k-N).array();
289  psi.col(k).array() += zi.array() * psi.col(k-N+1).array();
290  psi.col(k+1).array() = zr.array() * psi.col(k-N+1).array();
291  psi.col(k+1).array() -= zi.array() * psi.col(k-N).array();
292  double sqrt_1_N = sqrt(1./N);
293  psi.col(k) *= sqrt_1_N;
294  psi.col(k+1) *= sqrt_1_N;
295  k+=2;
296 
297  // Set psi_pq with q>0
298  // The rsq part of this calculation can be done in batch, which
299  // speeds things up a bit.
300  TMV_colRange(psi,k,k+N-1) = rsq.asDiagonal() * TMV_colRange(psi,k-2*N-1,k-N-2);
301  TMV_colRange(psi,k,k+N-1) -= (N-1.) * TMV_colRange(psi,k-2*N-1,k-N-2);
302  // The other calculation steps are different for each component:
303  for(int m=N-2,p=N-1,q=1;m>=0;--p,++q,m-=2) {
304  double pq = p*q;
305  if (m==0) {
306  psi.col(k) /= sqrt(pq);
307  if (q > 1) psi.col(k) -= sqrt(1.-(N-1.)/pq)*psi.col(k+2-4*N);
308  ++k;
309  } else {
310  TMV_colRange(psi,k,k+2) /= sqrt(pq);
311  if (q > 1)
312  TMV_colRange(psi,k,k+2) -=
313  sqrt(1.-(N-1.)/pq)*TMV_colRange(psi,k+2-4*N,k+4-4*N);
314  k+=2;
315  }
316  }
317  }
318  }
319 
320  void augmentPsi(DMatrix& psi, CDVectorView z, int order)
321  {
322  Assert(int(psi.TMV_rowsize()) >= (order+3)*(order+4)/2);
323  Assert(psi.TMV_colsize() == z.size());
324  Assert(order >= 1);
325  //Assert(psi.iscm());
326  //Assert(!psi.isconj());
327 
328  DVector rsq(z.size());
329  double* rsqit = TMV_ptr(rsq);
330  const std::complex<double>* zit = TMV_cptr(z);
331  const int zsize = z.size();
332  for(int i=0;i<zsize;++i) {
333  rsqit[i] = std::norm(zit[i]);
334  }
335 
336  DVector zr = z.TMV_realPart();
337  DVector zi = z.TMV_imagPart();
338  for(int N=order+1,k=N*(N+1)/2;N<=order+2;++N) {
339  psi.col(k).array() = zr.array() * psi.col(k-N).array();
340  psi.col(k).array() += zi.array() * psi.col(k-N+1).array();
341  psi.col(k+1).array() = zr.array() * psi.col(k-N+1).array();
342  psi.col(k+1).array() -= zi.array() * psi.col(k-N).array();
343  double sqrt_1_N = sqrt(1./N);
344  psi.col(k) *= sqrt_1_N;
345  psi.col(k+1) *= sqrt_1_N;
346  k+=2;
347 
348  TMV_colRange(psi,k,k+N-1) = rsq.asDiagonal() * TMV_colRange(psi,k-2*N-1,k-N-2);
349  TMV_colRange(psi,k,k+N-1) -=
350  (N-1.) * TMV_colRange(psi,k-2*N-1,k-N-2);
351  for(int m=N-2,p=N-1,q=1;m>=0;--p,++q,m-=2) {
352  double pq = p*q;
353  if (m==0) {
354  psi.col(k) /= sqrt(pq);
355  if (q > 1) psi.col(k) -= sqrt(1.-(N-1.)/pq)*psi.col(k+2-4*N);
356  ++k;
357  } else {
358  TMV_colRange(psi,k,k+2) /= sqrt(pq);
359  if (q > 1)
360  TMV_colRange(psi,k,k+2) -=
361  sqrt(1.-(N-1.)/pq)*TMV_colRange(psi,k+2-4*N,k+4-4*N);
362  k+=2;
363  }
364  }
365  }
366  }
367 
368 #endif
369 
370  void setupGx(DMatrix& Gx, int order1, int order2)
371  {
372  Assert(int(Gx.TMV_colsize()) == (order1+1)*(order1+2)/2);
373  Assert(int(Gx.TMV_rowsize()) == (order2+1)*(order2+2)/2);
374 
375  Gx.setZero();
376  for(int n=0,k=0;n<=order2;++n) for(int p=n,q=0;p>=q;--p,++q,++k) {
377  double dp = p;
378  double dq = q;
379  // d(psi(x,y))/dx = psi(x,y) Gx
380  // d/dx = 1/2(ap + aq - apt - aqt)
381  // Gx( p-1,q , pq ) = 1/2 sqrt(p)
382  // Gx( p,q-1 , pq ) = 1/2 sqrt(q)
383  // Gx( p+1,q , pq ) = -1/2 sqrt(p+1)
384  // Gx( p,q+1 , pq ) = -1/2 sqrt(q+1)
385  if (p==q) {
386  if (q>0 && n<=order1+1) Gx(k-n-2,k) = sqrt(dp)/2.;
387  if (n+1<=order1) Gx(k+n+1,k) = -sqrt(dp+1.)/2.;
388  } else if (p==q+1) {
389  if (n<=order1+1) Gx(k-n,k) = sqrt(dp);
390  if (q>0 && n<=order1+1) Gx(k-n-2,k) = sqrt(dq)/2.;
391  if (n+1<=order1) Gx(k+n+1,k) = -sqrt(dp+1.)/2.;
392  if (n+1<=order1) Gx(k+n+3,k) = -sqrt(dq+1.);
393  if (q>0 && n<=order1+1) Gx(k-n-1,k+1) = sqrt(dq)/2.;
394  if (n+1<=order1) Gx(k+n+2,k+1) = -sqrt(dp+1.)/2.;
395  } else {
396  if (n<=order1+1) Gx(k-n,k) = sqrt(dp)/2.;
397  if (q>0 && n<=order1+1) Gx(k-n-2,k) = sqrt(dq)/2.;
398  if (n+1<=order1) Gx(k+n+1,k) = -sqrt(dp+1.)/2.;
399  if (n+1<=order1) Gx(k+n+3,k) = -sqrt(dq+1.)/2.;
400  if (n<=order1+1) Gx(k-n+1,k+1) = sqrt(dp)/2.;
401  if (q>0 && n<=order1+1) Gx(k-n-1,k+1) = sqrt(dq)/2.;
402  if (n+1<=order1) Gx(k+n+2,k+1) = -sqrt(dp+1.)/2.;
403  if (n+1<=order1) Gx(k+n+4,k+1) = -sqrt(dq+1.)/2.;
404  }
405  if (p > q) ++k;
406  }
407  }
408 
409  void setupGy(DMatrix& Gy, int order1, int order2)
410  {
411  Assert(int(Gy.TMV_colsize()) == (order1+1)*(order1+2)/2);
412  Assert(int(Gy.TMV_rowsize()) == (order2+1)*(order2+2)/2);
413 
414  Gy.setZero();
415  for(int n=0,k=0;n<=order2;++n) for(int p=n,q=0;p>=q;--p,++q,++k) {
416  double dp = p;
417  double dq = q;
418  // d(psi(x,y))/dx = psi(x,y) Gx
419  // d/dy = 1/2 i (ap - aq + apt - aqt)
420  // Gy( p-1,q , pq ) = 1/2 i sqrt(p)
421  // Gy( p,q-1 , pq ) = -1/2 i sqrt(q)
422  // Gy( p+1,q , pq ) = 1/2 i sqrt(p+1)
423  // Gy( p,q+1 , pq ) = -1/2 i sqrt(q+1)
424  if (p==q) {
425  if (q>0 && n<=order1+1) Gy(k-n-1,k) = -sqrt(dp)/2.;
426  if (n+1<=order1) Gy(k+n+2,k) = sqrt(dp+1.)/2.;
427  } else if (p==q+1) {
428  if (q>0 && n<=order1+1) Gy(k-n-1,k) = -sqrt(dq)/2.;
429  if (n+1<=order1) Gy(k+n+2,k) = sqrt(dp+1.)/2.;
430  if (n<=order1+1) Gy(k-n,k+1) = -sqrt(dp);
431  if (q>0 && n<=order1+1) Gy(k-n-2,k+1) = sqrt(dq)/2.;
432  if (n+1<=order1) Gy(k+n+1,k+1) = -sqrt(dp+1.)/2.;
433  if (n+1<=order1) Gy(k+n+3,k+1) = sqrt(dq+1.);
434  } else {
435  if (n+1<=order1) Gy(k-n+1,k) = sqrt(dp)/2.;
436  if (q>0 && n<=order1+1) Gy(k-n-1,k) = -sqrt(dq)/2.;
437  if (n+1<=order1) Gy(k+n+2,k) = sqrt(dp+1.)/2.;
438  if (n+1<=order1) Gy(k+n+4,k) = -sqrt(dq+1.)/2.;
439  if (n<=order1+1) Gy(k-n,k+1) = -sqrt(dp)/2.;
440  if (q>0 && n<=order1+1) Gy(k-n-2,k+1) = sqrt(dq)/2.;
441  if (n+1<=order1) Gy(k+n+1,k+1) = -sqrt(dp+1.)/2.;
442  if (n+1<=order1) Gy(k+n+3,k+1) = sqrt(dq+1.)/2.;
443  }
444 
445  if (p > q) ++k;
446  }
447  }
448 
449  void setupGg1(DMatrix& Gg1, int order1, int order2)
450  {
451  Assert(int(Gg1.TMV_colsize()) == (order1+1)*(order1+2)/2);
452  Assert(int(Gg1.TMV_rowsize()) == (order2+1)*(order2+2)/2);
453 
454  Gg1.setZero();
455  for(int n=0,k=0;n<=order2;++n) for(int p=n,q=0;p>=q;--p,++q,++k) {
456  double dp = p;
457  double dq = q;
458  // d(psi(x,y))/dg1 = psi(x,y) Gg1
459  // Gg1 is G_g1
460  // d/dg1 = x d/dx - y d/dy
461  // = 1/2 (ap^2 + aq^2 - apt^2 - aqt^2)
462  // Gg1( p-2,q , pq ) = 1/2 sqrt(p(p-1))
463  // Gg1( p,q-2 , pq ) = 1/2 sqrt(q(q-1))
464  // Gg1( p+2,q , pq ) = -1/2 sqrt((p+1)(p+2))
465  // Gg1( p,q+2 , pq ) = -1/2 sqrt((q+1)(q+2))
466  if (p==q) {
467  if (q>1 && n<=order1+2) Gg1(k-2*n-3,k) = sqrt(dp*(dp-1.))/2.;
468  if (n+2<=order1) Gg1(k+2*n+3,k) = -sqrt((dp+1.)*(dp+2.))/2.;
469  } else if (p==q+1) {
470  if (p>1 && n<=order1+2) Gg1(k-2*n-1,k) = sqrt(dp*(dp-1.))/2.;
471  if (q>1 && n<=order1+2) Gg1(k-2*n-3,k) = sqrt(dq*(dq-1.))/2.;
472  if (n+2<=order1) Gg1(k+2*n+3,k) = -sqrt((dp+1.)*(dp+2.))/2.;
473  if (n+2<=order1) Gg1(k+2*n+5,k) = -sqrt((dq+1.)*(dq+2.))/2.;
474  if (p>1 && n<=order1+2) Gg1(k-2*n,k+1) = -sqrt(dp*(dp-1.))/2.;
475  if (q>1 && n<=order1+2) Gg1(k-2*n-2,k+1) = sqrt(dq*(dq-1.))/2.;
476  if (n+2<=order1) Gg1(k+2*n+4,k+1) = -sqrt((dp+1)*(dp+2.))/2.;
477  if (n+2<=order1) Gg1(k+2*n+6,k+1) = sqrt((dq+1)*(dq+2.))/2.;
478  } else if (p==q+2) {
479  if (n<=order1+2) Gg1(k-2*n+1,k) = sqrt(dp*(dp-1.));
480  if (q>1 && n<=order1+2) Gg1(k-2*n-3,k) = sqrt(dq*(dq-1.))/2.;
481  if (n+2<=order1) Gg1(k+2*n+3,k) = -sqrt((dp+1.)*(dp+2.))/2.;
482  if (n+2<=order1) Gg1(k+2*n+7,k) = -sqrt((dq+1.)*(dq+2.));
483  if (q>1 && n<=order1+2) Gg1(k-2*n-2,k+1) = sqrt(dq*(dq-1.))/2.;
484  if (n+2<=order1) Gg1(k+2*n+4,k+1) = -sqrt((dp+1.)*(dp+2.))/2.;
485  } else {
486  if (n<=order1+2) Gg1(k-2*n+1,k) = sqrt(dp*(dp-1.))/2.;
487  if (q>1 && n<=order1+2) Gg1(k-2*n-3,k) = sqrt(dq*(dq-1.))/2.;
488  if (n+2<=order1) Gg1(k+2*n+3,k) = -sqrt((dp+1.)*(dp+2.))/2.;
489  if (n+2<=order1) Gg1(k+2*n+7,k) = -sqrt((dq+1.)*(dq+2.))/2.;
490  if (n<=order1+2) Gg1(k-2*n+2,k+1) = sqrt(dp*(dp-1.))/2.;
491  if (q>1 && n<=order1+2) Gg1(k-2*n-2,k+1) = sqrt(dq*(dq-1.))/2.;
492  if (n+2<=order1) Gg1(k+2*n+4,k+1) = -sqrt((dp+1.)*(dp+2.))/2.;
493  if (n+2<=order1) Gg1(k+2*n+8,k+1) = -sqrt((dq+1.)*(dq+2.))/2.;
494  }
495 
496  if (p > q) ++k;
497  }
498  }
499 
500  void setupGg2(DMatrix& Gg2, int order1, int order2)
501  {
502  Assert(int(Gg2.TMV_colsize()) == (order1+1)*(order1+2)/2);
503  Assert(int(Gg2.TMV_rowsize()) == (order2+1)*(order2+2)/2);
504 
505  Gg2.setZero();
506  for(int n=0,k=0;n<=order2;++n) for(int p=n,q=0;p>=q;--p,++q,++k) {
507  double dp = p;
508  double dq = q;
509  // d(psi(x,y))/dg2 = psi(x,y) Gg2
510  // Gg2 is G_g2
511  // d/dg2 = y d/dx + x d/dy
512  // = 1/2 i (ap^2 - aq^2 + apt^2 - aqt^2)
513  // Gg2( p-2,q , pq ) = 1/2 i sqrt(p(p-1))
514  // Gg2( p,q-2 , pq ) = -1/2 i sqrt(q(q-1))
515  // Gg2( p+2,q , pq ) = 1/2 i sqrt((p+1)(p+2))
516  // Gg2( p,q+2 , pq ) = -1/2 i sqrt((q+1)(q+2))
517  if (p==q) {
518  if (q>1 && n<=order1+2) Gg2(k-2*n-2,k) = -sqrt(dp*(dp-1.))/2.;
519  if (n+2<=order1) Gg2(k+2*n+4,k) = sqrt((dp+1.)*(dp+2.))/2.;
520  } else if (p==q+1) {
521  if (p>1 && n<=order1+2) Gg2(k-2*n,k) = -sqrt(dp*(dp-1.))/2.;
522  if (q>1 && n<=order1+2) Gg2(k-2*n-2,k) = -sqrt(dq*(dq-1.))/2.;
523  if (n+2<=order1) Gg2(k+2*n+4,k) = sqrt((dp+1.)*(dp+2.))/2.;
524  if (n+2<=order1) Gg2(k+2*n+6,k) = sqrt((dq+1.)*(dq+2.))/2.;
525  if (p>1 && n<=order1+2) Gg2(k-2*n-1,k+1) = -sqrt(dp*(dp-1.))/2.;
526  if (q>1 && n<=order1+2) Gg2(k-2*n-3,k+1) = sqrt(dq*(dq-1.))/2.;
527  if (n+2<=order1) Gg2(k+2*n+3,k+1) = -sqrt((dp+1.)*(dp+2.))/2.;
528  if (n+2<=order1) Gg2(k+2*n+5,k+1) = sqrt((dq+1.)*(dq+2.))/2.;
529  } else if (p==q+2) {
530  if (q>1 && n<=order1+2) Gg2(k-2*n-2,k) = -sqrt(dq*(dq-1.))/2.;
531  if (n+2<=order1) Gg2(k+2*n+4,k) = sqrt((dp+1.)*(dp+2.))/2.;
532  if (n<=order2+2) Gg2(k-2*n+1,k+1) = -sqrt(dp*(dp-1.));
533  if (q>1 && n<=order1+2) Gg2(k-2*n-3,k+1) = sqrt(dq*(dq-1.))/2.;
534  if (n+2<=order1) Gg2(k+2*n+3,k+1) = -sqrt((dp+1.)*(dp+2.))/2.;
535  if (n+2<=order1) Gg2(k+2*n+7,k+1) = sqrt((dq+1.)*(dq+2.));
536  } else {
537  Gg2(k-2*n+2,k) = sqrt(dp*(dp-1.))/2.;
538  if (q>1 && n<=order1+2) Gg2(k-2*n-2,k) = -sqrt(dq*(dq-1.))/2.;
539  if (n+2<=order1) Gg2(k+2*n+4,k) = sqrt((dp+1.)*(dp+2.))/2.;
540  if (n+2<=order1) Gg2(k+2*n+8,k) = -sqrt((dq+1.)*(dq+2.))/2.;
541  if (n<=order2+2) Gg2(k-2*n+1,k+1) = -sqrt(dp*(dp-1.))/2.;
542  if (q>1 && n<=order1+2) Gg2(k-2*n-3,k+1) = sqrt(dq*(dq-1.))/2.;
543  if (n+2<=order1) Gg2(k+2*n+3,k+1) = -sqrt((dp+1.)*(dp+2.))/2.;
544  if (n+2<=order1) Gg2(k+2*n+7,k+1) = sqrt((dq+1.)*(dq+2.))/2.;
545  }
546 
547  if (p > q) ++k;
548  }
549  }
550 
551  void setupGmu(DMatrix& Gmu, int order1, int order2)
552  {
553  Assert(int(Gmu.TMV_colsize()) == (order1+1)*(order1+2)/2);
554  Assert(int(Gmu.TMV_rowsize()) == (order2+1)*(order2+2)/2);
555 
556  Gmu.setZero();
557  for(int n=0,k=0;n<=order2;++n) for(int p=n,q=0;p>=q;--p,++q,++k) {
558  double dp = p;
559  double dq = q;
560  // d(psi(x,y))/dmu = psi(x,y) Gmu
561  // d/dmu = x d/dx + y d/dy
562  // ap aq - apt aqt - 1
563  // Gmu( p-1,q-1 , pq ) = sqrt(pq)
564  // Gmu( p+1,q+1 , pq ) = -sqrt((p+1)(q+1))
565  // Gmu( pq , pq ) = -1
566  if (q>0 && n<=order1+2) Gmu(k-2*n-1,k) = sqrt(dp*dq);
567  if (n+2<=order1) Gmu(k+2*n+5,k) = -sqrt((dp+1.)*(dq+1.));
568  if (n<=order1) Gmu(k,k) = -1.;
569  if (p > q) {
570  if (q>0 && n<=order1+2) Gmu(k-2*n,k+1) = sqrt(dp*dq);
571  if (n+2<=order1) Gmu(k+2*n+6,k+1) = -sqrt((dp+1.)*(dq+1.));
572  if (n<=order1) Gmu(k+1,k+1) = -1.;
573  }
574 
575  if (p > q) ++k;
576  }
577  }
578 
579  void setupGth(DMatrix& Gth, int order1, int order2)
580  {
581  Assert(int(Gth.TMV_colsize()) == (order1+1)*(order1+2)/2);
582  Assert(int(Gth.TMV_rowsize()) == (order2+1)*(order2+2)/2);
583 
584  Gth.setZero();
585  for(int n=0,k=0;n<=order2;++n) for(int p=n,q=0;p>=q;--p,++q,++k) {
586  double dp = p;
587  double dq = q;
588  // d(psi(x,y))/dtheta = psi(x,y) Gth
589  // d/dtheta = x d/dy - y d/dx
590  // = m i
591  // Gth( pq , pq ) = m i
592  if (p>q && n<=order1) {
593  Gth(k+1,k) = (dp-dq);
594  Gth(k,k+1) = -(dp-dq);
595  }
596 
597  if (p > q) ++k;
598  }
599  }
600 
601 }}}}
void setupGg1(DMatrix &Gg1, int order1, int order2)
Definition: PsiHelper.cc:449
void setupGth(DMatrix &Gth, int order1, int order2)
Definition: PsiHelper.cc:579
T norm(const T &x)
Definition: Integrate.h:191
Eigen::Block< CDVector, Eigen::Dynamic, 1 > CDVectorView
Definition: MyMatrix.h:140
#define TMV_cptr(m)
Definition: MyMatrix.h:206
#define TMV_ptr(m)
Definition: MyMatrix.h:205
void makePsi(DMatrix &psi, CDVectorView z, int order, const DVector *coeff=0)
Definition: PsiHelper.cc:226
#define TMV_colRange(m, j1, j2)
Definition: MyMatrix.h:203
void setupGx(DMatrix &Gx, int order1, int order2)
Definition: PsiHelper.cc:370
void setupGy(DMatrix &Gy, int order1, int order2)
Definition: PsiHelper.cc:409
tuple m
Definition: lsstimport.py:48
void setupGg2(DMatrix &Gg2, int order1, int order2)
Definition: PsiHelper.cc:500
void augmentPsi(DMatrix &psi, CDVectorView z, int order)
Definition: PsiHelper.cc:320
double const INVSQRTPI
Definition: Angle.h:23
T real(const T &x)
Definition: Integrate.h:193
#define Assert(x)
Definition: dbg.h:73
void setupGmu(DMatrix &Gmu, int order1, int order2)
Definition: PsiHelper.cc:551