HepLib
Loading...
Searching...
No Matches
Integrates.cpp
Go to the documentation of this file.
1
6#include "SD.h"
7#include <math.h>
8#include <cmath>
9
10namespace HepLib::SD {
11
18 ex SecDec::ContinuousWRA(ex expr_in, int nc) {
19 auto expr = expr_in;
20 static symbol xwra("xwra");
21 expr = xyz_pow_simplify(expr);
22 exset pows_set;
23 expr.find(sqrt(w), pows_set);
24 expr.find(pow(w1,w2), pows_set);
25 exmap pow_map;
26 for(auto item : pows_set) {
27 if(item.has(WRA(w))) {
28 if(item.match(pow(w1,w2)) && !item.op(1).info(info_flags::integer)) {
29 pow_map[item] = exp(log(item.op(0)) * item.op(1));
30 } else if(item.match(sqrt(w))) {
31 pow_map[item] = exp(log(item.op(0))/2);
32 }
33 }
34 }
35
36 expr = expr.subs(pow_map);
37 expr = exp_simplify(expr);;
38
39 int log_id = 0;
40 exset logs_set;
41 expr.find(log(w), logs_set);
42 exmap log_map;
43 lst id_logz_lst;
44 for(auto item : logs_set){
45 if(item.has(WRA(w))) {
46 log_id++;
47 id_logz_lst.append(iWF(log_id,item.op(0)));
48 log_map[item] = iWF(log_id);
49 }
50 }
51 if(log_id<1) return expr.subs(WRA(w)==w);
52
53 expr = expr.subs(log_map).subs(WRA(w)==w);
54
55 exmap log_map2;
56 for(auto id_logz : id_logz_lst) {
57 auto id = id_logz.op(0);
58 auto zz = id_logz.op(1);
59 exset wra_set;
60 zz.find(WRA(w), wra_set);
61 if(wra_set.size()<1) return expr;
62 if(wra_set.size()!=1) {
63 cout << zz << endl;
64 throw Error("ContinuousWRA: too many WRA.");
65 }
66 zz = zz.subs(WRA(w)==xwra);
67 ex wra = (*(wra_set.begin())).op(0);
68
69 ex ret = log(zz.subs(xwra==wra));
70 if(nc<1) {
71 log_map2[iWF(id)] = ret;
72 break;
73 }
74 int total=0;
75 int ReIm[nc][2];
76 for(int k=0; k<=nc; k++) {
77 ex zzk;
78 if(k==0) zzk = NN(zz.subs(xwra==wra/(25*nc)));
79 else zzk = NN(zz.subs(xwra==k*wra/nc));
80 if(!is_a<numeric>(zzk)) throw Error("ContinuousWRA: zzk is not numeric: "+ex2str(zzk));
81 auto nzzk = ex_to<numeric>(zzk);
82 auto curR = real(nzzk);
83 auto curI = imag(nzzk);
84 if(is_zero(curI) && k==0 && curR<0) ReIm[total][1] = -1;
85 else if(is_zero(curR) || is_zero(curR)) continue;
86 else ReIm[total][1] = curI>0 ? 1 : -1;
87 ReIm[total][0] = curR>0 ? 1 : -1;
88 total++;
89 }
90
91 int cutN = 0;
92 for(int k=0; k<total-1; k++) {
93 if(ReIm[k][0]*ReIm[k+1][0]<0 && ReIm[k][1]*ReIm[k+1][1]<0) throw Error("ContinuousWRA: 1<->3 or 2<->4 happened.");
94 if(ReIm[k][0]<0 && ReIm[k+1][0]<0 && ReIm[k][1]*ReIm[k+1][1]<0) {
95 if(ReIm[k][1]>0) cutN++;
96 else cutN--;
97 }
98 }
99 if(cutN!=0) ret += I * cutN * 2 * Pi;
100 log_map2[iWF(id)] = ret;
101 }
102 expr = expr.subs(log_map2);
103 if(expr.has(iWF(w))) throw Error("ContinuousWRA: iWF(w) still exists in final result.");
104 return expr;
105 }
106
113 void SecDec::Integrates(const string & key, const string & pkey, int kid) {
114 if(MPDigits>0) {
115 auto pbit = mpfr::digits2bits(MPDigits);
116 if(mpfr::mpreal::get_default_prec()!=pbit) mpfr::mpreal::set_default_prec(pbit);
117 }
118 if(IsZero) return;
119
120 if(Verbose>0) cout << Color_HighLight << " Integrates @ " << now() << RESET << endl;
121
122 lst lstRE;
123 auto pid = getpid();
124 ostringstream fsofn, sofn, cmd;
125 if(key == "") {
126 sofn << pid << ".so";
127 fsofn << pid << "F.so";
128 } else {
129 set_precision(NNDigits);
130 sofn << key << ".so";
131 ostringstream garfn;
132 garfn << key << ".ci.gar";
133 archive ar;
134 ifstream in(garfn.str());
135 in >> ar;
136 in.close();
137 auto c = ar.unarchive_ex("c");
138 if(c!=19790923) throw Error("Integrates: *.ci.gar error!");
139 auto gl = ar.unarchive_ex("soLimit");
140 eps_lst = ex_to<lst>(ar.unarchive_ex("eps_lst"));
141 soLimit = ex2int(gl);
142 auto res = ar.unarchive_ex("res");
143 ciResult.clear();
144 for(auto item : ex_to<lst>(res)) ciResult.push_back(ex_to<lst>(item));
145 garfn.clear();
146 garfn.str("");
147 garfn << key;
148 if(pkey != "") garfn << "-" << pkey;
149 garfn << ".las.gar";
150 if(file_exists(garfn.str().c_str())) {
151 archive la_ar;
152 ifstream la_in(garfn.str());
153 la_in >> la_ar;
154 la_in.close();
155 auto la_c = la_ar.unarchive_ex("c");
156 auto la_res = la_ar.unarchive_ex("res");
157 if(la_c!=19790923) throw Error("Integrates: *.ci.gar error!");
158 for(auto item : ex_to<lst>(la_res)) {
159 LambdaMap[item.op(0)] = item.op(1);
160 }
161 }
162
163 if(kid>0) {
164 garfn.clear();
165 garfn.str("");
166 garfn << key;
167 if(pkey != "") garfn << "-" << pkey;
168 garfn << ".res.gar";
169 if(!file_exists(garfn.str().c_str())) {
170 throw Error("Integrates: File Not Found: " + garfn.str());
171 }
172
173 archive res_ar;
174 ifstream res_in(garfn.str());
175 res_in >> res_ar;
176 res_in.close();
177 auto res_c = res_ar.unarchive_ex("c");
178 auto relst = res_ar.unarchive_ex("relst");
179 if(res_c!=19790923) throw Error("*.res.gar error with kid!");
180 lstRE = ex_to<lst>(relst);
181 }
182 reset_precision();
183 }
184
185 void* main_module = dlopen(sofn.str().c_str(), RTLD_NOW);
186 if(main_module == nullptr) {
187 main_module = dlopen(("./"+sofn.str()).c_str(), RTLD_NOW);
188 if(main_module == nullptr) {
189 cout << "dlerror(): " << dlerror() << endl;
190 throw Error("Integrates: could not open main module!");
191 }
192 }
193
194 if(!Debug && key == "") {
195 if(file_exists(fsofn.str().c_str())) remove(fsofn.str().c_str());
196 remove(sofn.str().c_str());
197 }
198 vector<void*> ex_modules;
199 for(int n=1; true; n++) {
200 ostringstream ex_sofn;
201 if(key == "") {
202 ex_sofn << pid << "X" << n << ".so";
203 } else {
204 ex_sofn << key << "X" << n << ".so";
205 }
206 if(file_exists(ex_sofn.str().c_str())) {
207 void* module = dlopen(ex_sofn.str().c_str(), RTLD_NOW);
208 if(module == nullptr) {
209 module = dlopen(("./"+ex_sofn.str()).c_str(), RTLD_NOW);
210 if(module == nullptr) {
211 cout << "dlerror(): " << dlerror() << endl;
212 throw Error("Integrates: could not open ex-module!");
213 }
214 }
215 ex_modules.push_back(module);
216 if(!Debug && key == "") remove(ex_sofn.str().c_str());
217 } else break;
218 }
219
220 int npara = 0;
221 lst plRepl;
222 for(auto kv : Parameter) {
223 plRepl.append(PL(kv.first)==kv.second);
224 if(kv.first>npara) npara = kv.first;
225 }
226 plRepl.sort();
227 plRepl.unique();
228
229 int total = ciResult.size(), current = 0;
230 qREAL stot = sqrtq(total*1.Q);
231
232 ResultError = 0;
233 //----------------------------------------------------------------
234 for(auto &item : ciResult) {
235 current++;
236 if(kid>0 && current != kid) continue;
237 if(Verbose>0) {
238 cout << "\r \r";
239 cout << PRE << "\\--Integrating [" <<current<<"/"<<total<< "] " << flush;
240 }
241
242 unsigned int xsize = 0;
243 ex co;
244 exvector xs, fxs;
245 xsize = ex_to<numeric>(item.op(1)).to_int();
246 co = item.op(2).subs(plRepl).subs(iEpsilon==iEpsilonN);
247 if(co.has(WRA(w))) co = ContinuousWRA(co);
248
249 if(xsize<1) {
250 // { expr, xs.size(), kvf.op(0), -1}
251 ex exint = item.op(0).subs(plRepl).subs(iEpsilon==iEpsilonN);
252 if(exint.has(WRA(w))) exint = ContinuousWRA(exint);
253 ex res = VE(NN(exint),0);
254 if(Verbose>5) {
255 cout << "XDim=" << xsize << endl;
256 cout << Color_HighLight << " IRes = "<< HepLib::SD::VEResult(VESimplify(res)) << RESET << endl;
257 }
258 if(kid>0) {
259 lstRE[kid-1] = res*co;
260 break;
261 }
262 lstRE.append(res*co);
263 continue;
264 }
265
266 //{ idx, xs.size(), kvf.op(0), ft_n }
267 ex intid = item.op(0);
268 ex ftid = item.op(3);
269
270 if(co.is_zero()) {
271 lstRE.append(0);
272 continue;
273 }
274 co = collect_ex(co, eps);
275 if(co.is_zero()) {
276 lstRE.append(0);
277 continue;
278 }
279 if(co.has(PL(w))) throw Error("Integrates: PL found @ " + ex2str(co));
280 qREAL cmax = -1;
281 int reim = 0;
282 if(ReIm==3) reim = 3;
283
284 for(int si=co.ldegree(eps); si<=co.degree(eps); si++) {
285 auto tmp = co.coeff(eps, si);
286 if(tmp.has(eps)) throw Error("Integrates: eps found @ " + ex2str(tmp));
287 tmp = collect_ex(tmp, ep);
288 for(int i=tmp.ldegree(ep); i<=tmp.degree(ep); i++) {
289 auto ccRes = NN(tmp.coeff(ep, i)).expand();
290 lst css;
291 css.append(ccRes);
292 if(is_a<add>(ccRes)) {
293 for(auto item : ccRes) css.append(item);
294 }
295
296 exmap eps_map;
297 ex epn = ex(1)/111;
298 for(auto epi : eps_lst) {
299 eps_map[epi.op(0)] = epn;
300 epn = epn / 13;
301 }
302
303 for(int ci=0; ci<css.nops(); ci++) {
304 auto nt = NN(css.op(ci).subs(epz==1).subs(log(vs)==1).subs(vs==1).subs(nReplacement).subs(lst{CV(w1,w2)==w2}).subs(eps_map),30);
305 if(nt.has(ep)) throw Error("Integrates: ep found @ nt = "+ex2str(nt));
306
307 lst nt_lst;
308 if(!is_a<numeric>(nt)) {
309 auto cv_lst = collect_lst(nt,[](const ex &e)->bool{return Symbol::has(e);});
310 for(auto nti : cv_lst) {
311 auto nnt = nti.op(0);
312 if(!is_a<numeric>(nnt)) throw Error("Integrates: Not a number with nt = "+ex2str(nnt));
313 nt_lst.append(nnt);
314 }
315 } else nt_lst.append(nt);
316
317 for(auto nnt : nt_lst) {
318 if(ReIm!=3 && reim!=3) {
319 if(ex_to<numeric>(nnt).imag()==0) {
320 if(reim==2) reim = 3;
321 else reim = 1;
322 } else if(ex_to<numeric>(nnt).real()==0) {
323 if(reim==1) reim = 3;
324 else reim = 2;
325 } else {
326 reim = 3;
327 }
328 }
329 nnt = NN(abs(nnt)); // no PL here, then nReplacement
330
331 qREAL qnt = ex2q(nnt);
332 if(qnt > cmax) cmax = qnt;
333 }
334 }
335 }
336 }
337 if(cmax<=0) {
338 MapFunction logF([](const ex & e, MapFunction & self)->ex{
339 if(e.match(log(w)) && e.op(0)<1) {
340 return -log(1/e.op(0));
341 } else if(e.match(w1*log(w2))) {
342 ex cv = e.subs(w1*log(w2)==lst{w1,w2});
343 if(ex_to<numeric>(cv.op(0)).is_integer()) {
344 return self(log(pow(cv.op(1), cv.op(0))));
345 } else return e.map(self);
346 } else return e.map(self);
347 });
348 while(true) {
349 auto co2 = subs(exfactor(co), lst{exp(w1*log(w2)+w3)==pow(w2,w1)*exp(w3),exp(w1*log(w2))==pow(w2,w1)});
350 co2 = co2.subs(lst{log(w1)+log(w2)==log(w1*w2),log(w1)-log(w2)==log(w1/w2)});
351 co2 = logF(co2);
352 if(is_zero(co2-co)) break;
353 co = co2;
354 }
355 if(normal(co).is_zero()) {
356 lstRE.append(0);
357 continue;
358 }
359 throw Error("Integrates: cmax<=0 with co = "+ex2str(co));
360 }
361 if(reim!=3 && ReIm!=3) {
362 if(ReIm==2) {
363 if(reim==1) reim=2;
364 else if(reim==2) reim=1;
365 }
366 }
367
368 char d1[20], d2[20];
369 sprintf(d1, "%.6G", (double)(EpsAbs/cmax/stot));
370 sprintf(d2, "%.6G", (double)cmax);
371 if(Verbose>5) cout << "XDim=" << xsize << ", EpsAbs=" << d1 << "/" << d2 << endl;
372
373 auto las = LambdaMap[ftid];
374 bool hasF = (ftid>0);
375 if(hasF && las.is_zero()) throw Error("Integrates: lambda with the key(ft_n=" + ex2str(ftid) + ") is NOT found!");
376
377 if(hasF && !is_a<lst>(las)) {
378 if(!is_zero(las-ex(1979))) { // the convention for xPositive or explict real mode
379 throw Error("Integrates: something is wrong with the F-term @ ft_n = "+ex2str(ftid) + ", las=" + ex2str(las));
380 } else {
381 hasF = false;
382 }
383 }
384
385 IntegratorBase::SDD_Type fpD = nullptr;
386 IntegratorBase::SDQ_Type fpQ = nullptr;
387 IntegratorBase::SDMP_Type fpMP = nullptr;
388 IntegratorBase::FT_Type ftp = nullptr;
389 int idx = ex_to<numeric>(intid).to_int();
390 auto module = main_module;
391 if(idx>=soLimit) module = ex_modules[idx/soLimit-1];
392 ostringstream fname;
393 if(hasF) fname << "C";
394 fname << "SDD_" << idx;
395 fpD = (IntegratorBase::SDD_Type)dlsym(module, fname.str().c_str());
396
397 if(fpD==NULL) {
398 cout << "dlerror(): " << dlerror() << endl;
399 throw Error("Integrates: fpD==NULL");
400 }
401
402 fname.clear();
403 fname.str("");
404 if(hasF) fname << "C";
405 fname << "SDQ_" << idx;
406 fpQ = (IntegratorBase::SDQ_Type)dlsym(module, fname.str().c_str());
407 if(fpQ==NULL) {
408 cout << "dlerror(): " << dlerror() << endl;
409 throw Error("Integrates: fpQ==NULL");
410 }
411
412 fname.clear();
413 fname.str("");
414 if(hasF) fname << "C";
415 fname << "SDMP_" << idx;
416 fpMP = (IntegratorBase::SDMP_Type)dlsym(module, fname.str().c_str());
417 if(fpMP==NULL) {
418 cout << "dlerror(): " << dlerror() << endl;
419 throw Error("Integrates: fpMP==NULL");
420 }
421
422 if(is_a<lst>(las)) {
423 fname.clear();
424 fname.str("");
425 fname << "FT_" << idx;
426 ftp = (IntegratorBase::FT_Type)dlsym(module, fname.str().c_str());
427 if(ftp==NULL) {
428 cout << "dlerror(): " << dlerror() << endl;
429 throw Error("Integrates: ftp==NULL.");
430 }
431 }
432
433 dREAL dlas[las.nops()], dpl[npara];
434 qREAL qlas[las.nops()], qpl[npara];
435 mpREAL mplas[las.nops()], mppl[npara];
436 for(auto kv : Parameter) {
437 qpl[kv.first] = ex2q(kv.second);
438 dpl[kv.first] = (dREAL)qpl[kv.first];
439 mppl[kv.first] = qpl[kv.first];
440 }
441
442 if(Integrator==NULL) Integrator = new HCubature();
443 Integrator->ReIm = reim;
444 if(MPDigits>0) Integrator->MPDigits = MPDigits;
445 Integrator->IntegrandD = fpD;
446 Integrator->IntegrandQ = fpQ;
447 Integrator->IntegrandMP = fpMP;
448 Integrator->FT = ftp;
449 Integrator->dParameter = dpl;
450 Integrator->dLambda = dlas;
451 Integrator->qParameter = qpl;
452 Integrator->qLambda = qlas;
453 Integrator->mpParameter = mppl;
454 Integrator->mpLambda = mplas;
455 Integrator->XDim = xsize;
456
457 if(hasF) {
458 Integrator->EpsAbs = EpsAbs/cmax/stot/2;
459 Integrator->EpsRel = 0;
460
461 qREAL lamax = ex2q(las.op(las.nops()-1));
462 if(lamax > IntLaMax) lamax = IntLaMax;
463 int ctryR = 0, ctry = 0, ctryL = 0;
464 int smin = -1;
465 ex min_err, min_res;
466 size_t min_eval;
467 qREAL log_lamax = log10q(lamax);
468 qREAL log_lamin = log_lamax-1.Q;
469
470 ostringstream las_fn;
471 las_fn << key;
472 if(pkey != "") las_fn << "-" << pkey;
473 las_fn << "-" << current << ".las";
474 if(use_las && file_exists(las_fn.str().c_str())) {
475 std::ifstream las_ifs;
476 las_ifs.open(las_fn.str(), ios::in);
477 if (!las_ifs) throw Error("Integrates: failed to open *.las file!");
478 for(int i=0; i<las.nops()-1; i++) {
479 dREAL la_tmp;
480 las_ifs >> la_tmp;
481 dlas[i] = la_tmp;
482 qlas[i] = la_tmp;
483 mplas[i] = la_tmp;
484 }
485 las_ifs.close();
486 auto res = Integrator->Integrate();
487 auto res_tmp = res.subs(VE(w1, w2)==w2);
488 auto err = real_part(res_tmp);
489 if(err < imag_part(res_tmp)) err = imag_part(res_tmp);
490 min_err = err;
491 min_res = res;
492 } else {
493 // ---------------------------------------
494 auto prec = cout.precision();
495 cout.precision(5);
496 while(true) {
497 smin = -1;
498 min_err = 0;
499 ex lastResErr = 0;
500 bool err_break = false;
501 for(int s=0; s<=LambdaSplit; s++) {
502 if(Verbose>10 && s==0) {
503 if(ctryR>0 || ctry>0 || ctryL>0)
504 cout << " ------------------------------" << endl;
505 }
506 auto log_cla = (log_lamin + s * (log_lamax-log_lamin) / LambdaSplit);
507 auto cla = powq(10.Q, log_cla);
508 if(cla < 1E-10) throw Error("NIntegrate: too small lambda.");
509 for(int i=0; i<las.nops()-1; i++) {
510 qlas[i] = ex2q(las.op(i)) * cla;
511 mplas[i] = qlas[i];
512 dlas[i] = (dREAL)qlas[i];
513 }
514
515 auto res = Integrator->Integrate(CTryI);
516 if(Verbose>10) {
517 cout << "\r \r";
518 if(res.has(NaN)) cout << " λ=" << (double)cla << "/" << Integrator->NEval << ": " << NaN << endl;
519 else cout << " λ=" << (double)cla << "/" << Integrator->NEval << ": " << VEResult2(VESimplify(res)) << endl;
520 }
521
522 if(res.has(NaN) && s==0) continue;
523 else if(res.has(NaN)) break;
524 ex res_abs = NN(abs(res.subs(VE(w1,w2)==w1)));
525 if(lastResErr.is_zero()) lastResErr = res;
526 auto diff = VESimplify(lastResErr - res);
527 diff = diff.subs(VE(0,0)==0);
528 exset ves;
529 diff.find(VE(w0, w1), ves);
530 for(auto ve : ves) {
531 auto ve0 = abs(ve.op(0));
532 if(ve0>ve.op(1)) {
533 if(numeric("1.E10")*ve0<res_abs) continue; // avoid fluctuation aroud 0
534 if(numeric("1.E10")*ve0<q2ex(EpsAbs)) continue; // avoid fluctuation aroud 0
535 err_break = true;
536 break;
537 }
538 }
539 if(err_break) {
540 if(Verbose>10) cout << Color_HighLight << " Error Break ..." << RESET << endl;
541 break;
542 }
543 lastResErr = res;
544
545 auto res_tmp = res.subs(VE(w1, w2)==w2);
546 auto err = real_part(res_tmp);
547 if(err < imag_part(res_tmp)) err = imag_part(res_tmp);
548 if(smin<0 || err < min_err) {
549 min_err = err;
550 min_res = res;
551 min_eval = Integrator->NEval;
552 smin = s;
553 }
554 if(s>0 && min_err < q2ex(EpsAbs/cmax/stot)) {
555 // s>0 make sure at least 2 λs compatiable
556 if(Verbose>5) {
557 cout << Color_HighLight << " λ=" << (double)cla << "/" << min_eval << ": " << HepLib::SD::VEResult(VESimplify(min_res)) << RESET << endl;
558 }
559
560 smin = -2;
561 if(kid>0) lstRE[kid-1] = co * min_res;
562 else lstRE.append(co * min_res);
563 break;
564 }
565 if(err > 100 * min_err) break; // s>0 make sure at least 2 λs compatiable
566 }
567 if(smin == -2) break;
568 if(smin == -1) throw Error("Integrates: smin = -1, too small lambda (<1E-10)!");
569
570 if(smin <= 0) {
571 if((!err_break) && (ctryL >= CTryL || ctryR>0)) break;
572 log_lamax = log_lamin;
573 log_lamin -= 1.Q;
574 if(!err_break) ctryL++;
575 } else if(smin >= LambdaSplit) {
576 if(ctryR >= CTryR || ctryL>0) break;
577 log_lamin = log_lamax;
578 log_lamax += log10q(CTryRRatio);
579 ctryR++;
580 } else {
581 if(ctry >= CTryM) break;
582 auto la1 = log_lamin + (smin-1) * (log_lamax-log_lamin) / LambdaSplit;
583 auto la2 = log_lamin + (smin+1) * (log_lamax-log_lamin) / LambdaSplit;
584 log_lamin = la1;
585 log_lamax = la2;
586 ctry++;
587 }
588 }
589 cout.precision(prec);
590
591 if(smin == -2) {
592 if(kid>0) break;
593 continue;
594 }
595
596 auto log_cla = (log_lamin + smin * (log_lamax-log_lamin) / LambdaSplit);
597 auto cla = powq(10.Q, log_cla);
598 if(Verbose>5) cout << Color_HighLight << " Final λ = " << (double)cla << " / " << las.op(las.nops()-1) << RESET << endl;
599 for(int i=0; i<las.nops()-1; i++) {
600 qlas[i] = ex2q(las.op(i)) * cla;
601 dlas[i] = (dREAL)qlas[i];
602 mplas[i] = qlas[i];
603 }
604 // ---------------------------------------
605 }
606
607 // ---------------------------------------
608 // try HookeJeeves
609 // ---------------------------------------
610 if( use_ErrMin && (min_err > q2ex(1E5 * EpsAbs/cmax/stot)) ) {
611 ErrMin::lastResErr = min_res;
612 auto miner = new HookeJeeves();
613 ErrMin::miner = miner;
614 ErrMin::Integrator = Integrator;
615 dREAL oo[las.nops()-1], ip[las.nops()-1];
616 for(int i=0; i<las.nops()-1; i++) ip[i] = oo[i] = (dREAL)qlas[i];
617 ErrMin::lambda = oo;
618 ErrMin::err_max = 1E100;
619 auto oerrmin = ErrMin::err_min;
620 ErrMin::err_min = (dREAL)(oerrmin < 0 ? -oerrmin * ex2q(min_err) : oerrmin/cmax);
621 ErrMin::RunRND = 0;
622 miner->Minimize(las.nops()-1, ErrMin::IntError, ip);
623 delete miner;
624 ErrMin::err_min = oerrmin;
625 for(int i=0; i<las.nops()-1; i++) {
626 qlas[i] = ErrMin::lambda[i];
627 dlas[i] = (dREAL)qlas[i];
628 mplas[i] = qlas[i];
629 }
630
631 Integrator->dLambda = dlas;
632 Integrator->qLambda = qlas; // Integrator->Lambda changed in ErrMin
633 Integrator->mpLambda = mplas;
634 if(Verbose>5) {
635 cout << Color_HighLight << " Final λs: " << RESET;
636 for(int i=0; i<xsize; i++) {
637 char buffer[128];
638 quadmath_snprintf(buffer, sizeof buffer, "%.6QG", qlas[i]);
639 cout << buffer << " ";
640 }
641 cout << endl << " ------------------------------" << endl;
642 }
643 }
644 // ---------------------------------------
645
646 if(save_las) {
647 std::ofstream las_ofs;
648 las_ofs.open(las_fn.str(), ios::out);
649 if (las_ofs) {
650 for(int i=0; i<las.nops()-1; i++) {
651 dREAL la_tmp = (dREAL)qlas[i];
652 las_ofs << la_tmp << " ";
653 }
654 las_ofs << endl;
655 las_ofs.close();
656 }
657 }
658 }
659
660 Integrator->EpsAbs = EpsAbs/cmax/stot;
661 Integrator->EpsRel = 0;
662 auto res = Integrator->Integrate();
663 if(Verbose>5) {
664 cout << Color_HighLight << " IRes = "<< HepLib::SD::VEResult(VESimplify(res)) << RESET << endl;
665 }
666 if(res.has(NaN)) {
667 ResultError = NaN;
668 if(kid>0) lstRE[kid-1] = NaN;
669 else lstRE.append(NaN);
670 break;
671 } else {
672 if(kid>0) {
673 lstRE[kid-1] = co * res;
674 break;
675 } else lstRE.append(co * res);
676 }
677 }
678 if(Verbose>0 && Verbose<=5) cout << endl;
679 //----------------------------------------------------------------
680
681 if(use_dlclose) {
682 for(auto module : ex_modules) dlclose(module);
683 dlclose(main_module);
684 }
685
686 if(!ResultError.is_equal(NaN)) {
687 ResultError = 0;
688 for(auto item : lstRE) ResultError += item;
689 ResultError = VESimplify(ResultError);
690
691 if(key != "") {
692 ostringstream garfn;
693 garfn << key;
694 if(pkey != "") garfn << "-" << pkey;
695 garfn << ".res.gar";
696 archive ar;
697 ar.archive_ex(ResultError, "res");
698 ar.archive_ex(lstRE, "relst");
699 ar.archive_ex(19790923, "c");
700 ofstream out(garfn.str());
701 out << ar;
702 out.close();
703 }
704 }
705 }
706
713 void SecDec::ReIntegrates(const string & key, const string & pkey, qREAL err) {
714 if(IsZero) return;
715 if(Verbose>0) cout << Color_HighLight << " Integrates @ " << now() << RESET << endl;
716
717 lst lstRE;
718 auto pid = getpid();
719 ostringstream fsofn, sofn, cmd;
720 if(true) {
721 set_precision(NNDigits);
722 sofn << key << ".so";
723 ostringstream garfn;
724 garfn << key << ".ci.gar";
725 archive ar;
726 ifstream in(garfn.str());
727 in >> ar;
728 in.close();
729 auto c = ar.unarchive_ex("c");
730 if(c!=19790923) throw Error("Integrates: *.ci.gar error!");
731 auto gl = ar.unarchive_ex("soLimit");
732 eps_lst = ex_to<lst>(ar.unarchive_ex("eps_lst"));
733 soLimit = ex2int(gl);
734 auto res = ar.unarchive_ex("res");
735 ciResult.clear();
736 for(auto item : ex_to<lst>(res)) ciResult.push_back(ex_to<lst>(item));
737 garfn.clear();
738 garfn.str("");
739 garfn << key;
740 if(pkey != "") garfn << "-" << pkey;
741 garfn << ".las.gar";
742 if(file_exists(garfn.str().c_str())) {
743 archive la_ar;
744 ifstream la_in(garfn.str());
745 la_in >> la_ar;
746 la_in.close();
747 auto la_c = la_ar.unarchive_ex("c");
748 auto la_res = la_ar.unarchive_ex("res");
749 if(la_c!=19790923) throw Error("Integrates: *.ci.gar error!");
750 for(auto item : ex_to<lst>(la_res)) {
751 LambdaMap[item.op(0)] = item.op(1);
752 }
753 }
754
755 if(true) {
756 garfn.clear();
757 garfn.str("");
758 garfn << key;
759 if(pkey != "") garfn << "-" << pkey;
760 garfn << ".res.gar";
761 if(!file_exists(garfn.str().c_str())) {
762 throw Error("Integrates: File Not Found: " + garfn.str());
763 }
764
765 archive res_ar;
766 ifstream res_in(garfn.str());
767 res_in >> res_ar;
768 res_in.close();
769 auto res_c = res_ar.unarchive_ex("c");
770 auto relst = res_ar.unarchive_ex("relst");
771 if(res_c!=19790923) throw Error("*.res.gar error with ReIntegrates!");
772 lstRE = ex_to<lst>(relst);
773 }
774 reset_precision();
775
776 if(lstRE.nops()!=ciResult.size()) throw Error("Error: the size is not same between ciResult and lstRE!");
777 }
778
779 void* main_module = dlopen(sofn.str().c_str(), RTLD_NOW);
780 if(main_module == nullptr) {
781 main_module = dlopen(("./"+sofn.str()).c_str(), RTLD_NOW);
782 if(main_module == nullptr) {
783 cout << "dlerror(): " << dlerror() << endl;
784 throw Error("Integrates: could not open main module!");
785 }
786 }
787
788 vector<void*> ex_modules;
789 for(int n=1; true; n++) {
790 ostringstream ex_sofn;
791 ex_sofn << key << "X" << n << ".so";
792 if(file_exists(ex_sofn.str().c_str())) {
793 void* module = dlopen(ex_sofn.str().c_str(), RTLD_NOW);
794 if(module == nullptr) {
795 module = dlopen(("./"+ex_sofn.str()).c_str(), RTLD_NOW);
796 if(module == nullptr) {
797 cout << "dlerror(): " << dlerror() << endl;
798 throw Error("Integrates: could not open ex-module!");
799 }
800 }
801 ex_modules.push_back(module);
802 if(!Debug && key == "") remove(ex_sofn.str().c_str());
803 } else break;
804 }
805
806 int npara = 0;
807 lst plRepl;
808 for(auto kv : Parameter) {
809 plRepl.append(PL(kv.first)==kv.second);
810 if(kv.first>npara) npara = kv.first;
811 }
812 plRepl.sort();
813 plRepl.unique();
814
815 int total = ciResult.size(), current = 0;
816 qREAL stot = sqrtq(total*1.Q);
817
818 ResultError = 0;
819 //----------------------------------------------------------------
820 for(auto &item : ciResult) {
821 current++;
822 auto cmerr = ex2q(VEMaxErr(lstRE.op(current-1)));
823 if(cmerr < err) continue;
824 if(Verbose>10) {
825 char es[64];
826 quadmath_snprintf(es, sizeof es, "%.10QG", cmerr);
827 cout << PRE << "\\--Current Err: " << es << endl;
828 }
829 if(Verbose>0) {
830 cout << "\r \r";
831 cout << PRE << "\\--ReIntegrating [" <<current<<"/"<<total<< "] " << flush;
832 }
833
834 unsigned int xsize = 0;
835 ex co;
836 exvector xs, fxs;
837 xsize = ex_to<numeric>(item.op(1)).to_int();
838 co = item.op(2).subs(plRepl).subs(iEpsilon==iEpsilonN);
839 if(co.has(WRA(w))) co = ContinuousWRA(co);
840
841 if(xsize<1) {
842 // { expr, xs.size(), kvf.op(0), -1}
843 ex exint = item.op(0).subs(plRepl).subs(iEpsilon==iEpsilonN);
844 if(exint.has(WRA(w))) exint = ContinuousWRA(exint);
845 ex res = VE(NN(exint),0);
846 if(Verbose>5) {
847 cout << "XDim=" << xsize << endl;
848 cout << Color_HighLight << " IRes = "<< HepLib::SD::VEResult(VESimplify(res)) << RESET << endl;
849 }
850 lstRE[current-1] = res*co;
851 continue;
852 }
853
854 //{ idx, xs.size(), kvf.op(0), ft_n }
855 ex intid = item.op(0);
856 ex ftid = item.op(3);
857
858 if(co.is_zero()) {
859 lstRE[current-1] = 0;
860 continue;
861 }
862 co = collect_ex(co, eps);
863 if(co.is_zero()) {
864 lstRE[current-1] = 0;
865 continue;
866 }
867 if(co.has(PL(w))) throw Error("Integrates: PL found @ " + ex2str(co));
868 qREAL cmax = -1;
869 int reim = 0;
870 if(ReIm==3) reim = 3;
871
872 for(int si=co.ldegree(eps); si<=co.degree(eps); si++) {
873 auto tmp = co.coeff(eps, si);
874 if(tmp.has(eps)) throw Error("Integrates: eps found @ " + ex2str(tmp));
875 tmp = collect_ex(tmp, ep);
876 for(int i=tmp.ldegree(ep); i<=tmp.degree(ep); i++) {
877 auto ccRes = NN(tmp.coeff(ep, i)).expand();
878 lst css;
879 css.append(ccRes);
880 if(is_a<add>(ccRes)) {
881 for(auto item : ccRes) css.append(item);
882 }
883
884 exmap eps_map;
885 ex epn = ex(1)/111;
886 for(auto epi : eps_lst) {
887 eps_map[epi.op(0)] = epn;
888 epn = epn / 13;
889 }
890
891 for(int ci=0; ci<css.nops(); ci++) {
892 auto nt = NN(css.op(ci).subs(epz==1).subs(log(vs)==1).subs(vs==1).subs(nReplacement).subs(lst{CV(w1,w2)==w2}).subs(eps_map),30);
893 if(nt.has(ep)) throw Error("Integrates: ep found @ nt = "+ex2str(nt));
894
895 lst nt_lst;
896 if(!is_a<numeric>(nt)) {
897 auto cv_lst = collect_lst(nt,[](const ex &e)->bool{return Symbol::has(e);});
898 for(auto nti : cv_lst) {
899 auto nnt = nti.op(0);
900 if(!is_a<numeric>(nnt)) throw Error("Integrates: Not a number with nt = "+ex2str(nnt));
901 nt_lst.append(nnt);
902 }
903 } else nt_lst.append(nt);
904
905 for(auto nnt : nt_lst) {
906 if(ReIm!=3 && reim!=3) {
907 if(ex_to<numeric>(nnt).imag()==0) {
908 if(reim==2) reim = 3;
909 else reim = 1;
910 } else if(ex_to<numeric>(nnt).real()==0) {
911 if(reim==1) reim = 3;
912 else reim = 2;
913 } else {
914 reim = 3;
915 }
916 }
917 nnt = NN(abs(nnt)); // no PL here, then nReplacement
918
919 qREAL qnt = ex2q(nnt);
920 if(qnt > cmax) cmax = qnt;
921 }
922 }
923 }
924 }
925 if(cmax<=0) {
926 MapFunction logF([](const ex & e, MapFunction & self)->ex{
927 if(e.match(log(w)) && e.op(0)<1) {
928 return -log(1/e.op(0));
929 } else if(e.match(w1*log(w2))) {
930 ex cv = e.subs(w1*log(w2)==lst{w1,w2});
931 if(ex_to<numeric>(cv.op(0)).is_integer()) {
932 return self(log(pow(cv.op(1), cv.op(0))));
933 } else return e.map(self);
934 } else return e.map(self);
935 });
936 while(true) {
937 auto co2 = subs(exfactor(co),lst{exp(w1*log(w2)+w3)==pow(w2,w1)*exp(w3),exp(w1*log(w2))==pow(w2,w1)});
938 co2 = co2.subs(lst{log(w1)+log(w2)==log(w1*w2),log(w1)-log(w2)==log(w1/w2)});
939 co2 = logF(co2);
940 if(is_zero(co2-co)) break;
941 co = co2;
942 }
943 if(normal(co).is_zero()) {
944 lstRE[current-1] = 0;
945 continue;
946 }
947 throw Error("Integrates: cmax<=0 with co = "+ex2str(co));
948 }
949 if(reim!=3 && ReIm!=3) {
950 if(ReIm==2) {
951 if(reim==1) reim=2;
952 else if(reim==2) reim=1;
953 }
954 }
955
956 char d1[20], d2[20];
957 sprintf(d1, "%.6G", (double)(EpsAbs/cmax/stot));
958 sprintf(d2, "%.6G", (double)cmax);
959 if(Verbose>5) cout << "XDim=" << xsize << ", EpsAbs=" << d1 << "/" << d2 << endl;
960 auto las = LambdaMap[ftid];
961 bool hasF = (ftid>0);
962 if(hasF && las.is_zero()) throw Error("Integrates: lambda with the key(ft_n=" + ex2str(ftid) + ") is NOT found!");
963 if(hasF && !is_a<lst>(las)) {
964 if(!is_zero(las-ex(1979))) { // the convention for xPositive or explict real mode
965 throw Error("Integrates: something is wrong with the F-term @ ft_n = "+ex2str(ftid) + ", las=" + ex2str(las));
966 } else {
967 hasF = false;
968 }
969 }
970
971 IntegratorBase::SDD_Type fpD = nullptr;
972 IntegratorBase::SDQ_Type fpQ = nullptr;
973 IntegratorBase::SDMP_Type fpMP = nullptr;
974 IntegratorBase::FT_Type ftp = nullptr;
975 int idx = ex_to<numeric>(intid).to_int();
976 auto module = main_module;
977 if(idx>=soLimit) module = ex_modules[idx/soLimit-1];
978 ostringstream fname;
979 if(hasF) fname << "C";
980 fname << "SDD_" << idx;
981 fpD = (IntegratorBase::SDD_Type)dlsym(module, fname.str().c_str());
982 if(fpD==NULL) {
983 cout << "dlerror(): " << dlerror() << endl;
984 throw Error("Integrates: fpD==NULL");
985 }
986
987 fname.clear();
988 fname.str("");
989 if(hasF) fname << "C";
990 fname << "SDQ_" << idx;
991 fpQ = (IntegratorBase::SDQ_Type)dlsym(module, fname.str().c_str());
992 if(fpQ==NULL) {
993 cout << "dlerror(): " << dlerror() << endl;
994 throw Error("Integrates: fpQ==NULL");
995 }
996
997 fname.clear();
998 fname.str("");
999 if(hasF) fname << "C";
1000 fname << "SDMP_" << idx;
1001 fpMP = (IntegratorBase::SDMP_Type)dlsym(module, fname.str().c_str());
1002
1003 if(is_a<lst>(las)) {
1004 fname.clear();
1005 fname.str("");
1006 fname << "FT_" << idx;
1007 ftp = (IntegratorBase::FT_Type)dlsym(module, fname.str().c_str());
1008 if(ftp==NULL) {
1009 cout << "dlerror(): " << dlerror() << endl;
1010 throw Error("Integrates: ftp==NULL.");
1011 }
1012 }
1013
1014 dREAL dlas[las.nops()], dpl[npara];
1015 qREAL qlas[las.nops()], qpl[npara];
1016 mpREAL mplas[las.nops()], mppl[npara];
1017 for(auto kv : Parameter) {
1018 qpl[kv.first] = ex2q(kv.second);
1019 dpl[kv.first] = (dREAL)qpl[kv.first];
1020 mppl[kv.first] = qpl[kv.first];
1021 }
1022
1023 if(Integrator==NULL) Integrator = new HCubature();
1024 Integrator->ReIm = reim;
1025 if(MPDigits>0) Integrator->MPDigits = MPDigits;
1026 Integrator->IntegrandD = fpD;
1027 Integrator->IntegrandQ = fpQ;
1028 Integrator->IntegrandMP = fpMP;
1029 Integrator->FT = ftp;
1030 Integrator->dParameter = dpl;
1031 Integrator->dLambda = dlas;
1032 Integrator->qParameter = qpl;
1033 Integrator->qLambda = qlas;
1034 Integrator->mpParameter = mppl;
1035 Integrator->mpLambda = mplas;
1036 Integrator->XDim = xsize;
1037
1038 if(hasF) {
1039 qREAL lamax = ex2q(las.op(las.nops()-1));
1040 if(lamax > IntLaMax) lamax = IntLaMax;
1041
1042 int ctryR = 0, ctry = 0, ctryL = 0;
1043 int smin = -1;
1044 ex min_err, min_res;
1045 size_t min_eval;
1046 qREAL log_lamax = log10q(lamax);
1047 qREAL log_lamin = log_lamax-1.Q;
1048
1049 ostringstream las_fn;
1050 las_fn << key;
1051 if(pkey != "") las_fn << "-" << pkey;
1052 las_fn << "-" << current << ".las";
1053 if(use_las && file_exists(las_fn.str().c_str())) {
1054 std::ifstream las_ifs;
1055 las_ifs.open(las_fn.str(), ios::in);
1056 if (!las_ifs) throw Error("Integrates: failed to open *.las file!");
1057 for(int i=0; i<las.nops()-1; i++) {
1058 dREAL la_tmp;
1059 las_ifs >> la_tmp;
1060 qlas[i] = la_tmp;
1061 dlas[i] = (dREAL)qlas[i];
1062 mplas[i] = qlas[i];
1063 }
1064 las_ifs.close();
1065 auto res = Integrator->Integrate();
1066 auto res_tmp = res.subs(VE(w1, w2)==w2);
1067 auto err = real_part(res_tmp);
1068 if(err < imag_part(res_tmp)) err = imag_part(res_tmp);
1069 min_err = err;
1070 min_res = res;
1071 } else {
1072 // ---------------------------------------
1073 while(true) {
1074 smin = -1;
1075 min_err = 0;
1076 ex lastResErr = 0;
1077 bool err_break = false;
1078 for(int s=0; s<=LambdaSplit; s++) {
1079 if(Verbose>10 && s==0) {
1080 if(ctryR>0 || ctry>0 || ctryL>0)
1081 cout << " ------------------------------" << endl;
1082 }
1083 auto log_cla = (log_lamin + s * (log_lamax-log_lamin) / LambdaSplit);
1084 auto cla = powq(10.Q, log_cla);
1085 if(cla < 1E-10) throw Error("NIntegrate: too small lambda.");
1086 for(int i=0; i<las.nops()-1; i++) {
1087 qlas[i] = ex2q(las.op(i)) * cla;
1088 dlas[i] = (dREAL)qlas[i];
1089 mplas[i] = qlas[i];
1090 }
1091
1092 auto res = Integrator->Integrate(CTryI);
1093 if(Verbose>10) {
1094 cout << "\r \r";
1095 if(res.has(NaN)) cout << " λ=" << (double)cla << "/" << Integrator->NEval << ": " << NaN << endl;
1096 else cout << " λ=" << (double)cla << "/" << Integrator->NEval << ": " << VEResult2(VESimplify(res)) << endl;
1097 }
1098
1099 if(res.has(NaN) && s==0) continue;
1100 else if(res.has(NaN)) break;
1101 ex res_abs = NN(abs(res.subs(VE(w1,w2)==w1)));
1102 if(lastResErr.is_zero()) lastResErr = res;
1103 auto diff = VESimplify(lastResErr - res);
1104 diff = diff.subs(VE(0,0)==0);
1105 exset ves;
1106 diff.find(VE(w0, w1), ves);
1107 for(auto ve : ves) {
1108 auto ve0 = abs(ve.op(0));
1109 if(ve0>ve.op(1)) {
1110 if(numeric("1.E10")*ve0<res_abs) continue; // avoid fluctuation aroud 0
1111 if(numeric("1.E10")*ve0<q2ex(EpsAbs)) continue; // avoid fluctuation aroud 0
1112 err_break = true;
1113 break;
1114 }
1115 }
1116 if(err_break) {
1117 if(Verbose>10) cout << Color_HighLight << " Error Break ..." << RESET << endl;
1118 break;
1119 }
1120 lastResErr = res;
1121
1122 auto res_tmp = res.subs(VE(w1, w2)==w2);
1123 auto err = real_part(res_tmp);
1124 if(err < imag_part(res_tmp)) err = imag_part(res_tmp);
1125 if(smin<0 || err < min_err) {
1126 min_err = err;
1127 min_res = res;
1128 min_eval = Integrator->NEval;
1129 smin = s;
1130 }
1131 if(s>0 && min_err < q2ex(EpsAbs/cmax/stot)) {
1132 // s>0 make sure at least 2 λs compatiable
1133 if(Verbose>5) {
1134 cout << Color_HighLight << " λ=" << (double)cla << "/" << min_eval << ": " << HepLib::SD::VEResult(VESimplify(min_res)) << RESET << endl;
1135 }
1136
1137 smin = -2;
1138 lstRE[current-1] = co * min_res;
1139 break;
1140 }
1141 if(err > 100 * min_err) break; // s>0 make sure at least 2 λs compatiable
1142 }
1143 if(smin == -2) break;
1144 if(smin == -1) throw Error("Integrates: smin = -1, too small lambda (<1E-10)!");
1145
1146 if(smin <= 0) {
1147 if((!err_break) && (ctryL >= CTryL || ctryR>0)) break;
1148 log_lamax = log_lamin;
1149 log_lamin -= 1.Q;
1150 if(!err_break) ctryL++;
1151 } else if(smin >= LambdaSplit) {
1152 if(ctryR >= CTryR || ctryL>0) break;
1153 log_lamin = log_lamax;
1154 log_lamax += log10q(CTryRRatio);
1155 ctryR++;
1156 } else {
1157 if(ctry >= CTryM) break;
1158 auto la1 = log_lamin + (smin-1) * (log_lamax-log_lamin) / LambdaSplit;
1159 auto la2 = log_lamin + (smin+1) * (log_lamax-log_lamin) / LambdaSplit;
1160 log_lamin = la1;
1161 log_lamax = la2;
1162 ctry++;
1163 }
1164 }
1165
1166 if(smin == -2) continue;
1167
1168 auto log_cla = (log_lamin + smin * (log_lamax-log_lamin) / LambdaSplit);
1169 auto cla = powq(10.Q, log_cla);
1170 if(Verbose>5) cout << Color_HighLight << " Final λ = " << (double)cla << " / " << las.op(las.nops()-1) << RESET << endl;
1171 for(int i=0; i<las.nops()-1; i++) {
1172 qlas[i] = ex2q(las.op(i)) * cla;
1173 dlas[i] = (dREAL)qlas[i];
1174 mplas[i] = qlas[i];
1175 }
1176 // ---------------------------------------
1177 }
1178
1179 // ---------------------------------------
1180 // try HookeJeeves
1181 // ---------------------------------------
1182 if( use_ErrMin && (min_err > q2ex(1E5 * EpsAbs/cmax/stot)) ) {
1183 ErrMin::lastResErr = min_res;
1184 auto miner = new HookeJeeves();
1185 ErrMin::miner = miner;
1186 ErrMin::Integrator = Integrator;
1187 dREAL oo[las.nops()-1], ip[las.nops()-1];
1188 for(int i=0; i<las.nops()-1; i++) ip[i] = oo[i] = (dREAL)qlas[i];
1189 ErrMin::lambda = oo;
1190 ErrMin::err_max = 1E100;
1191 auto oerrmin = ErrMin::err_min;
1192 ErrMin::err_min = (dREAL)(oerrmin < 0 ? -oerrmin * ex2q(min_err) : oerrmin/cmax);
1193 ErrMin::RunRND = 0;
1194 miner->Minimize(las.nops()-1, ErrMin::IntError, ip);
1195 delete miner;
1196 ErrMin::err_min = oerrmin;
1197 for(int i=0; i<las.nops()-1; i++) {
1198 qlas[i] = ErrMin::lambda[i];
1199 dlas[i] = (dREAL)qlas[i];
1200 mplas[i] = qlas[i];
1201 }
1202
1203 Integrator->dLambda = dlas;
1204 Integrator->qLambda = qlas; // Integrator->Lambda changed in ErrMin
1205 Integrator->mpLambda = mplas;
1206 if(Verbose>5) {
1207 cout << Color_HighLight << " Final λs: " << RESET;
1208 for(int i=0; i<xsize; i++) {
1209 char buffer[128];
1210 quadmath_snprintf(buffer, sizeof buffer, "%.6QG", qlas[i]);
1211 cout << buffer << " ";
1212 }
1213 cout << endl << " ------------------------------" << endl;
1214 }
1215 }
1216 // ---------------------------------------
1217
1218 if(save_las) {
1219 std::ofstream las_ofs;
1220 las_ofs.open(las_fn.str(), ios::out);
1221 if (las_ofs) {
1222 for(int i=0; i<las.nops()-1; i++) {
1223 dREAL la_tmp = (dREAL)qlas[i];
1224 las_ofs << la_tmp << " ";
1225 }
1226 las_ofs << endl;
1227 las_ofs.close();
1228 }
1229 }
1230 }
1231
1232 auto res = Integrator->Integrate();
1233 if(Verbose>5) {
1234 cout << Color_HighLight << " IRes = "<< HepLib::SD::VEResult(VESimplify(res)) << RESET << endl;
1235 }
1236 if(res.has(NaN)) {
1237 ResultError = NaN;
1238 lstRE[current-1] = NaN;
1239 } else {
1240 lstRE[current-1] = co * res;
1241 }
1242 }
1243 if(Verbose>0 && Verbose<=5) cout << endl;
1244
1245 //----------------------------------------------------------------
1246
1247 if(use_dlclose) {
1248 for(auto module : ex_modules) dlclose(module);
1249 dlclose(main_module);
1250 }
1251
1252 if(!ResultError.is_equal(NaN)) {
1253 ResultError = 0;
1254 for(auto item : lstRE) ResultError += item;
1255 ResultError = VESimplify(ResultError);
1256 }
1257
1258 if(key != "") {
1259 ostringstream garfn;
1260 garfn << key;
1261 if(pkey != "") garfn << "-" << pkey;
1262 garfn << ".res.gar";
1263 archive ar;
1264 ar.archive_ex(ResultError, "res");
1265 ar.archive_ex(lstRE, "relst");
1266 ar.archive_ex(19790923, "c");
1267 ofstream out(garfn.str());
1268 out << ar;
1269 out.close();
1270 }
1271 }
1272
1273}
RESET
#define RESET
Definition BASIC.h:79
SD.h
SecDec header file.
HepLib::Error
class used to wrap error message
Definition BASIC.h:242
HepLib::MapFunction
class to wrap map_function of GiNaC
Definition BASIC.h:675
HepLib::MapFunction::subs
static ex subs(const ex &expr, const ex &pat, std::function< ex(const ex &)> f)
Definition Functions.cpp:42
HepLib::SD::ErrMin::lastResErr
static ex lastResErr
Definition SD.h:379
HepLib::SD::ErrMin::Integrator
static IntegratorBase * Integrator
Definition SD.h:370
HepLib::SD::ErrMin::RunRND
static size_t RunRND
Definition SD.h:375
HepLib::SD::ErrMin::err_max
static dREAL err_max
Definition SD.h:372
HepLib::SD::ErrMin::miner
static MinimizeBase * miner
Definition SD.h:376
HepLib::SD::ErrMin::IntError
static dREAL IntError(int nvars, dREAL *las, dREAL *n1, dREAL *n2)
Definition ErrMin.cpp:22
HepLib::SD::ErrMin::lambda
static dREAL * lambda
Definition SD.h:377
HepLib::SD::ErrMin::err_min
static dREAL err_min
Definition SD.h:373
HepLib::SD::HCubature
numerical integrator using HCubature
Definition SD.h:167
HepLib::SD::HookeJeeves
class to minimize a function using HookeJeeves
Definition SD.h:279
HepLib::SD::IntegratorBase::IntegrandMP
SDMP_Type IntegrandMP
Definition SD.h:145
HepLib::SD::IntegratorBase::dParameter
const dREAL * dParameter
Definition SD.h:147
HepLib::SD::IntegratorBase::FT_Type
qREAL(* FT_Type)(const qREAL xx[], const qREAL pl[])
Definition SD.h:139
HepLib::SD::IntegratorBase::SDD_Type
int(* SDD_Type)(const unsigned int xn, const dREAL x[], const unsigned int yn, dREAL y[], const dREAL pl[], const dREAL las[])
Definition SD.h:136
HepLib::SD::IntegratorBase::dLambda
const dREAL * dLambda
Definition SD.h:146
HepLib::SD::IntegratorBase::mpLambda
const mpREAL * mpLambda
Definition SD.h:150
HepLib::SD::IntegratorBase::SDQ_Type
int(* SDQ_Type)(const unsigned int xn, const qREAL x[], const unsigned int yn, qREAL y[], const qREAL pl[], const qREAL las[])
Definition SD.h:137
HepLib::SD::IntegratorBase::qParameter
const qREAL * qParameter
Definition SD.h:149
HepLib::SD::IntegratorBase::Integrate
virtual ex Integrate(size_t n=0)=0
HepLib::SD::IntegratorBase::SDMP_Type
int(* SDMP_Type)(const unsigned int xn, const mpREAL x[], const unsigned int yn, mpREAL y[], const mpREAL pl[], const mpREAL las[])
Definition SD.h:138
HepLib::SD::IntegratorBase::mpParameter
const mpREAL * mpParameter
Definition SD.h:151
HepLib::SD::IntegratorBase::qLambda
const qREAL * qLambda
Definition SD.h:148
HepLib::SD::SecDec::LambdaSplit
size_t LambdaSplit
Definition SD.h:451
HepLib::SD::SecDec::nReplacement
exmap nReplacement
Definition SD.h:424
HepLib::SD::SecDec::use_ErrMin
bool use_ErrMin
Definition SD.h:434
HepLib::SD::SecDec::CTryRRatio
dREAL CTryRRatio
Definition SD.h:457
HepLib::SD::SecDec::IntLaMax
qREAL IntLaMax
Definition SD.h:452
HepLib::SD::SecDec::Parameter
map< int, numeric > Parameter
Definition SD.h:443
HepLib::SD::SecDec::EpsAbs
qREAL EpsAbs
Definition SD.h:460
HepLib::SD::SecDec::save_las
bool save_las
Definition SD.h:436
HepLib::SD::SecDec::Integrator
IntegratorBase * Integrator
Definition SD.h:429
HepLib::SD::SecDec::use_dlclose
static bool use_dlclose
Definition SD.h:415
HepLib::SD::SecDec::Integrates
void Integrates(const string &key="", const string &pkey="", int kid=0)
Contours, note that here we need to provide the specific Parameter.
Definition Integrates.cpp:113
HepLib::SD::SecDec::ContinuousWRA
static ex ContinuousWRA(ex expr_in, int nc=15)
ContinuousWRA, note that here we need to provide the specific Parameter.
Definition Integrates.cpp:18
HepLib::SD::SecDec::ReIntegrates
void ReIntegrates(const string &key, const string &pkey, qREAL err)
Contours, note that here we need to provide the specific Parameter.
Definition Integrates.cpp:713
HepLib::SD::SecDec::CTryI
size_t CTryI
Definition SD.h:456
HepLib::Symbol::has
static bool has(const ex &e)
HepLib::SD
namespace for Numerical integration with Sector Decomposition method
Definition AsyMB.cpp:10
HepLib::SD::xyz_pow_simplify
ex xyz_pow_simplify(const ex expr_in)
Definition Helpers.cpp:349
HepLib::SD::exp_simplify
ex exp_simplify(const ex expr_in)
Definition Helpers.cpp:315
HepLib::SD::qREAL
__float128 qREAL
Definition SD.h:124
HepLib::SD::VEResult
ex VEResult(ex expr)
Definition Helpers.cpp:123
HepLib::SD::dREAL
long double dREAL
Definition SD.h:122
HepLib::SD::VESimplify
ex VESimplify(ex expr)
Definition Helpers.cpp:92
HepLib::SD::VEMaxErr
ex VEMaxErr(ex expr)
Definition Helpers.cpp:131
HepLib::SD::mpREAL
mpfr::mpreal mpREAL
Definition SD.h:126
HepLib::SD::VEResult2
ex VEResult2(ex expr)
Definition Helpers.cpp:127
HepLib::Color_HighLight
const char * Color_HighLight
Definition BASIC.cpp:248
HepLib::exfactor
ex exfactor(const ex &expr_in, int opt)
factorize a expression
Definition BASIC.cpp:1854
HepLib::ex2q
__float128 ex2q(ex num)
ex of numeric to __float128
Definition BASIC.cpp:879
HepLib::iEpsilon
const iSymbol iEpsilon
HepLib::reset_precision
void reset_precision()
Definition BASIC.cpp:2323
HepLib::q2ex
ex q2ex(__float128 num)
__float128 to ex
Definition BASIC.cpp:867
HepLib::NNDigits
int NNDigits
Definition Init.cpp:158
HepLib::NN
ex NN(ex expr, int digits)
the nuerical evaluation
Definition BASIC.cpp:1278
HepLib::w0
ex w0
Definition BASIC.h:499
HepLib::ep
const Symbol ep
HepLib::collect_ex
ex collect_ex(ex const &expr_in, std::function< bool(const ex &)> has_func, int opt)
the collect function like Mathematica
Definition BASIC.cpp:1202
HepLib::eps
const Symbol eps
HepLib::file_exists
bool file_exists(string fn)
Definition BASIC.h:289
HepLib::now
string now(bool use_date)
date/time string
Definition BASIC.cpp:525
HepLib::Debug
bool Debug
Definition Init.cpp:144
HepLib::w
ex w
Definition Init.cpp:93
HepLib::vs
const Symbol vs
HepLib::iEpsilonN
const ex iEpsilonN
Definition Init.cpp:141
HepLib::epz
const Symbol epz
HepLib::collect_lst
lst collect_lst(ex const &expr_in, std::function< bool(const ex &)> has_func, int opt)
the collect function like Mathematica, reture the lst { {c1,v1}, {c2,v2}, ... }
Definition BASIC.cpp:1222
HepLib::Verbose
int Verbose
Definition Init.cpp:142
HepLib::ex2str
string ex2str(const ex &expr)
convert ex to output string, the defalut printer format will be used
Definition BASIC.cpp:715
HepLib::w1
ex w1
Definition BASIC.h:499
HepLib::w3
ex w3
Definition BASIC.h:499
HepLib::set_precision
void set_precision(long prec, bool push)
Definition BASIC.cpp:2314
HepLib::w2
ex w2
Definition BASIC.h:499
HepLib::PRE
string PRE
Definition Init.cpp:143
HepLib::ex2int
int ex2int(ex num)
ex to integer
Definition BASIC.cpp:893
HepLib::NaN
const Symbol NaN
HepLib::subs
ex subs(const ex &e, init_list sl)
Definition BASIC.h:51
Generated by doxygen 1.9.8