(* ::Package:: *)

(* ::Subtitle:: *)
(*Pfaffians from MM *)


(* ::Text:: *)
(*load FiniteFlow*)


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


(* ::Text:: *)
(*copy pasted functions*)


numrep[params_]:=ReplaceAll[Thread[(#1->(1/RandomPrime[{10,10^5},#1]&)[Length[#1]]&)[params]]];
coeffarr[dir_,list_]:=With[{derdir=der[std,dir]},If[derdir\[Intersection]list=={},ConstantArray[0,Length/@{list,list}],Normal[CoefficientArrays[der[std,dir],list][[2]]]]];
der[js_,dir_Integer]:=Module[{len,m},len=Length[First[filter[_j][js]]];If[dir>len,Print["dir > len"],m=UnitVector[len,dir];js/. j[n__]:>j@@({n}+m)]];
filter[pattern_]:=Union[Flatten[Cases[#1,pattern,-1]]]&;
indrows[mat_List]:=(Transpose/*rowReduce/*Transpose/*Apply[List]/*Position[x_/;Length[x["NonzeroValues"]]===1]/*Flatten/*Rest)[mat];
rowReduce=RowReduce/*SparseArray/*Apply[List]/*DeleteCases[x_/;x["Density"]===0.`]/*SparseArray;


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


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


data = Import["macaulay_data.txt"] // StringReplace[{"(1)*<<"->"j[",">>"->"]"}] // ToExpression // Factor;


{m1, m2, ext, std} = data[[3]];


Dimensions /@ {m1, m2, ext, std} // Column


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


mons = Join[ext, std];


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


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


m1n = m1 // num;
MatrixRank[m1n]


(* ::Subsection:: *)
(*C1' and C2'*)


(* ::Text:: *)
(*18 is the Pfaffian index data[[1]]*)


c1 = coeffarr[18,ext];
c2 = coeffarr[18,std];
c12 = coeffarr[18,mons];


(* ::Text:: *)
(*C1' and C2'*)


c1 // Dimensions
c2 // Dimensions


(* ::Text:: *)
(*check*)


c12 . mons - der[std,18]
c12 . mons - c1 . ext - c2 . std


(* ::Section:: *)
(*Matrix equation*)


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


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


rows = indrows[m1n];
m1ind = m1[[rows]];


(* ::Text:: *)
(*wish to solve for c[i,j] in meq == 0*)


Clear[cij,meq]

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

meq = c1 - cij . m1ind // Flatten;
meq = (meq == 0  // Thread // DeleteCases[#, True]&);


(* ::Subsection:: *)
(*solve with FiniteFlow*)


(* ::Text:: *)
(*solution*)


Clear[sol]

sol = FFSparseSolve[
	meq,
	meq // filter[_c]
];


(* ::Text:: *)
(*update the matrix*)


Clear[csol]

csol = cij //. Dispatch[sol] // Factor;
csol // filter[_c]


(* ::Text:: *)
(*check that solution satisfies matrix equation (should get {0})*)


c1 - num[csol] . num[m1ind] // Flatten // DeleteDuplicates


(* ::Section:: *)
(*Pfaffian*)


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


m2ind = m2[[rows]];


(* ::Text:: *)
(*in direction r*)


Clear[pr]
pr = c2 - csol . m2ind // Factor;


(* ::Text:: *)
(*rank check*)


pr;
pr // num // MatrixRank


ptaka = data[[2,1]];
ptaka - pr /. {e -> 1/7, d -> -7/13, r -> 1/7}


Export["2L_0M_box_pfaffan_matrix_r.m", pr]
