LSSTApplications  20.0.0
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 {
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:
72  std::shared_ptr<Function1<ReturnT> > _functionPtr;
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 {
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  : _functionPtr(function.clone()),
124  _measurementList(measurementList),
125  _varianceList(varianceList),
126  _xPositionList(xPositionList),
127  _errorDef(errorDef) {}
128 
129 template <typename ReturnT>
130 MinimizerFunctionBase2<ReturnT>::MinimizerFunctionBase2(Function2<ReturnT> const &function,
131  std::vector<double> const &measurementList,
132  std::vector<double> const &varianceList,
133  std::vector<double> const &xPositionList,
134  std::vector<double> const &yPositionList,
135  double errorDef)
136  : _functionPtr(function.clone()),
137  _measurementList(measurementList),
138  _varianceList(varianceList),
139  _xPositionList(xPositionList),
140  _yPositionList(yPositionList),
141  _errorDef(errorDef) {}
142 
143 // Only method we need to set up; basically this is a chi^2 routine
144 template <typename ReturnT>
145 double MinimizerFunctionBase1<ReturnT>::operator()(const std::vector<double> &par) const {
146  // Initialize the function with the fit parameters
147  this->_functionPtr->setParameters(par);
148 
149  double chi2 = 0.0;
150  for (unsigned int i = 0; i < this->_measurementList.size(); i++) {
151  double resid = (*(this->_functionPtr))(this->_xPositionList[i]) - this->_measurementList[i];
152  chi2 += resid * resid / this->_varianceList[i];
153  }
154 
155  return chi2;
156 }
157 
158 template <typename ReturnT>
159 double MinimizerFunctionBase2<ReturnT>::operator()(const std::vector<double> &par) const {
160  // Initialize the function with the fit parameters
161  this->_functionPtr->setParameters(par);
162 
163  double chi2 = 0.0;
164  for (unsigned int i = 0; i < this->_measurementList.size(); i++) {
165  double resid = (*(this->_functionPtr))(this->_xPositionList[i], this->_yPositionList[i]) -
166  this->_measurementList[i];
167  chi2 += resid * resid / this->_varianceList[i];
168  }
169 
170  return chi2;
171 }
173 
174 template <typename ReturnT>
175 FitResults minimize(Function1<ReturnT> const &function, std::vector<double> const &initialParameterList,
176  std::vector<double> const &stepSizeList, std::vector<double> const &measurementList,
177  std::vector<double> const &varianceList, std::vector<double> const &xPositionList,
178  double errorDef) {
179  unsigned int const nParameters = function.getNParameters();
180  if (initialParameterList.size() != nParameters) {
181  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "initialParameterList is the wrong length");
182  }
183  if (stepSizeList.size() != nParameters) {
184  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "stepSizeList is the wrong length");
185  }
186  unsigned int const nMeasurements = measurementList.size();
187  if (varianceList.size() != nMeasurements) {
188  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "varianceList is the wrong length");
189  }
190  if (xPositionList.size() != nMeasurements) {
191  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "xPositionList is the wrong length");
192  }
193 
194  MinimizerFunctionBase1<ReturnT> minimizerFunc(function, measurementList, varianceList, xPositionList,
195  errorDef);
196 
197  ROOT::Minuit2::MnUserParameters fitPar;
198  std::vector<std::string> paramNames;
199  for (unsigned int i = 0; i < nParameters; ++i) {
200  paramNames.push_back((boost::format("p%d") % i).str());
201  fitPar.Add(paramNames[i].c_str(), initialParameterList[i], stepSizeList[i]);
202  }
203 
204  ROOT::Minuit2::MnMigrad migrad(minimizerFunc, fitPar);
205  ROOT::Minuit2::FunctionMinimum min = migrad();
206  ROOT::Minuit2::MnMinos minos(minimizerFunc, min);
207 
208  FitResults fitResults;
209  fitResults.chiSq = min.Fval();
210  fitResults.isValid = min.IsValid() && std::isfinite(fitResults.chiSq);
211  if (!fitResults.isValid) {
212  LOGL_WARN("afw.math.minimize", "Fit failed to converge");
213  }
214 
215  for (unsigned int i = 0; i < nParameters; ++i) {
216  fitResults.parameterList.push_back(min.UserState().Value(paramNames[i].c_str()));
217  if (fitResults.isValid) {
218  fitResults.parameterErrorList.push_back(minos(i));
219  } else {
220  double e = min.UserState().Error(paramNames[i].c_str());
221  std::pair<double, double> ep(-1 * e, e);
222  fitResults.parameterErrorList.push_back(ep);
223  }
224  }
225  return fitResults;
226 }
227 
228 template <typename ReturnT>
229 FitResults minimize(Function2<ReturnT> const &function, std::vector<double> const &initialParameterList,
230  std::vector<double> const &stepSizeList, std::vector<double> const &measurementList,
231  std::vector<double> const &varianceList, std::vector<double> const &xPositionList,
232  std::vector<double> const &yPositionList, double errorDef) {
233  unsigned int const nParameters = function.getNParameters();
234  if (initialParameterList.size() != nParameters) {
235  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "initialParameterList is the wrong length");
236  }
237  if (stepSizeList.size() != nParameters) {
238  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "stepSizeList is the wrong length");
239  }
240  unsigned int const nMeasurements = measurementList.size();
241  if (varianceList.size() != nMeasurements) {
242  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "varianceList is the wrong length");
243  }
244  if (xPositionList.size() != nMeasurements) {
245  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "xPositionList is the wrong length");
246  }
247  if (yPositionList.size() != nMeasurements) {
248  throw LSST_EXCEPT(pex::exceptions::InvalidParameterError, "yPositionList is the wrong length");
249  }
250 
251  MinimizerFunctionBase2<ReturnT> minimizerFunc(function, measurementList, varianceList, xPositionList,
252  yPositionList, errorDef);
253 
254  ROOT::Minuit2::MnUserParameters fitPar;
255  std::vector<std::string> paramNames;
256  for (unsigned int i = 0; i < nParameters; ++i) {
257  paramNames.push_back((boost::format("p%d") % i).str());
258  fitPar.Add(paramNames[i].c_str(), initialParameterList[i], stepSizeList[i]);
259  }
260 
261  ROOT::Minuit2::MnMigrad migrad(minimizerFunc, fitPar);
262  ROOT::Minuit2::FunctionMinimum min = migrad();
263  ROOT::Minuit2::MnMinos minos(minimizerFunc, min);
264 
265  FitResults fitResults;
266  fitResults.chiSq = min.Fval();
267  fitResults.isValid = min.IsValid() && std::isfinite(fitResults.chiSq);
268  if (!fitResults.isValid) {
269  LOGL_WARN("afw.math.minimize", "Fit failed to converge");
270  }
271 
272  for (unsigned int i = 0; i < nParameters; ++i) {
273  fitResults.parameterList.push_back(min.UserState().Value(paramNames[i].c_str()));
274  if (fitResults.isValid) {
275  fitResults.parameterErrorList.push_back(minos(i));
276  } else {
277  double e = min.UserState().Error(paramNames[i].c_str());
278  std::pair<double, double> ep(-1 * e, e);
279  fitResults.parameterErrorList.push_back(ep);
280  }
281  }
282  return fitResults;
283 }
284 
285 // Explicit instantiation
287 #define NL /* */
288 #define minimizeFuncs(ReturnT) \
289  template FitResults minimize(Function1<ReturnT> const &, std::vector<double> const &, \
290  std::vector<double> const &, std::vector<double> const &, \
291  std::vector<double> const &, std::vector<double> const &, double); \
292  NL template FitResults minimize(Function2<ReturnT> const &, std::vector<double> const &, \
293  std::vector<double> const &, std::vector<double> const &, \
294  std::vector<double> const &, std::vector<double> const &, \
295  std::vector<double> const &, double);
296 
297 minimizeFuncs(float) minimizeFuncs(double)
299 } // namespace math
300 } // namespace afw
301 } // namespace lsst
lsst::afw::math::FitResults
Results from minimizing a function.
Definition: minimize.h:44
std::shared_ptr
STL class.
std::pair< double, double >
lsst::afw::math::FitResults::isValid
bool isValid
true if the fit converged; false otherwise
Definition: minimize.h:46
std::vector< double >
std::vector::size
T size(T... args)
LOGL_WARN
#define LOGL_WARN(logger, message...)
Definition: Log.h:536
pex.config.history.format
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174
lsst::afw
Definition: imageAlgorithm.dox:1
minimize.h
lsst::afw::math::minimize
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:175
std::vector::push_back
T push_back(T... args)
lsst::afw::math::FitResults::parameterErrorList
std::vector< std::pair< double, double > > parameterErrorList
negative,positive (1 sigma?) error for each parameter
Definition: minimize.h:50
std::isfinite
T isfinite(T... args)
lsst::afw::math::FitResults::chiSq
double chiSq
chi squared; may be nan or infinite, but only if isValid false
Definition: minimize.h:47
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::afw::math::Function2
A Function taking two arguments.
Definition: Function.h:259
lsst::pex::exceptions::InvalidParameterError
Reports invalid arguments.
Definition: Runtime.h:66
min
int min
Definition: BoundedField.cc:103
lsst::afw::math::FitResults::parameterList
std::vector< double > parameterList
fit parameters
Definition: minimize.h:48
lsst::afw::math::Function1
A Function taking one argument.
Definition: Function.h:202
lsst::afw::image.slicing.clone
clone
Definition: slicing.py:257
Log.h
LSST DM logging module built on log4cxx.
Generated on Wed Jun 24 2020 18:10:02 for LSSTApplications by   doxygen 1.8.18