BeginPackage["IHX`",{"Local`","Profile`","Subspace`","Graph`"}]

DimB[8, {2, 1, 1, 1, 1, 1, 1}] = Pass
DimB[8, {2, 1, 1, 1, 1, 1}] = Pass
DimB[8, {1, 1, 1, 1, 1, 1, 1}] = Pass

SetupCompress[g_Graph] := ($CompressionKey = {
  g,
  Length[g]^(Range[Plus @@ Length /@ g]-1)
})

CompressGraph[g_Graph] := (Cases[g,_Integer,{2}]-1) . $CompressionKey[[2]]

UncompressGraph[c_Integer] := Block[
  {len=Length[$CompressionKey[[1]]], t=c, r},
  $CompressionKey[[1]] /. i_Integer :> (
    t=(t-(r=Mod[t,len]))/len;
    1+r
  )
]

IHX[ig_,opts___] := Profile[IHX,Module[
  {gs,rs,c,z,t,fig},
  fig=If[Head[ig]===List, First[ig], ig];
  Print["Type: ",(Plus @@ Head /@ fig) /. (_?EdgeQ) -> 0,"..."];
  t=TimeUsed[];
  {gs,rs}= Make[ig, opts, SolveIHX -> True];
  Print[
    (Plus @@ Head /@ fig) /. (_?EdgeQ) -> 0,"; ",
    "#graphs=",c=Length[gs],", ",
    "#rels=",AddVector[rs,0],", ",
    "time=",TimeUsed[]-t
  ];
  {gs,rs}
]]

Options[Make] := {
  TryUOVs -> False,
  SolveIHX -> True,
  AllGraphsGiven -> False,
  Relations -> {}
}

Make[0,opts___] := {{},EmptySubspace[1]}
Make[ig_Graph,opts___] := Make[{ig},opts]

Make[igs_List,opts___] := Profile[Make,Module[
  {rels,checked={},InCanonicalForm={},k,SerialNumber,mold,
    n,g,ng,bg,gI,gH,gX,l=0,AppendToAllGraphs,c,at,r,ndashed,pos,
    AllGraphs,
    solve = (SolveIHX /. {opts} /. Options[Make]),
    allgiven = (AllGraphsGiven /. {opts} /. Options[Make]),
    prerels = (Relations /. {opts} /. Options[Make]),
    TryUOVs = (TryUOVs /. {opts} /. Options[Make])
  },

  AppendToAllGraphs[g_] := Profile[AppendToAllGraphs,
    ng=CanonicalForm[g];
    If[!TryUOVs,ng=ReleaseHold[ng]];
    If[0===ng, 0,
      c=NumericalCoefficient[ng];
      bg=BareGraph[ng];
      If[l==0, SetupCompress[mold=bg]];
      at=SerialNumber[CompressGraph[bg]];
      If[Head[at]===SerialNumber,
        SerialNumber[CompressGraph[bg]]=++l;
        AllGraphs[l]=CompressGraph[bg];
        If[solve && MatchQ[ng,Hold[0*_]],
          AddVector[rels, BasisVector[l]]
        ];
        AppendTo[checked,Table[False,{ndashed}]];
        If[n>1 && l-n>1,
          checked[[l-n+1,First[Cases[ng,Marker[j_]:>inpos[j],Infinity]]]]=True
        ];
        c*BasisVector[l],
        If[at-n>1,
          checked[[
            at-n+1,
            First[Cases[ng,Marker[j_] :> inpos[j],Infinity]]
          ]]=True
        ];
        c*BasisVector[at]
      ] 
    ]
  ];

  n=1;
  ndashed=Count[First[igs],_?EdgeQ];
  Print[StringForm["Appending `` graphs...",Length[igs]]];
  AppendToAllGraphs /@ igs;
  If[solve,
    rels=EmptySubspace[If[allgiven, l, 200000]];
    Print[StringForm["Considering `` given relations...",Length[prerels]]];
    Scan[AddVector[rels,#]&, prerels];
    nrels=AddVector[rels,0],
    rels={}
  ];
  If[l>0,pos=Flatten[Position[mold,_?EdgeQ,{1}]]];

  While[n<=l,
    If[Mod[n,10]===0,
      Print[n,"(",l-n,"): ",
        If[solve,StringForm["``/",l-nrels],""],
        l,"."
      ]
    ];
    inpos[j_] := First[First[Position[pos,j]]];
    ndashed=Length[pos];

    Do[
      If[!checked[[1,i]],
        gI=UncompressGraph[AllGraphs[n]]~Append~Marker[pos[[i]]];
        gH=I2H[gI,pos[[i]]];
        gH=AppendToAllGraphs[gH];
        gX=I2X[gI,pos[[i]]];
        gX=AppendToAllGraphs[gX];
        If[solve,
          r=BasisVector[n]-gH+gX;
          Profile[AddVector,nrels=AddVector[rels,r];]
        ]
      ],
      {i,1,ndashed}
    ];
    checked=Rest[checked];
    ++n
  ];
  {AllGraphs /@ Range[l],rels}
]]

TagDashed[g_Graph] := Module[
  {i=0},
  g /. Dashed :> (
    d=ToExpression["Dashed"~StringJoin~ToString[++i]];
    SetAttributes[d,Orderless];
    EdgeQ[d]=True;
    EdgeQ[_d]=True;
    d
  )
]

ColoredUV[c_] := ToExpression["C"~StringJoin~ToString[c]]

Unknown /: c_*Unknown = Unknown
Unknown /: c_+Unknown = Unknown
TooHard /: c_*TooHard = TooHard
TooHard /: c_+TooHard = TooHard

BSeedGraph::dontknow = "I Don't know how to make a `1` `2` graph"

BSeedGraph::notpartition = "`1` is not a partition"

BSeedGraph[_,cols_List] := (
  Message[BSeedGraph::notpartition,cols];
  Unknown
) /; !OrderedQ[Reverse[cols]]

BSeedGraph[1,{}] := ThetaGraph[]

BSeedGraph[m_,{}] := Prism[m] /; m>=2

BSeedGraph[_,{1}] = 0

BSeedGraph[2,{3}] = 0

BSeedGraph[4,{3}] = 0

BSeedGraph[m_,{c_,1}] := Nest[Jump, Ring[{C2}~Join~Table[C1,{c}]], m-c-1] /;
  (m>=c+2 || (m==c+1 && EvenQ[m]))

BSeedGraph[m_,{c_,1}] := Join[
  ReplacePart[Ring[m-1] /. UV -> C1, Dashed[3,3m-2], 1],
  Graph[OV[1,3m-1,3m], C1[3m-2], C2[3m-2]]
] /; (m>=3 && m==c+1 && OddQ[m])

BSeedGraph[2,{1,1,1}]=Graph[OV[2,3,4],C1[1],C2[1],C3[1]]

BSeedGraph[m_,{1,1,1}] := Nest[
  Jump,
  Graph[
    Dashed[4,5],C2[4],C3[4],OV[1,2,3],OV[1,6,7],
    Dashed[5,8],Dashed[5,8],OV[6,7,9],C1[8]
  ],
  m-3
] /; m>=3

BSeedGraph[m_,{c_}] := 0 /; m<c

BSeedGraph[m_,{m_}] := 0 /; OddQ[m]

BSeedGraph[m_,{c_}] := Nest[Jump,Ring[Table[C1,{c}]],m-c] /;
  (m>=c && EvenQ[c])

BSeedGraph[m_,{c_}] := Nest[
  AddUnivalent[#,C1,2]&,
  Nest[Jump,ThetaGraph[C1],m-c-1],
  c-1
] /; (m>c+1 && OddQ[c])

BSeedGraph[m_,{c_}] := (BSeedGraph[m,{c}] = First[
  Make[
    Nest[Jump,Ring[Table[ColoredUV[1],{c}]],m-c] /. OV -> UOV,
    SolveIHX -> False
  ] /. UOV -> OV
]) /; (m==c+1 && OddQ[c])

BSeedGraph[m_,{m_,1}] = 0

BSeedGraph[m_,cols_List] := Module[
  {i,cl},
  cl=Flatten[Table[Table[ColoredUV[i],{cols[[i]]}],{i,1,Length[cols]}]];
  cl=RotateLeft[cl];
  {cl[[2]],cl[[-3]]}={cl[[-3]],cl[[2]]};
  If[cl[[1]]===cl[[2]] || cl[[-1]]===cl[[-2]],
    Message[BSeedGraph::dontknow, m, cols]; Unknown,
    Nest[Jump,Tree[cl],m+1-(Plus @@ cols)]
  ]
] /; (
  Length[cols]>=2 &&
  m>=(Plus @@ cols)-1 &&
  !MatchQ[cols,{_,1}] &&
  cols =!= {1,1,1} &&
  !(Length[cols]==2 && m==First[cols]+1 && EvenQ[m]) &&
  {m,cols} =!= {8,{4,3,2}} &&
  {m,cols} =!= {6,{4,3}}
)

BSeedGraph[m_,{c_,2}] := 0 /; (m==c+1 && EvenQ[m])

BSeedGraph[8,{4,3,2}] = Tree@{C1,C2,C1,C1,C1,C2,C3,C2,C3}
BSeedGraph[6,{4,3}] = Tree@{C1,C2,C1,C1,C2,C1,C2}

BSeedGraph[m_,cols_List] := (
  Message[BSeedGraph::dontknow, m, cols];
  Unknown
)

BBasis[m_,cols_List] := BBasis[m,cols] = Module[
  {n,diags,rels,basis,c,i,fname,s},
  If[(g=BSeedGraph[m,cols])===Unknown,Unknown,
    {diags,rels}=IHX[g];
    n=Length[diags];
    basis=Complement[
      Range[n],
      SubspaceBasis[rels] //. {v_Plus :> First[v], c_.*BasisVector[i_] :> i}
    ];
    basis=UncompressGraph /@ (diags[[basis]]);
    fname="BBasis/"~StringJoin~ToString[m]~StringJoin~(
      StringJoin @@ (ToString[StringForm["-``",#]]& /@ cols)
    );
    s=OpenWrite[fname];
    WriteString[s,"BBasis[",m,",",cols,"] =\n"];
    Write[s,basis];
    Close[s];
    Run["gzip -f", fname];
    basis
  ]
]

DimB[m_,cols_List] := DimB[m,Reverse[Sort[cols]]] /; !OrderedQ[Reverse[cols]]

DimB[m_,cols_List] := DimB[m,DeleteCases[cols,0]] /; MemberQ[cols,0]

DimB[1,{2}] = 1

DimB[m_,{2}] := DimB[m-1,{}] /; m>=2

DimB[m_,{1,1}] := DimB[m,{2}] /; m<=10

DimB[m_,cols_List] := (m-1)! /; (Union[cols]==={1} && Length[cols]===m+1)

DimB[m_,cols_List] := If[BBasis[m,cols]===Unknown,
  Unknown,
  Length[BBasis[m,cols]]
]

If[Head[DimBTableTimeLimit]===Symbol,DimBTableTimeLimit=5]

PrintDimBTable[m_] := Module[
  {dimb,e,part},
  dimb=DimB[m,{}]=TimeConstrained[DimB[m,{}],DimBTableTimeLimit,TooHard];
  Print["DimB[",m,",",{},"] = ",dimb,";"];
  Do[
    dimb=DimB[m,{e}]=TimeConstrained[DimB[m,{e}],DimBTableTimeLimit,TooHard];
    Print["DimB[",m,",",{e},"] = ",dimb,";"];
    For[part={e-1,1}, e>1 && part =!= {e}, part=NextPartition[part],
      dimb=DimB[m,part]=TimeConstrained[
        DimB[m,part],
        DimBTableTimeLimit,
        TooHard
      ];
      Print["DimB[",m,",",part,"] = ",dimb,";"];
    ];
    Print[" "],
    {e,1,m+1}
  ]
]

DimBTable[m_Integer] := Module[
  {dimb,e,part,table},
  table=DimBTable[beta[m,{}]~is~DimB[m,{}]];
  Do[
    dimb=DimB[m,{e}]=TimeConstrained[DimB[m,{e}],DimBTableTimeLimit,TooHard];
    AppendTo[table,beta[m,{e}]~is~dimb];
    For[part={e-1,1}, e>1 && part =!= {e}, part=NextPartition[part],
      dimb=DimB[m,part]=TimeConstrained[
        DimB[m,part],
        DimBTableTimeLimit,
        TooHard
      ];
      AppendTo[table,beta[m,part]~is~dimb]
    ],
    {e,1,m+1}
  ];
  table
]

DimBTable[m_Integer * AndUp] := DimBTable @@ (
  (##[[1,1]]~is~#[[2]])& /@ (
    Select[DownValues[DimB],TrueQ[#[[1,1,1]] >= m]&] /. DimB -> beta
  )
)

DimBTable[m_Integer * Known] := DimBTable @@ (
  (##[[1,1]]~is~#[[2]])& /@ (
    Select[DownValues[DimB],TrueQ[#[[1,1,1]] == m]&] /. DimB -> beta
  )
)

DimP[m_,Ihara] := {0,0,1,0,1,0,2,1,4,2,9,7,21}[[m]] /; m<=13

DimP[m_,n_] := DimP[m,n] = Module[
  {e,sum=0,part,k,h,l},
  Do[
    sum += DimB[m,{e}] * n;
    For[part={e-1,1}, part =!= {e}, part=NextPartition[part],
      sum += If[TrueQ[(l=Length[part])>n], 0, Expand[
        DimB[m,part] * Product[n+1-k, {k,l}] / (
          Times @@ (((List @@ Plus @@ h /@ part) /. _h -> 1)!)
        )
      ]]
    ],
    {e,2,m+1}
  ];
  sum
]

DimPT[m_,n_] := DimPT[m,n] = Module[
  {e,sum=0,part,k,h,l},
  Do[
    sum += DimB[m,{e}] * n;
    For[part={e-1,1}, part =!= {e}, part=NextPartition[part],
      sum += If[TrueQ[(l=Length[part])>n], 0, Expand[
        DimB[m,part] * Product[n+1-k, {k,l}] / (
          Times @@ (((List @@ Plus @@ h /@ part) /. _h -> 1)!)
        )
      ]]
    ],
    {e,m+1,m+1}
  ];
  sum
]

DimA[1,n_] := DimP[1,n]
DimA[m_,n_] := Module[
  {sum=0,part,h,k,l,dimp},
  dimp=Table[DimP[k,n],{k,m-1}];
  sum += DimP[m,n];
  For[part={m-1,1}, part =!= {m}, part=NextPartition[part],
    sum += (Times @@ h /@ part) /. {
      h[k_]^l_ :> Binomial[dimp[[k]]+l-1,l],
      h[k_] :> dimp[[k]]
    }
  ];
  sum
] /; m>=2

DimAr[1,n_] := DimP[1,n] - n
DimAr[m_,n_] := Module[
  {sum=0,part,h,k,l,dimp},
  dimp={DimP[1,n] - n}~Join~Table[DimP[k,n],{k,2,m-1}];
  sum += DimP[m,n];
  For[part={m-1,1}, part =!= {m}, part=NextPartition[part],
    sum += (Times @@ h /@ part) /. {
      h[k_]^l_ :> Binomial[dimp[[k]]+l-1,l],
      h[k_] :> dimp[[k]]
    }
  ];
  sum
] /; m>=2

(* DimAR --- all primitive knot-only invariants removed *)
DimAR[1,n_] := DimP[1,n] - n
DimAR[m_,n_] := Module[
  {sum=0,part,h,k,l,dimp},
  dimp=Table[DimP[k,n] - n*DimP[k,1], {k,m-1}];
  sum += DimP[m,n] - n*DimP[m,1];
  For[part={m-1,1}, part =!= {m}, part=NextPartition[part],
    sum += (Times @@ h /@ part) /. {
      h[k_]^l_ :> Binomial[dimp[[k]]+l-1,l],
      h[k_] :> dimp[[k]]
    }
  ];
  sum
] /; m>=2

(* DimATR - same as dimAR, but working only with Tree-level invariants *)
DimATR[1,n_] := DimPT[1,n] - n
DimATR[m_,n_] := Module[
  {sum=0,part,h,k,l,dimp},
  dimp=Table[DimPT[k,n] - n*DimPT[k,1], {k,m-1}];
  sum += DimPT[m,n] - n*DimPT[m,1];
  For[part={m-1,1}, part =!= {m}, part=NextPartition[part],
    sum += (Times @@ h /@ part) /. {
      h[k_]^l_ :> Binomial[dimp[[k]]+l-1,l],
      h[k_] :> dimp[[k]]
    }
  ];
  sum
] /; m>=2

EndPackage[]
