LSSTApplications  10.0+286,10.0+36,10.0+46,10.0-2-g4f67435,10.1+152,10.1+37,11.0,11.0+1,11.0-1-g47edd16,11.0-1-g60db491,11.0-1-g7418c06,11.0-2-g04d2804,11.0-2-g68503cd,11.0-2-g818369d,11.0-2-gb8b8ce7
LSSTDataManagementBasePackage
Classes | Typedefs | Enumerations | Functions
lsst::meas::algorithms::shapelet Namespace Reference

Classes

class  Position
 
class  Bounds
 
class  AssignableToBVec
 
class  BVec
 
class  ConvertibleString
 
class  ConfigFile
 
class  Ellipse
 
class  EllipseSolver3
 
class  FittedPsf
 
class  FittedPsfAtXY
 
class  BoundForm
 
class  Form
 
struct  RangeException
 
class  Function2D
 
class  Constant2D
 
class  Polynomial2D
 
class  Histogram
 
class  Legendre2D
 
class  NLSolver
 
struct  FileNotFoundError
 
struct  ParameterException
 
struct  ReadException
 
struct  WriteException
 
struct  ProcessingException
 
class  Pixel
 
class  PixelList
 
class  PotentialStar
 
class  SizeMagnitudeStarSelectorException
 
class  SizeMagnitudeStarSelectorAlgo
 
class  CrudeSolver
 
struct  PixelListSorter
 

Typedefs

typedef Eigen::VectorXd DVector
 
typedef Eigen::Block< DVector,
Eigen::Dynamic, 1 > 		DVectorView
 
typedef const DVectorView DConstVectorView
 
typedef Eigen::MatrixXd DMatrix
 
typedef Eigen::Block< DMatrix > DMatrixView
 
typedef const DMatrixView DConstMatrixView
 
typedef Eigen::Matrix< double, 2, 2 > DSmallMatrix22
 
typedef Eigen::RowVectorXd DRowVector
 
typedef Eigen::MatrixXd DBandMatrix
 
typedef Eigen::VectorXcd CDVector
 
typedef Eigen::Block< CDVector,
Eigen::Dynamic, 1 > 		CDVectorView
 
typedef const CDVectorView CDConstVectorView
 
typedef Eigen::MatrixXcd CDMatrix
 
typedef Eigen::Block< CDMatrix > CDMatrixView
 
typedef const CDMatrixView CDConstMatrixView
 
typedef Eigen::VectorXf FVector
 
typedef Eigen::Block< FVector,
Eigen::Dynamic, 1 > 		FVectorView
 
typedef const FVectorView FConstVectorView
 
typedef Eigen::MatrixXf FMatrix
 
typedef Eigen::Block< FMatrix > FMatrixView
 
typedef const FMatrixView FConstMatrixView
 
typedef Eigen::Matrix< float, 2, 2 > FSmallMatrix22
 
typedef Eigen::VectorXcf CFVector
 
typedef Eigen::Block< CFVector,
Eigen::Dynamic, 1 > 		CFVectorView
 
typedef const CFVectorView CFConstVectorView
 
typedef Eigen::MatrixXcf CFMatrix
 
typedef Eigen::Block< CFMatrix > CFMatrixView
 
typedef const CFMatrixView CFConstMatrixView
 
typedef DVector DDiagMatrix
 
typedef DMatrix DSymMatrix
 

Enumerations

enum  ExitCode {
  SUCCESS = 0, FAILURE, FAILURE_FILE_NOT_FOUND, FAILURE_PARAMETER_ERROR,
  FAILURE_READ_ERROR, FAILURE_WRITE_ERROR, FAILURE_PROCESSING_ERROR
}
 

Functions

double fact (int i)
 
double sqrtfact (int i)
 
double binom (int i, int j)
 
double sqrtn (int i)
 
std::complex< double > operator- (const std::complex< double > &z1, const Position &p2)
 
std::complex< double > operator/ (const Position &p1, double x)
 
std::ostream & operator<< (std::ostream &os, const Position &pos)
 
std::istream & operator>> (std::istream &os, Position &pos)
 
std::ostream & operator<< (std::ostream &fout, const Bounds &b)
 
std::istream & operator>> (std::istream &fin, Bounds &b)
 
std::ostream & operator<< (std::ostream &os, const BVec &b)
 
void calculateZTransform (std::complex< double > z, int order1, int order2, DMatrix &T)
 
void calculateZTransform (std::complex< double > z, int order, DMatrix &T)
 
void augmentZTransformCols (std::complex< double > z, int order, DMatrix &T)
 
void applyZ (std::complex< double > z, BVec &b)
 
void calculateMuTransform (double mu, int order1, int order2, DMatrix &D)
 
void calculateMuTransform (double mu, int order, DMatrix &D)
 
void augmentMuTransformRows (double mu, int order, DMatrix &D)
 
void augmentMuTransformCols (double mu, int order, DMatrix &D)
 
void applyMu (double mu, BVec &b)
 
void calculateThetaTransform (double theta, int order1, int order2, DBandMatrix &R)
 
void calculateThetaTransform (double theta, int order, DBandMatrix &R)
 
void applyTheta (double theta, BVec &b)
 
void calculateGTransform (std::complex< double > g, int order1, int order2, DMatrix &S)
 
void calculateGTransform (std::complex< double > g, int order, DMatrix &S)
 
void augmentGTransformCols (std::complex< double > g, int order, DMatrix &S)
 
void applyG (std::complex< double > g, BVec &b)
 
void calculatePsfConvolve (const BVec &bpsf, int order1, int order2, double sigma, DMatrix &C)
 
void calculatePsfConvolve (const BVec &bpsf, int order, double sigma, DMatrix &C)
 
void applyPsf (const BVec &bpsf, BVec &b)
 
std::ostream & operator<< (std::ostream &os, const ConfigFile &cf)
 
std::istream & operator>> (std::istream &is, ConfigFile &cf)
 
std::ostream & operator<< (std::ostream &os, const Ellipse &s)
 
std::complex< double > addShears (const std::complex< double > g1, const std::complex< double > g2)
 
void CalculateShift (std::complex< double > c1, std::complex< double > g1, std::complex< double > m1, std::complex< double > c2, std::complex< double > g2, std::complex< double > m2, std::complex< double > &c3, std::complex< double > &g3, std::complex< double > &m3)
 
void setupFloat (std::ostream &s, const Form &f)
 
void setupInt (std::ostream &s, const Form &f)
 
void setup (std::ostream &os, const BoundForm< double > &bf)
 
void setup (std::ostream &os, const BoundForm< long double > &bf)
 
void setup (std::ostream &os, const BoundForm< float > &bf)
 
void setup (std::ostream &os, const BoundForm< std::complex< double > > &bf)
 
void setup (std::ostream &os, const BoundForm< std::complex< long double > > &bf)
 
void setup (std::ostream &os, const BoundForm< std::complex< float > > &bf)
 
void setup (std::ostream &os, const BoundForm< int > &bf)
 
void setup (std::ostream &os, const BoundForm< short > &bf)
 
void setup (std::ostream &os, const BoundForm< long > &bf)
 
void setup (std::ostream &os, const BoundForm< unsigned int > &bf)
 
void setup (std::ostream &os, const BoundForm< unsigned short > &bf)
 
void setup (std::ostream &os, const BoundForm< unsigned long > &bf)
 
template<typename T >
void setup (std::ostream &os, const BoundForm< T > &bf)
 
template<typename T >
std::ostream & operator<< (std::ostream &os, const BoundForm< T > &bf)
 
std::ostream & operator<< (std::ostream &fout, const Function2D &f)
 
std::istream & operator>> (std::istream &fin, std::auto_ptr< Function2D > &f)
 
void PrintFlags (const std::vector< long > &flags, std::ostream &os)
 
const char * Text (const ExitCode &code)
 
int Status (ExitCode code, const ConfigFile &params)
 
void getSubPixList (PixelList &pix, const PixelList &allPix, std::complex< double > cen_offset, std::complex< double > shear, double aperture, long &flag)
 
void getSubPixList (PixelList &pix, const PixelList &allPix, std::complex< double > cen_offset, double aperture, long &flag)
 
void getSubPixList (PixelList &pix, const PixelList &allPix, double aperture, long &flag)
 
void makePsi (DMatrix &psi, CDVectorView z, int order, const DVector *coeff=0)
 
void makePsi (DVector &psi, std::complex< double > z, int order)
 
void augmentPsi (DMatrix &psi, CDVectorView z, int order)
 
void setupGx (DMatrix &Gx, int order1, int order2)
 
void setupGy (DMatrix &Gy, int order1, int order2)
 
void setupGg1 (DMatrix &Gg1, int order1, int order2)
 
void setupGg2 (DMatrix &Gg2, int order1, int order2)
 
void setupGmu (DMatrix &Gmu, int order1, int order2)
 
void setupGth (DMatrix &Gth, int order1, int order2)
 
int fitSize (const int xOrder, const int yOrder)
 
double absSq (const double &x)
 
double betai (double a, double b, double x)
 
bool Equivalent (double chisq1, double chisq2, int n1, int n2, double equivProb)
 
double betacf (double a, double b, double x)
 
double gammln (double x)
 

Typedef Documentation

typedef const CDMatrixView lsst::meas::algorithms::shapelet::CDConstMatrixView

Definition at line 144 of file MyMatrix.h.

typedef const CDVectorView lsst::meas::algorithms::shapelet::CDConstVectorView

Definition at line 141 of file MyMatrix.h.

typedef Eigen::MatrixXcd lsst::meas::algorithms::shapelet::CDMatrix

Definition at line 142 of file MyMatrix.h.

typedef Eigen::Block<CDMatrix> lsst::meas::algorithms::shapelet::CDMatrixView

Definition at line 143 of file MyMatrix.h.

typedef Eigen::VectorXcd lsst::meas::algorithms::shapelet::CDVector

Definition at line 139 of file MyMatrix.h.

typedef Eigen::Block<CDVector,Eigen::Dynamic,1> lsst::meas::algorithms::shapelet::CDVectorView

Definition at line 140 of file MyMatrix.h.

typedef const CFMatrixView lsst::meas::algorithms::shapelet::CFConstMatrixView

Definition at line 159 of file MyMatrix.h.

typedef const CFVectorView lsst::meas::algorithms::shapelet::CFConstVectorView

Definition at line 156 of file MyMatrix.h.

typedef Eigen::MatrixXcf lsst::meas::algorithms::shapelet::CFMatrix

Definition at line 157 of file MyMatrix.h.

typedef Eigen::Block<CFMatrix> lsst::meas::algorithms::shapelet::CFMatrixView

Definition at line 158 of file MyMatrix.h.

typedef Eigen::VectorXcf lsst::meas::algorithms::shapelet::CFVector

Definition at line 154 of file MyMatrix.h.

typedef Eigen::Block<CFVector,Eigen::Dynamic,1> lsst::meas::algorithms::shapelet::CFVectorView

Definition at line 155 of file MyMatrix.h.

typedef Eigen::MatrixXd lsst::meas::algorithms::shapelet::DBandMatrix

Definition at line 137 of file MyMatrix.h.

typedef const DMatrixView lsst::meas::algorithms::shapelet::DConstMatrixView

Definition at line 134 of file MyMatrix.h.

typedef const DVectorView lsst::meas::algorithms::shapelet::DConstVectorView

Definition at line 131 of file MyMatrix.h.

typedef DVector lsst::meas::algorithms::shapelet::DDiagMatrix

Definition at line 161 of file MyMatrix.h.

typedef Eigen::MatrixXd lsst::meas::algorithms::shapelet::DMatrix

Definition at line 132 of file MyMatrix.h.

typedef Eigen::Block<DMatrix> lsst::meas::algorithms::shapelet::DMatrixView

Definition at line 133 of file MyMatrix.h.

typedef Eigen::RowVectorXd lsst::meas::algorithms::shapelet::DRowVector

Definition at line 136 of file MyMatrix.h.

typedef Eigen::Matrix<double,2,2> lsst::meas::algorithms::shapelet::DSmallMatrix22

Definition at line 135 of file MyMatrix.h.

typedef DMatrix lsst::meas::algorithms::shapelet::DSymMatrix

Definition at line 162 of file MyMatrix.h.

typedef Eigen::VectorXd lsst::meas::algorithms::shapelet::DVector

Definition at line 129 of file MyMatrix.h.

typedef Eigen::Block<DVector,Eigen::Dynamic,1> lsst::meas::algorithms::shapelet::DVectorView

Definition at line 130 of file MyMatrix.h.

typedef const FMatrixView lsst::meas::algorithms::shapelet::FConstMatrixView

Definition at line 151 of file MyMatrix.h.

typedef const FVectorView lsst::meas::algorithms::shapelet::FConstVectorView

Definition at line 148 of file MyMatrix.h.

typedef Eigen::MatrixXf lsst::meas::algorithms::shapelet::FMatrix

Definition at line 149 of file MyMatrix.h.

typedef Eigen::Block<FMatrix> lsst::meas::algorithms::shapelet::FMatrixView

Definition at line 150 of file MyMatrix.h.

typedef Eigen::Matrix<float,2,2> lsst::meas::algorithms::shapelet::FSmallMatrix22

Definition at line 152 of file MyMatrix.h.

typedef Eigen::VectorXf lsst::meas::algorithms::shapelet::FVector

Definition at line 146 of file MyMatrix.h.

typedef Eigen::Block<FVector,Eigen::Dynamic,1> lsst::meas::algorithms::shapelet::FVectorView

Definition at line 147 of file MyMatrix.h.

Enumeration Type Documentation

enum lsst::meas::algorithms::shapelet::ExitCode
Enumerator
SUCCESS 
FAILURE 
FAILURE_FILE_NOT_FOUND 
FAILURE_PARAMETER_ERROR 
FAILURE_READ_ERROR 
FAILURE_WRITE_ERROR 
FAILURE_PROCESSING_ERROR 

Definition at line 139 of file Params.h.

139  {
140  SUCCESS = 0,
141  FAILURE,
142  FAILURE_FILE_NOT_FOUND,
143  FAILURE_PARAMETER_ERROR,
144  FAILURE_READ_ERROR,
145  FAILURE_WRITE_ERROR,
146  FAILURE_PROCESSING_ERROR
147  };

Function Documentation

double lsst::meas::algorithms::shapelet::absSq ( const double &  x)
inline

Definition at line 530 of file Function2D.cc.

530 { return x*x; }
x
double x
Definition: ChebyshevBoundedField.cc:305
std::complex<double> lsst::meas::algorithms::shapelet::addShears ( const std::complex< double >  g1,
const std::complex< double >  g2 
)
void lsst::meas::algorithms::shapelet::applyG ( std::complex< double >  g,
BVec &  b 
)

Definition at line 666 of file BVec.cc.

667  {
668  if (g != 0.0) {
669  DMatrix S(int(b.size()),int(b.size()));
670  calculateGTransform(g,b.getOrder(),S);
671  b.vec() = S * b.vec();
672  }
673  }
lsst::meas::algorithms::shapelet::calculateGTransform
void calculateGTransform(std::complex< double > g, int order1, int order2, DMatrix &S)
Definition: BVec.cc:468
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
void lsst::meas::algorithms::shapelet::applyMu ( double  mu,
BVec &  b 
)

Definition at line 410 of file BVec.cc.

411  {
412  if (mu != 0.0) {
413  DMatrix D(int(b.size()),int(b.size()));
414  calculateMuTransform(mu,b.getOrder(),D);
415  b.vec() = D * b.vec();
416  }
417  }
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
lsst::meas::algorithms::shapelet::calculateMuTransform
void calculateMuTransform(double mu, int order1, int order2, DMatrix &D)
Definition: BVec.cc:227
void lsst::meas::algorithms::shapelet::applyPsf ( const BVec &  bpsf,
BVec &  b 
)

Definition at line 895 of file BVec.cc.

896  {
897  DMatrix C(int(b.size()),int(b.size()));
898  calculatePsfConvolve(bpsf,b.getOrder(),b.getSigma(),C);
899  b.vec() = C * b.vec();
900  }
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
lsst::meas::algorithms::shapelet::calculatePsfConvolve
void calculatePsfConvolve(const BVec &bpsf, int order1, int order2, double sigma, DMatrix &C)
Definition: BVec.cc:675
void lsst::meas::algorithms::shapelet::applyTheta ( double  theta,
BVec &  b 
)

Definition at line 455 of file BVec.cc.

456  {
457  if (theta != 0.0) {
458 #ifdef USE_TMV
459  DBandMatrix R(int(b.size()),int(b.size()),1,1);
460 #else
461  DBandMatrix R(int(b.size()),int(b.size()));
462 #endif
463  calculateThetaTransform(theta,b.getOrder(),R);
464  b.vec() = R * b.vec();
465  }
466  }
lsst::meas::algorithms::shapelet::calculateThetaTransform
void calculateThetaTransform(double theta, int order1, int order2, DBandMatrix &R)
Definition: BVec.cc:419
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
void lsst::meas::algorithms::shapelet::applyZ ( std::complex< double >  z,
BVec &  b 
)

Definition at line 218 of file BVec.cc.

219  {
220  if (z != 0.0) {
221  DMatrix T(int(b.size()),int(b.size()));
222  calculateZTransform(z/b.getSigma(),b.getOrder(),T);
223  b.vec() = T * b.vec();
224  }
225  }
lsst::meas::algorithms::shapelet::calculateZTransform
void calculateZTransform(std::complex< double > z, int order1, int order2, DMatrix &T)
Definition: BVec.cc:84
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
void lsst::meas::algorithms::shapelet::augmentGTransformCols ( std::complex< double >  g,
int  order,
DMatrix &  S 
)

Definition at line 577 of file BVec.cc.

578  {
579  const int order1 = order;
580  const int order2 = order+2;
581  const int size1 = (order1+1)*(order1+2)/2;
582  //const int size2 = (order2+1)*(order2+2)/2;
583  Assert(int(S.TMV_colsize()) == size1);
584  Assert(int(S.TMV_rowsize()) == (order2+1)*(order2+2)/2);
585 
586  if (g == 0.0) return;
587 
588  double absg = std::abs(g);
589  double normg = std::norm(g);
590  std::vector<std::complex<double> > phase(order1+order2+1);
591  phase[0] = 1.;
592  std::complex<double> ph = -std::conj(g)/absg;
593  for(int i=1;i<=order1+order2;++i) phase[i] = phase[i-1]*ph;
594 
595  double te = absg;
596  double se = sqrt(1.-normg);
597 
598  std::vector<std::vector<double> > f(
599  order2+1,std::vector<double>(order1+1,0.));
600 
601  f[0][0] = sqrt(se); // f(0,0)
602  // only terms with p+s even are non-zero.
603  for(int p=2;p<=order2;p+=2) {
604  f[p][0] = f[p-2][0]*(-te/2.)*sqrtn(p*(p-1))/double(p/2); // f(p,0)
605  }
606 
607  for(int s=0;s<order1;++s) {
608  if (s%2==1) {
609  f[0][s+1] = sqrtn(s)*te*f[0][s-1]/sqrtn(s+1); // f(0,s+1)
610  }
611  for(int p=s%2+1;p<=order2;p+=2) {
612  double temp = sqrtn(p)*se*f[p-1][s];
613  if (s>0) temp += sqrtn(s)*te*f[p][s-1];
614  temp /= sqrtn(s+1);
615  f[p][s+1] = temp; // f(p,s+1)
616  }
617  }
618 
619  for(int n=order1+1,pq=size1;n<=order2;++n) {
620  for(int p=n,q=0;p>=q;--p,++q,++pq) {
621  const std::vector<double>& fp = f[p];
622  const std::vector<double>& fq = f[q];
623  double* Spq = TMV_ptr(S.col(pq));
624  double* Spq1 = p>q?TMV_ptr(S.col(pq+1)):0;
625  for(int nn=n%2,st=(nn==0?0:1);nn<=order1;nn+=2,st+=nn) {
626  for(int s=nn,t=0;s>=t;--s,++t,++st) {
627 
628  double s0 = fp[s] * fq[t];
629  int iphase = s-t-p+q;
630  std::complex<double> s1 = s0 *
631  (iphase >= 0 ?
632  phase[iphase/2] :
633  std::conj(phase[-iphase/2]));
634 
635  if (p==q) {
636  if (s==t) {
637  Spq[st] = std::real(s1);
638  } else {
639  Spq[st] = std::real(s1);
640  Spq[st+1] = std::imag(s1);
641  ++st;
642  }
643  } else if (s==t) {
644  Spq[st] = 2.*std::real(s1);
645  Spq1[st] = -2.*std::imag(s1);
646  } else {
647  s0 = fq[s] * fp[t];
648  iphase = s-t-q+p;
649  std::complex<double> s2 = s0 *
650  (iphase >= 0 ?
651  phase[iphase/2] :
652  std::conj(phase[-iphase/2]));
653  Spq[st]= std::real(s1) + std::real(s2);
654  Spq1[st] = -std::imag(s1) + std::imag(s2);
655  Spq[st+1] = std::imag(s1) + std::imag(s2);
656  Spq1[st+1] = std::real(s1) - std::real(s2);
657  ++st;
658  }
659  }
660  }
661  if (p!=q) ++pq;
662  }
663  }
664  }
lsst::afw::math::details::norm
T norm(const T &x)
Definition: Integrate.h:191
TMV_ptr
#define TMV_ptr(m)
Definition: MyMatrix.h:205
lsst::afw::math::details::real
T real(const T &x)
Definition: Integrate.h:193
Assert
#define Assert(x)
Definition: dbg.h:73
void lsst::meas::algorithms::shapelet::augmentMuTransformCols ( double  mu,
int  order,
DMatrix &  D 
)

Definition at line 316 of file BVec.cc.

317  {
318  const int order1 = order;
319  const int order2 = order+2;
320  const int size1 = (order1+1)*(order1+2)/2;
321  //const int size2 = (order2+1)*(order2+2)/2;
322  Assert(int(D.TMV_colsize()) >= size1);
323  Assert(int(D.TMV_rowsize()) >= (order2+1)*(order2+2)/2);
324 
325  if (mu == 0.0) return;
326 
327  double tmu = tanh(mu);
328  double smu = 1./cosh(mu);
329 
330  double Dqq = exp(-mu)*smu;
331  for(int q=order1/2;q>0;--q) Dqq *= tmu;
332  for(int n=order1+1,pq=size1;n<=order2;++n) {
333  for(int p=n,q=0;p>=q;--p,++q,++pq) {
334  double* Dpq = TMV_ptr(D.col(pq));
335  double* Dpq_n = TMV_ptr(D.col(pq-n));
336  double* Dpq_n_2 = q>0?TMV_ptr(D.col(pq-n-2)):0;
337  double* Dpq1 = p>q?TMV_ptr(D.col(pq+1)):0;
338  if (p==q) {
339  if (p > 0) Dqq *= tmu;
340  Dpq[0] = Dqq; // D(0,pq)
341  }
342  for(int m=1,st=1;m<=order1;++m) {
343  for(int s=m,t=0;s>=t;--s,++t,++st) {
344  if (p-q==s-t) {
345  double temp;
346  if (t == 0) {
347  temp = smu*sqrtn(p)*Dpq_n[st-m]/sqrtn(s);
348  } else {
349  temp = -tmu*sqrtn(s)*Dpq[st-2*m-1];
350  if (q > 0) {
351  temp += smu*sqrtn(q)*Dpq_n_2[st-m-2];
352  }
353  temp /= sqrtn(t);
354  }
355  Dpq[st] = temp;
356  if (s!=t) Dpq1[st+1] = temp;
357  }
358  if (s!=t) ++st;
359  }
360  }
361  if (p!=q) ++pq;
362  }
363  }
364  }
TMV_ptr
#define TMV_ptr(m)
Definition: MyMatrix.h:205
lsstimport.m
tuple m
Definition: lsstimport.py:48
Assert
#define Assert(x)
Definition: dbg.h:73
void lsst::meas::algorithms::shapelet::augmentMuTransformRows ( double  mu,
int  order,
DMatrix &  D 
)

Definition at line 366 of file BVec.cc.

367  {
368  const int order1 = order+2;
369  const int order2 = order;
370  //const int size1 = (order1+1)*(order1+2)/2;
371  const int size2 = (order2+1)*(order2+2)/2;
372  Assert(int(D.TMV_colsize()) == (order1+1)*(order1+2)/2);
373  Assert(int(D.TMV_rowsize()) == size2);
374 
375  if (mu == 0.0) return;
376 
377  double tmu = tanh(mu);
378  double smu = 1./cosh(mu);
379 
380  for(int n=0,pq=0;n<=order2;++n) {
381  for(int p=n,q=0;p>=q;--p,++q,++pq) {
382  double* Dpq = TMV_ptr(D.col(pq));
383  double* Dpq_n = TMV_ptr(D.col(pq-n));
384  double* Dpq_n_2 = q>0?TMV_ptr(D.col(pq-n-2)):0;
385  double* Dpq1 = p>q?TMV_ptr(D.col(pq+1)):0;
386  for(int m=order2+1,st=size2;m<=order1;++m) {
387  for(int s=m,t=0;s>=t;--s,++t,++st) {
388  if (p-q==s-t) {
389  double temp;
390  if (t == 0) {
391  temp = smu*sqrtn(p)*Dpq_n[st-m]/sqrtn(s);
392  } else {
393  temp = -tmu*sqrtn(s)*Dpq[st-2*m-1];
394  if (q > 0) {
395  temp += smu*sqrtn(q)*Dpq_n_2[st-m-2];
396  }
397  temp /= sqrtn(t);
398  }
399  Dpq[st] = temp;
400  if (s!=t) Dpq1[st+1] = temp;
401  }
402  if (s!=t) ++st;
403  }
404  }
405  if (p!=q) ++pq;
406  }
407  }
408  }
TMV_ptr
#define TMV_ptr(m)
Definition: MyMatrix.h:205
lsstimport.m
tuple m
Definition: lsstimport.py:48
Assert
#define Assert(x)
Definition: dbg.h:73
void lsst::meas::algorithms::shapelet::augmentPsi ( DMatrix &  psi,
CDVectorView  z,
int  order 
)

Definition at line 320 of file PsiHelper.cc.

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  }
lsst::afw::math::details::norm
T norm(const T &x)
Definition: Integrate.h:191
TMV_cptr
#define TMV_cptr(m)
Definition: MyMatrix.h:206
TMV_ptr
#define TMV_ptr(m)
Definition: MyMatrix.h:205
TMV_colRange
#define TMV_colRange(m, j1, j2)
Definition: MyMatrix.h:203
lsstimport.m
tuple m
Definition: lsstimport.py:48
Assert
#define Assert(x)
Definition: dbg.h:73
void lsst::meas::algorithms::shapelet::augmentZTransformCols ( std::complex< double >  z,
int  order,
DMatrix &  T 
)

Definition at line 156 of file BVec.cc.

157  {
158  const int order1 = order;
159  const int order2 = order+2;
160  const int size1 = (order1+1)*(order1+2)/2;
161  //const int size2 = (order2+1)*(order2+2)/2;
162  Assert(int(T.TMV_colsize()) >= size1);
163  Assert(int(T.TMV_rowsize()) >= (order2+1)*(order2+2)/2);
164 
165  if (z == 0.0) return;
166 
167  std::complex<double> zo2 = -z/2.;
168  std::vector<std::vector<std::complex<double> > > f(
169  order2+1,std::vector<std::complex<double> >(order1+1));
170 
171  f[0][0] = exp(-std::norm(z)/8.); // f(0,0)
172  for(int p=1;p<=order2;++p) {
173  f[p][0] = f[p-1][0]*(-zo2)/sqrtn(p); // f(p,0)
174  }
175  for(int s=0;s<order1;++s) {
176  f[0][s+1] = std::conj(zo2)*f[0][s]/sqrtn(s+1); // f(0,s+1)
177  for(int p=1;p<=order2;++p) {
178  f[p][s+1] = (sqrtn(p)*f[p-1][s] + std::conj(zo2)*f[p][s])/
179  sqrtn(s+1); // f(p,s+1)
180  }
181  }
182 
183  for(int n=order1+1,pq=size1;n<=order2;++n) {
184  for(int p=n,q=0;p>=q;--p,++q,++pq) {
185  const std::vector<std::complex<double> >& fp = f[p];
186  const std::vector<std::complex<double> >& fq = f[q];
187  double* Tpq = TMV_ptr(T.col(pq));
188  double* Tpq1 = Tpq + TMV_stepj(T);
189  for(int nn=0,st=0;nn<=order1;++nn) {
190  for(int s=nn,t=0;s>=t;--s,++t,++st) {
191  std::complex<double> t1 = fp[s] * std::conj(fq[t]);
192  if (p==q) {
193  if (s==t) {
194  Tpq[st] = std::real(t1);
195  } else {
196  Tpq[st] = std::real(t1);
197  Tpq[st+1] = std::imag(t1);
198  ++st;
199  }
200  } else if (s==t) {
201  Tpq[st] = 2.*std::real(t1);
202  Tpq1[st] = -2.*std::imag(t1);
203  } else {
204  std::complex<double> t2 = fq[s] * std::conj(fp[t]);
205  Tpq[st] = std::real(t1) + std::real(t2);
206  Tpq[st+1]= std::imag(t1) + std::imag(t2);
207  Tpq1[st] = -std::imag(t1) + std::imag(t2);
208  Tpq1[st+1] = std::real(t1) - std::real(t2);
209  ++st;
210  }
211  }
212  }
213  if (p!=q) ++pq;
214  }
215  }
216  }
lsst::afw::math::details::norm
T norm(const T &x)
Definition: Integrate.h:191
TMV_ptr
#define TMV_ptr(m)
Definition: MyMatrix.h:205
lsst::afw::math::details::real
T real(const T &x)
Definition: Integrate.h:193
Assert
#define Assert(x)
Definition: dbg.h:73
TMV_stepj
#define TMV_stepj(m)
Definition: MyMatrix.h:208
double lsst::meas::algorithms::shapelet::betacf ( double  a,
double  b,
double  x 
)
inline

Definition at line 650 of file Function2D.cc.

651  {
652  double qab = a+b;
653  double qap = a+1.0;
654  double qam = a-1.0;
655  double c = 1.0;
656  double d = 1.0 - qab*x/qap;
657  if (std::abs(d) < FPMIN) d = FPMIN;
658  d = 1.0/d;
659  double h = d;
660  int m;
661  for(m=1;m<=MAXIT;++m) {
662  int m2 = 2*m;
663  double aa = m*(b-m)*x/((qam+m2)*(a+m2));
664  d = 1.0+aa*d;
665  if (std::abs(d) < FPMIN) d = FPMIN;
666  c = 1.0+aa/c;
667  if (std::abs(c) < FPMIN) c = FPMIN;
668  d = 1.0/d;
669  h *= d*c;
670  aa = -(a+m)*(qab+m)*x/((a+m2)*(qap+m2));
671  d = 1.0+aa*d;
672  if (std::abs(d) < FPMIN) d = FPMIN;
673  c = 1.0+aa/c;
674  if (std::abs(c) < FPMIN) c = FPMIN;
675  d = 1.0/d;
676  double del = d*c;
677  h *= del;
678  if (std::abs(del-1.0) < EPS) break;
679  }
680  if (m>MAXIT) {
681  throw std::runtime_error(
682  "a or b too big in betacf, or MAXIT too small");
683  }
684  return h;
685  }
MAXIT
#define MAXIT
Definition: Function2D.cc:646
FPMIN
#define FPMIN
Definition: Function2D.cc:648
x
double x
Definition: ChebyshevBoundedField.cc:305
EPS
#define EPS
Definition: Function2D.cc:647
lsstimport.m
tuple m
Definition: lsstimport.py:48
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
double lsst::meas::algorithms::shapelet::betai ( double  a,
double  b,
double  x 
)
inline

Definition at line 704 of file Function2D.cc.

705  {
706  if (x<0.0 || x>1.0) throw std::runtime_error("Bad x in betai");
707  if (x==0.0) return 0.;
708  if (x==1.0) return 1.;
709  double bt = exp(gammln(a+b)-gammln(a)-gammln(b)+a*log(x)+b*log(1.0-x));
710  if (x < (a+1.0)/(a+b+2.0)) return bt*betacf(a,b,x)/a;
711  else return 1.0-bt*betacf(b,a,1.0-x)/b;
712  }
lsst::meas::algorithms::shapelet::betacf
double betacf(double a, double b, double x)
Definition: Function2D.cc:650
lsst.log.log.log
def log
Definition: log.py:85
x
double x
Definition: ChebyshevBoundedField.cc:305
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
double lsst::meas::algorithms::shapelet::binom ( int  i,
int  j 
)

Definition at line 85 of file BinomFact_omp.cc.

86  {
87  // return iCj, i!/(j!(i-j)!)
88  static std::vector<std::vector<double> > f(10);
89  static bool first=true;
90  if (first) {
91  f[0] = std::vector<double>(1,1.);
92  f[1] = std::vector<double>(2,1.);
93  for(int i1=2;i1<10;i1++) {
94  f[i1] = std::vector<double>(i1+1);
95  f[i1][0] = f[i1][i1] = 1.;
96  for(int j1=1;j1<i1;j1++) f[i1][j1] = f[i1-1][j1-1] + f[i1-1][j1];
97  }
98  first = false;
99  }
100  if (i>=(int)f.size()) {
101 #ifdef _OPENMP
102 #pragma omp critical (binom)
103 #endif
104  {
105  for(int i1=f.size();i1<=i;i1++) {
106  f.push_back(std::vector<double>(i1+1,1.));
107  for(int j1=1;j1<i1;j1++) f[i1][j1] = f[i1-1][j1-1] + f[i1-1][j1];
108  }
109  }
110  assert(i==(int)f.size()-1);
111  }
112  if (j<0 || j>i) return 0.;
113  assert(i<(int)f.size());
114  assert(j<(int)f[i].size());
115  return f[i][j];
116  }
void lsst::meas::algorithms::shapelet::calculateGTransform ( std::complex< double >  g,
int  order1,
int  order2,
DMatrix &  S 
)

Definition at line 468 of file BVec.cc.

470  {
471  const int size1 = (order1+1)*(order1+2)/2;
472  const int size2 = (order2+1)*(order2+2)/2;
473  Assert(int(S.TMV_colsize()) >= size1);
474  Assert(int(S.TMV_rowsize()) >= size2);
475  // S(st,pq) = f(p,s) f(q,t) (eta/|eta|)^(s-t-p+q)
476  // f(p,0) = sqrt(p!)/(p/2)! sqrt(sech(|eta|/2)) (-tanh(|eta|/2)/2)^p/2
477  // f(p,s+1) = (sqrt(p) sech(|eta|/2) f(p-1,s) +
478  // sqrt(s) tanh(|eta|/2) f(p,s-1))/sqrt(s+1)
479  //
480  // tanh(|eta|/2) = |g|
481  // sech(|eta|/2) = sqrt(1-|g|^2)
482  //
483  // Note: Like with the matrix in applyMu, this one is also fairly
484  // sparse. Could get a speedup by expoiting that, but currently don't.
485  // I'll wait until the speedup is found to be necessary.
486 
487  S.setZero();
488 
489  if (g == 0.0) {
490  S.TMV_diagpart(0,0,std::min(size1,size2)).TMV_setAllTo(1.);
491  return;
492  }
493 
494  double absg = std::abs(g);
495  double normg = std::norm(g);
496  std::vector<std::complex<double> > phase(order1+order2+1);
497  phase[0] = 1.;
498  std::complex<double> ph = -std::conj(g)/absg;
499  // I'm not sure why conj was needed here. Maybe there is an
500  // error in the phase indices below that this corrects.
501  for(int i=1;i<=order1+order2;++i) phase[i] = phase[i-1]*ph;
502 
503  double te = absg;
504  double se = sqrt(1.-normg);
505 
506  std::vector<std::vector<double> > f(
507  order2+1,std::vector<double>(order1+1,0.));
508 
509  f[0][0] = sqrt(se); // f(0,0)
510  // only terms with p+s even are non-zero.
511  for(int p=2;p<=order2;p+=2) {
512  f[p][0] = f[p-2][0]*(-te/2.)*sqrtn(p*(p-1))/double(p/2); // f(p,0)
513  }
514 
515  for(int s=0;s<order1;++s) {
516  if (s%2==1) {
517  f[0][s+1] = sqrtn(s)*te*f[0][s-1]/sqrtn(s+1); // f(0,s+1)
518  }
519  for(int p=s%2+1;p<=order2;p+=2) {
520  double temp = sqrtn(p)*se*f[p-1][s];
521  if (s>0) temp += sqrtn(s)*te*f[p][s-1];
522  temp /= sqrtn(s+1);
523  f[p][s+1] = temp; // f(p,s+1)
524  }
525  }
526 
527  for(int n=0,pq=0;n<=order2;++n) {
528  for(int p=n,q=0;p>=q;--p,++q,++pq) {
529  const std::vector<double>& fp = f[p];
530  const std::vector<double>& fq = f[q];
531  double* Spq = TMV_ptr(S.col(pq));
532  double* Spq1 = p>q?TMV_ptr(S.col(pq+1)):0;
533  for(int nn=n%2,st=(nn==0?0:1);nn<=order1;nn+=2,st+=nn) {
534  for(int s=nn,t=0;s>=t;--s,++t,++st) {
535 
536  double s0 = fp[s] * fq[t];
537 
538  int iphase = s-t-p+q;
539  std::complex<double> s1 = s0 *
540  (iphase >= 0 ?
541  phase[iphase/2] :
542  std::conj(phase[-iphase/2]));
543 
544  if (p==q) {
545  if (s==t) {
546  Spq[st] = std::real(s1);
547  } else {
548  Spq[st] = std::real(s1);
549  Spq[st+1] = std::imag(s1);
550  ++st;
551  }
552  } else if (s==t) {
553  // b_st = t_stpq b_pq + t_stqp b_qp
554  // In this case with s=t, t_stpq == t_stqp*
555  Spq[st] = 2.*std::real(s1);
556  Spq1[st] = -2.*std::imag(s1);
557  } else {
558  s0 = fq[s] * fp[t];
559  iphase = s-t-q+p;
560  std::complex<double> s2 = s0 *
561  (iphase >= 0 ?
562  phase[iphase/2] :
563  std::conj(phase[-iphase/2]));
564  Spq[st] = std::real(s1) + std::real(s2);
565  Spq[st+1] = std::imag(s1) + std::imag(s2);
566  Spq1[st] = -std::imag(s1) + std::imag(s2);
567  Spq1[st+1] = std::real(s1) - std::real(s2);
568  ++st;
569  }
570  }
571  }
572  if (p!=q) ++pq;
573  }
574  }
575  }
lsst::afw::math::details::norm
T norm(const T &x)
Definition: Integrate.h:191
TMV_ptr
#define TMV_ptr(m)
Definition: MyMatrix.h:205
lsst::afw::math::details::real
T real(const T &x)
Definition: Integrate.h:193
Assert
#define Assert(x)
Definition: dbg.h:73
void lsst::meas::algorithms::shapelet::calculateGTransform ( std::complex< double >  g,
int  order,
DMatrix &  S 
)
inline

Definition at line 127 of file BVec.h.

128  { calculateGTransform(g,order,order,S); }
lsst::meas::algorithms::shapelet::calculateGTransform
void calculateGTransform(std::complex< double > g, int order1, int order2, DMatrix &S)
Definition: BVec.cc:468
void lsst::meas::algorithms::shapelet::calculateMuTransform ( double  mu,
int  order1,
int  order2,
DMatrix &  D 
)

Definition at line 227 of file BVec.cc.

228  {
229  const int size1 = (order1+1)*(order1+2)/2;
230  const int size2 = (order2+1)*(order2+2)/2;
231  Assert(int(D.TMV_colsize()) >= size1);
232  Assert(int(D.TMV_rowsize()) >= size2);
233  // First I should point out an error in BJ02. It lists
234  // D(pq,00) = e^mu sech(mu) tanh(mu)^p delta_pq
235  // This is wrong.
236  // This is the formula for D(00,pq).
237  //
238  // Second, this formula has a different normalization than what we
239  // want. This gives the transformation I(x,y) -> I(exp(-mu)x,exp(-mu)y).
240  // However, we want the transformation for the same I, but measured in
241  // a different coordinate system. The difference is in the fact that
242  // the flux scales as sigma^2, so this gives another factor of exp(-2mu).
243  // So the formula we need is:
244  // D(00,pq) = e^-mu sech(mu) tanh(mu)^p delta_pq
245  //
246  // Third, the recursion suggested for calculating D(st,pq)
247  // is a bit cumbersome since it has a q+1 on the rhs, so the
248  // calculation of D(30,30) for example, would require knowledge
249  // of D(00,33) which is a higher order than we really need.
250  //
251  // An easier recursion is the following:
252  // D(00,pq) = e^-mu sech(mu) tanh(mu)^p delta_pq
253  // D(s0,pq) = 1/sqrt(s) sech(mu) sqrt(p) D(s-10,p-1q)
254  // D(st,pq) = 1/sqrt(t) [ sech(mu) sqrt(q) D(st-1,pq-1)
255  // - tanh(mu) sqrt(s) D(s-1t-1,pq) ]
256  // Note: Only s-t = p-q terms are nonzero.
257  // This means that this can be significantly sped up
258  // by forming smaller matrices for each m, and using
259  // permutations to rotate b so that these elements are
260  // continuous and then rotate back after multiplying.
261  // However, I'll wait to implement this until such speed
262  // is found to be necessary.
263 
264  D.setZero();
265 
266  if (mu == 0.0) {
267  D.TMV_diagpart(0,0,std::min(size1,size2)).TMV_setAllTo(1.);
268  return;
269  }
270 
271  double tmu = tanh(mu);
272  double smu = 1./cosh(mu);
273 
274  int minorder = std::min(order1,order2);
275  double Dqq = exp(-mu)*smu;
276  // Dqq = exp(-mu) sech(mu) (-tanh(mu))^q
277  // This variable keeps the latest Dqq value:
278  for(int n=0,pq=0;n<=order2;++n) {
279  for(int p=n,q=0;p>=q;--p,++q,++pq) {
280  double* Dpq = TMV_ptr(D.col(pq));
281  double* Dpq_n = TMV_ptr(D.col(pq-n));
282  double* Dpq_n_2 = q>0?TMV_ptr(D.col(pq-n-2)):0;
283  double* Dpq1 = p>q?TMV_ptr(D.col(pq+1)):0;
284  if (p==q) {
285  if (p > 0) Dqq *= tmu;
286  Dpq[0] = Dqq; // D(0,pq)
287  }
288  for(int m=1,st=1;m<=minorder;++m) {
289  for(int s=m,t=0;s>=t;--s,++t,++st) {
290  if (p-q==s-t) {
291  double temp;
292  if (t == 0) {
293  temp = smu*sqrtn(p)*Dpq_n[st-m]/sqrtn(s);
294  } else {
295  temp = -tmu*sqrtn(s)*Dpq[st-2*m-1];
296  if (q > 0) {
297  temp += smu*sqrtn(q)*Dpq_n_2[st-m-2];
298  }
299  temp /= sqrtn(t);
300  }
301  if (s == t) {
302  Dpq[st] = temp; // D(st,pq)
303  } else {
304  Dpq[st] = temp; // D(st,pq)
305  Dpq1[st+1] = temp; // D(st+1,pq+1)
306  }
307  }
308  if (s!=t) ++st;
309  }
310  }
311  if (p!=q) ++pq;
312  }
313  }
314  }
TMV_ptr
#define TMV_ptr(m)
Definition: MyMatrix.h:205
lsstimport.m
tuple m
Definition: lsstimport.py:48
Assert
#define Assert(x)
Definition: dbg.h:73
void lsst::meas::algorithms::shapelet::calculateMuTransform ( double  mu,
int  order,
DMatrix &  D 
)
inline

Definition at line 113 of file BVec.h.

114  { calculateMuTransform(mu,order,order,D); }
lsst::meas::algorithms::shapelet::calculateMuTransform
void calculateMuTransform(double mu, int order1, int order2, DMatrix &D)
Definition: BVec.cc:227
void lsst::meas::algorithms::shapelet::calculatePsfConvolve ( const BVec &  bpsf,
int  order1,
int  order2,
double  sigma,
DMatrix &  C 
)

Definition at line 675 of file BVec.cc.

677  {
678  //xdbg<<"Start calculatePsfConvolve\n";
679  //xdbg<<"bpsf = "<<bpsf<<std::endl;
680  //xdbg<<"order = "<<order<<std::endl;
681  //xdbg<<"sigma = "<<sigma<<std::endl;
682  // Just make the matrix which multiplies b to effect the convolution.
683  // b_obs = C * b_init
684  // However, we do not actually use it this way. Rather, we use it to
685  // switch the model from:
686  // I = Sum psi_pq b_obs_pq
687  // to
688  // I = Sum psi_pq C b_init_pq
689  // We use this to solve for the ML b_init.
690  const int order3 = bpsf.getOrder();
691  Assert(int(C.TMV_colsize()) >= (order1+1)*(order2+2)/2);
692  Assert(int(C.TMV_rowsize()) >= (order2+1)*(order2+2)/2);
693  Assert(int(bpsf.size()) == (order3+1)*(order3+2)/2);
694 
695  // C(st,pq) = 1/sqrt(pi) Sum_uv sqrt(p!u!q!v!/s!t!)/w!
696  // G(s,p,u) G(t,q,v) bpsf_uv
697  // The sum is only over terms for which p+u-s == q+v-t,
698  // and w = p+u-s = q+v-t >= 0
699  //
700  // G(0,p,u) = binom(p+u,u) (-A)^u B^p
701  // G(s+1,p,u) = A G(s,p-1,u) + B G(s,p,u-1)
702  // where A = sigma_init / sigma_obs, B = sigma_psf / sigma_obs
703  // D = A^2 = 1-B^2
704  //
705  // It is more efficient to combine the sqrt(p!u!/s!w!) into the G(s,p,u).
706  // Call the product H(s,p,u). We need to translate the above formulae:
707  // H(0,p,u) = sqrt(p!u!/(p+u)!) (p+u)!/(p!u!) (-A)^u B^p
708  // = sqrt((p+u)!/(p!u!)) (-A)^u B^p
709  // H(s+1,p,u) = sqrt(p!u!/(s+1)!(p+u-s-1)!) * [
710  // A H(s,p-1,u) / sqrt((p-1)!u!/s!(p+u-s-1)!) +
711  // B H(s,p,u-1) / sqrt(p!(u-1)!/s!(p+u-s-1)!) ]
712  // = A sqrt(p)/sqrt(s+1) H(s,p-1,u) +
713  // B sqrt(u)/sqrt(s+1) H(s,p,u-1)
714 
715  C.setZero();
716 
717  // sigma^2 = exp(mu)
718  // D = sigma_i^2 / (sigma_i^2 + sigma_psf^2)
719  // = 1 / (1 + (sigma_psf/sigma_i)^2)
720  double D = 1./(1.+pow(bpsf.getSigma()/sigma,2));
721  double A = sqrt(D);
722  double B = sqrt(1-D);
723 
724  std::vector<std::vector<std::vector<double> > > H(
725  order1+1,std::vector<std::vector<double> >(
726  order2+1,std::vector<double>(order3+1)));
727 
728  H[0][0][0] = 1.; // H[0](0,0)
729  for(int u=0;u<=order3;++u) {
730  if (u>0) H[0][0][u] = -A * H[0][0][u-1];
731  for(int p=1;p<=order2;++p)
732  H[0][p][u] = B*sqrtn(p+u)/sqrtn(p)*H[0][p-1][u]; // H[0](p,u)
733  }
734  for(int s=0;s<order1;++s) {
735  H[s+1][0][0] = 0.;
736  for(int p=1;p<=order2;++p)
737  H[s+1][p][0] = A*sqrtn(p)*H[s][p-1][0]/sqrtn(s+1);
738  for(int u=1;u<=order3;++u) {
739  H[s+1][0][u] = B*sqrtn(u)*H[s][0][u-1]/sqrtn(s+1);
740  for(int p=1;p<=order2;++p)
741  H[s+1][p][u] = (A*sqrtn(p)*H[s][p-1][u] +
742  B*sqrtn(u)*H[s][p][u-1])/ sqrtn(s+1);
743  }
744  }
745 
746  int pq = 0;
747  const double* bpsfv = TMV_cptr(bpsf.vec());
748  for(int n=0;n<=order2;++n) {
749  for(int p=n,q=0;p>=q;(p==q?++pq:pq+=2),--p,++q) {
750  //xdbg<<"Start column "<<pq<<" = "<<p<<','<<q<<std::endl;
751  int pmq = p-q;
752  double* Cpq = TMV_ptr(C.col(pq));
753  double* Cpq1 = p>q?TMV_ptr(C.col(pq+1)):0;
754  int st = 0;
755  for(int nn=0;nn<=order1;++nn) {
756  for(int s=nn,t=0;s>=t;(s==t?++st:st+=2),--s,++t) {
757  //xdbg<<"st = "<<st<<" = "<<s<<','<<t<<std::endl;
758  int smt = s-t;
759  double Cpqst = 0.;
760  double Cpq1st = 0.;
761  double Cpqst1 = 0.;
762  double Cpq1st1 = 0.;
763  const std::vector<double>& Hsp = H[s][p];
764  const std::vector<double>& Hsq = H[s][q];
765  const std::vector<double>& Htp = H[t][p];
766  const std::vector<double>& Htq = H[t][q];
767  int uv0 = 0;
768  int parity = (n+nn)%2;
769  for(int upv=0;upv<=order3;++upv,uv0+=upv) {
770  if (upv % 2 != parity) continue;
771  // There are three values of u-v that are worth considering:
772  // u-v = (s-t) - (p-q) >= 0
773  // u-v = (s-t) + (p-q) > 0
774  // u-v = (p-q) - (s-t) < 0
775 
776  int umv = smt-pmq;
777  if (umv >= 0 && umv <= upv) {
778  // First do terms with p>=q, u>=v (always keep s>=t)
779  // s-t = p-q + u-v
780  int u = (upv+umv)/2;
781  int v = (upv-umv)/2;
782 
783  int w = p+u-s;
784  Assert(q+v-t == w);
785  //Assert((w >= 0) == (umv >= 0));
786  if (w >= 0) {
787  int uv = uv0 + 2*v;
788  if (umv == 0) {
789  double temp = Hsp[u]*Htq[v]*bpsfv[uv];
790  if (s==t) {
791  Assert(p==q);
792  Cpqst += temp;
793  } else {
794  Cpqst += temp;
795  Cpq1st1 += temp;
796  }
797  } else {
798  Assert(s>t);
799  double tempr = Hsp[u]*Htq[v];
800  double tempi = tempr * bpsfv[uv+1];
801  tempr *= bpsfv[uv];
802  if (p==q) {
803  Cpqst += tempr;
804  Cpqst1 += tempi;
805  } else {
806  Assert(p>q);
807  Cpqst += tempr;
808  Cpqst1 += tempi;
809  Cpq1st -= tempi;
810  Cpq1st1 += tempr;
811  }
812  }
813  }
814  }
815 
816  umv = smt+pmq;
817  if (pmq != 0 && umv <= upv) {
818  // Next p<q, u>v. Implement by swapping p,q
819  // These terms account for the fact that
820  // b_init_qp = b_init_pq*
821  // s-t = q-p + u-v
822 
823  Assert(umv > 0);
824  int u = (upv+umv)/2;
825  int v = (upv-umv)/2;
826 
827  int w = q+u-s;
828  Assert(p+v-t == w);
829  //Assert((w >= 0) == (umv > 0));
830  if (w >= 0) {
831  int uv = uv0 + 2*v;
832  //Assert(w > 0);
833  //Assert((w > 0) == (umv > 0));
834  Assert(u>v);
835  double tempr = Hsq[u]*Htp[v];
836  double tempi = tempr * bpsfv[uv+1];
837  tempr *= bpsfv[uv];
838  if (smt==0) {
839  Cpqst += tempr;
840  Cpq1st += tempi;
841  } else {
842  Cpqst += tempr;
843  Cpqst1 += tempi;
844  Cpq1st += tempi;
845  Cpq1st1 -= tempr;
846  }
847  }
848  }
849 
850  umv = pmq-smt;
851  if (umv > 0 && umv <= upv) {
852  // Next p>q, u<v.
853  // These terms account for b_psf_vu = b_psf_uv*
854  int u = (upv+umv)/2;
855  int v = (upv-umv)/2;
856 
857  int w = p+v-s;
858  Assert(q+u-t == w);
859  if (w >= 0) {
860  int uv = uv0 + 2*v;
861  // s-t = p-q + v-u
862  Assert(p>q);
863  double tempr = Hsp[v]*Htq[u];
864  double tempi = -tempr * bpsfv[uv+1];
865  tempr *= bpsfv[uv];
866  if (smt==0) {
867  Cpqst += tempr;
868  Cpq1st -= tempi;
869  } else {
870  Cpqst += tempr;
871  Cpqst1 += tempi;
872  Cpq1st -= tempi;
873  Cpq1st1 += tempr;
874  }
875  }
876  }
877  }
878  Assert(uv0 == int(bpsf.size()));
879  //xdbg<<"Cpq[st] = "<<Cpqst<<std::endl;
880  Cpq[st] = Cpqst;
881  if (smt != 0) Cpq[st+1] = Cpqst1;
882  if (pmq != 0) {
883  Cpq1[st] = Cpq1st;
884  if (smt != 0) Cpq1[st+1] = Cpq1st1;
885  }
886  }
887  }
888  Assert(st == int(C.TMV_colsize()));
889  }
890  }
891  Assert(pq == int(C.TMV_rowsize()));
892  C /= afwGeom::SQRTPI;
893  }
TMV_cptr
#define TMV_cptr(m)
Definition: MyMatrix.h:206
TMV_ptr
#define TMV_ptr(m)
Definition: MyMatrix.h:205
sigma
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:43
w
double w
Definition: CoaddPsf.cc:57
lsst::afw::geom::SQRTPI
double const SQRTPI
Definition: Angle.h:22
Assert
#define Assert(x)
Definition: dbg.h:73
void lsst::meas::algorithms::shapelet::calculatePsfConvolve ( const BVec &  bpsf,
int  order,
double  sigma,
DMatrix &  C 
)
inline

Definition at line 134 of file BVec.h.

136  { calculatePsfConvolve(bpsf,order,order,sigma,C); }
sigma
afw::table::Key< double > sigma
Definition: GaussianPsf.cc:43
lsst::meas::algorithms::shapelet::calculatePsfConvolve
void calculatePsfConvolve(const BVec &bpsf, int order1, int order2, double sigma, DMatrix &C)
Definition: BVec.cc:675
void lsst::meas::algorithms::shapelet::CalculateShift ( std::complex< double >  c1,
std::complex< double >  g1,
std::complex< double >  m1,
std::complex< double >  c2,
std::complex< double >  g2,
std::complex< double >  m2,
std::complex< double > &  c3,
std::complex< double > &  g3,
std::complex< double > &  m3 
)

Definition at line 38 of file Ellipse.cc.

43  {
44  // The first set (c1,g1,m1) effect a transformation:
45  //
46  // z' = exp(-m1)/sqrt(1-|g1|^2) ( z-c1 - g1 (z-c1)* )
47  //
48  // Now execute a second transformation (c2,g2,m2):
49  //
50  // z'' = exp(-m2)/sqrt(1-|g2|^2) *
51  // [ ( exp(-m1)/sqrt(1-|g1|^2) ( z-c1 - g1 (z-c1)* ) - c2 ) -
52  // g2 * ( exp(-m1)/sqrt(1-|g1|^2) ( z-c1 - g1 (z-c1)* ) - c2 )* ]
53  //
54  // We need to figure out what this corresponds to in terms of an
55  // effective:
56  // z'' = exp(-m3)/sqrt(1-|g3|^2) ( z-c3 - g3 (z-c3)* )
57  //
58  // This is pretty complicated, but we can do it in steps.
59  //
60  // First, just take the terms with z or z*:
61  //
62  // z'' = exp(-m2)/sqrt(1-|g2|^2)/sqrt(1-|g1|^2) *
63  // ( exp(-m1) (z - g1 z*) - g2 exp(-m1*) (z* - g1* z) )
64  // = exp(-m2)/sqrt(1-|g2|^2)/sqrt(1-|g1|^2) *
65  // ( (exp(-m1) + exp(-m1*) g1* g2) z -
66  // (exp(-m1) g1 + exp(-m1*) g2) z* )
67  // = exp(-m1-m2)/sqrt(1-|g2|^2)/sqrt(1-|g1|^2)
68  // ( (1 + R g1* g2) z - (g1 + R g2) z* )
69  // where R == exp(2i imag(m1))
70  //
71  // So,
72  //
73  // g3 = (g1 + R g2) / (1 + R g1* g2)
74  // exp(-m3) = exp(-m1-m2) (1+R g1* g2) sqrt(1-|g3|^2)/
75  // sqrt(1-|g1|^2) / sqrt(1-|g2|^2)
76  //
77  // The g terms simplify (after a bit of algebra):
78  // exp(-m3) = exp(-m1-m2) (1+R g1* g2)/|1+R g1* g2|
79  //
80  // Finally, the translation terms are a bit messy.
81  // The translation terms in the equation for z'' are:
82  //
83  // exp(-m2)/sqrt(1-|g2|^2) *
84  // [ ( exp(-m1)/sqrt(1-|g1|^2) ( -c1 - g1 (-c1)* ) - c2 ) -
85  // g2 * ( exp(-m1)/sqrt(1-|g1|^2) ( -c1 - g1 (-c1)* ) - c2 )* ]
86  //
87  // This is some value, which we set equal to:
88  //
89  // exp(-m3)/sqrt(1-|g3|^2) (-c3 - g3(-c3)* )
90  //
91  // So solving for c3 - g3 c3* is straightforward.
92  // c3 - g3 c3* = X
93  // c3* - g3* c3 = X*
94  // g3 c3* - |g3|^2 c3 = g3 X*
95  // c3 - |g3|^2 c3 = X + g3 X*
96  // c3 = (X + g3 X*) / (1-|g3|^2)
97  //
98 
99  dbg<<"Start CalculateShift\n";
100  dbg<<"c1,g1,m1 = "<<c1<<" "<<g1<<" "<<m1<<std::endl;
101  dbg<<"c2,g2,m2 = "<<c2<<" "<<g2<<" "<<m2<<std::endl;
102 
103  std::complex<double> Rg2 = std::polar(1.,2*imag(m1)) * g2;
104  double normg1 = norm(g1);
105  double normg2 = norm(g2);
106  std::complex<double> denom = 1.+conj(g1)*Rg2;
107  g3 = (g1+Rg2) / denom;
108  dbg<<"g3 = "<<g3<<std::endl;
109 
110  //xdbg<<"Miralda-Escude formula gives g3 = "<<addShears(g1,g2)<<std::endl;
111 
112  m3 = m1 + m2 - std::complex<double>(0.,1.)*arg(denom);
113  dbg<<"m3 = "<<m3<<std::endl;
114 
115  std::complex<double> X = -c1 - g1 * conj(-c1);
116  X = exp(-m1)/sqrt(1.-normg1) * X - c2;
117  X = exp(-m2)/sqrt(1.-normg2) * (X - g2 * conj(X));
118  double normg3 = norm(g3);
119  X /= -exp(-m3)/sqrt(1.-normg3);
120  c3 = (X + g3*conj(X)) / (1.-normg3);
121  dbg<<"c3 = "<<c3<<std::endl;
122  }
lsst::afw::math::details::norm
T norm(const T &x)
Definition: Integrate.h:191
dbg
#define dbg
Definition: dbg.h:69
void lsst::meas::algorithms::shapelet::calculateThetaTransform ( double  theta,
int  order1,
int  order2,
DBandMatrix &  R 
)

Definition at line 419 of file BVec.cc.

421  {
422  const int size1 = (order1+1)*(order1+2)/2;
423  const int size2 = (order2+1)*(order2+2)/2;
424  Assert(int(R.TMV_colsize()) >= size1);
425  Assert(int(R.TMV_rowsize()) >= size2);
426 
427  R.setZero();
428 
429  if (theta == 0.0) {
430  R.TMV_diagpart(0,0,std::min(size1,size2)).TMV_setAllTo(1.);
431  return;
432  }
433 
434  int minorder = std::min(order1,order2);
435  std::vector<std::complex<double> > expimt(minorder+1);
436  expimt[0] = 1.;
437  if (minorder > 0) expimt[1] = std::polar(1.,theta);
438  for(int m=2;m<=minorder;++m) expimt[m] = expimt[m-1] * expimt[1];
439 
440  for(int n=0,pq=0;n<=minorder;++n) {
441  for(int p=n,q=0,m=n;p>=q;--p,++q,++pq,m-=2) {
442  if (m==0) {
443  R(pq,pq) = 1.;
444  } else {
445  R(pq,pq) = real(expimt[m]);
446  R(pq,pq+1) = -imag(expimt[m]);
447  R(pq+1,pq) = imag(expimt[m]);
448  R(pq+1,pq+1) = real(expimt[m]);
449  ++pq;
450  }
451  }
452  }
453  }
lsstimport.m
tuple m
Definition: lsstimport.py:48
lsst::afw::math::details::real
T real(const T &x)
Definition: Integrate.h:193
Assert
#define Assert(x)
Definition: dbg.h:73
void lsst::meas::algorithms::shapelet::calculateThetaTransform ( double  theta,
int  order,
DBandMatrix &  R 
)
inline

Definition at line 121 of file BVec.h.

122  { calculateThetaTransform(theta,order,order,R); }
lsst::meas::algorithms::shapelet::calculateThetaTransform
void calculateThetaTransform(double theta, int order1, int order2, DBandMatrix &R)
Definition: BVec.cc:419
void lsst::meas::algorithms::shapelet::calculateZTransform ( std::complex< double >  z,
int  order1,
int  order2,
DMatrix &  T 
)

Definition at line 84 of file BVec.cc.

86  {
87  const int size1 = (order1+1)*(order1+2)/2;
88  const int size2 = (order2+1)*(order2+2)/2;
89  Assert(int(T.TMV_colsize()) >= size1);
90  Assert(int(T.TMV_rowsize()) >= size2);
91 
92  T.setZero();
93 
94  if (z == 0.0) {
95  T.TMV_diagpart(0,0,std::min(size1,size2)).TMV_setAllTo(1.);
96  return;
97  }
98 
99  // T(st,pq) = f(p,s) f*(q,t)
100  // f(p,0) = (-z/2)^p / sqrt(p!) exp(-|z|^2/8)
101  // f(p,s+1) = (sqrt(p) f(p-1,s) + 1/2 z* f(p,s) )/sqrt(s+1)
102 
103  std::complex<double> zo2 = -z/2.;
104  std::vector<std::vector<std::complex<double> > > f(
105  order2+1,std::vector<std::complex<double> >(order1+1));
106 
107  f[0][0] = exp(-std::norm(z)/8.); // f(0,0)
108  for(int p=1;p<=order2;++p) {
109  f[p][0] = f[p-1][0]*(-zo2)/sqrtn(p); // f(p,0)
110  }
111  for(int s=0;s<order1;++s) {
112  f[0][s+1] = std::conj(zo2)*f[0][s]/sqrtn(s+1); // f(0,s+1)
113  for(int p=1;p<=order2;++p) {
114  f[p][s+1] = (sqrtn(p)*f[p-1][s] + std::conj(zo2)*f[p][s])/
115  sqrtn(s+1); // f(p,s+1)
116  }
117  }
118 
119  for(int n=0,pq=0;n<=order2;++n) {
120  for(int p=n,q=0;p>=q;--p,++q,++pq) {
121  double* Tpq = TMV_ptr(T.col(pq));
122  double* Tpq1 = Tpq + TMV_stepj(T);
123  const std::vector<std::complex<double> >& fp = f[p];
124  const std::vector<std::complex<double> >& fq = f[q];
125  for(int nn=0,st=0;nn<=order1;++nn) {
126  for(int s=nn,t=0;s>=t;--s,++t,++st) {
127  std::complex<double> t1 = fp[s] * std::conj(fq[t]);
128  if (p==q) {
129  if (s==t) {
130  Tpq[st] = std::real(t1); // T(st,pq)
131  } else {
132  Tpq[st] = std::real(t1);
133  Tpq[st+1] = std::imag(t1);
134  ++st;
135  }
136  } else if (s==t) {
137  // b_st = t_stpq b_pq + t_stqp b_qp
138  // In this case with s=t, t_stpq == t_stqp*
139  Tpq[st] = 2.*std::real(t1);
140  Tpq1[st] = -2.*std::imag(t1);
141  } else {
142  std::complex<double> t2 = fq[s] * std::conj(fp[t]);
143  Tpq[st] = std::real(t1) + std::real(t2);
144  Tpq[st+1]= std::imag(t1) + std::imag(t2);
145  Tpq1[st] = -std::imag(t1) + std::imag(t2);
146  Tpq1[st+1] = std::real(t1) - std::real(t2);
147  ++st;
148  }
149  }
150  }
151  if (p!=q) ++pq;
152  }
153  }
154  }
lsst::afw::math::details::norm
T norm(const T &x)
Definition: Integrate.h:191
TMV_ptr
#define TMV_ptr(m)
Definition: MyMatrix.h:205
lsst::afw::math::details::real
T real(const T &x)
Definition: Integrate.h:193
Assert
#define Assert(x)
Definition: dbg.h:73
TMV_stepj
#define TMV_stepj(m)
Definition: MyMatrix.h:208
void lsst::meas::algorithms::shapelet::calculateZTransform ( std::complex< double >  z,
int  order,
DMatrix &  T 
)
inline

Definition at line 107 of file BVec.h.

108  { calculateZTransform(z,order,order,T); }
lsst::meas::algorithms::shapelet::calculateZTransform
void calculateZTransform(std::complex< double > z, int order1, int order2, DMatrix &T)
Definition: BVec.cc:84
bool lsst::meas::algorithms::shapelet::Equivalent ( double  chisq1,
double  chisq2,
int  n1,
int  n2,
double  equivProb 
)
inline

Definition at line 584 of file Function2D.cc.

586  {
587  if (chisq2 <= chisq1) return true;
588  // should only happpen when essentially equal but have rounding errors
589  if (chisq1 <= 0.) return (chisq2 <= 0.);
590 
591  Assert(n1 < n2);
592  if (n1 <= 0) return (n2 <= 0);
593  double f = (chisq2-chisq1)/(n2-n1) / (chisq1/n1);
594 
595  double prob = betai(0.5*n2,0.5*n1,n2/(n2+n1*f));
596  // = probability that these chisq would happen by chance if equiv
597  // (technically if underlying chisq2 really smaller or = to chisq1)
598  // so equiv if prob is large, not equiv if prob is small.
599  return (prob > 1.-equivProb);
600  }
lsst::meas::algorithms::shapelet::betai
double betai(double a, double b, double x)
Definition: Function2D.cc:704
Assert
#define Assert(x)
Definition: dbg.h:73
double lsst::meas::algorithms::shapelet::fact ( int  i)

Definition at line 36 of file BinomFact_omp.cc.

37  {
38  // return i!
39  assert(i>=0);
40  static std::vector<double> f(10);
41  static bool first=true;
42 
43  if (first) {
44  f[0] = f[1] = 1.;
45  for(int j=2;j<10;j++) f[j] = f[j-1]*(double)j;
46  first = false;
47  }
48  if (i>=(int)f.size()) {
49 #ifdef _OPENMP
50 #pragma omp critical (fact)
51 #endif
52  {
53  for(int j=f.size();j<=i;j++)
54  f.push_back(f[j-1]*(double)j);
55  }
56  assert(i==(int)f.size()-1);
57  }
58  assert(i<(int)f.size());
59  return f[i];
60  }
int lsst::meas::algorithms::shapelet::fitSize ( const int  xOrder,
const int  yOrder 
)
inline

Definition at line 397 of file Function2D.cc.

398  {
399  int lowOrder = std::min(xOrder,yOrder);
400  int highOrder = std::max(xOrder,yOrder);
401  return (lowOrder+1)*(lowOrder+2)/2 + (lowOrder+1)*(highOrder-lowOrder);
402  }
double lsst::meas::algorithms::shapelet::gammln ( double  x)
inline

Definition at line 691 of file Function2D.cc.

692  {
693  const double cof[6]={76.18009172947146,-86.50532032941677,
694  24.01409824083091,-1.231739572450155,0.1208650973866179e-2,
695  -0.5395239384953e-5};
696  double temp = x+5.5;
697  temp -= (x+0.5)*log(temp);
698  double ser = 1.000000000190015;
699  double y=x;
700  for(int j=0;j<6;++j) ser += cof[j]/(y+=1.0);
701  return -temp+log(2.5066282746310005*ser/x);
702  }
y
int y
Definition: GaussianCentroid.cc:36
lsst.log.log.log
def log
Definition: log.py:85
x
double x
Definition: ChebyshevBoundedField.cc:305
void lsst::meas::algorithms::shapelet::getSubPixList ( PixelList &  pix,
const PixelList &  allPix,
std::complex< double >  cen_offset,
std::complex< double >  shear,
double  aperture,
long &  flag 
)

Definition at line 254 of file Pixel.cc.

258  {
259  // Select a subset of allPix that are within the given aperture
260  const int nTot = allPix.size();
261  xdbg<<"Start GetSubPixList\n";
262  xdbg<<"allPix has "<<nTot<<" objects\n";
263  xdbg<<"new aperture = "<<aperture<<std::endl;
264  xdbg<<"cen_offset = "<<cen_offset<<std::endl;
265  xdbg<<"shear = "<<shear<<std::endl;
266 
267  double normg = norm(shear);
268  double g1 = real(shear);
269  double g2 = imag(shear);
270  double apsq = aperture*aperture;
271 
272  // Do this next loop in two passes. First figure out which
273  // pixels we want to use. Then we can resize pix to the full size
274  // we will need, and go back through and enter the pixels.
275  // This saves us a lot of resizing calls in vector, which are
276  // both slow and can fragment the memory.
277  std::vector<bool> shouldUsePix(nTot,false);
278  int nPix = 0;
279 
280  double peak = 0.;
281  for(int i=0;i<nTot;++i) {
282  std::complex<double> z = allPix[i].getPos() - cen_offset;
283  double u = real(z);
284  double v = imag(z);
285  // (1 + |g|^2) (u^2+v^2) - 2g1 (u^2-v^2) - 2g2 (2uv)
286  // u,v are in arcsec
287  double usq = u*u;
288  double vsq = v*v;
289  double rsq = (1.+normg)*(usq+vsq) - 2.*g1*(usq-vsq) - 4.*g2*u*v;
290  rsq /= (1.-normg);
291  if (rsq <= apsq) {
292  //xdbg<<"u,v = "<<u<<','<<v<<" rsq = "<<rsq<<std::endl;
293  shouldUsePix[i] = true;
294  ++nPix;
295  if (allPix[i].getFlux() > peak) peak = allPix[i].getFlux();
296  }
297  }
298 
299  xdbg<<"npix = "<<nPix<<std::endl;
300  pix.resize(nPix);
301 
302  xdbg<<"pixlist size = "<<nPix<<" = "<<nPix*sizeof(Pixel)<<
303  " bytes = "<<nPix*sizeof(Pixel)/1024.<<" KB\n";
304 
305  int k=0;
306  xdbg<<"Bright pixels are:\n";
307  for(int i=0;i<nTot;++i) if(shouldUsePix[i]) {
308  Pixel p = allPix[i];
309  p.setPos(p.getPos() - cen_offset);
310  pix[k++] = p;
311  if (p.getFlux() > peak / 10.) {
312  xdbg<<p.getPos()<<" "<<p.getFlux()<<std::endl;
313  }
314  }
315  Assert(k == int(pix.size()));
316 
317  if (nPix < 10) flag |= LT10PIX;
318  }
lsst::afw::math::details::norm
T norm(const T &x)
Definition: Integrate.h:191
xdbg
#define xdbg
Definition: dbg.h:70
lsst::afw::math::details::real
T real(const T &x)
Definition: Integrate.h:193
Assert
#define Assert(x)
Definition: dbg.h:73
void lsst::meas::algorithms::shapelet::getSubPixList ( PixelList &  pix,
const PixelList &  allPix,
std::complex< double >  cen_offset,
double  aperture,
long &  flag 
)
inline

Definition at line 126 of file Pixel.h.

129  { getSubPixList(pix,allPix,cen_offset,0.,aperture,flag); }
lsst::meas::algorithms::shapelet::getSubPixList
void getSubPixList(PixelList &pix, const PixelList &allPix, std::complex< double > cen_offset, std::complex< double > shear, double aperture, long &flag)
Definition: Pixel.cc:254
void lsst::meas::algorithms::shapelet::getSubPixList ( PixelList &  pix,
const PixelList &  allPix,
double  aperture,
long &  flag 
)
inline

Definition at line 131 of file Pixel.h.

134  { getSubPixList(pix,allPix,0.,0.,aperture,flag); }
lsst::meas::algorithms::shapelet::getSubPixList
void getSubPixList(PixelList &pix, const PixelList &allPix, std::complex< double > cen_offset, std::complex< double > shear, double aperture, long &flag)
Definition: Pixel.cc:254
void lsst::meas::algorithms::shapelet::makePsi ( DMatrix &  psi,
CDVectorView  z,
int  order,
const DVector *  coeff = 0 
)

Definition at line 226 of file PsiHelper.cc.

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  }
lsst::afw::geom::INVSQRTPI
double const INVSQRTPI
Definition: Angle.h:23
lsst::afw::math::details::norm
T norm(const T &x)
Definition: Integrate.h:191
TMV_cptr
#define TMV_cptr(m)
Definition: MyMatrix.h:206
TMV_ptr
#define TMV_ptr(m)
Definition: MyMatrix.h:205
TMV_colRange
#define TMV_colRange(m, j1, j2)
Definition: MyMatrix.h:203
lsstimport.m
tuple m
Definition: lsstimport.py:48
Assert
#define Assert(x)
Definition: dbg.h:73
void lsst::meas::algorithms::shapelet::makePsi ( DVector &  psi,
std::complex< double >  z,
int  order 
)
std::complex<double> lsst::meas::algorithms::shapelet::operator- ( const std::complex< double > &  z1,
const Position &  p2 
)
inline

Definition at line 70 of file Bounds.h.

72  { return (p2-z1); }
std::complex<double> lsst::meas::algorithms::shapelet::operator/ ( const Position &  p1,
double  x 
)
inline

Definition at line 74 of file Bounds.h.

76  { Position p2 = p1; p2 /= x; return p2; }
x
double x
Definition: ChebyshevBoundedField.cc:305
std::ostream& lsst::meas::algorithms::shapelet::operator<< ( std::ostream &  os,
const Position &  pos 
)
inline

Definition at line 78 of file Bounds.h.

79  { pos.write(os); return os; }
std::ostream& lsst::meas::algorithms::shapelet::operator<< ( std::ostream &  os,
const BVec &  b 
)
inline

Definition at line 96 of file BVec.h.

97  { os << b.getOrder()<<" "<<b.getSigma()<<" "<<b.vec(); return os; }
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
std::ostream& lsst::meas::algorithms::shapelet::operator<< ( std::ostream &  fout,
const Function2D &  f 
)
inline

Definition at line 158 of file Function2D.h.

159  { f.write(fout); return fout; }
template<typename T >
std::ostream& lsst::meas::algorithms::shapelet::operator<< ( std::ostream &  os,
const BoundForm< T > &  bf 
)
inline

Definition at line 171 of file Form.h.

172  {
173  std::ostringstream s;
174  setup(s,bf);
175  s << bf.val;
176  if (bf.f._nTrail>0) s << std::string(bf.f._nTrail,bf.f._trailCh);
177  os << s.str();
178  return os;
179  }
lsst::meas::algorithms::shapelet::setup
void setup(std::ostream &os, const BoundForm< double > &bf)
Definition: Form.h:129
std::ostream& lsst::meas::algorithms::shapelet::operator<< ( std::ostream &  fout,
const Bounds &  b 
)
inline

Definition at line 177 of file Bounds.h.

178  { b.write(fout); return fout;}
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
std::ostream& lsst::meas::algorithms::shapelet::operator<< ( std::ostream &  os,
const Ellipse &  s 
)
inline

Definition at line 186 of file Ellipse.h.

187  { s.write(os); return os; }
std::ostream& lsst::meas::algorithms::shapelet::operator<< ( std::ostream &  os,
const ConfigFile &  cf 
)
inline

Definition at line 341 of file ConfigFile.h.

342  { cf.write(os); return os; }
std::istream& lsst::meas::algorithms::shapelet::operator>> ( std::istream &  os,
Position &  pos 
)
inline

Definition at line 81 of file Bounds.h.

82  { pos.read(os); return os; }
std::istream& lsst::meas::algorithms::shapelet::operator>> ( std::istream &  fin,
std::auto_ptr< Function2D > &  f 
)
inline

Definition at line 161 of file Function2D.h.

163  { f = Function2D::read(fin); return fin; }
std::istream& lsst::meas::algorithms::shapelet::operator>> ( std::istream &  fin,
Bounds &  b 
)
inline

Definition at line 180 of file Bounds.h.

181  { b.read(fin); return fin;}
b
afw::table::Key< double > b
Definition: DoubleGaussianPsf.cc:49
std::istream& lsst::meas::algorithms::shapelet::operator>> ( std::istream &  is,
ConfigFile &  cf 
)
inline

Definition at line 343 of file ConfigFile.h.

344  { cf.read(is); return is; }
void lsst::meas::algorithms::shapelet::PrintFlags ( const std::vector< long > &  flags,
std::ostream &  os 
)

Definition at line 33 of file Params.cc.

34  {
35  const int nObj = flags.size();
36 
37  std::vector<long> nFlagCount(NFLAGS,0);
38  long nNoFlag = 0;
39 
40  for(int i=0;i<nObj;++i) {
41  long flag = flags[i];
42  if (!flag) ++nNoFlag;
43  for(long flagNum = 0; flagNum < NFLAGS; ++flagNum) {
44  if (flag & (1 << flagNum)) ++nFlagCount[flagNum];
45  }
46  }
47  os<<" Total N = "<<nObj<<std::endl;
48  os<<" # with no flags = "<<nNoFlag<<std::endl;
49  for(long flagNum = 0; flagNum < NFLAGS; ++flagNum) {
50  if (nFlagCount[flagNum]) {
51  os<<" # with "<<flagName[flagNum]<<" = "<<
52  nFlagCount[flagNum]<<std::endl;
53  }
54  }
55  }
void lsst::meas::algorithms::shapelet::setup ( std::ostream &  os,
const BoundForm< double > &  bf 
)
inline

Definition at line 129 of file Form.h.

130  { setupFloat(os,bf.f); }
lsst::meas::algorithms::shapelet::setupFloat
void setupFloat(std::ostream &s, const Form &f)
Definition: Form.h:91
void lsst::meas::algorithms::shapelet::setup ( std::ostream &  os,
const BoundForm< long double > &  bf 
)
inline

Definition at line 132 of file Form.h.

133  { setupFloat(os,bf.f); }
lsst::meas::algorithms::shapelet::setupFloat
void setupFloat(std::ostream &s, const Form &f)
Definition: Form.h:91
void lsst::meas::algorithms::shapelet::setup ( std::ostream &  os,
const BoundForm< float > &  bf 
)
inline

Definition at line 135 of file Form.h.

136  { setupFloat(os,bf.f); }
lsst::meas::algorithms::shapelet::setupFloat
void setupFloat(std::ostream &s, const Form &f)
Definition: Form.h:91
void lsst::meas::algorithms::shapelet::setup ( std::ostream &  os,
const BoundForm< std::complex< double > > &  bf 
)
inline

Definition at line 138 of file Form.h.

139  { setupFloat(os,bf.f); }
lsst::meas::algorithms::shapelet::setupFloat
void setupFloat(std::ostream &s, const Form &f)
Definition: Form.h:91
void lsst::meas::algorithms::shapelet::setup ( std::ostream &  os,
const BoundForm< std::complex< long double > > &  bf 
)
inline

Definition at line 141 of file Form.h.

143  { setupFloat(os,bf.f); }
lsst::meas::algorithms::shapelet::setupFloat
void setupFloat(std::ostream &s, const Form &f)
Definition: Form.h:91
void lsst::meas::algorithms::shapelet::setup ( std::ostream &  os,
const BoundForm< std::complex< float > > &  bf 
)
inline

Definition at line 145 of file Form.h.

146  { setupFloat(os,bf.f); }
lsst::meas::algorithms::shapelet::setupFloat
void setupFloat(std::ostream &s, const Form &f)
Definition: Form.h:91
void lsst::meas::algorithms::shapelet::setup ( std::ostream &  os,
const BoundForm< int > &  bf 
)
inline

Definition at line 148 of file Form.h.

149  { setupInt(os,bf.f); }
lsst::meas::algorithms::shapelet::setupInt
void setupInt(std::ostream &s, const Form &f)
Definition: Form.h:112
void lsst::meas::algorithms::shapelet::setup ( std::ostream &  os,
const BoundForm< short > &  bf 
)
inline

Definition at line 151 of file Form.h.

152  { setupInt(os,bf.f); }
lsst::meas::algorithms::shapelet::setupInt
void setupInt(std::ostream &s, const Form &f)
Definition: Form.h:112
void lsst::meas::algorithms::shapelet::setup ( std::ostream &  os,
const BoundForm< long > &  bf 
)
inline

Definition at line 154 of file Form.h.

155  { setupInt(os,bf.f); }
lsst::meas::algorithms::shapelet::setupInt
void setupInt(std::ostream &s, const Form &f)
Definition: Form.h:112
void lsst::meas::algorithms::shapelet::setup ( std::ostream &  os,
const BoundForm< unsigned int > &  bf 
)
inline

Definition at line 157 of file Form.h.

158  { setupInt(os,bf.f); }
lsst::meas::algorithms::shapelet::setupInt
void setupInt(std::ostream &s, const Form &f)
Definition: Form.h:112
void lsst::meas::algorithms::shapelet::setup ( std::ostream &  os,
const BoundForm< unsigned short > &  bf 
)
inline

Definition at line 160 of file Form.h.

161  { setupInt(os,bf.f); }
lsst::meas::algorithms::shapelet::setupInt
void setupInt(std::ostream &s, const Form &f)
Definition: Form.h:112
void lsst::meas::algorithms::shapelet::setup ( std::ostream &  os,
const BoundForm< unsigned long > &  bf 
)
inline

Definition at line 163 of file Form.h.

164  { setupInt(os,bf.f); }
lsst::meas::algorithms::shapelet::setupInt
void setupInt(std::ostream &s, const Form &f)
Definition: Form.h:112
template<typename T >
void lsst::meas::algorithms::shapelet::setup ( std::ostream &  os,
const BoundForm< T > &  bf 
)
inline

Definition at line 167 of file Form.h.

168  { setupFloat(os,bf.f); }
lsst::meas::algorithms::shapelet::setupFloat
void setupFloat(std::ostream &s, const Form &f)
Definition: Form.h:91
void lsst::meas::algorithms::shapelet::setupFloat ( std::ostream &  s,
const Form &  f 
)
inline

Definition at line 91 of file Form.h.

92  {
93  s.precision(f._prc);
94  s.setf(f._fmt,std::ios_base::floatfield);
95  s.setf(f._just,std::ios_base::adjustfield);
96  if (f._wdt) s.width(f._wdt);
97  if (f._newFillCh) s.fill(f._newFillCh);
98  if (f._doUpper && f._fmt == std::ios_base::scientific) {
99  if (f._doUpper>0) s.setf(std::ios_base::uppercase);
100  else s.unsetf(std::ios_base::uppercase);
101  }
102  if (f._doPlus) {
103  if (f._doPlus>0) s.setf(std::ios_base::showpos);
104  else s.unsetf(std::ios_base::showpos);
105  }
106  if (f._doTrail) {
107  if (f._doTrail>0) s.setf(std::ios_base::showpoint);
108  else s.unsetf(std::ios_base::showpoint);
109  }
110  }
void lsst::meas::algorithms::shapelet::setupGg1 ( DMatrix &  Gg1,
int  order1,
int  order2 
)

Definition at line 449 of file PsiHelper.cc.

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  }
Assert
#define Assert(x)
Definition: dbg.h:73
void lsst::meas::algorithms::shapelet::setupGg2 ( DMatrix &  Gg2,
int  order1,
int  order2 
)

Definition at line 500 of file PsiHelper.cc.

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  }
Assert
#define Assert(x)
Definition: dbg.h:73
void lsst::meas::algorithms::shapelet::setupGmu ( DMatrix &  Gmu,
int  order1,
int  order2 
)

Definition at line 551 of file PsiHelper.cc.

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  }
Assert
#define Assert(x)
Definition: dbg.h:73
void lsst::meas::algorithms::shapelet::setupGth ( DMatrix &  Gth,
int  order1,
int  order2 
)

Definition at line 579 of file PsiHelper.cc.

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  }
Assert
#define Assert(x)
Definition: dbg.h:73
void lsst::meas::algorithms::shapelet::setupGx ( DMatrix &  Gx,
int  order1,
int  order2 
)

Definition at line 370 of file PsiHelper.cc.

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  }
Assert
#define Assert(x)
Definition: dbg.h:73
void lsst::meas::algorithms::shapelet::setupGy ( DMatrix &  Gy,
int  order1,
int  order2 
)

Definition at line 409 of file PsiHelper.cc.

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  }
Assert
#define Assert(x)
Definition: dbg.h:73
void lsst::meas::algorithms::shapelet::setupInt ( std::ostream &  s,
const Form &  f 
)
inline

Definition at line 112 of file Form.h.

113  {
114  s.setf(f._just,std::ios_base::adjustfield);
115  s.setf(f._base,std::ios_base::basefield);
116  if (f._wdt) s.width(f._wdt);
117  if (f._newFillCh) s.fill(f._newFillCh);
118  if (f._doUpper && f._base == std::ios_base::hex) {
119  if (f._doUpper>0) s.setf(std::ios_base::uppercase);
120  else s.unsetf(std::ios_base::uppercase);
121  }
122  if (f._doPlus) {
123  if (f._doPlus>0) s.setf(std::ios_base::showpos);
124  else s.unsetf(std::ios_base::showpos);
125  }
126  if (f._base != std::ios_base::dec) s.setf(std::ios_base::showbase);
127  }
double lsst::meas::algorithms::shapelet::sqrtfact ( int  i)

Definition at line 62 of file BinomFact_omp.cc.

63  {
64  // return sqrt(i!)
65  static std::vector<double> f(10);
66  static bool first=true;
67  if (first) {
68  f[0] = f[1] = 1.;
69  for(int j=2;j<10;j++) f[j] = f[j-1]*sqrt((double)j);
70  first = false;
71  }
72  if (i>=(int)f.size()) {
73 #ifdef _OPENMP
74 #pragma omp critical (sqrtfact)
75 #endif
76  {
77  for(int j=f.size();j<=i;j++)
78  f.push_back(f[j-1]*sqrt((double)j));
79  }
80  }
81  assert(i<(int)f.size());
82  return f[i];
83  }
double lsst::meas::algorithms::shapelet::sqrtn ( int  i)

Definition at line 118 of file BinomFact_omp.cc.

119  {
120  // return sqrt(i)
121  static std::vector<double> f(10);
122  static bool first=true;
123  if (first) {
124  for(int j=0;j<10;j++) f[j] = std::sqrt((double)j);
125  first = false;
126  }
127  if (i>=(int)f.size()) {
128 #ifdef _OPENMP
129 #pragma omp critical (sqrtn)
130 #endif
131  {
132  for(int j=f.size();j<=i;j++)
133  f.push_back(std::sqrt((double)j));
134  }
135  }
136  assert(i<(int)f.size());
137  return f[i];
138  }
int lsst::meas::algorithms::shapelet::Status ( ExitCode  code,
const ConfigFile &  params 
)

Definition at line 80 of file Params.cc.

81  {
82  switch (code) {
83  case SUCCESS :
84  return params.read("success_status",2);
85  case FAILURE :
86  return params.read("failure_status",4);
87  case FAILURE_FILE_NOT_FOUND :
88  return params.read("file_not_found_status",5);
89  case FAILURE_PARAMETER_ERROR :
90  return params.read("parameter_error_status",5);
91  case FAILURE_READ_ERROR :
92  return params.read("read_error_status",5);
93  case FAILURE_WRITE_ERROR :
94  return params.read("write_error_status",4);
95  case FAILURE_PROCESSING_ERROR :
96  return params.read("processing_error_status",4);
97  default :
98  return 0;
99  }
100  }
const char * lsst::meas::algorithms::shapelet::Text ( const ExitCode &  code)

Definition at line 58 of file Params.cc.

59  {
60  switch (code) {
61  case SUCCESS :
62  return "SUCCESS";
63  case FAILURE :
64  return "FAILURE";
65  case FAILURE_FILE_NOT_FOUND :
66  return "FAILURE_FILE_NOT_FOUND";
67  case FAILURE_PARAMETER_ERROR :
68  return "FAILURE_PARAMETER_ERROR";
69  case FAILURE_READ_ERROR :
70  return "FAILURE_READ_ERROR";
71  case FAILURE_WRITE_ERROR :
72  return "FAILURE_WRITE_ERROR";
73  case FAILURE_PROCESSING_ERROR :
74  return "FAILURE_PROCESSING_ERROR";
75  default :
76  return "UNKNOWN";
77  }
78  }
Generated on Thu Sep 24 2015 02:29:38 for LSSTApplications by   doxygen 1.8.5