HepLib
Loading...
Searching...
No Matches
UKIRA.cpp
Go to the documentation of this file.
1
6#include "IBP.h"
7#include <cmath>
8
9namespace HepLib {
10
14 string UKIRA::Fout(const ex & expr) {
15 if(!using_uw) {
16 ex f = expr;
17 if(is_a<lst>(f)) f = F(f);
18 else if(expr.match(F(w1,w2))) f = F(f.op(1));
19 string fstr = ex2str(f);
20 string_replace_all(fstr,"{","");
21 string_replace_all(fstr,"}","");
22 string_replace_all(fstr,"(","[");
23 string_replace_all(fstr,")","]");
24 return fstr;
25 } else {
26 ex idx;
27 if(expr.match(F(w1,w2))) idx = expr.op(1);
28 else if(expr.match(F(w))) idx = expr.op(0);
29 else idx = expr;
30 return to_string(i2w[idx]);
31 }
32 }
33
37 ex UKIRA::Fin(const string & expr) {
38 if(!using_uw) {
39 string fstr = expr;
40 string_replace_all(fstr,"[","("+to_string(ProblemNumber)+",{");
41 string_replace_all(fstr,"]","})");
42 return str2ex(fstr);
43 } else {
44 auto cpos = expr.find("*");
45 if(cpos==string::npos) {
46 if(expr=="0") return 0;
47 throw Error("UKIRA::Fin with 0 or * NOT Found.");
48 }
49 auto wstr = expr.substr(0,cpos);
50 unsigned long long weight = stoull(wstr,NULL,0);
51 auto oex = str2ex(expr.substr(cpos+1,string::npos));
52 return F(ProblemNumber, w2i[weight]) * oex;
53 }
54 }
55
59 void UKIRA::Export() {
60
61 if(Integral.nops()<1) return;
62
63 int pdim = Propagator.nops();
64 lst InExternal;
65 for(auto ii : Internal) InExternal.append(ii);
66 for(auto ii : External) InExternal.append(ii);
67
68 if(ISP.nops()<1) {
69 for(auto it : Internal) {
70 for(auto ii : InExternal) ISP.append(it*ii);
71 }
72 ISP.sort();
73 ISP.unique();
74 }
75
76 if(ISP.nops() > pdim) {
77 cout << "ISP = " << ISP << endl;
78 cout << "Propagator = " << Propagator << endl;
79 throw Error("UKIRA::Export: #(ISP) > #(Propagator).");
80 }
81
82 lst sp2s, s2sp, ss;
83 int _pic=0;
84 for(auto item : ISP) {
85 _pic++;
86 Symbol si("P"+to_string(_pic));
87 ss.append(si);
88 sp2s.append(w*item==w*si);
89 sp2s.append(item==si);
90 s2sp.append(si==item);
91 }
92
93 if(!has_w(Replacement)) {
94 lst repl;
95 for(auto item : Replacement) {
96 repl.append(item);
97 repl.append(w*item.op(0) == w*item.op(1));
98 }
99 Replacement = repl;
100 }
101
102 lst leqns;
103 for(int i=0; i<ISP.nops(); i++) { // note NOT pdim
104 auto eq = Propagator.op(i).expand().subs(iEpsilon==0); // drop iEpsilon
105 eq = eq.subs(sp2s);
106 eq = eq.subs(Replacement);
107 if(eq.has(iWF(w))) throw Error("UKIRA::Export, iWF used in eq.");
108 leqns.append(eq == iWF(i));
109 }
110 auto s2p = lsolve(leqns, ss);
111 if(s2p.nops() != ISP.nops()) throw Error("UKIRA::Export: lsolve failed.");
112
113 if(DSP.nops()<1) {
114 for(auto p1 : Internal)
115 for(auto p2 : InExternal)
116 DSP.append(lst{p1,p2});
117 }
118
119 ibps.remove_all(); // no need
120 lst nsa;
121 for(int i=0; i<pdim; i++) nsa.append(a(i));
122 for(auto sp : DSP) {
123 auto ilp = sp.op(0);
124 auto iep = sp.op(1);
125
126 ex ibp = 0;
127 symbol ss;
128 for(int i=0; i<pdim; i++) {
129 auto ns = nsa;
130 ns.let_op(i) = nsa.op(i) + 1;
131 auto dp = Propagator.op(i).subs(ilp==ss).diff(ss).subs(ss==ilp);
132 ibp -= (a(i)+Shift[i+1]) * F(ns) * dp;
133 }
134
135 ibp = ibp * iep;
136 ibp = ibp.expand();
137 ibp = ibp.subs(sp2s);
138 ibp = ibp.subs(Replacement);
139 ibp = ibp.subs(s2p);
140
141 ex res = 0;
142 for(int i=0; i<pdim; i++) {
143 auto ci = ibp.coeff(iWF(i), 1);
144 ci = MapFunction([i](const ex &e, MapFunction &self)->ex {
145 if(!e.has(F(w))) return e;
146 else if(e.match(F(w))) {
147 lst tmp = ex_to<lst>(e.op(0));
148 tmp.let_op(i) = tmp.op(i)-1;
149 return F(tmp);
150 } else return e.map(self);
151 })(ci);
152 res += ci;
153 }
154 res += ibp.subs(lst{iWF(w)==0});
155 auto cilp = iep.coeff(ilp);
156 if(!is_zero(cilp)) res += d*cilp*F(nsa);
157 ibps.append(res);
158 }
159
160 // seeds generation
161 int rmax = -1, smax = -1;
162 int rrmax[pdim], ssmax[pdim];
163 for(int i=0; i<pdim; i++) rrmax[i] = ssmax[i] = -1;
164 for(auto integral : Integral) {
165 if(integral.nops()!=pdim) throw Error("UKIRA::Export, integral dimension not match propagators.");
166 int rr = 0;
167 int ss = 0;
168 for(int i=0; i<pdim; i++) {
169 auto item = ex2int(integral.op(i));
170 if(item>0) {
171 if(rrmax[i]<item) rrmax[i] = item;
172 rr += item;
173 } else {
174 item = 0-item;
175 if(ssmax[i]<item) ssmax[i] = item;
176 ss += item;
177 }
178 }
179 if(rmax<rr) rmax = rr;
180 if(smax<ss) smax = ss;
181 }
182
183 if(seed_option == 0) {
184 int _rrmax = -1, _ssmax = -1;
185 for(int i=0; i<pdim; i++) {
186 if(_rrmax<rrmax[i]) _rrmax = rrmax[i];
187 if(_ssmax<ssmax[i]) _ssmax = ssmax[i];
188 }
189 for(int i=0; i<pdim; i++) {
190 rrmax[i] = _rrmax + rap;
191 ssmax[i] = _ssmax + sap;
192 }
193 rmax += ra;
194 smax += sa;
195 } else {
196 for(int i=0; i<pdim; i++) {
197 rrmax[i] += rap;
198 ssmax[i] += sap;
199 }
200 rmax += ra;
201 smax += sa;
202 }
203
204 // generate IBP equations
205 int as[pdim];
206 for(int i=0; i<pdim; i++) as[i] = -ssmax[i];
207 vector<vector<int>> asvec;
208 while(true) {
209 for(int i=0; i<pdim; i++) {
210 if(as[i]+1>rrmax[i]) {
211 if(i+1 == pdim) goto done;
212 as[i] = -ssmax[i];
213 } else {
214 as[i] = as[i] + 1;
215 break;
216 }
217 }
218 int _rmax = 0, _smax = 0;
219 for(int i=0; i<pdim; i++) {
220 if(as[i]>0) _rmax += as[i];
221 else _smax -= as[i];
222 }
223 if(_rmax>rmax || _smax>smax) continue;
224 vector<int> asv;
225 for(int i=0; i<pdim; i++) asv.push_back(as[i]);
226 asvec.push_back(asv);
227 }
228 done: ;
229
230 int nCut = Cut.nops();
231 bool hasCut = (nCut>1);
232 int iCut[nCut+1];
233 for(int i=0; i<nCut; i++) iCut[i] = ex_to<numeric>(Cut.op(i)).to_int();
234 auto eqns_result =
235 GiNaC_Parallel(asvec.size(), [this,&asvec,pdim,hasCut,nCut,&iCut](int idx)->ex {
236 auto as = asvec[idx];
237 exmap sol;
238 for(int i=0; i<pdim; i++) sol[a(i)]=as[i];
239 lst eqns;
240 for(auto item : ibps) {
241 auto ii = item.subs(sol);
242 if(ii.is_zero()) continue;
243 exset fs;
244 find(ii, F(w), fs);
245 if(hasCut) {
246 exmap repl;
247 for(auto fi : fs) {
248 lst ns = ex_to<lst>(fi.op(0));
249 for(int nc=0; nc<nCut; nc++) {
250 int j = iCut[nc]-1;
251 if(ns.op(j)<=0) {
252 repl[fi]=0;
253 break;
254 }
255 }
256 }
257 ii = ii.subs(repl);
258 }
259 if(ii.is_zero()) continue;
260
261 exmap repl;
262 for(auto li : Internal) {
263 for(auto fi : fs) {
264 lst ns = ex_to<lst>(fi.op(0));
265 bool has = false;
266 for(int i=0; i<pdim; i++) {
267 if(is_zero(Shift[i+1]) && ns.op(i)<=0) continue;
268 if(Propagator.op(i).has(li)) {
269 has = true;
270 break;
271 }
272 }
273 if(!has) repl[fi]=0;
274 }
275 }
276 ii = ii.subs(repl);
277 if(ii.is_zero()) continue;
278
279 eqns.append(ii);
280 }
281 return eqns;
282 }, "SEED");
283
284 lst eqns;
285 for(auto ilst : eqns_result) {
286 if(ilst.nops()<1) continue;
287 for(auto eqn : ilst) eqns.append(eqn);
288 }
289
290 if(using_uw) {
291 exset fs;
292 for(auto eqn : eqns) find(eqn, F(w), fs);
293 for(auto intg : Integral) fs.insert(F(intg));
294 for(auto intg : PIntegral) fs.insert(F(intg));
295 exvector intg_vec;
296 for(auto fi : fs) {
297 lst rs,ss;
298 int sid=0, rsum=0, ssum=0;
299 auto idx_lst = fi.op(0);
300 for(int i=0; i<idx_lst.nops(); i++) {
301 auto idx = ex_to<numeric>(idx_lst.op(i)).to_int();
302 if(idx!=0) sid += std::pow(2,idx_lst.nops()-i-1);
303 if(idx>0) {
304 rs.append(idx);
305 rsum += idx;
306 } else {
307 ss.append(idx);
308 ssum -= idx;
309 }
310 }
311
312 lst item;
313 if(sort_option==0) { // {r+s,r,s}
314 item.append(rsum+ssum);
315 item.append(rsum);
316 item.append(ssum);
317 for(auto ii : ss) item.append(-ii);
318 for(auto ii : rs) item.append(ii);
319 } else if(sort_option==1) { // {S,r,s,ss,rr}
320 item.append(rsum+ssum);
321 item.append(rsum);
322 item.append(ssum);
323 item.append(sid);
324 for(auto ii : ss) item.append(ii);
325 for(auto ii : rs) item.append(ii);
326 } else if(sort_option==2) { // {S,s,r,rr,ss}
327 item.append(rsum+ssum);
328 item.append(ssum);
329 item.append(rsum);
330 item.append(sid);
331 for(auto ii : rs) item.append(ii);
332 for(auto ii : ss) item.append(ii);
333 } else if(sort_option==-1) { // {S,r,s,-ss,rr}
334 item.append(rsum+ssum);
335 item.append(rsum);
336 item.append(ssum);
337 item.append(sid);
338 for(auto ii : ss) item.append(-ii);
339 for(auto ii : rs) item.append(ii);
340 } else if(sort_option==-2) { // {S,s,r,rr,-ss}
341 item.append(rsum+ssum);
342 item.append(ssum);
343 item.append(rsum);
344 item.append(sid);
345 for(auto ii : rs) item.append(ii);
346 for(auto ii : ss) item.append(-ii);
347 }
348 item.append(fi);
349 intg_vec.push_back(item);
350 }
351
352 sort_vec(intg_vec);
353 unsigned long long int64 = 100000000000000;
354 for(auto intg : intg_vec) {
355 int64++;
356 unsigned long long weight = int64;
357 auto idx = intg.op(intg.nops()-1).op(0);
358 for(int nc=0; nc<nCut; nc++) {
359 int j = iCut[nc]-1;
360 if(idx.op(j)>1) {
361 weight += 100000000000000;
362 break;
363 }
364 }
365 i2w[idx] = weight;
366 w2i[weight] = idx;
367 }
368 }
369
370 if(WorkingDir.length()<1) WorkingDir = to_string(getpid())+"IBP";
371 string job_dir = WorkingDir + "/" + to_string(ProblemNumber);
372 system(("rm -rf "+job_dir).c_str());
373 if(!dir_exists(job_dir)) system(("mkdir -p "+job_dir).c_str());
374
375 ostringstream oss;
376 for(auto eqn : eqns) {
377 if(is_zero(eqn)) continue;
378 auto cvs = collect_lst(eqn,F(w));
379 for(auto cv : cvs) {
380 oss << Fout(cv.op(1)) << " * (" << cv.op(0) << ")" << endl;
381 }
382 oss << endl;
383 }
384
385 ofstream eqn_out(job_dir+"/equations");
386 string ostr = oss.str();
387 string_replace_all(ostr, "{", "");
388 string_replace_all(ostr, "}", "");
389 eqn_out << ostr << endl;
390 eqn_out.close();
391
392 oss.str("");
393 oss.clear();
394 oss << "jobs:" << endl;
395 oss << " - reduce_user_defined_system:" << endl;
396 oss << " input_system: " << endl;
397 oss << " config: false" << endl;
398 oss << " files: [equations]" << endl;
399 if(PIntegral.nops()>0) {
400 ostringstream oss2;
401 int nn = PIntegral.nops();
402 for(int i=0; i<nn; i++) oss2 << Fout(PIntegral.op(i)) << endl;
403 ofstream pref_out(job_dir+"/preferred");
404 pref_out << oss2.str() << endl;
405 pref_out.close();
406 oss << " preferred_masters: preferred" << endl;
407 }
408 oss << " - kira2file:" << endl;
409 oss << " target:" << endl;
410 if(!using_uw) oss << " - [F,integrals]" << endl;
411 else oss << " - [Tuserweight,integrals]" << endl;
412 ofstream job_out(job_dir+"/jobs");
413 job_out << oss.str() << endl;
414 job_out.close();
415
416 oss.str("");
417 oss.clear();
418 for(auto integral : Integral) oss << Fout(integral) << endl;
419 ofstream intg_out(job_dir+"/integrals");
420 intg_out << oss.str() << endl;
421 intg_out.close();
422 }
423
427 void UKIRA::Run() {
428 string job_dir = WorkingDir + "/" + to_string(ProblemNumber);
429 ostringstream cmd;
430 cmd << "cd " << job_dir << " && kira " << KArgs << " --silent jobs >/dev/null 2>&1";
431 system(cmd.str().c_str());
432 }
433
437 void UKIRA::Import() {
438 string job_dir = WorkingDir + "/" + to_string(ProblemNumber);
439 ostringstream fn;
440 if(!using_uw) fn << job_dir << "/results/F/kira_integrals.kira";
441 else fn << job_dir << "/results/Tuserweight/kira_integrals.kira";
442 auto strvec = file2strvec(fn.str());
443
444 ex exL=0, exR=0;
445 map<ex,int,ex_is_less> flags;
446 lst exRs;
447 for(auto intg : Integral) flags[F(ProblemNumber,intg)] = 1;
448 Rules.remove_all();
449 for(auto line : strvec) {
450 if(line.size()==0) {
451 if(!is_zero(exL)) {
452 Rules.append(exL==exR);
453 flags[exL] = 0;
454 exRs.append(exR);
455 }
456 exL = exR = 0;
457 } else if(is_zero(exL)) {
458 exL -= Fin(line);
459 if(!exL.match(F(w1,w2))) {
460 cout << line << endl;
461 throw Error("UKIRA::Import error found.");
462 }
463 } else {
464 exR += Fin(line);
465 }
466 }
467 if(!is_zero(exL)) {
468 Rules.append(exL==exR);
469 flags[exL] = 0;
470 exRs.append(exR);
471 }
472 MIntegral.remove_all();
473 for(auto kv : flags) {
474 if(kv.second!=0) MIntegral.append(kv.first);
475 }
476 exset miset;
477 find(exRs,F(w1,w2),miset);
478 for(auto mi : miset) MIntegral.append(mi);
479 MIntegral.sort();
480 MIntegral.unique();
481
482 }
483
484}
a
int * a
Definition Functions.cpp:234
IBP.h
IBP header file.
HepLib::Error
class used to wrap error message
Definition BASIC.h:242
HepLib::IBP::Internal
lst Internal
Definition IBP.h:27
HepLib::IBP::ISP
lst ISP
Definition IBP.h:34
HepLib::IBP::Replacement
lst Replacement
Definition IBP.h:29
HepLib::IBP::DSP
lst DSP
Definition IBP.h:33
HepLib::IBP::Integral
lst Integral
Definition IBP.h:31
HepLib::IBP::External
lst External
Definition IBP.h:28
HepLib::IBP::Propagator
lst Propagator
Definition IBP.h:30
HepLib::IBP::ProblemNumber
int ProblemNumber
Definition IBP.h:39
HepLib::Symbol
class extended to GiNaC symbol class, represent a positive symbol
Definition BASIC.h:113
HepLib::UKIRA::using_uw
bool using_uw
Definition IBP.h:119
HepLib::UKIRA::Export
void Export() override
Export input data for KIRA.
Definition UKIRA.cpp:59
HepLib
HepLib namespace.
Definition BASIC.cpp:17
HepLib::iEpsilon
const iSymbol iEpsilon
HepLib::sp
ex sp(const ex &a, const ex &b)
translated the vector dot a.b to a*b, useful in SecDec
Definition Pair.cpp:237
HepLib::has_w
bool has_w(const ex &e)
Definition BASIC.cpp:2233
HepLib::str2ex
ex str2ex(const string &expr, symtab stab)
convert string to ex expression, using Parser internally
Definition BASIC.cpp:670
HepLib::w
ex w
Definition Init.cpp:93
HepLib::d
const Symbol d
HepLib::file2strvec
vector< string > file2strvec(string filename, bool skip_empty)
read file content to string vector
Definition BASIC.cpp:809
HepLib::as
const Symbol as
HepLib::collect_lst
lst collect_lst(ex const &expr_in, std::function< bool(const ex &)> has_func, int opt)
the collect function like Mathematica, reture the lst { {c1,v1}, {c2,v2}, ... }
Definition BASIC.cpp:1222
HepLib::sort_vec
void sort_vec(exvector &ivec, bool less=true)
sort the list in less order, or the reverse
Definition Sort.cpp:54
HepLib::string_replace_all
void string_replace_all(string &str, const string &from, const string &to)
Definition Functions.cpp:148
HepLib::ex2str
string ex2str(const ex &expr)
convert ex to output string, the defalut printer format will be used
Definition BASIC.cpp:715
HepLib::dir_exists
bool dir_exists(string dir)
Definition BASIC.h:297
HepLib::w1
ex w1
Definition BASIC.h:499
HepLib::GiNaC_Parallel
exvector GiNaC_Parallel(int ntotal, std::function< ex(int)> f, const string &key)
GiNaC Parallel Evaluation using fork.
Definition BASIC.cpp:260
HepLib::w2
ex w2
Definition BASIC.h:499
HepLib::ex2int
int ex2int(ex num)
ex to integer
Definition BASIC.cpp:893
Generated by doxygen 1.9.8