text:
"""
    M = trmatop(n, m)

Define the transposition operator `M: X -> X'` for all `n x m` matrices.
"""
function trmatop(n::Int,m::Int)
  function prod(x)
    X = reshape(x, n, m)
    return adjoint(X)[:]
  end
  function tprod(x)
    X = reshape(x, m, n)
    return adjoint(X)[:]
  end
  function ctprod(x)
    X = reshape(x, m, n)
    return adjoint(X)[:]
  end
  F1 = typeof(prod)
  F2 = typeof(tprod)
  F3 = typeof(ctprod)
  m == n ? sym = true : sym = false
  return LinearOperator{Int}(n * m, n * m, sym, sym, prod, tprod, ctprod)
end
trmatop(n::Int) = trmatop(n,n)
trmatop(dims::Tuple{Int,Int}) = trmatop(dims[1],dims[2])
"""
    M = trmatop(A)

Define the transposition operator `M: X -> X'` of all matrices of the size of `A`.
"""
trmatop(A) = trmatop(size(A))
"""
    L = lyapop(A; disc = false, her = false)

Define, for an `n x n` matrix `A`, the continuous Lyapunov operator `L:X -> AX+XA'`
if `disc = false` or the discrete Lyapunov operator `L:X -> AXA'-X` if `disc = true`.
If `her = false` the Lyapunov operator `L:X -> Y` maps general square matrices `X`
into general square matrices `Y`, and the associated matrix `M = Matrix(L)` is
``n^2 \\times n^2``.
If `her = true` the Lyapunov operator `L:X -> Y` maps symmetric/Hermitian matrices `X`
into symmetric/Hermitian matrices `Y`, and the associated matrix `M = Matrix(L)` is
``n(n+1)/2 \\times n(n+1)/2``.
For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov
operators. Linear Algebra and its Applications 312:35–71, 2000.
"""
function lyapop(A; disc = false, her = false)
  n = LinearAlgebra.checksquare(A)
  T = eltype(A)
  function prod(x)
    T1 = promote_type(T, eltype(x))
    if her
      X = vec2triu(convert(Vector{T1}, x),her = true)
      if disc
        return triu2vec(utqu(X,A') - X)
      else
        Y = A * X
        return triu2vec(Y + Y')
      end
    else
      X = reshape(convert(Vector{T1}, x), n, n)
      if disc
        Y = A*X*A' - X
      else
        Y = A*X + X*A'
      end
      return Y[:]
    end
  end
  function tprod(x)
    T1 = promote_type(T, eltype(x))
    if her
      X = vec2triu(convert(Vector{T1}, x),her = true)
      if disc
        return triu2vec(utqu(X,A) - X)
      else
        Y = X * A
        return triu2vec(Y + adjoint(Y))
      end
    else
      X = reshape(convert(Vector{T1}, x), n, n)
      if disc
         Y = adjoint(A)*X*A - X
       else
         Y = adjoint(A)*X + X*A
       end
       return Y[:]
    end
  end
  function ctprod(x)
    T1 = promote_type(T, eltype(x))
    if her
      X = vec2triu(convert(Vector{T1}, x),her = true)
      if disc
        return triu2vec(utqu(X,A) - X)
      else
        Y = X * A
        return triu2vec(Y + Y')
      end
    else
      X = reshape(convert(Vector{T1}, x), n, n)
      if disc
        return (A'*X*A - X )[:]
      else
        return (A'*X + X*A)[:]
      end
    end
  end
  F1 = typeof(prod)
  F2 = typeof(tprod)
  F3 = typeof(ctprod)
  her ? N = Int(n*(n+1)/2) : N = n*n
  return LinearOperator{T}(N, N, false, false, prod, tprod, ctprod)
end
"""
    L = lyapop(A, E; disc = false, her = false)

Define, for a pair `(A,E)` of `n x n` matrices, the continuous Lyapunov operator `L:X -> AXE'+EXA'`
if `disc = false` or the discrete Lyapunov operator `L:X -> AXA'-EXE'` if `disc = true`.
If `her = false` the Lyapunov operator `L:X -> Y` maps general square matrices `X`
into general square matrices `Y`, and the associated matrix `M = Matrix(L)` is
``n^2 \\times n^2``.
If `her = true` the Lyapunov operator `L:X -> Y` maps symmetric/Hermitian matrices `X`
into symmetric/Hermitian matrices `Y`, and the associated `M = Matrix(L)` is a
``n(n+1)/2 \\times n(n+1)/2``.
For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov
operators. Linear Algebra and its Applications 312:35–71, 2000.
"""
function lyapop(A, E; disc = false, her = false)
  n = LinearAlgebra.checksquare(A)
  if n != LinearAlgebra.checksquare(E)
    throw(DimensionMismatch("E must be a square matrix of dimension $n"))
  end
  T = promote_type(eltype(A), eltype(E))
  function prod(x)
    T1 = promote_type(T, eltype(x))
    if her
      X = vec2triu(convert(Vector{T1}, x),her = true)
      if disc
        return triu2vec(utqu(X,A') - utqu(X,E'))
      else
        Y = A * X * E'
        return triu2vec(Y + Y')
      end
    else
      X = reshape(convert(Vector{T1}, x), n, n)
      if disc
        Y = A*X*A' - E*X*E'
      else
        Y = A*X*E' + E*X*A'
      end
      return Y[:]
    end
  end
  function tprod(x)
    T1 = promote_type(T, eltype(x))
    if her
      X = vec2triu(convert(Vector{T1}, x),her = true)
      if disc
        return triu2vec(utqu(X,A) - utqu(X,E))
      else
        Y = E' * X * A
        return triu2vec(Y + Y')
      end
    else
      X = reshape(convert(Vector{T1}, x), n, n)
      if disc
        return (A'*X*A - E'*X*E )[:]
      else
        return (A'*X*E + E'*X*A)[:]
      end
    end
  end
  function ctprod(x)
    T1 = promote_type(T, eltype(x))
    if her
      X = vec2triu(convert(Vector{T1}, x),her = true)
      if disc
        return triu2vec(utqu(X,A) - utqu(X,E))
      else
        Y = E' * X * A
        return triu2vec(Y + Y')
      end
    else
      X = reshape(convert(Vector{T1}, x), n, n)
      if disc
        return (A'*X*A - E'*X*E )[:]
      else
        return (A'*X*E + E'*X*A)[:]
      end
    end
  end
  F1 = typeof(prod)
  F2 = typeof(tprod)
  F3 = typeof(ctprod)
  her ? N = Int(n*(n+1)/2) : N = n*n
  return LinearOperator{T}(N, N, false, false, prod, tprod, ctprod)
end
"""
    LINV = invlyapop(A; disc = false, her = false)

Define `LINV`, the inverse of the continuous Lyapunov operator `L:X -> AX+XA'` for `disc = false`
or the inverse of the discrete Lyapunov operator `L:X -> AXA'-X` for `disc = true`, where
`A` is an `n x n` matrix.
If `her = false` the inverse Lyapunov operator `LINV:Y -> X` maps general square matrices `Y`
into general square matrices `X`, and the associated matrix `M = Matrix(LINV)` is
``n^2 \\times n^2``.
If `her = true` the inverse Lyapunov operator `LINV:Y -> X` maps symmetric/Hermitian matrices `Y`
into symmetric/Hermitian matrices `X`, and the associated matrix `M = Matrix(LINV)` is
``n(n+1)/2 \\times n(n+1)/2``.
For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov
operators. Linear Algebra and its Applications 312:35–71, 2000.
"""
function invlyapop(A; disc = false, her = false)
   n = LinearAlgebra.checksquare(A)
   T = eltype(A)
   function prod(x)
    T1 = promote_type(T, eltype(x))
     try
       if her
         Y = vec2triu(convert(Vector{T1}, x),her = true)
         if disc
            return triu2vec(lyapd(A,-Y))
         else
             return triu2vec(lyapc(A,-Y))
         end
       else
         Y = reshape(convert(Vector{T1}, x), n, n)
         if disc
           return sylvd(-A,A',-Y)[:]
         else
           return sylvc(A,A',Y)[:]
         end
       end
     catch err
      #  if isnothing(findfirst("LAPACKException",string(err))) &&
      #     isnothing(findfirst("SingularException",string(err)))
      if findfirst("LAPACKException",string(err)) === nothing &&
          findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
     end
   end
   function tprod(x)
    T1 = promote_type(T, eltype(x))
     try
       if her
         Y = vec2triu(convert(Vector{T1}, x),her = true)
         if disc
           return triu2vec(lyapd(A',-Y))
         else
           return triu2vec(lyapc(A',-Y))
         end
       else
         Y = reshape(convert(Vector{T1}, x), n, n)
         if disc
           return sylvd(-A',A,-Y)[:]
         else
            return sylvc(A',A,Y)[:]
         end
       end
     catch err
      #  if isnothing(findfirst("LAPACKException",string(err))) &&
      #     isnothing(findfirst("SingularException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing &&
          findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
     end
   end
   function ctprod(x)
    T1 = promote_type(T, eltype(x))
     try
       if her
         Y = vec2triu(convert(Vector{T1}, x),her = true)
         if disc
           return triu2vec(lyapd(A',-Y))
         else
           return triu2vec(lyapc(A',-Y))
         end
       else
         Y = reshape(convert(Vector{T1}, x), n, n)
         if disc
           return sylvd(-A',A,-Y)[:]
         else
           return sylvc(A',A,Y)[:]
         end
       end
     catch err
      #  if isnothing(findfirst("LAPACKException",string(err))) &&
      #     isnothing(findfirst("SingularException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing &&
          findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
      end
     end
   end
   F1 = typeof(prod)
   F2 = typeof(tprod)
   F3 = typeof(ctprod)
   her ? N = Int(n*(n+1)/2) : N = n*n
   return LinearOperator{T}(N, N, false, false, prod, tprod, ctprod)
end
"""
    LINV = invlyapop(A, E; disc = false, her = false)

Define `LINV`, the inverse of the continuous Lyapunov operator `L:X -> AXE'+EXA'` for `disc = false`
or the inverse of the discrete Lyapunov operator `L:X -> AXA'-EXE'` for `disc = true`, where
`(A,E)` is a pair of `n x n` matrices.
If `her = false` the inverse Lyapunov operator `LINV:Y -> X` maps general square matrices `Y`
into general square matrices `X`, and the associated matrix `M = Matrix(LINV)` is
``n^2 \\times n^2``.
If `her = true` the inverse Lyapunov operator `LINV:Y -> X` maps symmetric/Hermitian matrices `Y`
into symmetric/Hermitian matrices `X`, and the associated matrix `M = Matrix(LINV)` is
``n(n+1)/2 \\times n(n+1)/2``.
For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov
operators. Linear Algebra and its Applications 312:35–71, 2000.
"""
function invlyapop(A, E; disc = false, her = false)
   n = LinearAlgebra.checksquare(A)
   if n != LinearAlgebra.checksquare(E)
     throw(DimensionMismatch("E must be a square matrix of dimension $n"))
   end
   T = promote_type(eltype(A), eltype(E))
   function prod(x)
    T1 = promote_type(T, eltype(x))
     try
       if her
         Y = vec2triu(convert(Vector{T1}, x),her = true)
         if disc
            return triu2vec(lyapd(A,E,-Y))
         else
            return triu2vec(lyapc(A,E,-Y))
         end
       else
         Y = reshape(convert(Vector{T1}, x), n, n)
         if disc
           return gsylv(-A,A',E,E',-Y)[:]
         else
           return gsylv(A,E',E,A',Y)[:]
         end
       end
     catch err
      #  if isnothing(findfirst("LAPACKException",string(err))) &&
      #     isnothing(findfirst("SingularException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing &&
          findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
     end
   end
   function tprod(x)
    T1 = promote_type(T, eltype(x))
     try
       if her
         Y = vec2triu(convert(Vector{T1}, x),her = true)
         if disc
            return triu2vec(lyapd(A',E',-Y))
         else
            return triu2vec(lyapc(A',E',-Y))
         end
       else
         Y = reshape(convert(Vector{T1}, x), n, n)
         if disc
            return gsylv(-A',A,E',E,-Y)[:]
         else
            return gsylv(A',E,E',A,Y)[:]
         end
       end
     catch err
      #  if isnothing(findfirst("LAPACKException",string(err))) &&
      #     isnothing(findfirst("SingularException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing &&
          findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
     end
   end
   function ctprod(x)
    T1 = promote_type(T, eltype(x))
     try
       if her
         Y = vec2triu(convert(Vector{T1}, x),her = true)
         if disc
           return triu2vec(lyapd(A',E',-Y))
         else
           return triu2vec(lyapc(A',E',-Y))
         end
       else
         Y = reshape(convert(Vector{T1}, x), n, n)
         if disc
           return gsylv(-A',A,E',E,-Y)[:]
         else
           return gsylv(A',E,E',A,Y)[:]
         end
       end
     catch err
      #  if isnothing(findfirst("LAPACKException",string(err))) &&
      #     isnothing(findfirst("SingularException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing &&
          findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
      end
     end
   end
   F1 = typeof(prod)
   F2 = typeof(tprod)
   F3 = typeof(ctprod)
   her ? N = Int(n*(n+1)/2) : N = n*n
   return LinearOperator{T}(N, N, false, false, prod, tprod, ctprod)
end
"""
    LINV = invlyapsop(A; disc = false, her = false)

Define `LINV`, the inverse of the continuous Lyapunov operator `L:X -> AX+XA'` for `disc = false`
or the inverse of the discrete Lyapunov operator `L:X -> AXA'-X` for `disc = true`, where
`A` is an `n x n` matrix in Schur form.
If `her = false` the inverse Lyapunov operator `LINV:Y -> X` maps general square matrices `Y`
into general square matrices `X`, and the associated matrix `M = Matrix(LINV)` is
``n^2 \\times n^2``.
If `her = true` the inverse Lyapunov operator `LINV:Y -> X` maps symmetric/Hermitian matrices `Y`
into symmetric/Hermitian matrices `X`, and the associated matrix `M = Matrix(LINV)` is
``n(n+1)/2 \\times n(n+1)/2``.
For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov
operators. Linear Algebra and its Applications 312:35–71, 2000.
"""
function invlyapsop(A; disc = false, her = false)
   n = LinearAlgebra.checksquare(A)
   T = eltype(A)
   if !(T <: BlasFloat)
      T = promote_type(Float64,T)
   end
   if eltype(A) !== T
      A = convert(Matrix{T},A)
   end

   # check A is in Schur form
   if !isschur(A)
       error("The matrix A must be in Schur form")
   end
   cmplx = T<:Complex
   function prod(x)
     T1 = promote_type(T, eltype(x))
     if T !== T1
        if cmplx
           A = convert(Matrix{T1},A)
        else
           T1r = real(T1)
           if T1r !== T
              A = convert(Matrix{T1r},A)
           end
        end
     end
     try
       if her
         Y = vec2triu(convert(Vector{T1}, -x),her = true)
         disc ? lyapds!(A,Y) : lyapcs!(A,Y)
         return triu2vec(Y)
       else
         Y = reshape(convert(Vector{T1}, -x), n, n)
         if disc
            sylvds!(-A,A,Y,adjB = true)
            return Y[:]
         else
            sylvcs!(A,A,Y,adjB = true)
            return -Y[:]
         end
       end
     catch err
      #  if isnothing(findfirst("LAPACKException",string(err))) &&
      #     isnothing(findfirst("SingularException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing &&
          findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
     end
   end
   function tprod(x)
     T1 = promote_type(T, eltype(x))
     if T !== T1
        if cmplx
           A = convert(Matrix{T1},A)
        else
           T1r = real(T1)
           if T1r !== T
              A = convert(Matrix{T1r},A)
           end
        end
     end
     try
      if her
        Y = vec2triu(convert(Vector{T1}, -x),her = true)
        disc ? lyapds!(A,Y,adj = true) : lyapcs!(A,Y,adj = true)
        return triu2vec(Y)
      else
        Y = reshape(convert(Vector{T1}, -x), n, n)
        if disc
           sylvds!(-A,A,Y,adjA = true)
           return Y[:]
        else
           sylvcs!(A,A,Y,adjA = true)
         #   realcase = eltype(A) <: AbstractFloat && eltype(Y) <: AbstractFloat
         #   realcase ? (TA,TB) = ('T','N') : (TA,TB) = ('C','N')
         #   Y, scale = LAPACK.trsyl!(TA, TB, A, A, Y)
         #   rmul!(Y, inv(-scale))
            return -Y[:]
        end
       end
     catch err
      #  if isnothing(findfirst("LAPACKException",string(err))) &&
      #     isnothing(findfirst("SingularException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing &&
          findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
     end
   end
   function ctprod(x)
     T1 = promote_type(T, eltype(x))
     if T !== T1
       if cmplx
          A = convert(Matrix{T1},A)
       else
          T1r = real(T1)
          if T1r !== T
             A = convert(Matrix{T1r},A)
          end
       end
     end
     try
      if her
        Y = vec2triu(convert(Vector{T1}, -x),her = true)
        disc ? lyapds!(A,Y,adj = true) : lyapcs!(A,Y,adj = true)
        return triu2vec(Y)
      else
        Y = reshape(convert(Vector{T1}, -x), n, n)
        if disc
           sylvds!(-A,A,Y,adjA = true)
           return Y[:]
        else
           sylvcs!(A,A,Y,adjA = true)
         #   realcase = eltype(A) <: AbstractFloat && eltype(Y) <: AbstractFloat
         #   realcase ? (TA,TB) = ('T','N') : (TA,TB) = ('C','N')
         #   Y, scale = LAPACK.trsyl!(TA, TB, A, A, Y)
         #   rmul!(Y, inv(-scale))
           return -Y[:]
        end
      end
    catch err
      #  if isnothing(findfirst("LAPACKException",string(err))) &&
      #     isnothing(findfirst("SingularException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing &&
          findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
      end
     end
   end
   F1 = typeof(prod)
   F2 = typeof(tprod)
   F3 = typeof(ctprod)
   her ? N = Int(n*(n+1)/2) : N = n*n
   return LinearOperator{T}(N, N, false, false, prod, tprod, ctprod)
end
invlyapsop(A :: Schur; disc = false, her = false) = invlyapsop(A.T,disc = disc,her = her)
invlyapsop(A :: Adjoint; disc = false, her = false) = invlyapsop(A.parent,disc = disc,her = her)'
"""
    LINV = invlyapsop(A, E; disc = false, her = false)

Define `LINV`, the inverse of the continuous Lyapunov operator `L:X -> AXE'+EXA'` for `disc = false`
or the inverse of the discrete Lyapunov operator `L:X -> AXA'-EXE'` for `disc = true`, where
`(A,E)` is a pair of `n x n` matrices in generalized Schur form.
If `her = false` the inverse Lyapunov operator `LINV:Y -> X` maps general square matrices `Y`
into general square matrices `X`, and the associated matrix `M = Matrix(LINV)` is
``n^2 \\times n^2``.
If `her = true` the inverse Lyapunov operator `LINV:Y -> X` maps symmetric/Hermitian matrices `Y`
into symmetric/Hermitian matrices `X`, and the associated matrix `M = Matrix(LINV)` is
``n(n+1)/2 \\times n(n+1)/2``.
For the definitions of the Lyapunov operators see:

M. Konstantinov, V. Mehrmann, P. Petkov. On properties of Sylvester and Lyapunov
operators. Linear Algebra and its Applications 312:35–71, 2000.
"""
function invlyapsop(A, E; disc = false, her = false)
   n = LinearAlgebra.checksquare(A)
   if n != LinearAlgebra.checksquare(E)
     throw(DimensionMismatch("E must be a square matrix of dimension $n"))
   end
   if isa(A,Adjoint) || isa(E,Adjoint)
     error("No calls with adjoint matrices are supported")
   end
   T = promote_type(eltype(A), eltype(E))
   if !(T <: BlasFloat)
      T = promote_type(Float64,T)
   end
   if eltype(A) !== T
      A = convert(Matrix{T},A)
   end
   if eltype(A) !== T
     A = convert(Matrix{T},A)
   end
   if eltype(E) !== T
     E = convert(Matrix{T},E)
   end
   cmplx = T<:Complex

   # check (A,E) is in generalized Schur form
   if !isschur(A,E)
       error("The matrix pair (A,E) must be in generalized Schur form")
   end
   function prod(x)
     T1 = promote_type(T, eltype(x))
     if T !== T1
       if cmplx
         A = convert(Matrix{T1},A)
         E = convert(Matrix{T1},E)
       else
         T1r = real(T1)
         if T1r !== T
            A = convert(Matrix{T1r},A)
            E = convert(Matrix{T1},E)
         end
       end
     end
     try
       if her
         Y = vec2triu(convert(Vector{T1}, -x),her = true)
         disc ? lyapds!(A,E,Y) : lyapcs!(A,E,Y)
         return triu2vec(Y)
       else
         Y = copy(reshape(convert(Vector{T1}, x), n, n))
         disc ? gsylvs!(A,A,-E,E,Y,adjBD = true) :
                gsylvs!(A,E,E,A,Y,adjBD = true,DBSchur = true)
         return Y[:]
       end
     catch err
      #  if isnothing(findfirst("LAPACKException",string(err))) &&
      #     isnothing(findfirst("SingularException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing &&
          findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
     end
   end
   function tprod(x)
    T1 = promote_type(T, eltype(x))
    if T !== T1
      if cmplx
         A = convert(Matrix{T1},A)
         E = convert(Matrix{T1},E)
      else
         T1r = real(T1)
         if T1r !== T
            A = convert(Matrix{T1r},A)
            E = convert(Matrix{T1},E)
         end
      end
    end
    try
       if her
         Y = vec2triu(convert(Vector{T1}, -x),her = true)
         disc ? lyapds!(A,E,Y,adj = true) : lyapcs!(A,E,Y,adj = true)
         return triu2vec(Y)
       else
         Y = copy(reshape(convert(Vector{T1}, x), n, n))
         disc ? gsylvs!(A,A,-E,E,Y,adjAC = true) :
                gsylvs!(A,E,E,A,Y,adjAC = true,DBSchur = true)
         return Y[:]
       end
     catch err
      #  if isnothing(findfirst("LAPACKException",string(err))) &&
      #     isnothing(findfirst("SingularException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing &&
          findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
     end
   end
   function ctprod(x)
    T1 = promote_type(T, eltype(x))
    if T !== T1
      if cmplx
         A = convert(Matrix{T1},A)
         E = convert(Matrix{T1},E)
      else
         T1r = real(T1)
         if T1r !== T
            A = convert(Matrix{T1r},A)
            E = convert(Matrix{T1},E)
         end
      end
    end
    try
       if her
         Y = vec2triu(convert(Vector{T1}, -x),her = true)
         disc ? lyapds!(A,E,Y,adj = true) : lyapcs!(A,E,Y,adj = true)
         return triu2vec(Y)
       else
         Y = copy(reshape(convert(Vector{T1}, x), n, n))
         disc ? gsylvs!(A,A,-E,E,Y,adjAC = true) :
                gsylvs!(A,E,E,A,Y,adjAC = true,DBSchur = true)
         return Y[:]
       end
     catch err
      #  if isnothing(findfirst("LAPACKException",string(err))) &&
      #     isnothing(findfirst("SingularException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing &&
          findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
      end
     end
   end
   F1 = typeof(prod)
   F2 = typeof(tprod)
   F3 = typeof(ctprod)
   her ? N = Int(n*(n+1)/2) : N = n*n
   return LinearOperator{T}(N, N, false, false, prod, tprod, ctprod)
end
invlyapsop(AE :: GeneralizedSchur; disc = false, her = false) = invlyapsop(AE.S, AE.T, disc = disc, her = her)
invlyapsop(A :: Adjoint, E :: Adjoint; disc = false, her = false) = invlyapsop(A.parent,E.parent,disc = disc,her = her)'
"""
    M = sylvop(A, B; disc = false)

Define the continuous Sylvester operator `M: X -> AX+XB` if `disc = false`
or the discrete Sylvester operator `M: X -> AXB+X` if `disc = true`, where `A` and `B` are square matrices.
"""
function sylvop(A, B; disc = false)
  m = LinearAlgebra.checksquare(A)
  n = LinearAlgebra.checksquare(B)
  T = promote_type(eltype(A), eltype(B))
  function prod(x)
    T1 = promote_type(T, eltype(x))
    X = reshape(convert(Vector{T1}, x), m, n)
    disc ? Y = A * X * B + X : Y = A * X + X * B
    return Y[:]
  end
  function tprod(x)
    T1 = promote_type(T, eltype(x))
    X = reshape(convert(Vector{T1}, x), m, n)
    disc ? Y = adjoint(A)*X*adjoint(B) + X : Y = adjoint(A)*X + X*adjoint(B)
    return Y[:]
  end
  function ctprod(x)
    T1 = promote_type(T, eltype(x))
    X = reshape(convert(Vector{T1}, x), m, n)
    disc ? Y = A'*X*B' + X : Y = A'*X + X*B'
    return Y[:]
  end
  F1 = typeof(prod)
  F2 = typeof(tprod)
  F3 = typeof(ctprod)
  return LinearOperator{T}(m * n, n * m, false, false, prod, tprod, ctprod)
end
"""
    MINV = invsylvop(A, B; disc = false)

Define MINV, the inverse of the continuous Sylvester operator  `M: X -> AX+XB` if `disc = false`
or of the discrete Sylvester operator `M: X -> AXB+X` if `disc = true`, where `A` and `B` are square matrices.
"""
function invsylvop(A, B; disc = false)
  m = LinearAlgebra.checksquare(A)
  n = LinearAlgebra.checksquare(B)
  T = promote_type(eltype(A), eltype(B))
  function prod(x)
    T1 = promote_type(T, eltype(x))
    C = reshape(convert(Vector{T1}, x), m, n)
    try
      if disc
        return sylvd(A,B,C)[:]
      else
        return sylvc(A,B,C)[:]
      end
    catch err
      #  if isnothing(findfirst("LAPACKException",string(err))) &&
      #     isnothing(findfirst("SingularException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing &&
          findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
    end
  end
  function tprod(x)
    T1 = promote_type(T, eltype(x))
    C = reshape(convert(Vector{T1}, x), m, n)
    try
      if disc
        return sylvd(A',B',C)[:]
      else
        return sylvc(A',B',C)[:]
     end
    catch err
      #  if isnothing(findfirst("LAPACKException",string(err))) &&
      #     isnothing(findfirst("SingularException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing &&
          findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
         throw("ME:SingularException: Singular operator")
       end
    end
  end
  function ctprod(x)
    T1 = promote_type(T, eltype(x))
    C = reshape(convert(Vector{T1}, x), m, n)
    try
      if disc
        return sylvd(A',B',C)[:]
      else
        return sylvc(A',B',C)[:]
     end
    catch err
      #  if isnothing(findfirst("LAPACKException",string(err))) &&
      #     isnothing(findfirst("SingularException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing &&
          findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
         throw("ME:SingularException: Singular operator")
       end
    end
  end
  F1 = typeof(prod)
  F2 = typeof(tprod)
  F3 = typeof(ctprod)
  return LinearOperator{T}(m * n, n * m, false, false, prod, tprod, ctprod)
end
"""
    MINV = invsylvsop(A, B; disc = false)

Define MINV, the inverse of the continuous Sylvester operator  `M: X -> AX+XB` if `disc = false`
or of the discrete Sylvester operator `M: X -> AXB+X` if `disc = true`, where `A` and `B` are square matrices in Schur forms.
"""
function invsylvsop(A, B; disc = false)
  m = LinearAlgebra.checksquare(A)
  n = LinearAlgebra.checksquare(B)
  T = promote_type(eltype(A), eltype(B))
  if !(T <: BlasFloat)
     T = promote_type(Float64,T)
  end
  adjA = isa(A,Adjoint)
  adjB = isa(B,Adjoint)
  if eltype(A) !== T
    adjA ? A = convert(Matrix{T},A.parent)'  : A = convert(Matrix{T},A)
  end
  if eltype(B) !== T
    adjB ? B = convert(Matrix{T},B.parent)' :  B = convert(Matrix{T},B)
  end
  cmplx = T<:Complex
  if adjA
     if !isschur(A.parent)
         error("A must be in Schur form")
     end
  else
     if !isschur(A)
        error("A must be in Schur form")
     end
  end
  if adjB
     if !isschur(B.parent)
         error("B must be in Schur form")
     end
  else
     if !isschur(B)
        error("B must be in Schur form")
     end
  end
  function prod(x)
    T1 = promote_type(T, eltype(x))
    C = copy(reshape(convert(Vector{T1}, x), m, n))
    if T !== T1
      if cmplx
         adjA ? A = convert(Matrix{T1},A.parent)' : A = convert(Matrix{T1},A)
         adjB ? B = convert(Matrix{T1},B.parent)' : B = convert(Matrix{T1},B)
      else
        T1r = real(T1)
        if T1r !== T
          adjA ? A = convert(Matrix{T1r},A.parent)' : A = convert(Matrix{T1r},A)
          adjB ? B = convert(Matrix{T1r},B.parent)' : B = convert(Matrix{T1r},B)
        end
      end
    end
    try
       if disc
          if !adjA & !adjB
             sylvds!(A, B, C, adjA = false, adjB = false)
          elseif !adjA & adjB
             sylvds!(A, B.parent, C, adjA = false, adjB = true)
          elseif adjA & !adjB
             sylvds!(A.parent, B, C, adjA = true, adjB = false)
          else
             sylvds!(A.parent, B.parent, C, adjA = true, adjB = true)
          end
          return C[:]
       else
         if !adjA & !adjB
            sylvcs!(A, B, C, adjA = false, adjB = false)
         elseif !adjA & adjB
            sylvcs!(A, B.parent, C, adjA = false, adjB = true)
         elseif adjA & !adjB
            sylvcs!(A.parent, B, C, adjA = true, adjB = false)
         else
            sylvcs!(A.parent, B.parent, C, adjA = true, adjB = true)
         end
         return C[:]
       end
    catch err
      # if isnothing(findfirst("LAPACKException",string(err)))
      if findfirst("LAPACKException",string(err)) === nothing
         rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
    end
  end
  function tprod(x)
    T1 = promote_type(T, eltype(x))
    C = copy(reshape(convert(Vector{T1}, x), m, n))
    if T !== T1
      if cmplx
         adjA ? A = convert(Matrix{T1},A.parent)' : A = convert(Matrix{T1},A)
         adjB ? B = convert(Matrix{T1},B.parent)' : B = convert(Matrix{T1},B)
      else
        T1r = real(T1)
        if T1r !== T
          adjA ? A = convert(Matrix{T1r},A.parent)' : A = convert(Matrix{T1r},A)
          adjB ? B = convert(Matrix{T1r},B.parent)' : B = convert(Matrix{T1r},B)
        end
      end
    end
    try
       if disc
          if !adjA & !adjB
             sylvds!(A, B, C, adjA = true, adjB = true)
          elseif !adjA & adjB
             sylvds!(A, B.parent, C, adjA = true, adjB = false)
          elseif adjA & !adjB
             sylvds!(A.parent, B, C, adjA = false, adjB = true)
          else
             sylvds!(A.parent, B.parent, C, adjA = false, adjB = false)
          end
          return C[:]
       else
         if !adjA & !adjB
            sylvcs!(A, B, C, adjA = true, adjB = true)
         elseif !adjA & adjB
            sylvcs!(A, B.parent, C, adjA = true, adjB = false)
         elseif adjA & !adjB
            sylvcs!(A.parent, B, C, adjA = false, adjB = true)
         else
            sylvcs!(A.parent, B.parent, C, adjA = false, adjB = false)
         end
         return C[:]
         end
    catch err
        # if isnothing(findfirst("LAPACKException",string(err)))
        if findfirst("LAPACKException",string(err)) === nothing
           rethrow()
        else
           throw("ME:SingularException: Singular operator")
        end
     end
  end
  function ctprod(x)
    T1 = promote_type(T, eltype(x))
    C = copy(reshape(convert(Vector{T1}, x), m, n))
    if T !== T1
      if cmplx
         adjA ? A = convert(Matrix{T1},A.parent)' : A = convert(Matrix{T1},A)
         adjB ? B = convert(Matrix{T1},B.parent)' : B = convert(Matrix{T1},B)
      else
        T1r = real(T1)
        if T1r !== T
          adjA ? A = convert(Matrix{T1r},A.parent)' : A = convert(Matrix{T1r},A)
          adjB ? B = convert(Matrix{T1r},B.parent)' : B = convert(Matrix{T1r},B)
        end
      end
    end
    try
       if disc
          if !adjA & !adjB
             sylvds!(A, B, C, adjA = true, adjB = true)
          elseif !adjA & adjB
             sylvds!(A, B.parent, C, adjA = true, adjB = false)
          elseif adjA & !adjB
             sylvds!(A.parent, B, C, adjA = false, adjB = true)
          else
             sylvds!(A.parent, B.parent, C, adjA = false, adjB = false)
          end
          return C[:]
       else
         if !adjA & !adjB
            sylvcs!(A, B, C, adjA = true, adjB = true)
         elseif !adjA & adjB
            sylvcs!(A, B.parent, C, adjA = true, adjB = false)
         elseif adjA & !adjB
            sylvcs!(A.parent, B, C, adjA = false, adjB = true)
         else
            sylvcs!(A.parent, B.parent, C, adjA = false, adjB = false)
         end
         return C[:]
       end
    catch err
        # if isnothing(findfirst("LAPACKException",string(err)))
        if findfirst("LAPACKException",string(err)) === nothing
           rethrow()
        else
           throw("ME:SingularException: Singular operator")
        end
     end
  end
  F1 = typeof(prod)
  F2 = typeof(tprod)
  F3 = typeof(ctprod)
  return LinearOperator{T}(m * n, n * m, false, false, prod, tprod, ctprod)
end
invsylvsop(A :: Schur, B :: Schur; disc = false) = invsylvsop(A.T,B.T,disc = disc)
"""
    M = sylvop(A, B, C, D)

Define the generalized Sylvester operator `M: X -> AXB+CXD`, where `(A,C)` and `(B,D)` a pairs of square matrices.
"""
function sylvop(A, B, C, D)
  m = LinearAlgebra.checksquare(A)
  n = LinearAlgebra.checksquare(B)
  if [m; n] != LinearAlgebra.checksquare(C,D)
     throw(DimensionMismatch("A, B, C and D have incompatible dimensions"))
  end
  T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D))
  function prod(x)
    T1 = promote_type(T, eltype(x))
    X = reshape(convert(Vector{T1}, x), m, n)
    return (A * X * B + C * X * D)[:]
  end
  function tprod(x)
    T1 = promote_type(T, eltype(x))
    X = reshape(convert(Vector{T1}, x), m, n)
    return (adjoint(A) * X * adjoint(B) + adjoint(C) * X * adjoint(D) )[:]
  end
  function ctprod(x)
    T1 = promote_type(T, eltype(x))
    X = reshape(convert(Vector{T1}, x), m, n)
    return (A' * X * B' + C' * X * D' )[:]
  end
  F1 = typeof(prod)
  F2 = typeof(tprod)
  F3 = typeof(ctprod)
  return LinearOperator{T}(m * n, n * m, false, false, prod, tprod, ctprod)
end
"""
    MINV = invsylvop(A, B, C, D)

Define MINV, the inverse of the generalized Sylvester operator `M: X -> AXB+CXD`,
where (A,C) and (B,D) a pairs of square matrices.
"""
function invsylvop(A, B, C, D)
  m = LinearAlgebra.checksquare(A)
  n = LinearAlgebra.checksquare(B)
  if [m; n] != LinearAlgebra.checksquare(C,D)
     throw(DimensionMismatch("A, B, C and D have incompatible dimensions"))
  end
  T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D))
  function prod(x)
    T1 = promote_type(T, eltype(x))
    E = reshape(convert(Vector{T1}, x), m, n)
    try
       return gsylv(A,B,C,D,E)[:]
    catch err
       # if isnothing(findfirst("SingularException",string(err)))
       if findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
    end
  end
  function tprod(x)
   T1 = promote_type(T, eltype(x))
   E = reshape(convert(Vector{T1}, x), m, n)
    try
       return gsylv(A',B',C',D',E)[:]
    catch err
       # if isnothing(findfirst("SingularException",string(err)))
       if findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
    end
  end
  function ctprod(x)
    T1 = promote_type(T, eltype(x))
    E = reshape(convert(Vector{T1}, x), m, n)
    try
       return gsylv(A',B',C',D',E)[:]
    catch err
       # if isnothing(findfirst("SingularException",string(err)))
       if findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
    end
  end
  F1 = typeof(prod)
  F2 = typeof(tprod)
  F3 = typeof(ctprod)
  return LinearOperator{T}(m * n, n * m, false, false, prod, tprod, ctprod)
end
"""
    MINV = invsylvsop(A, B, C, D; DBSchur = false)

Define MINV, the inverse of the generalized Sylvester operator `M: X -> AXB+CXD`,
with the pairs `(A,C)` and `(B,D)` in generalized Schur forms. If `DBSchur = true`,
the pair `(D,B)` is in generalized Schur form.
"""
function invsylvsop(A, B, C, D; DBSchur = false)
  m = LinearAlgebra.checksquare(A)
  n = LinearAlgebra.checksquare(B)
  if [m; n] != LinearAlgebra.checksquare(C,D)
     throw(DimensionMismatch("A, B, C and D have incompatible dimensions"))
  end
  T = promote_type(eltype(A),eltype(B),eltype(C),eltype(D))
  if !(T <: BlasFloat)
     T = promote_type(Float64,T)
  end
  adjA = isa(A,Adjoint)
  adjB = isa(B,Adjoint)
  adjC = isa(C,Adjoint)
  adjD = isa(D,Adjoint)
  if adjA !== adjC
     error("Only calls with pairs (A,C) or (A',C') are allowed")
  end
  if adjB !== adjD
     error("Only calls with pairs (B,D) or (B',D') are allowed")
  end
  if eltype(A) !== T
    adjA ? A = convert(Matrix{T},A.parent)'  : A = convert(Matrix{T},A)
  end
  if eltype(B) !== T
    adjB ? B = convert(Matrix{T},B.parent)' :  B = convert(Matrix{T},B)
  end
  if eltype(C) !== T
   adjC ? C = convert(Matrix{T},C.parent)'  : C = convert(Matrix{T},C)
  end
  if eltype(D) !== T
    adjD ? D = convert(Matrix{T},D.parent)' :  D = convert(Matrix{T},D)
  end
  cmplx = T<:Complex
  adjAC = adjA & adjC
  if adjAC
     if !isschur(A.parent,C.parent)
         error("The pair (A,C) must be in generalized Schur form")
     end
  else
     if !isschur(A,C)
        error("The pair (A,C) must be in generalized Schur form")
     end
  end
  adjBD = adjB & adjD
  if adjBD
     if DBSchur
       if !isschur(D.parent, B.parent)
           error("The pair (D,B) must be in generalized Schur form")
       end
     else
        if !isschur(B.parent, D.parent)
            error("The pair (B,D) must be in generalized Schur form")
        end
     end
  else
     if DBSchur
        if !isschur(D,B)
           error("The pair (D,B) must be in generalized Schur form")
        end
     else
        if !isschur(B,D)
           error("The pair (B,D) must be in generalized Schur form")
        end
     end
  end
  function prod(x)
    T1 = promote_type(T, eltype(x))
    Y = copy(reshape(convert(Vector{T1}, x), m, n))
    if T !== T1
      if cmplx
         adjA ? A = convert(Matrix{T1},A.parent)' : A = convert(Matrix{T1},A)
         adjB ? B = convert(Matrix{T1},B.parent)' : B = convert(Matrix{T1},B)
         adjC ? C = convert(Matrix{T1},C.parent)' : C = convert(Matrix{T1},C)
         adjD ? D = convert(Matrix{T1},D.parent)' : D = convert(Matrix{T1},D)
         else
        T1r = real(T1)
        if T1r !== T
          adjA ? A = convert(Matrix{T1r},A.parent)' : A = convert(Matrix{T1r},A)
          adjB ? B = convert(Matrix{T1r},B.parent)' : B = convert(Matrix{T1r},B)
          adjC ? C = convert(Matrix{T1r},C.parent)' : C = convert(Matrix{T1r},C)
          adjD ? D = convert(Matrix{T1r},D.parent)' : D = convert(Matrix{T1r},D)
        end
      end
    end
    try
       if !adjAC & !adjBD
          gsylvs!(A, B, C, D, Y, adjAC = false, adjBD = false, DBSchur = DBSchur)
       elseif !adjAC & adjBD
          gsylvs!(A, B.parent, C, D.parent, Y, adjAC = false, adjBD = true, DBSchur = DBSchur)
       elseif adjAC & !adjBD
          gsylvs!(A.parent, B, C.parent, D, Y, adjAC = true, adjBD = false, DBSchur = DBSchur)
       else
          gsylvs!(A.parent, B.parent, C.parent, D.parent, Y, adjAC = true, adjBD = true, DBSchur = DBSchur)
       end
       return Y[:]
    catch err
       # if isnothing(findfirst("SingularException",string(err)))
       if findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
    end
  end
  function tprod(x)
    T1 = promote_type(T, eltype(x))
    Y = copy(reshape(convert(Vector{T1}, x), m, n))
    if T !== T1
      if cmplx
         adjA ? A = convert(Matrix{T1},A.parent)' : A = convert(Matrix{T1},A)
         adjB ? B = convert(Matrix{T1},B.parent)' : B = convert(Matrix{T1},B)
         adjC ? C = convert(Matrix{T1},C.parent)' : C = convert(Matrix{T1},C)
         adjD ? D = convert(Matrix{T1},D.parent)' : D = convert(Matrix{T1},D)
         else
        T1r = real(T1)
        if T1r !== T
          adjA ? A = convert(Matrix{T1r},A.parent)' : A = convert(Matrix{T1r},A)
          adjB ? B = convert(Matrix{T1r},B.parent)' : B = convert(Matrix{T1r},B)
          adjC ? C = convert(Matrix{T1r},C.parent)' : C = convert(Matrix{T1r},C)
          adjD ? D = convert(Matrix{T1r},D.parent)' : D = convert(Matrix{T1r},D)
        end
      end
    end
    try
       if !adjAC & !adjBD
          gsylvs!(A, B, C, D, Y, adjAC = true, adjBD = true, DBSchur = DBSchur)
       elseif !adjAC & adjBD
          gsylvs!(A, B.parent, C, D.parent, Y, adjAC = true, adjBD = false, DBSchur = DBSchur)
       elseif adjAC & !adjBD
          gsylvs!(A.parent, B, C.parent, D, Y, adjAC = false, adjBD = true, DBSchur = DBSchur)
       else
          gsylvs!(A.parent, B.parent, C.parent, D.parent, Y, adjAC = false, adjBD = false, DBSchur = DBSchur)
       end
       return Y[:]
    catch err
       # if isnothing(findfirst("SingularException",string(err)))
       if findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
    end
  end
  function ctprod(x)
    T1 = promote_type(T, eltype(x))
    Y = copy(reshape(convert(Vector{T1}, x), m, n))
    if T !== T1
      if cmplx
         adjA ? A = convert(Matrix{T1},A.parent)' : A = convert(Matrix{T1},A)
         adjB ? B = convert(Matrix{T1},B.parent)' : B = convert(Matrix{T1},B)
         adjC ? C = convert(Matrix{T1},C.parent)' : C = convert(Matrix{T1},C)
         adjD ? D = convert(Matrix{T1},D.parent)' : D = convert(Matrix{T1},D)
         else
        T1r = real(T1)
        if T1r !== T
          adjA ? A = convert(Matrix{T1r},A.parent)' : A = convert(Matrix{T1r},A)
          adjB ? B = convert(Matrix{T1r},B.parent)' : B = convert(Matrix{T1r},B)
          adjC ? C = convert(Matrix{T1r},C.parent)' : C = convert(Matrix{T1r},C)
          adjD ? D = convert(Matrix{T1r},D.parent)' : D = convert(Matrix{T1r},D)
        end
      end
    end
    try
       if !adjAC & !adjBD
          gsylvs!(A, B, C, D, Y, adjAC = true, adjBD = true, DBSchur = DBSchur)
       elseif !adjAC & adjBD
          gsylvs!(A, B.parent, C, D.parent, Y, adjAC = true, adjBD = false, DBSchur = DBSchur)
       elseif adjAC & !adjBD
          gsylvs!(A.parent, B, C.parent, D, Y, adjAC = false, adjBD = true, DBSchur = DBSchur)
       else
          gsylvs!(A.parent, B.parent, C.parent, D.parent, Y, adjAC = false, adjBD = false, DBSchur = DBSchur)
       end
       return Y[:]
    catch err
       # if isnothing(findfirst("SingularException",string(err)))
       if findfirst("SingularException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
    end
  end
  F1 = typeof(prod)
  F2 = typeof(tprod)
  F3 = typeof(ctprod)
  return LinearOperator{T}(m * n, n * m, false, false, prod, tprod, ctprod)
end
invsylvsop(AC :: GeneralizedSchur, BD :: GeneralizedSchur) = invsylvsop(AC.S,BD.S,AC.T,BD.T)

"""
    M = sylvsysop(A, B, C, D)

Define the operator `M: (X,Y) -> (AX+YB, CX+YD )`,
where `(A,C)` and `(B,D)` a pairs of square matrices.
"""
function sylvsysop(A, B, C, D)
  m = LinearAlgebra.checksquare(A)
  n = LinearAlgebra.checksquare(B)
  T = promote_type(eltype(A), eltype(B))
  if [m; n] != LinearAlgebra.checksquare(C,D)
     throw(DimensionMismatch("A, B, C and D have incompatible dimensions"))
  end
  mn = m*n
  function prod(x)
    T1 = promote_type(T, eltype(x))
    X = reshape(convert(Vector{T1}, x[1:mn]), m, n)
    Y = reshape(convert(Vector{T1}, x[mn+1:2*mn]), m, n)
    return ([A * X + Y * B C * X + Y * D])[:]
  end
  function tprod(x)
    T1 = promote_type(T, eltype(x))
    X = reshape(convert(Vector{T1}, x[1:mn]), m, n)
    Y = reshape(convert(Vector{T1}, x[mn+1:2*mn]), m, n)
    return [adjoint(A) * X + adjoint(C) * Y  X * adjoint(B) + Y * adjoint(D)][:]
  end
  function ctprod(x)
    T1 = promote_type(T, eltype(x))
    X = reshape(convert(Vector{T1}, x[1:mn]), m, n)
    Y = reshape(convert(Vector{T1}, x[mn+1:2*mn]), m, n)
    return [A' * X + C' * Y  X * B' + Y * D'][:]
  end
  F1 = typeof(prod)
  F2 = typeof(tprod)
  F3 = typeof(ctprod)
  return LinearOperator{T}(2*mn, 2*mn, false, false, prod, tprod, ctprod)
end
"""
    MINV = invsylvsysop(A, B, C, D)

Define MINV, the inverse of the linear operator `M: (X,Y) -> (AX+YB, CX+YD )`,
where `(A,C)` and `(B,D)` a pairs of square matrices.
"""
function invsylvsysop(A, B, C, D)
  m = LinearAlgebra.checksquare(A)
  n = LinearAlgebra.checksquare(B)
  if [m; n] != LinearAlgebra.checksquare(C,D)
     throw(DimensionMismatch("A, B, C and D have incompatible dimensions"))
  end
  T = promote_type(eltype(A), eltype(B), eltype(C), eltype(D))
  mn = m*n
  function prod(x)
    T1 = promote_type(T, eltype(x))
    E = reshape(convert(Vector{T1}, x[1:mn]), m, n)
    F = reshape(convert(Vector{T1}, x[mn+1:2*mn]), m, n)
    try
       (X,Y) = sylvsys(A,B,E,C,D,F)
       return [X Y][:]
    catch err
       # if isnothing(findfirst("LAPACKException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
    end
  end
  function tprod(x)
    T1 = promote_type(T, eltype(x))
    E = reshape(convert(Vector{T1}, x[1:mn]), m, n)
    F = reshape(convert(Vector{T1}, x[mn+1:2*mn]), m, n)
    try
       (X,Y) = dsylvsys(A',B',E,C',D',F)[:]
       return [X Y][:]
    catch err
       # if isnothing(findfirst("LAPACKException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
    end
  end
  function ctprod(x)
    T1 = promote_type(T, eltype(x))
    E = reshape(convert(Vector{T1}, x[1:mn]), m, n)
    F = reshape(convert(Vector{T1}, x[mn+1:2*mn]), m, n)
    try
       (X,Y) = dsylvsys(A',B',E,C',D',F)[:]
       return [X Y][:]
    catch err
       # if isnothing(findfirst("LAPACKException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
    end
  end
  F1 = typeof(prod)
  F2 = typeof(tprod)
  F3 = typeof(ctprod)
  return LinearOperator{T}(2*mn, 2*mn, false, false, prod, tprod, ctprod)
end
"""
    MINV = invsylvsyssop(A, B, C, D)

Define MINV, the inverse of the linear operator `M: (X,Y) -> (AX+YB, CX+YD)`,
with the pairs `(A,C)` and `(B,D)` in generalized Schur forms.
"""
function invsylvsyssop(A, B, C, D)
  m = LinearAlgebra.checksquare(A)
  n = LinearAlgebra.checksquare(B)
  if [m; n] != LinearAlgebra.checksquare(C,D)
     throw(DimensionMismatch("A, B, C and D have incompatible dimensions"))
  end
  T = promote_type(eltype(A),eltype(B),eltype(C),eltype(D))
  if !(T <: BlasFloat)
     T = promote_type(Float64,T)
  end
  cmplx = T<:Complex
  if isa(A,Adjoint) || isa(B,Adjoint) || isa(C,Adjoint)  || isa(D,Adjoint)
     error("Only calls with (A, B, C, D) without adjoints are allowed")
  end
  if eltype(A) !== T
     A = convert(Matrix{T},A)
  end
  if eltype(B) !== T
     B = convert(Matrix{T},B)
  end
  if eltype(C) !== T
     C = convert(Matrix{T},C)
  end
  if eltype(D) !== T
     D = convert(Matrix{T},D)
  end
  if !isschur(A,C)
     error("The pair (A,C) must be in generalized Schur form")
  end
  if !isschur(B,D)
     error("The pair (B,D) must be in generalized Schur form")
  end
  cmplx ? TA = 'C' : TA = 'T'
  mn = m*n
  function prod(x)
    T1 = promote_type(T, eltype(x))
    E = reshape(convert(Vector{T1}, x[1:mn]), m, n)
    F = reshape(convert(Vector{T1}, x[mn+1:2*mn]), m, n)
    if T !== T1
      if cmplx
         A = convert(Matrix{T1},A)
         B = convert(Matrix{T1},B)
         C = convert(Matrix{T1},C)
         D = convert(Matrix{T1},D)
      else
        T1r = real(T1)
        if T1r !== T
         A = convert(Matrix{T1r},A)
         B = convert(Matrix{T1r},B)
         C = convert(Matrix{T1r},C)
         D = convert(Matrix{T1r},D)
        end
      end
    end
    try
      X, Y = sylvsyss!(A,B,E,C,D,F)
      return [X Y][:]
    catch err
       # if isnothing(findfirst("LAPACKException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
    end
  end
  function tprod(x)
    T1 = promote_type(T, eltype(x))
    E = reshape(convert(Vector{T1}, x[1:mn]), m, n)
    F = reshape(convert(Vector{T1}, x[mn+1:2*mn]), m, n)
    if T !== T1
      if cmplx
         A = convert(Matrix{T1},A)
         B = convert(Matrix{T1},B)
         C = convert(Matrix{T1},C)
         D = convert(Matrix{T1},D)
      else
        T1r = real(T1)
        if T1r !== T
         A = convert(Matrix{T1r},A)
         B = convert(Matrix{T1r},B)
         C = convert(Matrix{T1r},C)
         D = convert(Matrix{T1r},D)
        end
      end
    end
    try
      X, Y = dsylvsyss!(A,B,E,C,D,F)
      return [X Y][:]
    catch err
       # if isnothing(findfirst("LAPACKException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
    end
  end
  function ctprod(x)
    T1 = promote_type(T, eltype(x))
    E = reshape(convert(Vector{T1}, x[1:mn]), m, n)
    F = reshape(convert(Vector{T1}, x[mn+1:2*mn]), m, n)
    if T !== T1
      if cmplx
         A = convert(Matrix{T1},A)
         B = convert(Matrix{T1},B)
         C = convert(Matrix{T1},C)
         D = convert(Matrix{T1},D)
      else
        T1r = real(T1)
        if T1r !== T
         A = convert(Matrix{T1r},A)
         B = convert(Matrix{T1r},B)
         C = convert(Matrix{T1r},C)
         D = convert(Matrix{T1r},D)
        end
      end
    end
    try
      X, Y = dsylvsyss!(A,B,E,C,D,F)
      return [X Y][:]
    catch err
       # if isnothing(findfirst("LAPACKException",string(err)))
       if findfirst("LAPACKException",string(err)) === nothing
          rethrow()
       else
          throw("ME:SingularException: Singular operator")
       end
    end
  end
  F1 = typeof(prod)
  F2 = typeof(tprod)
  F3 = typeof(ctprod)
  return LinearOperator{T}(2*mn, 2*mn, false, false, prod, tprod, ctprod)
end
invsylvsyssop(AC :: GeneralizedSchur, BD :: GeneralizedSchur) = invsylvsyssop(AC.S,BD.S,AC.T,BD.T)


text:
@testset "eigenvalues/eigenvectors: $MatT" for MatT in (AcbMatrix, AcbRefMatrix)
    A = [
        0.6873474041954415 0.7282180564881044 0.07360652513458521
        0.000835810121029068 0.9256166870757694 0.5363310989411239
        0.07387174694790022 0.4050436025621329 0.20226010388885896
    ]
    B = [
        0.8982563031334123 0.3029712969740874 0.8585014523679579
        0.7583002736998279 0.8854763478184455 0.3031103325817668
        0.2319572749472405 0.5769840251057949 0.5119507333628952
    ]
    M = MatT(A + B * im, prec = 64)
    # M = MatT(rand(3, 3) + im * rand(3, 3), prec = 64)

    VecT = typeof(similar(M, 3))

    @testset "approx_eig_qr" begin
        λs_a_r, revs_a = Arblib.approx_eig_qr(M, side = :right)
        @test λs_a_r isa VecT
        @test revs_a isa MatT

        ε = 1e-10

        λs_a_l, revs_a = Arblib.approx_eig_qr(M, tol = Arblib.Mag(ε), side = :left)
        @test λs_a_l isa VecT
        @test revs_a isa MatT

        @test all(abs.(λs_a_r - λs_a_l) .< ε)

        λs_r = similar(M, size(M, 1))
        Arblib.approx_eig_qr!(λs_r, M)
        @test Arblib.is_zero(λs_r - λs_a_r, length(λs_r))
    end

    @testset "eig_simple" begin
        λs1, _ = Arblib.eig_simple_vdhoeven_mourrain(M, side = :right)
        λs2, _ = Arblib.eig_simple_vdhoeven_mourrain(M, side = :left)

        @test all(Arblib.containszero, λs1 - λs2)

        # λs1, _ = Arblib.eig_simple_rump(M, side = :right)
        # segfaults in acb_mat_solve at /workspace/srcdir/arb-2.18.1/acb_mat/solve.c:17
        # Issue #321 in Arblib (fixed by #330)
        # λs2, _ = Arblib.eig_simple_rump(M, side=:left)
        # @test all(Arblib.containszero, λs1 - λs2)

        λs1, _ = Arblib.eig_simple(M, side = :right)
        λs2, _ = Arblib.eig_simple(M, side = :left)
        @test all(Arblib.containszero, λs1 - λs2)

        λs = similar(M, size(M, 1))
        Arblib.eig_simple_vdhoeven_mourrain!(λs, M)
        @test all(Arblib.containszero, λs - λs1)

        λs = similar(M, size(M, 1))
        Arblib.eig_simple_rump!(λs, M)
        @test all(Arblib.containszero, λs - λs1)

        λs = similar(M, size(M, 1))
        Arblib.eig_simple!(λs, M)
        @test all(Arblib.containszero, λs - λs1)
    end

    N = similar(M)
    evs = [Arb(2.0), Arb(2.0), Arb(rand())]
    N[1, 1], N[2, 2], N[3, 3] = evs

    N = M * N * M^-1

    @testset "enclosures" begin
        ε = Arblib.Mag()
        tol = 1e-12
        λ_approx, R_approx = Arblib.approx_eig_qr(M, tol = tol)

        @test Arblib.eig_global_enclosure!(
            ε,
            M,
            λ_approx,
            R_approx;
            prec = precision(M),
        ) isa Arblib.Mag

        @test ε <= tol

        λs = similar(M, size(M, 1))
        Arblib.eig_simple!(λs, M, λ_approx, R_approx)

        for λa in λ_approx
            a_real = let x = real(λa)
                m = Arblib.midref(x)
                r = Arblib.radref(x)
                Arblib.set_interval!(x, m - (r + ε), m + (r + ε))
            end

            a_imag = let x = imag(λa)
                m = Arblib.midref(x)
                r = Arblib.radref(x)
                Arblib.set_interval!(x, m - (r + ε), m + (r + ε))
            end

            a = Acb(a_real, a_imag)
            @test any(Arblib.containszero(a - λ) for λ in λs)
        end

        @test_throws Arblib.EigenvalueComputationError Arblib.eig_simple(N)

        λ_approx, R_approx = Arblib.approx_eig_qr(N)
        v = sortperm(λ_approx, by = abs, rev = true)

        λ = Acb(prec = precision(N))
        R = similar(N, (3, 1))
        Arblib.eig_enclosure_rump!(λ, R, N, λ_approx[v[1]], R_approx[:, v[1:1]])
        @test !isfinite(λ)

        λ = Acb(prec = precision(N))
        R = similar(N, (3, 1))
        Arblib.eig_enclosure_rump!(λ, R, N, λ_approx[v[3]], R_approx[:, v[3:3]])
        @test isfinite(λ)
        @test Arblib.contains_zero(λ - evs[3])

        λ = Acb(prec = precision(N))
        R = similar(N, (3, 2))
        Arblib.eig_enclosure_rump!(λ, R, N, λ_approx[v[1]], R_approx[:, v[1:2]])
        @test isfinite(λ)
        @test Arblib.contains_zero(λ - 2)
        @test all(Arblib.contains_zero.(N * R - R * λ))

        λ = Acb(prec = precision(N))
        R = similar(N, (3, 2))
        J = similar(N, (2, 2))

        Arblib.eig_enclosure_rump!(λ, J, R, N, λ_approx[v[2]], R_approx[:, v[1:2]])
        @test Arblib.contains_zero(λ - 2)
        @test all(Arblib.contains_zero.(N * R - R * J))
    end

    @testset "eig_multiple" begin
        λs = similar(N, 3)
        @test Arblib.eig_multiple_rump(N) isa VecT
        Arblib.eig_multiple_rump!(λs, N)

        v = sortperm(λs, by = abs, rev = true)
        @test all(Arblib.contains_zero.(λs[v] - evs))

        @test Arblib.eig_multiple(N) isa VecT
        Arblib.eig_multiple!(λs, N)

        v = sortperm(λs, by = abs, rev = true)
        @test all(Arblib.contains_zero.(λs[v] - evs))
    end
end

ext: jl


text:
export Connect

function Connect(model;pre=[],post=[],condition="true",expr="",parallel=false)
  var=match(r"^.*(?==)",expr).match
  value=match(r"(?<==).*",expr).match
  if !contains(model.connection,var)
    N=model.groups[:N]
    model.connection*="$var=SharedArray(Float64,N,N)\n"
  end

   if pre==[] || post==[]
      model.connection*="$var=SharedArray(Float64,(N,N),init=S->S[localindexes(S)]=$value)\n"
  elseif pre=="all" && post=="all"
    model.connection*="for j=1:N\n for i=1:N\n if $condition\n $var[i,j]=$value\n end\n end\n end\n"
  elseif pre=="all" && post!="all"
    for postTemp in [post;]
      i=model.groups[postTemp][:localindex]
      ii="(i-$(i[1]-1))"
      condition=replace(condition,"i",ii)
      value=replace(value,"i",ii)
      model.connection*="for j=1:N\n for i=$i\n if $condition\n $var[i,j]=$value\n end\n end\n end\n"
    end
  elseif pre!="all" && post=="all"
    for preTemp in [pre;]
      j=model.groups[preTemp][:localindex]
      jj="(j-$(j[1]-1))"
      condition=replace(condition,"j",jj)
      value=replace(value,"j",jj)
      model.connection*="for j=$j\n for i=1:N\n if $condition\n $var[i,j]=$value\n end\n end\n end\n"
    end
  else
    for preTemp in [pre;],postTemp in [post;]
      i=model.groups[postTemp][:localindex]
      j=model.groups[preTemp][:localindex]
      ii="(i-$(i[1]-1))"
      jj="(j-$(j[1]-1))"
      condition=replace(condition,"i",ii);condition=replace(condition,"j",jj)
      value=replace(value,"i",ii);value=replace(value,"j",jj)
      model.connection*="for j=$j\n for i=$i\n if $condition\n $var[i,j]=$value\n end\n end\n end\n"
    end
  end
end # end of function Connect


text:
using LinearAlgebra

# -------- Wachter et al (2018):
struct GenericGenerator <: AbstractGradientBasedGenerator
    loss::Symbol # loss function
    complexity::Function # complexity function
    λ::AbstractFloat # strength of penalty
    ϵ::AbstractFloat # step size
    τ::AbstractFloat # tolerance for convergence
end

"""
    GenericGenerator(
        ;
        loss::Symbol=:logitbinarycrossentropy,
        complexity::Function=norm,
        λ::AbstractFloat=0.1,
        ϵ::AbstractFloat=0.1,
        τ::AbstractFloat=1e-5
    )

An outer constructor method that instantiates a generic generator.

# Examples
```julia-repl
generator = GenericGenerator()
```
"""
GenericGenerator(
    ;
    loss::Symbol=:logitbinarycrossentropy,
    complexity::Function=norm,
    λ::AbstractFloat=0.1,
    ϵ::AbstractFloat=0.1,
    τ::AbstractFloat=1e-5
) = GenericGenerator(loss, complexity, λ, ϵ, τ)


text:
"""
basic2d_test()

Test nsol with the simple 2D problem.

"""
function basic2d_test()
    x0 = ones(2, 1)
    fv = zeros(2, 1)
    jv = zeros(2, 2)
    jsv = zeros(Float32, 2, 2)
    # local convergence testing
    #
    # single vs double Jacobian
    # jfact=nofact should have no effect on the nonlinear iteration
    #
    nout = nsol(basic2d!, x0, fv, jv; rtol = 1.e-10, sham = 1)
    sout = nsol(basic2d!, x0, fv, jsv; sham = 1, jfact = nofact)
    dss = norm(nout.solution - sout.solution)
    hss = norm(nout.history - sout.history)
    singleok = (norm(dss) < 1.e-7) && (norm(hss) < 1.e-7)
    if ~singleok
        println("Single/Double test fails.")
    end
    #
    # chord vs Newton
    #
    cout = nsol(basic2d!, x0, fv, jv; solver = "chord")
    dsc = norm(nout.solution - cout.solution)
    lch = length(cout.history)
    lnh = length(nout.history)
    jevals = sum(cout.stats.ijac)
    chordok = (dsc < 1.e-6) && (lch == 17) && (lnh == 5) && (jevals == 1)
    if ~chordok
        println("Chord/Newton test fails")
    end
    #
    # Analytic vs finite-difference Jacobian
    #
    eout = nsol(basic2d!, x0, fv, jv, jbasic2d!; sham = 1)
    fdok =
        (norm(eout.history - nout.history) < 1.e-6) &&
        (norm(eout.solution - nout.solution) < 1.e-10)
    if ~fdok
        println("FD/Analytic test fails.")
    end
    #
    # Shamanskii
    #
    s1out = nsol(basic2d!, x0, fv, jv; sham = 2, rtol = 1.e-10)
    dout1 = norm(s1out.solution - nout.solution)
    jevals1 = sum(s1out.stats.ijac)
    s2out = nsol(basic2d!, x0, fv, jv; sham = 2, rtol = 1.e-10, resdec = 0.5)
    jevals2 = sum(s2out.stats.ijac)
    dout2 = norm(s2out.solution - nout.solution)
    shamok = (dout1 < 1.e-10) && (dout2 < 1.e-10) && (jevals1 == 4) && (jevals2 == 3)
    if ~shamok
        println("Shamanskii test fails.")
    end
    #
    # Global convergence
    #
    x0a = [2.0, 0.5]
    FS = zeros(2)
    FPS = zeros(2, 2)
    FPSS = zeros(Float32, 2, 2)
    nouta = nsol(simple!, x0a, FS, FPS; keepsolhist = true, sham = 1)
    noutb = nsol(simple!, x0a, FS, FPSS, jsimple!; keepsolhist = true, sham = 1)
    noutc = nsol(simple!, x0a, FS, FPSS, jsimple!; armmax = 0, sham = 1)
    iarm = nouta.stats.iarm
    iarm2 = noutc.idid
    armok = (iarm[2] == 2) && ~iarm2
    preok = (norm(noutb.solhist - nouta.solhist, Inf) < 1.e-6)
    solok = (norm(noutb.solution - nouta.solution, Inf) < 1.e-10)
    globok = armok && preok && solok
    if ~globok
        println("Global test fails.")
    end
    return chordok && singleok && fdok && shamok && globok
end

function basic2d!(FV, x)
    FV[1] = x[1] * x[1] - 2.0
    FV[2] = exp(x[1] - 1) + x[2] * x[2] - 2.0
    return FV
end

function jbasic2d!(JV, FV, x)
    JV[1, 1] = 2 * x[1]
    JV[1, 2] = 0.0
    JV[2, 1] = exp(x[1] - 1)
    JV[2, 2] = 2 * x[2]
    return JV
end


text:
VERSION >= v"0.4.0-dev+6641" && __precompile__()

module Faker

#export
#bothify
  #email
#state_abbr
#latitude
#street_name
#address
#street_address
#postcode
#longitude
#country
#secondary_address
#street_suffix
#city_prefix
#city_suffix
#building_number
#city
#state
#company
#company_suffix
#catch_phrase
#bs
#date_time_ad
#month
#am_pm
#iso8601
#date_time_this_century
#date_time
#month_name
#date_time_this_month
#date_time_this_decade
#day_of_week
#day_of_month
#unix_time
#date_time_this_year
#timezone
#time
#date_time_between
#geo_coordinate
#century
#year
#email
#free_email
#safe_email
#user_name
#domain_name
#domain_word
#domain_suffix
#ip_v4_address
#ip_v6_address
#uri
#url
#text
#sentence
#word
#paragraphs
#words
#paragraph
#sentences
#phone_number
#password
#locale
#md5
#sha1
#null_boolean
#sha256
#country_code
#language_code
#boolean

importall Base
import Base.Callable
import Base: getindex
#:  #isequal isfinite isinf, isinteger,
            # isless, isnan, isnull, isnumber, isreal, issubnormal,
            # typemin, typemax

# ---------------------------------------------------------------------------------------
include("creditcard.jl")
include("persona.jl")
include("direccion.jl")
include("trabajo.jl")
include("datatime.jl")
include("phone_number.jl")
include("color.jl")
include("base.jl")
include("company.jl")
include("file.jl")
include("internet.jl")
include("lorem.jl")
include("misc.jl")
include("profile.jl")
include("country.jl")
include("user_agent.jl")

end # module


text:
# Run the examples
using MatLang
examplePath = joinpath( dirname(dirname(@__DIR__)), "examples")
include(joinpath(examplePath, "Language_Fundamentals", "usage_Entering_Commands.jl"))
include(joinpath(examplePath, "Language_Fundamentals", "usage_Matrices_and_Arrays.jl"))
include(joinpath(examplePath, "Language_Fundamentals", "Data_Types", "usage_Numeric_Types.jl"))

max_stars_repo_name: juliamatlab/MatLang

max_issues_repo_name: juliamatlab/MatLang

max_forks_repo_name: juliamatlab/MatLang


text:
module ElectricalVehicleAndCharger

println("ElectricalVehicleAndCharger: Demonstrates the ability to change models from Julia.")

using Modia
using Modia.Electric

using Modia:@equation
using Modia:addEquation!
using Modia:deleteEquation!

# Desired:
#   using ModiaMath: plot
#
# In order that ModiaMath need not to be defined in the user environment, it is included via Modia:
using Modia.ModiaMath: plot

# -----------------------------------

function concatenateTimeSeries(result1, result2)
  result = result1
  for (name, values) in result
    result[name] = vcat(result1[name], result2[name])
  end
  return result
end

minute = 60
hour = 60*minute

@model Battery begin
  # Very simple model: Battery modelled as a capacitor.
  @extends Capacitor()
  @inherits C, v
  charge = Var()
@equations begin
  charge = C*v
  end
end

@model Charger begin
  # Very simple model: a constant voltage behind a resistor.
  v0=1000
  r=100
  V = ConstantVoltage(V=v0)
  R = Resistor(R=r)
  p = Pin()
  n = Pin()
  v = Var()
@equations begin
  connect(V.p, R.n)
  connect(V.n, n)
  connect(R.p, p)
  v = p.v - n.v
  end
end
result = simulate(Charger, 1*hour)

@model ElectricVehicle begin
  # Very simple model: Battery modelled as a capacitor and the load as a resistor.
  B = Battery(C=1000)
  load = Resistor(R=1)
  p = Pin()
  n = Pin()
@equations begin
  connect(B.p, load.p)
  connect(B.n, load.n)
  connect(B.n, n)
  connect(B.p, p)
  end
end
result = simulate(ElectricVehicle, 1*hour)
Bv = result["B.v"][end]  # Available voltage after driving used to intialize next model

Bv = 10.0

@model ElectricalVehicleWithCharger begin
  ev = ElectricVehicle(B = Battery(C=1000, v = Float(start=Bv)))
  c = if rand(1:2)==1; (println("Super Charger"); Charger()) else (println("Standard Charger"); Charger(v0=500, r=200)) end
 # y = Var(size=())
end

global now = 0
then = 1*hour

global result = simulate(ElectricalVehicleWithCharger, then)
now = then
Bv = result["ev.B.v"][end]

for i in 1:5
  global now
  global result
  println("\nConnecting Electric Vehicle to Charger")
  #addEquation!(ElectricalVehicleWithCharger, :(connect(ev.p, c.p)))
  #addEquation!(ElectricalVehicleWithCharger, :(connect(ev.n, c.n)))
  addEquation!(ElectricalVehicleWithCharger, @equation(connect(this.ev.p, this.c.p)))
  addEquation!(ElectricalVehicleWithCharger, @equation(connect(this.ev.n, this.c.n)))
  #addEquation!(ElectricalVehicleWithCharger, @equation(this.y = this.ev.C.v))

  then = now+15*minute
  resultCharging = simulate(ElectricalVehicleWithCharger, then, startTime=now, logStatistics=false)
  now = then

  result = concatenateTimeSeries(result, resultCharging)
  plot(result, ("ev.B.charge"), figure=17)
  Bv = result["ev.B.v"][end]

  println("\nDisconnect")
  deleteEquation!(ElectricalVehicleWithCharger, @equation(connect(this.ev.p, this.c.p)))
  deleteEquation!(ElectricalVehicleWithCharger, @equation(connect(this.ev.n, this.c.n)))

  then = now+1*hour
  resultDriving = simulate(ElectricalVehicleWithCharger, then, startTime=now, logStatistics=false)
  now = then

  result = concatenateTimeSeries(result, resultDriving)
  plot(result, ("ev.B.charge"), figure=17)
  Bv = result["ev.B.v"][end]
end

end


text:
using Documenter, AppliAR

makedocs(;
    modules=[AppliAR],
    format=Documenter.HTML(),
    pages=[
        "Accounts Receivable" => "index.md",
        "1 - Domain" => "chapter1.md",
        "2 - API" => "chapter2.md",
        "3 - Infrastructure" => "chapter3.md",
        "4 - Example" => "chapter4.md",
    ],
    repo="https://github.com/rbontekoe/AppliAR.jl/blob/{commit}{path}#L{line}",
    sitename="AppliAR.jl",
    authors="Rob Bontekoe",
    assets=String[],
)

deploydocs(;
    repo="github.com/rbontekoe/AppliAR.jl",
)


text:
using JSON3
using DFTK
using JLD2
using MPI
using Unitful
using UnitfulAtomic

function setup()
    if DFTK.mpi_nprocs() > 1
        DFTK.disable_threading()
    else
        DFTK.setup_threading()
    end
end

inputerror(parameters, p) = error("Unknown value to parameter $p: $(get(parameters, p, ""))")

function get_dftk_model(parameters, extra)
    # Parse psp and DFT functional parameters
    functionals = [:lda_xc_teter93]
    if lowercase(parameters[:xc]) == "lda"
        psp_functional = "lda"
        functionals = [:lda_xc_teter93]
    elseif lowercase(parameters[:xc]) == "pbe"
        psp_functional = "pbe"
        functionals = [:gga_x_pbe, :gga_c_pbe]
    else
        inputerror(parameters, "xc")
    end

    if !isnothing(parameters[:functionals])
        funs = parameters[:functionals]
        if funs isa AbstractArray
            functionals = Symbol.(funs)
        else
            functionals = [Symbol(funs)]
        end
    end

    if parameters[:pps] == "hgh"
        psp_family = "hgh"
        psp_core = :fullcore
    elseif parameters[:pps] == "hgh.k"
        psp_family = "hgh"
        psp_core = :semicore
    else
        inputerror(parameters, "pps")
    end

    # Parse smearing and temperature
    temperature = 0.0u"eV"
    smearing = nothing
    if !isnothing(parameters[:smearing])
        psmear = parameters[:smearing]
        if lowercase(psmear[1]) == "fermi-dirac"
            smearing = Smearing.FermiDirac()
        elseif lowercase(psmear[1]) == "gaussian"
            smearing = Smearing.Gaussian()
        elseif lowercase(psmear[1]) == "methfessel-paxton"
            smearing = Smearing.MethfesselPaxton(psmear[3])
        else
            inputerror(parameters, "smearing")
        end
        temperature = convert(Float64, psmear[2]) * u"eV"
    end

    if parameters[:charge] != 0.0
        error("Charged systems not yet supported in DFTK.")
    end

    # TODO This is the place where spin-polarization should be added
    #      once it is in DFTK.

    # Build DFTK atoms
    psploader(symbol) = load_psp(symbol, functional=psp_functional,
                                 family=psp_family, core=psp_core)
    atoms = [ElementPsp(element.symbol, psp=psploader(element.symbol)) => positions
             for (element, positions) in extra[:atoms]]

    model_DFT(extra[:lattice], atoms, functionals; temperature=temperature,
              smearing=smearing)
end

function get_dftk_basis(parameters, extra; model=get_dftk_model(parameters, extra))
    # Build kpoint mesh
    kpts = parameters[:kpts]
    kgrid = [1, 1, 1]
    kshift = [0, 0, 0]
    if kpts isa Number
        # DFTK uses k-point spacing whereas ASE uses k-point density
        kgrid = kgrid_size_from_minimal_spacing(model.lattice, 1 / kpts / u"Å")
    elseif length(kpts) == 3 && all(kpt isa Number for kpt in kpts)
        kgrid = kpts  # Just a plain MP grid
    elseif length(kpts) == 4 && all(kpt isa Number for kpt in kpts[1:3])
        kpts[4] != "gamma" && error("Unknown value to kpts: $kpts")
        kgrid = kpts[1:3]
        kshift = Int.(iseven.(kgrid)) .// 2  # Shift MP grid to always contain Gamma
    elseif kpts isa AbstractArray
        kgrid = nothing
        kshift = nothing
        kcoords = [Vec3(kpt...) for kpt in kpts]
        ksymops = [[(Mat3{Int}(I), Vec3(zeros(3)))] for _ in 1:length(kcoords)]
    end

    # Convert ecut to Hartree
    Ecut = austrip(parameters[:ecut] * u"eV")
    if isnothing(kgrid)
        PlaneWaveBasis(model, Ecut, kcoords, ksymops)
    else
        PlaneWaveBasis(model, Ecut, kgrid=kgrid, kshift=kshift)
    end
end

function get_dftk_mixing(parameters, extra; basis=get_dftk_basis(parameters, extra))
    if isnothing(parameters[:mixing])
        basis.model.temperature > 0 ? KerkerMixing() : SimpleMixing(α=0.8)
    else
        include_string(Main, parameters[:mixing], @__FILE__)
    end
end

function get_dftk_scfres(parameters, extra; basis=get_dftk_basis(parameters, extra))
    mixing = get_dftk_mixing(parameters, extra; basis=basis)

    n_bands = DFTK.default_n_bands(basis.model)
    if !isnothing(parameters[:nbands])
        n_bands = parameters[:nbands]
    end

    if mpi_master()
        println()
        println("temperature     = $(basis.model.temperature)")
        println("smearing        = $(basis.model.smearing)")
        println("Ecut            = $(basis.Ecut)")
        println("fft_size        = $(basis.fft_size)")
        println("kgrid           = $(basis.kgrid)")
        println("kshift          = $(basis.kshift)")
        println("n_irreducible_k = $(sum(length, basis.krange_allprocs))")
        println("n_electr        = $(basis.model.n_electrons)")
        println()
        println("mixing          = $mixing")
        println("n_bands         = $n_bands")
        println()
        flush(stdout)
    end

    callback = DFTK.ScfSaveCheckpoints(extra[:checkpointfile]) ∘ DFTK.ScfDefaultCallback()
    self_consistent_field(basis; tol=parameters[:scftol], callback=callback,
                          mixing=mixing, n_bands=n_bands)
end

function parse_json_array(data::AbstractDict)
    # Simplistic parser for ndarrays as they are stored in ASE json files
    shape, type, values = data["__ndarray__"]
    shape = reverse(Int.(shape))  # row-major -> column-major
    if type == "float64"
        reshape(Float64.(values), shape...)
    elseif type == "int64"
        reshape(Int64.(values), shape...)
    else
        error("Type not implemented: $type")
    end
end

function load_state(file)
    endswith(file, ".json") || error("State file should end in .json")

    # Read input data
    if DFTK.mpi_master()
        str = open(fp -> read(fp, String), file)
    else
        str = nothing
    end
    str = MPI.bcast(str, 0, MPI.COMM_WORLD)
    res = copy(JSON3.read(str))

    # Parse atoms json
    atoms_json = JSON3.read(res[:atoms])
    if length(atoms_json["ids"]) != 1
        @warn "Only parsing last atoms object in json. Ignoring all others"
    end
    atoms_id = string(atoms_json["ids"][end])

    cell        = atoms_json[atoms_id]["cell"]
    lattice_Ang = parse_json_array(get(cell, "array", cell))

    positions_Ang = parse_json_array(atoms_json[atoms_id]["positions"])
    numbers       = parse_json_array(atoms_json[atoms_id]["numbers"])

    atoms = []
    positions_fractional = Array.(eachcol(lattice_Ang \ positions_Ang))
    for (i, pos) in enumerate(positions_fractional)
        id = findfirst(el_pos -> el_pos[1].Z == numbers[i], atoms)
        if isnothing(id)
            push!(atoms, ElementCoulomb(numbers[i]) => [pos])
        else
            push!(atoms[id][2], pos)
        end
    end

    res[:extra] = Dict{Symbol, Any}()
    res[:extra][:lattice]    = lattice_Ang * austrip(1u"Å")
    res[:extra][:atoms]      = atoms
    res
end

function save_state(file, state)
    endswith(file, ".json") || error("State file should end in .json")

    # Transpose arrays inside results, such that they are written
    # in a predictable order to disk
    save_results = Dict{Symbol, Any}()
    for (key, val) in pairs(state[:results])
        if val isa AbstractArray
            if ndims(val) == 2
                save_results[key] = Vector.(eachrow(val))
                continue
            elseif ndims(val) > 2
                error("Arrays of dimension larger 2 not supported")
            end
        end
        save_results[key] = val
    end

    save_dict = Dict(
        :parameters => state[:parameters],
        :results    => save_results,
        :atoms      => state[:atoms],
        :scfres     => state[:scfres],
    )
    if DFTK.mpi_master()
        open(file, "w") do fp
            JSON3.write(fp, save_dict)
        end
    end
    save_dict
end

function run_calculation(properties::AbstractArray, statefile::AbstractString)
    state = load_state(statefile)
    @assert endswith(state[:scfres], ".scfres.jld2")
    prefix = state[:scfres][1:end-12]
    state[:extra][:checkpointfile] = prefix * ".checkpoint.jld2"

    if !isfile(state[:scfres])
        save_scfres(state[:scfres], get_dftk_scfres(state[:parameters], state[:extra]))
    end
    scfres = load_scfres(state[:scfres])

    state[:results][:energy] = ustrip(auconvert(u"eV", scfres.energies.total))

    # For now always compute forces, since it's pretty useful (and cheap)
    begin  # if "forces" in properties
        forces = compute_forces_cart(scfres)

        # DFTK has forces as Hartree over fractional coordinates
        # ASE wants forces as eV / Å
        n_atoms = sum(length, forces)
        cart_forces = zeros(eltype(scfres.basis), n_atoms, 3)
        for (i, atforce) in enumerate(Iterators.flatten(forces))
            cart_forces[i, :] = @. ustrip(auconvert(u"eV/Å", atforce))
        end
        state[:results][:forces] = cart_forces
    end
    save_state(statefile, state)
end

function main()
    setup()
    run_calculation(ARGS[1:end-1], ARGS[end])
end

(abspath(PROGRAM_FILE) == @__FILE__) && main()

max_stars_repo_name: mfherbst/ase-dftk

max_issues_repo_name: mfherbst/ase-dftk

max_forks_repo_name: mfherbst/asedftk


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -23.072345451695163 <= q <= 18.38163265505858)

                     add_NL_constraint(m, :(gelu(-0.6895881478464707 + 0.42824131807966515*gelu(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4])) + -0.8758217593243995*gelu(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4])) + -0.3819126296826152*gelu(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4])) + -0.7523924735752416*gelu(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))) + gelu(0.5903696565764314 + 0.6045010784464115*gelu(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4])) + 0.19992590084675577*gelu(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4])) + -0.13493555263327828*gelu(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4])) + 0.0007730145672497635*gelu(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))) + gelu(-0.987583954230804 + 0.6297203029880079*gelu(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4])) + -0.939181196020844*gelu(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4])) + 0.23747324715066886*gelu(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4])) + 0.2556344288953589*gelu(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))) + gelu(-0.44105426128556546 + -0.39452594471110514*gelu(-0.9657925733535908 + -0.20993433529810623*$(x[1]) + -0.9858887282433288*$(x[2]) + -0.7926811137744871*$(x[3]) + -0.8802726339564089*$(x[4])) + 0.16620534874676274*gelu(-0.29556778373435044 + 0.6695148008497993*$(x[1]) + 0.7620585252555783*$(x[2]) + 0.8362720647957991*$(x[3]) + 0.8751574333999192*$(x[4])) + 0.6245292013722641*gelu(-0.8627823044260601 + -0.3188988197195952*$(x[1]) + -0.18058802482760727*$(x[2]) + -0.7255181409249616*$(x[3]) + -0.8658172698743405*$(x[4])) + 0.9983459126930239*gelu(0.5506348394687617 + 0.9730675890304568*$(x[1]) + 0.0064031244486910666*$(x[2]) + -0.590531971380631*$(x[3]) + -0.964793293764588*$(x[4]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
if Compat.is_unix()
    "R eventloop hook on Unix system"
    function polled_events()
        event_callback()
        nothing
    end
end

"Event Callback: allows R to process Julia events when R is busy.
For example, writing output to STDOUT while running an expensive R command."
function event_callback()
    # dump printBuffer to STDOUT when available
    flush_print_buffer(STDOUT)
    nothing
end

# there is no use now, maybe useful for the future.
function interrupts_pending(s::Bool=true)
    if Compat.is_windows()
        unsafe_store!(cglobal((:UserBreak,libR),Cint), s?1:0)
    else
        unsafe_store!(cglobal((:R_interrupts_pending,libR),Cint), s?1:0)
    end
    nothing
end

# this shouldn't exist if we could hook into Julia eventloop.
function process_events()
    ##FIXME: a dirty fix to prevent segfault right after a sigint
    if unsafe_load(cglobal((:R_interrupts_pending,libR),Cint)) == 0
        if Compat.is_windows() || Compat.is_apple()
            ccall((:R_ProcessEvents, libR), Void, ())
        end
        if Compat.is_unix()
            what = ccall((:R_checkActivity,libR),Ptr{Void},(Cint,Cint),0,1)
            if what != C_NULL
                R_InputHandlers = unsafe_load(cglobal((:R_InputHandlers,libR),Ptr{Void}))
                ccall((:R_runHandlers,libR),Void,(Ptr{Void},Ptr{Void}),R_InputHandlers,what)
            end
        end
    end
    nothing
end

global timeout = nothing

function rgui_start(silent=false)
    global timeout
    if timeout == nothing
        timeout = Base.Timer(x -> process_events(), 0.05, 0.05)
        return true
    else
        silent || error("eventloop is already running.")
        return false
    end
end

function rgui_stop(silent=false)
    global timeout
    if timeout != nothing
        close(timeout)
        timeout = nothing
        return true
    else
        silent || error("eventloop is not running.")
        return false
    end
end

function set_hook(hookname, value)
    l = rparse("""setHook(hookname, function(...) foo)""")
    l[1][2] = hookname
    l[1][3][3] = value
    reval(l)
end

function rgui_init()
    f = rlang(rgui_start, true)
    set_hook("plot.new", f)
    set_hook("grid.newpage", f)
    set_hook(rlang(:packageEvent, "rgl", "onLoad"), f)
end


text:
module AlignedSpans

using Dates, Intervals, Onda
using TimeSpans: TimeSpans, start, stop, format_duration
using StructTypes, ArrowTypes

export SpanRoundingMode, RoundInward, RoundSpanDown, ConstantSamplesRoundingMode
export AlignedSpan, consecutive_subspans, n_samples

# Make our own method so we can add methods for Intervals without piracy
duration(span) = TimeSpans.duration(span)

#####
##### Types and rounding modes
#####

"""
    SpanRoundingMode(start::RoundingMode, stop::RoundingMode)

Creates a rounding object for [`AlignedSpan`](@ref) to indicate how the `AlignedSpan`'s
endpoints should be determined from a given `span`s endpoints'.
"""
struct SpanRoundingMode
    start::RoundingMode
    stop::RoundingMode
end

"""
    ConstantSamplesRoundingMode(start::RoundingMode)

Creates a rounding object for [`AlignedSpan`](@ref) to indicate the `AlignedSpan`
should be constructed by the `start` and `duration` of the `span`, without regard to its `stop`.

If two `span`s have the same duration, then the resulting `AlignedSpan`'s will have the same
number of samples when constructed with this rounding mode.

See also [`AlignedSpan(sample_rate, span, mode::ConstantSamplesRoundingMode)`](@ref).
"""
struct ConstantSamplesRoundingMode
    start::RoundingMode
end

"""
    RoundInward = SpanRoundingMode(RoundUp, RoundDown)

This is a rounding mode where both ends of the continuous time interval are rounded "inwards"
to construct the largest span of indices such that all samples are entirely contained within it.

## Example

Consider a signal with sample rate 1 Hz.

```
Index       1   2   3   4   5
Time (s)    0   1   2   3   4
```

Now, consider the time span 1.5s (inclusive) to 2.5s (exclusive).
Using brackets to highlight this span:

```
Index       1   2     3     4   5
Time (s)    0   1  [  2  )  3   4
```

In code, this span is described by

```jldoctest RoundInward
julia> using AlignedSpans, Dates, TimeSpans

julia> ts = TimeSpan(Millisecond(1500), Millisecond(2500))
TimeSpan(00:00:01.500000000, 00:00:02.500000000)
```

The only sample within the span is at index 3. And indeed,
```jldoctest RoundInward
julia> aligned = AlignedSpan(1, ts, RoundInward)
AlignedSpan(1.0, 3, 3)

julia> AlignedSpans.indices(aligned)
3:3
```
gives an `AlignedSpan` with indices `3:3`.
"""
const RoundInward = SpanRoundingMode(RoundUp, RoundDown)

"""
    RoundSpanDown = SpanRoundingMode(RoundDown, RoundDown)

This is a rounding mode where *both* ends of the continuous time interval are rounded
downwards.

## Example

Consider a signal with sample rate 1 Hz.

```
Index       1   2   3   4   5
Time (s)    0   1   2   3   4
```

Now, consider the time span 1.5s (inclusive) to 2.5s (exclusive).
Using brackets to highlight this span:

```
Index       1   2     3     4   5
Time (s)    0   1  [  2  )  3   4
```

In code, this span is described by

```jldoctest RoundSpanDown
julia> using AlignedSpans, Dates, TimeSpans

julia> ts = TimeSpan(Millisecond(1500), Millisecond(2500))
TimeSpan(00:00:01.500000000, 00:00:02.500000000)
```

If we round both ends of the interval down to the nearest sample,
the start of the interval becomes 1s, and the stop of the interval
becomes 2s. Thus, the associated samples are at indices `2:3`. And indeed,
```jldoctest RoundSpanDown
julia> aligned = AlignedSpan(1, ts, RoundSpanDown)
AlignedSpan(1.0, 2, 3)

julia> AlignedSpans.indices(aligned)
2:3
```
gives an `AlignedSpan` with indices `2:3`.
"""
const RoundSpanDown = SpanRoundingMode(RoundDown, RoundDown)

"""
    AlignedSpan(sample_rate::Number, first_index::Int, last_index::Int)

Construct an `AlignedSpan` directly from a `sample_rate` and indices.
"""
struct AlignedSpan
    sample_rate::Float64
    first_index::Int64
    last_index::Int64
    function AlignedSpan(sample_rate::Number, first_index::Int, last_index::Int)
        if last_index < first_index
            throw(ArgumentError("Cannot create `AlignedSpan` with right-endpoint (`last_index=$(last_index)`) strictly smaller than left endpoint (`first_index=$(first_index)`)"))
        end
        return new(convert(Float64, sample_rate), first_index, last_index)
    end
end

#####
##### Continuous -> discrete interface
#####

# Methods for these API functions are provided in `interop.jl`.

"""
    AlignedSpans.start_index_from_time(sample_rate, span, rounding_mode)

Returns the index of a sample object obtained by rounding the start of `span` according to `rounding_mode`.

See also [`AlignedSpan(sample_rate, span, mode::SpanRoundingMode)`](@ref) and
[`AlignedSpan(sample_rate, span, mode::ConstantSamplesRoundingMode)`](@ref).
"""
function start_index_from_time end

"""
    AlignedSpans.stop_index_from_time(sample_rate, span, rounding_mode)

Returns the index of a sample object obtained by rounding the stop of `span` according to `rounding_mode`.

See also [`AlignedSpan(sample_rate, span, mode::SpanRoundingMode)`](@ref).
"""
function stop_index_from_time end

#####
##### Continuous -> discrete conversions
#####

"""
    AlignedSpan(sample_rate, span, mode::SpanRoundingMode)

Creates an `AlignedSpan` by rounding the left endpoint according to `mode.start`,
and the right endpoint by `mode.stop`.

If `mode.start==RoundUp`, then the left index of the resulting span is guaranteed
to be inside `span`. This is accomplished by checking if the left endpoint of the span
is exclusive, and if so, incrementing the index after rounding when necessary.

Likewise, if `mode.start==RoundDown`, then the right index of the resulting span is guaranteed
to be inside `span`. This is accomplished by checking if the right endpoint of the span
is exclusive, and if so, decrementing the index after rounding when necessary.

Note: `span` may be of any type which which provides methods for `AlignedSpans.start_index_from_time` and `AlignedSpans.stop_index_from_time`.
"""
function AlignedSpan(sample_rate, span, mode::SpanRoundingMode)
    first_index = start_index_from_time(sample_rate, span, mode.start)
    last_index = stop_index_from_time(sample_rate, span, mode.stop)
    if last_index < first_index
        throw(ArgumentError("No samples lie within `span`"))
    end
    return AlignedSpan(sample_rate, first_index, last_index)
end

"""
    AlignedSpan(sample_rate, span, mode::ConstantSamplesRoundingMode)

Creates an `AlignedSpan` whose left endpoint is rounded according to `mode.start`,
and whose right endpoint is determined so by the left endpoint and the number of samples,
given by `AlignedSpans.n_samples(sample_rate, duration(span))`.

Interface: `span` may be of any type which which provides a method for [`AlignedSpans.start_index_from_time`](@ref) and `TimeSpans.duration`.

## More detailed information

This is designed so that if `AlignedSpan(sample_rate, span, mode::ConstantSamplesRoundingMode)` is applied to multiple `span`s, with the same `sample_rate`, and the same durations, then the resulting `AlignedSpan`'s will have the same number of samples.

For this reason, we ask for `TimeSpans.duration(span)` to be defined, rather than a `n_samples(span)` function: the idea is that we want to only using the duration and the starting time, rather than the *actual* number of samples in this particular `span`.

In contrast, `AlignedSpan(sample_rate, span, RoundInward)` provides an `AlignedSpan` which includes only (and exactly) the samples contained within `span`.

If one wants to create a collection of consecutive, non-overlapping, `AlignedSpans` each with the same number of samples, then use [`consecutive_subspans`](@ref) instead.
"""
function AlignedSpan(sample_rate, span, mode::ConstantSamplesRoundingMode)
    first_index = start_index_from_time(sample_rate, span, mode.start)
    n = n_samples(sample_rate, duration(span))
    last_index = first_index + n - 1
    return AlignedSpan(sample_rate, first_index, last_index)
end

include("time_index_conversions.jl")
include("interop.jl")
include("utilities.jl")

end


text:
#---------------------------------------------------------
# # [MRI SENSE](@id 10-mri-sense)
#---------------------------------------------------------

#=
This page illustrates the `mri_smap_fit` and `mri_spectra` methods
in the Julia package
[`ImagePhantoms`](https://github.com/JuliaImageRecon/ImagePhantoms.jl)
for performing MRI simulations with realistic sensitivity encoding (SENSE).

This page was generated from a single Julia file:
[10-mri-sense.jl](@__REPO_ROOT_URL__/10-mri-sense.jl).
=#

#md # In any such Julia documentation,
#md # you can access the source code
#md # using the "Edit on GitHub" link in the top right.

#md # The corresponding notebook can be viewed in
#md # [nbviewer](http://nbviewer.jupyter.org/) here:
#md # [`10-mri-sense.ipynb`](@__NBVIEWER_ROOT_URL__/10-mri-sense.ipynb),
#md # and opened in [binder](https://mybinder.org/) here:
#md # [`10-mri-sense.ipynb`](@__BINDER_ROOT_URL__/10-mri-sense.ipynb).

# ### Setup

# Packages needed here.

using ImagePhantoms: ellipse_parameters, SheppLoganBrainWeb, Ellipse
using ImagePhantoms: phantom, mri_smap_fit, mri_spectra
using FFTW: fft, fftshift
using ImageGeoms: embed
using LazyGrids: ndgrid
using MIRT: ir_mri_sensemap_sim
using MIRTjim: jim, prompt; jim(:prompt, true)
using Random: seed!
using Unitful: mm

# The following line is helpful when running this example.jl file as a script;
# this way it will prompt user to hit a key after each figure is displayed.

isinteractive() ? jim(:prompt, true) : prompt(:draw);

# ### Overview

#=
Modern MRI scanners use multiple receive coils
each of which has its own "sensitivity map" (or "coil profile").
Realistic MRI simulations should account for the effects of those
sensitivity maps analytically, rather than committing the "inverse crime"
of using rasterized phantoms and maps.

See the 2012 paper
[Guerquin-Kern et al.](http://doi.org/10.1109/TMI.2011.2174158)
that combines analytical k-space values of the phantom
with an analytical model for the sensitivity maps.
This package follows the recommended approach from that paper.
We used the `mri_smap_fit` function to fit each sensitivity map
with a modest number of complex exponential signals.
Then, instead of using the `spectrum` function
we use the `spectra` function
to generate simulated k-space data
from analytical phantoms (like ellipses).
=#

# Because FFTW.fft cannot handle units, this function is a work-around.
function myfft(x::AbstractArray{T}) where {T <: Number}
    u = oneunit(T)
    return fftshift(fft(fftshift(x) / u)) * u
end

# ### Phantom

# Image geometry:

fovs = (256mm, 250mm)
nx, ny = (128, 100) .* 2
dx, dy = fovs ./ (nx,ny)
x = (-(nx÷2):(nx÷2-1)) * dx
y = (-(ny÷2):(ny÷2-1)) * dy

# Define Shepp-Logan phantom object, with random complex phases
# to make it a bit more realistic.

oa = ellipse_parameters(SheppLoganBrainWeb() ; disjoint=true, fovs)
seed!(0)
oa[:,end] = [1; randn(ComplexF32, 9)] # random phases
oa = Ellipse(oa)
oversample = 3
image0 = phantom(x, y, oa, oversample)
cfun = z -> cat(dims = ndims(z)+1, real(z), imag(z))
jim(x, y, cfun(image0), "Digital phantom\n (real | imag)"; ncol=1)

#=
In practice, sensitivity maps are usually estimated
only over portion of the image array,
so we define a simple `mask` here
to exercise this issue.
=#

mask = trues(nx,ny)
mask[:,[1:2;end-2:end]] .= false
mask[[1:8;end-8:end],:] .= false
@assert mask .* image0 == image0
jim(x, y, mask, "mask")

# ### Sensitivity maps

# Here we use simulated sensitivity maps from MIRT:

ncoil = 2
smap = ir_mri_sensemap_sim(dims=(nx,ny), ncoil=ncoil, orbit_start=[0])
jim(x, y, cfun(smap), "Sensitivity maps raw")

#=
Typical sensitivity map estimation methods
normalize the maps
so that the square-root of the sum of squares (SSoS) is unity:
=#

ssos = sqrt.(sum(abs.(smap).^2, dims=ndims(smap))) # SSoS
ssos = selectdim(ssos, ndims(smap), 1)
jim(ssos, "SSoS for ncoil=$ncoil")

for ic=1:ncoil # normalize
    selectdim(smap, ndims(smap), ic) ./= ssos
end
smap .*= mask
stacker = x -> [(@view x[:,:,i]) for i=1:size(x,3)]
smaps = stacker(smap) # code hereafter expects vector of maps
jim(x, y, cfun(smaps), "Sensitivity maps (masked and normalized)")

# ### Sensitivity map fitting using complex exponentials

#=
The `mri_smap_fit` function fits each `smap`
with a linear combination of complex exponential signals.
(These signals are not orthogonal due to the `mask`.)
With frequencies `-7:7/N`, the maximum error is ≤ 0.2%.
=#

deltas = (dx, dy)
kmax = 7
fit = mri_smap_fit(smaps, embed; mask, kmax, deltas)
jim(
 jim(x, y, cfun(smaps), "Original maps"; prompt=false, clim=(-1,1)),
 jim(x, y, cfun(fit.smaps), "Fit maps"; prompt=false, clim=(-1,1)),
 jim(x, y, cfun(100 * (fit.smaps - smaps)), "error * 100"; prompt=false),
)

# The fit coefficients are smaller near `kmax`
# so probably `kmax` is large enough.

coefs = map(x -> reshape(x,15,15), fit.coefs)
jim(-kmax:kmax, -kmax:kmax, cfun(coefs), "coefs", prompt=false)

# ### Compare FFT with analytical spectra

# Frequency sample vectors:
fx = (-(nx÷2):(nx÷2-1)) / (nx*dx) # crucial to match `mri_smap_basis` internals!
fy = (-(ny÷2):(ny÷2-1)) / (ny*dy)
gx, gy = ndgrid(fx, fy);

# Analytical spectra computation for complex phantom using all smaps.
# Note the `fit` argument.
kspace1 = mri_spectra(vec(gx), vec(gy), oa, fit)
kspace1 = [reshape(k, nx, ny) for k in kspace1]
p1 = jim(fx, fy, cfun(kspace1), "Analytical")

# FFT spectra computation based on digital image and sensitivity maps:
image2 = [image0 .* s for s in smaps] # digital
kspace2 = myfft.(image2) * (dx * dy)
p2 = jim(fx, fy, cfun(kspace2), "FFT-based")

#
p3 = jim(fx, fy, cfun(kspace2 - kspace1), "Error")

# Zoom in to illustrate similarity:

xlims = (-1,1) .* (0.06/mm)
ylims = (-1,1) .* (0.06/mm)
jim(
 jim(fx, fy, real(kspace1[1]), "Analytical"; xlims, ylims),
 jim(fx, fy, real(kspace2[1]), "FFT-based"; xlims, ylims),
 jim(fx, fy, real(kspace2[1] - kspace1[1]), "Error"; xlims, ylims),
)

#=
In summary, the `mri_smap_fit` and `mri_spectra` methods here
reproduce the approach in the 2012 Guerquin-Kern paper, cited above,
enabling parallel MRI simulations that avoid an inverse crime.
=#