(*
   A 3-regular graph is given (usually the tetrahedron for starters).
   Three operations can take place: 
   	A vertex can expand into a triangular face.
	A triangular face adjacent to 3 vertices may contract into a vertex.
	An edge adjacent to 4 vertices may flip its orientation.
	
   Various rules can be defined governing when such operations occur 
   (see way below)

*)

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

ReplaceParts[expr_] := expr
ReplaceParts[expr_, {p_, new_}, rest___ ] := 
	ReplaceParts[ ReplacePart[expr, new, p], rest ]

(* choose a random element from a list *)
Unprotect[Random]
Random[{}]:={}
Random[l_List]:=l[[Random[Integer,{1,Length[l]}]]]
Protect[Random]
ChooseSome[l_List, n_Integer] /; Length[l] > n := 
	Drop[l, {Random[Integer, {1, Length[l]}]}]
ChooseSome[l_List, n_Integer] := l

(* divides l_ into Length[plist_] sublists by removing those 
	elements indexed by plist_ (by default) *)
SplitIntoPieces[l_List, plist_, offsets_:{1,-1}] :=
	If[Subtract@@# > 0, {}, Take[l,#]]& /@
		( (offsets+# &) /@ Partition[
		   Join[{1-First[offsets]},plist,{Length[l]-Last[offsets]}],
		2, 1] )
			
(* Random partition of n into p non-negative integers *)
RandomPartition[n_Integer,p_Integer] :=
	Length /@ SplitIntoPieces[#, First /@ Position[#,_Real]]& @
	    Sort[Join[Range[1,n],Table[0.5+Random[Integer,{1,n-1}],{p-1}]]]
	    
(* divides l_ randomly into sublists of minimum length *)
RandomPieces[l_List, n_Integer, minlen_:1] /; Length[l]>=minlen n := 
	SplitIntoPieces[ l, Rest[ FoldList[Plus, 0, 
		(minlen+#)& /@ Drop[RandomPartition[Length[l]-n minlen, n],-1]
	] ], {1,0} ]

RemoveDuplicates[{a___,b_,c___,d_,e___}, sameq_:SameQ] := 
	RemoveDuplicates[{a,b,c,e},sameq] /; sameq[b,d]
RemoveDuplicates[l_List, sameq_:SameQ] := l

Shrink[l_List, predq_:(OrderedQ[{#1,#2}]&)] := 
	Fold[ If[predq[Last[#1],#2] && predq[#2,Last[#1]], #1, Append[#1,#2]]&, 
		Take[#,1], Drop[#,1] ]& @ Sort[ l, predq] 

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

(* 
   Representation for a graph (needn't be 3-regular, but is here) is
    Graph[ Null | {
    		{ adjacent vertex index numbers (in clockwise order) },
		{ x, y } (both must always be real numbers)
	}, ... ]
   (must be careful of the Null entries!)
   
   vertices are represented by integers
   edges are represented by integer pairs
   	in unique canonical form, each is sorted (lowest number first)
	when direction is important, they are {from_, to_}
   faces are lists of integers, in clockwise order, and in
   	canonical form with the lowest vertex first
	(uniqueness for canonical form is unfortunately not yet
		guaranteed for big faces which contain their lowest
		entry more than once!)
   points are real number pairs
*)

(************* Sample 3-regular graphs ***************)

Tetrahedron = Graph[ { {2,3,4}, {0., 0.} },
		     { {1,4,3}, {N[Sqrt[3]/2], -.5} },
		     { {1,2,4}, {-N[Sqrt[3]/2], -.5} },
		     { {1,3,2}, {0., 1.} } ];
		     
Cube = Graph[ { {2,5,4}, {-2.,2.} },
	      { {1,3,6}, {2.,2.} },
	      { {2,4,7}, {2.,-2.} },
	      { {3,1,8}, {-2.,-2.} },
	      { {1,6,8}, {-1.,1.} },
	      { {5,2,7}, {1.,1.} },
	      { {6,3,8}, {1.,-1.} },
	      { {4,5,7}, {-1.,-1.} } ];

Prism = Graph[ { {2,3,4}, {1.,0.} },
	       { {1,5,3}, {3.,1.2} },
	       { {1,2,6}, {3.,-1.2} },
	       { {1,6,5}, {-1.,0.} },
	       { {2,4,6}, {-3.,1.2} },
	       { {3,5,4}, {-3.,-1.2} }  ];

Kite8 = Graph[	{{3, 2, 6}, {0.821177, 0.979568}}, 
		{{1, 5, 6}, {2.92976, 0.489963}}, 
 		{{1, 6, 7}, {-0.712016, 2.59254}}, 
		{{5, 7, 8}, {0.262078, 0.077012}}, 
 		{{4, 8, 2}, {1.71583, -1.61262}}, 
		{{1, 2, 3}, {1.71583, 2.59254}}, 
 		{{4, 3, 8}, {-1.92594, 0.489963}}, 
		{{4, 7, 5}, {-0.712016, -1.61262}}  ];

Asymm12A = Graph[	{{4, 10, 6}, {0.487894, 2.10071}}, 
			{{3, 7, 8}, {1.99066, -0.660749}}, 
 			{{2, 5, 10}, {0.762115, -0.989507}}, 
			{{1, 6, 7}, {3.10143, 3.91404}}, 
 			{{8, 12, 3}, {-1.06565, -4.73899}}, 
			{{1, 9, 4}, {-1.06565, 4.86515}}, 
 			{{2, 4, 8}, {4.95596, 0.0630824}}, 
			{{2, 7, 5}, {3.10143, -3.78788}}, 
 			{{6, 11, 12}, {-4.40738, 2.2002}}, 
			{{3, 11, 1}, {-0.080433, 0.50467}}, 
 			{{9, 10, 12}, {-1.90956, 0.0721451}}, 
			{{9, 11, 5}, {-4.40738, -2.07404}}   ];
			
Asymm12B = Graph[	{{10, 11, 12}, {1.59138, -0.0113995}}, 
			{{5, 10, 9}, {0.209812, 0.280379}}, 
  			{{11, 5, 6}, {2.20657, 2.60048}}, 
			{{12, 9, 10}, {-0.153754, -2.30078}}, 
 			{{3, 2, 7}, {0.575919, 1.1745}}, 
			{{3, 7, 8}, {-0.153754, 3.13921}}, 
			{{6, 5, 8}, {-0.411879, 1.40568}}, 
			{{6, 7, 9}, {-2.04659, 1.62972}}, 
  			{{4, 8, 2}, {-2.04659, -0.791298}}, 
			{{1, 4, 2}, {0.868398, -0.208094}}, 
 			{{1, 3, 12}, {3.25701, 0.419213}}, 
			{{1, 11, 4}, {2.20657, -1.76205}}    ];

(************* 2-space operations ************)

(* takes an ordered list of points, finds centroid of enclosed polygon*)
centroid[xylist:{_List,_List,_List}] := (Plus @@ xylist) / Length[xylist]
centroid[xylist_List] :=
	(( Plus @@ ( (Times @@ #)& /@ # ) ) / 
	 ( Plus @@ ( Last /@ # ) ) )& [
	({centroid[#], facearea[#]}& @ Append[#,{0.,0.}] )& /@
		Partition[Append[xylist,First[xylist]],2,1] ]
centroid[g_Graph, plist_List] := 
	centroid[Last /@ List@@g[[plist]]]
	
(* takes a list of points, finds the area of the enclosed polygon *)
(* positive area is returned for clockwise orientation *)
facearea[xylist_List] := 
	1/2 (Plus @@ ( First[#] . perp[Last[#]] & /@ 
		Partition[Append[xylist,First[xylist]],2,1] ) )
facearea[g_Graph, plist_List] := facearea[Last /@ List@@g[[plist]]]

(* return a point frac_ from from_ to to_ *)
midpoint[{from_List, to_List}, frac_:.5] := frac to + (1-frac) from
midpoint[g_Graph, {from_Integer, to_Integer}, frac_:.5] :=
	midpoint[ Last /@ List@@g[[{from,to}]], frac ]
	
(* returns a point for from_ were from-to rotated 90* cw about its midpoint *)
(* and expanded by frac_ *)
leftpoint[{from_List, to_List}, frac_:1.] := 
	midpoint[{from+ frac perp[to-from], to}]
leftpoint[g_Graph, {from_Integer, to_Integer}, frac_:.5] := 
	leftpoint[ Last /@ List@@g[[{from,to}]] ]

(* perpendicular vector to the left of vec *)
perp[vec_List] := {-1,1} Reverse[vec]

left
right
collinear
(* on which side of vector to-from is point_ ? *)
whichside[from_List, to_List, point_List] :=
	Sign[ perp[to-from] . (point-from) ] /.
		{ -1 -> right, 0 -> collinear, 1 -> left }
		
(* yields frac_ s.t. midpoint[{fromA, toA}, frac] is on line fromB -> toB *)
crossfrac[{fromA_List, toA_List}, {fromB_List, toB_List}] :=
	rectify[fromA-fromB, toA-fromA, perp[toB-fromB]]
crosspoint[{fromA_List, toA_List}, {fromB_List, toB_List}] :=
	midpoint[{fromA,toA},crossfrac[{fromA,toA},{fromB,toB}]]

(* solves (vec + s adj) . perp == 0, for real s *)
rectify[vec_List, adj_List, perp_List] := -((vec . perp)/(adj . perp))

(* fractions for v_ to a_, v_ to b_, v_ to c_ to make a nice triangle *)
trifracs[v_List, {a_List, b_List, c_List}] :=
	Select[ NestList[RotateLeft, {a,b,c}, 2],
	    (whichside[v, #[[1]], #[[2]]] === whichside[v, #[[1]], #[[3]]])& 
	] // If[ Length[#] == 0, {.44, .44, .44},
		With[ { sides = First /@ #, 
			mid = First[ Complement[{a,b,c}, First /@ # ] ] },
		With[ { frac = crossfrac[{v, mid}, sides] },
		    RotateRight[{.44, Min[.44, .2/frac], Min[.44, .2/frac]},
				First[First[Position[{a,b,c},mid]]]-1 ]
		] ]
	]& 
trifracs[g_Graph, v_Integer] :=
	trifracs[Last[g[[v]]],#]& @ (Last /@ List@@g[[ g[[v,1]] ]])

tripoints[g_Graph, p_Integer] :=
	MapThread[ midpoint[g,{p,#1},#2]&, { g[[p,1]], trifracs[g,p] } ]

(* returns good choices for the rotated edge *)
flippoints[g_Graph, {x_Integer, y_Integer}] :=
    With[ { mid = midpoint[g,{x,y}] }, 
	With[ { face = LeftFace[g, #] },
	    If[ facearea[g, face] < 0, leftpoint[g, #],
	    	With[ { cent = centroid[g, face],
			ab = Last[ g[[#]] ]& /@ { LeftVertex[g, #], 
				RightVertex[g, Reverse[#]] } },
		    midpoint[ {mid, cent}, 
		    	Min[.6, If[#<0,.1,.7 #]& @ crossfrac[{mid, cent}, ab] ]
		    ]
		]
	    ]
	]& /@ { {y,x}, {x,y} }
    ]


(*************** simple graph operations ***************)

(* rotates until the lowest number is first *)
canonical[ rep_List ] := 
	RotateLeft[ rep, Position[rep, Min[rep]] [[1,1]] -1 ]
(* rotates until given number is first *)
canonical[ rep_List, v_ ] := 
	RotateLeft[ rep, Position[rep, v] [[1,1]] -1 ]
(* rotates to make given number first in 2nd list, returns 1st list *)
canonical[ target_List, guide_List, v_ ] := 
	RotateLeft[target, Position[guide, v] [[1,1]] -1 ]


(* here's the (slower) canonical to use once canonical form for
	faces with multiple instances becomes important *)
canonicalU[ rep_List ] := 
	RotateLeft[ rep, Position[#, First[Sort[#]]] [[1,1]] -1 ]& @
		NestList[RotateLeft, rep, Length[rep]]


(* The number of vertices within graph-theoretic distance d from v *)
NeighborCount[g_Graph, v_, d_Integer] := Length[NeighborBall[g,v,d]]
NeighborBall[g_Graph, v_, d_Integer] := Last[NeighborList[g,v,d]]
NeighborCountList[g_Graph, v_, d_Integer] := Length /@ NeighborList[g,v,d]

(* 
	Returns a list of lists vertices within graph-theoretic 
		distance i from v, for i=0...d
*)
NeighborList[g_Graph, v_Integer, d_Integer] := NeighborList[g,{v},d]
NeighborList[g_Graph, vs_List, d_Integer] := 
	NestList[ Union[ #, Union @@ (First /@ List@@g[[#]]) ]&, vs, d]
	
Neighbors[g_Graph, v_Integer] := g[[v,1]]

(************* parts of the graph **************)

GraphVertices[g_Graph] := 
	DeleteCases[ Range[1,Length[g]] Cases[g, 
		vinfo_ :> If[vinfo===Null,0,1] ], 0]

(* returns a list of unique {v1, v2} pairs, with v1<v2 *)
GraphEdges[g_Graph] := Union @@
	DeleteCases[ MapIndexed[ Function[ {vinfo, v}, 
		If[vinfo===Null, Null, 
			Sort[Append[v, #]]& /@ First[vinfo] ] ], g ],
	Null]
	

(* faces give vertices in clockwise order, minimal number first *)
(* returns a list of unique {v1,...vn} (v1 <= vi), including external face *)
GraphFaces[g_Graph] := Union @@ (
	{ LeftFace[g, #], LeftFace[g, Reverse[#]] }& /@ GraphEdges[g] )

(* returns the face incident on the left to the edge, in canonical form *)
LeftFace[g_Graph, {from_Integer, to_Integer}] := 
    canonical[Reverse[Drop[
    	FixedPoint[ Append[#, LeftVertex[g, Take[#,-2]]]&,
		{from, to}, SameTest -> (First[#2]===Last[#2]&) ],1]]]
RightFace[g_Graph, {from_Integer, to_Integer}] := LeftFace[g, {to, from}]
	
(* coming into to_ from from_, result_ is the leftmost choice from there *)
LeftVertex[g_Graph, {from_Integer, to_Integer}] :=
	First[ RotateLeft[ canonical[ First[ g[[to]] ], from ] ] ]
RightVertex[g_Graph, {from_Integer, to_Integer}] :=
	First[ RotateRight[ canonical[ First[ g[[to]] ], from ] ] ]

(* the directed edge... *)
LeftEdge[g_Graph, {from_Integer, to_Integer}] := 
	{ to, LeftVertex[g,{from,to}] } 
RightEdge[g_Graph, {from_Integer, to_Integer}] :=
	{ to, RightVertex[g,{from,to}] } 

(************* graphics and display **************)

GraphEdgeGraphics[g_Graph, options___] := 
	Graphics[ Join[ {AbsoluteThickness[0.25]}, 
		Line[ Last /@ List@@g[[#]] ]& /@ GraphEdges[g] ], 
		PlotRange->All, AspectRatio -> 1, options]
	
GraphVertexGraphics[g_Graph, options___] :=
	Graphics[ Join[ {PointSize[.01]}, 
		List @@ (Point /@ Last /@ DeleteCases[g, Null]) ], 
		PlotRange->All, AspectRatio -> 1, options]

GraphVertexLabels[g_Graph, options___] := 
	Graphics[ Join[ {PointSize[.05]},
		List@@MapIndexed[ Function[{vinfo,v},
			If[ vinfo===Null, {}, 
				{ GrayLevel[1], Point[Last[vinfo]], 
		  		GrayLevel[0], Text[First[v],Last[vinfo]] } ]
		], g]  ], PlotRange->All, AspectRatio -> 1, options ]
	
GraphFaceGraphics[g_Graph, options___] :=
	Graphics[ Join[ {GrayLevel[.34]},
		Polygon[ Last /@ List@@g[[#]] ]& /@ GraphFaces[g] ],
		PlotRange->All, AspectRatio -> 1, options]
	
ShowLabeledGraph[g_Graph, options___] :=
	(Show[{GraphEdgeGraphics[g],GraphVertexLabels[g]}, options]; g)

ShowVertexGraph[g_Graph, options___] :=
	(Show[{GraphEdgeGraphics[g],GraphVertexGraphics[g]}, options]; g)

ShowGraph[g_Graph, options___] :=
	(Show[GraphEdgeGraphics[g], options]; g)

ShowGraphArray[{gs__Graph}, options___] := 
	(Show[GraphicsArray[GraphEdgeGraphics[#]& /@ {gs}], options]; {gs})
		
(************* operations on graph **************)

(* returns { list-of-available-positions, new-graph-structure } *)
NewVertices[g_Graph, n_Integer] :=
	If[ Length[#]>=n,
		{ First /@ Take[#,3], g },
		{ Join[ First /@ #, Length[g] + Range[1, n-Length[#]]],
			Join[ g, Graph @@ Table[ Null, {n-Length[#]} ] ] }
	]& @ Position[g, Null] 

(* 
	vertices adjacent to face or edge in clockwise order, 
	those nearest x_ first, then nearest y_, etc
*)

Circling[g_Graph, {x_Integer} ] := First[g[[x]]]
Circling[g_Graph, {x_Integer, y_Integer} ] :=
	DeleteCases[ Join[
		canonical[First[ g[[x]] ], y], canonical[First[ g[[y]] ], x]
	], x|y ] //. {a___, r_, b___, r_, c___} -> {a,r,b,c} 

Circling[g_Graph, {x_Integer, y_Integer, z_Integer} ] :=
	DeleteCases[ Join[  
		canonical[First[ g[[x]] ], z], 
		canonical[First[ g[[y]] ], x],
		canonical[First[ g[[z]] ], y]
	], x|y|z ] //. {a___, r_, b___, r_, c___} -> {a,r,b,c} 

	
ExpandVertex3[g_Graph, p_Integer] := ExpandVertex3[g,g,p]
ExpandVertex3[og_Graph, g_Graph, p_Integer] := 
    Replace[ { NewVertices[ReplacePart[g, Null, p], 3], 
    		g[[p,1]], tripoints[og, p] }, 
	    { { {p1_, p2_, p3_}, ng_ }, {a_,b_,c_}, {xy1_, xy2_, xy3_} } :>
		ReplaceParts[ ng, 
		    { p1, { {a, p2, p3}, xy1 } },
		    { p2, { {p1, b, p3}, xy2 } },
		    { p3, { {p1, p2, c}, xy3 } },
		    { a, g[[a]] /. p->p1 },
		    { b, g[[b]] /. p->p2 },
		    { c, g[[c]] /. p->p3 }
		]
    ]
    
FlipEdge3[g_Graph, {x_Integer, y_Integer}] := FlipEdge3[g,g,{x,y}]
FlipEdge3[og_Graph, g_Graph, {x_Integer, y_Integer}] :=
    Replace[ Circling[g, {x,y}], { {a_, b_, c_, d_} :> 
	Replace[ flippoints[og,{x,y}], { nx_, ny_ } :> 
	    ReplaceParts[g,
		{ x, { {a,y,d}, nx } }, 
		{ y, { {b,c,x}, ny } },
		{ b, g[[b]] /. x->y },
		{ d, g[[d]] /. y->x }
	    ]
	],
	_ -> g } (* not valid *)
    ]
    
ContractFace3[g_Graph, face_List] := ContractFace3[g, g, face]
ContractFace3[og_Graph, g_Graph, face_List] :=
    Replace[ { face, Circling[g, face] }, {
    	{ {x_, y_, z_}, {a_, b_, c_} } :> ReplaceParts[g,
		{x, { {a,b,c}, centroid[og,{x,y,z}] } },
		{y, Null},
		{z, Null},
		{ b, g[[b]] /. y->x },
		{ c, g[[c]] /. z->x }
	],
	_ -> g } (* not valid *)
    ]
(* beware of contracting an outer triangle i.e. the external face *)


RandomGraph3[g_Graph, flipsperexpand_:2] :=
    Nest[ FlipEdge3[#, Random[GraphEdges[#]]]&,  
	ExpandVertex3[g, Random[GraphVertices[g]]],
	flipsperexpand ]

RandomGraph3[n_Integer, flipsperexpand_:2] :=
    Nest[ RandomGraph3[#, flipsperexpand]&, Tetrahedron, n] 
		

(************* local topological rules *************)


(* edgecoloring is like Graph, but with vinfo containing only
	corrosponding colors for the edge *)
(* edge_ is Graph[{c1, c2, c3},...] with colors corrosponding to edges in g *)
(* edgefunc[g, {x,y}] yields a color *)
ApplyEdgeColoring[g_Graph, edgefunc_] :=
	MapIndexed[ Function[ {vinfo, v},
		If[ vinfo===Null, Null, 
			edgefunc[g, Sort[Append[v,#]]]& /@ First[vinfo]
		] ], g]

(* vertexcoloring is like Graph, but with vinfo containing only
	the color for the vertex *)
(* vertexfunc[g, x] yields a color *)
ApplyVertexColoring[g_Graph, vertexfunc_] :=
	MapIndexed[  Function[ {vinfo, v},
		If[vinfo===Null, Null, 
			vertexfunc[g, First[v]]] ], g]

(* facecoloring is like Graph, but with vinfo containing only
	corrosponding colors for the face to the left of each edge *)
(* facefunc[g, vlist] yields a color *)
ApplyFaceColoring[g_Graph, facefunc_] :=
	MapIndexed[ Function[ {vinfo, v},
		If[ vinfo===Null, Null, 
			facefunc[g, LeftFace[g,Append[v,#]] ]& /@ First[vinfo]
		] ], g]

  
(* func[g, v] is a boolean function *)
ApplyVertexSelector[g_Graph, func_] := Select[GraphVertices[g], func[g,#]&]
  
(* func[g, {x,y}] is a boolean function *)
ApplyEdgeSelector[g_Graph, func_] := Select[GraphEdges[g], func[g,#]&]

(* func[g, face] is a boolean function *)
ApplyFaceSelector[g_Graph, func_] := Select[GraphFaces[g], func[g,#]&]
  

(* ApplyOperations[] will return Graph[] if this would be exceeded *)
$MaxGraphSize = 150

(* apply the appropriate operation to each given element, provided that
	there is no contradiction with another pending operation *)
ApplyOperations3[g_Graph, vertices_List, edges_List, faces_List] :=
    With[ { vs = vertices, 
	    es = Select[Union[edges], NeighborCountList[g,#,1]==={2,6} &],
	    fs = Select[Union[faces], NeighborCountList[g,#,1]==={3,6} &] },
    With[ { 
	xvs = Complement[vs, Union[ Union @@ es, Union @@ fs ]],
	fes = Select[es, Intersection[#, vs] === {} &&
		Max[ With[{nb=NeighborBall[g,#,1]}, 
		    Length[Intersection[nb,#]]& /@ 
		    	Union[Complement[es,{#}],fs] ] ] < 2 &], 
	cfs = Select[fs, Intersection[#, vs] === {} &&
		Max[ With[{nb=NeighborBall[g,#,1]}, 
		    Length[Intersection[nb,#]]& /@ 
		    	Union[Complement[fs,{#}],es] ] ] < 2 &] },
			 
    	If[ Length[g] + 2 Length[xvs] - 2 Length[cfs] > $MaxGraphSize,
	    Graph[],
	    Fold[ExpandVertex3[g,##]&,
	    	Fold[FlipEdge3[g,##]&,
		   Fold[ContractFace3[g,##]&, g, cfs], fes], xvs]
	]
    ] ] 
 
  
(********************* high-level rule interface ******************)

  (* for expand-contract edges rules *)
  
  
(* { 1 - {4 - {4...10} } } *)
TotalisticVertexRuleCount[d_] := {1,1,7}[[1+d]]
TotalisticVertexRulePositions[0] = { {1} -> 1 }
TotalisticVertexRulePositions[1] = { {1, 4} -> 1 }
TotalisticVertexRulePositions[2] = { {1, 4, n_} :> n-3 }

(* { 2 - {4 - {4...6} , 5 - {6...10}, 6 - {6...14} } } *)
TotalisticEdgeRuleCount[d_] := {1,3,17}[[1+d]]
TotalisticEdgeRulePositions[0] = { {2} -> 1 }
TotalisticEdgeRulePositions[1] = { {2, n_} -> n-3 }
TotalisticEdgeRulePositions[2] = 
	{ {2, 4, n_} :> n-3, {2, 5, n_} :> n-3, {2, 6, n_} :> n+3 }
 
(* { 3 - {4 - {4} , 5 - {6...8}, 6 - {6...12} } } *)
TotalisticFaceRuleCount[d_] := {1,3,11}[[1+d]]
TotalisticFaceRulePositions[0] = { {3} -> 1 }
TotalisticFaceRulePositions[1] = { {3, n_} :> n-3 }
TotalisticFaceRulePositions[2] = 
	{ {3, 4, n_} :> n-3, {3, 5, n_} :> n-4, {3, 6, n_} :> n-1 }

  
(* 
   A rule may be defined governing when such operations occur:
   	It must assign to each edge a 0 ("contracting") or a 1 ("expanding").
   Then:
	A vertex expands if each incident edge is "expanding".
	A face contracts if its 3 edges are "shrinking".
	An edge flips if it is "shrinking" and 
		the 4 adjacent edges are "expanding".
*)


contract
expand
(* select operations using new mechanism, and thus correct in all cases *)

EdgeRuleSelect[e_Graph][g_Graph, v_Integer] := 
	e[[v]] === {expand, expand, expand}
EdgeRuleSelect[e_Graph][g_Graph, {x_Integer, y_Integer}] :=
	canonical[e[[x]], First[ g[[x]] ],y] === 
		{contract, expand, expand} === 
			canonical[e[[y]], First[ g[[y]] ], x]
EdgeRuleSelect[e_Graph][g_Graph, {x_Integer, y_Integer, z_Integer}] :=
	((First[canonical[e[[#1]], First[ g[[#1]] ], #2]]& @@ #)& /@
		NestList[RotateLeft, {x,y,z}, 2]) ===
			{contract, contract, contract}
EdgeRuleSelect[e_Graph][g_Graph, _] := False



TotalisticEdgeRule[r_List, d_Integer][g_Graph, {a_Integer,b_Integer}] :=
    r[[ NeighborCountList[g, {a,b}, d] /. 
    		TotalisticEdgeRulePositions[d] ]] /. {0->contract, 1->expand}
		
TotalisticEdgeRule[n_Integer, d_Integer /; d<=3] := 
	TotalisticEdgeRule[Reverse[
		IntegerDigits[n,2,TotalisticEdgeRuleCount[d]]],d]
TotalisticEdgeRule[n_Integer] := TotalisticEdgeRule[n, 2] 

(* s.t Cube, Prism, Tetra don't trivialize too quickly *)
(* 2048 posibilities *)
TERuleClassA[n_Integer] := Fold[2 #1 + #2 &, 0, Reverse @
	Fold[ Insert[#1, Last[#2], First[#2]]&, IntegerDigits[n,2,11], 
	{ {1,1}, {4,1}, {8,1}, {9,1}, {11,1}, {13,0} } ] ]
	
(* s.t Cube, Prism, Tetra are as interesting as possible *)
(* 512 posibilities *)  
TERuleClassB[n_Integer] := Fold[2 #1 + #2 &, 0, Reverse @
	Fold[ Insert[#1, Last[#2], First[#2]]&, IntegerDigits[n,2,9], 
	{ {1,1}, {4,1}, {8,1}, {9,1}, {11,1}, {13,0}, {14,0}, {15,1} } ] ]
	
(* s.t fewer rules, should be as interesting as above *)
(* 256 posibilities *)  
TERuleClassC[n_Integer] := Fold[2 #1 + #2 &, 0, Reverse @
	Fold[ Insert[#1, Last[#2], First[#2]]&, IntegerDigits[n,2,8], 
	{ {1,1}, {2,1}, {3,0}, {4,1}, {8,1}, {9,1}, {10,0}, {11,1}, {13,0} } ] ]

(* s.t fewer rules, should be as interesting as above *)
(* 128 posibilities *)  
TERuleClassD[n_Integer] := Fold[2 #1 + #2 &, 0, Reverse @
	Fold[ Insert[#1, Last[#2], First[#2]]&, IntegerDigits[n,2,7], 
	{ {1,1}, {2,1}, {3,0}, {4,1}, {5,0}, {6,0}, 
		{8,1}, {9,1}, {11,1}, {13,0} } ] ]


(* more general type of rule, using new mechanism directly *)
(* each vertex, edge, and face decides to do something depending upon
	the number of vertices within 0, 1, 2, ... d of the object *)

TotalisticRule[{vl_List, vd_Integer}, _, _][g_Graph, v_Integer] :=
    vl[[ NeighborCountList[g, v, vd] /. 
    	TotalisticVertexRulePositions[vd] ]] /. {0->False, 1->True}

TotalisticRule[_, {el_List, ed_Integer}, _][g_Graph, {x_Integer, y_Integer}] :=
    el[[ NeighborCountList[g, {x,y}, ed] /. 
    	TotalisticEdgeRulePositions[ed] ]] /. {0->False, 1->True}

TotalisticRule[_, _, {fl_List, fd_Integer}][g_Graph, 
		{x_Integer, y_Integer, z_Integer}] :=
    fl[[ NeighborCountList[g, {x,y,z}, fd] /. 
    	TotalisticFaceRulePositions[fd] ]] /. {0->False, 1->True}

TotalisticRule[_, _, _][g_Graph, _] := False

TotalisticRule[ {vn_Integer, vd_Integer /; vd<3}, 
		{en_Integer, ed_Integer /; ed<3},	
		{fn_Integer, fd_Integer /; fd<3} ] := 
    TotalisticRule[ 
	{ Reverse[IntegerDigits[vn,2,TotalisticVertexRuleCount[vd]]], vd },
	{ Reverse[IntegerDigits[en,2,TotalisticEdgeRuleCount[ed]]], ed }, 
	{ Reverse[IntegerDigits[fn,2,TotalisticFaceRuleCount[fd]]], fd }  ]

TotalisticRule[ v_Integer, e_, f_] := TotalistRule[{v,2},e,f]
TotalisticRule[ v_, e_Integer, f_] := TotalistRule[v,{e,2},f]
TotalisticRule[ v_, e_, f_Integer] := TotalistRule[v,e,{f,2}]

RandomTotalisticRule[{pv_Real,vd_Integer}, 
			{pe_Real,ed_Integer}, {pf_Real,fd_Integer}] :=
    TotalisticRule[ 
	{ Table[If[Random[]<pv,1,0], {TotalisticVertexRuleCount[vd]}], vd },
        { Table[If[Random[]<pe,1,0], {TotalisticEdgeRuleCount[ed]}], ed }, 
        { Table[If[Random[]<pf,1,0], {TotalisticFaceRuleCount[fd]}], fd } ]
RandomTotalisticRule[ f___, n_Real, l___] := 
	RandomTotalisticRule[f,{n,2},l]

  (* for face-parity vertex-edge-face selection rules *)
(* each edge, vertex, face decides to act depending upon the 
	parity of the number of sides of incident faces *)
   
 (* deterministic, numbered rules *) 
FaceParityRule[vl_List, _, _][g_Graph, v_Integer] :=
    vl[[ 1+ Count[ EvenQ[Length[LeftFace[g,{v,#}]]]& /@
    	First[g[[v]]], True ] ]] /. {0->False, 1->True}

FaceParityRule[_, el_List, _][g_Graph, {x_Integer, y_Integer}] :=
    el[[1 +
  	3 Count[ EvenQ[Length[LeftFace[g,#]]]& /@ { {x,y}, {y,x} }, True ] +
    	  Count[ EvenQ[Length[RightFace[g, LeftEdge[g,#] ]]]& /@ 
    		{ {x,y}, {y,x} }, True ] ]] /. {0->False, 1->True}

FaceParityRule[_, _, fl_List][g_Graph, {x_Integer, y_Integer, z_Integer}] :=
    fl[[ 1+ Count[ EvenQ[Length[LeftFace[g,#]]]& /@ 
  	Partition[{x,y,z,x},2,1], True ] ]] /. {0->False, 1->True}

FaceParityRule[_, _, _][g_Graph, _] := False

FaceParityRule[vn_Integer,en_Integer,fn_Integer] := 
	FaceParityRule[ Reverse[IntegerDigits[vn,2,4]], 
		Reverse[IntegerDigits[en,2,9]], 
		Reverse[IntegerDigits[fn,2,4]] ]
FaceParityRule[n_Integer] := 
	FaceParityRule[Take[#,{1,4}], Take[#, {5, -5}], Take[#, {-4,-1}]]& @
		Reverse[IntegerDigits[n,2,17]]

RandomFaceParityRule[pv_Real, pe_Real, pf_Real] :=
	FaceParityRule[ Table[If[Random[]<pv,1,0], {4}],
		Table[If[Random[]<pe,1,0], {9}],
		Table[If[Random[]<pf,1,0], {4}] ]

 (* random rules *)

RandomRule[vn_Real, _, _][g_Graph, v_Integer] := 
	Random[] < vn
RandomRule[_, en_Real, _][g_Graph, {x_Integer, y_Integer}] :=
	Random[] < en
RandomRule[_, _, fn_Real][g_Graph, {x_Integer, y_Integer, z_Integer}] :=
	Random[] < fn

(* evolution *)

GraphEvolveStep[r_TotalisticEdgeRule, g_Graph] :=
	With[ {e =  ApplyEdgeColoring[g, r]}, 
	    ApplyOperations3[g, 
	    	ApplyVertexSelector[g,EdgeRuleSelect[e]],
		ApplyEdgeSelector[g,EdgeRuleSelect[e]], 
		ApplyFaceSelector[g,EdgeRuleSelect[e]] ] ]

GraphEvolveStep[r:(_FaceParityRule|_TotalisticRule|_RandomRule), g_Graph] := 
	ApplyOperations3[g, ApplyVertexSelector[g,r],
		ApplyEdgeSelector[g,r], ApplyFaceSelector[g,r] ]

GraphEvolve[rule_, g_Graph, steps_Integer] :=
	Nest[GraphEvolveStep[rule, #]&, g, steps]

GraphEvolveList[rule_, g_Graph, steps_Integer] :=
	NestList[GraphEvolveStep[rule, #]&, g, steps]



(********** general-purpose graph stuff **********)

RandomEmbedding[g_Graph] := 
	{ First[#], {Random[Real],Random[Real]} }& /@ g

RelaxCentroidal[g_Graph, circle_:None] :=
    Module[ { faces, areas, external, todo, ng },
	faces = GraphFaces[g];
	areas = facearea[g, #]& /@ faces; 
	external = First[ Last[ Position[ areas, x_?Negative ] ] ];
	todo = GraphVertices[g]; ng=g; 
	Which[	circle===Old,
		With[ { cent = centroid[ng, faces[[external]]], 
			rad = N[Sqrt[-areas[[external]]/Pi]] },
		    ng=ReplaceParts[ng, ##]& @@ (
		    	{ {#,2}, cent + ((# (rad^2 / (# . #))&) @ 
		    		(Last[ng[[#]]]-cent)) }& /@ 
		    	faces[[external]] );
		    todo = Complement[todo, faces[[external]]]
		],
		circle===New,
		external = Length[faces];
		With[ { cent = centroid[ng, faces[[external]]], 
			rad = N[Sqrt[Abs[ areas[[external]] ]/Pi]],
			ext = If[areas[[external]]>0,Reverse[#],#]& @ 
					faces[[external]] },
		    ng=ReplaceParts[ng, ##]& @@ MapIndexed[
			{ {#1,2}, cent + 10 rad N[{Cos[#], Sin[#]}]& [
				2 Pi First[#2] / Length[ext] ] }&, ext ];
		    todo = Complement[todo, faces[[external]]]
		],
		circle===Incomplete,
		With[ { pieces = RandomPieces[
			Nest[RotateLeft, #, Random[Integer,{1,Length[#]}]]& @
			faces[[external]], 3 ] }, 
		],
		True, todo = Complement[todo, faces[[external]]]
	];
	While[ Length[todo] > 0, With[ { v = First[todo] }, 
	    ng=ReplacePart[ng, 
	 	centroid[ng, Join @@ (canonical[#,v]& /@ 
			Select[faces, MemberQ[#, v]&] )],
		{v, 2}
	    ];
	    todo=Rest[todo]
	] ];
	ng
    ]


(* example

Show[ GraphicsArray[#] ]& [ 
    Prepend[ GraphEdgeGraphics /@ Last[#], Graphics[ {
	Text["TE Rule (B): " <> ToString[First[#]],{0,0}]} ] ]& /@ 
    Fold[ With[
	{new=GraphEvolveList[TotalisticEdgeRule[TERuleClassB[#2]],Asymm12B,5]},
	If[MemberQ[Last/@#1, new], #1, Append[#1, {#2,new}]] ]&, 
    {}, Range[0,511] ]
]

Show[ GraphicsArray[#] ]& [ 
    Prepend[ GraphEdgeGraphics /@ Last[#], Graphics[ {
	Text[FontForm["TE Rule: ", {"Helvetica", 6}], {0,0}],
	Text[FontForm[First[#], {"Helvetica", 6}], {0, -1 }] },
	PlotRange->{ -2, 1} ] ]& /@ 
    Fold[ With[
	{new=GraphEvolveList[TotalisticEdgeRule[#2],Asymm12B,5]},
	If[MemberQ[Last/@#1, new], #1, Append[#1, {#2,new}]] ]&, 
    {}, Table[Random[Integer,{0,2^17-1}],{16}] ]
]


First /@ Fold[ With[ {new=GraphEvolveList[TotalisticEdgeRule[#2],Asymm12A,5]},
    If[MemberQ[Last/@#1, new], #1, Append[#1, {#2,ShowGraphArray[new,
	PlotLabel -> "Totalistic Edge Rule: " <> ToString[#2] ]}]] ]&, 
    {}, Table[Random[Integer,{0,2^17-1}],{160}] ]


Show[ GraphicsArray[#] ]& [ 
    Prepend[ GraphEdgeGraphics /@ Last[#], Graphics[ {
	Text[FontForm["TE Rule (B): ", {"Helvetica", 6}], {0,0}],
	Text[FontForm[First[#], {"Helvetica", 6}], {0, -1 }] },
	PlotRange->{ -2, 1} ] ]& /@ 
    Fold[ With[
	{new=GraphEvolveList[TotalisticEdgeRule[TERuleClassB[#2]],Prism,5]},
	If[MemberQ[Last/@#1, new], #1, Append[#1, {#2,new}]] ]&, 
    {}, Range[0,511] ]
]


Show[ GraphicsArray[#] ]& [ 
    Prepend[ GraphEdgeGraphics /@ Last[#], Graphics[ {
	Text[FontForm["TE Rule (A): ", {"Helvetica", 6}], {0,0}],
	Text[FontForm[First[#], {"Helvetica", 6}], {0, -1 }] },
	PlotRange->{ -2, 1} ] ]& /@ 
    Fold[ With[
	{new=GraphEvolveList[TotalisticEdgeRule[TERuleClassA[#2]],Prism,5]},
	If[MemberQ[Last/@#1, new], #1, Append[#1, {#2,new}]] ]&, 
    {}, Range[0,2047] ]
]


$MaxGraphSize = 1000
ShowGraphArray[GraphEvolveList[TotalisticEdgeRule[7,1], Tetrahedron, 8]]

*)