(* :Title: Symbolica *)

(* :Author: Frank Zizza
	    Willamette University
	    zizza@willamette.edu
*)

(* :Summary:
   A package for computer experiments in symbolic dynamics. 
*)

(* :Context: Graphs` *)

(* :Package Version: 1.2 *)

(* :Copyright: Copyright 1993 *)

(* :History:
  Created by Frank Zizza at Willamette University, Summer 1993
*)

(* :Keywords: 
  Graph, LabeledGraph, Adjacency matrix, State splitting, 
  Irreducible, Resolving, Essential, Zeta functions.
*)

(* :Source:
  Adler, Lind & Marcus "Symbolic Dynamics". 
*)


(* :Mathematica Version: 2.1*)

(*****************************************************************************)


BeginPackage["Symbolica`"]


(***** USAGE MESSAGES *****)


Edge::usage = 

	"Edge[a, b] is a directed edge in a graph.  The first \
	argument is the source vertex and the second argument \
	is the target vertex.  Vertices can be any expression.  \
	Edge[a, b, c] is a labeled edge in a labeled graph.  \
	The first two arguments are the source and target vertices \
	and the last argument is the label of the edge."



Graph::usage = 

	"Graph[vertexlist, edgelist] is the primitive graph object.  \
	Graph objects should be created with MakeGraph.";




LabeledGraph::usage = 

	"LabeledGraph[vertexlist, edgelist] is the primitive labeled 
	graph object.  LabeledGraph objects should be created with \
	MakeLabeledGraph.";	




GetVertices::usage = 

	"GetVertices[{Edge[_, _, ___] ..}] calculates the list of \
	vertices in a list of edges."



SubListQ::usage = 

	"SubListQ[list1, list2] is a Boolean function that returns \
	True if and only if list1 is a sublist of list2."



UniqueQ::usage = 

	"UniqueQ[list1] is a Boolean function that returns True \
	if and only if there are no repeated elements in list1."



NonNegativeIntegerQ::usage = 

	"NonNegativeIntegerQ[int] is a Boolean function that returns \
	True if and only if int is a non-negative integer."



SNNIMQ::usage = 

	"SNNIMQ[ m ] is a Boolean function that returns True if and \
	only if m is a square matrix with non-negative integer entries."




MakeGraph::usage = 

	"MakeGraph constructs a directed graph.  There are several \
	alternatives using MakeGraph to construct a Graph object.\n\n\
									\
	MakeGraph[ list of edges, list of vertices ] constructs a Graph \
	object with the given edges and vertices.  Edges must be entered \
	as Edge[ source vertex, target vertex ] and do not need to be \
	distinct.  Vertices can be any Mathematica expression.  Also, \
	the list of vertices may contain vertices not found in any of \
	the edges in the list of edges.\n\n\
									\
	MakeGraph[ list of edges ] constructs a Graph object exactly as \
	above except the list of vertices are computed from the arguments \
	of all Edge objects found in list of edges.\n\n\
									\
	MakeGraph[ A = square non-negative integer matrix ] constructs a \
	Graph object with  A[[i, j]] edges from vertex i to vertex j.  \
	In this construction, vertices are labeled 1,2,3,... unless the \
	option DefaultVertexName is set to something other than the \
	default, Identity.\n\n\
									\
	MakeGraph[ForbiddenWords[___List]] constructs a graph whose \
	associated shift space has as the complement of its language, \
	the arguments to ForbiddenWords.  The arguments to ForbiddenWords \
	must be lists of symbols. Example: The construction of the golden \
	ratio shift space results from the input \
	MakeGraph[ ForbiddenWords[ {1,1} ] ].  \
	The default set of available symbols are {0, 1}.  Other symbol \
	sets can be used by setting the option DefaultSymbols.  \
	ForbiddenWords automatically eliminates arguments (words) \
	that contain other arguments (words).\n\n\
								\
	Graph objects can be displayed with Show.";




DefaultVertexName::usage = 

	"DefaultVertexName is an option for MakeGraph used when \
	MakeGraph is supplied with the adjacency matrix of a graph.  \
	It specifies the function to use in a call to Array that generates \
	the list of vertex names in the resulting Graph object.  The \
	default for this option is DefaultVertexName -> Identity which \
	creates vertices of the form 1, 2, 3,..  Another sensible \
	possibility is DefaultVertexName -> v which would result in vertices \
	labeled  v[1], v[2], v[3],...";



AllWordsLength::usage = 

	"AllWordsLength[n] returns a list of all possible words \
	of length n from the default alphabet {0, 1}.  \
	AllWordsLength[n, DefaultSymbols -> {a, b, c, ..}] returns \
	the list of all possible words of length n from the alphabet \
	{a, b, c, ...}.  The default alphabet can be reset in the \
	options of MakeGraph."


DefaultSymbols::usage =

	"DefaultSymbols is an option for MakeGraph and \
	AllWordsLength that describes the alphabet from which \
	words can be created.  The default setting is \
	DefaultSymbols -> {0,1}."



Language::usage = 

	"Language[forbiddenwords, opts] returns the set of words that \
	do not appear in the list forbiddenwords.  The list of all words \
	is constructed by AllWordsLength."



ForbiddenWords::usage =

	"ForbiddenWords[{word1, word2, ...}] represents a list of \
	forbidden words in the language associated to a shift space.  \
	ForbiddenWords automatically reduces overlap in words.  For \
	example, ForbiddenWords[{1, 1}, {1, 1, 0}, {0, 0, 0}] reduces to \
	ForbiddenWords[{1, 1}, {0, 0, 0}]."




MakeLabeledGraph::usage = 

	"MakeLabeledGraph constructs a directed labeled graph.  There \
	are several alternatives using MakeLabeledGraph to construct a \
	LabeledGraph object.\n\n
									\
	MakeLabeledGraph[ list of labeled edges, list of vertices ] \
	constructs a LabeledGraph object with the given labeled edges and \
	vertices.  Edges must be entered as \
	Edge[ source vertex, target vertex, label ] and do not need to \
	be distinct.  Vertices and labels can be any Mathematica expression. \
	Also, the list of vertices may contain vertices not found in any of \
	the edges in the list of edges.\n\n
									\
	MakeLabeledGraph[ list of labeled edges ] constructs a LabeledGraph \
	object exactly as above except the list of vertices are computed \
	from the arguments of all Edge objects found in list of edges.\n\n\
									\
	MakeLabeledGraph[A = square matrix with list of symbols as entries] \
	constructs a LabeledGraph object with Length[A[[i, j]]] edges from \
	vertex i to vertex j labeled A[[i, j]].   In this construction, \
	vertices are labeled 1,2,3,... unless the option DefaultVertexName \
	is set to something other than the default, Identity. \n \n

	LabeledGraph objects can be displayed with Show.";




Edges::usage = 

	"Edges[graph or labeled graph] returns the list of edges of \
	the input graph object."



Vertices::usage = 

	"Vertices[graph or labeled graph] returns the list of vertices \
	of the input graph."



Graph::usage = 

	"Graph[edges, vertices] is the internal representation of \
	directed graphs.  Graph objects should be created using \
	MakeGraph."


LabeledGraph::usage = 

	"LabeledGraph[edges, vertices] is the internal representation of \
	directed labeled graphs.  LabeledGraph objects should be created \
	using MakeGraph."


GraphQ::usage = 

	"GraphQ[object] returns True if the object is a graph or a \
	labeled graph.  Otherwise it returns false."




Matrix::usage = 

	"Matrix[ Graph object ] calculates the adjacency matrix for the \
	given graph."





SymbolicMatrix::usage = 

	"SymbolicMatrix[ LabeledGraph object ] calculates the symbolic \
	adjacency matrix for the given labeled graph.  The matrix is \
	expressed as a square matrices with entries that are lists of the \
	symbols from the corresponding graph vertices."



OutSplit::usage = 

	"OutSplit[g] state splits the graph g.  User will be prompted for \
	the partitions of each of the outgoing edges from every state.  \
	Partitions of a set are entered using the indices indicating position \
	in the list of edges the list.  For example the partition {{a,c},{b}} \
	of {a,b,c} should be entered as {{1,3},{2}}."



InSplit::usage = 

	"InSplit[g] state splits the graph g.  User will be prompted for \
	the partitions of each of the incoming edges into every state.  \
	Partitions of a set are entered using the indices.  For example \
	the partition {{a,b},{c}} of {a,b,c} should be entered as \
	{{1,2},{3}}."




CompleteOutSplit::usage = 

	"CompleteOutSplit[g] state splits the graph g assuming the partition \
	for every set of outgoing edges is partitioned into singleton sets."



CompleteInSplit::usage = 

	"CompleteInSplit[g] state splits the graph g assuming the partition \
	for every set of incoming edges is partitioned into singleton sets."


DisjointUnion::usage = 

	"DisjointUnion[g1, g2] constructs a graph that represents the \
	disjoint union of the two graphs g1 and g2.  Indices 1 and 2 \
	are used as suffixes to disjointify the set of states and vertices."



GraphProduct::usage = 

	"GraphProduct[graph1, graph2] constructs the cartesian product of \
	the two graphs."



LabelProduct::usage = 

	"LabelProduct[graph1, graph2] constructs the fiber product of the \
	two graphs."



SubsetConstruction::usage = 

	"SubsetConstruction[g] constructs a resolving graph g' whose states \
	are the set of non-empty subsets of the states of the labeled graph \
	g with an edge labeled a from state S to state T if there are edges \
	labeled a from every state of g in S to states in g contained in T."


ResolvingQ::usage = 

	"ResolvingQ[g] determines if for every vertex v in \
	the graph g, the edges leaving v have unique labels."



EventuallyResolvingQ::usage = 

	"EventuallyResolvingQ[labeledgraph, D] determines if \
	labeledgraph is eventually resolving with resolving delay \
	D.  EventuallyResolvingQ[labeledgraph] computes \
	EventuallyResolvingQ[labeledgraph, D] for successive D \
	(0 <= D <= number of vertices of labeledgraph) until True \
	returns."



IrreducibleQ::usage = 

	"IrreducibleQ[matrix] tests for irreducible matrices."



StrandedVertices::usage = 

	"StrandedVertices[g] returns a list of all stranded vertices \
	of the graph g."



EssentialQ::usage = 

	"EssentialQ[g] determines if the graph g has any stranded vertices."



Essentialize::usage = 

	"Essentialize[g] removes stranded vertices until no stranded \
	vertices remain."



ExteriorPower::usage = 

	"ExteriorPower[ labeledgraph, k] calculates the kth exterior power \
	of the labeledgraph.  This construction is used in the calculation \
	of the zeta function corresponding to the shift space associated to \
	a labeled graph.  k must be less than or equal to the number of \
	vertices in labeled graph and labeled graph must be resolving."
	


ZetaFunction::usage = 

	"ZetaFunction[graph or labeled graph, t] calculates the Zeta \
	function for graphs and labeled graphs in terms of the variable \
	t.  For labeled graphs, ZetaFunction uses an algorithm that \
	involves the exterior power construction of a resolving labeled \
	graph."


NEntropy::usage = 

	"NEntropy[graph] computes the logarithm base 2 of the spectral \
	radius of the adjacency matrix of graph.  NEntropy of resolving
	labeled graphs is computed from the underlying unlabeled graph.  \
	One can calculate the NEntropy of eventually resolving graphs as \
	NEntropy[UnderlyingGraph[labeledgraph]]."


UnderlyingGraph::usage = 

	"UnderlyingGraph[labeledgraph] returns the unlabeled graph \
	in labeledgraph."


SpectralRadius::usage = 

	"SpectralRadius[matrix] returns the maximum of the moduli of \
	the eigenvalues of matrix."

(*****************************************************************************)

Begin["`Private`"]
(* EDGE FUNCTIONS *)

Edge::NoVertices := 

	"An Edge was supplied with no vertices."


Edge::NoTarget =

	"An Edge was supplied with only one vertex: Edge[`1`].  \
	2 or 3 arguments are expected."


Edge::TooManyVertices = 

	"An Edge was supplied with too many arguments, Edge[`1`].  \
	2 or 3 arguments are expected."


Edge[] := 

	(Message[Edge::NoVertices]; HoldForm[Edge[]])


Edge[a_] := 

	(Message[Edge::NoTarget, a]; HoldForm[Edge[a]])


Edge[a__/; Length[{a}] > 3] :=

	(Message[Edge::TooManyVertices, a]; HoldForm[Edge[a]])


Format[Edge[a_, b_]] := 

	StringForm["<`1`,`2`>", a, b]


Format[Edge[a_, b_, c_]] := 

	StringForm["<`1`,`2`, `3`>", a, b, c]


(*****************************************************************************)

(* INTERNAL UTILITIES *)


GetVertices::NotEdgeList = 

	"GetVertices was not passed a list of Edges.";


GetVertices[elist:{__Edge}] := 

	Union[ Flatten[ Cases[ elist, Edge[a_, b_, c___] -> {a, b} ], 1 ] ]


GetVertices[___] := 

	(Message[GetVertices::NotEdgeList]; Abort[])


SubListQ::NotTwoLists = 

	"SubListQ was not passed two Lists."; 


SubListQ[list1_List, list2_List] := 

	TrueQ[ Complement[list1, list2] === {} ]


SubListQ[___] := 

	(Message[SubListQ::NotTwoLists]; Abort[])


UniqueQ[list_List] := 

	Signature[list] =!= 0


NonNegativeIntegerQ[n_] := 

	And[ n > -1, IntegerQ[n] ]


SNNIMQ[A_] := 

	And[ MatrixQ[A, NonNegativeIntegerQ], MatchQ[Dimensions[A],{n_, n_}] ]



(*****************************************************************************)


(* CREATING GRAPHS *)


MakeGraph::VertexError1 = 

	"List of vertices has duplications.  The list of vertices must \
	be a list of distinct expressions."


MakeGraph::VertexError2 = 

	"Supplied list of edges (first argument) contains vertices not \
	found in the supplied list of vertices (second argument)."


MakeGraph::BadMatrix = 

	"MakeGraph was expecting an adjacency matrix of a graph, but \
	it was not passed a square non-negative integer matrix."


MakeGraph::UseMakeLabeledGraph = 

	"Use MakeLabeledGraph in this context."


MakeGraph::What = 

	"MakeGraph cannot understand the input `1`";


Options[MakeGraph] = 

	{ DefaultVertexName -> Identity, 
          DefaultSymbols ->  {0, 1} };



(***** CREATING GRAPHS FROM EDGE LISTS *****)


MakeGraph[elist:{Edge[_,_] ...}, vlist_List/;Not[UniqueQ[vlist]]] :=  

	(Message[MakeGraph::VertexError1]; Abort[])


MakeGraph[elist:{Edge[_,_] ...}, vlist_List] := 
       
	Graph[elist, vlist] /;
                SubListQ[ GetVertices[elist], vlist ] && UniqueQ[vlist]


MakeGraph[elist:{Edge[_,_] ...}, vlist_List] :=  
	
	(Message[MakeGraph::VertexError2]; Abort[]) /; 
		Not[ SubListQ[ GetVertices[elist], vlist ] ]


MakeGraph[elist:{Edge[_,_] ...}] := 

	Graph[ elist, GetVertices[elist] ]



(***** CREATING GRAPHS FROM ADJACENCY MATRICES *****)


MakeGraph[A_?SNNIMQ, opts___Rule] := 

	Module[{n = Length[A], vert, dvn, edges},
		dvn = DefaultVertexName /. {opts} /. Options[MakeGraph]; 
		vert = Array[ dvn, n];
		edges = 
		  Table[Edge[vert[[i]],vert[[j]]],{i,n},{j,n},{k,A[[i, j]]}];
                Graph[ Flatten[edges], vert ] ]


MakeGraph[A:{{___}..}, opts___Rule] := 
        (Message[MakeGraph::BadMatrix]; Abort[])



(***** CREATING GRAPHS FROM FORBIDDEN WORD LISTS *****)


AllWordsLength[n_, opts___Rule] := 
       
	Module[{symbols},
                symbols = DefaultSymbols /. {opts} /. Options[MakeGraph];
                Distribute[ Table[ symbols , {n} ], List ] ]


Language[forbiddenwords_, opts___Rule] :=
       	
	Fold[DeleteCases,
                AllWordsLength[Max[Length /@ forbiddenwords], opts],
                forbiddenwords /. {a__Integer} -> {___, a, ___}]


ConstructEdge[{a_, b__}, {b__, c_}] := Edge[{a, b}, {b, c}]


ConstructEdge[{a_, b__}, {c__, d_}] := Sequence[]


SetAttributes[ForbiddenWords, {Orderless}]


ForbiddenWords[w1_, w2_, rest___] := 

	ForbiddenWords[w1, rest] /; SubListQ[w1, w2]


MakeGraph[ForbiddenWords[words___], opts___Rule] :=
        
	With[{lang = Language[{words}, opts]}, 
        MakeGraph[ Distribute[{lang, lang}, List, List, List, ConstructEdge]]]



(***** UNUSUAL CONDITIONS *****)


MakeGraph[{}] := Graph[{}, {}]


MakeGraph[{}, verts_] := Graph[{}, verts]


MakeGraph[{Edge[_,_,_] ..}, ___] := 
	
	Message[MakeGraph::UseMakeLabeledGraph]


MakeGraph[X___] := 

	(Message[MakeGraph::What, X]; Abort[])


(****************************************************************************)

(***** CREATING LABELED GRAPHS *****)


MakeLabeledGraph::VertexError = 

	"Supplied list of edges (first argument) contains vertices not \
	found in the supplied list of vertices (second argument)."


MakeLabeledGraph::EdgeError = 

	"MakeLabeledGraph was supplied with an Edge with no label: `1`."


MakeLabeledGraph::BadMatrix = 

	"MakeGraph was expecting an adjacency matrix of a labeled graph, \
	but it was not passed a square matrix with lists of symbols for \
	entries."


MakeLabeledGraph::What = 

	"MakeLabeledGraph cannot understand the input `1`";


Options[MakeLabeledGraph] = 

	{ DefaultVertexName -> Identity };

(***** CREATING LABELED GRAPHS FROM EDGE AND VERTEX LISTS *****)


MakeLabeledGraph[elist:{Edge[_, _, _] ...}, vlist_List] := 
        
	LabeledGraph[elist, vlist] /;
                SubListQ[ GetVertices[elist], vlist ]


MakeLabeledGraph[elist:{Edge[_,_,_] ...}, vlist_List] :=  

	(Message[MakeLabeledGraph::VertexError]; Abort[])


MakeLabeledGraph[elist:{___, bad:Edge[_,_], ___}, vlist_List] :=  

	(Message[MakeLabeledGraph::EdgeError, bad]; Abort[])


(*****  CREATING LABELED GRAPHS FROM EDGE LISTS *****)


MakeLabeledGraph[elist:{Edge[_, _, _] ...}] := 

	LabeledGraph[ elist, GetVertices[elist] ]


MakeLabeledGraph[elist:{___, bad:Edge[_,_], ___}] :=  

	(Message[MakeLabeledGraph::EdgeError, bad]; Abort[])



(*****  CREATING LABELED GRAPHS FROM SYMBOLIC ADJACENCY MATRICES *****)


MakeLabeledGraph[A:{{{___}..}..}, opts___Rule] := 

	Module[{n = Length[A], vert, dvn, edges},

		dvn = DefaultVertexName /.{opts} /. Options[MakeLabeledGraph];

		vert = Array[ dvn, n];

		edges = Table[ 
				Edge[ vert[[i]], vert[[j]], A[[i, j, k]]], 
			      {i, n}, {j, n}, {k,Length[A[[i, j]]]}];

                LabeledGraph[ Flatten[edges], vert ] ] /; 
			MatchQ[Dimensions[A],{n_, n_, ___}]


MakeLabeledGraph[A:{{___}..}, opts___Rule] := 

	(Message[MakeLabeledGraph::BadMatrix]; Abort[])


(***** UNUSUAL CONDITIONS *****)


MakeLabeledGraph[{}] := LabeledGraph[{}, {}]


MakeLabeledGraph[{}, verts_] := LabeledGraph[{}, verts]


MakeLabeledGraph[X___] := 
	(Message[MakeLabeledGraph::What, X]; Abort[])



(*****************************************************************************)

(***** GRAPH UTILITIES *****)


Edges::NotGraphObject = 

	"Input of Edges was not a graph."


Edges[g_?GraphQ] := First[g]



Edges[___] := 

	(Message[Edges::NotGraphObject]; Abort[])


Vertices::NotGraphObject = 

	"Input to Vertices was not a graph or labeled graph."


Vertices[g_?GraphQ] := Last[g]


GraphQ[g_] := TrueQ[Head[g] == Graph || Head[g] == LabeledGraph]


Matrix::NotGraph = 

	"Argument to Matrix was not a graph."


Matrix::UsingSymbolicMatrix = 

	"SymbolicMatrix has been called to calculate the symbolic adjacency \
	matrix of the labeled graph input.  The result is a matrix of lists \
	where each list represents the labels of edges between the \
	corresponding pair of vertices."


Matrix[ Graph[edges_, vertices_] ] := 
	
	With[{n = Length[vertices]}, 
		Table[ 
			Count[edges, Edge[vertices[[i]], vertices[[j]]] ],
		      {i, n},{j, n}]]


Matrix[ LabeledGraph[arg___] ] := 

	(Message[Matrix::UsingSymbolicMatrix]; 
	SymbolicMatrix[ LabeledGraph[ arg ] ])


Matrix[___] := (Message[Matrix::NotGraph]; Abort[])


SymbolicMatrix::NotLabeledGraph = 

	"Argument to SymbolicMatrix was not a labeled graph."


SymbolicMatrix[ LabeledGraph[edges_, vert_] ] := 
	
	With[{n = Length[vert]},
		Table[ 
			Cases[edges, Edge[vert[[i]], vert[[j]], a_] -> a ],
		      {i, n},{j, n}] ]


SymbolicMatrix[___] := 

	(Message[SymbolicMatrix::NotLabeledGraph]; Abort[])



(***** GRAPHICS FOR GRAPHS AND LABELED GRAPHS *****)


Index[l_List] := 

	Thread[{l, Range[Length[l]]}]


Magnitude[X_] := Sqrt[Plus @@ (X^2)]


circ[{x0_, y0_}, h_] := 

	Module[{r = 0.35, arr} , 
	    arr = ((1 + 2 r h){x0, y0} + #)& /@ 
		(rot[ArcTan[x0,y0] - 1.65] /@ arrow); 
	    {
		Circle[(1 + r h) {x0 , y0}, r h],
		GrayLevel[1],
	     	Polygon[arr],	
	     	GrayLevel[0],
		Line[arr]
	     }
	    ]


circ[{x0_, y0_}, h_, A_] := 

	Module[{r = 0.35, arr} , 
	    arr = ((1 + 2 r h){x0, y0} + #)& /@ 
		(rot[ArcTan[x0,y0] - 1.65] /@ arrow); 
	    {
		Circle[(1 + r h) {x0 , y0}, r h],
		GrayLevel[1],
	     	Polygon[arr],	
	     	GrayLevel[0],
		Line[arr],
	     	Text[A,(1 + 2 r h) {x0 , y0}]
	     }
	    ]


rule1 = {Edge[a_, a_, A___], n_} -> circ[a, n/2.5, A]


rule2 = {Edge[a_, b_, C___], n_} -> arc[a, b, n/6, C]


arrowsize = 0.15;


arrow = 
	arrowsize * 
	N[{{-1/3,1/2},{-1/2,3^(1/2)/2},{1,0},{-1/2,-3^(1/2)/2},{-1/3, -1/2}}];


rot[t_][{x_,y_}] := {{Cos[t], -Sin[t]}, {Sin[t], Cos[t]}} . {x,y};


arc[p:{_,_}, q:{_,_}, h_?Positive] := 
	
	Module[{p1 = N[p],p2 = N[q],p3,eq,R,x0,y0,r,p0,v1,v2,t1,t2},
		
		p3 = (p1 + p2)/2 + 
			h {-1, 1} Reverse[p2 - p1]/Magnitude[p2 - p1];

		eq[pt_] := (Plus @@ ((pt - {x0, y0})^2)) == R^2;

		eqs = eq /@ {p1, p2, p3};

		{x0, y0, r} = 
			First[ Cases[ {x0, y0, R} /. 
			 Solve[eqs, {x0, y0, R}], {_, _, _?Positive}] ];

		p0 = {x0, y0}; v1 = p1 - p0; v2 = p2 - p0;

		t1 = ArcTan @@ v2;

		t2 = t1 + ArcCos[(v1 . v2)/(Magnitude[v1]*Magnitude[v2])];

		arr = (p3 + #)& /@ (rot[ArcTan @@ (p2 - p1)] /@ arrow);

		{
		Circle[p0, r, {t1, t2}], 
		GrayLevel[1],
		Polygon[ arr ],
		GrayLevel[0],
		Line[ arr ]
		}

		]


arc[p:{_,_}, q:{_,_}, h_?Positive, C_] := 
	
	Module[{p1 = N[p],p2 = N[q],p3,eq,R,x0,y0,r,p0,v1,v2,t1,t2, arr},
		
		p3 = (p1 + p2)/2 + 
			h {-1, 1} Reverse[p2 - p1]/Magnitude[p2 - p1];

		eq[pt_] := (Plus @@ ((pt - {x0, y0})^2)) == R^2;

		eqs = eq /@ {p1, p2, p3};

		{x0, y0, r} = First[ Cases[ {x0, y0, R} /. 
			Solve[eqs, {x0, y0, R}], {_, _, _?Positive}] ];

		p0 = {x0, y0};

		v1 = p1 - p0;

		v2 = p2 - p0;

		t1 = ArcTan @@ v2;

		t2 = t1 + ArcCos[(v1 . v2)/(Magnitude[v1]*Magnitude[v2])];

		arr = (p3 + #)& /@ (rot[ArcTan @@ (p2 - p1)] /@ arrow);

		{
		Circle[p0, r, {t1, t2}], 
		GrayLevel[1],
		Polygon[ arr ],
		GrayLevel[0],
		Line[ arr ],
		Text[ToString[C], p3]
		}

		]


show[edges_, v_, labs_, opts___Rule] := 

	Module[{e,n = Length[v], p, t = Function[{x,y},Text[ToString[x], y]]},
		
		p = Chop[Drop[#, -1]&[
			Table[ {Cos[theta], Sin[theta]}, 
			   {theta, N[Pi+0.001], N[-Pi+0.001], N[-2Pi/n]}]]];

		e = Table[ Cases[ edges, Edge[v[[i]], v[[j]], ___] ], 
			{i, n}, {j, n}];

		e = Map[ Index, e, {2} ];

		e = Flatten[e, 2];

		e = e /. Thread[v -> p];

		e = e /. rule1 /. rule2;

		Show[Graphics[
			{
			Thickness[0.003],
			e, 
			GrayLevel[1],
			Disk[#, 0.12]& /@ p,
			GrayLevel[0],
			Thickness[0.004],
			Circle[#, 0.12]& /@ p,
			MapThread[t, {labs, p}]
			}
		     ], opts,
		    AspectRatio -> 1, PlotRange -> {{-2, 2},{-2, 2}}]
		
		]



Graph /: Show[g_Graph, opts___Rule] := 
	Module[{g1 = MakeGraph[Matrix[g], DefaultVertexName -> "v"]},
	show[Edges[g1], Vertices[g1], Vertices[g], opts]]



LabeledGraph /: Show[g_LabeledGraph, opts___Rule] := 
	Module[{g1 = MakeLabeledGraph[SymbolicMatrix[g], 
		DefaultVertexName -> "v"]},
	show[Edges[g1], Vertices[g1], Vertices[g], opts]]



(****************************************************************************)

(***** State Splitting Algorithms *****)


OutSplit[g_?GraphQ] := 
	
	Module[{partition, edges, newvertices}, 
		partition = Map[PartitionOutgoingEdges[g, #]&, Vertices[g] ];
		newvertices = Apply[ Array, 
			{#[[1,1,1]], Length[#]}& /@ partition, {1}];
		edges = MapIndexed[IndexFirst, #]& /@ partition;
		edges = OutDuplicate[#, newvertices]& /@ Flatten[edges];
		Head[g][Map[CollectIndices, Flatten[edges], {2}],
			CollectIndices /@ Flatten[newvertices]]
		]


InSplit[g_?GraphQ] := 

	Module[{partition, edges, newvertices}, 
		partition = Map[PartitionIncommingEdges[g, #]&, Vertices[g] ];
		newvertices = Apply[ Array, 
			{#[[1,1,2]], Length[#]}& /@ partition, {1}];
		edges = MapIndexed[IndexSecond, #]& /@ partition;
		edges = InDuplicate[#, newvertices]& /@ Flatten[edges];
		Head[g][Map[CollectIndices, Flatten[edges], {2}], 
			CollectIndices /@ Flatten[newvertices]]
		]


CompleteOutSplit[g_?GraphQ] := 

	Module[{partition, edges, newvertices}, 
		partition =  Map[CompletePartitionOutgoingEdges[g, #]&, 
					Vertices[g] ];
		newvertices = Apply[ Array, 
			{#[[1,1,1]], Length[#]}& /@ partition, {1}];
		edges = MapIndexed[IndexFirst, #]& /@ partition;
		edges = OutDuplicate[#, newvertices]& /@ Flatten[edges];
		Head[g][ Map[CollectIndices, Flatten[edges], {2}], 
			CollectIndices /@ Flatten[newvertices]]
		]


CompleteInSplit[g_?GraphQ] := 

	Module[{partition, edges, newvertices}, 
		partition =  Map[CompletePartitionIncommingEdges[g, #]&, 
					Vertices[g] ];
		newvertices = Apply[ Array, 
			{#[[1,1,2]], Length[#]}& /@ partition, {1}];
		edges = MapIndexed[IndexSecond, #]& /@ partition;
		edges = InDuplicate[#, newvertices]& /@ Flatten[edges];
		Head[g][Map[CollectIndices, Flatten[edges], {2}], 
			CollectIndices /@ Flatten[newvertices]]
		]


PartitionEdges::NotPartition = 

	"PartitionEdges was not supplied with a partition of the set of \
	edges."
		

PartitionOutgoingEdges[g_?GraphQ, I_] := 

	Module[{edges, n, partition},
		edges = Cases[ Edges[g], Edge[I, __]]; 
		If[(n = Length[edges]) == 1, {edges}, 
		partition = Input[ StringJoin[
			"Partition the ", 
			ToString[n], 
			" edges:\n", 
			ToString[edges], "\n"] ];
		If[Sort[Flatten[partition]] == Range[n],
    		Map[ Part[edges, #]&, partition] ],
		Message[PartitionEdges::NotPartition]]]
    		

CompletePartitionOutgoingEdges[g_?GraphQ, I_] := 

	Module[{edges, partition},
		edges = Cases[ Edges[g], Edge[I, __]]; 
		Partition[edges,1]]
		    		

IndexFirst[{Edge[i_, j_, a___]}, {n_}] := {Edge[i[n], j, a]}


IndexFirst[list_List, {n_}] := 

	{IndexFirst[{First[list]}, {n}], IndexFirst[Rest[list], {n}]}


OutDuplicate[Edge[i_, j_, a___], newvertices_] :=
 
	Map[ Edge[i, #, a]&,  Cases[ Flatten[newvertices], j[_] ]  ]


PartitionIncommingEdges[g_?GraphQ, I_] := 

	Module[{edges, partition},
		edges = Cases[ Edges[g], Edge[_, I, a___]];
		If[Length[edges] === 1, {edges},
		partition = Input[ StringJoin[
			"Partition the edges:\n", 
			ToString[edges], "\n"] ];
    		Map[ Part[edges, #]&, partition]]]


CompletePartitionIncommingEdges[g_?GraphQ, I_, a___] := 

	Module[{edges, partition},
		edges = Cases[ Edges[g], Edge[_, I, ___]];
		Partition[edges,1]]


IndexSecond[{Edge[i_, j_, a___]}, {n_}] := {Edge[i, j[n], a]}


IndexSecond[list_List, {n_}] := 

	{IndexSecond[{First[list]}, {n}], IndexSecond[Rest[list], {n}]}


InDuplicate[Edge[i_, j_, a___], newvertices_] :=
 
	Map[ Edge[#, j, a]&,  Cases[ Flatten[newvertices], i[_] ]  ]

CollectIndices[a_[n__][m_]] := CollectIndices[a[n, m]]

CollectIndices[x_] := x

(****************************************************************************)

(***** GRAPH CONSTRUCTIONS *****)

DisjointUnion[g1_Graph, g2_Graph] := 

	MakeGraph[ 
	   Join[Map[#[1]&, Edges[g1], {2}],Map[#[2]&, Edges[g2], {2}]],
	   Join[Map[#[1]&, Vertices[g1], {1}], Map[#[2]&, Vertices[g2], {1}]]]


DisjointUnion[g1_LabeledGraph, g2_LabeledGraph] := 

	MakeLabeledGraph[ 
	    Join[Map[#[1]&, Edges[g1],{2}],Map[#[2]&, Edges[g2], {2}]],
	    Join[Map[#[1]&, Vertices[g1],{1}], Map[#[2]&, Vertices[g2], {1}]]]


GraphProduct[ g1_Graph, g2_Graph] := 

	Module[{edges}, 
		edges = CartesianProduct[Edges[g1], Edges[g2]];
		edges = Thread[#, Edge]& /@ edges;
		MakeGraph[edges, CartesianProduct[Vertices[g1],Vertices[g2]]]]


GraphProduct[ g1_LabeledGraph, g2_LabeledGraph] := 

	Module[{edges}, 
		edges = CartesianProduct[Edges[g1], Edges[g2]];
		edges = Thread[#, Edge]& /@ edges;
		MakeLabeledGraph[edges, 
			CartesianProduct[Vertices[g1], Vertices[g2]]]]


LabelProduct[ g1_LabeledGraph, g2_LabeledGraph] :=

	Module[{graph1, edges},
	graph1 = GraphProduct[g1, g2];
	edges = Cases[Edges[graph1],
		 Edge[x_, y_, OrderedPair[a_, a_]] -> Edge[x, y, a]];
	MakeLabeledGraph[ edges ]]


SubsetConstruction[g_LabeledGraph] := 

	Module[{vertices, alphabet, edges},
		vertices = ReducedPowerSet[ Vertices[g] ];
		alphabet = Alphabet[g];
		edges = CartesianProduct[vertices, alphabet];
		edges = edges /. OrderedPair[T_, a_] :> Edge[T, Ta[T,a,g], a];
		edges = DeleteCases[edges, Edge[_, set[], _]];
		MakeLabeledGraph[edges, vertices]]

SubsetConstruction[g_Graph] :=

	(Message[SubsetConstruction::NotLabeledGraph]; HoldForm[SubsetConstruction[g]])


SubsetConstruction::NotLabeledGraph = "The subset construction applies only to labeled graphs."


Ta[T_set, a_, g_LabeledGraph] :=
 
	Apply[set, 
		Union[ #[[2]]& /@ 
			Select[ Edges[g], 
				(MemberQ[T, #[[1]]] && #[[3]] === a)&] ] ]


Alphabet[g_LabeledGraph] := 

	Union[ Last /@ Edges[g] ] 


PowerSet[list_List] := 

	Distribute[({set[], set[#]}&) /@ list, List, List, List, Join]


Format[set[a___]] := {a}


ReducedPowerSet[G_List] := Complement[PowerSet[G], {set[]}]


CartesianProduct[A_List, B_List] := 

	Distribute[{A, B}, List, List, List, OrderedPair]


Format[ OrderedPair[a_, b_] ] := StringForm["(``,``)", a, b]

(****************************************************************************)

(***** BOOLEAN FUNCTIONS FOR GRAPHS AND LABELED GRAPHS *****)


PositiveQ[A_?MatrixQ] := 

	Apply[ And, Positive[A], {0,1}]



IrreducibleQ[A_?MatrixQ] := 

	PositiveQ[ Sum[ MatrixPower[A, i], {i, 1, Length[A]} ] ]


IrreducibleQ[g_Graph] := 

	IrreducibleQ[Matrix[g]]


IrreducibleQ[g_LabeledGraph] := 

	IrreducibleQ[UnderlyingGraph[g]]


UnderlyingGraph[g_LabeledGraph] := 
	
	Module[{edges},
		edges = Drop[#, -1]& /@ Edges[g];
		MakeGraph[edges, Vertices[g]]]
		

EdgesFrom[vertex_, g_?GraphQ] := 

	Cases[ Edges[g], Edge[vertex, __] ]
	

ResolvingQ[g_LabeledGraph] := 

	Apply[ And, 
		Map[ UniqueQ, 
			Map[ #[[3]]&, 
				Map[ EdgesFrom[#, g]&, 
					Vertices[g] ], {2}]]]


(*****************************************************************************)

(***** REMOVING STRANDED VERTICES *****)



StrandedQ[vertex_, g_?GraphQ] := 
	Module[{edges}, 
		edges  = Edges[g];
		Count[edges,Edge[vertex, __]] == 0 || 
		Count[edges,Edge[_, vertex, ___]] == 0]

StrandedVertices[g_?GraphQ] := 

	Select[ Vertices[g], StrandedQ[#, g]&]


RemoveVertex[g_?GraphQ, vertex_] := 
	
	Head[g][
	DeleteCases[ 
		Edges[g], Edge[vertex, ___] | Edge[_, vertex, ___] ] ,
	Complement[Vertices[g], {vertex}]]


EssentializeStep[g_?GraphQ] := 

	Fold[ RemoveVertex, g, StrandedVertices[g]]



Essentialize[g_?GraphQ] := 

	FixedPoint[ EssentializeStep, g]



EssentialQ[g_?GraphQ] := 

	StrandedVertices[g] === {}


(***************************************************************************)

(***** ZETA FUNCTIONS *****)


Indices[n_Integer, k_Integer] :=
 
	Module[{var, itt},

		var = Table[Unique[i], {k}];

		itt = {{var[[1]],1, n-k+1}}    ~Join~
		Rest[ Transpose[{var,RotateRight[var]+1,Range[n-k+1, n]}]];

		Flatten[Table[ Evaluate[var],Evaluate[Sequence @@ itt]], k-1]]
	


ExteriorVertices[vertices_, power_Integer] := 

	Part[vertices, #]& /@ Indices[Length[vertices], power]
	

Follower[G_, i_, a_] := 

	Cases[ Edges[G], Edge[i, j_, a] -> j ]


Follower[G_, list_List, a_] := 

	Flatten[ Map[Follower[G, #, a]&, list] /. {} -> Null, 1]


ExteriorEdge[G_, OrderedPair[list_, a_]] :=  
	
	Module[{vertex},
		vertex = Follower[G, list, a];
		If[Signature[vertex] =!= 0,
			Edge[list, Sort[vertex],
				{Signature[list]*Signature[vertex],  a}],
			Null]]

ExteriorPower::notdfa = 

	"The exterior power of a labeled graph is only defined \
	for resolving labeled graphs."

ExteriorPower::cnt = 

	"The k-th exterior power of a labeled graph is only defined \
	for k less than or equal to the number of vertices in the \
	labeled graph."


ExteriorPower[G_LabeledGraph, k_Integer] :=
	 Message[ExteriorPower::notdfa] /; !ResolvingQ[G]



ExteriorPower[G_LabeledGraph, k_Integer] :=
	Message[ExteriorPower::cnt] /; (k > Length[Vertices[G]])



ExteriorPower[G_LabeledGraph, k_Integer] :=
	Module[{vert, edges},
		vert = ExteriorVertices[Vertices[G], k];
		edges = CartesianProduct[vert, Alphabet[G]];
		edges = Map[ ExteriorEdge[G, #]&, edges ];
		edges = Select[edges, FreeQ[#, Null]&];
		MakeLabeledGraph[edges, vert] ] /; 
			k <= Length[Vertices[G]] && ResolvingQ[G]



Projection[M_] := Apply[ Plus, Map[ First, M, {3}], {2} ]


ZetaFunction::notdfa = 

	"This algorithm for the zeta function of a labeled graph only works \
	for resolving labeled graphs.  You may construct an equivalent \
	resolving labeled graph using SubsetConstruction and then calculate \
	the Zeta function of that labeled graph, however, the zeta function \
	algorithm is exponential in the number of vertices."


ZetaFunction[G_LabeledGraph, t_] := 
	
	Module[{n = Length[Vertices[G]], exteriorpowers, matrices, poly, exp},

	     If[ Not[ ResolvingQ[G] ], 
	   	
		(*then*)Message[ZetaFunction::notdfa],	
		
		(*else*)exteriorpowers = Table[ExteriorPower[G, k], {k,1,n}];
		
		matrices = Map[SymbolicMatrix, exteriorpowers];

		matrices = Map[Projection, matrices];

		poly =  Det[IdentityMatrix[Length[#]] - t*#]& /@ matrices;

		Product[ poly[[k]]^((-1)^k), {k, 1, n}]

	        ]

	       ]


ZetaFunction[G_Graph, t_] := 
	Module[{A, n},
		A = Matrix[G]; n = Length[A];
		1/Det[ IdentityMatrix[n] - t A ] ]



(***************************************************************************)

(*Eventually Resolving*)


PathsFrom[G_, I_, 1] := Map[EdgePath, Cases[Edges[G], Edge[I, __]]]

AddEdge[ep:EdgePath[___Edge, Edge[_, J_, ___]], e:Edge[J_, __]] := 
	Append[ep, e]

AddEdge[EdgePath[___Edge, Edge[_, J_, ___]], Edge[K_, __]] := 
	Sequence[]

PathsFrom[G_, I_, n_/;n>1] :=  
	Flatten[Outer[AddEdge, PathsFrom[G, I, n-1], Edges[G]]]

GetLabels[ep:EdgePath[___Edge]] := {First[ep], List @@ Last /@ ep}

LabelsFrom[G_, I_, n_/;n>1] := GetLabels /@ PathsFrom[G, I, n]

resolvingQ[G_LabeledGraph, I_, n_Integer/;n>1] := 
	Module[{ep, lab},
	ep = LabelsFrom[G, I, n];
	lab = Union[Last /@ ep];
	And @@ Apply[SameQ, Cases[ep, {_Edge, #}]& /@ lab, {1}]]

EventuallyResolvingQ[G_LabeledGraph, n_Integer/;n>0] := 
	UniqueQ[Edges[G]] && 
   	And @@ (resolvingQ[G,#,n+1]& /@ Vertices[G])

EventuallyResolvingQ[G_LabeledGraph, 0] := ResolvingQ[G]


EventuallyResolvingQ[G_LabeledGraph] := 
	Module[{i, n = Length[Vertices[G]]},
	For[i = 0, (i <= n), i++,
	 Print["Testing: Delay " <> ToString[i] <> "."];
	  If[ EventuallyResolvingQ[G, i], 
	    Print["Resolving delay is " <> ToString[i] <> "."]; Return[True],
	    Print["Resolving delay " <> ToString[i] <> " failed."]]]]



(***************************************************************************)


(*NEntropy*)


NEntropy::no = "NEntropy only works for labeled graphs that are at
	worst eventually resolving."


NEntropy[g_Graph] := 
	Log[2,SpectralRadius[Matrix[g]]]


NEntropy[g_LabeledGraph] := 
	Log[2,SpectralRadius[Matrix[UnderlyingGraph[g]]]]/;
		ResolvingQ[g]

NEntropy[g_LabeledGraph] := 
	Log[2,SpectralRadius[Matrix[UnderlyingGraph[g]]]]/;
		EventuallyResolvingQ[g]

NEntropy[g_LabeledGraph] := 
	Message[NEntropy::no, g]



(***************************************************************************)


(*Spectral Radius*)


SpectralRadius[m_?MatrixQ] :=  Max[Abs[Eigenvalues[N[m]]]]




(***************************************************************************)




End[]

EndPackage[]
