https://github.com/root-project/root
Tip revision: 0ce7176ca81ebb7ce89cb318dafffcd1456a2a55 authored by Fons Rademakers on 29 August 2013, 14:42:02 UTC
make version v5-34-10.
make version v5-34-10.
Tip revision: 0ce7176
stressMathCore.cxx
// @(#)root/test:$Id$
// Author: Lorenzo Moneta 06/2005
///////////////////////////////////////////////////////////////////////////////////
//
// MathCore Benchmark test suite
// ==============================
//
// This program performs tests :
// - mathematical functions in particular the statistical functions by estimating
// pdf, cdf and quantiles. cdf are estimated directly and compared with calculated integral from pdf
// - physics vectors (2D, 3D and 4D) including I/O for every type and for both double and Double32_t
// - SMatrix and SVectors including I/O for double and Double32_t types
// - I/O of complex objects which dictionary has been generated using CINT (default) or Reflex
// TrackD and TrackD32 which contain physics vectors of double and Double32_t
// TrackErrD and TrackErrD32 which contain physics vectors and an SMatrix of double and Double32_t
// VecTrackD which contains an std::vector<TrackD>
//
//
// the program cun run only in compiled mode.
// To run outside ROOT do:
//
// > cd $ROOTSYS/test
// > make stressMathMore
// > ./stressMathMore
//
// to run using REflex set before compiling the environment variable useReflex.
//
// > export useReflex=1
// > make stressMathMore
// > ./stressMathMore
//
// to run inside ROOT using ACliC
// for using CINT you need first to have the library libTrackMathCoreDict.so
// (type: make libTrackMathCoreDict.so to make it)
//
// root> gSystem->Load("libMathCore");
// root> gSystem->Load("libTree");
// root> gSystem->Load("libHist");
// root> .x stressMathCore.cxx+
//
// for using Reflex dictionaries you need first to have the library libTrackMathCoreRflx.so
// (type: make libTrackMathCoreRflx.so to make it)
//
// root> gSystem->Load("libMathCore");
// root> gSystem->Load("libTree");
// root> gSystem->Load("libHist");
// root> gSystem->Load("libReflex");
// root> gSystem->Load("libCintex");
// root> gSystem->SetIncludePath("-DUSE_REFLEX");
// root> .x stressMathCore.cxx+
//
//
#ifndef __CINT__
#include "Math/DistFuncMathCore.h"
#ifdef USE_MATHMORE
#include "Math/DistMathMore.h"
#endif
#include "Math/IParamFunction.h"
#include "Math/Integrator.h"
#include <iostream>
#include <iomanip>
#include <limits>
#include <cmath>
#include "TBenchmark.h"
#include "TROOT.h"
#include "TRandom3.h"
#include "TSystem.h"
#include "TTree.h"
#include "TFile.h"
#include "TF1.h"
#include "Math/Vector2D.h"
#include "Math/Vector3D.h"
#include "Math/Vector4D.h"
#include "Math/VectorUtil.h"
#include "Math/SVector.h"
#include "Math/SMatrix.h"
#include "TrackMathCore.h"
//#define USE_REFLEX
#ifdef USE_REFLEX
#include "Cintex/Cintex.h"
#include "Reflex/Reflex.h"
#endif
using namespace ROOT::Math;
#endif
//#define DEBUG
bool debug = true; // print out reason of test failures
bool debugTime = false; // print out separate timings for vectors
bool removeFiles = true; // remove Output root files
void PrintTest(std::string name) {
std::cout << std::left << std::setw(40) << name;
}
void PrintStatus(int iret) {
if (iret == 0)
std::cout <<"\t\t................ OK" << std::endl;
else
std::cout <<"\t\t............ FAILED " << std::endl;
}
int compare( std::string name, double v1, double v2, double scale = 2.0) {
// ntest = ntest + 1;
//std::cout << std::setw(50) << std::left << name << ":\t";
// numerical double limit for epsilon
double eps = scale* std::numeric_limits<double>::epsilon();
int iret = 0;
double delta = v2 - v1;
double d = 0;
if (delta < 0 ) delta = - delta;
if (v1 == 0 || v2 == 0) {
if (delta > eps ) {
iret = 1;
}
}
// skip case v1 or v2 is infinity
else {
d = v1;
if ( v1 < 0) d = -d;
// add also case when delta is small by default
if ( delta/d > eps && delta > eps )
iret = 1;
}
if (iret) {
if (debug) {
int pr = std::cout.precision (18);
std::cout << "\nDiscrepancy in " << name.c_str() << "() :\n " << v1 << " != " << v2 << " discr = " << int(delta/d/eps)
<< " (Allowed discrepancy is " << eps << ")\n\n";
std::cout.precision (pr);
//nfail = nfail + 1;
}
}
//else
// std::cout <<".";
return iret;
}
#ifndef __CINT__
// trait class for distinguishing the number of parameters for the various functions
template<class Func, unsigned int NPAR>
struct Evaluator {
static double F(Func f, double x, const double * ) {
return f(x);
}
};
template<class Func>
struct Evaluator<Func, 1> {
static double F(Func f, double x, const double * p) {
return f(x,p[0]);
}
};
template<class Func>
struct Evaluator<Func, 2> {
static double F(Func f, double x, const double * p) {
return f(x,p[0],p[1]);
}
};
template<class Func>
struct Evaluator<Func, 3> {
static double F(Func f, double x, const double * p) {
return f(x,p[0],p[1],p[2]);
}
};
// global test variable
int NFuncTest = 100;
// statistical function class
// template on the number of parameters
template<class Func, class FuncQ, int NPAR, int NPARQ=NPAR-1>
class StatFunction : public ROOT::Math::IParamFunction {
public:
StatFunction(Func pdf, Func cdf, FuncQ quant) : fPdf(pdf), fCdf(cdf), fQuant(quant)
{
fScale1 = 1.0E6; //scale for cdf test (integral)
fScale2 = 10; //scale for quantile test
for(int i = 0; i< NPAR; ++i) fParams[i]=0;
}
unsigned int NPar() const { return NPAR; }
const double * Parameters() const { return fParams; }
ROOT::Math::IGenFunction * Clone() const { return new StatFunction(fPdf,fCdf,fQuant); }
void SetParameters(const double * p) { std::copy(p,p+NPAR,fParams); }
void SetParameters(double p0) { *fParams = p0; }
void SetParameters(double p0, double p1) { *fParams = p0; *(fParams+1) = p1; }
void SetParameters(double p0, double p1, double p2) { *fParams = p0; *(fParams+1) = p1; *(fParams+2) = p2; }
static void SetNTest(int n) { NFuncTest = n; }
double Cdf(double x) const {
return Evaluator<Func,NPAR>::F(fCdf,x, fParams);
}
double Quantile(double x) const {
double z = Evaluator<FuncQ,NPARQ>::F(fQuant,x, fParams);
if ((NPAR - NPARQ) == 1)
z += fParams[NPAR-1]; // adjust the offset
return z;
}
// test cumulative function
int Test(double x1, double x2, double xl = 1, double xu = 0, bool cumul = false);
void ScaleTol1(double s) { fScale1 *= s; }
void ScaleTol2(double s) { fScale2 *= s; }
private:
double DoEvalPar(double x, const double * ) const {
// implement explicitly using cached parameter values
return Evaluator<Func,NPAR>::F(fPdf,x, fParams);
}
Func fPdf;
Func fCdf;
FuncQ fQuant;
double fParams[NPAR];
double fScale1;
double fScale2;
};
// test cdf at value f
template<class F1, class F2, int N1, int N2>
int StatFunction<F1,F2,N1,N2>::Test(double xmin, double xmax, double xlow, double xup, bool c) {
int iret = 0;
// scan all values from xmin to xmax
double dx = (xmax-xmin)/NFuncTest;
#ifdef USE_MATHMORE
for (int i = 0; i < NFuncTest; ++i) {
double v1 = xmin + dx*i; // value used for testing
double q1 = Cdf(v1);
//std::cout << "v1 " << v1 << " pdf " << (*this)(v1) << " cdf " << q1 << " quantile " << Quantile(q1) << std::endl;
// calculate integral of pdf
Integrator ig(IntegrationOneDim::ADAPTIVESINGULAR, 1.E-12,1.E-12,100000);
ig.SetFunction(*this);
double q2 = 0;
if (!c) {
// lower intergal (cdf)
if (xlow >= xup || xlow > xmin)
q2 = ig.IntegralLow(v1);
else
q2 = ig.Integral(xlow,v1);
// use a larger scale (integral error is 10-9)
iret |= compare("test _cdf", q1, q2, fScale1 );
// test the quantile
double v2 = Quantile(q1);
iret |= compare("test _quantile", v1, v2, fScale2 );
}
else {
// upper integral (cdf_c)
if (xlow >= xup || xup < xmax)
q2 = ig.IntegralUp(v1);
else
q2 = ig.Integral(v1,xup);
iret |= compare("test _cdf_c", q1, q2, fScale1);
double v2 = Quantile(q1);
iret |= compare("test _quantile_c", v1, v2, fScale2 );
}
if (iret) {
std::cout << "Failed test for x = " << v1 << " p = ";
for (int j = 0; j < N1; ++j) std::cout << fParams[j] << "\t";
std::cout << std::endl;
break;
}
}
#else
// use TF1 for the integral
// std::cout << "xlow-xuo " << xlow << " " << xup << std::endl;
// std::cout << "xmin-xmax " << xmin << " " << xmax << std::endl;
double x1,x2 = 0;
if (xlow >= xup) {
x1 = -100; x2 = 100;
}
else if (xup < xmax) {
x1 = xlow; x2 = 100;
}
else {
x1=xlow; x2 = xup;
}
//std::cout << "x1-x2 " << x1 << " " << x2 << std::endl;
TF1 * f = new TF1("ftemp",ParamFunctor(*this),x1,x2,0);
for (int i = 0; i < NFuncTest; ++i) {
double v1 = xmin + dx*i; // value used for testing
double q1 = Cdf(v1);
//std::cout << "i = " << i << " v1 = " << v1 << " pdf " << (*this)(v1) << " cdf " << q1 << std::endl;
double q2 = 0;
if (!c) {
q2 = f->Integral(x1,v1);
// use a larger scale (integral error is 10-9)
iret |= compare("test _cdf", q1, q2, fScale1 );
// test the quantile
double v2 = Quantile(q1);
iret |= compare("test _quantile", v1, v2, fScale2 );
}
else {
// upper integral (cdf_c)
q2 = f->Integral(v1,x2);
iret |= compare("test _cdf_c", q1, q2, fScale1);
double v2 = Quantile(q1);
iret |= compare("test _quantile_c", v1, v2, fScale2 );
}
}
delete f;
#endif
if (c || iret != 0) PrintStatus(iret);
return iret;
}
// typedef defining the functions
typedef double ( * F0) ( double);
typedef double ( * F1) ( double, double);
typedef double ( * F2) ( double, double, double);
typedef double ( * F3) ( double, double, double, double);
typedef StatFunction<F2,F2,2,2> Dist_beta;
typedef StatFunction<F2,F1,2> Dist_breitwigner;
typedef StatFunction<F2,F1,2> Dist_chisquared;
typedef StatFunction<F3,F2,3> Dist_fdistribution;
typedef StatFunction<F3,F2,3> Dist_gamma;
typedef StatFunction<F2,F1,2> Dist_gaussian;
typedef StatFunction<F3,F2,3> Dist_lognormal;
typedef StatFunction<F2,F1,2> Dist_tdistribution;
typedef StatFunction<F2,F1,2> Dist_exponential;
typedef StatFunction<F2,F1,2> Dist_landau;
typedef StatFunction<F3,F2,3> Dist_uniform;
#define CREATE_DIST(name) Dist_ ##name dist( name ## _pdf, name ## _cdf, name ##_quantile );
#define CREATE_DIST_C(name) Dist_ ##name distc( name ## _pdf, name ## _cdf_c, name ##_quantile_c );
// #define CREATE_DIST_C(name) Dist_ ##name distc( name ## _pdf, name ## _cdf_c, 0 );
// #else
// #endif
template<class Distribution>
int TestDist(Distribution & d, double x1, double x2) {
int ir = 0;
ir |= d.Test(x1,x2);
return ir;
}
int testStatFunctions(int nfunc = 100 ) {
// test statistical functions
int iret = 0;
NFuncTest = nfunc;
{
PrintTest("Beta distribution");
CREATE_DIST(beta);
dist.SetParameters( 2, 2);
iret |= dist.Test(0.01,0.99,0.,1.);
CREATE_DIST_C(beta);
distc.SetParameters( 2, 2);
iret |= distc.Test(0.01,0.99,0.,1.,true);
}
{
PrintTest("Gamma distribution");
#ifdef USE_MATHMORE // gamma_quantile is in mathmore
CREATE_DIST(gamma);
dist.SetParameters( 2, 1);
iret |= dist.Test(0.05,5, 0.,1.);
#endif
CREATE_DIST_C(gamma);
distc.SetParameters( 2, 1);
iret |= distc.Test(0.05,5, 0.,1.,true);
}
{
PrintTest("Chisquare distribution");
#ifdef USE_MATHMORE
CREATE_DIST(chisquared);
dist.SetParameters( 10, 0);
dist.ScaleTol2(10);
iret |= dist.Test(0.05,30, 0.,1.);
#endif
CREATE_DIST_C(chisquared);
distc.SetParameters( 10, 0);
distc.ScaleTol2(10000000); // t.b.c.
iret |= distc.Test(0.05,30, 0.,1.,true);
}
{
PrintTest("Normal distribution ");
CREATE_DIST(gaussian);
dist.SetParameters( 2, 1);
dist.ScaleTol2(100);
iret |= dist.Test(-3,5);
CREATE_DIST_C(gaussian);
distc.SetParameters( 1, 0);
distc.ScaleTol2(100);
iret |= distc.Test(-3,5,1,0,true);
}
{
PrintTest("BreitWigner distribution ");
CREATE_DIST(breitwigner);
dist.SetParameters( 1);
#ifndef USE_MATHMORE
dist.ScaleTol1(1E8);
#endif
dist.ScaleTol2(10);
iret |= dist.Test(-5,5);
CREATE_DIST_C(breitwigner);
distc.SetParameters( 1);
#ifndef USE_MATHMORE
distc.ScaleTol1(1E8);
#endif
distc.ScaleTol2(10);
iret |= distc.Test(-5,5,1,0,true);
}
{
PrintTest("F distribution ");
CREATE_DIST(fdistribution);
dist.SetParameters( 5, 4);
dist.ScaleTol1(1000000);
dist.ScaleTol2(10);
// if enlarge scale test fails
iret |= dist.Test(0.05,5,0,1);
CREATE_DIST_C(fdistribution);
distc.SetParameters( 5, 4);
#ifndef USE_MATHMORE
distc.ScaleTol1(100000000);
#endif
distc.ScaleTol2(10);
// if enlarge scale test fails
iret |= distc.Test(0.05,5,0,1,true);
}
#ifdef USE_MATHMORE // wait t quantile is in mathcore
{
PrintTest("t distribution ");
CREATE_DIST(tdistribution);
dist.SetParameters( 10 );
// dist.ScaleTol1(1000);
dist.ScaleTol2(5000);
iret |= dist.Test(-10,10);
CREATE_DIST_C(tdistribution);
distc.SetParameters( 10 );
distc.ScaleTol2(10000); // t.b.c.
iret |= distc.Test(-10,10,1,0,true);
}
#endif
{
PrintTest("lognormal distribution");
CREATE_DIST(lognormal);
dist.SetParameters(1,1 );
dist.ScaleTol1(1000);
iret |= dist.Test(0.01,5,0,1);
CREATE_DIST_C(lognormal);
distc.SetParameters(1,1 );
#ifndef USE_MATHMORE
distc.ScaleTol1(1000000);
#endif
distc.ScaleTol2(1000000); // t.b.c.
iret |= distc.Test(0.01,5,0,1,true);
}
{
PrintTest("Exponential distribution");
CREATE_DIST(exponential);
dist.SetParameters( 2);
dist.ScaleTol2(100);
iret |= dist.Test(0.,5.,0.,1.);
CREATE_DIST_C(exponential);
distc.SetParameters( 2);
distc.ScaleTol2(100);
iret |= distc.Test(0.,5.,0.,1.,true);
}
{
PrintTest("Landau distribution");
CREATE_DIST(landau);
dist.SetParameters( 2);
// Landau is not very precise (put prec at 10-6)
// as indicated in Landau paper (
dist.ScaleTol1(10000);
dist.ScaleTol2(1.E10);
iret |= dist.Test(-1,10,-10.,1.E10);
CREATE_DIST_C(landau);
distc.SetParameters( 2);
distc.ScaleTol1(10000);
distc.ScaleTol2(1.0E10);
iret |= distc.Test(-1,10,-10.,1.E10,true);
}
{
PrintTest("Uniform distribution");
CREATE_DIST(uniform);
dist.SetParameters( 1, 2);
iret |= dist.Test(1.,2.,1.,2.);
CREATE_DIST_C(uniform);
distc.SetParameters( 1, 2);
iret |= distc.Test(1.,2.,1.,2.,true);
}
return iret;
}
//*******************************************************************************************************************
// GenVector tests
//*******************************************************************************************************************
// trait for getting vector name
template<class V>
struct VecType {
static std::string name() { return "MathCoreVector";}
};
template<>
struct VecType<XYVector> {
static std::string name() { return "XYVector";}
static std::string name32() { return "ROOT::Math::DisplacementVector2D<ROOT::Math::Cartesian2D<Double32_t> >";}
};
template<>
struct VecType<Polar2DVector> {
static std::string name() { return "Polar2DVector";}
};
template<>
struct VecType<XYZVector> {
static std::string name() { return "XYZVector";}
static std::string name32() { return "ROOT::Math::DisplacementVector3D<ROOT::Math::Cartesian3D<Double32_t> >";}
};
template<>
struct VecType<Polar3DVector> {
static std::string name() { return "Polar3DVector";}
};
template<>
struct VecType<RhoEtaPhiVector> {
static std::string name() { return "RhoEtaPhiVector";}
};
template<>
struct VecType<RhoZPhiVector> {
static std::string name() { return "RhoZPhiVector";}
};
template<>
struct VecType<PxPyPzEVector> {
static std::string name() { return "PxPyPzEVector";}
static std::string name32() { return "ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<Double32_t> >";}
};
template<>
struct VecType<PtEtaPhiEVector> {
static std::string name() { return "PtEtaPhiEVector";}
};
template<>
struct VecType<PtEtaPhiMVector> {
static std::string name() { return "PtEtaPhiMVector";}
};
template<>
struct VecType<PxPyPzMVector> {
static std::string name() { return "PxPyPzMVector";}
};
template<>
struct VecType<SVector<double,3> > {
static std::string name() { return "SVector3";}
static std::string name32() { return "ROOT::Math::SVector<Double32_t,3>";}
};
template<>
struct VecType<SVector<double,4> > {
static std::string name() { return "SVector4";}
static std::string name32() { return "ROOT::Math::SVector<Double32_t,4>";}
};
template<>
struct VecType<TrackD> {
static std::string name() { return "TrackD";}
};
template<>
struct VecType<TrackD32> {
static std::string name() { return "TrackD32";}
};
template<>
struct VecType<TrackErrD> {
static std::string name() { return "TrackErrD";}
};
template<>
struct VecType<TrackErrD32> {
static std::string name() { return "TrackErrD32";}
};
template<>
struct VecType<VecTrack<TrackD> > {
static std::string name() { return "VecTrackD";}
};
template<>
struct VecType<VecTrack<TrackErrD> > {
static std::string name() { return "VecTrackErrD";}
};
// generic (2 dim)
template<class V, int Dim>
struct VecOp {
template<class It>
static V Create(It &x, It &y, It & , It& ) { return V(*x++,*y++); }
template<class It>
static void Set(V & v, It &x, It &y, It &, It&) { v.SetXY(*x++,*y++); }
static double Add(const V & v) { return v.X() + v.Y(); }
static double Delta(const V & v1, const V & v2) { double d = ROOT::Math::VectorUtil::DeltaPhi(v1,v2); return d*d; } // is v2-v1
};
// specialized for 3D
template<class V>
struct VecOp<V,3> {
template<class It>
static V Create(It &x, It& y, It& z , It& ) { return V(*x++,*y++,*z++); }
template<class It>
static void Set(V & v, It & x, It &y, It &z, It&) { v.SetXYZ(*x++,*y++,*z++); }
static V Create(double x, double y, double z , double ) { return V(x,y,z); }
static void Set(V & v, double x, double y, double z, double) { v.SetXYZ(x,y,z); }
static double Add(const V & v) { return v.X() + v.Y() + v.Z(); }
static double Delta(const V & v1, const V & v2) { return ROOT::Math::VectorUtil::DeltaR2(v1,v2); }
};
// specialized for 4D
template<class V>
struct VecOp<V,4> {
template<class It>
static V Create(It &x, It &y, It &z , It &t ) { return V(*x++,*y++,*z++,*t++);}
template<class It>
static void Set(V & v, It & x, It &y, It &z, It &t) { v.SetXYZT(*x++,*y++,*z++,*t++); }
static double Add(const V & v) { return v.X() + v.Y() + v.Z() + v.E(); }
static double Delta(const V & v1, const V & v2) {
return ROOT::Math::VectorUtil::DeltaR2(v1,v2) + ROOT::Math::VectorUtil::InvariantMass(v1,v2); }
};
// specialized for SVector<3>
template<>
struct VecOp<SVector<double,3>,3> {
typedef SVector<double,3> V;
template<class It>
static V Create(It &x, It &y, It &z , It & ) { return V(*x++,*y++,*z++);}
static double Add(const V & v) { return v(0) + v(1) + v(2); }
};
// specialized for SVector<4>
template<>
struct VecOp<SVector<double,4>,4> {
typedef SVector<double,4> V;
template<class It>
static V Create(It &x, It &y, It &z , It &t ) { return V(*x++,*y++,*z++,*t++);}
static double Add(const V & v) { return v(0) + v(1) + v(2) + v(3); }
};
// internal structure to measure the time
TStopwatch gTimer;
double gTotTime;
struct Timer {
Timer() {
gTimer.Start();
}
~Timer() {
gTimer.Stop();
gTotTime += Time();
if (debugTime) printTime();
}
void printTime( std::string s = "") {
int pr = std::cout.precision(8);
std::cout << s << "\t" << " time = " << Time() << "\t(sec)\t"
// << time.CpuTime()
<< std::endl;
std::cout.precision(pr);
}
double Time() { return gTimer.RealTime(); } // use real time
TStopwatch gTimer;
double gTotTime;
};
template<int Dim>
class VectorTest {
private:
// global data variables
std::vector<double> dataX;
std::vector<double> dataY;
std::vector<double> dataZ;
std::vector<double> dataE;
int nGen;
int n2Loop ;
double fSum; // total sum of x,y,z,t (for testing first addition)
public:
VectorTest(int n1, int n2=0) :
nGen(n1),
n2Loop(n2)
{
gTotTime = 0;
}
double TotalTime() const { return gTotTime; } // use real time
double Sum() const { return fSum; }
int check(std::string name, double s1, double s2, double scale=1) {
int iret = 0;
PrintTest(name);
iret |= compare(name,s1,s2,scale);
PrintStatus(iret);
return iret;
}
void print(std::string name) {
PrintTest(name);
std::cout <<"\t\t..............\n";
}
void genData() {
// generate for all 4 d data
TRandom3 r(111); // use a fixed seed to be able to reproduce tests
fSum = 0;
for (int i = 0; i < nGen ; ++ i) {
// generate a 4D vector and stores only the interested dimensions
double phi = r.Rndm()*3.1415926535897931;
double eta = r.Uniform(-5.,5.);
double pt = r.Exp(10.);
double m = r.Uniform(0,10.);
if ( i%50 == 0 )
m = r.BreitWigner(1.,0.01);
double E = sqrt( m*m + pt*pt*std::cosh(eta)*std::cosh(eta) );
// fill vectors
PtEtaPhiEVector q( pt, eta, phi, E);
dataX.push_back( q.x() );
dataY.push_back( q.y() );
fSum += q.x() + q.y();
if (Dim >= 3) {
dataZ.push_back( q.z() );
fSum += q.z();
}
if (Dim >=4 ) {
dataE.push_back( q.t() );
fSum += q.t();
}
}
assert( int(dataX.size()) == nGen);
assert( int(dataY.size()) == nGen);
if (Dim >= 3) assert( int(dataZ.size()) == nGen);
if (Dim >=4 ) assert( int(dataE.size()) == nGen);
// // // dataZ.resize(nGen);
// // // dataE.resize(nGen);
}
// gen data for a Ndim matrix or vector
void genDataN() {
// generate for all 4 d data
TRandom3 r(111); // use a fixed seed to be able to reproduce tests
fSum = 0;
dataX.reserve(nGen*Dim);
for (int i = 0; i < nGen*Dim ; ++ i) {
// generate random data between [0,1]
double x = r.Rndm();
fSum += x;
dataX.push_back( x );
}
}
typedef std::vector<double>::const_iterator DataIt;
// test methods
template <class V>
void testCreate( std::vector<V > & dataV) {
Timer tim;
DataIt x = dataX.begin();
DataIt y = dataY.begin();
DataIt z = dataZ.begin();
DataIt t = dataE.begin();
while (x != dataX.end() ) {
dataV.push_back(VecOp<V,Dim>::Create(x,y,z,t) );
assert(int(dataV.size()) <= nGen);
}
}
template <class V>
void testCreateAndSet( std::vector<V > & dataV) {
Timer tim;
DataIt x = dataX.begin();
DataIt y = dataY.begin();
DataIt z = dataZ.begin();
DataIt t = dataE.begin();
while (x != dataX.end() ) {
V v;
VecOp<V,Dim>::Set( v, x,y,z,t);
dataV.push_back(v);
assert(int(dataV.size()) <= nGen);
}
}
template <class V>
double testAddition( const std::vector<V > & dataV) {
V v0;
Timer t;
for (int i = 0; i < nGen; ++i) {
v0 += dataV[i];
}
return VecOp<V,Dim>::Add(v0);
}
template <class V>
double testOperations( const std::vector<V > & dataV) {
double tot = 0;
Timer t;
for (int i = 0; i < nGen-1; ++i) {
const V & v1 = dataV[i];
const V & v2 = dataV[i+1];
double a = v1.R();
double b = v2.mag2(); // mag2 is defined for all dimensions;
double c = 1./v1.Dot(v2);
V v3 = c * ( v1/a + v2/b );
tot += VecOp<V,Dim>::Add(v3);
}
return tot;
}
// mantain loop in gen otherwise is proportional to N**@
template <class V>
double testDelta( const std::vector<V > & dataV) {
double tot = 0;
Timer t;
for (int i = 0; i < nGen-1; ++i) {
const V & v1 = dataV[i];
const V & v2 = dataV[i+1];
tot += VecOp<V,Dim>::Delta(v1,v2);
}
return tot;
}
// template <class V>
// double testDotProduct( const std::vector<V *> & dataV) {
// //unsigned int n = std::min(n2Loop, dataV.size() );
// double tot = 0;
// V v0 = *(dataV[0]);
// Timer t;
// for (unsigned int i = 0; i < nGen-1; ++i) {
// V & v1 = *(dataV[i]);
// tot += v0.Dot(v1);
// }
// return tot;
// }
template <class V1, class V2>
void testConversion( std::vector<V1 > & dataV1, std::vector<V2 > & dataV2) {
Timer t;
for (int i = 0; i < nGen; ++i) {
dataV2.push_back( V2( dataV1[i] ) );
}
}
// rotation
template <class V, class R>
double testRotation( std::vector<V > & dataV ) {
double sum = 0;
double rotAngle = 1;
Timer t;
for (unsigned int i = 0; i < nGen; ++i) {
V & v1 = dataV[i];
V v2 = v1;
v2.Rotate(rotAngle);
sum += VecOp<V,Dim>::Add(v2);
}
return sum;
}
template<class V>
double testWrite(const std::vector<V> & dataV, std::string typeName="", bool compress = false) {
std::string fname = VecType<V>::name() + ".root";
// replace < character with _
TFile file(fname.c_str(),"RECREATE","",compress);
// create tree
std::string tree_name="Tree with" + VecType<V>::name();
TTree tree("VectorTree",tree_name.c_str());
V *v1 = new V();
//std::cout << "typeID written : " << typeid(*v1).name() << std::endl;
// need to add namespace to full type name
if (typeName == "") {
typeName = "ROOT::Math::" + VecType<V>::name();
}
//std::cout << typeName << std::endl;
TBranch * br = tree.Branch("Vector branch",typeName.c_str(),&v1);
if (br == 0) {
std::cout << "Error creating branch for" << typeName << "\n\t typeid is "
<< typeid(*v1).name() << std::endl;
return -1;
}
Timer timer;
for (int i = 0; i < nGen; ++i) {
*v1 = dataV[i];
tree.Fill();
}
#ifdef DEBUG
tree.Print(); // debug
#endif
file.Write();
file.Close();
return file.GetSize();
}
template<class V>
int testRead(std::vector<V> & dataV) {
dataV.clear();
dataV.reserve(nGen);
std::string fname = VecType<V>::name() + ".root";
TFile f1(fname.c_str());
if (f1.IsZombie() ) {
std::cout << " Error opening file " << fname << std::endl;
return -1;
}
// create tree
TTree *tree = dynamic_cast<TTree*>(f1.Get("VectorTree"));
if (tree == 0) {
std::cout << " Error reading file " << fname << std::endl;
return -1;
}
V *v1 = 0;
//std::cout << "reading typeID : " << typeid(*v1).name() << std::endl;
// cast to void * to avoid a warning
tree->SetBranchAddress("Vector branch",(void *) &v1);
Timer timer;
int n = (int) tree->GetEntries();
if (n != nGen) {
std::cout << "wrong tree entries from file" << fname << std::endl;
return -1;
}
for (int i = 0; i < n; ++i) {
tree->GetEntry(i);
dataV.push_back(*v1);
}
if (removeFiles) gSystem->Unlink(fname.c_str());
return 0;
}
// test of SVEctor's or SMatrix
template<class V>
void testCreateSV( std::vector<V> & dataV) {
Timer tim;
//DataIt x = dataX.begin();
double * x = &dataX.front();
double * end = x + dataX.size();
//SVector cannot be created from a generic iterator (should be fixed)
while (x != end) {
dataV.push_back( V(x,x+Dim) );
assert(int(dataV.size()) <= nGen);
x += Dim;
}
}
template <class V>
double testAdditionSV( const std::vector<V > & dataV) {
V v0;
Timer t;
for (int i = 0; i < nGen; ++i) {
v0 += dataV[i];
}
double tot = 0;
typedef typename V::const_iterator It;
for (It itr = v0.begin(); itr != v0.end(); ++itr)
tot += *itr;
return tot;
}
template <class V>
double testAdditionTR( const std::vector<V > & dataV) {
V v0;
Timer t;
for (int i = 0; i < nGen; ++i) {
v0 += dataV[i];
}
#ifdef DEBUG
v0.Print();
#endif
return v0.Sum();
}
};
//--------------------------------------------------------------------------------------
// test of all physics vector (GenVector's)
//--------------------------------------------------------------------------------------
template<class V1, class V2, int Dim>
int testVector(int ngen, bool testio=false) {
int iret = 0;
VectorTest<Dim> a(ngen);
a.genData();
std::vector<V1> v1;
std::vector<V2> v2;
v1.reserve(ngen);
v2.reserve(ngen);
double s1, s2 = 0;
double scale = 1;
double sref1, sref2 = 0;
a.testCreate(v1); iret |= a.check(VecType<V1>::name()+" creation",v1.size(),ngen);
s1 = a.testAddition(v1); iret |= a.check(VecType<V1>::name()+" addition",s1,a.Sum(),Dim*20);
sref1 = s1;
v1.clear();
assert(v1.size() == 0);
a.testCreateAndSet(v1); iret |= a.check(VecType<V1>::name()+" creation",v1.size(),ngen);
s2 = a.testAddition(v1); iret |= a.check(VecType<V1>::name()+" setting",s2,s1);
a.testConversion(v1,v2); iret |= a.check(VecType<V1>::name()+" -> " + VecType<V2>::name(),v1.size(),v2.size() );
scale = 1000;
if (Dim == 2) scale = 1.E12; // to be understood
if (Dim == 3) scale = 1.E4; // problem with RhoEtaPhiVector
s2 = a.testAddition(v2); iret |= a.check("Vector conversion",s2,s1,scale);
sref2 = s2;
s1 = a.testOperations(v1); a.print(VecType<V1>::name()+" operations");
scale = Dim*20;
if (Dim==3 && VecType<V2>::name() == "RhoEtaPhiVector") scale *= 10; // for problem with RhoEtaPhi
if (Dim==4 && VecType<V2>::name() == "PtEtaPhiMVector") scale *= 10;
#if defined (R__LINUX) && !defined(R__B64)
// problem of precision on linux 32
if (Dim ==4) scale = 1000000000;
#endif
// for problem with PtEtaPhiE
if (Dim==4 && VecType<V2>::name() == "PtEtaPhiEVector") scale = 0.01/(std::numeric_limits<double>::epsilon());
s2 = a.testOperations(v2); iret |= a.check(VecType<V2>::name()+" operations",s2,s1,scale);
s1 = a.testDelta(v1); a.print(VecType<V1>::name()+" delta values");
scale = Dim*16;
if (Dim==4) scale *= 100; // for problem with PtEtaPhiE
s2 = a.testDelta(v2); iret |= a.check(VecType<V2>::name()+" delta values",s2,s1,scale);
double fsize = 0;
int ir = 0;
if (!testio) return iret;
double estSize = ngen*8 * Dim + 10000; //add extra bytes
scale = 0.1 / std::numeric_limits<double>::epsilon();
fsize = a.testWrite(v1); iret |= a.check(VecType<V1>::name()+" write",fsize,estSize,scale);
ir = a.testRead(v1); iret |= a.check(VecType<V1>::name()+" read",ir,0);
s1 = a.testAddition(v1); iret |= a.check(VecType<V1>::name()+" after read",s1,sref1);
// test io vector 2
fsize = a.testWrite(v2); iret |= a.check(VecType<V2>::name()+" write",fsize,estSize,scale);
ir = a.testRead(v2); iret |= a.check(VecType<V2>::name()+" read",ir,0);
scale = 100; // gcc 4.3.2 gives and error for RhoEtaPhiVector for 32 bits
s2 = a.testAddition(v2); iret |= a.check(VecType<V2>::name()+" after read",s2,sref2,scale);
// test io of double 32 for vector 1
if (Dim==2) return iret;
estSize = ngen*4 * Dim + 10000; //add extra bytes
scale = 0.1 / std::numeric_limits<double>::epsilon();
std::string typeName = VecType<V1>::name32();
fsize = a.testWrite(v1,typeName); iret |= a.check(VecType<V1>::name() +"_D32 write",fsize,estSize,scale);
ir = a.testRead(v1); iret |= a.check(VecType<V1>::name()+"_D32 read",ir,0);
s1 = a.testAddition(v1); iret |= a.check(VecType<V1>::name()+"_D32 after read",s1,sref1,1.E9);
return iret;
}
//--------------------------------------------------------------------------------------
// test of Svector of dim 3 or 4
//--------------------------------------------------------------------------------------
template<int Dim>
int testVector34(int ngen, bool testio=false) {
int iret = 0;
VectorTest<Dim> a(ngen);
a.genData();
typedef SVector<double, Dim> SV;
std::vector<SV> v1;
v1.reserve(ngen);
double s1 = 0;
//double scale = 1;
double sref1 = 0;
std::string name = "SVector<double," + Util::ToString(Dim) + ">";
a.testCreate(v1); iret |= a.check(name+" creation",v1.size(),ngen);
s1 = a.testAddition(v1); iret |= a.check(name+" addition",s1,a.Sum(),Dim*20);
sref1 = s1;
// test the io
double fsize = 0;
int ir = 0;
if (!testio) return iret;
std::string typeName = "ROOT::Math::"+name;
double estSize = ngen*8 * Dim + 47000; // add extra bytes (why so much ? ) this is for ngen=10000
double scale = 0.1 / std::numeric_limits<double>::epsilon();
fsize = a.testWrite(v1,typeName); iret |= a.check(name+" write",fsize,estSize,scale);
ir = a.testRead(v1); iret |= a.check(name+" read",ir,0);
s1 = a.testAddition(v1); iret |= a.check(name+" after read",s1,sref1);
//std::cout << "File size = " << fsize << " estimated " << 8 * Dim * ngen << std::endl;
// test Double32
estSize = ngen*4 * Dim + 47000; //add extra bytes
scale = 0.1 / std::numeric_limits<double>::epsilon();
typeName = VecType<SV>::name32();
fsize = a.testWrite(v1,typeName); iret |= a.check(VecType<SV>::name() +"_D32 write",fsize,estSize,scale);
ir = a.testRead(v1); iret |= a.check(VecType<SV>::name()+"_D32 read",ir,0);
s1 = a.testAddition(v1); iret |= a.check(VecType<SV>::name()+"_D32 after read",s1,sref1,1.E9);
return iret;
}
//--------------------------------------------------------------------------------------
// test of generic Svector
//--------------------------------------------------------------------------------------
template<int Dim>
int testSVector(int ngen, bool testio=false) {
// test the matrix if D2 is not equal to 1
int iret = 0;
VectorTest<Dim> a(ngen);
a.genDataN();
typedef SVector<double, Dim> SV;
std::vector<SV> v1;
v1.reserve(ngen);
double s1 = 0;
//double scale = 1;
double sref1 = 0;
std::string name = "SVector<double," + Util::ToString(Dim) + ">";
a.testCreateSV(v1); iret |= a.check(name+" creation",v1.size(),ngen);
s1 = a.testAdditionSV(v1); iret |= a.check(name+" addition",s1,a.Sum(),Dim*4);
sref1 = s1;
// test the io
double fsize = 0;
int ir = 0;
if (!testio) return iret;
std::string typeName = "ROOT::Math::"+name;
double estSize = ngen*8 * Dim + 10000;
double scale = 0.1 / std::numeric_limits<double>::epsilon();
fsize = a.testWrite(v1,typeName); iret |= a.check(name+" write",fsize,estSize,scale);
ir = a.testRead(v1); iret |= a.check(name+" read",ir,0);
s1 = a.testAdditionSV(v1); iret |= a.check(name+" after read",s1,sref1);
return iret;
}
template<int D1, int D2>
struct RepStd {
typedef typename ROOT::Math::MatRepStd<double,D1,D2> R;
static std::string name() {
return "ROOT::Math::MatRepStd<double," + Util::ToString(D1) + "," + Util::ToString(D2) + "> ";
}
static std::string name32() {
return "ROOT::Math::MatRepStd<Double32_t," + Util::ToString(D1) + "," + Util::ToString(D2) + "> ";
}
static std::string sname() { return ""; }
};
template<int D1>
struct RepSym {
typedef typename ROOT::Math::MatRepSym<double,D1> R;
static std::string name() {
return "ROOT::Math::MatRepSym<double," + Util::ToString(D1) + "> ";
}
static std::string name32() {
return "ROOT::Math::MatRepSym<Double32_t," + Util::ToString(D1) + "> ";
}
static std::string sname() { return "MatRepSym"; }
};
//--------------------------------------------------------------------------------------
// test of generic SMatrix
//--------------------------------------------------------------------------------------
template<int D1, int D2,class Rep >
int testSMatrix(int ngen, bool testio=false) {
// test the matrix if D2 is not equal to 1
int iret = 0;
typedef typename Rep::R R;
typedef SMatrix<double, D1, D2,R> SM;
// in the case of sym matrices SM::kSize is different than R::kSize
// need to use the R::kSize for Dim
const int Dim = R::kSize;
VectorTest<Dim> a(ngen);
a.genDataN();
std::vector<SM> v1;
v1.reserve(ngen);
double s1 = 0;
//double scale = 1;
double sref1 = 0;
std::string name0 = "SMatrix<double," + Util::ToString(D1) + "," + Util::ToString(D2);
std::string name = name0 + "," + Rep::sname() + ">";
a.testCreateSV(v1); iret |= a.check(name+" creation",v1.size(),ngen);
s1 = a.testAdditionSV(v1); iret |= a.check(name+" addition",s1,a.Sum(),Dim*4);
sref1 = s1;
// test the io
if (!testio) return iret;
double fsize = 0;
int ir = 0;
// the full name is needed for sym matrices
std::string typeName = "ROOT::Math::"+name0 + "," + Rep::name() + ">";
double estSize = ngen*8 * Dim + 10000;
double scale = 0.1 / std::numeric_limits<double>::epsilon();
fsize = a.testWrite(v1,typeName); iret |= a.check(name+" write",fsize,estSize,scale);
ir = a.testRead(v1); iret |= a.check(name+" read",ir,0);
s1 = a.testAdditionSV(v1); iret |= a.check(name+" after read",s1,sref1);
//std::cout << "File size = " << fsize << " estimated " << 8 * Dim * ngen << std::endl;
// test storing as Double32_t
// dictionay exist only for square matrices between 3 and 6
if (D1 != D2) return iret;
if (D1 < 3 || D1 > 6) return iret;
double fsize32 = 0;
name0 = "SMatrix<Double32_t," + Util::ToString(D1) + "," + Util::ToString(D2);
name = name0 + "," + Rep::sname() + ">";
typeName = "ROOT::Math::"+name0+ "," + Rep::name32() + ">";
estSize = ngen* 4 * Dim + 10000;
scale = 0.1 / std::numeric_limits<double>::epsilon();
fsize32 = a.testWrite(v1,typeName); iret |= a.check(name+" write",fsize32,estSize,scale);
ir = a.testRead(v1); iret |= a.check(name+" read",ir,0);
//we read back float (scale errors then)
s1 = a.testAdditionSV(v1); iret |= a.check(name+" after read",s1,sref1,1.E9);
return iret;
}
//--------------------------------------------------------------------------------------
// test of a track an object containing vector and matrices
//--------------------------------------------------------------------------------------
template<class T>
int testTrack(int ngen) {
int iret = 0;
const int Dim = T::kSize;
std::string name = T::Type();
//std::cout << "Dim is " << Dim << std::endl;
VectorTest<Dim> a(ngen);
a.genDataN();
std::vector<T> v1;
v1.reserve(ngen);
double s1 = 0;
//double scale = 1;
double sref1 = 0;
a.testCreateSV(v1); iret |= a.check(name+" creation",v1.size(),ngen);
s1 = a.testAdditionTR(v1); iret |= a.check(name+" addition",s1,a.Sum(),Dim*4);
sref1 = s1;
//std::cout << " reference sum is " << sref1 << std::endl;
double fsize = 0;
int ir = 0;
// the full name is needed for sym matrices
std::string typeName = name;
int wsize = 8;
if (T::IsD32() ) wsize = 4;
double estSize = ngen*wsize * Dim + 10000;
double scale = 0.2 / std::numeric_limits<double>::epsilon();
fsize = a.testWrite(v1,typeName); iret |= a.check(name+" write",fsize,estSize,scale);
ir = a.testRead(v1); iret |= a.check(name+" read",ir,0);
scale = 1;
if (T::IsD32() ) scale = 1.E9; // use float numbers
s1 = a.testAdditionTR(v1); iret |= a.check(name+" after read",s1,sref1,scale);
if (T::IsD32() ) {
// check at double precision type must fail otherwise Double's are stored
bool pdebug = debug; debug = false;
int ret = compare(" ",s1,sref1);
iret |= a.check(name+" Double32 test",ret,1);
debug = pdebug;
}
return iret;
}
int testGenVectors(int ngen,bool io) {
int iret = 0;
std::cout <<"******************************************************************************\n";
std::cout << "\tTest of Physics Vector (GenVector package)\n";
std::cout <<"******************************************************************************\n";
iret |= testVector<XYVector, Polar2DVector, 2>(ngen,io);
iret |= testVector<XYZVector, Polar3DVector, 3>(ngen,io);
iret |= testVector<XYZVector, RhoEtaPhiVector, 3>(ngen,io);
iret |= testVector<XYZVector, RhoZPhiVector, 3>(ngen,io);
iret |= testVector<XYZTVector, PtEtaPhiEVector, 4>(ngen,io);
iret |= testVector<XYZTVector, PtEtaPhiMVector, 4>(ngen,io);
iret |= testVector<XYZTVector, PxPyPzMVector, 4>(ngen,io);
return iret;
}
int testSMatrix(int ngen,bool io) {
int iret = 0;
std::cout <<"******************************************************************************\n";
std::cout << "\tTest of SMatrix package\n";
std::cout <<"******************************************************************************\n";
iret |= testVector34<3>(ngen,io);
iret |= testVector34<4>(ngen,io);
// test now matrices and vectors
iret |= testSVector<6>(ngen,io);
// asymetric matrices we have only 4x3 and 3x4
iret |= testSMatrix<3,4,RepStd<3,4> >(ngen,io);
iret |= testSMatrix<4,3,RepStd<4,3> >(ngen,io);
iret |= testSMatrix<3,3,RepStd<3,3> >(ngen,io);
// sym matrix
iret |= testSMatrix<5,5,RepSym<5> >(ngen,io);
return iret;
}
int testCompositeObj(int ngen) {
int iret = 0;
std::cout <<"******************************************************************************\n";
std::cout << "\tTest of a Composite Object (containing Vector's and Matrices)\n";
std::cout <<"******************************************************************************\n";
#ifndef USE_REFLEX
std::cout << "Test Using CINT library\n\n";
// put path relative to LD_LIBRARY_PATH
iret = gSystem->Load("../test/libTrackMathCoreDict");
if (iret !=0) {
// if not assume running from top ROOT dir (case of roottest)
iret = gSystem->Load("test/libTrackMathCoreDict");
if (iret !=0) {
std::cerr <<"Error Loading libTrackMathCoreDict" << std::endl;
return iret;
}
}
#else
std::cout << "Test Using Reflex library\n\n";
#ifdef DEBUG
ROOT::Cintex::Cintex::SetDebug(1);
#endif
ROOT::Cintex::Cintex::Enable();
iret = gSystem->Load("../test/libTrackMathCoreRflx");
if (iret !=0) {
// if not assume running from top ROOT dir (case of roottest)
iret = gSystem->Load("test/libTrackMathCoreRflx");
if (iret !=0) {
std::cerr <<"Error Loading libTrackMathCoreRflx" << std::endl;
return iret;
}
}
#endif
iret |= testTrack<TrackD>(ngen);
iret |= testTrack<TrackD32>(ngen);
iret |= testTrack<TrackErrD>(ngen);
iret |= testTrack<TrackErrD32>(ngen);
iret |= testTrack<VecTrack<TrackD> >(ngen);
iret |= testTrack<VecTrack<TrackErrD> >(ngen);
return iret;
}
#endif // endif ifndef __CINT__
int stressMathCore(double nscale = 1) {
int iret = 0;
#ifdef __CINT__
std::cout << "Test must be run in compile mode - use ACLIC to compile!!" << std::endl;
gSystem->Load("libMathCore");
gSystem->Load("libTree");
gROOT->ProcessLine(".L stressMathCore.cxx++");
return stressMathCore();
#endif
// iret |= gSystem->Load("libMathCore");
// iret |= gSystem->Load("libMathMore");
// if (iret !=0) return iret;
TBenchmark bm;
bm.Start("stressMathCore");
const int ntest = 10000;
int n = int(nscale*ntest);
//std::cout << "StressMathCore: test number n = " << n << std::endl;
iret |= testStatFunctions(n/10);
bool io = true;
iret |= gSystem->Load("libSmatrix");
if (iret !=0) {
std::cerr <<"Error Loading libSmatrix" << std::endl;
io = false;
}
iret |= testGenVectors(n,io);
iret |= testSMatrix(n,io);
if (io) iret |= testCompositeObj(n);
bm.Stop("stressMathCore");
std::cout <<"******************************************************************************\n";
bm.Print("stressMathCore");
const double reftime = 7.1; // needs to be updated // ref time on pcbrun4
double rootmarks = 860 * reftime / bm.GetCpuTime("stressMathCore");
std::cout << " ROOTMARKS = " << rootmarks << " ROOT version: " << gROOT->GetVersion() << "\t"
<< gROOT->GetSvnBranch() << "@" << gROOT->GetSvnRevision() << std::endl;
std::cout <<"*******************************************************************************\n";
if (iret !=0) std::cerr << "stressMathCore Test Failed !!" << std::endl;
return iret;
}
int main(int argc,const char *argv[]) {
double nscale = 1;
if (argc > 1) {
nscale = atof(argv[1]);
//nscale = std::pow(10.0,double(scale));
}
return stressMathCore(nscale);
}
