Logo ROOT   6.10/09
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 
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 
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
double par[1]
Definition: unuranDistr.cxx:38
std::vector< MinuitParameter > fParameters
int FindIndex(const std::string &) const
void SetLowerLimit(unsigned int, double)
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21
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...
TRObject operator()(const T1 &t1) const
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
double result[121]
MnUserCovariance Int2extCovariance(const MnAlgebraicVector &, const MnAlgebraicSymMatrix &) const
const Int_t n
Definition: legend1.C:16
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...