How to call Mathematica ( or alternatives) function evaluation for specific pdfs

[RENATO_QUAGLIANI]

Dear experts,
I came across a problem using the 2F1 function to evaluate analytical integrals of a function which I currently have implemented as a RooGenericPdf.
After a bit of digging it turned out that inside Mathematica the integrals I want to compute are able to be evaluated without prompting infinities or nan errors from gsl.
I had a look around and it seems like that the hypergeometric Gaussian function from gsl doesn’t support some of the tricks to allow the function to be defined when X<-1.
Therefore I was planning to test whether cephes is able to do so but ideally I would like to use the Mathematica function call and the 2F1 implementation. I believe this would be doable if Mathematica would be open source but maybe some others had similar issue and have been able to circumvent it somehow or maybe there is some trick to use which I am not aware of.

Thanks in advance for any help

EDIT:
From what I read it seems possible to compile some C++ code which callback some mathematica code, however a guideline would be good to have ,also is that the linking possible on LCG software stack or do I need to install it as a standalone licensed software?

---

Please read tips for efficient and successful posting and posting code

ROOT Version: Not Provided
Platform: Not Provided
Compiler: Not Provided

---

[RENATO_QUAGLIANI]

To give a bit more context : I attach here the function i am trying to analytically integrate.
According to mathematica (not an expert) , the definite integral is defined whenever the exp of the last term in the 2F1 function is >=-1 which for real numbers Is always the case.

Therefore , the way i coded up this is to check for |z| < 1 to use the standard call of hyperg, while for z < -1 to call the transformed one ( taken from Calculating Hypergeometric Function 2F1 for |z|>1 | Physics Forums )

Double_t RooExpTurnOn::IntegralVal( Double_t xVal) const {
  Double_t Integrand = TMath::Exp( xVal * bComb)/bComb;
  Double_t A = getA();
  Double_t B = getB();
  Double_t C = getC();
  Double_t Z = getZ(xVal);
  Double_t Factor = 0;
  if( abs(Z) <1){
    Factor = -1. + ROOT::Math::hyperg( A, B, C, -Z);
  }else if( -Z < -1 ){
    Factor = -1. + hyperg_z_GT1( A, B, C, -Z);
  }else if( -Z > 1){
    abort();
  }
  return Integrand * Factor;
}
Double_t RooExpTurnOn::hyperg_z_GT1(double a, double b, double c, double z) const { 
  //calculates 2F1 for z < -1 
  double coef1,coef2; 
  coef1=ROOT::Math::tgamma(c)*ROOT::Math::tgamma(b-a)*pow(1.-z,-a)/(ROOT::Math::tgamma(b)*ROOT::Math::tgamma(c-a)); 
  coef2=ROOT::Math::tgamma(c)*ROOT::Math::tgamma(a-b)*pow(1.-z,-b)/(ROOT::Math::tgamma(a)*ROOT::Math::tgamma(c-b)); 
  return coef1*ROOT::Math::hyperg(a,c-b,a-b+1.,1./(1.-z))+coef2*ROOT::Math::hyperg(b,c-a,b-a+1.,1./(1.-z)); 
}

This seems to do the trick.
Being not an expert of GSL , neither a mathematician , i wonder if anyone have some comments, as far as i understood gsl::hyperg doesn’t support the case z < -1 , but from what i can see, it can be easily implemented.

[etejedor]

Hello,

@moneta or @jonas should be able to answer soon, since the end of the Christmas break is close!

[RENATO_QUAGLIANI]

Thanks @etejedor , i just want to give an heads-up from what i find out the last 2-3 days :
Acoording to this paper (if i understood correctly) , in what i posted before i am using only one possible case of transformation of argument for the function. I looked a bit around on other packages in R, Fortran etc… and i ended up concluding that one could code up in ROOT::Math::hypergeo a specific 2F1(a,b,c,z) call which branches into the different equations depending of Re(z) regions ( since the ROOT::Math::hyperg() uses only double one can take profit of the way this work in hypergeo from R. See figure 1. and equations 15.3.4-15.3.9 .

RJ-2015-022.pdf

Not sure if this sounds like a good solution and possibly provide a nice extension to avoid errors on domain calls from the ROOT::Math::hyperg call

[moneta]

Hi,

I need to investigate in detail, but it seems to be a possible improvement that can be added to the hypergeometric functions we have in ROOT. If you want to provide a PR for this it will be great.
One thing I would recommend is to test these new branch of the function (for |z| > 1 ) by comparing the results with Mathematica.

Lorenzo

[RENATO_QUAGLIANI]

Ok, i see if i find some time , for the use case i currently have, mathematica is able to compute things, and no assumptions are done behind the scenes. I think whatever i could manage to put in the Pull request, might break down in some corner cases and dispatch all the transformation might be unsuitable.
So i think it could become more a catch “the singularity” cases effort than anything else. Maybe there are some GSL developers which one could ask to ? In principle all possible branchings are implemented in the arb library and in R from the hypergeo package.