(******************************************************************************
 ** TestHull.sml
 ** sml
 **
 ** Guy Blelloch
 ** Code for testing and running statistics on the Convex Hull code
 ******************************************************************************)

structure TestHull =
    struct
	structure G3 = ExactGeometry3d
	structure R = G3.Number
	structure Vertex = Vertex(type point = G3.Point.t)
	structure Simplex = SimplexFromVertex(structure Vertex = Vertex)
	structure C2 = SimpCompFromSimplex(structure Simplex = Simplex )

	structure IHull = IncrementalHull(structure Geometry = G3
					  structure Complex = C2)

	structure BHull = BulldozerHull(structure Geometry = G3
					structure Complex = C2)

	(* Changing the following line to either IHull or BHull
	    switches between testing one or the other Hull *)
	structure Hull = IHull

	fun toPoint(x,y,z) =
	    G3.Point.fromSeq(#[R.fromReal(x),R.fromReal(y), R.fromReal(z)])

	fun generatePoints distfunc state 0 = nil
	  | generatePoints distfunc state npts =
	    let
		val (Maps.p3D(x,y,z),state) = distfunc(state)
		val pt = toPoint(x,y,z)
	    in
		pt::(generatePoints distfunc state (npts-1))
	    end

	val inSphere = RC4RandDist.UnitBall.nextSeed

	fun onSphere(state) =
	    let
		val (pt, state) = inSphere state
	    in
		(Maps.mapToUnitSphere(pt),state)
	    end

	fun poles(a,b) = GaryDist.poleHeavy (Maps.p3D(a,a,b))

	fun uniformProj(state) = 
	    let
		val (Maps.p3D(x,y,z),state) = inSphere state
	    in
		(Maps.p3D(x,y,x*x+y*y+z*z),state)
	    end

	fun kuzmin(state) =
	    let
		val (pt,state) = inSphere state
	    in
		(Maps.kuzmin(pt),state)
	    end
	    

    (* Print out the points and their radius *)
    fun pp nil = ()
      | pp ((p1,p2,p3)::ll) =
      let
	fun rts p = (Real.toString p) ^ "  " 
      in
	(print (rts p1); print (rts p2); print (rts p3);
	 print (rts (Math.sqrt(p1*p1 + p2*p2 + p3*p3)));
	 print "\n"; 
	 pp ll)
      end

    fun printhull(hull) =
      let
	val simps = C2.simplices hull 2;
	val _ = Sequence.map (fn s => print(Simplex.toString(s) ^ "\n")) simps
      in 
	true
      end

    fun timeOnce(pts) =
	let
	    val tt = Timer.startCPUTimer ()
	    val hh = Hull.hull pts
	    val {usr=uu,sys=ss,gc=gg} = Timer.checkCPUTimer tt
	    val _ = print("\nTime for " ^ (Int.toString (length pts)) ^
			  " points: (usr,sys,gc) = (" ^
			  (Time.toString uu) ^ "," ^
			  (Time.toString ss) ^ "," ^
			  (Time.toString gg) ^ "), total = " ^
			  (Time.toString (Time.+(uu,Time.+(ss,gg)))) ^ "\n")
	in
	    ()
	end

    fun timeHull nil = ()
      | timeHull (npt::npts) =
      let
	val ptsr = generatePoints onSphere (RC4byte.State.start "Hello World") npt
	val _ = timeOnce(ptsr)
      in
	timeHull npts
      end

    fun complexToMaple hull =
      let
	val fname = "testnew.mpl"

	fun simpToMaple(s,t) =
	  let
	    fun ptTost pt = 
		let
		    val #[p,q,r] = G3.Point.toSeq pt
		in
		    "[" ^ 
		    (R.toString p) ^ ", " ^
		    (R.toString q) ^ ", " ^
		    (R.toString r) ^ "]"
		end

	    fun rmtilde #"~" = #"-" 
	      | rmtilde a = a

	    val prn = Vector.map  
	      ((String.map rmtilde) o ptTost o Vertex.location)
	      (Simplex.vertices s)

	    fun pp i = Vector.sub(prn,i)
	  in
	    "polygon([" ^ pp(0) ^ 
	    ",\n\t " ^ pp(1) ^ ",\n\t " ^ pp(2) ^ "]),\n" ^ t
	  end

	fun select s =  (* Remove threefold representation of simplicies *)
	  let
	    val vv = Simplex.vertices s
	    val v0 = Vector.sub(vv,0)
	    val v1 = Vector.sub(vv,1)
	    val v2 = Vector.sub(vv,2)
	  in
	    (case Vertex.compare(v0,v1) of LESS => true | _ => false) andalso 
	    (case Vertex.compare(v0,v2) of LESS => true | _ => false)
	  end

	val pp = Vector.foldl simpToMaple "" 
	  (Sequence.filter select (C2.simplices hull 2))

	val _ = print ("Creating Maple file " ^ fname ^ "\n")

	val outfile = TextIO.openOut fname
	val _       = TextIO.output(outfile, "with(plottools):\n")
	val _       = TextIO.output(outfile, "ll := [\n")
	val _       = TextIO.output(outfile,  pp)
	val _       = TextIO.output(outfile, 
				    "point([0.0,0.0,0.0]) ]:\n")
	val _       = TextIO.output(outfile, 
				    "plots[display](ll,scaling=constrained);")
	val _       = TextIO.closeOut outfile
      in
	()
      end

    fun complexToMathematica hull =
      let
	val fname = "testguy.mta"

	fun simpToMaple(s,t) =
	  let
	    fun ptTost pt = 
		let
		    val #[p,q,r] = G3.Point.toSeq pt
		in
		    "{" ^ 
		    (R.toString p) ^ ", " ^
		    (R.toString q) ^ ", " ^
		    (R.toString r) ^ "}"
		end

	    fun rmtilde #"~" = #"-" 
	      | rmtilde a = a

	    val prn = Vector.map  
	      ((String.map rmtilde) o ptTost o Vertex.location)
	      (Simplex.vertices s)

	    fun pp i = Vector.sub(prn,i)
	  in
	    "Polygon[{" ^ pp(0) ^ 
	    ",\n\t " ^ pp(1) ^ ",\n\t " ^ pp(2) ^ "}],\n" ^ t
	  end

	fun select s =  (* Remove threefold representation of simplicies *)
	  let
	    val vv = Simplex.vertices s
	    val v0 = Vector.sub(vv,0)
	    val v1 = Vector.sub(vv,1)
	    val v2 = Vector.sub(vv,2)
	  in
	    (case Vertex.compare(v0,v1) of LESS => true | _ => false) andalso 
	    (case Vertex.compare(v0,v2) of LESS => true | _ => false)
	  end

	val pp = Vector.foldl simpToMaple ""
	  (Sequence.filter select (C2.simplices hull 2))

	val _ = print ("Creating Mta file " ^ fname ^ "\n")

	val outfile = TextIO.openOut fname
	val _       = TextIO.output(outfile, "Show[Graphics3D[{\n")
	val _       = TextIO.output(outfile, pp)

	val _       = TextIO.output(outfile, 
				    "Point[{0.0,0.0,0.0}] }]]")
	val _       = TextIO.closeOut outfile
      in
	()
      end

(*
    val hull = Hull.hull(pts)

    val _ = complexToMaple(hull)
    val _ = complexToMathematica(hull)
*)

    fun testHull distributionFunction n = 
	let
	    val state = (RC4byte.State.start "Hello World")
	    val pts = generatePoints distributionFunction state n
	    open Hull
	    val _ = clearCount()
	    val _ = Vcount.n := 0
	    val h = hull(pts)
	    (* val _ = complexToMathematica(h) *)
	    val ret = {dict_ops = (3*(!addCount)+ 3*(!remCount) + 
		                  ((!searchCount)+(!remCount))),
		       dict_cmps = !Vcount.n,
		       float_ops = 7* (!lineSideCount) + 20*(!makePlaneCount),
		       make_plane = (!makePlaneCount),
		       line_side = (!lineSideCount)}
	in
	    (n,(Sequence.length (C2.simplices h 2)),ret)
	end

    fun testloop distributionFunction nil = nil
      | testloop distributionFunction (n::r) = 
	(testHull distributionFunction n)::(testloop distributionFunction r)

    val fieldlen = 12
    fun padstr(str,n) = StringCvt.padLeft #" " n str
    fun padint(i,n) = padstr(Int.toString(i),n)


    fun printStatsR(nil) = ""
      | printStatsR(stat::rest) = 
	let
	    val (n,t,{dict_ops = d,
		      dict_cmps = c,
		      float_ops = f,
		      make_plane = m,
		      line_side = l}) = stat
	    val line = (padint(n,8) ^ padint(t,8) ^ 
			padint(d,fieldlen) ^ padint(c,fieldlen) ^ 
			padint(f,fieldlen) ^ padint(m,fieldlen) ^ padint(l,fieldlen) ^
			"\n")
	in
	    line ^ printStatsR(rest)
	end

    fun statsToString(l) =
        ((padstr("",8) ^ 
		padstr("Triangs",8) ^ 
		padstr("Dict Ops",fieldlen) ^ 
		padstr("Dict Comps",fieldlen) ^ 
		padstr("Float Ops",fieldlen) ^ 
		padstr("Make Plane",fieldlen) ^ 
		padstr("Line Side",fieldlen) ^ 
		"\n") ^ printStatsR(l))

    val lengths = [1000,2000,4000,8000,16000,32000,64000,128000]
    val longLengths = [1000,2000,4000,8000,16000,32000,64000,128000,256000,512000,1024000]
    val dists = [(onSphere,"Uniform on sphere"),
		 (poles(100.0,1.0),"Equator heavy on sphere (100.0,100.0,1.0)"),
		 (poles(1.0,100.0),"Poles heavy on sphere (1.0,1.0,100.0)"),
		 (inSphere,"Uniform in sphere"),
		 (uniformProj,
		  "Border heavy (in sphere translated to: (x, y, x^2 +y^2 + z^2))")]

    fun testAll(lengths,dists,filename) =
	let
	    val stream = TextIO.openOut(filename)
	    fun testrec(nil) = ()
	      | testrec((f,name)::r) = (TextIO.output(stream,"\n" ^ name ^ "\n"); 
					TextIO.output(stream,
						      statsToString(testloop f lengths));
					testrec(r))
	    val _ = testrec(dists)
	    val _ = TextIO.closeOut(stream)
	in
	    ()
	end

    (* A simple test *)
    fun test () = testHull onSphere 1000

  end (* TestHull *)