LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
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
38namespace lsst {
39namespace meas {
40namespace astrom {
41namespace sip {
42
64template <class FittingFunc>
66public:
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
78private:
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
108template <class FittingFunc>
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
165template <class FittingFunc>
167 return expandParams(_par);
168}
169
171template <class FittingFunc>
172Eigen::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
185template <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
202template <class FittingFunc>
204 return getChiSq() / (double)(_nData - _nPar);
205}
206
208template <class FittingFunc>
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
221template <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
234template <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
244template <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
264template <class FittingFunc>
265double 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
286template <class FittingFunc>
287double 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
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
afw::table::Key< afw::table::Array< VariancePixelT > > variance
double z
Definition Match.cc:44
int y
Definition SpanSet.cc:48
Fit an lsst::afw::math::Function1 object to a set of data points in two dimensions.
Eigen::MatrixXd getErrors()
Companion function to getParams(). Returns uncertainties in the parameters as a matrix.
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)
double getReducedChiSq()
Return a measure of the goodness of fit.
std::vector< double > residuals()
Return a vector of residuals of the fit (i.e the difference between the input z values,...
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)
Eigen::MatrixXd getParams()
Build up a triangular matrix of the parameters.
Reports errors that are due to events beyond the control of the program.
Definition Runtime.h:104
T push_back(T... args)
T reserve(T... args)
T size(T... args)
int exponent
ImageT val
Definition CR.cc:146
table::Key< table::Array< double > > coeff
Definition PsfexPsf.cc:362
table::Key< int > order