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