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

(* Temporarily exported: *)
$GraphLink; MLCanonicalForm; CC

(* Profile labels *)
FindNewPartition; CFArtificialSubdivision; CFIteration; RightPermute;

Graph::usage = "\n
	Graph[v1,v2,...] represents a graph whose vertices are v1,v2,..."

dK::usage = "\n
	dK is the Kontsevich differential on the space of Graphs."

ThetaGraph::usage = "\n
	ThetaGraph[l1,l2,l3] returns a theta-graph with univalent vertices \n
	specified by the lists l1-3 on its three arcs. If any of these\n
	lists is omitted, it is taken as empty and if any of the arguments\n
	is not a list, it is taken as a list of length 1."

Marker::usage = "\n
	Marker[i,j,...] is a fake vertex in a graph, always put last and\n
	ignored in CanonicalForm, used only to mark vertices i,j,..."

MarkerQ::usage = "\n
	MarkerQ[v] returns True if v is a marker, and False otherwise."

EdgeQ::usage = "\n
	EdgeQ[v] returns True if v is a unoriented edge, and False otherwise."

OddVertexQ::usage = "\n
	OddVertexQ[v] returns True if v is an odd vertex, (that is, if\n
	v's name begins with \"Odd\"), and False otherwise."

Laplacian::usage = "\n
	Laplacian[g] returns the Laplacian matrix of the graph g."

AdjacencyMatrix::usage = "\n
	AdjacencyMatrix[g] returns the adjacency matrix matrix of the graph g."

AdjacencyMatrix2Graph::usage = "\n
	AdjacencyMatrix2Graph[mat] constructs the graph whose adjacency \n
	matrix is mat, using only UOV vertices."

Prism::usage = "\n
	Prism[m] generates a prism of degree m."

Ring::usage = "\n
	Ring[m] generates a ring of degree m. Ring[list] does same but with\n
	the univalent vertices taken from the list list."

Tree::usage = "\n
	Tree[list] generates a trivalent tree whose univalent leafs are\n
	colored as in the list list."

Graph2CC::usage = "\n
	Graph2CC[g] converts the Graph g to the CC format. (No UVs\n
	implemented yet)! Graph2CC[expr] does same to all Graphs appearing\n
	in expr."

CC2Graph::usage = "\n
	CC2Graph[cc] converts the Chinese Character cc to the Graph format.\n
	(No UVs implemented yet)! CC2Graph[expr] does same to all Chinese\n
	Characters appearing in expr."

Phi::usage = "\n
	Phi[expr] computes Phi of all Chinese Characters appearing in expr,\n
	using the external C++ program PhiOf."

UV::usage = "\n
	UV[i] represents an unmarked univalent vertex connected to the #i\n
	vertex in a graph."

OV::usage = "\n
	OV[i,j,k] represents an oriented trivalent vertex connected to\n
	vertices #i, #j, and #k."

UOV::usage = "\n
	UOV[i,j,k] represents an unoriented trivalent vertex connected to\n
	vertices #i, #j, and #k."

OddUOV::usage = "\n
	OddUOV[i,j,k] represents an odd unoriented trivalent vertex\n
	connected to vertices #i, #j, and #k."

Dashed::usage = "\n
	Dashed[i,j] represents a dashed line connecting two vertices."

Arrow::usage = "\n
	Arrow[i,j] represents an oriented edge connecting two vertices."

BareVertices::usage = "\n
	BareVertices[g] return a list of the bare (unlabeled) vertices of g."

BareGraph::usage = "\n
	BareGraph[g] strips numerical coefficients from g and removes all\n
	markers."

Permute::usage = "\n
	Permute[p,g] permutes the vertices of g according to\n
	the permutation p. p and g must have the same length.\n
	Permute[g,p] uses the right action of the permutation\n
	group on graphs."

Invert::usage = "\n
	Invert[p] inverts the permutation represented by the list p."

CanonicalForm::usage = "\n
	CanonicalForm[g] brings a graph g to a canonical form. The\n
	arcs emenating from a vertex are always assumed ordered, and\n
	their order is changed by CanonicalForm only through Permute."

UseNauty::usage = "\n
	UseNauty is an option for CanonicalForm. If True, CanonicalForm\n
	calls nauty rather than finding the canonical form on its own.\n
	Results may differ!!!"

InCanonicalForm::usage ="\n
	InCanonicalForm is a list containing graphs known to be in \n
	canonical form. Graphs are added to it automatically whenever\n
	CanonicalForm takes more than $CanonicalFormLearnTime seconds\n
	to finish."

SortingPermutation::usage = "\n
	SortingPermutation[expr] returns the permutation whose (right)\n
	action on expr will sort it."

NumericalCoefficient::usage = "\n
	NumericalCoefficient[expr] returns the Numerical Coefficient of term."

I2H::usage = "\n
	I2H[g,i] takes the I around the i'th dashed line in g and turns it\n
	into an H."

I2X::usage = "\n
	I2X[g,i] takes the I around the i'th dashed line in g and turns it\n
	into an X."

Jump::usage = "\n
	Jump[g,i] replaces the i'th vertex of g (assuming it satisfies \n
	EdgeQ) by a blob. Jump[g] does same to the first dashed line in g."

AddUnivalent::usage = "\n
	AddUnivalent[g,uv,i] adds a univalent vertex of type uv along the\n
	i'th vertex of g (assuming it satisfies EdgeQ). AddUnivalent[g,uv]\n
	does same to the first dashed line in g."

$CanonicalFormLearnTime=1

Begin["`Private`"]

SetAttributes[Dashed,Orderless]

SetAttributes[UOV,Orderless]
SetAttributes[OddUOV,Orderless]

InCanonicalForm={}

MarkerQ[v_] := (Head[v]===Marker)

EdgeQ[_Dashed] = True
EdgeQ[Dashed] = True
EdgeQ[_] = False

OrientedVertexQ[_OV] = True
OrientedVertexQ[OV] = True
OrientedVertexQ[_Arrow] = True
OrientedVertexQ[Arrow] = True
OrientedVertexQ[_] = False

OddVertexQ[v_Symbol[___]] := StringMatchQ[ToString[v],"Odd*"]

NumericalCoefficient[Hold[expr_]] := NumericalCoefficient[expr]
NumericalCoefficient[c_?NumberQ]:=c
NumericalCoefficient[c_?NumberQ * _] := c
NumericalCoefficient[_]=1

Invert[p_List] := Last /@ Sort[Thread[{p,Range[Length[p]]}]]

SortingPermutation[expr_] := Profile[SortingPermutation,Last /@ Sort[
  Thread[{expr,Range[Length[expr]]}]
]]

Permute[p_List,g_Graph] := Map[Replace[#,i_ :> p[[i]]]&,g,{2}][[Invert[p]]]

Permute[c_?NumberQ*g_Graph,p_List] := c Permute[g,p]

Permute[g_Graph,p_List] := Profile[RightPermute,
  rule=MapThread[Rule,{p,Range[Length[p]]}];
  RPsign=Signature[Position[g,_?OddVertexQ,{1}] /. rule];
  RPng=Map[Replace[#,rule]&,g,{2}][[p]];
  RPng=RPng /. ov_?OrientedVertexQ :> (RPsign*=Signature[ov]; Sort[ov]);
  RPsign*RPng
]

BareVertices[g_Graph] := Apply[List, List @@ g, {1}]

BareGraph[_*h_Hold]=BareGraph[h]
BareGraph[Hold[(c_:1)*(g_Graph)]] := BareGraph[g]
BareGraph[g_] := DeleteCases[
  If[Head[g]===Times, Times @@ Cases[g,_Graph], g],
  _?MarkerQ
]

Options[CanonicalForm] := {UseNauty -> True}
$NautyCalls=0
CanonicalForm[0,opts___]=0
CanonicalForm[c_?NumberQ * g_Graph,opts___] := c*CanonicalForm[g,opts]

CanonicalForm[g_Graph,opts___] := Profile[CanonicalForm, If[
  UseNauty /. {opts} /. Options[CanonicalForm],
  Block[{ng,done=False},
    While[!done,
      If[Head[$GraphLink]=!=LinkObject,
        $GraphLink=LinkOpen["Graph",LinkMode -> Launch];
        Print[LinkRead[$GraphLink]]
      ];
      ++$NautyCalls; LinkWrite[$GraphLink,MLCanonicalForm[g]];
      ng=LinkReadHeld[$GraphLink];
      If[Mod[$NautyCalls,15000]=!=0 && First[LinkError[$GraphLink]]===0,
        done=True,
        LinkClose[$GraphLink];
        Unset[$GraphLink]
      ]
    ];
    If[ReleaseHold[ng]===0, ng, ReleaseHold[ng]]
  ],(
    Block[
      {
        nv=Length[g],
        van=False,
        VertexList,VertexTypes,part,ng,c,gg,am,mat,i,blobQ
      },
Print["Warning: unmaintained program section!!!"];
      If[MemberQ[InCanonicalForm, BareGraph[g]], g,
        VertexList={
          List @@ MarkerQ /@ g,
          List @@ Length /@ g,
          List @@ Head /@ g
        };
        blobQ=Table[
          If[EdgeQ[g[[i]]] || Head[g[[i]]] === Arrow,
            Position[
              Position[Sort /@ g,Sort[g[[i]]],{1},Heads->False],
              {i}
            ][[1,1]],
            1
          ],
          {i,1,nv}
        ];
Profile[AdjacencyMatrix,
        am=AdjacencyMatrix[g];
        mat=am=am.am;
        VertexList=Join[
          VertexList,
          Table[Transpose[mat=mat.am,{1,1}],{Floor[nv/4]}]
        ];
];
        AppendTo[VertexList,blobQ];
        VertexList=Transpose[VertexList];
        ng=Permute[g,SortingPermutation[VertexList]];
        VertexList=Sort[VertexList];
        VertexTypes=Union[VertexList];
        part=VertexList /. Thread[Rule[VertexTypes,Range[Length[VertexTypes]]]];
        $CanonicalFormStartTime=TimeUsed[];
        If[Last[part]<nv,ng=CFIteration[ng,part]];
        If[$CanonicalFormStartTime+$CanonicalFormLearnTime<TimeUsed[] && ng=!=0,
          gg=BareGraph[ng];
          If[!MemberQ[InCanonicalForm,gg],AppendTo[InCanonicalForm,gg]]
        ];
        If[van,
          If[Head[ng]===Times,ng=Last[ng]];
          ng={{0,ng}};
          ng=Hold @@ ng;
          ng=ReplacePart[ng,Times,{1,0}]
        ];
        ng
      ]
    ]
  )
]]

CFIteration[0,_]=0
CFIteration[c_?NumberQ * g_Graph, part_List] := c*CFIteration[g,part]
CFIteration[g_Graph,part_List] := Profile[CFIteration,Block[
  {VertexInfo,nlp,InfoTypes,npart},
  If[MemberQ[InCanonicalForm,BareGraph[g]], g,
Profile[FindNewPartition,
    rule=MapThread[Rule,{Range[nv],part}];
    VertexInfo=Thread[{
      part,
      BareVertices[g] /. rule
    }];
    nlp=Length[InfoTypes=Union[VertexInfo]];
];
    If[nlp==Last[part] && nlp<nv, CFArtificialSubdivision[g,part],
      npart=VertexInfo /. Thread[Rule[InfoTypes,Range[nlp]]];
      If[nlp<nv,
        CFIteration[Permute[g,SortingPermutation[npart]],Sort[npart]],
        Permute[g,SortingPermutation[npart]]
      ]
    ]
  ]
]]

CFArtificialSubdivision[g_Graph,part_List] :=
Profile[CFArtificialSubdivision,Block[
  {ming=1,BareMing,ord,NCMing,npart0,i=1,j,k,cg,pos,Barecg},
  While[Length[pos=Position[part,i]]==1,++i];
  npart0 = part /. (j_ :> j+1 /; j>=i);
  k=0;
  While[k<Length[pos] && (k==0 || !MemberQ[InCanonicalForm,BareMing]),
    npart=ReplacePart[npart0,i,pos[[++k]]];
    If[nlp+1<nv,
      cg=CFIteration[Permute[g,SortingPermutation[npart]],Sort[npart]],
      cg=Permute[g,SortingPermutation[npart]]
    ];
    Barecg=BareGraph[cg];
    Which[
    k==1,
      ming=cg;
      BareMing=Barecg,
    (ord=Order[Barecg,BareMing])==1,
      ming=cg;
      BareMing=Barecg,
    ord==0 && NumericalCoefficient[cg]!=NCMing,
      van=True
    ];
    NCMing=NumericalCoefficient[ming]
  ];
  ming
]]

I2H[g_Graph,i_Integer] := Block[
  {j,j1,j2,k,k1,k2,ng},
  {j,k}=List @@ g[[i]];
  If[j===k,0,
    {j1,j2}=List @@ DeleteCases[g[[j]],i];
    {k1,k2}=List @@ DeleteCases[g[[k]],i];
    Which[
    j1===k1 || j2===k2, 0,
    j1===k2 || j2===k1, g,
    True,
      ng=g;
      ng[[j]]=g[[j]] /. {j2->j1, j1->k1};
      ng[[k]]=g[[k]] /. {k1->k2, k2->j2};
      ng=ReplacePart[ng,g[[j2]] /. j->k,j2];
      ReplacePart[ng,g[[k1]] /. k->j,k1]
    ]
  ]
]

I2X[g_Graph,i_Integer] := Block[
  {j,j1,j2,k,k1,k2,ng},
  {j,k}=List @@ g[[i]];
  If[j===k,0,
    {j1,j2}=List @@ DeleteCases[g[[j]],i];
    {k1,k2}=List @@ DeleteCases[g[[k]],i];
    Which[
    j1===k2 || j2===k1, 0,
    j1===k1 || j2===k2, -g,
    True,
      ng=g;
      ng[[j]]=g[[j]] /. j1->k1;
      ng[[k]]=g[[k]] /. {k1->k2, k2->j1};
      ng=ReplacePart[ng,g[[j1]] /. j->k,j1];
      ReplacePart[ng,g[[k1]] /. k->j,k1]
    ]
  ]
]

Jump[g_Graph,i_Integer] := Block[
  {j,k,l=Length[g]},
  {j,k}=List @@ g[[i]];
  Join[
    ReplacePart[ReplacePart[g,Dashed[j,l+1],i],g[[k]] /. i -> l+2,k],
    Graph[OV[i,l+3,l+4],Dashed[k,l+5],Dashed[l+1,l+5],
      Dashed[l+1,l+5],OV[l+2,l+3,l+4]
    ]
  ]
]

Jump[g_Graph] := Jump[g,First[First[Position[g,_Dashed,{1}]]]]

AddUnivalent[g_Graph,uv_,i_Integer] := Block[
  {j,k,l=Length[g]},
  {j,k}=List @@ g[[i]];
  Join[
    ReplacePart[ReplacePart[g,Dashed[j,l+1],i],g[[k]] /. i -> l+2,k],
    Graph[OV[i,l+2,l+3],Dashed[k,l+1],uv[l+1]]
  ]
]

AddUnivalent[g_Graph,uv_] := AddUnivalent[
  g,
  uv,
  First[First[Position[g,_Dashed,{1}]]]
]

Graph2CC[g_Graph] := Block[
  {cc},
  cc=DeleteCases[g,_Dashed,{1}];
  cc=cc /. OV[i_,j_,k_] :> OV[i-1,j-1,k-1];
  CC[Count[cc,_OV]/2,0] @@ cc
]
Graph2CC[expr_] := expr /. g_Graph :> Graph2CC[g]

CC2Graph[CC[m_,_][vs___]] := Join[
  Graph @@ Table[Dashed @@ (3m+First /@ Position[{vs},i]),{i,0,3m-1}],
  Graph[vs] /. i_Integer :> i+1
]
CC2Graph[expr_] := expr /. CC[m_,e_][vs___] :> CC2Graph[CC[m,e][vs]]

Prism[m_Integer]:= CC2Graph[CC[m,0] @@ Join @@
  (
    (
      {
        OV[#~Mod~m,(#+1)~Mod~m,m+#~Mod~m],
        OV[2*m+#~Mod~m,m+#~Mod~m,2*m+(#+1)~Mod~m]
      }& /@ Range[m]
    ) /. OV[i_,j_,k_] :> OV[i+1,j+1,k+1]
  )
]

Ring[m_Integer] := Block[
  {m3=3*m,i},
  Join @@ Table[
    Graph[UV[i+2],Dashed[i+2,i+5],OV[i,i-2,i+1]],
    {i,0,3m-3,3}
  ] /. i_Integer :> 1+(i~Mod~m3)
]

Ring[list_List] := Block[
  {m3=3*Length[list],i},
  Join @@ Table[
    Graph[(list[[1+i/3]])[i+2],Dashed[i+2,i+5],OV[i,i-2,i+1]],
    {i,0,m3-3,3}
  ] /. i_Integer :> 1+(i~Mod~m3)
]

Tree[list_List] := ReplacePart[
  Join[
    Graph[list[[1]][3]],
    Join @@ Table[
      Graph[list[[i]][3i-3],OV[3i-5,3i-4,3i-2],Dashed[3i-3,3i]],
      {i,2,Length[list]-1}
    ]
  ],
  Last[list][3Length[list]-6],
  -1
]

ThetaGraph[l1_:{},l2_:{},l3_:{}] := Block[
  {g,ll1,ll2,ll3},
  ll1=If[Head[l1]===List,l1,{l1}];
  ll2=If[Head[l2]===List,l2,{l2}];
  ll3=If[Head[l3]===List,l3,{l3}];
  g=Graph[Dashed[4,5],Dashed[4,5],Dashed[4,5],OV[1,2,3],OV[1,2,3]];
  g=Fold[AddUnivalent[#1,#2,1]&,g,ll1];
  g=Fold[AddUnivalent[#1,#2,2]&,g,ll2];
  Fold[AddUnivalent[#1,#2,3]&,g,ll3]
]

Phi[s_] := s /. {CC[m_, e_][vs___] :> RunThrough["PhiOf", CC[m, e][vs]]}

Laplacian[g_Graph] := Block[
  {lap,i,j},
  lap=DiagonalMatrix[Length /@ List @@ g];
  Do[
    Do[
      --(lap[[i,g[[i,j]]]]),
      {j,1,Length[g[[i]]]}
    ],
    {i,1,Length[g]}
  ];
  lap
]

AdjacencyMatrix[g_Graph] := Block[
  {am,i,j},
  am=Table[0,{Length[g]},{Length[g]}];
  Do[
    Do[
      ++(am[[i,g[[i,j]]]]),
      {j,1,Length[g[[i]]]}
    ],
    {i,1,Length[g]}
  ];
  am
]

AdjacencyMatrix2Graph[mat_?MatrixQ] := Graph @@ Table[
  UOV @@ Flatten[Position[mat[[i]],1]],
  {i,1,Length[mat]}
]

End[]
EndPackage[]
