Logo ROOT   6.12/07
Reference Guide
MnMinos.cxx
Go to the documentation of this file.
1 // @(#)root/minuit2:$Id$
2 // Authors: M. Winkler, F. James, L. Moneta, A. Zsenei 2003-2005
3 
4 /**********************************************************************
5  * *
6  * Copyright (c) 2005 LCG ROOT Math team, CERN/PH-SFT *
7  * *
8  **********************************************************************/
9 
10 #include "Minuit2/MnMinos.h"
12 #include "Minuit2/FCNBase.h"
14 #include "Minuit2/MnCross.h"
15 #include "Minuit2/MinosError.h"
16 
17 //#define DEBUG
18 
19 #if defined(DEBUG) || defined(WARNINGMSG)
20 #include "Minuit2/MnPrint.h"
21 #endif
22 
23 
24 namespace ROOT {
25 
26  namespace Minuit2 {
27 
28 
29 MnMinos::MnMinos(const FCNBase& fcn, const FunctionMinimum& min, unsigned int stra ) :
30  fFCN(fcn),
31  fMinimum(min),
32  fStrategy(MnStrategy(stra))
33 {
34  // construct from FCN + Minimum
35  // check if Error definition has been changed, in case re-update errors
36  if (fcn.Up() != min.Up() ) {
37 #ifdef WARNINGMSG
38  MN_INFO_MSG("MnMinos UP value has changed, need to update FunctionMinimum class");
39 #endif
40  }
41 }
42 
43 MnMinos::MnMinos(const FCNBase& fcn, const FunctionMinimum& min, const MnStrategy& stra) :
44  fFCN(fcn),
45  fMinimum(min),
46  fStrategy(stra)
47 {
48  // construct from FCN + Minimum
49  // check if Error definition has been changed, in case re-update errors
50  if (fcn.Up() != min.Up() ) {
51 #ifdef WARNINGMSG
52  MN_INFO_MSG("MnMinos UP value has changed, need to update FunctionMinimum class");
53 #endif
54  }
55 }
56 
57 
58 std::pair<double,double> MnMinos::operator()(unsigned int par, unsigned int maxcalls, double toler) const {
59  // do Minos analysis given the parameter index returning a pair for (lower,upper) errors
60  MinosError mnerr = Minos(par, maxcalls,toler);
61  return mnerr();
62 }
63 
64 double MnMinos::Lower(unsigned int par, unsigned int maxcalls, double toler) const {
65  // get lower error for parameter par
67  double err = fMinimum.UserState().Error(par);
68 
69  MnCross aopt = Loval(par, maxcalls,toler);
70 
71  double lower = aopt.IsValid() ? -1.*err*(1.+ aopt.Value()) : (aopt.AtLimit() ? upar.Parameter(par).LowerLimit() : upar.Value(par));
72 
73  return lower;
74 }
75 
76 double MnMinos::Upper(unsigned int par, unsigned int maxcalls, double toler) const {
77  // upper error for parameter par
78  MnCross aopt = Upval(par, maxcalls,toler);
79 
81  double err = fMinimum.UserState().Error(par);
82 
83  double upper = aopt.IsValid() ? err*(1.+ aopt.Value()) : (aopt.AtLimit() ? upar.Parameter(par).UpperLimit() : upar.Value(par));
84 
85  return upper;
86 }
87 
88 MinosError MnMinos::Minos(unsigned int par, unsigned int maxcalls, double toler) const {
89  // do full minos error anlysis (lower + upper) for parameter par
90  assert(fMinimum.IsValid());
91  assert(!fMinimum.UserState().Parameter(par).IsFixed());
92  assert(!fMinimum.UserState().Parameter(par).IsConst());
93 
94  MnCross up = Upval(par, maxcalls,toler);
95 #ifdef DEBUG
96  std::cout << "Function calls to find upper error " << up.NFcn() << std::endl;
97 #endif
98 
99  MnCross lo = Loval(par, maxcalls,toler);
100 
101 #ifdef DEBUG
102  std::cout << "Function calls to find lower error " << lo.NFcn() << std::endl;
103 #endif
104 
105  return MinosError(par, fMinimum.UserState().Value(par), lo, up);
106 }
107 
108 
109 MnCross MnMinos::FindCrossValue(int direction, unsigned int par, unsigned int maxcalls, double toler) const {
110  // get crossing value in the parameter direction :
111  // direction = + 1 upper value
112  // direction = -1 lower value
113  // pass now tolerance used for Migrad minimizations
114 
115  assert(direction == 1 || direction == -1);
116 #ifdef DEBUG
117  if (direction == 1)
118  std::cout << "\n--------- MnMinos --------- \n Determination of positive Minos error for parameter "
119  << par << std::endl;
120  else
121  std::cout << "\n--------- MnMinos --------- \n Determination of positive Minos error for parameter "
122  << par << std::endl;
123 #endif
124 
125  assert(fMinimum.IsValid());
126  assert(!fMinimum.UserState().Parameter(par).IsFixed());
127  assert(!fMinimum.UserState().Parameter(par).IsConst());
128  if(maxcalls == 0) {
129  unsigned int nvar = fMinimum.UserState().VariableParameters();
130  maxcalls = 2*(nvar+1)*(200 + 100*nvar + 5*nvar*nvar);
131  }
132 
133  std::vector<unsigned int> para(1, par);
134 
136  double err = direction * upar.Error(par);
137  double val = upar.Value(par) + err;
138  std::vector<double> xmid(1, val);
139  std::vector<double> xdir(1, err);
140 
141  double up = fFCN.Up();
142  unsigned int ind = upar.IntOfExt(par);
143  // get error matrix (methods return a copy)
145  // get internal parameters
146  const MnAlgebraicVector & xt = fMinimum.Parameters().Vec();
147  //LM: change to use err**2 (m(i,i) instead of err as in F77 version
148  double xunit = sqrt(up/m(ind,ind));
149  // LM (29/04/08) bug: change should be done in internal variables
150  for(unsigned int i = 0; i < m.Nrow(); i++) {
151  if(i == ind) continue;
152  double xdev = xunit*m(ind,i);
153  double xnew = xt(i) + direction * xdev;
154 
155  // transform to external values
156  unsigned int ext = upar.ExtOfInt(i);
157 
158  double unew = upar.Int2ext(i, xnew);
159 
160 #ifdef DEBUG
161  std::cout << "Parameter " << ext << " is set from " << upar.Value(ext) << " to " << unew << std::endl;
162 #endif
163  upar.SetValue(ext, unew);
164  }
165 
166  upar.Fix(par);
167  upar.SetValue(par, val);
168 
169 #ifdef DEBUG
170  std::cout << "Parameter " << par << " is fixed and set from " << fMinimum.UserState().Value(par) << " to " << val << std::endl;
171 #endif
172 
173 
174  MnFunctionCross cross(fFCN, upar, fMinimum.Fval(), fStrategy);
175  MnCross aopt = cross(para, xmid, xdir, toler, maxcalls);
176 
177 
178 #ifdef DEBUG
179  std::cout<<"----- MnMinos: aopt found from MnFunctionCross = "<<aopt.Value()<<std::endl << std::endl;
180 #endif
181 
182 #ifdef WARNINGMSG
183  const char * par_name = upar.Name(par);
184  if(aopt.AtMaxFcn())
185  MN_INFO_VAL2("MnMinos maximum number of function calls exceeded for Parameter ",par_name);
186  if(aopt.NewMinimum())
187  MN_INFO_VAL2("MnMinos new Minimum found while looking for Parameter ",par_name);
188  if (direction ==1) {
189  if(aopt.AtLimit())
190  MN_INFO_VAL2("MnMinos Parameter is at Upper limit.",par_name);
191  if(!aopt.IsValid())
192  MN_INFO_VAL2("MnMinos could not find Upper Value for Parameter ",par_name);
193  }
194  else {
195  if(aopt.AtLimit())
196  MN_INFO_VAL2("MnMinos Parameter is at Lower limit.",par_name);
197  if(!aopt.IsValid())
198  MN_INFO_VAL2("MnMinos could not find Lower Value for Parameter ",par_name);
199  }
200 #endif
201 
202  return aopt;
203 }
204 
205 MnCross MnMinos::Upval(unsigned int par, unsigned int maxcalls, double toler) const {
206  // return crossing in the lower parameter direction
207  return FindCrossValue(1,par,maxcalls,toler);
208 }
209 
210 MnCross MnMinos::Loval(unsigned int par, unsigned int maxcalls, double toler) const {
211  // return crossing in the lower parameter direction
212  return FindCrossValue(-1,par,maxcalls,toler);
213 }
214 
215 // #ifdef DEBUG
216 // std::cout << "\n--------- MnMinos --------- \n Determination of negative Minos error for parameter "
217 // << par << std::endl;
218 // #endif
219 
220 // assert(fMinimum.IsValid());
221 // assert(!fMinimum.UserState().Parameter(par).IsFixed());
222 // assert(!fMinimum.UserState().Parameter(par).IsConst());
223 // if(maxcalls == 0) {
224 // unsigned int nvar = fMinimum.UserState().VariableParameters();
225 // maxcalls = 2*(nvar+1)*(200 + 100*nvar + 5*nvar*nvar);
226 // }
227 // std::vector<unsigned int> para(1, par);
228 
229 // MnUserParameterState upar = fMinimum.UserState();
230 // double err = upar.Error(par);
231 // double val = upar.Value(par) - err;
232 // std::vector<double> xmid(1, val);
233 // std::vector<double> xdir(1, -err);
234 
235 // double up = fFCN.Up();
236 // unsigned int ind = upar.IntOfExt(par);
237 // MnAlgebraicSymMatrix m = fMinimum.Error().Matrix();
238 // double xunit = sqrt(up/m(ind,ind));
239 // // get internal parameters
240 // const MnAlgebraicVector & xt = fMinimum.Parameters().Vec();
241 
242 // for(unsigned int i = 0; i < m.Nrow(); i++) {
243 // if(i == ind) continue;
244 // double xdev = xunit*m(ind,i);
245 
246 // double xnew = xt(i) - xdev;
247 
248 // // transform to external values
249 // double unew = upar.Int2ext(i, xnew);
250 
251 // unsigned int ext = upar.ExtOfInt(i);
252 
253 // #ifdef DEBUG
254 // std::cout << "Parameter " << ext << " is set from " << upar.Value(ext) << " to " << unew << std::endl;
255 // #endif
256 // upar.SetValue(ext, unew);
257 // }
258 
259 // upar.Fix(par);
260 // upar.SetValue(par, val);
261 
262 // #ifdef DEBUG
263 // std::cout << "Parameter " << par << " is fixed and set from " << fMinimum.UserState().Value(par) << " to " << val << std::endl;
264 // #endif
265 
266 // // double edmmax = 0.5*0.1*fFCN.Up()*1.e-3;
267 // double toler = 0.01;
268 // MnFunctionCross cross(fFCN, upar, fMinimum.Fval(), fStrategy);
269 
270 // MnCross aopt = cross(para, xmid, xdir, toler, maxcalls);
271 
272 // #ifdef DEBUG
273 // std::cout<<"----- MnMinos: aopt found from MnFunctionCross = "<<aopt.Value()<<std::endl << std::endl;
274 // #endif
275 
276 // #ifdef WARNINGMSG
277 // if(aopt.AtLimit())
278 // MN_INFO_VAL2("MnMinos Parameter is at Lower limit.",par);
279 // if(aopt.AtMaxFcn())
280 // MN_INFO_VAL2("MnMinos maximum number of function calls exceeded for Parameter ",par);
281 // if(aopt.NewMinimum())
282 // MN_INFO_VAL2("MnMinos new Minimum found while looking for Parameter ",par);
283 // if(!aopt.IsValid())
284 // MN_INFO_VAL2("MnMinos could not find Lower Value for Parameter ",par);
285 // #endif
286 
287 // return aopt;
288 
289 // }
290 
291 
292  } // namespace Minuit2
293 
294 } // namespace ROOT
virtual double Up() const =0
Error definition of the function.
MnMinos(const FCNBase &fcn, const FunctionMinimum &min, unsigned int stra=1)
construct from FCN + Minimum + strategy
Definition: MnMinos.cxx:29
#define MN_INFO_VAL2(loc, x)
Definition: MnPrint.h:130
unsigned int NFcn() const
Definition: MnCross.h:65
MnCross Upval(unsigned int, unsigned int maxcalls=0, double toler=0.1) const
Definition: MnMinos.cxx:205
auto * m
Definition: textangle.C:8
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
const MnAlgebraicVector & Vec() const
MinosError Minos(unsigned int, unsigned int maxcalls=0, double toler=0.1) const
ask for MinosError (Lower + Upper) can be printed via std::cout
Definition: MnMinos.cxx:88
bool AtLimit() const
Definition: MnCross.h:62
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
double Lower(unsigned int, unsigned int maxcalls=0, double toler=0.1) const
calculate one side (negative or positive Error) of the Parameter give as input (optionally) maxcalls ...
Definition: MnMinos.cxx:64
const MinuitParameter & Parameter(unsigned int i) const
MnCross Loval(unsigned int, unsigned int maxcalls=0, double toler=0.1) const
Definition: MnMinos.cxx:210
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
unsigned int ExtOfInt(unsigned int) const
double sqrt(double)
const FunctionMinimum & fMinimum
Definition: MnMinos.h:71
MnCross FindCrossValue(int dir, unsigned int, unsigned int maxcalls, double toler) const
internal method to get crossing value via MnFunctionCross
Definition: MnMinos.cxx:109
class holding the full result of the minimization; both internal and external (MnUserParameterState) ...
const MinimumParameters & Parameters() const
bool IsValid() const
Definition: MnCross.h:61
Interface (abstract class) defining the function to be minimized, which has to be implemented by the ...
Definition: FCNBase.h:47
double Upper(unsigned int, unsigned int maxcalls=0, double toler=0.1) const
Definition: MnMinos.cxx:76
MnStrategy fStrategy
Definition: MnMinos.h:72
unsigned int IntOfExt(unsigned int) const
unsigned int Nrow() const
Definition: LASymMatrix.h:239
Class holding the result of Minos (lower and upper values) for a specific parameter.
Definition: MinosError.h:25
MnAlgebraicSymMatrix Matrix() const
Definition: MinimumError.h:58
class which holds the external user and/or internal Minuit representation of the parameters and error...
double Value() const
Definition: MnCross.h:59
const char * Name(unsigned int) const
std::pair< double, double > operator()(unsigned int, unsigned int maxcalls=0, double toler=0.1) const
returns the negative (pair.first) and the positive (pair.second) Minos Error of the Parameter ...
Definition: MnMinos.cxx:58
const FCNBase & fFCN
Definition: MnMinos.h:70
double Int2ext(unsigned int, double) const
const MnUserParameterState & UserState() const
API class for defining three levels of strategies: low (0), medium (1), high (>=2); acts on: Migrad (...
Definition: MnStrategy.h:27
const MinimumError & Error() const