ROOT  6.06/09
Reference Guide
MnUserParameterState.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 
12 #include "Minuit2/MinimumState.h"
13 
14 #include "Minuit2/MnPrint.h"
15 
16 
17 namespace ROOT {
18 
19  namespace Minuit2 {
20 
21 
22 //
23 // construct from user parameters (befor minimization)
24 //
25 MnUserParameterState::MnUserParameterState(const std::vector<double>& par, const std::vector<double>& err) :
26  fValid(true), fCovarianceValid(false), fGCCValid(false), fCovStatus(-1), fFVal(0.), fEDM(0.), fNFcn(0), fParameters(MnUserParameters(par, err)), fCovariance(MnUserCovariance()), fGlobalCC(MnGlobalCorrelationCoeff()), fIntParameters(par), fIntCovariance(MnUserCovariance())
27  {}
28 
30  fValid(true), fCovarianceValid(false), fGCCValid(false), fCovStatus(-1), fFVal(0.), fEDM(0.), fNFcn(0), fParameters(par), fCovariance(MnUserCovariance()), fGlobalCC(MnGlobalCorrelationCoeff()), fIntParameters(std::vector<double>()), fIntCovariance(MnUserCovariance()) {
31  // construct from user parameters (befor minimization)
32 
33  for(std::vector<MinuitParameter>::const_iterator ipar = MinuitParameters().begin(); ipar != MinuitParameters().end(); ipar++) {
34  if((*ipar).IsConst() || (*ipar).IsFixed()) continue;
35  if((*ipar).HasLimits())
36  fIntParameters.push_back(Ext2int((*ipar).Number(), (*ipar).Value()));
37  else
38  fIntParameters.push_back((*ipar).Value());
39  }
40 }
41 
42 //
43 // construct from user parameters + errors (befor minimization)
44 //
45 MnUserParameterState::MnUserParameterState(const std::vector<double>& par, const std::vector<double>& cov, unsigned int nrow) :
46  fValid(true), fCovarianceValid(true), fGCCValid(false), fCovStatus(-1), fFVal(0.), fEDM(0.), fNFcn(0),
47  fParameters(MnUserParameters()), fCovariance(MnUserCovariance(cov, nrow)), fGlobalCC(MnGlobalCorrelationCoeff()), fIntParameters(par), fIntCovariance(MnUserCovariance(cov, nrow)) {
48  // construct from user parameters + errors (before minimization) using std::vector for parameter error and // an std::vector of size n*(n+1)/2 for the covariance matrix and n (rank of cov matrix)
49 
50  std::vector<double> err; err.reserve(par.size());
51  for(unsigned int i = 0; i < par.size(); i++) {
52  assert(fCovariance(i,i) > 0.);
53  err.push_back(sqrt(fCovariance(i,i)));
54  }
55  fParameters = MnUserParameters(par, err);
57 }
58 
59 MnUserParameterState::MnUserParameterState(const std::vector<double>& par, const MnUserCovariance& cov) :
60  fValid(true), fCovarianceValid(true), fGCCValid(false), fCovStatus(-1), fFVal(0.), fEDM(0.), fNFcn(0),
61  fParameters(MnUserParameters()), fCovariance(cov), fGlobalCC(MnGlobalCorrelationCoeff()), fIntParameters(par), fIntCovariance(cov) {
62  //construct from user parameters + errors (before minimization) using std::vector (params) and MnUserCovariance class
63 
64  std::vector<double> err; err.reserve(par.size());
65  for(unsigned int i = 0; i < par.size(); i++) {
66  assert(fCovariance(i,i) > 0.);
67  err.push_back(sqrt(fCovariance(i,i)));
68  }
69  fParameters = MnUserParameters(par, err);
71 }
72 
73 
75  fValid(true), fCovarianceValid(true), fGCCValid(false), fCovStatus(-1), fFVal(0.), fEDM(0.), fNFcn(0),
76  fParameters(par), fCovariance(cov), fGlobalCC(MnGlobalCorrelationCoeff()), fIntParameters(std::vector<double>()), fIntCovariance(cov) {
77  //construct from user parameters + errors (befor minimization) using
78  // MnUserParameters and MnUserCovariance objects
79 
80  fIntCovariance.Scale(0.5);
81  for(std::vector<MinuitParameter>::const_iterator ipar = MinuitParameters().begin(); ipar != MinuitParameters().end(); ipar++) {
82  if((*ipar).IsConst() || (*ipar).IsFixed()) continue;
83  if((*ipar).HasLimits())
84  fIntParameters.push_back(Ext2int((*ipar).Number(), (*ipar).Value()));
85  else
86  fIntParameters.push_back((*ipar).Value());
87  }
89  //
90  // need to Fix that in case of limited parameters
91  // fIntCovariance = MnUserCovariance();
92  //
93 }
94 
95 //
96 //
98  fValid(st.IsValid()), fCovarianceValid(false), fGCCValid(false), fCovStatus(-1),
99  fFVal(st.Fval()), fEDM(st.Edm()), fNFcn(st.NFcn()), fParameters(MnUserParameters()), fCovariance(MnUserCovariance()), fGlobalCC(MnGlobalCorrelationCoeff()), fIntParameters(std::vector<double>()), fIntCovariance(MnUserCovariance()) {
100  //
101  // construct from internal parameters (after minimization)
102  //
103  //std::cout << "build a MnUSerParameterState after minimization.." << std::endl;
104 
105  for(std::vector<MinuitParameter>::const_iterator ipar = trafo.Parameters().begin(); ipar != trafo.Parameters().end(); ipar++) {
106  if((*ipar).IsConst()) {
107  Add((*ipar).GetName(), (*ipar).Value());
108  } else if((*ipar).IsFixed()) {
109  Add((*ipar).GetName(), (*ipar).Value(), (*ipar).Error());
110  if((*ipar).HasLimits()) {
111  if((*ipar).HasLowerLimit() && (*ipar).HasUpperLimit())
112  SetLimits((*ipar).GetName(), (*ipar).LowerLimit(),(*ipar).UpperLimit());
113  else if((*ipar).HasLowerLimit() && !(*ipar).HasUpperLimit())
114  SetLowerLimit((*ipar).GetName(), (*ipar).LowerLimit());
115  else
116  SetUpperLimit((*ipar).GetName(), (*ipar).UpperLimit());
117  }
118  Fix((*ipar).GetName());
119  } else if((*ipar).HasLimits()) {
120  unsigned int i = trafo.IntOfExt((*ipar).Number());
121  double err = st.HasCovariance() ? sqrt(2.*up*st.Error().InvHessian()(i,i)) : st.Parameters().Dirin()(i);
122  Add((*ipar).GetName(), trafo.Int2ext(i, st.Vec()(i)), trafo.Int2extError(i, st.Vec()(i), err));
123  if((*ipar).HasLowerLimit() && (*ipar).HasUpperLimit())
124  SetLimits((*ipar).GetName(), (*ipar).LowerLimit(), (*ipar).UpperLimit());
125  else if((*ipar).HasLowerLimit() && !(*ipar).HasUpperLimit())
126  SetLowerLimit((*ipar).GetName(), (*ipar).LowerLimit());
127  else
128  SetUpperLimit((*ipar).GetName(), (*ipar).UpperLimit());
129  } else {
130  unsigned int i = trafo.IntOfExt((*ipar).Number());
131  double err = st.HasCovariance() ? sqrt(2.*up*st.Error().InvHessian()(i,i)) : st.Parameters().Dirin()(i);
132  Add((*ipar).GetName(), st.Vec()(i), err);
133  }
134  }
135 
136  // need to be set afterwards because becore the ::Add method set fCovarianceValid to false
137  fCovarianceValid = st.Error().IsValid();
138 
139  fCovStatus = -1; // when not available
140  //if (st.Error().HesseFailed() || st.Error().InvertFailed() ) fCovStatus = -1;
141  // when available
142  if (st.Error().IsAvailable() ) fCovStatus = 0;
143 
144  if(fCovarianceValid) {
145  fCovariance = trafo.Int2extCovariance(st.Vec(), st.Error().InvHessian());
146  fIntCovariance = MnUserCovariance(std::vector<double>(st.Error().InvHessian().Data(), st.Error().InvHessian().Data()+st.Error().InvHessian().size()), st.Error().InvHessian().Nrow());
147  fCovariance.Scale(2.*up);
150 
152 
153  fCovStatus = 1; // when is valid
154  }
155  if (st.Error().IsMadePosDef() ) fCovStatus = 2;
156  if (st.Error().IsAccurate() ) fCovStatus = 3;
157 
158 }
159 
161  // invert covariance matrix and return Hessian
162  // need to copy in a MnSymMatrix
164  std::copy(fCovariance.Data().begin(), fCovariance.Data().end(), mat.Data() );
165  int ifail = Invert(mat);
166  if(ifail != 0) {
167 #ifdef WARNINGMSG
168  MN_INFO_MSG("MnUserParameterState:Hessian inversion fails- return diagonal matrix.");
169 #endif
171  for(unsigned int i = 0; i < fCovariance.Nrow(); i++) {
172  tmp(i,i) = 1./fCovariance(i,i);
173  }
174  return tmp;
175  }
176 
177  MnUserCovariance hessian( mat.Data(), fCovariance.Nrow() );
178  return hessian;
179 }
180 
181 // facade: forward interface of MnUserParameters and MnUserTransformation
182 // via MnUserParameterState
183 
184 
185 const std::vector<MinuitParameter>& MnUserParameterState::MinuitParameters() const {
186  //access to parameters (row-wise)
187  return fParameters.Parameters();
188 }
189 
190 std::vector<double> MnUserParameterState::Params() const {
191  //access to parameters in column-wise representation
192  return fParameters.Params();
193 }
194 std::vector<double> MnUserParameterState::Errors() const {
195  //access to errors in column-wise representation
196  return fParameters.Errors();
197 }
198 
199 const MinuitParameter& MnUserParameterState::Parameter(unsigned int i) const {
200  //access to single Parameter i
201  return fParameters.Parameter(i);
202 }
203 
204 void MnUserParameterState::Add(const std::string & name, double val, double err) {
205  //add free Parameter
206  if ( fParameters.Add(name, val, err) ) {
207  fIntParameters.push_back(val);
208  fCovarianceValid = false;
209  fGCCValid = false;
210  fValid = true;
211  }
212  else {
213  // redefine an existing parameter
214  int i = Index(name);
215  SetValue(i,val);
216  if (Parameter(i).IsConst() ) {
217  std::string msg = "Cannot modify status of constant parameter " + name;
218  MN_INFO_MSG2("MnUserParameterState::Add",msg.c_str());
219  return;
220  }
221  SetError(i,err);
222  // release if it was fixed
223  if (Parameter(i).IsFixed() ) Release(i);
224  }
225 
226 }
227 
228 void MnUserParameterState::Add(const std::string & name, double val, double err, double low, double up) {
229  //add limited Parameter
230  if ( fParameters.Add(name, val, err, low, up) ) {
231  fCovarianceValid = false;
232  fIntParameters.push_back(Ext2int(Index(name), val));
233  fGCCValid = false;
234  fValid = true;
235  }
236  else { // Parameter already exist - just set values
237  int i = Index(name);
238  SetValue(i,val);
239  if (Parameter(i).IsConst() ) {
240  std::string msg = "Cannot modify status of constant parameter " + name;
241  MN_INFO_MSG2("MnUserParameterState::Add",msg.c_str());
242  return;
243  }
244  SetError(i,err);
245  SetLimits(i,low,up);
246  // release if it was fixed
247  if (Parameter(i).IsFixed() ) Release(i);
248  }
249 
250 
251 }
252 
253 void MnUserParameterState::Add(const std::string & name, double val) {
254  //add const Parameter
255  if ( fParameters.Add(name, val) )
256  fValid = true;
257  else
258  SetValue(name,val);
259 }
260 
261 //interaction via external number of Parameter
262 
263 void MnUserParameterState::Fix(unsigned int e) {
264  // fix parameter e (external index)
265  if(!Parameter(e).IsFixed() && !Parameter(e).IsConst()) {
266  unsigned int i = IntOfExt(e);
267  if(fCovarianceValid) {
270  }
271  fIntParameters.erase(fIntParameters.begin()+i, fIntParameters.begin()+i+1);
272  }
273  fParameters.Fix(e);
274  fGCCValid = false;
275 }
276 
277 void MnUserParameterState::Release(unsigned int e) {
278  // release parameter e (external index)
279  // no-op if parameter is const
280  if (Parameter(e).IsConst() ) return;
281  fParameters.Release(e);
282  fCovarianceValid = false;
283  fGCCValid = false;
284  unsigned int i = IntOfExt(e);
285  if(Parameter(e).HasLimits())
286  fIntParameters.insert(fIntParameters.begin()+i, Ext2int(e, Parameter(e).Value()));
287  else
288  fIntParameters.insert(fIntParameters.begin()+i, Parameter(e).Value());
289 }
290 
291 void MnUserParameterState::SetValue(unsigned int e, double val) {
292  // set value for parameter e ( external index )
293  fParameters.SetValue(e, val);
294  if(!Parameter(e).IsFixed() && !Parameter(e).IsConst()) {
295  unsigned int i = IntOfExt(e);
296  if(Parameter(e).HasLimits())
297  fIntParameters[i] = Ext2int(e, val);
298  else
299  fIntParameters[i] = val;
300  }
301 }
302 
303 void MnUserParameterState::SetError(unsigned int e, double val) {
304  // set error for parameter e (external index)
305  fParameters.SetError(e, val);
306 }
307 
308 void MnUserParameterState::SetLimits(unsigned int e, double low, double up) {
309  // set limits for parameter e (external index)
310  fParameters.SetLimits(e, low, up);
311  fCovarianceValid = false;
312  fGCCValid = false;
313  if(!Parameter(e).IsFixed() && !Parameter(e).IsConst()) {
314  unsigned int i = IntOfExt(e);
315  if(low < fIntParameters[i] && fIntParameters[i] < up)
317  else if (low >= fIntParameters[i] )
318  fIntParameters[i] = Ext2int(e, low + 0.1 * Parameter(e).Error() );
319  else
320  fIntParameters[i] = Ext2int(e, up - 0.1 * Parameter(e).Error() );
321  }
322 }
323 
324 void MnUserParameterState::SetUpperLimit(unsigned int e, double up) {
325  // set upper limit for parameter e (external index)
326  fParameters.SetUpperLimit(e, up);
327  fCovarianceValid = false;
328  fGCCValid = false;
329  if(!Parameter(e).IsFixed() && !Parameter(e).IsConst()) {
330  unsigned int i = IntOfExt(e);
331  if(fIntParameters[i] < up)
333  else
334  fIntParameters[i] = Ext2int(e, up - 0.1 * Parameter(e).Error() );
335  }
336 }
337 
338 void MnUserParameterState::SetLowerLimit(unsigned int e, double low) {
339  // set lower limit for parameter e (external index)
340  fParameters.SetLowerLimit(e, low);
341  fCovarianceValid = false;
342  fGCCValid = false;
343  if(!Parameter(e).IsFixed() && !Parameter(e).IsConst()) {
344  unsigned int i = IntOfExt(e);
345  if(low < fIntParameters[i])
347  else
348  fIntParameters[i] = Ext2int(e, low + 0.1 * Parameter(e).Error() );
349  }
350 }
351 
353  // remove limit for parameter e (external index)
355  fCovarianceValid = false;
356  fGCCValid = false;
357  if(!Parameter(e).IsFixed() && !Parameter(e).IsConst())
358  fIntParameters[IntOfExt(e)] = Value(e);
359 }
360 
361 double MnUserParameterState::Value(unsigned int i) const {
362  // get value for parameter e (external index)
363  return fParameters.Value(i);
364 }
365 double MnUserParameterState::Error(unsigned int i) const {
366  // get error for parameter e (external index)
367  return fParameters.Error(i);
368 }
369 
370 //interaction via name of Parameter
371 
372 void MnUserParameterState::Fix(const std::string & name) { Fix(Index(name));}
373 
374 void MnUserParameterState::Release(const std::string & name) {Release(Index(name));}
375 
376 void MnUserParameterState::SetValue(const std::string & name, double val) {SetValue(Index(name), val);}
377 
378 void MnUserParameterState::SetError(const std::string & name, double val) { SetError(Index(name), val);}
379 
380 void MnUserParameterState::SetLimits(const std::string & name, double low, double up) {SetLimits(Index(name), low, up);}
381 
382 void MnUserParameterState::SetUpperLimit(const std::string & name, double up) { SetUpperLimit(Index(name), up);}
383 
384 void MnUserParameterState::SetLowerLimit(const std::string & name, double low) {SetLowerLimit(Index(name), low);}
385 
386 void MnUserParameterState::RemoveLimits(const std::string & name) {RemoveLimits(Index(name));}
387 
388 double MnUserParameterState::Value(const std::string & name) const {return Value(Index(name));}
389 
390 double MnUserParameterState::Error(const std::string & name) const {return Error(Index(name));}
391 
392 
393 unsigned int MnUserParameterState::Index(const std::string & name) const {
394  //convert name into external number of Parameter
395  return fParameters.Index(name);
396 }
397 
398 const char* MnUserParameterState::Name(unsigned int i) const {
399  //convert external number into name of Parameter (API returing a const char *)
400  return fParameters.Name(i);
401 }
402 const std::string & MnUserParameterState::GetName(unsigned int i) const {
403  //convert external number into name of Parameter (new interface returning a string)
404  return fParameters.GetName(i);
405 }
406 
407 // transformation internal <-> external (forward to transformation class)
408 
409 double MnUserParameterState::Int2ext(unsigned int i, double val) const {
410  // internal to external value
411  return fParameters.Trafo().Int2ext(i, val);
412 }
413 double MnUserParameterState::Ext2int(unsigned int e, double val) const {
414  // external to internal value
415  return fParameters.Trafo().Ext2int(e, val);
416 }
417 unsigned int MnUserParameterState::IntOfExt(unsigned int ext) const {
418  // return internal index for external index ext
419  return fParameters.Trafo().IntOfExt(ext);
420 }
421 unsigned int MnUserParameterState::ExtOfInt(unsigned int internal) const {
422  // return external index for internal index internal
423  return fParameters.Trafo().ExtOfInt(internal);
424 }
426  // return number of variable parameters
428 }
430  // return global parameter precision
431  return fParameters.Precision();
432 }
433 
435  // set global parameter precision
437 }
438 
439  } // namespace Minuit2
440 
441 } // namespace ROOT
double Ext2int(unsigned int, double) const
const char * Name(unsigned int) const
double par[1]
Definition: unuranDistr.cxx:38
unsigned int Index(const std::string &) const
void SetLowerLimit(unsigned int, double)
int Invert(LASymMatrix &)
Definition: LaInverse.cxx:22
double Int2extError(unsigned int, double, double) const
const std::vector< ROOT::Minuit2::MinuitParameter > & MinuitParameters() const
facade: forward interface of MnUserParameters and MnUserTransformation
Namespace for new ROOT classes and functions.
Definition: ROOT.py:1
const std::string & GetName(unsigned int) const
#define assert(cond)
Definition: unittest.h:542
const MnMachinePrecision & Precision() const
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
unsigned int size() const
Definition: LASymMatrix.h:237
const std::vector< ROOT::Minuit2::MinuitParameter > & Parameters() const
access to parameters (row-wise)
double Value(unsigned int) const
std::vector< double > Params() const
unsigned int Nrow() const
Definition: LASymMatrix.h:239
STL namespace.
void Add(const std::string &name, double val, double err)
void Fix(unsigned int)
interaction via external number of Parameter
class for the individual Minuit Parameter with Name and number; contains the input numbers for the mi...
ClassImp(TIterator) Bool_t TIterator return false
Compare two iterator objects.
Definition: TIterator.cxx:20
determines the relative floating point arithmetic precision.
void SetError(unsigned int, double)
bool Add(const std::string &, double, double)
Add free Parameter Name, Value, Error.
const std::vector< MinuitParameter > & Parameters() const
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
const std::string & GetName(unsigned int) const
const char * Name(unsigned int) const
unsigned int IntOfExt(unsigned int) const
const MinimumError & Error() const
Definition: MinimumState.h:62
double sqrt(double)
std::vector< double > Errors() const
double Error(unsigned int) const
void SetLimits(unsigned int, double, double)
double Int2ext(unsigned int, double) const
const MinuitParameter & Parameter(unsigned int) const
access to single Parameter
void SetLimits(unsigned int, double, double)
double Int2ext(unsigned int, double) const
void SetValue(unsigned int, double)
#define MN_INFO_MSG2(loc, str)
Definition: MnPrint.h:124
unsigned int ExtOfInt(unsigned int) const
const MnAlgebraicVector & Vec() const
Definition: MinimumState.h:59
const MnAlgebraicVector & Dirin() const
class to reduce the covariance matrix when a parameter is fixed by removing the corresponding row and...
void SetUpperLimit(unsigned int, double)
const MnMachinePrecision & Precision() const
unsigned int IntOfExt(unsigned int) const
unsigned int ExtOfInt(unsigned int internal) const
void SetLowerLimit(unsigned int, double)
class dealing with the transformation between user specified parameters (external) and internal param...
API class for the user interaction with the parameters; serves as input to the minimizer as well as o...
const MnAlgebraicSymMatrix & InvHessian() const
Definition: MinimumError.h:60
MnUserCovariance Int2extCovariance(const MnAlgebraicVector &, const MnAlgebraicSymMatrix &) const
const MinuitParameter & Parameter(unsigned int i) const
MnUserParameterState()
default constructor (invalid state)
const MinimumParameters & Parameters() const
Definition: MinimumState.h:58
class for global correlation coefficient
std::vector< double > Params() const
access to parameters and errors in column-wise representation
#define name(a, b)
Definition: linkTestLib0.cpp:5
double Ext2int(unsigned int, double) const
void SetUpperLimit(unsigned int, double)
unsigned int Index(const std::string &) const
const double * Data() const
Definition: LASymMatrix.h:233
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
std::vector< double > Errors() const
const std::vector< double > & Data() const
const MnUserTransformation & Trafo() const
Class containing the covariance matrix data represented as a vector of size n*(n+1)/2 Used to hide in...