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