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`"];

Graph combinatorics

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}]

Lopsided dimensions

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[{δ=ReducedDimensionOfGenerator[g]},
ToNumberField0[RootReduce[If[OddQ[d],(1/δ),1]ReducedDimensionsByDepth[g][[d+1,k]]],NumberFieldGenerator[g]]
]

Graph planar algebra

Planar algebra operations

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

GPA element takes : a graph, whether or not the left-hand region is shaded, {a,b} meaning a strands going up (?) and b going down(?), and rules assigning any length a+b loop to a number.  Length (a+b) loops are denoted as a sequence of vertices and integers (the integer specifies which edge between the vertices to use).

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___}->ζ_):>({b,nb,z,a,na}->ζ((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___}->ζ_):>({c,nc,z,a,na,b,nb}->ζ((LopsidedDimension[g,a])/(LopsidedDimension[g,c])))]

RotateOneClick[GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients_List]/;up≥2∧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}->ζ_):>({e,ne,z,nz,A,b,nb,D}->ζ((LopsidedDimension[g,z])/(LopsidedDimension[g,b])))]

RotateTwoClicks[GPAElement[g_GradedBigraph,shading:(True|False),{down_Integer,up_Integer},coefficients_List]/;up≥3∧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}->ζ_):>({e,ne,f,nf,z,nz,A,b,nb,c,nc,D}->ζ((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∧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}):∞]:=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___}->ζ_):>({A,v,B}->ζ 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___}->ζ_):>({A,v,B}->ζ (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}->ζ_):>({b,nb,a,na,B}->ζ 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___}->ζ_):>({b,nb,B,a,na}->ζ (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 ∧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δ[k_]:=qIntegerδ[k]=Function[{δ},Evaluate[Collect[δ qIntegerδ[k-1][δ]-qIntegerδ[k-2][δ],δ]]]
qIntegerδ[2]=#&;
qIntegerδ[1]=1&;
qIntegerδ[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[{δ=DimensionOfGenerator[g]},ToNumberField0[If[Xor[EvenQ[k],shading],δ,(1/δ)]((qIntegerδ[k-1][δ])/(qIntegerδ[k][δ])),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[{δ=DimensionOfGenerator[g]},ToNumberField0[If[Xor[EvenQ[k],shading],δ,(1/δ)]((qIntegerδ[k-1][δ])/(qIntegerδ[k][δ])),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∧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]
]

Arithmetic

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

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)]]
]

Technical details

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]}]]

Matrices

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/:α_ GPAMatrix[g_GradedBigraph,shading:(True|False),bdy_,values_]:=
GPAMatrix[g,shading,bdy,(#[[1]]->α #[[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∨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]

Lowest weight vectors

Dimension here is really just a stand - in for distance from the trivalent vertex, but allows an easy generalisation to other graphs.

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,∞}]]]/.n_/;n≥3:>∞]&/@sortedConditions,{0,1,2,∞}][[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),∞];
(condition=Collect[#,variablePattern,Together];
If[condition=!=0,
targetVariable=First[SortBy[Cases[condition,variablePattern,{0,∞}],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,∞}],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,λ0_]:=
Module[{λ=If[NumericQ[λ0],ToNumberField[RootReduce[λ0],NumberFieldGenerator[k[[1]]]],λ0]},
EnlargeNumberField[k[[1]]][λ];
(DebugPrint["... preparing first cap"];GPACoefficients[CapAt[k,1]])~Join~(DebugPrint["... preparing second cap"];GPACoefficients[CapAt[k,2]])~Join~(DebugPrint["... preparing rotation"];GPACoefficients[RotateTwoClicks[k]-λ k])
]

LowestWeightEigenspace[g_GradedBigraph,λ_]:=
If[LoadLowestWeightEigenspaces[g],LowestWeightEigenspace[g,λ],
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,λ];
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,∞]];
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,λ]=Function[{cc},Evaluate[matrix/.c->cc]]);
DebugPrint["... finished LowestWeightEigenspace"];
PackageLowestWeightEigenspaces[g];
result
]
]

SparseMatrixSolve[conditions_,variablePattern_,variablePreferenceFunction_:(#&)]:=
Module[{variables=SortBy[Union[Cases[conditions,variablePattern,∞]],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,λ_]:=
If[LoadLowestWeightEigenspaces[g],LowestWeightEigenspace[g,λ],
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,λ];
DebugPrint["Beginning progressive solve..."];
solutions=ManyTermSparseMatrixSolve[conditions,_f,LoopPreferences[#[[1]]]&];
replacements=Thread[Union[Cases[solutions,_f,∞]]->Table[c[i],{i,1,Length[solutions]}].(Union[Cases[solutions,_f,∞]]/.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,∞]];
matrix=matrix/.Thread[remainingf->Table[c[i],{i,Length[solutions]+1,Length[solutions]+Length[remainingf]}]];
DebugPrint["But now it's gone..."];
];
result=(LowestWeightEigenspace[g,λ]=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[];