HepLib
Contours.cpp
Go to the documentation of this file.
1 
6 #include "SD.h"
7 #include <math.h>
8 #include <cmath>
9 
10 namespace HepLib::SD {
11 
17  void SecDec::Contours(const string & key, const string & pkey) {
18  if(IsZero) return;
19  if(Minimizer==NULL) Minimizer = new HookeJeeves();
20 
21  if(key != "") {
23  ostringstream garfn;
24  garfn << key << ".ci.gar";
25  archive ar;
26  ifstream in(garfn.str());
27  in >> ar;
28  in.close();
29  auto c = ar.unarchive_ex("c");
30  if(c!=19790923) throw Error("Contours: *.ci.gar error!");
31  FT_N_XN = ex_to<lst>(ar.unarchive_ex("ftnxn"));
33  }
34 
35  //change 2->1 from GiNaC 1.7.7
36  if(FT_N_XN.nops()<1) return;
37  if(Verbose > 0) cout << Color_HighLight << " Contours @ " << now() << RESET << endl;
38 
39  exvector ftnxn_vec;
40  //change 1->0 from GiNaC 1.7.7
41  for(int i=0; i<FT_N_XN.nops(); i++) ftnxn_vec.push_back(FT_N_XN.op(i));
42 
43  auto pid = getpid();
44  ostringstream cmd;
45  cmd << "mkdir -p " << pid;
46  if(!dir_exists(to_string(pid))) auto rc = system(cmd.str().c_str());
47 
48  ostringstream fsofn;
49  if(key != "") {
50  fsofn << key << "F.so";
51  } else {
52  fsofn << pid << "F.so";
53  }
54  void* module = dlopen(fsofn.str().c_str(), RTLD_NOW);
55  if (module == nullptr) {
56  module = dlopen(("./"+fsofn.str()).c_str(), RTLD_NOW);
57  if (module == nullptr) {
58  cerr << ErrColor << "Contours: could not open compiled module!" << RESET << endl;
59  cout << "dlerror(): " << dlerror() << endl;
60  exit(1);
61  }
62  }
63 
64  GiNaC_Parallel_NP["LAS"]=0;
65  GiNaC_Parallel_RM["LAS"] = !Debug;
66  auto res =
67  GiNaC_Parallel(ftnxn_vec.size(), [&ftnxn_vec,this,module](int idx)->ex {
68  // return lst{ ft_n, lst{lambda-i, lambda-max} }
69  // with I*[lambda-i]*lambda, lambda < lambda-max
70  // note that lambda sequence only matches to x sequence in F-term
71  if(Verbose>5) cout << " λ: " << idx << " / " << ftnxn_vec.size() << flush;
72  auto ftnxn = ftnxn_vec[idx];
73  int npara = -1;
74  for(auto kv : Parameter) if(npara<kv.first) npara = kv.first;
75  dREAL paras[npara+1];
76  lst plRepl;
77  for(auto kv : Parameter) {
78  paras[kv.first] = (dREAL)ex2q(kv.second);
79  plRepl.append(PL(kv.first)==kv.second);
80  }
81 
82  auto ft = ftnxn.op(0);
83  ft = ft.subs(plRepl);
84  if(xSign(ft)>0 || CTLaMax<0) {
85  if(Verbose>5) {
86  cout << "\r \r";
87  cout << " λ: xSign>0 or CTLaMax<0, back to REAL mode!" << endl;
88  }
89  return lst{ ftnxn.op(1), 1979 }; // ft_id, las
90  }
91 
92  int nvars = ex_to<numeric>(ftnxn.op(2)).to_int();
93  ostringstream fname;
94 
95  if(CTMaxF) {
96  MinimizeBase::FunctionType fp = NULL, fp2 = NULL;
97  fname.clear();
98  fname.str("");
99  fname << "minF_"<<ftnxn.op(1);
100  fp = (MinimizeBase::FunctionType)dlsym(module, fname.str().c_str());
101  fname.clear();
102  fname.str("");
103  fname << "minFM_"<<ftnxn.op(1);
104  fp2 = (MinimizeBase::FunctionType)dlsym(module, fname.str().c_str());
105  if(fp==NULL || fp2==NULL) {
106  cerr << ErrColor << "Contours: fp==NULL/fp2==NULL" << RESET << endl;
107  cout << "dlerror(): " << dlerror() << endl;
108  exit(1);
109  }
110 
111  dREAL maxf = Minimizer->FindMinimum(nvars, fp, paras);
112  if(maxf>0) {
113  if(Verbose>5) {
114  cout << "\r \r";
115  cout << " λ: F>0 is Found, back to REAL mode!" << endl;
116  }
117  return lst{ ftnxn.op(1), 1979 }; // ft_id, las
118  }
119 
120  maxf = Minimizer->FindMinimum(nvars, fp2, paras);
121  if(maxf>0) {
122  if(Verbose>5) {
123  cout << "\r \r";
124  cout << " λ: F<0 is Found, back to REAL mode!" << endl;
125  }
126  return lst{ ftnxn.op(1), 1979 }; // ft_id, las
127  }
128  }
129 
130  dREAL nlas[nvars];
131  dREAL max_df = -1, max_f;
132  for(int i=0; i<nvars; i++) {
133  MinimizeBase::FunctionType dfp = NULL;
134  fname.clear();
135  fname.str("");
136  fname << "dirC_"<<ftnxn.op(1)<<"_"<<i;
137  dfp = (MinimizeBase::FunctionType)dlsym(module, fname.str().c_str());
138  if(dfp==NULL) {
139  cerr << ErrColor << "Contours: dfp==NULL" << RESET << endl;
140  cout << "dlerror(): " << dlerror() << endl;
141  exit(1);
142  }
143 
144  dREAL maxdf = Minimizer->FindMinimum(nvars, dfp, paras);
145  maxdf = -maxdf;
146  nlas[i] = maxdf;
147  if(max_df<maxdf) max_df = maxdf;
148  }
149 
150  for(int i=0; i<nvars; i++) { // From FIESTA
151  if(nlas[i] < 1) nlas[i] = 1;
152  else nlas[i] = 1/nlas[i];
153  }
154 
155  lst las;
156  for(int i=0; i<nvars; i++) {
157  las.append(q2ex(nlas[i]));
158  }
159 
161  fname.clear();
162  fname.str("");
163  fname << "imgF_"<<ftnxn.op(1);
164  fp = (MinimizeBase::FunctionType)dlsym(module, fname.str().c_str());
165  if(fp==NULL) {
166  cerr << ErrColor << "Contours: fp==NULL" << RESET << endl;
167  cout << "dlerror(): " << dlerror() << endl;
168  exit(1);
169  }
170 
171  dREAL laBegin = 0, laEnd = CTLaMax, min;
172  dREAL UB[nvars+1];
173  for(int i=0; i<nvars+1; i++) UB[i] = 1;
174 
175  min = laEnd;
176  while(true) {
177  UB[nvars] = min;
178  dREAL res = Minimizer->FindMinimum(nvars+1, fp, paras, nlas, UB, NULL, NULL, true, CTTryPTS, CTSavePTS);
179  if(res < -1E-5) laEnd = min;
180  else laBegin = min;
181 
182  if(laEnd < 1E-10) {
183  cout << "\r \r";
184  cout << " λ: " << ErrColor << "too small lambda!" << RESET << endl;
185  break;
186  }
187 
188  if(laEnd-laBegin <= 0.01*laEnd) break;
189  min = (laBegin + laEnd) / 2.0;
190  }
191  min = laBegin;
192 
193  las.append(q2ex(min));
194  if(Verbose>5) {
195  cout << "\r \r";
196  cout << " λ: " << NN(las,3) << endl;
197  }
198 
199  return lst{ ftnxn.op(1), las }; // ft_id, las
200 
201  }, "LAS");
202 
203  if(use_dlclose) dlclose(module);
204 
205  ostringstream garfn;
206  if(key != "") {
207  garfn << key;
208  if(pkey != "") garfn << "-" << pkey;
209  garfn << ".las.gar";
210  lst gar_res;
211  for(auto &item : res) gar_res.append(item);
212  archive ar;
213  ar.archive_ex(gar_res, "res");
214  ar.archive_ex(19790923, "c");
215  ofstream out(garfn.str());
216  out << ar;
217  out.close();
218  } else {
219  for(auto &item : res) LambdaMap[item.op(0)] = item.op(1);
220  }
221 
222  cmd.clear();
223  cmd.str("");
224  cmd << "rm -rf " << pid;
225  if(!Debug) auto rc = system(cmd.str().c_str());
226  }
227 
228 }
#define RESET
Definition: BASIC.h:79
SecDec header file.
class used to wrap error message
Definition: BASIC.h:242
class to minimize a function using HookeJeeves
Definition: SD.h:280
virtual dREAL FindMinimum(int nvars, FunctionType func, dREAL *PL=NULL, dREAL *las=NULL, dREAL *UB=NULL, dREAL *LB=NULL, dREAL *IP=NULL, bool compare0=false, int TryPTS=0, int SavePTS=0)=0
dREAL(* FunctionType)(int nvars, dREAL *x, dREAL *pl, dREAL *las)
Definition: SD.h:270
int CTSavePTS
Definition: SD.h:450
int CTTryPTS
Definition: SD.h:449
dREAL CTLaMax
Definition: SD.h:448
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
bool CTMaxF
Definition: SD.h:447
bool IsZero
Definition: SD.h:433
namespace for Numerical integration with Sector Decomposition method
Definition: AsyMB.cpp:10
long double dREAL
Definition: SD.h:121
const char * Color_HighLight
Definition: BASIC.cpp:248
void reset_precision()
Definition: BASIC.cpp:2313
ex q2ex(__float128 num)
__float128 to ex
Definition: BASIC.cpp:867
int NNDigits
Definition: Init.cpp:155
ex NN(ex expr, int digits)
the nuerical evaluation
Definition: BASIC.cpp:1278
map< string, bool > GiNaC_Parallel_RM
Definition: Init.cpp:148
const char * ErrColor
Definition: BASIC.cpp:246
string now(bool use_date)
date/time string
Definition: BASIC.cpp:525
bool Debug
Definition: Init.cpp:141
map< string, int > GiNaC_Parallel_NP
Definition: Init.cpp:144
int Verbose
Definition: Init.cpp:139
bool dir_exists(string dir)
Definition: BASIC.h:297
void set_precision(long prec, bool push)
Definition: BASIC.cpp:2304
int xSign(ex const expr)
the always sign for expr
Definition: BASIC.cpp:1312
exvector GiNaC_Parallel(int ntotal, std::function< ex(int)> f, const string &key)
GiNaC Parallel Evaluation using fork.
Definition: BASIC.cpp:260