Amplitude library (`alibrary.m`)

Tired of relying on large monolitic programs to generate your Feynman amplitudes? What you need is a library: alibrary.

This file deals with Feynman diagrams and all the related things. Generic Mathematica utils go into utils.m. Tests go into atestsuite.m.

First, load the other libraries, utils.m and library.m.

If[MatchQ[$InputFileName, s_String /; StringEndsQ[s, "alibrary.m"]],
  $Apath = FileNameDrop[$InputFileName] /. "" -> ".";
  ,
  $Apath = ".";
]

Get[$Apath <> "/utils.m"];
Get[$Apath <> "/library.m"];

Diagram generation

We use diagrams generated by QGraf and expect them to come in the following format (defined in qgraf-stylefile):

Diagram[id, sym-factor, {in-field ...}, {out-field ...}, {propagator ...}, {vertex ...}]

where

id is an arbitrary identified of the diagram;
sym-factor is the symmetry factor and the sign of the diagram;
in- and out-fields: F["field", field-idx, vertex-idx, momentum];
propagators: P["field", from-field-idx, to-field-idx, from-vertex-idx, to-vertex-idx, momentum]
vertices: V[vertex-idx, "fields", field-idx-1, momentum-1, field-idx-2, momentum-2, ...]

All the information here is directly as QGraf provides it, just packaged into a Mathematica expression.

To generate diagrams conveniently, use Diagrams.

`DiagramId[]`

Because Diagram[] objects are not Associations, here are some functions that help access their parts by names, instead of numbers. (In many places we will access these parts by numbers still; arguably those places should be changed to use these functions).

DiagramId[Diagram[id_, factor_, ifld_, ofld_, props_, verts_]] := id

DiagramSymmetryFactor[Diagram[_, factor_, _, _, _, _]] := Abs[factor]

DiagramIncomingFields[Diagram[_, _, ifld_List, _, _, _]] := ifld

DiagramOutgoingFields[Diagram[_, _, _, ofld_List, _, _]] := ofld

DiagramPropagators[Diagram[_, _, _, _, props_List, _]] := props

DiagramVertices[Diagram[_, _, _, _, _, verts_List]] := verts

`DiagramClosedLoops[]`

Return the number of closed loops comprised only of fields that match the given (string) pattern. Only really works for fields that always come in a pair in each vertex (e.g. fermions).

DiagramClosedLoops[fieldpat_] := DiagramFieldLoops[#, fieldpat]&
DiagramClosedLoops[Diagram[_, _, i_List, o_List, p_List, _], fieldpat_] :=
Module[{edges, SS, II, EE},
  edges = Join[
    i // Cases[ F[fieldpat, fi_, vi_, _] :> SS[fi] <-> II[vi] ],
    o // Cases[ F[fieldpat, fi_, vi_, _] :> II[vi] <-> EE[fi] ],
    p // Cases[ P[fieldpat, fi1_, fi2_, vi1_, vi2_, _] :> II[vi1] <-> II[vi2] ]
  ];
  edges // ConnectedComponents // Select[FreeQ[_SS|_EE]] // Length
]

`Diagrams[]`

Return a list of diagrams produced by QGraf with the given incoming and outgoing fields and the given number of loops.

The QGraf binary, QGraf model file, QGraf options, and the massless particle pattern are specified via the global variables $QGraf, $QGrafModel, $QGrafOptions, and $MasslessFieldPattern. These can be overriden via QGraf, QGrafModel, QGrafOptions, and MasslessFieldPattern options, but the idea is that a model file has already provided the global variables instead.

Diagrams[fieldsi_List, fieldso_List, loops_Integer, OptionsPattern[]] :=
Module[{model, options, nomasspat, tmpdir, tmpoutput, momi, momo, i, result},
  model = Replace[OptionValue[QGrafModel], None -> $QGrafModel];
  options = Replace[OptionValue[QGrafOptions], None -> $QGrafOptions];
  qgraf = Replace[OptionValue[QGraf], None -> $QGraf];
  nomasspat = Replace[OptionValue[MasslessFieldPattern], None -> $MasslessFieldPattern];
  FailUnless[MatchQ[model, _String]];
  tmpdir = MkTemp["qgraf", ""];
  EnsureDirectory[tmpdir];
  CopyFile[$Apath <> "/qgraf-stylefile", tmpdir <> "/stylefile"];
  MkFile[tmpdir <> "/modelfile", model];
  momi = If[Length[fieldsi] == 1, {"q"}, Table["q" <> ToString[i], {i, Length[fieldsi]}]];
  momo = If[Length[fieldso] == 1, {"q"}, Table["p" <> ToString[i], {i, Length[fieldso]}]];
  MkFile[tmpdir <> "/qgraf.dat",
    "output='output.m';\n",
    "style='stylefile';\n",
    "model='modelfile';\n",
    "in=", MapThread[{#1, "[", #2, "]"}&, {fieldsi, momi}] // Riffle[#, ", "]&, ";\n",
    "out=", MapThread[{#1, "[", #2, "]"}&, {fieldso, momo}] // Riffle[#, ", "]&, ";\n",
    "loops=", loops, ";\n",
    "loop_momentum=l;\n",
    "options=", options, ";\n",
    If[Length[fieldsi] > 1,
      Table[
        If[MatchQ[fieldsi[[i]], nomasspat], {"false=plink[", 1-2*i, "];\n"}, ""],
        {i, Length[fieldsi]}],
      ""],
    If[Length[fieldso] > 1,
      Table[
        If[MatchQ[fieldso[[i]], nomasspat], {"false=plink[", -2*i, "];\n"}, ""],
        {i, Length[fieldso]}],
      ""]
  ];
  SafeRun[MkString["cd '", tmpdir, "' && '", qgraf, "'"]];
  result = SafeGet[tmpdir <> "/output.m"];
  EnsureNoDirectory[tmpdir];
  result
]
Options[Diagrams] = {QGraf -> None, QGrafModel -> None, QGrafOptions -> None, MasslessFieldPattern -> None};

By default look for QGraf in PATH.

If[Not[MatchQ[$QGraf, _String]], $QGraf = "qgraf"; ];

Diagrams to graphs

`DiagramGraphEdges[]`

Convert a diagram to a list of undirected edges.

DiagramGraphEdges[Diagram[_, _, i_List, o_List, p_List, _]] := Join[
  p /. P[f_, fi1_, fi2_, vi1_, vi2_, mom_] :> Internal[vi1] <-> Internal[vi2],
  i /. F[f_, fi_, vi_, mom_] :> Incoming[fi] <-> Internal[vi],
  o /. F[f_, fi_, vi_, mom_] :> Internal[vi] <-> Outgoing[fi]
]

`DiagramGraph[]`

Convert a diagram to an undirected Graph object.

DiagramGraph[d_Diagram] := DiagramGraphEdges[d] // Graph
DiagramGraph[CutDiagram[d1_Diagram, d2_Diagram]] := Module[{e1, e2, x},
  e1 = DiagramGraphEdges[d1];
  e2 = DiagramGraphEdges[d2];
  Join[
    e1 // DeleteCases[_ <-> _Outgoing],
    e2 /.
      Incoming -> Incoming2 /.
      Internal -> Internal2 /.
      Cases[e1, (i_ <-> e_Outgoing) :> ((x_ <-> e) :> Style[x <-> i, Dashed])]
  ]
]

`DiagramXGraph[]`

Convert a diagram to an [[XGraph]] object, with informational labels on each edge.

DiagramXGraph[Diagram[_, _, i_List, o_List, p_List, _]] := Join[
   i /. F[f_, fi_, vi_, mom_] :> {IN[fi] -> vi,
      Text[MkString[f, "(", mom, ")"] // StringReplace[" " -> ""]],
      Gray, Thin},
   o /. F[f_, fi_, vi_, mom_] :> {vi -> OO[fi],
      Text[MkString[f, "(", mom, ")"] // StringReplace[" " -> ""]],
      Gray, Thin},
   p /. P[f_, fi1_, fi2_, vi1_, vi2_, mom_] :> {vi1 -> vi2,
      Text[MkString[f, "(", mom, ")"] // StringReplace[" " -> ""]],
      Switch[f, q, Thickness[0.01], _, {}]}
   ] // XGraph

`DiagramToGraphviz[]`

Convert a diagram to a source code for a Graphviz digraph object.

DiagramToGraphviz[Diagram[id_, _, i_List, o_List, p_List, _]] := Module[{},
  MkString[
   "digraph {\n",
   " fontsize=12; margin=0;\n",
   " node [shape=circle width=0.1 color=black];\n",
   " edge [fontsize=8; colorscheme=paired12];\n",
   i /. F[f_, fi_, vi_, mom_] :> fi // Union // Map[{" ", #, " [width=0.05 color=gray];\n"} &],
   o /. F[f_, fi_, vi_, mom_] :> fi // Union // Map[{" ", #, " [width=0.05 color=gray];\n"} &],
   i /. F[f_, fi_, vi_, mom_] :> {
     " ", fi, " -> ", vi,
     " [label=\"", f, "(", mom // ToString // StringReplace[" " -> ""], ")\",color=", FieldGraphvizColor[f] - 1, "];\n"},
   p /. P[f_, fi1_, fi2_, vi1_, vi2_, mom_] :> {
     " ", vi1, " -> ", vi2,
     " [label=\"", f, "(", mom // ToString // StringReplace[" " -> ""], ")\",color=", FieldGraphvizColor[f], ",style=bold];\n"},
   o /. F[f_, fi_, vi_, mom_] :> {
     " ", vi, " -> ", fi,
     " [label=\"", f, "(", mom // ToString // StringReplace[" " -> ""], ")\",color=", FieldGraphvizColor[f] - 1, "];\n"},
   "}\n"
   ]
]
DiagramToGraphviz[CutDiagram[
  Diagram[id1_, _, i1_List, o1_List, p1_List, v1_List], Diagram[id2_, _, i2_List, o2_List, p2_List, v2_List]
]] := Module[{},
  MkString[
   "digraph {\n",
   " fontsize=12; margin=0; label_scheme=2;\n",
   " node [shape=circle width=0.1 color=black];\n",
   " edge [fontsize=8; colorscheme=paired12];\n",
   i1 /. F[f_, fi_, vi_, mom_] :> fi // Union // Map[{" ", #, "01 [fontsize=10 width=0.05 color=gray label=\"", #, "\"];\n"} &],
   i2 /. F[f_, fi_, vi_, mom_] :> fi // Union // Map[{" ", #, "02 [fontsize=10 width=0.05 color=gray label=\"", #, "'\"];\n"} &],
   v1 /. V[id_, ___] :> id // Union // Map[{" ", #, "01 [label=\"", #, "\"];\n"} &],
   v2 /. V[id_, ___] :> id // Union // Map[{" ", #, "02 [label=\"", #, "'\"];\n"} &],
   o1 /. F[f_, fi_, vi_, mom_] :> fi // Union // Map[{" \"|edgelabel|", -#, "00\" [fontsize=10 width=0.05 shape=square color=gray label=\"", #, "\"];\n"} &],
   i1 /. F[f_, fi_, vi_, mom_] :> {
     " ", fi, "01 -> ", vi, "01",
     "[label=\"", f, "(", mom // ToString // StringReplace[" " -> ""], ")\",color=", FieldGraphvizColor[f] - 1, "];\n"
   },
   i2 /. F[f_, fi_, vi_, mom_] :> {
     " ", fi, "02 -> ", vi, "02",
     "[label=\"", f, "(", mom // ToString // StringReplace[" " -> ""], ")\",color=", FieldGraphvizColor[f] - 1, "];\n"
   },
   p1 /. P[f_, fi1_, fi2_, vi1_, vi2_, mom_] :> {
     " ", vi1, "01 -> ", vi2, "01",
     " [label=\"", f, "(", mom // ToString // StringReplace[" " -> ""], ")\",color=", FieldGraphvizColor[f], ",style=bold];\n"
   },
   p2 /. P[f_, fi1_, fi2_, vi1_, vi2_, mom_] :> {
     " ", vi1, "02 -> ", vi2, "02",
     " [label=\"", f, "(", mom // ToString // StringReplace[" " -> ""], ")\",color=", FieldGraphvizColor[f], ",style=bold];\n"
   },
   o1 /. F[f_, fi_, vi_, mom_] :> {
     " ", vi, "01 -> \"|edgelabel|", -fi, "00\"",
     " [label=\"", f, "(", mom // ToString // StringReplace[" " -> ""], ")\",color=", FieldGraphvizColor[f] - 1, "];\n"
   },
   o2 /. F[f_, fi_, vi_, mom_] :> {
     " ", vi, "02 -> \"|edgelabel|", -fi, "00\"",
     " [label=\"", f, "(", mom // ToString // StringReplace[" " -> ""], ")\",color=", FieldGraphvizColor[f] - 1, "];\n"
   },
   "}\n"
   ]
]

`FieldGraphvizColor[]`

Default Graphviz style for all fields. Styles for specific fields should be defined by the model files.

FieldGraphvizColor[field_] = 12;

`DiagramViz[]`

Use Graphviz to convert a diagram into a Graphics object. Useful for visualization.

DiagramViz[d:(_Diagram|_CutDiagram)] := Module[{tmp, pdf, result},
  tmp = MkTemp["diavis", ".gv"];
  pdf = tmp <> ".pdf";
  Export[tmp, DiagramToGraphviz[d], "String"];
  Run["neato -Tpdf -o", pdf, tmp];
  result = Import[pdf];
  DeleteFile[{tmp, pdf}];
  result // First
]

`GraphLayoutVertexCoordinates[]`

Take a graph defined by edges (pairs of nodes) and produce coordinates for all vertices by calling Graphviz (sfpd).

Example:

In[]:= GraphLayoutVertexCoordinates[{{1,2},{1,3},{2,3}}]
Out[]= {1 -> {272.32, 243.6},
        2 -> {27, 214.2},
        3 -> {175.3, 18}}

GraphLayoutVertexCoordinates[edges:{{_,_} ...}] :=
Module[{tmp, out, result, VertexToName, NameToVertex, nameidx},
  tmp = MkTemp["graph", ".gv"];
  out = tmp <> ".json";
  nameidx = 0;
  VertexToName[v_] := VertexToName[v] = (
    NameToVertex[ToString[nameidx]] = v;
    ToString[nameidx++]
  );
  MkFile[tmp,
    "graph {\n",
    edges // MapReplace[
      {a_, b_} :> {" ", VertexToName[a], " -- ", VertexToName[b], ";\n"}
    ],
    "}\n"
  ];
  Run["cat " <> tmp];
  Run["sfdp -Tjson -o", out, tmp];
  result = Import[out];
  DeleteFile[{tmp, out}];
  "objects" /. result //
    Map[
      NameToVertex["name" /. #] -> ToExpression[MkString["{", "pos" /. #, "}"]]&
    ]
]

`NormalizeCoordinates[]`

Scale and rotate the coordinates (in the format returned by GraphLayoutVertexCoordinates) so that the left vertices end up on the left, the right on the right, and the distance between them is scale.

NormalizeCoordinates[coordmap_, left_List, right_List, scale_] :=
Module[{l, r, angle, rotationmx, coordmap3, scalefactor},
  FailUnless[AllTrue[left, KeyExistsQ[coordmap, #]&]];
  FailUnless[AllTrue[right, KeyExistsQ[coordmap, #]&]];
  l = left /. coordmap // Mean;
  r = right /. coordmap // Mean;
  angle = r - l // Apply[ArcTan] // N;
  rotationmx = {{Cos[angle], Sin[angle]}, {-Sin[angle], Cos[angle]}};
  coordmap2 = coordmap // Association // Map[rotationmx.(# - l)&];
  l = left /. coordmap2 // Map[First] // Min;
  r = right /. coordmap2 // Map[First] // Max;
  scalefactor = scale/{r-l, Max[coordmap2[[;;,2]]] - Min[coordmap2[[;;,2]]]} // N;
  coordmap2 //
    Map[scalefactor (# - {(l+r)/2, 0})&] //
    Normal
]

TikZ interface for diagrams

`DiagramToTikZ[]`

Convert a diagram into TikZ, somewhat badly. It is expected that the user will manually edit the TikZ source afterwards using e.g. TikZiT.

DiagramToTikZ[Diagram[id_, _, i_List, o_List, p_List, v_List], OptionsPattern[]] :=
Module[{pstyles, ni, no, coords, scale, external},
  pstyles = OptionValue[PropagatorStyles] //
    If[# =!= None, #, p /. P[f_, ___] :> FieldTikZStyle[f]]&;
  FailUnless[Length[p] === Length[pstyles]];
  ni = Length[i];
  no = Length[o];
  scale = Max[ni,no,4];
  scale = {2.0,1.5};
  external = Join[
      i /. F[f_, fi_, vi_, mom_] :> Rule[fi, scale {-1/2 - 1/2/scale[[1]], If[ni<=1,0,(fi+1)/2/(ni-1)+1/2]} // N],
      o /. F[f_, fi_, vi_, mom_] :> Rule[fi, scale {+1/2 + 1/2/scale[[1]], If[no<=1,0,(fi+2)/2/(no-1)+1/2]} // N]
  ] // Association;
  coords = Graph[
      v /. V[vi_, __] :> vi,
      p /. P[f_, fi1_, fi2_, vi1_, vi2_, mom_] :> vi1 <-> vi2
    ] //
    GraphEmbedding[#, "SpringEmbedding"]& //
    MapThread[Rule, {v /. V[vi_, __] :> vi, #}]& //
    NormalizeCoordinates[#,
      i /. F[f_, fi_, vi_, mom_] :> vi,
      o /. F[f_, fi_, vi_, mom_] :> vi,
      scale]& //
    Association //
    If[
      Plus@@(Join[i, o] /. F[f_, fi_, vi_, mom_] :> EuclideanDistance[#[vi]*{1,+1}, external[fi]]^2) >
      Plus@@(Join[i, o] /. F[f_, fi_, vi_, mom_] :> EuclideanDistance[#[vi]*{1,-1}, external[fi]]^2),
      Map[#*{1,+1}&, #],
      Map[#*{1,-1}&, #]
    ]& //
    Join[external, #]&;
  MkString[
    "\\begin{tikzpicture}\n",
    "\t\\begin{pgfonlayer}{nodelayer}\n",
    i /. F[f_, fi_, vi_, mom_] :> {
      "\t\t\\node [style=none] (", fi, ") at (",
        coords[fi] // Round[#, 0.25]& // Map[FormatFixed[2]] // Riffle[#, ","]&,
      ") {};\n"
    },
    o /. F[f_, fi_, vi_, mom_] :> {
      "\t\t\\node [style=none] (", fi, ") at (",
        coords[fi] // Round[#, 0.25]& // Map[FormatFixed[2]] // Riffle[#, ","]&,
      ") {};\n"
    },
    v /. V[vi_, __] :> {
      "\t\t\\node [style=dot] (", vi, ") at (", 
        coords[vi] // Round[#, 0.25]& // Map[FormatFixed[2]] // Riffle[#, ","]&,
      ") {};\n"
    },
    "\t\\end{pgfonlayer}\n",
    "\t\\begin{pgfonlayer}{edgelayer}\n",
    i /. F[f_, fi_, vi_, mom_] :> {
      "\t\t\\draw [style=incoming ", FieldTikZStyle[f], "] (", fi, ") to (", vi, ");\n"
    },
    o /. F[f_, fi_, vi_, mom_] :> {
      "\t\t\\draw [style=outgoing ", FieldTikZStyle[f], "] (", vi, ") to (", fi, ");\n"
    },
    {p, pstyles} // Transpose // MapReplace[
      {P[f_, fi1_, fi2_, vi1_, vi2_, mom_], s_} :>
        {"\t\t\\draw [style=", s, "] (", vi1, ") to (", vi2, ");\n"}
    ],
    "\t\\end{pgfonlayer}\n",
    "\\end{tikzpicture}\n"
  ]
]
Options[DiagramToTikZ] = {PropagatorStyles -> None};

`FieldTikZStyle[]`

Default TikZ style for all fields. Styles for specific fields should be defined by the model files.

FieldTikZStyle[field_] = "edge";

`MkTikZGraph[]`

Create a TikZ graph file from the given edges and labels.

MkTikZGraph[filename_String, dx_, dy_, vleft_, vright_, edges:{{_, _, _} ...}, labels:{{_, _, _, _}...}] :=
Module[{vertexcoords, outer, inner, sx, sy, x1, y1, x2, y2, x, y, bends, b, e, i},
  vertexcoords = edges[[;;,;;2]] // GraphLayoutVertexCoordinates//MapAt[ToExpression, #, {;;,1}]&;
  (* Shorten outer legs by a half. *)
  outer = edges[[;;,;;2]] // Flatten // Union // Select[Count[edges, {#, _, _}|{_, #, _}] === 1&];
  inner = outer // Map[Cases[edges, {#, v_, _}|{v_, #, _} :> v]& /* First];
  vertexcoords = Join[
    MapThread[#1->(#2+#3)/2&, {outer, outer // Map[Replace[vertexcoords]], inner // Map[Replace[vertexcoords]]}],
    vertexcoords // DeleteCases[Rule[Alternatives @@ outer, _]]
  ];
  (* Scale the diagram to desired size *)
  sx = vertexcoords // #[[;;,2,1]]& // dx/(1 + Max[#] - Min[#])&;
  sy = vertexcoords // #[[;;,2,2]]& // dy/(1 + Max[#] - Min[#])&;
  {x1, y1} = vleft // Replace[vertexcoords];
  {x2, y2} = vright // Replace[vertexcoords];
  If[x2 < x1, sx = -sx];
  If[y2 < y1, sy = -sy];
  vertexcoords = vertexcoords // MapReplace[Rule[v_, {x_, y_}] :> Rule[v, {(x - x1)*sx, (y - y1)*sy }]];
  bends = <|
    1 -> {""},
    2 -> {", bend right", ", bend left"},
    3 -> {", bend right=45", "", ", bend left=45"},
    4 -> {", bend right=75, looseness=1.25", ", bend right", ", bend left", ", bend left=75, looseness=1.25"},
    5 -> {", bend right=75, looseness=1.25", ", bend right=45", "", ", bend left=45", ", bend left=75, looseness=1.25"}
  |>;
  MkFile[filename,
    "\\begin{tikzpicture}\n",
    "\t\\begin{pgfonlayer}{nodelayer}\n",
    vertexcoords // MapReplace[(v_ -> {x_, y_}) :>
      {"\t\t\\node [style=", v /. { Alternatives @@ outer -> "none", _ -> "dot"}, "] (", v, ") at (",
        x // Round[#, 0.25]& // NumberForm[#, {Infinity,3}]&, ", ",
        y // Round[#, 0.25]& // NumberForm[#, {Infinity,3}]&, ") {};\n"}
    ],
    labels // MapReplace[{v1_, v2_, style_, text_} :> (
      {x, y} = ((v1 // Replace[vertexcoords]) + (v2 // Replace[vertexcoords]))/2;
      {"\t\t\\node [style=", style, "] (", v1, ":", v2, ") at (",
        x // Round[#, 0.125]& // NumberForm[#, {Infinity,3}]&, ", ",
        y // Round[#, 0.125]& // NumberForm[#, {Infinity,3}]&, ") {", text, "};\n"}
    )],
    "	\\end{pgfonlayer}\n",
    "	\\begin{pgfonlayer}{edgelayer}\n",
    edges //
      GroupBy[Sort[#[[;;2]]]&] //
      Values //
      Map[Function[{edgs},
        b = bends[Length[edgs]];
        Table[
          e = edgs[[i]];
          {"\t\t\\draw [style=", e[[3]], b[[If[Sort[e[[;;2]]]===e[[;;2]], i, Length[edgs]+1-i]]], If[e[[1]] === e[[2]], ", loop", ""], "] (", e[[1]], ") to (", e[[2]], ");\n"},
          {i, Length[edgs]}
        ]
      ]],
    "	\\end{pgfonlayer}\n",
    "\\end{tikzpicture}\n"
  ];
];

`MkIntegralTikZ[]`

Create a TikZ file for an integral with given indices. The indices should come in the same order as the propagators of the diagram. Only non-negative indices work at the moment.

MkIntegralTikZ[filename_String, Diagram[_, _, ifld_List, ofld_List, props_List, verts_List], indices_List] :=
Module[{toremove, shrinkgr, vimap},
  FailUnless[Length[props] === Length[indices]];
  toremove = indices // PositionIndex // Lookup[#, 0, {}]&;
  vimap = Graph[verts[[;;,1]], props[[toremove, 4;;5]] // Map[Apply[UndirectedEdge]]] //
    ConnectedComponents //
    Select[Length[#] > 1&] //
    Map[Sort] //
    Map[Alternatives @@ # -> First[#]&];
  field2style = {"t" -> "top", _ -> "edge"};
  ifield2style = {"t"|"T" -> "incoming top", _ -> "incoming"};
  ofield2style = {"t"|"T" -> "outgoing top", "H" -> "outgoing higgs", _ -> "outgoing"};
  MkTikZGraph[filename, 3.0, 2.0, -1, -2, Join[
      ifld /. F[f_, fi_, vi_, mom_] :> {fi, vi /. vimap, f /. ifield2style},
      ofld /. F[f_, fi_, vi_, mom_] :> {vi /. vimap, fi, f /. ofield2style},
      Transpose[{props, indices}] // Map[Replace[{
        {_P, 0} :> Nothing,
        {P[f_, fi1_, fi2_, vi1_, vi2_, mom_], 1} :> {vi1 /. vimap, vi2 /. vimap, f /. field2style},
        {P[f_, fi1_, fi2_, vi1_, vi2_, mom_], idx_} :> {vi1 /. vimap, vi2 /. vimap, (f /. field2style) <> ",style=edge dot" <> ToString[idx-1]}
      }]]
    ],
    {}
  ];
]

IBP Bases

`ExpandScalarProducts[]`

Expand sp[..., ...] and take out constant factors, so that only sp between momenta are left. To know which symbols are momenta, this function takes a list of symbols, or a pattern.

ExpandScalarProducts[momnames_List] := ExpandScalarProducts[Alternatives @@ momnames]
ExpandScalarProducts[mompattern_] := ReplaceAll[sp[a_, b_] :> (
  Expand[a*b] // ReplaceAll[{
    (l:mompattern) (k:mompattern) :> Sort[sp[l, k]],
    (l:mompattern)^2 :> sp[l, l]
  }]
)]

`BToDen[]`

Convert B notation back to a product of den.

BToDen[ex_, bases_List] := Module[{bid2basis},
  bid2basis = bases // GroupBy[#["id"]&] // Map[Only];
  ex /.
    B[bid_, idx__] :> (bid2basis[bid]["denominators"]^{idx} // Apply[Times]) /.
    den[p_]^n_?Negative :> sp[p, p]^(-n) /.
    den[p_, m_]^n_?Negative :> (sp[p, p] - m)^(-n) /.
    den[p_, m_, irr]^n_?Negative :> (sp[p, p] - m)^(-n)
]
BToDen[bases_List] := BToDen[#, bases] &

`ToB[]`

Convert an expression with sp and den to B notation in a given basis. This is the slow version of it; the faster one uses FORM: see RunThroughForm and FormCallToB.

ToB[basis_Association] := ToB[#, basis]&
ToB[ex_, basis_Association] := Module[{indices},
  ex //
    ExpandScalarProducts[Alternatives @@ Join[basis["externalmom"], basis["loopmom"]]] //
    ReplaceAll[basis["nummap"]] //
    ReplaceAll[basis["denmap"]] //
    ReplaceAll[basis["sprules"]] //
    (*Together //*)
    Bracket[#, _DEN, #&, (
      indices = Table[0, Length[basis["denominators"]]];
      # /. DEN[i_]^n_. :> (indices[[i]] += n; 1) // #* B[basis["id"], Sequence @@ indices] &
    ) &]&
]

`IBPRelations[]`

List IBP relations for a basis, return a list of equations with terms of the form B[id, n[1], n[2], ...] and coefficients that depend on n[i].

IBPRelations[basis_Association] := Module[{dens, i, l, v},
  dens = basis["denominators"];
  Table[
    Product[DEN[i]^(n[i]), {i, Length[dens]}]*(
        Sum[-2 n[i] D[dens[[i, 1]], l] sp[v, dens[[i, 1]]] DEN[i], {i, Length[dens]}] +
        If[v === l, d, 0]
      )
    ,
    {l, basis["loopmom"]},
    {v, Join[basis["externalmom"], basis["loopmom"]]}
  ] // Flatten // Map[ToB[basis]] // Bracket[#, _B, Together]&
]

`VariableDimensions[]`

Figure out the dimensionality of variables in an expression (or a list of expressions) of a given dimensionality.

E.g.:

In[1]:= VariableDimensions[m + m^2/q, 1]
Out[1]= {m -> 1, q -> 1}

VariableDimensions[{}, _] := {}
VariableDimensions[expression_, dimension_] :=
Module[{DimOf, DimOfSymbol, ex, eqns, solution},
  DimOf[ex_List] := (
    ex // Map[Sow[DimOf[#] == DimOf[ex[[1]]]]&];
    DimOf[ex[[1]]]
  );
  DimOf[ex_Plus] := (
    ex // Apply[List] // Map[Sow[DimOf[#] == DimOf[ex[[1]]]]&];
    DimOf[ex[[1]]]
  );
  DimOf[ex_Times] := ex // Apply[List] // Map[DimOf] // Apply[Plus];
  DimOf[ex_^n_] := DimOf[ex]*n;
  DimOf[ex_?NumberQ] := 0;
  DimOf[ex_Symbol] := DimOfSymbol[ex];
  DimOf[ex_] := Error["Don't know the dimension of: ", ex];
  eqns = Reap[Sow[DimOf[expression // DeleteCases[0]] == dimension];][[2,1]] // Apply[And];
  If[eqns === False, Error["The dimension of ", expression, " can't be ", dimension]];
  solution = eqns // Solve[#, # // CaseUnion[_DimOfSymbol]]&;
  If[solution === {}, Error["Inconsistent dimension in ", expression]];
  solution[[1]] /. DimOfSymbol[ex_] :> ex
]

`NormalizeDens[]`

Normalize the den[] expressions by dropping the leading signs of the momenta.

NormalizeDens[ex_] := ex /. den[p_, x___] :> den[DropLeadingSign[Expand[p]], x] /. den[p_, 0] :> den[p]

`CompleteIBPBasis[]`

Complete a list of denominators to a full IBP basis, and return an IBP basis object with all the information. This object is used throughout this library, and this is the main way to create it.

You are advised to NormalizeDens the denominators and sort them before calling this function.

Example:

In[]:= CompleteIBPBasis[1, {den[l1], den[q-l1]}, {l1, l2}, {q}, {sp[q,q]->qq}]
Out[]=
<|
 "id" -> 1,
 "loopmom" -> {l1, l2},
 "externalmom" -> {q},
 "sprules" -> {sp[q,q] -> qq},
 "invariants" -> {qq},
 "denominators" -> { den[l1], den[l1-q], den[l2,0,irr], den[l1+l2,0,irr], den[l2+q,0,irr] },
 "denmap" -> <| den[l1] -> DEN[1], ... |>,
 "nummap" -> <| sp[l1,l1] -> 1/DEN[1], ... |>
|>

CompleteIBPBasis[bid_, denominators_List, loopmom_List, extmom_List, sprules_List] :=
Module[{L, M, p, k, nums, vars, c, mx, candidatemoms, denadd, numadd, cadd, mxadd, dens, Complete, rels},
  L = loopmom // Map[DropLeadingSign] // Apply[Alternatives];
  M = Join[loopmom, extmom] // Map[DropLeadingSign] // Apply[Alternatives];
  dens = denominators // NormalizeDens;
  nums = dens /. den[p_] :> p^2 /. den[p_, m2_, ___] :> p^2 - m2 // Expand;
  nums = nums /. (l:M) (k:M) :> Sort[sp[l, k]] /. (l:M)^2 :> sp[l, l] /. sprules;
  vars = Tuples[{loopmom, Join[loopmom, extmom]}] //
    Map[DropLeadingSign /* Sort /* Apply[sp]] //
    Union;
  If[nums =!= {},
    {c, mx} = CoefficientArrays[nums, vars] // Normal;
    (* nums == c + mx.vars *)
    If[MatrixRank[mx] =!= Length[mx],
      Print["nums=",nums];
      Print["vars=",vars];
      Print[c];
      Print[mx];
      Print[MatrixRank[mx]];
      Error["CompleteIBPBasis: denominators ", denominators, " are already linearly dependent"]
    ];
    ,
    {c, mx} = {{}, {}};
  ];
  candidatemoms =
    Subsets[Join[loopmom, extmom], {1, Infinity}] //
    (* Higher counts almost always give more IBP terms; the
     * minimal amount is sufficient in our case (but not in
     * general). *)
    Select[Length[#] < 3&] //
    Join[#, Cases[#, {a_, b_} :> {a, -b}]]& //
    (*
    Join[#, Cases[#, {a_, b_, c_} :> {a, -b, c}]]& //
    Join[#, Cases[#, {a_, b_, c_} :> {a, b, -c}]]& //
    Join[#, Cases[#, {a_, b_, c_} :> {a, -b, -c}]]& //
    *)
    Map[Apply[Plus] /* DropLeadingSign] //
    Select[NotFreeQ[L]];
  Complete[dens_, mx_, c_, mini_] := Module[{i},
    If[Length[mx] === Length[vars],
      Sow[
        <|
          "id" -> bid,
          "loopmom" -> (loopmom // DropLeadingSign),
          "externalmom" -> (extmom // DropLeadingSign),
          "sprules" -> sprules,
          "denominators" -> dens,
          "denmap" -> (
            MapIndexed[#1 -> DEN @@ #2 &, dens] //
            DeleteCases[den[_, _, irr] -> _] //
            ReplaceAll[(den[p_, x___] -> y_) :> {den[p, x] -> y, den[-p, x] -> y}] //
            Flatten //
            Association
          ),
          "nummap" -> (
            Inverse[mx].(Map[1/DEN[#] &, Range[Length[mx]]] - c) //
            Bracket[#, _DEN, Factor]& //
            MapThread[Rule, {vars, #}]& //
            Join[#, # /. sp[a_, b_] :> sp[b, a]]& //
            Association
          ),
          "invariants" -> ({dens // Cases[den[_, m_, ___] :> m], sprules[[;;,2]]} // CaseUnion[_Symbol])
        |>
      ]
      ,
      For[i = mini, i < Length[candidatemoms], i++,
        numadd = Expand[candidatemoms[[i]]^2];
        numadd = numadd /. (l:M) (k:M) :> Sort[sp[l, k]] /. (l:M)^2 :> sp[l, l] /. sprules;
        {cadd, mxadd} = CoefficientArrays[numadd, vars] // Normal;
        If[MatrixRank[Append[mx, mxadd]] === Length[mx] + 1,
          Complete[
            Append[dens, den[candidatemoms[[i]], 0, irr]],
            Append[mx, mxadd],
            Append[c, cadd],
            i+1]
        ]
      ]
    ]
  ];
  Reap[Complete[dens, mx, c, 1]][[2]] // Only //
  MinimalBy[(
    rels = IBPRelations[#] // Bracket[#, _B, CO]& // Map[Terms /* Length];
    {Plus@@rels, Plus@@(rels*rels), #["denominators"][[;;,1]]//Map[Terms/*Length]//Apply[Plus]}
  )&]//First
]

`IBPBasisMapInvariants[]`

Apply a replacement map to the invariants of a basis.

IBPBasisMapInvariants[basis_Association, invmap_] :=
  basis //
    Append["denominators" -> (basis["denominators"] /. invmap)] //
    Append["sprules" -> (basis["sprules"] /. invmap)] //
    Append["invariants" -> (basis["invariants"] /. invmap // CaseUnion[_Symbol])]
IBPBasisMapInvariants[bases_List, invmap_] := bases // Map[IBPBasisMapInvariants[#, invmap]&]

`InvariantMapUnderMomentaPermutation[]`

Return a list of substitutions of the invariants of a basis corresponding to an external momenta permutation.

InvariantMapUnderMomentaPermutation[basis_Association, momperm_] :=
Module[{extmom, sprules, i, j, sps, v1, v2, vars, OLD, NEW, x},
  extmom = basis["externalmom"];
  sprules = basis["sprules"];
  sps = Table[
    Sort[sp[extmom[[i]], extmom[[j]]]]
    ,
    {i, 1, Length[extmom]},
    {j, 1, Length[extmom]}
  ] // Apply[Join] // Union;
  vars = basis["invariants"];
  v1 = sps /. sprules /. x:(Alternatives@@vars) :> OLD[x];
  v2 = sps /. momperm // Map[Sort] // ReplaceAll[sprules] // ReplaceAll[x:(Alternatives@@vars) :> NEW[x]];
  v1 - v2 //
    Map[#==0&] //
    Solve[#, vars // Map[OLD]]& //
    Only //
    ReplaceAll[(NEW|OLD)[x_] :> x] //
    DeleteCases[x_ -> x_]
]

`IBPBasisMassDimensions[]`

Return a map from the invariants of a basis to their mass dimensions.

IBPBasisMassDimensions[basis_Association] := IBPBasisMassDimensions[{basis}]
IBPBasisMassDimensions[bases_List] :=
  {
    bases[[;;, "sprules", ;;, 2]],
    bases[[;;, "denominators"]] // Map[Cases[den[_, m_, ___] :> m]]
  } //
  Flatten //
  DeleteCases[0] //
  VariableDimensions[#, 2]&

`IBPBasisCross[]`

Return a copy of the basis, but with external momenta swapped inside the denominator list. No change otherwise, i.e. the external scalar product rules remain the same.

IBPBasisCross[bid_, basis_Association, externalmommap_] :=
  Append[
    CompleteIBPBasis[bid, basis["denominators"] /. externalmommap, basis["loopmom"], basis["externalmom"], basis["sprules"]],
    "invariants" -> basis["invariants"]
  ]

`IBPBasisSameQ[]`

Check if two bases are identical, up to the order of keys.

IBPBasisSameQ[b1_, b2_] := Sort[Normal[b1]] === Sort[Normal[b2]]

Feynson interface for integral symmetries

By default look for Feynson in PATH.

If[Not[MatchQ[$Feynson, _String]], $Feynson = "feynson -q -j4"; ];

`ZeroSectors[]`

Calculate the zero sectors of a given basis. Return a list, each element being B[basis-id, (1|0), ...], listing the topmost zero sectors.

ZeroSectors[basis_Association] :=
  RunThrough[$Feynson <> " zero-sectors -sj3 -", {
      basis["denominators"] /. {
        den[p_] :> p^2,
        den[p_, m_] :> p^2-m,
        den[p_, m_, irr] :> p^2-m,
        den[p_, m_, cut] :> p^2-m-CUT
      },
      basis["denominators"] /. den[_, _, cut] -> 1 /. den[___] -> 0,
      basis["loopmom"],
      basis["sprules"] /. Rule->List /. sp -> Times
    }] //
    Map[B[basis["id"], Sequence@@SectorIdToIndices[#, Length[basis["denominators"]]]]&]
ZeroSectors[bases_List] := bases // Map[ZeroSectors] // Apply[Join]

`ZeroSectorPattern[]`

Return a pattern that matches zero intergals (in the B notation) of a given basis.

ZeroSectorPattern[basis_Association] := ZeroSectors[basis] //
  MapReplace[B[bid_, idx__] :> B[bid, {idx} /. 1 -> _ /. 0 -> _?NonPositive // Apply[Sequence]]] //
  Apply[Alternatives]
ZeroSectorPattern[bases_List] := bases // Map[ZeroSectorPattern] // Apply[Alternatives]

`SymmetryMaps[]`

Return a list of momenta maps, such that applying them to the list of feynman integral families makes symmetries and subtopology relations explicit. So, families that are symmetric will have identical sets of denominators after the maps are applied. Families that are symmetric to a subtopology of a bigger family will have a subsets of the denominators.

The families are defined by their set of den[]s.

The latter families are guaranteed to be mapped to the former ones.

You can use UniqueSupertopologyMapping to figure out the topmost supertopologies after this.

SymmetryMaps[families_List, loopmom_List, sprules_] :=
Module[{densets, uniqdensets, densetindices, uniqdensetmaps},
  densets = families // NormalizeDens // Map[
    CaseUnion[_den] /* Union /* Select[NotFreeQ[Alternatives@@loopmom]]
  ];
  densets = densets /. den[p_] :> p^2 /. den[p_, m_] :> p^2-m /. den[p_, m_, cut] :> p^2-m-CUT;
  RunThrough[$Feynson <> " symmetrize -", {densets, loopmom, sprules // Map[Apply[List]]}] // Map[Map[Apply[Rule]]]
];
SymmetryMaps[families_List, loopmom_List] := SymmetryMaps[families, loopmom, {}]

`IntegralFamilyMappingRules[]`

Determine if the latter families of integrals can be expressed in terms of the earlier ones. For each family (defined by a list of den[] expression) return either {} if it is unique, or {fam, {n1, n2, ...}}, meaning that any integral in the given family with indices {i_1,i_2,...} is equal to an integral in the family number fam with indices {i_n1,i_n2,...}.

IntegralFamilyMappingRules[densets_List, loopmom_List, sprules_List] := Module[{},
  RunThrough[$Feynson <> " -d mapping-rules -", {
    densets /.
      den[p_] :> p^2 /.
      den[p_, m_] :> p^2-m /.
      den[p_, m_, irr] :> p^2-m /.
      den[p_, m_, cut] :> p^2-m-CUT,
    loopmom, sprules // Map[Apply[List]]
  }]
]

`IntegralEqualitySets[]`

Take a list of integrals (in the B notation) and a list of bases, determine which integrals are equal, and return a list of index sets, with each integral in a set being equal to each other.

IntegralEqualitySets[integrals_List, bases_List] := Module[{bid2basis, densets},
  bid2basis = bases // GroupBy[#["id"]&] // Map[Only];
  densets = integrals //
    MapReplace[
      B[bid_, idx___] :> MapThread[Which[
        #2 === 0, Nothing,
        #2 === 1, #1,
        #2 > 0, #1 - MkExpression["POW", #2],
        #2 < 0, #1 - MkExpression["POWm", -#2],
        True, Error["Unrecognized index:", #2]
      ]&, {bid2basis[bid]["denominators"], {idx}}]
    ];
  densets = densets /. den[p_] :> p^2 /. den[p_, m_] :> p^2-m /. den[p_, m_, cut] :> p^2-m-CUT;
  IntegralFamilyMappingRules[
    densets,
    bases[[1, "loopmom"]],
    bases[[1, "sprules"]] /. sp[p1_,p2_] :> p1*p2 // Map[Apply[List]]
  ] //
    MapIndexed[Replace[#1, {} -> {First[#2], densets[[First[#2]]] // Length // Range}]&]//
    MapAt[Replace[Except[0] -> 1], #, {;;, 2, ;;}]& //
    PositionIndex //
    Values
]

`UniqueIntegralIndices[]`

Take a list of integrals (in the B notation) and a list of bases, and return the indices of unique integrals, discarding the symmetric duplicates.

UniqueIntegralIndices[integrals_List, bases_List] :=
  IntegralEqualitySets[integrals, bases] // Map[Sort] // Map[First]

`IntegralUnion[]`

Take a list of integrals (in the B notation) and a list of bases, and return a list of integrals with all symmetric duplicates removed.

IntegralUnion[integrals_List, bases_List] :=
  integrals[[UniqueIntegralIndices[integrals, bases]]]

`IntegralMapOntoSectors[]`

Try to map the given integrals onto their symmetric equivalents in the given sectors. Return a map from the given integrals to integrals in the specified sectors.

Negative indices are not supported here.

IntegralMapOntoSectors[sectors_List, integrals_List, bases_List] :=
Module[{bid2basis, fullfamidx, n, shortfams, intdensets, intidx, idx, bid},
  bid2basis = bases // GroupBy[#["id"]&] // Map[Only];
  fullfamidx = sectors // MapReplace[
    B[bid_, idx__] :> (n=0; {idx} // Map[If[#===0,0,n+=1;n]&])
  ];
  shortfams = sectors // MapReplace[
    B[bid_, idx__] :> (MapThread[If[#1===0,Nothing,#2]&, {{idx}, bid2basis[bid]["denominators"]}])
  ];
  {intdensets, intidx} = integrals // MapReplace[
      B[bid_, idx__] :> {
        MapThread[If[#===0, Nothing, #2]&, {{idx}, bid2basis[bid]["denominators"]}],
        {idx} // DeleteCases[0]
      }
    ] //
    Transpose;
  IntegralFamilyMappingRules[
    Join[shortfams, intdensets] /.
      den[p_] :> p^2 /.
      den[p_, m_] :> p^2-m /.
      den[p_, m_, irr] :> p^2-m /.
      den[p_, m_, cut] :> p^2-m-CUT,
    bases[[1, "loopmom"]],
    bases[[1, "sprules"]] /. sp[p1_,p2_] :> p1*p2 // Map[Apply[List]]
  ][[Length[sectors]+1;;]] //
    MapIndexed1[Function[{rule, i},
      If[And[rule =!= {}, rule[[1]] <= Length[sectors]],
        idx = rule[[2]];
        bid = sectors[[rule[[1]], 1]];
        integrals[[i]] -> B[bid, intidx[[i]] //
          Prepend[0] //
          #[[idx + 1]]& //
          Prepend[0] //
          #[[fullfamidx[[rule[[1]]]] + 1]]& //
          Apply[Sequence]
        ]
        ,
        Nothing
      ]
    ]]
]

IntegralMapOntoSector[sector_B, integrals_List, bases_List] :=
  IntegralMapOntoSectors[{sector}, integrals, bases]

`IntegralMapOntoFamily[]`

Try to map the given integrals onto a given basis using symmetries. Return a map from the given integrals to integrals of the given basis.

Negative indices are not supported here.

IntegralMapOntoFamily[basis_Association, integrals_List, bases_List] :=
  MapIntegralsOntoSectors[
    B[basis["bid"], basis["denominators"] // MapReplace[{den[__,irr]->0, _den->1}]],
    integrals,
    bases]

FORM interface for integrand transformation

If[Not[MatchQ[$FORM, _String]], $FORM = "tform -w4"; ];

`RunThroughForm[]`

Run expressions through FORM (using library.frm), running the specified code fragment on it. (The code is constructed with MkString).

RunThroughForm[{}, _] := {}
RunThroughForm[exprs_List, code_] :=
 Module[{tmpsrc, tmpdst, tmplog, result, toform, fromform, i, expridxs},
  tmpsrc = MkTemp["amp", ".frm"];
  tmpdst = tmpsrc <> ".m";
  tmplog = tmpsrc // StringReplace[".frm" ~~ EndOfString -> ".log"];
  {toform, fromform} = AmpFormIndexMaps[exprs];
  MkFile[tmpsrc,
    "#include ", $Apath, "/library.frm\n",
    (* This whole dance is needed to distribute expressions
     * across workers; tform would keep everything in a single
     * process if we where to just do assign our expression to
     * EXPR directly.
     *)
    "Table EXTBL(1:", Length[exprs], ");\n",
    Table[
      {"Fill EXTBL(", i, ") = (",
        (* Replacing delta() with d_() doesn't work because
         * d_()^2 is broken. Instead use (d_()), which works.
         * See: github.com/vermaseren/form/issues/341.
         *)
        exprs[[i]] // ReplaceAll[toform] // AmpToForm // (*StringReplace["delta("->"d_("]*)
          StringReplace["delta(" ~~ a:(Except[")"] ...) ~~ ")" -> "(d_(" ~~ a ~~ "))"],
        ");\n"}
      ,
      {i, Length[exprs]}
    ],
    (*"L EXPR = <EX(1)>+...+<EX(", Length[exprs], ")>;\n",*)
    "L EXPR = sum_(xidx,1,", Length[exprs], ",EX(xidx));\n",
    ".sort:init;\n",
    (* Now that EX(n) are in different workers, we can insert
     * their values, and start the computations.
     *)
    "id EX(x?) = EX(x)*EXTBL(x);\n",
    (*"cleartable EXTBL;\n"*)
    "#call input\n",
    code,
    "#call output(", tmpdst, ")\n",
    ".end\n"
  ];
  Print["RunThroughForm: calling ", $FORM, " ", tmpsrc];
  Run[MkString[
    "env ",
      "FORMTMP=\"${FORMTMP:-$TMP}\" ",
      "FORMTMPSORT=\"${FORMTMPSORT:-$TMP}\" ",
      "FORMPATH='", $Apath, "' ",
    $FORM, " -q -Z -M -l '", tmpsrc, "'"]];
  Print["RunThroughForm: reading result (", FileByteCount[tmpdst]//FormatBytes, ")"];
  result = SafeGet[tmpdst]//TM;
  DeleteFile[{tmpsrc, tmpdst, tmplog}];
  Print["RunThroughForm: transforming it back"];
  result /. fromform // AmpFromForm // Terms //
      Map[Replace[EX[i_]*ex_. :> {i, ex}]] //
     GroupBy[#, First -> (#[[2]] &)] & // Map[Apply[Plus]] //
   Lookup[#, Range[Length[exprs]], 0] &
  ]
RunThroughForm[code_] := RunThroughForm[#, code]&
RunThroughForm[exprs_, code_] := RunThroughForm[{exprs}, code] // Only

FormCall[procname_String] := {"#call ", procname, "\n"}
FormCall[procnames__] := Map[FormCall, {procnames}]

FormCallZeroSectors[zerobs_List] := zerobs //
  MapReplace[B[bid_, idx__] :> {
    "id B(", bid, MapIndexed[ReplaceAll[#1, {
      0 -> {", x", #2, "?neg0_"},
      1 -> {", x", #2, "?"}
    }]&, {idx}], ") = 0;\n"
  }]

FormCallReplace[rules__] := {rules} // MapReplace[
  (sp[p_] -> v_) :> {
    "id sp(", p // InputForm, ") = ", v // InputForm, ";\n",
    "id dot(", p // InputForm, ",", p // InputForm, ") = ", v // InputForm, ";\n"
  },
  (sp[p1_,p2_] -> v_) :> {
    "id sp(", p1 // InputForm, ",", p2 // InputForm, ") = ", v // InputForm, ";\n",
    "id sp(", p2 // InputForm, ",", p1 // InputForm, ") = ", v // InputForm, ";\n",
    "id dot(", p1 // InputForm, ",", p2 // InputForm, ") = ", v // InputForm, ";\n"
  }
]

`FormCallToB[]`

FORM function to convert the current expression into B notation. To be used with RunThroughForm.

FormCallToB[bases_List] := MkString[
    "#procedure toBID\n",
    "* Assume BID^n factor are already supplied.\n",
    "#endprocedure\n",
    "#procedure toDEN\n",
    "  ",
    bases // Map[Function[{basis},
      {
        "if (match(only, BID^", basis["id"], "));\n",
        basis["denmap"] // Normal //
          Map[{
            "    id ", #[[1]]//AmpToForm, " = ",
            #[[2]] /. basis["sprules"] /. DEN[n_] :> MkString["DEN", n] // AmpToForm,
            ";\n"
          }&] // Union,
        basis["nummap"] // Normal //
          Map[{
            "    id ", #[[1]] /. sp->(Dot/*Sort) //AmpToForm, " = ",
            #[[2]] /. basis["sprules"] /. DEN[n_] :> MkString["DEN", n] // AmpToForm,
            ";\n"
          }&] // Union,
        basis["sprules"] // Normal //
          Map[{
            "    id ", #[[1]] /. sp->(Dot/*Sort) // AmpToForm, " = ",
            #[[2]] /. DEN[n_] :> MkString["DEN", n] // AmpToForm,
            ";\n"
          }&] // Union
      }
    ]] // Riffle[#, "  else"]&,
    "  else;\n",
    "    exit \"ERROR: toDEN: got a term without a proper BID^n factor.\";\n",
    "  endif;\n",
    "#endprocedure\n",
    "#call toB(", Length[bases[[1, "denominators"]]], ", toBID, toDEN)\n"
  ]

Kira interface for IBP reduction

`KiraBasisName[]`

Kira sorts bases by name instead of adhering to the order of definition. We shall make sure that both the numerical and the lexicographic orders match, which will prevent Kira from messing it up.

KiraBasisName[bid_Integer] := MkString["b", IntegerDigits[bid, 10, 5]]
KiraBasisName[name_String] := name
KiraBasisName[DimShift[name_, n_]] := MkString[KiraBasisName[name], "_dimshift", n]

`MkKiraKinematicsYaml[]`

Create Kira’s kinematics.yaml config file.

MkKiraKinematicsYaml[filename_, extmom_List, sprules_List, variabledimensions_List, one_] :=
  MaybeMkFile[filename,
    "kinematics:\n",
    " incoming_momenta: [", extmom // Riffle[#, ", "]&, "]\n",
    " kinematic_invariants:\n",
    variabledimensions //
      ReplaceAll[(var_ -> dim_) :> {"  - [", var , ", ", dim, "]\n"}],
    " scalarproduct_rules:\n",
    sprules //
      ReplaceAll[sp -> (sp /* Sort)] //
      MapReplace[
      (sp[p_] -> v_) :> {"  - [[", p//InputForm, ",", p//InputForm, "], ", v//InputForm, "]\n"},
      (sp[p1_, p2_] -> v_) :> {"  - [[", p1//InputForm, ",", p2//InputForm, "], ", v//InputForm, "]\n"}
    ] // Union,
    If[one =!= None,
      {" symbol_to_replace_by_one: ", one, "\n"}
      ,
      {"# symbol_to_replace_by_one: ", variabledimensions[[1,1]], "\n"}
    ]
  ];

`MkKiraIntegralFamiliesYaml[]`

Create Kira’s integralfamilies.yaml config file.

MkKiraIntegralFamiliesYaml[filename_, bases_List] :=
Module[{loopmom, extmom, dens, basis},
  MaybeMkFile[filename,
    "integralfamilies:\n",
    Table[
      loopmom = basis["loopmom"];
      extmom = basis["externalmom"];
      dens = basis["denominators"];
      {
        "  - name: \"", KiraBasisName[basis["id"]], "\"\n",
        "    loop_momenta: [", Riffle[loopmom, ", "], "]\n",
        "    top_level_sectors: [b", dens // MapReplace[den[_, _, irr] -> 0, _den -> 1], "]\n",
        "    propagators:\n",
        dens // Map[Replace[{
          den[p_] | den[p_, 0, ___] :> {"      - [\"", CForm[p], "\", 0]\n"},
          den[p_, m_, ___] :> {"      - [\"", CForm[p], "\", \"", CForm[m /. sp[q] -> qq], "\"]\n"},
          d_ :> Error["MkKiraConfig: bad denominator form: ", d]
        }]],
        If[NotFreeQ[dens, cut],
          {"    cut_propagators: [",
          Riffle[Range[Length[dens]] // Select[MatchQ[dens[[#]], den[_, _, cut]]&], ", "],
          "]\n"
          }
          ,
          {}
        ]
      }
      ,
      {basis, bases}]
  ];
]

`MkKiraJobsYaml[]`

Create Kira’s jobs file.

MkKiraJobsYaml[filename_, bids_List, topsectors_, mode_] :=
Module[{bid, sector, r, s}, 
  MaybeMkFile[filename,
    "jobs:\n",
    " - reduce_sectors:\n",
    "    reduce:\n",
    Table[
        r = Max[sector["r"], 1 + Plus@@sector["idx"]];
        s = Max[sector["s"], 1];
        {"     - {topologies: [", KiraBasisName[bid], "], sectors: [b", sector["idx"], "], r: ", r, ", s: ", s, "}\n"}
        ,
        {bid, bids},
        {sector, topsectors[bid]}],
    "    select_integrals:\n",
    "     select_mandatory_list:\n",
    Table[
        {"      - [", KiraBasisName[bid], ", \"", KiraBasisName[bid], ".integrals\"]\n"}
        ,
        {bid, bids}],
    "#     select_mandatory_recursively:\n",
    Table[
        {
        "#      - {topologies: [", KiraBasisName[bid],
            "], sectors: [b", sector["idx"],
            "], r: ", sector["r"],
            ", s: ", sector["s"],
            ", d: ", sector["d"], "}\n"},
        {bid, bids},
        {sector, topsectors[bid]}],
    "    integral_ordering: 8\n",
    "    preferred_masters: \"preferred-masters\"\n",
    Switch[mode,
      "reduce", {
        "    run_symmetries: true\n",
        "    run_initiate: true\n",
        "#    run_triangular: true\n",
        "#    run_back_substitution: true\n",
        "    run_firefly: true\n",
        " - kira2math:\n",
        "    target:\n",
        Table[
          {"     - [", KiraBasisName[bid], ", \"", KiraBasisName[bid], ".integrals\"]\n"},
          {bid, bids}],
        Table[
            {"#     - {topologies: [", KiraBasisName[bid], "], sectors: [b", sector["idx"], "], r: ", sector["r"], ", s: ", sector["s"], ", d: ", sector["d"], "}\n"},
            {bid, bids},
            {sector, topsectors[bid]}],
        "    reconstruct_mass: false\n"
      },
      "prepare", {
        "    run_symmetries: true\n",
        "    run_initiate: true\n",
        "    run_triangular: false\n",
        "    run_back_substitution: false\n",
        "    run_firefly: false\n"
      },
      "prepare-input", {
        "    run_symmetries: true\n",
        "    run_initiate: input\n",
        "    run_triangular: false\n",
        "    run_back_substitution: false\n",
        "    run_firefly: false\n"
      }
    ]
  ];
]

`MkKiraFinishUDSYaml[]`

Create Kira’s jobs file for UD system reduction.

MkKiraFinishUDSYaml[filename_, bids_List] :=
Module[{bid, sector, r, s}, 
  MaybeMkFile[filename,
    "jobs:\n",
    " - reduce_user_defined_system:\n",
    "    input_system:\n",
    "      files:\n",
    Table[{"       - \"input_kira/", KiraBasisName[bid], "\"\n"}, {bid, bids}],
    "      config: false\n",
    "      otf: true\n",
    "    select_integrals:\n",
    "     select_mandatory_list:\n",
    Table[{"      - [", KiraBasisName[bid], ", \"", KiraBasisName[bid], ".integrals\"]\n"}, {bid, bids}],
    "    integral_ordering: 8\n",
    "    preferred_masters: \"preferred-masters\"\n",
    "    run_symmetries: false\n",
    "    run_initiate: false\n",
    "    run_triangular: false\n",
    "    run_back_substitution: false\n",
    "    run_firefly: true\n",
    " - kira2math:\n",
    "    target:\n",
    Table[
      {"     - [", KiraBasisName[bid], ", \"", KiraBasisName[bid], ".integrals\"]\n"},
      {bid, bids}],
    "    reconstruct_mass: false\n"
  ];
]

`MkKiraIntegrals[]`

Create Kira’s integral list files, one per basis.

MkKiraIntegrals[dirname_, blist_] := Module[{bid, idxlist},
  Do[
    idxlist = blist // CaseUnion[B[bid, idx__] :> {idx}];
    MaybeMkFile[dirname <> "/" <> KiraBasisName[bid] <> ".integrals",
      idxlist // Map[{KiraBasisName[bid], "[", Riffle[#, ","], "]\n"}&]
    ];
    ,
    {bid, blist // CaseUnion[B[bid_, ___] :> bid]}];
]

`MkKiraIntegralList[]`

Create Kira’s integral list file.

MkKiraIntegralList[filename_, blist_] :=
  MaybeMkFile[filename,
    blist //
      CaseUnion[B[bid_, idx__] :> {bid, {idx}}] //
      Map[Apply[{KiraBasisName[#1], "[", Riffle[#2, ","], "]\n"}&]]
  ]

`IndicesToSectorId[]`

R = denominator power sum Dots = denominator dot count T = denominator count S = numerator power sum

IndicesToSectorId[idx_List] := Plus @@ Table[If[idx[[i]] > 0, 2^(i-1), 0], {i, Length[idx]}]
SectorIdToIndices[sector_Integer, ndens_Integer] := IntegerDigits[sector, 2, ndens] // Reverse
IndicesToR[idx_List] := idx // Cases[n_ /; n > 0 :> n] // Apply[Plus]
IndicesToDots[idx_List] := idx // Cases[n_ /; n>1 :> n-1] // Apply[Plus]
IndicesToT[idx_List] := idx // Count[n_ /; n > 0]
IndicesToS[idx_List] := idx // Cases[n_ /; n < 0 :> -n] // Apply[Plus]

BToR[B[bid_, idx___]] := IndicesToR[{idx}]
BToR[DimShift[b_B, _]] := BToR[b]
BToDots[B[bid_, idx___]] := IndicesToDots[{idx}]
BToDots[DimShift[b_B, _]] := BToDots[b]
BToT[B[bid_, idx___]] := IndicesToT[{idx}]
BToT[DimShift[b_B, _]] := BToT[b]
BToS[B[bid_, idx___]] := IndicesToS[{idx}]
BToS[DimShift[b_B, _]] := BToS[b]
BToBID[B[bid_, idx___]] := bid
BToBID[DimShift[b_B, _]] := BToBID[b]
BToIndices[B[bid_, idx___]] := {idx}
BToIndices[DimShift[b_B, n_]] := BToIndices[b]
BToDimension[_B] := 4
BToDimension[DimShift[_B, n_]] := 4+n

`BToSector[]`

Convert an integral in B notation into the corner integral of the corresponding sector (i.e. normalize all indices to 1 or 0.

BToSector[B[bid_, idx___]] :=
  B[bid, {idx} // MapReplace[{ n_ /; n <= 0 -> 0, _ -> 1}] // Apply[Sequence]]

BToSector[DimShift[b_B, _]] := BToSector[b]

`BToSectorPattern[]`

Convert an integral in B notation into a pattern matching the corresponding sector.

BToSectorPattern[B[bid_, idx___]] :=
  B[bid, {idx} // MapReplace[{ n_ /; n <= 0 -> _?NonPositive, _ -> _?Positive}] // Apply[Sequence]]

BToSectorPattern[DimShift[b_B, _]] := BToSectorPattern[b]

BToSectorPattern[bs_List] := bs // Map[BToSectorPattern] // Apply[Alternatives]

`TopSectors[]`

Figure out a list of topmost sectors containing all the supplied integrals. Each returned sector is an Association with keys id, idx, r, s, and d. The integrals are specified by the lists of their indices.

TopSectors[idxlist_List] :=
Module[{sector2i, sector2r, sector2s, sector2d, o2sectors, int, sector, r, s, d, o, sectors, done, i, ss},
  sector2i = <||>;
  sector2r = <||>;
  sector2s = <||>;
  sector2d = <||>;
  o2sectors = <||>;
  Do[
      sector = IndicesToSectorId[int];
      sector2i[sector] = SectorIdToIndices[sector, Length[int]];
      r = IndicesToR[int];
      s = IndicesToS[int];
      d = IndicesToDots[int];
      o = {s, r};
      o2sectors[o] = {o2sectors[o] /. _Missing -> {}, sector};
      sector2r[sector] = Max[r, sector2r[sector] /. _Missing -> 0];
      sector2d[sector] = Max[d, sector2d[sector] /. _Missing -> 0];
      (* Note: s=0 makes Kira produce false masters. It’s not
       * clear if we should only fix s=0 case, or if we need
       * to add +1 to all s. Currently we’re doing the former
       * inside [[MkKiraJobsYaml]].
       *)
      sector2s[sector] = Max[s, sector2s[sector] /. _Missing -> 0];
      ,
      {int, idxlist}
  ];
  sectors = {};
  done = {};
  Do[
      Do[
          If[MemberQ[done, sector], Continue[]];
          i = FirstPosition[done, ss_ /; (BitAnd[ss, sector] === sector)];
          If[MatchQ[i, _Missing],
              AppendTo[done, sector];
              AppendTo[sectors, sector];
              ,
              i = i[[1]];
              sector2r[done[[i]]] = Max[sector2r[sector], sector2r[done[[i]]]];
              sector2d[done[[i]]] = Max[sector2d[sector], sector2d[done[[i]]]];
              sector2s[done[[i]]] = Max[sector2s[sector], sector2s[done[[i]]]];
              ];
          ,
          {sector, o2sectors[o] // Flatten // Union // Reverse}
      ];
      ,
      {o, o2sectors // Keys // Sort // Reverse}
  ];
  Table[
    <|"id" -> sector, "idx" -> sector2i[sector], "r" -> sector2r[sector], "s" -> sector2s[sector], "d" -> sector2d[sector]|>
    ,
    {sector, sectors}] //
    Sort
]

`TopSectorMap[]`

Return an association from basis ids to top sector lists as given by TopSectors.

TopSectorMap[blist_] :=
  blist //
  GroupBy[Replace[B[bid_, ___] :> bid]] //
  Map[(#[[;;,2;;]])& /* Map[Apply[List]] /* TopSectors]

`SectorUnion[]`

Return a list of sectors that represent a union of the given ones. In other words, drop the subsectors.

SectorUnion[sectors_List] := Module[{},
  sectors //
    Union // 
    GroupBy[#[[1]]&] //
    Normal //
    MapReplace[
      (bid_ -> seclist_) :> (
        seclist[[;;, 2;;]] //
          Map[Apply[List]] //
          TopSectors //
          Map[#["idx"]&] //
          Map[B[bid, Sequence @@ #]&]
      )
    ] //
    Apply[Join]
]

`MkKiraEquations[]`

Create a list of preferred masters in Kira syntax. The masters themselves can be both single integrals in the B notation, or linear combinations of them.

MkKiraEquations[filename_String, masters:{(_B|_DimShift|_amplitude) ...}] :=
Module[{COEF},
  MaybeMkFile[filename,
    masters //
      MapReplace[{
        B[bid_, idx___] :>
          {KiraBasisName[bid], "[", Riffle[{idx}, ","], "]\n"},
        DimShift[B[bid_, idx___], n_] :>
          {KiraBasisName[bid], "_dimshift", n, "[", Riffle[{idx}, ","], "]\n"},
        amplitude[idx___] :>
          {"amplitude[", Riffle[{idx}, ","], "]\n"}
      }]
  ]
]

MkKiraEquations[filename_String, expressions_List] :=
Module[{COEF},
  MaybeMkFile[filename,
    expressions //
      Bracket[#, _B|_DimShift|_amplitude, Factor /* COEF]& //
      Map[
        Terms /*
        MapReplace[{
            B[bid_, idx___]*COEF[co_] :>
              {KiraBasisName[bid], "[", Riffle[{idx}, ","], "]*(", co // InputForm // ToString // StringReplace[" "->""], ")\n"},
            DimShift[B[bid_, idx___], n_]*COEF[co_] :>
              {KiraBasisName[bid], "_dimshift", n, "[", Riffle[{idx}, ","], "]*(", co // InputForm // ToString // StringReplace[" "->""], ")\n"},
            amplitude[idx___]*COEF[co_] :>
              {"amplitude[", Riffle[{idx}, ","], "]*(", co // InputForm // ToString // StringReplace[" "->""], ")\n"}
          }]
      ] //
        Riffle[#, "\n"]&
  ]
]

`MkKiraConfig[]`

Create Kira’s configuration directory for reduction of the given integrals under the given bases.

MkKiraConfig[dirname_, bases_List, blist_, OptionsPattern[]] :=
Module[{bid, bids, bid2topsector, idxlist, massdims},
  EnsureDirectory[dirname];
  EnsureDirectory[dirname <> "/config"];
  bids = blist // CaseUnion[B[bid_, ___] :> bid];
  bid2topsectors = Table[
    bid -> (blist // CaseUnion[B[bid, idx__] :> {idx}] // TopSectors // Sort)
    ,
    {bid, bids}] // Association;
  massdims = IBPBasisMassDimensions[bases];
  FailUnless[Sort[massdims[[;;,1]]] === Sort[bases[[1, "invariants"]]]];
  massdims = bases[[1, "invariants"]] // Map[(# -> (# /. massdims))&];
  MkKiraKinematicsYaml[dirname <> "/config/kinematics.yaml",
    bases[[1,"externalmom"]], bases[[1,"sprules"]], massdims, OptionValue[ReplaceByOne]];
  MkKiraIntegralFamiliesYaml[dirname <> "/config/integralfamilies.yaml",(* bases];*)
    bases // Select[MemberQ[bids, #["id"]]&]];
  MkKiraJobsYaml[dirname <> "/reduce.yaml", bids, bid2topsectors, "reduce"];
  MkKiraJobsYaml[dirname <> "/prepare.yaml", bids, bid2topsectors, "prepare"];
  MkKiraJobsYaml[dirname <> "/prepare-uds.yaml", bids, bid2topsectors, "prepare-input"];
  MkKiraFinishUDSYaml[dirname <> "/finish-uds.yaml", bids]
  MkKiraIntegrals[dirname, blist];
  MkKiraEquations[dirname <> "/preferred-masters", 
    OptionValue[PreferredMasters] // Select[NotFreeQ[B[Alternatives @@ bids, ___]]]];
  MaybeMkFile[dirname <> "/Makefile",
    "THREADS?= 1\n",
    "RUN ?=\n",
    "\n",
    "reduce.done: ",
      "reduce.yaml ",
      Table[{KiraBasisName[bid], ".integrals "}, {bid, bids}],
      "config/integralfamilies.yaml ",
      "config/kinematics.yaml ",
      "preferred-masters\n",
    "\t${RUN} tempwrap ",
      "./ -i 'b*integrals' -i 'pref*' -i '*yaml' ",
      "./ -o 'results/b*' -o '*.log' ",
      "-- $${KIRA:-kira} reduce.yaml --parallel=${THREADS} --bunch_size=8\n",
    "\tdate >$@\n",
    "\n",
    "prepare-uds.done: ",
      "prepare-uds.yaml ",
      Table[{KiraBasisName[bid], ".integrals "}, {bid, bids}],
      "config/integralfamilies.yaml ",
      "config/kinematics.yaml\n",
    "\trm -rf input_kira/\n",
    "\t${RUN} tempwrap ",
      "./ -i 'b*integrals' -i 'pref*' -i '*yaml' ",
      "./ -o 'input_kira/*' -o '*.log' ",
      "-- $${KIRA:-kira} prepare-uds.yaml --parallel=${THREADS} --bunch_size=8 >$@.log 2>&1\n",
    "\tdate >$@\n",
    "\n",
    "finish-uds.done: ",
      "sed-arguments ",
      "prepare-uds.done ",
      "finish-uds.yaml ",
      Table[{KiraBasisName[bid], ".integrals "}, {bid, bids}],
      "\n",
    "\t${RUN} tempwrap ",
      "./ -i 'b*integrals' -i 'input_kira/*' -i 'pref*' -i 'finish-uds.yaml' ",
      "./ -o 'results/b*' -o '*.log' ",
      "$$(cat sed-arguments) ",
      "-- $${KIRA:-kira} finish-uds.yaml --parallel=${THREADS} --bunch_size=8 >$@.log 2>&1\n",
    "\tdate >$@\n"
  ];
]
Options[MkKiraConfig] = {PreferredMasters -> {}, ReplaceByOne -> None};

KiraAmplitudeTag[n_, nprops_] := IntegerDigits[n, 2, nprops] // Reverse // Apply[amplitude]

`MkKiraConfigUDS[]`

Create Kira’s configuration directory for reduction of the given integrals under the given bases.

MkKiraConfigUDS[dirname_, bases_List, blist_, OptionsPattern[]] :=
Module[{bid, bids, bid2topsector, idxlist, massdims, sector, nprops, r, s},
  EnsureDirectory[dirname];
  EnsureDirectory[dirname <> "/config"];
  bids = blist // CaseUnion[B[bid_, ___] :> bid];
  bid2topsectors = Table[
    idxlist = blist // CaseUnion[B[bid, idx__] :> {idx}];
    bid -> (idxlist // TopSectors // Sort)
    ,
    {bid, bids}] // Association;
  massdims = IBPBasisMassDimensions[bases];
  FailUnless[Sort[massdims[[;;,1]]] === Sort[bases[[1, "invariants"]]]];
  massdims = bases[[1, "invariants"]] // Map[(# -> (# /. massdims))&];
  MkKiraKinematicsYaml[dirname <> "/config/kinematics.yaml",
    bases[[1,"externalmom"]], bases[[1,"sprules"]], massdims, OptionValue[ReplaceByOne]];
  MkKiraIntegralFamiliesYaml[dirname <> "/config/integralfamilies.yaml",(* bases];*)
    bases // Select[MemberQ[bids, #["id"]]&]];

jobs:

reduce_sectors: reduce:
- {r: 7, s: 2} weight_notation: false generate_input: { level: 0}

  MaybeMkFile[dirname <> "/prepare.yaml",
    "jobs:\n",
    " - reduce_sectors:\n",
    "    reduce:\n",
    Table[
        r = Max[sector["r"], 1 + Plus@@sector["idx"]];
        s = Max[sector["s"], 1];
        {"     - {topologies: [", KiraBasisName[bid], "], sectors: [b", sector["idx"], "], r: ", r, ", s: ", s, "}\n"}
        ,
        {bid, bids},
        {sector, bid2topsectors[bid]}],
    "    generate_input: {level: 0}\n",
    "    weight_notation: false\n"
  ];
  nprops = Length[bases[[1,"denominators"]]];
  MkKiraEquations[dirname <> "/expressions",
    blist // MapIndexed1[#1-KiraAmplitudeTag[#2, nprops]&]];
  MaybeMkFile[dirname <> "/weights",
    Table[{KiraAmplitudeTag[i, nprops] // ToString // StringReplace[" "->""], "\n"}, {i, Length[blist]}]];
  MaybeMkFile[dirname <> "/solve.yaml",
    "jobs:\n",
    " - reduce_user_defined_system:\n",
    "    input_system:\n",
    "     config: false\n",
    "     otf: true\n",
    "     files:\n",
    "      - \"expressions\"\n",
    Table[
      {"      - \"input_kira/", KiraBasisName[bid], "\"\n"},
      {bid, bids}],
    "    select_integrals:\n",
    "     select_mandatory_list:\n",
    "      - [\"weights\"]\n",
    "#    iterative_reduction: masterwise\n",
    "    preferred_masters: \"preferred-masters\"\n",
    "    run_initiate: true\n",
    "#    run_triangular: true\n",
    "#    run_back_substitution: true\n",
    "    run_firefly: true\n",
    " - kira2math:\n",
    "    target:\n",
    "      - [\"weights\"]\n"
  ];
  MkKiraEquations[dirname <> "/preferred-masters", 
    OptionValue[PreferredMasters] // Select[NotFreeQ[B[Alternatives @@ bids, ___]]]];
  MaybeMkFile[dirname <> "/Makefile",
    "THREADS?= 1\n",
    "RUN ?=\n",
    "\n",
    "prepare.done: ",
      "prepare.yaml ",
      "preferred-masters ",
      "config/integralfamilies.yaml ",
      "config/kinematics.yaml\n",
    "\trm -rf input_kira/ results/\n",
    "\t${RUN} tempwrap ",
      "./ -i '*yaml' -i 'pref*' ",
      "./ -o 'input_kira/*' -o '*.log' -o '*id2int' ",
      "-- $${KIRA:-kira} prepare.yaml --parallel=${THREADS} --bunch_size=8\n",
    "\tdate >$@\n",
    "\n",
    "solve.done: ",
      "sed-arguments ",
      "prepare.done ",
      "solve.yaml ",
      "weights ",
      "\n",
    "\t${RUN} tempwrap ",
      "./ -i '*yaml' -i expressions -i weights -i 'pref*' -i 'input_kira/*' ",
      "./ -o 'results/*' -o '*.log' ",
      "$$(cat sed-arguments) ",
      "-- $${KIRA:-kira} solve.yaml --parallel=${THREADS} --bunch_size=8\n",
    "\tdate >$@\n"
  ];
  MaybeMkFile[dirname <> "/sed-arguments", ""];
]
Options[MkKiraConfigUDS] = {PreferredMasters -> {}, ReplaceByOne -> None};

`MkKiraConfigByBasis[]`

Populate a directory with Kira subdirectories for each basis, and write a Makefile that runs Kira for each of them. With this done, just running make should already reduce each basis (separately).

The intended usage however is to parallelize the computations on a cluster. To this end, first define THREADS environment variable, that specifies how many threads should each Kira use, and the RUN environment variable, that specifies the command prefix for running jobs on the cluster (for example, srun -c16 --mem=50G --tmp=10G); then run make -j99.

The idea behind THREADS and RUN variables is the same as with SecDecPrepare.

Note that unlike the direct usage, to discover symmetries between master integrals in different bases, a separate combined IBP run is needed -- this time covering only the master integrals.

MkKiraConfigByBasis[dirname_, bases_List, blist_, OptionsPattern[]] :=
Module[{bid, bids, bid2basis, name},
  bids = blist // CaseUnion[B[bid_, ___] :> bid];
  bid2basis = bases // GroupBy[#["id"]&] // Map[Only];
  Do[
    MkKiraConfig[
      dirname <> "/" <> KiraBasisName[bid],
      {bid2basis[bid]},
      blist // CaseUnion[B[bid,___]],
      PreferredMasters -> OptionValue[PreferredMasters],
      ReplaceByOne -> OptionValue[ReplaceByOne]
    ];
    ,
    {bid, bids}];
  MaybeMkFile[dirname <> "/Makefile",
    "THREADS?= 1\n",
    "RUN ?=\n",
    Table[
      name = KiraBasisName[bid];
      {
        "\n",
        "reduce: ", name, "/reduce.done\n",
        "\n",
        name, "/reduce.done: ",
          name, "/reduce.yaml ",
          name, "/config/integralfamilies.yaml ",
          name, "/config/kinematics.yaml ",
          name, "/preferred-masters ",
          name, "/", name, ".integrals\n",
        "\t${RUN} tempwrap ",
          name, "/ -i 'b*integrals' -i 'pref*' -i '*yaml' ",
          name, "/ -o 'results/b*' -o '*.log' ",
          "-- $${KIRA:-kira} reduce.yaml --parallel=${THREADS} --bunch_size=8 >$@.log 2>&1\n",
        "\tdate >$@\n",
        "\n",
        "prepare-uds: ", name, "/prepare-uds.done\n",
        "\n",
        name, "/prepare-uds.done: ",
          name, "/prepare-uds.yaml ",
          name, "/config/integralfamilies.yaml ",
          name, "/config/kinematics.yaml ",
          name, "/preferred-masters ",
          name, "/", name, ".integrals\n",
        "\trm -rf input_kira/\n",
        "\t${RUN} tempwrap ",
          name, "/ -i 'b*integrals' -i 'pref*' -i '*yaml' ",
          name, "/ -o 'input_kira/*' -o '*.log' ",
          "-- $${KIRA:-kira} prepare-uds.yaml --parallel=${THREADS} --bunch_size=8 >$@.log 2>&1\n",
        "\tdate >$@\n",
        "\n",
        "finish-uds: ", name, "/finish-uds.done\n",
        "\n",
        name, "/finish-uds.done: ",
          "sed-arguments ",
          name, "/prepare-uds.done ",
          name, "/finish-uds.yaml ",
          name, "/preferred-masters ",
          name, "/", name, ".integrals\n",
        "\t${RUN} tempwrap ",
          name, "/ -i 'b*integrals' -i 'input_kira/*' -i 'pref*' -i 'finish-uds.yaml' ",
          name, "/ -o 'results/b*' -o '*.log' ",
          "$$(cat sed-arguments) ",
          "-- $${KIRA:-kira} finish-uds.yaml --parallel=${THREADS} --bunch_size=8 >$@.log 2>&1\n",
        "\tdate >$@\n"
      },
      {bid, bids}]
  ];
]
Options[MkKiraConfigByBasis] = {PreferredMasters -> {}, ReplaceByOne -> None};

`KiraApplyResults[]`

Read the IBP tables from a Kira directory, apply them to a given expression.

KiraApplyResults[ex_, confdir_String, bases_List] :=
Module[{exx, idx, bids, bid, ibpmapfiles, bvar, table, bmap},
  exx = ex;
  bids = exx // CaseUnion[B[bid_, __] :> bid];
  Do[
    ibpmapfiles = MkString[confdir, "/results/", KiraBasisName[bid], "/kira_", KiraBasisName[bid], ".integrals.m"] // FileNames;
    If[ibpmapfiles === {}, Continue[]];
    Print["* Loading IBP tables for basis ", bid];
    FailUnless[Length[ibpmapfiles] === 1];
    table = LoadKiraMap[First[ibpmapfiles]] // Association // TM;
    (*BMapLoadFromKira[bmap, First[ibpmapfiles]];
    (*table = RunThrough["sed 's/b0*\\([0-9]*\\)\\[/B[\\1,/g' '" <> First[ibpmapfiles] <> "'", ""] // TM;
    BMapLoad[bmap, table] // TM;*)
    Print["Masters: "];
    BMapMasters[bmap] // Map[Print["- ", #]&];
    *)
    Print["* Applying the IBP tables for basis ", bid, "; expression size ", exx//ByteCount//FormatBytes];
    (*
    table = {};
    exx = exx // BMapApply[bmap] // TM;
    BMapClear[bmap];
    *)
    exx = exx /. table // TM;
    table = None;
    ,
    {bid, bids}
  ];
  exx
]
KiraApplyResults[confdir_String, bases_List] := KiraApplyResults[#, confdir, bases]&

`KiraApplyResultsInBulk[]`

Read the IBP tables from a Kira directory, apply them to a given expression.

KiraApplyResultsInBulk[ex_, confdir_String, bases_List] :=
Module[{exx, idx, bids, bid, ibpmapfiles, bvar, table, bmap},
  exx = ex;
  bids = exx // CaseUnion[B[bid_, __] :> bid];
  Do[
    ibpmapfiles = MkString[confdir, "/results/", KiraBasisName[bid], "/kira_", KiraBasisName[bid], ".integrals.m"] // FileNames;
    If[ibpmapfiles === {}, Continue[]];
    Print["* Loading IBP tables for basis ", bid];
    FailUnless[Length[ibpmapfiles] === 1];
    table = LoadKiraMap[First[ibpmapfiles]] // TM;
    BMapLoad[bmap, table] // TM;
    table = None;
    ,
    {bid, bids}
  ];
  Print["* Applying the IBP tables"];
  exx = ex // BMapApply[bmap];
  BMapClear[bmap];
  exx
]
KiraApplyResultsInBulk[confdir_String, bases_List] := KiraApplyResultsInBulk[#, confdir, bases]&

`LoadKiraMap[]`

Load an IBP map from a Kira output file.

LoadKiraMap[filename_] := Module[{table},
  FailUnless[FileExistsQ[filename]];
  table = RunThrough["sed -E -e 's/b0*([0-9]+)_dimshift([0-9]+)\\[/B[DimShift[\\1,\\2],/g' -e 's/b0*([0-9]*)\\[/B[\\1,/g' '" <> filename <> "'", ""];
  FailUnless[MatchQ[table, _List]];
  table /. B[DimShift[bid_, n_], idx___] :> DimShift[B[bid, idx], n]
]
LoadKiraMap[filenames_List] := filenames // Map[LoadKiraMap] // Apply[Join]

`LoadKiraMasters[]`

Extract the list of master integrals from Kira results.

LoadKiraMasters[confdir_String] :=
  FileNames[confdir <> "/results/b*/masters"] //
  Map[
    ReadString /*
    (StringSplit[#, "\n"]&) /*
    Map[
      StringReplace["[" -> ","] /*
      StringReplace["]" ~~ ___ -> "]"] /*
      StringReplace["b" ~~ x:(DigitCharacter ...) ~~ "_dimshift" ~~ n:(DigitCharacter ...) :>
        "B[DimShift[" <> x <> "," <> n <>"]"] /*
      StringReplace["b" -> "B["] /*
      ToExpression /*
      ReplaceAll[B[DimShift[bid_, n_], idx___] :> DimShift[B[bid, idx], n]]
    ]
  ] //
  Apply[Join] //
  Union

`KiraIBP[]`

Perform full IBP reduction of an expression using Kira, replacing all the B[...] with linear combinations of master integrals.

Note that usually it’s not the best idea to run Kira from Mathematica. It is possible though.

KiraIBP[ex_, bases_List, OptionsPattern[]] := Module[{blist, confdir, result},
  confdir = MkTempDirectory["kira", ""];
  MkKiraConfig[confdir, bases, ex // CaseUnion[_B],
    PreferredMasters -> OptionValue[PreferredMasters],
    ReplaceByOne -> OptionValue[ReplaceByOne]];
  If[Run[MkString["make -C '", confdir, "'"]] // TM // # =!= 0&,
    EnsureNoDirectory[confdir];
    Error["Failed to run kira in " <> confdir];
  ];
  result = KiraApplyResults[ex, confdir, bases];
  EnsureNoDirectory[confdir];
  result
]
KiraIBP[bases_List, OptionsPattern[]] := KiraIBP[#, bases, 
    PreferredMasters -> OptionValue[PreferredMasters],
    ReplaceByOne -> OptionValue[ReplaceByOne]]&;
Options[KiraIBP] = {PreferredMasters -> {}, ReplaceByOne -> None};

`MkKiraMastersConfig[]`

Create the Kira configuration directory for KiraMasters.

MkKiraMastersConfig[dirname_String, bases_List, maxdots_Integer, maxnums_Integer, intordering_Integer] :=
Module[{idx},
  EnsureDirectory[dirname];
  EnsureDirectory[dirname <> "/config"];
  MkKiraKinematicsYaml[dirname <> "/config/kinematics.yaml",
    bases[[1,"externalmom"]],
    bases[[1,"sprules"]],
    IBPBasisMassDimensions[bases],
    bases[[1, "invariants", 1]]];
  MkKiraIntegralFamiliesYaml[dirname <> "/config/integralfamilies.yaml", bases];
  MaybeMkFile[dirname <> "/domasters.yaml",
    "jobs:\n",
    " - reduce_sectors:\n",
    "    reduce:\n",
    Table[
      idx = basis["denominators"] // MapReplace[den[_,_,irr] -> 0, _ -> 1];
      {"     - {topologies: [", KiraBasisName[basis["id"]], "], sectors: [b",
        idx, "], r: ", Count[idx,1] + maxdots, ", s: ", maxnums, "}\n"}
      ,
      {basis, bases}],
    "    select_integrals:\n",
    "     select_mandatory_recursively:\n",
    Table[
      idx = basis["denominators"] // MapReplace[den[_,_,irr] -> 0, _ -> 1];
      {"      - {topologies: [", KiraBasisName[basis["id"]], "], sectors: [b",
        idx, "], r: ", Count[idx,1] + maxdots, ", s: ", maxnums, ", d: ", maxdots, "}\n"}
      ,
      {basis, bases}],
    "    integral_ordering: ", intordering, "\n",
    "    run_symmetries: true\n",
    "    run_initiate: true\n",
    "    run_triangular: false\n",
    "    run_back_substitution: false\n",
    "    run_firefly: false\n"
  ];
]

`KiraMasters[]`

Run Kira to determine a set of master integrals in a given basis. The masters are selected according to Kira's integral ordering setting (an integer from 1 to 8; see Kira's documentation for what these mean).

KiraMasters[bases_List, maxdots_Integer, maxnums_Integer, OptionsPattern[]] :=
Module[{dirname, result},
  dirname = MkTemp["kira", ".masters"];
  MkKiraMastersConfig[dirname, bases, maxdots, maxnums, OptionValue[KiraIntegralOrdering]];
  SafeRun[MkString["cd '", dirname, "' && ", OptionValue[Kira], " domasters.yaml"]];
  result = LoadKiraMasters[dirname];
  EnsureNoDirectory[dirname];
  result
]
Options[KiraMasters] = {Kira -> "kira", KiraIntegralOrdering -> 8};

FIRE & LiteRed interface for IBP reduction

TrailingIrr[denominators_] := denominators /. {___, trail:Longest[den[_, _, irr] ...]} :> Length[{trail}]
LeadingIrr[denominators_] := denominators /. {lead:Longest[den[_, _, irr] ...], ___} :> Length[{lead}]

`PrepareFireStart[]`

Use FIRE and LiteRed to prepare the basis definition and symmetry files that FIRE will need for the reduction.

PrepareFireStart[basis_Association, confdir_String] := PrepareFireStart[basis, confdir, {}]
PrepareFireStart[basis_Association, confdir_String, invariantrules_] :=
Module[{tmpdir, props},
  Print["* Preparing FIRE basis ", basis["id"]];
  FIREPATH = Environment["FIREPATH"];
  FailUnless[FileExistsQ[FIREPATH <> "/FIRE6.m"]];
  props = basis["denominators"] //
    ReplaceAll[invariantrules] //
    ReplaceAll[{sp[p_] :> p^2, sp[p_,q_] :> p*q}] //
    MapReplace[den[p_] :> p^2, den[p_, m_, ___] :> p^2-m];
  tmpdir = MkTemp["fire", ""];
  EnsureCleanDirectory[tmpdir];
  Print["Directory: ", tmpdir];
  RunMathProgram[
    "Off[FrontEndObject::notavail];\n",
    "SetDirectory[\"", FIREPATH, "\"];\n",
    "Get[\"FIRE6.m\"];\n",
    "Internal = ", basis["loopmom"]//InputForm, ";\n",
    "External = ", basis["externalmom"]//InputForm, ";\n",
    "Propagators = ", props // InputForm, ";\n",
    "Replacements = ", basis["sprules"] /. sp[p_] :> p^2 /. sp[p_,q_] :> p*q /. invariantrules // InputForm, ";\n",
    "RESTRICTIONS = ",
      basis["denominators"] // MapIndexed[
        If[MatchQ[#1, den[_, _, cut]],
          Table[If[#2 === {idx}, -1, 0], {idx, Length[basis["denominators"]]}],
          Nothing
        ]&
      ]//InputForm,
    ";\n",
    "Print[\"Internal: \", Internal];\n",
    "Print[\"External: \", External];\n",
    "Print[\"Propagators: \", Propagators//InputForm];\n",
    "Print[\"Replacements: \", Replacements//InputForm];\n",
    "Print[\"RESTRICTIONS: \", RESTRICTIONS];\n",
    "Print[\"* PrepareIBP[]\"];\n",
    "PrepareIBP[];\n",
    "Print[\"* Prepare[]\"];\n",
    "Prepare[AutoDetectRestrictions->False];\n",
    "Print[\"* SaveStart[]\"];\n",
    "SaveStart[\"", tmpdir, "/start\"];\n",
    "Print[\"* Done with SaveStart[]\"];\n"
  ];
  RunMathProgram[
    "Off[FrontEndObject::notavail];\n",
    "Off[DiskSave::dir];\n",
    "Off[DiskSave::overwrite];\n",
    "SetDirectory[\"", FIREPATH, "/extra/LiteRed/Setup\"];\n",
    "Get[\"LiteRed.m\"];\n",
    "SetDirectory[\"", FIREPATH, "\"];\n",
    "Get[\"FIRE6.m\"];\n",
    "Internal = ", basis["loopmom"]//InputForm, ";\n",
    "External = ", basis["externalmom"]//InputForm, ";\n",
    "Propagators = ", props // InputForm, ";\n",
    "Replacements = ", basis["sprules"] /. sp[p_] :> p^2 /. sp[p_,q_] :> p*q /. invariantrules // InputForm, ";\n",
    "RESTRICTIONS = ",
      basis["denominators"] // MapIndexed[
        If[MatchQ[#1, den[_, _, cut]],
          Table[If[#2 === {idx}, -1, 0], {idx, Length[basis["denominators"]]}],
          Nothing
        ]&
      ]//InputForm,
    ";\n",
    "CreateNewBasis[basisx, Directory->(\"", tmpdir, "/litered\")];\n",
    "GenerateIBP[basisx];\n",
    "Print[\"* AnalyzeSectors[]\"];\n",
    "AnalyzeSectors[basisx,\n",
    "  ", basis["denominators"] // MapReplace[den[_, _, irr] -> 0, den[___] -> _] // InputForm, ",\n",
    "  CutDs -> (", basis["denominators"] // MapReplace[den[_, _, cut] -> 1, den[___] -> 0] // InputForm, ")];\n",
    "Print[\"* FindSymmetries[]\"];\n",
    "FindSymmetries[basisx];\n",
    "Print[\"* DiskSave[]\"];\n",
    "DiskSave[basisx];\n",
    "Print[\"* Done with DiskSave[]\"];\n"
  ];
  EnsureDirectory[confdir];
  RunMathProgram[
    "Off[FrontEndObject::notavail];\n",
    "SetDirectory[\"", FIREPATH, "\"];\n",
    "Get[\"FIRE6.m\"];\n",
    "LoadStart[\"", tmpdir, "/start\"];\n",
    "TransformRules[\"", tmpdir, "/litered\", \"", ExpandFileName[confdir], "/b", basis["id"], ".lbases\", ", basis["id"], "];\n",
    "SaveSBases[\"", ExpandFileName[confdir], "/b", basis["id"], "\"];\n",
    "Print[\"* Done with SaveSBases[]\"];\n"
  ];
  EnsureNoDirectory[tmpdir];
  MkFile[confdir <> "/b" <> ToString[basis["id"]] <> ".pos",
    "|", LeadingIrr[basis["denominators"]] + 1, ",", Length[basis["denominators"]] - TrailingIrr[basis["denominators"]], "|"
  ];
  Print["* Done with everything"];
]

`LoadFireMasters[]`

Get the list of master integrals from a FIRE .tables file.

LoadFireMasters[filename_String] := filename // SafeGet // #[[2, ;;, 2]]& // MapReplace[{bid_, idx_List} :> B[bid, Sequence @@ idx]]

`LoadFireTables[]`

Load IBP tables from FIRE .tables file.

LoadFireTables[filename_, coeff_: Identity, JoinTerms_: True] := Module[{temp, GGG, data},
    data = SafeGet[filename];
    temp = {GGG[##[[1]]], {GGG[##[[1]]], ##[[2]]} & /@ ##[[2]]} & /@ data[[1]];
    Set[GGG[##[[1]]], G[##[[2, 1]], ##[[2, 2]]]] & /@ data[[2]];
    temp = temp;
    Clear[GGG];
    temp = DeleteCases[temp, {a_, {{a_, "1"}}}];
    temp = {##[[1]], {##[[1]], ToExpression[##[[2]]]} & /@ ##[[2]]} & /@ temp;
    temp = {##[[1]], {##[[1]], coeff[##[[2]]]} & /@ ##[[2]]} & /@ temp;
    If[JoinTerms,
        temp = {##[[1]], Times @@@ ##[[2]]} & /@ temp;
        temp = {##[[1]], Plus @@ ##[[2]]} & /@ temp;
    ];
    Rule @@@ temp // ReplaceAll[G[bid_, idx_List] :> B[bid, Sequence @@ idx]]
 ]

pySecDec interface for numerical evaluation of integrals

`SecDecIntegralName[]`

Convert an integral name (B notation) into a filename.

SecDecIntegralName[integral_B] := integral //
  ToString //
  StringReplace[" " -> ""] //
  StringReplace["," -> "_"] //
  StringReplace["[" -> ""] //
  StringReplace["]" -> ""] //
  StringReplace["-" -> "m"]
SecDecIntegralName[DimShift[integral_B, n_]] :=
  SecDecIntegralName[integral] <> "_dshift" <> If[Negative[n], "m", ""] <> ToString[Abs[n]]

`ToSympy[]`

Convert an expression into a string in Sympy format.

This might need additional work to correctly convert more forms.

ToSympy[ex_] := ex /. Pi -> pi //
  InputForm //
  ToString //
  StringReplace[#, WordBoundary ~~ "Sqrt[" -> "sqrt["]& //
  StringReplace[#, "[" -> "("]& //
  StringReplace[#, "]" -> ")"]& //
  StringReplace[#, "*^" -> "e"]&

`SecDecPrepare[]`

Prepare pySecDec files in a given directory for the given list of integrals. These can then be compiled manually by running make compile.

Each integral in the list can be:

B[...] * prefactor;
DimShift[B[...], n] * prefactor;
{integral, epsilon-order}.

One can alternatively use SecDecCompile to both prepare and compile.

The resulting directory consists of a Makefile and a set of python scripts. While compilation can be performed by just make compile, there are two provisions to parallelize the compilation: the THREADS environment variable which sets how many threads should comilation of a single integral use, and the RUN variable that is inserted as a prefix before the slow generation and compilation commands. The intent is for you to be able to set RUN to the command that executes its arguments on a cluster (e.g. srun -c 10 --mem=10G), and THREADS to the number of processes to use on the remote machine (e.g. 10). With these two set, make -j99 compile will parallelize the compilation across 99 cluster machines.

Of course, make -j99 compile will also work well on a local machine, if it has 99 processors and lots of free memory.

SecDecPrepare[basedir_String, bases_List, integrals_List] :=
Module[{name, basisid, bid2basis, indices, basis, integral, p, m, dim, order, prefactor},
  bid2basis = bases // GroupBy[#["id"]&] // Map[Only];
  EnsureDirectory[basedir];
  Do[
    {name, integral, dim, order, prefactor} = integral // Replace[{
        pre_. * b_B :> {SecDecIntegralName[b], b, 4, 0, pre},
        pre_. * int:DimShift[b_B, n_] :> {SecDecIntegralName[int], b, 4 + n, 0, pre},
        {pre_. * b_B, ord_} :> {SecDecIntegralName[b], b, 4, ord, pre},
        {pre_. * int:DimShift[b_B, n_], ord_} :> {SecDecIntegralName[int], b, 4+n, ord, pre}
    }];
    Print["* Making ", basedir, "/", name, ".*"];
    basisid = integral[[1]];
    indices = integral[[2;;]] // Apply[List];
    basis = bid2basis[basisid];
    MaybeMkFile[basedir <> "/" <> name <> ".compile.py",
      "#!/usr/bin/env python3\n",
      "import os\n",
      "import subprocess\n",
      "import tempfile\n",
      "import pySecDec as psd\n",
      "loopint = psd.loop_integral.LoopIntegralFromPropagators(\n",
      "  loop_momenta = ['", basis["loopmom"] // Riffle[#, "','"]&, "'],\n",
      "  external_momenta = ['", basis["externalmom"] // Riffle[#, "','"]&, "'],\n",
      "  regulator = 'eps',\n",
      "  propagators = [\n",
      basis["denominators"] /. {
        den[p_] :> {"    '(", p//ToSympy, ")^2'"},
        den[p_,m_,___] :> {"    '(", p//ToSympy, ")^2-", m//ToSympy, "'"}
      } // Riffle[#, ",\n"]&,
      "\n",
      "  ],\n",
      "  powerlist = [", indices // Riffle[#, ","]&, "],\n",
      "  dimensionality='", dim, "-2*eps',\n",
      "  replacement_rules = [\n  ",
      basis["sprules"] //
        ReplaceAll[sp -> (sp /* Sort)] //
        Union //
        MapReplace[
          (sp[p1_, p2_] -> v_) :> {"   ('", p1//InputForm, "*", p2//InputForm, "', '", v//InputForm, "')"}
        ] // Riffle[#, ",\n  "]&,
      "\n  ]\n",
      ")\n",
      "if __name__ == '__main__':\n",
      "    subprocess.check_call(['rm', '-rf', '", name, "', '", name, "_data'])\n",
      "    threads = int(os.environ.get('THREADS', '1'))\n",
      "    cwd = os.getcwd()\n",
      "    with tempfile.TemporaryDirectory(prefix='psd') as tmp:\n",
      "        os.chdir(tmp)\n",
      "        psd.loop_integral.loop_package(\n",
      "            name = '", name, "',\n",
      "            loop_integral = loopint,\n",
      "            real_parameters = [",
        basis["invariants"] // Map[{"'", #, "'"}&] // Riffle[#, ", "]&,
      "],\n",
      "            additional_prefactor = '",
        (* Convert from pySecDec to the physical normalization. *)
        (I*Pi^(dim/2-eps)/(2*Pi)^(dim-2*eps))^Length[basis["loopmom"]] * prefactor // Factor // ToSympy,
      "',\n",
      "            decomposition_method = '", If[Count[indices, Except[0]] > 2, "geometric_ku", "iterative"],"',\n",
      "            use_dreadnaut = True,\n",
      "            requested_order = ", order, ",\n",
      "            processes = threads,\n",
      "            form_optimization_level = 4,\n",
      "            form_work_space = '100M',\n",
      "            form_threads = 1,\n",
      "            contour_deformation = True\n",
      "        )\n",
      "        subprocess.check_call(['make', '-C', '", name, "', '-j', str(threads), '", name, "_pylink.so'])\n",
      "        subprocess.check_call(['mv', '", name, "/", name, "_pylink.so', cwd])\n",
      "        subprocess.check_call(['mv', '", name, "/", name, "_data', cwd])\n",
      "        subprocess.check_call(['rm', '-rf', '", name, "'])\n"
    ];
    MaybeMkFile[basedir <> "/" <> name <> ".integrate.py",
      "#!/usr/bin/env python3\n",
      "import os\n",
      "import sys\n",
      "import pySecDec as psd\n",
      "if __name__ == '__main__':\n",
      "    argmap = dict(arg.split('=') for arg in sys.argv[1:])\n",
      "    parameters = [",
        basis["invariants"] // Map[{"float(argmap['", #, "'])"}&] // Riffle[#, ", "]&,
      "]\n",
      "    threads = int(os.environ.get('THREADS', '1'))\n",
      "    libfile = './", name, "_pylink.so'\n",
      "    lib = psd.integral_interface.IntegralLibrary(libfile)\n",
      "    lib.use_Qmc(transform='korobov3', cputhreads=threads, verbosity=1)\n",
      "    int_wo_prefactor, prefactor, int_with_prefactor = lib(real_parameters=parameters, verbose=True, epsrel=1e-6, epsabs=1e-18, wall_clock_limit=2*60*60)\n",
      "    print(int_with_prefactor, file=sys.stderr)\n",
      "    result = '{%s,\\n%s}' % psd.integral_interface.series_to_mathematica(int_with_prefactor)\n",
      "    print(result)\n"
    ];
    ,
    {integral, integrals}];
  Print["* Making ", basedir, "/Makefile"];
  MaybeMkFile[basedir <> "/Makefile",
    "THREADS ?= 1\n",
    "RUN ?=\n",
    Table[
      name = integral // Replace[{
        pre_. * b_B :> SecDecIntegralName[b],
        pre_. * int:DimShift[b_B, n_] :> SecDecIntegralName[int],
        {pre_. * b_B, ord_} :> SecDecIntegralName[b],
        {pre_. * int:DimShift[b_B, n_], ord_} :> SecDecIntegralName[int]
      }];
      {
        "\n",
        "compile: ", name , "_pylink.so\n",
        "integrate: ", name , "_value.m\n",
        "\n",
        name, "_pylink.so: ", name, ".compile.py\n",
        "\t${RUN} python3 ", name, ".compile.py\n",
        "\n",
        (* This part is tricky: because `$RUN` might send the
        * command to a cluster, and there is some delay in the
        * propagation of data across the network filesystem,
        * we can not allow `*.integrate.py` scripts to directly
        * read `invariants.txt` or write to `*_value.m` files. So
        * we supply the parameters as command line arguments,
        * and expect the results on stdout. Also to make sure
        * that if the script fails, `*_value.m` file is removed
        * (or at least marked as outdated), so it would be rebuilt
        * the next time `make` is executed.
        *)
        name, "_value.m: ", name, "_pylink.so ", name, ".integrate.py invariants.txt\n",
        "\ttouch -t 199001010000 $@\n",
        "\t${RUN} python3 ", name, ".integrate.py $$(cat invariants.txt) >$@ 2>$@.log || rm -f $@\n"
      }
      ,
      {integral, integrals}]
  ];
]

ToPythonString[text_String] := StringJoin[
  "\"",
  StringReplace[text, {
    "\n" -> "\\n",
    "\r" -> "\\r",
    "\t" -> "\\t",
    "\\" -> "\\\\",
    "\"" -> "\\\""
  }],
  "\""]

`SecDecPrepareSum[]`

Prepare pySecDec files in a given directory for the given weighted sums of integrals (in B[...] or DimShift[B[...],...] syntax). The sums must be specified as Association maps from string (sum names) into expressions.

The resulting directory consists of a single Python script compile.py. When executed, this script will compile a pySecDec disteval integration library, and put it into the disteval subfolder, with disteval/sum.json being the main entry point for the evaluation.

By default compile.py will use all the locally available CPUs; this can be overridden by setting the THREADS environment variable.

After the compilation one can use the disteval CLI, or [SecDecIntegrateSum] to use the integration library to evaluate the specified sums.

Note that the loop integration in the integrals is assumed to come with physical normalization, i.e. $1/(2\pi)^d$ per loop, which is different from the default pySecDec normalization ( $1/(i\pi^{d/2})$ per loop). Additional prefactors can be passed via the IntegralPrefactors option as an Association mapping from integrals to prefactors.

There is a provision for automatically detecting Euclidean integrals, and compiling them without contour deformaion, thus making the evaluation faster: the InvariantSignMap option. If certain parameters of the integrals are known to always be positive (or negative), add an entry to InvariantSignMap mapping the parameter to +1 or -1. This information will then be used to disable contour deformation for integrals with polynomials that are always positive.

SecDecPrepareSum[basedir_String, bases_List, sums_Association, OptionsPattern[]] :=
Module[{sum2int2co, integrals, int2idx, name, basisid, bid2basis, indices, basis, integral, p, m, dim, order, prefactor, cache, intfunnames, sumnames},
  sum2int2co = sums // Map[BracketAssociation[_DimShift|_B]];
  integrals = sum2int2co // Map[Keys] // Values // Apply[Join] // Union;
  int2idx = integrals // PositionIndex // Map[First];
  FailUnless[Length[sums] > 0];
  bid2basis = bases // GroupBy[#["id"]&] // Map[Only];
  EnsureDirectory[basedir];
  cache = <||>;
  intfunnames = integrals // MapIndexed1[MkString["int_", #2-1, "_", SecDecIntegralName[#1]]&];
  MaybeMkFile[basedir <> "/compile.py",
    "#!/usr/bin/env python3\n",
    "import multiprocessing\n",
    "import os\n",
    "import subprocess\n",
    "import sys\n",
    "import tempfile\n",
    "import sympy as sp\n",
    "import pySecDec as psd\n",
    "\n",
    "signmap = {\n",
    Table[{"  sp.Symbol('x", i-1, "'): 1,\n"}, {i, bases // Map[#["denominators"]& /* Length] // Max}],
    If[OptionValue[InvariantSignMap] === None,
      "",
      OptionValue[InvariantSignMap] // Normal // MapReplace[(k_ -> v_) :> {"  sp.Symbol('", k, "'): ", v // InputForm, ",\n"}]
    ],
    "}\n",
    "\n",
    "def needs_cdeform(name, Fpoly):\n",
    "  poly = sp.sympify(str(Fpoly))\n",
    "  terms = sp.Add.make_args(poly.expand())\n",
    "  res = not all((t.subs(signmap) > 0) == True for t in terms)\n",
    "  print(f'Needs-Contour-Deformation({name}, {poly}) = {res}')\n",
    "  return res\n",
    "\n",
    "def call(f):\n",
    "  return f()\n",
    Table[
      {
      "\n",
      "b", basis["id"], "_propagators = [\n",
      basis["denominators"] /. {
        den[p_] :> {"  '(", p//ToSympy, ")^2'"},
        den[p_,m_,___] :> {"  '(", p//ToSympy, ")^2-", m//ToSympy, "'"}
      } // Riffle[#, ",\n"]&,
      "\n]\n",
      "\n",
      "b", basis["id"], "_replacement_rules = ",
      If[KeyExistsQ[cache, basis["sprules"]],
        {cache[basis["sprules"]], "\n"}
        ,
        cache[basis["sprules"]] = MkString["b", basis["id"], "_replacement_rules"];
        {
        "[\n",
        basis["sprules"] //
          ReplaceAll[sp -> (sp /* Sort)] //
          Union //
          MapReplace[
            (sp[p1_, p2_] -> v_) :> {"  ('", p1//InputForm, "*", p2//InputForm, "', '", v//ToSympy, "')"}
          ] // Riffle[#, ",\n"]&,
        "\n]\n"
        }
      ]
      }
      ,
      {basis, bases}],
    Table[
      integral = integrals[[idx]];
      prefactor = Lookup[OptionValue[IntegralPrefactors], integral, 1];
      indices = BToIndices[integral];
      basis = bid2basis[BToBID[integral]];
      dim = BToDimension[integral];
      {
      "\n",
      "def ", intfunnames[[idx]], "():\n",
      "  li = psd.loop_integral.LoopIntegralFromPropagators(\n",
      "    loop_momenta = ['", basis["loopmom"] // Riffle[#, "','"]&, "'],\n",
      "    external_momenta = ['", basis["externalmom"] // Riffle[#, "','"]&, "'],\n",
      "    propagators = b", basis["id"], "_propagators,\n",
      "    powerlist = [", indices // Riffle[#, ","]&, "],\n",
      "    replacement_rules = b", basis["id"], "_replacement_rules,\n",
      "    dimensionality = '", dim, "-2*eps',\n",
      "    regulators = ['eps']\n",
      "  )\n",
      "  return psd.LoopPackage(\n",
      "    name = '", SecDecIntegralName[integral], "',\n",
      "    loop_integral = li,\n",
      "    real_parameters = [",
        basis["invariants"] // Map[{"'", #, "'"}&] // Riffle[#, ", "]&,
      "],\n",
      "    additional_prefactor = '",
        (* Convert from pySecDec to the physical normalization. *)
        (I*Pi^(dim/2-eps)/(2*Pi)^(dim-2*eps))^Length[basis["loopmom"]] * prefactor // Factor // ToSympy,
      "',\n",
      "    decomposition_method = '", If[Count[indices, Except[0]] > 2, "geometric_ku", "iterative"],"',\n",
      "    form_optimization_level = 4,\n",
      "    form_work_space = '100M',\n",
      "    form_threads = 1,\n",
      "    contour_deformation = needs_cdeform('", SecDecIntegralName[integral], "', li.F)\n",
      "  )\n"
      }
      ,
      {idx, Length[integrals]}
    ],
    "\n",
    "if __name__ == '__main__':\n",
    "  # Always start in the directory this file resides in\n",
    "  thisdir = os.path.dirname(sys.argv[0])\n",
    "  if thisdir: os.chdir(os.path.dirname(sys.argv[0]))\n",
    "\n",
    "  make_integrals = [\n",
    Table[{"    ", name, ",\n"}, {name, intfunnames}],
    "  ]\n",
    "\n",
    "  coefficients = {\n",
    sum2int2co //
      Normal //
      MapReplace[(sumname_ -> coefficients_) :> {
        "    ", sumname // ToPythonString, ": {\n",
        coefficients //
          Normal //
          MapReplace[(int_ -> coefficient_) :> {
            "      ", int2idx[int]-1, ": ", coefficient // ToSympy // ToPythonString, ",\n"
          }],
        "    },\n"
      }],
    "  }\n",
    "  try:\n",
    "    nthreads = int(os.environ['THREADS'])\n",
    "  except KeyError:\n",
    "    try:\n",
    "      nthreads = len(os.sched_getaffinity(0))\n",
    "    except AttributeError:\n",
    "      nthreads = os.cpu_count()\n",
    "  if nthreads > 1:\n",
    "    with multiprocessing.Pool(nthreads) as pool:\n",
    "      integrals = pool.map(call, make_integrals)\n",
    "  else:\n",
    "    integrals = [f() for f in make_integrals]\n",
    "  subprocess.check_call(['rm', '-rf', 'disteval'])\n",
    "  cwd = os.getcwd()\n",
    "  with tempfile.TemporaryDirectory(prefix='psd') as tmp:\n",
    "    os.chdir(tmp)\n",
    "    psd.sum_package('sum',\n",
    "      integrals,\n",
    "      coefficients = coefficients,\n",
    "      regulators = ['eps'],\n",
    "      requested_orders = [", OptionValue[Order], "],\n",
    "      real_parameters = [",
      basis["invariants"] // Map[{"'", #, "'"}&] // Riffle[#, ", "]&,
    "],\n",
    "      processes = nthreads,\n",
    "    )\n",
    "    subprocess.check_call(['make', '-C', 'sum', '-j', str(nthreads), 'disteval.done'])\n",
    "    subprocess.check_call(['cp', '-a', 'sum/disteval', cwd])\n",
    "    subprocess.check_call(['rm', '-rf', 'sum'])\n"
  ];
]
Options[SecDecPrepareSum] = {Order -> 0, IntegralPrefactors -> <||>, InvariantSignMap -> None};

`SecDecIntegrateSum[]`

Evaluate the sums contained in a given pySecDec integration library. The library is assume to have been prepared by SecDecPrepareSum, and compiled.

SecDecIntegrateSum[path_String, invariantmap_List, OptionsPattern[]] :=
Module[{tmp, result, regulators},
  tmp = MkTemp["pysecdec", ".output.json"];
  SafeRun["echo python3 -m pySecDec.disteval '", path, "' "
    "--epsrel=", OptionValue[EpsRel] // N // CForm, " ",
    "--epsabs=", OptionValue[EpsAbs] // N // CForm, " ",
    "--format=json ",
    "--lattice-candidates=11 ",
    invariantmap // Sort // MapReplace[(k_ -> v_) :> {k, "=", v//N//CForm}] // Riffle[#, " "]&, " ",
    ">'", tmp, "' || rm -f '", tmp, "'"
  ];
  SafeRun["python3 -m pySecDec.disteval '", path, "' "
    "--epsrel=", OptionValue[EpsRel] // N // CForm, " ",
    "--epsabs=", OptionValue[EpsAbs] // N // CForm, " ",
    "--format=json ",
    "--lattice-candidates=11 ",
    invariantmap // Sort // MapReplace[(k_ -> v_) :> {k, "=", v//N//CForm}] // Riffle[#, " "]&, " ",
    ">'", tmp, "' || rm -f '", tmp, "'"
  ];
  result = Import[tmp];
  If[tmp === $Failed, Error["Could not open ", tmp, "; pySecDec has failed"]];
  DeleteFile[tmp];
  result = result // Association;
  regulators = result["regulators"] // Map[ToExpression];
  <|
    "sums" -> (
      result["sums"] // Association // Map[MapReplace[
        {regpow_List, {re_, im_}, {dre_, dim_}} :>
        {regulators^regpow // Apply[Times], Complex[re, im], Complex[dre, dim]}
      ]]
    ),
    "integrals" -> (
      result["integrals"] // Association // Map[MapReplace[
        {regpow_List, {re_, im_}, {dre_, dim_}} :>
        {regulators^regpow // Apply[Times], Complex[re, im], Complex[dre, dim]}
      ]]
    )
  |>
]
Options[SecDecIntegrateSum] = {
    EpsRel -> 10.^-3,
    EpsAbs -> 10.^-8
};

`SecDecLeadingOrders[]`

Return the leading expansion orders in 'eps' of a list of integrals computed via pySecDec.

SecDecLeadingOrders[base_List, integrals_List, OptionsPattern[]] :=
Module[{bid2basis, tmpscript, tmpresult},
  bid2basis = bases // GroupBy[#["id"]&] // Map[Only];
  tmpscript = MkTemp["psdlo", ".py"];
  tmpresult = tmpscript <> ".m";
  MaybeMkFile[tmpscript,
    "#!/usr/bin/env python3\n",
    "import os\n",
    "import subprocess\n",
    "import tempfile\n",
    "import json\n",
    "import pySecDec as psd\n",
    "if __name__ == '__main__':\n",
    "    terms = []\n",
    Table[
      {name, integral, dim, order, prefactor} = integrals[[i]] // Replace[{
          pre_. * b_B :> {SecDecIntegralName[b], b, 4, 0, pre},
          pre_. * int:DimShift[b_B, n_] :> {SecDecIntegralName[int], b, 4 + n, 0, pre},
          {pre_. * b_B, ord_} :> {SecDecIntegralName[b], b, 4, ord, pre},
          {pre_. * int:DimShift[b_B, n_], ord_} :> {SecDecIntegralName[int], b, 4+n, ord, pre}
      }];
      basisid = integral[[1]];
      indices = integral[[2;;]] // Apply[List];
      basis = bid2basis[basisid];
      {
      "    loopint = psd.loop_integral.LoopIntegralFromPropagators(\n",
      "      loop_momenta = ['", basis["loopmom"] // Riffle[#, "','"]&, "'],\n",
      "      external_momenta = ['", basis["externalmom"] // Riffle[#, "','"]&, "'],\n",
      "      regulator = 'eps',\n",
      "      propagators = [\n",
      basis["denominators"] /. {
        den[p_] :> {"        '(", p//CForm, ")^2'"},
        den[p_,m_,___] :> {"        '(", p//CForm, ")^2-", m//CForm, "'"}
      } // Riffle[#, ",\n"]&,
      "\n",
      "      ],\n",
      "      powerlist = [", indices // Riffle[#, ","]&, "],\n",
      "      dimensionality='", dim, "-2*eps',\n",
      "      replacement_rules = [\n  ",
      basis["sprules"] //
        ReplaceAll[sp -> (sp /* Sort)] //
        Union //
        MapReplace[
          (sp[p1_, p2_] -> v_) :> {"      ('", p1//InputForm, "*", p2//InputForm, "', '", v//InputForm, "')"}
        ] // Riffle[#, ",\n  "]&,
      "\n      ]\n",
      "    )\n",
      "    terms.append(psd.LoopPackage(\n",
      "      name = 'integral", i, "',\n",
      "      loop_integral = loopint,\n",
      "      real_parameters = [",
        basis["invariants"] // Map[{"'", #, "'"}&] // Riffle[#, ", "]&,
      "],\n",
      "      additional_prefactor = '",
        (* Convert from pySecDec to the physical normalization. *)
        (I*Pi^(dim/2-eps)/(2*Pi)^(dim-2*eps))^Length[basis["loopmom"]] * prefactor // Factor // ToSympy,
      "',\n",
      "      decomposition_method = '", If[Count[indices, Except[0]] > 2, "geometric_ku", "iterative"],"',\n",
      "      use_dreadnaut = True,\n",
      "      requested_order = ", order, ",\n",
      "      processes = 1,\n",
      "      form_optimization_level = 4,\n",
      "      form_work_space = '100M',\n",
      "      form_threads = 1,\n",
      "      contour_deformation = False\n",
      "    ))\n"
      }
      ,
      {i, Length[integrals]}],
    "    with tempfile.TemporaryDirectory(prefix='psd') as tmp:\n",
    "        os.chdir(tmp)\n",
    "        psd.sum_package(\n",
    "            'sum',\n",
    "            terms,\n",
    "            regulators=['eps'],\n",
    "            requested_orders=[", OptionValue[MaxOrder], "],\n",
    "            real_parameters=[",
      bases[[;;,"invariants"]] // Apply[Join] // Union // Map[{"'", #, "'"}&] // Riffle[#, ", "]&,
    "],\n",
    "            processes = ", OptionValue[Threads], "\n",
    "        )\n",
    "        orders = []\n",
    "        for i in range(", Length[integrals], "):\n",
    "            with open(f'sum/integral{i+1}/disteval/integral{i+1}.json', 'r') as f:\n",
    "                desk = json.load(f)\n",
    "            if desk['orders']:\n",
    "                orders.append(desk['lowest_orders'][0] + desk['prefactor_lowest_orders'][0])\n",
    "            else:\n",
    "                orders.append(", OptionValue[MaxOrder], ")\n",
    "        with open('", tmpresult, "', 'w') as f:\n",
    "            f.write('{' + ','.join(map(str, orders)) + '}')\n"
  ];
  If[Run["python3 '" <> tmpscript <> "' >/dev/null"] =!= 0,
    Error["pySecDec script failed\n"];
  ];
  results = Get[tmpresult];
  DeleteFile[{tmpscript, tmpresult}];
  results
]
Options[SecDecLeadingOrders] = {Threads -> 2, MaxOrder -> 1};

`MkSecDecCoefficientFile[]`

Write a single coefficient file in pySecDec syntax. It is expected that overall regulator factors are factorized.

MkSecDecCoefficientFile[filename_String, value_, regulators_List] :=
Module[{regpattern, num, den, reg},
  regpattern = regulators // Apply[Alternatives];
  nums = dens = regs = {};
  value //
    Factors //
    MapReplace[{
      x_Integer :> (nums = {nums, x}),
      x_Rational :> (nums = {nums, x//Numerator}; dens = {dens, x//Denominator}),
      x:(regpattern^n_.) :> (regs = {regs, x}),
      x:(_^n_ /; n > 0) :> (nums = {nums, x}),
      x:(_^n_ /; n < 0) :> (dens = {dens, 1/x}),
      x:(_) :> (nums = {nums, x})
    }];
  If[value === 0, regs = regulators // Map[#^500000000&]];
  MkFile[filename,
    "numerator = 1\n",
    nums // Flatten // MapReplace[{
      x_^n_ :> {"*(", x // InputForm, ")^(", n, ")\n"},
      x_ :> {"*(", x // InputForm, ")\n"}
    }],
    ";\n",
    "denominator = 1\n",
    dens // Flatten // MapReplace[{
      x_^n_ :> {"*(", x // InputForm, ")^(", n, ")\n"},
      x_ :> {"*(", x // InputForm, ")\n"}
    }],
    ";\n",
    "regulator_factor = 1\n",
    regs // Flatten // MapReplace[{
      x_Symbol^n_ :> {"*", x // InputForm, "^(", n, ")\n"},
      x_Symbol :> {"*", x // InputForm, "\n"},
      x_^n_ :> {"*(", x // InputForm, ")^(", n, ")\n"},
      x_ :> {"*(", x // InputForm, ")\n"}
    }],
    ";\n"
  ];
];

`MkSecDecCoefficientDir[]`

Populate a directory with coefficient files in pySecDec syntax.

MkSecDecCoefficientDir[coefdir_String, integrals_List, coefficients_List, regulators_List] :=
Module[{name, basisid, bid2basis, indices, basis, integral, coeff, p, m, dim, order, prefactor, c},
  FailUnless[Length[coefficients] > 0];
  FailUnless[Length[coefficients[[1]]] === Length[integrals]];
  EnsureDirectory[coefdir];
  Do[
    integral = integrals[[idx]];
    name = SecDecIntegralName[integral];
    Do[
      c = coefficients[[cidx, idx]];
      If[c =!= 0,
        MkSecDecCoefficientFile[MkString[coefdir, "/", name, "_coefficient", cidx-1, ".txt"], c, regulators];
      ];
      ,
      {cidx, Length[coefficients]}];
    ,
    {idx, Length[integrals]}];
]

`SecDecCompile[]`

Compile integration libraries for a list of integrals belonging to a given set of bases in a given directory.

Note that by design if this function is re-run with the same (or slightly different) parameters, then no (or very little) recompilation will take place.

SecDecCompile[basedir_String, bases_List, integrals_List, jobs_:1] := Module[{},
  SecDecPrepare[basedir, bases, integrals];
  SafeRun["make -j", jobs, " -C '", basedir, "' compile"];
]

MkInvariantsTxt[filename_, invariantmap_] :=
  MaybeMkFile[filename,
    invariantmap // Sort // MapReplace[(k_ -> v_) :> {k, "=", v//N//CForm}] // Riffle[#, " "]&
  ]

MkSedArguments[filename_, invmap_] := Module[{},
  MaybeMkFile[filename, invmap //
    MapReplace[{
      (x_Symbol -> value_Integer /; value >= 0) :> { "-e s,", x, ",", value // InputForm, ",g" },
      (x_Symbol -> value_) :> { "-e s,", x, ",(", value // InputForm, "),g" }
    }] //
    Riffle[#, " "]&,
    "\n"
  ];
]

`SecDecIntegrate[]`

Integrate a set of integrals at a given phase-space point given by a map of invariant values. The direcory should have already been prepared with SecDecPrepare, and optionally compiled with SecDecCompile.

Here the same setup as with SecDecPrepare is used: you can set the THREADS and RUN environment variables to parallelize the integration.

As with SecDecCompile this function can be re-run twice, and the second time it will return immediately because all the integration is already saved to disk. Only if the phase-space point has changed will the new integration be performed.

SecDecIntegrate[basedir_String, integrals_List, invariantmap_List, jobs_:1] :=
Module[{invstring,filenames},
  MkInvariantsTxt[basedir <> "/invariants.txt", invariantmap];
  filenames = integrals // Map[MkString[SecDecIntegralName[#], "_value.m"]&];
  SafeRun["make -j", jobs, " -C '", basedir, "' ", filenames // Riffle[#, " "]&];
  filenames // Map[SafeGet[basedir <> "/" <> #]&]
]
SecDecIntegrate[basedir_String, integral_, variables_List] :=
  SecDecIntegrate[basedir, {integral}, variables] // Only

Dimensional recurrence

DenToNumerator[den[p_]] := sp[p,p]
DenToNumerator[den[p_, m_, ___]] := sp[p,p] - m

ExpandSP[sp[p1_]] := ExpandSP[sp[p1,p1]]
ExpandSP[sp[p1_,p2_]] := p1*p2//Expand//Terms//MapReplace[
  a_Symbol*b_Symbol*c_. :> Sort[sp[a,b]]*c,
  a_Symbol^2*c_.:> sp[a,a]*c
]//Apply[Plus]
ExpandSP[ex_] := ex /. s_sp :> ExpandSP[s]

`RaisingDRR[]`

“Raising” dimensional recurrence: expresses a given integral in d-2 space-time dimensions as a linear combination of integrals in d space-time dimensions.

Note that Minkowski metrics is assumed. Multiply by (-1)^L to get Euclidean.

Also be careful about the normalization: the integrals are assumed to use loop integration measure of $d^d l/(i\pi^{d/2})$ , which is not the physical $d^d l/(2\pi)^d$ . To obtain relations between physically normalized integrals, multiply the result of this function by $(4\pi)^L$ .

Finally, you might want to cleanup the output of this function using ZeroSectors, because some of the integrals will come out to be scaleless.

RaisingDRR[ex_, basis_Association] := Module[{loopmom, i, j, k, mx, op, OP, bid, ii, n, idx, result},
  loopmom = basis["loopmom"];
  mx = Table[
    If[i===j, 1, 1/2]
    Sum[
      OP[k] D[
        basis["denominators"][[k]] // DenToNumerator // ExpandSP,
        Sort[sp[loopmom[[i]], loopmom[[j]]]]
      ],
      {k, Length[basis["denominators"]]}]
    ,
    {i, Length[loopmom]},
    {j, Length[loopmom]}];
  op = Det[mx] // Bracket[#, _OP, Factor]&;
  bid = basis["id"];
  result = op * ex // Bracket[#, _B|_OP, #&, ReplaceRepeated[#,
    OP[n_]^k_. B[bid, idx__] :> (ii = {idx}; ii[[n]] += k; Pochhammer[ii[[n]]-k,k] B[bid, Sequence@@ii])
  ]&]&;
  If[NotFreeQ[result, _OP], Error["Failed to replace all OPs"]];
  (-1)^Length[loopmom] result
]

`LoweringDRR[]`

“Lowering” dimensional recurrence: expresses a given integral in d+2 space-time dimensions as a linear combination of integrals in d space-time dimensions.

Note that Minkowski metrics is assumed. Multiply by (-1)^L to get Euclidean.

Same normalization issue as in RaisingDRR; divide by $(4\pi)^L$ to get the relation between physically normalized integrals.

LoweringDRR[ex_, basis_Association] := Module[{extmom, loopmom, op, OP, i, k, n, bid, ii, idx, result},
  extmom = basis["externalmom"];
  loopmom = basis["loopmom"];
  op = Det[GramMatrix[Join[loopmom, extmom]] /. basis["sprules"]] /.
    basis["nummap"] /.
    basis["sprules"] /.
    DEN[n_] :> 1/OP[n] //
    Bracket[#, _B, Together]&;
  bid = basis["id"];
  result = op * ex // Bracket[#, _B|_OP, #&, ReplaceRepeated[#,
    OP[n_]^k_. B[bid, idx__] :> (ii = {idx}; ii[[n]] -= k; B[bid, Sequence@@ii])
  ]&]&;
  If[NotFreeQ[result, _OP], Error["Failed to replace all As"]];
  (
    (-2)^Length[loopmom]
    1/Det[GramMatrix[extmom] /. basis["sprules"]]
    1/Pochhammer[d-Length[extmom]-Length[loopmom]+1, Length[loopmom]]
    result
  )
]

Integral differentiation

`BasisExternalInvariantSymbols[]`

All symbols on the right-hand side of scalar product rules. Presumably names of Mandelstam variables and the like.

BasisExternalInvariantSymbols[basis_Association] := basis[["sprules",;;,2]] // CaseUnion[_Symbol]

`BasisExternalScalarProducts[]`

All distinct sp[p1, p2], for p1 and p2 being external momenta of a basis.

BasisExternalScalarProducts[basis_Association] := Module[{p1, p2},
  splist = Table[
    Sort[sp[p1,p2]],
    {p1, basis["externalmom"]},
    {p2, basis["externalmom"]}
  ] // Apply[Join] // Union
]

`GramMatrix[]`

Gram matrix, |pi * pj|, for a given list of pi.

GramMatrix[vectors_List] := Module[{p,q}, Table[Sort[sp[p,q]], {p, vectors}, {q, vectors}]]

`BasisGramMatrix[]`

Gram matrix, |pi * pj|, for a given basis.

BasisGramMatrix[basis_Association] := GramMatrix[basis["externalmom"]] /. basis["sprules"] // Together

`BasisInvGramMatrix[]`

The inverse of Gram matrix for a given basis.

BasisInvGramMatrix[basis_Association] := Inverse[BasisGramMatrix[basis]] // Together

`BasisInvGramMatrixMap[]`

The inverse of Gram matrix for a given basis, represented as a nested Association. This is so it could be indexed by momenta as BasisInvGramMatrixMap[basis][p1,p2].

BasisInvGramMatrixMap[basis_Association] := BasisInvGramMatrixMap[basis] = Module[{igm, emom, i, j},
  igm = BasisInvGramMatrix[basis];
  emom = basis["externalmom"];
  Table[
    emom[[i]] -> Association[
      Table[emom[[j]] -> igm[[i,j]], {j, Length[emom]}]
    ]
    ,
    {i, Length[emom]}
  ] // Association
]

ClearAll[BDiffByMomentum, BDiffByMass, BDiffBySP, BDiffByInv, BDiff];
BDiffByMomentum[basis_Association, indices_List, p_Symbol, pmul_Symbol] := Module[{dens, ddens},
  dens = {basis["denominators"], indices} // Transpose // DeleteCases[{_, 0}];
  ddens = dens // MapReplace[
    {d:den[mom_, ___], n_} :> -n d^(n+1) 2 Sort[sp[mom, pmul]] D[mom, p]
  ];
  Sum[
    Product[If[i === k, ddens[[i]], dens[[i,1]]^dens[[i,2]]], {i, Length[dens]}]
    ,
    {k, Length[dens]}]
]
BDiffByMass[basis_Association, indices_List, mass_Symbol] := Module[{dens, ddens},
  dens = {basis["denominators"], indices} // Transpose // DeleteCases[{_, 0}];
  ddens = dens // MapReplace[
    {d:den[mom_], n_} :> -n d^(n+1) D[ExpandSP[sp[mom, mom]], mass]0,
    {d:den[mom_, m_, ___], n_} :> -n d^(n+1) D[ExpandSP[sp[mom, mom]]-m, mass]
  ];
  Sum[
    Product[If[i === k, ddens[[i]], dens[[i,1]]^dens[[i,2]]], {i, Length[dens]}]
    ,
    {k, Length[dens]}]//ToB[basis]
]
BDiffBySP[basis_Association, indices_List, s:sp[p1_Symbol, p1_Symbol]] := BDiffBySP[basis, indices, s] =
Module[{igm, pi},
  igm = BasisInvGramMatrixMap[basis];
  1/2 Sum[igm[pi,p1] BDiffByMomentum[basis, indices, p1, pi], {pi, basis["externalmom"]}]//ToB[basis]
]
BDiffBySP[basis_Association, indices_List, s:sp[p1_Symbol, p2_Symbol]] := BDiffBySP[basis, indices, s] =
Module[{igm, pi},
  igm = BasisInvGramMatrixMap[basis];
  Sum[igm[pi,p2] BDiffByMomentum[basis, indices, p1, pi], {pi, basis["externalmom"]}]//ToB[basis]
]
BDiffByInv[basis_Association, indices_List, inv_Symbol] := Module[{invlist, splist, ds, s},
  FailUnless[Length[basis["denominators"]] === Length[indices]];
  invlist = BasisExternalInvariantSymbols[basis];
  splist = BasisExternalScalarProducts[basis];
  (BDiffByMass[basis, indices, inv]) + Sum[
    ds = D[s/.basis["sprules"], inv];
    If[ds === 0, 0, ds*BDiffBySP[basis, indices, s]],
    {s, splist}]
]

`BDiff[]`

Diffferentiate an integral in the B notation by a scalar product of external momenta, or an invariant of the basis.

BDiff[ex_, basis_Association, s_sp] := ex /. B[basis["id"], idx__] :> BDiffBySP[basis, {idx}, s]
BDiff[ex_, basis_Association, inv_Symbol] := ex /. B[basis["id"], idx__] :> BDiffByInv[basis, {idx}, inv]
BDiff[basis_Association, inv_Symbol] := BDiff[#, basis, inv]&

Amplitudes

ClearAll[Amplitude];

`Amplitude[]`

Convert a diagram into an amplitude.

Only the default cases are here; all actual Feynman rules are supplied by the model files.

Amplitude[
  dia:Diagram[id_, factor_, ifields_List, ofields_List, propagators_List, vertices_List]] := Flatten[{
    factor,
    propagators // Map[Amplitude],
    vertices // Map[Amplitude]
  }] // Apply[Times]
Amplitude[P[field_, fi1_, fi2_, _, _, p_]] :=
  Error["No Feynman rules for the propagator of ", field]
Amplitude[V[_, fields_, ___]] :=
  Error["No Feynman rules for the vertex ", fields]
Amplitude[x__] := Error["Don't know an amplitude for: ", x]

ClearAll[CutAmplitudeGlue];

`CutAmplitudeGlue[]`

Take two diagrams with the same final states and return an amplitude factor (i.e. delta functions) that join the outgoing fields. Only the default cases are here; all actual Feynman rules are supplied by the model files.

CutAmplitudeGlue[f1_, f2_] := Error["Can't connect these cut lines: ", f1, " and ", f2]
CutAmplitudeGlue[
  dia1:Diagram[id1_, factor1_, ifields1_List, ofields1_List, propagators1_List, vertices1_List],
  dia2:Diagram[id2_, factor2_, ifields2_List, ofields2_List, propagators2_List, vertices2_List]
] := (
  FailUnless[ofields1[[;;,1]] === ofields2[[;;,1]]];
  MapThread[CutAmplitudeGlue, {ofields1, ofields2}] // Apply[Times]
)

`DiagramFinalStateSum[]`

Compute the sum over the final state particle states (i.e. the spin sum and the polarization sum) of a product of two diagrams, the second being complex-conjugated.

DiagramFinalStateSum[
  dia1:Diagram[id1_, factor1_, ifields1_List, ofields1_List, propagators1_List, vertices1_List],
  dia2:Diagram[id2_, factor2_, ifields2_List, ofields2_List, propagators2_List, vertices2_List]
] := (
  FailUnless[ofields1[[;;,1]] === ofields2[[;;,1]]];
  MapThread[CutAmplitudeGlue, {ofields1, ofields2}] // Apply[Times]
)

Amplitude library (alibrary.m)

Contents

Diagram generation

DiagramId[]

DiagramClosedLoops[]

Diagrams[]

Diagrams to graphs

DiagramGraphEdges[]

DiagramGraph[]

DiagramXGraph[]

DiagramToGraphviz[]

FieldGraphvizColor[]

DiagramViz[]

GraphLayoutVertexCoordinates[]

NormalizeCoordinates[]

TikZ interface for diagrams

DiagramToTikZ[]

FieldTikZStyle[]

MkTikZGraph[]

MkIntegralTikZ[]

IBP Bases

ExpandScalarProducts[]

BToDen[]

ToB[]

IBPRelations[]

VariableDimensions[]

NormalizeDens[]

CompleteIBPBasis[]

IBPBasisMapInvariants[]

InvariantMapUnderMomentaPermutation[]

IBPBasisMassDimensions[]

IBPBasisCross[]

IBPBasisSameQ[]

Feynson interface for integral symmetries

ZeroSectors[]

ZeroSectorPattern[]

SymmetryMaps[]

IntegralFamilyMappingRules[]

IntegralEqualitySets[]

UniqueIntegralIndices[]

IntegralUnion[]

IntegralMapOntoSectors[]

IntegralMapOntoFamily[]