Usage in Python

HepLib provides a wrapper class HepLibW as an interface to other languages, the SWIG interace file HepLib.i is also installed <INSTALL_PATH>/include/HepLib.i. We will just show the usage within python, but it also applys to other languages with the help of SWIG.

Compile the python module

One can use the following makefile to generate the python module files: _HepLib.so and HepLib.py.

python 3.x(>3.8)

all: _HepLib.so
HLS=HepLib_SWIG
flatns=""

uname := $(shell uname -s)
ifeq ($(uname),Darwin)
  flatns = "-flat_namespace"
endif

$(HLS).cpp:
        swig -python -c++ -o $(HLS).cpp $$(heplib-config --prefix)/include/HepLib.i
$(HLS).o : $(HLS).cpp
        heplib++ -fPIC -c $(HLS).cpp $$(python3-config --cflags)
_HepLib.so : $(HLS).o
        heplib++ -shared $(flatns) $(HLS).o -o _HepLib.so $$(python3-config --ldflags --embed)
clean:
        rm -f $(HLS)* _HepLib.so HepLib.py

python 3.8

Python version for various C++ codes

Note: the various .cpp files are introduced in Reference: https://doi.org/10.1016/j.cpc.2021.107982

• 0.py (python version for 0.cpp)

#!/usr/bin/env python3

# python version for 0.cpp

from HepLib import *

x = symbol("x")
y = symbol("y")
z = symbol("z")
r = expr("2/3")
print("x ->", x, ", r ->", r)

e1 = r*x+2*y+pow(y,10)
print("e1 ->", e1)
e2 = (x+1)/(x-1)
print("e2 ->", e2)
e3 = sin(x+2*y)+3*z+41
print("e3 ->", e3)
e4 = e3+e2/exp(e1)
print("e4 ->", e4)
print()

x = symbol("x")
y = Symbol("y")
print("x ->", x, ", y->", y)
e1 = conjugate(x)
print("conjugate(x) ->", e1)
e2 = conjugate(y)
print("conjugate(y) ->", e2)
print()

• 1.py (python version for 1.cpp)

#!/usr/bin/env python3

# python version for 0.cpp

from HepLib import *

expr_str = "WF(1)+x(1)^2+sin(5)+power(a,n)"
e1 = expr(expr_str)
print(e1)
a = Symbol("a")
n = Symbol("n")
e2 = WF(expr(1))+pow(x(1),2)+sin(expr(5))+pow(a,n)
print(e1-e2)

print()

• 2.py (python version for2.cpp)

#!/usr/bin/env python3

# python version for 0.cpp

from HepLib import *

x = Symbol("x")
y = Symbol("y")
z = Symbol("z")
n = Symbol("n")
l1 = lst([x, y, x+z])
e1 = pow(sin(x),n)
print("l1 ->", l1)
print("e1 ->", e1)
tot = l1.nops()
print("l1.nops() ->", tot)
item1 = l1.op(0)
print("1st item of l1 ->", item1)
l1.let_op(2, e1)
print("updated l1 ->", l1)
tot = e1.nops()
print("e1.nops() ->", tot)
item2 = e1.op(1)
print("2nd item of e1 ->", item2)
print()

• 3.py (python version for3.cpp)

#!/usr/bin/env python3

# python version for 0.cpp

from HepLib import *

x = Symbol("x")
y = Symbol("y")
e0 = expr("x^4+x^3+x^2+x")
print("e0 ->", e0)
e1 = e0.subs([pow(x,w)>>pow(y,w+2)])
print("e1 ->",e1)

class mapClass(MapFunction):
    def map(self, e):
        if(e.match(pow(x,w)) and e.op(1).info("even")):
            return pow(y,e.op(1)+2)
        else:
            return e.map(self)

e2 = mapClass()(e0)
print("e2 ->",e2)

print()

• 4.py (python version for4.cpp)

#!/usr/bin/env python3

# python version for 0.cpp

from HepLib import *

x = Symbol("x")
y = Symbol("y")

data = exvec()
total = 100
for i in range(total):
    data.push_back(sin(exp(x+y*i)))

class ParFunRun(ParFun):
    def __call__(self, idx):
        ret = data[idx]
        ret = series(ret,x,5)
        #print("idx:", idx)
        return ret

f = ParFunRun()
set_Verbose(100)
set_Parallel_Process(4)
ret = Parallel(total, f)
#print(ret)

• 5.py (python version for5.cpp)

#!/usr/bin/env python3

# python version for 0.cpp

from HepLib import *

me = Symbol("me")
mm = Symbol("mm")
e = Symbol("e")
mu = Index("mu")
nu = Index("nu")
p = Vector("p")
P = Vector("P")
k = Vector("k")
K = Vector("K")
q = Vector("q")

letSP(p,me*me)
letSP(P,me*me)
letSP(k,mm*mm)
letSP(K,mm*mm)

def gpm(p, m):
    return GAS(p)+m*GAS(1)

tr1 = TR( gpm(P,-me)*GAS(mu)*gpm(p,me)*GAS(nu) );
tr2 = TR( gpm(k,mm)*GAS(mu)*gpm(K,-mm)*GAS(nu) );
res =  pow(e,4) / (4*pow(SP(q),2)) * tr1 * tr2;
set_form_using_dim4(True)
res = form(res)
res = factor(res);
print(res.subs(me>>0));
print()

set_form_using_su3(True)
a = IndexCA("a")
i = IndexCF("i")
j = IndexCF("j");
tr = TTR([a,a]);
print("tr1 =", form(tr))
tr = SUNT(a,i,j) * SUNT(a,j,i);
print("tr2 =", form(tr))
tr = SUNT([a,a],i,i)
print("tr3 =", form(tr))

print()

• 6.py (python version for6.cpp)

#!/usr/bin/env python3

from HepLib import *

A = Symbol("A")
e = Symbol("e")
ebar = Symbol("ebar")
mu = Symbol("mu")
mubar = Symbol("mubar")

p = Vector("p")
P = Vector("P")
k = Vector("k")
K = Vector("K")

me = Symbol("me")
mm = Symbol("mm")

proc = Process()
proc.Model = """
    [e, ebar, -]
    [mu, mubar, -]
    [A, A, +]
    [ebar, e, A]
    [mubar, mu, A]
"""
proc.In = "e[p],ebar[P]"
proc.Out = "mubar[k],mu[K]"
proc.Options = "onshell"
proc.Loops = 0

st = {}
st["p"] = p
st["P"] = P
st["k"] = k
st["K"] = K

amps = proc.Amplitudes(st)

set_InOutTeX(-1,"$e^-(p)$")
set_InOutTeX(-3,"$e^+(P)$")
set_InOutTeX(-2,"$\\mu^+(k)$")
set_InOutTeX(-4,"$\\mu^-(K)$")

set_LineTeX(Symbol("e"),"fermion, edge label=$e$")
set_LineTeX(Symbol("ebar"),"anti fermion, edge label=$e$")
set_LineTeX(Symbol("mu"),"fermion, edge label=$\\mu$")
set_LineTeX(Symbol("mbar"),"anti fermion, edge label=$\\mu$")
set_LineTeX(Symbol("A"),"photon, edge label=$\\gamma$")

Process.DrawPDF(amps, "amps.pdf")

class ClassFR(MapFunction):
    def __init__(self):
        MapFunction.__init__(self)

    def map(self, e):
        if(isFunction(e,"OutField") or isFunction(e,"InField")):
            return expr(1)
        elif(isFunction(e, "Propagator")):
            fi1 = e.op(0).op(1)
            fi2 = e.op(1).op(1)
            mom = e.op(2)
            if(e.op(0).op(0)==A):
                return (-I) * SP(LI(fi1),LI(fi2)) / SP(mom); # Feynman Gauge
            elif(e.op(0).op(0)==ebar):
                return I * Matrix(GAS(mom)+GAS(1)*me, DI(fi1),DI(fi2)) / (SP(mom)-me*me)
            elif(e.op(0).op(0)==mubar):
                return I * Matrix(GAS(mom)+GAS(1)*mm, DI(fi1),DI(fi2)) / (SP(mom)-mm*mm)
        elif(isFunction(e, "Vertex")):
            fi1 = e.op(0).op(1)
            fi2 = e.op(1).op(1)
            fi3 = e.op(2).op(1)
            if(e.op(0).op(0)==ebar):
                return I*Symbol("e")*Matrix(GAS(LI(fi3)),DI(fi1),DI(fi2))
            elif(e.op(0).op(0)==mubar):
                return I*Symbol("e")*Matrix(GAS(LI(fi3)),DI(fi1),DI(fi2))
        else:
            return e.map(self)

amps_FR = ClassFR()(amps[0])

print("amps_FR: ")
print(amps_FR);

ampL = amps_FR
ampR = IndexL2R(conjugate(ampL));
def SS1(p,m,i):
    return Matrix(GAS(p)+m*GAS(1),DI(i),RDI(i))
def SS2(p,m,i):
    return Matrix(GAS(p)+m*GAS(1),RDI(i),DI(i))

M2 = ampL * ampR * SS1(p,me,-1) * SS2(P,-me,-3) * SS1(k,-mm,-2) * SS2(K,mm,-4);
M2 = MatrixContract(M2);

print("M2: ")
print(M2)
print()

set_form_using_dim4(True)
letSP(p,me*me)
letSP(P,me*me)
letSP(k,mm*mm)
letSP(K,mm*mm)

res = form(M2);
print("Final M2:")
print(factor(res.subs(me>>0)))

• 7.py (python version for7.cpp)

#!/usr/bin/env python3

from HepLib import *

p = Vector("p")
q1 = Vector("q1")
expr = expr(1)/SP(q1) * expr(1)/(2*SP(p,q1)-SP(q1)) * expr(1)/(2*SP(p,q1)+SP(q1))
r = Apart(expr,[q1],[p]);
r1 = ApartIR2ex(r);
r2 = ApartIR2F(r);

print(r)
print()
print(r1)
print()
print(r2)
print()

r3 = Apart(expr, [SP(q1),SP(p,q1)])
print(r3)
print()

• 8.py (python version for8.cpp)

#!/usr/bin/env python3

import os
from HepLib import *

q1 = Symbol("q1")
p = Symbol("p")
m = Symbol("m")

fire = FIRE()
fire.Internal = exvec([ q1 ])
fire.External = exvec([ p ])
fire.Replacements = exvec([ p*p >> m*m ])

fire.Propagators = exvec([ q1*q1, 2*p*q1-q1*q1 ])
fire.Integrals = exvec([ expr("{2,1}") ])
fire.WorkingDir = "IBPdir";
fire.Reduce()

print("Reduced Rules:")
print(fire.Rules)
print("Master Integrals:")
print(fire.MIntegrals)

os.system("rm -rf IBPdir")

• 9.py (python version for9.cpp)

#!/usr/bin/env python3

import os
from HepLib import *

k = Symbol("k")
r = Symbol("r")
q = Symbol("q")
p1 = Symbol("p1")
p2 = Symbol("p2")
s = Symbol("s")

fp = FeynmanParameter()
fp.LoopMomenta = exvec([k,r,q]);
fp.Propagators= exvec([ -pow(k,2),-pow(k+p1+p2,2),-pow(-k+r,2),-pow(p1+r,2),-pow(k-q,2),-pow(p1+q,2) ]);
fp.Exponents = exvec(expr("{1+3*ep,1,1,1,1,1}"))
fp.lReplacements[p1*p1] = expr(0)
fp.lReplacements[p2*p2] = expr(0)
fp.lReplacements[p2*p1] = s/2
fp.lReplacements[s] = expr(-1)
fp.Prefactor = pow(I*pow(Pi,2-ep),-3) * pow(tgamma(1-ep),3)
work = SecDec()
set_Verbose(100)
work.Evaluate(fp)
print("Final Result & Error:")
print(work.VE)