Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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#include "Minuit2/MnPrint.h"
14
15namespace ROOT {
16
17namespace Minuit2 {
18
19//
20// construct from user parameters (before minimization)
21//
22MnUserParameterState::MnUserParameterState(const std::vector<double> &par, const std::vector<double> &err)
23 : fValid(true), fCovarianceValid(false), fGCCValid(false), fCovStatus(-1), fFVal(0.), fEDM(0.), fNFcn(0),
24 fParameters(MnUserParameters(par, err)), fCovariance(MnUserCovariance()), fGlobalCC(MnGlobalCorrelationCoeff()),
25 fIntParameters(par), fIntCovariance(MnUserCovariance())
26{
27}
28
30 : fValid(true), fCovarianceValid(false), fGCCValid(false), fCovStatus(-1), fFVal(0.), fEDM(0.), fNFcn(0),
31 fParameters(par), fCovariance(MnUserCovariance()), fGlobalCC(MnGlobalCorrelationCoeff()),
32 fIntParameters(std::vector<double>()), fIntCovariance(MnUserCovariance())
33{
34 // construct from user parameters (before minimization)
35
36 for (std::vector<MinuitParameter>::const_iterator ipar = MinuitParameters().begin();
37 ipar != MinuitParameters().end(); ++ipar) {
38 if ((*ipar).IsConst() || (*ipar).IsFixed())
39 continue;
40 if ((*ipar).HasLimits())
41 fIntParameters.push_back(Ext2int((*ipar).Number(), (*ipar).Value()));
42 else
43 fIntParameters.push_back((*ipar).Value());
44 }
45}
46
47//
48// construct from user parameters + errors (before minimization)
49//
50MnUserParameterState::MnUserParameterState(const std::vector<double> &par, const std::vector<double> &cov,
51 unsigned int nrow)
52 : fValid(true), fCovarianceValid(true), fGCCValid(false), fCovStatus(-1), fFVal(0.), fEDM(0.), fNFcn(0),
53 fParameters(MnUserParameters()), fCovariance(MnUserCovariance(cov, nrow)), fGlobalCC(MnGlobalCorrelationCoeff()),
54 fIntParameters(par), fIntCovariance(MnUserCovariance(cov, nrow))
55{
56 // construct from user parameters + errors (before minimization) using std::vector for parameter error and // an
57 // std::vector of size n*(n+1)/2 for the covariance matrix and n (rank of cov matrix)
58
59 std::vector<double> err;
60 err.reserve(par.size());
61 for (unsigned int i = 0; i < par.size(); i++) {
62 assert(fCovariance(i, i) > 0.);
63 err.push_back(std::sqrt(fCovariance(i, i)));
64 }
65 fParameters = MnUserParameters(par, err);
66 assert(fCovariance.Nrow() == VariableParameters());
67}
68
69MnUserParameterState::MnUserParameterState(const std::vector<double> &par, const MnUserCovariance &cov)
70 : fValid(true), fCovarianceValid(true), fGCCValid(false), fCovStatus(-1), fFVal(0.), fEDM(0.), fNFcn(0),
71 fParameters(MnUserParameters()), fCovariance(cov), fGlobalCC(MnGlobalCorrelationCoeff()), fIntParameters(par),
72 fIntCovariance(cov)
73{
74 // construct from user parameters + errors (before minimization) using std::vector (params) and MnUserCovariance
75 // class
76
77 std::vector<double> err;
78 err.reserve(par.size());
79 for (unsigned int i = 0; i < par.size(); i++) {
80 assert(fCovariance(i, i) > 0.);
81 err.push_back(std::sqrt(fCovariance(i, i)));
82 }
83 fParameters = MnUserParameters(par, err);
84 assert(fCovariance.Nrow() == VariableParameters());
85}
86
88 : fValid(true), fCovarianceValid(true), fGCCValid(false), fCovStatus(-1), fFVal(0.), fEDM(0.), fNFcn(0),
89 fParameters(par), fCovariance(cov), fGlobalCC(MnGlobalCorrelationCoeff()), fIntParameters(std::vector<double>()),
90 fIntCovariance(cov)
91{
92 // construct from user parameters + errors (before minimization) using
93 // MnUserParameters and MnUserCovariance objects
94
96 for (std::vector<MinuitParameter>::const_iterator ipar = MinuitParameters().begin();
97 ipar != MinuitParameters().end(); ++ipar) {
98 if ((*ipar).IsConst() || (*ipar).IsFixed())
99 continue;
100 if ((*ipar).HasLimits())
101 fIntParameters.push_back(Ext2int((*ipar).Number(), (*ipar).Value()));
102 else
103 fIntParameters.push_back((*ipar).Value());
104 }
105 assert(fCovariance.Nrow() == VariableParameters());
106 //
107 // need to Fix that in case of limited parameters
108 // fIntCovariance = MnUserCovariance();
109 //
110}
111
112//
113//
115 : fValid(st.IsValid()), fCovarianceValid(false), fGCCValid(false), fCovStatus(-1), fFVal(st.Fval()), fEDM(st.Edm()),
116 fNFcn(st.NFcn()), fParameters(MnUserParameters()), fCovariance(MnUserCovariance()),
117 fGlobalCC(MnGlobalCorrelationCoeff()), fIntParameters(std::vector<double>()), fIntCovariance(MnUserCovariance())
118{
119 //
120 // construct from internal parameters (after minimization)
121 //
122 // std::cout << "build a MnUserParameterState after minimization.." << std::endl;
123
124 for (std::vector<MinuitParameter>::const_iterator ipar = trafo.Parameters().begin();
125 ipar != trafo.Parameters().end(); ++ipar) {
126 if ((*ipar).IsConst()) {
127 Add((*ipar).GetName(), (*ipar).Value());
128 } else if ((*ipar).IsFixed()) {
129 Add((*ipar).GetName(), (*ipar).Value(), (*ipar).Error());
130 if ((*ipar).HasLimits()) {
131 if ((*ipar).HasLowerLimit() && (*ipar).HasUpperLimit())
132 SetLimits((*ipar).GetName(), (*ipar).LowerLimit(), (*ipar).UpperLimit());
133 else if ((*ipar).HasLowerLimit() && !(*ipar).HasUpperLimit())
134 SetLowerLimit((*ipar).GetName(), (*ipar).LowerLimit());
135 else
136 SetUpperLimit((*ipar).GetName(), (*ipar).UpperLimit());
137 }
138 Fix((*ipar).GetName());
139 } else if ((*ipar).HasLimits()) {
140 unsigned int i = trafo.IntOfExt((*ipar).Number());
141 double err =
142 st.Error().IsValid() ? std::sqrt(2. * up * st.Error().InvHessian()(i, i)) : st.Parameters().Dirin()(i);
143 Add((*ipar).GetName(), trafo.Int2ext(i, st.Vec()(i)), trafo.Int2extError(i, st.Vec()(i), err));
144 if ((*ipar).HasLowerLimit() && (*ipar).HasUpperLimit())
145 SetLimits((*ipar).GetName(), (*ipar).LowerLimit(), (*ipar).UpperLimit());
146 else if ((*ipar).HasLowerLimit() && !(*ipar).HasUpperLimit())
147 SetLowerLimit((*ipar).GetName(), (*ipar).LowerLimit());
148 else
149 SetUpperLimit((*ipar).GetName(), (*ipar).UpperLimit());
150 } else {
151 unsigned int i = trafo.IntOfExt((*ipar).Number());
152 double err =
153 st.Error().IsValid() ? std::sqrt(2. * up * st.Error().InvHessian()(i, i)) : st.Parameters().Dirin()(i);
154 Add((*ipar).GetName(), st.Vec()(i), err);
155 }
156 }
157
158 // need to be set afterwards because becore the ::Add method set fCovarianceValid to false
160
161 fCovStatus = -1; // when not available
162 // if (st.Error().HesseFailed() || st.Error().InvertFailed() ) fCovStatus = -1;
163 // when available
164 if (st.Error().IsAvailable())
165 fCovStatus = 0;
166
167 if (fCovarianceValid) {
168 fCovariance = trafo.Int2extCovariance(st.Vec(), st.Error().InvHessian());
170 MnUserCovariance(std::vector<double>(st.Error().InvHessian().Data(),
171 st.Error().InvHessian().Data() + st.Error().InvHessian().size()),
172 st.Error().InvHessian().Nrow());
173 fCovariance.Scale(2. * up);
176
177 assert(fCovariance.Nrow() == VariableParameters());
178
179 if (!st.Error().IsNotPosDef())
180 fCovStatus = 1; // when is valid and not NotPosDef
181 }
182 if (st.Error().IsMadePosDef())
183 fCovStatus = 2;
184 if (st.Error().IsAccurate())
185 fCovStatus = 3;
186
187}
188
190{
191 // invert covariance matrix and return Hessian
192 // need to copy in a MnSymMatrix
193 MnPrint print("MnUserParameterState::Hessian");
194
196 std::copy(fCovariance.Data().begin(), fCovariance.Data().end(), mat.Data());
197 int ifail = Invert(mat);
198 if (ifail != 0) {
199 print.Warn("Inversion failed; return diagonal matrix");
200
202 for (unsigned int i = 0; i < fCovariance.Nrow(); i++) {
203 tmp(i, i) = 1. / fCovariance(i, i);
204 }
205 return tmp;
206 }
207
208 MnUserCovariance hessian(mat.Data(), fCovariance.Nrow());
209 return hessian;
210}
211
212// facade: forward interface of MnUserParameters and MnUserTransformation
213// via MnUserParameterState
214
215const std::vector<MinuitParameter> &MnUserParameterState::MinuitParameters() const
216{
217 // access to parameters (row-wise)
218 return fParameters.Parameters();
219}
220
221std::vector<double> MnUserParameterState::Params() const
222{
223 // access to parameters in column-wise representation
224 return fParameters.Params();
225}
226std::vector<double> MnUserParameterState::Errors() const
227{
228 // access to errors in column-wise representation
229 return fParameters.Errors();
230}
231
233{
234 // access to single Parameter i
235 return fParameters.Parameter(i);
236}
237
238void MnUserParameterState::Add(const std::string &name, double val, double err)
239{
240 MnPrint print("MnUserParameterState::Add");
241
242 // add free Parameter
243 if (fParameters.Add(name, val, err)) {
244 fIntParameters.push_back(val);
245 fCovarianceValid = false;
246 fGCCValid = false;
247 fValid = true;
248 } else {
249 // redefine an existing parameter
250 int i = Index(name);
251 SetValue(i, val);
252 if (Parameter(i).IsConst()) {
253 print.Warn("Cannot modify status of constant parameter", name);
254 return;
255 }
256 SetError(i, err);
257 // release if it was fixed
258 if (Parameter(i).IsFixed())
259 Release(i);
260 }
261}
262
263void MnUserParameterState::Add(const std::string &name, double val, double err, double low, double up)
264{
265 MnPrint print("MnUserParameterState::Add");
266
267 // add limited Parameter
268 if (fParameters.Add(name, val, err, low, up)) {
269 fCovarianceValid = false;
270 fIntParameters.push_back(Ext2int(Index(name), val));
271 fGCCValid = false;
272 fValid = true;
273 } else { // Parameter already exist - just set values
274 int i = Index(name);
275 SetValue(i, val);
276 if (Parameter(i).IsConst()) {
277 print.Warn("Cannot modify status of constant parameter", name);
278 return;
279 }
280 SetError(i, err);
281 SetLimits(i, low, up);
282 // release if it was fixed
283 if (Parameter(i).IsFixed())
284 Release(i);
285 }
286}
287
288void MnUserParameterState::Add(const std::string &name, double val)
289{
290 // add const Parameter
291 if (fParameters.Add(name, val))
292 fValid = true;
293 else
294 SetValue(name, val);
295}
296
297// interaction via external number of Parameter
298
299void MnUserParameterState::Fix(unsigned int e)
300{
301 // fix parameter e (external index)
302 if (!Parameter(e).IsFixed() && !Parameter(e).IsConst()) {
303 unsigned int i = IntOfExt(e);
304 if (fCovarianceValid) {
305 // when covariance is valid remove fixed parameter row/column in covariance matrix
306 // use squeeze to remove first from Hessian and obtain then new covariance
309 }
310 fIntParameters.erase(fIntParameters.begin() + i, fIntParameters.begin() + i + 1);
311 }
313 fGCCValid = false;
314}
315
317{
318 // release parameter e (external index)
319 // no-op if parameter is const or if it is not fixed
320 if (Parameter(e).IsConst() || !Parameter(e).IsFixed())
321 return;
323 fCovarianceValid = false;
324 fGCCValid = false;
325 unsigned int i = IntOfExt(e);
326 if (Parameter(e).HasLimits())
327 fIntParameters.insert(fIntParameters.begin() + i, Ext2int(e, Parameter(e).Value()));
328 else
329 fIntParameters.insert(fIntParameters.begin() + i, Parameter(e).Value());
330}
331
332void MnUserParameterState::SetValue(unsigned int e, double val)
333{
334 // set value for parameter e ( external index )
335 fParameters.SetValue(e, val);
336 if (!Parameter(e).IsFixed() && !Parameter(e).IsConst()) {
337 unsigned int i = IntOfExt(e);
338 if (Parameter(e).HasLimits())
339 fIntParameters[i] = Ext2int(e, val);
340 else
341 fIntParameters[i] = val;
342 }
343}
344
345void MnUserParameterState::SetError(unsigned int e, double val)
346{
347 // set error for parameter e (external index)
348 fParameters.SetError(e, val);
349}
350
351void MnUserParameterState::SetLimits(unsigned int e, double low, double up)
352{
353 // set limits for parameter e (external index)
354 fParameters.SetLimits(e, low, up);
355 fCovarianceValid = false;
356 fGCCValid = false;
357 if (!Parameter(e).IsFixed() && !Parameter(e).IsConst()) {
358 unsigned int i = IntOfExt(e);
359 if (low < fIntParameters[i] && fIntParameters[i] < up)
361 else if (low >= fIntParameters[i])
362 fIntParameters[i] = Ext2int(e, low + 0.1 * Parameter(e).Error());
363 else
364 fIntParameters[i] = Ext2int(e, up - 0.1 * Parameter(e).Error());
365 }
366}
367
368void MnUserParameterState::SetUpperLimit(unsigned int e, double up)
369{
370 // set upper limit for parameter e (external index)
372 fCovarianceValid = false;
373 fGCCValid = false;
374 if (!Parameter(e).IsFixed() && !Parameter(e).IsConst()) {
375 unsigned int i = IntOfExt(e);
376 if (fIntParameters[i] < up)
378 else
379 fIntParameters[i] = Ext2int(e, up - 0.1 * Parameter(e).Error());
380 }
381}
382
383void MnUserParameterState::SetLowerLimit(unsigned int e, double low)
384{
385 // set lower limit for parameter e (external index)
387 fCovarianceValid = false;
388 fGCCValid = false;
389 if (!Parameter(e).IsFixed() && !Parameter(e).IsConst()) {
390 unsigned int i = IntOfExt(e);
391 if (low < fIntParameters[i])
393 else
394 fIntParameters[i] = Ext2int(e, low + 0.1 * Parameter(e).Error());
395 }
396}
397
399{
400 // remove limit for parameter e (external index)
402 fCovarianceValid = false;
403 fGCCValid = false;
404 if (!Parameter(e).IsFixed() && !Parameter(e).IsConst())
406}
407
408double MnUserParameterState::Value(unsigned int i) const
409{
410 // get value for parameter e (external index)
411 return fParameters.Value(i);
412}
413double MnUserParameterState::Error(unsigned int i) const
414{
415 // get error for parameter e (external index)
416 return fParameters.Error(i);
417}
418
419// interaction via name of Parameter
420
421void MnUserParameterState::Fix(const std::string &name)
422{
423 Fix(Index(name));
424}
425
426void MnUserParameterState::Release(const std::string &name)
427{
429}
430
431void MnUserParameterState::SetValue(const std::string &name, double val)
432{
433 SetValue(Index(name), val);
434}
435
436void MnUserParameterState::SetError(const std::string &name, double val)
437{
438 SetError(Index(name), val);
439}
440
441void MnUserParameterState::SetLimits(const std::string &name, double low, double up)
442{
443 SetLimits(Index(name), low, up);
444}
445
446void MnUserParameterState::SetUpperLimit(const std::string &name, double up)
447{
449}
450
451void MnUserParameterState::SetLowerLimit(const std::string &name, double low)
452{
453 SetLowerLimit(Index(name), low);
454}
455
457{
459}
460
461double MnUserParameterState::Value(const std::string &name) const
462{
463 return Value(Index(name));
464}
465
466double MnUserParameterState::Error(const std::string &name) const
467{
468 return Error(Index(name));
469}
470
471unsigned int MnUserParameterState::Index(const std::string &name) const
472{
473 // convert name into external number of Parameter
474 return fParameters.Index(name);
475}
476
477const char *MnUserParameterState::Name(unsigned int i) const
478{
479 // convert external number into name of Parameter (API returning a const char *)
480 return fParameters.Name(i);
481}
482const std::string &MnUserParameterState::GetName(unsigned int i) const
483{
484 // convert external number into name of Parameter (new interface returning a string)
485 return fParameters.GetName(i);
486}
487
488// transformation internal <-> external (forward to transformation class)
489
490double MnUserParameterState::Int2ext(unsigned int i, double val) const
491{
492 // internal to external value
493 return fParameters.Trafo().Int2ext(i, val);
494}
495double MnUserParameterState::Ext2int(unsigned int e, double val) const
496{
497 // external to internal value
498 return fParameters.Trafo().Ext2int(e, val);
499}
500unsigned int MnUserParameterState::IntOfExt(unsigned int ext) const
501{
502 // return internal index for external index ext
503 return fParameters.Trafo().IntOfExt(ext);
504}
505unsigned int MnUserParameterState::ExtOfInt(unsigned int internal) const
506{
507 // return external index for internal index internal
508 return fParameters.Trafo().ExtOfInt(internal);
509}
511{
512 // return number of variable parameters
514}
516{
517 // return global parameter precision
518 return fParameters.Precision();
519}
520
522{
523 // set global parameter precision
525}
526
527} // namespace Minuit2
528
529} // namespace ROOT
#define e(i)
Definition RSha256.hxx:103
char name[80]
Definition TGX11.cxx:110
Class describing a symmetric matrix of size n.
Definition LASymMatrix.h:45
const double * Data() const
unsigned int Nrow() const
unsigned int size() const
const MnAlgebraicSymMatrix & InvHessian() const
const MnAlgebraicVector & Dirin() const
MinimumState keeps the information (position, Gradient, 2nd deriv, etc) after one minimization step (...
const MinimumError & Error() const
const MnAlgebraicVector & Vec() const
const MinimumParameters & Parameters() const
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
Sets the relative floating point (double) arithmetic precision.
void Warn(const Ts &... args)
Definition MnPrint.h:135
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
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
CPyCppyy::Parameter Parameter
int Invert(LASymMatrix &)
Definition LaInverse.cxx:21
This file contains a specialised ROOT message handler to test for diagnostic in unit tests.