(* ::Package:: *)

(* ::Subtitle:: *)
(*Restriction of Pfaffian matrices coming from the GKZ system*)


SetDirectory @ NotebookDirectory[];
Get["../utils.m"];


(* ::Section:: *)
(*Pfaffians*)


(* ::Text:: *)
(*Pfaffian matrices from the Macaulay matrix method*)


{P@6, P@7, P@8, P@10} = Import["pfaffians_std_basis.m"];


(* ::Text:: *)
(*Seek x6, x7, x8 --> 1, or equivalently r6, r7, r8 --> 0*)


x2r = ReplaceAll[{x6 -> 1 + r6, x7 -> 1 + r7, x8 -> 1 + r8, x10 -> r10}]


(* ::Text:: *)
(*Rewrite Pfaffians in r-variables*)


P@6  = P@6  // x2r // Factor;
P@7  = P@7  // x2r // Factor;
P@8  = P@8  // x2r // Factor;
P@10 = P@10 // x2r // Factor;


(* ::Section:: *)
(*Residue matrices*)


(* ::Text:: *)
(*Residue matrix for P6 as r6, r7 --> 0*)


R6[1] = P[6] // RightComposition[
	mseries1[r6, -1],
	Part[#, 2, 1]&,
	Normal,
	mseries1[r7, 0],
	Part[#, 2, 1]&
]


(* ::Text:: *)
(*Residue matrix for P7 as r6, r7 --> 0*)


R7[1] = P[7] // RightComposition[ 
	mseries1[r7, -1],
	Part[#, 2, 1]&,
	Normal,
	mseries1[r6, 0],
	Part[#, 2, 1]&
]


(* ::Section:: *)
(* Normal form w.r.t. the expansion order {r6, r7, r8}*)


(* ::Text:: *)
(*Shearing transformation w.r.t. variable r8 (cf. Moser reduction)*)


D8[1] = Join[ConstantArray[1 / r8, 1], {}] // PadRight[#, 9, 1]& // DiagonalMatrix


(* ::Subsection:: *)
(*Bring P[8] to normal form*)


(* ::Text:: *)
(*Expand P[8] up to O(r6^0) + O(r7^0)*)


A8[0] = P[8] // RightComposition[
	mseries1[r6, 0],
	Part[#, 2, 1]&,
	Normal,
	mseries1[r7, 0],
	Part[#, 2, 1]&,
	Normal,
	Factor
];


(* ::Text:: *)
(*Moser reduction via shearing.*)
(*A8[1] has a simple pole in r8, whereas P[8] has higher-order poles*)


A8[1] = A8[0] // gauge[D8[1], r8] // Factor;


(* ::Text:: *)
(*Residue matrix for A8[1]*)


R8[1] = A8[1] // RightComposition[
	mseries1[r8, 0],
	Part[#, 2, 1]&,
	Normal,
	Factor
];


(* ::Text:: *)
(*New residue matrix for P[6]*)


R6[2] = R6[1] // RightComposition[
	gauge[D8[1], r6],
	Factor,
	mseries1[r8, 0],
	Part[#, 2, 1]&,
	Normal
];


(* ::Text:: *)
(*New residue matrix for P[7]*)


R7[2] = R7[1] // RightComposition[
	gauge[D8[1], r7],
	Factor,
	mseries1[r8, 0],
	Part[#, 2, 1]&,
	Normal
];


(* ::Subsection:: *)
(*Bring P[10] to normal form*)


(* ::Text:: *)
(*Expand P[10] in {r6, r7, r8} and apply the gauge transformation with D8[1]*)


P10 = P[10] // RightComposition[
    mseries1[r6, 0],
    Part[#, 2, 1]&,
    Normal,
    mseries1[r7, 0],
    Part[#, 2, 1]&,
    Normal,
    Factor,
    gauge[D8[1], r10],
    Factor,
    mseries1[r8, 0],
    Part[#, 2, 1]&,
    Normal,
    Factor,
    Identity
];


(* ::Text:: *)
(*P[10] is now finite in the limit*)


P10 /. {r6 | r7 | r8 -> 0};


(* ::Section:: *)
(*Restriction*)


(* ::Subsection:: *)
(*Matrices*)


(* ::Text:: *)
(*B-matrix (see eq. (5.11) in the restriction paper)*)


B = Import["B_matrix.m"] // x2r // ReplaceAll[r6 | r7 | r8 -> 0] // Factor;


(* ::Text:: *)
(*R-matrix (see eq. (5.11) in the restriction paper)*)


R = Join[R6[2], R7[2], R8[1]] // rowReduce;


(* ::Text:: *)
(*M-matrix by stacking B on top of R*)


M = Join[B, R];


(* ::Text:: *)
(*Check that it has full rank*)


M // MatrixRank


(* ::Text:: *)
(*Lee-Pomeransky prefactors and factors of s from rescaling integration variables zi --> zi/s*)


LPfacs = DiagonalMatrix[{
	(s^(-4 - \[Epsilon] - 9 \[Delta] \[Epsilon]) Gamma[2 - \[Epsilon]])/( Gamma[\[Delta] \[Epsilon]]^6 Gamma[-3 \[Epsilon] - 9 \[Delta] \[Epsilon]] Gamma[2 + \[Delta] \[Epsilon]]^3), 
	(s^(-4 - \[Epsilon] - 9 \[Delta] \[Epsilon]) Gamma[2 - \[Epsilon]])/(Gamma[\[Delta] \[Epsilon]]^6 Gamma[-3 \[Epsilon] -9 \[Delta] \[Epsilon]] Gamma[2 + \[Delta] \[Epsilon]]^3), 
	(s^(-3 - \[Epsilon] - 9 \[Delta] \[Epsilon]) Gamma[2 - \[Epsilon]])/( Gamma[\[Delta] \[Epsilon]]^4 Gamma[1 - 3 \[Epsilon] - 9 \[Delta] \[Epsilon]] Gamma[1 + \[Delta] \[Epsilon]]^5)
}];


(* ::Text:: *)
(*Extra factors of s to render the master integrals dimensionless*)


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


(* ::Text:: *)
(*Combine all prefactors (see eq. (5.8))*)


facs = LPfacs . G0


(* ::Subsection:: *)
(*Restriction*)


(* ::Text:: *)
(*Equation (3.14).*)
(*We then include prefactors and send the regular \[Delta] --> 0*)


Q10 = (D[M,r10] + M . P10) . Inverse[M] // RightComposition[
	#[[;;3 , ;;3]]&,
	Factor,
	conj[facs // Inverse],
	FunctionExpand,
	ReplaceAll[\[Delta] -> 0],
	Factor,
	Identity
];


(* ::Text:: *)
(*This result agrees with an independent calculation using LiteRed [1310.1145]*)


Q10 // MatrixForm
