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