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