(* ::Package:: *)

(* ::Title:: *)
(*Pfaffian system from Macaulay matrix*)


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


(* ::Text:: *)
(*x-directions corresponding to {x6,x7,x8,x10} (called {z6,z7,z8,z10} in the paper)*)


dirs = {1,2,3,4}


(* ::Text:: *)
(*dir = x10  for tests*)


dir = dirs[[4]]


(* ::Section:: *)
(*Macaulay data*)


(* ::Subsection:: *)
(*import*)


(* ::Text:: *)
(*import Macalay data computed via asir*)


data = mmimport["macaulay_matrix.txt"];
{Mext, Mstd, ext, std} = data;


(* ::Text:: *)
(*the GKZ system has rank 9 = size of std*)


Dimensions /@ {Mext, Mstd, ext, std} // Column


(* ::Text:: *)
(*all monomials*)


mons = Join[ext, std];


(* ::Subsection:: *)
(*numerics*)


(* ::Text:: *)
(*numerical values for parameters to speed up checks*)


params = Mext // Variables
num = numrep[params]


(* ::Text:: *)
(*find number of independent rows in Mext*)


Mextnum = Mext // num;
MatrixRank[Mextnum]


(* ::Subsection:: *)
(*C-matrices*)


(* ::Text:: *)
(*d_xi \[Bullet] std =  Cmons \[Bullet] mons = Cext \[Bullet] ext + Cstd \[Bullet] std*)


Cext[dir_]  := coeffarr[dir,ext];
Cstd[dir_]  := coeffarr[dir,std];
Cmons[dir_] := coeffarr[dir,mons];


(* ::Text:: *)
(*dimensions*)


Cext[dir] // Dimensions
Cstd[dir] // Dimensions


(* ::Text:: *)
(*check that deriv[x8] = d_8 \[Bullet] std*)


Cmons[dir] . mons - der[std,dir]
Cmons[dir] . mons - Cext[dir] . ext - Cstd[dir] . std


(* ::Section:: *)
(*Solve matrix equation*)


(* ::Subsection:: *)
(*setup*)


(* ::Text:: *)
(*independent rows*)


rows = indrows[Mextnum]
Mextind = Mext[[rows]];


(* ::Text:: *)
(*wish to solve for c[i,j]*)


cij = Table[c[i, j], {i, std // Length}, {j, Mextind // Dimensions // First}];


(* ::Text:: *)
(*matrix equatin meq == 0*)


Clear[meq]
meq[dir_] := meq[dir] = Cext[dir] - cij . Mextind // DeleteCases[#,0]& // Flatten;


(* ::Subsection:: *)
(*solve*)


(* ::Text:: *)
(*solve matrix equation for c's with FiniteFlow [1905.08019]*)


Clear[sol,csol]

sol[dir_] := sol[dir] = FFSparseSolve[
	meq[dir] == 0 // Thread,
	meq[dir] // filter[_c]
]

csol[dir_] := csol[dir] = cij //. Dispatch[sol@dir] // Factor;


(* ::Text:: *)
(*update the csol matrix and check that no c's are left*)


csol[dir];
%[[All,2]] // filter[_c]


(* ::Text:: *)
(*check that solution satisfies matrix equation*)


Cext[dir] - num[csol@dir] . num[Mextind] // Flatten // DeleteDuplicates


(* ::Section:: *)
(*Pfaffian system in std basis*)


(* ::Text:: *)
(*independent rows for Mstd*)


Mstdind = Mstd[[rows]];


(* ::Text:: *)
(*Pfaffian matrices*)


P[dir_] := P[dir] = Cstd[dir] - csol[dir] . Mstdind // Factor;


(* ::Text:: *)
(*Compute for all directions and save to file*)


Export["pfaffians_std_basis.m", P /@ dirs]


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


P[1]
P[1] // Dimensions



