Logo ROOT   6.18/05
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
19namespace ROOT {
20
21 namespace Minuit2 {
22
23
24class MnParStr {
25
26public:
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
37private:
38 const std::string & fName;
39};
40
41
42MnUserTransformation::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
63std::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
95double 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
109double 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
162double 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
177double 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
213unsigned 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
224std::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
234std::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
264bool 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
275bool 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
286bool 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
297void 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
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
317void 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
324void 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
330void 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
337void 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
343void 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
355void 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
362double 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
368double 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
376void MnUserTransformation::Fix(const std::string & name) {
377 // fix parameter
378 Fix(Index(name));
379}
380
381void MnUserTransformation::Release(const std::string & name) {
382 // release parameter
384}
385
386void MnUserTransformation::SetValue(const std::string & name, double val) {
387 // set value for parameter
388 SetValue(Index(name), val);
389}
390
391void MnUserTransformation::SetError(const std::string & name, double err) {
392 // set error
393 SetError(Index(name), err);
394}
395
396void MnUserTransformation::SetLimits(const std::string & name, double low, double up) {
397 // set lower/upper limits
398 SetLimits(Index(name), low, up);
399}
400
401void MnUserTransformation::SetUpperLimit(const std::string & name, double up) {
402 // set upper limit
404}
405
406void MnUserTransformation::SetLowerLimit(const std::string & name, double lo) {
407 // set lower limit
409}
410
411void MnUserTransformation::RemoveLimits(const std::string & name) {
412 // remove limits
414}
415
416double MnUserTransformation::Value(const std::string & name) const {
417 // get parameter value
418 return Value(Index(name));
419}
420
421double MnUserTransformation::Error(const std::string & name) const {
422 // get parameter error
423 return Error(Index(name));
424}
425
426unsigned 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
435int 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
444const 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
450const 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
char name[80]
Definition: TGX11.cxx:109
TRObject operator()(const T1 &t1) const
Class describing a symmetric matrix of size n.
Definition: LASymMatrix.h:51
unsigned int Nrow() const
Definition: LASymMatrix.h:239
unsigned int size() const
Definition: LAVector.h:198
class for the individual Minuit Parameter with Name and number; contains the input numbers for the mi...
determines the relative floating point arithmetic precision.
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
void SetUpperLimit(unsigned int, double)
bool Add(const std::string &, double, double)
double DInt2Ext(unsigned int, double) const
void SetLowerLimit(unsigned int, double)
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
std::vector< double > Errors() 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...
double DInt2Ext(double Value, double Upper, double Lower) const
double Int2ext(double Value, double Upper, double Lower) const
double Ext2int(double Value, double Upper, double Lower, const MnMachinePrecision &) const
Transformation from external to internal Parameter based on sqrt(1 + x**2)
double Ext2int(double Value, double Lower, const MnMachinePrecision &) const
Transformation from external to internal Parameter based on sqrt(1 + x**2)
double DInt2Ext(double Value, double Upper) const
double Ext2int(double Value, double Upper, const MnMachinePrecision &) const
double Int2ext(double Value, double Upper) const
const Int_t n
Definition: legend1.C:16
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Namespace for new ROOT classes and functions.
Definition: StringConv.hxx:21