# Computation of the Hausdorff Distance between Two Compact Convex Sets
## Appendix A
### Julia Computer Code
bold = bnew  end  end  end  return anew + hnew  end  
\"\"\"Finds the farthest point on A from B by Frank-Wolfe.\"\"\" function FrankWolfe(SuppA, PB, x0) (zold, xnew) = (copy(x0), similar(x0)) for iter = 1:100  xnew = SuppA(xold - PB(xold))  if norm(mew - xold) < 1.0e-10  break  else  xold.= xnew  end end  far = norm(xnew - PB(xnew))  return (far, xnew)  end \"\"\"Finds the farthest point on A from B by projected gradient ascent.\"\"\" function farthest(PA, PB, x0) (xold, xnew) = (copy(x0),copy(x0)) for iter = 1:100  xnew = PA(2 xold - PB(xold))  if norm(mew - xold) < 1.0e-10  break  else  xold.= xnew  end  end  far = norm(xnew - PB(xnew))  return (far, xnew)  end  
\"\"\"Finds the farthest point on A from B by homotopy.\"\"\" function farthest_homotopy(PA, PB, SA, CenterA, CenterB, x0, n, method)  x = BallProjection(x0) (far, homotopy_points, conv) = (0.0, 10, 1.0e-10)  for iter = 0:homotopy_points  if iter == 0 # ball to ball  d = norm(CenterA - CenterB) (far, x) = (d, (1 + 1 / d) * CenterA - CenterB / d)  elseif iter == homotopy_points # d(A, B)  if method == \"proj grad\"  (far, x) = farthest(PA, PB, x)  elseif method == \"Frank-Wolfe\"  (far, x) = [PERSON](SA, PB, x)  end  else # intermediate case  t = iter / homotopy_points  PMB(z) = MinkowskiNear(PB, z, t, conv)  if method == \"proj grad\"  PMA(z) = MinkowskiNear(PA, z, t, conv)  (far, x) = fartherthest(PMA, PMB, x)  elseif method == \"Frank-Wolfe\"  SN(z) = SA(t * z) + BallSupp((1 - t) * z)  (far, x) = FrankWolfe(SM, PMB, x)  end end end  return(far, x)  end  \"\"\"Orchestrates Hausdorff distance estimation.\"\"\" function master(ProjA, ProjB, SuppA, SuppB, CenterA, CenterB, method, homotopy, n, trials, io)  #  (count, tries, optimum, obj) = (0, 100, 0.0, zeros(trials)) x0 = zeros(n)  PA(z) = RecenterProjection(ProjA, z, CenterA)  PB(z) = RecenterProjection(ProjB, z, CenterB)  for trial = :trials  success = false  for i = 1:tries # find a point in A \ B  x0 = PA(randn(n))  if norm(PB(x0) - x0) > 1.0e-10  success = true  break  end  end if homotopy # solve for d(A, B) (objA, xA) = farthest_homotopy(PA, PB, SuppA, CenterA, CenterB, x0, n, method) else  if method == \"proj grad\" && success # solve for d(A, B)  (objA, xA) = farthest(PA, PB, x0)  elseif method == \"Frank-Wolfe\" && success  (objA, xA) = FrankWolfe(SuppA, PB, x0)  else  objA = 0.0  end  end success = false  for i = 1:tries # find point in B \ A  x0 = PB(randn(n))  if norm(PA(x0) - x0) > 1.0e-10  success = true  break  end  end if homotopy # solve for d(B, A)  (objB, xB) = farthest_homotopy(PB, PA, SuppB, CenterB, CenterA, x0, n, method)  else  if method == \"proj grad\" && success  (objB, xB,) = farthest(PB, PA, x0)  elseif method == \"Frank-Wolfe\" && success  (objB, xB) = FrankWolfe(SuppB, PA, x0)  else  objB = 0.0  end  end  end obj[trial] = max(objA, objB) # Hausdorff distance  if obj[trial] > optimum + 10.0e-10 # update count of maximum distance  count = 1  optimum = obj[trial]  elseif obj[trial] > optimum - 10.0e-8  count = count + 1  end  end (avg, stdev) = (mean(obj), std(obj))  if stdev < 1.0e-10 stdev = 0.0 end  return (fraction = count / trials, optimum, avg, stdev)  end  outfile = \"Hausdorff.out\";  io = open(outfile, \"w\");  trials = 100;points = 100000 println(io,\"Set Pair\",\" & \",\"p\",\" & \",\"Method\", & \",\"Monotopy\",\" & \",  \"Maximum\",\" & \",\"Mean\",\" & \",\"Std\",\" & \",\"Seconds\",\" \") for n in [2, 3, 10, 100, 1000] for i = 1:2 if i == 1 CenterA = ones(n) CenterB = zeros(n) ProjA = BoxProjection ProjB = BallAndOrhathatDepjection SuppB = BoxSupp SuppB = BallAndOrhatNumberSupp title = \"dH(Box, ball and orthant)\" elseif i == 2 CenterA = zeros(n) CenterB = ones(n) ProjA = SimplexProjection ProjB = LIBallProjection SuppB = SimplexSupp SuppB = LIBallSupp title = \"dH(simplex, LI ball)\" end # (method, homotopy) = (\"proj grad\", false) Random.seed1(1234) time = Gelapsed (fraction, optimum, avg, stdev) = master(ProjA, ProjB, SuppA, SuppB, CenterA, CenterB, method, homotopy, n, trials, io) println(io,title,\" & \",n,\" & \",method,\" & \",homotopy,\" & \", round(optimum, sigdigits=5),\" & \",round(avg, sigdigits=5),\" & \",  round(stdev, sigdigits=5),\" & \",round(time/trials, sigdigits=3),\" \")
method, homotopy) = (\"proj grad\", true); Random.seed1(1234) time = Gelapsed (fraction, optimum, avg, stdev) = master(ProjA, ProjB, SuppA, SuppB, CenterA, CenterB, method, homotopy, n, trials, io) println(io,title,\" & \",n,\" & \",method,\" & \",homotopy,\" & \", round(optimum, sigdigits=5),\" & \",round(avg, sigdigits=5),\" & \", round(stdev, sigdigits=5),\" & \",round(time/trials,sigdigits=3),\" \")
method, homotopy) = (\"Frank-Wolfe\", false); Random.seed1(1234) time = Gelapsed (fraction, optimum, avg, stdev) = master(ProjA, ProjB, SuppA, SuppB, CenterA, CenterB, method, homotopy, n, trials, io) println(io,title,\" & \",n,\" & \",method,\" & \",homotopy,\" & \", round(optimum, sigdigits=5),\" & \",round(avg, sigdigits=5),\" & \", round(stdev, sigdigits=5),\" & \",round(time/trials,sigdigits=3),\" \")
method, homotopy) = (\"Frank-Wolfe\", true); Random.seed1(1234) time = Gelapsed (fraction, optimum, avg, stdev) = master(ProjA, ProjB, SuppA, SuppB, CenterA, CenterB, method, homotopy, n, trials, io) println(io,title,\" & \",n,\" & \",method,\" & \",homotopy,\" & \", round(optimum, sigdigits=5),\" & \",round(avg, sigdigits=5),\" & \", round(stdev, sigdigits=5),\" & \",round(time/trials,sigdigits=3),\" \")
method, homotopy) = (\"Point cloud\", false); Random.seed1(1234) points = 10000 A = zeros(n, points) B = zeros(n, points) if i == 1 (a, b) = (ones(n), 2 * ones(n)) for j = 1:points  AL(: j) = RandomBox(a, b) B[:, j] = RandomBallOrhatant(1.0, n)