HepLib
SD.h
Go to the documentation of this file.
1 
6 #pragma once
7 
8 extern "C" {
9  #include <quadmath.h>
10 }
11 
12 #include <dlfcn.h>
13 #include <string>
14 #include <signal.h>
15 #include <sys/syscall.h>
16 #include <sys/wait.h>
17 #include <sstream>
18 #include <ios>
19 #include <regex>
20 #include <complex>
21 #include "mpreal.h"
22 #include "BASIC.h"
23 
27 namespace HepLib::SD {
28 
29  using namespace HepLib;
30 
31  /*-----------------------------------------------------*/
32  // Global Functions
33  /*-----------------------------------------------------*/
34  exvector get_xy_from(ex pol);
35  exvector get_x_from(ex pol);
36  exvector get_y_from(ex pol);
37  exvector get_z_from(ex pol);
38  exvector get_pl_from(ex pol);
39  int epsRank(ex,ex);
40  int vsRank(ex);
41  int x_free_index(ex expr);
42  int y_free_index(ex expr);
43  ex Factor(const ex expr);
44  ex FactorOutX(const ex expr);
45  ex exp_simplify(const ex);
46  ex pow_simplify(const ex);
47  ex xyz_pow_simplify(const ex expr);
48 
49  /*-----------------------------------------------------*/
50  // Customized GiNaC Function
51  /*-----------------------------------------------------*/
52  DECLARE_FUNCTION_1P(fabs)
53  DECLARE_FUNCTION_1P(PL)
54  DECLARE_FUNCTION_1P(CT)
59  DECLARE_FUNCTION_2P(CV) // not used internally, for user use only
60  DECLARE_FUNCTION_1P(WRA) // for wick rotation
61  extern int VEO_Digits;
62 
70  lst Exponent;
71  exmap lReplacement;
72  exmap tReplacement;
73  exmap nReplacement;
74  ex Prefactor = 1;
75  bool isQuasi = false;
76  bool isAsy = false;
77  };
78 
82  struct XIntegrand {
83  lst Function;
84  lst Exponent;
85  exmap nReplacement;
86  lst Deltas;
87  bool isAsy = false;
88  };
89 
93  class SecDecBase {
94  public:
95  virtual ~SecDecBase() { }
96  virtual vector<exmap> x2y(const ex &xpol) =0;
97  vector<exmap> x2y(const lst &xpols);
99  static bool VerifySD(vector<exmap> vmap, bool quick = true);
100  static ex XMonomials(const ex & expr);
101  };
102 
106  class SecDecG : public SecDecBase {
107  public:
108  vector<exmap> x2y(const ex &xpol) override;
109  static vector<vector<int>> RunQHull(const matrix &pts);
110  private:
111  vector<matrix> ZeroFaces(const matrix &pts);
112  matrix NormalVectors(const vector<matrix> &zfs);
113  matrix DualCone(const matrix &pts);
114  vector<matrix> SimplexCones(matrix pts);
115  };
116 
117  /*-----------------------------------------------------*/
118  // Integrator Classes
119  /*-----------------------------------------------------*/
120 
121  typedef long double dREAL;
122  typedef complex<dREAL> dCOMPLEX;
123  typedef __float128 qREAL;
124  typedef __complex128 qCOMPLEX;
125  typedef mpfr::mpreal mpREAL;
126  typedef complex<mpREAL> mpCOMPLEX;
127 
132  public:
133  virtual ~IntegratorBase() { }
134 
135  typedef int (*SDD_Type) (const unsigned int xn, const dREAL x[], const unsigned int yn, dREAL y[], const dREAL pl[], const dREAL las[]);
136  typedef int (*SDQ_Type) (const unsigned int xn, const qREAL x[], const unsigned int yn, qREAL y[], const qREAL pl[], const qREAL las[]);
137  typedef int (*SDMP_Type) (const unsigned int xn, const mpREAL x[], const unsigned int yn, mpREAL y[], const mpREAL pl[], const mpREAL las[]);
138  typedef qREAL (*FT_Type) (const qREAL xx[], const qREAL pl[]);
139  virtual ex Integrate(size_t n=0) =0;
140 
141  FT_Type FT = NULL;
142  SDD_Type IntegrandD = NULL;
143  SDQ_Type IntegrandQ = NULL;
144  SDMP_Type IntegrandMP = NULL;
145  const dREAL* dLambda;
147  const qREAL* qLambda;
149  const mpREAL* mpLambda;
151  int XDim;
152 
153  qREAL EpsAbs = 1E-5;
154  qREAL EpsRel = 0;
155  int ReIm = 3; // 1-Re, 2-Im, 3-ReIm
156  int MPDigits = 64;
157  size_t NEval = 0;
158  protected:
159  time_t StartTimer; // used internally
160  size_t RunTime = 0;
161  };
162 
166  class HCubature : public IntegratorBase {
167  public:
168  static int Wrapper(unsigned int xdim, size_t npts, const qREAL *x, void *fdata, unsigned int ydim, qREAL *y);
169 
170  typedef void (*PrintHookerType) (qREAL*, qREAL*, size_t *, void *);
171 
172  virtual ex Integrate(size_t n=0) override;
173  static void DefaultPrintHooker(qREAL*, qREAL*, size_t *, void*);
174  PrintHookerType PrintHooker = DefaultPrintHooker;
175  bool use_last = false;
176 
177  int DQMP = 0;
178  int QXDim = 0;
179  int MPXDim = 0;
180  qREAL QXLimit = 1E-6Q;
181  qREAL MPXLimit = 1E-8Q;
182  qREAL QFLimit = -1;
183  qREAL MPFLimit = -1;
184  size_t RunMAX = 100;
185  size_t RunPTS = 100000;
186  size_t MinPTS = 0;
187  size_t MaxPTS; // used internally
188  unsigned int Threads = 0;
189  private:
190  int NANMax = 250;
191  int nNAN = 0;
192  size_t lastNRUN = 0;
193  qREAL LastResult[2];
194  qREAL LastAbsErr[2];
195  int lastnNAN = 0;
196  int LastState = 0;
197 
198  int inDQMP(qREAL const *x);
199  };
200 
204  class HCubatureMP : public IntegratorBase {
205  public:
206  static int Wrapper(unsigned int xdim, size_t npts, const mpREAL *x, void *fdata, unsigned int ydim, mpREAL *y);
207  typedef void (*PrintHookerType) (mpREAL*, mpREAL*, size_t *, void *);
208  virtual ex Integrate(size_t n=0) override;
209  static void DefaultPrintHooker(mpREAL*, mpREAL*, size_t *, void*);
210  PrintHookerType PrintHooker = DefaultPrintHooker;
211  size_t RunMAX = 100;
212  size_t RunPTS = 100000;
213  size_t MinPTS = 0;
214  size_t MaxPTS; // used internally
215  unsigned int Threads = 0;
216  private:
217  ex mp2ex(const mpREAL & num);
218  int NANMax = 250;
219  int nNAN = 0;
220  bool use_last = false;
221  size_t lastNRUN = 0;
222  mpREAL LastResult[2];
223  mpREAL LastAbsErr[2];
224  int lastnNAN = 0;
225  int LastState = 0;
226 
227  int inDQMP(qREAL const *x);
228  };
229 
233  class TanhSinhMP : public IntegratorBase {
234  public:
235  static int Wrapper(unsigned ydim, mpREAL *y, mpREAL *e, unsigned xdim, const mpREAL *x, void *fdata);
236  typedef void (*PrintHookerType) (mpREAL *, mpREAL *, size_t *, void *);
237  virtual ex Integrate(size_t n=0) override;
238  static void DefaultPrintHooker(mpREAL *, mpREAL *, size_t *, void *);
239  PrintHookerType PrintHooker = DefaultPrintHooker;
240  TanhSinhMP(size_t k=10);
241  private:
242  ex mp2ex(const mpREAL & num);
243  size_t K = 10;
244  };
245 
249  class QuadMP : public IntegratorBase {
250  public:
251  static int Wrapper(unsigned yn, mpREAL *y, mpREAL *e, unsigned xdim, const mpREAL *x, void *fdata);
252  typedef void (*PrintHookerType) (mpREAL*, mpREAL*, size_t *, void *);
253  virtual ex Integrate(size_t n=0) override;
254  static void DefaultPrintHooker(mpREAL *, mpREAL *, size_t *, void *);
255  PrintHookerType PrintHooker = DefaultPrintHooker;
256  QuadMP() { }
257  QuadMP(size_t m) : mGK(m) { }
258  size_t nGK = 100;
259  size_t mGK = 10;
260  private:
261  ex mp2ex(const mpREAL & num);
262  };
263 
264  typedef long double dREAL;
268  class MinimizeBase {
269  public:
270  typedef dREAL (*FunctionType)(int nvars, dREAL* x, dREAL* pl, dREAL *las);
271  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;
272  dREAL ZeroValue = -1E-20;
273  virtual void Minimize(int nvars, FunctionType func, dREAL *ip)=0;
274  virtual void ForceStop()=0;
275  };
276 
280  class HookeJeeves : public MinimizeBase {
281  public:
282  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) override;
283  bool Exit = false;
284  virtual void Minimize(int nvars, FunctionType func, dREAL *ip) override;
285  virtual void ForceStop() override;
286 
287  private:
288  dREAL best_nearby(dREAL* delta, dREAL* point, dREAL prevbest, int nvars);
289  int hooke(int nvars, dREAL* startpt, dREAL* endpt, dREAL rho, dREAL epsilon, int itermax);
290  dREAL ObjectWrapper(int nvars, dREAL* x);
291  FunctionType ObjectFunction;
292  dREAL UpperBound[50];
293  dREAL LowerBound[50];
294  dREAL *PL;
295  dREAL *LAS;
296  };
297 
301  class CppFormat : public print_csrc_cl_N {
302  GINAC_DECLARE_PRINT_CONTEXT(CppFormat, print_csrc_cl_N)
303  public:
304  CppFormat(ostream &os, const string & s = "L", unsigned opt = 0);
305  string suffix;
306  string MQuote = "\"";
307 
308  template<class T> const CppFormat & operator << (const T & v) const {
309  s << v;
310  return *this;
311  };
312  const CppFormat & operator << (const basic & v) const;
313  const CppFormat & operator << (const ex & v) const;
314  const CppFormat & operator << (const lst & v) const;
315  const CppFormat & operator<<(std::ostream& (*v)(std::ostream&)) const;
316 
317  #ifndef DOXYGEN_SKIP
318  class _init {
319  public: _init();
320  };
321  #endif
322  private:
323  #ifndef DOXYGEN_SKIP
324  static _init CppFormat_init;
325  #endif
326  static void print_integer(const CppFormat & c, const cln::cl_I & x);
327  static void print_real(const CppFormat & c, const cln::cl_R & x);
328  static void print_numeric(const numeric & p, const CppFormat & c, unsigned level);
329  };
330 
331  class ExFormat : public print_dflt {
332  GINAC_DECLARE_PRINT_CONTEXT(ExFormat, print_dflt)
333  public:
334  ExFormat(ostream &os, const string & s = "L", unsigned opt = 0);
335  string suffix;
336  string type = "ex";
337  string MQuote = "\"";
338 
339  template<class T> const ExFormat & operator << (const T & v) const {
340  s << v;
341  return *this;
342  };
343  const ExFormat & operator << (const basic & v) const;
344  const ExFormat & operator << (const ex & v) const;
345  const ExFormat & operator << (const lst & v) const;
346  const ExFormat & operator<<(std::ostream& (*v)(std::ostream&)) const;
347 
348  class _init {
349  public: _init();
350  };
351  private:
352  static _init ExFormat_init;
353  static void print_integer(const ExFormat & c, const cln::cl_I & x);
354  static void print_real(const ExFormat & c, const cln::cl_R & x);
355  static void print_numeric(const numeric & p, const ExFormat & c, unsigned level);
356  };
357 
358  /*-----------------------------------------------------*/
359  // VE
360  /*-----------------------------------------------------*/
361  ex VESimplify(ex expr);
362  ex VEResult(ex expr);
363  ex VEResult2(ex expr); // keep two digits in error
364  ex VEMaxErr(ex expr);
365 
369  class ErrMin {
370  public:
372  static qREAL *paras;
373  static dREAL err_max;
374  static dREAL err_min;
375  static size_t MaxRND;
376  static size_t RunRND;
378  static dREAL *lambda;
379  static dREAL hjRHO;
380  static ex lastResErr;
381  static dREAL IntError(int nvars, dREAL *las, dREAL *n1, dREAL *n2);
382  };
383 
387  class ChengWu {
388  public:
389  static bool isProjective(const ex fe, const ex delta);
390  static void Projectivize(ex &fe, const ex delta, const ex xsum=0);
391  static void Scalelize(ex &fe, const lst xs, const ex cy);
392  static void Scalelize(ex &fe, const ex xi, const ex cy);
393  static exvector Binarize(ex const fe, ex const eqn);
394  static void Binarize(ex const fe, ex const eqn, exvector & ovec);
395  static bool isLinearizable(const ex ft, const ex delta, lst & xcs);
396  static void Linearize(const lst xcs, ex & fe, ex & ft);
397  static bool isPartilizable(const ex ft, const ex delta, lst &xcs, int mode=0);
398  static void Partilize(const lst xcs, const lst delta, const ex fe, exvector & ret_lst);
399 
400  static exvector Evaluate(const ex & fe);
401  static exvector WickRotation(const exvector & fe_vec);
402  static exvector Apply(const exvector & fe_vec, const ex & ft=0);
403  inline static exvector Apply(const ex & fe, const ex & ft=0) {
404  exvector fe_vec;
405  fe_vec.push_back(fe);
406  return Apply(fe_vec, ft);
407  }
408  };
409 
413  class SecDec {
414 
415  public:
416  static bool use_dlclose;
417  static string cpp;
418 
419  lst eps_lst = { lst{eps,0}, lst{ep,0} }; // { {epi, epiN}, ... }
420  int vsN = 0;
421  int PoleRequested = -5;
422  bool vs_before_ep = false;
423  bool use_XMonomials = true;
424  bool disable_Contour = false;
426  exvector FunExp; // each item : { {f1,f2,...}, {n1,n2,...}, { delta_list1, delta_list2 } }
427  exvector Integrands;
428  exvector expResult;
430  IntegratorBase *Integrator = NULL;
431  MinimizeBase *Minimizer = NULL;
433  bool IsZero = false;
434  bool CheckEnd = false;
435  bool use_ErrMin = false;
436  bool use_las = false;
437  bool save_las = false;
438  int IBF = 0; // 0 - not use IBF
439  bool use_Normalizes = true;
440  bool use_XReOrders = false;
441  int MPDigits = 0; // digits in mpREAL for MP
442  lst BisectionPoints = lst { ex(1)/13, ex(1)/19, ex(1)/29, ex(1)/59, ex(1)/41, ex(1)/37, ex(1)/43, ex(1)/53 };
443 
444  map<int, numeric> Parameter; // used Contours and Integrates, use PL in Prepares part
445 
446  // used in Contours
447  bool CTMaxF = true;
448  dREAL CTLaMax = 10; // CTLaMax<0 for explict REAL mode
449  int CTTryPTS = 3;
450  int CTSavePTS = 3;
451 
452  size_t LambdaSplit = 5;
453  qREAL IntLaMax = 50;
454  int CTryM = 1; // try lambda in Middle
455  int CTryL = 1; // try lambda in Left
456  int CTryR = 1; // try lambda in Right
457  size_t CTryI = 0; // integrator limit in CTry
458  dREAL CTryRRatio = 1.5;
459  int soLimit = 10000;
460 
461  qREAL EpsAbs = 1E-4;
462  int ReIm = 3; // 1-Re, 2-Im, 3-ReIm
463 
464  void Initialize(FeynmanParameter fpi);
465  void Initialize(XIntegrand xint);
466  void BiSection(ex xi, ex x0);
467  void Normalizes();
468  void Scalelesses();
469  void SDPrepares();
470  void EpsExpands();
471  void RemoveDeltas();
472  void XReOrders();
473  void XTogethers();
474  void XExpands();
475  void KillPowers(int bits=1+2);
476  bool IsBad(ex f, vector<exmap> vmap);
477  exvector AutoEnd(ex po_ex);
478  void CIPrepares(const string & key = "");
479  void Contours(const string & key = "", const string & pkey = "");
480  void Integrates(const string & key="", const string & pkey="", int kid=0);
481  void ReIntegrates(const string & key, const string & pkey, qREAL err);
482  void Evaluate(FeynmanParameter fpi, const string & key = "");
483  void Evaluate(XIntegrand xint, const string & key = "");
484  void Evaluate(const exvector & FunExp, const string & key = "");
485  void MB();
486  void XEnd();
487  void ChengWu(const ex & ft=0);
488 
489  static bool VerifySD(vector<exmap> vmap, bool quick = true);
490  static ex XRefined(ex const & ft);
491  static lst XRefined_lst(ex const & ft);
492  static ex PrefactorFIESTA(int nLoop);
493  ex VEResult();
494  void VEPrint(bool endlQ=true);
495  static ex PExpand(ex xpol, bool delta=true);
496  static int PRank(matrix m);
497  static ex ContinuousWRA(ex expr_in, int nc=15);
498 
499  ~SecDec();
500 
501  private:
502  exvector DS(const ex po_ex);
503  lst Normalize(const ex &input);
504  void DoAsy();
505  bool KillPowerD(ex fe, int kpi);
506  bool KillPower(ex fe, int kpi, int bits);
507 
508  void CompileMatDet();
509  vector<lst> ciResult;
510  lst FT_N_XN; // list of { ft, n, xn }
511  exmap LambdaMap;
512 
513  };
514 
518  class cseParser {
519  public:
520  ex Parse(ex expr, bool reset=true);
521  string oc = "o";
522  int on();
523  vector<pair<int, ex>> os();
524  private:
525  map<ex, ex, ex_is_less> ex_var_map;
526  int no = 0;
527  vector<pair<int, ex>> o_ex_vec;
528  map<int, int> used;
529  };
530 
531  class cse_Parser {
532  public:
533  ex Parse(ex expr) { return Parse(expr, true); }
534  string v = "v";
535  int vn();
536  const vector<pair<int,ex>> & vs();
537  private:
538  ex Parse(ex expr, bool reset);
539  map<ex, int, ex_is_less> exn;
540  int no = 0;
541  exvector exv;
542  map<int, int> used;
543  vector<pair<int,ex>> on_ex_vec;
544  };
545 
546 }
547 
548 
Basic header file.
class to manipulate with Cheng-Wu theorem
Definition: SD.h:387
static exvector Apply(const ex &fe, const ex &ft=0)
Definition: SD.h:403
class to export GiNaC expression to cpp format
Definition: SD.h:301
string suffix
Definition: SD.h:305
ErrMin with HookeJeeves.
Definition: SD.h:369
static ex lastResErr
Definition: SD.h:380
static IntegratorBase * Integrator
Definition: SD.h:371
static size_t RunRND
Definition: SD.h:376
static dREAL err_max
Definition: SD.h:373
static MinimizeBase * miner
Definition: SD.h:377
static qREAL * paras
Definition: SD.h:372
static dREAL * lambda
Definition: SD.h:378
static size_t MaxRND
Definition: SD.h:375
static dREAL hjRHO
Definition: SD.h:379
static dREAL err_min
Definition: SD.h:374
string suffix
Definition: SD.h:335
numerical integrator using HCubatureMP
Definition: SD.h:204
numerical integrator using HCubature
Definition: SD.h:166
size_t MaxPTS
Definition: SD.h:187
class to minimize a function using HookeJeeves
Definition: SD.h:280
base for numerical integrator
Definition: SD.h:131
const dREAL * dParameter
Definition: SD.h:146
const dREAL * dLambda
Definition: SD.h:145
const mpREAL * mpLambda
Definition: SD.h:149
const qREAL * qParameter
Definition: SD.h:148
virtual ~IntegratorBase()
Definition: SD.h:133
virtual ex Integrate(size_t n=0)=0
const mpREAL * mpParameter
Definition: SD.h:150
const qREAL * qLambda
Definition: SD.h:147
base for class to minimize a function
Definition: SD.h:268
virtual void Minimize(int nvars, FunctionType func, dREAL *ip)=0
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
virtual void ForceStop()=0
numerical integrator using TanhSinhMP
Definition: SD.h:249
QuadMP(size_t m)
Definition: SD.h:257
base class of SecDec
Definition: SD.h:93
virtual ~SecDecBase()
Definition: SD.h:95
bool use_XMonomials
Definition: SD.h:98
virtual vector< exmap > x2y(const ex &xpol)=0
SecDec by geometric method.
Definition: SD.h:106
SecDec the main class to use Sector Decompostion method.
Definition: SD.h:413
exvector expResult
Definition: SD.h:428
exmap nReplacement
Definition: SD.h:425
map< int, numeric > Parameter
Definition: SD.h:444
static bool use_dlclose
Definition: SD.h:416
ex ResultError
Definition: SD.h:432
static string cpp
Definition: SD.h:417
exvector FunExp
Definition: SD.h:426
exvector Integrands
Definition: SD.h:427
numerical integrator using TanhSinhMP
Definition: SD.h:233
class for Common SubExpression Parser
Definition: SD.h:518
ex Parse(ex expr)
Definition: SD.h:533
namespace for Numerical integration with Sector Decomposition method
Definition: AsyMB.cpp:10
ex FactorOutX(const ex expr)
Definition: Helpers.cpp:278
exvector get_z_from(ex pol)
Definition: Helpers.cpp:59
exvector get_xy_from(ex pol)
Definition: Helpers.cpp:14
ex xyz_pow_simplify(const ex expr_in)
Definition: Helpers.cpp:336
DECLARE_FUNCTION_2P(CV) DECLARE_FUNCTION_1P(WRA) extern int VEO_Digits
ex exp_simplify(const ex expr_in)
Definition: Helpers.cpp:302
int x_free_index(ex expr)
Definition: Helpers.cpp:165
complex< dREAL > dCOMPLEX
Definition: SD.h:122
int vsRank(ex expr_in)
Definition: Helpers.cpp:211
__float128 qREAL
Definition: SD.h:123
ex VEResult(ex expr)
Definition: Helpers.cpp:123
long double dREAL
Definition: SD.h:121
ex VESimplify(ex expr)
Definition: Helpers.cpp:92
ex VEMaxErr(ex expr)
Definition: Helpers.cpp:131
complex< mpREAL > mpCOMPLEX
Definition: SD.h:126
exvector get_pl_from(ex pol)
Definition: Helpers.cpp:74
exvector get_y_from(ex pol)
Definition: Helpers.cpp:44
exvector get_x_from(ex pol)
Definition: Helpers.cpp:29
int y_free_index(ex expr)
Definition: Helpers.cpp:180
__complex128 qCOMPLEX
Definition: SD.h:124
mpfr::mpreal mpREAL
Definition: SD.h:125
int epsRank(ex expr_in, ex epi)
Definition: Helpers.cpp:195
ex Factor(const ex expr_in)
Definition: Helpers.cpp:244
ex VEResult2(ex expr)
Definition: Helpers.cpp:127
ex pow_simplify(const ex expr_in)
Definition: Helpers.cpp:317
HepLib namespace.
Definition: BASIC.cpp:17
const Symbol ep
const Symbol eps
const Symbol vs
bool IsZero(const ex &e)
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