gromov_hausdorff.lean

/-
Copyright (c) 2018 Sébastien Gouëzel. All rights reserved.
Released under Apache 2.0 license as described in the file LICENSE.
Author: Sébastien Gouëzel

The Gromov-Hausdorff distance on nonempty compact spaces.

We introduces the space of all nonempty compact spaces, up to isometry,
called `GH_space`, and endow it with a metric space structure. The distance,
known as the Gromov-Hausdorff distance, is defined as follows: given two
nonempty compact spaces X and Y, their distance is the minimum Hausdorff distance
between all possible isometric embeddings of X and Y in all metric spaces.
To define properly the Gromov-Hausdorff space, we consider the non-empty
compact subsets of ℓ^∞(ℝ) up to isometry, which is a well-defined type,
and define the distance as the infimum of the Hausdorff distance over all
embeddings in ℓ^∞(ℝ). We prove that this coincides with the previous description,
as all separable metric spaces embed isometrically into ℓ^∞(ℝ), through an
embedding called the Kuratowski embedding.
To prove that we have a distance, we should show that if spaces can be coupled
to be arbitrarily close, then they are isometric. More generally, the Gromov-Hausdorff
distance is realized, i.e., there is a coupling for which the Hausdorff distance
is exactly the Gromov-Hausdorff distance. This follows from a compactness
argument, essentially following from Arzela-Ascoli.
-/

import topology.bounded_continuous_function topology.metric_space.hausdorff_distance
       tactic.tidy

noncomputable theory
local attribute [instance, priority 0] classical.prop_decidable
universes u v w

open classical lattice set function topological_space filter metric quotient bounded_continuous_function
open sum (inl inr)
set_option class.instance_max_depth 60

@[reducible] def l_infty_ℝ : Type := bounded_continuous_function ℕ ℝ

local attribute [instance] metric_space_sum

namespace Gromov_Hausdorff

section Gromov_Hausdorff_realized
/- Preliminary 1 : construction of the optimal coupling.
This section contains preliminaries to show that the Gromov-Hausdorff distance
is realized. For this, we consider candidate distances on the disjoint union
α ⊕ β of two compact nonempty spaces, almost realizing the Gromov-Hausdorff
distance, and show that they form a compact family by applying Arzela-Ascoli
theorem. The existence of a minimizer follows. -/

section definitions
variables (α : Type u) (β : Type v)
  [metric_space α] [compact_space α] [nonempty α]
  [metric_space β] [compact_space β] [nonempty β]

@[reducible] private def prod_space_fun : Type* := ((α ⊕ β) × (α ⊕ β)) → ℝ
@[reducible] private def Cb : Type* := bounded_continuous_function ((α ⊕ β) × (α ⊕ β)) ℝ

private def max_var : ℝ := 2 * diam (univ : set α) + 1 + 2 * diam (univ : set β)

private lemma one_le_max_var : 1 ≤ max_var α β := calc
  (1 : real) = 2 * 0 + 1 + 2 * 0 : by simp
  ... ≤ 2 * diam (univ : set α) + 1 + 2 * diam (univ : set β) :
  begin
    apply_rules [add_le_add, mul_le_mul_of_nonneg_left, diam_nonneg, diam_nonneg],
    norm_num, norm_num, norm_num
  end

/-- The set of functions on α ⊕ β that are candidates distances to realize the
minimum of the Hausdorff distances between α and β in a coupling -/
private def candidates : set (prod_space_fun α β) :=
  {f | (((((∀x y : α, f (sum.inl x, sum.inl y) = dist x y)
    ∧ (∀x y : β, f (sum.inr x, sum.inr y) = dist x y))
    ∧ (∀x y, f (x, y) = f (y, x)))
    ∧ (∀x y z, f (x, z) ≤ f (x, y) + f (y, z)))
    ∧ (∀x, f (x, x) = 0))
    ∧ (∀x y, f (x, y) ≤ max_var α β) }

/-- Version of the set of candidates in bounded_continuous_functions, to apply
Arzela-Ascoli -/
private def candidates_b : set (Cb α β) := {f : Cb α β | f.val ∈ candidates α β}

end definitions --section

section constructions

variables {α : Type u} {β : Type v}
variables [metric_space α] [compact_space α] [nonempty α] [metric_space β] [compact_space β] [nonempty β]
variables {f : prod_space_fun α β} {x y z t : α ⊕ β}
local attribute [instance, priority 0] inhabited_of_nonempty'

private lemma max_var_bound : dist x y ≤ max_var α β := calc
  dist x y ≤ diam (univ : set (α ⊕ β)) :
    dist_le_diam_of_mem (bounded_of_compact compact_univ) (mem_univ _) (mem_univ _)
  ... = diam (inl '' (univ : set α) ∪ inr '' (univ : set β)) :
    by apply congr_arg; ext x y z; cases x; simp [mem_univ, mem_range_self]
  ... ≤ diam (inl '' (univ : set α)) + dist (inl (default α)) (inr (default β)) + diam (inr '' (univ : set β)) :
    diam_union (mem_image_of_mem _ (mem_univ _)) (mem_image_of_mem _ (mem_univ _))
  ... = diam (univ : set α) + (dist (default α) (default α) + 1 + dist (default β) (default β)) + diam (univ : set β) :
    by rw [diam_image_of_isometry isometry_on_inl, diam_image_of_isometry isometry_on_inr];
      [refl, by apply_instance, by apply_instance]
  ... = 1 * diam (univ : set α) + 1 + 1 * diam (univ : set β) : by simp
  ... ≤ 2 * diam (univ : set α) + 1 + 2 * diam (univ : set β) :
  begin
    apply_rules [add_le_add, mul_le_mul_of_nonneg_right, diam_nonneg, diam_nonneg, le_refl],
    norm_num, norm_num
  end

private lemma candidates_symm (fA : f ∈ candidates α β) : f (x, y) = f (y ,x) :=
by rw candidates at fA; exact fA.1.1.1.2 x y

private lemma candidates_triangle (fA : f ∈ candidates α β) : f (x, z) ≤ f (x, y) + f (y, z) :=
by rw candidates at fA; exact fA.1.1.2 x y z

private lemma candidates_refl (fA : f ∈ candidates α β) : f (x, x) = 0 :=
by rw candidates at fA; exact fA.1.2 x

private lemma candidates_nonneg (fA : f ∈ candidates α β) : 0 ≤ f (x, y) :=
begin
  have : 0 ≤ 2 * f (x, y) := calc
    0 = f (x, x) : (candidates_refl fA).symm
    ... ≤ f (x, y) + f (y, x) : candidates_triangle fA
   ... = f (x, y) + f (x, y) : by rw [candidates_symm fA]
    ... = 2 * f (x, y) : by ring,
  by linarith
end

private lemma candidates_dist_inl (fA : f ∈ candidates α β) (x y: α) : f (inl x, inl y) = dist x y :=
by rw candidates at fA; exact fA.1.1.1.1.1 x y

private lemma candidates_dist_inr (fA : f ∈ candidates α β) (x y : β) : f (inr x, inr y) = dist x y :=
by rw candidates at fA; exact fA.1.1.1.1.2 x y

private lemma candidates_le_max_var (fA : f ∈ candidates α β) : f (x, y) ≤ max_var α β :=
by rw candidates at fA; exact fA.2 x y

/-- candidates are bounded by max_var α β -/
private lemma candidates_dist_bound  (fA : f ∈ candidates α β) :
  ∀ {x y : α ⊕ β}, f (x, y) ≤ max_var α β * dist x y
| (inl x) (inl y) := calc
    f (inl x, inl y) = dist x y : candidates_dist_inl fA x y
    ... = dist (inl x) (inl y) : begin rw @sum.dist_eq α β; refl end
    ... = 1 * dist (inl x) (inl y) : by simp
    ... ≤ max_var α β * dist (inl x) (inl y) :
      mul_le_mul_of_nonneg_right (one_le_max_var α β) dist_nonneg
| (inl x) (inr y) := calc
    f (inl x, inr y) ≤ max_var α β : candidates_le_max_var fA
    ... = max_var α β * 1 : by simp
    ... ≤ max_var α β * dist (inl x) (inr y) :
      mul_le_mul_of_nonneg_left sum.one_dist_le (le_trans (zero_le_one) (one_le_max_var α β))
| (inr x) (inl y) := calc
    f (inr x, inl y) ≤ max_var α β : candidates_le_max_var fA
    ... = max_var α β * 1 : by simp
    ... ≤ max_var α β * dist (inl x) (inr y) :
      mul_le_mul_of_nonneg_left sum.one_dist_le (le_trans (zero_le_one) (one_le_max_var α β))
| (inr x) (inr y) := calc
    f (inr x, inr y) = dist x y : candidates_dist_inr fA x y
    ... = dist (inr x) (inr y) : by rw @sum.dist_eq α β; refl
    ... = 1 * dist (inr x) (inr y) : by simp
    ... ≤ max_var α β * dist (inr x) (inr y) :
      mul_le_mul_of_nonneg_right (one_le_max_var α β) dist_nonneg

/-- Technical lemma to prove that candidates are Lipschitz -/
private lemma candidates_lipschitz_aux (fA : f ∈ candidates α β) : f (x, y) - f (z, t) ≤ 2 * max_var α β * dist (x, y) (z, t) :=
calc
  f (x, y) - f(z, t) ≤ f (x, t) + f (t, y) - f (z, t) : add_le_add_right (candidates_triangle fA) _
  ... ≤ (f (x, z) + f (z, t) + f(t, y)) - f (z, t) :
    add_le_add_right (add_le_add_right (candidates_triangle fA) _ ) _
  ... = f (x, z) + f (t, y) : by simp
  ... ≤ max_var α β * dist x z + max_var α β * dist t y :
    add_le_add (candidates_dist_bound fA) (candidates_dist_bound fA)
  ... ≤ max_var α β * max (dist x z) (dist t y) + max_var α β * max (dist x z) (dist t y) :
  begin
    apply add_le_add,
    apply mul_le_mul_of_nonneg_left (le_max_left (dist x z) (dist t y)) (le_trans zero_le_one (one_le_max_var α β)),
    apply mul_le_mul_of_nonneg_left (le_max_right (dist x z) (dist t y)) (le_trans zero_le_one (one_le_max_var α β)),
  end
  ... = 2 * max_var α β * max (dist x z) (dist y t) :
    by simp [dist_comm]; ring
  ... = 2 * max_var α β * dist (x, y) (z, t) : by refl

/-- Candidates are Lipschitz -/
private lemma candidates_lipschitz (fA : f ∈ candidates α β) (p q : (α ⊕ β) × (α ⊕ β)) :
  dist (f p) (f q) ≤ 2 * max_var α β * dist p q :=
begin
  rcases p with ⟨x, y⟩,
  rcases q with ⟨z, t⟩,
  rw real.dist_eq,
  apply abs_le_of_le_of_neg_le,
  { exact candidates_lipschitz_aux fA },
  { have : -(f (x, y) - f (z, t)) = f (z, t) - f (x, y) := by ring,
    rw [this, dist_comm],
    exact candidates_lipschitz_aux fA }
end

/-- candidates give rise to elements of bounded_continuous_functions -/
private def candidates_b_of_candidates (f : prod_space_fun α β) (fA : f ∈ candidates α β) : Cb α β :=
bounded_continuous_function.mk_of_compact f (continuous_of_lipschitz (candidates_lipschitz fA))

private lemma candidates_b_of_candidates_mem (f : prod_space_fun α β) (fA : f ∈ candidates α β) :
  candidates_b_of_candidates f fA ∈ candidates_b α β := fA

/-- The distance on α ⊕ β is a candidate -/
private lemma dist_mem_candidates : (λp : (α ⊕ β) × (α ⊕ β), dist p.1 p.2) ∈ candidates α β :=
begin
  simp only [candidates, dist_comm, forall_const, and_true, add_comm, eq_self_iff_true,
             and_self, sum.forall, set.mem_set_of_eq, dist_self],
  repeat { split
    <|> exact (λa y z, dist_triangle_left _ _ _)
    <|> exact (λx y, by refl)
    <|> exact (λx y, max_var_bound) }
end

private def candidates_b_dist (α : Type u) (β : Type v) [metric_space α] [compact_space α] [inhabited α]
  [metric_space β] [compact_space β] [inhabited β] : Cb α β := candidates_b_of_candidates _ dist_mem_candidates

private lemma candidates_b_dist_mem_candidates_b : candidates_b_dist α β ∈ candidates_b α β :=
candidates_b_of_candidates_mem _ _

private lemma candidates_b_ne_empty : candidates_b α β ≠ ∅ :=
ne_empty_of_mem candidates_b_dist_mem_candidates_b

/-- To apply Arzela-Ascoli, we need to check that the set of candidates is closed and equicontinuous.
Equicontinuity follows from the Lipschitz control, we check closedness -/
private lemma closed_candidates_b : is_closed (candidates_b α β) :=
begin
  have I1 : ∀x y, is_closed {f : Cb α β | f (inl x, inl y) = dist x y} :=
    λx y, is_closed_eq continuous_evalx continuous_const,
  have I2 : ∀x y, is_closed {f : Cb α β | f (inr x, inr y) = dist x y } :=
    λx y, is_closed_eq continuous_evalx continuous_const,
  have I3 : ∀x y, is_closed {f : Cb α β | f (x, y) = f (y, x)} :=
    λx y, is_closed_eq continuous_evalx continuous_evalx,
  have I4 : ∀x y z, is_closed {f : Cb α β | f (x, z) ≤ f (x, y) + f (y, z)} :=
    λx y z, is_closed_le continuous_evalx (continuous_add continuous_evalx continuous_evalx),
  have I5 : ∀x, is_closed {f : Cb α β | f (x, x) = 0} :=
    λx, is_closed_eq continuous_evalx continuous_const,
  have I6 : ∀x y, is_closed {f : Cb α β | f (x, y) ≤ max_var α β} :=
    λx y, is_closed_le continuous_evalx continuous_const,
  have : candidates_b α β = (⋂x y, {f : Cb α β | f ((@inl α β x), (@inl α β y)) = dist x y})
               ∩ (⋂x y, {f : Cb α β | f ((@inr α β x), (@inr α β y)) = dist x y})
               ∩ (⋂x y, {f : Cb α β | f (x, y) = f (y, x)})
               ∩ (⋂x y z, {f : Cb α β | f (x, z) ≤ f (x, y) + f (y, z)})
               ∩ (⋂x, {f : Cb α β | f (x, x) = 0})
               ∩ (⋂x y, {f : Cb α β | f (x, y) ≤ max_var α β}) :=
    begin ext, unfold candidates_b, unfold candidates, simp [-sum.forall], refl end,
  rw this,
  repeat { apply is_closed_inter _ _
       <|> apply is_closed_Inter _
       <|> apply I1 _ _
       <|> apply I2 _ _
       <|> apply I3 _ _
       <|> apply I4 _ _ _
       <|> apply I5 _
       <|> apply I6 _ _
       <|> assume x },
end

/-- Compactness of candidates (in bounded_continuous_functions) follows -/
private lemma compact_candidates_b : compact (candidates_b α β) :=
begin
  refine arzela_ascoli₂ (Icc 0 (max_var α β)) compact_Icc (candidates_b α β) closed_candidates_b _ _,
  { rintros f ⟨x1, x2⟩ hf,
    simp only [set.mem_Icc],
    exact ⟨candidates_nonneg hf, candidates_le_max_var hf⟩ },
  { refine equicontinuous_of_continuity_modulus (λt, 2 * max_var α β * t) _ _ _,
    { have : tendsto (λ (t : ℝ), 2 * max_var α β * t) (nhds 0) (nhds (2 * max_var α β * 0)) :=
        tendsto_mul tendsto_const_nhds tendsto_id,
      simpa using this },
    { assume x y f hf,
      exact candidates_lipschitz hf _ _ }}
end

/-- We will then choose the candidate minimizing the Hausdorff distance. Except that we are not
in a metric space setting, so we need to define our custom version of Hausdorff distance,
called HD, and prove its basic properties. -/
private def HD (f : Cb α β) := max (supr (λx:α, infi (λy:β, f (inl x, inr y))))
                                   (supr (λy:β, infi (λx:α, f (inl x, inr y))))

/- We will show that HD is continuous on bounded_continuous_functions, to deduce that its
minimum on the compact set candidates_b is attained. Since it is defined in terms of
infimum and supremum on ℝ, which is only conditionnally complete, we will need all the time
to check that the defining sets are bounded below or above. This is done in the next few
technical lemmas -/

private lemma HD_below_aux1 {f : Cb α β} (C : ℝ) {x : α} : bdd_below (range (λ (y : β), f (inl x, inr y) + C)) :=
let ⟨cf, hcf⟩ := (real.bounded_iff_bdd_below_bdd_above.1 bounded_range).1 in
⟨cf + C, forall_range_iff.2 (λi, add_le_add_right ((λx, hcf (f x) (mem_range_self _)) _) _)⟩

private lemma HD_bound_aux1 (f : Cb α β) (C : ℝ) : bdd_above (range (λ (x : α), infi (λy:β, f (inl x, inr y) + C))) :=
begin
  rcases (real.bounded_iff_bdd_below_bdd_above.1 bounded_range).2 with ⟨Cf, hCf⟩,
  refine ⟨Cf + C, forall_range_iff.2 (λx, _)⟩,
  calc infi (λy:β, f (inl x, inr y) + C) ≤ f (inl x, inr (default β)) + C :
    cinfi_le (HD_below_aux1 C)
    ... ≤ Cf + C : add_le_add ((λx, hCf (f x) (mem_range_self _)) _) (le_refl _)
end

private lemma HD_below_aux2 {f : Cb α β} (C : ℝ) {y : β} : bdd_below (range (λ (x : α), f (inl x, inr y) + C)) :=
let ⟨cf, hcf⟩ := (real.bounded_iff_bdd_below_bdd_above.1 bounded_range).1 in
⟨cf + C, forall_range_iff.2 (λi, add_le_add_right ((λx, hcf (f x) (mem_range_self _)) _) _)⟩

private lemma HD_bound_aux2 (f : Cb α β) (C : ℝ) : bdd_above (range (λ (y : β), infi (λx:α, f (inl x, inr y) + C))) :=
begin
  rcases (real.bounded_iff_bdd_below_bdd_above.1 bounded_range).2 with ⟨Cf, hCf⟩,
  refine ⟨Cf + C, forall_range_iff.2 (λy, _)⟩,
  calc infi (λx:α, f (inl x, inr y) + C) ≤ f (inl (default α), inr y) + C :
    cinfi_le (HD_below_aux2 C)
  ... ≤ Cf + C : add_le_add ((λx, hCf (f x) (mem_range_self _)) _) (le_refl _)
end

/-- Explicit bound on HD (dist). This means that when looking for minimizers it will
be sufficient to look for functions with HD(f) bounded by this bound. -/
private lemma HD_candidates_b_dist_le : HD (candidates_b_dist α β) ≤ diam (univ : set α) + 1 + diam (univ : set β) :=
begin
  refine max_le (csupr_le (λx, _)) (csupr_le (λy, _)),
  { have A : infi (λy:β, candidates_b_dist α β (inl x, inr y)) ≤ candidates_b_dist α β (inl x, inr (default β)) :=
      cinfi_le (by simpa using HD_below_aux1 0),
    have B : dist (inl x) (inr (default β)) ≤ diam (univ : set α) + 1 + diam (univ : set β) := calc
      dist (inl x) (inr (default β)) = dist x (default α) + 1 + dist (default β) (default β) : rfl
      ... ≤ diam (univ : set α) + 1 + diam (univ : set β) :
      begin
        apply add_le_add (add_le_add _ (le_refl _)),
        exact dist_le_diam_of_mem (bounded_of_compact (compact_univ)) (mem_univ _) (mem_univ _),
        exact dist_le_diam_of_mem (bounded_of_compact (compact_univ)) (mem_univ _) (mem_univ _)
      end,
    exact le_trans A B },
  { have A : infi (λx:α, candidates_b_dist α β (inl x, inr y)) ≤ candidates_b_dist α β (inl (default α), inr y) :=
      cinfi_le (by simpa using HD_below_aux2 0),
    have B : dist (inl (default α)) (inr y) ≤ diam (univ : set α) + 1 + diam (univ : set β) := calc
      dist (inl (default α)) (inr y) = dist (default α) (default α) + 1 + dist (default β) y : rfl
      ... ≤ diam (univ : set α) + 1 + diam (univ : set β) :
      begin
        apply add_le_add (add_le_add _ (le_refl _)),
        exact dist_le_diam_of_mem (bounded_of_compact (compact_univ)) (mem_univ _) (mem_univ _),
        exact dist_le_diam_of_mem (bounded_of_compact (compact_univ)) (mem_univ _) (mem_univ _)
      end,
    exact le_trans A B },
end

/- To check that HD is continuous, we check that it is Lipschitz. As HD is a max, we
prove separately inequalities controlling the two terms (relying too heavily on copy-paste...) -/
private lemma HD_lipschitz_aux1 (f g : Cb α β) :
  supr (λx:α, infi (λy:β, f (inl x, inr y))) ≤ supr (λx:α, infi (λy:β, g (inl x, inr y))) + dist f g :=
begin
  rcases (real.bounded_iff_bdd_below_bdd_above.1 bounded_range).1 with ⟨cg, hcg⟩,
  have Hcg : ∀x, cg ≤ g x := λx, hcg (g x) (mem_range_self _),
  rcases (real.bounded_iff_bdd_below_bdd_above.1 bounded_range).1 with ⟨cf, hcf⟩,
  have Hcf : ∀x, cf ≤ f x := λx, hcf (f x) (mem_range_self _),

  -- prove the inequality but with `dist f g` inside, by using inequalities comparing
  -- supr to supr and infi to infi
  have Z : supr (λx:α, infi (λy:β, f (inl x, inr y))) ≤ supr (λx:α, infi (λy:β, g (inl x, inr y) + dist f g)) :=
    csupr_le_csupr (HD_bound_aux1 _ (dist f g))
      (λx, cinfi_le_cinfi ⟨cf, forall_range_iff.2(λi, Hcf _)⟩ (λy, coe_le_coe_add_dist)),
  -- move the `dist f g` out of the infimum and the supremum, arguing that continuous monotone maps
  -- (here the addition of `dist f g`) preserve infimum and supremum
  have E1 : ∀x, infi (λy:β, g (inl x, inr y)) + dist f g =
             infi ((λz, z + dist f g) ∘ (λy:β, (g (inl x, inr y)))),
  { assume x,
    refine cinfi_of_cinfi_of_monotone_of_continuous (_ : continuous (λ (z : ℝ), z + dist f g)) _ _,
    { exact continuous_add continuous_id continuous_const },
    { assume x y hx, simpa },
    { show bdd_below (range (λ (y : β), g (inl x, inr y))),
        from ⟨cg, forall_range_iff.2(λi, Hcg _)⟩ }},
  have E2 : supr (λx:α, infi (λy:β, g (inl x, inr y))) + dist f g =
         supr ((λz, z + dist f g) ∘ (λx:α, infi (λy:β, g (inl x, inr y)))),
  { refine csupr_of_csupr_of_monotone_of_continuous (_ : continuous (λ (z : ℝ), z + dist f g)) _ _,
    { exact continuous_add continuous_id continuous_const },
    { assume x y hx, simpa },
    { by simpa using HD_bound_aux1 _ 0 }},
  -- deduce the result from the above two steps
  simp only [add_comm] at Z,
  simpa [E2, E1, function.comp]
end

private lemma HD_lipschitz_aux2 (f g : Cb α β) :
  supr (λy:β, infi (λx:α, f (inl x, inr y))) ≤ supr (λy:β, infi (λx:α, g (inl x, inr y))) + dist f g :=
begin
  rcases (real.bounded_iff_bdd_below_bdd_above.1 bounded_range).1 with ⟨cg, hcg⟩,
  have Hcg : ∀x, cg ≤ g x := λx, hcg (g x) (mem_range_self _),
  rcases (real.bounded_iff_bdd_below_bdd_above.1 bounded_range).1 with ⟨cf, hcf⟩,
  have Hcf : ∀x, cf ≤ f x := λx, hcf (f x) (mem_range_self _),

  -- prove the inequality but with `dist f g` inside, by using inequalities comparing
  -- supr to supr and infi to infi
  have Z : supr (λy:β, infi (λx:α, f (inl x, inr y))) ≤ supr (λy:β, infi (λx:α, g (inl x, inr y) + dist f g)) :=
    csupr_le_csupr (HD_bound_aux2 _ (dist f g))
      (λy, cinfi_le_cinfi  ⟨cf, forall_range_iff.2(λi, Hcf _)⟩ (λy, coe_le_coe_add_dist)),
  -- move the `dist f g` out of the infimum and the supremum, arguing that continuous monotone maps
  -- (here the addition of `dist f g`) preserve infimum and supremum
  have E1 : ∀y, infi (λx:α, g (inl x, inr y)) + dist f g =
             infi ((λz, z + dist f g) ∘ (λx:α, (g (inl x, inr y)))),
  { assume y,
    refine cinfi_of_cinfi_of_monotone_of_continuous (_ : continuous (λ (z : ℝ), z + dist f g)) _ _,
    { exact continuous_add continuous_id continuous_const },
    { assume x y hx, simpa },
    { show bdd_below (range (λx:α, g (inl x, inr y))),
        from ⟨cg, forall_range_iff.2(λi, Hcg _)⟩ }},
  have E2 : supr (λy:β, infi (λx:α, g (inl x, inr y))) + dist f g =
         supr ((λz, z + dist f g) ∘ (λy:β, infi (λx:α, g (inl x, inr y)))),
  { refine csupr_of_csupr_of_monotone_of_continuous (_ : continuous (λ (z : ℝ), z + dist f g)) _ _,
    { exact continuous_add continuous_id continuous_const },
    { assume x y hx, simpa },
    { by simpa using HD_bound_aux2 _ 0 }},
  -- deduce the result from the above two steps
  simp only [add_comm] at Z,
  simpa [E2, E1, function.comp]
end

private lemma HD_lipschitz_aux3 (f g : Cb α β) : HD f ≤ HD g + dist f g :=
max_le (le_trans (HD_lipschitz_aux1 f g) (add_le_add_right (le_max_left _ _) _))
       (le_trans (HD_lipschitz_aux2 f g) (add_le_add_right (le_max_right _ _) _))

/-- Conclude that HD, being Lipschitz, is continuous -/
private lemma HD_continuous : continuous (HD : Cb α β → ℝ) :=
uniform_continuous.continuous $ uniform_continuous_of_le_add 1 $
λf g, begin simp, exact HD_lipschitz_aux3 _ _ end

end constructions --section

section consequences
variables (α : Type u) (β : Type v) [metric_space α] [compact_space α] [nonempty α] [metric_space β] [compact_space β] [nonempty β]

/- Now that we have proved that the set of candidates is compact, and that HD is continuous,
we can finally select a candidate minimizing HD. This will be the candidate realizing the
optimal coupling. -/
private lemma exists_minimizer : ∃f ∈ candidates_b α β, ∀g ∈ candidates_b α β, HD f ≤ HD g :=
exists_forall_le_of_compact_of_continuous _ HD_continuous _ compact_candidates_b candidates_b_ne_empty

private definition optimal_GH_dist : Cb α β := classical.some (exists_minimizer α β)

private lemma optimal_GH_dist_mem_candidates_b : optimal_GH_dist α β ∈ candidates_b α β :=
by cases (classical.some_spec (exists_minimizer α β)); assumption

private lemma HD_optimal_GH_dist_le (g : Cb α β) (hg : g ∈ candidates_b α β) : HD (optimal_GH_dist α β) ≤ HD g :=
let ⟨Z1, Z2⟩ := classical.some_spec (exists_minimizer α β) in Z2 g hg

/-- With the optimal candidate, construct a premetric space structure on α ⊕ β, on which the
predistance is given by the candidate. Then, we will identify points at 0 predistance
to obtain a genuine metric space -/
def premetric_optimal_GH_dist : premetric_space (α ⊕ β) :=
{ dist := λp q, optimal_GH_dist α β (p, q),
  dist_self := λx, candidates_refl (optimal_GH_dist_mem_candidates_b α β),
  dist_comm := λx y, candidates_symm (optimal_GH_dist_mem_candidates_b α β),
  dist_triangle := λx y z, candidates_triangle (optimal_GH_dist_mem_candidates_b α β) }

local attribute [instance] premetric_optimal_GH_dist

/-- A metric space which realizes the optimal coupling between α and β -/
@[reducible] definition optimal_GH_coupling : Type* :=
metric_quot (α ⊕ β)

instance : metric_space (optimal_GH_coupling α β) := by apply_instance

private lemma optimal_GH_dist.dist_eq (p q : α ⊕ β) : dist ⟦p⟧ ⟦q⟧ = (optimal_GH_dist α β).val (p, q) := rfl

/-- Injection of α in the optimal coupling between α and β -/
def optimal_GH_injl (x : α) : optimal_GH_coupling α β := ⟦inl x⟧

/-- The injection of α in the optimal coupling between α and β is an isometry. -/
lemma isometry_optimal_GH_injl : isometry_on (optimal_GH_injl α β) univ :=
begin
  assume x y hx hy,
  change dist ⟦inl x⟧ ⟦inl y⟧ = dist x y,
  rw [optimal_GH_dist.dist_eq α β],
  exact candidates_dist_inl (optimal_GH_dist_mem_candidates_b α β) _ _,
end

/-- Injection of β  in the optimal coupling between α and β -/
def optimal_GH_injr (y : β) : optimal_GH_coupling α β := ⟦inr y⟧

/-- The injection of β in the optimal coupling between α and β is an isometry. -/
lemma isometry_optimal_GH_injr : isometry_on (optimal_GH_injr α β) univ :=
begin
  assume x y hx hy,
  change dist ⟦inr x⟧ ⟦inr y⟧ = dist x y,
  rw [optimal_GH_dist.dist_eq α β],
  exact candidates_dist_inr (optimal_GH_dist_mem_candidates_b α β) _ _,
end

/-- The optimal coupling between two compact spaces α and β is still a compact space -/
instance compact_space_optimal_GH_coupling : compact_space (optimal_GH_coupling α β) :=
⟨begin
  have : (univ : set (optimal_GH_coupling α β)) = (optimal_GH_injl α β '' univ) ∪ (optimal_GH_injr α β '' univ),
  { refine subset.antisymm (λxc hxc, _) (subset_univ _),
    rcases quotient.exists_rep xc with ⟨x, hx⟩,
    cases x; rw ← hx,
    { have : ⟦inl x⟧ = optimal_GH_injl α β x := rfl,
      rw this,
      exact mem_union_left _ (mem_image_of_mem _ (mem_univ _)) },
    { have : ⟦inr x⟧ = optimal_GH_injr α β x := rfl,
      rw this,
      exact mem_union_right _ (mem_image_of_mem _ (mem_univ _)) }},
  rw this,
  exact compact_union_of_compact
    (compact_image (compact_univ) (continuous_of_isometry (isometry_optimal_GH_injl α β)))
    (compact_image (compact_univ) (continuous_of_isometry (isometry_optimal_GH_injr α β)))
end⟩

/-- For any candidate f, HD(f) is larger than or equal to the Hausdorff distance in the
optimal coupling. This follows from the fact that HD of the optimal candidate is exactly
the Hausdorff distance in the optimal coupling, although we only prove here the inequality
we need. -/
private lemma Hausdorff_dist_optimal_le_HD {f} (h : f ∈ candidates_b α β) :
  Hausdorff_dist (range (optimal_GH_injl α β)) (range (optimal_GH_injr α β)) ≤ HD f :=
begin
  refine le_trans (le_of_forall_le_of_dense (λr hr, _)) (HD_optimal_GH_dist_le α β f h),
  have A : ∀ x ∈ range (optimal_GH_injl α β), ∃ y ∈ range (optimal_GH_injr α β), dist x y ≤ r,
  { assume x hx,
    rcases mem_range.1 hx with ⟨z, hz⟩,
    rw ← hz,
    have I1 : supr (λx:α, infi (λy:β, optimal_GH_dist α β (inl x, inr y))) < r :=
      lt_of_le_of_lt (le_max_left _ _) hr,
    have I2 : infi (λy:β, optimal_GH_dist α β (inl z, inr y)) ≤ supr (λx:α, infi (λy:β, optimal_GH_dist α β (inl x, inr y))) :=
      le_cSup (by simpa using HD_bound_aux1 _ 0) (mem_range_self _),
    have I : infi (λy:β, optimal_GH_dist α β (inl z, inr y)) < r := lt_of_le_of_lt I2 I1,
    rcases exists_lt_of_cInf_lt (by simpa) I with ⟨r', r'range, hr'⟩,
    rcases mem_range.1 r'range with ⟨z', hz'⟩,
    existsi [optimal_GH_injr α β z', mem_range_self _],
    have : (optimal_GH_dist α β) (inl z, inr z') ≤ r := begin rw hz', exact le_of_lt hr' end,
    exact this },
  refine Hausdorff_dist.le_of_mem_dist_le_of_mem_dist_le _ _ _,
  { rcases exists_mem_of_nonempty α with ⟨xα, _⟩,
    have : optimal_GH_injl α β xα ∈ range (optimal_GH_injl α β) := mem_range_self _,
    rcases A _ this with ⟨y, yrange, hy⟩,
    exact le_trans dist_nonneg hy },
  { exact A },
  { assume y hy,
    rcases mem_range.1 hy with ⟨z, hz⟩,
    rw ← hz,
    have I1 : supr (λy:β, infi (λx:α, optimal_GH_dist α β (inl x, inr y))) < r :=
      lt_of_le_of_lt (le_max_right _ _) hr,
    have I2 : infi (λx:α, optimal_GH_dist α β (inl x, inr z)) ≤ supr (λy:β, infi (λx:α, optimal_GH_dist α β (inl x, inr y))) :=
      le_cSup (by simpa using HD_bound_aux2 _ 0) (mem_range_self _),
    have I : infi (λx:α, optimal_GH_dist α β (inl x, inr z)) < r := lt_of_le_of_lt I2 I1,
    rcases exists_lt_of_cInf_lt (by simpa) I with ⟨r', r'range, hr'⟩,
    rcases mem_range.1 r'range with ⟨z', hz'⟩,
    existsi [optimal_GH_injl α β z', mem_range_self _],
    have : (optimal_GH_dist α β) (inl z', inr z) ≤ r := begin rw hz', exact le_of_lt hr' end,
    rw dist_comm,
    exact this }
end

end consequences
/- We are done with the construction of the optimal coupling, prerequisiste number 1-/
end Gromov_Hausdorff_realized --section


section Kuratowski_embedding
variables {α : Type u} {f g : l_infty_ℝ} {n : nat} {C : real} [metric_space α] (x : ℕ → α) (a b : α)

/-- A metric space can be embedded in `l^∞_ℝ` via the distances to points in
a fixed countable set, if this set is dense. This map is given in the next definition,
without density assumptions. -/
def embedding_of_subset : l_infty_ℝ :=
of_normed_group_discrete (λn, dist a (x n) - dist (x 0) (x n)) (dist a (x 0)) (λ_, abs_dist_sub_le _ _ _)

lemma embedding_of_subset_coe : embedding_of_subset x a n = dist a (x n) - dist (x 0) (x n) := rfl

/-- The embedding map is always a semi-contraction. -/
lemma embedding_of_subset_dist_le (a b : α) : dist (embedding_of_subset x a) (embedding_of_subset x b) ≤ dist a b :=
begin
  refine (dist_le dist_nonneg).2 (λn, _),
  have A : dist a (x n) + (dist (x 0) (x n) + (-dist b (x n) + -dist (x 0) (x n))) = dist a (x n) - dist b (x n) := by ring,
  simp only [embedding_of_subset_coe, real.dist_eq, A, add_comm, neg_add_rev, _root_.neg_neg, sub_eq_add_neg, add_left_comm],
  exact abs_dist_sub_le _ _ _
end

/-- When the reference set is dense, the embedding map is an isometry on its image. -/
lemma embedding_of_subset_isometry (H : closure (range x) = univ) : isometry_on (embedding_of_subset x) univ :=
begin
  assume a b _ _,
  refine le_antisymm (embedding_of_subset_dist_le x a b) (real.le_of_forall_epsilon_le (λe epos, _)),
  /- First step: find n with dist a (x n) < e -/
  have A : a ∈ closure (range x) := by have B := mem_univ a; rwa [← H] at B,
  rcases mem_closure_iff'.1 A (e/2) (half_pos epos) with ⟨d, ⟨drange, hd⟩⟩,
  cases drange with n dn,
  rw [← dn] at hd,
  /- Second step: use the norm control at index n to conclude -/
  have C : dist b (x n) - dist a (x n) = embedding_of_subset x b n - embedding_of_subset x a n :=
    by simp [embedding_of_subset_coe]; ring,
  have := calc
    dist a b ≤ dist a (x n) + dist (x n) b : dist_triangle _ _ _
    ...    = 2 * dist a (x n) + (dist b (x n) - dist a (x n)) : by simp [dist_comm]; ring
    ...    ≤ 2 * dist a (x n) + abs (dist b (x n) - dist a (x n)) :
      by apply_rules [add_le_add_left, le_abs_self]
    ...    ≤ 2 * (e/2) + abs (embedding_of_subset x b n - embedding_of_subset x a n) :
      begin rw [C], apply_rules [add_le_add, mul_le_mul_of_nonneg_left, le_of_lt hd, le_refl], norm_num end
    ...    ≤ 2 * (e/2) + dist (embedding_of_subset x b) (embedding_of_subset x a) :
      begin rw [← coe_diff], apply add_le_add_left, rw [coe_diff, ←real.dist_eq], apply dist_coe_le_dist end
    ...    = dist (embedding_of_subset x b) (embedding_of_subset x a) + e : by ring,
  simpa [dist_comm] using this
end

/-- Every separable metric space embeds isometrically in l_infty_ℝ. -/
theorem exists_isometric_embedding (α : Type u) [metric_space α] [separable_space α] :
  ∃(f : α → l_infty_ℝ), isometry_on f univ :=
begin
  by_cases h : (univ : set α) = ∅,
  { simp [isometry_on_empty, h] },
  { /- We construct a map x : ℕ → α with dense image -/
    rcases exists_mem_of_ne_empty h with basepoint,
    haveI : inhabited α := ⟨basepoint⟩,
    have : ∃s:set α, countable s ∧ closure s = univ := separable_space.exists_countable_closure_eq_univ _,
    rcases this with ⟨S, ⟨S_countable, S_dense⟩⟩,
    rcases countable_iff_exists_surjective.1 S_countable with ⟨x, x_range⟩,
    have : closure (range x) = univ := univ_subset_iff.1 (by rw [← S_dense]; apply closure_mono; assumption),
    /- Use embedding_of_subset to construct the desired isometry -/
    exact ⟨embedding_of_subset x, embedding_of_subset_isometry x this⟩ }
end

/-- The Kuratowski embedding is an isometric embedding of a separable metric space in ℓ^∞(ℝ) -/
def Kuratowski_embedding (α : Type u) [metric_space α] [separable_space α] : α → l_infty_ℝ :=
  classical.some (exists_isometric_embedding α)

/-- The Kuratowski embedding is an isometry -/
lemma Kuratowski_embedding_isometry {α : Type u} [metric_space α] [separable_space α] :
  isometry_on (Kuratowski_embedding α) univ :=
classical.some_spec (exists_isometric_embedding α)

/-- Version of the Kuratowski embedding for nonempty compacts -/
def nonempty_compacts.Kuratowski_embedding (α : Type u) [metric_space α] [compact_space α] [i : nonempty α] :
  nonempty_compacts l_infty_ℝ :=
⟨range (Kuratowski_embedding α),
begin
  split,
  { rcases exists_mem_of_nonempty α with ⟨x, hx⟩,
    have A : Kuratowski_embedding α x ∈ range (Kuratowski_embedding α) := ⟨x, by simp⟩,
    apply ne_empty_of_mem A },
  { rw ← image_univ,
    exact compact_image compact_univ (continuous_of_isometry Kuratowski_embedding_isometry) },
end⟩

end Kuratowski_embedding

section GH_space
/- In this section, we define the Gromov-Hausdorff space, denoted `GH_space` as the quotient
of nonempty compact subsets of ℓ^∞(ℝ) identifying isometric sets.
Using the Kuratwoski embedding, we get a canonical map `to_GH_space` mapping any nonempty
compact type to GH_space. -/

/-- Equivalence relation identifying two nonempty compact sets which are isometric -/
private definition isometry_rel : nonempty_compacts l_infty_ℝ → nonempty_compacts l_infty_ℝ → Prop :=
  λx y, ∃Φ : l_infty_ℝ → l_infty_ℝ, isometry_on Φ x.val ∧ Φ '' x.val = y.val

/-- This is indeed an equivalence relation -/
private lemma is_equivalence_isometry_rel : equivalence isometry_rel :=
begin
  repeat { split },
  { exact λx, ⟨id, ⟨λx y hx hy, by simp, by ext; simp⟩⟩ },
  { rintros x y ⟨Φ, ⟨isomΦ, imΦ⟩⟩,
    existsi (inv_fun_on Φ x.val),
    split,
    { rw ← imΦ,
      apply isometry_on_inv_fun_on isomΦ },
    { rw ← imΦ,
      apply inv_fun_on_image (inj_on_of_isometry_on isomΦ) (subset.refl _) }},
  { rintros x y z ⟨Φ, ⟨isomΦ, imΦ⟩⟩  ⟨Ψ, ⟨isomΨ, imΨ⟩⟩,
    existsi (Ψ ∘ Φ),
    split,
    { apply isometry_on_comp _ isomΨ isomΦ,
      rw [maps_to', imΦ] },
    { rw [image_comp, imΦ, imΨ] }}
end

/-- setoid instance identifying two isometric nonempty compact subspaces of ℓ^∞(ℝ) -/
instance isometry_rel.setoid : setoid (nonempty_compacts l_infty_ℝ) :=
setoid.mk isometry_rel is_equivalence_isometry_rel

/-- The Gromov-Hausdorff space -/
definition GH_space : Type := quotient (isometry_rel.setoid)

/-- Map any nonempty compact type to GH_space -/
definition to_GH_space (α : Type u) [metric_space α] [compact_space α] [nonempty α] : GH_space :=
  ⟦nonempty_compacts.Kuratowski_embedding α⟧

definition rep_GH_space (p : GH_space) : Type := (some (quot.exists_rep p)).val

lemma eq_to_GH_space_iff {α : Type u} [metric_space α] [compact_space α] [nonempty α] {p : nonempty_compacts l_infty_ℝ} :
  ⟦p⟧ = to_GH_space α ↔ ∃Ψ : α → l_infty_ℝ, isometry_on Ψ univ ∧ range Ψ = p.val :=
begin
  simp only [to_GH_space, quotient.eq],
  split,
  { assume h,
    have h' := setoid.symm h,
    change ∃Φ : l_infty_ℝ → l_infty_ℝ, isometry_on Φ (range (Kuratowski_embedding α))
            ∧ Φ '' (range (Kuratowski_embedding α)) = p.val at h',
    rcases h' with ⟨Φ, ⟨Φisom, Φrange⟩⟩,
    existsi (Φ ∘ (Kuratowski_embedding α)),
    split,
    { exact isometry_on_comp (maps_to_range _ _) Φisom Kuratowski_embedding_isometry },
    { by rwa range_comp }},
  { rintros ⟨Ψ, ⟨isomΨ, rangeΨ⟩⟩,
    letI : inhabited α := classical.inhabited_of_nonempty (by assumption),
    change (∃Φ : l_infty_ℝ → l_infty_ℝ, isometry_on Φ p.val ∧ Φ '' p.val = range (Kuratowski_embedding α)),
    existsi (Kuratowski_embedding α) ∘ (inv_fun_on Ψ univ),
    split,
    { apply isometry_on_comp (maps_to_univ _ _) Kuratowski_embedding_isometry,
      rw [← rangeΨ, ← image_univ],
      apply isometry_on_inv_fun_on isomΨ },
    { rw [← rangeΨ, ← image_univ, image_comp, inv_fun_on_image, image_univ],
      { exact inj_on_of_isometry_on isomΨ },
      {exact subset.refl _ }}}
end

lemma eq_to_GH_space {p : nonempty_compacts l_infty_ℝ} : ⟦p⟧ = to_GH_space p.val :=
eq_to_GH_space_iff.2 ⟨((λx, x) : p.val → l_infty_ℝ), isometry_on_subtype.val, subtype_val_range⟩

section
local attribute [reducible] rep_GH_space

instance rep_GH_space_metric_space {p : GH_space} : metric_space (rep_GH_space p) :=
by apply_instance

instance rep_GH_space_compact_space {p : GH_space} : compact_space (rep_GH_space p) :=
by apply_instance

instance rep_GH_space_nonempty {p : GH_space} : nonempty (rep_GH_space p) :=
by apply_instance
end

lemma to_GH_space_rep_GH_space (p : GH_space) : to_GH_space (rep_GH_space p) = p :=
begin
  change to_GH_space ((some (quot.exists_rep p)).val) = p,
  rw ← eq_to_GH_space,
  exact some_spec ((quot.exists_rep p))
end

/-- Two nonempty compact spaces have the same image in GH_space if and only if they are isometric -/
lemma to_GH_space_eq_to_GH_space_iff_isometric {α : Type u} [metric_space α] [compact_space α] [nonempty α]
  {β : Type u} [metric_space β] [compact_space β] [nonempty β] :
  to_GH_space α = to_GH_space β ↔ ∃Φ : α → β, isometry_on Φ univ ∧ range Φ = univ :=
⟨begin
  assume h,
  rcases eq_to_GH_space_iff.1 h with ⟨Ψ, ⟨Ψisom, Ψrange⟩⟩,
  let Φ := Kuratowski_embedding α,
  have Ψrange' : range Ψ = range Φ := Ψrange,
  letI : inhabited β := inhabited_of_nonempty',
  existsi ((inv_fun_on Ψ univ) ∘ Φ),
  split,
  { apply isometry_on_comp (maps_to_range _ univ) _ Kuratowski_embedding_isometry,
    rw [← Ψrange', ← image_univ],
    exact isometry_on_inv_fun_on Ψisom },
  { rw [← image_univ, image_comp, image_univ, ← Ψrange', ← image_univ],
    apply inv_fun_on_image (inj_on_of_isometry_on Ψisom) (subset_univ _) }
end,
begin
  rintros ⟨Φ, ⟨Φisom, Φrange⟩⟩,
  apply eq.symm,
  erw [eq_to_GH_space_iff],
  let Ψ := Kuratowski_embedding β,
  existsi (Ψ ∘ Φ),
  split,
  { exact isometry_on_comp (maps_to_univ _ _) Kuratowski_embedding_isometry Φisom },
  { rw [← image_univ, image_comp, image_univ, Φrange, image_univ], refl }
end⟩

/-- Distance on GH_space : the distance between two nonempty compact spaces is the infimum
Hausdorff distance between isometric copies of the two spaces in a metric space. For the definition,
we only consider embeddings in ℓ^∞(ℝ), but we will prove below that it works for all spaces. -/
instance : has_dist (GH_space) :=
{ dist := λx y, Inf ((λp : nonempty_compacts l_infty_ℝ × nonempty_compacts l_infty_ℝ, Hausdorff_dist p.1.val p.2.val) ''
                      (set.prod {a | ⟦a⟧ = x} {b | ⟦b⟧ = y})) }

def GH_dist (α : Type u) (β : Type v) [metric_space α] [nonempty α] [compact_space α]
  [metric_space β] [nonempty β] [compact_space β] : ℝ := dist (to_GH_space α) (to_GH_space β)

lemma dist_GH_dist (p q : GH_space) : dist p q = GH_dist (rep_GH_space p) (rep_GH_space q) :=
by rw [GH_dist, to_GH_space_rep_GH_space, to_GH_space_rep_GH_space]

/-- The Gromov-Hausdorff distance between two spaces is bounded by the Hausdorff distance
of isometric copies of the spaces, in any metric space. -/
theorem GH_dist_le_Hausdorff_dist {α : Type u} [metric_space α] [compact_space α] [nonempty α]
  {β : Type v} [metric_space β] [compact_space β] [nonempty β]
  {γ : Type w} [metric_space γ] {Φ : α → γ} {Ψ : β → γ} (ha : isometry_on Φ univ) (hb : isometry_on Ψ univ) :
  GH_dist α β ≤ Hausdorff_dist (range Φ) (range Ψ) :=
begin
  /- For the proof, we want to embed γ in ℓ^∞(ℝ), to say that the Hausdorff distance is realized
  in ℓ^∞(ℝ) and therefore bounded below by the Gromov-Hausdorff-distance. However, γ is not
  separable in general. We restrict to the union of the images of α and β in γ, which is
  separable and therefore embeddable in ℓ^∞(ℝ). -/
  rcases exists_mem_of_nonempty α with ⟨xα, _⟩,
  letI : inhabited α := ⟨xα⟩,
  letI : inhabited β := classical.inhabited_of_nonempty (by assumption),
  let s : set γ := (range Φ) ∪ (range Ψ),
  let Φ' : α → subtype s := λy, ⟨Φ y, mem_union_left _ (mem_range_self _)⟩,
  let Ψ' : β → subtype s := λy, ⟨Ψ y, mem_union_right _ (mem_range_self _)⟩,
  have : isometry_on Φ' univ := λx y hx hy, ha hx hy,
  have : isometry_on Ψ' univ := λx y hx hy, hb hx hy,
  have : compact s,
  { apply compact_union_of_compact,
    { rw ← image_univ,
      apply compact_image (compact_univ) (continuous_of_isometry ha) },
    { rw ← image_univ,
      apply compact_image (compact_univ) (continuous_of_isometry hb) }},
  letI : metric_space (subtype s) := by apply_instance,
  haveI : compact_space (subtype s) := ⟨compact_iff_compact_univ.1 ‹compact s›⟩,
  haveI : nonempty (subtype s) := ⟨Φ' xα⟩,
  have ΦΦ' : Φ = subtype.val ∘ Φ' := by funext; refl,
  have ΨΨ' : Ψ = subtype.val ∘ Ψ' := by funext; refl,
  have : Hausdorff_dist (range Φ) (range Ψ) = Hausdorff_dist (range Φ') (range Ψ'),
  { rw [ΦΦ', ΨΨ', range_comp, range_comp],
    exact Hausdorff_dist.image (subset_univ _) (subset_univ _) (isometry_on_subtype.val) },
  rw this,
  -- Embed s in ℓ^∞(ℝ) through its Kuratowski embedding
  let F := Kuratowski_embedding (subtype s),
  have : Hausdorff_dist (F '' (range Φ')) (F '' (range Ψ')) = Hausdorff_dist (range Φ') (range Ψ') :=
    Hausdorff_dist.image (subset_univ _) (subset_univ _) Kuratowski_embedding_isometry,
  rw ← this,
  -- Let A and B be the images of α and β under this embedding. They are in ℓ^∞(ℝ), and
  -- their Hausdorff distance is the same as in the original space.
  let A : nonempty_compacts l_infty_ℝ := ⟨F '' (range Φ'), ⟨by simp, begin
      rw [← range_comp, ← image_univ],
      exact compact_image compact_univ
        (continuous.comp (continuous_of_isometry ‹isometry_on Φ' univ›) (continuous_of_isometry Kuratowski_embedding_isometry)),
    end⟩⟩,
  let B : nonempty_compacts l_infty_ℝ := ⟨F '' (range Ψ'), ⟨by simp, begin
      rw [← range_comp, ← image_univ],
      exact compact_image compact_univ
        (continuous.comp (continuous_of_isometry ‹isometry_on Ψ' univ›) (continuous_of_isometry Kuratowski_embedding_isometry)),
    end⟩⟩,
  have Aα : ⟦A⟧ = to_GH_space α,
  { rw [to_GH_space, quotient.eq],
    change ∃G : l_infty_ℝ → l_infty_ℝ, isometry_on G (F '' (range Φ')) ∧ G '' (F '' (range Φ')) = range (Kuratowski_embedding α),
    existsi (Kuratowski_embedding α) ∘ (inv_fun_on (F ∘ Φ') (univ : set α)),
    split,
    { apply isometry_on_comp (maps_to_univ _ _) Kuratowski_embedding_isometry,
      rw [← range_comp, ← image_univ],
      exact isometry_on_inv_fun_on (isometry_on_comp (maps_to_univ _ _) Kuratowski_embedding_isometry ‹isometry_on Φ' univ›) },
    { rw [← range_comp, image_comp, ← image_univ, inv_fun_on_image _ (subset.refl _), image_univ],
      exact inj_on_of_isometry_on (isometry_on_comp (maps_to_univ _ _) Kuratowski_embedding_isometry ‹isometry_on Φ' univ›) }},
  have Bβ : ⟦B⟧ = to_GH_space β,
  { rw [to_GH_space, quotient.eq],
    change ∃G : l_infty_ℝ → l_infty_ℝ, isometry_on G (F '' (range Ψ')) ∧ G '' (F '' (range Ψ')) = range (Kuratowski_embedding β),
    existsi (Kuratowski_embedding β) ∘ (inv_fun_on (F ∘ Ψ') (univ : set β)),
    split,
    { apply isometry_on_comp (maps_to_univ _ _) Kuratowski_embedding_isometry,
      rw [← range_comp, ← image_univ],
      exact isometry_on_inv_fun_on (isometry_on_comp (maps_to_univ _ _) Kuratowski_embedding_isometry ‹isometry_on Ψ' univ›) },
    { rw [← range_comp, image_comp, ← image_univ, inv_fun_on_image _ (subset.refl _), image_univ],
      exact inj_on_of_isometry_on (isometry_on_comp (maps_to_univ _ _) Kuratowski_embedding_isometry ‹isometry_on Ψ' univ›) }},
  refine cInf_le ⟨0, begin simp, assume t _ _ _ _ ht, rw ← ht, exact Hausdorff_dist.nonneg end⟩ _,
  apply (mem_image _ _ _).2,
  existsi (⟨A, B⟩ : nonempty_compacts l_infty_ℝ × nonempty_compacts l_infty_ℝ),
  simp [Aα, Bβ],
end

local attribute [instance, priority 0] inhabited_of_nonempty'

/-- The optimal coupling constructed above realizes exactly the Gromov-Hausdorff distance,
essentially by design. -/
lemma Hausdorff_dist_optimal {α : Type u} [metric_space α] [compact_space α] [nonempty α]
  {β : Type v} [metric_space β] [compact_space β] [nonempty β] :
  Hausdorff_dist (range (optimal_GH_injl α β)) (range (optimal_GH_injr α β)) = GH_dist α β :=
begin
  /- we only need to check the inequality ≤, as the other one follows from the previous lemma.
     As the Gromov-Hausdorff distance is an infimum, we need to check that the Hausdorff distance
     in the optimal coupling is smaller than the Hausdorff distance of any coupling.
     First, we check this for couplings which already have small Hausdorff distance: in this
     case, the induced "distance" on α ⊕ β belongs to the candidates family introduced in the
     definition of the optimal coupling, and the conclusion follows from the optimality
     of the optimal coupling within this family.
  -/
  have A : ∀p q : nonempty_compacts (l_infty_ℝ), ⟦p⟧ = to_GH_space α → ⟦q⟧ = to_GH_space β →
        Hausdorff_dist (p.val) (q.val) < diam (univ : set α) + 1 + diam (univ : set β) →
        Hausdorff_dist (range (optimal_GH_injl α β)) (range (optimal_GH_injr α β)) ≤ Hausdorff_dist (p.val) (q.val),
  { assume p q hp hq bound,
    rcases eq_to_GH_space_iff.1 hp with ⟨Φ, ⟨Φisom, Φrange⟩⟩,
    rcases eq_to_GH_space_iff.1 hq with ⟨Ψ, ⟨Ψisom, Ψrange⟩⟩,
    have I : diam (range Φ ∪ range Ψ) ≤ 2 * diam (univ : set α) + 1 + 2 * diam (univ : set β),
    { rcases exists_mem_of_nonempty α with ⟨xα, _⟩,
      have : ∃y ∈ range Ψ, dist (Φ xα) y < diam (univ : set α) + 1 + diam (univ : set β),
      { rw Ψrange,
        have : Φ xα ∈ p.val := begin rw ← Φrange, exact mem_range_self _ end,
        exact exists_dist_lt_of_Hausdorff_dist_lt this bound
          (Hausdorff_edist_ne_top_of_bounded_of_ne_empty p.2.1 q.2.1 (bounded_of_compact p.2.2) (bounded_of_compact q.2.2)) },
      rcases this with ⟨y, hy, dy⟩,
      rcases mem_range.1 hy with ⟨z, hzy⟩,
      rw ← hzy at dy,
      have DΦ : diam (range Φ) = diam (univ : set α) :=
        begin rw [← image_univ], apply diam_image_of_isometry Φisom end,
      have DΨ : diam (range Ψ) = diam (univ : set β) :=
        begin rw [← image_univ], apply diam_image_of_isometry Ψisom end,
      calc
        diam (range Φ ∪ range Ψ) ≤ diam (range Φ) + dist (Φ xα) (Ψ z) + diam (range Ψ) :
          diam_union (mem_range_self _) (mem_range_self _)
        ... ≤ diam (univ : set α) + (diam (univ : set α) + 1 + diam (univ : set β)) + diam (univ : set β) :
          by rw [DΦ, DΨ]; apply add_le_add (add_le_add (le_refl _) (le_of_lt dy)) (le_refl _)
        ... = 2 * diam (univ : set α) + 1 + 2 * diam (univ : set β) : by ring },

    let f : α ⊕ β → l_infty_ℝ := λx, match x with | inl y := Φ y | inr z := Ψ z end,
    let F : (α ⊕ β) × (α ⊕ β) → ℝ := λp, dist (f p.1) (f p.2),
    -- check that the induced "distance" is a candidate
    have Fgood : F ∈ candidates α β,
    { simp only [candidates, forall_const, and_true, add_comm, eq_self_iff_true, dist_eq_zero,
                 and_self, set.mem_set_of_eq],
      repeat {split},
      { exact λx y, calc
        F (inl x, inl y) = dist (Φ x) (Φ y) : rfl
        ... = dist x y : Φisom (mem_univ _) (mem_univ _) },
      { exact λx y, calc
        F (inr x, inr y) = dist (Ψ x) (Ψ y) : rfl
        ... = dist x y : Ψisom (mem_univ _) (mem_univ _) },
      { exact λx y, dist_comm _ _ },
      { exact λx y z, dist_triangle _ _ _ },
      { exact λx y, calc
        F (x, y) ≤ diam (range Φ ∪ range Ψ) :
        begin
          have A : ∀z : α ⊕ β, f z ∈ range Φ ∪ range Ψ,
          { assume z,
            cases z,
            { apply mem_union_left, apply mem_range_self },
            { apply mem_union_right, apply mem_range_self }},
          refine dist_le_diam_of_mem _ (A _) (A _),
          rw [Φrange, Ψrange],
          exact bounded_of_compact (compact_union_of_compact p.2.2 q.2.2),
        end
        ... ≤ 2 * diam (univ : set α) + 1 + 2 * diam (univ : set β) : I }},
    let Fb := candidates_b_of_candidates F Fgood,
    have : Hausdorff_dist (range (optimal_GH_injl α β)) (range (optimal_GH_injr α β)) ≤ HD Fb :=
      Hausdorff_dist_optimal_le_HD _ _ (candidates_b_of_candidates_mem F Fgood),
    refine le_trans this (le_of_forall_le_of_dense (λr hr, _)),
    have I1 : ∀x : α, infi (λy:β, Fb (inl x, inr y)) ≤ r,
    { assume x,
      have : f (inl x) ∈ p.val := by rw [← Φrange]; apply mem_range_self,
      rcases exists_dist_lt_of_Hausdorff_dist_lt this hr
        (Hausdorff_edist_ne_top_of_bounded_of_ne_empty p.2.1 q.2.1 (bounded_of_compact p.2.2) (bounded_of_compact q.2.2))
        with ⟨z, zq, hz⟩,
      have : z ∈ range Ψ := by rwa [← Ψrange] at zq,
      rcases mem_range.1 this with ⟨y, hy⟩,
      calc infi (λy:β, Fb (inl x, inr y)) ≤ Fb (inl x, inr y) :
          cinfi_le (by simpa using HD_below_aux1 0)
        ... = dist (Φ x) (Ψ y) : rfl
        ... = dist (f (inl x)) z : by rw hy
        ... ≤ r : le_of_lt hz },
    have I2 : ∀y : β, infi (λx:α, Fb (inl x, inr y)) ≤ r,
    { assume y,
      have : f (inr y) ∈ q.val := by rw [← Ψrange]; apply mem_range_self,
      rcases exists_dist_lt_of_Hausdorff_dist_lt' this hr
        (Hausdorff_edist_ne_top_of_bounded_of_ne_empty p.2.1 q.2.1 (bounded_of_compact p.2.2) (bounded_of_compact q.2.2))
        with ⟨z, zq, hz⟩,
      have : z ∈ range Φ := by rwa [← Φrange] at zq,
      rcases mem_range.1 this with ⟨x, hx⟩,
      calc infi (λx:α, Fb (inl x, inr y)) ≤ Fb (inl x, inr y) :
          cinfi_le (by simpa using HD_below_aux2 0)
        ... = dist (Φ x) (Ψ y) : rfl
        ... = dist z (f (inr y)) : by rw hx
        ... ≤ r : le_of_lt hz },
    simp [HD, csupr_le I1, csupr_le I2] },

  /- Get the same inequality for any coupling. If the coupling is quite good, the desired
  inequality has been proved above. If it is bad, then the inequality is obvious. -/
  have B : ∀p q : nonempty_compacts (l_infty_ℝ), ⟦p⟧ = to_GH_space α → ⟦q⟧ = to_GH_space β →
        Hausdorff_dist (range (optimal_GH_injl α β)) (range (optimal_GH_injr α β)) ≤ Hausdorff_dist (p.val) (q.val),
  { assume p q hp hq,
    by_cases h : Hausdorff_dist (p.val) (q.val) < diam (univ : set α) + 1 + diam (univ : set β),
    { exact A p q hp hq h },
    { calc Hausdorff_dist (range (optimal_GH_injl α β)) (range (optimal_GH_injr α β)) ≤ HD (candidates_b_dist α β) :
             Hausdorff_dist_optimal_le_HD _ _ (candidates_b_dist_mem_candidates_b)
           ... ≤ diam (univ : set α) + 1 + diam (univ : set β) : HD_candidates_b_dist_le
           ... ≤ Hausdorff_dist (p.val) (q.val) : not_lt.1 h }},
  refine le_antisymm _ _,
  { apply le_cInf,
    { rw ne_empty_iff_exists_mem,
      simp only [set.mem_image, nonempty_of_inhabited, set.mem_set_of_eq, prod.exists],
      existsi [Hausdorff_dist (nonempty_compacts.Kuratowski_embedding α).val (nonempty_compacts.Kuratowski_embedding β).val,
               nonempty_compacts.Kuratowski_embedding α, nonempty_compacts.Kuratowski_embedding β],
      simp [to_GH_space, -quotient.eq] },
    { rintro b ⟨⟨p, q⟩, ⟨hp, hq⟩, rfl⟩,
      exact B p q hp hq }},
  { exact GH_dist_le_Hausdorff_dist (isometry_optimal_GH_injl α β) (isometry_optimal_GH_injr α β) }
end

/-- The Gromov-Hausdorff distance can also be realized by a coupling in ℓ^∞(ℝ), by embedding
the optimal coupling through its Kuratowski embedding. -/
theorem GH_dist_eq_Hausdorff_dist (α : Type u) [metric_space α] [compact_space α] [nonempty α]
  (β : Type v) [metric_space β] [compact_space β] [nonempty β] :
  ∃Φ : α → l_infty_ℝ, ∃Ψ : β → l_infty_ℝ, isometry_on Φ univ ∧ isometry_on Ψ univ ∧
  GH_dist α β = Hausdorff_dist (range Φ) (range Ψ) :=
begin
  let F := Kuratowski_embedding (optimal_GH_coupling α β),
  let Φ := F ∘ optimal_GH_injl α β,
  let Ψ := F ∘ optimal_GH_injr α β,
  refine ⟨Φ, Ψ, _, _, _⟩,
  { exact isometry_on_comp (maps_to_univ _ _) Kuratowski_embedding_isometry (isometry_optimal_GH_injl α β) },
  { exact isometry_on_comp (maps_to_univ _ _) Kuratowski_embedding_isometry (isometry_optimal_GH_injr α β) },
  { rw [← image_univ, ← image_univ, image_comp F, image_univ, image_comp F (optimal_GH_injr α β),
      image_univ, ← Hausdorff_dist_optimal],
    exact (Hausdorff_dist.image (subset_univ _) (subset_univ _) Kuratowski_embedding_isometry).symm },
end

-- without the next two lines, { exact closed_of_compact (range Φ) hΦ } in the next
-- proof is very slow, as the t2_space instance is very hard to find
local attribute [instance, priority 0] orderable_topology.t2_space
local attribute [instance, priority 0] ordered_topology.to_t2_space

/-- The Gromov-Hausdorff distance defines a genuine distance on the Gromov-Hausdorff space. -/
instance GH_space_metric_space : metric_space GH_space :=
{ dist_self := λx, begin
    rcases exists_rep x with ⟨y, hy⟩,
    refine le_antisymm _ _,
    { apply cInf_le,
      { exact ⟨0, by rintro b ⟨⟨u, v⟩, ⟨hu, hv⟩, rfl⟩; exact Hausdorff_dist.nonneg ⟩},
      { simp, existsi [y, y], simpa }},
    { apply le_cInf,
      { simp only [set.image_eq_empty, ne.def],
        apply ne_empty_iff_exists_mem.2,
        existsi (⟨y, y⟩ : nonempty_compacts l_infty_ℝ × nonempty_compacts l_infty_ℝ),
        simpa },
      { rintro b ⟨⟨u, v⟩, ⟨hu, hv⟩, rfl⟩, exact Hausdorff_dist.nonneg }},
  end,
  dist_comm := λx y, begin
    have A : (λ (p : nonempty_compacts l_infty_ℝ × nonempty_compacts l_infty_ℝ),
                 Hausdorff_dist ((p.fst).val) ((p.snd).val)) '' (set.prod {a | ⟦a⟧ = x} {b | ⟦b⟧ = y})
           = ((λ (p : nonempty_compacts l_infty_ℝ × nonempty_compacts l_infty_ℝ),
                 Hausdorff_dist ((p.fst).val) ((p.snd).val)) ∘ prod.swap) '' (set.prod {a | ⟦a⟧ = x} {b | ⟦b⟧ = y}) :=
      by congr; funext; simp; rw Hausdorff_dist.comm,
    simp only [dist, A, image_comp, prod.swap, image_swap_prod],
  end,
  eq_of_dist_eq_zero := λx y hxy, begin
    /- To show that two spaces at zero distance are isometric, we argue that the distance
    is realized by some coupling. In this coupling, the two spaces are at zero Hausdorff distance,
    i.e., they coincide. Therefore, the original spaces are isometric. -/
    rcases GH_dist_eq_Hausdorff_dist (rep_GH_space x) (rep_GH_space y) with ⟨Φ, Ψ, Φisom, Ψisom, DΦΨ⟩,
    rw [← dist_GH_dist, hxy] at DΦΨ,
    have : range Φ = range Ψ,
    { have hΦ : compact (range Φ) :=
        by rw [← image_univ]; exact compact_image compact_univ (continuous_of_isometry Φisom),
      have hΨ : compact (range Ψ) :=
        by rw [← image_univ]; exact compact_image compact_univ (continuous_of_isometry Ψisom),
      apply (Hausdorff_dist.zero_iff_eq_of_closed _ _ _).1 (DΦΨ.symm),
      { exact closed_of_compact (range Φ) hΦ },
      { exact closed_of_compact (range Ψ) hΨ },
      { exact Hausdorff_edist_ne_top_of_bounded_of_ne_empty (by simp [-nonempty_subtype])
          (by simp [-nonempty_subtype]) (bounded_of_compact hΦ) (bounded_of_compact hΨ) }},
    let F := (inv_fun_on Ψ univ) ∘ Φ,
    rw [← to_GH_space_rep_GH_space x, ← to_GH_space_rep_GH_space y],
    refine to_GH_space_eq_to_GH_space_iff_isometric.2 ⟨F, ⟨_, _⟩⟩,
    { apply isometry_on_comp (maps_to_range _ univ) _ Φisom,
      rw [‹range Φ = range Ψ›, ← image_univ],
      exact isometry_on_inv_fun_on Ψisom },
    { rw [← image_univ, image_comp, image_univ, ‹range Φ = range Ψ›, ← image_univ],
      exact inv_fun_on_image (inj_on_of_isometry_on Ψisom) (subset_univ _) },
  end,
  dist_triangle := λx y z, begin
    /- To show the triangular inequality between X, Y and Z, realize an optimal coupling
    between X and Y in a space γ1, and an optimal coupling between Y and Z in a space γ2. Then,
    glue these metric spaces along Y. We get a new space γ in which X and Y are optimally coupled,
    as well as Y and Z. Apply the triangle inequality for the Hausdorff distance in γ to conclude. -/
    let X := rep_GH_space x,
    let Y := rep_GH_space y,
    let Z := rep_GH_space z,
    let γ1 := optimal_GH_coupling X Y,
    let γ2 := optimal_GH_coupling Y Z,
    let Φ : Y → γ1 := optimal_GH_injr X Y,
    have hΦ : isometry_on Φ univ := isometry_optimal_GH_injr X Y,
    let Ψ : Y → γ2 := optimal_GH_injl Y Z,
    have hΨ : isometry_on Ψ univ := isometry_optimal_GH_injl Y Z,
    let γ := glue_space hΦ hΨ,
    letI : metric_space γ := metric.metric_space_glue_space hΦ hΨ,
    have Comm : (to_glue_l hΦ hΨ) ∘ (optimal_GH_injr X Y) = (to_glue_r hΦ hΨ) ∘ (optimal_GH_injl Y Z) :=
      to_glue_commute hΦ hΨ,
    calc dist x z = dist (to_GH_space X) (to_GH_space Z) :
        by rw [to_GH_space_rep_GH_space x, to_GH_space_rep_GH_space z]
      ... ≤ Hausdorff_dist (range ((to_glue_l hΦ hΨ) ∘ (optimal_GH_injl X Y)))
                       (range ((to_glue_r hΦ hΨ) ∘ (optimal_GH_injr Y Z))) :
        GH_dist_le_Hausdorff_dist
          (isometry_on_comp (maps_to_univ _ _) (to_glue_l_isometry hΦ hΨ) (isometry_optimal_GH_injl X Y))
          (isometry_on_comp (maps_to_univ _ _) (to_glue_r_isometry hΦ hΨ) (isometry_optimal_GH_injr Y Z))
      ... ≤ Hausdorff_dist (range ((to_glue_l hΦ hΨ) ∘ (optimal_GH_injl X Y)))
                           (range ((to_glue_l hΦ hΨ) ∘ (optimal_GH_injr X Y)))
          + Hausdorff_dist (range ((to_glue_l hΦ hΨ) ∘ (optimal_GH_injr X Y)))
                           (range ((to_glue_r hΦ hΨ) ∘ (optimal_GH_injr Y Z))) :
        begin
          refine Hausdorff_dist.triangle (Hausdorff_edist_ne_top_of_bounded_of_ne_empty
            (by simp [-nonempty_subtype]) (by simp [-nonempty_subtype]) _ _),
          { rw [← image_univ],
            exact bounded_of_compact (compact_image compact_univ (continuous_of_isometry
              (isometry_on_comp (maps_to_univ _ _) (to_glue_l_isometry hΦ hΨ) (isometry_optimal_GH_injl X Y)))) },
          { rw [← image_univ],
            exact bounded_of_compact (compact_image compact_univ (continuous_of_isometry
              (isometry_on_comp (maps_to_univ _ _) (to_glue_l_isometry hΦ hΨ) (isometry_optimal_GH_injr X Y)))) }
        end
      ... = Hausdorff_dist ((to_glue_l hΦ hΨ) '' (range (optimal_GH_injl X Y)))
                           ((to_glue_l hΦ hΨ) '' (range (optimal_GH_injr X Y)))
          + Hausdorff_dist ((to_glue_r hΦ hΨ) '' (range (optimal_GH_injl Y Z)))
                           ((to_glue_r hΦ hΨ) '' (range (optimal_GH_injr Y Z))) :
        by simp only [eq.symm range_comp, Comm, eq_self_iff_true, add_right_inj]
      ... = Hausdorff_dist (range (optimal_GH_injl X Y))
                           (range (optimal_GH_injr X Y))
          + Hausdorff_dist (range (optimal_GH_injl Y Z))
                           (range (optimal_GH_injr Y Z)) :
        by rw [Hausdorff_dist.image (subset_univ _) (subset_univ _) (to_glue_l_isometry hΦ hΨ),
               Hausdorff_dist.image (subset_univ _) (subset_univ _) (to_glue_r_isometry hΦ hΨ)]
      ... = dist (to_GH_space X) (to_GH_space Y) + dist (to_GH_space Y) (to_GH_space Z) :
        by rw [Hausdorff_dist_optimal, Hausdorff_dist_optimal, GH_dist, GH_dist]
      ... = dist x y + dist y z:
        by rw [to_GH_space_rep_GH_space x, to_GH_space_rep_GH_space y, to_GH_space_rep_GH_space z]
  end }

end GH_space --section

section
/- In this section, we show that if two metric spaces are isometric up to ε2, then their Gromov-Hausdorff distance
is bounded by ε2 / 2. More generally, if there are subsets which are ε1-dense and ε3-dense in two spaces, and
isometric up to ε2, then the Gromov-Hausdorff distance between the spaces is bounded by ε1 + ε2/2 + ε3.
For this, we construct a suitable coupling between the two spaces, by gluing them (approximately) along
the two matching subsets. -/

variables {α : Type u} [metric_space α] [compact_space α] [nonempty α]
          {β : Type v} [metric_space β] [compact_space β] [nonempty β]

/-- If there are subsets which are ε1-dense and ε3-dense in two spaces, and
isometric up to ε2, then the Gromov-Hausdorff distance between the spaces is bounded by ε1 + ε2/2 + ε3. -/
theorem GH_dist_le_of_approx_subsets {s : set α} (Φ : s → β) {ε1 ε2 ε3 : ℝ}
  (hs : ∀x : α, ∃y ∈ s, dist x y ≤ ε1) (hs' : ∀x : β, ∃y : s, dist x (Φ y) ≤ ε3)
  (H : ∀x y : s, abs (dist x y - dist (Φ x) (Φ y)) ≤ ε2) :
  GH_dist α β ≤ ε1 + ε2 / 2 + ε3 :=
begin
  refine real.le_of_forall_epsilon_le (λδ δ0, _),
  rcases exists_mem_of_nonempty α with ⟨xα, _⟩,
  rcases hs xα with ⟨xs, hxs, Dxs⟩,
  haveI : nonempty s := ⟨⟨xs, hxs⟩⟩,
  have sne : s ≠ ∅ := ne_empty_of_mem hxs,
  letI : nonempty (subtype s) := ⟨⟨xs, hxs⟩⟩,
  have : 0 ≤ ε2 := le_trans (abs_nonneg _) (H ⟨xs, hxs⟩ ⟨xs, hxs⟩),
  have : ∀ p q : s, abs (dist p q - dist (Φ p) (Φ q)) ≤ 2 * (ε2/2 + δ) := λp q, calc
    abs (dist p q - dist (Φ p) (Φ q)) ≤ ε2 : H p q
    ... ≤ 2 * (ε2/2 + δ) : by linarith,
    -- glue α and β along the almost matching subsets
  letI : metric_space (α ⊕ β) := glue_metric_approx (@subtype.val α s) (λx, Φ x) (ε2/2 + δ) (by linarith) this,
  let Fl := @sum.inl α β,
  let Fr := @sum.inr α β,
  have Il : isometry_on Fl univ := λx y hx hy, rfl,
  have Ir : isometry_on Fr univ := λx y hx hy, rfl,
  /- The proof goes as follows : the GH_dist is bounded by the Hausdorff distance of the images in the
  coupling, which is bounded (using the triangular inequality) by the sum of the Hausdorff distances
  of α and s (in the coupling or, equivalently in the original space), of s and Φ s, and of Φ s and β
  (in the coupling or, equivalently, in the original space). The first term is bounded by ε1, by density.
  The third one is bounded by ε3. And the middle one is bounded by ε2/2 as in the coupling the points x and Φ x
  are at distance ε2/2 by construction of the coupling (in fact ε2/2 + δ where δ is an arbitrarily small positive
  constant where positivity is used to ensure that the coupling is really a metric space and not a premetric
  space on α ⊕ β. -/
  have : GH_dist α β ≤ Hausdorff_dist (range Fl) (range Fr) :=
    GH_dist_le_Hausdorff_dist Il Ir,
  have : Hausdorff_dist (range Fl) (range Fr) ≤ Hausdorff_dist (range Fl) (Fl '' s)
                                              + Hausdorff_dist (Fl '' s) (range Fr),
  { have B : bounded (range Fl) := bounded_of_compact (compact_range (continuous_of_isometry Il)),
    exact Hausdorff_dist.triangle (Hausdorff_edist_ne_top_of_bounded_of_ne_empty (by simpa) (by simpa)
      B (bounded.subset (image_subset_range _ _) B)) },
  have : Hausdorff_dist (Fl '' s) (range Fr) ≤ Hausdorff_dist (Fl '' s) (Fr '' (range Φ))
                                             + Hausdorff_dist (Fr '' (range Φ)) (range Fr),
  { have B : bounded (range Fr) := bounded_of_compact (compact_range (continuous_of_isometry Ir)),
    exact Hausdorff_dist.triangle' (Hausdorff_edist_ne_top_of_bounded_of_ne_empty
      (by simpa [-nonempty_subtype]) (by simpa) (bounded.subset (image_subset_range _ _) B) B) },
  have : Hausdorff_dist (range Fl) (Fl '' s) ≤ ε1,
  { rw [← image_univ, Hausdorff_dist.image (subset_univ _) (subset_univ _) Il],
    have : 0 ≤ ε1 := le_trans dist_nonneg Dxs,
    refine Hausdorff_dist.le_of_mem_dist_le_of_mem_dist_le this (λx hx, hs x)
      (λx hx, ⟨x, mem_univ _, by simpa⟩) },
  have : Hausdorff_dist (Fl '' s) (Fr '' (range Φ)) ≤ ε2/2 + δ,
  { refine Hausdorff_dist.le_of_mem_dist_le_of_mem_dist_le (by linarith) _ _,
    { assume x' hx',
      rcases (set.mem_image _ _ _).1 hx' with ⟨x, ⟨x_in_s, xx'⟩⟩,
      rw ← xx',
      use [Fr (Φ ⟨x, x_in_s⟩), mem_image_of_mem Fr (mem_range_self _)],
      exact le_of_eq (glue_dist_glued_points (@subtype.val α s) Φ (ε2/2 + δ) ⟨x, x_in_s⟩) },
    { assume x' hx',
      rcases (set.mem_image _ _ _).1 hx' with ⟨y, ⟨y_in_s', yx'⟩⟩,
      rcases mem_range.1 y_in_s' with ⟨x, xy⟩,
      use [Fl x, mem_image_of_mem _ x.2],
      rw [← yx', ← xy, dist_comm],
      exact le_of_eq (glue_dist_glued_points (@subtype.val α s) Φ (ε2/2 + δ) x) }},
  have : Hausdorff_dist (Fr '' (range Φ)) (range Fr) ≤ ε3,
  { rw [← @image_univ _ _ Fr, Hausdorff_dist.image (subset_univ _) (subset_univ _) Ir],
    rcases exists_mem_of_nonempty β with ⟨xβ, _⟩,
    rcases hs' xβ with ⟨xs', Dxs'⟩,
    have : 0 ≤ ε3 := le_trans dist_nonneg Dxs',
    refine Hausdorff_dist.le_of_mem_dist_le_of_mem_dist_le this (λx hx, ⟨x, mem_univ _, by simpa⟩) (λx _, _),
    rcases hs' x with ⟨y, Dy⟩,
    exact ⟨Φ y, mem_range_self _, Dy⟩ },
  linarith
end
end

lemma dependent_helper {F G : Σn:ℕ, (fin n → fin n → ℤ)} (h : F = G) {i j: ℕ}
  {iF : i < F.1} {iG : i < G.1} {jF : j < F.1} {jG : j < G.1} :
  F.2 ⟨i, iF⟩ ⟨j, jF⟩ = G.2 ⟨i, iG⟩ ⟨j, jG⟩ :=
by subst h; refl

/-- The Gromov-Hausdorff space is second countable. -/
lemma second_countable : second_countable_topology GH_space :=
begin
  refine second_countable_of_countable_discretization (λδ δpos, _),
  let ε := (2/5) * δ,
  have εpos : 0 < ε := mul_pos (by norm_num) δpos,
  have : ∀p:GH_space, ∃s : set (rep_GH_space p), finite s ∧ (univ ⊆ (⋃x∈s, ball x ε)) :=
    λp, by simpa using finite_cover_balls_of_compact (@compact_univ (rep_GH_space p) _ _ _) εpos,
  -- for each p, s p is a finite ε-dense subset of p (or rather the metric space
  -- rep_GH_space p representing p)
  choose s hs using this,
  have : ∀p:GH_space, ∀t:set (rep_GH_space p), finite t → ∃n:ℕ, ∃e:equiv t (fin n), true,
  { assume p t ht,
    letI : fintype t := finite.fintype ht,
    rcases fintype.exists_equiv_fin t with ⟨n, hn⟩,
    rcases exists_mem_of_nonempty (↥t ≃ fin n) with ⟨e, _⟩,
    exact ⟨n, e, trivial⟩ },
  choose N e hne using this,
  -- cardinality of a nice finite subset of rep_GH_space p, called N p
  let N := λp:GH_space, N p (s p) (hs p).1,
  -- equiv from s p, a nice finite subset of rep_GH_space p, to fin (N p), called E p
  let E := λp:GH_space, e p (s p) (hs p).1,
  -- associate to p ∈ GH_space a function F recording the distances of points
  -- in the ε-dense set s p.
  let F : GH_space → Σn:ℕ, (fin n → fin n → ℤ) :=
    λp, ⟨N p, λa b, floor (ε⁻¹ * dist ((E p).inv_fun a) ((E p).inv_fun b))⟩,
  refine ⟨_, by apply_instance, F, λp q hpq, _⟩,
  /- As the target space of F is countable, it suffices to show that two points
  p and q with F p = F q are at distance ≤ δ.
  For this, we construct a map Φ from s p ⊆ rep_GH_space p (representing p)
  to rep_GH_space q (representing q) which is almost an isometry on s p, and
  with image s q. For this, we compose the identification of s p with fin (N p)
  and the inverse of the identification of s q with fin (N q). Together with
  the fact that N p = N q, this constructs Ψ between s p and s q, and then
  composing with the canonical inclusion we get Φ. -/
  have Npq : N p = N q := (sigma.mk.inj_iff.1 hpq).1,
  let Ψ : s p → s q := λx, (E q).inv_fun (fin.cast Npq ((E p).to_fun x)),
  let Φ : s p → rep_GH_space q := λx, Ψ x,
  -- Use the almost isometry Φ to show that rep_GH_space p and rep_GH_space q
  -- are within controlled Gromov-Hausdorff distance.
  have M : GH_dist (rep_GH_space p) (rep_GH_space q) ≤ ε + ε/2 + ε,
  { refine GH_dist_le_of_approx_subsets Φ  _ _ _,
    show ∀x : rep_GH_space p, ∃ (y : rep_GH_space p) (H : y ∈ s p), dist x y ≤ ε,
    { -- by construction, s p is ε-dense
      assume x,
      have : x ∈ ⋃y∈(s p), ball y ε := (hs p).2 (mem_univ _),
      rcases mem_bUnion_iff.1 this with ⟨y, ⟨ys, hy⟩⟩,
      exact ⟨y, ys, le_of_lt hy⟩ },
    show ∀x : rep_GH_space q, ∃ (z : s p), dist x (Φ z) ≤ ε,
    { -- by construction, s q is ε-dense, and it is the range of Φ
      assume x,
      have : x ∈ ⋃y∈(s q), ball y ε := (hs q).2 (mem_univ _),
      rcases mem_bUnion_iff.1 this with ⟨y, ⟨ys, hy⟩⟩,
      let i := ((E q).to_fun ⟨y, ys⟩).1,
      let hi := ((E q).to_fun ⟨y, ys⟩).2,
      have ihi_eq : (⟨i, hi⟩ : fin (N q)) = (E q).to_fun ⟨y, ys⟩ := by rw fin.ext_iff,
      have hiq : i < N q := hi,
      have hip : i < N p, { rwa Npq.symm at hiq },
      let z := (E p).inv_fun ⟨i, hip⟩,
      use z,
      have C1 : (E p).to_fun z = ⟨i, hip⟩ := (E p).right_inv ⟨i, hip⟩,
      have C2 : fin.cast Npq ⟨i, hip⟩ = ⟨i, hi⟩ := rfl,
      have C3 : (E q).inv_fun ⟨i, hi⟩ = ⟨y, ys⟩ := by rw ihi_eq; exact (E q).left_inv ⟨y, ys⟩,
      have : Φ z = y :=
        by simp only [Φ, Ψ]; rw [C1, C2, C3]; refl,
      rw this,
      exact le_of_lt hy },
    show ∀x y : s p, abs (dist x y - dist (Φ x) (Φ y)) ≤ ε,
    { /- the distance between x and y is encoded in F p, and the distance between
      Φ x and Φ y (two points of s q) is encoded in F q, all this up to ε.
      As F p = F q, the distances are almost equal. -/
      assume x y,
      have : dist (Φ x) (Φ y) = dist (Ψ x) (Ψ y) := rfl,
      rw this,
      -- introduce i, that codes both x and Φ x in fin (N p) = fin (N q)
      let i := ((E p).to_fun x).1,
      have hip : i < N p := ((E p).to_fun x).2,
      have hiq : i < N q, by rwa Npq at hip,
      have i' : i = ((E q).to_fun (Ψ x)).1, by simp [Ψ, (E q).right_inv _]; refl,
      -- introduce j, that codes both y and Φ y in fin (N p) = fin (N q)
      let j := ((E p).to_fun y).1,
      have hjp : j < N p := ((E p).to_fun y).2,
      have hjq : j < N q, by rwa Npq at hjp,
      have j' : j = ((E q).to_fun (Ψ y)).1, by simp [Ψ, (E q).right_inv _]; refl,
      -- Express dist x y in terms of F p
      have : (F p).2 ((E p).to_fun x) ((E p).to_fun y) = floor (ε⁻¹ * dist x y),
        by simp only [F, (E p).left_inv _],
      have Ap : (F p).2 ⟨i, hip⟩ ⟨j, hjp⟩ = floor (ε⁻¹ * dist x y),
        by rw ← this; congr; apply (fin.ext_iff _ _).2; refl,
      -- Express dist (Φ x) (Φ y) in terms of F q
      have : (F q).2 ((E q).to_fun (Ψ x)) ((E q).to_fun (Ψ y)) = floor (ε⁻¹ * dist (Ψ x) (Ψ y)),
        by simp only [F, (E q).left_inv _],
      have Aq : (F q).2 ⟨i, hiq⟩ ⟨j, hjq⟩ = floor (ε⁻¹ * dist (Ψ x) (Ψ y)),
        by rw ← this; congr; apply (fin.ext_iff _ _).2; [exact i', exact j'],
      -- use the equality between F p and F q to deduce that the distances have equal
      -- integer parts
      have : (F p).2 ⟨i, hip⟩ ⟨j, hjp⟩ = (F q).2 ⟨i, hiq⟩ ⟨j, hjq⟩ :=
        dependent_helper hpq,
      rw [Ap, Aq] at this,
      -- deduce that the distances coincide up to ε, by a straightforward computation
      -- that should be automated
      have I := calc
        abs (ε⁻¹) * abs (dist x y - dist (Ψ x) (Ψ y)) =
          abs (ε⁻¹ * (dist x y - dist (Ψ x) (Ψ y))) : (abs_mul _ _).symm
        ... = abs ((ε⁻¹ * dist x y) - (ε⁻¹ * dist (Ψ x) (Ψ y))) : by congr; ring
        ... ≤ 1 : le_of_lt (real.abs_diff_lt_one_of_floor_eq_floor this),
       calc
         abs (dist x y - dist (Ψ x) (Ψ y)) = (ε * ε⁻¹) * abs (dist x y - dist (Ψ x) (Ψ y)) :
           by rw [mul_inv_cancel (ne_of_gt εpos), one_mul]
        ... = ε * (abs (ε⁻¹) * abs (dist x y - dist (Ψ x) (Ψ y))) :
          by rw [abs_of_nonneg (le_of_lt (inv_pos εpos)), mul_assoc]
        ... ≤ ε * 1 : mul_le_mul_of_nonneg_left I (le_of_lt εpos)
        ... = ε : mul_one _ }},
  calc dist p q = GH_dist (rep_GH_space p) (rep_GH_space q) : dist_GH_dist p q
    ... ≤ ε + ε/2 + ε : M
    ... = δ : by simp [ε]; ring
end

end Gromov_Hausdorff --namespace