Back to home page

EIC code displayed by LXR

 
 

    


Warning, file /include/root/TFormula.h was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 // @(#)root/hist:$Id$
0002 // Author: Maciej Zimnoch   30/09/2013
0003 
0004 /*************************************************************************
0005  * Copyright (C) 1995-2013, Rene Brun and Fons Rademakers.               *
0006  * All rights reserved.                                                  *
0007  *                                                                       *
0008  * For the licensing terms see $ROOTSYS/LICENSE.                         *
0009  * For the list of contributors see $ROOTSYS/README/CREDITS.             *
0010  *************************************************************************/
0011 #ifndef ROOT_TFormula
0012 #define ROOT_TFormula
0013 
0014 
0015 #include "TNamed.h"
0016 #include "TBits.h"
0017 #include "TInterpreter.h"
0018 #include "TMath.h"
0019 #include <Math/Types.h>
0020 
0021 #include <atomic>
0022 #include <cassert>
0023 #include <list>
0024 #include <map>
0025 #include <string>
0026 #include <vector>
0027 
0028 class TMethodCall;
0029 
0030 
0031 class TFormulaFunction
0032 {
0033 public:
0034    TString  fName;
0035    TString  fBody;
0036    Int_t    fNargs;
0037    Bool_t   fFound;
0038    Bool_t   fFuncCall;
0039    const char *  GetName() const    { return fName.Data(); }
0040    const char *  GetBody() const    { return fBody.Data(); }
0041    Int_t    GetNargs() const   { return fNargs;}
0042    Bool_t   IsFuncCall() const { return fFuncCall;}
0043    TFormulaFunction(){}
0044    TFormulaFunction(const TString &name, const TString &body, int numArgs)
0045       : fName(name),fBody(body),fNargs(numArgs),fFound(false),fFuncCall(true) {}
0046    TFormulaFunction(const TString& name)
0047    : fName(name),fBody(""),fNargs(0),fFound(false),fFuncCall(false){}
0048    Bool_t operator<(const TFormulaFunction &rhv) const
0049    {
0050       // order by length - first the longer ones to avoid replacing wrong functions
0051       if ( fName.Length() < rhv.fName.Length() )
0052          return true;
0053       else if ( fName.Length() > rhv.fName.Length() )
0054          return false;
0055       // case of equal length
0056       return fName < rhv.fName && fBody < rhv.fBody;
0057    }
0058    Bool_t operator==(const TFormulaFunction &rhv) const
0059    {
0060       return fName == rhv.fName && fBody == rhv.fBody && fNargs == rhv.fNargs;
0061    }
0062 };
0063 
0064 class TFormulaVariable
0065 {
0066 public:
0067    TString fName;
0068    Double_t fValue;
0069    Int_t fArrayPos;
0070    Bool_t fFound;
0071    const char * GetName() const     { return fName.Data(); }
0072    Double_t GetInitialValue() const    { return fValue; }
0073    Int_t    GetArrayPos() const { return fArrayPos; }
0074    TFormulaVariable():fName(""),fValue(-1),fArrayPos(-1),fFound(false){}
0075    TFormulaVariable(const TString &name, Double_t value, Int_t pos)
0076    : fName(name), fValue(value), fArrayPos(pos),fFound(false) {}
0077    Bool_t operator<(const TFormulaVariable &rhv) const
0078    {
0079       return fName < rhv.fName;
0080    }
0081 };
0082 
0083 struct TFormulaParamOrder {
0084    bool operator() (const TString& a, const TString& b) const;
0085 };
0086 
0087 
0088 class TFormula : public TNamed
0089 {
0090 private:
0091 
0092 // All data members are transient apart from the string defining the formula and the parameter values
0093    TString           fClingInput;                  ///<! Input function passed to Cling
0094    std::vector<Double_t>  fClingVariables;         ///<! Cached variables
0095    std::vector<Double_t>  fClingParameters;        ///<  Parameter values
0096    Bool_t            fReadyToExecute;              ///<! Transient to force initialization
0097    std::atomic<Bool_t>  fClingInitialized;         ///<! Transient to force re-initialization
0098    Bool_t            fAllParametersSetted;         ///<  Flag to control if all parameters are setted
0099    Bool_t            fLazyInitialization = kFALSE; ///<! Transient flag to control lazy initialization (needed for reading from files)
0100    std::unique_ptr<TMethodCall> fMethod;           ///<! Pointer to methodcall
0101    TString           fClingName;                   ///<! Unique name passed to Cling to define the function ( double clingName(double*x, double*p) )
0102    std::string       fSavedInputFormula;           ///<! Unique name used to defined the function and used in the global map (need to be saved in case of lazy initialization)
0103 
0104    using CallFuncSignature = TInterpreter::CallFuncIFacePtr_t::Generic_t;
0105    std::string       fGradGenerationInput;         ///<! Input query to clad to generate a gradient
0106    std::string       fHessGenerationInput;         ///<! Input query to clad to generate a hessian
0107    CallFuncSignature fFuncPtr = nullptr;           ///<! Function pointer, owned by the JIT.
0108    CallFuncSignature fGradFuncPtr = nullptr;       ///<! Function pointer, owned by the JIT.
0109    CallFuncSignature fHessFuncPtr = nullptr;       ///<! Function pointer, owned by the JIT.
0110    void *   fLambdaPtr = nullptr;                  ///<! Pointer to the lambda function
0111    static bool       fIsCladRuntimeIncluded;
0112 
0113    void     InputFormulaIntoCling();
0114    Bool_t   PrepareEvalMethod();
0115    void     FillDefaults();
0116    void     HandlePolN(TString &formula);
0117    void     HandleParametrizedFunctions(TString &formula);
0118    void HandleParamRanges(TString &formula);
0119    void HandleFunctionArguments(TString &formula);
0120    void     HandleExponentiation(TString &formula);
0121    void     HandleLinear(TString &formula);
0122    Bool_t   InitLambdaExpression(const char * formula);
0123    static Bool_t   IsDefaultVariableName(const TString &name);
0124    void ReplaceAllNames(TString &formula, std::map<TString, TString> &substitutions);
0125    void FillParametrizedFunctions(std::map<std::pair<TString, Int_t>, std::pair<TString, TString>> &functions);
0126    void FillVecFunctionsShurtCuts();
0127    void ReInitializeEvalMethod();
0128    std::string GetGradientFuncName() const {
0129       return std::string(GetUniqueFuncName().Data()) + "_grad_1";
0130    }
0131    std::string GetHessianFuncName() const {
0132       return std::string(GetUniqueFuncName().Data()) + "_hessian_1";
0133    }
0134    bool HasGradientGenerationFailed() const {
0135       return !fGradFuncPtr && !fGradGenerationInput.empty();
0136    }
0137    bool HasHessianGenerationFailed() const {
0138       return !fHessFuncPtr && !fHessGenerationInput.empty();
0139    }
0140 
0141 protected:
0142 
0143    std::list<TFormulaFunction>         fFuncs;              ///<!
0144    std::map<TString,TFormulaVariable>  fVars;               ///<!  List of  variable names
0145    std::map<TString,Int_t,TFormulaParamOrder>   fParams;    ///<|| List of  parameter names
0146    std::map<TString,Double_t>          fConsts;             ///<!
0147    std::map<TString,TString>           fFunctionsShortcuts; ///<!
0148    TString                             fFormula;            ///<   String representing the formula expression
0149    Int_t                               fNdim;               ///<   Dimension - needed for lambda expressions
0150    Int_t                               fNpar;               ///<!  Number of parameter (transient since we save the vector)
0151    Int_t                               fNumber;             ///<   Number used to identify pre-defined functions (gaus, expo,..)
0152    std::vector<TObject*>               fLinearParts;        ///<   Vector of linear functions
0153    Bool_t                              fVectorized = false; ///<   Whether we should use vectorized or regular variables
0154    // (we default to false since a lot of functions still cannot be expressed in vectorized form)
0155 
0156    static Bool_t IsOperator(const char c);
0157    static Bool_t IsBracket(const char c);
0158    static Bool_t IsFunctionNameChar(const char c);
0159    static Bool_t IsHexadecimal(const TString & formula, int ipos);
0160    static Bool_t IsAParameterName(const TString & formula, int ipos);
0161    void   ExtractFunctors(TString &formula);
0162    void   PreProcessFormula(TString &formula);
0163    void   ProcessFormula(TString &formula);
0164    Bool_t PrepareFormula(TString &formula);
0165    void   ReplaceParamName(TString &formula, const TString & oldname, const TString & name);
0166    void   DoAddParameter(const TString &name, Double_t value, bool processFormula);
0167    void   DoSetParameters(const Double_t * p, Int_t size);
0168    void   SetPredefinedParamNames();
0169 
0170    Double_t       DoEval(const Double_t * x, const Double_t * p = nullptr) const;
0171 #ifdef R__HAS_VECCORE
0172    ROOT::Double_v DoEvalVec(const ROOT::Double_v *x, const Double_t *p = nullptr) const;
0173 #endif
0174 
0175 public:
0176 
0177    enum EStatusBits {
0178       kNotGlobal     = BIT(10),    ///< Don't store in gROOT->GetListOfFunction (it should be protected)
0179       kNormalized    = BIT(14),    ///< Set to true if the TFormula (ex gausn) is normalized
0180       kLinear        = BIT(16),    ///< Set to true if the TFormula is for linear fitting
0181       kLambda        = BIT(17)     ///< Set to true if TFormula has been build with a lambda
0182    };
0183    using CladStorage = std::vector<Double_t>;
0184 
0185                   TFormula();
0186           ~TFormula() override;
0187    TFormula&      operator=(const TFormula &rhs);
0188    TFormula(const char *name, const char * formula = "", bool addToGlobList = true, bool vectorize = false);
0189    TFormula(const char *name, const char * formula, int ndim, int npar, bool addToGlobList = true);
0190                   TFormula(const TFormula &formula);
0191    //               TFormula(const char *name, Int_t nparams, Int_t ndims);
0192 
0193    void           AddParameter(const TString &name, Double_t value = 0) { DoAddParameter(name,value,true); }
0194    void           AddVariable(const TString &name, Double_t value = 0);
0195    void           AddVariables(const TString *vars, const Int_t size);
0196    Int_t          Compile(const char *expression="");
0197    void   Copy(TObject &f1) const override;
0198    void   Clear(Option_t * option="") override;
0199    template <typename... Args>
0200    Double_t       Eval(Args... args) const;
0201    Double_t       EvalPar(const Double_t *x, const Double_t *params = nullptr) const;
0202 
0203    /// Generate gradient computation routine with respect to the parameters.
0204    /// \returns true if a gradient was generated and GradientPar can be called.
0205    bool GenerateGradientPar();
0206 
0207    /// Generate hessian computation routine with respect to the parameters.
0208    /// \returns true if a hessian was generated and HessianPar can be called.
0209    bool GenerateHessianPar();
0210 
0211    /// Compute the gradient employing automatic differentiation.
0212    ///
0213    /// \param[in] x - The given variables, if nullptr the already stored
0214    ///                variables are used.
0215    /// \param[out] result - The result of the computation wrt each direction.
0216    void GradientPar(const Double_t *x, TFormula::CladStorage& result);
0217 
0218    void GradientPar(const Double_t *x, Double_t *result);
0219 
0220    /// Compute the gradient employing automatic differentiation.
0221    ///
0222    /// \param[in] x - The given variables, if nullptr the already stored
0223    ///                variables are used.
0224    /// \param[out] result - The 2D hessian matrix flattened to form a vector
0225    ///                      in row-major order.
0226    void HessianPar(const Double_t *x, TFormula::CladStorage& result);
0227 
0228    void HessianPar(const Double_t *x, Double_t *result);
0229 
0230    // query if TFormula provides gradient computation using AD (CLAD)
0231    bool HasGeneratedGradient() const {
0232       return fGradFuncPtr != nullptr;
0233    }
0234 
0235    // query if TFormula provides hessian computation using AD (CLAD)
0236    bool HasGeneratedHessian() const {
0237       return fHessFuncPtr != nullptr;
0238    }
0239 
0240    static Bool_t IsScientificNotation(const TString & formula, int ipos);
0241 
0242    // template <class T>
0243    // T Eval(T x, T y = 0, T z = 0, T t = 0) const;
0244    template <class T>
0245    T EvalPar(const T *x, const Double_t *params = nullptr) const {
0246       return  EvalParVec(x, params);
0247    }
0248 #ifdef R__HAS_VECCORE
0249    ROOT::Double_v EvalParVec(const ROOT::Double_v *x, const Double_t *params = nullptr) const;
0250 #endif
0251    TString        GetExpFormula(Option_t *option = "") const;
0252    TString        GetGradientFormula() const;
0253    TString        GetHessianFormula() const;
0254    TString        GetUniqueFuncName() const {
0255       assert(fClingName.Length() && "TFormula is not initialized yet!");
0256       return fClingName;
0257    }
0258 
0259    const TObject *GetLinearPart(Int_t i) const;
0260    Int_t          GetNdim() const {return fNdim;}
0261    Int_t          GetNpar() const {return fNpar;}
0262    Int_t          GetNumber() const { return fNumber; }
0263    const char *   GetParName(Int_t ipar) const;
0264    Int_t          GetParNumber(const char * name) const;
0265    Double_t       GetParameter(const char * name) const;
0266    Double_t       GetParameter(Int_t param) const;
0267    Double_t*      GetParameters() const;
0268    void           GetParameters(Double_t *params) const;
0269    Double_t       GetVariable(const char *name) const;
0270    Int_t          GetVarNumber(const char *name) const;
0271    TString        GetVarName(Int_t ivar) const;
0272    Bool_t         IsValid() const { return fReadyToExecute && fClingInitialized; }
0273    Bool_t IsVectorized() const { return fVectorized; }
0274    Bool_t         IsLinear() const { return TestBit(kLinear); }
0275    void           Print(Option_t *option = "") const override;
0276    void           SetName(const char* name) override;
0277    void           SetParameter(const char* name, Double_t value);
0278    void           SetParameter(Int_t param, Double_t value);
0279    void           SetParameters(const Double_t *params);
0280    //void           SetParameters(const pair<TString,Double_t> *params, const Int_t size);
0281    template <typename... Args>
0282    void           SetParameters(Double_t arg1, Args &&... args);
0283    void           SetParName(Int_t ipar, const char *name);
0284    template <typename... Args>
0285    void           SetParNames(Args &&... args);
0286    void           SetVariable(const TString &name, Double_t value);
0287    void           SetVariables(const std::pair<TString,Double_t> *vars, const Int_t size);
0288    void SetVectorized(Bool_t vectorized);
0289 
0290    ClassDefOverride(TFormula,14)
0291 };
0292 
0293 ////////////////////////////////////////////////////////////////////////////////
0294 /// Set a list of parameters.
0295 /// The order is by default the alphabetic order given to the parameters,
0296 /// apart if the users has defined explicitly the parameter names.
0297 /// NaN values will be skipped, meaning that the corresponding parameters will not be changed.
0298 
0299 template <typename... Args>
0300 void TFormula::SetParameters(Double_t arg1, Args &&...args)
0301 {
0302    int i = 0;
0303    for (double val : {arg1, static_cast<Double_t>(args)...}) {
0304       if(!TMath::IsNaN(val)) SetParameter(i++, val);
0305    }
0306 }
0307 
0308 ////////////////////////////////////////////////////////////////////////////////
0309 /// Set parameter names.
0310 /// Empty strings will be skipped, meaning that the corresponding name will not be changed.
0311 template <typename... Args>
0312 void TFormula::SetParNames(Args &&...args)
0313 {
0314    int i = 0;
0315    for (auto name : {static_cast<std::string const&>(args)...}) {
0316       if(!name.empty()) SetParName(i++, name.c_str());
0317    }
0318 }
0319 
0320 ////////////////////////////////////////////////////////////////////////////////
0321 /// Set first 1, 2, 3 or 4 variables (e.g. x, y, z and t)
0322 /// and evaluate formula.
0323 
0324 template <typename... Args>
0325 Double_t TFormula::Eval(Args... args) const
0326 {
0327    if (sizeof...(args) > 4) {
0328       Error("Eval", "Eval() only support setting up to 4 variables");
0329    }
0330    double xxx[] = {static_cast<Double_t>(args)...};
0331    return EvalPar(xxx, nullptr);
0332 }
0333 
0334 #endif