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