HepLib
QuadMP.cpp
Go to the documentation of this file.
1 
6 #include <math.h>
7 #include <complex>
8 extern "C" {
9 #include <quadmath.h>
10 }
11 #include "mpreal.h"
12 #include "Lib3_GK.h"
13 #include "SD.h"
14 
15 /* error return codes */
16 #define SUCCESS 0
17 #define FAILURE 1
18 
19 using namespace std;
20 namespace {
21  typedef mpfr::mpreal mpREAL;
22  typedef const mpREAL & mpREAL_t;
23  typedef complex<mpREAL> mpCOMPLEX;
24  typedef std::function<int(unsigned yn, mpREAL *y, mpREAL *e, const mpREAL & x, void *fdata)> f1Type;
25  typedef const f1Type & f1Type_t;
26  typedef std::function<int(unsigned yn, mpREAL *y, mpREAL *e, unsigned xdim, const mpREAL *x, void *fdata)> fnType;
27  typedef const fnType & fnType_t;
28  typedef void (*PrintHookerType) (mpREAL *, mpREAL *, size_t *, void *);
29 }
30 
31 extern mpREAL mpPi;
32 extern mpREAL mpEuler;
33 extern mpCOMPLEX mpiEpsilon;
34 
35 namespace {
36  size_t _nGK_ = 10000;
37  size_t _mGK_ = 10; // sets (2m+1)-point Gauss-Kronrod
38 
39  int QuadPack1(unsigned yn, mpREAL *oval, mpREAL *oerr, f1Type_t f, mpREAL_t epsabs, PrintHookerType PrintHooker, void *fdata) {
40  GK gk(_nGK_, _mGK_);
41  Function F(yn,f,fdata);
42  try {
43  return gk.QAG(F, 0, 1, epsabs, 0, oval, oerr, PrintHooker);
44  } catch (const char* reason) {
45  cout << reason << endl;
46  throw reason;
47  }
48  return SUCCESS;
49  }
50 
51  int QuadPackN(unsigned yn, mpREAL *oval, mpREAL *oerr, unsigned xdim, fnType_t f, mpREAL_t eps, PrintHookerType PrintHooker, void *fdata) {
52  if(xdim==1) {
53  auto f1 = [f,eps](unsigned yn, mpREAL *y, mpREAL *e, mpREAL_t x, void *fdata)->int {
54  mpREAL xs[1];
55  xs[0] = x;
56  return f(yn,y,e,1,xs,fdata);
57  };
58  return QuadPack1(yn, oval, oerr, f1, eps, PrintHooker, fdata);
59  }
60 
61  auto f1 = [f,eps,xdim](unsigned yn, mpREAL *y, mpREAL *e, mpREAL_t x, void *fdata)->int {
62  auto f2 =[f,x](unsigned yn, mpREAL *y1, mpREAL *e1, unsigned xdim1, const mpREAL *xs1, void *fdata1)->int {
63  mpREAL xs[xdim1+1];
64  if(true) { // insert first
65  xs[0] = x;
66  for(int i=0; i<xdim1; i++) xs[i+1] = xs1[i];
67  } else { // insert last
68  for(int i=0; i<xdim1; i++) xs[i] = xs1[i];
69  xs[xdim1] = x;
70  }
71  return f(yn,y1,e1,xdim1+1,xs,fdata1);
72  };
73  return QuadPackN(yn,y,e,xdim-1,f2,eps/xdim,NULL,fdata);
74  };
75  return QuadPack1(yn, oval, oerr, f1, eps, PrintHooker, fdata);
76  }
77 
78 }
79 
80 namespace HepLib::SD {
81 
82 /*-----------------------------------------------------*/
83 // QuadMP Classes
84 /*-----------------------------------------------------*/
85 
86 ex QuadMP::mp2ex(const mpREAL & num) {
87  ostringstream oss;
88  oss.precision(MPDigits);
89  oss << num;
90  string sn = oss.str();
91  numeric ret(sn.c_str());
92  return ret;
93 }
94 
95 int QuadMP::Wrapper(unsigned ydim, mpREAL *y, mpREAL *e, unsigned xdim, const mpREAL *x, void *fdata) {
96  auto self = (QuadMP*)fdata;
97  self->IntegrandMP(xdim, x, ydim, y, self->mpParameter, self->mpLambda);
98  // Check NaN/Inf
99  bool ok = true;
100  for(int j=0; j<ydim; j++) {
101  if(!isfinite(y[j])) { ok = false; break; }
102  }
103  if(!ok) {
104  mpfr::mpreal::set_default_prec(mpfr::digits2bits(self->MPDigits*10));
105  self->IntegrandMP(xdim, x, ydim, y, self->mpParameter, self->mpLambda);
106  mpfr::mpreal::set_default_prec(mpfr::digits2bits(self->MPDigits));
107  }
108  for(int j=0; j<ydim; j++) e[j] = 0; // assume error can be ignored
109  return SUCCESS;
110 }
111 
112 void QuadMP::DefaultPrintHooker(mpREAL* result, mpREAL* epsabs, size_t * nrun, void *fdata) {
113  auto self = (QuadMP*)fdata;
114  if(*nrun == self->nGK + 1979) return;
115  if(self->RunTime>0) {
116  auto cur_timer = time(NULL);
117  auto used_time = difftime(cur_timer,self->StartTimer);
118  if(used_time>self->RunTime) {
119  self->NEval = *nrun;
120  *nrun = self->nGK + 1979;
121  if(Verbose>10) cout << WarnColor << " Exit with Run out of Time: " << used_time << RESET << endl;
122  return;
123  }
124  }
125  if(Verbose>10) {
126  auto r0 = result[0];
127  auto r1 = result[1];
128  auto e0 = epsabs[0].toString(3);
129  auto e1 = epsabs[1].toString(3);
130  cout << " L: " << (*nrun) << ", ";
131  if(self->ReIm==3 || self->ReIm==1) cout << "[" << r0 << ", " << e0 << "]";
132  if(self->ReIm==3 || self->ReIm==2) cout << "+I*[" << r1 << ", " << e1 << "]";
133  cout << endl;
134  }
135  self->NEval = *nrun;
136 
137  bool rExit = (epsabs[0] < self->EpsAbs+1E-50Q) || (epsabs[0] < fabs(result[0])*self->EpsRel+1E-50Q);
138  bool iExit = (epsabs[1] < self->EpsAbs+1E-50Q) || (epsabs[1] < fabs(result[1])*self->EpsRel+1E-50Q);
139  if(rExit && iExit) {
140  *nrun = self->nGK + 1979;
141  return;
142  }
143 
144  auto pid = getpid();
145  ostringstream fn;
146  fn << pid << ".int.done";
147  if(file_exists(fn.str().c_str())) {
148  self->NEval = *nrun;
149  *nrun = self->nGK + 1979;
150  ostringstream cmd;
151  cmd << "rm " << fn.str();
152  auto rc = system(cmd.str().c_str());
153  if(Verbose>10) cout << " Exit: " << fn.str() << endl;
154  }
155 }
156 
157 ex QuadMP::Integrate(size_t n) {
158  if(mpfr_buildopt_tls_p()<=0) throw Error("Integrate: mpfr_buildopt_tls_p()<=0.");
159  mpfr_free_cache();
160  mpfr::mpreal::set_default_prec(mpfr::digits2bits(MPDigits));
161  mpPi = mpfr::const_pi();
162  mpEuler = mpfr::const_euler();
163  mpiEpsilon = complex<mpREAL>(0,mpfr::machine_epsilon()*100);
164 
165  unsigned int xdim = XDim;
166  unsigned int ydim = 2;
167  mpREAL result[ydim], estabs[ydim];
168 
169  NEval = 0;
170  _mGK_ = mGK;
171  _nGK_ = n==0 ? nGK : n;
172  StartTimer = time(NULL);
173  auto nok = QuadPackN(ydim, result, estabs, xdim, Wrapper, EpsAbs, n==0 ? PrintHooker : NULL, this);
174 
175  if(nok) {
176  mpREAL abs_res = sqrt(result[0]*result[0]+result[1]*result[1]);
177  mpREAL abs_est = sqrt(estabs[0]*estabs[0]+estabs[1]*estabs[1]);
178  mpREAL mpfr_eps = 10*mpfr::machine_epsilon();
179  if( (abs_res < mpfr_eps) && (abs_est < mpfr_eps) ) {
180  cout << ErrColor << "QuadMP Failed with 0 result returned!" << RESET << endl;
181  return NaN;
182  }
183  }
184 
185  ex FResult = 0;
186  if(isnan(result[0]) || isnan(result[1])) FResult += NaN;
187  else {
188  try{
189  FResult += VE(mp2ex(result[0]), mp2ex(estabs[0]));
190  FResult += VE(mp2ex(result[1]), mp2ex(estabs[1])) * I;
191  } catch(...) {
192  FResult += NaN;
193  }
194  }
195 
196  return FResult;
197 }
198 
199 }
#define RESET
Definition: BASIC.h:79
complex< mpREAL > mpCOMPLEX
Definition: HCubature.cpp:20
mpfr::mpreal mpREAL
Definition: HCubature.cpp:19
bool file_exists(const char *fn)
Definition: Process.cpp:9
mpREAL mpEuler
#define SUCCESS
Definition: QuadMP.cpp:16
mpCOMPLEX mpiEpsilon
mpREAL mpPi
SecDec header file.
class used to wrap error message
Definition: BASIC.h:242
numerical integrator using TanhSinhMP
Definition: SD.h:249
namespace for Numerical integration with Sector Decomposition method
Definition: AsyMB.cpp:10
mpfr::mpreal mpREAL
Definition: SD.h:125
const char * ErrColor
Definition: BASIC.cpp:246
const Symbol eps
const char * WarnColor
Definition: BASIC.cpp:247
int Verbose
Definition: Init.cpp:139
const Symbol NaN