LSSTApplications  20.0.0
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 #ifndef LEAST_SQ_FITTER_2D
26 #define LEAST_SQ_FITTER_2D
27 
28 #include <cstdio>
29 #include <memory>
30 #include <vector>
31 
32 #include "Eigen/Core"
33 #include "Eigen/SVD"
34 
37 
38 namespace lsst {
39 namespace meas {
40 namespace astrom {
41 namespace sip {
42 
64 template <class FittingFunc>
66 public:
68  const std::vector<double> &s, int order);
69 
70  Eigen::MatrixXd getParams();
71  Eigen::MatrixXd getErrors();
72  double valueAt(double x, double y);
74 
75  double getChiSq();
76  double getReducedChiSq();
77 
78 private:
79  void initFunctions();
80 
81  Eigen::MatrixXd expandParams(Eigen::VectorXd const &input) const;
82 
83  double func2d(double x, double y, int term);
84  double func1d(double value, int exponent);
85 
86  std::vector<double> _x, _y, _z, _s;
87  int _order; // Degree of polynomial to fit, e.g 4=> cubic
88  int _nPar; // Number of parameters in fitting eqn, e.g x^2, xy, y^2, x^3,
89  int _nData; // Number of data points, == _x.size()
90 
91  Eigen::JacobiSVD<Eigen::MatrixXd> _svd;
92  Eigen::VectorXd _par;
93 
95 };
96 
97 // The .cc part
98 
108 template <class FittingFunc>
110  const std::vector<double> &z, const std::vector<double> &s,
111  int order)
112  : _x(x), _y(y), _z(z), _s(s), _order(order), _nPar(0), _par(1) {
113  //_nPar, the number of terms to fix (x^2, xy, y^2 etc.) is \Sigma^(order+1) 1
114  _nPar = 0;
115  for (int i = 0; i < order; ++i) {
116  _nPar += i + 1;
117  }
118 
119  // Check input vectors are the same size
120  _nData = _x.size();
121  if (_nData != static_cast<int>(_y.size())) {
122  throw LSST_EXCEPT(pex::exceptions::RuntimeError, "x and y vectors of different lengths");
123  }
124  if (_nData != static_cast<int>(_s.size())) {
125  throw LSST_EXCEPT(pex::exceptions::RuntimeError, "x and s vectors of different lengths");
126  }
127  if (_nData != static_cast<int>(_z.size())) {
128  throw LSST_EXCEPT(pex::exceptions::RuntimeError, "x and z vectors of different lengths");
129  }
130 
131  for (int i = 0; i < _nData; ++i) {
132  if (_s[i] == 0.0) {
133  std::string msg = "Illegal zero value for fit weight encountered.";
135  }
136  }
137 
138  if (_nData < _order) {
139  throw LSST_EXCEPT(pex::exceptions::RuntimeError, "Fewer data points than parameters");
140  }
141 
142  initFunctions();
143 
144  Eigen::MatrixXd design(_nData, _nPar);
145  Eigen::VectorXd rhs(_nData);
146  for (int i = 0; i < _nData; ++i) {
147  rhs[i] = z[i] / s[i];
148  for (int j = 0; j < _nPar; ++j) {
149  design(i, j) = func2d(_x[i], _y[i], j) / _s[i];
150  }
151  }
152  _svd.compute(design, Eigen::ComputeThinU | Eigen::ComputeThinV);
153  _par = _svd.solve(rhs);
154 }
155 
165 template <class FittingFunc>
167  return expandParams(_par);
168 }
169 
171 template <class FittingFunc>
172 Eigen::MatrixXd LeastSqFitter2d<FittingFunc>::expandParams(Eigen::VectorXd const &input) const {
173  Eigen::MatrixXd out = Eigen::MatrixXd::Zero(_order, _order);
174  int count = 0;
175  for (int i = 0; i < _order; ++i) {
176  for (int j = 0; j < _order - i; ++j) {
177  out(i, j) = input[count++];
178  }
179  }
180  return out;
181 }
182 
185 template <class FittingFunc>
187  double chisq = 0;
188  for (int i = 0; i < _nData; ++i) {
189  double val = _z[i] - valueAt(_x[i], _y[i]);
190  val /= _s[i];
191  chisq += pow(val, 2);
192  }
193 
194  return chisq;
195 }
196 
202 template <class FittingFunc>
204  return getChiSq() / (double)(_nData - _nPar);
205 }
206 
208 template <class FittingFunc>
209 double LeastSqFitter2d<FittingFunc>::valueAt(double x, double y) {
210  double val = 0;
211 
212  // Sum the values of the different orders to get the value of the fitted function
213  for (int i = 0; i < _nPar; ++i) {
214  val += _par[i] * func2d(x, y, i);
215  }
216  return val;
217 }
218 
221 template <class FittingFunc>
224  out.reserve(_nData);
225 
226  for (int i = 0; i < _nData; ++i) {
227  out.push_back(_z[i] - valueAt(_x[i], _y[i]));
228  }
229 
230  return out;
231 }
232 
234 template <class FittingFunc>
236  Eigen::ArrayXd variance(_nPar);
237  for (int i = 0; i < _nPar; ++i) {
238  variance[i] = _svd.matrixV().row(i).dot(
239  (_svd.singularValues().array().inverse().square() * _svd.matrixV().col(i).array()).matrix());
240  }
241  return expandParams(variance.sqrt().matrix());
242 }
243 
244 template <class FittingFunc>
246  // Initialise the array of functions. _funcArray[i] is a object of type math::Function1 of order i
247  _funcArray.reserve(_order);
248 
250  coeff.reserve(_order);
251 
252  coeff.push_back(1);
253  for (int i = 0; i < _order; ++i) {
254  std::shared_ptr<FittingFunc> p(new FittingFunc(coeff));
255  _funcArray.push_back(p);
256 
257  coeff[i] = 0;
258  coeff.push_back(1); // coeff now looks like [0,0,...,0,1]
259  }
260 }
261 
262 // The ith term in the fitting polynomial is of the form x^a * y^b. This function figures
263 // out the value of a and b, then calculates the value of the ith term at the given x and y
264 template <class FittingFunc>
265 double LeastSqFitter2d<FittingFunc>::func2d(double x, double y, int term) {
266  int yexp = 0; // y exponent
267  int xexp = 0; // x exponent
268 
269  for (int i = 0; i < term; ++i) {
270  yexp = (yexp + 1) % (_order - xexp);
271  if (yexp == 0) {
272  xexp++;
273  }
274  }
275 
276  double xcomp = func1d(x, xexp); // x component of polynomial
277  double ycomp = func1d(y, yexp); // y component
278 
279 #if 0 // A useful debugging printf statement
280  printf("The %i(th) function: x^%i * y^%i = %.1f * %.1f\n", term, xexp, yexp, xcomp, ycomp);
281 #endif
282 
283  return xcomp * ycomp;
284 }
285 
286 template <class FittingFunc>
287 double LeastSqFitter2d<FittingFunc>::func1d(double value, int exponent) {
288  return (*_funcArray[exponent])(value);
289 }
290 
291 } // namespace sip
292 } // namespace astrom
293 } // namespace meas
294 } // namespace lsst
295 
296 #endif
y
int y
Definition: SpanSet.cc:49
std::string
STL class.
std::shared_ptr< FittingFunc >
lsst::meas::astrom::sip::LeastSqFitter2d::getChiSq
double getChiSq()
Return a measure of the goodness of fit.
Definition: LeastSqFitter2d.h:186
lsst::meas::astrom::sip::LeastSqFitter2d::getReducedChiSq
double getReducedChiSq()
Return a measure of the goodness of fit.
Definition: LeastSqFitter2d.h:203
std::vector::reserve
T reserve(T... args)
std::vector< double >
std::vector::size
T size(T... args)
FunctionLibrary.h
lsst::meas::astrom::sip::LeastSqFitter2d::residuals
std::vector< double > residuals()
Return a vector of residuals of the fit (i.e the difference between the input z values,...
Definition: LeastSqFitter2d.h:222
val
ImageT val
Definition: CR.cc:146
lsst::meas::astrom::sip::LeastSqFitter2d::valueAt
double valueAt(double x, double y)
Return the value of the best fit function at a given position (x,y)
Definition: LeastSqFitter2d.h:209
lsst::meas::astrom::sip::LeastSqFitter2d
Fit an lsst::afw::math::Function1 object to a set of data points in two dimensions.
Definition: LeastSqFitter2d.h:65
std::vector::push_back
T push_back(T... args)
coeff
table::Key< table::Array< double > > coeff
Definition: PsfexPsf.cc:362
lsst::meas::astrom::sip::LeastSqFitter2d::getErrors
Eigen::MatrixXd getErrors()
Companion function to getParams(). Returns uncertainties in the parameters as a matrix.
Definition: LeastSqFitter2d.h:235
z
double z
Definition: Match.cc:44
x
double x
Definition: ChebyshevBoundedField.cc:277
lsst::meas::astrom::sip::LeastSqFitter2d::getParams
Eigen::MatrixXd getParams()
Build up a triangular matrix of the parameters.
Definition: LeastSqFitter2d.h:166
lsst
A base class for image defects.
Definition: imageAlgorithm.dox:1
LSST_EXCEPT
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
Runtime.h
lsst::meas::astrom::sip::LeastSqFitter2d::LeastSqFitter2d
LeastSqFitter2d(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &z, const std::vector< double > &s, int order)
Fit a 2d polynomial to a set of data points z(x, y)
Definition: LeastSqFitter2d.h:109
variance
afw::table::Key< afw::table::Array< VariancePixelT > > variance
Definition: HeavyFootprint.cc:218
exponent
int exponent
Definition: orientation.cc:41
lsst::pex::exceptions::RuntimeError
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104