(* ::Package:: *)

(************************************************************************)
(* This file was generated automatically by the Mathematica front end.  *)
(* It contains Initialization cells from a Notebook file, which         *)
(* typically will have the same name as this file except ending in      *)
(* ".nb" instead of ".m".                                               *)
(*                                                                      *)
(* This file is intended to be loaded into the Mathematica kernel using *)
(* the package loading commands Get or Needs.  Doing so is equivalent   *)
(* to using the Evaluate Initialization Cells menu command in the front *)
(* end.                                                                 *)
(*                                                                      *)
(* DO NOT EDIT THIS FILE.  This entire file is regenerated              *)
(* automatically each time the parent Notebook file is saved in the     *)
(* Mathematica front end.  Any changes you make to this file will be    *)
(* overwritten.                                                         *)
(************************************************************************)



BeginPackage["FusionAtlas`GraphPlanarAlgebra`",{"FusionAtlas`","FusionAtlas`Bigraphs`","FusionAtlas`PackageData`","FusionAtlas`Debugging`"}];


Vertex;


GPAElement;


GPAMatrix;


NumberFieldGenerator;EnlargeNumberField;


LopsidedDimension;


EmptyGPAElement;


GPACoefficients;


GPAMultiply;


RotateOneClick;RotateTwoClicks;


CapTopLeft;CapAt;


AddStrandOnRight; AddCupOnTop;


JonesWenzlIdempotent;


LowestWeightEigenvectorConditions;


LowestWeightEigenspace;


Begin["`Private`"];


EvenVertices[g_GradedBigraph]:=EvenVertices[g]=Flatten[Table[Vertex[d,k],{d,0,GraphDepth[g],2},{k,1,RankAtDepth[g,d]}]]
OddVertices[g_GradedBigraph]:=OddVertices[g]=Flatten[Table[Vertex[d,k],{d,1,GraphDepth[g],2},{k,1,RankAtDepth[g,d]}]]


EdgesAdjacentTo[g_GradedBigraph,Vertex[d_,k_]]:=EdgesAdjacentTo[g,Vertex[d,k]]=If[d<GraphDepth[g],Inner[List,Table[Vertex[d+1,l],{l,1,RankAtDepth[g,d+1]}],g[[d+1]][[All,k]],List],{}]~Join~If[d>0,Inner[List,Table[Vertex[d-1,l],{l,1,RankAtDepth[g,d-1]}],g[[d]][[k]],List],{}]/.{x_,n_Integer}:>Sequence@@Table[{x,l},{l,1,n}]


GraphPaths[g_GradedBigraph,shading:(True|False),0]:=GraphPaths[g,shading,0]={#}&/@If[shading,OddVertices[g],EvenVertices[g]]


GraphPaths[g_GradedBigraph,shading:(True|False),length_Integer]/;length>0:=GraphPaths[g,shading,length]=Flatten[Function[{path},path~Join~Reverse[#]&/@EdgesAdjacentTo[g,Last[path]]]/@GraphPaths[g,shading,length-1],1]


GraphPathsFrom[g_GradedBigraph,length_Integer,initial_Vertex]/;length>0:=Flatten[Function[{path},path~Join~Reverse[#]&/@EdgesAdjacentTo[g,Last[path]]]/@GraphPathsFrom[g,length-1,initial],1]


GraphPathsFrom[g_GradedBigraph,0,initial_Vertex]:={{initial}}


GraphPathsTo[g_GradedBigraph,length_Integer,final_Vertex]:=Reverse/@GraphPathsFrom[g,length,final]


GraphPathsBetween[g_GradedBigraph,length_Integer,initial_Vertex,final_Vertex]:=GraphPathsBetween[g,length,initial,final]=Cases[GraphPathsFrom[g,length,initial],{___,final}]


GraphLoops[g_GradedBigraph,shading:(True|False),length_Integer]/;length>=2:=GraphLoops[g,shading,length]=Cases[GraphPaths[g,shading,length],{a_,z___,a_}:>{a,z}]


NumberFieldGenerator[g_GradedBigraph]:=NumberFieldGenerator[g]=RootReduce[Sqrt[2-ReducedDimensionOfGenerator[g]^2]]


ToNumberField0[x_Real?NumericQ,y_?NumericQ]:=x
ToNumberField0[x_Complex?NumericQ,y_?NumericQ]:=x
ToNumberField0[x_?NumericQ,y_?NumericQ]:=ToNumberField[RootReduce[x],y]
ToNumberField0[x_,y_]:=x


EnlargeNumberField[g_GradedBigraph][x_?NumericQ]:=Module[{},
(*Off[ToNumberField::"nnfel"];*)
Switch[ToNumberField0[x,NumberFieldGenerator[g]],
_ToNumberField,
NumberFieldGenerator[g]=ToNumberField[{RootReduce[x],NumberFieldGenerator[g]}][[1,1]];
DebugPrint["Enlarging number field for ",GraphToString[g],": ",NumberFieldGenerator[g]],
_,0
];
(*On[ToNumberField::"nnfel"];*)
NumberFieldGenerator[g]
]


LopsidedDimension[g_GradedBigraph,Vertex[d_,k_]]:=LopsidedDimension[g,Vertex[d,k]]=
With[{\[Delta]=ReducedDimensionOfGenerator[g]},
ToNumberField0[RootReduce[If[OddQ[d],1/\[Delta],1]ReducedDimensionsByDepth[g][[d+1,k]]],NumberFieldGenerator[g]]
]


EmptyGPAElement[g_GradedBigraph,shading:(True|False)]:=GPAElement[g,shading,{0,0},{#}->1&/@If[shading,OddVertices[g],EvenVertices[g]]]


RotateOneClick[gpa:GPAElement[_GradedBigraph,_Boolean,{0,0},_List]]:=gpa
RotateTwoClicks[gpa:GPAElement[_GradedBigraph,_Boolean,{0,0},_List]]:=gpa



RotateOneClick[GPAElement[g_GradedBigraph,shading:(True|False),{0,up_Integer},coefficients_List]/;up>=2]:=GPAElement[g,!shading,{0,up},coefficients/.({a_Vertex,na_Integer,b_Vertex,nb_Integer,z___}->\[Zeta]_):>({b,nb,z,a,na}->\[Zeta] LopsidedDimension[g,a]/LopsidedDimension[g,b])]


RotateTwoClicks[GPAElement[g_GradedBigraph,shading:(True|False),{0,up_Integer},coefficients_List]/;up>=4]:=GPAElement[g,shading,{0,up},coefficients/.({a_Vertex,na_Integer,b_Vertex,nb_Integer,c_Vertex,nc_Integer,z___}->\[Zeta]_):>({c,nc,z,a,na,b,nb}->\[Zeta] LopsidedDimension[g,a]/LopsidedDimension[g,c])]


RotateOneClick[GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients_List]/;up>=2\[And]down>=1]:=GPAElement[g,!shading,{down,up},coefficients/.({z_Vertex,nz_Integer,A:RepeatedSequence[_Vertex,_Integer,{up-2}],b_Vertex,nb_Integer,D:RepeatedSequence[_Vertex,_Integer,{up-1}],e_Vertex,ne_Integer}->\[Zeta]_):>({e,ne,z,nz,A,b,nb,D}->\[Zeta] LopsidedDimension[g,z]/LopsidedDimension[g,b])]


RotateTwoClicks[GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients_List]/;up>=3\[And]down>=2]:=GPAElement[g,shading,{down,up},coefficients/.({z_Vertex,nz_Integer,A:RepeatedSequence[_Vertex,_Integer,{up-3}],b_Vertex,nb_Integer,c_Vertex,nc_Integer,D:RepeatedSequence[_Vertex,_Integer,{up-2}],e_Vertex,ne_Integer,f_Vertex,nf_Integer}->\[Zeta]_):>({e,ne,f,nf,z,nz,A,b,nb,c,nc,D}->\[Zeta] (LopsidedDimension[g,z]LopsidedDimension[g,f])/(LopsidedDimension[g,b]LopsidedDimension[g,c]))]


CapTopLeft[GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients_List]]:=
With[{dispatch=Dispatch[coefficients]},
ExpandCoefficientRules[
GPAElement[g,shading,{down,up-2},{{a_Vertex,na_Integer,z___}:>1/LopsidedDimension[g,a] Plus@@({a,#[[2]],#[[1]],#[[2]],a,na,z}/.dispatch&/@EdgesAdjacentTo[g,a])}]
]
]


CapAt[GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients_List],k_Integer]/;1<=k<up:=
With[{dispatch=Dispatch[coefficients]},
ExpandCoefficientRules[
GPAElement[g,shading,{down,up-2},{{A:RepeatedSequence[_Vertex,_Integer,{k-1}],a_Vertex,na_Integer,z___}:>1/LopsidedDimension[g,a] Plus@@({A,a,#[[2]],#[[1]],#[[2]],a,na,z}/.dispatch&/@EdgesAdjacentTo[g,a])}]
]
]


PartialTrace[GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer}/;down>=2\[And]up>=2,coefficients_List]]:=
ExpandCoefficientRules[
GPAElement[g,shading,{down-1,up-1},{{A:RepeatedSequence[_Vertex,_Integer,{up-1}],a_Vertex,na_Integer,B:RepeatedSequence[_Vertex,_Integer,{down-2}]}:>Plus@@((LopsidedDimension[g,#[[1]]]/LopsidedDimension[g,a] ({A,a,#[[2]],#[[1]],#[[2]],a,na,B}/.coefficients))&/@EdgesAdjacentTo[g,a])}]
]


RepeatedSequence[A__,repeat:(_Integer|{_Integer}|{_Integer,_Integer}):\[Infinity]]:=Repeated[PatternSequence[A],repeat]


TurnUpBottomRightCorner[GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients_]/;down>=1]:=
GPAElement[g,shading,{down-1,up+1},
coefficients/.({A:RepeatedSequence[_Vertex,_Integer,{up}],v_Vertex,B___}->\[Zeta]_):>({A,v,B}->\[Zeta] LopsidedDimension[g,v])
]


TurnDownTopRightCorner[GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients_]/;up>=1]:=
GPAElement[g,shading,{down+1,up-1},coefficients/.({A:RepeatedSequence[_Vertex,_Integer,{up-1}],v_Vertex,B___}->\[Zeta]_):>({A,v,B}->\[Zeta] 1/LopsidedDimension[g,v])
]


TurnUpBottomLeftCorner[GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients_]/;down>=1]:=
GPAElement[g,!shading,{down-1,up+1},
coefficients/.({a_Vertex,na_Integer,B___,b_Vertex,nb_Integer}->\[Zeta]_):>({b,nb,a,na,B}->\[Zeta] LopsidedDimension[g,a])
]


TurnDownTopLeftCorner[GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients_]/;up>=1]:=
GPAElement[g,!shading,{down+1,up-1},
coefficients/.({a_Vertex,na_Integer,b_Vertex,nb_Integer,B___}->\[Zeta]_):>({b,nb,B,a,na}->\[Zeta] 1/LopsidedDimension[g,b])]


AddStrandOnRight[GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients_]/;down>=1]:=ExpandCoefficientRules[GPAElement[g,shading,{down+1,up+1},{{A:RepeatedSequence[_Vertex,_Integer,{up}],a_Vertex,na_Integer,b_Vertex,na_Integer,a_Vertex,nc_Integer,B:RepeatedSequence[_Vertex,_Integer,{down-1}]}:>({A,a,nc,B}/.coefficients)}]]


AddStrandOnRight[GPAElement[g_GradedBigraph,shading:(True|False),{0,up_Integer},coefficients_]/;up>=1]:=ExpandCoefficientRules[GPAElement[g,shading,{1,up+1},{{a_Vertex,nc_Integer,A:RepeatedSequence[_Vertex,_Integer,{up-1}],a_Vertex,na_Integer,b_Vertex,na_Integer}:>{a,nc,A}/.coefficients}]]


AddStrandOnRight[GPAElement[g_GradedBigraph,shading:(True|False),{0,0},coefficients_]]:=
ExpandCoefficientRules[GPAElement[g,shading,{1,1},{{a_Vertex,na_Integer,b_Vertex,na_Integer}:>({a}/.coefficients)}]]


AddCupOnTop[GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients_],cupPosition_Integer]/;(up+down>=1 \[And]1<= cupPosition<=up+1):=ExpandCoefficientRules[GPAElement[g,shading,{down,up+2},{{A:RepeatedSequence[_Vertex,_Integer,{cupPosition-1}],a_Vertex,na_Integer,b_Vertex,na_Integer,a_Vertex,nc_Integer,B:RepeatedSequence[_Vertex,_Integer,{up+down-cupPosition}]}:>LopsidedDimension[g,b]* ({A,a,nc,B}/.coefficients)}]]


AddCupOnTop[GPAElement[g_GradedBigraph,shading:(True|False),{0,0},coefficients_],1]:=ExpandCoefficientRules[GPAElement[g,shading,{2,0},{{a_Vertex,na_Integer,b_Vertex,na_Integer}:>LopsidedDimension[g,b]*({a}/.coefficients)}]]


GPAMultiply[GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,middle_Integer},coefficients1_],GPAElement[g_GradedBigraph,shading:(True|False),{middle_Integer,up_Integer},coefficients2_]]:=ExpandCoefficientRules[
GPAElement[g,shading,{down,up},{{A:RepeatedSequence[_Vertex,_Integer,{up}],B:RepeatedSequence[_Vertex,_Integer,{down}]}:>Plus@@(({A}~Join~Most[#]/.coefficients2)(Most[Reverse[#]]~Join~{B}/.coefficients1)&/@GraphPathsBetween[g,middle,{B}[[1]],{A}[[1]]])}
]
]


qInteger\[Delta][k_]:=qInteger\[Delta][k]=Function[{\[Delta]},Evaluate[Collect[\[Delta] qInteger\[Delta][k-1][\[Delta]]-qInteger\[Delta][k-2][\[Delta]],\[Delta]]]]
qInteger\[Delta][2]=#&;
qInteger\[Delta][1]=1&;
qInteger\[Delta][0]=0&;


(*
WenzlRecursion[gpa_GPAElement,x_]:=
AddStrandOnRight[gpa]-x GPAMultiply[TurnDownTopRightCorner[gpa],TurnUpBottomRightCorner[gpa]]
*)


(*
JonesWenzlIdempotent[g_GradedBigraph,0,shading:(True|False)]:=EmptyGPAElement[g,shading];
JonesWenzlIdempotent[g_GradedBigraph,1,shading:(True|False)]:=JonesWenzlIdempotent[g,1,shading]=WenzlRecursion[JonesWenzlIdempotent[g,0,shading],0];
*)


(*
JonesWenzlIdempotent[g_GradedBigraph,k_,shading:(True|False)]/;k>1:=
If[LoadIdempotents[g],JonesWenzlIdempotent[g,k,shading],
Module[{result},
JonesWenzlIdempotent[g,k,shading]=result=WenzlRecursion[JonesWenzlIdempotent[g,k-1,shading],With[{\[Delta]=DimensionOfGenerator[g]},ToNumberField0[If[Xor[EvenQ[k],shading],\[Delta],1/\[Delta]] qInteger\[Delta][k-1][\[Delta]]/qInteger\[Delta][k][\[Delta]],NumberFieldGenerator[g]]]];
PackageIdempotents[g];
result
]
]
*)


WenzlRecursion[gpa_GPAMatrix,x_]:=
(DebugPrint["Beginning Wenzl recursion calculation."];
(DebugPrint["   adding strand on right"];AddStrandOnRight[gpa])-x (DebugPrint["   turning down top right"];TurnDownTopRightCorner[gpa]).(DebugPrint["   turning up bottom right"];TurnUpBottomRightCorner[gpa]))


JonesWenzlIdempotent[g_GradedBigraph,0,shading:(True|False)]:=EmptyGPAMatrix[g,shading];
JonesWenzlIdempotent[g_GradedBigraph,1,shading:(True|False)]:=JonesWenzlIdempotent[g,1,shading]=AddStrandOnRight[JonesWenzlIdempotent[g,0,shading]];


JonesWenzlIdempotent[g_GradedBigraph,k_,shading:(True|False)]/;k>1:=
If[LoadIdempotents[g],JonesWenzlIdempotent[g,k,shading],
Module[{result},
JonesWenzlIdempotent[g,k,shading]=result=WenzlRecursion[JonesWenzlIdempotent[g,k-1,shading],With[{\[Delta]=DimensionOfGenerator[g]},ToNumberField0[If[Xor[EvenQ[k],shading],\[Delta],1/\[Delta]] qInteger\[Delta][k-1][\[Delta]]/qInteger\[Delta][k][\[Delta]],NumberFieldGenerator[g]]]];
PackageIdempotents[g];
result
]
]


LoadIdempotents[g_GradedBigraph]:=Module[{},
Off[Get::noopen,Needs::nocont];
Needs["FusionAtlas`Data`"<>GraphToString[g]<>"`Idempotents`"];
On[Get::noopen,Needs::nocont];
LoadIdempotents[g]=False;
True
]


GraphPlanarAlgebraIdempotent[g_GradedBigraph,Vertex[0,1]]:=EmptyGPAElement[g,False]


GraphPlanarAlgebraIdempotent[g_GradedBigraph,Vertex[depth_,k_]]/;depth>0\[And]Length[Cases[Transpose[g[[depth]]],UnitVector[RankAtDepth[g,depth],k]]]>0:=GraphPlanarAlgebraIdempotent[g,Vertex[depth,k]]=
Module[{l,x},
l=FindFirst[Transpose[g[[depth]]],UnitVector[RankAtDepth[g,depth],k]];
(*Work out the lopsided version.*)
x=If[depth>1,RootReduce[g[[depth-1,l]].ReducedDimensionsByDepth[g][[depth-2+1]]/ReducedDimensionsByDepth[g][[depth-1+1,l]]],0];
WenzlRecursion[GraphPlanarAlgebraIdempotent[g,Vertex[depth-1,l]],x]
]


GPAElement/:z_ GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients_List]:=GPAElement[g,shading,{down,up},coefficients/.(loop_List->\[Zeta]_):>(loop->z \[Zeta])]


GPAElement/:GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients1_List]+GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients2_List]:=
With[{dispatch1=Dispatch[coefficients1],dispatch2=Dispatch[coefficients2]},
ExpandCoefficientRules[GPAElement[g,shading,{down,up},{loop_List:>(loop/.dispatch1)+(loop/.dispatch2)}]]
]


GPAElement/:GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients1_List]+GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients2_List]:=
With[{dispatch1=Dispatch[coefficients1],dispatch2=Dispatch[coefficients2]},
GPAElement[g,shading,{down,up},Thread[GraphLoops[g,shading,down+up]->(GraphLoops[g,shading,down+up]/.dispatch1)+(GraphLoops[g,shading,down+up]/.dispatch2)]]
]


ExpandCoefficientRules[GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients_List]]:=Module[{result},
expansionCounter=0;
DebugPrint[".... expanding coefficients"];
result=GPAElement[g,shading,{down,up},
(If[Mod[++expansionCounter,10000]==0,DebugPrint[expansionCounter]];#->Replace[(#/.coefficients),_List->0])&/@GraphLoops[g,shading,down+up]];
DebugPrint[".... finished expanding coefficients"];
result
]


GPACoefficients[GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients_List]]:=Last/@coefficients


VariableGPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},f_]/;up+down>=2:=
ExpandCoefficientRules[GPAElement[g,shading,{down,up},{loop:{PatternSequence[_Vertex,_Integer]..}:>f[loop]}]]


GPAMatrix[GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},values_]]:=
Module[{initialVertices,finalVertices,pathsList,matrices},
initialVertices=If[shading,OddVertices[g],EvenVertices[g]];
finalVertices=If[Xor[shading,OddQ[down]],OddVertices[g],EvenVertices[g]];
pathsList=DeleteCases[Flatten[Outer[{{#1,#2},GraphPathsBetween[g,up,#1,#2],GraphPathsBetween[g,down,#1,#2]}&,initialVertices,finalVertices],1],{{_Vertex,_Vertex},{}}];
matrices=pathsList/.{{initial_Vertex,final_Vertex},upPaths_List,downPaths_List}:>({initial,final}->(Outer[#1~Join~Most[Reverse[Most[#2]]]&,upPaths,downPaths,1]/.values));
GPAMatrix[g,shading,{down,up},matrices]
]


GPAMatrix/:GPAMatrix[g_GradedBigraph,shading:(True|False),{down_,middle_},values1_].GPAMatrix[g_GradedBigraph,shading:(True|False),{middle_,up_},values2_]:=
GPAMatrix[g,shading,{down,up},(#[[1,1]]->#[[2,2]].#[[1,2]])&/@Transpose[{values1,values2}]]


GPAMatrix/:GPAMatrix[g_GradedBigraph,shading:(True|False),bdy_,values1_]+GPAMatrix[g_GradedBigraph,shading:(True|False),bdy_,values2_]:=
GPAMatrix[g,shading,bdy,(#[[1,1]]->#[[1,2]]+#[[2,2]])&/@Transpose[{values1,values2}]]


GPAMatrix/:\[Alpha]_ GPAMatrix[g_GradedBigraph,shading:(True|False),bdy_,values_]:=
GPAMatrix[g,shading,bdy,(#[[1]]->\[Alpha] #[[2]])&/@values]


EmptyGPAMatrix[g_GradedBigraph,shading:(True|False)]:=
GPAMatrix[g,shading,{0,0},If[shading,OddVertices[g],EvenVertices[g]]/.v_Vertex:>({v,v}->{{1}})]


TurnDownTopRightCorner[GPAMatrix[g_GradedBigraph,shading:(True|False),{down_,up_},coefficients_]]:=
Module[{initialVertices,finalVertices,pathsList,matrices,findMatrixEntry},
initialVertices=If[shading,OddVertices[g],EvenVertices[g]];
finalVertices=If[Xor[shading,OddQ[down+1]],OddVertices[g],EvenVertices[g]];

pathsList=DeleteCases[Flatten[Outer[{{#1,#2},GraphPathsBetween[g,up-1,#1,#2],GraphPathsBetween[g,down+1,#1,#2]}&,initialVertices,finalVertices],1],{{_Vertex,_Vertex},{}}];

findMatrixEntry[upPath_,downPath_]:=Module[{s,f,i,j},
s=downPath[[1]];
f=downPath[[-3]];
i=FindFirst[GraphPathsBetween[g,up,s,f],upPath~Join~{downPath[[-2]],downPath[[-3]]}];
j=FindFirst[GraphPathsBetween[g,down,s,f],Drop[downPath,-2]];
({s,f}/.coefficients)[[i,j]]
];

matrices=pathsList/.{{initial_Vertex,final_Vertex},upPaths_List,downPaths_List}:>({initial,final}->Outer[( LopsidedDimension[g,#2[[-1]]]^-1 findMatrixEntry[#1,#2])&,upPaths,downPaths,1]);

GPAMatrix[g,shading,{down+1,up-1},matrices]
]


TurnUpBottomRightCorner[GPAMatrix[g_GradedBigraph,shading:(True|False),{down_,up_},coefficients_]]:=
Module[{initialVertices,finalVertices,pathsList,matrices,findMatrixEntry},
initialVertices=If[shading,OddVertices[g],EvenVertices[g]];
finalVertices=If[Xor[shading,OddQ[up+1]],OddVertices[g],EvenVertices[g]];

pathsList=DeleteCases[Flatten[Outer[{{#1,#2},GraphPathsBetween[g,up+1,#1,#2],GraphPathsBetween[g,down-1,#1,#2]}&,initialVertices,finalVertices],1],{{_Vertex,_Vertex},{}}];

findMatrixEntry[upPath_,downPath_]:=Module[{s,f,i,j},
s=upPath[[1]];
f=upPath[[-3]];
i=FindFirst[GraphPathsBetween[g,up,s,f],Drop[upPath,-2]];
j=FindFirst[GraphPathsBetween[g,down,s,f],downPath~Join~{upPath[[-2]],upPath[[-3]]}];
({s,f}/.coefficients)[[i,j]]
];

matrices=pathsList/.{{initial_Vertex,final_Vertex},upPaths_List,downPaths_List}:>({initial,final}->Outer[(LopsidedDimension[g,#1[[-3]]]findMatrixEntry[#1,#2])&,upPaths,downPaths,1]);

GPAMatrix[g,shading,{down-1,up+1},matrices]
]


AddStrandOnRight[GPAMatrix[g_GradedBigraph,shading:(True|False),{down_,up_},coefficients_]]:=
Module[{initialVertices,finalVertices,pathsList,matrices,findMatrixEntry},
initialVertices=If[shading,OddVertices[g],EvenVertices[g]];
finalVertices=If[Xor[shading,OddQ[down+1]],OddVertices[g],EvenVertices[g]];

pathsList=DeleteCases[Flatten[Outer[{{#1,#2},GraphPathsBetween[g,up+1,#1,#2],GraphPathsBetween[g,down+1,#1,#2]}&,initialVertices,finalVertices],1],{{_Vertex,_Vertex},{}}];

findMatrixEntry[upPath_,downPath_]:=Module[{s,f,i,j},
If[downPath[[-3]]=!=upPath[[-3]],0,
s=downPath[[1]];
f=downPath[[-3]];
i=FindFirst[GraphPathsBetween[g,up,s,f],Drop[upPath,-2]];
j=FindFirst[GraphPathsBetween[g,down,s,f],Drop[downPath,-2]];
({s,f}/.coefficients)[[i,j]]
]
];

matrices=pathsList/.{{initial_Vertex,final_Vertex},upPaths_List,downPaths_List}:>({initial,final}->Outer[findMatrixEntry[#1,#2]&,upPaths,downPaths,1]);

GPAMatrix[g,shading,{down+1,up+1},matrices]
]


GPAElement[GPAMatrix[g_GradedBigraph,shading:(True|False),{down_,up_},values_]]:=
GPAElement[g,shading,{down,up},Flatten[values/.({initial_Vertex,final_Vertex}->matrix_):>
With[{
downPaths=GraphPathsBetween[g,down,initial,final],
upPaths=GraphPathsBetween[g,up,initial,final]
},
If[Length[downPaths]==0\[Or]Length[upPaths]==0,{},
MapThread[Rule,{Outer[#1~Join~Most[Reverse[Most[#2]]]&,upPaths,downPaths,1],matrix},2]]
]
]]


GPACoefficients[GPAMatrix[g_GradedBigraph,shading:(True|False),bdy_,matrices_List]]:=Flatten[Last/@matrices]


SolveTwoTermEquations[conditions_,variablePattern_,variablePreferenceFunction_:(#&)]:=
Module[{constantConditions,oneVariableConditions,twoVariableConditions,sortedConditions,condition,targetVariable,solutions},
sortedConditions=conditions;
While[
DebugPrint["... sorting conditions"];
{constantConditions,oneVariableConditions,twoVariableConditions,sortedConditions}=Flatten/@Reap[Sow[#,Length[Union[Cases[#,variablePattern,{0,\[Infinity]}]]]/.n_/;n>=3:>\[Infinity]]&/@sortedConditions,{0,1,2,\[Infinity]}][[2]];

If[Length[DeleteCases[constantConditions,0]]>0,
DebugPrint["... no solution! This should be impossible; 0 always works."];
DebugPrint[constantConditions]
];
DebugPrint["... found ",Length[oneVariableConditions]," and ", Length[twoVariableConditions]," one- and two-variable conditions."];
Length[oneVariableConditions]+Length[twoVariableConditions]>0,
Cases[Collect[oneVariableConditions,variablePattern,RootReduce[Together[#]]&],x:variablePattern:>(x=0),\[Infinity]];
(condition=Collect[#,variablePattern,Together];
If[condition=!=0,
targetVariable=First[SortBy[Cases[condition,variablePattern,{0,\[Infinity]}],variablePreferenceFunction]];
solutions=Solve[condition==0,targetVariable];
If[Length[solutions]==1,
Collect[solutions[[1]],variablePattern,Together]/.(x_->y_):>(x=y),
Print["Something went wrong:",condition," ",solutions];
];
];)&/@twoVariableConditions;
];
sortedConditions
]


ProgressiveSolve[conditions_,variablePattern_,variablePreferenceFunction_:(#&)]:=
Module[{sortedConditions,n,counter=0,condition,targetVariable,solutions},
sortedConditions=conditions;
sortedConditions=SolveTwoTermEquations[sortedConditions,variablePattern];
(
++counter;
If[Mod[counter,20]==0,DebugPrint[counter,"/",Length[sortedConditions]]];
condition=Collect[#,variablePattern,Together];
If[condition=!=0,
targetVariable=First[SortBy[Cases[condition,variablePattern,{0,\[Infinity]}],variablePreferenceFunction]];
solutions=Solve[condition==0,targetVariable];
If[Length[solutions]==1,
Collect[solutions[[1]],variablePattern,Together]/.(x_->y_):>(x=y),
Print["Something went wrong:",condition," ",solutions];
];
];
)&/@sortedConditions;
{}
]


VertexPreferences[g_GradedBigraph][Vertex[d_,k_]]:=-LopsidedDimension[g,Vertex[d,k]]+0.1 k+0.01d


LoopPreferences[g_][loop_]:=Table[1.+.01k,{k,1,Length[loop]/2}].(VertexPreferences[g]/@loop[[;;;;2]])


LowestWeightEigenvectorConditions[k_GPAElement,\[Lambda]0_]:=
Module[{\[Lambda]=If[NumericQ[\[Lambda]0],ToNumberField[RootReduce[\[Lambda]0],NumberFieldGenerator[k[[1]]]],\[Lambda]0]},
EnlargeNumberField[k[[1]]][\[Lambda]];
(DebugPrint["... preparing first cap"];GPACoefficients[CapAt[k,1]])~Join~(DebugPrint["... preparing second cap"];GPACoefficients[CapAt[k,2]])~Join~(DebugPrint["... preparing rotation"];GPACoefficients[RotateTwoClicks[k]-\[Lambda] k])
]


LowestWeightEigenspace[g_GradedBigraph,\[Lambda]_]:=
If[LoadLowestWeightEigenspaces[g],LowestWeightEigenspace[g,\[Lambda]],
Module[{result,f,kernel,conditions,solution,variables,matrix,c},
DebugPrint["Preparing test kernel..."];
kernel=VariableGPAElement[g,False,{DepthOfBranchPoint[g]+1,DepthOfBranchPoint[g]+1},f];
DebugPrint["Preparing low weight eigenspace conditions..."];
conditions=LowestWeightEigenvectorConditions[kernel,\[Lambda]];
DebugPrint["Beginning progressive solve..."];
solution=ProgressiveSolve[conditions,_f,LoopPreferences[#[[1]]]&];
DebugPrint["Unfolding definitions for f..."];
DownValues[f]=DownValues[f]/.(a_:>b_):>(a:>Evaluate[b]);
DebugPrint["Applying solution rules to test kernel..."];
kernel=kernel/.(x_->y_):>(x->(Collect[y/.solution,_f,Together]));
DebugPrint["Finding remaining variables..."];
variables=Union[Cases[kernel,_f,\[Infinity]]];
DebugPrint["Replacing variables..."];
matrix=GPAMatrix[kernel]/.Thread[variables->Table[c[i],{i,1,Length[variables]}]];
(*DebugPrint["Preparing basis."];
Table[GPAMatrix[kernel]/.Thread[variables->UnitVector[Length[variables],i]],{i,1,Length[variables]}]*)
result=(LowestWeightEigenspace[g,\[Lambda]]=Function[{cc},Evaluate[matrix/.c->cc]]);
DebugPrint["... finished LowestWeightEigenspace"];
PackageLowestWeightEigenspaces[g];
result
]
]


SparseMatrixSolve[conditions_,variablePattern_,variablePreferenceFunction_:(#&)]:=
Module[{variables=SortBy[Union[Cases[conditions,variablePattern,\[Infinity]]],variablePreferenceFunction],matrix,nullspace},
(*matrix=SparseArray[{i_,j_}:>Coefficient[conditions[[i]],variables[[j]]],{Length[conditions],Length[variables]}];*)
DebugPrint["preparing SparseArray"];
matrix=CoefficientArrays[conditions,variables][[2]];
DebugPrint["finding nullspace..."];
nullspace=NullSpace[matrix,Method->"OneStepRowReduction"];
(* TODO, prepare a list of rules *)
Thread[variables->#]&/@nullspace
]


ManyTermSparseMatrixSolve[conditions_,variablePattern_,variablePreferenceFunction_:(#&)]:=
SparseMatrixSolve[SolveTwoTermEquations[conditions,variablePattern,variablePreferenceFunction],variablePattern,variablePreferenceFunction]


CollectGPAMatrix[GPAMatrix[g_,s_,b_,matrices_],p_]:=GPAMatrix[g,s,b,matrices/.((v_->m_):>(v->Collect[m,p,Together]))]


LowestWeightEigenspace[g_GradedBigraph,\[Lambda]_]:=
If[LoadLowestWeightEigenspaces[g],LowestWeightEigenspace[g,\[Lambda]],
Module[{result,f,kernel,conditions,solutions,variables,matrix,c,replacements},
DebugPrint["Preparing test kernel..."];
kernel=VariableGPAElement[g,False,{DepthOfBranchPoint[g]+1,DepthOfBranchPoint[g]+1},f];
DebugPrint["Preparing low weight eigenspace conditions..."];
conditions=LowestWeightEigenvectorConditions[kernel,\[Lambda]];
DebugPrint["Beginning progressive solve..."];
solutions=ManyTermSparseMatrixSolve[conditions,_f,LoopPreferences[#[[1]]]&];
replacements=Thread[Union[Cases[solutions,_f,\[Infinity]]]->Table[c[i],{i,1,Length[solutions]}].(Union[Cases[solutions,_f,\[Infinity]]]/.solutions)];
matrix=GPAMatrix[kernel]/.replacements;
matrix=CollectGPAMatrix[matrix,_c];
If[!FreeQ[matrix,f],
DebugPrint["Oh dear --- the variable hasn't disappeared yet..."];
remainingf=Union[Cases[matrix,_f,\[Infinity]]];
matrix=matrix/.Thread[remainingf->Table[c[i],{i,Length[solutions]+1,Length[solutions]+Length[remainingf]}]];
DebugPrint["But now it's gone..."];
];
result=(LowestWeightEigenspace[g,\[Lambda]]=Function[{cc},Evaluate[matrix/.c->cc]]);
DebugPrint["... finished LowestWeightEigenspace"];
PackageLowestWeightEigenspaces[g];
result
]
]


LoadLowestWeightEigenspaces[g_GradedBigraph]:=Module[{},
Off[Get::noopen,Needs::nocont];
Needs["FusionAtlas`Data`"<>GraphToString[g]<>"`LowestWeightEigenspace`"];
On[Get::noopen,Needs::nocont];
LoadLowestWeightEigenspaces[g]=False;
True
]


End[];


EndPackage[];
