(* ::Package:: *)

(* ::Title:: *)
(*From integrals to D-module basis*)


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


(* ::Section:: *)
(*GKZ data*)


(* ::Text:: *)
(*GKZ A-matrix*)


myA = {
	{1,1,1,1,1,1,1,1,1,1},
	{1,1,1,0,0,0,0,0,1,0},
	{1,0,0,1,1,1,0,0,1,1},
	{0,0,0,1,0,0,1,1,0,1},
	{0,1,0,0,1,0,1,0,1,0},
	{0,0,1,0,0,1,0,1,0,1}
};


(* ::Text:: *)
(*Nullspace of the A-matrix (needed for R-vectors)*)


ns = myA // NullSpace;
% // Map@MatrixForm


(* ::Text:: *)
(*Simplex*)


sigma = {1,2,3,4,5,9}


(* ::Text:: *)
(*Exponent parameters*)


beta = Range[6] // Map@b // Map@ind2sym


(* ::Section:: *)
(*Integral basis*)


(* ::Text:: *)
(*Integral notation: j[d0/2, propagators] where d = d0-2*eps*)
(*We have two sunrises and one N-box.*)


ints = {j[4,0,2,2,0,2], j[4,2,2,0,2,0], j[4,1,1,1,1,1]}


(* ::Text:: *)
(*Corresponding R-vectors*)


rvecs = ints /. j[d0_,n__] :> rvec[myA, {d0/2,n}, ns]


(* ::Text:: *)
(*asir script to find D-module basis*)


script = asirscript[myA, beta, sigma, rvecs]


(* ::Text:: *)
(*We import the result of this script.*)
(*Then we replace b-exponents with \[Epsilon] and \[Delta] and rewrite to j-notation given by j[a,b,c,d] = dx6^a * dx7^b * dx8^c * dx10^d.*)


basis = Import["d_module_basis.txt"] // ToExpression
basis = basis // niceop[#, {dx6, dx7, dx8, dx10}]&


(* ::Section:: *)
(*B-matrix for restriction*)


(* ::Text:: *)
(*Standard monomials*)


std = {j[0,0,2,0], j[1,0,0,1], j[0,1,0,1], j[0,0,0,2], j[1,0,0,0], j[0,1,0,0], j[0,0,1,0], j[0,0,0,1], j[0,0,0,0]}


(* ::Text:: *)
(*The B-matrix in equation (5.12)*)


B = CoefficientArrays[basis, std] [[2]] // Normal // Factor;
B // MatrixForm


Export["B_matrix.m", B]
