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