Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
MnUserTransformation.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/MnMatrix.h"
13
14#include <algorithm>
15#include <cstdio>
16#include <string>
17#include <sstream>
18
19namespace ROOT {
20
21namespace Minuit2 {
22
23class MnParStr {
24
25public:
26 MnParStr(const std::string &name) : fName(name) {}
27
29
30 bool operator()(const MinuitParameter &par) const
31 {
32 // return (strcmp(par.Name(), fName) == 0);
33 return par.GetName() == fName;
34 }
35
36private:
37 const std::string &fName;
38};
39
40MnUserTransformation::MnUserTransformation(const std::vector<double> &par, const std::vector<double> &err)
41 : fPrecision(MnMachinePrecision()), fParameters(std::vector<MinuitParameter>()),
42 fExtOfInt(std::vector<unsigned int>()), fDoubleLimTrafo(SinParameterTransformation()),
43 fUpperLimTrafo(SqrtUpParameterTransformation()), fLowerLimTrafo(SqrtLowParameterTransformation()),
44 fCache(std::vector<double>())
45{
46 // constructor from a vector of parameter values and a vector of errors (step sizes)
47 // class has as data member the transformation objects (all of the types),
48 // the std::vector of MinuitParameter objects and the vector with the index conversions from
49 // internal to external (fExtOfInt)
50
51 fParameters.reserve(par.size());
52 fExtOfInt.reserve(par.size());
53 fCache.reserve(par.size());
54
55 std::string parName;
56 for (unsigned int i = 0; i < par.size(); i++) {
57 std::ostringstream buf;
58 buf << "p" << i;
59 parName = buf.str();
60 Add(parName, par[i], err[i]);
61 }
62}
63
64//#ifdef MINUIT2_THREAD_SAFE
65// this if a thread-safe implementation needed if want to share transformation object between the threads
66std::vector<double> MnUserTransformation::operator()(const MnAlgebraicVector &pstates) const
67{
68 // transform an internal Minuit vector of internal values in a std::vector of external values
69 // fixed parameters will have their fixed values
70 unsigned int n = pstates.size();
71 // need to initialize to the stored (initial values) parameter values for the fixed ones
72 std::vector<double> pcache(fCache);
73 for (unsigned int i = 0; i < n; i++) {
74 if (fParameters[fExtOfInt[i]].HasLimits()) {
75 pcache[fExtOfInt[i]] = Int2ext(i, pstates(i));
76 } else {
77 pcache[fExtOfInt[i]] = pstates(i);
78 }
79 }
80 return pcache;
81}
82
83// #else
84// const std::vector<double> & MnUserTransformation::operator()(const MnAlgebraicVector& pstates) const {
85// // transform an internal Minuit vector of internal values in a std::vector of external values
86// // std::vector<double> Cache(pstates.size() );
87// for(unsigned int i = 0; i < pstates.size(); i++) {
88// if(fParameters[fExtOfInt[i]].HasLimits()) {
89// fCache[fExtOfInt[i]] = Int2ext(i, pstates(i));
90// } else {
91// fCache[fExtOfInt[i]] = pstates(i);
92// }
93// }
94
95// return fCache;
96// }
97// #endif
98
99double MnUserTransformation::Int2ext(unsigned int i, double val) const
100{
101 // return external value from internal value for parameter i
102 if (fParameters[fExtOfInt[i]].HasLimits()) {
103 if (fParameters[fExtOfInt[i]].HasUpperLimit() && fParameters[fExtOfInt[i]].HasLowerLimit())
104 return fDoubleLimTrafo.Int2ext(val, fParameters[fExtOfInt[i]].UpperLimit(),
105 fParameters[fExtOfInt[i]].LowerLimit());
106 else if (fParameters[fExtOfInt[i]].HasUpperLimit() && !fParameters[fExtOfInt[i]].HasLowerLimit())
107 return fUpperLimTrafo.Int2ext(val, fParameters[fExtOfInt[i]].UpperLimit());
108 else
109 return fLowerLimTrafo.Int2ext(val, fParameters[fExtOfInt[i]].LowerLimit());
110 }
111
112 return val;
113}
114
115double MnUserTransformation::Int2extError(unsigned int i, double val, double err) const
116{
117 // return external error from internal error for parameter i
118
119 // err = sigma Value == std::sqrt(cov(i,i))
120 double dx = err;
121
122 if (fParameters[fExtOfInt[i]].HasLimits()) {
123 double ui = Int2ext(i, val);
124 double du1 = Int2ext(i, val + dx) - ui;
125 double du2 = Int2ext(i, val - dx) - ui;
126 if (fParameters[fExtOfInt[i]].HasUpperLimit() && fParameters[fExtOfInt[i]].HasLowerLimit()) {
127 // double al = fParameters[fExtOfInt[i]].Lower();
128 // double ba = fParameters[fExtOfInt[i]].Upper() - al;
129 // double du1 = al + 0.5*(sin(val + dx) + 1.)*ba - ui;
130 // double du2 = al + 0.5*(sin(val - dx) + 1.)*ba - ui;
131 // if(dx > 1.) du1 = ba;
132 if (dx > 1.)
133 du1 = fParameters[fExtOfInt[i]].UpperLimit() - fParameters[fExtOfInt[i]].LowerLimit();
134 dx = 0.5 * (std::fabs(du1) + std::fabs(du2));
135 } else {
136 dx = 0.5 * (std::fabs(du1) + std::fabs(du2));
137 }
138 }
139
140 return dx;
141}
142
145{
146 // return the external covariance matrix from the internal error matrix and the internal parameter value
147 // the vector of internal parameter is needed for the derivatives (Jacobian of the transformation)
148 // Vext(i,j) = Vint(i,j) * dPext(i)/dPint(i) * dPext(j)/dPint(j)
149
150 MnUserCovariance result(cov.Nrow());
151 for (unsigned int i = 0; i < vec.size(); i++) {
152 double dxdi = 1.;
153 if (fParameters[fExtOfInt[i]].HasLimits()) {
154 // dxdi = 0.5*std::fabs((fParameters[fExtOfInt[i]].Upper() -
155 // fParameters[fExtOfInt[i]].Lower())*cos(vec(i)));
156 dxdi = DInt2Ext(i, vec(i));
157 }
158 for (unsigned int j = i; j < vec.size(); j++) {
159 double dxdj = 1.;
160 if (fParameters[fExtOfInt[j]].HasLimits()) {
161 // dxdj = 0.5*std::fabs((fParameters[fExtOfInt[j]].Upper() -
162 // fParameters[fExtOfInt[j]].Lower())*cos(vec(j)));
163 dxdj = DInt2Ext(j, vec(j));
164 }
165 result(i, j) = dxdi * cov(i, j) * dxdj;
166 }
167 // double diag = Int2extError(i, vec(i), std::sqrt(cov(i,i)));
168 // result(i,i) = diag*diag;
169 }
170
171 return result;
172}
173
174double MnUserTransformation::Ext2int(unsigned int i, double val) const
175{
176 // return the internal value for parameter i with external value val
177
178 if (fParameters[i].HasLimits()) {
179 if (fParameters[i].HasUpperLimit() && fParameters[i].HasLowerLimit())
180 return fDoubleLimTrafo.Ext2int(val, fParameters[i].UpperLimit(), fParameters[i].LowerLimit(), Precision());
181 else if (fParameters[i].HasUpperLimit() && !fParameters[i].HasLowerLimit())
182 return fUpperLimTrafo.Ext2int(val, fParameters[i].UpperLimit(), Precision());
183 else
184 return fLowerLimTrafo.Ext2int(val, fParameters[i].LowerLimit(), Precision());
185 }
186
187 return val;
188}
189
190double MnUserTransformation::DInt2Ext(unsigned int i, double val) const
191{
192 // return the derivative of the int->ext transformation: dPext(i) / dPint(i)
193 // for the parameter i with value val
194
195 double dd = 1.;
196 if (fParameters[fExtOfInt[i]].HasLimits()) {
197 if (fParameters[fExtOfInt[i]].HasUpperLimit() && fParameters[fExtOfInt[i]].HasLowerLimit())
198 // dd = 0.5*std::fabs((fParameters[fExtOfInt[i]].Upper() -
199 // fParameters[fExtOfInt[i]].Lower())*cos(vec(i)));
200 dd = fDoubleLimTrafo.DInt2Ext(val, fParameters[fExtOfInt[i]].UpperLimit(),
201 fParameters[fExtOfInt[i]].LowerLimit());
202 else if (fParameters[fExtOfInt[i]].HasUpperLimit() && !fParameters[fExtOfInt[i]].HasLowerLimit())
203 dd = fUpperLimTrafo.DInt2Ext(val, fParameters[fExtOfInt[i]].UpperLimit());
204 else
205 dd = fLowerLimTrafo.DInt2Ext(val, fParameters[fExtOfInt[i]].LowerLimit());
206 }
207
208 return dd;
209}
210
211/*
212 double MnUserTransformation::dExt2Int(unsigned int, double) const {
213 double dd = 1.;
214
215 if(fParameters[fExtOfInt[i]].HasLimits()) {
216 if(fParameters[fExtOfInt[i]].HasUpperLimit() && fParameters[fExtOfInt[i]].HasLowerLimit())
217 // dd = 0.5*std::fabs((fParameters[fExtOfInt[i]].Upper() -
218 fParameters[fExtOfInt[i]].Lower())*cos(vec(i))); dd = fDoubleLimTrafo.dExt2Int(val,
219 fParameters[fExtOfInt[i]].UpperLimit(), fParameters[fExtOfInt[i]].LowerLimit()); else
220 if(fParameters[fExtOfInt[i]].HasUpperLimit() && !fParameters[fExtOfInt[i]].HasLowerLimit()) dd =
221 fUpperLimTrafo.dExt2Int(val, fParameters[fExtOfInt[i]].UpperLimit()); else dd = fLowerLimTrafo.dExtInt(val,
222 fParameters[fExtOfInt[i]].LowerLimit());
223 }
224
225 return dd;
226 }
227 */
228
229unsigned int MnUserTransformation::IntOfExt(unsigned int ext) const
230{
231 // return internal index given external one ext
232 assert(ext < fParameters.size());
233 assert(!fParameters[ext].IsFixed());
234 assert(!fParameters[ext].IsConst());
235 std::vector<unsigned int>::const_iterator iind = std::find(fExtOfInt.begin(), fExtOfInt.end(), ext);
236 assert(iind != fExtOfInt.end());
237
238 return (iind - fExtOfInt.begin());
239}
240
241std::vector<double> MnUserTransformation::Params() const
242{
243 // return std::vector of double with parameter values
244 unsigned int n = fParameters.size();
245 std::vector<double> result(n);
246 for (unsigned int i = 0; i < n; ++i)
247 result[i] = fParameters[i].Value();
248
249 return result;
250}
251
252std::vector<double> MnUserTransformation::Errors() const
253{
254 // return std::vector of double with parameter errors
255 std::vector<double> result;
256 result.reserve(fParameters.size());
257 for (std::vector<MinuitParameter>::const_iterator ipar = Parameters().begin(); ipar != Parameters().end(); ++ipar)
258 result.push_back((*ipar).Error());
259
260 return result;
261}
262
264{
265 // return the MinuitParameter object for index n (external)
266 assert(n < fParameters.size());
267 return fParameters[n];
268}
269
270// bool MnUserTransformation::Remove(const std::string & name) {
271// // remove parameter with name
272// // useful if want to re-define a parameter
273// // return false if parameter does not exist
274// std::vector<MinuitParameter>::iterator itr = std::find_if(fParameters.begin(), fParameters.end(), MnParStr(name)
275// ); if (itr == fParameters.end() ) return false; int n = itr - fParameters.begin(); if (n < 0 || n >=
276// fParameters.size() ) return false; fParameters.erase(itr); fCache.erase( fExtOfInt.begin() + n);
277// std::vector<unsigned int>::iterator iind = std::find(fExtOfInt.begin(), fExtOfInt.end(), n);
278// if (iind != fExtOfInt.end()) fExtOfInt.erase(iind);
279// }
280
281bool MnUserTransformation::Add(const std::string &name, double val, double err)
282{
283 // add a new unlimited parameter giving name, value and err (step size)
284 // return false if parameter already exists
285 if (std::find_if(fParameters.begin(), fParameters.end(), MnParStr(name)) != fParameters.end())
286 return false;
287 fExtOfInt.push_back(fParameters.size());
288 fCache.push_back(val);
289 fParameters.push_back(MinuitParameter(fParameters.size(), name, val, err));
290 return true;
291}
292
293bool MnUserTransformation::Add(const std::string &name, double val, double err, double low, double up)
294{
295 // add a new limited parameter giving name, value, err (step size) and lower/upper limits
296 // return false if parameter already exists
297 if (std::find_if(fParameters.begin(), fParameters.end(), MnParStr(name)) != fParameters.end())
298 return false;
299 fExtOfInt.push_back(fParameters.size());
300 fCache.push_back(val);
301 fParameters.push_back(MinuitParameter(fParameters.size(), name, val, err, low, up));
302 return true;
303}
304
305bool MnUserTransformation::Add(const std::string &name, double val)
306{
307 // add a new constant parameter giving name and value
308 // return false if parameter already exists
309 if (std::find_if(fParameters.begin(), fParameters.end(), MnParStr(name)) != fParameters.end())
310 return false;
311 fCache.push_back(val);
312 // costant parameter - do not add in list of internals (fExtOfInt)
313 fParameters.push_back(MinuitParameter(fParameters.size(), name, val));
314 return true;
315}
316
317void MnUserTransformation::Fix(unsigned int n)
318{
319 // fix parameter n (external index)
320 assert(n < fParameters.size());
321 std::vector<unsigned int>::iterator iind = std::find(fExtOfInt.begin(), fExtOfInt.end(), n);
322 if (iind != fExtOfInt.end())
323 fExtOfInt.erase(iind, iind + 1);
324 fParameters[n].Fix();
325}
326
328{
329 // release parameter n (external index)
330 assert(n < fParameters.size());
331 std::vector<unsigned int>::const_iterator iind = std::find(fExtOfInt.begin(), fExtOfInt.end(), n);
332 if (iind == fExtOfInt.end()) {
333 fExtOfInt.push_back(n);
334 std::sort(fExtOfInt.begin(), fExtOfInt.end());
335 }
336 fParameters[n].Release();
337}
338
339void MnUserTransformation::SetValue(unsigned int n, double val)
340{
341 // set value for parameter n (external index)
342 assert(n < fParameters.size());
343 fParameters[n].SetValue(val);
344 fCache[n] = val;
345}
346
347void MnUserTransformation::SetError(unsigned int n, double err)
348{
349 // set error for parameter n (external index)
350 assert(n < fParameters.size());
351 fParameters[n].SetError(err);
352}
353
354void MnUserTransformation::SetLimits(unsigned int n, double low, double up)
355{
356 // set limits (lower/upper) for parameter n (external index)
357 assert(n < fParameters.size());
358 assert(low != up);
359 fParameters[n].SetLimits(low, up);
360}
361
362void MnUserTransformation::SetUpperLimit(unsigned int n, double up)
363{
364 // set upper limit for parameter n (external index)
365 assert(n < fParameters.size());
366 fParameters[n].SetUpperLimit(up);
367}
368
369void MnUserTransformation::SetLowerLimit(unsigned int n, double lo)
370{
371 // set lower limit for parameter n (external index)
372 assert(n < fParameters.size());
373 fParameters[n].SetLowerLimit(lo);
374}
375
377{
378 // remove limits for parameter n (external index)
379 assert(n < fParameters.size());
380 fParameters[n].RemoveLimits();
381}
382
383void MnUserTransformation::SetName(unsigned int n, const std::string &name)
384{
385 // set name for parameter n (external index)
386 assert(n < fParameters.size());
387 fParameters[n].SetName(name);
388}
389
390double MnUserTransformation::Value(unsigned int n) const
391{
392 // get value for parameter n (external index)
393 assert(n < fParameters.size());
394 return fParameters[n].Value();
395}
396
397double MnUserTransformation::Error(unsigned int n) const
398{
399 // get error for parameter n (external index)
400 assert(n < fParameters.size());
401 return fParameters[n].Error();
402}
403
404// interface by parameter name
405
406void MnUserTransformation::Fix(const std::string &name)
407{
408 // fix parameter
409 Fix(Index(name));
410}
411
412void MnUserTransformation::Release(const std::string &name)
413{
414 // release parameter
416}
417
418void MnUserTransformation::SetValue(const std::string &name, double val)
419{
420 // set value for parameter
421 SetValue(Index(name), val);
422}
423
424void MnUserTransformation::SetError(const std::string &name, double err)
425{
426 // set error
427 SetError(Index(name), err);
428}
429
430void MnUserTransformation::SetLimits(const std::string &name, double low, double up)
431{
432 // set lower/upper limits
433 SetLimits(Index(name), low, up);
434}
435
436void MnUserTransformation::SetUpperLimit(const std::string &name, double up)
437{
438 // set upper limit
440}
441
442void MnUserTransformation::SetLowerLimit(const std::string &name, double lo)
443{
444 // set lower limit
446}
447
449{
450 // remove limits
452}
453
454double MnUserTransformation::Value(const std::string &name) const
455{
456 // get parameter value
457 return Value(Index(name));
458}
459
460double MnUserTransformation::Error(const std::string &name) const
461{
462 // get parameter error
463 return Error(Index(name));
464}
465
466unsigned int MnUserTransformation::Index(const std::string &name) const
467{
468 // get index (external) corresponding to name
469 std::vector<MinuitParameter>::const_iterator ipar =
470 std::find_if(fParameters.begin(), fParameters.end(), MnParStr(name));
471 assert(ipar != fParameters.end());
472 // return (ipar - fParameters.begin());
473 return (*ipar).Number();
474}
475
476int MnUserTransformation::FindIndex(const std::string &name) const
477{
478 // find index (external) corresponding to name - return -1 if not found
479 std::vector<MinuitParameter>::const_iterator ipar =
480 std::find_if(fParameters.begin(), fParameters.end(), MnParStr(name));
481 if (ipar == fParameters.end())
482 return -1;
483 return (*ipar).Number();
484}
485
486const std::string &MnUserTransformation::GetName(unsigned int n) const
487{
488 // get name corresponding to index (external)
489 assert(n < fParameters.size());
490 return fParameters[n].GetName();
491}
492
493const char *MnUserTransformation::Name(unsigned int n) const
494{
495 // get name corresponding to index (external)
496 return GetName(n).c_str();
497}
498
499} // namespace Minuit2
500
501} // namespace ROOT
double
char name[80]
Definition TGX11.cxx:110
Class describing a symmetric matrix of size n.
Definition LASymMatrix.h:45
unsigned int Nrow() const
unsigned int size() const
Definition LAVector.h:227
class for the individual Minuit Parameter with Name and number; contains the input numbers for the mi...
const std::string & GetName() const
Sets the relative floating point (double) arithmetic precision.
MnParStr(const std::string &name)
bool operator()(const MinuitParameter &par) const
Class containing the covariance matrix data represented as a vector of size n*(n+1)/2 Used to hide in...
double Ext2int(unsigned int, double) const
MnUserCovariance Int2extCovariance(const MnAlgebraicVector &, const MnAlgebraicSymMatrix &) const
bool Add(const std::string &, double, double)
double DInt2Ext(unsigned int, double) const
void SetName(unsigned int, const std::string &)
std::vector< MinuitParameter > fParameters
std::vector< double > operator()(const MnAlgebraicVector &) const
const std::vector< MinuitParameter > & Parameters() const
unsigned int IntOfExt(unsigned int) const
const char * Name(unsigned int) const
std::vector< unsigned int > fExtOfInt
std::vector< double > Params() const
access to parameters and errors in column-wise representation
SqrtLowParameterTransformation fLowerLimTrafo
int FindIndex(const std::string &) const
double Int2ext(unsigned int, double) const
unsigned int Index(const std::string &) const
SinParameterTransformation fDoubleLimTrafo
const std::string & GetName(unsigned int) const
void SetLimits(unsigned int, double, double)
double Int2extError(unsigned int, double, double) const
const MnMachinePrecision & Precision() const
forwarded interface
SqrtUpParameterTransformation fUpperLimTrafo
const MinuitParameter & Parameter(unsigned int) const
class for the transformation for double-limited parameter Using a sin function one goes from a double...
long double DInt2Ext(long double Value, long double Upper, long double Lower) const
long double Int2ext(long double Value, long double Upper, long double Lower) const
long double Ext2int(long double Value, long double Upper, long double Lower, const MnMachinePrecision &) const
Transformation from external to internal Parameter based on sqrt(1 + x**2)
long double Ext2int(long double Value, long double Lower, const MnMachinePrecision &) const
long double DInt2Ext(long double Value, long double Lower) const
long double Int2ext(long double Value, long double Lower) const
Transformation from external to internal Parameter based on sqrt(1 + x**2)
long double Int2ext(long double Value, long double Upper) const
long double DInt2Ext(long double Value, long double Upper) const
long double Ext2int(long double Value, long double Upper, const MnMachinePrecision &) const
const Int_t n
Definition legend1.C:16
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...