HepLib
ApartIBP.cpp
Go to the documentation of this file.
1 
6 #include "HEP.h"
7 
8 namespace HepLib {
9 
10  namespace {
11 
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());
15  GiNaC_Parallel(air_vec.size(), [&air_vec,&save_dir](int idx) {
16  garWrite(air_vec[idx], save_dir+"/AIR2F/air-"+to_string(idx)+".gar");
17  return 0;
18  }, "AIR2F_AIR");
19 
20  GiNaC_Parallel(ibp_vec.size(), [&ibp_vec,&save_dir](int idx) {
21  garWrite(ibp_vec[idx]->TO(), save_dir+"/AIR2F/ibp-"+to_string(idx)+".gar");
22  return 0;
23  }, "AIR2F_IBP");
24 
25  ostringstream oss;
26  oss << air_vec.size() << " " << ibp_vec.size() << " " << IntFs.nops() << endl;
27  if(IntFs.nops()>0) {
28  oss << IntFs.op(0).op(1).nops() << endl;
29  for(auto const & f : IntFs) {
30  oss << f.op(0);
31  for(auto const & n : f.op(1)) oss << " " << n;
32  oss << endl;
33  }
34  }
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());
38  ofs.close();
39  }
40 
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;
45 
46  IntFs.remove_all();
47  if(nf>0) {
48  int nn;
49  ifs >> nn;
50  for(int i=0; i<nf; i++) {
51  int pn, ni;
52  ifs >> pn;
53  lst ns;
54  for(int j=0; j<nn; j++) {
55  ifs >> ni;
56  ns.append(ni);
57  }
58  IntFs.append(F(pn, ns));
59  }
60  }
61 
62  for(int idx=0; idx<nair; idx++) {
63  air_vec[idx] = garRead(save_dir+"/AIR2F/air-"+to_string(idx)+".gar");
64  }
65 
66  ibp_vec.resize(nibp);
67  for(int idx=0; idx<nibp; idx++) {
68  IBP* ibp;
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();
73  else ibp = new IBP();
74  ex ibp_from = garRead(save_dir+"/AIR2F/ibp-"+to_string(idx)+".gar");
75  ibp->FROM(ibp_from);
76  ibp_vec[idx] = ibp;
77  }
78  }
79 
80  }
81 
87  ex F2ex(const ex & expr_in) {
88  ex ret = expr_in;
89  ret = MapFunction([](const ex & e, MapFunction &self)->ex{
90  if(!e.has(F(w1,w2))) return e;
91  else if(e.match(F(w1, w2))) {
92  auto ps = e.op(0);
93  auto ns = e.op(1);
94  ex res = 1;
95  for(int i=0; i<ps.nops(); i++) res *= pow(ps.op(i), ex(0)-ns.op(i));
96  return res;
97  } else return e.map(self);
98  })(ret);
99  return ret;
100  }
101 
109  void ApartIBP(exvector &air_vec, AIOption aio) {
110  if(aio.smap.size()<1) aio.init_smap();
111  int IBPmethod = aio.IBPmethod;
112  int rc;
113 
114  lst lmom = ex_to<lst>(aio.Internal);
115  lst emom = ex_to<lst>(aio.External);
116 
117  string wdir;
118  if(aio.SaveDir != "") {
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";
122  } else {
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";
127  }
128 
129  lst IntFs;
130  vector<IBP*> ibp_vec;
131  if(aio.SaveDir != "" && file_exists(aio.SaveDir+"/AIR2F.gar")) {
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; // update working directory
135  if(Verbose > 1) cout << " @ " << now(false) << endl;
136  goto AIR2F_Done;
137  }
138 
139  if(aio.SaveDir != "") {
140  if(file_exists(aio.SaveDir+"/AP.gar")) {
141  if(Verbose > 1) cout << PRE << "\\--Reading AP.gar" << flush;
142  garRead(air_vec, aio.SaveDir+"/AP.gar");
143  if(Verbose > 1) cout << " @ " << now(false) << endl;
144  goto Apart_Done;
145  } else rc = system(("mkdir -p "+aio.SaveDir).c_str());
146  }
147 
148  if(true) {
149  int av_size = air_vec.size();
150  air_vec = GiNaC_Parallel(av_size, [air_vec,lmom] (int idx) {
151  return collect_lst(air_vec[idx], lmom, o_flint);
152  }, "ApPre");
153 
154  exset vset;
155  for(int i=0; i<av_size; i++) {
156  for(auto cv : ex_to<lst>(air_vec[i])) vset.insert(cv.op(1));
157  }
158  exvector vvec(vset.begin(), vset.end());
159  vset.clear();
160 
161  exmap v2api;
162  for(int i=0; i<vvec.size(); i++) v2api[vvec[i]] = i;
163  for(int i=0; i<av_size; i++) {
164  lst av_item;
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)]});
168  cvs.remove_all();
169  air_vec[i] = av_item;
170  }
171  v2api.clear();
172 
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);
178  return air;
179  }, "Apart");
180  vvec.clear();
181 
182  exset ap_set;
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());
185  ap_set.clear();
186  exmap ap_rules;
187  if(aio.ap_rules) ap_rules = ApartRules(ap_ir_vec); // including ApartIRC
188  ap_ir_vec.clear();
189 
190  if(false) { // Parallel Version
191  ap_vec = GiNaC_Parallel(ap_vec.size(), [&aio,&ap_vec,&ap_rules] (int idx) {
192  auto const & cvs = ap_vec[idx];
193  ex res = 0;
194  if(aio.ap_rules) {
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;
199  }
200  } else {
201  for(auto cv : cvs) res += cv.op(0) * ApartIRC(cv.op(1));
202  }
203  return res;
204  }, "ApRule");
205  } else {
206  int vn = ap_vec.size();
207  for(int idx=0; idx<vn; idx++) {
208  lst cvs = ex_to<lst>(ap_vec[idx]);
209  ex res = 0;
210  if(aio.ap_rules) {
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;
215  }
216  } else {
217  for(auto cv : cvs) res += cv.op(0) * ApartIRC(cv.op(1));
218  }
219  ap_vec[idx] = res;
220  }
221  }
222  ap_rules.clear();
223 
224  if(GiNaC_Parallel_NP.find("ApPost")==GiNaC_Parallel_NP.end() && CpuCores()>8) GiNaC_Parallel_NP["ApPost"] = 8;
225  if(GiNaC_Parallel_NB.find("ApPost")==GiNaC_Parallel_NB.end() && CpuCores()>100) GiNaC_Parallel_NB["ApPost"] = 100;
226  air_vec = GiNaC_Parallel(av_size, [&air_vec,&ap_vec] (int idx) {
227  lst cvs = ex_to<lst>(air_vec[idx]);
228  ex res = 0;
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];
232  }
233  res = collect_ex(res, ApartIR(w1,w2), o_flint);
234  return res; // air_vec updated to ApartIR
235  }, "ApPost");
236  ap_vec.clear();
237 
238  if(aio.SaveDir != "") {
239  garWrite(air_vec,aio.SaveDir+"/AP.gar");
240  garRead(air_vec,aio.SaveDir+"/AP.gar");
241  }
242  }
243  Apart_Done: ;
244 
245  if(true) {
246 
247  exvector AIR;
248  if(true) {
249  auto ret = GiNaC_Parallel(air_vec.size(), [&air_vec](int idx)->ex {
250  auto air = air_vec[idx];
251  exset airs;
252  find(air, ApartIR(w1,w2), airs);
253  lst ret;
254  for(auto item : airs) ret.append(item);
255  return ret;
256  }, "ApIRC");
257  exset intg;
258  for(auto airs : ret) for(auto air : ex_to<lst>(airs)) intg.insert(air);
259  AIR = exvector(intg.begin(), intg.end());
260  }
261 
262  for(auto sp : aio.CSP) SP_map.erase(sp);
263  // from here, Vector will be replaced by its name Symbol
264 
265  lst repls;
266  auto sps = sp_map();
267  for(auto kv : sps) {
268  repls.append(w*kv.first == w*kv.second);
269  repls.append(kv.first == kv.second);
270  }
271 
272  lst loops, exts; // to match FIRE notation, not Vector, just Symbol
273  for(auto li : lmom) {
274  if(is_a<Vector>(li)) loops.append(ex_to<Vector>(li).name);
275  else loops.append(li);
276  }
277  for(auto li : emom) {
278  if(is_a<Vector>(li)) exts.append(ex_to<Vector>(li).name);
279  else exts.append(li);
280  }
281 
282  if(Verbose>0) cout << PRE << "\\--Prepare " << WHITE << "IBP" << RESET << " reduction @ " << now(false) << flush;
283 
284  exmap AIR2F;
285  std::map<ex, IBP*, ex_is_less> p2IBP;
286  int pn=1;
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;
292  }
293  auto const & ir = AIR[i];
294  auto mat = ex_to<matrix>(ir.op(0));
295  auto vars = ex_to<lst>(ir.op(1));
296  lst pns;
297  int nrow = mat.rows();
298  int den_tot = 0;
299  for(int c=0; c<mat.cols(); c++) {
300  ex pc = 0;
301  for(int r=0; r<nrow-2; r++) pc += mat(r,c) * vars.op(r);
302  pc += mat(nrow-2,c);
303  pc = SP2sp(pc);
304  ex nc = ex(0)-mat(nrow-1,c);
305  int ncn;
306  if(nc>0) ncn = -1;
307  else ncn = 1;
308  if(ncn==-1) den_tot++;
309  pns.append(lst{ ncn, pc, nc }); // note the convension, ncn just for sorting
310  }
311  bool pn_sector = false;
312  if(aio.pn_sector>0 && den_tot>=aio.pn_sector) pn_sector = true;
313  if(!pn_sector) { // back to original format
314  for(int i=0; i<pns.nops(); i++) pns.let_op(i) = lst{ pns.op(i).op(1), pns.op(i).op(2) };
315  }
316  sort_lst(pns);
317  if(pn_sector) { // back to original format
318  for(int i=0; i<pns.nops(); i++) pns.let_op(i) = lst{ pns.op(i).op(1), pns.op(i).op(2) };
319  }
320 
321  int nCut = aio.Cut.nops();
322  if(nCut>0) {
323  ex cuts = aio.Cut;
324  cuts = cuts.subs(SP_map,nopat);
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 });
327  }
328 
329  lst props, ns;
330  for(auto item : pns) {
331  props.append(item.op(0));
332  ns.append(item.op(1));
333  }
334 
335  ex key = props;
336  if(pn_sector) {
337  lst nss;
338  for(int i=0; i<ns.nops(); i++) nss.append(ns.op(i)>0 ? 1 : 0);
339  key = lst{props,nss};
340  }
341 
342  auto kv = p2IBP.find(key);
343  if(kv==p2IBP.end()) {
344  IBP* ibp;
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();
349  else {
350  ibp = new IBP();
351  IBPmethod = 0;
352  }
353 
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);
365  ibp->DSP.append(sp);
366  }
367  }
368  if(pn_sector) {
369  lst sector;
370  for(auto const & item : ns) {
371  if(item>0) sector.append(1);
372  else sector.append(0);
373  }
374  ibp->SECTOR = sector;
375  }
376  ibp->WorkingDir = wdir;
377  ibp->ProblemNumber = pn;
378  pn++;
379  if(nCut>0) {
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);
382  }
383  ibp_vec.push_back(ibp);
384  ibp->Integral.append(ns);
385  AIR2F[AIR[i]] = F(ibp->ProblemNumber, ns);
386  } else {
387  IBP* ibp = kv->second;
388  ibp->Integral.append(ns);
389  AIR2F[AIR[i]] = F(ibp->ProblemNumber, ns);
390  }
391  }
392  if(Verbose>0) cout << endl;
393 
394  if(Verbose>0) cout << PRE << "\\--Total Ints/Pros: " << WHITE << ntot << "/" << ibp_vec.size() << RESET << " @ " << now(false) << endl;
395 
396  if(true) {
397  //vector<IBP*> ibp_vec2;
398  //for(auto ibp : ibp_vec) ibp_vec2.push_back(ibp);
399  auto int_fr = FindRules(ibp_vec, false, aio.UF);
400  IntFs = int_fr.second;
401  if(GiNaC_Parallel_NP.find("AIR2F")==GiNaC_Parallel_NP.end() && CpuCores()>8) GiNaC_Parallel_NP["AIR2F"] = 8;
402  if(GiNaC_Parallel_NB.find("AIR2F")==GiNaC_Parallel_NB.end() && CpuCores()>100) GiNaC_Parallel_NB["AIR2F"] = 100;
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);
408  return air;
409  }, "AIR2F");
410  if(aio.SaveDir != "") AIR2F_Save(aio.SaveDir, air_vec, IntFs, ibp_vec);
411  }
412  }
413  AIR2F_Done: ;
414 
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));
421  lst ps, ns;
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));
426  }
427  return F(ps,ns);
428  } else return e.map(self);
429  });
430 
431  if(IBPmethod==0) { // no IBP reduction
432  air_vec = GiNaC_Parallel(air_vec.size(), [&air_vec,&_F2ex](int idx)->ex {
433  auto res = air_vec[idx];
434  return _F2ex(res);
435  }, "F2F");
436  for(auto fp : ibp_vec) delete fp;
437  if(aio.SaveDir == "") rc = system(("rm -rf "+wdir).c_str());
438  return;
439  }
440 
441  exmap ibpRules; // IBP rules for problem pn
442  if(aio.SaveDir != "" && file_exists(aio.SaveDir+"/Rules.gar")) {
443  goto Rules_Done;
444  }
445  if(true) {
446  vector<IBP*> ibp_vec_re;
447  if(true) {
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));
452  }
453 
454  int nints = 0;
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);
460  }
461 
462  if(Verbose>0) cout << PRE << "\\--Refined Ints/Pros: " << WHITE << nints << "/" << ibp_vec_re.size() << RESET << " @ " << now(false) << endl;
463  }
464 
465  if(IBPmethod==1) {
466  //if(GiNaC_Parallel_NB.find("Expo")==GiNaC_Parallel_NB.end()) GiNaC_Parallel_NB["Expo"] = 1;
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 };
470  return ret;
471  }, "Expo");
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));
475  }
476 
477  int nproc = aio.NIBP;
478  if(nproc<1) nproc = 8;
479  int cproc = 0;
480  if(nproc<1) nproc = 1;
481  size_t nibp = ibp_vec_re.size();
482 
483  //#define using_openMP
484  #ifdef using_openMP
485  #pragma omp parallel for num_threads(nproc) schedule(dynamic, 1)
486  for(int pi=0; pi<nibp; pi++) {
487  if(Verbose>1) {
488  #pragma omp critical
489  {
490  cout << "\r \r" << flush;
491  cout << PRE << "\\--" << WHITE << "FIRE" << RESET << " Reduction [" << (++cproc) << "/" << nibp << "] " << flush;
492  }
493  }
494  ibp_vec_re[pi]->Run();
495  }
496  if(Verbose>1 && nibp>0) cout << "@" << now(false) << endl;
497  #else
498  if(nproc>1) {
499  if(GiNaC_Parallel_NP.find("FIRE")==GiNaC_Parallel_NB.end()) GiNaC_Parallel_NP["FIRE"] = nproc;
500  if(GiNaC_Parallel_NB.find("FIRE")==GiNaC_Parallel_NB.end()) GiNaC_Parallel_NB["FIRE"] = 1;
501  GiNaC_Parallel(nibp, [&ibp_vec_re](int idx)->ex {
502  ibp_vec_re[idx]->Run();
503  return 0;
504  }, "FIRE");
505  } else {
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();
509  }
510  if(Verbose>1 && nibp>0) cout << "@" << now(false) << endl;
511  }
512  #endif
513 
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();
518  }, "Impo");
519  for(int i=0; i<ibp_vec_re.size(); i++) ibp_vec_re[i]->FROM(ret[i]);
520  } else {
521  cproc = 0;
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;
524  item->Import();
525  }
526  if(Verbose>1 && ibp_vec_re.size()>0) cout << "@" << now(false) << endl;
527  }
528  //IBP::ReShare(ibp_vec_re);
529 
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());
534  }
535 
536  // Find Rules in MIs
537  exmap miRules = FindRules(ibp_vec_re, true, aio.UF).first;
538  if(true) { // scope for ret
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);
542  lst res;
543  for(auto ri : rules) res.append(lst {
544  ri.op(0),
545  collect_ex(ri.op(1).subs(miRules,nopat),F(w1,w2),o_flint)
546  });
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 });
550  }
551  auto pn = ibp_vec_re[idx]->ProblemNumber;
552  if(aio.SaveDir != "") {
553  garWrite(aio.SaveDir+"/Rules/"+to_string(pn)+".gar", res);
554  return 0;
555  } else return res;
556  }, "FR2MI");
557  if(aio.SaveDir != "") {
558  garWrite(aio.SaveDir+"/Rules.gar", 1);
559  } else {
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);
562  }
563  }
564  }
565  }
566  Rules_Done: ;
567 
568  if(GiNaC_Parallel_NP.find("F2MI")==GiNaC_Parallel_NP.end() && CpuCores()>16) GiNaC_Parallel_NP["F2MI"] = 16;
569  air_vec =
570  GiNaC_Parallel(air_vec.size(), [&air_vec,&ibpRules,&_F2ex,&aio](int idx)->ex {
571  ex res = air_vec[idx];
572  exmap rules;
573  if(aio.SaveDir != "") {
574  exset fs;
575  find(res, F(w1,w2), fs);
576  exset pns;
577  for(auto fi : fs) pns.insert(fi.op(0));
578  for(auto pn : pns) {
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);
581  }
582  } else rules = ibpRules;
583  res = res.subs(rules,nopat);
584  if(aio.pat.nops()>0) {
585  auto cvs = collect_lst(res, aio.pat);
586  res = 0;
587  for(auto cv : cvs) {
588  auto c = cv.op(0);
589  auto v = cv.op(1);
590  if(aio.cv!=nullptr) {
591  auto _cv = aio.cv(c,v);
592  c = _cv.op(0);
593  v = _cv.op(1);
594  }
595  res += c * v;
596  }
597  }
598  return _F2ex(res);
599  }, "F2MI");
600 
601  for(auto fp : ibp_vec) delete fp;
602  }
603 
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) {
616 
617  AIOption aio;
618  aio.IBPmethod = IBPmethod;
619  //aio.pn_sector = 4;
620  aio.Internal = loops;
621  aio.External = exts;
622  aio.Cut = cut_props;
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));
627  }
628  aio.CSP.sort();
629  aio.CSP.unique();
630  }
631  for(auto li : loops) aio.smap[SP(li)] = 1;
632  aio.UF = uf;
633  ApartIBP(air_vec, aio);
634  }
635 
636 }
#define WHITE
Definition: BASIC.h:87
#define RESET
Definition: BASIC.h:79
HEP header file.
IBP base class for IBP reduction.
Definition: IBP.h:24
class to wrap map_function of GiNaC
Definition: BASIC.h:632
HepLib namespace.
Definition: BASIC.cpp:17
ex SP2sp(const ex &exin)
convert SP(a,b) to sp(a,b)
Definition: Pair.cpp:302
ex sp(const ex &a, const ex &b)
translated the vector dot a.b to a*b, useful in SecDec
Definition: Pair.cpp:237
exmap sp_map()
the SP_map with SP(a,b) replaced to sp(a,b)
Definition: Pair.cpp:318
const int o_flint
Definition: Init.cpp:110
exmap SP_map
Definition: Init.cpp:179
ex ApartIRC(const ex &expr_in)
complete the ApartIR elements
Definition: Apart.cpp:914
void ApartIBP(exvector &air_vec, AIOption aio)
perform IBP reduction on the Aparted input
Definition: ApartIBP.cpp:109
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
void garRead(const string &garfn, map< string, ex > &resMap)
garRead from file, and output in a map
Definition: BASIC.cpp:591
map< string, int > GiNaC_Parallel_NB
Definition: Init.cpp:147
bool file_exists(string fn)
Definition: BASIC.h:289
pair< exmap, lst > FindRules(vector< IBP * > fs, bool mi, std::function< lst(const IBP &, const ex &)> uf)
Find Rules for Integral or Master Integral.
Definition: IBP.cpp:566
string now(bool use_date)
date/time string
Definition: BASIC.cpp:525
ex F2ex(const ex &expr_in)
convert F(ps, ns) to normal ex, ns is like FIRE convention
Definition: ApartIBP.cpp:87
ex w
Definition: Init.cpp:90
int CpuCores()
return the cpu cores using OpenMP
Definition: BASIC.cpp:1820
map< string, int > GiNaC_Parallel_NP
Definition: Init.cpp:144
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
Definition: ApartIBP.cpp:614
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
exmap ApartRules(const exvector &airs, bool irc)
Definition: Apart.cpp:156
int Verbose
Definition: Init.cpp:139
void sort_lst(lst &ilst, bool less=true)
sort the list in less order, or the reverse
Definition: Sort.cpp:79
void garWrite(const string &garfn, const map< string, ex > &resMap)
garWrite to write the string-key map to the archive
Definition: BASIC.cpp:639
ex w1
Definition: BASIC.h:499
unsigned nopat
Definition: Init.cpp:88
exvector GiNaC_Parallel(int ntotal, std::function< ex(int)> f, const string &key)
GiNaC Parallel Evaluation using fork.
Definition: BASIC.cpp:260
ex w2
Definition: BASIC.h:499
string PRE
Definition: Init.cpp:140
ex SP(const ex &a, bool use_map=false)
Definition: Pair.cpp:166
std::function< lst(const IBP &, const ex &)> UF
Definition: HEP.h:538
exmap smap
Definition: HEP.h:527
lst Internal
Definition: HEP.h:524
string SaveDir
Definition: HEP.h:537
int IBPmethod
Definition: HEP.h:523
lst External
Definition: HEP.h:525
bool ap_rules
Definition: HEP.h:522
void init_smap()
Definition: HEP.h:539