Theory Cblinfun_Code

section ‹‹Cblinfun_Code› -- Support for code generation›

text ‹This theory provides support for code generation involving on complex vector spaces and
  bounded operators (e.g., types ‹cblinfun› and ‹ell2›).
  To fully support code generation, in addition to importing this theory,
  one need to activate support for code generation (import theory ‹Jordan_Normal_Form.Matrix_Impl›)
  and for real and complex numbers (import theory ‹Real_Impl.Real_Impl› for support of reals of the 
  form ‹a + b * sqrt c› or ‹Algebraic_Numbers.Real_Factorization› (much slower) for support of algebraic reals;
  support of complex numbers comes "for free").

  The builtin support for real and complex numbers (in ‹Complex_Main›) is not sufficient because it
  does not support the computation of square-roots which are used in the setup below.

  It is also recommended to import ‹HOL-Library.Code_Target_Numeral› for faster support of nats 
  and integers.›

theory Cblinfun_Code
  imports
    Cblinfun_Matrix Containers.Set_Impl Jordan_Normal_Form.Matrix_Kernel
begin

no_notation "Lattice.meet" (infixl "⊓ı" 70)
no_notation "Lattice.join" (infixl "⊔ı" 65)
hide_const (open) Coset.kernel
hide_const (open) Matrix_Kernel.kernel
hide_const (open) Order.bottom Order.top

unbundle lattice_syntax
unbundle jnf_notation
unbundle cblinfun_notation



subsection ‹Code equations for cblinfun operators›

text ‹In this subsection, we define the code for all operations involving only 
  operators (no combinations of operators/vectors/subspaces)›


text ‹The following lemma registers cblinfun as an abstract datatype with 
  constructor const‹cblinfun_of_mat›.
  That means that in generated code, all cblinfun operators will be represented
  as term‹cblinfun_of_mat X› where X is a matrix.
  In code equations for operations involving operators (e.g., +), we
  can then write the equation directly in terms of matrices
  by writing, e.g., term‹mat_of_cblinfun (A+B)› in the lhs,
  and in the rhs we define the matrix that corresponds to the sum of A,B.
  In the rhs, we can access the matrices corresponding to A,B by
  writing term‹mat_of_cblinfun B›.
  (See, e.g., lemma @{thm [source] mat_of_cblinfun_plus}.)

  See cite‹"code-generation-tutorial"› for more information on 
  @{theory_text ‹[code abstype]›}.›

declare mat_of_cblinfun_inverse [code abstype]

declare mat_of_cblinfun_plus[code]
  ― ‹Code equation for addition of cblinfuns›

declare mat_of_cblinfun_id[code]
  ― ‹Code equation for computing the identity operator›

declare mat_of_cblinfun_1[code]
  ― ‹Code equation for computing the one-dimensional identity›

declare mat_of_cblinfun_zero[code]
  ― ‹Code equation for computing the zero operator›


declare mat_of_cblinfun_uminus[code]
  ― ‹Code equation for computing the unary minus on cblinfun's›


declare mat_of_cblinfun_minus[code]
  ― ‹Code equation for computing the difference of cblinfun's›


declare mat_of_cblinfun_classical_operator[code]
  ― ‹Code equation for computing the "classical operator"›

declare mat_of_cblinfun_explicit_cblinfun[code]
  ― ‹Code equation for computing the const‹explicit_cblinfun››

declare mat_of_cblinfun_compose[code]
  ― ‹Code equation for computing the composition/product of cblinfun's›

declare mat_of_cblinfun_scaleC[code]
  ― ‹Code equation for multiplication with complex scalar›

declare mat_of_cblinfun_scaleR[code]
  ― ‹Code equation for multiplication with real scalar›

declare mat_of_cblinfun_adj[code]
  ― ‹Code equation for computing the adj›

text ‹This instantiation defines a code equation for equality tests for cblinfun.›
instantiation cblinfun :: (onb_enum,onb_enum) equal begin
definition [code]: "equal_cblinfun M N ⟷ mat_of_cblinfun M = mat_of_cblinfun N" 
  for M N :: "'a ⇒CL 'b"
instance 
  apply intro_classes
  unfolding equal_cblinfun_def 
  using mat_of_cblinfun_inj injD by fastforce
end

subsection ‹Vectors›

text ‹In this section, we define code for operations on vectors. As with operators above,
  we do this by using an isomorphism between finite vectors
  (i.e., types T of sort ‹complex_vector›) and the type typ‹complex vec› from
  session‹Jordan_Normal_Form›. We have developed such an isomorphism in 
  theory ‹Cblinfun_Matrix› for 
  any type T of sort ‹onb_enum› (i.e., any type with a finite canonical orthonormal basis)
  as was done above for bounded operators.
  Unfortunately, we cannot declare code equations for a type class, 
  code equations must be related to a specific type constructor.
  So we give code definition only for vectors of type typ‹'a ell2› (where typ‹'a›
  must be of sort ‹enum› to make make sure that typ‹'a ell2› is finite dimensional).
  
  The isomorphism between typ‹'a ell2› is given by the constants ‹ell2_of_vec›
  and ‹vec_of_ell2› which are copies of the more general const‹basis_enum_of_vec›
  and const‹vec_of_basis_enum› but with a more restricted type to be usable in our code equations.
›

definition ell2_of_vec :: "complex vec ⇒ 'a::enum ell2" where "ell2_of_vec = basis_enum_of_vec"
definition vec_of_ell2 :: "'a::enum ell2 ⇒ complex vec" where "vec_of_ell2 = vec_of_basis_enum"

text ‹The following theorem registers the isomorphism ‹ell2_of_vec›/‹vec_of_ell2›
  for code generation. From now on,
  code for operations on typ‹_ ell2› can be expressed by declarations such
  as term‹vec_of_ell2 (f a b) = g (vec_of_ell2 a) (vec_of_ell2 b)›
  if the operation f on typ‹_ ell2› corresponds to the operation g on
  typ‹complex vec›.›

lemma vec_of_ell2_inverse [code abstype]:
  "ell2_of_vec (vec_of_ell2 B) = B" 
  unfolding ell2_of_vec_def vec_of_ell2_def
  by (rule vec_of_basis_enum_inverse)

text ‹This instantiation defines a code equation for equality tests for ell2.›
instantiation ell2 :: (enum) equal begin
definition [code]: "equal_ell2 M N ⟷ vec_of_ell2 M = vec_of_ell2 N" 
  for M N :: "'a::enum ell2"
instance 
  apply intro_classes
  unfolding equal_ell2_def
  by (metis vec_of_ell2_inverse)
end

lemma vec_of_ell2_carrier_vec[simp]: ‹vec_of_ell2 v ∈ carrier_vec CARD('a)› for v :: ‹'a::enum ell2›
  using vec_of_basis_enum_carrier_vec[of v]
  apply (simp only: canonical_basis_length canonical_basis_length_ell2)
  by (simp add: vec_of_ell2_def)

lemma vec_of_ell2_zero[code]:
  ― ‹Code equation for computing the zero vector›
  "vec_of_ell2 (0::'a::enum ell2) = zero_vec (CARD('a))"
  by (simp add: vec_of_ell2_def vec_of_basis_enum_zero)

lemma vec_of_ell2_ket[code]:
  ― ‹Code equation for computing a standard basis vector›
  "vec_of_ell2 (ket i) = unit_vec (CARD('a)) (enum_idx i)" 
  for i::"'a::enum"
  using vec_of_ell2_def vec_of_basis_enum_ket by metis

lemma vec_of_ell2_scaleC[code]: 
  ― ‹Code equation for multiplying a vector with a complex scalar›
  "vec_of_ell2 (scaleC a ψ) = smult_vec a (vec_of_ell2 ψ)"
  for ψ :: "'a::enum ell2"
  by (simp add: vec_of_ell2_def vec_of_basis_enum_scaleC)

lemma vec_of_ell2_scaleR[code]: 
  ― ‹Code equation for multiplying a vector with a real scalar›
  "vec_of_ell2 (scaleR a ψ) = smult_vec (complex_of_real a) (vec_of_ell2 ψ)"
  for ψ :: "'a::enum ell2"
  by (simp add: vec_of_ell2_def vec_of_basis_enum_scaleR)

lemma ell2_of_vec_plus[code]:
  ― ‹Code equation for adding vectors›
  "vec_of_ell2 (x + y) =  (vec_of_ell2 x) + (vec_of_ell2 y)" for x y :: "'a::enum ell2"
  by (simp add: vec_of_ell2_def vec_of_basis_enum_add) 

lemma ell2_of_vec_minus[code]:
  ― ‹Code equation for subtracting vectors›
  "vec_of_ell2 (x - y) =  (vec_of_ell2 x) - (vec_of_ell2 y)" for x y :: "'a::enum ell2"
  by (simp add: vec_of_ell2_def vec_of_basis_enum_minus)

lemma ell2_of_vec_uminus[code]:
  ― ‹Code equation for negating a vector›
  "vec_of_ell2 (- y) =  - (vec_of_ell2 y)" for y :: "'a::enum ell2"
  by (simp add: vec_of_ell2_def vec_of_basis_enum_uminus)

lemma cinner_ell2_code [code]: "cinner ψ φ = cscalar_prod (vec_of_ell2 φ) (vec_of_ell2 ψ)"
  ― ‹Code equation for the inner product of vectors›
  by (simp add: cscalar_prod_vec_of_basis_enum vec_of_ell2_def)

lemma norm_ell2_code [code]: 
  ― ‹Code equation for the norm of a vector›
  "norm ψ = norm_vec (vec_of_ell2 ψ)"
  by (simp add: norm_vec_of_basis_enum vec_of_ell2_def)

lemma times_ell2_code[code]: 
  ― ‹Code equation for the product in the algebra of one-dimensional vectors›
  fixes ψ φ :: "'a::{CARD_1,enum} ell2"
  shows "vec_of_ell2 (ψ * φ)
   = vec_of_list [vec_index (vec_of_ell2 ψ) 0 * vec_index (vec_of_ell2 φ) 0]"
  by (simp add: vec_of_ell2_def vec_of_basis_enum_times)

lemma divide_ell2_code[code]: 
  ― ‹Code equation for the product in the algebra of one-dimensional vectors›
  fixes ψ φ :: "'a::{CARD_1,enum} ell2"
  shows "vec_of_ell2 (ψ / φ)
   = vec_of_list [vec_index (vec_of_ell2 ψ) 0 / vec_index (vec_of_ell2 φ) 0]"
  by (simp add: vec_of_ell2_def vec_of_basis_enum_divide)

lemma inverse_ell2_code[code]: 
  ― ‹Code equation for the product in the algebra of one-dimensional vectors›
  fixes ψ :: "'a::{CARD_1,enum} ell2"
  shows "vec_of_ell2 (inverse ψ)
   = vec_of_list [inverse (vec_index (vec_of_ell2 ψ) 0)]"
  by (simp add: vec_of_ell2_def vec_of_basis_enum_to_inverse)

lemma one_ell2_code[code]: 
  ― ‹Code equation for the unit in the algebra of one-dimensional vectors›
  "vec_of_ell2 (1 :: 'a::{CARD_1,enum} ell2) = vec_of_list [1]"
  by (simp add: vec_of_ell2_def vec_of_basis_enum_1) 

subsection ‹Vector/Matrix›

text ‹We proceed to give code equations for operations involving both
  operators (cblinfun) and vectors. As explained above, we have to restrict
  the equations to vectors of type typ‹'a ell2› even though the theory is available
  for any type of class class‹onb_enum›. As a consequence, we run into an
  additional technicality now. For example, to define a code equation for applying
  an operator to a vector, we might try to give the following lemma:

⬚‹lemma cblinfun_apply_ell2[code]:
  "vec_of_ell2 (M *V x) = (mult_mat_vec (mat_of_cblinfun M) (vec_of_ell2 x))"
  by (simp add: mat_of_cblinfun_cblinfun_apply vec_of_ell2_def)›

  Unfortunately, this does not work, Isabelle produces the warning 
  "Projection as head in equation", most likely due to the fact that
  the type of term‹(*V)› in the equation is less general than the type of 
  term‹(*V)› (it is restricted to @{type ell2}). We overcome this problem
  by defining a constant ‹cblinfun_apply_ell2› which is equal to term‹(*V)›
  but has a more restricted type. We then instruct the code generation 
  to replace occurrences of term‹(*V)› by ‹cblinfun_apply_ell2› (where possible),
  and we add code generation for ‹cblinfun_apply_ell2› instead of term‹(*V)›.
›

definition cblinfun_apply_ell2 :: "'a ell2 ⇒CL 'b ell2 ⇒ 'a ell2 ⇒ 'b ell2" 
  where [code del, code_abbrev]: "cblinfun_apply_ell2 = (*V)"
    ― ‹@{attribute code_abbrev} instructs the code generation to replace the
     rhs term‹(*V)› by the lhs term‹cblinfun_apply_ell2› before starting 
     the actual code generation.›

lemma cblinfun_apply_ell2[code]:
  ― ‹Code equation for term‹cblinfun_apply_ell2›, i.e., for applying an operator
     to an type‹ell2› vector›
  "vec_of_ell2 (cblinfun_apply_ell2 M x) = (mult_mat_vec (mat_of_cblinfun M) (vec_of_ell2 x))"
  by (simp add: cblinfun_apply_ell2_def mat_of_cblinfun_cblinfun_apply vec_of_ell2_def)

text ‹For the constant term‹vector_to_cblinfun› (canonical isomorphism from
  vectors to operators), we have the same problem and define a constant
  ‹vector_to_cblinfun_code› with more restricted type›

definition vector_to_cblinfun_code :: "'a ell2 ⇒ 'b::one_dim ⇒CL 'a ell2" where
  [code del,code_abbrev]: "vector_to_cblinfun_code = vector_to_cblinfun"
  ― ‹@{attribute code_abbrev} instructs the code generation to replace the
     rhs term‹vector_to_cblinfun› by the lhs term‹vector_to_cblinfun_code›
     before starting the actual code generation.›

lemma vector_to_cblinfun_code[code]: 
  ― ‹Code equation for translating a vector into an operation (single-column matrix)›
  "mat_of_cblinfun (vector_to_cblinfun_code ψ) = mat_of_cols (CARD('a)) [vec_of_ell2 ψ]"
  for ψ::"'a::enum ell2"
  by (simp add: mat_of_cblinfun_vector_to_cblinfun  vec_of_ell2_def vector_to_cblinfun_code_def)

definition butterfly_code :: ‹'a ell2 ⇒ 'b ell2 ⇒ 'b ell2 ⇒CL 'a ell2› 
  where [code del, code_abbrev]: ‹butterfly_code = butterfly› 
lemma butterfly_code[code]: ‹mat_of_cblinfun (butterfly_code s t)
   = mat_of_cols (CARD('a)) [vec_of_ell2 s] * mat_of_rows (CARD('b)) [map_vec cnj (vec_of_ell2 t)]›
  for s :: ‹'a::enum ell2› and t :: ‹'b::enum ell2›
  by (simp add: butterfly_code_def butterfly_def vector_to_cblinfun_code mat_of_cblinfun_compose
      mat_of_cblinfun_adj mat_adjoint_def map_map_vec_cols
      flip: vector_to_cblinfun_code_def map_vec_conjugate[abs_def])

subsection ‹Subspaces›

text ‹In this section, we define code equations for handling subspaces, i.e.,
values of type typ‹'a ccsubspace›. We choose to computationally represent
a subspace by a list of vectors that span the subspace. That is,
if term‹vecs› are vectors (type typ‹complex vec›), ‹SPAN vecs› is defined to be their
span. Then the code generation can simply represent all subspaces in this form, and 
we need to define the operations on subspaces in terms of list of vectors 
(e.g., the closed union of two subspaces would be computed as the concatenation 
of the two lists, to give one of the simplest examples).

To support this, ‹SPAN› is declared as a "‹code_datatype›".
(Not as an abstract datatype like term‹cblinfun_of_mat›/term‹mat_of_cblinfun›
because that would require ‹SPAN› to be injective.)

Then all code equations for different operations need to be formulated as
functions of values of the form ‹SPAN x›. (E.g., ‹SPAN x + SPAN y = SPAN (…)›.)›

definition [code del]: "SPAN x = (let n = length (canonical_basis :: 'a::onb_enum list) in
    ccspan (basis_enum_of_vec ` Set.filter (λv. dim_vec v = n) (set x)) :: 'a ccsubspace)"
  ― ‹The SPAN of vectors x, as a type‹ccsubspace›.
      We filter out vectors of the wrong dimension because ‹SPAN› needs to have
      well-defined behavior even in cases that would not actually occur in an execution.›
code_datatype SPAN

text ‹We first declare code equations for term‹Proj›, i.e., for
turning a subspace into a projector. This means, we would need a code equation
of the form ‹mat_of_cblinfun (Proj (SPAN S)) = …›. However, this equation is
not accepted by the code generation for reasons we do not understand. But
if we define an auxiliary constant ‹mat_of_cblinfun_Proj_code› that stands for 
‹mat_of_cblinfun (Proj _)›, define a code equation for ‹mat_of_cblinfun_Proj_code›,
and then define a code equation for ‹mat_of_cblinfun (Proj S)› in terms of 
‹mat_of_cblinfun_Proj_code›, Isabelle accepts the code equations.›

definition "mat_of_cblinfun_Proj_code S = mat_of_cblinfun (Proj S)"
declare mat_of_cblinfun_Proj_code_def[symmetric, code]

lemma mat_of_cblinfun_Proj_code_code[code]: 
  ― ‹Code equation for computing a projector onto a set S of vectors.
     We first make the vectors S into an orthonormal basis using
     the Gram-Schmidt procedure and then compute the projector
     as the sum of the "butterflies" ‹x * x*› of the vectors ‹x∈S›
     (done by term‹mk_projector›).›
  "mat_of_cblinfun_Proj_code (SPAN S :: 'a::onb_enum ccsubspace) = 
    (let d = length (canonical_basis :: 'a list) in mk_projector d (filter (λv. dim_vec v = d) S))"
proof -
  have *: "map_option vec_of_basis_enum (if dim_vec x = length (canonical_basis :: 'a list) then Some (basis_enum_of_vec x :: 'a) else None)
      = (if dim_vec x = length (canonical_basis :: 'a list) then Some x else None)" for x
    by auto
  show ?thesis
    unfolding SPAN_def mat_of_cblinfun_Proj_code_def
    using mat_of_cblinfun_Proj_ccspan[where S = 
        "map basis_enum_of_vec (filter (λv. dim_vec v = (length (canonical_basis :: 'a list))) S) :: 'a list"]
    apply (simp only: Let_def map_filter_map_filter filter_set image_set map_map_filter o_def)
    unfolding *
    by (simp add: map_filter_map_filter[symmetric])
qed

lemma top_ccsubspace_code[code]: 
  ― ‹Code equation for term‹top›, the subspace containing everything.
      Top is represented as the span of the standard basis vectors.›
  "(top::'a ccsubspace) =
      (let n = length (canonical_basis :: 'a::onb_enum list) in SPAN (unit_vecs n))"
  unfolding SPAN_def
  apply (simp only: index_unit_vec Let_def map_filter_map_filter filter_set image_set map_map_filter 
      map_filter_map o_def unit_vecs_def)
  apply (simp add: basis_enum_of_vec_unit_vec)
  apply (subst nth_image)
  by auto

lemma bot_as_span[code]: 
  ― ‹Code equation for term‹bot›, the subspace containing everything.
      Top is represented as the span of the standard basis vectors.›
  "(bot::'a::onb_enum ccsubspace) = SPAN []"
  unfolding SPAN_def by (auto simp: Set.filter_def)


lemma sup_spans[code]:
  ― ‹Code equation for the join (lub) of two subspaces (union of the generating lists)›
  "SPAN A ⊔ SPAN B = SPAN (A @ B)"
  unfolding SPAN_def 
  by (auto simp: ccspan_union image_Un filter_Un Let_def)

text ‹We do not need an equation for term‹(+)› because term‹(+)›
is defined in terms of term‹(⊔)› (for type‹ccsubspace›), thus the code generation automatically
computes term‹(+)› in terms of the code for term‹(⊔)››

definition [code del,code_abbrev]: "Span_code (S::'a::enum ell2 set) = (ccspan S)"
  ― ‹A copy of term‹ccspan› with restricted type. For analogous reasons as
     term‹cblinfun_apply_ell2›, see there for explanations›

lemma span_Set_Monad[code]: "Span_code (Set_Monad l) = (SPAN (map vec_of_ell2 l))"
  ― ‹Code equation for the span of a finite set. (term‹Set_Monad› is a datatype
     constructor that represents sets as lists in the computation.)›
  apply (simp add: Span_code_def SPAN_def Let_def)
  apply (subst Set_filter_unchanged)
   apply (auto simp add: vec_of_ell2_def)[1]
  by (metis (no_types, lifting) ell2_of_vec_def image_image map_idI set_map vec_of_ell2_inverse)

text ‹This instantiation defines a code equation for equality tests for type‹ccsubspace›.
      The actual code for equality tests is given below (lemma ‹equal_ccsubspace_code›).›
instantiation ccsubspace :: (onb_enum) equal begin
definition [code del]: "equal_ccsubspace (A::'a ccsubspace) B = (A=B)"
instance apply intro_classes unfolding equal_ccsubspace_def by simp
end

lemma leq_ccsubspace_code[code]:
  ― ‹Code equation for deciding inclusion of one space in another.
     Uses the constant term‹is_subspace_of_vec_list› which implements the actual
     computation by checking for each generator of A whether it is in the
     span of B (by orthogonal projection onto an orthonormal basis of B
     which is computed using Gram-Schmidt).›
  "SPAN A ≤ (SPAN B :: 'a::onb_enum ccsubspace)
      ⟷ (let d = length (canonical_basis :: 'a list) in
          is_subspace_of_vec_list d
          (filter (λv. dim_vec v = d) A)
          (filter (λv. dim_vec v = d) B))"
proof -
  define d A' B' where "d = length (canonical_basis :: 'a list)"
    and "A' = filter (λv. dim_vec v = d) A"
    and "B' = filter (λv. dim_vec v = d) B"

  show ?thesis
    unfolding SPAN_def d_def[symmetric] filter_set Let_def
      A'_def[symmetric] B'_def[symmetric] image_set
    apply (subst ccspan_leq_using_vec)
    unfolding d_def[symmetric] map_map o_def
    apply (subst map_cong[where xs=A', OF refl])
     apply (rule basis_enum_of_vec_inverse)
     apply (simp add: A'_def d_def)
    apply (subst map_cong[where xs=B', OF refl])
     apply (rule basis_enum_of_vec_inverse)
    by (simp_all add: B'_def d_def)
qed

lemma equal_ccsubspace_code[code]:
  ― ‹Code equation for equality test. By checking mutual inclusion
      (for which we have code by the preceding code equation).›
  "HOL.equal (A::_ ccsubspace) B = (A≤B ∧ B≤A)"
  unfolding equal_ccsubspace_def by auto

lemma cblinfun_image_code[code]:
  ― ‹Code equation for applying an operator term‹A› to a subspace. 
      Simply by multiplying each generator with term‹A››
  "A *S SPAN S = (let d = length (canonical_basis :: 'a list) in
         SPAN (map (mult_mat_vec (mat_of_cblinfun A))
               (filter (λv. dim_vec v = d) S)))"
  for A::"'a::onb_enum ⇒CL'b::onb_enum"
proof -
  define dA dB S'
    where "dA = length (canonical_basis :: 'a list)"
      and "dB = length (canonical_basis :: 'b list)"
      and "S' = filter (λv. dim_vec v = dA) S"

  have "cblinfun_image A (SPAN S) = A *S ccspan (set (map basis_enum_of_vec S'))"
    unfolding SPAN_def dA_def[symmetric] Let_def S'_def filter_set
    by simp
  also have "… = ccspan ((λx. basis_enum_of_vec 
            (mat_of_cblinfun A *v vec_of_basis_enum (basis_enum_of_vec x :: 'a))) ` set S')"
    apply (subst cblinfun_image_ccspan_using_vec)
    by (simp add: image_image)
  also have "… = ccspan ((λx. basis_enum_of_vec (mat_of_cblinfun A *v x)) ` set S')"
    apply (subst image_cong[OF refl])
     apply (subst basis_enum_of_vec_inverse)
    by (auto simp add: S'_def dA_def)
  also have "… = SPAN (map (mult_mat_vec (mat_of_cblinfun A)) S')"
    unfolding SPAN_def dB_def[symmetric] Let_def filter_set 
    apply (subst filter_True)
    by (simp_all add: dB_def mat_of_cblinfun_def image_image)

  finally show ?thesis
    unfolding dA_def[symmetric] S'_def[symmetric] Let_def
    by simp
qed

definition [code del, code_abbrev]: "range_cblinfun_code A = A *S top"
  ― ‹A new constant for the special case of applying an operator to the subspace term‹top›
  (i.e., for computing the range of the operator). We do this to be able to give
  more specialized code for this specific situation. (The generic code for
  term‹(*S)› would work but is less efficient because it involves repeated matrix 
  multiplications. @{attribute code_abbrev} makes sure occurrences of term‹A *S top›
  are replaced before starting the actual code generation.›

lemma range_cblinfun_code[code]: 
  ― ‹Code equation for computing the range of an operator term‹A›.
      Returns the columns of the matrix representation of term‹A›.›
  fixes A :: "'a::onb_enum ⇒CL 'b::onb_enum"
  shows "range_cblinfun_code A = SPAN (cols (mat_of_cblinfun A))"
proof -
  define dA dB
    where "dA = length (canonical_basis :: 'a list)"
      and "dB = length (canonical_basis :: 'b list)"
  have carrier_A: "mat_of_cblinfun A ∈ carrier_mat dB dA"
    unfolding mat_of_cblinfun_def dA_def dB_def by simp

  have "range_cblinfun_code A = A *S SPAN (unit_vecs dA)"
    unfolding range_cblinfun_code_def
    by (metis dA_def top_ccsubspace_code)
  also have "… = SPAN (map (λi. mat_of_cblinfun A *v unit_vec dA i) [0..<dA])"
    unfolding cblinfun_image_code dA_def[symmetric] Let_def
    apply (subst filter_True)
     apply (meson carrier_vecD subset_code(1) unit_vecs_carrier)
    by (simp add: unit_vecs_def o_def)
  also have "… = SPAN (map (λx. mat_of_cblinfun A *v col (1m dA) x) [0..<dA])"
    apply (subst map_cong[OF refl])
    by auto
  also have "… = SPAN (map (col (mat_of_cblinfun A * 1m dA)) [0..<dA])"
    apply (subst map_cong[OF refl])
     apply (subst col_mult2[symmetric])
        apply (rule carrier_A)
    by auto
  also have "… = SPAN (cols (mat_of_cblinfun A))"
    unfolding cols_def dA_def[symmetric]
    apply (subst right_mult_one_mat[OF carrier_A])
    using carrier_A by blast
  finally show ?thesis
    by -
qed


lemma uminus_Span_code[code]: "- X = range_cblinfun_code (id_cblinfun - Proj X)"
  ― ‹Code equation for the orthogonal complement of a subspace term‹X›. 
      Computed as the range of one minus the projector on term‹X››
  unfolding range_cblinfun_code_def
  by (metis Proj_ortho_compl Proj_range)

lemma kernel_code[code]: 
  ― ‹Computes the kernel of an operator term‹A›.
      This is implemented using the existing functions 
      for transforming a matrix into row echelon form (term‹gauss_jordan_single›)
      and for computing a basis of the kernel of such a matrix
      (term‹find_base_vectors›)›
  "kernel A = SPAN (find_base_vectors (gauss_jordan_single (mat_of_cblinfun A)))" 
  for A::"('a::onb_enum,'b::onb_enum) cblinfun"
proof -
  define dA dB Am Ag base
    where "dA = length (canonical_basis :: 'a list)"
      and "dB = length (canonical_basis :: 'b list)"
      and "Am = mat_of_cblinfun A"
      and "Ag = gauss_jordan_single Am"
      and "base = find_base_vectors Ag"

  interpret complex_vec_space dA.

  have Am_carrier: "Am ∈ carrier_mat dB dA"
    unfolding Am_def mat_of_cblinfun_def dA_def dB_def by simp

  have row_echelon: "row_echelon_form Ag"
    unfolding Ag_def
    using Am_carrier refl by (rule gauss_jordan_single)

  have Ag_carrier: "Ag ∈ carrier_mat dB dA"
    unfolding Ag_def
    using Am_carrier refl by (rule gauss_jordan_single(2))

  have base_carrier: "set base ⊆ carrier_vec dA"
    unfolding base_def
    using find_base_vectors(1)[OF row_echelon Ag_carrier]
    using Ag_carrier mat_kernel_def by blast

  interpret k: kernel dB dA Ag
    apply standard using Ag_carrier by simp

  have basis_base: "kernel.basis dA Ag (set base)"
    using row_echelon Ag_carrier unfolding base_def
    by (rule find_base_vectors(3))


  have "space_as_set (SPAN base)
       = space_as_set (ccspan (basis_enum_of_vec ` set base :: 'a set))"
    unfolding SPAN_def dA_def[symmetric] Let_def filter_set
    apply (subst filter_True)
    using base_carrier by auto

  also have "… = cspan (basis_enum_of_vec ` set base)"
    apply transfer apply (subst closure_finite_cspan)
    by simp_all

  also have "… = basis_enum_of_vec ` span (set base)"
    apply (subst basis_enum_of_vec_span)
    using base_carrier dA_def by auto

  also have "… = basis_enum_of_vec ` mat_kernel Ag"
    using basis_base k.Ker.basis_def k.span_same by auto

  also have "… = basis_enum_of_vec ` {v ∈ carrier_vec dA. Ag *v v = 0v dB}"
    apply (rule arg_cong[where f="λx. basis_enum_of_vec ` x"])
    unfolding mat_kernel_def using Ag_carrier
    by simp

  also have "… = basis_enum_of_vec ` {v ∈ carrier_vec dA. Am *v v = 0v dB}"
    using gauss_jordan_single(1)[OF Am_carrier Ag_def[symmetric]]
    by auto

  also have "… = {w. A *V w = 0}"
  proof -
    have "basis_enum_of_vec ` {v ∈ carrier_vec dA. Am *v v = 0v dB}
        = basis_enum_of_vec ` {v ∈ carrier_vec dA. A *V basis_enum_of_vec v = 0}"
      apply (rule arg_cong[where f="λt. basis_enum_of_vec ` t"])
      apply (rule Collect_cong)
      apply (simp add: Am_def)
      by (metis Am_carrier Am_def carrier_matD(2) carrier_vecD dB_def mat_carrier 
          mat_of_cblinfun_def mat_of_cblinfun_cblinfun_apply vec_of_basis_enum_inverse 
          basis_enum_of_vec_inverse vec_of_basis_enum_zero)
    also have "… = {w ∈ basis_enum_of_vec ` carrier_vec dA. A *V w = 0}"
      apply (subst Compr_image_eq[symmetric])
      by simp
    also have "… = {w. A *V w = 0}"
      apply auto
      by (metis (no_types, lifting) Am_carrier Am_def carrier_matD(2) carrier_vec_dim_vec dim_vec_of_basis_enum' image_iff mat_carrier mat_of_cblinfun_def vec_of_basis_enum_inverse)
    finally show ?thesis
      by -
  qed

  also have "… = space_as_set (kernel A)"
    apply transfer by auto

  finally have "SPAN base = kernel A"
    by (simp add: space_as_set_inject)
  then show ?thesis
    by (simp add: base_def Ag_def Am_def)
qed

lemma inf_ccsubspace_code[code]: 
  ― ‹Code equation for intersection of subspaces.
     Reduced to orthogonal complement and sum of subspaces
     for which we already have code equations.›
  "(A::'a::onb_enum ccsubspace) ⊓ B = - (- A ⊔ - B)"
  by (subst ortho_involution[symmetric], subst compl_inf, simp)

lemma Sup_ccsubspace_code[code]:
  ― ‹Supremum (sum) of a set of subspaces. Implemented
     by repeated pairwise sum.›
  "Sup (Set_Monad l :: 'a::onb_enum ccsubspace set) = fold sup l bot"
  unfolding Set_Monad_def
  by (simp add: Sup_set_fold)


lemma Inf_ccsubspace_code[code]: 
  ― ‹Infimum (intersection) of a set of subspaces. 
      Implemented by the orthogonal complement of the supremum.›
  "Inf (Set_Monad l :: 'a::onb_enum ccsubspace set)
  = - Sup (Set_Monad (map uminus l))"
  unfolding Set_Monad_def
  apply (induction l)
  by auto

subsection ‹Miscellanea›

text ‹This is a hack to circumvent a bug in the code generation. The automatically
  generated code for the class class‹uniformity› has a type that is different from
  what the generated code later assumes, leading to compilation errors (in ML at least)
  in any expression involving typ‹_ ell2› (even if the constant const‹uniformity› is
  not actually used).
  
  The fragment below circumvents this by forcing Isabelle to use the right type.
  (The logically useless fragment "‹let x = ((=)::'a⇒_⇒_)›" achieves this.)›
lemma uniformity_ell2_code[code]: "(uniformity :: ('a ell2 * _) filter) = Filter.abstract_filter (%_.
    Code.abort STR ''no uniformity'' (%_. 
    let x = ((=)::'a⇒_⇒_) in uniformity))"
  by simp

text ‹Code equation for term‹UNIV›. 
  It is now implemented via type class class‹enum› 
  (which provides a list of all values).›
declare [[code drop: UNIV]]
declare enum_class.UNIV_enum[code]

text ‹Setup for code generation involving sets of type‹ell2›/type‹ccsubspace›.
  This configures to use lists for representing sets in code.›
derive (eq) ceq ccsubspace
derive (no) ccompare ccsubspace
derive (monad) set_impl ccsubspace
derive (eq) ceq ell2
derive (no) ccompare ell2
derive (monad) set_impl ell2


unbundle no_lattice_syntax
unbundle no_jnf_notation
unbundle no_cblinfun_notation

end