Needs["LinAlg`"]

Print["Loading CCReduce..."]

$MaxDegree=Infinity

Clear[MultiLinear]
MultiLinear[op_]:=(op[lft___,0,rgt___]=0;
        op[lft___,x_?NumberQ * y_,rgt___]:=x*op[lft,y,rgt];
        op[lft___,x_Symbol * y_,rgt___]:=x*op[lft,y,rgt];
        op[lft___,x_+y_,rgt___]:=op[lft,x,rgt]+op[lft,y,rgt])

SortingRuleForType[type_] := SortingRuleForType[type] = Module[
  {Cs,i},
  Print["Now sorting ",type," ..."];
  Cs=Table[{-Coefficient[type,ColoredUV[i]],ColoredUV[i]},{i,1,16}];
  Cs=Sort[DeleteCases[Cs,{0,_}]];
  Thread[Rule[Last /@ Cs, ColoredUV /@ Range[Length[Cs]]]]
]
SortingBackRuleForType[type_] := SortingBackRuleForType[type] =
  Reverse /@ SortingRuleForType[type]

Type2Class[type_] := Type2Class[type] = Module[
  {OVs,Cs,rtype,n,k},
  rtype=type-(OVs=Coefficient[type,OV])*OV;
  For[Cs={};n=1, rtype=!=0 && n<=16, ++n,
    AppendTo[Cs,k=Coefficient[rtype,ColoredUV[n]]];
    rtype-=k*ColoredUV[n]
  ];
  {(OVs+Plus @@ Cs)/2, Cs}
]

Type2FileName[type_] := Type2FileName[type] = Module[
  {m,cols},
  {m,cols}=Type2Class[type];
  ToString[m]~StringJoin~(
    StringJoin @@ (ToString[StringForm["-``",#]]& /@ cols)
  )
]

ToCCNames[expr_] := Expand[CCReduce[expr] /. g_Graph :> NameOfGraph[g]]

$MissingNames={}
NameOfGraph[g_Graph] := Block[
  {type,rule},
  rule=SortingRuleForType[type=GraphType[g]];
  type=type /. rule;
  If[Head[BasisForType[type]]===BasisForType, (
    If[!MemberQ[$MissingNames,type],
      AppendTo[$MissingNames,type];
      Print["No names for type ",type]
    ];
    g
  ),
    NameOfGraph[g] = Dot[
      (g /. rule)~InTermsOf~CCReduce[Last /@ BasisForType[type]],
      First /@ BasisForType[type]
    ] /. (Reverse /@ rule)
  ]
]

CCReduce[expr_] := Profile[CCReduce,
  Expand[(expr /. g_Graph :> Times @@ SplitGraph[g])
    /. g_Graph :> Block[
      {c,ng,type,rule},
      {c,ng}=CanonicalForm[g /. (rule=SortingRuleForType[GraphType[g]])];
      If[c===0, 0,
        type=GraphType[ng];
        If[Head[GraphsOfType[type]] =!= List,Module[{fname,graphs,rels,s},
          fname="CCReduceData/"~StringJoin~Type2FileName[type];
          If[MemberQ[FileNames["*","CCReduceData"],fname],(
            Print["Loading ",fname," ..."];
            {graphs,rels}=Get[fname];
            s=NewSubspace[Length[graphs]];
            Scan[AddVector[s,#]&, rels];
          ),
            {graphs,s}=IHX[ng,TryUOVs -> True];
            graphs=UncompressGraph /@ graphs;
            Print["Writing ",fname," ..."];
            Put[{graphs,SubspaceBasis[s]},fname]
          ];
          {GraphsOfType[type],RelationsForType[type]}={graphs,s}
        ]];
        at=Position[GraphsOfType[type],ng,{1},Heads->False];
        If[at==={},
          Print["**** CannotCCReduce !!! ****"]; CannotCCReduce,
          at=at[[1,1]];
          c=Expand[c*ReduceMod[RelationsForType[type],BasisVector[at]]];
          c /. BasisVector[at_] :> GraphsOfType[type][[at]] /. (
            Reverse /@ rule
          )
        ]
      ]
    ]
  ]
]

Fuse[g_Graph,l1_List,l2_List] := Profile[Fuse,Module[
  {ng,i,del,redirect,rg},
  ng=g; redirect={};
  rg[i_] := ng[[i,1]] //. redirect;
  Do[
    v1=l1[[i]]; v2=l2[[i]];
    If[Length[ng[[rg[v1]]]]==1 || Length[ng[[rg[v2]]]]==1,
      AppendTo[redirect,v1->rg[v2]]; AppendTo[redirect,v2->rg[v1]],
      ng=ReplacePart[ng,Dashed[rg[v1],rg[v2]],v1];
      AppendTo[redirect,v2->v1]
    ],
    {i,Length[l1]}
  ];
  ng=Delete[ng,List /@ (del=First /@ redirect)] //. redirect;
  ng /. Thread[Rule[Complement[Range[Length[g]],del],Range[Length[ng]]]]
]]

Fuse[g1_Graph,g2_Graph,l1_List,l2_List] := Fuse[
  MergeGraphs[GraphList[g1,g2]],
  l1, Length[g1]+l2
]

MultiLinear[LeftStableFusion]
LeftStableFusion[p_Times,g_Graph] := LeftStableFusion[MergeGraphs[p],g]
LeftStableFusion[g_Graph,p_Times] := LeftStableFusion[g,MergeGraphs[p]]
LeftStableFusion[g1_Graph,g2_Graph,vts__List] := Profile[LeftStableFusion,
Module[
  {p,i,f,h},
  p=Map[
    {
      Flatten[First /@ Position[g1,First[#]]],
      Flatten[First /@ Position[g2,Last[#]]]
    }&,
    {vts}
  ];
  If[And @@ Apply[LessEqual, Map[Length,p,{2}], {1}],
    p=Table[
      {p[[i,1]],#}& /@ Permutations[Length[p[[i,1]]],p[[i,2]]],
      {i,Length[p]}
    ];
    p=Flatten[Outer @@ ({f}~Join~Apply[h,p,{1}])];
    p=p /. x_f :> (f @@ MapThread[Join,List @@ x]);
    p=Plus @@ (p /. f[x_,y_] :> CCReduce[Fuse[g1,g2,x,y]]),
    0
  ]
]]

BCH8=X + Y + b[X, Y]/2 + b[X, b[X, Y]]/12 - b[X, b[X, b[X, b[X, Y]]]]/720 + 
  b[X, b[X, b[X, b[X, b[X, b[X, Y]]]]]]/30240 - 
  b[X, b[X, b[X, b[X, b[X, b[Y, b[X, Y]]]]]]]/20160 - 
  b[X, b[X, b[X, b[X, b[Y, b[X, Y]]]]]]/10080 + 
  b[X, b[X, b[X, b[X, b[Y, b[Y, b[X, Y]]]]]]]/20160 - 
  b[X, b[X, b[X, b[Y, b[X, Y]]]]]/1440 + 
  b[X, b[X, b[X, b[Y, b[X, b[X, b[X, Y]]]]]]]/12096 - 
  b[X, b[X, b[X, b[Y, b[X, b[Y, b[X, Y]]]]]]]/5040 - 
  b[X, b[X, b[X, b[Y, b[Y, b[X, Y]]]]]]/7560 - 
  (5*b[X, b[X, b[X, b[Y, b[Y, b[Y, b[X, Y]]]]]]])/24192 - 
  b[X, b[X, b[Y, b[X, Y]]]]/120 + 
  b[X, b[X, b[Y, b[X, b[X, b[X, Y]]]]]]/2520 + 
  b[X, b[X, b[Y, b[X, b[X, b[Y, b[X, Y]]]]]]]/3360 + 
  b[X, b[X, b[Y, b[X, b[Y, b[X, Y]]]]]]/3360 + 
  b[X, b[X, b[Y, b[X, b[Y, b[Y, b[X, Y]]]]]]]/2520 - 
  b[X, b[X, b[Y, b[Y, b[X, Y]]]]]/720 - 
  b[X, b[X, b[Y, b[Y, b[Y, b[X, Y]]]]]]/3360 + 
  b[X, b[X, b[Y, b[Y, b[Y, b[Y, b[X, Y]]]]]]]/20160 - 
  b[X, b[Y, b[X, Y]]]/24 + b[X, b[Y, b[X, b[X, b[X, Y]]]]]/720 - 
  b[X, b[Y, b[X, b[X, b[X, b[X, b[X, Y]]]]]]]/20160 - 
  b[X, b[Y, b[X, b[X, b[X, b[Y, b[X, Y]]]]]]]/6720 + 
  b[X, b[Y, b[X, b[X, b[Y, b[X, Y]]]]]]/1680 + 
  b[X, b[Y, b[X, b[X, b[Y, b[Y, b[X, Y]]]]]]]/15120 + 
  b[X, b[Y, b[X, b[Y, b[X, Y]]]]]/240 - 
  b[X, b[Y, b[X, b[Y, b[X, b[X, b[X, Y]]]]]]]/10080 - 
  b[X, b[Y, b[X, b[Y, b[X, b[Y, b[X, Y]]]]]]]/2016 + 
  b[X, b[Y, b[X, b[Y, b[Y, b[X, Y]]]]]]/1260 - 
  b[X, b[Y, b[X, b[Y, b[Y, b[Y, b[X, Y]]]]]]]/20160 - 
  b[X, b[Y, b[Y, b[X, Y]]]]/360 + 
  b[X, b[Y, b[Y, b[X, b[X, b[Y, b[X, Y]]]]]]]/20160 - 
  b[X, b[Y, b[Y, b[X, b[Y, b[Y, b[X, Y]]]]]]]/10080 + 
  b[X, b[Y, b[Y, b[Y, b[X, Y]]]]]/1440 + 
  b[X, b[Y, b[Y, b[Y, b[Y, b[X, Y]]]]]]/10080 - 
  b[X, b[Y, b[Y, b[Y, b[Y, b[Y, b[X, Y]]]]]]]/60480 - b[Y, b[X, Y]]/12 + 
  b[Y, b[X, b[X, b[X, Y]]]]/360 - 
  b[Y, b[X, b[X, b[X, b[X, b[X, Y]]]]]]/10080 - 
  b[Y, b[X, b[X, b[X, b[Y, b[X, Y]]]]]]/3360 + 
  b[Y, b[X, b[X, b[Y, b[Y, b[X, Y]]]]]]/7560 + 
  b[Y, b[X, b[Y, b[X, Y]]]]/120 - 
  b[Y, b[X, b[Y, b[X, b[X, b[X, Y]]]]]]/5040 - 
  b[Y, b[X, b[Y, b[X, b[Y, b[X, Y]]]]]]/1008 - 
  b[Y, b[X, b[Y, b[Y, b[Y, b[X, Y]]]]]]/10080 + 
  b[Y, b[Y, b[X, b[X, b[Y, b[X, Y]]]]]]/10080 - 
  b[Y, b[Y, b[X, b[Y, b[Y, b[X, Y]]]]]]/5040 + 
  b[Y, b[Y, b[Y, b[X, Y]]]]/720 - b[Y, b[Y, b[Y, b[Y, b[Y, b[X, Y]]]]]]/30240

BCH[md_Integer,xd_Integer,yd_Integer] := BCH[md,xd,yd] = Dot[
  (List @@ BCH8),
  Block[
    {d,d1,d2},
    ((List @@ BCH8) //. {
      X:>d[xd], Y:> d[yd],
      b[d[d1_],d[d2_]] :> d[d1+d2],
      c_?NumberQ*d[d1_] :> d[d1]
    }) /. d[d1_] :> If[d1>md,0,1]
  ]
]

ExpBCH[md_Integer,xd_Integer,yd_Integer] := ExpBCH[md,xd,yd] = Module[
  {exp,d,d1,d2,n},
  exp=Expand[Sum[(BCH8-X-Y)^n/n!,{n,0,md-xd-yd}]];
  Dot[
    (List @@ exp),
    ((List @@ exp) //. {
      X:>d[xd], Y:> d[yd],
      b[d[d1_],d[d2_]] :> d[d1+d2],
      d[d1_]*d[d2_] :> d[d1+d2],
      d[d1_]^n_ :> d[n*d1],
      c_?NumberQ*d[d1_] :> d[d1]
    }) /. d[d1_] :> If[d1>md,0,1]
  ]
]

ExpBCH[0,_]=ExpBCH[_,0]=1
ExpBCH[xd_Integer,yd_Integer] := ExpBCH[xd,yd] = Module[
  {exp,n,mon},
  Print["Computing ExpBCH[",xd,",",yd,"]..."];
  exp=Expand[Sum[(BCH8-X-Y)^n/n!,{n,0,Min[xd,yd]}]];
  Select[exp,(
    mon=#//.b[x_,y_]:>x*y;
    Exponent[mon,X]<=xd && Exponent[mon,Y]<=yd
  )&]
]

Lie2Graph[z_Symbol,b[x_Symbol,y_Symbol]] := Graph[z[4],x[4],y[4],OV[1,2,3]]
Lie2Graph[z_Symbol,b[x_Symbol,y_b]] := Module[{T},
  Fuse[Graph[z[4],x[4],T[4],OV[1,2,3]], Lie2Graph[T,y], {3}, {1}]
]
Lie2Graph[z_Symbol,b[x_b,y_Symbol]] := Module[{T},
  Fuse[Graph[z[4],T[4],y[4],OV[1,2,3]], Lie2Graph[T,x], {2}, {1}]
]
Lie2Graph[z_Symbol,b[x_b,y_b]] := Module[{Tx,Ty},
  Fuse[
    Fuse[Graph[z[4],Tx[4],Ty[4],OV[1,2,3]], Lie2Graph[Tx,x], {2}, {1}],
    Lie2Graph[Ty,y],
    {3}, {1}
  ]
]
Lie2Graph[z_Symbol,expr_] := expr /. x_b :> Lie2Graph[z,x]

FusionBCH[xd_Integer,yd_Integer,z_Symbol] := FusionBCH[xd,yd,z] =
  Graph[]-1+(
    MergeGraphs[Lie2Graph[z,ExpBCH[xd,yd]]] /.
      g_Graph :> (Times @@ CanonicalForm[g])
  )

AMergeColors[g_Graph,x_Symbol,y_Symbol,z_Symbol] := Profile[AMergeColors,
  CCReduce[LeftStableFusion[
    FusionBCH[Count[g,_x],Count[g,_y],z],
    g, {X,x}, {Y,y}
  ] /. {x->z, y->z}]
]
AMergeColors[expr_,x_Symbol,y_Symbol,z_Symbol] := Expand[
  MergeGraphs[Expand[expr]] /. g_Graph :> AMergeColors[g,x,y,z]
]

ExpBCH[md_Integer,xd_Integer,yd_Integer] := ExpBCH[md,xd,yd] = Module[
  {exp,d,d1,d2,n},
  exp=Expand[Sum[(BCH8-X-Y)^n/n!,{n,0,md-xd-yd}]];
  Dot[
    (List @@ exp),
    ((List @@ exp) //. {
      X:>d[xd], Y:> d[yd],
      b[d[d1_],d[d2_]] :> d[d1+d2],
      d[d1_]*d[d2_] :> d[d1+d2],
      d[d1_]^n_ :> d[n*d1],
      c_?NumberQ*d[d1_] :> d[d1]
    }) /. d[d1_] :> If[d1>md,0,1]
  ]
]

ModByDegree[deg_,expr_] := Module[
  {h,d2},
  d2=2deg;
  Expand[expr] /. g_Graph :> h^Count[g, _[_] | _[_,_,_]]*g
    /. h^(p_:1) :> If[p<d2,1,0]
]

LoopNumber[g_Graph] := (Count[g,_OV]+2-Count[g,_[_]])/2

A3Multiply[x_,y_] := Profile[A3Multiply,Module[
  {res},
  res=ModByDegree[$MaxDegree+1, x*(y /. {C1->C4, C2->C5, C3->C6})];
  res=AMergeColors[res,C1,C4,C1];
  res=AMergeColors[res,C2,C5,C2];
  AMergeColors[res,C3,C6,C3]
]]

A3Multiply[expr_] := Expand[
  expr /. prod_p :> Fold[CCReduce[A3Multiply[#1,#2]]&,First[prod],Rest[prod]]
]

A4Multiply[x_,y_] := Module[
  {res},
  res=ModByDegree[$MaxDegree, x*(y /. {C1->C5, C2->C6, C3->C7, C4->C8})];
  res=AMergeColors[g,C1,C5,C1];
  res=AMergeColors[g,C2,C6,C2];
  res=AMergeColors[g,C3,C7,C3];
  AMergeColors[g,C4,C8,C4]
]

A4Multiply[expr_] := Expand[
  expr /. prod_p :> Fold[CCReduce[A4Multiply[#1,#2]]&,First[prod],Rest[prod]]
]

BExp[expr_] := Module[
  {term=1,sum=0,n=0},
  While[term=!=0,
    sum+=Expand[term/n!];
    ++n;
    term=ModByDegree[$MaxDegree+1,Expand[term*expr]]
  ];
  sum
]

A3Exp[expr_] := Module[
  {term=1,sum=0,n=0},
  While[term=!=0,
    sum+=Expand[term/n!];
    ++n;
    term=CCReduce[term~A3Multiply~expr]
  ];
  sum
]

BLog[1]=0
BLog[1+expr_] := Module[{term,s,i},
  For[
    i=1; term=ModByDegree[$MaxDegree+1,Expand[expr]] ; s=0,
    term=!=0,
    ++i; term=ModByDegree[$MaxDegree+1,Expand[term*(-expr)]],
    s+=Expand[term/i]
  ];
  s
]

A3Invert[1]=1
A3Invert[1+expr_] := Module[{term,s},
  For[
    term=Expand[-expr] ; s=1,
    !term===0,
    term=term~A3Multiply~(-expr),
    s+=term
  ];
  s
]

A4Invert[1]=1
A4Invert[1+expr_] := Module[{term,s},
  For[
    term=Expand[-expr] ; s=1,
    !term===0,
    term=term~A4Multiply~(-expr),
    s+=term
  ];
  s
]

BPentagon[phi_] := Module[
  {inv},
  phi ~A4Multiply~ BAct[1,{2,3},4][phi] ~A4Multiply~ BAct[2,3,4][phi]
    ~A4Multiply~ BAct[1,2,{3,4}][inv=A3Invert[phi]] ~A4Multiply~
    BAct[{1,2},3,4][inv]
]

BHexagonP[r_,phi_] := (
  BAct[{1,2},3][A3Invert[r]] ~A3Multiply~ phi ~A3Multiply~ BAct[2,3][r]
  ~A3Multiply~ BAct[1,3,2][A3Invert[phi]] ~A3Multiply~ BAct[1,3][r]
  ~A3Multiply~ BAct[3,1,2][phi]
)

BHexagonM[r_,phi_] := Module[
  {inv},
  BAct[{1,2},3][r] ~A3Multiply~ phi ~A3Multiply~ BAct[2,3][inv=A3Invert[r]]
  ~A3Multiply~ BAct[1,3,2][A3Invert[phi]] ~A3Multiply~ BAct[1,3][inv]
  ~A3Multiply~ BAct[3,1,2][phi]
]

ModifiedB[n_] := BernoulliB[n]/(2n*n!)

WheelingElement[k_] := WheelingElement[k] = MergeGraphs[Expand[Product[
  Sum[
    (ModifiedB[2n]Ring[2n])^p/p!,
    {p,0,Floor[k/(2n)]}
  ],
  {n,1,Floor[k/2]}
]]]-1+Graph[]
Wheel[g_Graph,c_] := LeftStableFusion[WheelingElement[Count[g,_c]], g, {UV,c}]
Wheel[expr_,c_] := Expand[MergeGraphs[expr] /. g_Graph :> Wheel[g,c]]

DeWheelingElement[k_] := DeWheelingElement[k] = MergeGraphs[Expand[Product[
  Sum[
    (-ModifiedB[2n]Ring[2n])^p/p!,
    {p,0,Floor[k/(2n)]}
  ],
  {n,1,Floor[k/2]}
]]]-1+Graph[]
DeWheel[g_Graph,c_] := LeftStableFusion[DeWheelingElement[Count[g,_c]],g,{UV,c}]
DeWheel[expr_,c_] := Expand[MergeGraphs[expr] /. g_Graph :> DeWheel[g,c]]

CD2CDP[cd_CD] := Module[
  {m=Length[cd]/2,i,beg,end},
  beg=Table[Position[cd,i],{i,1,m}];
  end=Reverse /@ beg;
  CDP @@ (Range[2m] /. Thread[Rule[Flatten[beg],Flatten[end]]])
]

CD2Graph[expr_] := expr /. cd_CD :> 
  Graph @@ MapIndexed[(ColoredUV@@#2)[#1]&, CD2CDP[cd]]

CD2CC[expr_] := Expand[expr /. cd_CD :> Fold[
  CCReduce[AMergeColors[#1,C1,#2,C1]]&,
  CD2Graph[cd],
  ColoredUV /@ Range[2,Length[cd]]
]]

BAct[spec___][expr_] := Module[{rule,i},
  rule=Replace[#,c_Integer:>{c}]& /@ {spec};
  rule=Table[
    ColoredUV[i][j_] -> Plus @@ (ColoredUV[#][j]& /@ rule[[i]]),
    {i,Length[rule]}
  ];
  CCReduce[
    expr /. rule /. p[lft___,0,rgt___] -> 0 /. g_Graph :> Distribute[g]
  ]
]

(* Need to Do:

Apply to Phi.

*)

