FusionAtlasPaths={"~/projects/fusionatlas/package/"};
$Path=$Path~Join~FusionAtlasPaths;
<<FusionAtlas`

PrepareFusionEquations[g_GradedGraph,Z_]/;RankAtDepth[g,1]==1:=
Module[{},
PrepareAssociativityEquations[g,Z]~Join~PreparePrincipalGraphEquations[g,Z]
]

PrepareFusionEquations[h0_BigraphWithDuals, h1_BigraphWithDuals,Z_,extensions:(True|False):False]/;RankAtDepth[h0,1]==1∧RankAtDepth[h1,1]==1:=
Module[{g,h,equations},
{h[0],h[1]}={h0,h1};
{g[0],g[1]}=If[extensions,
If[EvenQ[GraphDepth[#]],BigraphWithDuals[Append[#[[1]],{Table[1,{RankAtDepth[#[[1]],GraphDepth[#]]}]}],#[[2]]],
BigraphWithDuals[Append[Append[#[[1]],{Table[1,{RankAtDepth[#[[1]],GraphDepth[#]]}]}],{{1}}],Append[#[[2]],{1}]]
]&/@{h0,h1},
{h0,h1}];
ApplyDepthEquations[g[0][[1]],g[1][[1]],Z];
equations=DeleteCases[(PrepareAssociativityEquations[g[0][[1]],g[1][[1]],Z]~Join~PreparePrincipalGraphEquations[g[0][[1]],0,Z]~Join~PreparePrincipalGraphEquations[g[1][[1]],1,Z]~Join~PrepareDualityEquations[g[0],g[1],Z])(*/.{Z[{_,_,d1_,_},{_,_,d2_,_},{_,_,d3_,_}]/;(d3<Abs[d1-d2])∨(d3>d1+d2)->0}*),True](*~Join~PrepareDepthEquations[g[0][[1]],g[1][[1]],Z]*);
If[extensions,
DeleteCases[equations,eq_/;!FreeQ[eq,{L_,_,d_,1}/;d>GraphDepth[h[L]]]],
equations
]
]

Clear[FindPartialFusionRules,FindFusionRules]

FindPartialFusionRules[g_GradedGraph,Z_]/;RankAtDepth[g,1]==1:=
Module[{multiplicities,equations,d1,d2,d3,k1,k2,k3},
equations=SolveAllSolvableEquations[PrepareFusionEquations[g,Z],Z];
multiplicities=Flatten[Table[Z[{d1,k1},{d2,k2},{d3,k3}],
{d1,0,GraphDepth[g]},{k1,1,RankAtDepth[g,d1]},{d2,0,GraphDepth[g]},{k2,1,RankAtDepth[g,d2]},{d3,0,GraphDepth[g]},{k3,1,RankAtDepth[g,d3]}],{{1,2},{3,4},{5,6}}];
PartialFusionRules[g,{{0,0,0}->multiplicities},equations]
]

FindPartialFusionRules[g0_BigraphWithDuals, g1_BigraphWithDuals,Z_,extensions:(True|False):False]/;RankAtDepth[g0,1]==1∧RankAtDepth[g1,1]==1:=
Module[{multiplicities,equations,ZZ,d1,d2,d3,k1,k2,k3,g,L,M,R},
{g[0],g[1]}={g0,g1};
equations=SolveAllSolvableEquations[PrepareFusionEquations[g0,g1,Z,extensions],Z];
ZZ[L_,M_,R_]:=ZZ[L,M,R]=Flatten[Table[Z[{L,M,d1,k1},{M,R,d2,k2},{L,R,d3,k3}],{d1,Mod[M+L,2],GraphDepth[g[L]],2},{k1,1,RankAtDepth[g[L],d1]},{d2,Mod[M+R,2],GraphDepth[g[M]],2},{k2,1,RankAtDepth[g[M],d2]},{d3,Mod[L+R,2],GraphDepth[g[L]],2},{k3,1,RankAtDepth[g[L],d3]}],{{1,2},{3,4},{5,6}}];
multiplicities=Flatten[Table[{L,M,R}->ZZ[L,M,R],{L,0,1},{M,0,1},{R,0,1}]];
PartialFusionRules[{g0,g1},multiplicities,equations]
]

SolvePartialFusionRules[PartialFusionRules[g_,multiplicities_,equations_],Z_]:=
FusionRules[g,multiplicities]/.Solve[equations,Union[Cases[equations,_Z,∞]]]

FindFusionRules[g_GradedGraph,Z_]:=
Module[{partialFusionRules=FindPartialFusionRules[g,Z],possibleDimensionList,possibleValues,vars,dim,solutions,numberCases},
DebugPrint["Completed partial fusion rules."];
dim[d_,k_]:=DimensionsByDepth[g][[d+1,k]];
vars=Union[Cases[partialFusionRules,_Z,∞]];
possibleDimensionList=(vars/.(a_Z:>{a,Range[0,Floor[(dim@@a[[1]])*(dim@@a[[2]])/(dim@@a[[3]])]]}));
numberCases=Times@@(1+Length/@(Last/@possibleDimensionList));
DebugPrint["Trying ",numberCases," solutions."];
If[numberCases>1000000,Abort[]];
possibleValues=Flatten[Outer[List,Sequence @@(Last/@possibleDimensionList)],Length[possibleDimensionList]-1];
solutions=Select[possibleValues,(And@@partialFusionRules[[3]]/.Thread[vars->#])&];
FusionRules[g,#]&/@(partialFusionRules[[2]]/.Thread[vars->#]&/@solutions)
]

FindFusionRules[g0_BigraphWithDuals, g1_BigraphWithDuals,Z_]:=
Module[{partialFusionRules=FindPartialFusionRules[g0,g1,Z],possibleDimensionList,possibleValues,vars,dim,solutions,numberCases,g,conditions},
{g[0],g[1]}={g0,g1};
DebugPrint["Completed partial fusion rules."];
If[MemberQ[partialFusionRules[[3]],False],Return[{}]];
dim[a_,b_,d_,k_]:=ReducedDimensionsByDepth[g[a][[1]]][[d+1,k]];
vars=Union[Cases[partialFusionRules,_Z,∞]];
conditions=And@@(partialFusionRules[[3]]~Join~Thread[Flatten[partialFusionRules[[2,All,2]]]≥0]);
possibleDimensionList=(vars/.(a_Z:>{a,Range[0,Floor[(dim@@a[[1]])*(dim@@a[[2]])/(dim@@a[[3]])]]}));
numberCases=Times@@(1+Length/@(Last/@possibleDimensionList));
DebugPrint["Trying ",numberCases," solutions."];
If[numberCases>1000000,
DebugPrint[partialFusionRules];
DebugPrint[vars];
DebugPrint[conditions];
DebugPrint[possibleDimensionList];
];
If[Length[vars]==0,
solutions={{}},
possibleValues=Flatten[Outer[List,Sequence @@(Last/@possibleDimensionList)],Length[possibleDimensionList]-1];
solutions=Select[possibleValues,(conditions/.Thread[vars->#])&];
];
FusionRules[{g0,g1},#]&/@(partialFusionRules[[2]]/.Thread[vars->#]&/@solutions)
]

FindFusionRules[g_GradedGraph]:=FindFusionRules[g]=Module[{Z},FindFusionRules[g,Z]]

FindFusionRules[g0_BigraphWithDuals, g1_BigraphWithDuals]:=FindFusionRules[g0, g1]=Module[{Z},FindFusionRules[g0,g1,Z]]

PrepareAssociativityEquations[g_GradedGraph,Z_]:=Module[{ZZ,d1,d2,d3,k1,k2,k3},
ZZ=Flatten[Table[Z[{d1,k1},{d2,k2},{d3,k3}],
{d1,0,GraphDepth[g]},{k1,1,RankAtDepth[g,d1]},{d2,0,GraphDepth[g]},{k2,1,RankAtDepth[g,d2]},{d3,0,GraphDepth[g]},{k3,1,RankAtDepth[g,d3]}],{{1,2},{3,4},{5,6}}];

DeleteCases[Union[#==0&/@Flatten[ZZ.ZZ-Transpose[ZZ.Transpose[ZZ,{2,1,3}],{2,3,1,4}]]],True]
]

ApplyDepthEquations[g0_GradedBigraph,g1_GradedBigraph,Z_]:=
Module[{g},
{g[0],g[1]}={g0,g1};
Cases[Flatten[Table[Z[{L,M,d1,k1},{M,R,d2,k2},{L,R,d3,k3}],
{L,0,1},{M,0,1},{R,0,1},
{d1,Mod[M+L,2],GraphDepth[g[L]],2},{k1,1,RankAtDepth[g[L],d1]},{d2,Mod[M+R,2],GraphDepth[g[M]],2},{k2,1,RankAtDepth[g[M],d2]},{d3,Mod[L+R,2],GraphDepth[g[L]],2},{k3,1,RankAtDepth[g[L],d3]}]],(z:(Z[{_,_,d1_,_},{_,_,d2_,_},{_,_,d3_,_}]/;(d3<Abs[d1-d2])∨(d3>d1+d2))):>(z=0)];
]

PrepareAssociativityEquations[g0_GradedBigraph,g1_GradedBigraph,Z_]:=Module[{g,ZZ,associativityEquations,d1,d2,d3,k1,k2,k3,A,B,C,D},
DebugPrint["beginning preparation of associativity equations"];
{g[0],g[1]}={g0,g1};
ZZ[L_,M_,R_]:=ZZ[L,M,R]=Flatten[Table[Z[{L,M,d1,k1},{M,R,d2,k2},{L,R,d3,k3}],{d1,Mod[M+L,2],GraphDepth[g[L]],2},{k1,1,RankAtDepth[g[L],d1]},{d2,Mod[M+R,2],GraphDepth[g[M]],2},{k2,1,RankAtDepth[g[M],d2]},{d3,Mod[L+R,2],GraphDepth[g[L]],2},{k3,1,RankAtDepth[g[L],d3]}],{{1,2},{3,4},{5,6}}];
associativityEquations[A_,B_,C_,D_]:=
DeleteCases[Union[#==0&/@Flatten[ZZ[A,B,C].ZZ[A,C,D]-Transpose[ZZ[B,C,D].Transpose[ZZ[A,B,D],{2,1,3}],{2,3,1,4}]]],True];
Flatten[Table[DebugPrint[A,B,C,D];associativityEquations[A,B,C,D],{A,0,1},{B,0,1},{C,0,1},{D,0,1}]]
]

PreparePrincipalGraphEquations[g_GradedGraph,Z_]:=
Module[{tensorWithX,tensorWith1,d1,k1,d2,k2},

tensorWithX[d1_,k1_]:=
Flatten[Table[
Z[{d1,k1},{1,1},{d2,k2}]-
Switch[d2-d1,
0,g[[1,d1+1,k2,k1]],
1,g[[2,d1+1,k2,k1]],
-1,g[[2,d1,k1,k2]],
_,0]==0,
{d2,0,GraphDepth[g]},{k2,1,RankAtDepth[g,d2]}]];

tensorWith1[d1_,k1_]:=
Flatten[Table[{
Z[{d1,k1},{0,1},{d2,k2}]-If[d1==d2∧k1==k2,1,0]==0,
Z[{0,1},{d1,k1},{d2,k2}]-If[d1==d2∧k1==k2,1,0]==0
},
{d2,0,GraphDepth[g],1},{k2,1,RankAtDepth[g,d2]}]];


Flatten[Table[tensorWithX[d1,k1]~Join~tensorWith1[d1,k1],{d1,0,GraphDepth[g]},{k1,1,RankAtDepth[g,d1]}]]
]

PreparePrincipalGraphEquations[g_GradedBigraph,L:(0|1),Z_]:=
Module[{tensorWithX,tensorWith1,d1,k1,d2,k2,R},
DebugPrint["beginning preparation of principal graph ",L, " equations"];
tensorWithX[L,R_,d1_,k1_]:=
Flatten[Table[
Z[{L,R,d1,k1},{R,1-R,1,1},{L,1-R,d2,k2}]-
Switch[d2-d1,
1,g[[d1+1,k2,k1]],
-1,g[[d1,k1,k2]],
_,0]==0,
{d2,1-Mod[d1,2],GraphDepth[g],2},{k2,1,RankAtDepth[g,d2]}]];

tensorWith1[L,R_,d1_,k1_]:=
Flatten[Table[{
Z[{L,R,d1,k1},{R,R,0,1},{L,R,d2,k2}]-If[d1==d2∧k1==k2,1,0]==0,
Z[{L,L,0,1},{L,R,d1,k1},{L,R,d2,k2}]-If[d1==d2∧k1==k2,1,0]==0
},
{d2,Mod[d1,2],GraphDepth[g],2},{k2,1,RankAtDepth[g,d2]}]];


Flatten[Table[tensorWithX[L,R,d1,k1]~Join~tensorWith1[L,R,d1,k1],{R,0,1},{d1,0,GraphDepth[g]},{k1,1,RankAtDepth[g,d1]}]]
]

PrepareDualityEquations[g0_BigraphWithDuals, g1_BigraphWithDuals,Z_]:=
Module[{g,dual,dualityEquations,d1,k1,d2,k2,d3,k3,A,B,C},
DebugPrint["beginning preparation of duality equations"];

{g[0],g[1]}={g0,g1};

dual[{A_,A_,d_/;EvenQ[d],k_}]:={A,A,d,g[A][[2,d/2+1,k]]};
dual[{A_,B_,d_/;OddQ[d],k_}]:={B,A,d,k};

dualityEquations[L_,M_,R_]:=Flatten[Table[Z[{L,M,d1,k1},{M,R,d2,k2},{L,R,d3,k3}]-Z[dual[{L,R,d3,k3}],{L,M,d1,k1},dual[{M,R,d2,k2}]]==0,{d1,Mod[M+L,2],GraphDepth[g[L]],2},{k1,1,RankAtDepth[g[L],d1]},{d2,Mod[M+R,2],GraphDepth[g[M]],2},{k2,1,RankAtDepth[g[M],d2]},{d3,Mod[L+R,2],GraphDepth[g[L]],2},{k3,1,RankAtDepth[g[L],d3]}],{{1,2},{3,4},{5,6}}];

Flatten[Table[dualityEquations[A,B,C],{A,0,1},{B,0,1},{C,0,1}]]
]

SolvablePatterns[Z_]:=SolvablePatterns[Z]={{1,_Z==0,SolveSolvableEquations[Z]},{1,-_Z==0,SolveSolvableEquations[Z]},{1,_Integer+_Z==0,SolveSolvableEquations[Z]},{1,_Integer-_Z==0,SolveSolvableEquations[Z]},{2,a_Z+b_Z==0,AllVariablesAreZero[Z]},{2,-a_Z-b_Z==0,AllVariablesAreZero[Z]},{2,a_Z-b_Z==0,SolveSolvableEquations[Z]},{∞,e_/;(SimplifyingEquations=True;Length[CoefficientArrays[e[[1]],Union[Cases[e,_Z,∞]]]]==2),SolveSolvableEquations[Z]}};

SolveSolvableEquations[Z_][F_]:=
Module[{solutions},
solutions=Solve[F,Union[Cases[F,_Z,∞]]];
Switch[Length[solutions],
1,
DebugPrint["solved for ",Length[solutions[[1]]], " variables."];
solutions[[1]]/.(a_->b_):>(a=b);
True,
0,
False,
_,DebugPrint["oh dear, multiple solutions to linear equations!"];True
]
]

AllVariablesAreZero[Z_][eq_]:=(Cases[eq,_Z,∞]/.a_Z:>(a=0);True)

SolveSomeSolvableEquations[eq_,Z_]:=Module[{equationsWithKVariables,targetEquations,i,numberOfVariables,pattern,solveMethod},
equationsWithKVariables[k_]:=equationsWithKVariables[k]=Cases[eq,e_/;Length[Union[Cases[e,_Z,∞]]]≤k];
For[i=1,i≤Length[SolvablePatterns[Z]],++i,
DebugPrint["Trying pattern ",i];
{numberOfVariables,pattern,solveMethod}=SolvablePatterns[Z][[i]];
targetEquations=Cases[equationsWithKVariables[numberOfVariables],pattern];
If[Length[targetEquations]>0,
If[solveMethod[targetEquations],
DebugPrint["found ", Length[targetEquations], " equations"];
Return[Length[targetEquations]],
DebugPrint["found a contradiction"];
Return[-1]
]
];
];

]

SolveAllSolvableEquations[eq0_,Z_]:=
Module[{eq=eq0,r=0},
SimplifyingEquations=False;
While[(!MemberQ[eq,False])∧(r=SolveSomeSolvableEquations[eq,Z])>0,
If[SimplifyingEquations,eq=Collect[eq,_Z,Together]];
eq=DeleteCases[Union[eq],True]
];
If[r==-1,PrependTo[eq,False]];
DeleteCases[Union[eq],True]
]

EvenPart[FusionRules[{g0_,g1_},multiplicities_]]:=FusionRules[EvenPart[g0[[1]]],{{0,0,0}->({0,0,0}/.multiplicities)}]

FlipFusionRules[FusionRules[{g0_,g1_},multiplicities_]]:=FusionRules[{g1,g0},multiplicities/.({a_,b_,c_}->T_):>((1-{a,b,c})->T)]