HepLib
HCubature.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 "SD.h"
13 
14 using namespace std;
15 typedef __float128 qREAL;
16 typedef __complex128 qCOMPLEX;
17 typedef long double dREAL;
18 typedef complex<dREAL> dCOMPLEX;
19 typedef mpfr::mpreal mpREAL;
20 typedef complex<mpREAL> mpCOMPLEX;
21 
22 extern const qCOMPLEX qiEpsilon;
23 extern mpREAL mpPi;
24 extern mpREAL mpEuler;
25 extern mpCOMPLEX mpiEpsilon;
26 
27 #include "Lib3_HCubature.h"
28 namespace HepLib::SD {
29 
30 /*-----------------------------------------------------*/
31 // HCubature Classes
32 /*-----------------------------------------------------*/
33 
34  int HCubature::Wrapper(unsigned int xdim, size_t npts, const qREAL *x, void *fdata, unsigned int ydim, qREAL *y) {
35  auto self = (HCubature*)fdata;
36  bool NaNQ = false;
37 
38  unsigned int nthreads = self->Threads>0 ? self->Threads : omp_get_num_procs();
39  #pragma omp parallel for num_threads(nthreads) schedule(dynamic, 1)
40  for(int i=0; i<npts; i++) {
41  mpfr_free_cache();
42  mpfr::mpreal::set_default_prec(mpfr::digits2bits(self->MPDigits));
43  int iDQMP = self->inDQMP(x+i*xdim);
44  if( (self->IntegrandMP!=NULL) && (self->DQMP>2 || iDQMP>2) ) {
45  mpREAL mpx[xdim], mpy[ydim];
46  for(int j=0; j<xdim; j++) mpx[j] = x[i*xdim+j];
47  self->IntegrandMP(xdim, mpx, ydim, mpy, self->mpParameter, self->mpLambda);
48  for(int j=0; j<ydim; j++) y[i*ydim+j] = mpy[j].toFloat128();
49  } else if(self->DQMP>1 || iDQMP>1) {
50  self->IntegrandQ(xdim, x+i*xdim, ydim, y+i*ydim, self->qParameter, self->qLambda);
51  bool ok = true;
52  for(int j=0; j<ydim; j++) {
53  qREAL ytmp = y[i*ydim+j];
54  if(isnanq(ytmp) || isinfq(ytmp)) { ok = false; break; }
55  }
56 
57  if(!ok && (self->IntegrandMP!=NULL)) {
58  mpREAL mpx[xdim], mpy[ydim];
59  for(int j=0; j<xdim; j++) mpx[j] = x[i*xdim+j];
60  self->IntegrandMP(xdim, mpx, ydim, mpy, self->mpParameter, self->mpLambda);
61  for(int j=0; j<ydim; j++) y[i*ydim+j] = mpy[j].toFloat128();
62  }
63  } else {
64  dREAL dx[xdim], dy[ydim];
65  for(int j=0; j<xdim; j++) dx[j] = (dREAL)x[i*xdim+j];
66  self->IntegrandD(xdim, dx, ydim, dy, self->dParameter, self->dLambda);
67  for(int j=0; j<ydim; j++) y[i*ydim+j] = dy[j];
68  bool ok = true;
69  for(int j=0; j<ydim; j++) {
70  qREAL ytmp = y[i*ydim+j];
71  if(isnanq(ytmp) || isinfq(ytmp)) { ok = false; break; }
72  }
73 
74  if(!ok) self->IntegrandQ(xdim, x+i*xdim, ydim, y+i*ydim, self->qParameter, self->qLambda);
75 
76  ok = true;
77  for(int j=0; j<ydim; j++) {
78  qREAL ytmp = y[i*ydim+j];
79  if(isnanq(ytmp) || isinfq(ytmp)) { ok = false; break; }
80  }
81  if(!ok && (self->IntegrandMP!=NULL)) {
82  mpREAL mpx[xdim], mpy[ydim];
83  for(int j=0; j<xdim; j++) mpx[j] = x[i*xdim+j];
84  self->IntegrandMP(xdim, mpx, ydim, mpy, self->mpParameter, self->mpLambda);
85  for(int j=0; j<ydim; j++) y[i*ydim+j] = mpy[j].toFloat128();
86  }
87  }
88 
89  // Final Check NaN/Inf
90  bool ok = true;
91  for(int j=0; j<ydim; j++) {
92  qREAL ytmp = y[i*ydim+j];
93  if(isnanq(ytmp) || isinfq(ytmp)) { ok = false; break; }
94  }
95  if(!ok && (self->IntegrandMP!=NULL)) {
96  mpfr_free_cache();
97  mpfr::mpreal::set_default_prec(mpfr::digits2bits(self->MPDigits*10));
98  mpREAL mpx[xdim], mpy[ydim];
99  for(int j=0; j<xdim; j++) mpx[j] = x[i*xdim+j];
100  self->IntegrandMP(xdim, mpx, ydim, mpy, self->mpParameter, self->mpLambda);
101  for(int j=0; j<ydim; j++) y[i*ydim+j] = mpy[j].toFloat128();
102  mpfr::mpreal::set_default_prec(mpfr::digits2bits(self->MPDigits));
103  }
104 
105  // final check
106  for(int j=0; j<ydim; j++) {
107  qREAL ytmp = y[i*ydim+j];
108  if(isnanq(ytmp) || isinfq(ytmp)) {
109  #pragma omp atomic
110  self->nNAN++;
111  if(self->nNAN > self->NANMax) { NaNQ = true; break; }
112  else y[i*ydim+j] = 0;
113  }
114  }
115  if(self->ReIm == 1) y[i*ydim+1] = 0;
116  else if(self->ReIm == 2) y[i*ydim+0] = 0;
117  mpfr_free_cache();
118  }
119 
120  return NaNQ ? 1 : 0;
121  }
122 
123  void HCubature::DefaultPrintHooker(qREAL *result, qREAL *epsabs, size_t *nrun, void *fdata) {
124  auto self = (HCubature*)fdata;
125  if(*nrun == self->MaxPTS + 1979) return;
126  if(self->RunTime>0) {
127  auto cur_timer = time(NULL);
128  auto used_time = difftime(cur_timer,self->StartTimer);
129  if(used_time>self->RunTime) {
130  self->NEval = *nrun;
131  *nrun = self->MaxPTS + 1979;
132  if(Verbose>10) cout << WarnColor << " Exit with Run out of Time: " << used_time << RESET << endl;
133  return;
134  }
135  }
136 
137  if((*nrun-self->NEval) >= self->RunPTS) {
138  if(Verbose>10) {
139  char r0[64], r1[64], e0[32], e1[32];
140  quadmath_snprintf(r0, sizeof r0, "%.10QG", result[0]);
141  quadmath_snprintf(r1, sizeof r1, "%.10QG", result[1]);
142  quadmath_snprintf(e0, sizeof e0, "%.5QG", epsabs[0]);
143  quadmath_snprintf(e1, sizeof e1, "%.5QG", epsabs[1]);
144  cout << " N: " << (*nrun) << ", ";
145  if(self->ReIm==3 || self->ReIm==1) cout << "[" << r0 << ", " << e0 << "]";
146  if(self->ReIm==3 || self->ReIm==2) cout << "+I*[" << r1 << ", " << e1 << "]";
147  cout << endl;
148  }
149  self->NEval = *nrun;
150  }
151 
152  if((isnanq(result[0]) || isnanq(result[1]) || isnanq(epsabs[0]) || isnanq(epsabs[1])) || (isinfq(result[0]) || isinfq(result[1]) || isinfq(epsabs[0]) || isinfq(epsabs[1]))) {
153  self->NEval = *nrun;
154  *nrun = self->MaxPTS + 1979;
155  if(self->LastState>0) self->LastState = -1;
156  if(Verbose>10) cout << ErrColor << " Exit with NaN, LastN=" << self->lastNRUN << RESET << endl;
157  return;
158  }
159 
160  if(epsabs[0] > 1E30*self->EpsAbs || epsabs[1] > 1E30*self->EpsAbs) {
161  self->NEval = *nrun;
162  *nrun = self->MaxPTS + 1979;
163  if(self->LastState>0) self->LastState = -1;
164  if(Verbose>10) cout << WarnColor << " Exit with EpsAbs, LastN=" << self->lastNRUN << RESET << endl;
165  return;
166  }
167 
168  if((self->LastState == 0) || (epsabs[0]<=2*self->LastAbsErr[0] && epsabs[1]<=2*self->LastAbsErr[1])) {
169  self->LastResult[0] = result[0];
170  self->LastResult[1] = result[1];
171  self->LastAbsErr[0] = epsabs[0];
172  self->LastAbsErr[1] = epsabs[1];
173  self->LastState = 1;
174  self->lastNRUN = *nrun;
175  self->lastnNAN = self->nNAN;
176  }
177 
178  bool rExit = (epsabs[0] < self->EpsAbs+1E-50Q) || (epsabs[0] < fabsq(result[0])*self->EpsRel+1E-50Q);
179  bool iExit = (epsabs[1] < self->EpsAbs+1E-50Q) || (epsabs[1] < fabsq(result[1])*self->EpsRel+1E-50Q);
180  if(rExit && iExit && (*nrun)>self->MinPTS) {
181  self->NEval = *nrun;
182  *nrun = self->MaxPTS + 1979;
183  return;
184  }
185 
186  auto pid = getpid();
187  ostringstream fn;
188  fn << pid << ".int.done";
189  if(file_exists(fn.str().c_str())) {
190  self->NEval = *nrun;
191  *nrun = self->MaxPTS + 1979;
192  ostringstream cmd;
193  cmd << "rm " << fn.str();
194  auto rc = system(cmd.str().c_str());
195  if(Verbose>10) cout << " Exit: " << fn.str() << endl;
196  }
197  }
198 
199  ex HCubature::Integrate(size_t tn) {
200  if(mpfr_buildopt_tls_p()<=0) throw Error("Integrate: mpfr_buildopt_tls_p()<=0.");
201  mpfr_free_cache();
202  mpfr::mpreal::set_default_prec(mpfr::digits2bits(MPDigits));
203  mpPi = mpfr::const_pi();
204  mpEuler = mpfr::const_euler();
205  mpiEpsilon = complex<mpREAL>(0,mpfr::machine_epsilon()*100);
206 
207  unsigned int xdim = XDim;
208  unsigned int ydim = 2;
209  qREAL result[ydim], estabs[ydim];
210 
211  qREAL xmin[xdim], xmax[xdim];
212  for(int i=0; i<xdim; i++) {
213  xmin[i] = 0.0Q;
214  xmax[i] = 1.0Q;
215  }
216  LastState = 0;
217  NEval = 0;
218  nNAN = 0;
219 
220  size_t _MinPTS_, _RunPTS_;
221  if(tn==0) {
222  _RunPTS_ = RunPTS;
223  MaxPTS = RunPTS * RunMAX;
224  _MinPTS_ = MinPTS>0 ? MinPTS : _RunPTS_/10;
225  } else {
226  MaxPTS = tn;
227  if(MaxPTS<10000) MaxPTS = 10000;
228  _RunPTS_ = MaxPTS/5;
229  _MinPTS_ = MinPTS>0 ? MinPTS : _RunPTS_/10;
230  }
231  StartTimer = time(NULL);
232 
233  int nok = hcubature_v(ydim, Wrapper, this, xdim, xmin, xmax, _MinPTS_, _RunPTS_, MaxPTS, EpsAbs, EpsRel, result, estabs, tn==0 ? PrintHooker : NULL);
234 
235  if(nok) {
236  if( (cabsq(result[0]+result[1]*1.Qi) < FLT128_EPSILON) && (cabsq(estabs[0]+estabs[1]*1.Qi) < FLT128_EPSILON) ) {
237  cout << ErrColor << "HCubature Failed with 0 result returned!" << RESET << endl;
238  return NaN;
239  }
240  }
241 
242  if(LastState==-1 && use_last) {
243  result[0] = LastResult[0];
244  result[1] = LastResult[1];
245  estabs[0] = LastAbsErr[0];
246  estabs[1] = LastAbsErr[1];
247  NEval = lastNRUN;
248  nNAN = lastnNAN;
249  }
250 
251  ex FResult = 0;
252  if(isnanq(result[0]) || isnanq(result[1])) FResult += NaN;
253  else {
254  try{
255  FResult += VE(q2ex(result[0]), q2ex(estabs[0]));
256  FResult += VE(q2ex(result[1]), q2ex(estabs[1])) * I;
257  } catch(...) {
258  FResult += NaN;
259  }
260  }
261 
262  // Check nNAN / NEval
263  if(nNAN * 1000 > NEval && NEval>0) {
264  cout << ErrColor << "NAN=" << nNAN << " v.s. RUN=" << NEval << RESET << endl;
265  }
266 
267  return FResult;
268  }
269 
270  int HCubature::inDQMP(qREAL const *x) {
271  unsigned int xdim = XDim;
272 
273  if(xdim<=MPXDim) return 3;
274  qREAL xmin = 100;
275  for(int i=0; i<xdim; i++) {
276  if(x[i] < xmin) xmin = x[i];
277  }
278  if(xmin < MPXLimit) return 3;
279 
280  if(FT!=NULL) {
281  qREAL ft = 1E50;
282  static FT_Type last_ft = NULL;
283  static qREAL ft0;
284  if(last_ft!=FT) {
285  qREAL x0[xdim];
286  for(int i=0; i<xdim; i++) x0[i]=0.521Q;
287  qREAL ft0 = fabsq(FT(x0, qParameter));
288  if(ft0<1E-50) ft0 = 1;
289  last_ft = FT;
290  }
291  ft = fabsq(FT(x, qParameter));
292  ft = ft/ft0;
293  if(ft<MPFLimit) return 3;
294  else if(ft<QFLimit) return 2;
295  }
296 
297  if(xdim <= QXDim || xmin < QXLimit) return 2;
298 
299  return 1;
300  }
301 
302 }
#define RESET
Definition: BASIC.h:79
complex< mpREAL > mpCOMPLEX
Definition: HCubature.cpp:20
complex< dREAL > dCOMPLEX
Definition: HCubature.cpp:18
mpfr::mpreal mpREAL
Definition: HCubature.cpp:19
__complex128 qCOMPLEX
Definition: HCubature.cpp:16
long double dREAL
Definition: HCubature.cpp:17
__float128 qREAL
Definition: HCubature.cpp:15
const qCOMPLEX qiEpsilon
mpREAL mpEuler
mpCOMPLEX mpiEpsilon
mpREAL mpPi
bool file_exists(const char *fn)
Definition: Process.cpp:9
SecDec header file.
class used to wrap error message
Definition: BASIC.h:242
numerical integrator using HCubature
Definition: SD.h:166
namespace for Numerical integration with Sector Decomposition method
Definition: AsyMB.cpp:10
__float128 qREAL
Definition: SD.h:123
long double dREAL
Definition: SD.h:121
mpfr::mpreal mpREAL
Definition: SD.h:125
ex q2ex(__float128 num)
__float128 to ex
Definition: BASIC.cpp:867
const char * ErrColor
Definition: BASIC.cpp:246
const char * WarnColor
Definition: BASIC.cpp:247
int Verbose
Definition: Init.cpp:139
const Symbol NaN