LSSTApplications  20.0.0
LSSTDataManagementBasePackage
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 
38 namespace lsst {
39 namespace meas {
40 namespace astrom {
41 namespace sip {
42 
63 template <class FittingFunc>
65 public:
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 
78 private:
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 
103 template <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 
138 template <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 
148 template <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 
159 template <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 
172 template <class FittingFunc>
174  FittingFunc f = getBestFitFunction();
175 
176  return f(x);
177 }
178 
181 template <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 
198 template <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 
217 template <class FittingFunc>
219  return getChiSq() / (double)(_nData - _order);
220 }
221 
224 template <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 
240 template <class FittingFunc>
241 double 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
y
int y
Definition: SpanSet.cc:49
lsst::meas::astrom::sip::LeastSqFitter1d::LeastSqFitter1d
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)
Definition: LeastSqFitter1d.h:104
std::shared_ptr< FittingFunc >
std::vector::reserve
T reserve(T... args)
std::vector< double >
std::vector::size
T size(T... args)
lsst::meas::astrom::sip::LeastSqFitter1d::residuals
std::vector< double > residuals()
Return a vector of residuals of the fit (i.e the difference between the input y values,...
Definition: LeastSqFitter1d.h:182
FunctionLibrary.h
lsst::meas::astrom::sip::LeastSqFitter1d::getErrors
Eigen::VectorXd getErrors()
Return the 1 sigma uncertainties in the best fit parameters as an Eigen::Matrix.
Definition: LeastSqFitter1d.h:149
val
ImageT val
Definition: CR.cc:146
lsst::meas::astrom::sip::LeastSqFitter1d
Fit an lsst::afw::math::Function1 object to a set of data points in one dimension.
Definition: LeastSqFitter1d.h:64
lsst::meas::astrom::sip::LeastSqFitter1d::getChiSq
double getChiSq()
Return a measure of the goodness of fit.
Definition: LeastSqFitter1d.h:199
std::vector::push_back
T push_back(T... args)
coeff
table::Key< table::Array< double > > coeff
Definition: PsfexPsf.cc:362
x
double x
Definition: ChebyshevBoundedField.cc:277
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
lsst::meas::astrom::sip::LeastSqFitter1d::valueAt
double valueAt(double x)
Calculate the value of the function at a given point.
Definition: LeastSqFitter1d.h:173
Runtime.h
variance
afw::table::Key< afw::table::Array< VariancePixelT > > variance
Definition: HeavyFootprint.cc:218
exponent
int exponent
Definition: orientation.cc:41
lsst::meas::astrom::sip::LeastSqFitter1d::getParams
Eigen::VectorXd getParams()
Return the best fit parameters as an Eigen::Matrix.
Definition: LeastSqFitter1d.h:139
lsst::meas::astrom::sip::LeastSqFitter1d::getReducedChiSq
double getReducedChiSq()
Return a measure of the goodness of fit.
Definition: LeastSqFitter1d.h:218
lsst::meas::astrom::sip::LeastSqFitter1d::getBestFitFunction
FittingFunc getBestFitFunction()
Return the best fit polynomial as a lsst::afw::math::Function1 object.
Definition: LeastSqFitter1d.h:160
lsst::pex::exceptions::RuntimeError
Reports errors that are due to events beyond the control of the program.
Definition: Runtime.h:104