(* ::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.                                                         *)
(************************************************************************)



SetAttributes[{P,S}, Orderless]; dot /: dot[_]^k_ /; k>=2 := 0;
(\[Sigma]_S)[i_] := \[Sigma][i] = First@Cases[\[Sigma], P[i,j_] :> j];


EC[\[Lambda]_List] := Module[{\[Rho], ec=\[Lambda]}, (* "Finding Equivalence Classes" *)
Do[ \[Rho] = First /@ Position[ec, i];
ec=Append[Delete[ec, List /@ \[Rho]], Union@@(ec[[\[Rho]]])],
{i, Union @@ \[Lambda]} ]; ec];
EC[\[Lambda]__S] := EC[Join[\[Lambda]] /. S|P->List];
ECP[\[Lambda]__] := Union@@Replace[EC[\[Lambda]], c_ :> ((# -> First[c])& /@ c), {1}] ;


VCLaw[\[Beta]_S, \[Mu]_S, \[Tau]_S] := VCLaw[\[Beta],\[Mu],\[Tau]]=Module[{ins, outs, p, \[Chi]s, h, dec} ,
ins = First /@ Join[EC[\[Beta],\[Mu]], EC[\[Mu],\[Tau]]];
outs=First/@EC[\[Beta],\[Tau]]; p=ECP[\[Beta],\[Mu],\[Tau]];
\[Chi]s = Times @@ (h/@ Join[ins, outs] /. p);
\[Chi]s *= PowerExpand[(Times @@ (h /@ (Last /@ p)))^(-1/2)];
dec = \[Chi]s /. h[i_]^\[Chi]_. :> (2dot[i])^((2-\[Chi])/2);
dec *= Product[If[i==(i/.p), 1, dot[i]+dot[i/.p]], {i,outs}];
{Expand[dec], Table[dot[i]->dot[i/.p], {i, Union[ins]}]}];
VC[\[Beta]_S, \[Mu]_S, \[Tau]_S, dots_] := Module[{dec, dotslaw},
{dec, dotslaw} = VCLaw[\[Beta],\[Mu],\[Tau]]; Expand[dec*(dots /.dotslaw)] ];


m0[i_,j_][\[Sigma]_S] := m0[i,j][\[Sigma]] = If[\[Sigma][i]==j, DeleteCases[\[Sigma], P[i,j]],
Append[DeleteCases[\[Sigma], P[i,_]|P[_,j]], P[\[Sigma][i],\[Sigma][j]]] ];
m[i_,j_][\[Sigma]_S] := m0[i,j][\[Sigma]]*If[\[Sigma][i]==j, {q,q^-1}, {1}];
m[i_,j_][q^k_. \[Sigma]_S] := q^k m[i,j][\[Sigma]];


m[i_,j_][Cob[\[Beta]_S, \[Tau]_S, dots_]] := Module[{p, ijdot, np, ndots, x},
p=ECP[\[Beta],\[Tau]]; ijdot=dot@Min[i,j]; np=ECP[m0[i,j][\[Beta]],m0[i,j][\[Tau]]];
ndots = Which[\[Beta][i]==j && \[Tau][i]==j, ({
 {ijdot, 0},
 {1, ijdot}
}),
\[Beta][i]==j && \[Tau][i]!=j, {{1, ijdot}},
\[Beta][i]!=j && \[Tau][i]==j, {{ijdot},{1}},\[Beta][i]!=j && \[Tau][i]!=j, {{If[(i/.p)!=(j/.p), 1, dot[\[Beta][i]]+dot[\[Tau][i]]]}} ];
ndots=Expand[dots*ndots] /. dot[k_] :> dot[k /. {i->\[Beta][i], j->\[Beta][j]} /. {i->\[Tau][i], j->\[Tau][j]} /. np ];
If[\[Beta][i]==j && \[Tau][i]==j, Coefficient[ndots /. ijdot -> x, x], ndots]];


m[i_,j_][Kom[\[CapitalOmega]_, d_]] := Kom[
Flatten /@ Map[m[i,j], \[CapitalOmega], {2}],
Table[If[Length@\[CapitalOmega][[k]]==0 || Length@\[CapitalOmega][[k+1]]==0, 0,
Table[m[i,j][Cob[\[CapitalOmega][[k,b]], \[CapitalOmega][[k+1, a]], d[[k,a,b]]] /. q->1],
{a, Length@\[CapitalOmega][[k+1]]}, {b, Length@\[CapitalOmega][[k]]}
] // ArrayFlatten  ],
{k,Length@d}]  ];


(Kom[\[CapitalOmega]_, d_]//Cob[q^p1_. \[Beta]_, q^p2_. \[Tau]_, 1]) := Module[{L, \[Rho], \[Delta], k},
L=Length[\[CapitalOmega]]; \[Rho][k_] := \[Rho][k]=Length[\[CapitalOmega][[k]]]; \[Rho][0]=\[Rho][L+1]=0;
Kom[
MapThread[Join, List @@@ {
Append[\[CapitalOmega] /. \[Sigma]_S :> q^p1 Join[\[Beta],\[Sigma]], {}],
Prepend[\[CapitalOmega] /. \[Sigma]_S :> q^p2 Join[\[Tau],\[Sigma]], {}] }],
Table[
If[\[Rho][k]+\[Rho][k-1]==0 || \[Rho][k+1]+\[Rho][k]==0, 0,
\[Delta]=Table[0, {\[Rho][k+1]+\[Rho][k]}, {\[Rho][k]+\[Rho][k-1]}];
If[\[Rho][k]\[Rho][k+1]!=0, \[Delta][[1;;\[Rho][k+1], 1;;\[Rho][k]]]=d[[k]]];
If[\[Rho][k]!=0, \[Delta][[\[Rho][k+1]+1;;\[Rho][k+1]+\[Rho][k], 1;;\[Rho][k]]]=(-1)^k IdentityMatrix[\[Rho][k]]];
If[\[Rho][k-1]\[Rho][k]!=0, \[Delta][[\[Rho][k+1]+1;;\[Rho][k+1]+\[Rho][k], \[Rho][k]+1;;\[Rho][k]+\[Rho][k-1]]]=d[[k-1]]];
\[Delta]
], {k,L}   ] ] ]


Contract[kom_Kom] := Module[{\[CapitalOmega], d, L, \[Rho], k, done, a, b, \[Phi], \[Gamma]\[Delta]},
{\[CapitalOmega], d} = List @@ kom; L=Length@d;\[Rho][k_] := Length@\[CapitalOmega][[k]];
For[k=1, k<=L, ++k,
done=False; While[!done, done=True;
For[a=1, a<=\[Rho][k+1], ++a, For[b=1, b<=\[Rho][k], ++b,
If[NumberQ[\[Phi]=d[[k,a,b]]] && \[Phi]!=0 && \[CapitalOmega][[k+1,a]]==\[CapitalOmega][[k,b]],
done=False;
If[\[Rho][k]<=1 || \[Rho][k+1]<=1, d[[k]]=0,
\[Gamma]\[Delta] = Table[
VC[\[CapitalOmega][[k,n]]/.q->1, \[CapitalOmega][[k+1,a]]/.q -> 1,\[CapitalOmega][[k+1,m]]/.q->1, d[[k,a,n]]d[[k,m,b]]],  {m, \[Rho][k+1]}, {n, \[Rho][k]} ];
d[[k]] = Expand@Drop[d[[k]]-\[Phi]^-1 \[Gamma]\[Delta], {a}, {b}]];
\[CapitalOmega][[k]]=Drop[\[CapitalOmega][[k]], {b}]; \[CapitalOmega][[k+1]]=Drop[\[CapitalOmega][[k+1]], {a}];
If[k>1 && \[Rho][k-1]>0, d[[k-1]]= Drop[d[[k-1]], {b}]];
If[k<L && \[Rho][k+2]>0, d[[k+1]]= Drop[d[[k+1]], {}, {a}]];
If[a<=\[Rho][k+1], --a]; b=\[Rho][k] ]  ]] ] ];
Kom[\[CapitalOmega], d] ];


Kom[] = Kom[{{S[]}},{}];
Cob@Xp[i_,j_,k_,l_]:=Cob[q S[P[-i,j],P[k,-l]],q^2 S[P[-i,-l],P[j,k]],1];
Cob@Xm[i_,j_,k_,l_]:=Cob[q^-2 S[P[-i,-j],P[k,l]],q^-1 S[P[-i,l],P[-j,k]],1];
Cob[x_X] := Cob[If[PositiveQ[x],Xp@@x, Xm@@x]];


KhComplex[L_] := Module[
{pd=PD[L], kom=Kom[], inside={}, pos},
While[Length[pd]>0,
pos = Last[Ordering[(Length[(List @@ #)\[Intersection]inside])& /@ pd]];
kom=kom//Cob[pd[[pos]]];
(kom = Contract[kom // m[#, -#]])& /@ ((List@@pd[[pos]])\[Intersection]inside);
inside = inside\[Union](List @@ pd[[pos]]); pd=Drop[pd, {pos}] ];
kom ];
KhPoly[L_] := Expand[t^(-Length@Select[PD@L, NegativeQ]+Range[0,Crossings[L]]).(List @@ Plus @@@ First @ KhComplex[L]) /. S[] -> 1]



