#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;
ok &= testVec1D("3+[0]",constant_function, 3.333);
ok &= testVec1D("3",constant_function, 3.333);
ok &= testVec1D("sin(x)",std::sin,1);
ok &= testVec1D("cos(x)",std::cos,1);
ok &= testVec1D("exp(x)",std::exp,1);
ok &= testVec1D("log(x)",std::log,2);
ok &= testVec1D("log10(x)",TMath::Log10,2);
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);
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);
ok &= testVec2D("min(x,y)",TMath::Min,2,3);
ok &= testVec2D("max(x,y)",TMath::Max,2,3);
ok &= testVec2D("atan2(x,y)",TMath::ATan2,2,3);
return ok;
}