Logo ROOT  
Reference Guide
FlexibleInterpVar.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id: cranmer $
2// Author: Kyle Cranmer, Akira Shibata
3// Author: Giovanni Petrucciani (UCSD) (log-interpolation)
4/*************************************************************************
5 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12////////////////////////////////////////////////////////////////////////////////
13
14/** \class RooStats::HistFactory::FlexibleInterpVar
15 * \ingroup HistFactory
16 */
17
18#include "Riostream.h"
19#include <math.h>
20#include "TMath.h"
21
22#include "RooAbsReal.h"
23#include "RooRealVar.h"
24#include "RooArgList.h"
25#include "RooMsgService.h"
26#include "RooTrace.h"
27
29
30using namespace std;
31
33
34using namespace RooStats;
35using namespace HistFactory;
36
37////////////////////////////////////////////////////////////////////////////////
38/// Default constructor
39
40FlexibleInterpVar::FlexibleInterpVar()
41{
42 _nominal = 0;
44 _logInit = false ;
46}
47
48
49////////////////////////////////////////////////////////////////////////////////
50
51FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title,
52 const RooArgList& paramList,
53 double argNominal, std::vector<double> lowVec, std::vector<double> highVec) :
54 RooAbsReal(name, title),
55 _paramList("paramList","List of paramficients",this),
56 _nominal(argNominal), _low(lowVec), _high(highVec), _interpBoundary(1.)
57{
58 _logInit = false ;
59
60 for (auto param : paramList) {
61 if (!dynamic_cast<RooAbsReal*>(param)) {
62 coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") ERROR: paramficient " << param->GetName()
63 << " is not of type RooAbsReal" << endl ;
64 R__ASSERT(0) ;
65 }
66 _paramList.add(*param) ;
67 _interpCode.push_back(0); // default code
68 }
69 if (int(_low.size() ) != _paramList.getSize() || _low.size() != _high.size()) {
70 coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") invalid input low/high vectors " << endl;
71 R__ASSERT(int(_low.size() ) == _paramList.getSize());
72 R__ASSERT(_low.size() == _high.size());
73 }
74
76}
77
78////////////////////////////////////////////////////////////////////////////////
79
80FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title,
81 const RooArgList& paramList,
82 double argNominal, const RooArgList& lowList, const RooArgList& highList) :
83 RooAbsReal(name, title),
84 _paramList("paramList","List of paramficients",this),
85 _nominal(argNominal), _interpBoundary(1.)
86{
87 for (auto const *val : static_range_cast<RooAbsReal *>(lowList)){
88 _low.push_back(val->getVal()) ;
89 }
90
91 for (auto const *val : static_range_cast<RooAbsReal *>(highList)) {
92 _high.push_back(val->getVal()) ;
93 }
94
95 _logInit = false ;
96
97 for (auto param : paramList) {
98 if (!dynamic_cast<RooAbsReal*>(param)) {
99 coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") ERROR: paramficient " << param->GetName()
100 << " is not of type RooAbsReal" << endl ;
101 R__ASSERT(0) ;
102 }
103 _paramList.add(*param) ;
104 _interpCode.push_back(0); // default code
105 }
106 if (int(_low.size() ) != _paramList.getSize() || _low.size() != _high.size()) {
107 coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") invalid input low/high lists " << endl;
108 R__ASSERT(int(_low.size() ) == _paramList.getSize());
109 R__ASSERT(_low.size() == _high.size());
110 }
111
113}
114
115
116
117
118////////////////////////////////////////////////////////////////////////////////
119
120FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title,
121 const RooArgList& paramList,
122 double argNominal, vector<double> lowVec, vector<double> highVec,
123 vector<int> code) :
124 RooAbsReal(name, title),
125 _paramList("paramList","List of paramficients",this),
126 _nominal(argNominal), _low(lowVec), _high(highVec), _interpCode(code), _interpBoundary(1.)
127{
128 _logInit = false ;
129
130 for (auto param : paramList) {
131 if (!dynamic_cast<RooAbsReal*>(param)) {
132 coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") ERROR: paramficient " << param->GetName()
133 << " is not of type RooAbsReal" << endl ;
134 // use R__ASSERT which remains also in release mode
135 R__ASSERT(0) ;
136 }
137 _paramList.add(*param) ;
138 }
139
140 if (int(_low.size() ) != _paramList.getSize() || _low.size() != _high.size() || _low.size() != _interpCode.size()) {
141 coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") invalid input vectors " << endl;
142 R__ASSERT(int(_low.size() ) == _paramList.getSize());
143 R__ASSERT(_low.size() == _high.size());
144 R__ASSERT(_low.size() == _interpCode.size());
145 }
146
148}
149
150////////////////////////////////////////////////////////////////////////////////
151/// Constructor of flat polynomial function
152
153FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title) :
154 RooAbsReal(name, title),
155 _paramList("paramList","List of coefficients",this),
156 _nominal(0), _interpBoundary(1.)
157{
158 _logInit = false ;
160}
161
162////////////////////////////////////////////////////////////////////////////////
163
165 RooAbsReal(other, name),
166 _paramList("paramList",this,other._paramList),
167 _nominal(other._nominal), _low(other._low), _high(other._high), _interpCode(other._interpCode), _interpBoundary(other._interpBoundary)
168
169{
170 // Copy constructor
171 _logInit = false ;
173
174}
175
176
177////////////////////////////////////////////////////////////////////////////////
178/// Destructor
179
181{
183}
184
185
186////////////////////////////////////////////////////////////////////////////////
187
189 int index = _paramList.index(&param);
190 if(index<0){
191 coutE(InputArguments) << "FlexibleInterpVar::setInterpCode ERROR: " << param.GetName()
192 << " is not in list" << endl ;
193 } else {
194 coutW(InputArguments) << "FlexibleInterpVar::setInterpCode : " << param.GetName()
195 << " is now " << code << endl ;
196 _interpCode.at(index) = code;
197 }
198 // GHL: Adding suggestion by Swagato:
199 _logInit = false ;
201}
202
203////////////////////////////////////////////////////////////////////////////////
204
206 for(unsigned int i=0; i<_interpCode.size(); ++i){
207 _interpCode.at(i) = code;
208 }
209 // GHL: Adding suggestion by Swagato:
210 _logInit = false ;
212
213}
214
215////////////////////////////////////////////////////////////////////////////////
216
217void FlexibleInterpVar::setNominal(double newNominal){
218 coutW(InputArguments) << "FlexibleInterpVar::setNominal : nominal is now " << newNominal << endl ;
219 _nominal = newNominal;
220
221 _logInit = false ;
222
224}
225
226////////////////////////////////////////////////////////////////////////////////
227
228void FlexibleInterpVar::setLow(RooAbsReal& param, double newLow){
229 int index = _paramList.index(&param);
230 if(index<0){
231 coutE(InputArguments) << "FlexibleInterpVar::setLow ERROR: " << param.GetName()
232 << " is not in list" << endl ;
233 } else {
234 coutW(InputArguments) << "FlexibleInterpVar::setLow : " << param.GetName()
235 << " is now " << newLow << endl ;
236 _low.at(index) = newLow;
237 }
238
239 _logInit = false ;
240
242}
243
244////////////////////////////////////////////////////////////////////////////////
245
246void FlexibleInterpVar::setHigh(RooAbsReal& param, double newHigh){
247 int index = _paramList.index(&param);
248 if(index<0){
249 coutE(InputArguments) << "FlexibleInterpVar::setHigh ERROR: " << param.GetName()
250 << " is not in list" << endl ;
251 } else {
252 coutW(InputArguments) << "FlexibleInterpVar::setHigh : " << param.GetName()
253 << " is now " << newHigh << endl ;
254 _high.at(index) = newHigh;
255 }
256
257 _logInit = false ;
259}
260
261////////////////////////////////////////////////////////////////////////////////
262
264 for(unsigned int i=0; i<_interpCode.size(); ++i){
265 coutI(InputArguments) <<"interp code for " << _paramList.at(i)->GetName() << " = " << _interpCode.at(i) <<endl;
266 // GHL: Adding suggestion by Swagato:
267 if( _low.at(i) <= 0.001 ) coutE(InputArguments) << GetName() << ", " << _paramList.at(i)->GetName() << ": low value = " << _low.at(i) << endl;
268 if( _high.at(i) <= 0.001 ) coutE(InputArguments) << GetName() << ", " << _paramList.at(i)->GetName() << ": high value = " << _high.at(i) << endl;
269 }
270
271}
272
273////////////////////////////////////////////////////////////////////////////////
274
275double FlexibleInterpVar::PolyInterpValue(int i, double x) const {
276 // code for polynomial interpolation used when interpCode=4
277
278 double boundary = _interpBoundary;
279
280 double x0 = boundary;
281
282
283 // cache the polynomial coefficient values
284 // which do not depend on x but on the boundaries values
285 if (!_logInit) {
286
287 _logInit=true ;
288
289 unsigned int n = _low.size();
290 assert(n == _high.size() );
291
292 _polCoeff.resize(n*6) ;
293
294 for (unsigned int j = 0; j < n ; j++) {
295
296 // location of the 6 coefficient for the j-th variable
297 double * coeff = &_polCoeff[j * 6];
298
299 // GHL: Swagato's suggestions
300 double pow_up = std::pow(_high[j]/_nominal, x0);
301 double pow_down = std::pow(_low[j]/_nominal, x0);
302 double logHi = std::log(_high[j]) ;
303 double logLo = std::log(_low[j] );
304 double pow_up_log = _high[j] <= 0.0 ? 0.0 : pow_up * logHi;
305 double pow_down_log = _low[j] <= 0.0 ? 0.0 : -pow_down * logLo;
306 double pow_up_log2 = _high[j] <= 0.0 ? 0.0 : pow_up_log * logHi;
307 double pow_down_log2= _low[j] <= 0.0 ? 0.0 : -pow_down_log* logLo;
308
309 double S0 = (pow_up+pow_down)/2;
310 double A0 = (pow_up-pow_down)/2;
311 double S1 = (pow_up_log+pow_down_log)/2;
312 double A1 = (pow_up_log-pow_down_log)/2;
313 double S2 = (pow_up_log2+pow_down_log2)/2;
314 double A2 = (pow_up_log2-pow_down_log2)/2;
315
316 //fcns+der+2nd_der are eq at bd
317
318 // cache coefficient of the polynomial
319 coeff[0] = 1./(8*x0) *( 15*A0 - 7*x0*S1 + x0*x0*A2);
320 coeff[1] = 1./(8*x0*x0) *(-24 + 24*S0 - 9*x0*A1 + x0*x0*S2);
321 coeff[2] = 1./(4*pow(x0, 3))*( - 5*A0 + 5*x0*S1 - x0*x0*A2);
322 coeff[3] = 1./(4*pow(x0, 4))*( 12 - 12*S0 + 7*x0*A1 - x0*x0*S2);
323 coeff[4] = 1./(8*pow(x0, 5))*( + 3*A0 - 3*x0*S1 + x0*x0*A2);
324 coeff[5] = 1./(8*pow(x0, 6))*( -8 + 8*S0 - 5*x0*A1 + x0*x0*S2);
325
326 }
327
328 }
329
330 // GHL: Swagato's suggestions
331 // if( _low[i] == 0 ) _low[i] = 0.0001;
332 // if( _high[i] == 0 ) _high[i] = 0.0001;
333
334 // get pointer to location of coefficients in the vector
335
336 assert(int(_polCoeff.size()) > i );
337 const double * coefficients = &_polCoeff.front() + 6*i;
338
339 double a = coefficients[0];
340 double b = coefficients[1];
341 double c = coefficients[2];
342 double d = coefficients[3];
343 double e = coefficients[4];
344 double f = coefficients[5];
345
346
347 // evaluate the 6-th degree polynomial using Horner's method
348 double value = 1. + x * (a + x * ( b + x * ( c + x * ( d + x * ( e + x * f ) ) ) ) );
349 return value;
350}
351
352////////////////////////////////////////////////////////////////////////////////
353/// Const getters
354
356double FlexibleInterpVar::nominal() const { return _nominal; }
357const std::vector<double>& FlexibleInterpVar::low() const { return _low; }
358const std::vector<double>& FlexibleInterpVar::high() const { return _high; }
359
360////////////////////////////////////////////////////////////////////////////////
361/// Calculate and return value of polynomial
362
364{
365 double total(_nominal) ;
366 int i=0;
367
368 for (auto arg : _paramList) {
369 auto param = static_cast<const RooAbsReal*>(arg);
370
371 Int_t icode = _interpCode[i] ;
372
373 switch(icode) {
374
375 case 0: {
376 // piece-wise linear
377 if(param->getVal()>0)
378 total += param->getVal()*(_high[i] - _nominal );
379 else
380 total += param->getVal()*(_nominal - _low[i]);
381 break ;
382 }
383 case 1: {
384 // pice-wise log
385 if(param->getVal()>=0)
386 total *= pow(_high[i]/_nominal, +param->getVal());
387 else
388 total *= pow(_low[i]/_nominal, -param->getVal());
389 break ;
390 }
391 case 2: {
392 // parabolic with linear
393 double a = 0.5*(_high[i]+_low[i])-_nominal;
394 double b = 0.5*(_high[i]-_low[i]);
395 double c = 0;
396 if(param->getVal()>1 ){
397 total += (2*a+b)*(param->getVal()-1)+_high[i]-_nominal;
398 } else if(param->getVal()<-1 ) {
399 total += -1*(2*a-b)*(param->getVal()+1)+_low[i]-_nominal;
400 } else {
401 total += a*pow(param->getVal(),2) + b*param->getVal()+c;
402 }
403 break ;
404 }
405 case 3: {
406 //parabolic version of log-normal
407 double a = 0.5*(_high[i]+_low[i])-_nominal;
408 double b = 0.5*(_high[i]-_low[i]);
409 double c = 0;
410 if(param->getVal()>1 ){
411 total += (2*a+b)*(param->getVal()-1)+_high[i]-_nominal;
412 } else if(param->getVal()<-1 ) {
413 total += -1*(2*a-b)*(param->getVal()+1)+_low[i]-_nominal;
414 } else {
415 total += a*pow(param->getVal(),2) + b*param->getVal()+c;
416 }
417 break ;
418 }
419
420 case 4: {
421 double boundary = _interpBoundary;
422 double x = param->getVal();
423 //std::cout << icode << " param " << param->GetName() << " " << param->getVal() << " boundary " << boundary << std::endl;
424
425 if(x >= boundary)
426 {
427 total *= std::pow(_high[i]/_nominal, +param->getVal());
428 }
429 else if (x <= -boundary)
430 {
431 total *= std::pow(_low[i]/_nominal, -param->getVal());
432 }
433 else if (x != 0)
434 {
435 total *= PolyInterpValue(i, x);
436 }
437 break ;
438 }
439 default: {
440 coutE(InputArguments) << "FlexibleInterpVar::evaluate ERROR: " << param->GetName()
441 << " with unknown interpolation code" << endl ;
442 }
443 }
444 ++i;
445 }
446
447 if(total<=0) {
449 }
450
451 return total;
452}
453
454void FlexibleInterpVar::printMultiline(ostream& os, Int_t contents,
455 bool verbose, TString indent) const
456{
458 os << indent << "--- FlexibleInterpVar ---" << endl;
460}
461
463{
464 for (int i=0;i<(int)_low.size();i++) {
465 auto& param = static_cast<RooAbsReal&>(_paramList[i]);
466 os << setw(36) << param.GetName()<<": "<<setw(7) << _low[i]<<" "<<setw(7) << _high[i]
467 <<endl;
468 }
469}
470
471
#define d(i)
Definition: RSha256.hxx:102
#define f(i)
Definition: RSha256.hxx:104
#define S0(x)
Definition: RSha256.hxx:88
#define S1(x)
Definition: RSha256.hxx:89
#define c(i)
Definition: RSha256.hxx:101
#define e(i)
Definition: RSha256.hxx:103
#define coutI(a)
Definition: RooMsgService.h:34
#define coutW(a)
Definition: RooMsgService.h:36
#define coutE(a)
Definition: RooMsgService.h:37
#define TRACE_DESTROY
Definition: RooTrace.h:24
#define TRACE_CREATE
Definition: RooTrace.h:23
#define ClassImp(name)
Definition: Rtypes.h:375
static void indent(ostringstream &buf, int indent_level)
#define R__ASSERT(e)
Definition: TError.h:118
static unsigned int total
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 b
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
char name[80]
Definition: TGX11.cxx:110
void setValueDirty()
Mark the element dirty. This forces a re-evaluation when a value is requested.
Definition: RooAbsArg.h:490
Int_t getSize() const
Return the number of elements in the collection.
Int_t index(const RooAbsArg *arg) const
Returns index of given arg, or -1 if arg is not in the collection.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:62
void printMultiline(std::ostream &os, Int_t contents, bool verbose=false, TString indent="") const override
Structure printing.
Definition: RooAbsReal.cxx:464
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:22
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:110
bool add(const RooAbsArg &var, bool valueServer, bool shapeServer, bool silent)
Overloaded RooCollection_t::add() method insert object into set and registers object as server to own...
std::vector< double > _polCoeff
! cached polynomial coefficients
bool _logInit
! flag used for caching polynomial coefficients
void printMultiline(std::ostream &os, Int_t contents, bool verbose=false, TString indent="") const override
Interface for detailed printing of object.
void setInterpCode(RooAbsReal &param, int code)
const std::vector< double > & high() const
void setLow(RooAbsReal &param, double newLow)
const std::vector< double > & low() const
void setHigh(RooAbsReal &param, double newHigh)
double evaluate() const override
Calculate and return value of polynomial.
double PolyInterpValue(int i, double x) const
virtual void printFlexibleInterpVars(std::ostream &os) const
const RooListProxy & variables() const
Const getters.
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
Basic string class.
Definition: TString.h:136
RVec< PromoteTypes< T0, T1 > > pow(const T0 &x, const RVec< T1 > &v)
Definition: RVec.hxx:1783
RVec< PromoteType< T > > log(const RVec< T > &v)
Definition: RVec.hxx:1778
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
@ HistFactory
Definition: RooGlobalFunc.h:63
@ InputArguments
Definition: RooGlobalFunc.h:62
Namespace for the RooStats classes.
Definition: Asimov.h:19
static T Min()
Returns maximum representation for type T.
Definition: TMath.h:923
TArc a
Definition: textangle.C:12