Theory EdmondsKarp_Impl

section ‹Implementation of the Edmonds-Karp Algorithm›
theory EdmondsKarp_Impl
imports 
  EdmondsKarp_Algo
  Augmenting_Path_BFS
  Refine_Imperative_HOL.IICF
begin

  text ‹We now implement the Edmonds-Karp algorithm.
    Note that, during the implementation, we explicitly write down the 
    whole refined algorithm several times. As refinement is modular, most 
    of these copies could be avoided--- we inserted them deliberately for
    documentation purposes.
    ›

  subsection ‹Refinement to Residual Graph›
    text ‹As a first step towards implementation, we refine the algorithm
      to work directly on residual graphs. For this, we first have to 
      establish a relation between flows in a network and residual graphs.
      ›
    
  subsubsection ‹Refinement of Operations›
  context Network 
  begin
    text ‹We define the relation between residual graphs and flows›
    definition "cfi_rel  br flow_of_cf (RGraph c s t)"

    text ‹It can also be characterized the other way round, i.e., 
      mapping flows to residual graphs:›
    lemma cfi_rel_alt: "cfi_rel = {(cf,f). cf = residualGraph c f  NFlow c s t f}"
      unfolding cfi_rel_def br_def
      by (auto 
          simp: NFlow.is_RGraph RGraph.is_NFlow 
          simp: RPreGraph.rg_fo_inv[OF RGraph.this_loc_rpg]
          simp: NPreflow.fo_rg_inv[OF NFlow.axioms(1)])

    
    text ‹Initially, the residual graph for the zero flow equals the original network›
    lemma residualGraph_zero_flow: "residualGraph c (λ_. 0) = c" 
      unfolding residualGraph_def by (auto intro!: ext)
    lemma flow_of_c: "flow_of_cf c = (λ_. 0)"
      by (auto simp add: flow_of_cf_def[abs_def])

    text ‹The residual capacity is naturally defined on residual graphs›
    definition "resCap_cf cf p  Min {cf e | e. eset p}"
    lemma (in NFlow) resCap_cf_refine: "resCap_cf cf p = resCap p"
      unfolding resCap_cf_def resCap_def ..

    text ‹Augmentation can be done by @{const Graph.augment_cf}.› 

    
    lemma (in NFlow) augment_cf_refine_aux: (* For snippet *)
      assumes AUG: "isAugmentingPath p"
      shows "residualGraph c (augment (augmentingFlow p)) (u,v) = (
        if (u,v)set p then (residualGraph c f (u,v) - resCap p)
        else if (v,u)set p then (residualGraph c f (u,v) + resCap p)
        else residualGraph c f (u,v))"
      using augment_alt[OF AUG] by (auto simp: Graph.augment_cf_def)

    lemma augment_cf_refine:
      assumes R: "(cf,f)cfi_rel"
      assumes AUG: "NPreflow.isAugmentingPath c s t f p"
      shows "(Graph.augment_cf cf (set p) (resCap_cf cf p), 
          NFlow.augment_with_path c f p)  cfi_rel"
    proof -    
      from R have [simp]: "cf = residualGraph c f" and "NFlow c s t f"
        by (auto simp: cfi_rel_alt br_def)
      then interpret f: NFlow c s t f by simp
      
      show ?thesis 
        unfolding f.augment_with_path_def
      proof (simp add: cfi_rel_alt; safe intro!: ext)
        fix u v
        show "Graph.augment_cf f.cf (set p) (resCap_cf f.cf p) (u,v) 
              = residualGraph c (f.augment (f.augmentingFlow p)) (u,v)"
          unfolding f.augment_cf_refine_aux[OF AUG]
          unfolding f.cf.augment_cf_def
          by (auto simp: f.resCap_cf_refine)
      qed (rule f.augment_pres_nflow[OF AUG])
    qed  

    text ‹We rephrase the specification of shortest augmenting path to
      take a residual graph as parameter›
    (* TODO: This actually rephrases the spec to make it look more similar to 
      what BFS does later. This rephrasing does not belong here, but where we 
      implement it with BFS. *)
    definition "find_shortest_augmenting_spec_cf cf  
      assert (RGraph c s t cf) 
      SPEC (λ
        None  ¬Graph.connected cf s t 
      | Some p  Graph.isShortestPath cf s p t)"

    lemma (in RGraph) find_shortest_augmenting_spec_cf_refine: 
       "find_shortest_augmenting_spec_cf cf 
       find_shortest_augmenting_spec (flow_of_cf cf)"
      unfolding f_def[symmetric]
      unfolding find_shortest_augmenting_spec_cf_def 
        and find_shortest_augmenting_spec_def
      by (auto 
        simp: pw_le_iff refine_pw_simps 
        simp: this_loc rg_is_cf
        simp: f.isAugmentingPath_def Graph.connected_def Graph.isSimplePath_def 
        dest: cf.shortestPath_is_path
        split: option.split)
      
    text ‹This leads to the following refined algorithm›  
    definition "edka2  do {
      let cf = c;

      (cf,_)  whileT 
        (λ(cf,brk). ¬brk) 
        (λ(cf,_). do {
          assert (RGraph c s t cf);
          p  find_shortest_augmenting_spec_cf cf;
          case p of 
            None  return (cf,True)
          | Some p  do {
              assert (p[]);
              assert (Graph.isShortestPath cf s p t);
              let cf = Graph.augment_cf cf (set p) (resCap_cf cf p);
              assert (RGraph c s t cf);
              return (cf, False)
            }  
        })
        (cf,False);
      assert (RGraph c s t cf);
      let f = flow_of_cf cf;  
      return f
    }"

    lemma edka2_refine: "edka2  Id edka"
    proof -
      have [refine_dref_RELATES]: "RELATES cfi_rel" by (simp add: RELATES_def)

      show ?thesis
        unfolding edka2_def edka_def
        (*apply (rewrite in "let f' = NFlow.augmentingFlow c _ _ in _" Let_def)
        apply (rewrite in "let f = flow_of_cf _ in _" Let_def)*)
        apply (refine_rcg)
        apply refine_dref_type
        apply vc_solve

        ― ‹Solve some left-over verification conditions one by one›
        apply (drule NFlow.is_RGraph; 
            auto simp: cfi_rel_def br_def residualGraph_zero_flow flow_of_c; 
            fail)
        apply (auto simp: cfi_rel_def br_def; fail)
        using RGraph.find_shortest_augmenting_spec_cf_refine
        apply (auto simp: cfi_rel_def br_def; fail)
        apply (auto simp: cfi_rel_def br_def simp: RPreGraph.rg_fo_inv[OF RGraph.this_loc_rpg]; fail)
        apply (drule (1) augment_cf_refine; simp add: cfi_rel_def br_def; fail)
        apply (simp add: augment_cf_refine; fail)
        apply (auto simp: cfi_rel_def br_def; fail)
        apply (auto simp: cfi_rel_def br_def; fail)
        done

    qed    

    subsection ‹Implementation of Bottleneck Computation and Augmentation›  
    text ‹We will access the capacities in the residual graph
      only by a get-operation, which asserts that the edges are valid›
    
    abbreviation (input) valid_edge :: "edge  bool" where
      "valid_edge  λ(u,v). uV  vV"

    definition cf_get 
      :: "'capacity graph  edge  'capacity nres" 
      where "cf_get cf e  ASSERT (valid_edge e)  RETURN (cf e)"  
    definition cf_set 
      :: "'capacity graph  edge  'capacity  'capacity graph nres"
      where "cf_set cf e cap  ASSERT (valid_edge e)  RETURN (cf(e:=cap))"  

    definition resCap_cf_impl :: "'capacity graph  path  'capacity nres" 
    where "resCap_cf_impl cf p  
      case p of
        []  RETURN (0::'capacity)
      | (e#p)  do {
          cap  cf_get cf e;
          ASSERT (distinct p);
          nfoldli 
            p (λ_. True)
            (λe cap. do {
              cape  cf_get cf e;
              RETURN (min cape cap)
            }) 
            cap
        }"

    lemma (in RGraph) resCap_cf_impl_refine:   
      assumes AUG: "cf.isSimplePath s p t"
      shows "resCap_cf_impl cf p  SPEC (λr. r = resCap_cf cf p)"
    proof -
      (* TODO: Can we exploit Min.set_eq_fold *)

      note [simp del] = Min_insert
      note [simp] = Min_insert[symmetric]
      from AUG[THEN cf.isSPath_distinct]
      have "distinct p" .
      moreover from AUG cf.isPath_edgeset have "set p  cf.E"
        by (auto simp: cf.isSimplePath_def)
      hence "set p  Collect valid_edge"  
        using cf.E_ss_VxV by simp
      moreover from AUG have "p[]" by (auto simp: s_not_t) 
        then obtain e p' where "p=e#p'" by (auto simp: neq_Nil_conv)
      ultimately show ?thesis  
        unfolding resCap_cf_impl_def resCap_cf_def cf_get_def
        apply (simp only: list.case)
        apply (refine_vcg nfoldli_rule[where 
            I = "λl l' cap. 
              cap = Min (cf`insert e (set l)) 
             set (l@l')  Collect valid_edge"])
        apply (auto intro!: arg_cong[where f=Min])
        done
    qed    

    definition (in Graph) 
      "augment_edge e cap  (c(
                  e := c e - cap, 
        prod.swap e := c (prod.swap e) + cap))"

    (* TODO: This would be much simpler to prove if we had a characterization 
      of simple-path only depending on p. *)    
    lemma (in Graph) augment_cf_inductive:
      fixes e cap
      defines "c'  augment_edge e cap"
      assumes P: "isSimplePath s (e#p) t"
      shows "augment_cf (insert e (set p)) cap = Graph.augment_cf c' (set p) cap"
      and "s'. Graph.isSimplePath c' s' p t"  
    proof -
      obtain u v where [simp]: "e=(u,v)" by (cases e)

      from isSPath_no_selfloop[OF P] have [simp]: "u. (u,u)set p" "uv" by auto

      from isSPath_nt_parallel[OF P] have [simp]: "(v,u)set p" by auto
      from isSPath_distinct[OF P] have [simp]: "(u,v)set p" by auto

      show "augment_cf (insert e (set p)) cap = Graph.augment_cf c' (set p) cap"
        apply (rule ext)  
        unfolding Graph.augment_cf_def c'_def Graph.augment_edge_def
        by auto

      have "Graph.isSimplePath c' v p t"  
        unfolding Graph.isSimplePath_def
        apply rule
        apply (rule transfer_path)
        unfolding Graph.E_def
        apply (auto simp: c'_def Graph.augment_edge_def) []
        using P apply (auto simp: isSimplePath_def) []
        using P apply (auto simp: isSimplePath_def) []
        done
      thus "s'. Graph.isSimplePath c' s' p t" .. 

    qed    
        
    definition "augment_edge_impl cf e cap  do {
      v  cf_get cf e; cf  cf_set cf e (v-cap);
      let e = prod.swap e;
      v  cf_get cf e; cf  cf_set cf e (v+cap);
      RETURN cf
    }"

    lemma augment_edge_impl_refine: 
      assumes "valid_edge e" "u. e(u,u)"
      shows "augment_edge_impl cf e cap 
           (spec r. r = Graph.augment_edge cf e cap)"
      using assms
      unfolding augment_edge_impl_def Graph.augment_edge_def 
      unfolding cf_get_def cf_set_def
      apply refine_vcg
      apply auto
      done
      
    definition augment_cf_impl 
      :: "'capacity graph  path  'capacity  'capacity graph nres" 
      where
      "augment_cf_impl cf p x  do {
        (recT D. λ
          ([],cf)  return cf
        | (e#p,cf)  do {
            cf  augment_edge_impl cf e x;
            D (p,cf)
          }  
        ) (p,cf)
      }"

    text ‹Deriving the corresponding recursion equations›  
    lemma augment_cf_impl_simps[simp]: 
      "augment_cf_impl cf [] x = return cf"
      "augment_cf_impl cf (e#p) x = do { 
        cf  augment_edge_impl cf e x; 
        augment_cf_impl cf p x}"
      apply (simp add: augment_cf_impl_def)
      apply (subst RECT_unfold, refine_mono)
      apply simp
      
      apply (simp add: augment_cf_impl_def)
      apply (subst RECT_unfold, refine_mono)
      apply simp
      done      

    lemma augment_cf_impl_aux:    
      assumes "eset p. valid_edge e"
      assumes "s. Graph.isSimplePath cf s p t"
      shows "augment_cf_impl cf p x  RETURN (Graph.augment_cf cf (set p) x)"
      using assms
      apply (induction p arbitrary: cf)
      apply (simp add: Graph.augment_cf_empty)

      apply clarsimp
      apply (subst Graph.augment_cf_inductive, assumption)

      apply (refine_vcg augment_edge_impl_refine[THEN order_trans])
      apply simp
      apply simp
      apply (auto dest: Graph.isSPath_no_selfloop) []
      apply (rule order_trans, rprems)
        apply (drule Graph.augment_cf_inductive(2)[where cap=x]; simp)
        apply simp
      done  
      
    lemma (in RGraph) augment_cf_impl_refine:     
      assumes "Graph.isSimplePath cf s p t"
      shows "augment_cf_impl cf p x  RETURN (Graph.augment_cf cf (set p) x)"
      apply (rule augment_cf_impl_aux)
      using assms cf.E_ss_VxV apply (auto simp: cf.isSimplePath_def dest!: cf.isPath_edgeset) []
      using assms by blast
      
    text ‹Finally, we arrive at the algorithm where augmentation is 
      implemented algorithmically: ›  
    definition "edka3  do {
      let cf = c;

      (cf,_)  whileT 
        (λ(cf,brk). ¬brk) 
        (λ(cf,_). do {
          assert (RGraph c s t cf);
          p  find_shortest_augmenting_spec_cf cf;
          case p of 
            None  return (cf,True)
          | Some p  do {
              assert (p[]);
              assert (Graph.isShortestPath cf s p t);
              bn  resCap_cf_impl cf p;
              cf  augment_cf_impl cf p bn;
              assert (RGraph c s t cf);
              return (cf, False)
            }  
        })
        (cf,False);
      assert (RGraph c s t cf);
      let f = flow_of_cf cf;  
      return f
    }"

    lemma edka3_refine: "edka3  Id edka2"
      unfolding edka3_def edka2_def
      apply (rewrite in "let cf = Graph.augment_cf _ _ _ in _" Let_def)
      apply refine_rcg
      apply refine_dref_type
      apply (vc_solve)
      apply (drule Graph.shortestPath_is_simple)
      apply (frule (1) RGraph.resCap_cf_impl_refine)
      apply (frule (1) RGraph.augment_cf_impl_refine)
      apply (auto simp: pw_le_iff refine_pw_simps)
      done
      
    
    subsection ‹Refinement to use BFS›

    text ‹We refine the Edmonds-Karp algorithm to use breadth first search (BFS)›
    definition "edka4  do {
      let cf = c;

      (cf,_)  whileT 
        (λ(cf,brk). ¬brk) 
        (λ(cf,_). do {
          assert (RGraph c s t cf);
          p  Graph.bfs cf s t;
          case p of 
            None  return (cf,True)
          | Some p  do {
              assert (p[]);
              assert (Graph.isShortestPath cf s p t);
              bn  resCap_cf_impl cf p;
              cf  augment_cf_impl cf p bn;
              assert (RGraph c s t cf);
              return (cf, False)
            }  
        })
        (cf,False);
      assert (RGraph c s t cf);
      let f = flow_of_cf cf;  
      return f
    }"

    text ‹A shortest path can be obtained by BFS›  
    lemma bfs_refines_shortest_augmenting_spec: 
      "Graph.bfs cf s t  find_shortest_augmenting_spec_cf cf"
      unfolding find_shortest_augmenting_spec_cf_def
      apply (rule le_ASSERTI)
      apply (rule order_trans)
      apply (rule Graph.bfs_correct)
      apply (simp add: RPreGraph.resV_netV[OF RGraph.this_loc_rpg] s_node)
      apply (simp add: RPreGraph.resV_netV[OF RGraph.this_loc_rpg])
      apply (simp)
      done

    lemma edka4_refine: "edka4  Id edka3"
      unfolding edka4_def edka3_def
      apply refine_rcg
      apply refine_dref_type
      apply (vc_solve simp: bfs_refines_shortest_augmenting_spec)
      done


    subsection ‹Implementing the Successor Function for BFS›  

    text ‹We implement the successor function in two steps.
      The first step shows how to obtain the successor function by
      filtering the list of adjacent nodes. This step contains the idea   
      of the implementation. The second step is purely technical, and makes 
      explicit the recursion of the filter function as a recursion combinator
      in the monad. This is required for the Sepref tool.
      ›

    text ‹Note: We use @{term filter_rev} here, as it is tail-recursive, 
      and we are not interested in the order of successors.›
    definition "rg_succ am cf u   
      filter_rev (λv. cf (u,v) > 0) (am u)"
  
    lemma (in RGraph) rg_succ_ref1: "is_adj_map am 
       (rg_succ am cf u, Graph.E cf``{u})  Idlist_set_rel"
      unfolding Graph.E_def
      apply (clarsimp simp: list_set_rel_def br_def rg_succ_def filter_rev_alt; 
        intro conjI)
      using cfE_ss_invE resE_nonNegative 
      apply (auto 
        simp: is_adj_map_def less_le Graph.E_def 
        simp del: cf.zero_cap_simp zero_cap_simp) []
      apply (auto simp: is_adj_map_def) []
      done

    definition ps_get_op :: "_  node  node list nres" 
      where "ps_get_op am u  assert (uV)  return (am u)"

    definition monadic_filter_rev_aux 
      :: "'a list  ('a  bool nres)  'a list  'a list nres"
    where
      "monadic_filter_rev_aux a P l  (recT D. (λ(l,a). case l of
        []  return a 
      | (v#l)  do {
          c  P v;
          let a = (if c then v#a else a);
          D (l,a)
        }
      )) (l,a)"

    lemma monadic_filter_rev_aux_rule:
      assumes "x. xset l  P x  SPEC (λr. r=Q x)"
      shows "monadic_filter_rev_aux a P l  SPEC (λr. r=filter_rev_aux a Q l)"
      using assms
      apply (induction l arbitrary: a)

      apply (unfold monadic_filter_rev_aux_def) []
      apply (subst RECT_unfold, refine_mono)
      apply (fold monadic_filter_rev_aux_def) []
      apply simp

      apply (unfold monadic_filter_rev_aux_def) []
      apply (subst RECT_unfold, refine_mono)
      apply (fold monadic_filter_rev_aux_def) []
      apply (auto simp: pw_le_iff refine_pw_simps)
      done

    definition "monadic_filter_rev = monadic_filter_rev_aux []"

    lemma monadic_filter_rev_rule:
      assumes "x. xset l  P x  (spec r. r=Q x)"
      shows "monadic_filter_rev P l  (spec r. r=filter_rev Q l)"
      using monadic_filter_rev_aux_rule[where a="[]"] assms
      by (auto simp: monadic_filter_rev_def filter_rev_def)

    definition "rg_succ2 am cf u  do {
      l  ps_get_op am u;
      monadic_filter_rev (λv. do {
        x  cf_get cf (u,v);
        return (x>0)
      }) l
    }"

    lemma (in RGraph) rg_succ_ref2: 
      assumes PS: "is_adj_map am" and V: "uV"
      shows "rg_succ2 am cf u  return (rg_succ am cf u)"
    proof -
      have "vset (am u). valid_edge (u,v)"
        using PS V
        by (auto simp: is_adj_map_def Graph.V_def)
      
      thus ?thesis  
        unfolding rg_succ2_def rg_succ_def ps_get_op_def cf_get_def
        apply (refine_vcg monadic_filter_rev_rule[
            where Q="(λv. 0 < cf (u, v))", THEN order_trans])
        by (vc_solve simp: V)
    qed    

    lemma (in RGraph) rg_succ_ref:
      assumes A: "is_adj_map am"
      assumes B: "uV"
      shows "rg_succ2 am cf u  SPEC (λl. (l,cf.E``{u})  Idlist_set_rel)"
      using rg_succ_ref1[OF A, of u] rg_succ_ref2[OF A B]
      by (auto simp: pw_le_iff refine_pw_simps)


    subsection ‹Adding Tabulation of Input›  
    text ‹
      Next, we add functions that will be refined to tabulate the input of 
      the algorithm, i.e., the network's capacity matrix and adjacency map,
      into efficient representations. 
      The capacity matrix is tabulated to give the initial residual graph,
      and the adjacency map is tabulated for faster access.

      Note, on the abstract level, the tabulation functions are just identity,
      and merely serve as marker constants for implementation.
      ›
    definition init_cf :: "'capacity graph nres" 
      ― ‹Initialization of residual graph from network›
      where "init_cf  RETURN c"
    definition init_ps :: "(node  node list)  _" 
      ― ‹Initialization of adjacency map›
      where "init_ps am  ASSERT (is_adj_map am)  RETURN am"

    definition compute_rflow :: "'capacity graph  'capacity flow nres" 
      ― ‹Extraction of result flow from residual graph›
      where
      "compute_rflow cf  ASSERT (RGraph c s t cf)  RETURN (flow_of_cf cf)"

    definition "bfs2_op am cf  Graph.bfs2 cf (rg_succ2 am cf) s t"

    text ‹We split the algorithm into a tabulation function, and the 
      running of the actual algorithm:›
    definition "edka5_tabulate am  do {
      cf  init_cf;
      am  init_ps am;
      return (cf,am)
    }"

    definition "edka5_run cf am  do {
      (cf,_)  whileT 
        (λ(cf,brk). ¬brk) 
        (λ(cf,_). do {
          assert (RGraph c s t cf);
          p  bfs2_op am cf;
          case p of 
            None  return (cf,True)
          | Some p  do {
              assert (p[]);
              assert (Graph.isShortestPath cf s p t);
              bn  resCap_cf_impl cf p;
              cf  augment_cf_impl cf p bn;
              assert (RGraph c s t cf);
              return (cf, False)
            }  
        })
        (cf,False);
      f  compute_rflow cf;  
      return f
    }"

    definition "edka5 am  do {
      (cf,am)  edka5_tabulate am;
      edka5_run cf am
    }"

    lemma edka5_refine: "is_adj_map am  edka5 am  Id edka4"
      unfolding edka5_def edka5_tabulate_def edka5_run_def
        edka4_def init_cf_def compute_rflow_def
        init_ps_def Let_def nres_monad_laws bfs2_op_def
      apply refine_rcg
      apply refine_dref_type
      apply (vc_solve simp: )
      apply (rule refine_IdD)
      apply (rule Graph.bfs2_refine)
      apply (simp add: RPreGraph.resV_netV[OF RGraph.this_loc_rpg])
      apply (simp add: RGraph.rg_succ_ref)
      done

  end    

  subsection ‹Imperative Implementation›  
  text ‹In this section we provide an efficient imperative implementation,
    using the Sepref tool. It is mostly technical, setting up the mappings
    from abstract to concrete data structures, and then refining the algorithm,
    function by function.  
    ›

  text ‹
    This is also the point where we have to choose the implementation of 
    capacities. Up to here, they have been a polymorphic type with a
    typeclass constraint of being a linearly ordered integral domain.
    Here, we switch to @{typ [source] capacity_impl} (@{typ capacity_impl}).
    ›
  locale Network_Impl = Network c s t for c :: "capacity_impl graph" and s t

  text ‹Moreover, we assume that the nodes are natural numbers less 
    than some number @{term N}, which will become an additional parameter 
    of our algorithm. ›
  locale Edka_Impl = Network_Impl +
    fixes N :: nat
    assumes V_ss: "V{0..<N}"
  begin  
    lemma this_loc: "Edka_Impl c s t N" by unfold_locales

    lemma E_ss: "E  {0..<N}×{0..<N}" using E_ss_VxV V_ss by auto

    lemma mtx_nonzero_iff[simp]: "mtx_nonzero c = E" unfolding E_def by (auto simp: mtx_nonzero_def)

    lemma mtx_nonzeroN: "mtx_nonzero c  {0..<N}×{0..<N}" using E_ss by simp

    lemma [simp]: "vV  v<N" using V_ss by auto


    text ‹Declare some variables to Sepref. ›
    lemmas [id_rules] = 
      itypeI[Pure.of N "TYPE(nat)"]  
      itypeI[Pure.of s "TYPE(node)"]  
      itypeI[Pure.of t "TYPE(node)"]  
      itypeI[Pure.of c "TYPE(capacity_impl graph)"]  
    text ‹Instruct Sepref to not refine these parameters. This is expressed
      by using identity as refinement relation.›
    lemmas [sepref_import_param] = 
      IdI[of N]
      IdI[of s]
      IdI[of t]
      (*IdI[of c]*)

    lemma [sepref_fr_rules]: "(uncurry0 (return c),uncurry0 (return c))unit_assnk a pure (nat_rel×rnat_rel  int_rel)"
      apply sepref_to_hoare by sep_auto


    subsubsection ‹Implementation of Adjacency Map by Array›  
    definition "is_am am psi 
       Al. psi a l 
          * (length l = N  (i<N. l!i = am i) 
               (iN. am i = []))"
  
    lemma is_am_precise[safe_constraint_rules]: "precise (is_am)"
      apply rule
      unfolding is_am_def
      apply clarsimp
      apply (rename_tac l l')
      apply prec_extract_eqs
      apply (rule ext)
      apply (rename_tac i)
      apply (case_tac "i<length l'")
      apply fastforce+
      done

    sepref_decl_intf i_ps is "nat  nat list" 

    definition (in -) "ps_get_imp psi u  Array.nth psi u"

    lemma [def_pat_rules]: "Network.ps_get_op$c  UNPROTECT ps_get_op" by simp
    sepref_register "PR_CONST ps_get_op" :: "i_ps  node  node list nres"

    lemma ps_get_op_refine[sepref_fr_rules]: 
      "(uncurry ps_get_imp, uncurry (PR_CONST ps_get_op)) 
         is_amk *a (pure Id)k a list_assn (pure Id)"
      unfolding list_assn_pure_conv
      apply sepref_to_hoare
      using V_ss
      by (sep_auto 
            simp: is_am_def pure_def ps_get_imp_def 
            simp: ps_get_op_def refine_pw_simps)

    lemma is_pred_succ_no_node: "is_adj_map a; uV  a u = []"
      unfolding is_adj_map_def V_def
      by auto

    lemma [sepref_fr_rules]: "(Array.make N, PR_CONST init_ps) 
       (pure Id)k a is_am" 
      apply sepref_to_hoare
      using V_ss
      by (sep_auto simp: init_ps_def refine_pw_simps is_am_def pure_def
        intro: is_pred_succ_no_node)

    lemma [def_pat_rules]: "Network.init_ps$c  UNPROTECT init_ps" by simp
    sepref_register "PR_CONST init_ps" :: "(node  node list)  i_ps nres"

    subsubsection ‹Implementation of Capacity Matrix by Array›  
    lemma [def_pat_rules]: "Network.cf_get$c  UNPROTECT cf_get" by simp
    lemma [def_pat_rules]: "Network.cf_set$c  UNPROTECT cf_set" by simp

    sepref_register 
      "PR_CONST cf_get" :: "capacity_impl i_mtx  edge  capacity_impl nres"
    sepref_register 
      "PR_CONST cf_set" :: "capacity_impl i_mtx  edge  capacity_impl 
         capacity_impl i_mtx nres"

    text ‹We have to link the matrix implementation, which encodes the bound, 
      to the abstract assertion of the bound›

    sepref_definition cf_get_impl is "uncurry (PR_CONST cf_get)" :: "(asmtx_assn N id_assn)k *a (prod_assn id_assn id_assn)k a id_assn"
      unfolding PR_CONST_def cf_get_def[abs_def]
      by sepref
    lemmas [sepref_fr_rules] = cf_get_impl.refine
    lemmas [sepref_opt_simps] = cf_get_impl_def

    sepref_definition cf_set_impl is "uncurry2 (PR_CONST cf_set)" :: "(asmtx_assn N id_assn)d *a (prod_assn id_assn id_assn)k *a id_assnk a asmtx_assn N id_assn"
      unfolding PR_CONST_def cf_set_def[abs_def]
      by sepref
    lemmas [sepref_fr_rules] = cf_set_impl.refine
    lemmas [sepref_opt_simps] = cf_set_impl_def


    sepref_thm init_cf_impl is "uncurry0 (PR_CONST init_cf)" :: "unit_assnk a asmtx_assn N id_assn"
      unfolding PR_CONST_def init_cf_def 
      using E_ss
      apply (rewrite op_mtx_new_def[of c, symmetric])
      apply (rewrite amtx_fold_custom_new[of N N])
      by sepref

    concrete_definition (in -) init_cf_impl uses Edka_Impl.init_cf_impl.refine_raw is "(uncurry0 ?f,_)_" 
    prepare_code_thms (in -) init_cf_impl_def
    lemmas [sepref_fr_rules] = init_cf_impl.refine[OF this_loc]  

    (* TODO: Use sepref to synthesize the get-operations! *)
    lemma amtx_cnv: "amtx_assn N M id_assn = IICF_Array_Matrix.is_amtx N M" 
      by (simp add: amtx_assn_def)

    (*

    lemma init_cf_imp_refine[sepref_fr_rules]: 
      "(uncurry0 (mtx_new N c), uncurry0 (PR_CONST init_cf)) 
        ∈ (pure unit_rel)ka (asmtx_assn N id_assn)"
      unfolding asmtx_cnv
      apply sepref_to_hoare
      using E_ss
      by (sep_auto simp: init_cf_def)
    *)  

    lemma [def_pat_rules]: "Network.init_cf$c  UNPROTECT init_cf" by simp
    sepref_register "PR_CONST init_cf" :: "capacity_impl i_mtx nres"

    subsubsection ‹Representing Result Flow as Residual Graph›
    definition (in Network_Impl) "is_rflow N f cfi 
       Acf. asmtx_assn N id_assn cf cfi * (RGraph c s t cf  f = flow_of_cf cf)"
    lemma is_rflow_precise[safe_constraint_rules]: "precise (is_rflow N)"
      apply rule
      unfolding is_rflow_def
      apply (clarsimp simp: amtx_assn_def)
      apply prec_extract_eqs
      apply simp
      done

    sepref_decl_intf i_rflow is "nat×nat  int"

    lemma [sepref_fr_rules]: 
      "(λcfi. return cfi, PR_CONST compute_rflow)  (asmtx_assn N id_assn)d a is_rflow N"
      unfolding amtx_cnv
      apply sepref_to_hoare
      apply (sep_auto simp: amtx_cnv compute_rflow_def is_rflow_def refine_pw_simps hn_ctxt_def)
      done

    lemma [def_pat_rules]: 
      "Network.compute_rflow$c$s$t  UNPROTECT compute_rflow" by simp
    sepref_register 
      "PR_CONST compute_rflow" :: "capacity_impl i_mtx  i_rflow nres"


    subsubsection ‹Implementation of Functions›  

    schematic_goal rg_succ2_impl:
      fixes am :: "node  node list" and cf :: "capacity_impl graph"
      notes [id_rules] = 
        itypeI[Pure.of u "TYPE(node)"]
        itypeI[Pure.of am "TYPE(i_ps)"]
        itypeI[Pure.of cf "TYPE(capacity_impl i_mtx)"]
      notes [sepref_import_param] = IdI[of N]
      notes [sepref_fr_rules] = HOL_list_empty_hnr
      shows "hn_refine (hn_ctxt is_am am psi * hn_ctxt (asmtx_assn N id_assn) cf cfi * hn_val nat_rel u ui) (?c::?'c Heap)  ?R (rg_succ2 am cf u)"
      unfolding rg_succ2_def APP_def monadic_filter_rev_def monadic_filter_rev_aux_def
      (* TODO: Make setting up combinators for sepref simpler, then we do not need to unfold! *)
      using [[id_debug, goals_limit = 1]]
      by sepref
    concrete_definition (in -) succ_imp uses Edka_Impl.rg_succ2_impl
    prepare_code_thms (in -) succ_imp_def

    lemma succ_imp_refine[sepref_fr_rules]: 
      "(uncurry2 (succ_imp N), uncurry2 (PR_CONST rg_succ2)) 
         is_amk *a (asmtx_assn N id_assn)k *a (pure Id)k a list_assn (pure Id)"
      apply rule
      using succ_imp.refine[OF this_loc]            
      by (auto simp: hn_ctxt_def mult_ac split: prod.split)

    lemma [def_pat_rules]: "Network.rg_succ2$c  UNPROTECT rg_succ2" by simp
    sepref_register 
      "PR_CONST rg_succ2" :: "i_ps  capacity_impl i_mtx  node  node list nres"

    
    lemma [sepref_import_param]: "(min,min)IdIdId" by simp


    abbreviation "is_path  list_assn (prod_assn (pure Id) (pure Id))"

    schematic_goal resCap_imp_impl:
      fixes am :: "node  node list" and cf :: "capacity_impl graph" and p pi
      notes [id_rules] = 
        itypeI[Pure.of p "TYPE(edge list)"]
        itypeI[Pure.of cf "TYPE(capacity_impl i_mtx)"]
      notes [sepref_import_param] = IdI[of N]
      shows "hn_refine 
        (hn_ctxt (asmtx_assn N id_assn) cf cfi * hn_ctxt is_path p pi) 
        (?c::?'c Heap)  ?R 
        (resCap_cf_impl cf p)"
      unfolding resCap_cf_impl_def APP_def
      using [[id_debug, goals_limit = 1]]
      by sepref
    concrete_definition (in -) resCap_imp uses Edka_Impl.resCap_imp_impl
    prepare_code_thms (in -) resCap_imp_def

    lemma resCap_impl_refine[sepref_fr_rules]: 
      "(uncurry (resCap_imp N), uncurry (PR_CONST resCap_cf_impl)) 
         (asmtx_assn N id_assn)k *a (is_path)k a (pure Id)"
      apply sepref_to_hnr
      apply (rule hn_refine_preI)
      apply (clarsimp 
        simp: uncurry_def list_assn_pure_conv hn_ctxt_def 
        split: prod.split)
      apply (clarsimp simp: pure_def)
      apply (rule hn_refine_cons[OF _ resCap_imp.refine[OF this_loc] _])
      apply (simp add: list_assn_pure_conv hn_ctxt_def)
      apply (simp add: pure_def)
      apply (sep_auto simp add: hn_ctxt_def pure_def intro!: enttI)
      apply (simp add: pure_def)
      done

    lemma [def_pat_rules]: 
      "Network.resCap_cf_impl$c  UNPROTECT resCap_cf_impl" 
      by simp
    sepref_register "PR_CONST resCap_cf_impl" 
      :: "capacity_impl i_mtx  path  capacity_impl nres"
    
    sepref_thm augment_imp is "uncurry2 (PR_CONST augment_cf_impl)" :: "((asmtx_assn N id_assn)d *a (is_path)k *a (pure Id)k a asmtx_assn N id_assn)"
      unfolding augment_cf_impl_def[abs_def] augment_edge_impl_def PR_CONST_def
      using [[id_debug, goals_limit = 1]]
      by sepref 
    concrete_definition (in -) augment_imp uses Edka_Impl.augment_imp.refine_raw is "(uncurry2 ?f,_)_"
    prepare_code_thms (in -) augment_imp_def
    lemma augment_impl_refine[sepref_fr_rules]: 
      "(uncurry2 (augment_imp N), uncurry2 (PR_CONST augment_cf_impl)) 
         (asmtx_assn N id_assn)d *a (is_path)k *a (pure Id)k a asmtx_assn N id_assn"
      using augment_imp.refine[OF this_loc] .

    lemma [def_pat_rules]: 
      "Network.augment_cf_impl$c  UNPROTECT augment_cf_impl" 
      by simp
    sepref_register "PR_CONST augment_cf_impl" 
      :: "capacity_impl i_mtx  path  capacity_impl  capacity_impl i_mtx nres"

    sublocale bfs: Impl_Succ 
      "snd" 
      "TYPE(i_ps × capacity_impl i_mtx)" 
      "PR_CONST (λ(am,cf). rg_succ2 am cf)" 
      "prod_assn is_am (asmtx_assn N id_assn)" 
      "λ(am,cf). succ_imp N am cf"
      unfolding APP_def
      apply unfold_locales
      apply (simp add: fold_partial_uncurry)
      apply (rule hfref_cons[OF succ_imp_refine[unfolded PR_CONST_def]])
      by auto
      
    definition (in -) "bfsi' N s t psi cfi 
       bfs_impl (λ(am, cf). succ_imp N am cf) (psi,cfi) s t"

    lemma [sepref_fr_rules]: 
      "(uncurry (bfsi' N s t),uncurry (PR_CONST bfs2_op)) 
         is_amk *a (asmtx_assn N id_assn)k a option_assn is_path"
      unfolding bfsi'_def[abs_def] bfs2_op_def[abs_def] 
      using bfs.bfs_impl_fr_rule unfolding bfs.op_bfs_def[abs_def]
      apply (clarsimp simp: hfref_def all_to_meta)
      apply (rule hn_refine_cons[rotated])
      apply rprems
      apply (sep_auto simp: pure_def intro!: enttI)
      apply (sep_auto simp: pure_def)
      apply (sep_auto simp: pure_def)
      done

    lemma [def_pat_rules]: "Network.bfs2_op$c$s$t  UNPROTECT bfs2_op" by simp
    sepref_register "PR_CONST bfs2_op" 
      :: "i_ps  capacity_impl i_mtx  path option nres"  


    schematic_goal edka_imp_tabulate_impl:
      notes [sepref_opt_simps] = heap_WHILET_def
      fixes am :: "node  node list" and cf :: "capacity_impl graph"
      notes [id_rules] = 
        itypeI[Pure.of am "TYPE(node  node list)"]
      notes [sepref_import_param] = IdI[of am]
      shows "hn_refine (emp) (?c::?'c Heap)  ?R (edka5_tabulate am)"
      unfolding edka5_tabulate_def
      using [[id_debug, goals_limit = 1]]
      by sepref

    concrete_definition (in -) edka_imp_tabulate 
      uses Edka_Impl.edka_imp_tabulate_impl
    prepare_code_thms (in -) edka_imp_tabulate_def

    lemma edka_imp_tabulate_refine[sepref_fr_rules]: 
      "(edka_imp_tabulate c N, PR_CONST edka5_tabulate) 
       (pure Id)k a prod_assn (asmtx_assn N id_assn) is_am"
      apply (rule)
      apply (rule hn_refine_preI)
      apply (clarsimp 
        simp: uncurry_def list_assn_pure_conv hn_ctxt_def 
        split: prod.split)
      apply (rule hn_refine_cons[OF _ edka_imp_tabulate.refine[OF this_loc]])
      apply (sep_auto simp: hn_ctxt_def pure_def)+
      done

    lemma [def_pat_rules]: 
      "Network.edka5_tabulate$c  UNPROTECT edka5_tabulate" 
      by simp
    sepref_register "PR_CONST edka5_tabulate"
      :: "(node  node list)  (capacity_impl i_mtx × i_ps) nres"


    schematic_goal edka_imp_run_impl:
      notes [sepref_opt_simps] = heap_WHILET_def
      fixes am :: "node  node list" and cf :: "capacity_impl graph"
      notes [id_rules] = 
        itypeI[Pure.of cf "TYPE(capacity_impl i_mtx)"]
        itypeI[Pure.of am "TYPE(i_ps)"]
      shows "hn_refine 
        (hn_ctxt (asmtx_assn N id_assn) cf cfi * hn_ctxt is_am am psi) 
        (?c::?'c Heap)  ?R  
        (edka5_run cf am)"
      unfolding edka5_run_def
      using [[id_debug, goals_limit = 1]]
      by sepref

    concrete_definition (in -) edka_imp_run uses Edka_Impl.edka_imp_run_impl
    prepare_code_thms (in -) edka_imp_run_def

    thm edka_imp_run_def
    lemma edka_imp_run_refine[sepref_fr_rules]: 
      "(uncurry (edka_imp_run s t N), uncurry (PR_CONST edka5_run)) 
         (asmtx_assn N id_assn)d *a (is_am)k a is_rflow N"
      apply rule
      apply (clarsimp 
        simp: uncurry_def list_assn_pure_conv hn_ctxt_def 
        split: prod.split)
      apply (rule hn_refine_cons[OF _ edka_imp_run.refine[OF this_loc] _])
      apply (sep_auto simp: hn_ctxt_def)+
      done

    lemma [def_pat_rules]: 
      "Network.edka5_run$c$s$t  UNPROTECT edka5_run" 
      by simp
    sepref_register "PR_CONST edka5_run" 
      :: "capacity_impl i_mtx  i_ps  i_rflow nres"


    schematic_goal edka_imp_impl:
      notes [sepref_opt_simps] = heap_WHILET_def
      fixes am :: "node  node list" and cf :: "capacity_impl graph"
      notes [id_rules] = 
        itypeI[Pure.of am "TYPE(node  node list)"]
      notes [sepref_import_param] = IdI[of am]
      shows "hn_refine (emp) (?c::?'c Heap)  ?R (edka5 am)"
      unfolding edka5_def
      using [[id_debug, goals_limit = 1]]
      by sepref

    concrete_definition (in -) edka_imp uses Edka_Impl.edka_imp_impl
    prepare_code_thms (in -) edka_imp_def
    lemmas edka_imp_refine = edka_imp.refine[OF this_loc]

    thm pat_rules TrueI def_pat_rules
    

  end

  export_code edka_imp checking SML_imp

  subsection ‹Correctness Theorem for Implementation›
  text ‹We combine all refinement steps to derive a correctness 
    theorem for the implementation›
  context Network_Impl begin
    theorem edka_imp_correct: 
      assumes VN: "Graph.V c  {0..<N}"
      assumes ABS_PS: "is_adj_map am"
      shows "
        <emp> 
          edka_imp c s t N am 
        <λfi. Af. is_rflow N f fi * (isMaxFlow f)>t"
    proof -
      interpret Edka_Impl by unfold_locales fact

      note edka5_refine[OF ABS_PS]
      also note edka4_refine                 
      also note edka3_refine    
      also note edka2_refine 
      also note edka_refine
      also note edka_partial_refine
      also note fofu_partial_correct
      finally have "edka5 am  SPEC isMaxFlow" .
      from hn_refine_ref[OF this edka_imp_refine]
      show ?thesis 
        by (simp add: hn_refine_def)
    qed
  end    
end