Theory Quick_Sort_Average_Case

(*
  File:     Quick_Sort_Average_Case.thy
  Author:   Manuel Eberl <manuel@pruvisto.org>

  Definition and average-case analysis of the standard deterministic QuickSort algorithm
*)
section ‹Average case analysis of deterministic QuickSort›
theory Quick_Sort_Average_Case
  imports Randomised_Quick_Sort
begin
  
subsection ‹Definition of deterministic QuickSort›
  
text ‹
  This is the functional description of the standard variant of deterministic QuickSort that 
  always chooses the first list element as the pivot as given by Hoare in 1962~cite"hoare". 
  For a list that is already sorted, this leads to $n(n-1)$ 
  comparisons, but as is well known, the average case is not that bad.
›
fun quicksort :: "('a × 'a) set  'a list  'a list" where
  "quicksort _ [] = []"
| "quicksort R (x # xs) = 
     quicksort R (filter (λy. (y,x)  R) xs) @ [x] @ quicksort R (filter (λy. (y,x)  R) xs)"

text ‹
  We can easily show that this QuickSort is correct:
›
theorem mset_quicksort [simp]: "mset (quicksort R xs) = mset xs"
  by (induction R xs rule: quicksort.induct) (simp_all)

corollary set_quicksort [simp]: "set (quicksort R xs) = set xs"
  by (induction R xs rule: quicksort.induct) auto

theorem sorted_wrt_quicksort: 
  assumes "trans R" and "total_on (set xs) R" and "x. x  set xs  (x, x)  R"
  shows   "sorted_wrt R (quicksort R xs)"
using assms
proof (induction R xs rule: quicksort.induct)
  case (2 R x xs)
  have total: "(a, b)  R" if "(b, a)  R" "a  set (x#xs)" "b  set (x#xs)" for a b
    using "2.prems" that unfolding total_on_def by (cases "a = b") auto
    
  have *: "sorted_wrt R (quicksort R (filter (λy. (y,x)  R) xs))"
          "sorted_wrt R (quicksort R (filter (λy. (y,x)  R) xs))"
    by ((rule 2 total_on_subset[OF total_on (set (x#xs)) R]) | force)+
  show ?case
    by (auto intro!: sorted_wrt_append sorted_wrt.intros trans R * 
             intro: transD[OF trans R] dest!: total simp: total_on_def)
qed auto

corollary sorted_wrt_quicksort':
  assumes "linorder_on A R" and "set xs  A"
  shows   "sorted_wrt R (quicksort R xs)"
  by (rule sorted_wrt_quicksort)
     (insert assms, auto simp: linorder_on_def refl_on_def dest: total_on_subset)

text ‹
  We now define another version of QuickSort that is identical to the previous one but also 
  counts the number of comparisons that were made.
›
fun quicksort' :: "('a × 'a) set  'a list  'a list × nat" where
  "quicksort' _ [] = ([], 0)"
| "quicksort' R (x # xs) = (
     let (ls, rs)  = partition (λy. (y,x)  R) xs;
         (ls', n1) = quicksort' R ls;
         (rs', n2) = quicksort' R rs
     in
         (ls' @ [x] @ rs', length xs + n1 + n2))"

text ‹
  For convenience, we also define a function that computes only the number of comparisons that 
  were made and not the result list.
›
fun qs_cost :: "('a × 'a) set  'a list  nat" where
  "qs_cost _ [] = 0"
| "qs_cost R (x # xs) = 
     length xs + qs_cost R (filter (λy. (y,x)R) xs) + qs_cost R (filter (λy. (y,x)R) xs)"


text ‹
  It is obvious that the original QuickSort and the cost function are the projections 
  of the cost-counting QuickSort.
›  
lemma fst_quicksort' [simp]: "fst (quicksort' R xs) = quicksort R xs"
  by (induction R xs rule: quicksort.induct) (simp_all add: case_prod_unfold Let_def o_def)

lemma snd_quicksort' [simp]: "snd (quicksort' R xs) = qs_cost R xs"
  by (induction R xs rule: quicksort.induct) (simp_all add: case_prod_unfold Let_def o_def)

    
subsection ‹Analysis›

text ‹
  We will reduce the average-case analysis to showing that it is essentially equivalent to 
  the randomised QuickSort we analysed earlier. Similar, but more direct analyses are given 
  by Hoare~cite"hoare" and Sedgewick~cite"sedgewick". 

  The proof is relatively straightforward -- but still a bit messy. We show that the cost 
  distribution of QuickSort run on a random permutation of a set of size $n$ is exactly the same 
  as that of randomised QuickSort being run on any fixed list of size $n$ (which we analysed 
  before):  
›
theorem qs_cost_average_conv_rqs_cost:
  assumes "finite A" and "linorder_on B R" and "A  B"
  shows   "map_pmf (qs_cost R) (pmf_of_set (permutations_of_set A)) = rqs_cost (card A)"
using assms(1,3)
proof (induction A rule: finite_psubset_induct)
  case (psubset A)
  show ?case
  proof (cases "A = {}")
    case True
    thus ?thesis by (simp add: pmf_of_set_singleton)
  next
    case False
    note A = finite A A  {}
    define n where "n = card A - 1"
    from A have "pmf_of_set (permutations_of_set A) = 
      do {x  pmf_of_set A; xs  pmf_of_set (permutations_of_set (A - {x})); return_pmf (x#xs)}"
      by (rule random_permutation_of_set)
    also have "map_pmf (qs_cost R)  =
                 do {
                   x  pmf_of_set A;
                   xs  pmf_of_set (permutations_of_set (A - {x}));
                   return_pmf (length xs + qs_cost R [yxs. (y,x)R] + qs_cost R [yxs. (y,x)R])
                 }" by (simp add: map_bind_pmf)
    also have " = map_pmf (λm. n + m) (
          do {
            x  pmf_of_set A;
            xs  pmf_of_set (permutations_of_set (A - {x}));
            return_pmf (qs_cost R [yxs. (y,x)R] + qs_cost R [yxs. (y,x)R])
          })" (is "_ = map_pmf _ ?X") using A unfolding n_def map_bind_pmf
      by (intro bind_pmf_cong map_pmf_cong refl) (auto simp: length_finite_permutations_of_set)
    also have "?X = do {
                      x  pmf_of_set A;
                      (ls,rs)  map_pmf (partition (λy. (y,x)R)) 
                                   (pmf_of_set (permutations_of_set (A - {x})));
                      return_pmf (qs_cost R ls + qs_cost R rs)
                    }" by (simp add: bind_map_pmf o_def)
    also have " = do {
                      x  pmf_of_set A;
                      (n1, n2)  pair_pmf 
                        (rqs_cost (linorder_rank R A x)) (rqs_cost (n - linorder_rank R A x));
                      return_pmf (n1 + n2)}"
    proof (intro bind_pmf_cong refl, goal_cases)
      case (1 x)
      have "map_pmf (partition (λy. (y,x)R)) (pmf_of_set (permutations_of_set (A - {x})))
               (λ(ls, rs). return_pmf (qs_cost R ls + qs_cost R rs)) = 
            map_pmf (λ(n1, n2). n1 + n2) (pair_pmf
              (map_pmf (qs_cost R) (pmf_of_set (permutations_of_set {xa  A - {x}. (xa, x)  R})))
              (map_pmf (qs_cost R) (pmf_of_set (permutations_of_set {xa  A - {x}. (xa, x)  R}))))"
        (is "_ = map_pmf _ (pair_pmf ?X ?Y)")
        by (subst partition_random_permutations)
           (simp_all add: map_pmf_def case_prod_unfold bind_return_pmf bind_assoc_pmf pair_pmf_def A)
      also {
        have "{xa  A - {x}. (xa, x)  R}  A - {x}" by blast
        also have "  A" using 1 A by auto
        finally have subset: "{xa  A - {x}. (xa, x)  R}  A" .
        also have "  B" by fact
        finally have "?X = rqs_cost (card {xa  A - {x}. (xa, x)  R})" using subset
          by (intro psubset.IH) auto
        also have "card {xa  A - {x}. (xa, x)  R} = linorder_rank R A x"
          by (simp add: linorder_rank_def)
        finally have "?X = rqs_cost " .
      }
      also {
        have "{xa  A - {x}. (xa, x)  R}  A - {x}" by blast
        also have "  A" using 1 A by auto
        finally have subset: "{xa  A - {x}. (xa, x)  R}  A" .
        also have "  B" by fact
        finally have "?Y = rqs_cost (card {xa  A - {x}. (xa, x)  R})" using subset
          by (intro psubset.IH) auto
        also {
          have "card ({yA-{x}. (y,x)R}  {yA-{x}. (y,x)R}) = 
                  linorder_rank R A x + card {xa  A - {x}. (xa, x)  R}"
            unfolding linorder_rank_def using A by (intro card_Un_disjoint) auto
          also have "{yA-{x}. (y,x)R}  {yA-{x}. (y,x)R} = A - {x}" by blast
          also have "card  = n" using A 1 by (simp add: n_def)
          finally have "card {xa  A - {x}. (xa, x)  R} = n - linorder_rank R A x" by simp
        }
        finally have "?Y = rqs_cost (n - linorder_rank R A x)" .
      }
      finally show ?case by (simp add: case_prod_unfold map_pmf_def)
    qed
    also have " = do {
                      i  map_pmf (linorder_rank R A) (pmf_of_set A);
                      (n1, n2)  pair_pmf (rqs_cost i) (rqs_cost (n - i));
                      return_pmf (n1 + n2)
                    }" by (simp add: bind_map_pmf)
    also have "map_pmf (linorder_rank R A) (pmf_of_set A) = pmf_of_set {..<card A}"
      by (intro map_pmf_of_set_bij_betw bij_betw_linorder_rank[OF assms(2)] A psubset.prems)
    also from A have "card A > 0" by (intro Nat.gr0I) auto
    hence "{..<card A} = {..n}" by (auto simp: n_def)
    also have "map_pmf (λm. n + m) (
                 do {
                      i  pmf_of_set {..n};
                      (n1, n2)  pair_pmf (rqs_cost i) (rqs_cost (n - i));
                      return_pmf (n1 + n2)
                    }) = rqs_cost (Suc n)"
      by (simp add: pair_pmf_def map_bind_pmf case_prod_unfold
                    bind_assoc_pmf bind_return_pmf add_ac)
    also from A have "card A > 0" by (intro Nat.gr0I) auto
    hence "Suc n = card A" by (simp add: n_def)
    finally show ?thesis .
  qed
qed

text ‹
  We therefore have the same expectation as well. (Note that we showed 
  @{thm rqs_cost_exp_eq [no_vars]} and @{thm rqs_cost_exp_asymp_equiv [no_vars]} before.
›
corollary expectation_qs_cost: 
  assumes "finite A" and "linorder_on B R" and "A  B"
  defines "random_list  pmf_of_set (permutations_of_set A)"
  shows   "measure_pmf.expectation (map_pmf (qs_cost R) random_list) real = 
             rqs_cost_exp (card A)"
  unfolding random_list_def
  by (subst qs_cost_average_conv_rqs_cost[OF assms(1-3)]) (simp add: expectation_rqs_cost)

end