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