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