(* ::Package:: *)

(* ::Subtitle:: *)
(*Full solution vector*)


SetDirectory @ NotebookDirectory[];
Import["../utils.m"]
Off[Set::nosym]


(* ::Section:: *)
(*Extra utils*)


(* ::Text:: *)
(*Import AMFlow[2201.11669] results and map phase-space points to z-variables*)


import[file_] := Import[file] // RightComposition[
	ToExpression,
	Chop,
	ReplaceAll[eps -> \[Epsilon]],
	MapAt[
		Append[#, {z1 -> M/(-s) , z2 -> (-t)/(-s)} /. #]&,
		1
	],
	MapAt[
		Flatten,
		1
	],
	ReplaceAll // Map
];


(* ::Text:: *)
(*Compare \[Epsilon] coefficients between two results*)


compare[res1_List, res2_List] := Module[{coeffs,coeff1,coeff2},
	coeffs = Table[
		Coefficient[#, \[Epsilon], i],
		{i,-2,4}
	]&;
	
	coeff1 = res1 // Map@coeffs;
	coeff2 = res2 // Map@coeffs;
	
	(coeff1+1)/(coeff2+1) // Simplify
]


(* ::Section:: *)
(*Import*)


(* ::Text:: *)
(*8 AMFlow results for massive and masless master integrals.*)
(*pspt[n] contains the phase-space point for (s, t, M=m^2).*)
(*(s, t) points are the same for the the massive and massless results, so we only save one list pspt[n].*)


Do[
	{pspt[n], massiveNum[n]}  = import["numerical_master_integrals/amflow_massive_pspt"  <> ToString[n] <> ".txt"];
	{_ ,      masslessNum[n]} = import["numerical_master_integrals/amflow_massless_pspt" <> ToString[n] <> ".txt"],
	{n,8}
]


(* ::Text:: *)
(*Example*)


pspt[1]
Print[""]

massiveNum[1]
Print[""]

masslessNum[1]


(* ::Text:: *)
(*Master integrals*)


massiveMIs  = Import["massive_master_integrals.m"]
masslessMIs = Import["massless_master_integrals.m"]


(* ::Section:: *)
(*Solution vectors*)


(* ::Text:: *)
(*Gauge transformation G0 from eq. (5.18) of the restriction paper.*)
(*This was used to rescale massive master integrals to become unitless.*)


G0 = {(-s)^(1-\[Epsilon]), (-s)^-\[Epsilon], (-s)^-\[Epsilon], (-s)^(-1-\[Epsilon]), (-s)^(-2-\[Epsilon])} // DiagonalMatrix


(* ::Text:: *)
(*Solution vector for eigenvalue \[Lambda]1 = 0.*)
(*"box" refers to the massless box.*)


I\[Lambda]1 = Import["holomorphic_solution_vector.m"] /. j[\[Nu]__] :> j[box, \[Nu]]


(* ::Text:: *)
(*Solution vector for eigenvalue \[Lambda]2 = -\[Epsilon].*)


I\[Lambda]2 = Import["log_solution_vector.m"] // Factor


(* ::Text:: *)
(*Full solution vector expressed in (s, t, M=m^2)*)


Isol = G0 . (I\[Lambda]1 + I\[Lambda]2) //. {z1 -> M/(-s), z2 -> t/s, (M/(-s))^p_ * (-s)^(q_:1) :> M^p *(-s)^(q-p)}


(* ::Section:: *)
(*Fix boundary constants C[1], C[2]*)


(* ::Subsection:: *)
(*C[1]*)


(* ::Text:: *)
(*Exact result for the massive 1-loop tadpole j[0,0,1,0], where M = m^2.*)


tadpole = - Gamma[\[Epsilon]-1] * M^(1-\[Epsilon])


(* ::Text:: *)
(*Match with first element of the full solution vector *)


c1 = Solve[tadpole == (Isol // First), C[1]] // Flatten // ReplaceAll


(* ::Subsection:: *)
(*C[2]*)


(* ::Text:: *)
(*Rule to expand boundary constant C[2] in \[Epsilon]*)


c2Expand = C[2] -> Sum[c[n] \[Epsilon]^n , {n,-2,4}] // ReplaceAll 


(* ::Text:: *)
(*Evaluate full solution vector at one phase space point and expand C[2] to get equations for c2[n].*)
(*Note that Isol depends on the massless integrals because of I\[Lambda]1, so we evaluate those with masslessNum[1].*)


IsolNum = Isol // RightComposition[
	pspt[1],
	masslessNum[1],
	N,
	c1,
	c2Expand,
	Series[#, {\[Epsilon],0,4}]&,
	Normal,
	Collect[#, \[Epsilon], FullSimplify]&
];


(* ::Text:: *)
(*We match with the numerical evaluation of the 5th massive master integral (the Bhabha box)*)


c2 = IsolNum - massiveMIs // RightComposition[
	Part[#, 5]&,
	massiveNum[1],
	CoefficientArrays[\[Epsilon]^2 * # // Expand, \[Epsilon]]&,
	Flatten,
	Equal[#, 0]&,
	Solve,
	Flatten,
	ReplaceAll
]


(* ::Text:: *)
(*So we have found C[2]*)


C[2] // c2Expand // c2


(* ::Section:: *)
(*Numerical checks*)


(* ::Text:: *)
(*Redefine Isol to include results for c1 and c2*)


Isol = Isol // c2Expand // c1 // c2 // Factor;


(* ::Text:: *)
(*Our result evaluated at the nth phase-space point*)


myRes[n_] := Isol // RightComposition[
	pspt[n],
	masslessNum[n],
	N,
	Series[#, {\[Epsilon],0,4}]&,
	Normal,
	Collect[#, \[Epsilon], FullSimplify]&
]


(* ::Text:: *)
(*Compare our \[Epsilon] coefficients with numerical evaluations at other phase-space points.*)
(*The (i,j)th entry of the matrices below represents the ith MI and the jth \[Epsilon] coefficient.*)
(*We take ratios between our results and those of AMFlow, so the results should be close to 1.*)
(*The triangle (4th row) and the box (5th row), the only integrals to receive logarithmic corrections, agree well with the numerical results.*)


Do[
	"Phase space point " <> ToString[n] // Print;
	pspt[n] // Print;
	res1 = myRes[n];
	res2 = massiveMIs // massiveNum[n];
	compare[res1, res2] // MatrixForm // Print,
	{n,2,8}
]



