25 #ifndef LEAST_SQ_FITTER_2D
26 #define LEAST_SQ_FITTER_2D
64 template <
class FittingFunc>
81 Eigen::MatrixXd expandParams(Eigen::VectorXd
const &input)
const;
83 double func2d(
double x,
double y,
int term);
84 double func1d(
double value,
int exponent);
91 Eigen::JacobiSVD<Eigen::MatrixXd> _svd;
108 template <
class FittingFunc>
112 : _x(
x), _y(
y), _z(
z), _s(s), _order(order), _nPar(0), _par(1) {
115 for (
int i = 0; i < order; ++i) {
121 if (_nData !=
static_cast<int>(_y.
size())) {
124 if (_nData !=
static_cast<int>(_s.
size())) {
127 if (_nData !=
static_cast<int>(_z.
size())) {
131 for (
int i = 0; i < _nData; ++i) {
133 std::string msg =
"Illegal zero value for fit weight encountered.";
138 if (_nData < _order) {
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];
152 _svd.compute(design, Eigen::ComputeThinU | Eigen::ComputeThinV);
153 _par = _svd.solve(rhs);
165 template <
class FittingFunc>
167 return expandParams(_par);
171 template <
class FittingFunc>
173 Eigen::MatrixXd out = Eigen::MatrixXd::Zero(_order, _order);
175 for (
int i = 0; i < _order; ++i) {
176 for (
int j = 0; j < _order - i; ++j) {
177 out(i, j) = input[count++];
185 template <
class FittingFunc>
188 for (
int i = 0; i < _nData; ++i) {
189 double val = _z[i] - valueAt(_x[i], _y[i]);
191 chisq += pow(
val, 2);
202 template <
class FittingFunc>
204 return getChiSq() / (double)(_nData - _nPar);
208 template <
class FittingFunc>
213 for (
int i = 0; i < _nPar; ++i) {
214 val += _par[i] * func2d(
x,
y, i);
221 template <
class FittingFunc>
226 for (
int i = 0; i < _nData; ++i) {
227 out.
push_back(_z[i] - valueAt(_x[i], _y[i]));
234 template <
class FittingFunc>
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());
241 return expandParams(
variance.sqrt().matrix());
244 template <
class FittingFunc>
247 _funcArray.reserve(_order);
250 coeff.reserve(_order);
253 for (
int i = 0; i < _order; ++i) {
255 _funcArray.push_back(p);
264 template <
class FittingFunc>
265 double LeastSqFitter2d<FittingFunc>::func2d(
double x,
double y,
int term) {
269 for (
int i = 0; i < term; ++i) {
270 yexp = (yexp + 1) % (_order - xexp);
276 double xcomp = func1d(
x, xexp);
277 double ycomp = func1d(
y, yexp);
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);
283 return xcomp * ycomp;
286 template <
class FittingFunc>
287 double LeastSqFitter2d<FittingFunc>::func1d(
double value,
int exponent) {
288 return (*_funcArray[
exponent])(value);