HepLib
Loading...
Searching...
No Matches
InitEval.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
13 if(SecDec!=NULL) delete SecDec;
14 if(Integrator!=NULL) delete Integrator;
15 if(Minimizer!=NULL) delete Minimizer;
16 SecDec = NULL;
17 Integrator = NULL;
18 Minimizer = NULL;
19 }
20
21 void Replacement2(exmap &repl) {
22 auto tmp = repl;
23 for(auto &kv : repl) {
24 kv.second = Symbol::set_all(kv.second.subs(tmp));
25 }
26 }
27
29
30 if(fp.Propagator.nops()<1) {
31 FunExp.clear();
32 if(fp.LoopMomenta.nops()>0 || fp.tLoopMomenta.nops()>0) {
33 IsZero = true;
34 return;
35 } else {
36 FunExp.push_back(lst{lst{Symbol::set_all(fp.Prefactor)}, lst{1}});
37 }
38 return;
39 }
40
41 if(fp.Propagator.nops() != fp.Exponent.nops()) {
42 throw Error("Initialize: the length of Propagator and Exponent are NOT equal.");
43 }
44
45 if(Symbol::set_all(fp.Prefactor).is_zero()) {
46 IsZero = true;
47 return;
48 }
49
50 IsZero = false;
51
52 bool wFound = false;
53 for(auto kv : fp.lReplacement) {
54 if(has_w(kv.first)) {
55 wFound = true;
56 break;
57 }
58 }
59 if(!wFound) {
60 auto repl = fp.lReplacement;
61 for(auto kv : repl) fp.lReplacement[w*kv.first] = w*kv.second;
62 }
63 wFound = false;
64 for(auto kv : fp.tReplacement) {
65 if(has_w(kv.first)) {
66 wFound = true;
67 break;
68 }
69 }
70 if(!wFound) {
71 auto repl = fp.tReplacement;
72 for(auto kv : repl) fp.tReplacement[w*kv.first] = w*kv.second;
73 }
74 wFound = false;
75 for(auto kv : fp.nReplacement) {
76 if(has_w(kv.first)) {
77 wFound = true;
78 break;
79 }
80 }
81 if(!wFound) {
82 auto repl = fp.nReplacement;
83 for(auto kv : repl) fp.nReplacement[w*kv.first] = w*kv.second;
84 }
85
86
87 for(auto kv: fp.lReplacement) {
88 if((lst{kv.first, kv.second}).has(iEpsilon)) {
89 throw Error("Initialize: (lst{kv.first, kv.second}).has(iEpsilon) @1");
90 }
91 }
92 for(auto kv: fp.tReplacement) {
93 if((lst{kv.first, kv.second}).has(iEpsilon)) {
94 throw Error("Initialize: (lst{kv.first, kv.second}).has(iEpsilon) @2");
95 }
96 }
97 for(auto kv: fp.nReplacement) {
98 if((lst{kv.first, kv.second}).has(iEpsilon)) {
99 throw Error("Initialize: (lst{kv.first, kv.second}).has(iEpsilon) @3");
100 }
101 }
102
103 auto ps = Symbol::set_all(fp.Propagator);
104 auto ns = Symbol::set_all(fp.Exponent);
105
106 auto ls = fp.LoopMomenta;
107 auto tls = fp.tLoopMomenta;
108
109 for(auto item : ls) {
110 if(!is_a<symbol>(item)) throw Error("SecDec::Initialize failed, NOT a symbol: "+ex2str(item));
111 }
112 for(auto item : tls) {
113 if(!is_a<symbol>(item)) throw Error("SecDec::Initialize failed, NOT a symbol: "+ex2str(item));
114 }
115
119
120 auto lsubs = fp.lReplacement;
121 auto tsubs = fp.tReplacement;
122 auto nsubs = fp.nReplacement;
124
125 if(Verbose > 0) cout << Color_HighLight << " Initialize @ " << now() << RESET << endl;
126
127 ex asgn = 1;
128 ex a = 0;
129 ex ad = (4-2*ep)*ls.nops();
130 if(fp.isQuasi) ad += (3-2*ep)*tls.nops();
131 else ad += (2-2*ep)*tls.nops();
132 int xn = ps.nops();
133 ex rem = 0;
134 exmap xtNeg;
135
136 ex pre = Symbol::set_all(fp.Prefactor); // come from below
137 for(int i=0; i<ps.nops(); i++) {
138 bool ltQ = false; {
139 auto tps = Symbol::set_all(ps.op(i).expand().subs(lsubs).subs(tsubs));
140 for(auto lsi : ls) {
141 if(tps.has(lsi)) {
142 ltQ = true;
143 break;
144 }
145 }
146
147 if(!ltQ) {
148 for(auto lsi : tls) {
149 if(tps.has(lsi)) {
150 ltQ = true;
151 break;
152 }
153 }
154 }
155 }
156
157 if(ltQ) a += ns.op(i);
158 ex sgn = 0;
159
160 if(!ltQ) {
161 pre = pre * pow(Symbol::set_all(ps.op(i).expand().subs(lsubs).subs(tsubs)), ex(0)-ns.op(i));
162 ns.let_op(i) = 0;
163 ps.let_op(i) = 1;
164 continue;
165 } else if(ns.op(i).info(info_flags::negint) || is_zero(ns.op(i))) {
166 xtNeg[x(i)]=0;
167 if(is_zero(ns.op(i))) continue;
168 }
169
170 auto p = ps.op(i).expand().subs(lsubs).subs(tsubs).subs(nsubs);
171 p = Symbol::set_all(p);
172 p = p.subs(lsubs).subs(tsubs).subs(nsubs);
173 p = Symbol::set_all(p);
174
175 // check loop^2
176 for(auto m : ls) {
177 if(!is_a<numeric>(p.coeff(m,2))) {
178 cout << ErrColor << "not numeric: " << p.coeff(m,2) << endl;
179 cout << "nsubs = " << nsubs << RESET << endl;
180 exit(1);
181 }
182 numeric nm = ex_to<numeric>(p.coeff(m,2));
183 if(nm.is_zero()) continue;
184 sgn = nm>0 ? -1 : 1;
185 break;
186 }
187 // check iEpsilon
188 if(sgn.is_zero()) {
189 if(!is_a<numeric>(p.coeff(iEpsilon))) {
190 throw Error("Initialize: (!is_a<numeric>(p.coeff(iEpsilon)))");
191 }
192 numeric nm = ex_to<numeric>(p.coeff(iEpsilon));
193 if(!nm.is_zero()) sgn = nm>0 ? -1 : 1;
194 }
195 // check tloop^2
196 if(sgn.is_zero()) {
197 for(auto m : tls) {
198 if(!is_a<numeric>(p.coeff(m,2))) {
199 throw Error("Initialize: NOT numeric: " + ex2str(p.coeff(m,2)));
200 }
201 numeric nm = ex_to<numeric>(p.coeff(m,2));
202 if(nm.is_zero()) continue;
203 sgn = nm>0 ? -1 : 1;
204 break;
205 }
206 }
207 // others
208 if(sgn.is_zero()) {
209 sgn = 1;
210 if(is_a<numeric>(p) && ex_to<numeric>(p)>0) sgn = -1;
211 if(Verbose>0 && (!is_a<numeric>(ns.op(i)) || ns.op(i)>0)) {
212 cout << WarnColor << " - Warning: Can NOT determine the iEpsilon sign." << RESET << endl;
213 cout << WarnColor << " - " << p << " from " << ps.op(i) << RESET << endl;
214 }
215 }
216
217 p = (ps.op(i)*sgn).subs(iEpsilon==0);
218 if(sgn==-1) asgn *= exp(I * Pi * ns.op(i));
219 rem += x(i) * p;
220 }
221
222 rem = Symbol::set_all(rem.expand());
223 lst uList1, uList2;
224 ex u=1, cu=1;
225
226 // Loop
227 if(ls.nops()>0) {
228 u=1;
229 for(int i=0; i<ls.nops(); i++) {
230 auto t2 = rem.coeff(ls.op(i),2);
231 auto t1 = rem.coeff(ls.op(i),1);
232 auto t0 = rem.coeff(ls.op(i),0);
233 u *= (-t2);
234 if(t2==0) {
235 IsZero = true;
236 return;
237 }
238 rem = expand(t0 - pow(t1,2)/(4*t2));
239 }
240 rem = normal(Symbol::set_all(rem.subs(lsubs).subs(lsubs)));
241 u = normal(Symbol::set_all(u.subs(lsubs)));
242 for(auto m: tls) {
243 if(u.has(m)) {
244 cerr << ErrColor << "Initialize: u.has(m), " << u << ", " << m << RESET << endl;
245 exit(1);
246 }
247 }
248
249 cu *= u;
250 auto u_nd = numer_denom(u);
251 ex usgn = FactorOutX(u_nd.op(1)).subs(xtNeg).subs(x(w)==ex(1)/2).subs(nsubs);
252 if(is_zero(usgn)) usgn = FactorOutX(u_nd.op(1)).subs(xtNeg).subs(x(w)==ex(1)/3).subs(nsubs);
253 if(!is_a<numeric>(usgn) || is_zero(usgn)) throw Error("Initialize: usgn is zero or non-numeric.");
254 usgn = normal(usgn)>0 ? 1 : -1;
255
256 uList1.append(usgn*u_nd.op(0));
257 uList2.append(-(4-2*ep)/2);
258 if((usgn*u_nd.op(0)) != 1) {
259 uList1.append(usgn*u_nd.op(1));
260 uList2.append((4-2*ep)/2);
261 }
262 } else {
263 rem = normal(Symbol::set_all(rem.subs(lsubs).subs(lsubs)));
264 }
265
266 // t-Loop
267 if(tls.nops()>0) {
268 u=1;
269 for(int i=0; i<tls.nops(); i++) {
270 auto t2 = rem.coeff(tls.op(i),2);
271 auto t1 = rem.coeff(tls.op(i),1);
272 auto t0 = rem.coeff(tls.op(i),0);
273 u *= (-t2);
274 if(t2.is_zero()) {
275 IsZero = true;
276 return;
277 }
278 rem = expand(t0 - pow(t1,2)/(4*t2));
279 }
280 rem = normal(Symbol::set_all(rem.subs(tsubs)));
281 u = normal(Symbol::set_all(u.subs(lsubs)));
282 for(auto m: tls) {
283 if(u.has(m)) {
284 cerr << ErrColor << "Initialize: u.has(m), " << u << ", " << m << RESET << endl;
285 exit(1);
286 }
287 }
288
289 cu *= u;
290 auto u_nd = numer_denom(u);
291 ex usgn = FactorOutX(u_nd.op(1)).subs(xtNeg).subs(x(w)==ex(1)/2).subs(nsubs);
292 if(is_zero(usgn)) usgn = FactorOutX(u_nd.op(1)).subs(xtNeg).subs(x(w)==ex(1)/3).subs(nsubs);
293 if(!is_a<numeric>(usgn) || is_zero(usgn)) throw Error("Initialize: usgn is zero or non-numeric.");
294 usgn = normal(usgn)>0 ? 1 : -1;
295
296 uList1.append(usgn*u_nd.op(0));
297 if(fp.isQuasi) uList2.append(-(3-2*ep)/2);
298 else uList2.append(-(2-2*ep)/2);
299 if(usgn*u_nd.op(1) != 1) {
300 uList1.append(usgn*u_nd.op(1));
301 if(fp.isQuasi) uList2.append((3-2*ep)/2);
302 else uList2.append((2-2*ep)/2);
303 }
304 }
305
306 u = normal(cu);
307 auto u_nd = numer_denom(u);
308 rem = normal(rem * u);
309 auto rem_nd = numer_denom(rem);
310
311 ex usgn = FactorOutX(u_nd.op(1)).subs(xtNeg).subs(x(w)==ex(1)/2).subs(nsubs);
312 if(!is_a<numeric>(usgn) || is_zero(usgn)) throw Error("Initialize: usgn is zero or non-numeric.");
313 usgn = normal(usgn)>0 ? 1 : -1;
314 ex fsgn = FactorOutX(rem_nd.op(1)).subs(xtNeg).subs(x(w)==ex(1)/2).subs(nsubs);
315 if(!is_a<numeric>(fsgn) || is_zero(fsgn)) throw Error("Initialize: fsgn is zero or non-numeric.");
316 fsgn = normal(fsgn)>0 ? 1 : -1;
317
318 lst fList1, fList2;
319 fList1.append(usgn*u_nd.op(0));
320 fList2.append(a-ad/2);
321 fList1.append(fsgn*rem_nd.op(0));
322 fList2.append(-a+ad/2);
323 if(usgn*u_nd.op(1) != 1) {
324 fList1.append(usgn*u_nd.op(1));
325 fList2.append(-a+ad/2);
326 }
327 if(fsgn*rem_nd.op(1) != 1) {
328 fList1.append(fsgn*rem_nd.op(1));
329 fList2.append(a-ad/2);
330 }
331
332 for(int i=0; i<uList1.nops(); i++) {
333 fList1.append(uList1[i]);
334 fList2.append(uList2[i]);
335 }
336
337 exvector ret;
338 ret.push_back(lst{fList1, fList2});
339
340 // negative index
341 for(int i=0; i<xn; i++) {
342 if(ns.op(i).info(info_flags::negint)) {
343 for(int j=0; j<ex(0)-ns.op(i); j++) {
344 exvector nret;
345 for(auto fe : ret) {
346 auto plst = ex_to<lst>(fe.op(0));
347 auto nlst = ex_to<lst>(fe.op(1));
348
349 ex nxi=0;
350 for(int ij=0; ij<plst.nops(); ij++) {
351 auto ldeg = expand_ex(plst.op(ij),x(w)).ldegree(x(i));
352 if(ldeg>0) {
353 plst.let_op(ij) = collect_common_factors(plst.op(ij) / pow(x(i),ldeg));
354 nxi += ldeg * nlst.op(ij);
355 }
356 }
357 nxi = normal(nxi);
358 if(!is_zero(nxi)) {
359 plst.append(x(i));
360 nlst.append(nxi);
361 }
362
363 for(int ij=0; ij<nlst.nops(); ij++) { // note the "-" sign
364 auto dtmp = ex(0)-nlst.op(ij) * diff_ex(plst.op(ij),x(i));
365 if(is_zero(dtmp)) continue;
366 auto plst2 = plst;
367 auto nlst2 = nlst;
368 if(is_zero(nlst.op(ij)-1)) {
369 plst2.let_op(ij) = dtmp;
370 } else {
371 nlst2.let_op(ij) = nlst.op(ij)-1;
372 int nn = plst.nops();
373 if(!is_zero(nlst.op(nn-1)-1)) {
374 plst2.append(dtmp);
375 nlst2.append(1);
376 } else plst2.let_op(nn-1) = plst.op(nn-1) * dtmp;
377 }
378 nret.push_back(lst{plst2, nlst2});
379 }
380 }
381 ret = exvector(std::move(nret));
382 }
383
384 for(auto &fe : ret) {
385 ex xpn = 0;
386 int nps = fe.op(0).nops();
387 for(int k=0; k<nps; k++) {
388 auto tmp = fe.op(0).op(k);
389 auto ldeg = expand_ex(tmp,x(w)).ldegree(x(i));
390 if(ldeg>0) {
391 xpn += ldeg * fe.op(1).op(k);
392 tmp = collect_common_factors(tmp / pow(x(i),ldeg));
393 }
394 tmp = tmp.subs(x(i)==0);
395 let_op(fe,0,k,tmp);
396 }
397
398 xpn = normal(xpn);
399 if(!is_zero(xpn)) {
400 if(is_a<numeric>(xpn) && xpn<0) {
401 cout << "xpn=" << xpn << " :> " << i << endl;
402 cout << fe << endl;
403 throw Error("Initialize: xpn < 0 found.");
404 }
405 fe = 0;
406 }
407 }
408 }}
409
410 // simplification
411 // ex pre = fp.Prefactor; // moved to above
412 pre *= asgn * pow(I,ls.nops()+(fp.isQuasi ? tls.nops() : 0)) * pow(Pi, ad/2) * tgamma(a-ad/2);
413 for(int i=0; i<ns.nops(); i++) {
414 if(is_a<numeric>(ns.op(i)) && ns.op(i)<=0) continue;
415 pre /= tgamma(ns.op(i));
416 }
417 if(tls.nops()>0 && (!fp.isQuasi)) pre *= exp(I * Pi * tls.nops()*(2-2*ep)/2);
418
419 ex xpre = 1;
420 for(int i=0; i<ns.nops(); i++) {
421 if(is_a<numeric>(ns.op(i)) && ns.op(i)<=1) continue;
422 else {
423 for(auto &fe : ret) {
424 if(is_zero(fe)) continue;
425 let_op_append(fe, 0, x(i));
426 let_op_append(fe, 1, ns.op(i)-1);
427 }
428 }
429 }
430
431 for(auto &fe : ret) {
432 if(is_zero(fe)) continue;
433 let_op_append(fe, 0, pre);
434 let_op_append(fe, 1, 1);
435 if(xpre != 1) {
436 let_op_append(fe, 0, xpre);
437 let_op_append(fe, 1, 1);
438 }
439 auto nnn = fe.op(0).nops();
440 for(int j=0; j<nnn; j++) {
441 fe.let_op(0).let_op(j) = collect_common_factors(fe.op(0).op(j));
442 fe.let_op(1).let_op(j) = collect_common_factors(fe.op(1).op(j));
443 }
444 }
445
446 lst delta;
447 for(int i=0; i<ns.nops(); i++) {
448 if(is_a<numeric>(ns.op(i)) && ns.op(i)<=0) continue;
449 delta.append(x(i));
450 }
451 FunExp.clear();
452 for(auto &fe : ret) {
453 if(is_zero(fe)) continue;
454 let_op_append(fe, lst{delta});
455 FunExp.push_back(fe);
456 }
457 Normalizes();
458 exmap eps_map;
459 ex epn = ex(1)/111;
460 for(auto epi : eps_lst) {
461 eps_map[epi.op(0)] = epn;
462 epn = epn / 13;
463 }
464 // final check non-f term positive
465 for(auto fe : FunExp) {
466 auto fs = fe.op(0);
467 auto ns = fe.op(1);
468 for(int i=0; i<fs.nops(); i++) {
469 if(i==1 || ns.op(i).info(info_flags::integer)) continue;
470 auto nv = Symbol::set_all(fs.op(i)).subs(nReplacement).subs(lst{CV(w1,w2)==w2}).subs(eps_map);
471 if(!xPositive(nv)) {
472 cout << "fs = " << fs << endl << "ns = " << ns << endl;
473 throw Error("Initialize: non-positive u-term found.");
474 }
475 }
476 }
477
478 if(fp.isAsy) DoAsy();
479 XReOrders();
480 Normalizes();
481 }
482
484 IsZero = false;
487
488 for(int di=0; di<xint.Deltas.nops(); di++) {
489 auto delta = xint.Deltas.op(di);
490 if(!is_a<lst>(delta) || delta.nops()<1) {
491 cout << ErrColor << "Deltas is NOT valide: " << xint.Deltas << RESET << endl;
492 exit(1);
493 }
494 }
495
496 FunExp.clear();
497 FunExp.shrink_to_fit();
498 if(xint.Deltas.nops()>0) FunExp.push_back(lst{xint.Function, xint.Exponent, xint.Deltas});
499 else FunExp.push_back(lst{xint.Function, xint.Exponent});
500
501 Normalizes();
502 if(xint.isAsy) DoAsy();
503 XReOrders();
504 Normalizes();
505 }
506
507 void SecDec::Evaluate(FeynmanParameter fp, const string & key) {
508
509 if(Verbose>1) cout << endl << " Starting @ " << now() << endl;
510
511 Initialize(fp);
512 MB();
513 if(FunExp.size()<1) return;
514 Scalelesses();
515 ChengWu();
516 RemoveDeltas();
517 KillPowers();
518 SDPrepares();
519 EpsExpands();
520 CIPrepares(key);
521 auto pps = GiNaC_Parallel_Process;
523 Contours(key);
524 Integrates(key);
526 if(Verbose>1) cout << " Finished @ " << now() << endl << endl;
527 }
528
529 void SecDec::Evaluate(XIntegrand xint, const string & key) {
530
531 if(Verbose>1) cout << endl << " Starting @ " << now() << endl;
532
533 Initialize(xint);
534 MB();
535 if(FunExp.size()<1) return;
536 Scalelesses();
537 ChengWu();
538 RemoveDeltas();
539 KillPowers();
540 SDPrepares();
541 EpsExpands();
542 CIPrepares(key);
543 auto pps = GiNaC_Parallel_Process;
545 Contours(key);
546 Integrates(key);
548 if(Verbose>1) cout << " Finished @ " << now() << endl << endl;
549 }
550
551 void SecDec::Evaluate(const exvector & funexp, const string & key) {
552
553 if(Verbose>1) cout << endl << " Starting @ " << now() << endl;
554
555 FunExp = funexp;
556 MB();
557 if(FunExp.size()<1) return;
558 Scalelesses();
559 ChengWu();
560 RemoveDeltas();
561 KillPowers();
562 SDPrepares();
563 EpsExpands();
564 CIPrepares(key);
565 auto pps = GiNaC_Parallel_Process;
567 Contours(key);
568 Integrates(key);
570 if(Verbose>1) cout << " Finished @ " << now() << endl << endl;
571 }
572
573}
int * a
#define RESET
Definition BASIC.h:79
SecDec header file.
class used to wrap error message
Definition BASIC.h:242
class to manipulate with Cheng-Wu theorem
Definition SD.h:387
SecDec the main class to use Sector Decompostion method.
Definition SD.h:413
void KillPowers(int bits=1+2)
Kill Powers will call KillPowerD or KillPower.
exmap nReplacement
Definition SD.h:425
SecDecBase * SecDec
Definition SD.h:429
MinimizeBase * Minimizer
Definition SD.h:431
void Contours(const string &key="", const string &pkey="")
Contours, note that here we need to provide the specific Parameter.
Definition Contours.cpp:17
IntegratorBase * Integrator
Definition SD.h:430
void Initialize(FeynmanParameter fpi)
Definition InitEval.cpp:28
exvector FunExp
Definition SD.h:426
void Integrates(const string &key="", const string &pkey="", int kid=0)
Contours, note that here we need to provide the specific Parameter.
void Evaluate(FeynmanParameter fpi, const string &key="")
Definition InitEval.cpp:507
void CIPrepares(const string &key="")
Prepare for the Contours and Integrates calls .so will be generated, for more detailed exported funct...
static ex set_all(const ex &expr)
Definition BASIC.cpp:147
namespace for Numerical integration with Sector Decomposition method
Definition AsyMB.cpp:10
ex FactorOutX(const ex expr)
Definition Helpers.cpp:278
void Replacement2(exmap &repl)
Definition InitEval.cpp:21
ex expand_ex(ex const &expr_in, std::function< bool(const ex &)> has_func)
the expand like Mathematica
Definition BASIC.cpp:1188
const char * Color_HighLight
Definition BASIC.cpp:248
const iSymbol iEpsilon
int GiNaC_Parallel_Process
Definition Init.cpp:146
void let_op(ex &ex_in, int index1, int index2, const ex item)
update index1-th.index2-th of expression with item
Definition BASIC.cpp:1560
bool xPositive(ex const expr)
check the expr is xPositive, i.e., each x-monomial item is postive
Definition BASIC.cpp:1290
const char * ErrColor
Definition BASIC.cpp:246
const Symbol ep
bool has_w(const ex &e)
Definition BASIC.cpp:2233
const char * WarnColor
Definition BASIC.cpp:247
string now(bool use_date)
date/time string
Definition BASIC.cpp:525
ex w
Definition Init.cpp:93
ex diff_ex(ex const expr, ex const xp, unsigned nth, bool expand)
the differential like Mathematica
Definition BASIC.cpp:1063
int Verbose
Definition Init.cpp:142
void let_op_append(ex &ex_in, const ex item)
append item into expression
Definition BASIC.cpp:1324
string ex2str(const ex &expr)
convert ex to output string, the defalut printer format will be used
Definition BASIC.cpp:715
ex w1
Definition BASIC.h:499
ex w2
Definition BASIC.h:499
ex subs(const ex &e, init_list sl)
Definition BASIC.h:51
wrap parameters for loop integrals
Definition SD.h:66
wrap parameters for generic parameter integrals
Definition SD.h:82
exmap nReplacement
Definition SD.h:85