Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 
10 #include "Minuit2/MnUserParameterState.h"
11 #include "Minuit2/MnCovarianceSqueeze.h"
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 
29 MnUserParameterState::MnUserParameterState(const MnUserParameters& par) :
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);
56  assert(fCovariance.Nrow() == VariableParameters());
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);
70  assert(fCovariance.Nrow() == VariableParameters());
71 }
72 
73 
74 MnUserParameterState::MnUserParameterState(const MnUserParameters& par, const MnUserCovariance& cov) :
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  }
88  assert(fCovariance.Nrow() == VariableParameters());
89  //
90  // need to Fix that in case of limited parameters
91  // fIntCovariance = MnUserCovariance();
92  //
93 }
94 
95 //
96 //
97 MnUserParameterState::MnUserParameterState(const MinimumState& st, double up, const MnUserTransformation& trafo) :
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.Error().IsValid() ? 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.Error().IsValid() ? 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);
148  fGlobalCC = MnGlobalCorrelationCoeff(st.Error().InvHessian());
149  fGCCValid = fGlobalCC.IsValid();
150 
151  assert(fCovariance.Nrow() == VariableParameters());
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 
160 MnUserCovariance MnUserParameterState::Hessian() const {
161  // invert covariance matrix and return Hessian
162  // need to copy in a MnSymMatrix
163  MnAlgebraicSymMatrix mat(fCovariance.Nrow() );
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
170  MnUserCovariance tmp(fCovariance.Nrow());
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) {
268  fCovariance = MnCovarianceSqueeze()(fCovariance, i);
269  fIntCovariance = MnCovarianceSqueeze()(fIntCovariance, i);
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)
316  fIntParameters[i] = Ext2int(e, fIntParameters[i]);
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)
332  fIntParameters[i] = Ext2int(e, fIntParameters[i]);
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])
346  fIntParameters[i] = Ext2int(e, fIntParameters[i]);
347  else
348  fIntParameters[i] = Ext2int(e, low + 0.1 * Parameter(e).Error() );
349  }
350 }
351 
352 void MnUserParameterState::RemoveLimits(unsigned int e) {
353  // remove limit for parameter e (external index)
354  fParameters.RemoveLimits(e);
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 }
425 unsigned int MnUserParameterState::VariableParameters() const {
426  // return number of variable parameters
427  return fParameters.Trafo().VariableParameters();
428 }
429 const MnMachinePrecision& MnUserParameterState::Precision() const {
430  // return global parameter precision
431  return fParameters.Precision();
432 }
433 
434 void MnUserParameterState::SetPrecision(double eps) {
435  // set global parameter precision
436  fParameters.SetPrecision(eps);
437 }
438 
439  } // namespace Minuit2
440 
441 } // namespace ROOT