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
Function2D.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 #include <algorithm>
25 #include <complex>
26 #include <iostream>
27 #include <vector>
28 #include <stdexcept>
29 #include <cmath>
30 
33 
34 namespace lsst {
35 namespace meas {
36 namespace algorithms {
37 namespace shapelet {
38 
39  Constant2D::Constant2D(std::istream& fin) : Function2D()
40  {
41  if(!(fin >> (*_coeffs)(0,0))) throw std::runtime_error("reading constant");
42  }
43 
44  void Constant2D::write(std::ostream& fout) const
45  {
46  int oldPrec = fout.precision(6);
47  std::ios_base::fmtflags oldf =
48  fout.setf(std::ios_base::scientific,std::ios_base::floatfield);
49  fout << "C " << (*_coeffs)(0,0) << std::endl;
50  if (!fout) throw std::runtime_error("writing (constant) function");
51  fout.precision(oldPrec);
52  fout.flags(oldf);
53  }
54 
56  {
57  const Constant2D* crhs = dynamic_cast<const Constant2D*>(&rhs);
58  Assert(crhs);
59  (*_coeffs)(0,0) += (*crhs->_coeffs)(0,0);
60  }
61 
62 
63  void Polynomial2D::setFunction(int xo, int yo, const DVector& fVect)
64  {
65  if(_xOrder != xo || _yOrder != yo) {
66  _xOrder = xo; _yOrder = yo;
67  _coeffs.reset(new DMatrix(xo+1,yo+1));
68  _coeffs->setZero();
69  }
70  int k=0;
71 
72  int maxOrder = std::max(xo,yo);
73  for(int m=0;m<=maxOrder;++m) {
74  int i0 = std::min(xo,m);
75  int len = std::min(yo,i0)+1;
76 #ifdef USE_TMV
77  DVectorView coeffDiag = _coeffs->subVector(i0,m-i0,-1,1,len);
78  tmv::ConstVectorView<double> subf = fVect.subVector(k,k+len);
79  coeffDiag = subf;
80 #else
81  for(int i=0;i<len;++i)
82  (*_coeffs)(i0-i,m-i0+i) = fVect(k+i);
83 #endif
84  k += len;
85  }
86 
87  Assert(k==(int)fVect.size());
88  }
89 
90  Polynomial2D::Polynomial2D(std::istream& fin) :
91  Function2D()
92  {
93  // Order of parameters: (example is for xorder = 2, yorder = 3
94  // xorder(2) yorder(3) a00 a10 a01 a20 a11 a02 a21 a12 a03
95  // where f = a00 + a10 x + a01 y + a20 x^2 + a11 xy + a02 y^2
96  // + a21 x^2 y + a12 xy^2 + a03 y^3
97  // Note that total order doesn't go past the max of xorder and yorder.
98  // Also note that a30 is not listed since xorder is only 2.
99  // Note that aij are complex numbers so each is listed as
100  // real_part imag_part.
101  int xo,yo;
102  if (!(fin >> xo >> yo >> _scale))
103  throw std::runtime_error("reading xorder,yorder,scale");
104  _xOrder = xo;
105  _yOrder = yo;
106  _coeffs.reset(new DMatrix(xo+1,yo+1));
107  _coeffs->setZero();
108  int maxOrder = std::max(xo,yo);
109  for(int m=0;m<=maxOrder;++m) {
110  int i0 = std::min(xo,m);
111  int len = std::min(yo,i0)+1;
112 #ifdef USE_TMV
113  DVectorView coeffDiag = _coeffs->subVector(i0,m-i0,-1,1,len);
114  for(int i=0;i<len;++i) fin >> coeffDiag(i);
115 #else
116  for(int i=0;i<len;++i) fin >> (*_coeffs)(i0-i,m-i0+i);
117 #endif
118  }
119  if (!fin) throw std::runtime_error("reading (polynomial)");
120  }
121 
122  void Polynomial2D::write(std::ostream& fout) const
123  {
124  int oldPrec = fout.precision(6);
125  std::ios_base::fmtflags oldf =
126  fout.setf(std::ios_base::scientific,std::ios_base::floatfield);
127  int maxOrder = std::max(_xOrder,_yOrder);
128  if (maxOrder == 0) {
129  fout << "C " << (*_coeffs)(0,0) << std::endl;
130  } else {
131  fout << "P " << _xOrder << ' ' << _yOrder << ' ' << _scale << ' ';
132  for(int m=0;m<=maxOrder;++m) {
133  int i0 = std::min(_xOrder,m);
134  int len = std::min(_yOrder,i0)+1;
135 #ifdef USE_TMV
136  DVectorView coeffDiag = _coeffs->subVector(i0,m-i0,-1,1,len);
137  for(int i=0;i<len;++i) fout << coeffDiag(i) << ' ';
138 #else
139  for(int i=0;i<len;++i) fout << (*_coeffs)(i0-i,m-i0+i);
140 #endif
141  }
142  }
143  fout << std::endl;
144  if (!fout) throw std::runtime_error("writing (polynomial) function");
145  fout.flags(oldf);
146  fout.precision(oldPrec);
147  }
148 
149  void Polynomial2D::addLinear(double a, double b, double c)
150  {
151  (*_coeffs)(0,0) += a;
152  (*_coeffs)(1,0) += b*_scale;
153  (*_coeffs)(0,1) += c*_scale;
154  }
155 
157  double a, double b, double c, double d, double e, double f)
158  {
159  // F(x,y) = Sum_i,j a(i,j) x^i y^j
160  // F'(x,y) = F(a+bx+cy,d+ex+fy)
161  // = Sum_i,j a(i,j) (a+bx+cy)^i (d+ex+fy)^j
162  // = Sum_i,j a(i,j) (Sum_kl iCk kCl a^i-k (bx)^k-l (cy)^l) *
163  // Sum_mn jCm mCn d^j-m (ex)^m-n (fy)^n
164  int maxOrder = std::max(_xOrder,_yOrder);
165  std::vector<double> scaleToThe(maxOrder+1);
166  scaleToThe[0] = 1.0; scaleToThe[1] = _scale;
167  for(int i=2;i<=maxOrder;++i) scaleToThe[i] = scaleToThe[i-1]*_scale;
168  std::vector<double> aToThe(maxOrder+1);
169  std::vector<double> bToThe(maxOrder+1);
170  std::vector<double> cToThe(maxOrder+1);
171  std::vector<double> dToThe(maxOrder+1);
172  std::vector<double> eToThe(maxOrder+1);
173  std::vector<double> fToThe(maxOrder+1);
174  aToThe[0] = 1.; aToThe[1] = a;
175  bToThe[0] = 1.; bToThe[1] = b;
176  cToThe[0] = 1.; cToThe[1] = c;
177  dToThe[0] = 1.; dToThe[1] = d;
178  eToThe[0] = 1.; eToThe[1] = e;
179  fToThe[0] = 1.; fToThe[1] = f;
180  for(int i=2;i<=maxOrder;++i) aToThe[i] = aToThe[i-1]*a;
181  for(int i=2;i<=maxOrder;++i) bToThe[i] = bToThe[i-1]*b;
182  for(int i=2;i<=maxOrder;++i) cToThe[i] = cToThe[i-1]*c;
183  for(int i=2;i<=maxOrder;++i) dToThe[i] = dToThe[i-1]*d;
184  for(int i=2;i<=maxOrder;++i) eToThe[i] = eToThe[i-1]*e;
185  for(int i=2;i<=maxOrder;++i) fToThe[i] = fToThe[i-1]*f;
186  _xOrder = maxOrder; _yOrder = maxOrder;
187 
188  DMatrix binom(maxOrder+1,maxOrder+1);
189  binom(0,0) = 1.0;
190  for(int n=1;n<=maxOrder;++n) {
191  binom(n,0) = 1.0;
192  binom(n,n) = 1.0;
193  for(int m=1;m<n;++m) {
194  binom(n,m) = binom(n-1,m-1) + binom(n-1,m);
195  }
196  }
197 
198  std::auto_ptr<DMatrix> oldCoeffs = _coeffs;
199  _coeffs.reset(new DMatrix(_xOrder+1,_yOrder+1));
200  _coeffs->setZero();
201  for(int i=0;i<=_xOrder;++i) for(int j=0;j<=_yOrder&&i+j<=maxOrder;++j) {
202  for(int k=0;k<=i;++k) for(int l=0;l<=k;++l) {
203  for(int m=0;m<=j;++m) for(int n=0;n<=m;++n) {
204  (*_coeffs)(k-l+m-n,l+n) +=
205  binom(i,k)*binom(k,l)*binom(j,m)*binom(m,n)*
206  aToThe[i-k]*bToThe[k-l]*cToThe[l]*
207  dToThe[j-m]*eToThe[m-n]*fToThe[n]*
208  (*oldCoeffs)(i,j)/scaleToThe[i+j-k-m];
209  }
210  }
211  }
212  }
213 
215  {
216  const Polynomial2D* prhs = dynamic_cast<const Polynomial2D*>(&rhs);
217  Assert(prhs);
218  Assert(_scale == prhs->_scale);
219  if (_xOrder == prhs->_xOrder && _yOrder == prhs->_yOrder) {
220  *_coeffs += *prhs->_coeffs;
221  } else {
222  int newXOrder = std::max(_xOrder,prhs->_xOrder);
223  int newYOrder = std::max(_yOrder,prhs->_yOrder);
224  std::auto_ptr<DMatrix > newCoeffs(
225  new DMatrix(newXOrder+1,newYOrder+1));
226  newCoeffs->setZero();
227  newCoeffs->TMV_subMatrix(0,_xOrder+1,0,_yOrder+1) = *_coeffs;
228  newCoeffs->TMV_subMatrix(0,prhs->_xOrder+1,0,prhs->_yOrder+1) +=
229  *prhs->_coeffs;
230  _coeffs = newCoeffs;
231  _xOrder = newXOrder;
232  _yOrder = newYOrder;
233  }
234  }
235 
237  const Polynomial2D& f, const Polynomial2D& g)
238  {
239  // h(x,y) = f(x,y) * g(x,y)
240  // = (Sum_ij f_ij x^i y^j) (Sum_mn g_mn x^m y^n)
241  // = Sum_ijmj f_ij g_mn x^(i+m) y^(j+n)
242  Assert(_scale == f._scale);
243  Assert(_scale == g._scale);
244  int newXOrder = f.getXOrder() + g.getXOrder();
245  int newYOrder = f.getYOrder() + g.getYOrder();
246  if (_xOrder != newXOrder || _yOrder != newYOrder) {
247  _xOrder = newXOrder;
248  _yOrder = newYOrder;
249  _coeffs.reset(new DMatrix(_xOrder+1,_yOrder+1));
250  }
251  _coeffs->setZero();
252  for(int i=0;i<=f.getXOrder();++i) for(int j=0;j<=f.getYOrder();++j) {
253  for(int m=0;m<=g.getXOrder();++m) for(int n=0;n<=g.getYOrder();++n) {
254  Assert (i+m <= _xOrder && j+n <= _yOrder);
255  (*_coeffs)(i+m,j+n) += (*f._coeffs)(i,j) * (*g._coeffs)(m,n);
256  }
257  }
258  }
259 
260  std::auto_ptr<Function2D> Polynomial2D::dFdX() const
261  {
262  if (_xOrder == 0) {
263  return std::auto_ptr<Function2D>(new Constant2D());
264  }
265  if (_xOrder == 1 && _yOrder == 0) {
266  return std::auto_ptr<Function2D>(
267  new Constant2D((*_coeffs)(1,0)));
268  }
269 
270  int newXOrder = _xOrder-1;
271  int newYOrder = _xOrder > _yOrder ? _yOrder : _yOrder-1;
272 
273  std::auto_ptr<Polynomial2D> temp(
274  new Polynomial2D(newXOrder,newYOrder));
275 
276  int maxOrder = std::max(newXOrder,newYOrder);
277  for(int i=newXOrder;i>=0;--i) {
278  for(int j=std::min(maxOrder-i,newYOrder);j>=0;--j) {
279  Assert(i+1<=_xOrder);
280  Assert(j<=_yOrder);
281  Assert(i+1+j<=std::max(_xOrder,_yOrder));
282  (*temp->_coeffs)(i,j) = (*_coeffs)(i+1,j)*(i+1.)/_scale;
283  }
284  }
285  return std::auto_ptr<Function2D>(temp);
286  }
287 
288  std::auto_ptr<Function2D> Polynomial2D::dFdY() const
289  {
290  if (_yOrder == 0) {
291  return std::auto_ptr<Function2D>(new Constant2D());
292  }
293  if (_yOrder == 1 && _xOrder == 0) {
294  return std::auto_ptr<Function2D>(
295  new Constant2D((*_coeffs)(0,1)));
296  }
297 
298  int newXOrder = _yOrder > _xOrder ? _xOrder : _xOrder-1;
299  int newYOrder = _yOrder-1;
300 
301  std::auto_ptr<Polynomial2D> temp(
302  new Polynomial2D(newXOrder,newYOrder));
303 
304  int maxOrder = std::max(newXOrder,newYOrder);
305  for(int i=newXOrder;i>=0;--i)
306  for(int j=std::min(maxOrder-i,newYOrder);j>=0;--j) {
307  Assert(i<=_xOrder);
308  Assert(j+1<=_yOrder);
309  Assert(i+j+1<=std::max(_xOrder,_yOrder));
310  (*temp->_coeffs)(i,j) = (*_coeffs)(i,j+1)*(j+1.)/_scale;
311  }
312  return std::auto_ptr<Function2D>(temp);
313  }
314 
315  std::auto_ptr<Function2D> Function2D::conj() const
316  {
317  std::auto_ptr<Function2D> temp = copy();
318  TMV_conjugateSelf(*(temp->_coeffs));
319  return temp;
320  }
321 
322  double Function2D::operator()(double x,double y) const
323  {
324  DVector px = definePX(_xOrder,x);
325  DVector py = definePY(_yOrder,y);
326  double result = EIGEN_ToScalar(EIGEN_Transpose(px) * (*_coeffs) * py);
327  return result;
328  }
329 
330  std::auto_ptr<Function2D> Function2D::read(std::istream& fin)
331  {
332  char fc,tc;
333 
334  fin >> fc >> tc;
335  if (tc != 'D' && tc != 'C') fin.putback(tc);
336  std::auto_ptr<Function2D> ret;
337  switch(fc) {
338  case 'C' : ret.reset(new Constant2D(fin));
339  break;
340  case 'P' : ret.reset(new Polynomial2D(fin));
341  break;
342 #ifdef LEGENDRE2D_H
343  case 'L' : ret.reset(new Legendre2D(fin));
344  break;
345 #endif
346 #ifdef CHEBY2D_H
347  case 'X' : ret.reset(new Cheby2D(fin));
348  break;
349 #endif
350  default: throw std::runtime_error("invalid type");
351  }
352  return ret;
353  }
354 
356  double a, double b, double c,
357  const Function2D& f, const Function2D& g)
358  {
359  if(dynamic_cast<Constant2D*>(this)) {
360  Assert(dynamic_cast<const Constant2D*>(&f));
361  Assert(dynamic_cast<const Constant2D*>(&g));
362  }
363  if(dynamic_cast<Polynomial2D*>(this)) {
364  Assert(dynamic_cast<const Polynomial2D*>(&f));
365  Assert(dynamic_cast<const Polynomial2D*>(&g));
366  }
367 #ifdef LEGENDRE2D_H
368  if(dynamic_cast<Legendre2D*>(this)) {
369  Assert(dynamic_cast<const Legendre2D*>(&f));
370  Assert(dynamic_cast<const Legendre2D*>(&g));
371  Assert(dynamic_cast<const Legendre2D*>(&f)->getBounds() ==
372  dynamic_cast<const Legendre2D*>(&g)->getBounds());
373  }
374 #endif
375 #ifdef CHEBY2D_H
376  if(dynamic_cast<Cheby2D*>(this)) {
377  Assert(dynamic_cast<const Cheby2D*>(&f));
378  Assert(dynamic_cast<const Cheby2D*>(&g));
379  Assert(dynamic_cast<const Cheby2D*>(&f)->getBounds() ==
380  dynamic_cast<const Cheby2D*>(&g)->getBounds());
381  }
382 #endif
383  Assert(f.getXOrder() == g.getXOrder());
384  Assert(f.getYOrder() == g.getYOrder());
385 
386  if (_xOrder != f.getXOrder() || _yOrder != f.getYOrder()) {
387  _xOrder = f.getXOrder();
388  _yOrder = f.getYOrder();
389  _coeffs.reset(new DMatrix(_xOrder+1,_yOrder+1));
390  _coeffs->setZero();
391  } else _coeffs->setZero();
392  for(int i=0;i<=_xOrder;++i) for(int j=0;j<=_yOrder;++j) {
393  (*_coeffs)(i,j) = a + b*f.getCoeffs()(i,j) + c*g.getCoeffs()(i,j);
394  }
395  }
396 
397  inline int fitSize(const int xOrder, const int yOrder)
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  }
403 
405  int xOrder, int yOrder,
406  const std::vector<Position>& pos, const std::vector<double>& vals,
407  const std::vector<bool>& shouldUse, DVector *f,
408  const std::vector<double>* sigList,
409  int *dof, DVector *diff, DMatrix* cov)
410  {
411  // f(x,y) = Sum_pq k_pq px_p py_q
412  // = P . K (where each is a vector in pq)
413  // chisq = Sum_n ((v_n - f(x,y))/s_n)^2
414  // minchisq => 0 = Sum_n (v_n - f(x,y)) P_pq/s_n^2
415  // => [Sum_n P_pq P/s_n^2].K = [Sum_n v_n P_pq/s_n^2]
416  // Using ^ to represent an outer product, this can be written:
417  //
418  // [Sum_n (P/s_n)^(P/s_n)] . K = [Sum_n (v_n/s_n) (P/s_n)]
419  //
420  // Or if P' is a matrix with n index for rows, and pq index for columns,
421  // with each element P'(n,pq) = P_n,pq/s_n
422  // and v' is a vector with v'(n) = v_n/s_n
423  //
424  // Then, this can be written:
425  //
426  // P' K = v'
427  //
428  // The solution to this equation, then, gives our answer for K.
429  //
430  Assert(pos.size() == vals.size());
431  Assert(shouldUse.size() == vals.size());
432  Assert((sigList==0) || (sigList->size() == vals.size()));
433  xdbg<<"Start SimpleFit: size = "<<pos.size()<<std::endl;
434  xdbg<<"order = "<<xOrder<<','<<yOrder<<std::endl;
435 
436  int highOrder = std::max(xOrder,yOrder);
437  int size = fitSize(xOrder,yOrder);
438  xdbg<<"size = "<<size<<std::endl;
439 
440  const int nVals = vals.size();
441 
442  Assert(int(f->size()) == size);
443  Assert(!diff || int(diff->size()) == nVals);
444  Assert(!cov ||
445  (int(cov->TMV_colsize()) == size && int(cov->TMV_rowsize()) == size));
446 
447  int nUse = std::count(shouldUse.begin(),shouldUse.end(),true);
448  xdbg<<"nuse = "<<nUse<<std::endl;
449 
450  DMatrix P(nUse,size);
451  P.setZero();
452  DVector V(nUse);
453 
454  int ii=0;
455  for(int i=0;i<nVals;++i) if (shouldUse[i]) {
456  if (sigList) {
457  Assert((*sigList)[i] > 0.);
458  V(ii) = vals[i]/(*sigList)[i];
459  } else {
460  V(ii) = vals[i];
461  }
462 
463  DVector px = definePX(xOrder,pos[i].getX());
464  DVector py = definePY(yOrder,pos[i].getY());
465  int pq=0;
466  for(int pplusq=0;pplusq<=highOrder;++pplusq) {
467  for(int p=std::min(pplusq,xOrder),q=pplusq-p;
468  q<=std::min(pplusq,yOrder);--p,++q,++pq) {
469  Assert(p<int(px.size()));
470  Assert(q<int(py.size()));
471  Assert(pq<int(P.TMV_rowsize()));
472  P(ii,pq) = px[p]*py[q];
473  }
474  }
475  Assert(pq == int(P.TMV_rowsize()));
476  if (sigList) P.row(ii) /= (*sigList)[i];
477  ++ii;
478  }
479  Assert(ii==nUse);
480 
481  xdbg<<"Done make V,P\n";
482  xdbg<<"V = "<<EIGEN_Transpose(V)<<std::endl;
483  //P.divideUsing(tmv::QR);
484  //P.saveDiv();
485  //*f = V/P;
486  TMV_QR(P);
487  TMV_QR_Solve(P,*f,V);
488  xdbg<<"*f = "<<EIGEN_Transpose(*f)<<std::endl;
489 
490  if (diff) {
491  int k=0;
492  for(int i=0;i<nVals;++i) {
493  if (shouldUse[i]) {
494  (*diff)(i) =
495  V(k) - EIGEN_ToScalar(P.row(k) * (*f));
496  ++k;
497  } else {
498  (*diff)(i) = 0.;
499  }
500  }
501  }
502  if (dof) {
503  *dof = P.TMV_colsize() - P.TMV_rowsize();
504  if (*dof < 0) *dof = 0;
505  }
506  if (cov) {
507  //P.makeInverseATA(*cov);
508  TMV_QR_InverseATA(P,*cov);
509  }
510  xdbg<<"Done simple fit\n";
511  }
512 
514  int order, const std::vector<Position>& pos,
515  const std::vector<double>& vals, const std::vector<bool>& shouldUse,
516  const std::vector<double>* sig,
517  double *chisqOut, int *dofOut, DMatrix* cov)
518  {
519  DVector fVect(fitSize(order,order));
520  if (chisqOut) {
521  DVector diff(vals.size());
522  doSimpleFit(order,order,pos,vals,shouldUse,&fVect,sig,dofOut,&diff,cov);
523  *chisqOut = diff.TMV_normSq();
524  } else {
525  doSimpleFit(order,order,pos,vals,shouldUse,&fVect,sig);
526  }
527  setFunction(order,order,fVect);
528  }
529 
530  inline double absSq(const double& x) { return x*x; }
531 
532  //inline double absSq(const std::complex<double>& x) { return std::norm(x); }
533 
535  int order,double nSig,
536  const std::vector<Position>& pos, const std::vector<double>& vals,
537  std::vector<bool> *shouldUse,
538  const std::vector<double>* sig, double *chisqOut, int *dofOut,
539  DMatrix *cov)
540  {
541  xdbg<<"start outlier fit\n";
542  const int nVals = vals.size();
543  bool isDone=false;
544  DVector fVect(fitSize(order,order));
545  int dof;
546  double chisq=0.;
547  double nSigSq = nSig*nSig;
548  while (!isDone) {
549  DVector diff(nVals);
550  xdbg<<"before dosimple\n";
551  doSimpleFit(order,order,pos,vals,*shouldUse,&fVect,sig,&dof,&diff,cov);
552  xdbg<<"after dosimple\n";
553  // Caclulate chisq, keeping the vector diffsq for later when
554  // looking for outliers
555  chisq = diff.TMV_normSq();
556  // If sigmas are given, then chisq should be 1, since diff is
557  // already normalized by sig_i. But if not, then this effectively
558  // assumes that all the given errors are off by some uniform factor.
559 
560  // Look for outliers, setting isDone = false if any are found
561  isDone = true;
562  if (dof <= 0) break;
563  double thresh=nSigSq*chisq/dof;
564  xdbg<<"chisq = "<<chisq<<", thresh = "<<thresh<<std::endl;
565  for(int i=0;i<nVals;++i) if( (*shouldUse)[i]) {
566  xdbg<<"pos ="<<pos[i]<<", v = "<<vals[i]<<
567  " diff = "<<diff[i]<<std::endl;
568  if (absSq(diff[i]) > thresh) {
569  isDone = false;
570  (*shouldUse)[i] = false;
571  xdbg<<i<<" ";
572  }
573  }
574  if (!isDone) xdbg<<" are outliers\n";
575  }
576  setFunction(order,order,fVect);
577  if (chisqOut) *chisqOut = chisq;
578  if (dofOut) *dofOut = dof;
579  xdbg<<"done outlier fit\n";
580  }
581 
582  inline double betai(double a,double b,double x);
583 
584  inline bool Equivalent(double chisq1,double chisq2, int n1, int n2,
585  double equivProb)
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  }
601 
603  int maxOrder, double equivProb,
604  const std::vector<Position>& pos,const std::vector<double>& vals,
605  const std::vector<bool>& shouldUse, const std::vector<double>* sig,
606  double *chisqOut, int *dofOut, DMatrix* cov)
607  {
608  xdbg<<"Start OrderFit\n";
609  DVector fVectMax(fitSize(maxOrder,maxOrder));
610  DVector diff(vals.size());
611  int dofmax;
612  doSimpleFit(maxOrder,maxOrder,pos,vals,shouldUse,
613  &fVectMax,sig,&dofmax,&diff,cov);
614  double chisqmax = diff.TMV_normSq();
615  xdbg<<"chisq,dof(n="<<maxOrder<<") = "<<chisqmax<<','<<dofmax<<std::endl;
616  int tryOrder;
617  double chisq=-1.;
618  int dof=-1;
619  std::auto_ptr<DVector > fVect(0);
620  for(tryOrder=0;tryOrder<maxOrder;++tryOrder) {
621  fVect.reset(new DVector(fitSize(tryOrder,tryOrder)));
622  doSimpleFit(tryOrder,tryOrder,pos,vals,shouldUse,
623  fVect.get(),sig,&dof,&diff,cov);
624  chisq = diff.TMV_normSq();
625  xdbg<<"chisq,dof(n="<<tryOrder<<") = "<<chisq<<','<<dof<<".... ";
626  if (Equivalent(chisqmax,chisq,dofmax,dof,equivProb)) {
627  xdbg<<"equiv\n";
628  break;
629  }
630  xdbg<<"not equiv\n";
631  }
632  if (tryOrder == maxOrder) {
633  setFunction(tryOrder,tryOrder,fVectMax);
634  if(chisqOut) *chisqOut = chisqmax;
635  if(dofOut) *dofOut = dofmax;
636  } else {
637  Assert(fVect.get());
638  setFunction(tryOrder,tryOrder,*fVect);
639  if(chisqOut) *chisqOut = chisq;
640  if(dofOut) *dofOut = dof;
641  }
642  }
643 
644  // These are numerical recipes routines:
645 
646 #define MAXIT 100
647 #define EPS 3.0e-7
648 #define FPMIN 1.0e-30
649 
650  inline double betacf(double a, double b, double x)
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  }
686 
687 #undef MAXIT
688 #undef EPS
689 #undef FPMIN
690 
691  inline double gammln(double x)
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  }
703 
704  inline double betai(double a,double b,double x)
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  }
713 
714 }}}}
int y
Eigen::Block< DVector, Eigen::Dynamic, 1 > DVectorView
Definition: MyMatrix.h:130
#define MAXIT
Definition: Function2D.cc:646
virtual DVector definePY(int order, double y) const =0
virtual DVector definePX(int order, double x) const =0
#define TMV_conjugateSelf(m)
Definition: MyMatrix.h:210
virtual std::auto_ptr< Function2D > conj() const
Definition: Function2D.cc:315
#define TMV_QR_Solve(m, x, b)
Definition: MyMatrix.h:227
#define EIGEN_ToScalar(m)
Definition: MyMatrix.h:212
virtual void setFunction(int _xOrder, int _yOrder, const DVector &fVect)=0
double betacf(double a, double b, double x)
Definition: Function2D.cc:650
virtual void outlierFit(int order, double nsig, const std::vector< Position > &pos, const std::vector< double > &v, std::vector< bool > *shouldUse, const std::vector< double > *sigList=0, double *chisqout=0, int *dofout=0, DMatrix *cov=0)
Definition: Function2D.cc:534
#define FPMIN
Definition: Function2D.cc:648
virtual void write(std::ostream &fout) const
Definition: Function2D.cc:44
def log
Definition: log.py:85
virtual void linearPreTransform(double a, double b, double c, double d, double e, double f)
Definition: Function2D.cc:156
int fitSize(const int xOrder, const int yOrder)
Definition: Function2D.cc:397
virtual std::auto_ptr< Function2D > dFdY() const
Definition: Function2D.cc:288
double betai(double a, double b, double x)
Definition: Function2D.cc:704
#define TMV_QR(m)
Definition: MyMatrix.h:222
virtual void operator+=(const Function2D &rhs)
Definition: Function2D.cc:214
virtual void write(std::ostream &fout) const
Definition: Function2D.cc:122
virtual std::auto_ptr< Function2D > dFdX() const
Definition: Function2D.cc:260
#define EPS
Definition: Function2D.cc:647
virtual void addLinear(double a, double b, double c)
Definition: Function2D.cc:149
int x
double chisq
#define TMV_QR_InverseATA(m, cov)
Definition: MyMatrix.h:229
void makeProductOf(const Polynomial2D &f, const Polynomial2D &g)
Definition: Function2D.cc:236
tuple m
Definition: lsstimport.py:48
virtual void linearTransform(double a, double b, double c, const Function2D &f, const Function2D &g)
Definition: Function2D.cc:355
virtual double operator()(double x, double y) const
Definition: Function2D.cc:322
afw::table::Key< double > b
virtual void simpleFit(int order, const std::vector< Position > &pos, const std::vector< double > &v, const std::vector< bool > &shouldUse, const std::vector< double > *sigList=0, double *chisqOut=0, int *dofOut=0, DMatrix *cov=0)
Definition: Function2D.cc:513
virtual void setFunction(int xOrder, int yOrder, const DVector &fVect)
Definition: Function2D.cc:63
void doSimpleFit(int xOrder, int yOrder, const std::vector< Position > &pos, const std::vector< double > &v, const std::vector< bool > &shouldUse, DVector *f, const std::vector< double > *sigList=0, int *dof=0, DVector *diff=0, DMatrix *cov=0)
Definition: Function2D.cc:404
#define xdbg
Definition: dbg.h:70
virtual std::auto_ptr< Function2D > copy() const =0
double absSq(const double &x)
Definition: Function2D.cc:530
virtual void operator+=(const Function2D &rhs)
Definition: Function2D.cc:55
#define EIGEN_Transpose(m)
Definition: MyMatrix.h:211
#define Assert(x)
Definition: dbg.h:73
static std::auto_ptr< Function2D > read(std::istream &fin)
Definition: Function2D.cc:330
virtual void orderFit(int maxOrder, double equivProb, const std::vector< Position > &pos, const std::vector< double > &v, const std::vector< bool > &shouldUse, const std::vector< double > *sigList=0, double *chisqout=0, int *dofout=0, DMatrix *cov=0)
Definition: Function2D.cc:602
bool Equivalent(double chisq1, double chisq2, int n1, int n2, double equivProb)
Definition: Function2D.cc:584