Logo ROOT   6.18/05
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
13
14#include "Minuit2/MnPrint.h"
15
16
17namespace ROOT {
18
19 namespace Minuit2 {
20
21
22//
23// construct from user parameters (befor minimization)
24//
25MnUserParameterState::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//
45MnUserParameterState::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
59MnUserParameterState::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
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
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//
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
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
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
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
185const std::vector<MinuitParameter>& MnUserParameterState::MinuitParameters() const {
186 //access to parameters (row-wise)
187 return fParameters.Parameters();
188}
189
190std::vector<double> MnUserParameterState::Params() const {
191 //access to parameters in column-wise representation
192 return fParameters.Params();
193}
194std::vector<double> MnUserParameterState::Errors() const {
195 //access to errors in column-wise representation
196 return fParameters.Errors();
197}
198
200 //access to single Parameter i
201 return fParameters.Parameter(i);
202}
203
204void 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
228void 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
253void 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
263void 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 }
274 fGCCValid = false;
275}
276
278 // release parameter e (external index)
279 // no-op if parameter is const
280 if (Parameter(e).IsConst() ) return;
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
291void 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
303void MnUserParameterState::SetError(unsigned int e, double val) {
304 // set error for parameter e (external index)
305 fParameters.SetError(e, val);
306}
307
308void 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
324void MnUserParameterState::SetUpperLimit(unsigned int e, double up) {
325 // set upper limit for parameter e (external index)
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
338void MnUserParameterState::SetLowerLimit(unsigned int e, double low) {
339 // set lower limit for parameter e (external index)
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())
359}
360
361double MnUserParameterState::Value(unsigned int i) const {
362 // get value for parameter e (external index)
363 return fParameters.Value(i);
364}
365double 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
372void MnUserParameterState::Fix(const std::string & name) { Fix(Index(name));}
373
375
376void MnUserParameterState::SetValue(const std::string & name, double val) {SetValue(Index(name), val);}
377
378void MnUserParameterState::SetError(const std::string & name, double val) { SetError(Index(name), val);}
379
380void MnUserParameterState::SetLimits(const std::string & name, double low, double up) {SetLimits(Index(name), low, up);}
381
382void MnUserParameterState::SetUpperLimit(const std::string & name, double up) { SetUpperLimit(Index(name), up);}
383
384void MnUserParameterState::SetLowerLimit(const std::string & name, double low) {SetLowerLimit(Index(name), low);}
385
387
388double MnUserParameterState::Value(const std::string & name) const {return Value(Index(name));}
389
390double MnUserParameterState::Error(const std::string & name) const {return Error(Index(name));}
391
392
393unsigned int MnUserParameterState::Index(const std::string & name) const {
394 //convert name into external number of Parameter
395 return fParameters.Index(name);
396}
397
398const 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}
402const 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
409double MnUserParameterState::Int2ext(unsigned int i, double val) const {
410 // internal to external value
411 return fParameters.Trafo().Int2ext(i, val);
412}
413double MnUserParameterState::Ext2int(unsigned int e, double val) const {
414 // external to internal value
415 return fParameters.Trafo().Ext2int(e, val);
416}
417unsigned int MnUserParameterState::IntOfExt(unsigned int ext) const {
418 // return internal index for external index ext
419 return fParameters.Trafo().IntOfExt(ext);
420}
421unsigned 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
#define MN_INFO_MSG2(loc, str)
Definition: MnPrint.h:124
#define MN_INFO_MSG(str)
Definition: MnPrint.h:110
#define e(i)
Definition: RSha256.hxx:103
char name[80]
Definition: TGX11.cxx:109
double sqrt(double)
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
const double * Data() const
Definition: LASymMatrix.h:233
unsigned int Nrow() const
Definition: LASymMatrix.h:239
unsigned int size() const
Definition: LASymMatrix.h:237
const MnAlgebraicSymMatrix & InvHessian() const
Definition: MinimumError.h:60
const MnAlgebraicVector & Dirin() const
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
Definition: MinimumState.h:29
const MinimumError & Error() const
Definition: MinimumState.h:62
const MnAlgebraicVector & Vec() const
Definition: MinimumState.h:59
const MinimumParameters & Parameters() const
Definition: MinimumState.h:58
class for the individual Minuit Parameter with Name and number; contains the input numbers for the mi...
class to reduce the covariance matrix when a parameter is fixed by removing the corresponding row and...
class for global correlation coefficient
determines the relative floating point arithmetic precision.
Class containing the covariance matrix data represented as a vector of size n*(n+1)/2 Used to hide in...
const std::vector< double > & Data() const
const MnMachinePrecision & Precision() const
void SetLimits(unsigned int, double, double)
unsigned int Index(const std::string &) const
const std::string & GetName(unsigned int) const
double Int2ext(unsigned int, double) const
const MinuitParameter & Parameter(unsigned int i) const
double Ext2int(unsigned int, double) const
void Add(const std::string &name, double val, double err)
unsigned int ExtOfInt(unsigned int) const
const char * Name(unsigned int) const
MnUserParameterState()
default constructor (invalid state)
const std::vector< ROOT::Minuit2::MinuitParameter > & MinuitParameters() const
facade: forward interface of MnUserParameters and MnUserTransformation
unsigned int IntOfExt(unsigned int) const
void SetUpperLimit(unsigned int, double)
std::vector< double > Errors() const
void SetLowerLimit(unsigned int, double)
std::vector< double > Params() const
API class for the user interaction with the parameters; serves as input to the minimizer as well as o...
double Error(unsigned int) const
std::vector< double > Params() const
access to parameters and errors in column-wise representation
const char * Name(unsigned int) const
const MinuitParameter & Parameter(unsigned int) const
access to single Parameter
unsigned int Index(const std::string &) const
double Value(unsigned int) const
const MnMachinePrecision & Precision() const
void Fix(unsigned int)
interaction via external number of Parameter
void SetLowerLimit(unsigned int, double)
void SetError(unsigned int, double)
void SetValue(unsigned int, double)
const std::vector< ROOT::Minuit2::MinuitParameter > & Parameters() const
access to parameters (row-wise)
const MnUserTransformation & Trafo() const
std::vector< double > Errors() const
const std::string & GetName(unsigned int) const
void SetUpperLimit(unsigned int, double)
bool Add(const std::string &, double, double)
Add free Parameter Name, Value, Error.
void SetLimits(unsigned int, double, double)
class dealing with the transformation between user specified parameters (external) and internal param...
double Ext2int(unsigned int, double) const
MnUserCovariance Int2extCovariance(const MnAlgebraicVector &, const MnAlgebraicSymMatrix &) const
unsigned int ExtOfInt(unsigned int internal) const
const std::vector< MinuitParameter > & Parameters() const
unsigned int IntOfExt(unsigned int) const
double Int2ext(unsigned int, double) const
double Int2extError(unsigned int, double, double) const
int Invert(LASymMatrix &)
Definition: LaInverse.cxx:22
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21