12 void AIR2F_Save(
string save_dir, exvector air_vec,
const lst & IntFs, vector<IBP*> &ibp_vec) {
13 auto rc = system((
"mkdir -p "+save_dir+
"/AIR2F").c_str());
14 rc = system((
"rm -f "+save_dir+
"/AIR2F/*.gar > /dev/null").c_str());
16 garWrite(air_vec[idx], save_dir+
"/AIR2F/air-"+to_string(idx)+
".gar");
21 garWrite(ibp_vec[idx]->TO(), save_dir+
"/AIR2F/ibp-"+to_string(idx)+
".gar");
26 oss << air_vec.size() <<
" " << ibp_vec.size() <<
" " << IntFs.nops() << endl;
28 oss << IntFs.op(0).op(1).nops() << endl;
29 for(
auto const & f : IntFs) {
31 for(
auto const & n : f.op(1)) oss <<
" " << n;
35 auto oss_str = oss.str();
36 fstream ofs(save_dir+
"/AIR2F.gar", fstream::out);
37 ofs.write(oss_str.c_str(), oss_str.size());
41 void AIR2F_Get(
string save_dir, exvector &air_vec, lst &IntFs, vector<IBP*> &ibp_vec,
int IBPmethod) {
42 fstream ifs(save_dir+
"/AIR2F.gar", fstream::in);
43 size_t nair, nibp, nf;
44 ifs >> nair >> nibp >> nf;
50 for(
int i=0; i<nf; i++) {
54 for(
int j=0; j<nn; j++) {
58 IntFs.append(F(pn, ns));
62 for(
int idx=0; idx<nair; idx++) {
63 air_vec[idx] =
garRead(save_dir+
"/AIR2F/air-"+to_string(idx)+
".gar");
67 for(
int idx=0; idx<nibp; idx++) {
69 if(IBPmethod==0) ibp =
new IBP();
70 else if(IBPmethod==1) ibp =
new FIRE();
71 else if(IBPmethod==2) ibp =
new KIRA();
72 else if(IBPmethod==3) ibp =
new UKIRA();
74 ex ibp_from =
garRead(save_dir+
"/AIR2F/ibp-"+to_string(idx)+
".gar");
87 ex
F2ex(
const ex & expr_in) {
90 if(!e.has(F(
w1,
w2)))
return e;
91 else if(e.match(F(
w1,
w2))) {
95 for(
int i=0; i<ps.nops(); i++) res *= pow(ps.op(i), ex(0)-ns.op(i));
97 }
else return e.map(
self);
114 lst lmom = ex_to<lst>(aio.
Internal);
115 lst emom = ex_to<lst>(aio.
External);
119 if(IBPmethod==1) wdir = aio.
SaveDir +
"/FIRE";
120 else if(IBPmethod==2) wdir = aio.
SaveDir +
"/KIRA";
121 else if(IBPmethod==3) wdir = aio.
SaveDir +
"/UKIRA";
123 wdir = to_string(getpid());
124 if(IBPmethod==1) wdir = wdir +
"_FIRE";
125 else if(IBPmethod==2) wdir = wdir +
"_KIRA";
126 else if(IBPmethod==3) wdir = wdir +
"_UKIRA";
130 vector<IBP*> ibp_vec;
132 if(
Verbose > 1) cout <<
PRE <<
"\\--Reading AIR2F" << flush;
133 AIR2F_Get(aio.
SaveDir, air_vec, IntFs, ibp_vec, IBPmethod);
134 for(
auto ibp : ibp_vec) ibp->WorkingDir = wdir;
135 if(
Verbose > 1) cout <<
" @ " <<
now(
false) << endl;
141 if(
Verbose > 1) cout <<
PRE <<
"\\--Reading AP.gar" << flush;
143 if(
Verbose > 1) cout <<
" @ " <<
now(
false) << endl;
145 }
else rc = system((
"mkdir -p "+aio.
SaveDir).c_str());
149 int av_size = air_vec.size();
155 for(
int i=0; i<av_size; i++) {
156 for(
auto cv : ex_to<lst>(air_vec[i])) vset.insert(cv.op(1));
158 exvector vvec(vset.begin(), vset.end());
162 for(
int i=0; i<vvec.size(); i++) v2api[vvec[i]] = i;
163 for(
int i=0; i<av_size; i++) {
165 lst cvs = ex_to<lst>(air_vec[i]);
166 air_vec[i] = av_item;
167 for(
auto & cv : cvs) av_item.append(lst{cv.op(0), v2api[cv.op(1)]});
169 air_vec[i] = av_item;
173 auto ap_vec =
GiNaC_Parallel(vvec.size(), [&vvec,lmom,emom,aio] (
int idx) {
174 auto air = vvec[idx];
175 air = Apart(air,lmom,emom,aio.smap);
176 air = air.subs(ApartIR(1,w)==aio.apart1);
177 air = collect_lst(air, ApartIR(w1,w2), o_flintf);
183 for(
auto cvs : ap_vec)
for(
auto cv : cvs)
if(is_a<matrix>(cv.op(1).op(0))) ap_set.insert(cv.op(1));
184 exvector ap_ir_vec(ap_set.begin(), ap_set.end());
191 ap_vec =
GiNaC_Parallel(ap_vec.size(), [&aio,&ap_vec,&ap_rules] (
int idx) {
192 auto const & cvs = ap_vec[idx];
195 for(auto const & cv : cvs) {
196 auto fi = ap_rules.find(cv.op(1));
197 if(fi==ap_rules.end()) res += cv.op(0) * cv.op(1);
198 else res += cv.op(0) * fi->second;
201 for(auto cv : cvs) res += cv.op(0) * ApartIRC(cv.op(1));
206 int vn = ap_vec.size();
207 for(
int idx=0; idx<vn; idx++) {
208 lst cvs = ex_to<lst>(ap_vec[idx]);
211 for(
auto const & cv : cvs) {
212 auto fi = ap_rules.find(cv.op(1));
213 if(fi==ap_rules.end()) res += cv.op(0) * cv.op(1);
214 else res += cv.op(0) * fi->second;
217 for(
auto cv : cvs) res += cv.op(0) *
ApartIRC(cv.op(1));
227 lst cvs = ex_to<lst>(air_vec[idx]);
229 for(
auto const & cv : cvs) {
230 int idx = ex_to<numeric>(cv.op(1)).to_int();
231 res += cv.op(0) * ap_vec[idx];
249 auto ret =
GiNaC_Parallel(air_vec.size(), [&air_vec](
int idx)->ex {
250 auto air = air_vec[idx];
252 find(air, ApartIR(w1,w2), airs);
254 for(auto item : airs) ret.append(item);
258 for(
auto airs : ret)
for(
auto air : ex_to<lst>(airs)) intg.insert(air);
259 AIR = exvector(intg.begin(), intg.end());
268 repls.append(
w*kv.first ==
w*kv.second);
269 repls.append(kv.first == kv.second);
273 for(
auto li : lmom) {
274 if(is_a<Vector>(li)) loops.append(ex_to<Vector>(li).name);
275 else loops.append(li);
277 for(
auto li : emom) {
278 if(is_a<Vector>(li)) exts.append(ex_to<Vector>(li).name);
279 else exts.append(li);
282 if(
Verbose>0) cout <<
PRE <<
"\\--Prepare " <<
WHITE <<
"IBP" <<
RESET <<
" reduction @ " <<
now(
false) << flush;
285 std::map<ex, IBP*, ex_is_less> p2IBP;
287 int ntot = AIR.size();
288 for(
int i=0; i<ntot; i++) {
289 if(
Verbose>0 && (((i+1)%1000)==0 || i+1==ntot)) {
290 cout <<
"\r \r" << flush;
291 cout <<
PRE <<
"\\--Prepare " <<
WHITE <<
"IBP" <<
RESET <<
" reduction [" << (i+1) <<
"/" << ntot <<
"] @ " <<
now(
false) << flush;
293 auto const & ir = AIR[i];
294 auto mat = ex_to<matrix>(ir.op(0));
295 auto vars = ex_to<lst>(ir.op(1));
297 int nrow = mat.rows();
299 for(
int c=0; c<mat.cols(); c++) {
301 for(
int r=0; r<nrow-2; r++) pc += mat(r,c) * vars.op(r);
304 ex nc = ex(0)-mat(nrow-1,c);
308 if(ncn==-1) den_tot++;
309 pns.append(lst{ ncn, pc, nc });
311 bool pn_sector =
false;
312 if(aio.pn_sector>0 && den_tot>=aio.pn_sector) pn_sector =
true;
314 for(
int i=0; i<pns.nops(); i++) pns.let_op(i) = lst{ pns.op(i).op(1), pns.op(i).op(2) };
318 for(
int i=0; i<pns.nops(); i++) pns.let_op(i) = lst{ pns.op(i).op(1), pns.op(i).op(2) };
321 int nCut = aio.Cut.nops();
325 if(aio.CutFirst)
for(
auto cut : cuts) pns.prepend(lst{
SP2sp(cut), 1 });
326 else for(
auto cut : cuts) pns.append(lst{
SP2sp(cut), 1 });
330 for(
auto item : pns) {
331 props.append(item.op(0));
332 ns.append(item.op(1));
338 for(
int i=0; i<ns.nops(); i++) nss.append(ns.op(i)>0 ? 1 : 0);
339 key = lst{props,nss};
342 auto kv = p2IBP.find(key);
343 if(kv==p2IBP.end()) {
345 if(IBPmethod==0) ibp =
new IBP();
346 else if(IBPmethod==1) ibp =
new FIRE();
347 else if(IBPmethod==2) ibp =
new KIRA();
348 else if(IBPmethod==3) ibp =
new UKIRA();
354 p2IBP.insert(make_pair(key,ibp));
355 ibp->Propagator = props;
356 ibp->Internal = loops;
357 ibp->External = exts;
358 ibp->Replacement = repls;
359 if(aio.ISP.nops()>0)
for(
auto item : aio.ISP) ibp->ISP.append(
SP2sp(item));
360 if(aio.DSP.nops()>0) {
361 for(
auto item : aio.DSP) {
362 lst
sp = ex_to<lst>(item);
363 if(is_a<Vector>(
sp.op(0)))
sp.let_op(0) = (ex_to<Vector>(
sp.op(0)).name);
364 if(is_a<Vector>(
sp.op(1)))
sp.let_op(1) = (ex_to<Vector>(
sp.op(1)).name);
370 for(
auto const & item : ns) {
371 if(item>0) sector.append(1);
372 else sector.append(0);
374 ibp->SECTOR = sector;
376 ibp->WorkingDir = wdir;
377 ibp->ProblemNumber = pn;
380 if(aio.CutFirst)
for(
int i=0; i<nCut; i++) ibp->Cut.append(i+1);
381 else for(
int i=0; i<nCut; i++) ibp->Cut.append(nCut-i);
383 ibp_vec.push_back(ibp);
384 ibp->Integral.append(ns);
385 AIR2F[AIR[i]] = F(ibp->ProblemNumber, ns);
387 IBP* ibp = kv->second;
388 ibp->Integral.append(ns);
389 AIR2F[AIR[i]] = F(ibp->ProblemNumber, ns);
394 if(
Verbose>0) cout <<
PRE <<
"\\--Total Ints/Pros: " <<
WHITE << ntot <<
"/" << ibp_vec.size() <<
RESET <<
" @ " <<
now(
false) << endl;
399 auto int_fr =
FindRules(ibp_vec,
false, aio.UF);
400 IntFs = int_fr.second;
403 air_vec =
GiNaC_Parallel(air_vec.size(), [&air_vec,&AIR2F,&int_fr] (
int idx) {
404 auto air = air_vec[idx];
405 air = air.subs(AIR2F,nopat);
406 air = air.subs(int_fr.first,nopat);
407 air = collect_ex(air, F(w1,w2), o_flint);
410 if(aio.SaveDir !=
"") AIR2F_Save(aio.SaveDir, air_vec, IntFs, ibp_vec);
415 MapFunction _F2ex([&ibp_vec,aio](
const ex &e, MapFunction &
self)->ex {
416 if(!e.has(F(
w1,
w2)))
return e;
417 else if(e.match(F(
w1,
w2))) {
418 int pn = ex_to<numeric>(e.op(0)).to_int();
419 auto pso = ex_to<lst>(ibp_vec[pn-1]->Propagator);
420 auto nso = ex_to<lst>(e.op(1));
422 for(
int i=0; i<pso.nops(); i++) {
423 if(!aio.keep0F && nso.op(i).is_zero())
continue;
424 ps.append(pso.op(i));
425 ns.append(nso.op(i));
428 }
else return e.map(
self);
432 air_vec =
GiNaC_Parallel(air_vec.size(), [&air_vec,&_F2ex](
int idx)->ex {
433 auto res = air_vec[idx];
436 for(
auto fp : ibp_vec)
delete fp;
437 if(aio.SaveDir ==
"") rc = system((
"rm -rf "+wdir).c_str());
442 if(aio.SaveDir !=
"" &&
file_exists(aio.SaveDir+
"/Rules.gar")) {
446 vector<IBP*> ibp_vec_re;
448 map<int,lst> pn_ints_map;
449 for(
auto item : IntFs) {
450 int pn = ex_to<numeric>(item.op(0)).to_int();
451 pn_ints_map[pn].append(item.op(1));
455 for(
auto pi : pn_ints_map) {
456 auto ibp = ibp_vec[pi.first-1];
457 ibp->Integral = pi.second;
458 nints += ibp->Integral.nops();
459 ibp_vec_re.push_back(ibp);
462 if(
Verbose>0) cout <<
PRE <<
"\\--Refined Ints/Pros: " <<
WHITE << nints <<
"/" << ibp_vec_re.size() <<
RESET <<
" @ " <<
now(
false) << endl;
467 auto pRes =
GiNaC_Parallel(ibp_vec_re.size(), [&ibp_vec_re](
int idx)->ex {
468 ibp_vec_re[idx]->Export();
469 auto ret = lst{ ibp_vec_re[idx]->IsAlwaysZero ? 1 : 0, ibp_vec_re[idx]->Rules };
472 for(
int i=0; i<ibp_vec_re.size(); i++) {
473 ibp_vec_re[i]->IsAlwaysZero = (pRes[i].op(0)==1 ? true :
false);
474 ibp_vec_re[i]->Rules = ex_to<lst>(pRes[i].op(1));
477 int nproc = aio.NIBP;
478 if(nproc<1) nproc = 8;
480 if(nproc<1) nproc = 1;
481 size_t nibp = ibp_vec_re.size();
485 #pragma omp parallel for num_threads(nproc) schedule(dynamic, 1)
486 for(
int pi=0; pi<nibp; pi++) {
490 cout <<
"\r \r" << flush;
491 cout <<
PRE <<
"\\--" <<
WHITE <<
"FIRE" <<
RESET <<
" Reduction [" << (++cproc) <<
"/" << nibp <<
"] " << flush;
494 ibp_vec_re[pi]->Run();
496 if(
Verbose>1 && nibp>0) cout <<
"@" <<
now(
false) << endl;
502 ibp_vec_re[idx]->Run();
506 for(
int pi=0; pi<nibp; pi++) {
507 if(
Verbose>1) cout <<
"\r \r" <<
PRE <<
"\\--" <<
WHITE <<
"FIRE" <<
RESET <<
" Reduction [" << (++cproc) <<
"/" << nibp <<
"] " << flush;
508 ibp_vec_re[pi]->Run();
510 if(
Verbose>1 && nibp>0) cout <<
"@" <<
now(
false) << endl;
514 if(ibp_vec_re.size()>100) {
515 auto ret =
GiNaC_Parallel(ibp_vec_re.size(), [&ibp_vec_re,wdir](
int idx)->ex {
516 ibp_vec_re[idx]->Import();
517 return ibp_vec_re[idx]->TO();
519 for(
int i=0; i<ibp_vec_re.size(); i++) ibp_vec_re[i]->FROM(ret[i]);
522 for(
auto item : ibp_vec_re) {
523 if(
Verbose>1) cout <<
"\r \r" <<
PRE <<
"\\--" <<
WHITE <<
"FIRE" <<
RESET <<
" Import [" << (++cproc) <<
"/" << ibp_vec_re.size() <<
"] " << flush;
526 if(
Verbose>1 && ibp_vec_re.size()>0) cout <<
"@" <<
now(
false) << endl;
530 if(aio.SaveDir ==
"") rc = system((
"rm -rf "+wdir).c_str());
531 }
else if(IBPmethod==2 || IBPmethod==3) {
532 for(
auto ibp : ibp_vec_re) ibp->Reduce();
533 if(aio.SaveDir ==
"") rc = system((
"rm -rf "+wdir).c_str());
537 exmap miRules =
FindRules(ibp_vec_re,
true, aio.UF).first;
539 if(aio.SaveDir !=
"") rc = system((
"mkdir -p "+aio.SaveDir+
"/Rules").c_str());
540 auto rules_vec =
GiNaC_Parallel(ibp_vec_re.size(), [&ibp_vec_re,&miRules,&aio](
int idx)->ex {
541 lst rules = ex_to<lst>(ibp_vec_re[idx]->Rules);
543 for(auto ri : rules) res.append(lst {
545 collect_ex(ri.op(1).subs(miRules,nopat),F(w1,w2),o_flint)
547 for(
auto mi : ibp_vec_re[idx]->MIntegral) {
548 auto fi = miRules.find(mi);
549 if(fi!=miRules.end()) res.append(lst{ mi, fi->second });
551 auto pn = ibp_vec_re[idx]->ProblemNumber;
552 if(aio.SaveDir !=
"") {
553 garWrite(aio.SaveDir+
"/Rules/"+to_string(pn)+
".gar", res);
557 if(aio.SaveDir !=
"") {
558 garWrite(aio.SaveDir+
"/Rules.gar", 1);
560 for(
auto rs : rules_vec) {
561 for(
auto ri : rs)
if(ri.op(0)!=ri.op(1)) ibpRules[ri.op(0)] = ri.op(1);
570 GiNaC_Parallel(air_vec.size(), [&air_vec,&ibpRules,&_F2ex,&aio](
int idx)->ex {
571 ex res = air_vec[idx];
573 if(aio.SaveDir !=
"") {
575 find(res, F(w1,w2), fs);
577 for(auto fi : fs) pns.insert(fi.op(0));
579 auto rs = ex_to<lst>(garRead(aio.SaveDir+
"/Rules/"+ex2str(pn)+
".gar"));
580 for(auto ri : rs) if(ri.op(0)!=ri.op(1)) rules[ri.op(0)] = ri.op(1);
582 }
else rules = ibpRules;
583 res = res.subs(rules,
nopat);
584 if(aio.pat.nops()>0) {
585 auto cvs = collect_lst(res, aio.pat);
590 if(aio.cv!=nullptr) {
591 auto _cv = aio.cv(c,v);
601 for(
auto fp : ibp_vec)
delete fp;
614 void ApartIBP(exvector &air_vec,
int IBPmethod,
const lst & loops,
const lst & exts,
const lst & cut_props,
615 std::function<lst(
const IBP &,
const ex &)> uf) {
623 if(cut_props.nops()>0) {
624 for(
auto p1 : loops) {
625 for(
auto p2 : loops) aio.
CSP.append(
SP(p1,p2));
626 for(
auto p2 : exts) aio.
CSP.append(
SP(p1,p2));
631 for(
auto li : loops) aio.
smap[
SP(li)] = 1;
IBP base class for IBP reduction.
class to wrap map_function of GiNaC
ex SP2sp(const ex &exin)
convert SP(a,b) to sp(a,b)
ex sp(const ex &a, const ex &b)
translated the vector dot a.b to a*b, useful in SecDec
exmap sp_map()
the SP_map with SP(a,b) replaced to sp(a,b)
ex ApartIRC(const ex &expr_in)
complete the ApartIR elements
void ApartIBP(exvector &air_vec, AIOption aio)
perform IBP reduction on the Aparted input
ex collect_ex(ex const &expr_in, std::function< bool(const ex &)> has_func, int opt)
the collect function like Mathematica
void garRead(const string &garfn, map< string, ex > &resMap)
garRead from file, and output in a map
map< string, int > GiNaC_Parallel_NB
bool file_exists(string fn)
pair< exmap, lst > FindRules(vector< IBP * > fs, bool mi, std::function< lst(const IBP &, const ex &)> uf)
Find Rules for Integral or Master Integral.
string now(bool use_date)
date/time string
ex F2ex(const ex &expr_in)
convert F(ps, ns) to normal ex, ns is like FIRE convention
int CpuCores()
return the cpu cores using OpenMP
map< string, int > GiNaC_Parallel_NP
void ApartIBP(exvector &air_vec, int IBPmethod, const lst &loops, const lst &exts, const lst &cut_props, std::function< lst(const IBP &, const ex &)> uf)
perform IBP reduction on the Aparted input
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}, ... }
exmap ApartRules(const exvector &airs, bool irc)
void sort_lst(lst &ilst, bool less=true)
sort the list in less order, or the reverse
void garWrite(const string &garfn, const map< string, ex > &resMap)
garWrite to write the string-key map to the archive
exvector GiNaC_Parallel(int ntotal, std::function< ex(int)> f, const string &key)
GiNaC Parallel Evaluation using fork.
ex SP(const ex &a, bool use_map=false)
std::function< lst(const IBP &, const ex &)> UF