root/test/TFormulaVecTests.h at 096f2dfa7d105e1817cecaf26439a845f46159e2 · guitargeek/root · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
#include <TF1.h>
#include <TF2.h>
#include <iostream>
#include <TString.h>

bool verbose  = true;
//#define VEC_IN_CTOR

bool CheckValues(const TString & testName, double r, double r0) {

   bool ret =  TMath::AreEqualAbs(r,r0,1.E-12);
   if (!ret) {
      std::cout << "Test failed for " << testName << "  : " <<  r << "   " << r0 <<  std::endl;
   }
   else
      if (verbose) std::cout << testName << "\t ok\n";

   return ret;
}


typedef double ( * FreeFunc1D) (double );
bool testVec1D(TF1 * f1, const TString & formula, FreeFunc1D func, double x ) {

   auto y = f1->Eval(x);
   auto y0 = func(x);

   bool ret = CheckValues(formula, y, y0);

   // check passing double_v interface
#ifdef R__HAS_VECCORE
   ROOT::Double_v vx = x;
   ROOT::Double_v vy = f1->EvalPar(&vx, nullptr);
   double y2 = vecCore::Get(vy,0);
   ret &= CheckValues(formula+TString("_v"), y2, y0);
#endif

   return ret;
}
bool testVec1D(const TString & formula, FreeFunc1D func, double x = 1.) {

   // test first by vectorizing in the constructor
   auto f1 = new TF1("f",formula,0,1,"VEC");
   bool ret = testVec1D(f1,formula,func,x);
   // test by vectorizing afterwards calling SetVectorized
   auto f2 = new TF1("f2",formula);
   f2->SetVectorized(true);
   //std::cout << "set vectprization" << std::endl;
   //f2->Print("V");
   ret &=  testVec1D(f2,formula,func,x);
   // test by removing vectorization
   f2->SetVectorized(false);
   ret &=  testVec1D(f1,formula,func,x);
   return ret;
}


typedef double ( * FreeFunc2D) (double, double );
bool testVec2D(TF2 * f1, const TString & formula, FreeFunc2D func, double x, double y ) {

   auto r = f1->Eval(x,y);
   auto r0 = func(x,y);

   bool ret = CheckValues(formula,r,r0);

   // check passing double_v interface
#ifdef R__HAS_VECCORE
   ROOT::Double_v vx[2] = { x, y};
   ROOT::Double_v vy = f1->EvalPar(vx, nullptr);
   double r2 = vecCore::Get(vy,0);
   ret &= CheckValues(formula+TString("_v"), r2, r0);
#endif

   return ret;
}

bool testVec2D(const TString & formula, FreeFunc2D func, double x = 1., double y = 1.) {

   auto f1 = new TF2("f",formula,0,1,0,1,"VEC");
   bool ret = testVec2D(f1, formula, func, x, y);

   auto f2 = new TF2("f",formula);
   f2->SetVectorized(true);
   ret &=  testVec2D(f2,formula,func,x,y);
   return ret;
}


double constant_function(double ) { return 3; }


bool testVecFormula() {

   bool ok = true;

   // Note that we have to disable several tests when using the Vc backend
   // because of linker errors. This is because Vc is always built as a static
   // library, and linked in MathCore, but the interpreter cannot lookup
   // missing symbols for Vc math functions if the formula requires them and
   // they are not already linked into MathCore by chance.
   // See also:
   // https://its.cern.ch/jira/browse/ROOT-10614
   // https://github.com/root-project/root/pull/20769

   ok &= testVec1D("3+[0]",constant_function, 3.333);
   ok &= testVec1D("3",constant_function, 3.333);
#ifndef R__HAS_VC
   ok &= testVec1D("sin(x)",std::sin,1);
   ok &= testVec1D("cos(x)",std::cos,1);
#endif
   ok &= testVec1D("exp(x)",std::exp,1);
   ok &= testVec1D("log(x)",std::log,2);
   ok &= testVec1D("log10(x)",TMath::Log10,2);
#ifndef R__HAS_VC
   ok &= testVec1D("tan(x)",std::tan,1);
   ok &= testVec1D("sinh(x)",std::sinh,1);
   ok &= testVec1D("cosh(x)",std::cosh,1);
   ok &= testVec1D("tanh(x)",std::tanh,1);
   ok &= testVec1D("asin(x)",std::asin,.1);
   ok &= testVec1D("acos(x)",std::acos,.1);
#endif
   ok &= testVec1D("atan(x)",std::atan,.1);
   ok &= testVec1D("sqrt(x)",std::sqrt,2);
   ok &= testVec1D("abs(x)",std::abs,-1);
   ok &= testVec2D("pow(x,y)",std::pow,2,3);
#ifndef R__HAS_VC
   ok &= testVec2D("min(x,y)",TMath::Min,2,3);
   ok &= testVec2D("max(x,y)",TMath::Max,2,3);
#endif
   ok &= testVec2D("atan2(x,y)",TMath::ATan2,2,3);

   return ok;
}