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
LeastSqFitter2d.h
Go to the documentation of this file.
1 // -*- LSST-C++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010 LSST Corporation.
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 
25 
26 #ifndef LEAST_SQ_FITTER_2D
27 #define LEAST_SQ_FITTER_2D
28 
29 
30 #include <cstdio>
31 #include <vector>
32 
33 #include "boost/shared_ptr.hpp"
34 #include "Eigen/Core"
35 #include "Eigen/SVD"
36 
38 #include "lsst/pex/logging/Trace.h"
40 
41 namespace except = lsst::pex::exceptions;
42 namespace pexLogging = lsst::pex::logging;
43 
44 
45 namespace lsst {
46 namespace meas {
47 namespace astrom {
48 namespace sip {
49 
71 template <class FittingFunc>class LeastSqFitter2d {
72 public:
73  LeastSqFitter2d(const std::vector<double> &x, const std::vector<double> &y, const std::vector<double> &z,
74  const std::vector<double> &s, int order);
75 
76  Eigen::MatrixXd getParams();
77  Eigen::MatrixXd getErrors();
78  double valueAt(double x, double y);
79  std::vector<double> residuals();
80 
81  double getChiSq();
82  double getReducedChiSq();
83 
84 private:
85  void initFunctions();
86 
87  Eigen::MatrixXd expandParams(Eigen::VectorXd const & input) const;
88 
89  double func2d(double x, double y, int term);
90  double func1d(double value, int exponent);
91 
92  std::vector<double> _x, _y, _z, _s;
93  int _order; //Degree of polynomial to fit, e.g 4=> cubic
94  int _nPar; //Number of parameters in fitting eqn, e.g x^2, xy, y^2, x^3,
95  int _nData; //Number of data points, == _x.size()
96 
97  Eigen::JacobiSVD<Eigen::MatrixXd> _svd;
98  Eigen::VectorXd _par;
99 
100  std::vector<boost::shared_ptr<FittingFunc> > _funcArray;
101 
102 
103 };
104 
105 
106 
107 
108 //The .cc part
109 namespace sip = lsst::meas::astrom::sip;
110 namespace math = lsst::afw::math;
111 
121 template<class FittingFunc> LeastSqFitter2d<FittingFunc>::LeastSqFitter2d(const std::vector<double> &x,
122  const std::vector<double> &y, const std::vector<double> &z, const std::vector<double> &s, int order) :
123  _x(x), _y(y), _z(z), _s(s), _order(order), _nPar(0), _par(1) {
124 
125  //_nPar, the number of terms to fix (x^2, xy, y^2 etc.) is \Sigma^(order+1) 1
126  _nPar = 0;
127  for (int i = 0; i < order; ++i) {
128  _nPar += i + 1;
129  }
130 
131  //Check input vectors are the same size
132  _nData = _x.size();
133  if (_nData != static_cast<int>(_y.size())) {
134  throw LSST_EXCEPT(except::RuntimeError, "x and y vectors of different lengths");
135  }
136  if (_nData != static_cast<int>(_s.size())) {
137  throw LSST_EXCEPT(except::RuntimeError, "x and s vectors of different lengths");
138  }
139  if (_nData != static_cast<int>(_z.size())) {
140  throw LSST_EXCEPT(except::RuntimeError, "x and z vectors of different lengths");
141  }
142 
143  for (int i = 0; i < _nData; ++i) {
144  if ( _s[i] == 0.0 ) {
145  std::string msg = "Illegal zero value for fit weight encountered.";
146  throw LSST_EXCEPT(except::RuntimeError, msg);
147  }
148  }
149 
150  if (_nData < _order) {
151  throw LSST_EXCEPT(except::RuntimeError, "Fewer data points than parameters");
152  }
153 
154  initFunctions();
155 
156  Eigen::MatrixXd design(_nData, _nPar);
157  Eigen::VectorXd rhs(_nData);
158  for (int i = 0; i < _nData; ++i) {
159  rhs[i] = z[i] / s[i];
160  for (int j = 0; j < _nPar; ++j) {
161  design(i, j) = func2d(_x[i], _y[i], j) / _s[i];
162  }
163  }
164  _svd.compute(design, Eigen::ComputeThinU | Eigen::ComputeThinV);
165  _par = _svd.solve(rhs);
166 }
167 
168 
178 template<class FittingFunc> Eigen::MatrixXd LeastSqFitter2d<FittingFunc>::getParams() {
179  return expandParams(_par);
180 }
181 
183 template<class FittingFunc>
184 Eigen::MatrixXd LeastSqFitter2d<FittingFunc>::expandParams(Eigen::VectorXd const & input) const {
185  Eigen::MatrixXd out = Eigen::MatrixXd::Zero(_order, _order);
186  int count = 0;
187  for (int i = 0; i < _order; ++i) {
188  for (int j = 0; j < _order - i; ++j) {
189  out(i, j) = input[count++];
190  }
191  }
192  return out;
193 }
194 
197 template<class FittingFunc> double LeastSqFitter2d<FittingFunc>::getChiSq() {
198 
199  double chisq = 0;
200  for (int i = 0; i < _nData; ++i) {
201  double val = _z[i] - valueAt(_x[i], _y[i]);
202  val /= _s[i];
203  chisq += pow(val, 2);
204  }
205 
206  return chisq;
207 }
208 
209 
215 template<class FittingFunc> double LeastSqFitter2d<FittingFunc>::getReducedChiSq() {
216  return getChiSq()/(double) (_nData - _nPar);
217 }
218 
219 
220 
222 template<class FittingFunc> double LeastSqFitter2d<FittingFunc>::valueAt(double x, double y){
223  double val = 0;
224 
225  //Sum the values of the different orders to get the value of the fitted function
226  for (int i = 0; i < _nPar; ++i) {
227  val += _par[i] * func2d(x, y, i);
228  }
229  return val;
230 }
231 
234 template<class FittingFunc> std::vector<double> LeastSqFitter2d<FittingFunc>::residuals(){
235 
236  std::vector<double> out;
237  out.reserve(_nData);
238 
239  for (int i = 0; i < _nData; ++i) {
240  out.push_back(_z[i] - valueAt(_x[i], _y[i]));
241  }
242 
243  return out;
244 }
245 
246 
248 template<class FittingFunc>
250  Eigen::ArrayXd variance(_nPar);
251  for (int i = 0; i < _nPar; ++i) {
252  variance[i] = _svd.matrixV().row(i).dot(
253  (_svd.singularValues().array().inverse().square() * _svd.matrixV().col(i).array()).matrix()
254  );
255  }
256  return expandParams(variance.sqrt().matrix());
257 }
258 
259 
260 template<class FittingFunc> void LeastSqFitter2d<FittingFunc>::initFunctions() {
261  //Initialise the array of functions. _funcArray[i] is a object of type math::Function1 of order i
262  _funcArray.reserve(_order);
263 
264  std::vector<double> coeff;
265  coeff.reserve( _order);
266 
267  coeff.push_back(1);
268  for (int i = 0; i < _order; ++i) {
269  boost::shared_ptr<FittingFunc> p(new FittingFunc(coeff));
270  _funcArray.push_back(p);
271 
272  coeff[i] = 0;
273  coeff.push_back(1); //coeff now looks like [0,0,...,0,1]
274  }
275 }
276 
277 
278 //The ith term in the fitting polynomial is of the form x^a * y^b. This function figures
279 //out the value of a and b, then calculates the value of the ith term at the given x and y
280 template<class FittingFunc> double LeastSqFitter2d<FittingFunc>::func2d(double x, double y, int term) {
281 
282  int yexp = 0; //y exponent
283  int xexp = 0; //x exponent
284 
285  for (int i = 0; i<term; ++i) {
286  yexp = (yexp + 1) % (_order - xexp);
287  if ( yexp == 0){
288  xexp++;
289  }
290  }
291 
292  double xcomp = func1d(x, xexp); //x component of polynomial
293  double ycomp = func1d(y, yexp); //y component
294 
295  #if 0 //A useful debugging printf statement
296  printf("The %i(th) function: x^%i * y^%i = %.1f * %.1f\n", term, xexp, yexp, xcomp, ycomp);
297  #endif
298 
299  return xcomp*ycomp;
300 }
301 
302 template<class FittingFunc> double LeastSqFitter2d<FittingFunc>::func1d(double value, int exponent) {
303  return (*_funcArray[exponent])(value);
304 }
305 
306 
307 
308 }}}}
309 
310 #endif
int y
Fit an lsst::afw::math::Function1 object to a set of data points in two dimensions.
Eigen::MatrixXd expandParams(Eigen::VectorXd const &input) const
Turn a flattened parameter-like vector into a triangular matrix.
double func2d(double x, double y, int term)
double getReducedChiSq()
Return a measure of the goodness of fit. Where is the number of data points, and is the number of ...
definition of the Trace messaging facilities
Eigen::JacobiSVD< Eigen::MatrixXd > _svd
Define a collection of useful Functions.
double getChiSq()
Return a measure of the goodness of fit. .
double valueAt(double x, double y)
Return the value of the best fit function at a given position (x,y)
int x
double chisq
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
double func1d(double value, int exponent)
Eigen::MatrixXd getErrors()
Companion function to getParams(). Returns uncertainties in the parameters as a matrix.
LeastSqFitter2d(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &z, const std::vector< double > &s, int order)
ImageT val
Definition: CR.cc:154
Eigen::JacobiSVD< Eigen::MatrixXd > _svd
std::vector< boost::shared_ptr< FittingFunc > > _funcArray