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
LeastSqFitter1d.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_1D
26#define LEAST_SQ_FITTER_1D
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
63template <class FittingFunc>
65public:
67 int order);
68
69 Eigen::VectorXd getParams();
70 Eigen::VectorXd getErrors();
71 FittingFunc getBestFitFunction();
72 double valueAt(double x);
74
75 double getChiSq();
76 double getReducedChiSq();
77
78private:
79 void initFunctions();
80
81 double func1d(double value, int exponent);
82
83 std::vector<double> _x, _y, _s;
84 int _order; // Degree of polynomial to fit, e.g 4=> cubic
85 int _nData; // Number of data points, == _x.size()
86
87 Eigen::JacobiSVD<Eigen::MatrixXd> _svd;
88 Eigen::VectorXd _par;
89
91};
92
93// The .cc part
94
103template <class FittingFunc>
105 const std::vector<double> &s, int order)
106 : _x(x), _y(y), _s(s), _order(order) {
107 if (order == 0) {
108 throw LSST_EXCEPT(pex::exceptions::RuntimeError, "Fit order must be >= 1");
109 }
110
111 _nData = _x.size();
112 if (_nData != static_cast<int>(_y.size())) {
113 throw LSST_EXCEPT(pex::exceptions::RuntimeError, "x and y vectors of different lengths");
114 }
115 if (_nData != static_cast<int>(_s.size())) {
116 throw LSST_EXCEPT(pex::exceptions::RuntimeError, "x and s vectors of different lengths");
117 }
118
119 if (_nData < _order) {
120 throw LSST_EXCEPT(pex::exceptions::RuntimeError, "Fewer data points than parameters");
121 }
122
123 initFunctions();
124
125 Eigen::MatrixXd design(_nData, _order);
126 Eigen::VectorXd rhs(_nData);
127 for (int i = 0; i < _nData; ++i) {
128 rhs[i] = y[i] / _s[i];
129 for (int j = 0; j < _order; ++j) {
130 design(i, j) = func1d(_x[i], j) / _s[i];
131 }
132 }
133 _svd.compute(design, Eigen::ComputeThinU | Eigen::ComputeThinV);
134 _par = _svd.solve(rhs);
135}
136
138template <class FittingFunc>
140 Eigen::VectorXd vec = Eigen::VectorXd::Zero(_order);
141 for (int i = 0; i < _order; ++i) {
142 vec(i) = _par(i);
143 }
144 return vec;
145}
146
148template <class FittingFunc>
150 Eigen::ArrayXd variance(_order);
151 for (int i = 0; i < _order; ++i) {
152 variance[i] = _svd.matrixV().row(i).dot(
153 (_svd.singularValues().array().inverse().square() * _svd.matrixV().col(i).array()).matrix());
154 }
155 return variance.sqrt().matrix();
156}
157
159template <class FittingFunc>
161 // FittingFunc and LeastSqFitter disagree on the definition of order of a function.
162 // LSF says that a linear function is order 2 (two coefficients), FF says only 1
163 FittingFunc func(_order - 1);
164
165 for (int i = 0; i < _order; ++i) {
166 func.setParameter(i, _par(i));
167 }
168 return func;
169}
170
172template <class FittingFunc>
174 FittingFunc f = getBestFitFunction();
175
176 return f(x);
177}
178
181template <class FittingFunc>
184 out.reserve(_nData);
185
186 FittingFunc f = getBestFitFunction();
187
188 for (int i = 0; i < _nData; ++i) {
189 out.push_back(_y[i] - f(_x[i]));
190 }
191
192 return out;
193}
194
198template <class FittingFunc>
200 FittingFunc f = getBestFitFunction();
201
202 double chisq = 0;
203 for (int i = 0; i < _nData; ++i) {
204 double val = _y[i] - f(_x[i]);
205 val /= _s[i];
206 chisq += pow(val, 2);
207 }
208
209 return chisq;
210}
211
217template <class FittingFunc>
219 return getChiSq() / (double)(_nData - _order);
220}
221
224template <class FittingFunc>
226 _funcArray.reserve(_order);
227
229 coeff.reserve(_order);
230
231 coeff.push_back(1.0);
232 for (int i = 0; i < _order; ++i) {
233 std::shared_ptr<FittingFunc> p(new FittingFunc(coeff));
234 _funcArray.push_back(p);
235 coeff[i] = 0.0;
236 coeff.push_back(1.0); // coeff now looks like [0,0,...,0,1]
237 }
238}
239
240template <class FittingFunc>
241double LeastSqFitter1d<FittingFunc>::func1d(double value, int exponent) {
242 return (*_funcArray[exponent])(value);
243}
244
245} // namespace sip
246} // namespace astrom
247} // namespace meas
248} // namespace lsst
249
250#endif
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
afw::table::Key< afw::table::Array< VariancePixelT > > variance
int y
Definition SpanSet.cc:48
Fit an lsst::afw::math::Function1 object to a set of data points in one dimension.
std::vector< double > residuals()
Return a vector of residuals of the fit (i.e the difference between the input y values,...
double getChiSq()
Return a measure of the goodness of fit.
Eigen::VectorXd getParams()
Return the best fit parameters as an Eigen::Matrix.
double getReducedChiSq()
Return a measure of the goodness of fit.
Eigen::VectorXd getErrors()
Return the 1 sigma uncertainties in the best fit parameters as an Eigen::Matrix.
FittingFunc getBestFitFunction()
Return the best fit polynomial as a lsst::afw::math::Function1 object.
LeastSqFitter1d(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &s, int order)
Fit a 1d polynomial to a set of data points z(x, y)
double valueAt(double x)
Calculate the value of the function at a given point.
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