(* ::Package:: *)

Quit


(* ::Subsection::Closed:: *)
(*Setup*)


(*Boundary condition, scaleless integrals vanish:*)
F[a1_,a2_,a3_,a4_,a5_] /; TrueQ[Sum[Sign[Sign[ai-1/2]+1],{ai,{a1,a2,a3,a4,a5}}]<3] := 0;
F[a1_,a2_,a3_,a4_,a5_] /; TrueQ[Sum[Sign[Sign[ai-1/2]+1],{ai,{a1,a2,a3,a4}}]<=3 && a5<=0] := 0;
F[a1_,a2_,a3_,a4_,a5_] /; TrueQ[(a1<=0&&a2<=0) || (a1<=0&&a3<=0) || (a2<=0&&a4<=0) || (a3<=0&&a4<=0)] := 0;
(*Symmetry conditions:*)
F[a1_,a2_,a3_,a4_,a5_] /; TrueQ[(Abs[a3]+Abs[a4] < Abs[a1]+Abs[a2]) || (Abs[a3]+Abs[a4] == Abs[a1]+Abs[a2] && a3<a1)] := F[a3,a4,a1,a2,a5];
F[a1_,a2_,a3_,a4_,a5_] /; TrueQ[(Abs[a2]+Abs[a4] < Abs[a1]+Abs[a3]) || (Abs[a2]+Abs[a4] == Abs[a1]+Abs[a3] && a2<a1)] := F[a2,a1,a4,a3,a5];


(*Verify relations between denominators and scalar products*)
{D1==S[l1,l1],D2==S[l1,l1]+2S[l1,q]+q2,D3==S[l1,l1]+2S[l1,l2]+S[l2,l2],D4==S[l1,l1]+2S[l1,l2]+2S[l1,q]+2S[l2,q]+S[l2,l2]+q2,D5==S[l2,l2]}/.{
S[l1,l1]->D1,S[l1,q]->(D2-D1-q2)/2,S[l1,l2]->(D3-D1-D5)/2,S[l2,l2]->D5,S[l2,q]->(D1-D2-D3+D4)/2,S[q,q]->q2}


(*Generate general IBPs*)
Attributes[S] = {Orderless};
IBP[l_,p_,{a1_,a2_,a3_,a4_,a5_}] := Simplify[Expand[Switch[l,
    (*Derivative wrt l1*)
    l1,If[TrueQ[p==l1],d*F[a1,a2,a3,a4,a5],0] - 
    a1*(2*S[l1,p])*F[a1+1,a2,a3,a4,a5] - 
    a2*(2*S[l1,p] + 2*S[q,p])*F[a1,a2+1,a3,a4,a5] - 
    a3*(2*S[l1,p] + 2*S[l2,p])*F[a1,a2,a3+1,a4,a5] - 
    a4*(2*S[l1,p] + 2*S[l2,p] + 2*S[q,p])*F[a1,a2,a3,a4+1,a5],
    (*Derivative wrt l2*)
    l2,If[TrueQ[p==l2],d*F[a1,a2,a3,a4,a5],0] - 
    a3*(2*S[l1,p] + 2*S[l2,p])*F[a1,a2,a3+1,a4,a5] - 
    a4*(2*S[l1,p] + 2*S[l2,p] + 2*S[q,p])*F[a1,a2,a3,a4+1,a5] - 
    a5*(2*S[l2,p])*F[a1,a2,a3,a4,a5+1]
    ] /. {
  S[l1,l1]->D1,S[l1,q]->(D2-D1-q2)/2,S[l1,l2]->(D3-D1-D5)/2,
  S[l2,l2]->D5,S[l2,q]->(D1-D2-D3+D4)/2,S[q,q]->q2}] /. {
  D1*F[i1_,i2_,i3_,i4_,i5_]->F[i1-1,i2,i3,i4,i5],
  D2*F[i1_,i2_,i3_,i4_,i5_]->F[i1,i2-1,i3,i4,i5],
  D3*F[i1_,i2_,i3_,i4_,i5_]->F[i1,i2,i3-1,i4,i5],
  D4*F[i1_,i2_,i3_,i4_,i5_]->F[i1,i2,i3,i4-1,i5],
  D5*F[i1_,i2_,i3_,i4_,i5_]->F[i1,i2,i3,i4,i5-1]}];

IBP[1,{a1_,a2_,a3_,a4_,a5_}] = IBP[l1,l1,{a1,a2,a3,a4,a5}];
IBP[2,{a1_,a2_,a3_,a4_,a5_}] = IBP[l1,l2,{a1,a2,a3,a4,a5}];
IBP[3,{a1_,a2_,a3_,a4_,a5_}] = IBP[l1,q,{a1,a2,a3,a4,a5}];
IBP[4,{a1_,a2_,a3_,a4_,a5_}] = IBP[l2,l1,{a1,a2,a3,a4,a5}];
IBP[5,{a1_,a2_,a3_,a4_,a5_}] = IBP[l2,l2,{a1,a2,a3,a4,a5}];
IBP[6,{a1_,a2_,a3_,a4_,a5_}] = IBP[l2,q,{a1,a2,a3,a4,a5}];


(* Process an individual IBP: If IBP is independent, solve for the least simple integral and add to known expressions *)
PrintInfo::usage="Print info about processes IBP.";
Options[ProcessIBP] = {PrintInfo ->True};

ProcessIBP[i_,{a1_,a2_,a3_,a4_,a5_},OptionsPattern[]] := Module[{ibp,ints,numdens,numpows,int,res},
  ibp = Collect[IBP[i,{a1,a2,a3,a4,a5}],F[__],Simplify];
  (*Stop if IBP is redundant*)
     If[TrueQ[ibp==0],
        If[TrueQ[OptionValue[PrintInfo]],Print["IBP is not independent."]];
        Return[Null]];
  (*Find the most complicated integral in the IBP*)
     (*unordered list of integrals that appear, need rescaling with a because otherwise integrals with coefficient 1 are missed*)
     ints = Union[Cases[ibp /. F[x__]->a*F[x], (x___ : 1) p__F :> List@@ p, {2}]];
     (*number of denominators for each integral*)
     numdens = Map[Total[Sign[Sign[Sequence[#]-{1/2,1/2,1/2,1/2,1/2}]+{1,1,1,1,1}]]&,ints,{1}];
     (*remove integrals with non-maximal number of denominators*)
     ints = Delete[ints,Position[numdens,n_ /; n<Max[numdens]]];
     (*number of powers for each integral*)
     numpows = Map[Total[Abs[Sequence[#]]]&,ints,{1}];
     (*remove integrals with non-maximal number of powers*)
     ints = Delete[ints,Position[numpows,n_ /; n<Max[numpows]]];
     (*For the toy example no further ordering is imposed, pick the first integral of the remaining ones*)
     int = F@@(ints[[1]]);
  (*Solve the IBP for this integral*)
     res = Collect[(int /. {Solve[ibp==0, int][[1]]})[[1]],F[__],Simplify];
     If[TrueQ[OptionValue[PrintInfo]],Print[{int,res}]];
  (*Apparently needs some trickery to add this to the function definition*)
     ReleaseHold[{int,res} /. {a_,b_}:>Hold[a=b]];
  ];


(* Process all IBPs for a given set of indices *)
ProcessIBPs[{a1_,a2_,a3_,a4_,a5_},opts___] := (
ProcessIBP[1,{a1,a2,a3,a4,a5},opts];
ProcessIBP[2,{a1,a2,a3,a4,a5},opts];
ProcessIBP[3,{a1,a2,a3,a4,a5},opts];
ProcessIBP[4,{a1,a2,a3,a4,a5},opts];
ProcessIBP[5,{a1,a2,a3,a4,a5},opts];
ProcessIBP[6,{a1,a2,a3,a4,a5},opts];
);


(*Reduce a sector up to maxpows extra powers*)
ReduceSector[dir_List,maxpows_,opts___] := Module[{a0,i,auxlist,ProcessIBPsaux},
   If[Length[dir] =!= 5,Print["Incompatible input!"]; Return[Null]];
   If[maxpows<0,Return[Null]];
   (*Index vector with no extra powers, i.e. the smallest Sum[Abs[ai]] in the sector*)
      a0=Table[Sign[Sign[dir[[j]]-1/2]+1],{j,1,5}];
   (*Start from lowest number of extra powers and increase until the cutoff is reached*)
      Print["Processing ",a0,"."];
      ProcessIBPs[a0,opts];
      i=1;
      While[i<=maxpows,
         Print["Processing ",a0," plus ",i," extra powers."];
         With[{aux1=Table[ToExpression["j"<>ToString[k]],{k,1,i}],
               aux2=Sequence@@Table[{ToExpression["j"<>ToString[k]],1,5},{k,1,i}]},
            auxlist = ProcessIBPsaux/@(Flatten[Table[a0+Sum[dir[[j]]*Normal[SparseArray[j->1,5]],{j,aux1}],aux2],i-1]);
            auxlist /. {ProcessIBPsaux[x_] :> ProcessIBPs[x,opts]}
         ];
         i++;
      ];
];


(* ::Subsection::Closed:: *)
(*Reduction*)


(*Start in lowest non-trivial (i.e. not scaleless) sector with direction {1,-1,-1,1,1}*)
(*Start with lowest possible powers of indices*)
ProcessIBP[1,{1,0,0,1,1}]


(*Definition of the integral has been added to boundary and symmetry conditions*)
?F


(*Do this for all IBPs with {1,0,0,1,1}*)
ProcessIBPs[{1,0,0,1,1}]


(*Now add one extra power*)
ProcessIBPs[{2,0,0,1,1}]


(*Consider all integrals in the sector with one extra power*)
ReduceSector[{1,-1,-1,1,1},1]


(*Consider all integrals in the sector with up to 3 extra powers*)
ReduceSector[{1,-1,-1,1,1},3,PrintInfo->False]


(*We find one master integral in the sector*)
F[4,-2,0,1,1]//Simplify
F[1,-1,-2,2,2]//Simplify
F[2,-1,-1,2,2]//Simplify


(*Move to next complicated sectors*)
ReduceSector[{1,1,1,1,-1},3,PrintInfo->False]


F[3,2,1,1,-1]//Simplify
F[3,1,3,1,-1]//Simplify


ReduceSector[{1,1,1,-1,1},3,PrintInfo->False]


F[3,2,1,-1,1]//Simplify
F[3,1,3,-1,1]//Simplify
F[1,1,3,-1,3]//Simplify


ReduceSector[{1,1,-1,1,1},3,PrintInfo->False]


(*Most complicated sector*)
ReduceSector[{1,1,1,1,1},3,PrintInfo->False]


F[2,2,1,2,2]


(* ::Subsection::Closed:: *)
(*Runtime comparison*)


(*Consider F[2,2,1,2,2]*)


(*Naive strategy (however already contains elements of Laporta): start in sector of desired integral*)


Timing[ReduceSector[{1,1,1,1,1},3,PrintInfo->False]]


F[2,2,1,2,2]//Simplify


Timing[ReduceSector[{-1,1,1,1,1},3,PrintInfo->False]]


F[2,2,1,2,2]//Simplify


Timing[ReduceSector[{1,1,1,1,-1},3,PrintInfo->False]]


F[2,2,1,2,2]//Simplify


Timing[ReduceSector[{1,-1,-1,1,1},3,PrintInfo->False]]


F[2,2,1,2,2]//Simplify


(*Laporta strategy*)


Timing[
ReduceSector[{1,-1,-1,1,1},3,PrintInfo->False];
ReduceSector[{1,1,1,1,-1},3,PrintInfo->False];
ReduceSector[{1,1,1,-1,1},3,PrintInfo->False];
ReduceSector[{1,1,-1,1,1},3,PrintInfo->False];
ReduceSector[{1,1,1,1,1},3,PrintInfo->False]]


F[2,2,1,2,2]