(* ::Package:: *)


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


(* ::Section:: *)
(*Restriction unequal to equal mass sunrise*)


(* ::Text:: *)
(*3-mass sunrise: z1 = M1/s , z2 = M2/s , z3 = M3/s*)


{P1, P2, P3} = Import["unequal_mass_pfaffians.m"];


(* ::Text:: *)
(*Equal mass limit z1, z2, z3 -> z, and in new variables x1, x2 -> 0*)


shiftz = ReplaceAll[{z1 -> z, z2 -> z + x1, z3 -> z + x2}];


Pz = P1 + P2 + P3 // shiftz // Factor;
Px1 = P2 // shiftz // Factor;
Px2 = P3 // shiftz // Factor;


(* ::Text:: *)
(*All Pfaffian matrices are finite in the x1, x2 -> 0 limit*)


{Pz, Px1, Px2} // mseries1[x1, 0]


{Pz, Px1, Px2} // mseries1[x2, 0]


(* ::Text:: *)
(*Master integrals for the unequal mass system*)


mastersuneqm = Import["unequal_mass_master_integrals.m"]


(* ::Text:: *)
(*Master integrals for the equal mass system*)


masterseqm = Import["equal_mass_master_integrals.m"]


B = CoefficientArrays[masterseqm, mastersuneqm] [[2]] // Normal;
B // MatrixForm


(* ::Text:: *)
(*The symmetry relations encoded in R-matrix*)


R = {
    {1, -1, 0, 0, 0, 0, 0},
    {1, 0, -1, 0, 0, 0, 0},
    {0, 0, 0, 0, 1, -1, 0},
    {0, 0, 0, 0, 1, 0, -1},
    Nothing
} // rowReduce;


(* ::Text:: *)
(*M-matrix*)


M = Join[B, R]


(* ::Text:: *)
(*Restriction*)


Q = (D[M, z] + M . (Pz /. x1 | x2 -> 0)) . Inverse[M] // RightComposition[
	#[[;;3 , ;;3]]&,
	Factor,
	Identity
];


Qanswer = Import["equal_mass_pfaffian_answer.m"]


Q - Qanswer // Simplify


(* ::Section::Closed:: *)
(*Derive the symmetry relations using the secondary-like equation*)


(* ::Text:: *)
(*Derive the secondary-like equation*)


sys = Pz // RightComposition[
    ReplaceAll[x1 | x2 -> 0],
    Factor,
    {-#, #}&,
    Transpose[#, {1, 3, 2}]&,
    Outer[Times, #[[1]], #[[1]] // Length // IdentityMatrix, 2]
    +
    Outer[Times, #[[2]] // Length // IdentityMatrix, #[[2]], 2]&,
    Flatten[#, {{1, 4}, {2, 3}}]&,
    Identity
];


(* ::Text:: *)
(*The same system generated explicitely*)
Y = Array[y, Pz // Dimensions] // Map2[ind2sym]
ys = Y // Flatten;
sysY = Transpose[-Pz /. x1 | x2 -> 0] . Y + Y . Transpose[Pz /. x1 | x2 -> 0] // CoefficientArrays[#, ys]& // Last // Flatten[#, {{1, 2}, {3}}]& // Normal;
sys - sysY // Simplify // Flatten // ContainsOnly[{0}]


(* ::Text:: *)
(*Export to Maple*)


MaplePut[sys, "secondary_equation.mma", "sys"]


(* ::Text:: *)
(*At this point you should run the `secondary.mpl` script*)


(* ::Text:: *)
(*Import from Maple*)


sol = MapleGet["secondary_equation_solutions.mma"] // Transpose // Map[Partition[#, Pz // Length]&];


(* ::Text:: *)
(*Check against R-matrix (the specific eigenspaces might be different for you)*)


Rmaple = Join[
    sol[[2]] + IdentityMatrix[7] // getNullSpace,
    sol[[3]] // getNullSpace,
    sol[[5]] - IdentityMatrix[7] // getNullSpace,
    sol[[1]] + IdentityMatrix[7] // getNullSpace,
    {}
] // rowReduce // Normal;

Rmaple - R // Normal // Factor // Flatten // ContainsOnly[{0}]
