functor GHSurfaceSimplify
    (structure Complex : SIMPLICIAL_COMPLEX
     structure Metric : ERROR_METRIC
     sharing type Metric.Geometry.Point.point = Complex.Vertex.point
	     )
 : SURFACE_SIMPLIFY where Complex = Complex

=

struct

    structure Complex = Complex

    (* Provides forEachSimplex traversal primitive *)
    structure Traverser = TraverserFromComplex(structure Complex = Complex)

    type 'a complex = 'a Complex.complex

    (* raised for conditions which should never occur
       (failed patterns matches or invariants *)
    exception Strange

	(* local alias *)
    structure Vertex = Complex.Vertex
    structure Number = Metric.Geometry.Number
    structure Point = Metric.Geometry.Point

    structure VertexKey : TOTAL_ORD =
	VertexKey(structure Vertex = Complex.Vertex)

    structure VertexTable : TABLE =
	RedBlackTable(structure Index = VertexKey)

    structure Contraction =
	struct
	    type contraction = (Vertex.vertex * Vertex.vertex * Point.point * Number.t)
	    fun getError(_, _, _, e) = e

	    structure ContractionKey : TOTAL_ORD = 
		struct
		    type t = contraction
		    fun compare((_, _, _, e1), (_, _, _, e2)) = Number.compare(e1, e2)
		end

		(* convenience alias *)
	    structure Vec = Point.Vec

	    fun toString((v1, v2, v', cost) : contraction) =
		    "#" ^ Vertex.toString v1 ^ " #" ^ Vertex.toString v2 ^ " "
		    ^ Point.toString v' ^ " " ^ Number.toString cost
		
	    fun pointMid(a, b) =
		let
		    val (a', b') = (Point.toVec(a), Point.toVec(b))
		in
		    Point.fromVec (Vec.scale(Vec.+(a', b'), Number.fromReal(0.5)))
		end

	    fun optimize(#[a, b], quadrics : Metric.metric VertexTable.table) =
		let
		    fun getq x = case VertexTable.find(x, quadrics) of
			NONE => raise Strange
		      | SOME(q) => q
		    val Qa = getq a
		    val Qb = getq b
		    val Q = Metric.plus(Qa, Qb)

		    fun fallback(i,j) = 
			(* Should also allow placement on a line;
			   current code corresponds to Garlands 
			   MX_PLACE_ENDORMID *)
			let
			    val ipos = Vertex.location i
			    val jpos = Vertex.location j
			    val mid = pointMid(ipos, jpos)
			    val ei = Metric.evaluate(Q, ipos)
			    val ej = Metric.evaluate(Q, jpos)
			    val emid = Metric.evaluate(Q, mid)
			    open Number
			    val (best_end, best_e) = if ei < ej then (ipos, ei)
				else (jpos, ej)
			    val (best, error) = if emid < best_e then (mid, emid)
				else (best_end, best_e)
			in
			    (best, error)
			end
		    
		    val (new_pos, error) = 
			case Metric.optimal(Q) of
			    NONE => fallback(a, b)
			  | SOME(v) => (v, Metric.evaluate(Q, v))

		    val pair = (a, b, new_pos, error)
		in
		    pair
		end
	end


    structure ContractionQueue : PRIORITY_QUEUE =
	LeftistHeap (structure Priority = Contraction.ContractionKey)

    type 'a contractibleComplex =
	'a complex *  Metric.metric VertexTable.table * 
	ContractionQueue.pqueue * 
	(('a option -> 'a option) * ('a option -> bool))

    type number = Metric.Geometry.Number.t

    fun makeTable(C : 'a complex, travFuns) =
	let
	    val vertices = Complex.vertices C
	    val tbl = Sequence.foldl
		(fn (vertex, t) => VertexTable.insert((vertex, Metric.zero), t))
		VertexTable.empty
		vertices
	    fun initFace(s, tbl) =
		let
		    val vertices = Sequence.map 
			Complex.Vertex.location
			(Complex.Simplex.vertices s)
		    (* Calculate metric for this face *)
		    val q = Metric.init(vertices)

		    fun addQ qo = 
			(case qo of 
			     NONE => raise Strange
			   | SOME(q') => SOME(Metric.plus(q, q')))

		    (* Distribute metric to vertices *)
		    val #[v1, v2, v3] = Complex.Simplex.vertices s
		    val tbl = VertexTable.update addQ (v1, tbl)
		    val tbl = VertexTable.update addQ (v2, tbl)
		    val tbl = VertexTable.update addQ (v3, tbl)
		in
		    tbl
		end
	in
	    Traverser.forEachSimplex (2, initFace, travFuns) (C, tbl) 
	end


    fun makeQueue(c : 'a complex, table, travFuns) =
	let
	    fun insertPair(s, queue) =
		let
		    val vertices = Complex.Simplex.vertices s
		in
		    ContractionQueue.insert(Contraction.optimize(vertices, table), queue)
		end
	in
	    Traverser.forEachSimplex (1, insertPair, travFuns) (c, ContractionQueue.empty)
	end

    fun initialize travFuns (C : 'a complex) : 'a contractibleComplex =
	let
	    val _ = if Complex.dim C <> 2 then raise Strange else ()
	    val table = makeTable(C, travFuns)
	    val queue = makeQueue(C, table, travFuns)
	in
	    (C, table, queue, travFuns)
	end

    structure Simplex = Complex.Simplex

    fun mergeVertices(c, v1 : Vertex.vertex, v2 : Vertex.vertex,
		       v' : Vertex.vertex) =
	let
	    val faces = Sequence.flatten #[(Complex.grep c (2, Simplex.vertex v1)), 
					   (Complex.grep c (2, Simplex.vertex v2))]
		
		
	    (* c' is complex without any involved faces *)
	    val c' = Sequence.foldl (fn (s, c) => Complex.rem(c, s)) c faces

	
	    fun veq(v1, v2) = case (Vertex.compare(v1, v2)) of
		EQUAL => true
	      | _ => false
		    
	    fun replace' v = if veq(v, v1) orelse veq(v, v2) then v' else v
	    fun replace(s) = Sequence.map replace' s
	
	    (* Works strictly with 2d simplices *)
	    (* Find degenerate triangles *)
	    fun isDup(v : Vertex.vertex Sequence.seq) =
		let
		    val #[v1, v2, v3] = v
		in
		    veq(v1,v2) orelse veq(v1,v3) orelse veq(v2, v3)
		end
	    fun partition f s =
		(* Should call Sequence.partition, but it doesn't
		   exist *)
		let
		    val kept = Sequence.filter f s
		    val dropped = Sequence.filter (fn x => not (f x)) s
		in
		    (kept, dropped)
		end


	    val faces = Sequence.map (fn s  => (replace (Simplex.vertices s))) faces
	    val (dead, changed) = partition isDup faces
	    val c' = Sequence.foldl (fn (s, c) => Complex.add(c, Simplex.fromSeq s)) c' changed
		(* might be cheaper to pull out list of new edges here, rather than grep in
		   update_pairs *)
	in
	    (c', Sequence.map (Simplex.fromSeq) changed, Sequence.map Simplex.fromSeq dead)
	end


    fun mergeAndUpdatePairs((v1, v2, v', error), complex, metrics, queue) =
	let
	    val (complex', changed, dead) = mergeVertices(complex, v1, v2, v')
	    val newEdges = Complex.grep complex' (1, Complex.Simplex.vertex v')
	    fun insertPair(s, queue) =
		ContractionQueue.insert(
			   Contraction.optimize(Complex.Simplex.vertices s, metrics), 
			   queue)
	in
	    (complex', Sequence.foldl insertPair queue newEdges, changed, dead)
	end

    fun writeContractionState(outfile, (c, metrics, _, _)) =
	let
	    fun writeMatrix v =
		let
		    val SOME(metric) = VertexTable.find(v, metrics)
		    val str = Metric.toString metric
		in
		    TextIO.output(outfile, 
				  "Matrix #" ^ (Simplex.Vertex.toString v) ^ " " ^ str ^ "\n")
		end
	in
	    Sequence.app writeMatrix (Complex.vertices c);
	    TextIO.flushOut(outfile)
	end
	    



    type contractionRecord = 
	(Complex.Simplex.Vertex.vertex list * Complex.Simplex.Vertex.vertex * number)

    fun contract' (c, metrics, queue, travFuns) =
	let

	    (* We may well have junk in the queue, so make sure all
	       vertices are actually in the complex *)
	    fun firstRealContraction(queue) =
		let
		    val (con as (v1, v2, _, _), queue') = ContractionQueue.deleteMin(queue)
		    fun isV(v) =
			Sequence.length(Complex.grep c (2, Complex.Simplex.vertex v)) > 0
		in
		    if isV(v1) andalso isV(v2) then
			(con, queue')
		    else
			firstRealContraction(queue')
		end
	    val ((v1, v2, loc, error), queue') = firstRealContraction queue
	    val v' = Vertex.new loc
	    val (metrics', SOME(M1)) = VertexTable.deleteReturn(v1, metrics)
	    val (metrics', SOME(M2)) = VertexTable.deleteReturn(v2, metrics')
	    val metrics' = VertexTable.insert((v', Metric.plus(M1, M2)), metrics')
(*	    val (complex', queue', changed, dead) = mergeAndUpdatePairs(
				(v1, v2, v', error), c, metrics', queue')
		                handle Complex.Book => contract'(c, metrics, queue', travFuns) *)
				     (* Until complex code allows books to be formed,
				        discard any contractions that would create one *)
	    fun extract(c, q, changed, dead) = ((c, metrics', q, travFuns), ([v1, v2], v', error), (changed, dead))
	in
	    extract(mergeAndUpdatePairs((v1, v2, v', error), c, metrics', queue'))
	     handle Complex.Book => contract'(c, metrics, queue', travFuns)
(*
	    ((complex', metrics', queue', travFuns), ([v1, v2], v', error), (changed, dead))
*)
	end

    fun contract(c) =
	let
	    val (result, _, _) = contract'(c)
	in
	    result
	end

    fun result (c, _, _, _) = c

    fun error (_, _, queue, _) =
	let
	    val nextPair = ContractionQueue.findMin queue
	in
	    Contraction.getError nextPair
	end

end