Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TFormula.cxx
Go to the documentation of this file.
1// @(#)root/hist:$Id$
2// Author: Maciej Zimnoch 30/09/2013
3
4/*************************************************************************
5 * Copyright (C) 1995-2013, 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#include "TROOT.h"
13#include "TBuffer.h"
14#include "TMethod.h"
15#include "TMath.h"
16#include "TF1.h"
17#include "TMethodCall.h"
18#include <TBenchmark.h>
19#include "TError.h"
20#include "TInterpreter.h"
21#include "TInterpreterValue.h"
22#include "TFormula.h"
23#include "TRegexp.h"
24#include <array>
25#include <cassert>
26#include <iostream>
27#include <unordered_map>
28#include <functional>
29#include <set>
30
31using namespace std;
32
33#ifdef WIN32
34#pragma optimize("",off)
35#endif
36#include "v5/TFormula.h"
37
39
40/** \class TFormula TFormula.h "inc/TFormula.h"
41 \ingroup Hist
42 The Formula class
43
44 This is a new version of the TFormula class based on Cling.
45 This class is not 100% backward compatible with the old TFormula class, which is still available in ROOT as
46 `ROOT::v5::TFormula`. Some of the TFormula member functions available in version 5, such as
47 `Analyze` and `AnalyzeFunction` are not available in the new TFormula.
48 On the other hand formula expressions which were valid in version 5 are still valid in TFormula version 6
49
50 This class has been implemented during Google Summer of Code 2013 by Maciej Zimnoch.
51
52 ### Example of valid expressions:
53
54 - `sin(x)/x`
55 - `[0]*sin(x) + [1]*exp(-[2]*x)`
56 - `x + y**2`
57 - `x^2 + y^2`
58 - `[0]*pow([1],4)`
59 - `2*pi*sqrt(x/y)`
60 - `gaus(0)*expo(3) + ypol3(5)*x`
61 - `gausn(0)*expo(3) + ypol3(5)*x`
62 - `gaus(x, [0..2]) + expo(y, [3..4])`
63
64 In the last examples above:
65
66 - `gaus(0)` is a substitute for `[0]*exp(-0.5*((x-[1])/[2])**2)`
67 and (0) means start numbering parameters at 0
68 - `gausn(0)` is a substitute for `[0]*exp(-0.5*((x-[1])/[2])**2)/(sqrt(2*pi)*[2]))`
69 and (0) means start numbering parameters at 0
70 - `expo(3)` is a substitute for `exp([3]+[4]*x)`
71 - `pol3(5)` is a substitute for `par[5]+par[6]*x+par[7]*x**2+par[8]*x**3`
72 (`PolN` stands for Polynomial of degree N)
73 - `gaus(x, [0..2])` is a more explicit way of writing `gaus(0)`
74 - `expo(y, [3..4])` is a substitute for `exp([3]+[4]*y)`
75
76 `TMath` functions can be part of the expression, eg:
77
78 - `TMath::Landau(x)*sin(x)`
79 - `TMath::Erf(x)`
80
81 Formula may contain constants, eg:
82
83 - `sqrt2`
84 - `e`
85 - `pi`
86 - `ln10`
87 - `infinity`
88
89 and more.
90
91 Formulas may also contain other user-defined ROOT functions defined with a
92 TFormula, eg, where `f1` is defined on one x-dimension and 2 parameters:
93
94 - `f1(x, [omega], [phi])`
95 - `f1([0..1])`
96 - `f1([1], [0])`
97 - `f1(y)`
98
99 To replace only parameter names, the dimension variable can be dropped.
100 Alternatively, to change only the dimension variable, the parameters can be
101 dropped. Note that if a parameter is dropped or keeps its old name, its old
102 value will be copied to the new function. The syntax used in the examples
103 above also applies to the predefined parametrized functions like `gaus` and
104 `expo`.
105
106 Comparisons operators are also supported `(&amp;&amp;, ||, ==, &lt;=, &gt;=, !)`
107
108 Examples:
109
110 `sin(x*(x&lt;0.5 || x&gt;1))`
111
112 If the result of a comparison is TRUE, the result is 1, otherwise 0.
113
114 Already predefined names can be given. For example, if the formula
115
116 `TFormula old("old",sin(x*(x&lt;0.5 || x&gt;1)))`
117
118 one can assign a name to the formula. By default the name of the object = title = formula itself.
119
120 `TFormula new("new","x*old")`
121
122 is equivalent to:
123
124 `TFormula new("new","x*sin(x*(x&lt;0.5 || x&gt;1))")`
125
126 The class supports unlimited number of variables and parameters.
127 By default the names which can be used for the variables are `x,y,z,t` or
128 `x[0],x[1],x[2],x[3],....x[N]` for N-dimensional formulas.
129
130 This class is not anymore the base class for the function classes `TF1`, but it has now
131 a data member of TF1 which can be accessed via `TF1::GetFormula`.
132
133 ### An expanded note on variables and parameters
134
135 In a TFormula, a variable is a defined by a name `x`, `y`, `z` or `t` or an
136 index like `x[0]`, `x[1]`, `x[2]`; that is `x[N]` where N is an integer.
137
138 ```
139 TFormula("", "x[0] * x[1] + 10")
140 ```
141
142 Parameters are similar and can take any name. It is specified using brackets
143 e.g. `[expected_mass]` or `[0]`.
144
145 ```
146 TFormula("", "exp([expected_mass])-1")
147 ```
148
149 Variables and parameters can be combined in the same TFormula. Here we consider
150 a very simple case where we have an exponential decay after some time t and a
151 number of events with timestamps for which we want to evaluate this function.
152
153 ```
154 TFormula tf ("", "[0]*exp(-[1]*t)");
155 tf.SetParameter(0, 1);
156 tf.SetParameter(1, 0.5);
157
158 for (auto & event : events) {
159 tf.Eval(event.t);
160 }
161 ```
162
163 The distinction between variables and parameters arose from the TFormula's
164 application in fitting. There parameters are fitted to the data provided
165 through variables. In other applications this distinction can go away.
166
167 Parameter values can be provided dynamically using `TFormula::EvalPar`
168 instead of `TFormula::Eval`. In this way parameters can be used identically
169 to variables. See below for an example that uses only parameters to model a
170 function.
171
172 ```
173 Int_t params[2] = {1, 2}; // {vel_x, vel_y}
174 TFormula tf ("", "[vel_x]/sqrt(([vel_x + vel_y])**2)");
175
176 tf.EvalPar(nullptr, params);
177 ```
178
179 ### A note on operators
180
181 All operators of C/C++ are allowed in a TFormula with a few caveats.
182
183 The operators `|`, `&`, `%` can be used but will raise an error if used in
184 conjunction with a variable or a parameter. Variables and parameters are treated
185 as doubles internally for which these operators are not defined.
186 This means the following command will run successfully
187 ```root -l -q -e TFormula("", "x+(10%3)").Eval(0)```
188 but not
189 ```root -l -q -e TFormula("", "x%10").Eval(0)```.
190
191 The operator `^` is defined to mean exponentiation instead of the C/C++
192 interpretaion xor. `**` is added, also meaning exponentiation.
193
194 The operators `++` and `@` are added, and are shorthand for the a linear
195 function. That means the expression `x@2` will be expanded to
196 ```[n]*x + [n+1]*2``` where n is the first previously unused parameter number.
197
198 \class TFormulaFunction
199 Helper class for TFormula
200
201 \class TFormulaVariable
202 Another helper class for TFormula
203
204 \class TFormulaParamOrder
205 Functor defining the parameter order
206*/
207
208// prefix used for function name passed to Cling
209static const TString gNamePrefix = "TFormula__";
210
211// static map of function pointers and expressions
212//static std::unordered_map<std::string, TInterpreter::CallFuncIFacePtr_t::Generic_t> gClingFunctions = std::unordered_map<TString, TInterpreter::CallFuncIFacePtr_t::Generic_t>();
213static std::unordered_map<std::string, void *> gClingFunctions = std::unordered_map<std::string, void * >();
214
215static void R__v5TFormulaUpdater(Int_t nobjects, TObject **from, TObject **to)
216{
217 auto **fromv5 = (ROOT::v5::TFormula **)from;
218 auto **target = (TFormula **)to;
219
220 for (int i = 0; i < nobjects; ++i) {
221 if (fromv5[i] && target[i]) {
222 TFormula fnew(fromv5[i]->GetName(), fromv5[i]->GetExpFormula());
223 *(target[i]) = fnew;
224 target[i]->SetParameters(fromv5[i]->GetParameters());
225 }
226 }
227}
228
229using TFormulaUpdater_t = void (*)(Int_t nobjects, TObject **from, TObject **to);
231
233
234////////////////////////////////////////////////////////////////////////////////
236{
237 // operator ":" must be handled separately
238 const static std::set<char> ops {'+','^','-','/','*','<','>','|','&','!','=','?','%'};
239 return ops.end() != ops.find(c);
240}
241
242////////////////////////////////////////////////////////////////////////////////
244{
245 // Note that square brackets do not count as brackets here!!!
246 char brackets[] = { ')','(','{','}'};
247 Int_t bracketsLen = sizeof(brackets)/sizeof(char);
248 for(Int_t i = 0; i < bracketsLen; ++i)
249 if(brackets[i] == c)
250 return true;
251 return false;
252}
253
254////////////////////////////////////////////////////////////////////////////////
256{
257 return !IsBracket(c) && !IsOperator(c) && c != ',' && c != ' ';
258}
259
260////////////////////////////////////////////////////////////////////////////////
262{
263 return name == "x" || name == "z" || name == "y" || name == "t";
264}
265
266////////////////////////////////////////////////////////////////////////////////
268{
269 // check if the character at position i is part of a scientific notation
270 if ( (formula[i] == 'e' || formula[i] == 'E') && (i > 0 && i < formula.Length()-1) ) {
271 // handle cases: 2e+3 2e-3 2e3 and 2.e+3
272 if ( (isdigit(formula[i-1]) || formula[i-1] == '.') && ( isdigit(formula[i+1]) || formula[i+1] == '+' || formula[i+1] == '-' ) )
273 return true;
274 }
275 return false;
276}
277
278////////////////////////////////////////////////////////////////////////////////
280{
281 // check if the character at position i is part of a scientific notation
282 if ( (formula[i] == 'x' || formula[i] == 'X') && (i > 0 && i < formula.Length()-1) && formula[i-1] == '0') {
283 if (isdigit(formula[i+1]) )
284 return true;
285 static char hex_values[12] = { 'a','A', 'b','B','c','C','d','D','e','E','f','F'};
286 for (int jjj = 0; jjj < 12; ++jjj) {
287 if (formula[i+1] == hex_values[jjj])
288 return true;
289 }
290 }
291 // else
292 // return false;
293 // // handle cases: 2e+3 2e-3 2e3 and 2.e+3
294 // if ( (isdigit(formula[i-1]) || formula[i-1] == '.') && ( isdigit(formula[i+1]) || formula[i+1] == '+' || formula[i+1] == '-' ) )
295 // return true;
296 // }
297 return false;
298}
299////////////////////////////////////////////////////////////////////////////
300// check is given position is in a parameter name i.e. within "[ ]"
301////
302Bool_t TFormula::IsAParameterName(const TString & formula, int pos) {
303
304 Bool_t foundOpenParenthesis = false;
305 if (pos == 0 || pos == formula.Length()-1) return false;
306 for (int i = pos-1; i >=0; i--) {
307 if (formula[i] == ']' ) return false;
308 if (formula[i] == '[' ) {
309 foundOpenParenthesis = true;
310 break;
311 }
312 }
313 if (!foundOpenParenthesis ) return false;
314
315 // search after the position
316 for (int i = pos+1; i < formula.Length(); i++) {
317 if (formula[i] == ']' ) return true;
318 }
319 return false;
320}
321
322
323////////////////////////////////////////////////////////////////////////////////
325 // implement comparison used to set parameter orders in TFormula
326 // want p2 to be before p10
327
328 // Returns true if (a < b), meaning a comes before b, and false if (a >= b)
329
330 TRegexp numericPattern("p?[0-9]+");
331 Ssiz_t len; // buffer to store length of regex match
332
333 int patternStart = numericPattern.Index(a, &len);
334 bool aNumeric = (patternStart == 0 && len == a.Length());
335
336 patternStart = numericPattern.Index(b, &len);
337 bool bNumeric = (patternStart == 0 && len == b.Length());
338
339 if (aNumeric && !bNumeric)
340 return true; // assume a (numeric) is always before b (not numeric)
341 else if (!aNumeric && bNumeric)
342 return false; // b comes before a
343 else if (!aNumeric && !bNumeric)
344 return a < b;
345 else {
346 int aInt = (a[0] == 'p') ? TString(a(1, a.Length())).Atoi() : a.Atoi();
347 int bInt = (b[0] == 'p') ? TString(b(1, b.Length())).Atoi() : b.Atoi();
348 return aInt < bInt;
349 }
350
351}
352
353////////////////////////////////////////////////////////////////////////////////
354void TFormula::ReplaceAllNames(TString &formula, map<TString, TString> &substitutions)
355{
356 /// Apply the name substitutions to the formula, doing all replacements in one pass
357
358 for (int i = 0; i < formula.Length(); i++) {
359 // start of name
360 // (a little subtle, since we want to match names like "{V0}" and "[0]")
361 if (isalpha(formula[i]) || formula[i] == '{' || formula[i] == '[') {
362 int j; // index to end of name
363 for (j = i + 1;
364 j < formula.Length() && (IsFunctionNameChar(formula[j]) // square brackets are function name chars
365 || (formula[i] == '{' && formula[j] == '}'));
366 j++)
367 ;
368 TString name = (TString)formula(i, j - i);
369
370 // std::cout << "Looking for name: " << name << std::endl;
371
372 // if we find the name, do the substitution
373 if (substitutions.find(name) != substitutions.end()) {
374 formula.Replace(i, name.Length(), "(" + substitutions[name] + ")");
375 i += substitutions[name].Length() + 2 - 1; // +2 for parentheses
376 // std::cout << "made substitution: " << name << " to " << substitutions[name] << std::endl;
377 } else if (isalpha(formula[i])) {
378 // if formula[i] is alpha, can skip to end of candidate name, otherwise, we'll just
379 // move one character ahead and try again
380 i += name.Length() - 1;
381 }
382 }
383 }
384}
385
386////////////////////////////////////////////////////////////////////////////////
388{
389 fName = "";
390 fTitle = "";
391 fClingInput = "";
392 fReadyToExecute = false;
393 fClingInitialized = false;
394 fAllParametersSetted = false;
395 fNdim = 0;
396 fNpar = 0;
397 fNumber = 0;
398 fClingName = "";
399 fFormula = "";
400 fLambdaPtr = nullptr;
401}
402
403////////////////////////////////////////////////////////////////////////////////
404static bool IsReservedName(const char* name){
405 if (strlen(name)!=1) return false;
406 for (auto const & specialName : {"x","y","z","t"}){
407 if (strcmp(name,specialName)==0) return true;
408 }
409 return false;
410}
411
412////////////////////////////////////////////////////////////////////////////////
414{
415
416 // N.B. a memory leak may happen if user set bit after constructing the object,
417 // Setting of bit should be done only internally
420 gROOT->GetListOfFunctions()->Remove(this);
421 }
422
423 int nLinParts = fLinearParts.size();
424 if (nLinParts > 0) {
425 for (int i = 0; i < nLinParts; ++i) delete fLinearParts[i];
426 }
427}
428
429////////////////////////////////////////////////////////////////////////////////
430TFormula::TFormula(const char *name, const char *formula, bool addToGlobList, bool vectorize) :
431 TNamed(name,formula),
432 fClingInput(formula),fFormula(formula)
433{
434 fReadyToExecute = false;
435 fClingInitialized = false;
436 fNdim = 0;
437 fNpar = 0;
438 fNumber = 0;
439 fLambdaPtr = nullptr;
440 fVectorized = vectorize;
441#ifndef R__HAS_VECCORE
442 fVectorized = false;
443#endif
444
445 FillDefaults();
446
447
448 //fName = gNamePrefix + name; // is this needed
449
450 // do not process null formulas.
451 if (!fFormula.IsNull() ) {
453
454 bool ok = PrepareFormula(fFormula);
455 // if the formula has been correctly initialized add to the list of global functions
456 if (ok) {
457 if (addToGlobList && gROOT) {
458 TFormula *old = 0;
460 old = dynamic_cast<TFormula *>(gROOT->GetListOfFunctions()->FindObject(name));
461 if (old)
462 gROOT->GetListOfFunctions()->Remove(old);
463 if (IsReservedName(name))
464 Error("TFormula", "The name %s is reserved as a TFormula variable name.\n", name);
465 else
466 gROOT->GetListOfFunctions()->Add(this);
467 }
468 SetBit(kNotGlobal,!addToGlobList);
469 }
470 }
471}
472
473////////////////////////////////////////////////////////////////////////////////
474/// Constructor from a full compile-able C++ expression
475
476TFormula::TFormula(const char *name, const char *formula, int ndim, int npar, bool addToGlobList) :
477 TNamed(name,formula),
478 fClingInput(formula),fFormula(formula)
479{
480 fReadyToExecute = false;
481 fClingInitialized = false;
482 fNpar = 0;
483 fNumber = 0;
484 fLambdaPtr = nullptr;
485 fFuncPtr = nullptr;
486 fGradFuncPtr = nullptr;
487
488
489 fNdim = ndim;
490 for (int i = 0; i < npar; ++i) {
491 DoAddParameter(TString::Format("p%d",i), 0, false);
492 }
494 assert (fNpar == npar);
495
496 bool ret = InitLambdaExpression(formula);
497
498 if (ret) {
499
501
502 fReadyToExecute = true;
503
504 if (addToGlobList && gROOT) {
505 TFormula *old = 0;
507 old = dynamic_cast<TFormula*> ( gROOT->GetListOfFunctions()->FindObject(name) );
508 if (old)
509 gROOT->GetListOfFunctions()->Remove(old);
510 if (IsReservedName(name))
511 Error("TFormula","The name %s is reserved as a TFormula variable name.\n",name);
512 else
513 gROOT->GetListOfFunctions()->Add(this);
514 }
515 SetBit(kNotGlobal,!addToGlobList);
516 }
517 else
518 Error("TFormula","Syntax error in building the lambda expression %s", formula );
519}
520
521////////////////////////////////////////////////////////////////////////////////
523 TNamed(formula.GetName(),formula.GetTitle())
524{
525 formula.Copy(*this);
526
529 TFormula *old = (TFormula*)gROOT->GetListOfFunctions()->FindObject(formula.GetName());
530 if (old)
531 gROOT->GetListOfFunctions()->Remove(old);
532
533 if (IsReservedName(formula.GetName())) {
534 Error("TFormula","The name %s is reserved as a TFormula variable name.\n",formula.GetName());
535 } else
536 gROOT->GetListOfFunctions()->Add(this);
537 }
538
539}
540
541////////////////////////////////////////////////////////////////////////////////
542/// = operator.
543
545{
546
547 if (this != &rhs) {
548 rhs.Copy(*this);
549 }
550 return *this;
551}
552
553////////////////////////////////////////////////////////////////////////////////
555
556 std::string lambdaExpression = formula;
557
558 // check if formula exist already in the map
559 {
561
562 auto funcit = gClingFunctions.find(lambdaExpression);
563 if (funcit != gClingFunctions.end() ) {
564 fLambdaPtr = funcit->second;
565 fClingInitialized = true;
566 return true;
567 }
568 }
569
570 // to be sure the interpreter is initialized
573
574 // set the cling name using hash of the static formulae map
575 auto hasher = gClingFunctions.hash_function();
576 TString lambdaName = TString::Format("lambda__id%zu", hasher(lambdaExpression) );
577
578 //lambdaExpression = TString::Format("[&](double * x, double *){ return %s ;}",formula);
579 //TString lambdaName = TString::Format("mylambda_%s",GetName() );
580 TString lineExpr = TString::Format("std::function<double(double*,double*)> %s = %s ;",lambdaName.Data(), lambdaExpression.c_str() );
581 gInterpreter->ProcessLine(lineExpr);
582 fLambdaPtr = (void*) gInterpreter->ProcessLine(TString(lambdaName)+TString(";")); // add ; to avoid printing
583 if (fLambdaPtr != nullptr) {
585 gClingFunctions.insert ( std::make_pair ( lambdaExpression, fLambdaPtr) );
586 fClingInitialized = true;
587 return true;
588 }
589 fClingInitialized = false;
590 return false;
591}
592
593////////////////////////////////////////////////////////////////////////////////
594/// Compile the given expression with Cling
595/// backward compatibility method to be used in combination with the empty constructor
596/// if no expression is given , the current stored formula (retrieved with GetExpFormula()) or the title is used.
597/// return 0 if the formula compilation is successful
598
599Int_t TFormula::Compile(const char *expression)
600{
601 TString formula = expression;
602 if (formula.IsNull() ) {
603 formula = fFormula;
604 if (formula.IsNull() ) formula = GetTitle();
605 }
606
607 if (formula.IsNull() ) return -1;
608
609 // do not re-process if it was done before
610 if (IsValid() && formula == fFormula ) return 0;
611
612 // clear if a formula was already existing
613 if (!fFormula.IsNull() ) Clear();
614
615 fFormula = formula;
616
618 bool ret = InitLambdaExpression(fFormula);
619 return (ret) ? 0 : 1;
620 }
621
622 if (fVars.empty() ) FillDefaults();
623 // prepare the formula for Cling
624 //printf("compile: processing formula %s\n",fFormula.Data() );
626 // pass formula in CLing
627 bool ret = PrepareFormula(fFormula);
628
629 return (ret) ? 0 : 1;
630}
631
632////////////////////////////////////////////////////////////////////////////////
633void TFormula::Copy(TObject &obj) const
634{
635 TNamed::Copy(obj);
636 // need to copy also cling parameters
637 TFormula & fnew = dynamic_cast<TFormula&>(obj);
638
641
642 fnew.fFuncs = fFuncs;
643 fnew.fVars = fVars;
644 fnew.fParams = fParams;
645 fnew.fConsts = fConsts;
647 fnew.fFormula = fFormula;
648 fnew.fNdim = fNdim;
649 fnew.fNpar = fNpar;
650 fnew.fNumber = fNumber;
653 // copy Linear parts (it is a vector of TFormula pointers) needs to be copied one by one
654 // looping at all the elements
655 // delete first previous elements
656 int nLinParts = fnew.fLinearParts.size();
657 if (nLinParts > 0) {
658 for (int i = 0; i < nLinParts; ++i) delete fnew.fLinearParts[i];
659 fnew.fLinearParts.clear();
660 }
661 // old size that needs to be copied
662 nLinParts = fLinearParts.size();
663 if (nLinParts > 0) {
664 fnew.fLinearParts.reserve(nLinParts);
665 for (int i = 0; i < nLinParts; ++i) {
666 TFormula * linearNew = new TFormula();
667 TFormula * linearOld = (TFormula*) fLinearParts[i];
668 if (linearOld) {
669 linearOld->Copy(*linearNew);
670 fnew.fLinearParts.push_back(linearNew);
671 }
672 else
673 Warning("Copy","Function %s - expr %s has a dummy linear part %d",GetName(),GetExpFormula().Data(),i);
674 }
675 }
676
681 fnew.fClingName = fClingName;
684
685 // case of function based on a C++ expression (lambda's) which is ready to be compiled
687
688 bool ret = fnew.InitLambdaExpression(fnew.fFormula);
689 if (ret) {
691 fnew.fReadyToExecute = true;
692 }
693 else {
694 Error("TFormula","Syntax error in building the lambda expression %s", fFormula.Data() );
695 fnew.fReadyToExecute = false;
696 }
697 }
698
699 // use copy-constructor of TMethodCall
700 // if c++-14 could use std::make_unique
701 TMethodCall *m = (fMethod) ? new TMethodCall(*fMethod) : nullptr;
702 fnew.fMethod.reset(m);
703 TMethodCall *gm = (fGradMethod) ? new TMethodCall(*fGradMethod) : nullptr;
704 fnew.fGradMethod.reset(gm);
705
706 fnew.fFuncPtr = fFuncPtr;
709
710}
711
712////////////////////////////////////////////////////////////////////////////////
713/// Clear the formula setting expression to empty and reset the variables and
714/// parameters containers.
715
717{
718 fNdim = 0;
719 fNpar = 0;
720 fNumber = 0;
721 fFormula = "";
722 fClingName = "";
723
724 fMethod.reset();
725 fGradMethod.reset();
726
727 fClingVariables.clear();
728 fClingParameters.clear();
729 fReadyToExecute = false;
730 fClingInitialized = false;
731 fAllParametersSetted = false;
732 fFuncs.clear();
733 fVars.clear();
734 fParams.clear();
735 fConsts.clear();
736 fFunctionsShortcuts.clear();
737
738 // delete linear parts
739 int nLinParts = fLinearParts.size();
740 if (nLinParts > 0) {
741 for (int i = 0; i < nLinParts; ++i) delete fLinearParts[i];
742 }
743 fLinearParts.clear();
744
745}
746
747// Returns nullptr on failure.
748static std::unique_ptr<TMethodCall>
749prepareMethod(bool HasParameters, bool HasVariables, const char* FuncName,
750 bool IsVectorized, bool IsGradient = false) {
751 std::unique_ptr<TMethodCall> Method = std::unique_ptr<TMethodCall>(new TMethodCall());
752
753 TString prototypeArguments = "";
754 if (HasVariables || HasParameters) {
755 if (IsVectorized)
756 prototypeArguments.Append("ROOT::Double_v*");
757 else
758 prototypeArguments.Append("Double_t*");
759 }
760 auto AddDoublePtrParam = [&prototypeArguments] () {
761 prototypeArguments.Append(",");
762 prototypeArguments.Append("Double_t*");
763 };
764 if (HasParameters)
765 AddDoublePtrParam();
766
767 // We need an extra Double_t* for the gradient return result.
768 if (IsGradient)
769 AddDoublePtrParam();
770
771 // Initialize the method call using real function name (cling name) defined
772 // by ProcessFormula
773 Method->InitWithPrototype(FuncName, prototypeArguments);
774 if (!Method->IsValid()) {
775 Error("prepareMethod",
776 "Can't compile function %s prototype with arguments %s", FuncName,
777 prototypeArguments.Data());
778 return nullptr;
779 }
780
781 return Method;
782}
783
785 if (!method) return nullptr;
786 CallFunc_t *callfunc = method->GetCallFunc();
787
788 if (!gCling->CallFunc_IsValid(callfunc)) {
789 Error("prepareFuncPtr", "Callfunc retuned from Cling is not valid");
790 return nullptr;
791 }
792
794 = gCling->CallFunc_IFacePtr(callfunc).fGeneric;
795 if (!Result) {
796 Error("prepareFuncPtr", "Compiled function pointer is null");
797 return nullptr;
798 }
799 return Result;
800}
801
802////////////////////////////////////////////////////////////////////////////////
803/// Sets TMethodCall to function inside Cling environment.
804/// TFormula uses it to execute function.
805/// After call, TFormula should be ready to evaluate formula.
806/// Returns false on failure.
807
809{
810 if (!fMethod) {
811 Bool_t hasParameters = (fNpar > 0);
812 Bool_t hasVariables = (fNdim > 0);
813 fMethod = prepareMethod(hasParameters, hasVariables, fClingName,fVectorized);
814 if (!fMethod) return false;
816 }
817 return fFuncPtr;
818}
819
820////////////////////////////////////////////////////////////////////////////////
821/// Inputs formula, transfered to C++ code into Cling
822
824{
825
827 // make sure the interpreter is initialized
830
831 // Trigger autoloading / autoparsing (ROOT-9840):
832 TString triggerAutoparsing = "namespace ROOT_TFormula_triggerAutoParse {\n"; triggerAutoparsing += fClingInput + "\n}";
833 gCling->ProcessLine(triggerAutoparsing);
834
835 // add pragma for optimization of the formula
836 fClingInput = TString("#pragma cling optimize(2)\n") + fClingInput;
837
838 // Now that all libraries and headers are loaded, Declare() a performant version
839 // of the same code:
842 if (!fClingInitialized) Error("InputFormulaIntoCling","Error compiling formula expression in Cling");
843 }
844}
845
846////////////////////////////////////////////////////////////////////////////////
847/// Fill structures with default variables, constants and function shortcuts
848
850{
851 const TString defvars[] = { "x","y","z","t"};
852 const pair<TString, Double_t> defconsts[] = {{"pi", TMath::Pi()},
853 {"sqrt2", TMath::Sqrt2()},
854 {"infinity", TMath::Infinity()},
855 {"e", TMath::E()},
856 {"ln10", TMath::Ln10()},
857 {"loge", TMath::LogE()},
858 {"c", TMath::C()},
859 {"g", TMath::G()},
860 {"h", TMath::H()},
861 {"k", TMath::K()},
862 {"sigma", TMath::Sigma()},
863 {"r", TMath::R()},
864 {"eg", TMath::EulerGamma()},
865 {"true", 1},
866 {"false", 0}};
867 // const pair<TString,Double_t> defconsts[] = { {"pi",TMath::Pi()}, {"sqrt2",TMath::Sqrt2()},
868 // {"infinity",TMath::Infinity()}, {"ln10",TMath::Ln10()},
869 // {"loge",TMath::LogE()}, {"true",1},{"false",0} };
870 const pair<TString,TString> funShortcuts[] =
871 { {"sin","TMath::Sin" },
872 {"cos","TMath::Cos" }, {"exp","TMath::Exp"}, {"log","TMath::Log"}, {"log10","TMath::Log10"},
873 {"tan","TMath::Tan"}, {"sinh","TMath::SinH"}, {"cosh","TMath::CosH"},
874 {"tanh","TMath::TanH"}, {"asin","TMath::ASin"}, {"acos","TMath::ACos"},
875 {"atan","TMath::ATan"}, {"atan2","TMath::ATan2"}, {"sqrt","TMath::Sqrt"},
876 {"ceil","TMath::Ceil"}, {"floor","TMath::Floor"}, {"pow","TMath::Power"},
877 {"binomial","TMath::Binomial"},{"abs","TMath::Abs"},
878 {"min","TMath::Min"},{"max","TMath::Max"},{"sign","TMath::Sign" },
879 {"sq","TMath::Sq"}
880 };
881
882 std::vector<TString> defvars2(10);
883 for (int i = 0; i < 9; ++i)
884 defvars2[i] = TString::Format("x[%d]",i);
885
886 for (const auto &var : defvars) {
887 int pos = fVars.size();
888 fVars[var] = TFormulaVariable(var, 0, pos);
889 fClingVariables.push_back(0);
890 }
891 // add also the variables defined like x[0],x[1],x[2],...
892 // support up to x[9] - if needed extend that to higher value
893 // const int maxdim = 10;
894 // for (int i = 0; i < maxdim; ++i) {
895 // TString xvar = TString::Format("x[%d]",i);
896 // fVars[xvar] = TFormulaVariable(xvar,0,i);
897 // fClingVariables.push_back(0);
898 // }
899
900 for (auto con : defconsts) {
901 fConsts[con.first] = con.second;
902 }
903 if (fVectorized) {
905 } else {
906 for (auto fun : funShortcuts) {
907 fFunctionsShortcuts[fun.first] = fun.second;
908 }
909 }
910}
911
912////////////////////////////////////////////////////////////////////////////////
913/// Fill the shortcuts for vectorized functions
914/// We will replace for example sin with vecCore::Mat::Sin
915///
916
918#ifdef R__HAS_VECCORE
919 const pair<TString,TString> vecFunShortcuts[] =
920 { {"sin","vecCore::math::Sin" },
921 {"cos","vecCore::math::Cos" }, {"exp","vecCore::math::Exp"}, {"log","vecCore::math::Log"}, {"log10","vecCore::math::Log10"},
922 {"tan","vecCore::math::Tan"},
923 //{"sinh","vecCore::math::Sinh"}, {"cosh","vecCore::math::Cosh"},{"tanh","vecCore::math::Tanh"},
924 {"asin","vecCore::math::ASin"},
925 {"acos","TMath::Pi()/2-vecCore::math::ASin"},
926 {"atan","vecCore::math::ATan"},
927 {"atan2","vecCore::math::ATan2"}, {"sqrt","vecCore::math::Sqrt"},
928 {"ceil","vecCore::math::Ceil"}, {"floor","vecCore::math::Floor"}, {"pow","vecCore::math::Pow"},
929 {"cbrt","vecCore::math::Cbrt"},{"abs","vecCore::math::Abs"},
930 {"min","vecCore::math::Min"},{"max","vecCore::math::Max"},{"sign","vecCore::math::Sign" }
931 //{"sq","TMath::Sq"}, {"binomial","TMath::Binomial"} // this last two functions will not work in vectorized mode
932 };
933 // replace in the data member maps fFunctionsShortcuts
934 for (auto fun : vecFunShortcuts) {
935 fFunctionsShortcuts[fun.first] = fun.second;
936 }
937#endif
938 // do nothing in case Veccore is not enabled
939}
940
941
942////////////////////////////////////////////////////////////////////////////////
943/// Handling polN
944/// If before 'pol' exist any name, this name will be treated as variable used in polynomial
945/// eg.
946/// varpol2(5) will be replaced with: [5] + [6]*var + [7]*var^2
947/// Empty name is treated like variable x.
948
950{
951 Int_t polPos = formula.Index("pol");
952 while (polPos != kNPOS && !IsAParameterName(formula, polPos)) {
953
954 Bool_t defaultVariable = false;
955 TString variable;
956 Int_t openingBracketPos = formula.Index('(', polPos);
957 Bool_t defaultCounter = openingBracketPos == kNPOS;
958 Bool_t defaultDegree = true;
959 Int_t degree, counter;
960 TString sdegree;
961 if (!defaultCounter) {
962 // verify first of opening parenthesis belongs to pol expression
963 // character between 'pol' and '(' must all be digits
964 sdegree = formula(polPos + 3, openingBracketPos - polPos - 3);
965 if (!sdegree.IsDigit())
966 defaultCounter = true;
967 }
968 if (!defaultCounter) {
969 degree = sdegree.Atoi();
970 counter = TString(formula(openingBracketPos + 1, formula.Index(')', polPos) - openingBracketPos)).Atoi();
971 } else {
972 Int_t temp = polPos + 3;
973 while (temp < formula.Length() && isdigit(formula[temp])) {
974 defaultDegree = false;
975 temp++;
976 }
977 degree = TString(formula(polPos + 3, temp - polPos - 3)).Atoi();
978 counter = 0;
979 }
980
981 TString replacement = TString::Format("[%d]", counter);
982 if (polPos - 1 < 0 || !IsFunctionNameChar(formula[polPos - 1]) || formula[polPos - 1] == ':') {
983 variable = "x";
984 defaultVariable = true;
985 } else {
986 Int_t tmp = polPos - 1;
987 while (tmp >= 0 && IsFunctionNameChar(formula[tmp]) && formula[tmp] != ':') {
988 tmp--;
989 }
990 variable = formula(tmp + 1, polPos - (tmp + 1));
991 }
992 Int_t param = counter + 1;
993 Int_t tmp = 1;
994 while (tmp <= degree) {
995 if (tmp > 1)
996 replacement.Append(TString::Format("+[%d]*%s^%d", param, variable.Data(), tmp));
997 else
998 replacement.Append(TString::Format("+[%d]*%s", param, variable.Data()));
999 param++;
1000 tmp++;
1001 }
1002 // add parenthesis before and after
1003 if (degree > 0) {
1004 replacement.Insert(0, '(');
1005 replacement.Append(')');
1006 }
1007 TString pattern;
1008 if (defaultCounter && !defaultDegree) {
1009 pattern = TString::Format("%spol%d", (defaultVariable ? "" : variable.Data()), degree);
1010 } else if (defaultCounter && defaultDegree) {
1011 pattern = TString::Format("%spol", (defaultVariable ? "" : variable.Data()));
1012 } else {
1013 pattern = TString::Format("%spol%d(%d)", (defaultVariable ? "" : variable.Data()), degree, counter);
1014 }
1015
1016 if (!formula.Contains(pattern)) {
1017 Error("HandlePolN", "Error handling polynomial function - expression is %s - trying to replace %s with %s ",
1018 formula.Data(), pattern.Data(), replacement.Data());
1019 break;
1020 }
1021 if (formula == pattern) {
1022 // case of single polynomial
1023 SetBit(kLinear, 1);
1024 fNumber = 300 + degree;
1025 }
1026 formula.ReplaceAll(pattern, replacement);
1027 polPos = formula.Index("pol");
1028 }
1029}
1030
1031////////////////////////////////////////////////////////////////////////////////
1032/// Handling parametrized functions
1033/// Function can be normalized, and have different variable then x.
1034/// Variables should be placed in brackets after function name.
1035/// No brackets are treated like [x].
1036/// Normalized function has char 'n' after name, eg.
1037/// gausn[var](0) will be replaced with [0]*exp(-0.5*((var-[1])/[2])^2)/(sqrt(2*pi)*[2])
1038///
1039/// Adding function is easy, just follow these rules, and add to
1040/// `TFormula::FillParametrizedFunctions` defined further below:
1041///
1042/// - Key for function map is pair of name and dimension of function
1043/// - value of key is a pair function body and normalized function body
1044/// - {Vn} is a place where variable appear, n represents n-th variable from variable list.
1045/// Count starts from 0.
1046/// - [num] stands for parameter number.
1047/// If user pass to function argument 5, num will stand for (5 + num) parameter.
1048///
1049
1051{
1052 // define all parametrized functions
1053 map< pair<TString,Int_t> ,pair<TString,TString> > functions;
1054 FillParametrizedFunctions(functions);
1055
1056 map<TString,Int_t> functionsNumbers;
1057 functionsNumbers["gaus"] = 100;
1058 functionsNumbers["bigaus"] = 102;
1059 functionsNumbers["landau"] = 400;
1060 functionsNumbers["expo"] = 200;
1061 functionsNumbers["crystalball"] = 500;
1062
1063 // replace old names xygaus -> gaus[x,y]
1064 formula.ReplaceAll("xyzgaus","gaus[x,y,z]");
1065 formula.ReplaceAll("xygaus","gaus[x,y]");
1066 formula.ReplaceAll("xgaus","gaus[x]");
1067 formula.ReplaceAll("ygaus","gaus[y]");
1068 formula.ReplaceAll("zgaus","gaus[z]");
1069 formula.ReplaceAll("xexpo","expo[x]");
1070 formula.ReplaceAll("yexpo","expo[y]");
1071 formula.ReplaceAll("zexpo","expo[z]");
1072 formula.ReplaceAll("xylandau","landau[x,y]");
1073 formula.ReplaceAll("xyexpo","expo[x,y]");
1074 // at the moment pre-defined functions have no more than 3 dimensions
1075 const char * defaultVariableNames[] = { "x","y","z"};
1076
1077 for (map<pair<TString, Int_t>, pair<TString, TString>>::iterator it = functions.begin(); it != functions.end();
1078 ++it) {
1079
1080 TString funName = it->first.first;
1081 Int_t funDim = it->first.second;
1082 Int_t funPos = formula.Index(funName);
1083
1084 // std::cout << formula << " ---- " << funName << " " << funPos << std::endl;
1085 while (funPos != kNPOS && !IsAParameterName(formula, funPos)) {
1086
1087 // should also check that function is not something else (e.g. exponential - parse the expo)
1088 Int_t lastFunPos = funPos + funName.Length();
1089
1090 // check that first and last character is not a special character
1091 Int_t iposBefore = funPos - 1;
1092 // std::cout << "looping on funpos is " << funPos << " formula is " << formula << " function " << funName <<
1093 // std::endl;
1094 if (iposBefore >= 0) {
1095 assert(iposBefore < formula.Length());
1096 //if (isalpha(formula[iposBefore])) {
1097 if (IsFunctionNameChar(formula[iposBefore])) {
1098 // std::cout << "previous character for function " << funName << " is " << formula[iposBefore] << "- skip
1099 // " << std::endl;
1100 funPos = formula.Index(funName, lastFunPos);
1101 continue;
1102 }
1103 }
1104
1105 Bool_t isNormalized = false;
1106 if (lastFunPos < formula.Length()) {
1107 // check if function is normalized by looking at "n" character after function name (e.g. gausn)
1108 isNormalized = (formula[lastFunPos] == 'n');
1109 if (isNormalized)
1110 lastFunPos += 1;
1111 if (lastFunPos < formula.Length()) {
1112 char c = formula[lastFunPos];
1113 // check if also last character is not alphanumeric or is not an operator and not a parenthesis ( or [.
1114 // Parenthesis [] are used to express the variables
1115 if (isalnum(c) || (!IsOperator(c) && c != '(' && c != ')' && c != '[' && c != ']')) {
1116 // std::cout << "last character for function " << funName << " is " << c << " skip .." << std::endl;
1117 funPos = formula.Index(funName, lastFunPos);
1118 continue;
1119 }
1120 }
1121 }
1122
1123 if (isNormalized) {
1124 SetBit(kNormalized, 1);
1125 }
1126 std::vector<TString> variables;
1127 Int_t dim = 0;
1128 TString varList = "";
1129 Bool_t defaultVariables = false;
1130
1131 // check if function has specified the [...] e.g. gaus[x,y]
1132 Int_t openingBracketPos = funPos + funName.Length() + (isNormalized ? 1 : 0);
1133 Int_t closingBracketPos = kNPOS;
1134 if (openingBracketPos > formula.Length() || formula[openingBracketPos] != '[') {
1135 dim = funDim;
1136 variables.resize(dim);
1137 for (Int_t idim = 0; idim < dim; ++idim)
1138 variables[idim] = defaultVariableNames[idim];
1139 defaultVariables = true;
1140 } else {
1141 // in case of [..] found, assume they specify all the variables. Use it to get function dimension
1142 closingBracketPos = formula.Index(']', openingBracketPos);
1143 varList = formula(openingBracketPos + 1, closingBracketPos - openingBracketPos - 1);
1144 dim = varList.CountChar(',') + 1;
1145 variables.resize(dim);
1146 Int_t Nvar = 0;
1147 TString varName = "";
1148 for (Int_t i = 0; i < varList.Length(); ++i) {
1149 if (IsFunctionNameChar(varList[i])) {
1150 varName.Append(varList[i]);
1151 }
1152 if (varList[i] == ',') {
1153 variables[Nvar] = varName;
1154 varName = "";
1155 Nvar++;
1156 }
1157 }
1158 if (varName != "") // we will miss last variable
1159 {
1160 variables[Nvar] = varName;
1161 }
1162 }
1163 // check if dimension obtained from [...] is compatible with what is defined in existing pre-defined functions
1164 // std::cout << " Found dim = " << dim << " and function dimension is " << funDim << std::endl;
1165 if (dim != funDim) {
1166 pair<TString, Int_t> key = make_pair(funName, dim);
1167 if (functions.find(key) == functions.end()) {
1168 Error("PreProcessFormula", "Dimension of function %s is detected to be of dimension %d and is not "
1169 "compatible with existing pre-defined function which has dim %d",
1170 funName.Data(), dim, funDim);
1171 return;
1172 }
1173 // skip the particular function found - we might find later on the corresponding pre-defined function
1174 funPos = formula.Index(funName, lastFunPos);
1175 continue;
1176 }
1177 // look now for the (..) brackets to get the parameter counter (e.g. gaus(0) + gaus(3) )
1178 // need to start for a position
1179 Int_t openingParenthesisPos = (closingBracketPos == kNPOS) ? openingBracketPos : closingBracketPos + 1;
1180 bool defaultCounter = (openingParenthesisPos > formula.Length() || formula[openingParenthesisPos] != '(');
1181
1182 // Int_t openingParenthesisPos = formula.Index('(',funPos);
1183 // Bool_t defaultCounter = (openingParenthesisPos == kNPOS);
1184 Int_t counter;
1185 if (defaultCounter) {
1186 counter = 0;
1187 } else {
1188 // Check whether this is just a number in parentheses. If not, leave
1189 // it to `HandleFunctionArguments` to be parsed
1190
1191 TRegexp counterPattern("([0-9]+)");
1192 Ssiz_t len;
1193 if (counterPattern.Index(formula, &len, openingParenthesisPos) == -1) {
1194 funPos = formula.Index(funName, funPos + 1);
1195 continue;
1196 } else {
1197 counter =
1198 TString(formula(openingParenthesisPos + 1, formula.Index(')', funPos) - openingParenthesisPos - 1))
1199 .Atoi();
1200 }
1201 }
1202 // std::cout << "openingParenthesisPos " << openingParenthesisPos << " counter is " << counter << std::endl;
1203
1204 TString body = (isNormalized ? it->second.second : it->second.first);
1205 if (isNormalized && body == "") {
1206 Error("PreprocessFormula", "%d dimension function %s has no normalized form.", it->first.second,
1207 funName.Data());
1208 break;
1209 }
1210 for (int i = 0; i < body.Length(); ++i) {
1211 if (body[i] == '{') {
1212 // replace {Vn} with variable names
1213 i += 2; // skip '{' and 'V'
1214 Int_t num = TString(body(i, body.Index('}', i) - i)).Atoi();
1215 TString variable = variables[num];
1216 TString pattern = TString::Format("{V%d}", num);
1217 i -= 2; // restore original position
1218 body.Replace(i, pattern.Length(), variable, variable.Length());
1219 i += variable.Length() - 1; // update i to reflect change in body string
1220 } else if (body[i] == '[') {
1221 // update parameter counters in case of many functions (e.g. gaus(0)+gaus(3) )
1222 Int_t tmp = i;
1223 while (tmp < body.Length() && body[tmp] != ']') {
1224 tmp++;
1225 }
1226 Int_t num = TString(body(i + 1, tmp - 1 - i)).Atoi();
1227 num += counter;
1228 TString replacement = TString::Format("%d", num);
1229
1230 body.Replace(i + 1, tmp - 1 - i, replacement, replacement.Length());
1231 i += replacement.Length() + 1;
1232 }
1233 }
1234 TString pattern;
1235 if (defaultCounter && defaultVariables) {
1236 pattern = TString::Format("%s%s", funName.Data(), (isNormalized ? "n" : ""));
1237 }
1238 if (!defaultCounter && defaultVariables) {
1239 pattern = TString::Format("%s%s(%d)", funName.Data(), (isNormalized ? "n" : ""), counter);
1240 }
1241 if (defaultCounter && !defaultVariables) {
1242 pattern = TString::Format("%s%s[%s]", funName.Data(), (isNormalized ? "n" : ""), varList.Data());
1243 }
1244 if (!defaultCounter && !defaultVariables) {
1245 pattern =
1246 TString::Format("%s%s[%s](%d)", funName.Data(), (isNormalized ? "n" : ""), varList.Data(), counter);
1247 }
1248 TString replacement = body;
1249
1250 // set the number (only in case a function exists without anything else
1251 if (fNumber == 0 && formula.Length() <= (pattern.Length() - funPos) + 1) { // leave 1 extra
1252 fNumber = functionsNumbers[funName] + 10 * (dim - 1);
1253 }
1254
1255 // std::cout << " replace " << pattern << " with " << replacement << std::endl;
1256
1257 formula.Replace(funPos, pattern.Length(), replacement, replacement.Length());
1258
1259 funPos = formula.Index(funName);
1260 }
1261 // std::cout << " End loop of " << funName << " formula is now " << formula << std::endl;
1262 }
1263}
1264
1265////////////////////////////////////////////////////////////////////////////////
1266/// Handling parameter ranges, in the form of [1..5]
1268{
1269 TRegexp rangePattern("\\[[0-9]+\\.\\.[0-9]+\\]");
1270 Ssiz_t len;
1271 int matchIdx = 0;
1272 while ((matchIdx = rangePattern.Index(formula, &len, matchIdx)) != -1) {
1273 int startIdx = matchIdx + 1;
1274 int endIdx = formula.Index("..", startIdx) + 2; // +2 for ".."
1275 int startCnt = TString(formula(startIdx, formula.Length())).Atoi();
1276 int endCnt = TString(formula(endIdx, formula.Length())).Atoi();
1277
1278 if (endCnt <= startCnt)
1279 Error("HandleParamRanges", "End parameter (%d) <= start parameter (%d) in parameter range", endCnt, startCnt);
1280
1281 TString newString = "[";
1282 for (int cnt = startCnt; cnt < endCnt; cnt++)
1283 newString += TString::Format("%d],[", cnt);
1284 newString += TString::Format("%d]", endCnt);
1285
1286 // std::cout << "newString generated by HandleParamRanges is " << newString << std::endl;
1287 formula.Replace(matchIdx, formula.Index("]", matchIdx) + 1 - matchIdx, newString);
1288
1289 matchIdx += newString.Length();
1290 }
1291
1292 // std::cout << "final formula is now " << formula << std::endl;
1293}
1294
1295////////////////////////////////////////////////////////////////////////////////
1296/// Handling user functions (and parametrized functions)
1297/// to take variables and optionally parameters as arguments
1299{
1300 // std::cout << "calling `HandleFunctionArguments` on " << formula << std::endl;
1301
1302 // Define parametrized functions, in case we need to use them
1303 std::map<std::pair<TString, Int_t>, std::pair<TString, TString>> parFunctions;
1304 FillParametrizedFunctions(parFunctions);
1305
1306 // loop through characters
1307 for (Int_t i = 0; i < formula.Length(); ++i) {
1308 // List of things to ignore (copied from `TFormula::ExtractFunctors`)
1309
1310 // ignore things that start with square brackets
1311 if (formula[i] == '[') {
1312 while (formula[i] != ']')
1313 i++;
1314 continue;
1315 }
1316 // ignore strings
1317 if (formula[i] == '\"') {
1318 do
1319 i++;
1320 while (formula[i] != '\"');
1321 continue;
1322 }
1323 // ignore numbers (scientific notation)
1324 if (IsScientificNotation(formula, i))
1325 continue;
1326 // ignore x in hexadecimal number
1327 if (IsHexadecimal(formula, i)) {
1328 while (!IsOperator(formula[i]) && i < formula.Length())
1329 i++;
1330 continue;
1331 }
1332
1333 // investigate possible start of function name
1334 if (isalpha(formula[i]) && !IsOperator(formula[i])) {
1335 // std::cout << "character : " << i << " " << formula[i] << " is not an operator and is alpha" << std::endl;
1336
1337 int j; // index to end of name
1338 for (j = i; j < formula.Length() && IsFunctionNameChar(formula[j]); j++)
1339 ;
1340 TString name = (TString)formula(i, j - i);
1341 // std::cout << "parsed name " << name << std::endl;
1342
1343 // Count arguments (careful about parentheses depth)
1344 // Make list of indices where each argument is separated
1345 int nArguments = 1;
1346 int depth = 1;
1347 std::vector<int> argSeparators;
1348 argSeparators.push_back(j); // opening parenthesis
1349 int k; // index for end of closing parenthesis
1350 for (k = j + 1; depth >= 1 && k < formula.Length(); k++) {
1351 if (formula[k] == ',' && depth == 1) {
1352 nArguments++;
1353 argSeparators.push_back(k);
1354 } else if (formula[k] == '(')
1355 depth++;
1356 else if (formula[k] == ')')
1357 depth--;
1358 }
1359 argSeparators.push_back(k - 1); // closing parenthesis
1360
1361 // retrieve `f` (code copied from ExtractFunctors)
1362 TObject *obj = 0;
1363 {
1365 obj = gROOT->GetListOfFunctions()->FindObject(name);
1366 }
1367 TFormula *f = dynamic_cast<TFormula *>(obj);
1368 if (!f) {
1369 // maybe object is a TF1
1370 TF1 *f1 = dynamic_cast<TF1 *>(obj);
1371 if (f1)
1372 f = f1->GetFormula();
1373 }
1374 // `f` should be found by now, if it was a user-defined function.
1375 // The other possibility we need to consider is that this is a
1376 // parametrized function (else case below)
1377
1378 bool nameRecognized = (f != nullptr);
1379
1380 // Get ndim, npar, and replacementFormula of function
1381 int ndim = 0;
1382 int npar = 0;
1383 TString replacementFormula;
1384 if (f) {
1385 ndim = f->GetNdim();
1386 npar = f->GetNpar();
1387 replacementFormula = f->GetExpFormula();
1388 } else {
1389 // otherwise, try to match default parametrized functions
1390
1391 for (auto keyval : parFunctions) {
1392 // (name, ndim)
1393 pair<TString, Int_t> name_ndim = keyval.first;
1394 // (formula without normalization, formula with normalization)
1395 pair<TString, TString> formulaPair = keyval.second;
1396
1397 // match names like gaus, gausn, breitwigner
1398 if (name == name_ndim.first)
1399 replacementFormula = formulaPair.first;
1400 else if (name == name_ndim.first + "n" && formulaPair.second != "")
1401 replacementFormula = formulaPair.second;
1402 else
1403 continue;
1404
1405 // set ndim
1406 ndim = name_ndim.second;
1407
1408 // go through replacementFormula to find the number of parameters
1409 npar = 0;
1410 int idx = 0;
1411 while ((idx = replacementFormula.Index('[', idx)) != kNPOS) {
1412 npar = max(npar, 1 + TString(replacementFormula(idx + 1, replacementFormula.Length())).Atoi());
1413 idx = replacementFormula.Index(']', idx);
1414 if (idx == kNPOS)
1415 Error("HandleFunctionArguments", "Square brackets not matching in formula %s",
1416 (const char *)replacementFormula);
1417 }
1418 // npar should be set correctly now
1419
1420 // break if number of arguments is good (note: `gaus`, has two
1421 // definitions with different numbers of arguments, but it works
1422 // out so that it should be unambiguous)
1423 if (nArguments == ndim + npar || nArguments == npar || nArguments == ndim) {
1424 nameRecognized = true;
1425 break;
1426 }
1427 }
1428 }
1429 if (nameRecognized && ndim > 4)
1430 Error("HandleFunctionArguments", "Number of dimensions %d greater than 4. Cannot parse formula.", ndim);
1431
1432 // if we have "recognizedName(...", then apply substitutions
1433 if (nameRecognized && j < formula.Length() && formula[j] == '(') {
1434 // std::cout << "naive replacement formula: " << replacementFormula << std::endl;
1435 // std::cout << "formula: " << formula << std::endl;
1436
1437 // map to rename each argument in `replacementFormula`
1438 map<TString, TString> argSubstitutions;
1439
1440 const char *defaultVariableNames[] = {"x", "y", "z", "t"};
1441
1442 // check nArguments and add to argSubstitutions map as appropriate
1443 bool canReplace = false;
1444 if (nArguments == ndim + npar) {
1445 // loop through all variables and parameters, filling in argSubstitutions
1446 for (int argNr = 0; argNr < nArguments; argNr++) {
1447
1448 // Get new name (for either variable or parameter)
1449 TString newName =
1450 TString(formula(argSeparators[argNr] + 1, argSeparators[argNr + 1] - argSeparators[argNr] - 1));
1451 PreProcessFormula(newName); // so that nesting works
1452
1453 // Get old name(s)
1454 // and add to argSubstitutions map as appropriate
1455 if (argNr < ndim) { // variable
1456 TString oldName = (f) ? TString::Format("x[%d]", argNr) : TString::Format("{V%d}", argNr);
1457 argSubstitutions[oldName] = newName;
1458
1459 if (f)
1460 argSubstitutions[defaultVariableNames[argNr]] = newName;
1461
1462 } else { // parameter
1463 int parNr = argNr - ndim;
1464 TString oldName =
1465 (f) ? TString::Format("[%s]", f->GetParName(parNr)) : TString::Format("[%d]", parNr);
1466 argSubstitutions[oldName] = newName;
1467
1468 // If the name stays the same, keep the old value of the parameter
1469 if (f && oldName == newName)
1470 DoAddParameter(f->GetParName(parNr), f->GetParameter(parNr), false);
1471 }
1472 }
1473
1474 canReplace = true;
1475 } else if (nArguments == npar) {
1476 // Try to assume variables are implicit (need all arguments to be
1477 // parameters)
1478
1479 // loop to check if all arguments are parameters
1480 bool varsImplicit = true;
1481 for (int argNr = 0; argNr < nArguments && varsImplicit; argNr++) {
1482 int openIdx = argSeparators[argNr] + 1;
1483 int closeIdx = argSeparators[argNr + 1] - 1;
1484
1485 // check brackets on either end
1486 if (formula[openIdx] != '[' || formula[closeIdx] != ']' || closeIdx <= openIdx + 1)
1487 varsImplicit = false;
1488
1489 // check that the middle is a single function-name
1490 for (int idx = openIdx + 1; idx < closeIdx && varsImplicit; idx++)
1491 if (!IsFunctionNameChar(formula[idx]))
1492 varsImplicit = false;
1493
1494 if (!varsImplicit)
1495 Warning("HandleFunctionArguments",
1496 "Argument %d is not a parameter. Cannot assume variables are implicit.", argNr);
1497 }
1498
1499 // loop to replace parameter names
1500 if (varsImplicit) {
1501 // if parametrized function, still need to replace parameter names
1502 if (!f) {
1503 for (int dim = 0; dim < ndim; dim++) {
1504 argSubstitutions[TString::Format("{V%d}", dim)] = defaultVariableNames[dim];
1505 }
1506 }
1507
1508 for (int argNr = 0; argNr < nArguments; argNr++) {
1509 TString oldName =
1510 (f) ? TString::Format("[%s]", f->GetParName(argNr)) : TString::Format("[%d]", argNr);
1511 TString newName =
1512 TString(formula(argSeparators[argNr] + 1, argSeparators[argNr + 1] - argSeparators[argNr] - 1));
1513
1514 // preprocess the formula so that nesting works
1515 PreProcessFormula(newName);
1516 argSubstitutions[oldName] = newName;
1517
1518 // If the name stays the same, keep the old value of the parameter
1519 if (f && oldName == newName)
1520 DoAddParameter(f->GetParName(argNr), f->GetParameter(argNr), false);
1521 }
1522
1523 canReplace = true;
1524 }
1525 }
1526 if (!canReplace && nArguments == ndim) {
1527 // Treat parameters as implicit
1528
1529 // loop to replace variable names
1530 for (int argNr = 0; argNr < nArguments; argNr++) {
1531 TString oldName = (f) ? TString::Format("x[%d]", argNr) : TString::Format("{V%d}", argNr);
1532 TString newName =
1533 TString(formula(argSeparators[argNr] + 1, argSeparators[argNr + 1] - argSeparators[argNr] - 1));
1534
1535 // preprocess so nesting works
1536 PreProcessFormula(newName);
1537 argSubstitutions[oldName] = newName;
1538
1539 if (f) // x, y, z are not used in parametrized function definitions
1540 argSubstitutions[defaultVariableNames[argNr]] = newName;
1541 }
1542
1543 if (f) {
1544 // keep old values of the parameters
1545 for (int parNr = 0; parNr < npar; parNr++)
1546 DoAddParameter(f->GetParName(parNr), f->GetParameter(parNr), false);
1547 }
1548
1549 canReplace = true;
1550 }
1551
1552 if (canReplace)
1553 ReplaceAllNames(replacementFormula, argSubstitutions);
1554 // std::cout << "after replacement, replacementFormula is " << replacementFormula << std::endl;
1555
1556 if (canReplace) {
1557 // std::cout << "about to replace position " << i << " length " << k-i << " in formula : " << formula <<
1558 // std::endl;
1559 formula.Replace(i, k - i, replacementFormula);
1560 i += replacementFormula.Length() - 1; // skip to end of replacement
1561 // std::cout << "new formula is : " << formula << std::endl;
1562 } else {
1563 Warning("HandleFunctionArguments", "Unable to make replacement. Number of parameters doesn't work : "
1564 "%d arguments, %d dimensions, %d parameters",
1565 nArguments, ndim, npar);
1566 i = j;
1567 }
1568
1569 } else {
1570 i = j; // skip to end of candidate "name"
1571 }
1572 }
1573 }
1574
1575}
1576
1577////////////////////////////////////////////////////////////////////////////////
1578/// Handling exponentiation
1579/// Can handle multiple carets, eg.
1580/// 2^3^4 will be treated like 2^(3^4)
1581
1583{
1584 Int_t caretPos = formula.Last('^');
1585 while (caretPos != kNPOS && !IsAParameterName(formula, caretPos)) {
1586
1587 TString right, left;
1588 Int_t temp = caretPos;
1589 temp--;
1590 // get the expression in ( ) which has the operator^ applied
1591 if (formula[temp] == ')') {
1592 Int_t depth = 1;
1593 temp--;
1594 while (depth != 0 && temp > 0) {
1595 if (formula[temp] == ')')
1596 depth++;
1597 if (formula[temp] == '(')
1598 depth--;
1599 temp--;
1600 }
1601 if (depth == 0)
1602 temp++;
1603 }
1604 // this in case of someting like sin(x+2)^2
1605 do {
1606 temp--; // go down one
1607 // handle scientific notation cases (1.e-2 ^ 3 )
1608 if (temp >= 2 && IsScientificNotation(formula, temp - 1))
1609 temp -= 3;
1610 } while (temp >= 0 && !IsOperator(formula[temp]) && !IsBracket(formula[temp]));
1611
1612 assert(temp + 1 >= 0);
1613 Int_t leftPos = temp + 1;
1614 left = formula(leftPos, caretPos - leftPos);
1615 // std::cout << "left to replace is " << left << std::endl;
1616
1617 // look now at the expression after the ^ operator
1618 temp = caretPos;
1619 temp++;
1620 if (temp >= formula.Length()) {
1621 Error("HandleExponentiation", "Invalid position of operator ^");
1622 return;
1623 }
1624 if (formula[temp] == '(') {
1625 Int_t depth = 1;
1626 temp++;
1627 while (depth != 0 && temp < formula.Length()) {
1628 if (formula[temp] == ')')
1629 depth--;
1630 if (formula[temp] == '(')
1631 depth++;
1632 temp++;
1633 }
1634 temp--;
1635 } else {
1636 // handle case first character is operator - or + continue
1637 if (formula[temp] == '-' || formula[temp] == '+')
1638 temp++;
1639 // handle cases x^-2 or x^+2
1640 // need to handle also cases x^sin(x+y)
1641 Int_t depth = 0;
1642 // stop right expression if is an operator or if is a ")" from a zero depth
1643 while (temp < formula.Length() && ((depth > 0) || !IsOperator(formula[temp]))) {
1644 temp++;
1645 // handle scientific notation cases (1.e-2 ^ 3 )
1646 if (temp >= 2 && IsScientificNotation(formula, temp))
1647 temp += 2;
1648 // for internal parenthesis
1649 if (temp < formula.Length() && formula[temp] == '(')
1650 depth++;
1651 if (temp < formula.Length() && formula[temp] == ')') {
1652 if (depth > 0)
1653 depth--;
1654 else
1655 break; // case of end of a previously started expression e.g. sin(x^2)
1656 }
1657 }
1658 }
1659 right = formula(caretPos + 1, (temp - 1) - caretPos);
1660 // std::cout << "right to replace is " << right << std::endl;
1661
1662 TString pattern = TString::Format("%s^%s", left.Data(), right.Data());
1663 TString replacement = TString::Format("pow(%s,%s)", left.Data(), right.Data());
1664
1665 // std::cout << "pattern : " << pattern << std::endl;
1666 // std::cout << "replacement : " << replacement << std::endl;
1667 formula.Replace(leftPos, pattern.Length(), replacement, replacement.Length());
1668
1669 caretPos = formula.Last('^');
1670 }
1671}
1672
1673////////////////////////////////////////////////////////////////////////////////
1674/// Handle linear functions defined with the operator ++.
1675
1677{
1678 // Handle Linear functions identified with "@" operator
1679 Int_t linPos = formula.Index("@");
1680 if (linPos == kNPOS ) return; // function is not linear
1681 Int_t nofLinParts = formula.CountChar((int)'@');
1682 assert(nofLinParts > 0);
1683 fLinearParts.reserve(nofLinParts + 1);
1684 Int_t Nlinear = 0;
1685 bool first = true;
1686 while (linPos != kNPOS && !IsAParameterName(formula, linPos)) {
1687 SetBit(kLinear, 1);
1688 // analyze left part only the first time
1689 Int_t temp = 0;
1690 TString left;
1691 if (first) {
1692 temp = linPos - 1;
1693 while (temp >= 0 && formula[temp] != '@') {
1694 temp--;
1695 }
1696 left = formula(temp + 1, linPos - (temp + 1));
1697 }
1698 temp = linPos + 1;
1699 while (temp < formula.Length() && formula[temp] != '@') {
1700 temp++;
1701 }
1702 TString right = formula(linPos + 1, temp - (linPos + 1));
1703
1704 TString pattern =
1705 (first) ? TString::Format("%s@%s", left.Data(), right.Data()) : TString::Format("@%s", right.Data());
1706 TString replacement =
1707 (first) ? TString::Format("([%d]*(%s))+([%d]*(%s))", Nlinear, left.Data(), Nlinear + 1, right.Data())
1708 : TString::Format("+([%d]*(%s))", Nlinear, right.Data());
1709 Nlinear += (first) ? 2 : 1;
1710
1711 formula.ReplaceAll(pattern, replacement);
1712 if (first) {
1713 TFormula *lin1 = new TFormula("__linear1", left, false);
1714 fLinearParts.push_back(lin1);
1715 }
1716 TFormula *lin2 = new TFormula("__linear2", right, false);
1717 fLinearParts.push_back(lin2);
1718
1719 linPos = formula.Index("@");
1720 first = false;
1721 }
1722}
1723
1724////////////////////////////////////////////////////////////////////////////////
1725/// Preprocessing of formula
1726/// Replace all ** by ^, and removes spaces.
1727/// Handle also parametrized functions like polN,gaus,expo,landau
1728/// and exponentiation.
1729/// Similar functionality should be added here.
1730
1732{
1733 formula.ReplaceAll("**","^");
1734 formula.ReplaceAll("++","@"); // for linear functions
1735 formula.ReplaceAll(" ","");
1736 HandlePolN(formula);
1738 HandleParamRanges(formula);
1739 HandleFunctionArguments(formula);
1740 HandleExponentiation(formula);
1741 // "++" wil be dealt with Handle Linear
1742 HandleLinear(formula);
1743 // special case for "--" and "++"
1744 // ("++" needs to be written with whitespace that is removed before but then we re-add it again
1745 formula.ReplaceAll("--","- -");
1746 formula.ReplaceAll("++","+ +");
1747}
1748
1749////////////////////////////////////////////////////////////////////////////////
1750/// prepare the formula to be executed
1751/// normally is called with fFormula
1752
1754{
1755 fFuncs.clear();
1756 fReadyToExecute = false;
1757 ExtractFunctors(formula);
1758
1759 // update the expression with the new formula
1760 fFormula = formula;
1761 // save formula to parse variable and parameters for Cling
1762 fClingInput = formula;
1763 // replace all { and }
1764 fFormula.ReplaceAll("{","");
1765 fFormula.ReplaceAll("}","");
1766
1767 // std::cout << "functors are extracted formula is " << std::endl;
1768 // std::cout << fFormula << std::endl << std::endl;
1769
1770 fFuncs.sort();
1771 fFuncs.unique();
1772
1773 // use inputFormula for Cling
1775
1776 // for pre-defined functions (need after processing)
1777 if (fNumber != 0) SetPredefinedParamNames();
1778
1780}
1781
1782////////////////////////////////////////////////////////////////////////////////
1783/// Extracts functors from formula, and put them in fFuncs.
1784/// Simple grammar:
1785/// - <function> := name(arg1,arg2...)
1786/// - <variable> := name
1787/// - <parameter> := [number]
1788/// - <name> := String containing lower and upper letters, numbers, underscores
1789/// - <number> := Integer number
1790/// Operators are omitted.
1791
1793{
1794 // std::cout << "Commencing ExtractFunctors on " << formula << std::endl;
1795
1796 TString name = "";
1797 TString body = "";
1798 // printf("formula is : %s \n",formula.Data() );
1799 for (Int_t i = 0; i < formula.Length(); ++i) {
1800
1801 // std::cout << "loop on character : " << i << " " << formula[i] << std::endl;
1802 // case of parameters
1803 if (formula[i] == '[') {
1804 Int_t tmp = i;
1805 i++;
1806 TString param = "";
1807 while (i < formula.Length() && formula[i] != ']') {
1808 param.Append(formula[i++]);
1809 }
1810 i++;
1811 // rename parameter name XX to pXX
1812 // std::cout << "examine parameters " << param << std::endl;
1813 int paramIndex = -1;
1814 if (param.IsDigit()) {
1815 paramIndex = param.Atoi();
1816 param.Insert(0, 'p'); // needed for the replacement
1817 if (paramIndex >= fNpar || fParams.find(param) == fParams.end()) {
1818 // add all parameters up to given index found
1819 for (int idx = 0; idx <= paramIndex; ++idx) {
1820 TString pname = TString::Format("p%d", idx);
1821 if (fParams.find(pname) == fParams.end())
1822 DoAddParameter(pname, 0, false);
1823 }
1824 }
1825 } else {
1826 // handle whitespace characters in parname
1827 param.ReplaceAll("\\s", " ");
1828
1829 // only add if parameter does not already exist, because maybe
1830 // `HandleFunctionArguments` already assigned a default value to the
1831 // parameter
1832 if (fParams.find(param) == fParams.end() || GetParNumber(param) < 0 ||
1833 (unsigned)GetParNumber(param) >= fClingParameters.size()) {
1834 // std::cout << "Setting parameter " << param << " to 0" << std::endl;
1835 DoAddParameter(param, 0, false);
1836 }
1837 }
1838 TString replacement = TString::Format("{[%s]}", param.Data());
1839 formula.Replace(tmp, i - tmp, replacement, replacement.Length());
1840 fFuncs.push_back(TFormulaFunction(param));
1841 // we need to change index i after replacing since string length changes
1842 // and we need to re-calculate i position
1843 int deltai = replacement.Length() - (i-tmp);
1844 i += deltai;
1845 // printf("found parameter %s \n",param.Data() );
1846 continue;
1847 }
1848 // case of strings
1849 if (formula[i] == '\"') {
1850 // look for next instance of "\"
1851 do {
1852 i++;
1853 } while (formula[i] != '\"');
1854 }
1855 // case of e or E for numbers in exponential notaton (e.g. 2.2e-3)
1856 if (IsScientificNotation(formula, i))
1857 continue;
1858 // case of x for hexadecimal numbers
1859 if (IsHexadecimal(formula, i)) {
1860 // find position of operator
1861 // do not check cases if character is not only a to f, but accept anything
1862 while (!IsOperator(formula[i]) && i < formula.Length()) {
1863 i++;
1864 }
1865 continue;
1866 }
1867
1868 // std::cout << "investigating character : " << i << " " << formula[i] << " of formula " << formula <<
1869 // std::endl;
1870 // look for variable and function names. They start in C++ with alphanumeric characters
1871 if (isalpha(formula[i]) &&
1872 !IsOperator(formula[i])) // not really needed to check if operator (if isalpha is not an operator)
1873 {
1874 // std::cout << "character : " << i << " " << formula[i] << " is not an operator and is alpha " <<
1875 // std::endl;
1876
1877 while (i < formula.Length() && IsFunctionNameChar(formula[i])) {
1878 // need special case for separating operator ":" from scope operator "::"
1879 if (formula[i] == ':' && ((i + 1) < formula.Length())) {
1880 if (formula[i + 1] == ':') {
1881 // case of :: (scopeOperator)
1882 name.Append("::");
1883 i += 2;
1884 continue;
1885 } else
1886 break;
1887 }
1888
1889 name.Append(formula[i++]);
1890 }
1891 // printf(" build a name %s \n",name.Data() );
1892 if (formula[i] == '(') {
1893 i++;
1894 if (formula[i] == ')') {
1895 fFuncs.push_back(TFormulaFunction(name, body, 0));
1896 name = body = "";
1897 continue;
1898 }
1899 Int_t depth = 1;
1900 Int_t args = 1; // we will miss first argument
1901 while (depth != 0 && i < formula.Length()) {
1902 switch (formula[i]) {
1903 case '(': depth++; break;
1904 case ')': depth--; break;
1905 case ',':
1906 if (depth == 1)
1907 args++;
1908 break;
1909 }
1910 if (depth != 0) // we don't want last ')' inside body
1911 {
1912 body.Append(formula[i++]);
1913 }
1914 }
1915 Int_t originalBodyLen = body.Length();
1916 ExtractFunctors(body);
1917 formula.Replace(i - originalBodyLen, originalBodyLen, body, body.Length());
1918 i += body.Length() - originalBodyLen;
1919 fFuncs.push_back(TFormulaFunction(name, body, args));
1920 } else {
1921
1922 // std::cout << "check if character : " << i << " " << formula[i] << " from name " << name << " is a
1923 // function " << std::endl;
1924
1925 // check if function is provided by gROOT
1926 TObject *obj = 0;
1927 // exclude case function name is x,y,z,t
1928 if (!IsReservedName(name))
1929 {
1931 obj = gROOT->GetListOfFunctions()->FindObject(name);
1932 }
1933 TFormula *f = dynamic_cast<TFormula *>(obj);
1934 if (!f) {
1935 // maybe object is a TF1
1936 TF1 *f1 = dynamic_cast<TF1 *>(obj);
1937 if (f1)
1938 f = f1->GetFormula();
1939 }
1940 if (f) {
1941 // Replacing user formula the old way (as opposed to 'HandleFunctionArguments')
1942 // Note this is only for replacing functions that do
1943 // not specify variables and/or parameters in brackets
1944 // (the other case is done by `HandleFunctionArguments`)
1945
1946 TString replacementFormula = f->GetExpFormula();
1947
1948 // analyze expression string
1949 // std::cout << "formula to replace for " << f->GetName() << " is " << replacementFormula <<
1950 // std::endl;
1951 PreProcessFormula(replacementFormula);
1952 // we need to define different parameters if we use the unnamed default parameters ([0])
1953 // I need to replace all the terms in the functor for backward compatibility of the case
1954 // f1("[0]*x") f2("[0]*x") f1+f2 - it is weird but it is better to support
1955 // std::cout << "current number of parameter is " << fNpar << std::endl;
1956 int nparOffset = 0;
1957 // if (fParams.find("0") != fParams.end() ) {
1958 // do in any case if parameters are existing
1959 std::vector<TString> newNames;
1960 if (fNpar > 0) {
1961 nparOffset = fNpar;
1962 newNames.resize(f->GetNpar());
1963 // start from higher number to avoid overlap
1964 for (int jpar = f->GetNpar() - 1; jpar >= 0; --jpar) {
1965 // parameters name have a "p" added in front
1966 TString pj = TString(f->GetParName(jpar));
1967 if (pj[0] == 'p' && TString(pj(1, pj.Length())).IsDigit()) {
1968 TString oldName = TString::Format("[%s]", f->GetParName(jpar));
1969 TString newName = TString::Format("[p%d]", nparOffset + jpar);
1970 // std::cout << "replace - parameter " << f->GetParName(jpar) << " with " << newName <<
1971 // std::endl;
1972 replacementFormula.ReplaceAll(oldName, newName);
1973 newNames[jpar] = newName;
1974 } else
1975 newNames[jpar] = f->GetParName(jpar);
1976 }
1977 // std::cout << "after replacing params " << replacementFormula << std::endl;
1978 }
1979 ExtractFunctors(replacementFormula);
1980 // std::cout << "after re-extracting functors " << replacementFormula << std::endl;
1981
1982 // set parameter value from replacement formula
1983 for (int jpar = 0; jpar < f->GetNpar(); ++jpar) {
1984 if (nparOffset > 0) {
1985 // parameter have an offset- so take this into account
1986 assert((int)newNames.size() == f->GetNpar());
1987 SetParameter(newNames[jpar], f->GetParameter(jpar));
1988 } else
1989 // names are the same between current formula and replaced one
1990 SetParameter(f->GetParName(jpar), f->GetParameter(jpar));
1991 }
1992 // need to add parenthesis at begin and end of replacementFormula
1993 replacementFormula.Insert(0, '(');
1994 replacementFormula.Insert(replacementFormula.Length(), ')');
1995 formula.Replace(i - name.Length(), name.Length(), replacementFormula, replacementFormula.Length());
1996 // move forward the index i of the main loop
1997 i += replacementFormula.Length() - name.Length();
1998
1999 // we have extracted all the functor for "fname"
2000 // std::cout << "We have extracted all the functors for fname" << std::endl;
2001 // std::cout << " i = " << i << " f[i] = " << formula[i] << " - " << formula << std::endl;
2002 name = "";
2003
2004 continue;
2005 }
2006
2007 // add now functor in
2008 TString replacement = TString::Format("{%s}", name.Data());
2009 formula.Replace(i - name.Length(), name.Length(), replacement, replacement.Length());
2010 i += 2;
2011 fFuncs.push_back(TFormulaFunction(name));
2012 }
2013 }
2014 name = body = "";
2015 }
2016}
2017
2018////////////////////////////////////////////////////////////////////////////////
2019/// Iterates through functors in fFuncs and performs the appropriate action.
2020/// If functor has 0 arguments (has only name) can be:
2021/// - variable
2022/// * will be replaced with x[num], where x is an array containing value of this variable under num.
2023/// - pre-defined formula
2024/// * will be replaced with formulas body
2025/// - constant
2026/// * will be replaced with constant value
2027/// - parameter
2028/// * will be replaced with p[num], where p is an array containing value of this parameter under num.
2029/// If has arguments it can be :
2030/// - function shortcut, eg. sin
2031/// * will be replaced with fullname of function, eg. sin -> TMath::Sin
2032/// - function from cling environment, eg. TMath::BreitWigner(x,y,z)
2033/// * first check if function exists, and has same number of arguments, then accept it and set as found.
2034/// If all functors after iteration are matched with corresponding action,
2035/// it inputs C++ code of formula into cling, and sets flag that formula is ready to evaluate.
2036
2038{
2039 // std::cout << "Begin: formula is " << formula << " list of functors " << fFuncs.size() << std::endl;
2040
2041 for (list<TFormulaFunction>::iterator funcsIt = fFuncs.begin(); funcsIt != fFuncs.end(); ++funcsIt) {
2042 TFormulaFunction &fun = *funcsIt;
2043
2044 // std::cout << "fun is " << fun.GetName() << std::endl;
2045
2046 if (fun.fFound)
2047 continue;
2048 if (fun.IsFuncCall()) {
2049 // replace with pre-defined functions
2050 map<TString, TString>::iterator it = fFunctionsShortcuts.find(fun.GetName());
2051 if (it != fFunctionsShortcuts.end()) {
2052 TString shortcut = it->first;
2053 TString full = it->second;
2054 // std::cout << " functor " << fun.GetName() << " found - replace " << shortcut << " with " << full << " in
2055 // " << formula << std::endl;
2056 // replace all functors
2057 Ssiz_t index = formula.Index(shortcut, 0);
2058 while (index != kNPOS) {
2059 // check that function is not in a namespace and is not in other characters
2060 // std::cout << "analyzing " << shortcut << " in " << formula << std::endl;
2061 Ssiz_t i2 = index + shortcut.Length();
2062 if ((index > 0) && (isalpha(formula[index - 1]) || formula[index - 1] == ':')) {
2063 index = formula.Index(shortcut, i2);
2064 continue;
2065 }
2066 if (i2 < formula.Length() && formula[i2] != '(') {
2067 index = formula.Index(shortcut, i2);
2068 continue;
2069 }
2070 // now replace the string
2071 formula.Replace(index, shortcut.Length(), full);
2072 Ssiz_t inext = index + full.Length();
2073 index = formula.Index(shortcut, inext);
2074 fun.fFound = true;
2075 }
2076 }
2077 // for functions we can live it to cling to decide if it is a valid function or NOT
2078 // We don't need to retrieve this information from the ROOT interpreter
2079 // we assume that the function is then found and all the following code does not need to be there
2080#ifdef TFORMULA_CHECK_FUNCTIONS
2081
2082 if (fun.fName.Contains("::")) // add support for nested namespaces
2083 {
2084 // look for last occurence of "::"
2085 std::string name(fun.fName.Data());
2086 size_t index = name.rfind("::");
2087 assert(index != std::string::npos);
2088 TString className = fun.fName(0, fun.fName(0, index).Length());
2089 TString functionName = fun.fName(index + 2, fun.fName.Length());
2090
2091 Bool_t silent = true;
2092 TClass *tclass = TClass::GetClass(className, silent);
2093 // std::cout << "looking for class " << className << std::endl;
2094 const TList *methodList = tclass->GetListOfAllPublicMethods();
2095 TIter next(methodList);
2096 TMethod *p;
2097 while ((p = (TMethod *)next())) {
2098 if (strcmp(p->GetName(), functionName.Data()) == 0 &&
2099 (fun.GetNargs() <= p->GetNargs() && fun.GetNargs() >= p->GetNargs() - p->GetNargsOpt())) {
2100 fun.fFound = true;
2101 break;
2102 }
2103 }
2104 }
2105 if (!fun.fFound) {
2106 // try to look into all the global functions in gROOT
2107 TFunction *f;
2108 {
2110 f = (TFunction *)gROOT->GetListOfGlobalFunctions(true)->FindObject(fun.fName);
2111 }
2112 // if found a function with matching arguments
2113 if (f && fun.GetNargs() <= f->GetNargs() && fun.GetNargs() >= f->GetNargs() - f->GetNargsOpt()) {
2114 fun.fFound = true;
2115 }
2116 }
2117
2118 if (!fun.fFound) {
2119 // ignore not found functions
2120 if (gDebug)
2121 Info("TFormula", "Could not find %s function with %d argument(s)", fun.GetName(), fun.GetNargs());
2122 fun.fFound = false;
2123 }
2124#endif
2125 } else {
2126 TFormula *old = 0;
2127 {
2129 old = (TFormula *)gROOT->GetListOfFunctions()->FindObject(gNamePrefix + fun.fName);
2130 }
2131 if (old) {
2132 // we should not go here (this analysis is done before in ExtractFunctors)
2133 assert(false);
2134 fun.fFound = true;
2135 TString pattern = TString::Format("{%s}", fun.GetName());
2136 TString replacement = old->GetExpFormula();
2137 PreProcessFormula(replacement);
2138 ExtractFunctors(replacement);
2139 formula.ReplaceAll(pattern, replacement);
2140 continue;
2141 }
2142 // looking for default variables defined in fVars
2143
2144 map<TString, TFormulaVariable>::iterator varsIt = fVars.find(fun.GetName());
2145 if (varsIt != fVars.end()) {
2146
2147 TString name = (*varsIt).second.GetName();
2148 Double_t value = (*varsIt).second.fValue;
2149
2150 AddVariable(name, value); // this set the cling variable
2151 if (!fVars[name].fFound) {
2152
2153 fVars[name].fFound = true;
2154 int varDim = (*varsIt).second.fArrayPos; // variable dimensions (0 for x, 1 for y, 2, for z)
2155 if (varDim >= fNdim) {
2156 fNdim = varDim + 1;
2157
2158 // we need to be sure that all other variables are added with position less
2159 for (auto &v : fVars) {
2160 if (v.second.fArrayPos < varDim && !v.second.fFound) {
2161 AddVariable(v.first, v.second.fValue);
2162 v.second.fFound = true;
2163 }
2164 }
2165 }
2166 }
2167 // remove the "{.. }" added around the variable
2168 TString pattern = TString::Format("{%s}", name.Data());
2169 TString replacement = TString::Format("x[%d]", (*varsIt).second.fArrayPos);
2170 formula.ReplaceAll(pattern, replacement);
2171
2172 // std::cout << "Found an observable for " << fun.GetName() << std::endl;
2173
2174 fun.fFound = true;
2175 continue;
2176 }
2177 // check for observables defined as x[0],x[1],....
2178 // maybe could use a regular expression here
2179 // only in case match with defined variables is not successful
2180 TString funname = fun.GetName();
2181 if (funname.Contains("x[") && funname.Contains("]")) {
2182 TString sdigit = funname(2, funname.Index("]"));
2183 int digit = sdigit.Atoi();
2184 if (digit >= fNdim) {
2185 fNdim = digit + 1;
2186 // we need to add the variables in fVars all of them before x[n]
2187 for (int j = 0; j < fNdim; ++j) {
2188 TString vname = TString::Format("x[%d]", j);
2189 if (fVars.find(vname) == fVars.end()) {
2190 fVars[vname] = TFormulaVariable(vname, 0, j);
2191 fVars[vname].fFound = true;
2192 AddVariable(vname, 0.);
2193 }
2194 }
2195 }
2196 // std::cout << "Found matching observable for " << funname << std::endl;
2197 fun.fFound = true;
2198 // remove the "{.. }" added around the variable
2199 TString pattern = TString::Format("{%s}", funname.Data());
2200 formula.ReplaceAll(pattern, funname);
2201 continue;
2202 }
2203 //}
2204
2205 auto paramsIt = fParams.find(fun.GetName());
2206 if (paramsIt != fParams.end()) {
2207 // TString name = (*paramsIt).second.GetName();
2208 TString pattern = TString::Format("{[%s]}", fun.GetName());
2209 // std::cout << "pattern is " << pattern << std::endl;
2210 if (formula.Index(pattern) != kNPOS) {
2211 // TString replacement = TString::Format("p[%d]",(*paramsIt).second.fArrayPos);
2212 TString replacement = TString::Format("p[%d]", (*paramsIt).second);
2213 // std::cout << "replace pattern " << pattern << " with " << replacement << std::endl;
2214 formula.ReplaceAll(pattern, replacement);
2215 }
2216 fun.fFound = true;
2217 continue;
2218 } else {
2219 // std::cout << "functor " << fun.GetName() << " is not a parameter " << std::endl;
2220 }
2221
2222 // looking for constants (needs to be done after looking at the parameters)
2223 map<TString, Double_t>::iterator constIt = fConsts.find(fun.GetName());
2224 if (constIt != fConsts.end()) {
2225 TString pattern = TString::Format("{%s}", fun.GetName());
2226 TString value = TString::Format("%lf", (*constIt).second);
2227 formula.ReplaceAll(pattern, value);
2228 fun.fFound = true;
2229 // std::cout << "constant with name " << fun.GetName() << " is found " << std::endl;
2230 continue;
2231 }
2232
2233 fun.fFound = false;
2234 }
2235 }
2236 // std::cout << "End: formula is " << formula << std::endl;
2237
2238 // ignore case of functors have been matched - try to pass it to Cling
2239 if (!fReadyToExecute) {
2240 fReadyToExecute = true;
2241 Bool_t hasVariables = (fNdim > 0);
2242 Bool_t hasParameters = (fNpar > 0);
2243 if (!hasParameters) {
2244 fAllParametersSetted = true;
2245 }
2246 // assume a function without variables is always 1-dimensional ???
2247 // if (hasParameters && !hasVariables) {
2248 // fNdim = 1;
2249 // AddVariable("x", 0);
2250 // hasVariables = true;
2251 // }
2252 // does not make sense to vectorize function which is of FNDim=0
2253 if (!hasVariables) fVectorized=false;
2254 // when there are no variables but only parameter we still need to ad
2255 //Bool_t hasBoth = hasVariables && hasParameters;
2256 Bool_t inputIntoCling = (formula.Length() > 0);
2257 if (inputIntoCling) {
2258 // save copy of inputFormula in a std::strig for the unordered map
2259 // and also formula is same as FClingInput typically and it will be modified
2260 std::string inputFormula(formula.Data());
2261
2262 // The name we really use for the unordered map will have a flag that
2263 // says whether the formula is vectorized
2264 std::string inputFormulaVecFlag = inputFormula;
2265 if (fVectorized)
2266 inputFormulaVecFlag += " (vectorized)";
2267
2268 TString argType = fVectorized ? "ROOT::Double_v" : "Double_t";
2269
2270 // valid input formula - try to put into Cling (in case of no variables but only parameter we need to add the standard signature)
2271 TString argumentsPrototype = TString::Format("%s%s%s", ( (hasVariables || hasParameters) ? (argType + " *x").Data() : ""),
2272 (hasParameters ? "," : ""), (hasParameters ? "Double_t *p" : ""));
2273
2274 // set the name for Cling using the hash_function
2276
2277 // check if formula exist already in the map
2279
2280 // std::cout << "gClingFunctions list" << std::endl;
2281 // for (auto thing : gClingFunctions)
2282 // std::cout << "gClingFunctions : " << thing.first << std::endl;
2283
2284 auto funcit = gClingFunctions.find(inputFormulaVecFlag);
2285
2286 if (funcit != gClingFunctions.end()) {
2287 fFuncPtr = (TFormula::CallFuncSignature)funcit->second;
2288 fClingInitialized = true;
2289 inputIntoCling = false;
2290 }
2291
2292
2293
2294 // set the cling name using hash of the static formulae map
2295 auto hasher = gClingFunctions.hash_function();
2296 fClingName = TString::Format("%s__id%zu", gNamePrefix.Data(), hasher(inputFormulaVecFlag));
2297
2298 fClingInput = TString::Format("%s %s(%s){ return %s ; }", argType.Data(), fClingName.Data(),
2299 argumentsPrototype.Data(), inputFormula.c_str());
2300
2301
2302 // std::cout << "Input Formula " << inputFormula << " \t vec formula : " << inputFormulaVecFlag << std::endl;
2303 // std::cout << "Cling functions existing " << std::endl;
2304 // for (auto & ff : gClingFunctions)
2305 // std::cout << ff.first << std::endl;
2306 // std::cout << "\n";
2307 // std::cout << fClingName << std::endl;
2308
2309 // this is not needed (maybe can be re-added in case of recompilation of identical expressions
2310 // // check in case of a change if need to re-initialize
2311 // if (fClingInitialized) {
2312 // if (oldClingInput == fClingInput)
2313 // inputIntoCling = false;
2314 // else
2315 // fClingInitialized = false;
2316 // }
2317
2318 if (inputIntoCling) {
2319 if (!fLazyInitialization) {
2321 if (fClingInitialized) {
2322 // if Cling has been successfully initialized
2323 // put function ptr in the static map
2325 gClingFunctions.insert(std::make_pair(inputFormulaVecFlag, (void *)fFuncPtr));
2326 }
2327 }
2328 if (!fClingInitialized) {
2329 // needed in case of lazy initialization of failure compiling the expression
2330 fSavedInputFormula = inputFormulaVecFlag;
2331 }
2332
2333 } else {
2334 fAllParametersSetted = true;
2335 fClingInitialized = true;
2336 }
2337 }
2338 }
2339
2340 // In case of a Cling Error check components which are not found in Cling
2341 // check that all formula components are matched otherwise emit an error
2343 //Bool_t allFunctorsMatched = false;
2344 for (list<TFormulaFunction>::iterator it = fFuncs.begin(); it != fFuncs.end(); ++it) {
2345 // functions are now by default always not checked
2346 if (!it->fFound && !it->IsFuncCall()) {
2347 //allFunctorsMatched = false;
2348 if (it->GetNargs() == 0)
2349 Error("ProcessFormula", "\"%s\" has not been matched in the formula expression", it->GetName());
2350 else
2351 Error("ProcessFormula", "Could not find %s function with %d argument(s)", it->GetName(), it->GetNargs());
2352 }
2353 }
2354 Error("ProcessFormula","Formula \"%s\" is invalid !", GetExpFormula().Data() );
2355 fReadyToExecute = false;
2356 }
2357
2358 // clean up un-used default variables in case formula is valid
2359 //if (fClingInitialized && fReadyToExecute) {
2360 //don't check fClingInitialized in case of lazy execution
2361 if (fReadyToExecute) {
2362 auto itvar = fVars.begin();
2363 // need this loop because after erase change iterators
2364 do {
2365 if (!itvar->second.fFound) {
2366 // std::cout << "Erase variable " << itvar->first << std::endl;
2367 itvar = fVars.erase(itvar);
2368 } else
2369 itvar++;
2370 } while (itvar != fVars.end());
2371 }
2372}
2373
2374////////////////////////////////////////////////////////////////////////////////
2375/// Fill map with parametrized functions
2376
2377void TFormula::FillParametrizedFunctions(map<pair<TString, Int_t>, pair<TString, TString>> &functions)
2378{
2379 // map< pair<TString,Int_t> ,pair<TString,TString> > functions;
2380 functions.insert(
2381 make_pair(make_pair("gaus", 1), make_pair("[0]*exp(-0.5*(({V0}-[1])/[2])*(({V0}-[1])/[2]))",
2382 "[0]*exp(-0.5*(({V0}-[1])/[2])*(({V0}-[1])/[2]))/(sqrt(2*pi)*[2])")));
2383 functions.insert(make_pair(make_pair("landau", 1), make_pair("[0]*TMath::Landau({V0},[1],[2],false)",
2384 "[0]*TMath::Landau({V0},[1],[2],true)")));
2385 functions.insert(make_pair(make_pair("expo", 1), make_pair("exp([0]+[1]*{V0})", "")));
2386 functions.insert(
2387 make_pair(make_pair("crystalball", 1), make_pair("[0]*ROOT::Math::crystalball_function({V0},[3],[4],[2],[1])",
2388 "[0]*ROOT::Math::crystalball_pdf({V0},[3],[4],[2],[1])")));
2389 functions.insert(
2390 make_pair(make_pair("breitwigner", 1), make_pair("[0]*ROOT::Math::breitwigner_pdf({V0},[2],[1])",
2391 "[0]*ROOT::Math::breitwigner_pdf({V0},[2],[4],[1])")));
2392 // chebyshev polynomial
2393 functions.insert(make_pair(make_pair("cheb0", 1), make_pair("ROOT::Math::Chebyshev0({V0},[0])", "")));
2394 functions.insert(make_pair(make_pair("cheb1", 1), make_pair("ROOT::Math::Chebyshev1({V0},[0],[1])", "")));
2395 functions.insert(make_pair(make_pair("cheb2", 1), make_pair("ROOT::Math::Chebyshev2({V0},[0],[1],[2])", "")));
2396 functions.insert(make_pair(make_pair("cheb3", 1), make_pair("ROOT::Math::Chebyshev3({V0},[0],[1],[2],[3])", "")));
2397 functions.insert(
2398 make_pair(make_pair("cheb4", 1), make_pair("ROOT::Math::Chebyshev4({V0},[0],[1],[2],[3],[4])", "")));
2399 functions.insert(
2400 make_pair(make_pair("cheb5", 1), make_pair("ROOT::Math::Chebyshev5({V0},[0],[1],[2],[3],[4],[5])", "")));
2401 functions.insert(
2402 make_pair(make_pair("cheb6", 1), make_pair("ROOT::Math::Chebyshev6({V0},[0],[1],[2],[3],[4],[5],[6])", "")));
2403 functions.insert(
2404 make_pair(make_pair("cheb7", 1), make_pair("ROOT::Math::Chebyshev7({V0},[0],[1],[2],[3],[4],[5],[6],[7])", "")));
2405 functions.insert(make_pair(make_pair("cheb8", 1),
2406 make_pair("ROOT::Math::Chebyshev8({V0},[0],[1],[2],[3],[4],[5],[6],[7],[8])", "")));
2407 functions.insert(make_pair(make_pair("cheb9", 1),
2408 make_pair("ROOT::Math::Chebyshev9({V0},[0],[1],[2],[3],[4],[5],[6],[7],[8],[9])", "")));
2409 functions.insert(
2410 make_pair(make_pair("cheb10", 1),
2411 make_pair("ROOT::Math::Chebyshev10({V0},[0],[1],[2],[3],[4],[5],[6],[7],[8],[9],[10])", "")));
2412 // 2-dimensional functions
2413 functions.insert(
2414 make_pair(make_pair("gaus", 2), make_pair("[0]*exp(-0.5*(({V0}-[1])/[2])^2 - 0.5*(({V1}-[3])/[4])^2)", "")));
2415 functions.insert(
2416 make_pair(make_pair("landau", 2),
2417 make_pair("[0]*TMath::Landau({V0},[1],[2],false)*TMath::Landau({V1},[3],[4],false)", "")));
2418 functions.insert(make_pair(make_pair("expo", 2), make_pair("exp([0]+[1]*{V0})", "exp([0]+[1]*{V0}+[2]*{V1})")));
2419 // 3-dimensional function
2420 functions.insert(
2421 make_pair(make_pair("gaus", 3), make_pair("[0]*exp(-0.5*(({V0}-[1])/[2])^2 - 0.5*(({V1}-[3])/[4])^2 - 0.5*(({V2}-[5])/[6])^2)", "")));
2422 // gaussian with correlations
2423 functions.insert(
2424 make_pair(make_pair("bigaus", 2), make_pair("[0]*ROOT::Math::bigaussian_pdf({V0},{V1},[2],[4],[5],[1],[3])",
2425 "[0]*ROOT::Math::bigaussian_pdf({V0},{V1},[2],[4],[5],[1],[3])")));
2426}
2427
2428////////////////////////////////////////////////////////////////////////////////
2429/// Set parameter names only in case of pre-defined functions.
2430
2432
2433 if (fNumber == 0) return;
2434
2435 if (fNumber == 100) { // Gaussian
2436 SetParName(0,"Constant");
2437 SetParName(1,"Mean");
2438 SetParName(2,"Sigma");
2439 return;
2440 }
2441 if (fNumber == 110) {
2442 SetParName(0,"Constant");
2443 SetParName(1,"MeanX");
2444 SetParName(2,"SigmaX");
2445 SetParName(3,"MeanY");
2446 SetParName(4,"SigmaY");
2447 return;
2448 }
2449 if (fNumber == 120) {
2450 SetParName(0,"Constant");
2451 SetParName(1,"MeanX");
2452 SetParName(2,"SigmaX");
2453 SetParName(3,"MeanY");
2454 SetParName(4,"SigmaY");
2455 SetParName(5,"MeanZ");
2456 SetParName(6,"SigmaZ");
2457 return;
2458 }
2459 if (fNumber == 112) { // bigaus
2460 SetParName(0,"Constant");
2461 SetParName(1,"MeanX");
2462 SetParName(2,"SigmaX");
2463 SetParName(3,"MeanY");
2464 SetParName(4,"SigmaY");
2465 SetParName(5,"Rho");
2466 return;
2467 }
2468 if (fNumber == 200) { // exponential
2469 SetParName(0,"Constant");
2470 SetParName(1,"Slope");
2471 return;
2472 }
2473 if (fNumber == 400) { // landau
2474 SetParName(0,"Constant");
2475 SetParName(1,"MPV");
2476 SetParName(2,"Sigma");
2477 return;
2478 }
2479 if (fNumber == 500) { // crystal-ball
2480 SetParName(0,"Constant");
2481 SetParName(1,"Mean");
2482 SetParName(2,"Sigma");
2483 SetParName(3,"Alpha");
2484 SetParName(4,"N");
2485 return;
2486 }
2487 if (fNumber == 600) { // breit-wigner
2488 SetParName(0,"Constant");
2489 SetParName(1,"Mean");
2490 SetParName(2,"Gamma");
2491 return;
2492 }
2493 // if formula is a polynomial (or chebyshev), set parameter names
2494 // not needed anymore (p0 is assigned by default)
2495 // if (fNumber == (300+fNpar-1) ) {
2496 // for (int i = 0; i < fNpar; i++) SetParName(i,TString::Format("p%d",i));
2497 // return;
2498 // }
2499
2500 // // general case if parameters are digits (XX) change to pXX
2501 // auto paramMap = fParams; // need to copy the map because SetParName is going to modify it
2502 // for ( auto & p : paramMap) {
2503 // if (p.first.IsDigit() )
2504 // SetParName(p.second,TString::Format("p%s",p.first.Data()));
2505 // }
2506
2507 return;
2508}
2509
2510////////////////////////////////////////////////////////////////////////////////
2511/// Return linear part.
2512
2514{
2515 if (!fLinearParts.empty()) {
2516 int n = fLinearParts.size();
2517 if (i < 0 || i >= n ) {
2518 Error("GetLinearPart","Formula %s has only %d linear parts - requested %d",GetName(),n,i);
2519 return nullptr;
2520 }
2521 return fLinearParts[i];
2522 }
2523 return nullptr;
2524}
2525
2526////////////////////////////////////////////////////////////////////////////////
2527/// Adds variable to known variables, and reprocess formula.
2528
2529void TFormula::AddVariable(const TString &name, double value)
2530{
2531 if (fVars.find(name) != fVars.end()) {
2532 TFormulaVariable &var = fVars[name];
2533 var.fValue = value;
2534
2535 // If the position is not defined in the Cling vectors, make space for it
2536 // but normally is variable is defined in fVars a slot should be also present in fClingVariables
2537 if (var.fArrayPos < 0) {
2538 var.fArrayPos = fVars.size();
2539 }
2540 if (var.fArrayPos >= (int)fClingVariables.size()) {
2541 fClingVariables.resize(var.fArrayPos + 1);
2542 }
2543 fClingVariables[var.fArrayPos] = value;
2544 } else {
2545 TFormulaVariable var(name, value, fVars.size());
2546 fVars[name] = var;
2547 fClingVariables.push_back(value);
2548 if (!fFormula.IsNull()) {
2549 // printf("process formula again - %s \n",fClingInput.Data() );
2551 }
2552 }
2553}
2554
2555////////////////////////////////////////////////////////////////////////////////
2556/// Adds multiple variables.
2557/// First argument is an array of pairs<TString,Double>, where
2558/// first argument is name of variable,
2559/// second argument represents value.
2560/// size - number of variables passed in first argument
2561
2562void TFormula::AddVariables(const TString *vars, const Int_t size)
2563{
2564 Bool_t anyNewVar = false;
2565 for (Int_t i = 0; i < size; ++i) {
2566
2567 const TString &vname = vars[i];
2568
2569 TFormulaVariable &var = fVars[vname];
2570 if (var.fArrayPos < 0) {
2571
2572 var.fName = vname;
2573 var.fArrayPos = fVars.size();
2574 anyNewVar = true;
2575 var.fValue = 0;
2576 if (var.fArrayPos >= (int)fClingVariables.capacity()) {
2577 Int_t multiplier = 2;
2578 if (fFuncs.size() > 100) {
2579 multiplier = TMath::Floor(TMath::Log10(fFuncs.size()) * 10);
2580 }
2581 fClingVariables.reserve(multiplier * fClingVariables.capacity());
2582 }
2583 fClingVariables.push_back(0.0);
2584 }
2585 // else
2586 // {
2587 // var.fValue = v.second;
2588 // fClingVariables[var.fArrayPos] = v.second;
2589 // }
2590 }
2591 if (anyNewVar && !fFormula.IsNull()) {
2593 }
2594}
2595
2596////////////////////////////////////////////////////////////////////////////////
2597/// Set the name of the formula. We need to allow the list of function to
2598/// properly handle the hashes.
2599
2600void TFormula::SetName(const char* name)
2601{
2602 if (IsReservedName(name)) {
2603 Error("SetName", "The name \'%s\' is reserved as a TFormula variable name.\n"
2604 "\tThis function will not be renamed.",
2605 name);
2606 } else {
2607 // Here we need to remove and re-add to keep the hashes consistent with
2608 // the underlying names.
2609 auto listOfFunctions = gROOT->GetListOfFunctions();
2610 TObject* thisAsFunctionInList = nullptr;
2612 if (listOfFunctions){
2613 thisAsFunctionInList = listOfFunctions->FindObject(this);
2614 if (thisAsFunctionInList) listOfFunctions->Remove(thisAsFunctionInList);
2615 }
2617 if (thisAsFunctionInList) listOfFunctions->Add(thisAsFunctionInList);
2618 }
2619}
2620
2621////////////////////////////////////////////////////////////////////////////////
2622///
2623/// Sets multiple variables.
2624/// First argument is an array of pairs<TString,Double>, where
2625/// first argument is name of variable,
2626/// second argument represents value.
2627/// size - number of variables passed in first argument
2628
2629void TFormula::SetVariables(const pair<TString,Double_t> *vars, const Int_t size)
2630{
2631 for(Int_t i = 0; i < size; ++i)
2632 {
2633 auto &v = vars[i];
2634 if (fVars.find(v.first) != fVars.end()) {
2635 fVars[v.first].fValue = v.second;
2636 fClingVariables[fVars[v.first].fArrayPos] = v.second;
2637 } else {
2638 Error("SetVariables", "Variable %s is not defined.", v.first.Data());
2639 }
2640 }
2641}
2642
2643////////////////////////////////////////////////////////////////////////////////
2644/// Returns variable value.
2645
2647{
2648 const auto nameIt = fVars.find(name);
2649 if (fVars.end() == nameIt) {
2650 Error("GetVariable", "Variable %s is not defined.", name);
2651 return -1;
2652 }
2653 return nameIt->second.fValue;
2654}
2655
2656////////////////////////////////////////////////////////////////////////////////
2657/// Returns variable number (positon in array) given its name.
2658
2660{
2661 const auto nameIt = fVars.find(name);
2662 if (fVars.end() == nameIt) {
2663 Error("GetVarNumber", "Variable %s is not defined.", name);
2664 return -1;
2665 }
2666 return nameIt->second.fArrayPos;
2667}
2668
2669////////////////////////////////////////////////////////////////////////////////
2670/// Returns variable name given its position in the array.
2671
2673{
2674 if (ivar < 0 || ivar >= fNdim) return "";
2675
2676 // need to loop on the map to find corresponding variable
2677 for ( auto & v : fVars) {
2678 if (v.second.fArrayPos == ivar) return v.first;
2679 }
2680 Error("GetVarName","Variable with index %d not found !!",ivar);
2681 //return TString::Format("x%d",ivar);
2682 return "";
2683}
2684
2685////////////////////////////////////////////////////////////////////////////////
2686/// Sets variable value.
2687
2689{
2690 if (fVars.find(name) == fVars.end()) {
2691 Error("SetVariable", "Variable %s is not defined.", name.Data());
2692 return;
2693 }
2694 fVars[name].fValue = value;
2695 fClingVariables[fVars[name].fArrayPos] = value;
2696}
2697
2698////////////////////////////////////////////////////////////////////////////////
2699/// Adds parameter to known parameters.
2700/// User should use SetParameter, because parameters are added during initialization part,
2701/// and after that adding new will be pointless.
2702
2703void TFormula::DoAddParameter(const TString &name, Double_t value, Bool_t processFormula)
2704{
2705 //std::cout << "adding parameter " << name << std::endl;
2706
2707 // if parameter is already defined in fParams - just set the new value
2708 if(fParams.find(name) != fParams.end() )
2709 {
2710 int ipos = fParams[name];
2711 // TFormulaVariable & par = fParams[name];
2712 // par.fValue = value;
2713 if (ipos < 0) {
2714 ipos = fParams.size();
2715 fParams[name] = ipos;
2716 }
2717 //
2718 if (ipos >= (int)fClingParameters.size()) {
2719 if (ipos >= (int)fClingParameters.capacity())
2720 fClingParameters.reserve(TMath::Max(int(fParams.size()), ipos + 1));
2721 fClingParameters.insert(fClingParameters.end(), ipos + 1 - fClingParameters.size(), 0.0);
2722 }
2723 fClingParameters[ipos] = value;
2724 } else {
2725 // new parameter defined
2726 fNpar++;
2727 // TFormulaVariable(name,value,fParams.size());
2728 int pos = fParams.size();
2729 // fParams.insert(std::make_pair<TString,TFormulaVariable>(name,TFormulaVariable(name,value,pos)));
2730 auto ret = fParams.insert(std::make_pair(name, pos));
2731 // map returns a std::pair<iterator, bool>
2732 // use map the order for default position of parameters in the vector
2733 // (i.e use the alphabetic order)
2734 if (ret.second) {
2735 // a new element is inserted
2736 if (ret.first == fParams.begin())
2737 pos = 0;
2738 else {
2739 auto previous = (ret.first);
2740 --previous;
2741 pos = previous->second + 1;
2742 }
2743
2744 if (pos < (int)fClingParameters.size())
2745 fClingParameters.insert(fClingParameters.begin() + pos, value);
2746 else {
2747 // this should not happen
2748 if (pos > (int)fClingParameters.size())
2749 Warning("inserting parameter %s at pos %d when vector size is %d \n", name.Data(), pos,
2750 (int)fClingParameters.size());
2751
2752 if (pos >= (int)fClingParameters.capacity())
2753 fClingParameters.reserve(TMath::Max(int(fParams.size()), pos + 1));
2754 fClingParameters.insert(fClingParameters.end(), pos + 1 - fClingParameters.size(), 0.0);
2755 fClingParameters[pos] = value;
2756 }
2757
2758 // need to adjust all other positions
2759 for (auto it = ret.first; it != fParams.end(); ++it) {
2760 it->second = pos;
2761 pos++;
2762 }
2763
2764 // for (auto & p : fParams)
2765 // std::cout << "Parameter " << p.first << " position " << p.second << " value " <<
2766 // fClingParameters[p.second] << std::endl;
2767 // printf("inserted parameters size params %d size cling %d \n",fParams.size(), fClingParameters.size() );
2768 }
2769 if (processFormula) {
2770 // replace first in input parameter name with [name]
2773 }
2774 }
2775}
2776
2777////////////////////////////////////////////////////////////////////////////////
2778/// Return parameter index given a name (return -1 for not existing parameters)
2779/// non need to print an error
2780
2781Int_t TFormula::GetParNumber(const char * name) const {
2782 auto it = fParams.find(name);
2783 if (it == fParams.end()) {
2784 return -1;
2785 }
2786 return it->second;
2787
2788}
2789
2790////////////////////////////////////////////////////////////////////////////////
2791/// Returns parameter value given by string.
2792
2794{
2795 const int i = GetParNumber(name);
2796 if (i == -1) {
2797 Error("GetParameter","Parameter %s is not defined.",name);
2798 return TMath::QuietNaN();
2799 }
2800
2801 return GetParameter( i );
2802}
2803
2804////////////////////////////////////////////////////////////////////////////////
2805/// Return parameter value given by integer.
2806
2808{
2809 //TString name = TString::Format("%d",param);
2810 if(param >=0 && param < (int) fClingParameters.size())
2811 return fClingParameters[param];
2812 Error("GetParameter","wrong index used - use GetParameter(name)");
2813 return TMath::QuietNaN();
2814}
2815
2816////////////////////////////////////////////////////////////////////////////////
2817/// Return parameter name given by integer.
2818
2819const char * TFormula::GetParName(Int_t ipar) const
2820{
2821 if (ipar < 0 || ipar >= fNpar) return "";
2822
2823 // need to loop on the map to find corresponding parameter
2824 for ( auto & p : fParams) {
2825 if (p.second == ipar) return p.first.Data();
2826 }
2827 Error("GetParName","Parameter with index %d not found !!",ipar);
2828 //return TString::Format("p%d",ipar);
2829 return "";
2830}
2831
2832////////////////////////////////////////////////////////////////////////////////
2834{
2835 if(!fClingParameters.empty())
2836 return const_cast<Double_t*>(&fClingParameters[0]);
2837 return 0;
2838}
2839
2841{
2842 for (Int_t i = 0; i < fNpar; ++i) {
2843 if (Int_t(fClingParameters.size()) > i)
2844 params[i] = fClingParameters[i];
2845 else
2846 params[i] = -1;
2847 }
2848}
2849
2850////////////////////////////////////////////////////////////////////////////////
2851/// Sets parameter value.
2852
2853void TFormula::SetParameter(const char *name, Double_t value)
2854{
2855 SetParameter( GetParNumber(name), value);
2856
2857 // do we need this ???
2858#ifdef OLDPARAMS
2859 if (fParams.find(name) == fParams.end()) {
2860 Error("SetParameter", "Parameter %s is not defined.", name.Data());
2861 return;
2862 }
2863 fParams[name].fValue = value;
2864 fParams[name].fFound = true;
2865 fClingParameters[fParams[name].fArrayPos] = value;
2866 fAllParametersSetted = true;
2867 for (map<TString, TFormulaVariable>::iterator it = fParams.begin(); it != fParams.end(); ++it) {
2868 if (!it->second.fFound) {
2869 fAllParametersSetted = false;
2870 break;
2871 }
2872 }
2873#endif
2874}
2875
2876#ifdef OLDPARAMS
2877
2878////////////////////////////////////////////////////////////////////////////////
2879/// Set multiple parameters.
2880/// First argument is an array of pairs<TString,Double>, where
2881/// first argument is name of parameter,
2882/// second argument represents value.
2883/// size - number of params passed in first argument
2884
2885void TFormula::SetParameters(const pair<TString,Double_t> *params,const Int_t size)
2886{
2887 for(Int_t i = 0 ; i < size ; ++i)
2888 {
2889 pair<TString, Double_t> p = params[i];
2890 if (fParams.find(p.first) == fParams.end()) {
2891 Error("SetParameters", "Parameter %s is not defined", p.first.Data());
2892 continue;
2893 }
2894 fParams[p.first].fValue = p.second;
2895 fParams[p.first].fFound = true;
2896 fClingParameters[fParams[p.first].fArrayPos] = p.second;
2897 }
2898 fAllParametersSetted = true;
2899 for (map<TString, TFormulaVariable>::iterator it = fParams.begin(); it != fParams.end(); ++it) {
2900 if (!it->second.fFound) {
2901 fAllParametersSetted = false;
2902 break;
2903 }
2904 }
2905}
2906#endif
2907
2908////////////////////////////////////////////////////////////////////////////////
2910{
2911 if(!params || size < 0 || size > fNpar) return;
2912 // reset vector of cling parameters
2913 if (size != (int) fClingParameters.size() ) {
2914 Warning("SetParameters","size is not same of cling parameter size %d - %d",size,int(fClingParameters.size()) );
2915 for (Int_t i = 0; i < size; ++i) {
2916 TString name = TString::Format("%d", i);
2917 SetParameter(name, params[i]);
2918 }
2919 return;
2920 }
2921 fAllParametersSetted = true;
2922 std::copy(params, params+size, fClingParameters.begin() );
2923}
2924
2925////////////////////////////////////////////////////////////////////////////////
2926/// Set a vector of parameters value.
2927/// Order in the vector is by default the alphabetic order given to the parameters
2928/// apart if the users has defined explicitly the parameter names
2929
2931{
2932 DoSetParameters(params,fNpar);
2933}
2934
2935////////////////////////////////////////////////////////////////////////////////
2936/// Set a list of parameters.
2937/// The order is by default the alphabetic order given to the parameters
2938/// apart if the users has defined explicitly the parameter names
2939
2941 Double_t p7, Double_t p8, Double_t p9, Double_t p10)
2942{
2943 if(fNpar >= 1) SetParameter(0,p0);
2944 if(fNpar >= 2) SetParameter(1,p1);
2945 if(fNpar >= 3) SetParameter(2,p2);
2946 if(fNpar >= 4) SetParameter(3,p3);
2947 if(fNpar >= 5) SetParameter(4,p4);
2948 if(fNpar >= 6) SetParameter(5,p5);
2949 if(fNpar >= 7) SetParameter(6,p6);
2950 if(fNpar >= 8) SetParameter(7,p7);
2951 if(fNpar >= 9) SetParameter(8,p8);
2952 if(fNpar >= 10) SetParameter(9,p9);
2953 if(fNpar >= 11) SetParameter(10,p10);
2954}
2955
2956////////////////////////////////////////////////////////////////////////////////
2957/// Set a parameter given a parameter index
2958/// The parameter index is by default the alphabetic order given to the parameters
2959/// apart if the users has defined explicitly the parameter names
2960
2962{
2963 if (param < 0 || param >= fNpar) return;
2964 assert(int(fClingParameters.size()) == fNpar);
2965 fClingParameters[param] = value;
2966 // TString name = TString::Format("%d",param);
2967 // SetParameter(name,value);
2968}
2969
2970////////////////////////////////////////////////////////////////////////////////
2971void TFormula::SetParNames(const char *name0, const char *name1, const char *name2, const char *name3,
2972 const char *name4, const char *name5, const char *name6, const char *name7,
2973 const char *name8, const char *name9, const char *name10)
2974{
2975 if (fNpar >= 1)
2976 SetParName(0, name0);
2977 if (fNpar >= 2)
2978 SetParName(1, name1);
2979 if (fNpar >= 3)
2980 SetParName(2, name2);
2981 if (fNpar >= 4)
2982 SetParName(3, name3);
2983 if (fNpar >= 5)
2984 SetParName(4, name4);
2985 if (fNpar >= 6)
2986 SetParName(5, name5);
2987 if (fNpar >= 7)
2988 SetParName(6, name6);
2989 if (fNpar >= 8)
2990 SetParName(7, name7);
2991 if (fNpar >= 9)
2992 SetParName(8, name8);
2993 if (fNpar >= 10)
2994 SetParName(9, name9);
2995 if (fNpar >= 11)
2996 SetParName(10, name10);
2997}
2998
2999////////////////////////////////////////////////////////////////////////////////
3000void TFormula::SetParName(Int_t ipar, const char * name)
3001{
3002
3003 if (ipar < 0 || ipar > fNpar) {
3004 Error("SetParName","Wrong Parameter index %d ",ipar);
3005 return;
3006 }
3007 TString oldName;
3008 // find parameter with given index
3009 for ( auto &it : fParams) {
3010 if (it.second == ipar) {
3011 oldName = it.first;
3012 fParams.erase(oldName);
3013 fParams.insert(std::make_pair(name, ipar) );
3014 break;
3015 }
3016 }
3017 if (oldName.IsNull() ) {
3018 Error("SetParName","Parameter %d is not existing.",ipar);
3019 return;
3020 }
3021
3022 //replace also parameter name in formula expression in case is not a lambda
3024
3025}
3026
3027////////////////////////////////////////////////////////////////////////////////
3028/// Replace in Formula expression the parameter name.
3029
3030void TFormula::ReplaceParamName(TString & formula, const TString & oldName, const TString & name){
3031 if (!formula.IsNull() ) {
3032 bool found = false;
3033 for(list<TFormulaFunction>::iterator it = fFuncs.begin(); it != fFuncs.end(); ++it)
3034 {
3035 if (oldName == it->GetName()) {
3036 found = true;
3037 it->fName = name;
3038 break;
3039 }
3040 }
3041 if (!found) {
3042 Error("SetParName", "Parameter %s is not defined.", oldName.Data());
3043 return;
3044 }
3045 // change whitespace to \s to avoid problems in parsing
3046 TString newName = name;
3047 newName.ReplaceAll(" ", "\\s");
3048 TString pattern = TString::Format("[%s]", oldName.Data());
3049 TString replacement = TString::Format("[%s]", newName.Data());
3050 formula.ReplaceAll(pattern, replacement);
3051 }
3052
3053}
3054
3055////////////////////////////////////////////////////////////////////////////////
3057{
3058#ifdef R__HAS_VECCORE
3059 if (fNdim == 0) {
3060 Info("SetVectorized","Cannot vectorized a function of zero dimension");
3061 return;
3062 }
3063 if (vectorized != fVectorized) {
3064 if (!fFormula)
3065 Error("SetVectorized", "Cannot set vectorized to %d -- Formula is missing", vectorized);
3066
3067 fVectorized = vectorized;
3068 // no need to JIT a new signature in case of zero dimension
3069 //if (fNdim== 0) return;
3070 fClingInitialized = false;
3071 fReadyToExecute = false;
3072 fClingName = "";
3074
3075 fMethod.reset();
3076 // should I add fGradMethod.reset() ?
3077
3078 FillVecFunctionsShurtCuts(); // to replace with the right vectorized signature (e.g. sin -> vecCore::math::Sin)
3081 }
3082#else
3083 if (vectorized)
3084 Warning("SetVectorized", "Cannot set vectorized -- try building with option -Dbuiltin_veccore=On");
3085#endif
3086}
3087
3088////////////////////////////////////////////////////////////////////////////////
3089Double_t TFormula::EvalPar(const Double_t *x,const Double_t *params) const
3090{
3091 if (!fVectorized)
3092 return DoEval(x, params);
3093
3094#ifdef R__HAS_VECCORE
3095
3096 if (fNdim == 0 || !x) {
3097 ROOT::Double_v ret = DoEvalVec(nullptr, params);
3098 return vecCore::Get( ret, 0 );
3099 }
3100
3101 // otherwise, regular Double_t inputs on a vectorized function
3102
3103 // convert our input into vectors then convert back
3104 if (gDebug)
3105 Info("EvalPar", "Function is vectorized - converting Double_t into ROOT::Double_v and back");
3106
3107 if (fNdim < 5) {
3108 const int maxDim = 4;
3109 std::array<ROOT::Double_v, maxDim> xvec;
3110 for (int i = 0; i < fNdim; i++)
3111 xvec[i] = x[i];
3112
3113 ROOT::Double_v ans = DoEvalVec(xvec.data(), params);
3114 return vecCore::Get(ans, 0);
3115 }
3116 // allocating a vector is much slower (we do only for dim > 4)
3117 std::vector<ROOT::Double_v> xvec(fNdim);
3118 for (int i = 0; i < fNdim; i++)
3119 xvec[i] = x[i];
3120
3121 ROOT::Double_v ans = DoEvalVec(xvec.data(), params);
3122 return vecCore::Get(ans, 0);
3123
3124#else
3125 // this should never happen, because fVectorized can only be set true with
3126 // R__HAS_VECCORE, but just in case:
3127 Error("EvalPar", "Formula is vectorized (even though VECCORE is disabled!)");
3128 return TMath::QuietNaN();
3129#endif
3130}
3131
3133
3134static bool functionExists(const string &Name) {
3135 return gInterpreter->GetFunction(/*cl*/0, Name.c_str());
3136}
3137
3138/// returns true on success.
3140{
3141 // We already have generated the gradient.
3142 if (fGradMethod)
3143 return true;
3144
3146 // FIXME: Move this elsewhere
3149 gInterpreter->Declare("#include <Math/CladDerivator.h>\n#pragma clad OFF");
3150 }
3151
3152 // Check if the gradient request was made as part of another TFormula.
3153 // This can happen when we create multiple TFormula objects with the same
3154 // formula. In that case, the hasher will give identical id and we can
3155 // reuse the already generated gradient function.
3157 std::string GradReqFuncName = GetGradientFuncName() + "_req";
3158 // We want to call clad::differentiate(TFormula_id);
3159 fGradGenerationInput = std::string("#pragma cling optimize(2)\n") +
3160 "#pragma clad ON\n" +
3161 "void " + GradReqFuncName + "() {\n" +
3162 "clad::gradient(" + std::string(fClingName.Data()) + ");\n }\n" +
3163 "#pragma clad OFF";
3164
3165 if (!gInterpreter->Declare(fGradGenerationInput.c_str()))
3166 return false;
3167 }
3168
3169 Bool_t hasParameters = (fNpar > 0);
3170 Bool_t hasVariables = (fNdim > 0);
3171 std::string GradFuncName = GetGradientFuncName();
3172 fGradMethod = prepareMethod(hasParameters, hasVariables,
3173 GradFuncName.c_str(),
3174 fVectorized, /*IsGradient*/ true);
3176 return true;
3177 }
3178 return false;
3179}
3180
3182{
3183 if (DoEval(x) == TMath::QuietNaN())
3184 return;
3185
3186 if (!fClingInitialized) {
3187 Error("GradientPar", "Could not initialize the formula!");
3188 return;
3189 }
3190
3191 if (!GenerateGradientPar()) {
3192 Error("GradientPar", "Could not generate a gradient for the formula %s!",
3193 fClingName.Data());
3194 return;
3195 }
3196
3197 if ((int)result.size() < fNpar) {
3198 Warning("GradientPar",
3199 "The size of gradient result is %zu but %d is required. Resizing.",
3200 result.size(), fNpar);
3201 result.resize(fNpar);
3202 }
3203 GradientPar(x, result.data());
3204}
3205
3207{
3208 void* args[3];
3209 const double * vars = (x) ? x : fClingVariables.data();
3210 args[0] = &vars;
3211 if (fNpar <= 0) {
3212 // __attribute__((used)) extern "C" void __cf_0(void* obj, int nargs, void** args, void* ret)
3213 // {
3214 // if (ret) {
3215 // new (ret) (double) (((double (&)(double*))TFormula____id)(*(double**)args[0]));
3216 // return;
3217 // } else {
3218 // ((double (&)(double*))TFormula____id)(*(double**)args[0]);
3219 // return;
3220 // }
3221 // }
3222 args[1] = &result;
3223 (*fGradFuncPtr)(0, 2, args, /*ret*/nullptr); // We do not use ret in a return-void func.
3224 } else {
3225 // __attribute__((used)) extern "C" void __cf_0(void* obj, int nargs, void** args, void* ret)
3226 // {
3227 // ((void (&)(double*, double*,
3228 // double*))TFormula____id_grad)(*(double**)args[0], *(double**)args[1],
3229 // *(double**)args[2]);
3230 // return;
3231 // }
3232 const double *pars = fClingParameters.data();
3233 args[1] = &pars;
3234 args[2] = &result;
3235 (*fGradFuncPtr)(0, 3, args, /*ret*/nullptr); // We do not use ret in a return-void func.
3236 }
3237}
3238
3239////////////////////////////////////////////////////////////////////////////////
3240#ifdef R__HAS_VECCORE
3241// ROOT::Double_v TFormula::Eval(ROOT::Double_v x, ROOT::Double_v y, ROOT::Double_v z, ROOT::Double_v t) const
3242// {
3243// ROOT::Double_v xxx[] = {x, y, z, t};
3244// return EvalPar(xxx, nullptr);
3245// }
3246
3247ROOT::Double_v TFormula::EvalParVec(const ROOT::Double_v *x, const Double_t *params) const
3248{
3249 if (fVectorized)
3250 return DoEvalVec(x, params);
3251
3252 if (fNdim == 0 || !x)
3253 return DoEval(nullptr, params); // automatic conversion to vectorized
3254
3255 // otherwise, trying to input vectors into a scalar function
3256
3257 if (gDebug)
3258 Info("EvalPar", "Function is not vectorized - converting ROOT::Double_v into Double_t and back");
3259
3260 const int vecSize = vecCore::VectorSize<ROOT::Double_v>();
3261 std::vector<Double_t> xscalars(vecSize*fNdim);
3262
3263 for (int i = 0; i < vecSize; i++)
3264 for (int j = 0; j < fNdim; j++)
3265 xscalars[i*fNdim+j] = vecCore::Get(x[j],i);
3266
3267 ROOT::Double_v answers(0.);
3268 for (int i = 0; i < vecSize; i++)
3269 vecCore::Set(answers, i, DoEval(&xscalars[i*fNdim], params));
3270
3271 return answers;
3272}
3273#endif
3274
3275////////////////////////////////////////////////////////////////////////////////
3276/// Sets first 4 variables (e.g. x, y, z, t) and evaluate formula.
3277
3279{
3280 double xxx[4] = {x,y,z,t};
3281 return EvalPar(xxx, nullptr); // takes care of case where formula is vectorized
3282}
3283
3284////////////////////////////////////////////////////////////////////////////////
3285/// Sets first 3 variables (e.g. x, y, z) and evaluate formula.
3286
3288{
3289 double xxx[3] = {x,y,z};
3290 return EvalPar(xxx, nullptr);
3291}
3292
3293////////////////////////////////////////////////////////////////////////////////
3294/// Sets first 2 variables (e.g. x and y) and evaluate formula.
3295
3297{
3298 double xxx[2] = {x,y};
3299 return EvalPar(xxx, nullptr);
3300}
3301
3302////////////////////////////////////////////////////////////////////////////////
3303/// Sets first variable (e.g. x) and evaluate formula.
3304
3306{
3307 double * xxx = &x;
3308 return EvalPar(xxx, nullptr);
3309}
3310
3311////////////////////////////////////////////////////////////////////////////////
3312/// Evaluate formula.
3313/// If formula is not ready to execute(missing parameters/variables),
3314/// print these which are not known.
3315/// If parameter has default value, and has not been set, appropriate warning is shown.
3316
3317Double_t TFormula::DoEval(const double * x, const double * params) const
3318{
3319 if(!fReadyToExecute)
3320 {
3321 Error("Eval", "Formula is invalid and not ready to execute ");
3322 for (auto it = fFuncs.begin(); it != fFuncs.end(); ++it) {
3323 TFormulaFunction fun = *it;
3324 if (!fun.fFound) {
3325 printf("%s is unknown.\n", fun.GetName());
3326 }
3327 }
3328 return TMath::QuietNaN();
3329 }
3330
3331 // Lazy initialization is set and needed when reading from a file
3333 // try recompiling the formula. We need to lock because this is not anymore thread safe
3335 // check again in case another thread has initialized the formula (see ROOT-10994)
3336 if (!fClingInitialized) {
3337 auto thisFormula = const_cast<TFormula*>(this);
3338 thisFormula->ReInitializeEvalMethod();
3339 }
3340 if (!fClingInitialized) {
3341 Error("DoEval", "Formula has error and it is not properly initialized ");
3342 return TMath::QuietNaN();
3343 }
3344 }
3345
3346 if (fLambdaPtr && TestBit(TFormula::kLambda)) {// case of lambda functions
3347 std::function<double(double *, double *)> & fptr = * ( (std::function<double(double *, double *)> *) fLambdaPtr);
3348 assert(x);
3349 //double * v = (x) ? const_cast<double*>(x) : const_cast<double*>(fClingVariables.data());
3350 double * v = const_cast<double*>(x);
3351 double * p = (params) ? const_cast<double*>(params) : const_cast<double*>(fClingParameters.data());
3352 return fptr(v, p);
3353 }
3354
3355
3356 Double_t result = 0;
3357 void* args[2];
3358 double * vars = (x) ? const_cast<double*>(x) : const_cast<double*>(fClingVariables.data());
3359 args[0] = &vars;
3360 if (fNpar <= 0) {
3361 (*fFuncPtr)(0, 1, args, &result);
3362 } else {
3363 double *pars = (params) ? const_cast<double *>(params) : const_cast<double *>(fClingParameters.data());
3364 args[1] = &pars;
3365 (*fFuncPtr)(0, 2, args, &result);
3366 }
3367 return result;
3368}
3369
3370////////////////////////////////////////////////////////////////////////////////
3371// Copied from DoEval, but this is the vectorized version
3372#ifdef R__HAS_VECCORE
3373ROOT::Double_v TFormula::DoEvalVec(const ROOT::Double_v *x, const double *params) const
3374{
3375 if (!fReadyToExecute) {
3376 Error("Eval", "Formula is invalid and not ready to execute ");
3377 for (auto it = fFuncs.begin(); it != fFuncs.end(); ++it) {
3378 TFormulaFunction fun = *it;
3379 if (!fun.fFound) {
3380 printf("%s is unknown.\n", fun.GetName());
3381 }
3382 }
3383 return TMath::QuietNaN();
3384 }
3385 // todo maybe save lambda ptr stuff for later
3386
3388 // try recompiling the formula. We need to lock because this is not anymore thread safe
3390 // check again in case another thread has initialized the formula (see ROOT-10994)
3391 if (!fClingInitialized) {
3392 auto thisFormula = const_cast<TFormula*>(this);
3393 thisFormula->ReInitializeEvalMethod();
3394 }
3395 if (!fClingInitialized) {
3396 Error("DoEval", "Formula has error and it is not properly initialized ");
3398 return res;
3399 }
3400 }
3401
3402 ROOT::Double_v result = 0;
3403 void *args[2];
3404
3405 ROOT::Double_v *vars = const_cast<ROOT::Double_v *>(x);
3406 args[0] = &vars;
3407 if (fNpar <= 0) {
3408 (*fFuncPtr)(0, 1, args, &result);
3409 }else {
3410 double *pars = (params) ? const_cast<double *>(params) : const_cast<double *>(fClingParameters.data());
3411 args[1] = &pars;
3412 (*fFuncPtr)(0, 2, args, &result);
3413 }
3414 return result;
3415}
3416#endif // R__HAS_VECCORE
3417
3418
3419//////////////////////////////////////////////////////////////////////////////
3420/// Re-initialize eval method
3421///
3422/// This function is called by DoEval and DoEvalVector in case of a previous failure
3423/// or in case of reading from a file
3424////////////////////////////////////////////////////////////////////////////////
3426
3427
3428 if (TestBit(TFormula::kLambda) ) {
3429 Info("ReInitializeEvalMethod","compile now lambda expression function using Cling");
3431 fLazyInitialization = false;
3432 return;
3433 }
3434 fMethod.reset();
3435
3436 if (!fLazyInitialization) Warning("ReInitializeEvalMethod", "Formula is NOT properly initialized - try calling again TFormula::PrepareEvalMethod");
3437 //else Info("ReInitializeEvalMethod", "Compile now the formula expression using Cling");
3438
3439 // check first if formula exists in the global map
3440 {
3441
3443
3444 // std::cout << "gClingFunctions list" << std::endl;
3445 // for (auto thing : gClingFunctions)
3446 // std::cout << "gClingFunctions : " << thing.first << std::endl;
3447
3448 auto funcit = gClingFunctions.find(fSavedInputFormula);
3449
3450 if (funcit != gClingFunctions.end()) {
3451 fFuncPtr = (TFormula::CallFuncSignature)funcit->second;
3452 fClingInitialized = true;
3453 fLazyInitialization = false;
3454 return;
3455 }
3456 }
3457 // compile now formula using cling
3459 if (fClingInitialized && !fLazyInitialization) Info("ReInitializeEvalMethod", "Formula is now properly initialized !!");
3460 fLazyInitialization = false;
3461
3462 // add function pointer in global map
3463 if (fClingInitialized) {
3464 R__ASSERT(!fSavedInputFormula.empty());
3465 // if Cling has been successfully initialized
3466 // put function ptr in the static map
3468 gClingFunctions.insert(std::make_pair(fSavedInputFormula, (void *)fFuncPtr));
3469 }
3470
3471
3472 return;
3473}
3474
3475////////////////////////////////////////////////////////////////////////////////
3476/// Return the expression formula.
3477///
3478/// - If option = "P" replace the parameter names with their values
3479/// - If option = "CLING" return the actual expression used to build the function passed to cling
3480/// - If option = "CLINGP" replace in the CLING expression the parameter with their values
3481
3483{
3484 TString opt(option);
3485 if (opt.IsNull() || TestBit(TFormula::kLambda) ) return fFormula;
3486 opt.ToUpper();
3487
3488 // if (opt.Contains("N") ) {
3489 // TString formula = fFormula;
3490 // ReplaceParName(formula, ....)
3491 // }
3492
3493 if (opt.Contains("CLING") ) {
3494 std::string clingFunc = fClingInput.Data();
3495 std::size_t found = clingFunc.find("return");
3496 std::size_t found2 = clingFunc.rfind(";");
3497 if (found == std::string::npos || found2 == std::string::npos) {
3498 Error("GetExpFormula","Invalid Cling expression - return default formula expression");
3499 return fFormula;
3500 }
3501 TString clingFormula = fClingInput(found+7,found2-found-7);
3502 // to be implemented
3503 if (!opt.Contains("P")) return clingFormula;
3504 // replace all "p[" with "[parname"
3505 int i = 0;
3506 while (i < clingFormula.Length()-2 ) {
3507 // look for p[number
3508 if (clingFormula[i] == 'p' && clingFormula[i+1] == '[' && isdigit(clingFormula[i+2]) ) {
3509 int j = i+3;
3510 while ( isdigit(clingFormula[j]) ) { j++;}
3511 if (clingFormula[j] != ']') {
3512 Error("GetExpFormula","Parameters not found - invalid expression - return default cling formula");
3513 return clingFormula;
3514 }
3515 TString parNumbName = clingFormula(i+2,j-i-2);
3516 int parNumber = parNumbName.Atoi();
3517 assert(parNumber < fNpar);
3518 TString replacement = TString::Format("%f",GetParameter(parNumber));
3519 clingFormula.Replace(i,j-i+1, replacement );
3520 i += replacement.Length();
3521 }
3522 i++;
3523 }
3524 return clingFormula;
3525 }
3526 if (opt.Contains("P") ) {
3527 // replace parameter names with their values
3528 TString expFormula = fFormula;
3529 int i = 0;
3530 while (i < expFormula.Length()-2 ) {
3531 // look for [parName]
3532 if (expFormula[i] == '[') {
3533 int j = i+1;
3534 while ( expFormula[j] != ']' ) { j++;}
3535 if (expFormula[j] != ']') {
3536 Error("GetExpFormula","Parameter names not found - invalid expression - return default formula");
3537 return expFormula;
3538 }
3539 TString parName = expFormula(i+1,j-i-1);
3540 TString replacement = TString::Format("%g",GetParameter(parName));
3541 expFormula.Replace(i,j-i+1, replacement );
3542 i += replacement.Length();
3543 }
3544 i++;
3545 }
3546 return expFormula;
3547 }
3548 Warning("GetExpFormula","Invalid option - return default formula expression");
3549 return fFormula;
3550}
3551
3553 std::unique_ptr<TInterpreterValue> v = gInterpreter->MakeInterpreterValue();
3554 gInterpreter->Evaluate(GetGradientFuncName().c_str(), *v);
3555 return v->ToString();
3556}
3557
3558////////////////////////////////////////////////////////////////////////////////
3559/// Print the formula and its attributes.
3560
3561void TFormula::Print(Option_t *option) const
3562{
3563 printf(" %20s : %s Ndim= %d, Npar= %d, Number= %d \n",GetName(),GetTitle(), fNdim,fNpar,fNumber);
3564 printf(" Formula expression: \n");
3565 printf("\t%s \n",fFormula.Data() );
3566 TString opt(option);
3567 opt.ToUpper();
3568 // do an evaluation as a cross-check
3569 //if (fReadyToExecute) Eval();
3570
3571 if (opt.Contains("V") ) {
3572 if (fNdim > 0 && !TestBit(TFormula::kLambda)) {
3573 printf("List of Variables: \n");
3574 assert(int(fClingVariables.size()) >= fNdim);
3575 for ( int ivar = 0; ivar < fNdim ; ++ivar) {
3576 printf("Var%4d %20s = %10f \n",ivar,GetVarName(ivar).Data(), fClingVariables[ivar]);
3577 }
3578 }
3579 if (fNpar > 0) {
3580 printf("List of Parameters: \n");
3581 if ( int(fClingParameters.size()) < fNpar)
3582 Error("Print","Number of stored parameters in vector %zu in map %zu is different than fNpar %d",fClingParameters.size(), fParams.size(), fNpar);
3583 assert(int(fClingParameters.size()) >= fNpar);
3584 // print with order passed to Cling function
3585 for ( int ipar = 0; ipar < fNpar ; ++ipar) {
3586 printf("Par%4d %20s = %10f \n",ipar,GetParName(ipar), fClingParameters[ipar] );
3587 }
3588 }
3589 printf("Expression passed to Cling:\n");
3590 printf("\t%s\n",fClingInput.Data() );
3591 if (fGradFuncPtr) {
3592 printf("Generated Gradient:\n");
3593 printf("%s\n", fGradGenerationInput.c_str());
3594 printf("%s\n", GetGradientFormula().Data());
3595 }
3596 }
3597 if(!fReadyToExecute)
3598 {
3599 Warning("Print", "Formula is not ready to execute. Missing parameters/variables");
3600 for (list<TFormulaFunction>::const_iterator it = fFuncs.begin(); it != fFuncs.end(); ++it) {
3601 TFormulaFunction fun = *it;
3602 if (!fun.fFound) {
3603 printf("%s is unknown.\n", fun.GetName());
3604 }
3605 }
3606 }
3607 if (!fAllParametersSetted) {
3608 // we can skip this
3609 // Info("Print","Not all parameters are set.");
3610 // for(map<TString,TFormulaVariable>::const_iterator it = fParams.begin(); it != fParams.end(); ++it)
3611 // {
3612 // pair<TString,TFormulaVariable> param = *it;
3613 // if(!param.second.fFound)
3614 // {
3615 // printf("%s has default value %lf\n",param.first.Data(),param.second.GetInitialValue());
3616 // }
3617 // }
3618 }
3619}
3620
3621////////////////////////////////////////////////////////////////////////////////
3622/// Stream a class object.
3623
3624void TFormula::Streamer(TBuffer &b)
3625{
3626 if (b.IsReading() ) {
3627 UInt_t R__s, R__c;
3628 Version_t v = b.ReadVersion(&R__s, &R__c);
3629 //std::cout << "version " << v << std::endl;
3630 if (v <= 8 && v > 3 && v != 6) {
3631 // old TFormula class
3633 // read old TFormula class
3634 fold->Streamer(b, v, R__s, R__c, TFormula::Class());
3635 //std::cout << "read old tformula class " << std::endl;
3636 TFormula fnew(fold->GetName(), fold->GetExpFormula() );
3637
3638 *this = fnew;
3639
3640 // printf("copying content in a new TFormula \n");
3641 SetParameters(fold->GetParameters() );
3642 if (!fReadyToExecute ) {
3643 Error("Streamer","Old formula read from file is NOT valid");
3644 Print("v");
3645 }
3646 delete fold;
3647 return;
3648 }
3649 else if (v > 8) {
3650 // new TFormula class
3651 b.ReadClassBuffer(TFormula::Class(), this, v, R__s, R__c);
3652
3653 //std::cout << "reading npar = " << GetNpar() << std::endl;
3654
3655 // initialize the formula
3656 // need to set size of fClingVariables which is transient
3657 //fClingVariables.resize(fNdim);
3658
3659 // case of formula contains only parameters
3660 if (fFormula.IsNull() ) return;
3661
3662
3663 // store parameter values, names and order
3664 std::vector<double> parValues = fClingParameters;
3665 auto paramMap = fParams;
3666 fNpar = fParams.size();
3667
3668 fLazyInitialization = true; // when reading we initialize the formula later to avoid problem of recursive Jitting
3669
3670 if (!TestBit(TFormula::kLambda) ) {
3671
3672 // save dimension read from the file (stored for V >=12)
3673 // and we check after initializing if it is the same
3674 int ndim = fNdim;
3675 fNdim = 0;
3676
3677 //std::cout << "Streamer::Reading preprocess the formula " << fFormula << " ndim = " << fNdim << " npar = " << fNpar << std::endl;
3678 // for ( auto &p : fParams)
3679 // std::cout << "parameter " << p.first << " index " << p.second << std::endl;
3680
3681 fClingParameters.clear(); // need to be reset before re-initializing it
3682
3683 FillDefaults();
3684
3685
3687
3688 //std::cout << "Streamer::after pre-process the formula " << fFormula << " ndim = " << fNdim << " npar = " << fNpar << std::endl;
3689
3691
3692 //std::cout << "Streamer::after prepared " << fClingInput << " ndim = " << fNdim << " npar = " << fNpar << std::endl;
3693
3694
3695 // restore parameter values
3696 if (fNpar != (int) parValues.size() ) {
3697 Error("Streamer","number of parameters computed (%d) is not same as the stored parameters (%d)",fNpar,int(parValues.size()) );
3698 Print("v");
3699 }
3700 if (v > 11 && fNdim != ndim) {
3701 Error("Streamer","number of dimension computed (%d) is not same as the stored value (%d)",fNdim, ndim );
3702 Print("v");
3703 }
3704 }
3705 else {
3706 // we also delay the initializtion of lamda expressions
3707 if (!fLazyInitialization) {
3708 bool ret = InitLambdaExpression(fFormula);
3709 if (ret) {
3710 fClingInitialized = true;
3711 }
3712 }else {
3713 fReadyToExecute = true;
3714 }
3715 }
3716 assert(fNpar == (int) parValues.size() );
3717 std::copy( parValues.begin(), parValues.end(), fClingParameters.begin() );
3718 // restore parameter names and order
3719 if (fParams.size() != paramMap.size() ) {
3720 Warning("Streamer","number of parameters list found (%zu) is not same as the stored one (%zu) - use re-created list",fParams.size(),paramMap.size()) ;
3721 //Print("v");
3722 }
3723 else
3724 //assert(fParams.size() == paramMap.size() );
3725 fParams = paramMap;
3726
3727 // input formula into Cling
3728 // need to replace in cling the name of the pointer of this object
3729 // TString oldClingName = fClingName;
3730 // fClingName.Replace(fClingName.Index("_0x")+1,fClingName.Length(), TString::Format("%p",this) );
3731 // fClingInput.ReplaceAll(oldClingName, fClingName);
3732 // InputFormulaIntoCling();
3733
3734 if (!TestBit(kNotGlobal)) {
3736 gROOT->GetListOfFunctions()->Add(this);
3737 }
3738 if (!fReadyToExecute ) {
3739 Error("Streamer","Formula read from file is NOT ready to execute");
3740 Print("v");
3741 }
3742 //std::cout << "reading 2 npar = " << GetNpar() << std::endl;
3743
3744 return;
3745 }
3746 else {
3747 Error("Streamer","Reading version %d is not supported",v);
3748 return;
3749 }
3750 }
3751 else {
3752 // case of writing
3753 b.WriteClassBuffer(TFormula::Class(), this);
3754 // std::cout << "writing npar = " << GetNpar() << std::endl;
3755 }
3756}
double
#define b(i)
Definition RSha256.hxx:100
#define f(i)
Definition RSha256.hxx:104
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
const Ssiz_t kNPOS
Definition RtypesCore.h:115
int Int_t
Definition RtypesCore.h:45
short Version_t
Definition RtypesCore.h:65
double Double_t
Definition RtypesCore.h:59
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:364
#define R__ASSERT(e)
Definition TError.h:120
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition TError.cxx:187
void GetParameters(TFitEditor::FuncParams_t &pars, TF1 *func)
Stores the parameters of the given function into pars.
static std::unordered_map< std::string, void * > gClingFunctions
Definition TFormula.cxx:213
void(*)(Int_t nobjects, TObject **from, TObject **to) TFormulaUpdater_t
Definition TFormula.cxx:229
static const TString gNamePrefix
Definition TFormula.cxx:209
static bool IsReservedName(const char *name)
Definition TFormula.cxx:404
bool R__SetClonesArrayTFormulaUpdater(TFormulaUpdater_t func)
static std::unique_ptr< TMethodCall > prepareMethod(bool HasParameters, bool HasVariables, const char *FuncName, bool IsVectorized, bool IsGradient=false)
Definition TFormula.cxx:749
static bool functionExists(const string &Name)
int R__RegisterTFormulaUpdaterTrigger
Definition TFormula.cxx:232
static TInterpreter::CallFuncIFacePtr_t::Generic_t prepareFuncPtr(TMethodCall *method)
Definition TFormula.cxx:784
static void R__v5TFormulaUpdater(Int_t nobjects, TObject **from, TObject **to)
Definition TFormula.cxx:215
char name[80]
Definition TGX11.cxx:110
R__EXTERN TInterpreter * gCling
#define gInterpreter
Int_t gDebug
Definition TROOT.cxx:590
R__EXTERN TVirtualMutex * gROOTMutex
Definition TROOT.h:63
#define gROOT
Definition TROOT.h:406
typedef void((*Func_t)())
#define R__LOCKGUARD(mutex)
The FORMULA class (ROOT version 5)
Definition TFormula.h:65
virtual Double_t * GetParameters() const
Definition TFormula.h:243
virtual TString GetExpFormula(Option_t *option="") const
Reconstruct the formula expression from the internal TFormula member variables.
void Streamer(TBuffer &b, const TClass *onfile_class)
Stream a class object.
Buffer base class used for serializing objects.
Definition TBuffer.h:43
TClass instances represent classes, structs and namespaces in the ROOT type system.
Definition TClass.h:80
const TList * GetListOfAllPublicMethods(Bool_t load=kTRUE)
Returns a list of all public methods of this class and its base classes.
Definition TClass.cxx:3822
static TClass * GetClass(const char *name, Bool_t load=kTRUE, Bool_t silent=kFALSE)
Static method returning pointer to TClass of the specified class name.
Definition TClass.cxx:2957
1-Dim function class
Definition TF1.h:213
virtual TFormula * GetFormula()
Definition TF1.h:453
Helper class for TFormula.
Definition TFormula.h:30
const char * GetName() const
Definition TFormula.h:37
Int_t GetNargs() const
Definition TFormula.h:39
TString fName
Definition TFormula.h:32
Bool_t IsFuncCall() const
Definition TFormula.h:40
Another helper class for TFormula.
Definition TFormula.h:63
TString fName
Definition TFormula.h:65
Double_t fValue
Definition TFormula.h:66
The Formula class.
Definition TFormula.h:87
Double_t GetParameter(const char *name) const
Returns parameter value given by string.
Bool_t PrepareFormula(TString &formula)
prepare the formula to be executed normally is called with fFormula
std::atomic< Bool_t > fClingInitialized
trasient to force initialization
Definition TFormula.h:96
void FillDefaults()
Fill structures with default variables, constants and function shortcuts.
Definition TFormula.cxx:849
const TObject * GetLinearPart(Int_t i) const
Return linear part.
void InputFormulaIntoCling()
Inputs formula, transfered to C++ code into Cling.
Definition TFormula.cxx:823
virtual ~TFormula()
Definition TFormula.cxx:413
void DoAddParameter(const TString &name, Double_t value, bool processFormula)
Adds parameter to known parameters.
void HandleLinear(TString &formula)
Handle linear functions defined with the operator ++.
TString fClingName
pointer to a methodcall
Definition TFormula.h:101
TString GetVarName(Int_t ivar) const
Returns variable name given its position in the array.
void SetVariable(const TString &name, Double_t value)
Sets variable value.
std::vector< Double_t > GradientStorage
Definition TFormula.h:177
Int_t GetParNumber(const char *name) const
Return parameter index given a name (return -1 for not existing parameters) non need to print an erro...
void DoSetParameters(const Double_t *p, Int_t size)
TString GetGradientFormula() const
Double_t Eval(Double_t x) const
Sets first variable (e.g. x) and evaluate formula.
void HandleParametrizedFunctions(TString &formula)
Handling parametrized functions Function can be normalized, and have different variable then x.
TInterpreter::CallFuncIFacePtr_t::Generic_t CallFuncSignature
unique name used to defined the function and used in the global map (need to be saved in case of lazy...
Definition TFormula.h:104
Double_t * GetParameters() const
void SetParName(Int_t ipar, const char *name)
std::string fGradGenerationInput
Definition TFormula.h:105
Double_t EvalPar(const Double_t *x, const Double_t *params=0) const
static Bool_t IsAParameterName(const TString &formula, int ipos)
Definition TFormula.cxx:302
bool HasGradientGenerationFailed() const
Definition TFormula.h:130
void ReplaceParamName(TString &formula, const TString &oldname, const TString &name)
Replace in Formula expression the parameter name.
void SetPredefinedParamNames()
Set parameter names only in case of pre-defined functions.
std::map< TString, Double_t > fConsts
Definition TFormula.h:139
std::map< TString, TString > fFunctionsShortcuts
Definition TFormula.h:140
void HandleParamRanges(TString &formula)
Handling parameter ranges, in the form of [1..5].
std::vector< TObject * > fLinearParts
Definition TFormula.h:145
static Bool_t IsBracket(const char c)
Definition TFormula.cxx:243
Bool_t fVectorized
Definition TFormula.h:146
TString fClingInput
Definition TFormula.h:92
Bool_t PrepareEvalMethod()
Sets TMethodCall to function inside Cling environment.
Definition TFormula.cxx:808
Int_t fNumber
Number of parameter (transient since we save the vector)
Definition TFormula.h:144
std::string GetGradientFuncName() const
Definition TFormula.h:126
static bool fIsCladRuntimeIncluded
pointer to the lambda function
Definition TFormula.h:109
Double_t GetVariable(const char *name) const
Returns variable value.
std::list< TFormulaFunction > fFuncs
Definition TFormula.h:136
Bool_t fAllParametersSetted
transient to force re-initialization
Definition TFormula.h:97
void ProcessFormula(TString &formula)
Iterates through functors in fFuncs and performs the appropriate action.
static Bool_t IsOperator(const char c)
Definition TFormula.cxx:235
void SetVectorized(Bool_t vectorized)
void FillVecFunctionsShurtCuts()
Fill the shortcuts for vectorized functions We will replace for example sin with vecCore::Mat::Sin.
Definition TFormula.cxx:917
const char * GetParName(Int_t ipar) const
Return parameter name given by integer.
std::map< TString, TFormulaVariable > fVars
Definition TFormula.h:137
CallFuncSignature fFuncPtr
input query to clad to generate a gradient
Definition TFormula.h:106
void HandlePolN(TString &formula)
Handling polN If before 'pol' exist any name, this name will be treated as variable used in polynomia...
Definition TFormula.cxx:949
static Bool_t IsHexadecimal(const TString &formula, int ipos)
Definition TFormula.cxx:279
void ExtractFunctors(TString &formula)
Extracts functors from formula, and put them in fFuncs.
Bool_t InitLambdaExpression(const char *formula)
Definition TFormula.cxx:554
void SetParameters(const Double_t *params)
Set a vector of parameters value.
void Print(Option_t *option="") const
Print the formula and its attributes.
std::string fSavedInputFormula
unique name passed to Cling to define the function ( double clingName(double*x, double*p) )
Definition TFormula.h:102
void GradientPar(const Double_t *x, TFormula::GradientStorage &result)
Compute the gradient employing automatic differentiation.
void SetName(const char *name)
Set the name of the formula.
Bool_t IsValid() const
Definition TFormula.h:241
void ReplaceAllNames(TString &formula, std::map< TString, TString > &substitutions)
Definition TFormula.cxx:354
static Bool_t IsDefaultVariableName(const TString &name)
Definition TFormula.cxx:261
void FillParametrizedFunctions(std::map< std::pair< TString, Int_t >, std::pair< TString, TString > > &functions)
Fill map with parametrized functions.
void SetParNames(const char *name0="p0", const char *name1="p1", const char *name2="p2", const char *name3="p3", const char *name4="p4", const char *name5="p5", const char *name6="p6", const char *name7="p7", const char *name8="p8", const char *name9="p9", const char *name10="p10")
void AddVariables(const TString *vars, const Int_t size)
Adds multiple variables.
Int_t GetVarNumber(const char *name) const
Returns variable number (positon in array) given its name.
std::vector< Double_t > fClingVariables
input function passed to Cling
Definition TFormula.h:93
std::unique_ptr< TMethodCall > fMethod
transient flag to control lazy initialization (needed for reading from files)
Definition TFormula.h:99
TString fFormula
Definition TFormula.h:141
static Bool_t IsScientificNotation(const TString &formula, int ipos)
Definition TFormula.cxx:267
CallFuncSignature fGradFuncPtr
function pointer, owned by the JIT.
Definition TFormula.h:107
std::vector< Double_t > fClingParameters
cached variables
Definition TFormula.h:94
void SetParameter(const char *name, Double_t value)
Sets parameter value.
std::unique_ptr< TMethodCall > fGradMethod
pointer to methodcall
Definition TFormula.h:100
void PreProcessFormula(TString &formula)
Preprocessing of formula Replace all ** by ^, and removes spaces.
TString GetExpFormula(Option_t *option="") const
Return the expression formula.
Int_t fNdim
Definition TFormula.h:142
void SetVariables(const std::pair< TString, Double_t > *vars, const Int_t size)
Sets multiple variables.
void ReInitializeEvalMethod()
Re-initialize eval method.
@ kNotGlobal
Definition TFormula.h:172
@ kNormalized
Definition TFormula.h:173
void * fLambdaPtr
function pointer, owned by the JIT.
Definition TFormula.h:108
TFormula & operator=(const TFormula &rhs)
= operator.
Definition TFormula.cxx:544
Int_t Compile(const char *expression="")
Compile the given expression with Cling backward compatibility method to be used in combination with ...
Definition TFormula.cxx:599
Int_t fNpar
Definition TFormula.h:143
void HandleFunctionArguments(TString &formula)
Handling user functions (and parametrized functions) to take variables and optionally parameters as a...
std::map< TString, Int_t, TFormulaParamOrder > fParams
list of variable names
Definition TFormula.h:138
void AddVariable(const TString &name, Double_t value=0)
Adds variable to known variables, and reprocess formula.
Bool_t fLazyInitialization
Definition TFormula.h:98
virtual void Clear(Option_t *option="")
Clear the formula setting expression to empty and reset the variables and parameters containers.
Definition TFormula.cxx:716
void HandleExponentiation(TString &formula)
Handling exponentiation Can handle multiple carets, eg.
static Bool_t IsFunctionNameChar(const char c)
Definition TFormula.cxx:255
virtual void Copy(TObject &f1) const
Copy this to obj.
Definition TFormula.cxx:633
Double_t DoEval(const Double_t *x, const Double_t *p=nullptr) const
Evaluate formula.
Bool_t fReadyToExecute
Definition TFormula.h:95
bool GenerateGradientPar()
Generate gradient computation routine with respect to the parameters.
Global functions class (global functions are obtained from CINT).
Definition TFunction.h:30
Int_t GetNargsOpt() const
Number of function optional (default) arguments.
Int_t GetNargs() const
Number of function arguments.
virtual Bool_t Declare(const char *code)=0
virtual Bool_t CallFunc_IsValid(CallFunc_t *) const
virtual Long_t ProcessLine(const char *line, EErrorCode *error=0)=0
virtual CallFuncIFacePtr_t CallFunc_IFacePtr(CallFunc_t *) const
A doubly linked list.
Definition TList.h:44
Method or function calling interface.
Definition TMethodCall.h:37
CallFunc_t * GetCallFunc() const
Definition TMethodCall.h:93
Each ROOT class (see TClass) has a linked list of methods.
Definition TMethod.h:38
The TNamed class is the base class for all named ROOT classes.
Definition TNamed.h:29
virtual void Copy(TObject &named) const
Copy this to obj.
Definition TNamed.cxx:94
TString fTitle
Definition TNamed.h:33
TString fName
Definition TNamed.h:32
virtual void SetName(const char *name)
Set the name of the TNamed.
Definition TNamed.cxx:140
virtual const char * GetTitle() const
Returns title of object.
Definition TNamed.h:48
virtual const char * GetName() const
Returns name of object.
Definition TNamed.h:47
Mother of all ROOT objects.
Definition TObject.h:37
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition TObject.h:187
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:879
virtual TObject * FindObject(const char *name) const
Must be redefined in derived classes.
Definition TObject.cxx:323
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:696
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:893
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition TObject.cxx:867
Regular expression class.
Definition TRegexp.h:31
Ssiz_t Index(const TString &str, Ssiz_t *len, Ssiz_t start=0) const
Find the first occurrence of the regexp in string and return the position, or -1 if there is no match...
Definition TRegexp.cxx:209
Basic string class.
Definition TString.h:136
Ssiz_t Length() const
Definition TString.h:410
TString & Insert(Ssiz_t pos, const char *s)
Definition TString.h:649
Int_t Atoi() const
Return integer value of string.
Definition TString.cxx:1941
TString & Replace(Ssiz_t pos, Ssiz_t n, const char *s)
Definition TString.h:682
const char * Data() const
Definition TString.h:369
Bool_t IsDigit() const
Returns true if all characters in string are digits (0-9) or white spaces, i.e.
Definition TString.cxx:1783
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition TString.h:692
Ssiz_t Last(char c) const
Find last occurrence of a character c.
Definition TString.cxx:912
void ToUpper()
Change string to upper case.
Definition TString.cxx:1158
Bool_t IsNull() const
Definition TString.h:407
Int_t CountChar(Int_t c) const
Return number of times character c occurs in the string.
Definition TString.cxx:496
TString & Append(const char *cs)
Definition TString.h:564
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString.
Definition TString.cxx:2331
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition TString.h:624
Ssiz_t Index(const char *pat, Ssiz_t i=0, ECaseCompare cmp=kExact) const
Definition TString.h:639
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
TF1 * f1
Definition legend1.C:11
Double_t Double_v
Definition Types.h:51
TROOT * GetROOT()
Definition TROOT.cxx:465
constexpr Double_t G()
Gravitational constant in: .
Definition TMath.h:138
constexpr Double_t C()
Velocity of light in .
Definition TMath.h:117
Short_t Max(Short_t a, Short_t b)
Definition TMathBase.h:212
Double_t QuietNaN()
Returns a quiet NaN as defined by IEEE 754
Definition TMath.h:901
Double_t Floor(Double_t x)
Definition TMath.h:703
constexpr Double_t K()
Boltzmann's constant in .
Definition TMath.h:253
constexpr Double_t Sqrt2()
Definition TMath.h:88
constexpr Double_t E()
Base of natural log:
Definition TMath.h:96
constexpr Double_t Sigma()
Stefan-Boltzmann constant in .
Definition TMath.h:277
constexpr Double_t H()
Planck's constant in .
Definition TMath.h:192
constexpr Double_t LogE()
Base-10 log of e (to convert ln to log)
Definition TMath.h:110
constexpr Double_t Ln10()
Natural log of 10 (to convert log to ln)
Definition TMath.h:103
constexpr Double_t EulerGamma()
Euler-Mascheroni Constant.
Definition TMath.h:339
constexpr Double_t Pi()
Definition TMath.h:37
constexpr Double_t R()
Universal gas constant ( ) in
Definition TMath.h:309
Double_t Log10(Double_t x)
Definition TMath.h:764
Double_t Infinity()
Returns an infinity as defined by the IEEE standard.
Definition TMath.h:914
Definition first.py:1
bool operator()(const TString &a, const TString &b) const
Definition TFormula.cxx:324
void(* Generic_t)(void *, int, void **, void *)
auto * m
Definition textangle.C:8