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