ROOT Mathematical Libraries

1. ROOT Mathematical Libraries Lorenzo Moneta CERN/PH-SFT

2. Root Math Work Package • Main responsibilities for this work package: • Basic mathematical functions • Numerical algorithms • Random numbers • Linear algebra • Physics and geometry vectors (3D and 4D) • Fitting and minimization • Histograms (math and statistical part) • Statistics (confidence levels, multivariate analysis) • People contributing: • A. Kreshuk, L. Moneta (CERN) • E. Offermann • J. Leydold, Vienna (Unuran) • M. Fischler, J. Marraffino, W. Brown, FNAL • others: W. Verkerke (RooFit), TMVA team, etc..

3. Outline • New core math library • merge Math in libCore with current libMathCore • numerical algorithms • functor classes (changes in TF1) • Recent developments of ROOT Math Libraries : • Random numbers improvements • Physics and vector package (GenVector) • SMatrix package • MathMore • Fitting and Minimization • plans for new fitting classes • fitting GUI (new fit panel) • Other recent developments: • Histogram comparison (Chi2 test) • FFT , Unuran, TMVA • Future plans

4. New Structure of ROOT Math Libs not yet released with dependency no dependency

5. New Core Math Library • Restructuring of ROOT libraries: new core Math lib with • Math classes from libCore: • TRandom classes, TComplex , TMath (excluding TMathBase) • Classes from current version of MathCore • basic mathematical and statistical functions • remove duplications with TMath • physics vector ( 3D and LorentzVector + Rotation and Boost classes ) • Implementation of basic numerical algorithms • from TF1 : derivation, numerical integration, Brent minimization • Interfaces for functions and algorithms • needed for using different implementations (with plug-in manager) • Functor classes • use of C++ callable objects in TF1 and in the algorithms

6. Classes/Functions size of Library (KB) size of Library and Dictionary (KB) TMath 109 240 TRandom, 1,2,3 55 150 TComplex 4 70 ROOT::Math functions 16 150 Physics Vector 116 (~2000) Numerical algorithms 100 300 Total for libMath 400 900 (~3MB) Library size • Estimation for the size of the library (on Linux) • Dictionary for physics vectors will be in a • separate library ? • move them in libPhysics with TLorentzVector ? ~ 1Mb

7. Math Dictionaries • Template classes (like GenVector and SMatrix) the dictionary is provided for the most used types • double, float and Double32_t • all types of 3D and Lorentz Vectors (12 in total) • for all the main vector operations • for SVector and SMatrix classes up to dimension 7 • dictionary is the dominant part of the library : • 2 Mb on Linux of a 2.3 Mb library for current MathCore • libSMatrix ( ~ 3 Mb) contains only the dictionary. • do not use if you don’t need I/O or to run interactively

8. MathMore Library • C++ interface to GSL algorithms and functions • requires a GSL (version >= 1.8) already installed • Numerical algorithms for 1D functions: • Numerical Derivation • central evaluation (5 points rule) and forward/backward • Numerical Integration • adaptive integration for finite and infinite intervals • Root Finders • bracketing and polishing algorithms using derivatives • Minimization • Golden section and Brent algorithm • Interpolation • linear, polynomial, cubic and Akima spline • Chebyshev polynomials (for function approximation) • Random Numbers generation and distributions

9. Numerical Algorithms • Collect in the new MathCore basic numerical algorithms present in various places in ROOT • start with algorithms from TF1 (derivation, integration, etc..) • provide a coherent interface for these algorithms • maintain the current methods in TF1 for user convenience and backward compatibility • design based on classes already developed in MathMore • Derivator, Integrator, RootFinder, etc.. • provides a basic implementation using code from TF1 • use plug-in manager for alternative implementation (from GSL) • Algorithms could be used directly (with-out the need of creating a TF1 objects) • user would need to provide a C++ callable object (functor)

10. Example: Integration • Entry point is a Integrator class which will be in MathCore • can be used directly of via TF1 • Integrator class internally uses the chosen implementation: • Gauss Integrator : simple algorithm from MathCore (used now in TF1::Integral) • GSL Integrator based on QAGS (from MathMore)

11. Functor class • Use a Functor class to wrap any C++ callable object which has a right signature • Can work with: • free C functions • C++ classes (or structs) implementing operator() • member function of a class • The requirement is that they must have the right signature • Example: a function of type : • double f ( double )

12. Examples of using Functors double func(double x) { ............} Functor f1(func); f1(x); // evaluate functor at value x • Free Functions • Classes implementing operator() • Member functions structMyFunction { doubleoperator()(double x) {.....} }; Functor f2( MyFunction() ); f2(x); // evaluate functor at value x classMyClass { ...... doubleEval(double x) {.....} }; MyClass * p = new MyClass(); Functor f3( p, &MyClass::Eval ); f3(x); // evaluate functor at value x

13. Advantages of Functor classes • They are very convenient for users • easier to defining a functor than implementing a class following an abstract interface • Have value semantics (like shared pointers) • Very flexible • users can customize the callable objects using its state • this is not possible when using a free function • Example of usage in C++: • STL algorithm are based on similar concept • very powerful and flexible • boost::function (will be in next C++ standard)

14. Functor Usage in TF1 • Use a Functor class in TF1 instead of using a pointer to a free function (fFunction) • Functor class with needed signature(ROOT::Math::ParamFunctor) • Use this class as a data member of TF1 • Have template constructors for callable objects and for member functions // template constructors from general callable object template <typename Func> TF1(const char * name, Func f, Double_t xmin, Double_t xmax, Int_t npar) {..... fFunction = ROOT::Math::ParamFunctor(f); } // template constructors from a member function template <class PtrObj, typename MemFn> TF1(const char * name,const PtrObj & p, MemFn fn, Double_t xmin,Double_t xmax.. {..... fFunction = ROOT::Math::ParamFunctor(p,memFn); }

15. Example of Functor Usage • working example: • possible now with the current TF1 only by using globals objects doubleMyFunc (double *x, double *p ) { return TMath::Gaus(x[0],p[0],p[1] ); } structMyDerivFunc { MyDerivFunc(TF1 * f): fFunc(f) {} doubleoperator()(double *x,double *) const { return fFunc->Derivative(*x); } TF1 * fFunc; }; structMyIntegFunc { MyIntegFunc(TF1 * f): fFunc(f) {} doubleoperator()(double *x,double *) const { double a = fFunc->GetXmin(); return fFunc->Integral(a, *x); } TF1 * fFunc; }; voidtestTF1Functor() { double xmin = -10; double xmax = 10; TF1 * f1 = new TF1("f1",MyFunc,xmin,xmax,2); f1->SetParameters(0.,1.); f1->SetMaximum(3); f1->SetMinimum(-1); f1->Draw(); // create derivatives function using // the parameter of the parent function TF1 * f2 = new TF1("f2",MyDerivFunc(f1), xmin, xmax, 0); f2->SetLineColor(kBlue); f2->Draw("same"); // create integral function TF1 * f3 = new TF1("f3",MyIntegFunc(f1), xmin, xmax, 0); f3->SetLineColor(kRed); f3->Draw("same"); }

16. Functors in ROOT • Different signatures can be easily matched using some Functor adapters • using something similar to std::bind2nd functions taking two parameters can be adapted to function with one parameter • one could project a multidimensional function in a 1D function • Numerical Algorithm (Integration, Root Finder classes, etc..) will work also accepting a Functor • use a template method relying on the defined functor signature • Better to define also in TF1 an operator() for TF1::Eval • it would avoid need of an additional adapter • Expect to make the TF1 changes for the June release

17. Improvements in Pseudo-Random Number Generators • Mersenne-Twister generator (TRandom3) is now used for gRandom • fast and excellent pseudo-random quality • very long period, ~106000 • replace obsolete TRandom2 with TausWorthe generator from L’Ecuyer (fast and good quality, but not so long period: 1026) • add RanLux generator (TRandom1) • use a better linear congruential for TRandom • old one had seeding problems and a not uniform coverage • need to maintain for backward compatibility a generator based • on a state of only 32 bits (very short period : 231 ~ 109) • must NOT be used in any statistics application • we should rename classes: • more meaningful names and use M-T for the base class

18. Further TRandom Improvements • New faster algorithm for generating Gaussian number • acceptance-complement ratio method (W.Hörmann,G.Derflinger) from UNU.RAN • faster than polar method (Box-Muller) which requires 2 numbers + sin, sqrt and log calculations • one of the fastest existing methods • New Poisson algorithm • previous one was buggy for large mu values • exact algorithm (rejection from a Lorentzian distribution) • if μ constant not very efficient • better using algorithm from UNU.RAN • Improve performances of TRandom::Landau

19. Random Number Performances • generators: • Gauss random numbers Tests run on lxplus slc4 dual-core Intel 64 ns/call

20. Physics and Geometry Vectors • Classes for 3D and 4D vectors with their operations and transformations (rotations) • functionality as in CLHEP Vector and Geometry packages • Work done in collaboration with Fermilab computing group (M. Fischler, W. Brown and J. Marraffino) • Main features of the new classes: • generic scalar contained type • i.e. single or double precision • generic coordinate system concept • i.e. cartesian, polar and cylindrical • Used now by CMS and LHCb • Classes do not inherit from TObject and cannot be used in ROOT collections • Plan in the future to re-implement TLorentzVector using new classes

21. SMatrix Package • Package initially developed by T. Glebe for HeraB • Matrix and vector classes of arbitrary type • For fixed (not dynamic) matrix and vector sizes : • size must be knows at compile time • SMatrix< double, 2 , 5> • SVector< double, 5 > • Complementary and NOT a replacement of TMatrix • Optimized for small matrix sizes (dim <= 6): • use expression templates to avoid temporaries • Support for symmetric matrices (thanks to J.Palacios, LHCb) • storage of only n*(n+1)/2 elements • Support for basic operations and matrix inversion • not full linear algebra functionality • Used by LHCb, CMS and now ATLAS

22. Matrix Operations Performances • Comparison ROOT (TMatrix/SMatrix) and CLHEP (HepMatrix) • lxplus ( new Intel dual-core 64 bits) running slc4 with gcc 3.4.6

23. Kalman Filter Tests • CPU performances in the Kalman filter update equations • addition,multiplication and inversion • test varying matrix sizes • 2 <N1,N2 ≤ 6 • 6 < max(N1,N2) ≤ 10 • Useful exercise also for TMatrix • achieved substantial improvements since last workshop

24. Fitting and Minimization • New C++ version of Minuit (Minuit2) since ROOT v5.08 • Same basic functionality as in old version • Migrad, Simplex, Minos algorithms • Extended functionality: • single side parameter limits • added Fumili method for Chi2 and likelihood fits • implemented a ROOT fitter interface (TVirtualFitter) • TFitterMinuit and TFitterFumili • validated with extensive testing • OO package for generic function minimization • easy to extend by inserting new minimization algorithms • plan to add constrained minimization

25. Proposal for new ROOT Fitter • Proposal for new fitting classes • a mini fitting framework (like a simplified RooFit) • idea is to improve and extend functionality of TVirtualFitter with new set of fitting classes • Require a modular design: • easy possibility to extend and add new complex functionality : • perform parallel fits (use in a multi-threads environment) • add new minimization algorithms • implement parameter constraints • easy maintainability in the long term 1

26. Problems with TVirtualFitter • class designed for TMinuit, difficult to adapt for other minimizers • i.e. TVirtualFitter::ExecuteCommand • no separation Minimization-Fitting • it is more an interface for Minimization • users needs to provide the function to be minimized • must be a free function (signature needed by TMinuit): • func(int &, double *, double &f, double *par, int iflag) • there is no possibility to pass an object as function to be minimized • currently a global static instance of a TVirtualFitter is needed

27. Fitter Design Proposal

28. New Fitter Characteristics • Decoupling of Fitter from the various data sources • coupling only at the level of the FitData classes • tune the fit data according to the source • optimize memory vs CPU performances • Have an abstract interface for the Minimizer • instantiate the Minimizer classes via the plug-in manager • user can deal directly with minimizer interface • Have minimal Function interfaces • decouple Fitter from complex function objects like TF1 • describe only the Math functionality • evaluation, derivative, possibly integral (for the pdf) • state with parameters (for the model functions)

29. Fitter class • Fitter class glue together data and the model function • Fitter::Fit( IParamFunction & , const FitData & ) • create a concrete objective function • a Least Square or Likelihood function using a reference to the data and the model function • create a concrete Minimizer class with chosen configuration • use plug-in manager to load the corresponding library • TMinuit, Minuit2, MathMore (GSL Minimizer), etc... • according to the chosen implementation type(Minuit, Minuit2,Fumili, etc..) and configuration • find the minimum of the objective function • perform optionally error analysis (MINOS) • fill and return the FitResult class • parameter values, errors, error matrix, etc...

30. Current Status • Have already a prototype for Least Square fits working with a Minuit, Minuit2 and GSL implementations • Have in CVS next development cycle (after June release) • use then for FitPanel • re-implement TH1::Fit • re-implement TVirtualFitter methods to maintain backward compatibility // fit inputs TH1 * h1 = ..... TF1 * func = ..... ROOT::Fit::BinData d; // fill the data set from the histogram ROOT::Fit::FillData(d,h1); // create wrapped parametric function ROOT::Math::WrappedTF1 f(*func); // create fitter ROOT::Fit::Fitter // set minimizer and configuration fitter.Config().SetMinimizer(“Minuit2”); //perform the fit bool ret = fitter.Fit(d,f); //retrieve optionally the fit results if (ret) fitter.Result().Print();

31. Possible Extensions • Provide set of pre-defined functions (pdf) like in RooFit • have a catalog of the most used model functions • providing analytical implementations for the gradient, integral , etc.. • Provide eventually possibility to compose functions: • additions :h(x) = f(x) + g(x) • multiplications: h(x) = f(x)g(x) • composition: h(x) = g ( f(x) ) h(x,y) = f(x) g(y) • convolution:h(x) = ∫f(x-t)g(t)dt • Have a more complete type of fitting methods • Support for non-trivial constraints • using a minimizer which supports them

32. New GUI for Fitting • developed a new Fit panel (see Ilka presentation) • panel for fitting ROOT objects (TH1, TGraph etc...)

33. UNURAN • new package for generating non-uniform pseudo-random numbers • developed ROOT interface to a C package developed by J. Leydold et al. (Vienna) • various methods for 1D, multi-dimensional discrete and empirical distributions (from a set of data) • methods also for standard distributions (like Gauss, Poisson, ..) • see J. Leydold presentation on Wednesday • provides efficient, exact methods and works for non truncated functions • TF1::GetRandom() requires a range and it is approximate

34. μs/call Unuran 2.6 FOAM 2.55 TF1::GetRandom 0.25 Unuran Performances • Very efficient when distribution parameters are constant • Multi-dimensional distributions

35. FFT • Included in ROOT a common base class (TVirtualFFT) • add a functions to use it from TH1 (TH1::FFT) • Implemented an interface to the popular FFTW package (see www.fftw.org) • support for one and multi-dimensional transforms • support for complex and real transformations • TFFTComplex for complex input/ • complex output transforms • TFFTRealComplexfor real input/ • complex output • TFFTComplexRealfor complex input/real output • TFFTRealfor real input/output

36. Histogram Comparison • Improvements in Chi2 test for comparing histograms • algorithm from N. Gagunashvili and implemented in C++ by D. Haertl • add possibility to use weighted histograms • comparison of histograms with different scales • produce normalized residuals

37. Statistical Tools • TMVA: new package for multivariate analysis • Provides various methods for signal/background discrimination: • see presentation on Wednesday from H. Voss • Plan for a re-organization of statistical tools in ROOT • driven by requirements from LHC experiments • see Wednesday presentations: • RooStat (from K.Cranmer) • a new statistical tools based on RooFit • MadTools (from W. Quayle) • library with statistical methods needed for LHC analysis

38. Future Plans • MathCore • complete restructuring and work next months on the integration with ROOT analysis objects • remove some duplication (obsolete algorithms or functions) • MathMore • complete in MathMore the GSL wrapper • quasi-random numbers, differential equations, multi-dimensional algorithms (minimization, integration, etc..) • Complete design and developments of new fitting classes and then start integrating in ROOT • easier to use various fitting and minimization methods

39. Documentation • In addition to THtml provide also online doc based on Doxygen for the new classes of MathCore, MathMore, SMatrix and Minuit2 • provided for every new ROOT release, for example • http://www.cern.ch/mathlibs/sw/5_15_04/SMatrix/html/index.html • Written a new Math chapter in the ROOT User Guide (chapter 13, 227-248) • describe random numbers, MathCore (GenVector), mathematical functions, SMatrix • see ftp://root.cern.ch/root/doc/chapter13.pdf • Separate docs exist for other packages • Minuit2, RooFit, TMVA

40. References • MathCore online doc:http://www.cern.ch/mathlibs/sw/MathCore/html/index.html • MathMore online doc:http://www.cern.ch/mathlibs/sw/MathMore/html/index.html • SMatrix online doc: http://www.cern.ch/mathlibs/sw/SMatrix/html/index.html • Minuit2 online doc: http://www.cern.ch/mathlibs/sw/Minuit2/html/index.html • RooFit homepage: http://roofit.sourceforge.net/ • TMVA homepage: http://tmva.sourceforge.net/ • FFTW homepage: http://www.fftw.org/ • Histogram comparison paper: http://arxiv.org/abs/physics/0605123 • SPlot paper: http://arxiv.org/abs/physics/0402083 • UNURAN homepage: http://statmath.wu-wien.ac.at/unuran/