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
minimize.cc
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/*
26 * Definition of member functions for minuit adaptors
27 *
28 * This file is meant to be included by lsst/afw/math/minimize.h
29 */
30
31#include <string>
32#include <vector>
33#include <memory>
34
35#include "Minuit2/FunctionMinimum.h"
36#include "Minuit2/MnMigrad.h"
37#include "Minuit2/MnMinos.h"
38
39#include "lsst/log/Log.h"
41
42namespace lsst {
43namespace afw {
44namespace math {
45
46namespace {
47/*
48 * Minuit wrapper for a function(x)
49 */
50template <typename ReturnT>
51class MinimizerFunctionBase1 : public ROOT::Minuit2::FCNBase {
52public:
53 explicit MinimizerFunctionBase1(Function1<ReturnT> const &function,
54 std::vector<double> const &measurementList,
55 std::vector<double> const &varianceList,
56 std::vector<double> const &xPositionList, double errorDef);
57 MinimizerFunctionBase1(MinimizerFunctionBase1 const &) = default;
58 MinimizerFunctionBase1(MinimizerFunctionBase1 &&) = default;
59 MinimizerFunctionBase1 &operator=(MinimizerFunctionBase1 const &) = default;
60 MinimizerFunctionBase1 &operator=(MinimizerFunctionBase1 &&) = default;
61 ~MinimizerFunctionBase1() override = default;
62 // Required by ROOT::Minuit2::FCNBase
63 double Up() const override { return _errorDef; }
64 double operator()(const std::vector<double> &) const override;
65
66#if 0 // not used
67 inline std::vector<double> getMeasurements() const {return _measurementList;}
68 inline std::vector<double> getVariances() const {return _varianceList;}
69 inline std::vector<double> getPositions() const {return _xPositionList;}
70 inline void setErrorDef(double def) {_errorDef=def;}
71#endif
72private:
74 std::vector<double> _measurementList;
75 std::vector<double> _varianceList;
76 std::vector<double> _xPositionList;
77 double _errorDef;
78};
79
80/*
81 * Minuit wrapper for a function(x, y)
82 */
83template <typename ReturnT>
84class MinimizerFunctionBase2 : public ROOT::Minuit2::FCNBase {
85public:
86 explicit MinimizerFunctionBase2(Function2<ReturnT> const &function,
87 std::vector<double> const &measurementList,
88 std::vector<double> const &varianceList,
89 std::vector<double> const &xPositionList,
90 std::vector<double> const &yPositionList, double errorDef);
91 MinimizerFunctionBase2(MinimizerFunctionBase2 const &) = default;
92 MinimizerFunctionBase2(MinimizerFunctionBase2 &&) = default;
93 MinimizerFunctionBase2 &operator=(MinimizerFunctionBase2 const &) = default;
94 MinimizerFunctionBase2 &operator=(MinimizerFunctionBase2 &&) = default;
95 ~MinimizerFunctionBase2() override = default;
96 // Required by ROOT::Minuit2::FCNBase
97 double Up() const override { return _errorDef; }
98 double operator()(const std::vector<double> &par) const override;
99
100#if 0 // not used
101 inline std::vector<double> getMeasurements() const {return _measurementList;}
102 inline std::vector<double> getVariances() const {return _varianceList;}
103 inline std::vector<double> getPosition1() const {return _xPositionList;}
104 inline std::vector<double> getPosition2() const {return _yPositionList;}
105 inline void setErrorDef(double def) {_errorDef=def;}
106#endif
107private:
109 std::vector<double> _measurementList;
110 std::vector<double> _varianceList;
111 std::vector<double> _xPositionList;
112 std::vector<double> _yPositionList;
113 double _errorDef;
114};
115} // namespace
116
118template <typename ReturnT>
119MinimizerFunctionBase1<ReturnT>::MinimizerFunctionBase1(Function1<ReturnT> const &function,
120 std::vector<double> const &measurementList,
121 std::vector<double> const &varianceList,
122 std::vector<double> const &xPositionList,
123 double errorDef)
124 : _functionPtr(function.clone()),
125 _measurementList(measurementList),
126 _varianceList(varianceList),
127 _xPositionList(xPositionList),
128 _errorDef(errorDef) {}
129
130template <typename ReturnT>
131MinimizerFunctionBase2<ReturnT>::MinimizerFunctionBase2(Function2<ReturnT> const &function,
132 std::vector<double> const &measurementList,
133 std::vector<double> const &varianceList,
134 std::vector<double> const &xPositionList,
135 std::vector<double> const &yPositionList,
136 double errorDef)
137 : _functionPtr(function.clone()),
138 _measurementList(measurementList),
139 _varianceList(varianceList),
140 _xPositionList(xPositionList),
141 _yPositionList(yPositionList),
142 _errorDef(errorDef) {}
143
144// Only method we need to set up; basically this is a chi^2 routine
145template <typename ReturnT>
146double MinimizerFunctionBase1<ReturnT>::operator()(const std::vector<double> &par) const {
147 // Initialize the function with the fit parameters
148 this->_functionPtr->setParameters(par);
149
150 double chi2 = 0.0;
151 for (unsigned int i = 0; i < this->_measurementList.size(); i++) {
152 double resid = (*(this->_functionPtr))(this->_xPositionList[i]) - this->_measurementList[i];
153 chi2 += resid * resid / this->_varianceList[i];
154 }
155
156 return chi2;
157}
158
159template <typename ReturnT>
160double MinimizerFunctionBase2<ReturnT>::operator()(const std::vector<double> &par) const {
161 // Initialize the function with the fit parameters
162 this->_functionPtr->setParameters(par);
163
164 double chi2 = 0.0;
165 for (unsigned int i = 0; i < this->_measurementList.size(); i++) {
166 double resid = (*(this->_functionPtr))(this->_xPositionList[i], this->_yPositionList[i]) -
167 this->_measurementList[i];
168 chi2 += resid * resid / this->_varianceList[i];
169 }
170
171 return chi2;
172}
174
175template <typename ReturnT>
176FitResults minimize(Function1<ReturnT> const &function, std::vector<double> const &initialParameterList,
177 std::vector<double> const &stepSizeList, std::vector<double> const &measurementList,
178 std::vector<double> const &varianceList, std::vector<double> const &xPositionList,
179 double errorDef) {
180 unsigned int const nParameters = function.getNParameters();
181 if (initialParameterList.size() != nParameters) {
182 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "initialParameterList is the wrong length");
183 }
184 if (stepSizeList.size() != nParameters) {
185 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "stepSizeList is the wrong length");
186 }
187 unsigned int const nMeasurements = measurementList.size();
188 if (varianceList.size() != nMeasurements) {
189 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "varianceList is the wrong length");
190 }
191 if (xPositionList.size() != nMeasurements) {
192 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "xPositionList is the wrong length");
193 }
194
195 MinimizerFunctionBase1<ReturnT> minimizerFunc(function, measurementList, varianceList, xPositionList,
196 errorDef);
197
198 ROOT::Minuit2::MnUserParameters fitPar;
199 std::vector<std::string> paramNames;
200 for (unsigned int i = 0; i < nParameters; ++i) {
201 paramNames.push_back((boost::format("p%d") % i).str());
202 fitPar.Add(paramNames[i].c_str(), initialParameterList[i], stepSizeList[i]);
203 }
204
205 ROOT::Minuit2::MnMigrad migrad(minimizerFunc, fitPar);
206 ROOT::Minuit2::FunctionMinimum min = migrad();
207 ROOT::Minuit2::MnMinos minos(minimizerFunc, min);
208
209 FitResults fitResults;
210 fitResults.chiSq = min.Fval();
211 fitResults.isValid = min.IsValid() && std::isfinite(fitResults.chiSq);
212 if (!fitResults.isValid) {
213 LOGL_WARN("lsst.afw.math.minimize", "Fit failed to converge");
214 }
215
216 for (unsigned int i = 0; i < nParameters; ++i) {
217 fitResults.parameterList.push_back(min.UserState().Value(paramNames[i].c_str()));
218 if (fitResults.isValid) {
219 fitResults.parameterErrorList.push_back(minos(i));
220 } else {
221 double e = min.UserState().Error(paramNames[i].c_str());
222 std::pair<double, double> ep(-1 * e, e);
223 fitResults.parameterErrorList.push_back(ep);
224 }
225 }
226 return fitResults;
227}
228
229template <typename ReturnT>
230FitResults minimize(Function2<ReturnT> const &function, std::vector<double> const &initialParameterList,
231 std::vector<double> const &stepSizeList, std::vector<double> const &measurementList,
232 std::vector<double> const &varianceList, std::vector<double> const &xPositionList,
233 std::vector<double> const &yPositionList, double errorDef) {
234 unsigned int const nParameters = function.getNParameters();
235 if (initialParameterList.size() != nParameters) {
236 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "initialParameterList is the wrong length");
237 }
238 if (stepSizeList.size() != nParameters) {
239 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "stepSizeList is the wrong length");
240 }
241 unsigned int const nMeasurements = measurementList.size();
242 if (varianceList.size() != nMeasurements) {
243 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "varianceList is the wrong length");
244 }
245 if (xPositionList.size() != nMeasurements) {
246 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "xPositionList is the wrong length");
247 }
248 if (yPositionList.size() != nMeasurements) {
249 throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "yPositionList is the wrong length");
250 }
251
252 MinimizerFunctionBase2<ReturnT> minimizerFunc(function, measurementList, varianceList, xPositionList,
253 yPositionList, errorDef);
254
255 ROOT::Minuit2::MnUserParameters fitPar;
256 std::vector<std::string> paramNames;
257 for (unsigned int i = 0; i < nParameters; ++i) {
258 paramNames.push_back((boost::format("p%d") % i).str());
259 fitPar.Add(paramNames[i].c_str(), initialParameterList[i], stepSizeList[i]);
260 }
261
262 ROOT::Minuit2::MnMigrad migrad(minimizerFunc, fitPar);
263 ROOT::Minuit2::FunctionMinimum min = migrad();
264 ROOT::Minuit2::MnMinos minos(minimizerFunc, min);
265
266 FitResults fitResults;
267 fitResults.chiSq = min.Fval();
268 fitResults.isValid = min.IsValid() && std::isfinite(fitResults.chiSq);
269 if (!fitResults.isValid) {
270 LOGL_WARN("lsst.afw.math.minimize", "Fit failed to converge");
271 }
272
273 for (unsigned int i = 0; i < nParameters; ++i) {
274 fitResults.parameterList.push_back(min.UserState().Value(paramNames[i].c_str()));
275 if (fitResults.isValid) {
276 fitResults.parameterErrorList.push_back(minos(i));
277 } else {
278 double e = min.UserState().Error(paramNames[i].c_str());
279 std::pair<double, double> ep(-1 * e, e);
280 fitResults.parameterErrorList.push_back(ep);
281 }
282 }
283 return fitResults;
284}
285
286// Explicit instantiation
288#define NL /* */
289#define minimizeFuncs(ReturnT) \
290 template FitResults minimize(Function1<ReturnT> const &, std::vector<double> const &, \
291 std::vector<double> const &, std::vector<double> const &, \
292 std::vector<double> const &, std::vector<double> const &, double); \
293 NL template FitResults minimize(Function2<ReturnT> const &, std::vector<double> const &, \
294 std::vector<double> const &, std::vector<double> const &, \
295 std::vector<double> const &, std::vector<double> const &, \
296 std::vector<double> const &, double);
297
298minimizeFuncs(float) minimizeFuncs(double)
300} // namespace math
301} // namespace afw
302} // namespace lsst
int min
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
LSST DM logging module built on log4cxx.
#define LOGL_WARN(logger, message...)
Log a warn-level message using a varargs/printf style interface.
Definition Log.h:547
A Function taking one argument.
Definition Function.h:202
A Function taking two arguments.
Definition Function.h:259
unsigned int getNParameters() const noexcept
Return the number of function parameters.
Definition Function.h:112
Reports invalid arguments.
Definition Runtime.h:66
T isfinite(T... args)
FitResults minimize(lsst::afw::math::Function1< ReturnT > const &function, std::vector< double > const &initialParameterList, std::vector< double > const &stepSizeList, std::vector< double > const &measurementList, std::vector< double > const &varianceList, std::vector< double > const &xPositionList, double errorDef)
Find the minimum of a function(x)
Definition minimize.cc:176
T push_back(T... args)
T size(T... args)
Results from minimizing a function.
Definition minimize.h:44
std::vector< double > parameterList
fit parameters
Definition minimize.h:48
bool isValid
true if the fit converged; false otherwise
Definition minimize.h:46
double chiSq
chi squared; may be nan or infinite, but only if isValid false
Definition minimize.h:47
std::vector< std::pair< double, double > > parameterErrorList
negative,positive (1 sigma?) error for each parameter
Definition minimize.h:50