text:
# Continuous Lyapunov equations
"""
    X = lyapc(A, C)

Compute `X`, the symmetric or hermitian solution of the continuous Lyapunov equation

      AX + XA' + C = 0,

where `A` is a square real or complex matrix and `C` is a symmetric or hermitian
matrix. `A` must not have two eigenvalues `α` and `β` such that `α+β = 0`.

The following particular cases are also adressed:

    X = lyapc(α*I,C)  or  X = lyapc(α,C)

Solve the matrix equation `(α+α')X + C = 0`.

    x = lyapc(α,γ)

Solve the equation `(α+α')x + γ = 0`.

# Example
```jldoctest
julia> A = [3. 4.; 5. 6.]
2×2 Array{Float64,2}:
 3.0  4.0
 5.0  6.0

julia> C = [1. 1.; 1. 2.]
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  2.0

julia> X = lyapc(A, C)
2×2 Array{Float64,2}:
  0.5  -0.5
 -0.5   0.25

julia> A*X + X*A' + C
2×2 Array{Float64,2}:
 -8.88178e-16   2.22045e-16
  2.22045e-16  -4.44089e-16
```
"""
function lyapc(A::AbstractMatrix, C::AbstractMatrix)
   """
   The Bartels-Steward Schur form based method is employed.

   Reference:
   R. H. Bartels and G. W. Stewart. Algorithm 432: Solution of the matrix equation AX+XB=C.
   Comm. ACM, 15:820–826, 1972.
   """
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && ishermitian(C)) ||
      throw(DimensionMismatch("C must be a symmetric/hermitian matrix of dimension $n"))

   adj = isa(A,Adjoint)

   T2 = promote_type(eltype(A), eltype(C))
   T2 <: BlasFloat  || (T2 = promote_type(Float64,T2))
   eltype(A) == T2 || (adj ? A = convert(Matrix{T2},A.parent)' : A = convert(Matrix{T2},A))
   eltype(C) == T2 || (C = convert(Matrix{T2},C))

   # Reduce A to Schur form and transform C
   if adj
      AS, Q = schur(A.parent)
   else
      AS, Q = schur(A)
   end

   #X = Q'*C*Q
   X = utqu(C,Q)
   lyapcs!(AS, X, adj = adj)
   #X <- Q*X*Q'
   utqu!(X,Q')
end
# (α+α')X  + C = 0
lyapc(A::UniformScaling, C::AbstractMatrix) = -C/(A+A')
lyapc(A::Union{Real,Complex}, C::AbstractMatrix) = real(A) == 0 ? throw(SingularException(1)) : -C/(A+A')
# (α+α')x + γ = 0
lyapc(A::Union{Real,Complex}, C::Union{Real,Complex})  = real(A) == 0 ? throw(SingularException(1)) : -C/(A+A')
"""
    X = lyapc(A, E, C)

Compute `X`, the symmetric or hermitian solution of the
generalized continuous Lyapunov equation

     AXE' + EXA' + C = 0,

where `A` and `E` are square real or complex matrices and `C` is a symmetric or
hermitian matrix. The pencil `A-λE` must not have two eigenvalues `α` and `β` such that `α+β = 0`.

The following particular cases are also adressed:

    X = lyapc(A,β*I,C)  or  X = lyapc(A,β,C)

Solve the matrix equation `AXβ' + βXA' + C = 0`.

    X = lyapc(α*I,E,C)  or  X = lyapc(α,E,C)

Solve the matrix equation `αXE' + EXα' + C = 0`.

    X = lyapc(α*I,β*I,C)  or  X = lyapc(α,β,C)

Solve the matrix equation `(αβ'+α'β)X + C = 0`.

    x = lyapc(α,β,γ)

Solve the equation `(αβ'+α'β)x + γ = 0`.

# Example
```jldoctest
julia> A = [3. 4.; 5. 6.]
2×2 Array{Float64,2}:
 3.0  4.0
 5.0  6.0

julia> E = [ 1. 2.; 0. 1.]
2×2 Array{Float64,2}:
 1.0  2.0
 0.0  1.0

julia> C = [1. 1.; 1. 2.]
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  2.0

julia> X = lyapc(A, E, C)
2×2 Array{Float64,2}:
 -2.5   2.5
  2.5  -2.25

julia> A*X*E' + E*X*A' + C
2×2 Array{Float64,2}:
 -5.32907e-15  -2.66454e-15
 -4.44089e-15   0.0
```
"""
function lyapc(A::AbstractMatrix, E::AbstractMatrix, C::AbstractMatrix)
   """
   The extension of the Bartels-Steward method based on the generalized Schur form
   is employed.

   Reference:
   T. Penzl. Numerical solution of generalized Lyapunov equations.
   Adv. Comput. Math., 8:33–48, 1998.
   """
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && ishermitian(C)) ||
      throw(DimensionMismatch("C must be a symmetric/hermitian matrix of dimension $n"))
   isequal(E,I) && size(E,1) == n && (return lyapc(A, C))
   LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a square matrix of dimension $n"))

   adjA = isa(A,Adjoint)
   adjE = isa(E,Adjoint)
   if adjA && !adjE
      A = copy(A)
      adjA = false
   elseif !adjA && adjE
      E = copy(E)
      adjE = false
   end

   adj = adjA & adjE

   T2 = promote_type(eltype(A), eltype(E), eltype(C))
   T2 <: BlasFloat  || (T2 = promote_type(Float64,T2))
   eltype(A) == T2 || (adj ? A = convert(Matrix{T2},A.parent)' : A = convert(Matrix{T2},A))
   eltype(E) == T2 || (adj ? E = convert(Matrix{T2},E.parent)' : E = convert(Matrix{T2},E))
   eltype(C) == T2 || (C = convert(Matrix{T2},C))

   # Reduce (A,E) to generalized Schur form and transform C
   # (as,es) = (q'*A*z, q'*E*z)
   if adj
      as, es, q, z = schur(A.parent,E.parent)
   else
      as, es, q, z = schur(A,E)
   end
   if adj
      #x = z'*C*z
      x = utqu(C,z)
      lyapcs!(as,es,x,adj = true)
      #x = q*x*q'
      utqu!(x,q')
   else
      #x = q'*C*q
      x = utqu(C,q)
      lyapcs!(as,es,x)
      #x = z*x*z'
      utqu!(x,z')
   end
end
# Aβ'X + XA'β + C = 0
lyapc(A::AbstractMatrix, E::Union{Real,Complex,UniformScaling}, C::AbstractMatrix) =
      isa(A,Adjoint) ? lyapc((A.parent*E')',C) : lyapc(A*E',C)
# αXE' + EXα' + C = 0
lyapc(A::Union{Real,Complex,UniformScaling}, E::AbstractMatrix, C::AbstractMatrix) =
      isa(E,Adjoint) ? lyapc((E.parent*A')',C) : lyapc(A'*E,C)
# (αβ'+α'β)X + C = 0
lyapc(A::UniformScaling, E::UniformScaling, C::AbstractMatrix) = -(A*E'+A'*E)\C
lyapc(A::Union{Real,Complex}, E::Union{Real,Complex}, C::AbstractMatrix) = real(A*E') == 0 ? throw(SingularException(1)) : -C/(A*E'+A'*E)
# (αβ'+α'β)X  + γ = 0
lyapc(A::Union{Real,Complex}, E::Union{Real,Complex}, C::Union{Real,Complex}) = real(A*E') == 0 ? throw(SingularException(1)) : -C/(A*E'+A'*E)
# Discrete Lyapunov equations
"""
    X = lyapd(A, C)

Compute `X`, the symmetric or hermitian solution
of the discrete Lyapunov equation

       AXA' - X + C = 0,

where `A` is a square real or complex matrix and `C` is a symmetric or hermitian
matrix. `A` must not have two eigenvalues `α` and `β` such that `αβ = 1`.
The following particular cases are also adressed:

    X = lyapd(α*I,C)  or  X = lyapd(α,C)

Solve the matrix equation `(αα'-1)X + C = 0`.

    x = lyapd(α,γ)

Solve the equation `(αα'-1)x + γ = 0`.

# Example
```jldoctest
julia> A = [3. 4.; 5. 6.]
2×2 Array{Float64,2}:
 3.0  4.0
 5.0  6.0

julia> C = [1. 1.; 1. 2.]
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  2.0

julia> X = lyapd(A, C)
2×2 Array{Float64,2}:
  0.2375  -0.2125
 -0.2125   0.1375

julia> A*X*A' - X + C
2×2 Array{Float64,2}:
 5.55112e-16  6.66134e-16
 2.22045e-16  4.44089e-16
```
"""
function lyapd(A::AbstractMatrix, C::AbstractMatrix)
   """
   The discrete analog of the Bartels-Steward method based on the Schur form
   is employed.

   Reference:
   G. Kitagawa. An Algorithm for solving the matrix equation X = F X F' + S,
   International Journal of Control, 25:745-753, 1977.
   """
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && ishermitian(C)) ||
       throw(DimensionMismatch("C must be a symmetric/hermitian matrix of dimension $n"))

   adj = isa(A,Adjoint)

   T2 = promote_type(eltype(A), eltype(C))
   T2 <: BlasFloat  || (T2 = promote_type(Float64,T2))
   eltype(A) == T2 || (adj ? A = convert(Matrix{T2},A.parent)' : A = convert(Matrix{T2},A))
   eltype(C) == T2 || (C = convert(Matrix{T2},C))

   # Reduce A to Schur form and transform C
   if adj
      AS, Q = schur(A.parent)
   else
      AS, Q = schur(A)
   end
   #X = Q'*C*Q
   X = utqu(C,Q)
   lyapds!(AS, X, adj = adj)
   #X <- Q*X*Q'
   utqu!(X,Q')
end
# (αα'-1)X + C = 0
lyapd(A::UniformScaling, C::AbstractMatrix) = (I-A'*A)\C
lyapd(A::Union{Real,Complex}, C::AbstractMatrix) = A*A' == 1 ? throw(SingularException(1)) : C/(1-A'*A)
# (αα'-1)x + γ = 0
lyapd(A::Union{Real,Complex}, C::Union{Real,Complex}) = A*A' == 1 ? throw(SingularException(1)) : C/(one(C)-A'*A)
"""
    X = lyapd(A, E, C)

Compute `X`, the symmetric or hermitian solution
of the generalized discrete Lyapunov equation

         AXA' - EXE' + C = 0,

where `A` and `E` are square real or complex matrices and `C` is a symmetric
or hermitian matrix. The pencil `A-λE` must not have two eigenvalues `α` and `β`
such that `αβ = 1`.
The following particular cases are also adressed:

    X = lyapd(A,β*I,C)  or  X = lyapd(A,β,C)

Solve the matrix equation `AXA' - βXβ' + C = 0`.

    X = lyapd(α*I,E,C)  or  X = lyapd(α,E,C)

Solve the matrix equation `αXα' - EXE' + C = 0`.

    X = lyapd(α*I,β*I,C)  or  X = lyapd(α,β,C)

Solve the matrix equation `(αα'-ββ')X + C = 0`.

    x = lyapd(α,β,γ)

Solve the equation `(αα'-ββ')x + γ = 0`.

# Example
```jldoctest
julia> A = [3. 4.; 5. 6.]
2×2 Array{Float64,2}:
 3.0  4.0
 5.0  6.0

julia> E = [ 1. 2.; 0. -1.]
2×2 Array{Float64,2}:
 1.0   2.0
 0.0  -1.0

julia> C = [1. 1.; 1. 2.]
2×2 Array{Float64,2}:
 1.0  1.0
 1.0  2.0

julia> X = lyapd(A, E, C)
2×2 Array{Float64,2}:
  1.775  -1.225
 -1.225   0.775

julia> A*X*A' - E*X*E' + C
2×2 Array{Float64,2}:
 -2.22045e-16  -4.44089e-16
 -1.33227e-15   1.11022e-15
```
"""
function lyapd(A::AbstractMatrix, E::AbstractMatrix, C::AbstractMatrix)
   """
   The extension of the Bartels-Steward method based on the generalized Schur form
   is employed.

   Reference:
   T. Penzl. Numerical solution of generalized Lyapunov equations.
   Adv. Comput. Math., 8:33–48, 1998.
   """
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && ishermitian(C)) ||
      throw(DimensionMismatch("C must be a symmetric/hermitian matrix of dimension $n"))
   isequal(E,I) && size(E,1) == n && (return lyapd(A, C))
   LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a square matrix of dimension $n"))

   adjA = isa(A,Adjoint)
   adjE = isa(E,Adjoint)
   if adjA && !adjE
      A = copy(A)
      adjA = false
   elseif !adjA && adjE
      E = copy(E)
      adjE = false
   end

   adj = adjA & adjE

   T2 = promote_type(eltype(A), eltype(E), eltype(C))
   T2 <: BlasFloat  || (T2 = promote_type(Float64,T2))
   eltype(A) == T2 || (adj ? A = convert(Matrix{T2},A.parent)' : A = convert(Matrix{T2},A))
   eltype(E) == T2 || (adj ? E = convert(Matrix{T2},E.parent)' : E = convert(Matrix{T2},E))
   eltype(C) == T2 || (C = convert(Matrix{T2},C))

   # Reduce (A,E) to generalized Schur form and transform C
   # (as,es) = (q'*A*z, q'*E*z)
   if adj
      as, es, q, z = schur(A.parent,E.parent)
   else
      as, es, q, z = schur(A,E)
   end
   if adj
      #x = z'*C*z
      x = utqu(C,z)
      lyapds!(as,es,x,adj = true)
      #x = q*x*q'
      utqu!(x,q')
   else
      #x = q'*C*q
      x = utqu(C,q)
      lyapds!(as,es,x)
      #x = z*x*z'
      utqu!(x,z')
   end
end
# AXA' - Xββ' + C = 0
lyapd(A::AbstractMatrix, E::Union{Real,Complex,UniformScaling}, C::AbstractMatrix) =
      isa(A,Adjoint) ? lyapd((A.parent/E)',C/(E*E')) : lyapd(A/E,C/(E*E'))
# αXα' - EXE' + C = 0
lyapd(A::Union{Real,Complex,UniformScaling}, E::AbstractMatrix, C::AbstractMatrix) =
      isa(E,Adjoint) ? lyapd((E.parent/A)',C/(-A*A')) : lyapd(E/A,C/(-A*A'))
# (α'α-β'β)X + C = 0
lyapd(A::UniformScaling, E::UniformScaling, C::AbstractMatrix) = C/(E'*E-A'*A)
lyapd(A::Union{Real,Complex}, E::Union{Real,Complex}, C::AbstractMatrix) = A*A'== E*E' ? throw(SingularException(1)) : C/(E'*E-A'*A)
# (α'α-β'β)X + γ  = 0
lyapd(A::Union{Real,Complex}, E::Union{Real,Complex}, C::Union{Real,Complex}) = A*A'== E*E' ? throw(SingularException(1)) : C/(E'*E-A'*A)
"""
    lyapcs!(A,C;adj = false)

Solve the continuous Lyapunov matrix equation

                op(A)X + Xop(A)' + C = 0,

where `op(A) = A` if `adj = false` and `op(A) = A'` if `adj = true`.
`A` is a square real matrix in a real Schur form, or a square complex matrix in a
complex Schur form and `C` is a symmetric or hermitian matrix.
`A` must not have two eigenvalues `α` and `β` such that `α+β = 0`.
`C` contains on output the solution `X`.
"""
function lyapcs!(A::T1, C::Union{T1,T2}; adj = false) where
   {T1<:Union{Matrix{Float32},Matrix{Float64}}, T2<:Union{Matrix{Complex{Float64}},Matrix{Complex{Float32}}} }
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && ishermitian(C)) ||
      throw(DimensionMismatch("C must be a $n x $n symmetric/hermitian matrix"))

   T = eltype(A)
   TR = real(eltype(C))
   T == TR || error("TypeError: for real part of C expected Type{$T}, got Type{$TR}")
   ZERO = zero(T)
   # determine the structure of the real Schur form
   ba = fill(1,n)
   p = 1
   if n > 1
      d = [diag(A,-1);zeros(1)]
      i = 1
      p = 0
      while i <= n
         p += 1
         if d[i] != ZERO
            ba[p] = 2
            i += 1
         end
         i += 1
      end
   end

   W = Array{eltype(A),2}(I,2,2)
   if adj
      """
      The (K,L)th block of X is determined starting from
      upper-left corner column by column by

      A(K,K)'*X(K,L) + X(K,L)*A(L,L) = -C(K,L) - R(K,L),

      where
               K-1                    L-1
      R(K,L) = SUM [A(I,K)'*X(I,L)] + SUM [X(K,J)*A(J,L)].
               I=1                    J=1
      """
      j = 1
      for ll = 1:p
          dl = ba[ll]
          l = j:j+dl-1
          i = j
          for kk = ll:p
              dk = ba[kk]
              k = i:i+dk-1
              y = C[k,l]
              if kk > ll
                 ia = j:i-1
                 y += A[ia,k]'*C[ia,l]
              end
              Z = (kron(W[1:dl,1:dl],transpose(A[k,k]))+kron(transpose(A[l,l]),W[1:dk,1:dk]))\(-y[:])
              isfinite(maximum(abs.(Z))) ? C[k,l] = Z : throw("ME:SingularException: A has eigenvalues α and β such that α+β ≈ 0")
              if i == j && dk == 2
                 temp = C[k,l]
                 C[k,l] = (temp'+temp)/2
              else
                 C[l,k] = C[k,l]'
              end
              i += dk
           end
           j += dl
           if j <= n
              for jr = j:n
                 for ir = jr:n
                     for lll = l
                         C[ir,jr] += C[ir,lll]*A[lll,jr] + (A[lll,ir]*C[jr,lll])'
                         C[jr,ir] = C[ir,jr]'
                      end
                  end
               end
           end
      end
   else
      """
      The (K,L)th block of X is determined starting from
      bottom-right corner column by column by

      A(K,K)*X(K,L) + X(K,L)*A(L,L)' = -C(K,L) - R(K,L),

      where
                N                     N
      R(K,L) = SUM [A(K,I)*X(I,L)] + SUM [X(K,J)*A(L,J)'].
              I=K+1                 J=L+1
      """
      j = n
      for ll = p:-1:1
          dl = ba[ll]
          l = j-dl+1:j
          i = j
          for kk = ll:-1:1
              dk = ba[kk]
              i1 = i-dk+1
              k = i1:i
              y = C[k,l]
              if kk < ll
                 ia = i+1:j
                 y += A[k,ia]*C[ia,l]
              end
              Z = (kron(W[1:dl,1:dl],A[k,k])+kron(A[l,l],W[1:dk,1:dk]))\(-y[:])
              isfinite(maximum(abs.(Z))) ? C[k,l] = Z : throw("ME:SingularException: A has eigenvalues α and β such that α+β ≈ 0")
              if i == j && dl == 2
                 temp = C[k,l]
                 C[k,l] = (temp'+temp)/2
              else
                 C[l,k] = C[k,l]'
               end
              i = i-dk
           end
           j = j-dl
           if j >= 0
              for jr = 1:j
                 for ir = 1:jr
                     for lll = l
                         C[ir,jr] += C[ir,lll]*A[jr,lll] + A[ir,lll]*C[lll,jr]
                         C[jr,ir] = C[ir,jr]'
                      end
                  end
               end
           end
       end
   end
end
function lyapcs!(A::T1, C::T1; adj = false) where T1<:Union{Matrix{Complex{Float64}},Matrix{Complex{Float32}}}
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && ishermitian(C)) ||
      throw(DimensionMismatch("C must be a $n x $n symmetric/hermitian matrix"))

   if adj
      """
      The (K,L)th element of X is determined starting from
      upper-left corner column by column by

      A(K,K)'*X(K,L) + X(K,L)*A(L,L) = -C(K,L) - R(K,L),

      where
               K-1                    L-1
      R(K,L) = SUM [A(I,K)'*X(I,L)] + SUM [X(K,J)*A(J,L)].
               I=1                    J=1
      """
      for l = 1:n
          for k = l:n
              y = C[k,l]
              for ia = l:k-1
                  y += A[ia,k]'*C[ia,l]
              end
              Z = -y/(A[k,k]'+A[l,l])
              isfinite(Z) ? C[k,l] = Z : throw("ME:SingularException: A has eigenvalues α and β such that α+β ≈ 0")
              if k != l
                 C[l,k] = C[k,l]'
              end
           end
           for jr = l+1:n
               for ir = jr:n
                   C[ir,jr] += C[ir,l]*A[l,jr] + (A[l,ir]*C[jr,l])'
                   if jr != ir
                      C[jr,ir] = C[ir,jr]'
                   end
               end
           end
      end
   else
      """
      The (K,L)th element of X is determined starting from
      bottom-right corner column by column by

      A(K,K)*X(K,L) + X(K,L)*A(L,L)' = -C(K,L) - R(K,L),

      where
                N                     N
      R(K,L) = SUM [A(K,I)*X(I,L)] + SUM [X(K,J)*A(L,J)'].
              I=K+1                 J=L+1
      """
      for l = n:-1:1
          for k = l:-1:1
              y = C[k,l]
              for ia = k+1:l
                 y += A[k,ia]*C[ia,l]
              end
              Z = -y/(A[k,k]+A[l,l]')
              isfinite(Z) ? C[k,l] = Z : throw("ME:SingularException: A has eigenvalues α and β such that α+β ≈ 0")
              if k != l
                 C[l,k] = C[k,l]'
               end
           end
           for jr = 1:l-1
               for ir = 1:jr
                  C[ir,jr] += C[ir,l]*A[jr,l]' + A[ir,l]*C[l,jr]
                  if ir != jr
                     C[jr,ir] = C[ir,jr]'
                  end
               end
           end
       end
   end
end
"""
    lyapcs!(A, E, C; adj = false)

Solve the generalized continuous Lyapunov matrix equation

                op(A)Xop(E)' + op(E)Xop(A)' + C = 0,

where `op(A) = A` and `op(E) = E` if `adj = false` and `op(A) = A'` and
`op(E) = E'` if `adj = true`. The pair `(A,E)` is in a generalized real or
complex Schur form and `C` is a symmetric or hermitian matrix.
The pencil `A-λE` must not have two eigenvalues `α` and `β` such that `α+β = 0`.
The computed symmetric or hermitian solution `X` is contained in `C`.
"""
function lyapcs!(A::T1, E::Union{T1,UniformScaling{Bool}}, C::Union{T1,T2}; adj = false) where
      {T1<:Union{Matrix{Float32},Matrix{Float64}}, T2<:Union{Matrix{Complex{Float64}},Matrix{Complex{Float32}}} }
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && ishermitian(C)) ||
      throw(DimensionMismatch("C must be a $n x $n symmetric/hermitian matrix"))
   (typeof(E) == UniformScaling{Bool} || (isequal(E,I) && size(E,1) == n)) && (lyapcs!(A, C, adj = adj); return)
   LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a $n x $n matrix or I"))

   T = eltype(A)
   TR = real(eltype(C))
   T == TR || error("TypeError: for real part of C expected Type{$T}, got Type{$TR}")

   # determine the structure of the generalized real Schur form
   ba = fill(1,n)
   p = 1
   if n > 1
      d = [diag(A,-1);zeros(1)]
      i = 1
      p = 0
      while i <= n
         p += 1
         if d[i] != 0
            ba[p] = 2
            i += 1
         end
         i += 1
      end
   end

   W = Array{eltype(C),2}(undef,n,2)
   if adj
      """
      The (K,L)th block of X is determined starting from the
      upper-left corner column by column by

      A(K,K)'*X(K,L)*E(L,L) + E(K,K)'*X(K,L)*A(L,L) = -C(K,L) - R(K,L),

      where
                K           L-1
      R(K,L) = SUM {A(I,K)'*SUM [X(I,J)*E(J,L)]} +
               I=1          J=1

                K           L-1
               SUM {E(I,K)'*SUM [X(I,J)*A(J,L)]} +
               I=1          J=1

               K-1
              {SUM [A(I,K)'*X(I,L)]}*E(L,L) +
               I=1

               K-1
              {SUM [E(I,K)'*X(I,L)]}*A(L,L).
               I=1
      """
      i = 1
      for kk = 1:p
          dk = ba[kk]
          dkk = 1:dk
          k = i:i+dk-1
          j = 1
          for ll = 1:kk
             j1 = j+ba[ll]-1
             l = j:j1
             y = C[k,l]
             if kk > 1
                 ir = 1:i-1
                 C[l,k] = C[l,ir]*A[ir,k]
                 W[l,dkk] = C[l,ir]*E[ir,k]
                 ic = 1:j1
                 y += C[ic,k]'*E[ic,l] + W[ic,dkk]'*A[ic,l]
             end
             Z = ((kron(E[l,l],A[k,k])+kron(A[l,l],E[k,k]))')\(-y[:])
             isfinite(maximum(abs.(Z))) ? C[k,l] = Z : throw("ME:SingularException: A-λE has eigenvalues α and β such that α+β ≈ 0")
             if i == j
                if dk == 2
                   temp = C[k,l]
                   C[k,l] = Hermitian((temp'+temp)/2)
                else
                   C[k,l] = real(C[k,l])
                end
             end
             j += ba[ll]
             if j <= i
                C[l,k] += C[k,l]'*A[k,k]
                W[l,dkk] += C[k,l]'*E[k,k]
             end
          end
          if kk > 1
             ir = 1:i-1
             C[ir,k] = C[k,ir]'
          end
          i += dk
      end
   else
      """
      The (K,L)th block of X is determined starting from
      bottom-right corner column by column by

      A(K,K)*X(K,L)*E(L,L)' + E(K,K)*X(K,L)*A(L,L)' = -C(K,L) - R(K,L),

      where

                N            N
      R(K,L) = SUM {A(K,I)* SUM [X(I,J)*E(L,J)']} +
               I=K         J=L+1

                N            N
               SUM {E(K,I)* SUM [X(I,J)*A(L,J)']} +
               I=K         J=L+1

                  N
               { SUM [A(K,J)*X(J,L)]}*E(L,L)' +
                J=K+1

                  N
               { SUM [E(K,J)*X(J,L)]}*A(L,L)'
                J=K+1
      """
      j = n
      for ll = p:-1:1
        dl = ba[ll]
        l = j-dl+1:j
        dll = 1:dl
        i = n
        for kk = p:-1:ll
            dk = ba[kk]
            i1 = i-dk+1
            k = i1:i
            y = C[l,k]
            if ll < p
               ir = j+1:n
               C[k,l] = C[k,ir]*A[l,ir]'
               W[k,dll] = C[k,ir]*E[l,ir]'
               ic = i1:n
               y += (E[k,ic]*C[ic,l] + A[k,ic]*W[ic,dll])'
            end
            Z = (kron(E[k,k],A[l,l])+kron(A[k,k],E[l,l]))\(-y[:])
            isfinite(maximum(abs.(Z))) ? C[l,k] = Z : throw("ME:SingularException: A-λE has eigenvalues α and β such that α+β ≈ 0")
            if i == j
               if dk == 2
                  temp = C[l,k]
                  C[l,k] = Hermitian((temp'+temp)/2)
               else
                  C[l,k] = real(C[l,k])
               end
            end
            i = i-dk
            if i >= j
               C[k,l] += (A[l,l]*C[l,k])'
               W[k,dll] += (E[l,l]*C[l,k])'
            else
               break
            end
        end
        if ll < p
           ir = j+1:n
           C[ir,l] = C[l,ir]'
        end
        j = j-dl
      end
   end
end
function lyapcs!(A::T1, E::Union{T1,UniformScaling{Bool}}, C::T1; adj = false) where T1<:Union{Matrix{Complex{Float64}},Matrix{Complex{Float32}}}
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && ishermitian(C)) ||
      throw(DimensionMismatch("C must be a $n x $n symmetric/hermitian matrix"))
   (typeof(E) == UniformScaling{Bool} || (isequal(E,I) && size(E,1) == n)) && (lyapcs!(A, C, adj = adj); return)
   LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a $n x $n matrix or I"))

   W = Array{Complex{Float64},1}(undef,n)
   # Compute the hermitian solution
   if adj
      for k = 1:n
         for l = 1:k
            y = C[k,l]
            if k > 1
               C[l,k] = C[l,1]*A[1,k]
               W[l] = C[l,1]*E[1,k]
               for ir = 2:k-1
                  C[l,k] +=  C[l,ir]*A[ir,k]
                  W[l] += C[l,ir]*E[ir,k]
               end
               for ic = 1:l
                   y += C[ic,k]'*E[ic,l] + W[ic]'*A[ic,l]
                end
            end
            Z = -y/(A[k,k]'*E[l,l]+E[k,k]'*A[l,l])
            isfinite(Z) ? C[k,l] = Z : throw("ME:SingularException: A-λE has eigenvalues α and β such that α+β ≈ 0")
            if k == l
               C[k,l] = real(C[k,l])
            end
            if l < k
               C[l,k] += C[k,l]'*A[k,k]
               W[l] += C[k,l]'*E[k,k]
            end
         end
         for ir = 1:k-1
             C[ir,k] = C[k,ir]'
         end
      end
   else
      for l = n:-1:1
        for k = n:-1:l
            y = C[l,k]
            if l < n
               C[k,l] = C[k,l+1]*A[l,l+1]'
               W[k] = C[k,l+1]*E[l,l+1]'
               for ir = l+2:n
                  C[k,l] += C[k,ir]*A[l,ir]'
                  W[k] += C[k,ir]*E[l,ir]'
               end
               for ic = k:n
                   y += (E[k,ic]*C[ic,l] + A[k,ic]*W[ic])'
               end
            end
            Z = -y/(E[k,k]'*A[l,l]+A[k,k]'*E[l,l])
            isfinite(Z) ? C[l,k] = Z : throw("ME:SingularException: A-λE has eigenvalues α and β such that α+β ≈ 0")
            if k == l
               C[k,l] = real(C[k,l])
            end
            if k > l
               C[k,l] += (A[l,l]*C[l,k])'
               W[k] += (E[l,l]*C[l,k])'
            end
        end
        if l < n
           for ir = l+1:n
               C[ir,l] = C[l,ir]'
           end
        end
      end
   end
end

"""
    lyapds!(A, C; adj = false)

Solve the discrete Lyapunov matrix equation

                op(A)Xop(A)' - X + C = 0,

where `op(A) = A` if `adj = false` and `op(A) = A'` if `adj = true`.
`A` is in a real or complex Schur form and `C` is a symmetric or hermitian matrix.
`A` must not have two eigenvalues `α` and `β` such that `αβ = 1`.
The computed symmetric or hermitian solution `X` is contained in `C`.
"""
function lyapds!(A::T1, C::Union{T1,T2}; adj = false) where
   {T1<:Union{Matrix{Float32},Matrix{Float64}}, T2<:Union{Matrix{Complex{Float64}},Matrix{Complex{Float32}}} }
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && ishermitian(C)) ||
      throw(DimensionMismatch("C must be a $n x $n symmetric/hermitian matrix"))

   T = eltype(A)
   TR = real(eltype(C))
   T == TR || error("TypeError: for real part of C expected Type{$T}, got Type{$TR}")

   # determine the structure of the real Schur form
   ba = fill(1,n)
   p = 1
   if n > 1
      d = [diag(A,-1);zeros(1)]
      i = 1
      p = 0
      while i <= n
         p += 1
         if d[i] != 0
            ba[p] = 2
            i += 1
         end
         i += 1
      end
   end

   if adj
      """
      The (K,L)th block of X is determined starting from the
      upper-left corner column by column by

      A(K,K)'*X(K,L)*A(L,L) - X(K,L) = -C(K,L) - R(K,L),

      where
                K           L-1
      R(K,L) = SUM {A(I,K)'*SUM [X(I,J)*A(J,L)]} +
               I=1          J=1

                K-1
               {SUM [A(I,K)'*X(I,L)]}*A(L,L).
                I=1
      """
      i = 1
      for kk = 1:p
          dk = ba[kk]
          k = i:i+dk-1
          j = 1
          for ll = 1:kk
              j1 = j+ba[ll]-1
              l = j:j1
              y = C[k,l]
              if kk > 1
                 ir = 1:i-1
                 C[l,k] = C[l,ir]*A[ir,k]
                 ic = 1:j1
                 y += C[ic,k]'*A[ic,l]
              end
              Z = (I-kron(A[l,l]',A[k,k]'))\y[:]
              isfinite(maximum(abs.(Z))) ? C[k,l] = Z : throw("ME:SingularException: A has eigenvalues α and β such that αβ ≈ 1")
              if i == j
                 if dk == 2
                    temp = C[k,l]
                    C[k,l] = Hermitian((temp'+temp)/2)
                 else
                    C[k,l] = real(C[k,l])
                 end
              end
              j += ba[ll]
              if j <= i
                 C[l,k] += C[k,l]'*A[k,k]
              end
          end
          if kk > 1
             ir = 1:i-1
             C[ir,k] = C[k,ir]'
          end
          i += dk
      end
   else
      """
      The (K,L)th block of X is determined starting from
      bottom-right corner column by column by

      A(K,K)*X(K,L)*A(L,L)' - X(K,L) = -C(K,L) - R(K,L),

      where

                N            N
      R(K,L) = SUM {A(K,I)* SUM [X(I,J)*A(L,J)']} +
               I=K         J=L+1

                  N
               { SUM [A(K,J)*X(J,L)]}*A(L,L)'
                J=K+1
      """
      j = n
      for ll = p:-1:1
          dl = ba[ll]
          l = j-dl+1:j
          i = n
          for kk = p:-1:ll
              i1 = i-ba[kk]+1
              k = i1:i
              y = C[l,k]
              if ll < p
                 ir = j+1:n
                 C[k,l] = C[k,ir]*A[l,ir]'
                 ic = i1:n
                 y += (A[k,ic]*C[ic,l])'
              end
              Z = (I-kron(A[k,k],A[l,l]))\y[:]
              isfinite(maximum(abs.(Z))) ? C[l,k] = Z : throw("ME:SingularException: A has eigenvalues α and β such that αβ ≈ 1")
              if i == j
                 if dl == 2
                    temp = C[l,k]
                    C[l,k] = Hermitian((temp'+temp)/2)
                 else
                    C[l,k] = real(C[l,k])
                 end
              end
              i = i-ba[kk]
              if i >= j
                 C[k,l] += (A[l,l]*C[l,k])'
              else
                 break
              end
          end
          if ll < p
             ir = i+2:n
             C[ir,l] = C[l,ir]'
          end
          j = j-dl
      end
   end
end

function lyapds!(A::T1, C::T1; adj = false) where T1<:Union{Matrix{Complex{Float64}},Matrix{Complex{Float32}}}
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && ishermitian(C)) ||
      throw(DimensionMismatch("C must be a $n x $n symmetric/hermitian matrix"))

   # Compute the hermitian solution
   if adj
      for k = 1:n
         for l = 1:k
            y = C[k,l]
            if k > 1
               C[l,k] = C[l,1]*A[1,k]
               for ir = 2:k-1
                  C[l,k] +=  C[l,ir]*A[ir,k]
               end
               for ic = 1:l
                   y += C[ic,k]'*A[ic,l]
                end
            end
            Z = y/(I-A[k,k]'*A[l,l])
            isfinite(Z) ? C[k,l] = Z : throw("ME:SingularException: A has eigenvalues α and β such that αβ ≈ 1")
            if k == l
               C[k,l] = real(C[k,l])
            end
            if l < k
               C[l,k] += C[k,l]'*A[k,k]
            end
         end
         for ir = 1:k-1
             C[ir,k] = C[k,ir]'
         end
      end
   else
      for l = n:-1:1
        for k = n:-1:l
            y = C[l,k]
            if l < n
               C[k,l] = C[k,l+1]*A[l,l+1]'
               for ir = l+2:n
                  C[k,l] += C[k,ir]*A[l,ir]'
               end
               for ic = k:n
                   y += (A[k,ic]*C[ic,l])'
               end
            end
            Z = y/(I-A[k,k]'*A[l,l])
            isfinite(Z) ? C[l,k] = Z : throw("ME:SingularException: A has eigenvalues α and β such that αβ ≈ 1")
            if k == l
               C[k,l] = real(C[k,l])
            end
            if k > l
               C[k,l] += (A[l,l]*C[l,k])'
            end
        end
        if l < n
           for ir = l+1:n
               C[ir,l] = C[l,ir]'
           end
        end
      end
   end
end
"""
    lyapds!(A, E, C; adj = false)

Solve the generalized discrete Lyapunov matrix equation

                op(A)Xop(A)' - op(E)Xop(E)' + C = 0,

where `op(A) = A` and `op(E) = E` if `adj = false` and `op(A) = A'` and
`op(E) = E'` if `adj = true`. The pair `(A,E)` is in a generalized real or
complex Schur form and `C` is a symmetric or hermitian matrix.
The pencil `A-λE` must not have two eigenvalues `α` and `β` such that `αβ = 1`.
The computed symmetric or hermitian solution `X` is contained in `C`.
"""
function lyapds!(A::T1, E::Union{T1,UniformScaling{Bool}}, C::Union{T1,T2}; adj = false) where
   {T1<:Union{Matrix{Float32},Matrix{Float64}}, T2<:Union{Matrix{Complex{Float64}},Matrix{Complex{Float32}}} }
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && ishermitian(C)) ||
      throw(DimensionMismatch("C must be a $n x $n symmetric/hermitian matrix"))
   (typeof(E) == UniformScaling{Bool} || (isequal(E,I) && size(E,1) == n)) && (lyapds!(A, C, adj = adj); return)
   LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a $n x $n matrix or I"))

   T = eltype(A)
   TR = real(eltype(C))
   T == TR || error("TypeError: for real part of C expected Type{$T}, got Type{$TR}")

   # determine the structure of the real Schur form
   ba = fill(1,n)
   p = 1
   if n > 1
      d = [diag(A,-1);zeros(1)]
      i = 1
      p = 0
      while i <= n
         p += 1
         if d[i] != 0
            ba[p] = 2
            i += 1
         end
         i += 1
      end
   end

   W = Array{eltype(C),2}(undef,n,2)
   if adj
      """
      The (K,L)th block of X is determined starting from the
      upper-left corner column by column by

      A(K,K)'*X(K,L)*A(L,L) - E(K,K)'*X(K,L)*E(L,L) = -C(K,L) - R(K,L),

      where
                K           L-1
      R(K,L) = SUM {A(I,K)'*SUM [X(I,J)*A(J,L)]} -
               I=1          J=1

                K           L-1
               SUM {E(I,K)'*SUM [X(I,J)*E(J,L)]} +
               I=1          J=1

                K-1
               {SUM [A(I,K)'*X(I,L)]}*A(L,L) -
                I=1

                K-1
               {SUM [E(I,K)'*X(I,L)]}*E(L,L).
                I=1
      """
      i = 1
      for kk = 1:p
          dk = ba[kk]
          dkk = 1:dk
          k = i:i+dk-1
          j = 1
          for ll = 1:kk
             j1 = j+ba[ll]-1
             l = j:j1
             y = C[k,l]
             if kk > 1
                 ir = 1:i-1
                 C[l,k] = C[l,ir]*A[ir,k]
                 W[l,dkk] = C[l,ir]*E[ir,k]
                 ic = 1:j1
                 y += C[ic,k]'*A[ic,l] -W[ic,dkk]'*E[ic,l]
             end
             Z = (kron(E[l,l]',E[k,k]')-kron(A[l,l]',A[k,k]'))\y[:]
             isfinite(maximum(abs.(Z))) ? C[k,l] = Z : throw("ME:SingularException: A-λE has eigenvalues α and β such that αβ ≈ 1")
             if i == j
                if dk == 2
                   temp = C[k,l]
                   C[k,l] = Hermitian((temp'+temp)/2)
                else
                   C[k,l] = real(C[k,l])
                end
             end
             j += ba[ll]
             if j <= i
                 C[l,k] += C[k,l]'*A[k,k]
                 W[l,dkk] += C[k,l]'*E[k,k]
             end
          end
          if kk > 1
             ir = 1:i-1
             C[ir,k] = C[k,ir]'
          end
          i += dk
      end
   else
      """
      The (K,L)th block of X is determined starting from
      bottom-right corner column by column by

      A(K,K)*X(K,L)*A(L,L)' - E(K,K)*X(K,L)*E(L,L)' = -C(K,L) - R(K,L),

      where

                N            N
      R(K,L) = SUM {A(K,I)* SUM [X(I,J)*A(L,J)']} -
               I=K         J=L+1

                N            N
               SUM {E(K,I)* SUM [X(I,J)*E(L,J)']} +
               I=K         J=L+1

                  N
               { SUM [A(K,J)*X(J,L)]}*A(L,L)' -
                J=K+1

                  N
               { SUM [E(K,J)*X(J,L)]}*E(L,L)'
                J=K+1
   """
      j = n
      for ll = p:-1:1
        dl = ba[ll]
        l = j-dl+1:j
        dll = 1:dl
        i = n
        for kk = p:-1:ll
            i1 = i-ba[kk]+1
            k = i1:i
            y = C[l,k]
            if ll < p
               ir = j+1:n
               C[k,l] = C[k,ir]*A[l,ir]'
               W[k,dll] = C[k,ir]*E[l,ir]'
               ic = i1:n
               y += (A[k,ic]*C[ic,l] - E[k,ic]*W[ic,dll])'
            end
            Z = (kron(E[k,k],E[l,l])-kron(A[k,k],A[l,l]))\y[:]
            isfinite(maximum(abs.(Z))) ? C[l,k] = Z : throw("ME:SingularException: A-λE has eigenvalues α and β such that αβ ≈ 1")
            if i == j
               if dl == 2
                  temp = C[l,k]
                  C[l,k] = Hermitian((temp'+temp)/2)
               else
                  C[l,k] = real(C[l,k])
               end
            end
            i = i-ba[kk]
            if i >= j
               C[k,l] += (A[l,l]*C[l,k])'
               W[k,dll] += (E[l,l]*C[l,k])'
            else
               break
            end
        end
        if ll < p
           ir = i+2:n
           C[ir,l] = C[l,ir]'
        end
        j = j-dl
      end
   end
end

function lyapds!(A::T1, E::Union{T1,UniformScaling{Bool}}, C::T1; adj = false) where T1<:Union{Matrix{Complex{Float64}},Matrix{Complex{Float32}}}
   n = LinearAlgebra.checksquare(A)
   (LinearAlgebra.checksquare(C) == n && ishermitian(C)) ||
      throw(DimensionMismatch("C must be a $n x $n hermitian matrix"))
   (typeof(E) == UniformScaling{Bool} || (isequal(E,I) && size(E,1) == n)) && (lyapds!(A, C, adj = adj); return)
   LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a $n x $n matrix or I"))

   W = Array{Complex{Float64},1}(undef,n)
   # Compute the hermitian solution
   if adj
      for k = 1:n
         for l = 1:k
            y = C[k,l]
            if k > 1
               C[l,k] = C[l,1]*A[1,k]
               W[l] = C[l,1]*E[1,k]
               for ir = 2:k-1
                  C[l,k] +=  C[l,ir]*A[ir,k]
                  W[l] += C[l,ir]*E[ir,k]
               end
               for ic = 1:l
                   y += C[ic,k]'*A[ic,l] - W[ic]'*E[ic,l]
                end
            end
            Z = y/(E[k,k]'*E[l,l]-A[k,k]'*A[l,l])
            isfinite(Z) ? C[k,l] = Z : throw("ME:SingularException: A-λE has eigenvalues α and β such that αβ ≈ 1")
            if k == l
               C[k,l] = real(C[k,l])
            end
            if l < k
               C[l,k] += C[k,l]'*A[k,k]
               W[l] += C[k,l]'*E[k,k]
            end
         end
         for ir = 1:k-1
             C[ir,k] = C[k,ir]'
         end
      end
   else
      for l = n:-1:1
        for k = n:-1:l
            y = C[l,k]
            if l < n
               C[k,l] = C[k,l+1]*A[l,l+1]'
               W[k] = C[k,l+1]*E[l,l+1]'
               for ir = l+2:n
                  C[k,l] += C[k,ir]*A[l,ir]'
                  W[k] += C[k,ir]*E[l,ir]'
               end
               for ic = k:n
                   y += (A[k,ic]*C[ic,l] - E[k,ic]*W[ic])'
               end
            end
            Z = y/(E[k,k]'*E[l,l]-A[k,k]'*A[l,l])
            isfinite(Z) ? C[l,k] = Z : throw("ME:SingularException: A-λE has eigenvalues α and β such that αβ ≈ 1")
            if k == l
               C[k,l] = real(C[k,l])
            end
            if k > l
               C[k,l] += (A[l,l]*C[l,k])'
               W[k] += (E[l,l]*C[l,k])'
            end
        end
        if l < n
           for ir = l+1:n
               C[ir,l] = C[l,ir]'
           end
        end
      end
   end
end


text:
using Documenter, WaveFD

makedocs(sitename = "WaveFD", modules=[WaveFD])

deploydocs(
    repo = "github.com/ChevronETC/WaveFD.jl.git",
)

ext: jl

max_stars_repo_name: ChevronETC/WaveFD

max_issues_repo_name: ChevronETC/WaveFD

max_forks_repo_name: ChevronETC/WaveFD


text:
"""
This file has been adapted from an example found at https://github.com/GenXProject/GenX
"""

cd(dirname(@__FILE__))
settings_path = joinpath(pwd(), "Settings")
#=
environment_path = "../../../package_activate.jl"
include(environment_path) #Run this line to activate the Julia virtual environment for GenX; skip it, if the appropriate package versions are installed
=#
### Set relevant directory paths
src_path = "$GENX_REPO_PATH"

inpath = pwd()

### Load GenX
println("Loading packages")
push!(LOAD_PATH, src_path)
using GenX
using YAML

genx_settings = joinpath(settings_path, "genx_settings.yml") #Settings YAML file path
mysetup = YAML.load(open(genx_settings)) # mysetup dictionary stores settings and GenX-specific parameters

### Cluster time series inputs if necessary and if specified by the user
TDRpath = joinpath(inpath, mysetup["TimeDomainReductionFolder"])
if mysetup["TimeDomainReduction"] == 1
    if (!isfile(TDRpath*"/Load_data.csv")) || (!isfile(TDRpath*"/Generators_variability.csv")) || (!isfile(TDRpath*"/Fuels_data.csv"))
        println("Clustering Time Series Data...")
        cluster_inputs(inpath, settings_path, mysetup)
    else
        println("Time Series Data Already Clustered.")
    end
end

### Configure solver
println("Configuring Solver")
OPTIMIZER = configure_solver(mysetup["Solver"], settings_path)

#### Running a case

### Load inputs
println("Loading Inputs")
myinputs = Dict() # myinputs dictionary will store read-in data and computed parameters
myinputs = load_inputs(mysetup, inpath)

### Generate model
println("Generating the Optimization Model")
EP = generate_model(mysetup, myinputs, OPTIMIZER)

### Solve model
println("Solving Model")
EP, solve_time = solve_model(EP, mysetup)
myinputs["solve_time"] = solve_time # Store the model solve time in myinputs

### Write output
# Run MGA if the MGA flag is set to 1 else only save the least cost solution
println("Writing Output")
outpath = "$inpath/Results"
write_outputs(EP, outpath, mysetup, myinputs)
if mysetup["ModelingToGenerateAlternatives"] == 1
    println("Starting Model to Generate Alternatives (MGA) Iterations")
    mga(EP, inpath, mysetup, myinputs, outpath)
end

max_stars_repo_name: danielolsen/pygenx

max_issues_repo_name: danielolsen/pygenx

max_forks_repo_name: danielolsen/pygenx


text:
using Documenter, Air

makedocs(;
         modules=[Air],
         format=Documenter.HTML(),
         pages=[
             "Home" => "index.md",
             "Persistent Arrays" => "parray.md",
             "Persistent Dictionaries" => "pdict.md",
             "Persistent Sets" => "pset.md",
             "Persistent Weighted Dictionaries" => "pwdict.md",
             "Persistent Weighted Sets" => "pwset.md",
             "Persistent Lazy Dictionaries" => "lazydict.md"],
         repo="https://github.com/noahbenson/Air.jl/blob/{commit}{path}#L{line}",
         sitename="Air.jl",
         authors="Noah C. Benson",
         #assets=String[],
)

deploydocs(;
           repo="github.com/noahbenson/Air.jl",
           deploy_config=Documenter.GitHubActions(),
           push_preview=false
)


text:
using LightOSM, GeoInterface, Plots, DataFrames
using ArchGDAL: createmultilinestring

g = graph_from_download(
    :place_name, place_name="moabit, berlin germany",
    network_type=:bike
)
# reverse coordinates for plotting
reverse!.(g.node_coordinates)

g_simple, weights, node_gdf, edge_gdf = simplify_graph(g)

# join edges in mulitlinestring for faster plotting
all_edges = createmultilinestring(coordinates.(edge_gdf.geom))

# node validation

# nodes from original graph
plot(all_edges, color=:black, size=(1200,800))
scatter!(first.(g.node_coordinates), last.(g.node_coordinates), color=:red)

# nodes from simplified graph
plot(all_edges, color=:black, size=(1200,800))
scatter!(node_gdf.geom, color=:green)

# edge validation

function highway_gdf(osmg::OSMGraph)
    function _geometrize_way(way)
        createlinestring(map(id -> coordinates(osmg.nodes[id]), way.nodes))
    end
    geom = map(way -> _geometrize_way(way), values(osmg.highways))
    return DataFrame(; id = collect(keys(osmg.highways)), geom)
end


text:
module FiVI

using POMDPs
using POMDPPolicies
using POMDPModelTools
using LinearAlgebra
using FiniteHorizonPOMDPs
using Distributions
import POMDPs: Solver, solve
import Base: ==, hash, convert, getindex
import FiniteHorizonPOMDPs: InStageDistribution, FixedHorizonPOMDPWrapper

export
    FiVISolver,
    AlphaVec,
    solve

include("solver.jl")

end


text:
using ConvexBodyProximityQueries
using Test
using LinearAlgebra
using StaticArrays

@testset "ConvexBodyProximityQueries" begin
    include("support.jl")
    include("simplex.jl")
    include("proximity.jl")
end


text:
@testset "Opening models from BIGG" begin

    # E. coli models - realistic size models
    iJO1366_xml = download_data_file(
        "http://bigg.ucsd.edu/static/models/iJO1366.xml",
        joinpath("data", "iJO1366.xml"),
        "d6d9ec61ef6f155db5bb2f49549119dc13b96f6098b403ef82ea4240b27232eb",
    )
    sbmlmodel = load_model(iJO1366_xml)
    @test sbmlmodel isa SBMLModel
    @test n_reactions(sbmlmodel) == 2583

    iJO1366_mat = download_data_file(
        "http://bigg.ucsd.edu/static/models/iJO1366.mat",
        joinpath("data", "iJO1366.mat"),
        "b5cfe21b6369a00e45d600b783f89521f5cc953e25ee52c5f1d0a3f83743be30",
    )
    matlabmodel = load_model(iJO1366_mat)
    @test matlabmodel isa MATModel
    @test n_reactions(matlabmodel) == 2583

    iJO1366_json = download_data_file(
        "http://bigg.ucsd.edu/static/models/iJO1366.json",
        joinpath("data", "iJO1366.json"),
        "9376a93f62ad430719f23e612154dd94c67e0d7c9545ed9d17a4d0c347672313",
    )
    jsonmodel = load_model(iJO1366_json)
    @test jsonmodel isa JSONModel
    @test n_reactions(jsonmodel) == 2583

    @test Set(lowercase.(reactions(sbmlmodel))) ==
          Set("r_" .* lowercase.(reactions(matlabmodel)))
    @test Set(lowercase.(reactions(sbmlmodel))) ==
          Set("r_" .* lowercase.(reactions(jsonmodel)))
end


text:
# Julia wrapper for header: /Users/bieler/Downloads/SDL2_ttf-2.0.14/libSDL2_ttf.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0
include("SDL_ttf_h.jl")

function TTF_Linked_Version()
    ccall((:TTF_Linked_Version, libsdl2_ttf), Ptr{SDL_version}, ())
end

function TTF_ByteSwappedUNICODE(swapped::Cint)
    ccall((:TTF_ByteSwappedUNICODE, libsdl2_ttf), Cvoid, (Cint,), swapped)
end

function TTF_Init()
    ccall((:TTF_Init, libsdl2_ttf), Cint, ())
end

function TTF_OpenFont(file, ptsize)
    ccall((:TTF_OpenFont, libsdl2_ttf), Ptr{TTF_Font}, (Cstring, Cint), file, ptsize)
end

function TTF_OpenFontIndex(file, ptsize::Cint, index::Clong)
    ccall((:TTF_OpenFontIndex, libsdl2_ttf), Ptr{TTF_Font}, (Cstring, Cint, Clong), file, ptsize, index)
end

function TTF_OpenFontRW(src, freesrc::Cint, ptsize::Cint)
    ccall((:TTF_OpenFontRW, libsdl2_ttf), Ptr{TTF_Font}, (Ptr{RWops}, Cint, Cint), src, freesrc, ptsize)
end

function TTF_OpenFontIndexRW(src, freesrc::Cint, ptsize::Cint, index::Clong)
    ccall((:TTF_OpenFontIndexRW, libsdl2_ttf), Ptr{TTF_Font}, (Ptr{RWops}, Cint, Cint, Clong), src, freesrc, ptsize, index)
end

function TTF_GetFontStyle(font)
    ccall((:TTF_GetFontStyle, libsdl2_ttf), Cint, (Ptr{TTF_Font},), font)
end

function TTF_SetFontStyle(font, style::Cint)
    ccall((:TTF_SetFontStyle, libsdl2_ttf), Cvoid, (Ptr{TTF_Font}, Cint), font, style)
end

function TTF_GetFontOutline(font)
    ccall((:TTF_GetFontOutline, libsdl2_ttf), Cint, (Ptr{TTF_Font},), font)
end

function TTF_SetFontOutline(font, outline::Cint)
    ccall((:TTF_SetFontOutline, libsdl2_ttf), Cvoid, (Ptr{TTF_Font}, Cint), font, outline)
end

function TTF_GetFontHinting(font)
    ccall((:TTF_GetFontHinting, libsdl2_ttf), Cint, (Ptr{TTF_Font},), font)
end

function TTF_SetFontHinting(font, hinting::Cint)
    ccall((:TTF_SetFontHinting, libsdl2_ttf), Cvoid, (Ptr{TTF_Font}, Cint), font, hinting)
end

function TTF_FontHeight(font)
    ccall((:TTF_FontHeight, libsdl2_ttf), Cint, (Ptr{TTF_Font},), font)
end

function TTF_FontAscent(font)
    ccall((:TTF_FontAscent, libsdl2_ttf), Cint, (Ptr{TTF_Font},), font)
end

function TTF_FontDescent(font)
    ccall((:TTF_FontDescent, libsdl2_ttf), Cint, (Ptr{TTF_Font},), font)
end

function TTF_FontLineSkip(font)
    ccall((:TTF_FontLineSkip, libsdl2_ttf), Cint, (Ptr{TTF_Font},), font)
end

function TTF_GetFontKerning(font)
    ccall((:TTF_GetFontKerning, libsdl2_ttf), Cint, (Ptr{TTF_Font},), font)
end

function TTF_SetFontKerning(font, allowed::Cint)
    ccall((:TTF_SetFontKerning, libsdl2_ttf), Cvoid, (Ptr{TTF_Font}, Cint), font, allowed)
end

function TTF_FontFaces(font)
    ccall((:TTF_FontFaces, libsdl2_ttf), Clong, (Ptr{TTF_Font},), font)
end

function TTF_FontFaceIsFixedWidth(font)
    ccall((:TTF_FontFaceIsFixedWidth, libsdl2_ttf), Cint, (Ptr{TTF_Font},), font)
end

function TTF_FontFaceFamilyName(font)
    ccall((:TTF_FontFaceFamilyName, libsdl2_ttf), Cstring, (Ptr{TTF_Font},), font)
end

function TTF_FontFaceStyleName(font)
    ccall((:TTF_FontFaceStyleName, libsdl2_ttf), Cstring, (Ptr{TTF_Font},), font)
end

function TTF_GlyphIsProvided(font, ch::Uint16)
    ccall((:TTF_GlyphIsProvided, libsdl2_ttf), Cint, (Ptr{TTF_Font}, Uint16), font, ch)
end

function TTF_GlyphMetrics(font, ch::Uint16, minx, maxx, miny, maxy, advance)
    ccall((:TTF_GlyphMetrics, libsdl2_ttf), Cint, (Ptr{TTF_Font}, Uint16, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), font, ch, minx, maxx, miny, maxy, advance)
end

function TTF_SizeText(font, text, w, h)
    ccall((:TTF_SizeText, libsdl2_ttf), Cint, (Ptr{TTF_Font}, Cstring, Ptr{Cint}, Ptr{Cint}), font, text, w, h)
end

function TTF_SizeUTF8(font, text, w, h)
    ccall((:TTF_SizeUTF8, libsdl2_ttf), Cint, (Ptr{TTF_Font}, Cstring, Ptr{Cint}, Ptr{Cint}), font, text, w, h)
end

function TTF_SizeUNICODE(font, text, w, h)
    ccall((:TTF_SizeUNICODE, libsdl2_ttf), Cint, (Ptr{TTF_Font}, Ptr{Uint16}, Ptr{Cint}, Ptr{Cint}), font, text, w, h)
end

function TTF_RenderText_Solid(font, text, fg::Color)
    ccall((:TTF_RenderText_Solid, libsdl2_ttf), Ptr{Surface}, (Ptr{TTF_Font}, Cstring, Color), font, text, fg)
end

function TTF_RenderUTF8_Solid(font, text, fg::Color)
    ccall((:TTF_RenderUTF8_Solid, libsdl2_ttf), Ptr{Surface}, (Ptr{TTF_Font}, Cstring, Color), font, text, fg)
end

function TTF_RenderUNICODE_Solid(font, text, fg::Color)
    ccall((:TTF_RenderUNICODE_Solid, libsdl2_ttf), Ptr{Surface}, (Ptr{TTF_Font}, Ptr{Uint16}, Color), font, text, fg)
end

function TTF_RenderGlyph_Solid(font, ch::Uint16, fg::Color)
    ccall((:TTF_RenderGlyph_Solid, libsdl2_ttf), Ptr{Surface}, (Ptr{TTF_Font}, Uint16, Color), font, ch, fg)
end

function TTF_RenderText_Shaded(font, text, fg::Color, bg::Color)
    ccall((:TTF_RenderText_Shaded, libsdl2_ttf), Ptr{Surface}, (Ptr{TTF_Font}, Cstring, Color, Color), font, text, fg, bg)
end

function TTF_RenderUTF8_Shaded(font, text, fg::Color, bg::Color)
    ccall((:TTF_RenderUTF8_Shaded, libsdl2_ttf), Ptr{Surface}, (Ptr{TTF_Font}, Cstring, Color, Color), font, text, fg, bg)
end

function TTF_RenderUNICODE_Shaded(font, text, fg::Color, bg::Color)
    ccall((:TTF_RenderUNICODE_Shaded, libsdl2_ttf), Ptr{Surface}, (Ptr{TTF_Font}, Ptr{Uint16}, Color, Color), font, text, fg, bg)
end

function TTF_RenderGlyph_Shaded(font, ch::Uint16, fg::Color, bg::Color)
    ccall((:TTF_RenderGlyph_Shaded, libsdl2_ttf), Ptr{Surface}, (Ptr{TTF_Font}, Uint16, Color, Color), font, ch, fg, bg)
end

function TTF_RenderText_Blended(font, text, fg::Color)
    ccall((:TTF_RenderText_Blended, libsdl2_ttf), Ptr{Surface}, (Ptr{TTF_Font}, Cstring, Color), font, text, fg)
end

function TTF_RenderUTF8_Blended(font, text, fg::Color)
    ccall((:TTF_RenderUTF8_Blended, libsdl2_ttf), Ptr{Surface}, (Ptr{TTF_Font}, Cstring, Color), font, text, fg)
end

function TTF_RenderUNICODE_Blended(font, text, fg::Color)
    ccall((:TTF_RenderUNICODE_Blended, libsdl2_ttf), Ptr{Surface}, (Ptr{TTF_Font}, Ptr{Uint16}, Color), font, text, fg)
end

function TTF_RenderText_Blended_Wrapped(font, text, fg::Color, wrapLength::Uint32)
    ccall((:TTF_RenderText_Blended_Wrapped, libsdl2_ttf), Ptr{Surface}, (Ptr{TTF_Font}, Cstring, Color, Uint32), font, text, fg, wrapLength)
end

function TTF_RenderUTF8_Blended_Wrapped(font, text, fg::Color, wrapLength::Uint32)
    ccall((:TTF_RenderUTF8_Blended_Wrapped, libsdl2_ttf), Ptr{Surface}, (Ptr{TTF_Font}, Cstring, Color, Uint32), font, text, fg, wrapLength)
end

function TTF_RenderUNICODE_Blended_Wrapped(font, text, fg::Color, wrapLength::Uint32)
    ccall((:TTF_RenderUNICODE_Blended_Wrapped, libsdl2_ttf), Ptr{Surface}, (Ptr{TTF_Font}, Ptr{Uint16}, Color, Uint32), font, text, fg, wrapLength)
end

function TTF_RenderGlyph_Blended(font, ch::Uint16, fg::Color)
    ccall((:TTF_RenderGlyph_Blended, libsdl2_ttf), Ptr{Surface}, (Ptr{TTF_Font}, Uint16, Color), font, ch, fg)
end

function TTF_CloseFont(font)
    ccall((:TTF_CloseFont, libsdl2_ttf), Cvoid, (Ptr{TTF_Font},), font)
end

function TTF_Quit()
    ccall((:TTF_Quit, libsdl2_ttf), Cvoid, ())
end

function TTF_WasInit()
    ccall((:TTF_WasInit, libsdl2_ttf), Cint, ())
end

function TTF_GetFontKerningSize(font, prev_index::Cint, index::Cint)
    ccall((:TTF_GetFontKerningSize, libsdl2_ttf), Cint, (Ptr{TTF_Font}, Cint, Cint), font, prev_index, index)
end

function TTF_GetFontKerningSizeGlyphs(font, previous_ch::Uint16, ch::Uint16)
    ccall((:TTF_GetFontKerningSizeGlyphs, libsdl2_ttf), Cint, (Ptr{TTF_Font}, Uint16, Uint16), font, previous_ch, ch)
end

max_stars_repo_name: Kyjor/julia-game

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: Kyjor/julia-game

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: Kyjor/julia-game

max_forks_repo_licenses: MIT Unlicense


text:
## ----------------------------------------------------------------------------
# Functional interface
function _open_src_libs(cmd = "code")
    for path in bibtex_paths()
        !isfile(path) && continue
        run(`$(cmd) $(path)`; wait = false)
    end
end

function _open_book_root(cmd = "code")
    broot = bookdir()
    if isdir(broot)
        run(`$(cmd) $(broot)`; wait = false)
    end
end


text:
@enum CONTEXT_TYPE CPU=1 GPU=2 CPU_PINNED=3

"""
    Context(dev_type, dev_id)

A context describes the device type and id on which computation should be carried on.
"""
immutable Context
  device_type :: CONTEXT_TYPE
  device_id   :: Int
end
Context(dev_type :: Union{CONTEXT_TYPE, Int}, dev_id :: Int = 0) =
    Context(convert(CONTEXT_TYPE, dev_type), dev_id)

function Base.show(io :: IO, ctx :: Context)
  print(io, "$(ctx.device_type)$(ctx.device_id)")
end

"""
    cpu(dev_id)

Get a CPU context with a specific id. `cpu()` is usually the default context for many
operations when no context is specified.

# Arguments
* `dev_id::Int = 0`: the CPU id.
"""
function cpu(dev_id::Int=0)
  return Context(CPU, dev_id)
end

"""
    gpu(dev_id)

Get a GPU context with a specific id. The K GPUs on a node is typically numbered as 0,...,K-1.

# Arguments
* `dev_id :: Int = 0` the GPU device id.
"""
function gpu(dev_id::Int=0)
  return Context(GPU, dev_id)
end


text:
"""
    plot(path::RegularizationPath, args...; <keyword arguments>)

Plots a `RegularizationPath` fitted with the `Lasso` package.

The minimum AICc segment is represented by a solid vertical line and the CVmin
and CV1se cross-validation selected segments by dashed vertical lines.

By default it shows nonzero coefficients at the AICc in color and the rest grayed out.

LassoPlot uses Plots.jl, so you can set any supported backend before plotting,
and add any features of the plot afterwards.

# Example:
```julia
    using Lasso, LassoPath, Plots
    path = fit(LassoPath, X, y, dist, link)
    plot(path)
```
# Arguments
- `args...` additional arguments passed to Plots.plot()

# Keywords
- `x=:segment` one of (:segment, :λ, :logλ)
- `varnames=nothing` specify variable names
- `select=MinAICc()` Path segment selector (e.g. MinBIC(), MinCVmse(path, 5), MinCV1se(path, 5)) for which coefficients
    will be shown in color. The rest are grayed out.
- `showselectors=[MinAICc(), MinCVmse(path, 5), MinCV1se(path, 5)]` shown vertical lines
- `selectedvars=[]` Subset of the variables to present, or empty vector for all
- `kwargs...` additional keyword arguments passed along to fit(GammaLassoPath,...)
"""
function Plots.plot(path::RegularizationPath, args...;
    x=:segment, varnames=nothing, selectedvars=[], select=MinAICc(), showselectors=[MinAICc(), MinCVmse(path, 10), MinCV1se(path, 10)])
    β=coef(path)
    if hasintercept(path)
        β = β[2:end,:]
    end

    (p,nλ)=size(β)

    if varnames==nothing
        varnames=[Symbol("x$i") for i=1:p]
    end

    indata=DataFrame()
    if x==:λ
        indata[x]=path.λ
    elseif x==:logλ
        indata[x]=log.(path.λ)
    else
        x=:segment
        indata[x]=1:nλ
    end
    outdata = deepcopy(indata)

    # automatic selectors
    # xintercept = Float64[]
    dashed_vlines=Float64[]
    solid_vlines=Float64[]

    selectors = union([select], showselectors)
    ixselect = 0

    for s in selectors
        ixshown = segselect(path, s)
        if select == s
            push!(solid_vlines,indata[ixshown,x])
            ixselect = ixshown
        else
            push!(dashed_vlines,indata[ixshown,x])
        end
    end

    if length(selectedvars) == 0
        if isa(select, AllSeg)
            selectedvars = 1:p
        else
            selectedvars = findall(!iszero, β[:,ixselect])
        end
    end

    # colored paths
    for j in selectedvars
        indata[varnames[j]]=Vector(β[j,:])
    end

    # grayed out paths
    for j in setdiff(1:p,selectedvars)
        outdata[varnames[j]]=Vector(β[j,:])
    end

    inmdframe=melt(indata,x)
    outmdframe=melt(outdata,x)
    rename!(inmdframe,:value=>:coefficients)
    rename!(outmdframe,:value=>:coefficients)
    inmdframe = inmdframe[convert(BitArray,map(b->!isnan(b),inmdframe[:coefficients])),:]
    outmdframe = outmdframe[convert(BitArray,map(b->!isnan(b),outmdframe[:coefficients])),:]

    p = plot(xlabel=string(x), ylabel="Coefficient", args...)
    if size(inmdframe,1) > 0
      @df inmdframe plot!(cols(x), :coefficients, group=:variable)
    end
    if size(outmdframe,1) > 0
      @df outmdframe plot!(cols(x), :coefficients, group=:variable, palette=:grays)
    end
    if length(dashed_vlines) > 0
        vline!(dashed_vlines, line = (:dash, 0.5, 2, :black), label="")
    end
    if length(solid_vlines) > 0
        vline!(solid_vlines, line = (:solid, 0.5, 2, :black), label="")
    end

    p
end


text:
############
# dspev!
############

using LinearAlgebra

if isdefined(ProximalOperators, :dspev!) && isdefined(ProximalOperators, :dspevV!)

a = [1.0,2.0,3.0,5.0,6.0,9.0]
W_ref = [0.0,0.6992647456322766,14.300735254367698]
Z_ref = [0.9486832980505137	0.17781910596911388	-0.26149639682478454;
  0.0	-0.8269242138935418	-0.5623133863572413;
  -0.3162277660168381	0.5334573179073402	-0.7844891904743537]
A_ref = [1.0 2.0 3.0; 2.0 5.0 6.0; 3.0 6.0 9.0]

a_copy = copy(a); W1, Z1 = ProximalOperators.dspev!(:V,:L,a_copy)
a_copy = copy(a); W2, Z2 = ProximalOperators.dspevV!(:L,a_copy)

A1 = Z1*Diagonal(W1)*Z1'
A2 = Z2*Diagonal(W2)*Z2'

@test all((W1-W_ref)./(1 .+ abs.(W_ref)) .<= 1e-8)
@test all((A1-A_ref)./(1 .+ abs.(A_ref)) .<= 1e-8)
@test all((W2-W_ref)./(1 .+ abs.(W_ref)) .<= 1e-8)
@test all((A2-A_ref)./(1 .+ abs.(A_ref)) .<= 1e-8)

end


text:
module B

using A

export BType

struct BType
	x::AType
end

struct SameNameType
    y::Int
end

end


text:
"""
Represents any kind of wrapper structure that was generated from a Vulkan structure.
`D` is a `Bool` parameter indicating whether the structure has specific dependencies or not.
"""
abstract type VulkanStruct{D} end

Base.broadcastable(x::VulkanStruct) = Ref(x)

"""
Represents a structure that will never be requested by API functions.
"""
abstract type ReturnedOnly <: VulkanStruct{false} end

const FunctionPtr = Union{Ptr{Cvoid}, Base.CFunction}
const OptionalPtr{T} = Union{T, Ptr{Cvoid}}


text:
using LinearAlgebra
using CompScienceMeshes
using SauterSchwab3D

const pI = point(1,0,0)
const pII = point(0,0,0)
const pIII = point(0.5,0.7,0)
const pIV = point(0.8,-0.3,0)

const P = simplex(pI,pIII,pII)
const Q = simplex(pI,pIV,pII)

const sing = SauterSchwab3D.singularity_detection(P,Q)

Accuracy2 = 30
ce_ref = CommonEdge4D(sing,SauterSchwab3D._legendre(Accuracy2,0.0,1.0))

function integrand(x,y)
      return ((x-pI)'*(y-pII))*exp(-im*1*norm(x-y))/(4pi*norm(x-y))
end

function INTEGRAND(u,v)
   if u[2]< 0.0 || u[2]>1-u[1] || u[1]<0.0 || u[1]>1.0
      println(u)
   end
   if v[2]< 0.0 || v[2]>1-v[1] || v[1]<0.0 || v[1]>1.0
      println(v)
   end
   n1 = neighborhood(P,u)
   n2 = neighborhood(Q,v)
   x = cartesian(n1)
   y = cartesian(n2)
   output = integrand(x,y)*jacobian(n1)*jacobian(n2)
   return(output)
end

print("Ref: ")
ref = sauterschwab_parameterized(INTEGRAND, ce_ref)
println(ref)
println()

res_tp =[]
res_sp =[]
res_gm = []
n1 = []
n2 = []
n3 = []
for i in 2:1:20
   Accuracy = i
   ce = CommonEdge4D(sing,SauterSchwab3D._legendre(Accuracy,0.0,1.0))

   int_tp = sauterschwab_parameterized(INTEGRAND, ce)

   num_pts = 5*length(ce.qps)^4
   push!(n1,num_pts)

   push!(res_tp,int_tp)
end

for i in 2:1:7
   Accuracy = i
   ce_s = CommonEdge4D_S(sing,(SauterSchwab3D._legendre(Accuracy,0.0,1.0),
                          SauterSchwab3D._shunnham2D(Accuracy)))

   int_sp = sauterschwab_parameterized(INTEGRAND, ce_s)

   num_pts = length(ce_s.qps[1])^2*length(ce_s.qps[2]) + 4*length(ce_s.qps[2])^2
   push!(n2,num_pts)
   push!(res_sp,int_sp)
end

for i in 2:1:12
   Accuracy = i
   ce_gm = CommonEdge4D_S(sing,(SauterSchwab3D._legendre(Accuracy,0.0,1.0),
                          SauterSchwab3D._grundmannMoeller2D(2*Accuracy-1)))

   int_gm = sauterschwab_parameterized(INTEGRAND, ce_gm)

   num_pts = length(ce_gm.qps[1])^2*length(ce_gm.qps[2]) + 4*length(ce_gm.qps[2])^2
   push!(n3,num_pts)
   push!(res_gm,int_gm)
end

err_tp = norm.(res_tp.-ref)/norm(ref)
err_sp = norm.(res_sp.-ref)/norm(ref)
err_gm = norm.(res_gm.-ref)/norm(ref)

println(err_tp)
println(err_sp)
println(err_gm)

using Plots
plot( yaxis=:log,xaxis=:log, fontfamily="Times")
plot!(n1,err_tp, label="Gauss Tensor-Product",markershape=:circle)
plot!(n2,err_sp, label="Simplex Tensor-Product",markershape=:rect)
#plot!(n3,err_gm, label="Simplex-Product GM",markershape=:x)
plot!(xlims=(3e1,1e5),ylims=(1e-7,1e-1))
plot!(xlabel="Quad. points/Func. evals.", ylabel="Rel. Error.", title="Common Edge 4D",legend=:bottomleft)

using BenchmarkTools

ref = -0.0028991793711368725 + 0.0012868587156426375im

ce = CommonEdge4D(sing,SauterSchwab3D._legendre(5,0.0,1.0))

int_tp = sauterschwab_parameterized(INTEGRAND, ce)

num_pts = 5*length(ce.qps)^4
println("#Pts: ",num_pts)
err_tp = norm.(int_tp.-ref)/norm(ref)
println("Rel Err: ",err_tp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, ce) end

ce_s = CommonEdge4D_S(sing,(SauterSchwab3D._legendre(5,0.0,1.0),
SauterSchwab3D._shunnham2D(5)))

int_sp = sauterschwab_parameterized(INTEGRAND, ce_s)

num_pts = length(ce_s.qps[1])^2*length(ce_s.qps[2]) + 4*length(ce_s.qps[2])^2

println("#Pts: ",num_pts)
err_sp = norm.(int_sp.-ref)/norm(ref)
println("Rel Err: ",err_sp)
@time   for i in 1:100  sauterschwab_parameterized(INTEGRAND, ce_s) end

#=

ce_gm = CommonEdge4D_S((SauterSchwab3D._legendre(5,0.0,1.0),
SauterSchwab3D._grundmannMoeller2D(2*5-1)))

int_gm = sauterschwab_parameterized(INTEGRAND, ce_gm)

num_pts = length(ce_gm.qps[1])^2*length(ce_gm.qps[2]) + 4*length(ce_gm.qps[2])^2

println("#Pts: ",num_pts)
err_gm = norm.(int_gm.-ref)/norm(ref)
println("Rel Err: ",err_gm)
=#