HepLib
TanhSinhMP.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 <omp.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 ydim, mpREAL *y, mpREAL *e, const mpREAL & x, void *fdata)> f1Type;
25  typedef const f1Type & f1Type_t;
26  typedef std::function<int(unsigned ydim, 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 extern mpREAL mpPi;
31 extern mpREAL mpEuler;
32 
33 namespace {
34  int CPUCORES = 8;
35  size_t kmax = 10;
36  // int f from a to b, parallel version
37  int TanhSinh1(mpREAL *oval, mpREAL *oerr, f1Type_t f, unsigned ydim, mpREAL_t epsrel, PrintHookerType PrintHooker, void *fdata) {
38  mpREAL a=0, b=1; // integration domain (a,b)
39  mpREAL c = (a+b)/2; // gamma
40  mpREAL d = (b-a)/2; // sigma
41  const mpREAL tol = 10*epsrel;
42  mpREAL s[ydim], e_s[ydim], val[ydim], err[ydim], v[ydim], h = 2;
43  if(f(ydim, s, e_s, c, fdata)) return FAILURE;
44  size_t k = 0;
45  while(k<=kmax) {
46  mpREAL p[ydim], fp[ydim], fm[ydim], q, t, eh;
47  for(int j=0; j<ydim; j++) p[j] = fp[j] = fm[j] = 0;
48  h /= 2;
49  eh = exp(h);
50  t = eh;
51  if (k > 0) eh *= eh;
52  bool parallel = true;
53  if(omp_in_parallel() && ( !omp_get_nested() || omp_get_active_level()>1 )) parallel = false;
54  if(parallel) { // Parallel
55  while(true) {
56  int total = CPUCORES/2;
57  mpREAL xs[total], ws[total], fps[total*ydim], e_fps[total*ydim], fms[total*ydim], e_fms[total*ydim];
58  int RC1[total], RC2[total];
59  for(int i=0; i<total; i++) {
60  mpREAL u = exp(1/t-t);
61  mpREAL r = 2*u/(1+u);
62  ws[i] = (t+1/t)*r/(1+u);
63  xs[i] = d*r;
64  t *= eh;
65  RC1[i] = RC2[i] = 0;
66  }
67  auto prec = mpfr::mpreal::get_default_prec();
68  auto rnd = mpfr::mpreal::get_default_rnd();
69  #pragma omp parallel for
70  for(int ii=0; ii<2*total; ii++) {
71  int i = ii/2, i2 = ii%2;
72  mpfr::mpreal::set_default_prec(prec);
73  mpfr::mpreal::set_default_rnd(rnd);
74  mpREAL_t x = xs[i];
75  if(i2==0) { if (a+x > a) RC1[i] = f(ydim,fps+ydim*i,e_fps+ydim*i,a+x,fdata); }
76  else if(i2==1) { if (b-x < b) RC2[i] = f(ydim,fms+ydim*i,e_fms+ydim*i,b-x,fdata); }
77  mpfr_free_cache();
78  }
79  bool ok;
80  for(int i=0; i<total; i++) {
81  if(RC1[i]!=0 || RC2[i]!=0) return FAILURE;
82  mpREAL_t x = xs[i];
83  ok = true;
84  for(int j=0; j<ydim; j++) {
85  if (a+x > a) if (isfinite(fps[i])) {
86  fp[j] = fps[i*ydim+j];
87  if(e_s[j] < e_fps[i*ydim+j]) e_s[j] = e_fps[i*ydim+j];
88  }
89  if (b-x < b) if (isfinite(fms[i])) {
90  fm[j] = fms[i*ydim+j];
91  if(e_s[j] < e_fms[i*ydim+j]) e_s[j] = e_fms[i*ydim+j];
92  }
93  q = ws[i]*(fp[j]+fm[j]);
94  p[j] += q;
95  ok = ok && (fabs(q)<=epsrel*fabs(p[j]));
96  }
97  if(ok) break;
98  }
99  if(ok) break;
100  }
101  } else { // Non-Parallel
102  while(true) {
103  mpREAL u = exp(1/t-t);
104  mpREAL r = 2*u/(1+u);
105  mpREAL w = (t+1/t)*r/(1+u);
106  mpREAL x = d*r;
107  t *= eh;
108  mpREAL y[ydim], e_y[ydim];
109  if (a+x > a) {
110  if(f(ydim,y,e_y,a+x,fdata)) return FAILURE;
111  for(int j=0; j<ydim; j++) if (isfinite(y[j])) {
112  fp[j] = y[j];
113  if(e_s[j] < e_y[j]) e_s[j] = e_y[j];
114  }
115  }
116  if (b-x < b) {
117  if(f(ydim,y,e_y,b-x,fdata)) return FAILURE;
118  for(int j=0; j<ydim; j++) if (isfinite(y[j])) {
119  fm[j] = y[j];
120  if(e_s[j] < e_y[j]) e_s[j] = e_y[j];
121  }
122  }
123  bool ok = true;
124  for(int j=0; j<ydim; j++) {
125  q = w*(fp[j]+fm[j]);
126  p[j] += q;
127  ok = ok && (fabs(q)<=epsrel*fabs(p[j]));
128  }
129  if(ok) break;
130  }
131  }
132  for(int j=0; j<ydim; j++) {
133  v[j] = s[j] - p[j];
134  s[j] += p[j];
135  val[j] = d*h*s[j];
136  err[j] = fabs(d*h*v[j]) + (b-a)*e_s[j];
137  }
138  if(k > 10) {
139  cout << RED << "Large Level: " << RESET << k << endl;
140  for(int j=0; j<ydim; j++) {
141  cout << val[j] << " +- " << err[j].toString(3) << endl;
142  }
143  cout << endl;
144  }
145  size_t kk = k;
146  if(PrintHooker) PrintHooker(val, err, &kk, fdata);
147  if(kk!=k) {
148  for(int j=0; j<ydim; j++) {
149  oval[j] = val[j];
150  oerr[j] = err[j];
151  }
152  return SUCCESS;
153  }
154  bool ok = true;
155  for(int j=0; j<ydim; j++) {
156  ok = ok && (fabs(v[j])<=tol*fabs(s[j]));
157  if(!ok) break;
158  }
159  if(ok) break;
160  ++k;
161  }
162  for(int j=0; j<ydim; j++) {
163  oval[j] = val[j];
164  oerr[j] = err[j];
165  }
166  return SUCCESS;
167  }
168 
169  int TanhSinhN(mpREAL *val, mpREAL *err, fnType_t f, unsigned xdim, unsigned ydim, mpREAL_t eps, PrintHookerType PrintHooker, void *fdata) {
170  if(xdim==1) {
171  auto f1 = [f](unsigned ydim, mpREAL *y, mpREAL *e, mpREAL_t x, void *fdata)->int {
172  mpREAL xs[1];
173  xs[0] = x;
174  if(f(ydim,y,e,1,xs,fdata)) return FAILURE;
175  return SUCCESS;
176  };
177  return TanhSinh1(val,err,f1,ydim,eps,PrintHooker,fdata);
178  } else {
179  auto f1 = [f,xdim,eps](unsigned ydim, mpREAL *y, mpREAL *e, mpREAL_t x, void *fdata)->int {
180  auto f2 =[f,x](unsigned ydim, mpREAL *y1, mpREAL *e1, unsigned xdim1, const mpREAL *xs1, void *fdata)->int {
181  mpREAL xs[xdim1+1];
182  if(true) { // insert first
183  xs[0] = x;
184  for(int i=1; i<=xdim1; i++) xs[i] = xs1[i-1];
185  } else { // insert last
186  for(int i=0; i<xdim1; i++) xs[i] = xs1[i];
187  xs[xdim1] = x;
188  }
189  if(f(ydim,y1,e1,xdim1+1,xs,fdata)) return FAILURE;
190  return SUCCESS;
191  };
192  return TanhSinhN(y,e,f2,xdim-1,ydim,eps/xdim,NULL,fdata);
193  };
194  return TanhSinh1(val,err,f1,ydim,eps,PrintHooker,fdata);
195  }
196  }
197 }
198 
199 namespace HepLib::SD {
200 
201 /*-----------------------------------------------------*/
202 // TanhSinhMP Classes
203 /*-----------------------------------------------------*/
204 
205 TanhSinhMP::TanhSinhMP(size_t k) { K = k; }
206 
207 ex TanhSinhMP::mp2ex(mpREAL_t num) {
208  ostringstream oss;
209  oss.precision(MPDigits);
210  oss << num;
211  string sn = oss.str();
212  numeric ret(sn.c_str());
213  return ret;
214 }
215 
216 int TanhSinhMP::Wrapper(unsigned ydim, mpREAL *y, mpREAL *e, unsigned xdim, const mpREAL *x, void *fdata) {
217  for(int j=0; j<ydim; j++) e[j] = 0; // assume error can be ignored
218  auto self = (TanhSinhMP*)fdata;
219  self->IntegrandMP(xdim, x, ydim, y, self->mpParameter, self->mpLambda);
220  bool ok = true;
221  for(int j=0; j<ydim; j++) {
222  if(!isfinite(y[j])) { ok = false; break; }
223  }
224  if(!ok) {
225  mpfr::mpreal::set_default_prec(mpfr::digits2bits(self->MPDigits*10));
226  self->IntegrandMP(xdim, x, ydim, y, self->mpParameter, self->mpLambda);
227  mpfr::mpreal::set_default_prec(mpfr::digits2bits(self->MPDigits));
228  }
229  return SUCCESS;
230 }
231 
232 void TanhSinhMP::DefaultPrintHooker(mpREAL* result, mpREAL* epsabs, size_t * nrun, void *fdata) {
233  auto self = (TanhSinhMP*)fdata;
234  if(*nrun == self->K + 1979) return;
235  if(self->RunTime>0) {
236  auto cur_timer = time(NULL);
237  auto used_time = difftime(cur_timer,self->StartTimer);
238  if(used_time>self->RunTime) {
239  self->NEval = *nrun;
240  *nrun = self->K + 1979;
241  if(Verbose>10) cout << WarnColor << " Exit with Run out of Time: " << used_time << RESET << endl;
242  return;
243  }
244  }
245  if(Verbose>10 && self->K>0) {
246  auto r0 = result[0];
247  auto r1 = result[1];
248  auto e0 = epsabs[0].toString(3);
249  auto e1 = epsabs[1].toString(3);
250  cout << " L: " << (*nrun) << ", ";
251  if(self->ReIm==3 || self->ReIm==1) cout << "[" << r0 << ", " << e0 << "]";
252  if(self->ReIm==3 || self->ReIm==2) cout << "+I*[" << r1 << ", " << e1 << "]";
253  cout << endl;
254  }
255  self->NEval = *nrun;
256 
257  bool rExit = (epsabs[0] < self->EpsAbs+1E-50Q) || (epsabs[0] < fabs(result[0])*self->EpsRel+1E-50Q);
258  bool iExit = (epsabs[1] < self->EpsAbs+1E-50Q) || (epsabs[1] < fabs(result[1])*self->EpsRel+1E-50Q);
259  if(rExit && iExit) {
260  *nrun = self->K + 1979;
261  return;
262  }
263 
264  auto pid = getpid();
265  ostringstream fn;
266  fn << pid << ".int.done";
267  if(file_exists(fn.str().c_str())) {
268  self->NEval = *nrun;
269  *nrun = self->K + 1979;
270  ostringstream cmd;
271  cmd << "rm " << fn.str();
272  system(cmd.str().c_str());
273  if(Verbose>10) cout << " Exit: " << fn.str() << endl;
274  }
275 }
276 
277 ex TanhSinhMP::Integrate(size_t n) {
278  kmax = n==0 ? K : n;
279  CPUCORES = omp_get_num_procs();
280  if(mpfr_buildopt_tls_p()<=0) throw Error("Integrate: mpfr_buildopt_tls_p()<=0.");
281  mpfr_free_cache();
282  mpfr::mpreal::set_default_prec(mpfr::digits2bits(MPDigits));
283  mpPi = mpfr::const_pi();
284  mpEuler = mpfr::const_euler();
285 
286  unsigned int xdim = XDim;
287  unsigned int ydim = 2;
288  mpREAL result[ydim], estabs[ydim];
289  NEval = 0;
290  mpREAL epsrel = 1E-2;
291  int nok = TanhSinhN(result, estabs, Wrapper, xdim, ydim, epsrel, NULL, this);
292  epsrel = min(EpsAbs/(fabs(result[0])+EpsAbs), EpsAbs/(fabs(result[1])+EpsAbs));
293  if(epsrel < 1E-2) {
294  StartTimer = time(NULL);
295  nok = TanhSinhN(result, estabs, Wrapper, xdim, ydim, epsrel, n==0 ? PrintHooker : NULL, this);
296  }
297 
298  if(nok) {
299  mpREAL abs_res = sqrt(result[0]*result[0]+result[1]*result[1]);
300  mpREAL abs_est = sqrt(estabs[0]*estabs[0]+estabs[1]*estabs[1]);
301  mpREAL mpfr_eps = 10*mpfr::machine_epsilon();
302  if( (abs_res < mpfr_eps) && (abs_est < mpfr_eps) ) {
303  cout << ErrColor << "TanhSinhMP Failed with 0 result returned!" << RESET << endl;
304  return NaN;
305  }
306  }
307 
308  ex FResult = 0;
309  if(!isfinite(result[0]) || !isfinite(result[1])) FResult += NaN;
310  else {
311  try{
312  FResult += VE(mp2ex(result[0]), mp2ex(estabs[0]));
313  FResult += VE(mp2ex(result[1]), mp2ex(estabs[1])) * I;
314  } catch(...) {
315  FResult += NaN;
316  }
317  }
318 
319  return FResult;
320 }
321 
322 }
int * a
Definition: Functions.cpp:234
#define RED
Definition: BASIC.h:81
#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
SecDec header file.
#define FAILURE
Definition: TanhSinhMP.cpp:17
mpREAL mpEuler
#define SUCCESS
Definition: TanhSinhMP.cpp:16
mpREAL mpPi
class used to wrap error message
Definition: BASIC.h:242
numerical integrator using TanhSinhMP
Definition: SD.h:233
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
ex w
Definition: Init.cpp:90
const Symbol d
int Verbose
Definition: Init.cpp:139
const Symbol NaN