text:
"""
    arec(A, G, Q = 0; as = false, rtol::Real = nϵ) -> (X, EVALS, Z)

Compute `X`, the hermitian/symmetric stabilizing solution (if `as = false`) or
anti-stabilizing solution (if `as = true`) of the continuous-time
algebraic Riccati equation

     A'X + XA - XGX + Q = 0,

where `G` and `Q` are hermitian/symmetric matrices or uniform scaling operators.
Scalar-valued `G` and `Q` are interpreted as appropriately sized uniform scaling operators `G*I` and `Q*I`.

By default, the lower bound for the 1-norm reciprocal condition number `rtol` is `n*ϵ`, where `n` is the order of `A`
and `ϵ` is the _machine epsilon_ of the element type of `A`.

`EVALS` is a vector containing the (stable or anti-stable) eigenvalues of `A-GX`.
`Z = [ U; V ]` is an orthogonal basis for the stable/anti-stable deflating subspace such that `X = V/U`.

`Reference:`
Laub, A.J., A Schur Method for Solving Algebraic Riccati equations.
IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.

# Example
```jldoctest
julia> using LinearAlgebra

julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
 -6.0  -2.0   1.0
  5.0   1.0  -1.0
 -4.0  -2.0  -1.0

julia> G = [1. 0. 0.; 0. 5. 0.; 0. 0. 10.]
3×3 Array{Float64,2}:
 1.0  0.0   0.0
 0.0  5.0   0.0
 0.0  0.0  10.0

julia> X, CLSEIG = arec(A,G,2I);

julia> X
3×3 Array{Float64,2}:
  0.459589   0.333603   -0.144406
  0.333603   0.65916    -0.0999216
 -0.144406  -0.0999216   0.340483

julia> A'*X+X*A-X*G*X+2I
3×3 Array{Float64,2}:
  2.22045e-16  4.44089e-16  -1.77636e-15
  4.44089e-16  6.66134e-16   1.11022e-16
 -1.77636e-15  1.11022e-16  -1.33227e-15

julia> CLSEIG
3-element Array{Complex{Float64},1}:
 -4.411547592296008 + 2.4222082620381102im
 -4.411547592296008 - 2.4222082620381102im
 -4.337128244724371 + 0.0im

julia> eigvals(A-G*X)
3-element Array{Complex{Float64},1}:
 -4.4115475922960075 - 2.4222082620381076im
 -4.4115475922960075 + 2.4222082620381076im
  -4.337128244724374 + 0.0im
```
"""
function arec(A::AbstractMatrix, G::Union{AbstractMatrix,UniformScaling,Real,Complex}, Q::Union{AbstractMatrix,UniformScaling,Real,Complex} = zero(eltype(A));
              as = false, rtol::Real = size(A,1)*eps(real(float(one(eltype(A))))))
    n = LinearAlgebra.checksquare(A)
    T = promote_type( eltype(A), eltype(G), eltype(Q) )
    if typeof(G) <: AbstractArray
       (LinearAlgebra.checksquare(G) == n && ishermitian(G)) ||
          throw(DimensionMismatch("G must be a symmetric/hermitian matrix of dimension $n"))
    else
       G = G*I
       iszero(imag(G.λ)) || throw("G must be a symmetric/hermitian matrix")
    end
    if typeof(Q) <: AbstractArray
       (LinearAlgebra.checksquare(Q) == n && ishermitian(Q)) ||
          throw(DimensionMismatch("Q must be a symmetric/hermitian matrix of dimension $n"))
    else
       Q = Q*I
       iszero(imag(Q.λ)) || throw("Q must be a symmetric/hermitian matrix")
    end
    T <: BlasFloat || (T = promote_type(Float64,T))
    TR = real(T)
    epsm = eps(TR)
    eltype(A) == T || (A = convert(Matrix{T},A))
    eltype(G) == T || (typeof(G) <: AbstractMatrix ? G = convert(Matrix{T},G) : G = convert(T,G.λ)*I)
    eltype(Q) == T || (typeof(Q) <: AbstractMatrix ? Q = convert(Matrix{T},Q) : Q = convert(T,Q.λ)*I)

    n == 0 && (return  zeros(T,0,0), zeros(T,0), zeros(T,m,0) )

    S = schur([A  -G; -Q  -copy(A')])

    as ? select = real(S.values) .> 0 : select = real(S.values) .< 0
    n == length(filter(y-> y == true,select)) || error("The Hamiltonian matrix is not dichotomic")
    ordschur!(S, select)

    n2 = n+n
    ix = 1:n
    F = _LUwithRicTest(S.Z[ix, ix],rtol)
    x = S.Z[n+1:n2, ix]/F
    return  (x+x')/2, S.values[ix], S.Z[:,ix]
end
function _LUwithRicTest(Z11::AbstractArray,rtol::Real)
   try
      F = LinearAlgebra.lu(Z11)
      Z11norm = opnorm(Z11,1)
      Z11norm > 2*rtol ? (rcond = LAPACK.gecon!('1',F.factors,Z11norm)) : (rcond = zero(eltype(Z11)))
      rcond <= rtol ? error("no finite solution exists for the Riccati equation") : (return  F)
    catch
      error("no finite solution exists for the Riccati equation")
   end
end
"""
    arec(A, B, R, Q, S; as = false, rtol::Real = nϵ, orth = false) -> (X, EVALS, F, Z)

Computes `X`, the hermitian/symmetric stabilizing solution (if `as = false`) or
anti-stabilizing solution (if `as = true`) of the continuous-time
algebraic Riccati equation

     A'X + XA - (XB+S)R^(-1)(B'X+S') + Q = 0,

where `R` and `Q` are hermitian/symmetric matrices or uniform scaling operators such that `R` is nonsingular.
Scalar-valued `R` and `Q` are interpreted as appropriately sized uniform scaling operators `R*I` and `Q*I`.
`S`, if not specified, is set to `S = zeros(size(B))`.

By default, the lower bound for the 1-norm reciprocal condition number `rtol` is `n*ϵ`, where `n` is the order of `A`
and `ϵ` is the _machine epsilon_ of the element type of `A`.

`EVALS` is a vector containing the (stable or anti-stable) eigenvalues of `A-BF`.
`F` is the stabilizing or anti-stabilizing gain matrix `F = R^(-1)(B'X+S')`.
`Z = [ U; V; W ]` is a basis for the relevant stable/anti-stable deflating subspace such that `X = V/U` and `F = -W/U`.
An orthogonal basis Z can be determined, with an increased computational cost, by setting `orth = true`.

`Reference:`
Laub, A.J., A Schur Method for Solving Algebraic Riccati equations.
IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.

# Example
```jldoctest
julia> using LinearAlgebra

julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
 -6.0  -2.0   1.0
  5.0   1.0  -1.0
 -4.0  -2.0  -1.0

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

julia> R = [1. 0.; 0. 5.]
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  5.0

julia> X, CLSEIG, F = arec(A,B,R,2I);

julia> X
3×3 Array{Float64,2}:
  0.522588   0.303007  -0.327227
  0.303007   0.650895  -0.132608
 -0.327227  -0.132608   0.629825

julia> A'*X+X*A-X*B*inv(R)*B'*X+2I
3×3 Array{Float64,2}:
 -2.66454e-15  -1.55431e-15   8.88178e-16
 -1.55431e-15   2.22045e-15  -2.9976e-15
  9.99201e-16  -2.9976e-15    4.44089e-16

julia> CLSEIG
3-element Array{Complex{Float64},1}:
   -4.37703628399912 + 2.8107164873731247im
   -4.37703628399912 - 2.8107164873731247im
 -1.8663764577096091 + 0.0im

julia> eigvals(A-B*F)
3-element Array{Complex{Float64},1}:
  -4.377036283999118 - 2.8107164873731234im
  -4.377036283999118 + 2.8107164873731234im
 -1.8663764577096063 + 0.0im
```
"""
function arec(A::AbstractMatrix, B::AbstractVecOrMat, R::Union{AbstractMatrix,UniformScaling,Real,Complex},
   Q::Union{AbstractMatrix,UniformScaling,Real,Complex}, S::AbstractVecOrMat = zeros(eltype(B),size(B));
   as = false, rtol::Real = size(A,1)*eps(real(float(one(eltype(A))))), orth = false)
   if orth
      return garec(A, I, B, 0, R, Q, S; as = as, rtol = rtol)
   else
      return arec(A, B, 0, R, Q, S; as = as, rtol = rtol)
   end
end
"""
    arec(A, B, G, R, Q, S; as = false, rtol::Real = nϵ, orth = false) -> (X, EVALS, F, Z)

Computes `X`, the hermitian/symmetric stabilizing solution (if `as = false`) or
anti-stabilizing solution (if `as = true`) of the continuous-time
algebraic Riccati equation

     A'X + XA - XGX - (XB+S)R^(-1)(B'X+S') + Q = 0,

where `G`, `R` and `Q` are hermitian/symmetric matrices or uniform scaling operators such that `R` is nonsingular.
Scalar-valued `G`, `R` and `Q` are interpreted as appropriately sized uniform scaling operators `G*I`, `R*I` and `Q*I`.

By default, the lower bound for the 1-norm reciprocal condition number `rtol` is `n*ϵ`, where `n` is the order of `A`
and `ϵ` is the _machine epsilon_ of the element type of `A`.

`EVALS` is a vector containing the (stable or anti-stable) eigenvalues of `A-BF-GX`.
`F` is the stabilizing or anti-stabilizing gain matrix `F = R^(-1)(B'X+S')`.
`Z = [ U; V; W ]` is a basis for the relevant stable/anti-stable deflating subspace such that `X = V/U` and `F = -W/U`.
An orthogonal basis Z can be determined, with an increased computational cost, by setting `orth = true`.

`Reference:`
Laub, A.J., A Schur Method for Solving Algebraic Riccati equations.
IEEE Trans. Auto. Contr., AC-24, pp. 913-921, 1979.
"""
function arec(A::AbstractMatrix, B::AbstractVecOrMat, G::Union{AbstractMatrix,UniformScaling,Real,Complex},
              R::Union{AbstractMatrix,UniformScaling,Real,Complex}, Q::Union{AbstractMatrix,UniformScaling,Real,Complex},
              S::AbstractVecOrMat; as = false, rtol::Real = size(A,1)*eps(real(float(one(eltype(A))))), orth = false)
    orth && (return garec(A, I, B, G, R, Q, S; as = as, rtol = rtol))

    n = LinearAlgebra.checksquare(A)
    T = promote_type( eltype(A), eltype(B), eltype(G), eltype(Q), eltype(R), eltype(S) )
    typeof(B) <: AbstractVector ? (nb, m) = (length(B), 1) : (nb, m) = size(B)
    n == nb || throw(DimensionMismatch("B must be a matrix with row dimension $n or a vector of length $n"))
    if typeof(G) <: AbstractArray
      (LinearAlgebra.checksquare(G) == n && ishermitian(G)) ||
          throw(DimensionMismatch("G must be a symmetric/hermitian matrix of dimension $n"))
   else
     G = G*I
     iszero(imag(G.λ)) || throw("G must be a symmetric/hermitian matrix")
   end
   if typeof(R) <: AbstractArray
      (LinearAlgebra.checksquare(R) == m && ishermitian(R)) ||
         throw(DimensionMismatch("R must be a symmetric/hermitian matrix of dimension $m"))
   else
      R = R*I
      iszero(imag(R.λ)) || throw("R must be a symmetric/hermitian matrix")
   end
   if typeof(Q) <: AbstractArray
     (LinearAlgebra.checksquare(Q) == n && ishermitian(Q)) ||
        throw(DimensionMismatch("Q must be a symmetric/hermitian matrix of dimension $n"))
   else
     Q = Q*I
     iszero(imag(Q.λ)) || throw("Q must be a symmetric/hermitian matrix")
   end
   typeof(S) <: AbstractVector ? (ns, ms) = (length(S), 1) : (ns, ms) = size(S)
   (n == ns && m == ms) ||
      throw(DimensionMismatch("S must be a $n x $m matrix or a vector of length $n"))
   T <: BlasFloat || (T = promote_type(Float64,T))
   TR = real(T)
   epsm = eps(TR)
   eltype(A) == T || (A = convert(Matrix{T},A))
   eltype(B) == T || (typeof(B) <: AbstractVector ? B = convert(Vector{T},B) : B = convert(Matrix{T},B))
   if typeof(G) <: AbstractArray
       (LinearAlgebra.checksquare(G) == n && ishermitian(G)) ||
          throw(DimensionMismatch("G must be a symmetric/hermitian matrix of dimension $n"))
   else
      G = G*I
      iszero(imag(G.λ)) || throw("G must be a symmetric/hermitian matrix")
   end
   if typeof(R) <: AbstractArray
      (LinearAlgebra.checksquare(R) == m && ishermitian(R)) ||
         throw(DimensionMismatch("R must be a symmetric/hermitian matrix of dimension $m"))
   else
      R = R*I
      iszero(imag(R.λ)) || throw("R must be a symmetric/hermitian matrix")
   end
   if typeof(Q) <: AbstractArray
     (LinearAlgebra.checksquare(Q) == n && ishermitian(Q)) ||
        throw(DimensionMismatch("Q must be a symmetric/hermitian matrix of dimension $n"))
   else
      Q = Q*I
      iszero(imag(Q.λ)) || throw("Q must be a symmetric/hermitian matrix")
   end
   if eltype(S) != T
      if typeof(S) <: AbstractVector
         S = convert(Vector{T},S)
      else
         S = convert(Matrix{T},S)
      end
   end

   n == 0 && (return  zeros(T,0,0), zeros(T,0), zeros(T,m,0), zeros(T,m,0) )

   S0flag = iszero(S)
   typeof(R) <: UniformScaling && (R = Matrix{T}(R,m,m))
   SR = schur(R)
   D = real(diag(SR.T))
   Da = abs.(D)
   minDa, = findmin(Da)
   maxDa, = findmax(Da)
   minDa <= epsm*maxDa && error("R must be non-singular")
   if minDa > sqrt(epsm)*maxDa && maxDa > 100*eps(max(opnorm(A,1),opnorm(G,1),opnorm(Q,1)))
      #Dinv = diagm(0 => 1 ./ D)
      Dinv = Diagonal(1 ./ D)
      Bu = B*SR.Z
      #G = G + Bu*Dinv*Bu'
      #G = utqu(Dinv,Bu')
      G += utqu(Dinv,Bu')
      if S0flag
         sol = arec(A,G,Q; as = as, rtol = rtol)
         w2 = SR.Z*Dinv*Bu'
         f = w2*sol[1]
         z = [sol[3]; w2*(sol[3])[n+1:end,:]]
      else
         Su = S*SR.Z
         #Q -= Su*Dinv*Su'
         Q -= utqu(Dinv,Su')
         sol = arec(A-Bu*Dinv*Su',G,Q; as = as, rtol = rtol)
         w1 = SR.Z*Dinv*Su'
         w2 = SR.Z*Dinv*Bu'
         f = w1+w2*sol[1]
         #f = SR.Z*Dinv*(Bu'*sol[1]+Su')
         z = [sol[3]; [w1 w2]*sol[3] ]
      end
      return sol[1], sol[2], f, z
   else
      # use implicit form
      @warn "R nearly singular: using the orthogonal reduction method"
      return garec(A, I, B, G, R, Q, S; as = as, rtol = rtol)
   end
end
"""
    garec(A, E, G, Q = 0; as = false, rtol::Real = nϵ) -> (X, EVALS, Z)

Compute `X`, the hermitian/symmetric stabilizing solution (if `as = false`) or
anti-stabilizing solution (if `as = true`) of the generalized continuous-time
algebraic Riccati equation

    A'XE + E'XA - E'XGXE + Q = 0,

where `G` and `Q` are hermitian/symmetric matrices or uniform scaling operators and `E` is a nonsingular matrix.
Scalar-valued `G` and `Q` are interpreted as appropriately sized uniform scaling operators `G*I` and `Q*I`.

By default, the lower bound for the 1-norm reciprocal condition number `rtol` is `n*ϵ`, where `n` is the order of `A`
and `ϵ` is the _machine epsilon_ of the element type of `A`.

`EVALS` is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair `(A-GXE,E)`.
`Z = [ U; V ]` is an orthogonal basis for the stable/anti-stable deflating subspace such that `X = V/(EU)`.

`Reference:`
W.F. Arnold, III and A.J. Laub,
Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations,
Proc. IEEE, 72:1746-1754, 1984.
"""
function garec(A::AbstractMatrix, E::Union{AbstractMatrix,UniformScaling}, G::Union{AbstractMatrix,UniformScaling,Real,Complex},
               Q::Union{AbstractMatrix,UniformScaling,Real,Complex} = zero(eltype(A));
               as = false, rtol::Real = size(A,1)*eps(real(float(one(eltype(A))))))
    n = LinearAlgebra.checksquare(A)
    T = promote_type( eltype(A), eltype(G), eltype(Q) )
    eident = (E == I)
    if !eident
       LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a $n x $n matrix or I"))
       eident = isequal(E,I)
       eident ? (E = I) : (T = promote_type(T,eltype(E)))
    end
    if typeof(G) <: AbstractArray
       (LinearAlgebra.checksquare(G) == n && ishermitian(G)) ||
          throw(DimensionMismatch("G must be a symmetric/hermitian matrix of dimension $n"))
    else
      G = G*I
      iszero(imag(G.λ)) || throw("G must be a symmetric/hermitian matrix")
    end
    if typeof(Q) <: AbstractArray
       (LinearAlgebra.checksquare(Q) == n && ishermitian(Q)) ||
          throw(DimensionMismatch("Q must be a symmetric/hermitian matrix of dimension $n"))
    else
       Q = Q*I
       iszero(imag(Q.λ)) || throw("Q must be a symmetric/hermitian matrix")
    end
    T <: BlasFloat || (T = promote_type(Float64,T))
    TR = real(T)
    epsm = eps(TR)
    eltype(A) == T || (A = convert(Matrix{T},A))
    eident || eltype(E) == T || (E = convert(Matrix{T},E))
    eltype(G) == T || (typeof(G) <: AbstractMatrix ? G = convert(Matrix{T},G) : G = convert(T,G.λ)*I)
    eltype(Q) == T || (typeof(Q) <: AbstractMatrix ? Q = convert(Matrix{T},Q) : Q = convert(T,Q.λ)*I)

    n == 0 && (return  zeros(T,0,0), zeros(T,0), zeros(T,m,0) )

    if !eident
       Et = LinearAlgebra.LAPACK.getrf!(copy(E))
       LinearAlgebra.LAPACK.gecon!('1',Et[1],opnorm(E,1))  < epsm && error("E must be non-singular")
    end
    #  Method:  A stable/anti-stable deflating subspace Z1 = [Z11; Z21] of the pencil
    #       L -s P := [  A  -G ]  - s [ E  0  ]
    #                 [ -Q  -A']      [ 0  E' ]
    #  is determined and the solution X is computed as X = Z21*inv(E*Z11).
    L = [ A -G; -Q -A']
    P = [ E zeros(T,n,n); zeros(T,n,n) E']
    LPS = schur(L,P)
    as ? select = real.(LPS.α ./ LPS.β) .> 0 : select = real.(LPS.α ./ LPS.β) .< 0
    n == length(filter(y-> y == true,select)) ||
       error("The Hamiltonian/skew-Hamiltonian pencil is not dichotomic")
    ordschur!(LPS, select)
    i1 = 1:n
    i2 = n+1:2n
    F = _LUwithRicTest(LPS.Z[i1, i1],rtol)
    if eident
       x = LPS.Z[i2,i1]/F
    else
       x = LPS.Z[i2,i1]/(E*LPS.Z[i1,i1])
    end

    return  (x+x')/2, LPS.values[i1], LPS.Z[:,i1]
end
"""
    garec(A, E, B, R, Q, S; as = false, rtol::Real = nϵ) -> (X, EVALS, F, Z)

Compute `X`, the hermitian/symmetric stabilizing solution (if `as = false`) or
anti-stabilizing solution (if `as = true`) of the generalized continuous-time
algebraic Riccati equation

    A'XE + E'XA - (E'XB+S)R^(-1)(B'XE+S') + Q = 0,

where `R` and `Q` are hermitian/symmetric matrices such that `R` is nonsingular, and
`E` is a nonsingular matrix.
Scalar-valued `R` and `Q` are interpreted as appropriately sized uniform scaling operators `R*I` and `Q*I`.
`S`, if not specified, is set to `S = zeros(size(B))`.

By default, the lower bound for the 1-norm reciprocal condition number `rtol` is `n*ϵ`, where `n` is the order of `A`
and `ϵ` is the _machine epsilon_ of the element type of `A`.

`EVALS` is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair `(A-BF,E)`.
`F` is the stabilizing/anti-stabilizing gain matrix `F = R^(-1)(B'XE+S')`.
`Z = [ U; V; W ]` is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that `X = V/(EU)` and `F = -W/U`.

`Reference:`
W.F. Arnold, III and A.J. Laub,
Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations,
Proc. IEEE, 72:1746-1754, 1984.

# Example
```jldoctest
julia> using LinearAlgebra

julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
 -6.0  -2.0   1.0
  5.0   1.0  -1.0
 -4.0  -2.0  -1.0

julia> E = [10. 3. 0.; 0. 5. -1.; 0. 0. 10.]
3×3 Array{Float64,2}:
 10.0  3.0   0.0
  0.0  5.0  -1.0
  0.0  0.0  10.0

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

julia> R = [1. 0.; 0. 5.]
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  5.0

julia> X, CLSEIG, F = garec(A,E,B,R,2I);

julia> X
3×3 Array{Float64,2}:
  0.0502214   0.0284089   -0.0303703
  0.0284089   0.111219    -0.00259162
 -0.0303703  -0.00259162   0.0618395

julia> A'*X*E+E'*X*A-E'*X*B*inv(R)*B'*X*E+2I
3×3 Array{Float64,2}:
  1.55431e-15  -1.9984e-15   -3.33067e-15
 -1.77636e-15   1.33227e-15  -3.33067e-15
 -2.88658e-15  -3.21965e-15   1.11022e-15

julia> CLSEIG
3-element Array{Complex{Float64},1}:
  -0.6184265391601464 + 0.2913286844595737im
  -0.6184265391601464 - 0.2913286844595737im
 -0.21613059964451786 + 0.0im

julia> eigvals(A-B*F,E)
3-element Array{Complex{Float64},1}:
 -0.6184265391601462 - 0.29132868445957383im
 -0.6184265391601462 + 0.2913286844595739im
  -0.216130599644518 + 0.0im
```
"""
function garec(A::AbstractMatrix, E::Union{AbstractMatrix,UniformScaling}, B::AbstractVecOrMat, R::Union{AbstractMatrix,UniformScaling,Real,Complex},
   Q::Union{AbstractMatrix,UniformScaling,Real,Complex}, S::AbstractVecOrMat = zeros(eltype(B),size(B));
   as = false, rtol::Real = size(A,1)*eps(real(float(one(eltype(A))))))
   garec(A, E, B, 0, R, Q, S; as = as, rtol = rtol)
end
"""
    garec(A, E, B, G, R, Q, S; as = false, rtol::Real = nϵ) -> (X, EVALS, F, Z)

Compute `X`, the hermitian/symmetric stabilizing solution (if `as = false`) or
anti-stabilizing solution (if `as = true`) of the generalized continuous-time
algebraic Riccati equation

    A'XE + E'XA - E'XGXE - (E'XB+S)R^(-1)(B'XE+S') + Q = 0,

where `G`, `Q` and `R` are hermitian/symmetric matrices such that `R` is nonsingular, and
`E` is a nonsingular matrix.
Scalar-valued `G`, `R` and `Q` are interpreted as appropriately sized uniform scaling operators `G*I`, `R*I` and `Q*I`.

By default, the lower bound for the 1-norm reciprocal condition number `rtol` is `n*ϵ`, where `n` is the order of `A`
and `ϵ` is the _machine epsilon_ of the element type of `A`.

`EVALS` is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair `(A-BF-GXE,E)`.
`F` is the stabilizing/anti-stabilizing gain matrix `F = R^(-1)(B'XE+S')`.
`Z = [ U; V; W ]` is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that `X = V/(EU)` and `F = -W/U`.

`Reference:`
W.F. Arnold, III and A.J. Laub,
Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations,
Proc. IEEE, 72:1746-1754, 1984.
"""
function garec(A::AbstractMatrix, E::Union{AbstractMatrix,UniformScaling}, B::AbstractVecOrMat,
               G::Union{AbstractMatrix,UniformScaling,Real,Complex}, R::Union{AbstractMatrix,UniformScaling,Real,Complex},
               Q::Union{AbstractMatrix,UniformScaling,Real,Complex}, S::AbstractVecOrMat;
               as = false, rtol::Real = size(A,1)*eps(real(float(one(eltype(A))))))
    n = LinearAlgebra.checksquare(A)
    T = promote_type( eltype(A), eltype(B), eltype(G), eltype(Q), eltype(R), eltype(S) )
    typeof(B) <: AbstractVector ? (nb, m) = (length(B), 1) : (nb, m) = size(B)
    n == nb || throw(DimensionMismatch("B must be a matrix with row dimension $n or a vector of length $n"))
    eident = (E == I)
    if !eident
       LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a $n x $n matrix or I"))
       eident = isequal(E,I)
       eident ? (E = I) : (T = promote_type(T,eltype(E)))
    end
    if typeof(G) <: AbstractArray
       (LinearAlgebra.checksquare(G) == n && ishermitian(G)) ||
           throw(DimensionMismatch("G must be a symmetric/hermitian matrix of dimension $n"))
    else
      G = G*I
      iszero(imag(G.λ)) || throw("G must be a symmetric/hermitian matrix")
    end
    if typeof(R) <: AbstractArray
       (LinearAlgebra.checksquare(R) == m && ishermitian(R)) ||
          throw(DimensionMismatch("R must be a symmetric/hermitian matrix of dimension $m"))
    else
       R = R*I
       iszero(imag(R.λ)) || throw("R must be a symmetric/hermitian matrix")
    end
    if typeof(Q) <: AbstractArray
      (LinearAlgebra.checksquare(Q) == n && ishermitian(Q)) ||
         throw(DimensionMismatch("Q must be a symmetric/hermitian matrix of dimension $n"))
    else
      Q = Q*I
      iszero(imag(Q.λ)) || throw("Q must be a symmetric/hermitian matrix")
    end
    typeof(S) <: AbstractVector ? (ns, ms) = (length(S), 1) : (ns, ms) = size(S)
    (n == ns && m == ms) || throw(DimensionMismatch("S must be a $n x $m matrix or a vector of length $n"))
    T <: BlasFloat || (T = promote_type(Float64,T))
    TR = real(T)
    epsm = eps(TR)
    eltype(A) == T || (A = convert(Matrix{T},A))
    eident || eltype(E) == T || (E = convert(Matrix{T},E))
    if !eident
       Et = LinearAlgebra.LAPACK.getrf!(copy(E))
       LinearAlgebra.LAPACK.gecon!('1',Et[1],opnorm(E,1))  < epsm && error("E must be non-singular")
    end
    eltype(B) == T || (typeof(B) <: AbstractVector ? B = convert(Vector{T},B) : B = convert(Matrix{T},B))
    eltype(G) == T || (typeof(G) <: AbstractMatrix ? G = convert(Matrix{T},G) : G = convert(T,G.λ)*I)
    eltype(Q) == T || (typeof(Q) <: AbstractMatrix ? Q = convert(Matrix{T},Q) : Q = convert(T,Q.λ)*I)
    eltype(R) == T || (typeof(R) <: AbstractMatrix ? R = convert(Matrix{T},R) : R = convert(T,R.λ)*I)
    eltype(S) == T || (typeof(S) <: AbstractVector ? S = convert(Vector{T},S) : S = convert(Matrix{T},S))

    n == 0 && (return  zeros(T,0,0), zeros(T,0), zeros(T,m,0), zeros(T,m,0) )

    cond(R)*epsm < 1 || error("R must be non-singular")

    #  Method:  A stable/ant-stable deflating subspace Z1 = [Z11; Z21; Z31] of the pencil
    #               [  A  -G    B ]      [ E  0  0 ]
    #      L -s P = [ -Q  -A'  -S ]  - s [ 0  E' 0 ]
    #               [  S'  B'   R ]      [ 0  0  0 ]
    # is determined and the solution X and feedback F are computed as
    #          X = Z21*inv(E*Z11),   F = -Z31*inv(Z11).

    #deflate m simple infinite eigenvalues
    n2 = n+n;
    W = qr(Matrix([S; B; R]));
    cond(W.R) * epsm  < 1 || error("The extended Hamiltonian/skew-Hamiltonian pencil is not regular")

    #z = W.Q[:,m+1:m+n2]
    z = W.Q*[fill(false,m,n2); I ]

    iric = 1:n2
    i1 = 1:n
    i2 = n+1:n2
    i3 = n2+1:n2+m
    L11 = [ A -G B; -Q -A' -S]*z
    P11 = [ E*z[i1,:]; E'*z[i2,:] ]
    LPS = schur(L11,P11)
    as ? select = real.(LPS.α ./ LPS.β) .> 0 : select = real.(LPS.α ./ LPS.β) .< 0
    n == length(filter(y-> y == true,select)) ||
         error("The extended Hamiltonian/skew-Hamiltonian pencil is not dichotomic")
    ordschur!(LPS, select)

    z[:,i1] = z[:,iric]*LPS.Z[:,i1];

    F = _LUwithRicTest(z[i1,i1],rtol)
    if eident
       x = z[n+1:end,i1]/F
       f = -x[n+1:end,:]
       x = x[i1,:]
    else
       f = -z[i3,i1]/F
       x = z[i2,i1]/(E*z[i1,i1])
    end

    return  (x+x')/2, LPS.values[i1] , f, z[:,i1]
end

"""
    ared(A, B, R, Q, S; as = false, rtol::Real = nϵ) -> (X, EVALS, F, Z)

Compute `X`, the hermitian/symmetric stabilizing solution (if `as = false`) or
anti-stabilizing solution (if `as = true`) of the discrete-time algebraic Riccati equation

    A'XA - X - (A'XB+S)(R+B'XB)^(-1)(B'XA+S') + Q = 0,

where `R` and `Q` are hermitian/symmetric matrices.
Scalar-valued `R` and `Q` are interpreted as appropriately sized uniform scaling operators `R*I` and `Q*I`.
`S`, if not specified, is set to `S = zeros(size(B))`.

By default, the lower bound for the 1-norm reciprocal condition number `rtol` is `n*ϵ`, where `n` is the order of `A`
and `ϵ` is the _machine epsilon_ of the element type of `A`.

`EVALS` is a vector containing the (stable) eigenvalues of `A-BF`.
`F` is the stabilizing gain matrix `F = (R+B'XB)^(-1)(B'XA+S')`.
`Z = [ U; V; W ]` is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that `X = V/(EU)` and `F = -W/U`.

`Reference:`
W.F. Arnold, III and A.J. Laub,
Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations,
Proc. IEEE, 72:1746-1754, 1984.

# Example
```jldoctest
julia> using LinearAlgebra

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

julia> B = [ 0.; sqrt(2.) ]
2-element Array{Float64,1}:
 0.0
 1.4142135623730951

julia> R = 1.
1.0

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

julia> X, CLSEIG, F = ared(A,B,R,Q);

julia> X
2×2 Array{Float64,2}:
  1.0  -1.0
 -1.0   1.5

julia> A'*X*A-X-A'*X*B*inv(R+B'*X*B)*B'*X*A+Q
2×2 Array{Float64,2}:
  0.0          -3.33067e-16
 -3.33067e-16   8.88178e-16

julia> CLSEIG
2-element Array{Complex{Float64},1}:
 0.4999999999999998 - 0.0im
               -0.0 - 0.0im

julia> eigvals(A-B*F)
2-element Array{Float64,1}:
 -2.7755575615628914e-16
  0.5
```
"""
function ared(A::AbstractMatrix, B::AbstractVecOrMat, R::Union{AbstractMatrix,UniformScaling,Real,Complex},
              Q::Union{AbstractMatrix,UniformScaling,Real,Complex}, S::AbstractVecOrMat = zeros(eltype(B),size(B));
              as = false, rtol::Real = size(A,1)*eps(real(float(one(eltype(A))))))
    gared(A, I, B, R, Q, S; as = as, rtol = rtol)
end
"""
    gared(A, E, B, R, Q, S; as = false, rtol::Real = nϵ) -> (X, EVALS, F, Z)

Compute `X`, the hermitian/symmetric stabilizing solution (if `as = false`) or
anti-stabilizing solution (if `as = true`) of the generalized discrete-time
algebraic Riccati equation

    A'XA - E'XE - (A'XB+S)(R+B'XB)^(-1)(B'XA+S') + Q = 0,

where `R` and `Q` are hermitian/symmetric matrices, and `E` ist non-singular.
Scalar-valued `R` and `Q` are interpreted as appropriately sized uniform scaling operators `R*I` and `Q*I`.
`S`, if not specified, is set to `S = zeros(size(B))`.

By default, the lower bound for the 1-norm reciprocal condition number `rtol` is `n*ϵ`, where `n` is the order of `A`
and `ϵ` is the _machine epsilon_ of the element type of `A`.

`EVALS` is a vector containing the (stable or anti-stable) generalized eigenvalues of the pair `(A-BF,E)`.
`F` is the stabilizing/anti-stabilizing gain matrix `F = (R+B'XB)^(-1)(B'XA+S')`.
`Z = [ U; V; W ]` is an orthogonal basis for the relevant stable/anti-stable deflating subspace such that `X = V/(EU)` and `F = -W/U`.

`Reference:`
W.F. Arnold, III and A.J. Laub,
Generalized Eigenproblem Algorithms and Software for Algebraic Riccati Equations,
Proc. IEEE, 72:1746-1754, 1984.

# Example
```jldoctest
julia> using LinearAlgebra

julia> A = [-6. -2. 1.; 5. 1. -1; -4. -2. -1.]
3×3 Array{Float64,2}:
 -6.0  -2.0   1.0
  5.0   1.0  -1.0
 -4.0  -2.0  -1.0

julia> E = [10. 3. 0.; 0. 5. -1.; 0. 0. 10.]
3×3 Array{Float64,2}:
 10.0  3.0   0.0
  0.0  5.0  -1.0
  0.0  0.0  10.0

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

julia> R = [1. 0.; 0. 5.]
2×2 Array{Float64,2}:
 1.0  0.0
 0.0  5.0

julia> X, CLSEIG, F = gared(A,E,B,R,2I);

julia> X
3×3 Array{Float64,2}:
  0.065865   -0.0147205  -0.0100407
 -0.0147205   0.0885939   0.0101422
 -0.0100407   0.0101422   0.0234425

julia> A'*X*A-E'*X*E-A'*X*B*inv(R+B'*X*B)*B'*X*A+2I
3×3 Array{Float64,2}:
 -1.33227e-15  -2.48412e-15   1.38778e-16
 -2.498e-15    -4.44089e-16  -6.50521e-16
  1.80411e-16  -5.91541e-16  -1.33227e-15

julia> CLSEIG
3-element Array{Complex{Float64},1}:
  -0.084235615751339 - 0.0im
  -0.190533552034239 - 0.0im
 -0.5238922629921539 - 0.0im

julia> eigvals(A-B*F,E)
3-element Array{Float64,1}:
 -0.5238922629921539
 -0.19053355203423886
 -0.08423561575133902
```
"""
function gared(A::AbstractMatrix, E::Union{AbstractMatrix,UniformScaling}, B::AbstractVecOrMat,
               R::Union{AbstractMatrix,UniformScaling,Real,Complex}, Q::Union{AbstractMatrix,UniformScaling,Real,Complex},
               S::AbstractVecOrMat = zeros(eltype(B),size(B)); as = false, rtol::Real = size(A,1)*eps(real(float(one(eltype(A))))))
    n = LinearAlgebra.checksquare(A)
    T = promote_type( eltype(A), eltype(B), eltype(R), eltype(Q), eltype(S) )
    typeof(B) <: AbstractVector ? (nb, m) = (length(B), 1) : (nb, m) = size(B)
    n == nb || throw(DimensionMismatch("B must be a matrix with row dimension $n or a vector of length $n"))
    if typeof(E) <: UniformScaling{Bool}
       eident = true
       E = I
    else
       LinearAlgebra.checksquare(E) == n || throw(DimensionMismatch("E must be a $n x $n matrix or I"))
       eident = isequal(E,I)
       eident ? E = I : T = promote_type(T,eltype(E))
     end
    if typeof(R) <: AbstractArray
       (LinearAlgebra.checksquare(R) == m && ishermitian(R)) ||
         throw(DimensionMismatch("R must be a symmetric/hermitian matrix of dimension $m"))
    else
      R = R*I
      iszero(imag(R.λ)) || throw("R must be a symmetric/hermitian matrix")
    end
    if typeof(Q) <: AbstractArray
       (LinearAlgebra.checksquare(Q) == n && ishermitian(Q)) ||
          throw(DimensionMismatch("Q must be a symmetric/hermitian matrix of dimension $n"))
    else
       Q = Q*I
       iszero(imag(Q.λ)) || throw("Q must be a symmetric/hermitian matrix")
    end
    typeof(S) <: AbstractVector ? (ns, ms) = (length(S), 1) : (ns, ms) = size(S)
    (n == ns && m == ms) || throw(DimensionMismatch("S must be a $n x $m matrix or a vector of length $n"))
    T <: BlasFloat || (T = promote_type(Float64,T))
    TR = real(T)
    epsm = eps(TR)
    eltype(A) == T || (A = convert(Matrix{T},A))
    eident || eltype(E) == T || (E = convert(Matrix{T},E))
    if !eident
       Et = LinearAlgebra.LAPACK.getrf!(copy(E))
       LinearAlgebra.LAPACK.gecon!('1',Et[1],opnorm(E,1))  < epsm && error("E must be non-singular")
    end
    eltype(B) == T || (typeof(B) <: AbstractVector ? B = convert(Vector{T},B) : B = convert(Matrix{T},B))
    eltype(Q) == T || (typeof(Q) <: AbstractMatrix ? Q = convert(Matrix{T},Q) : Q = convert(T,Q.λ)*I)
    eltype(R) == T || (typeof(R) <: AbstractMatrix ? R = convert(Matrix{T},R) : R = convert(T,R.λ)*I)
    eltype(S) == T || (typeof(S) <: AbstractVector ? S = convert(Vector{T},S) : S = convert(Matrix{T},S))

    n == 0 && (return  zeros(T,0,0), zeros(T,0), zeros(T,m,0), zeros(T,m,0) )

    if !eident
      Et = LinearAlgebra.LAPACK.getrf!(copy(E))
      LinearAlgebra.LAPACK.gecon!('1',Et[1],opnorm(E,1))  < epsm && error("E must be non-singular")
    end
    #  Method:  A stable deflating subspace Z1 = [Z11; Z21; Z31] of the pencil
    #                   [  A   0    B ]      [ E  0  0 ]
    #          L -z P = [ -Q   E'  -S ]  - z [ 0  A' 0 ]
    #                   [ S'   0    R ]      [ 0 -B' 0 ]
    #  is computed and the solution X and feedback F are computed as
    #          X = Z21*inv(E*Z11),   F = Z31*inv(Z11).
    n2 = n+n;
    F = qr([A'; -B'])
    L2 = F.Q'*[-Q  E' -S; copy(S') zeros(T,m,n) R]
    P2 = [zeros(T,n,n) F.R zeros(T,n,m)]

    G = qr(L2[n+1:n+m,:]')
    cond(G.R) * epsm  < 1 || error("The extended symplectic pencil is not regular")
    z = (G.Q*I)[:,[m+1:m+n2; 1:m]]

    i1 = 1:n
    i2 = n+1:n2
    i3 = n2+1:n2+m
    iric = 1:n2

    L1 = [ A zeros(T,n,n) B; L2[i1,:]]*z
    P1 = [ E zeros(T,n,n+m); P2]*z

    as ? PLS = schur(L1[iric,iric],P1[iric,iric]) : PLS = schur(P1[iric,iric],L1[iric,iric])
    select = abs.(PLS.α) .> abs.(PLS.β)

    n == length(filter(y-> y == true,select)) || error("The extended simplectic pencil is not dichotomic")

    ordschur!(PLS, select)
    z[:,i1]= z[:,iric]*PLS.Z[:,i1]

    F = _LUwithRicTest(z[i1,i1],rtol)
    if eident
       x = z[n+1:end,i1]/F
       f = -x[n+1:end,:]
       x = x[i1,:]
    else
       f = -z[i3,i1]/F
       x = z[i2,i1]/(E*z[i1,i1])
    end

    as ? iev = i2 : iev = i1
    clseig = PLS.β[iev] ./ PLS.α[iev]
    if as && T <: Complex
      clseig =  conj(clseig)
    end
    return  (x+x')/2, clseig, f, z[:,i1]
end

max_stars_repo_name: andreasvarga/MatrixEquations

max_issues_repo_name: andreasvarga/MatrixEquations

max_forks_repo_name: andreasvarga/MatrixEquations


text:
"""
    UInt256

256 bit unsigned integer, only a few operations are implemented. Holds the hash
value of `Block` and `Transaction`.

# Example

```julia
double_sha256(block)
double_sha256(tx)

zero(UInt256)
UInt256("234")
BTCParser.to_unsigned((0x00, 0x01, 0x02, 0x03,
                       0x00, 0x01, 0x02, 0x03,
                       0x00, 0x01, 0x02, 0x03,
                       0x00, 0x01, 0x02, 0x03,
                       0x00, 0x01, 0x02, 0x03,
                       0x00, 0x01, 0x02, 0x03,
                       0x00, 0x01, 0x02, 0x03,
                       0x00, 0x01, 0x02, 0x03))
```
"""
primitive type UInt256 <: Unsigned 256 end

# promotion_rule(::Type{UInt256}, ::Type{UInt}) = UInt256

# print does not work, have to implement
# ndigits0z(x::UInt256)

Base.show(io::IOStream, x::UInt256) = show(io, "0x" * string(x, base = 16))

# Base.show(io::IOStream, x::Array{UInt256}) = show(io, reinterpret(UInt8, x))

# Base.print(io::IOStream, x::UInt256)        = print(io, string(x, base = 16))
# Base.print(io::IOStream, x::Array{UInt256}) = show(io, x)

# Base.display(io::IOStream, x::UInt256)        = show(x)
# Base.display(io::IOStream, x::Array{UInt256}) = print(x)

# TODO: pad is currently being ignored
function Base.string(x::UInt256; base = 10, pad::Int = 0, neg::Bool = false)

    # The below only works if base is a power of 2, for now restrict it to 16
    # only...
    @assert base == 16

    x |> to_byte_tuple |>
        y -> map(i -> string(y[i], base = base, pad = 2), 32:-1:1) |>
        z -> string(z...)
    # tmp1 = convert(Ptr{UInt8}, pointer_from_objref(x))
    # tmp3 = map(i -> string(unsafe_load(tmp1, i), base = base, pad = pad), 32:-1:1)
    # string(tmp3...)
end

# TODO: I don't think there is a way to do this with bitwise operators in julia
function Base.bswap(x::UInt256)
    reinterpret(UInt256, reinterpret(UInt8, UInt256[x])[end:-1:1])[1]
end

function Base.zero(x::UInt256)
    reinterpret(UInt256, zeros(UInt8, 32))[1]
end

function Base.zero(::Type{UInt256})
    reinterpret(UInt256, zeros(UInt8, 32))[1]
end

# this one takes ~ 2x as long
# function Base.bswap(x::UInt256)
#     convert(UInt256, convert(NTuple{32, UInt8}, x)[end:-1:1])
# end

# Required to use this as keys in Dict
# TODO: this is probably not good:
Base.hash(x::UInt256, h::UInt) = Base.hash(x % UInt64, h)
Base.hash(x::UInt256) = Base.hash(x, zero(UInt))

# Base.rem(x::UInt256, ::Type{T}) where T = reinterpret(T,  Int256[x])[1]
# TODO: big/little endian
function Base.rem(x::UInt256, ::Type{T}) where {T <: Integer}
    Ref(x) |>
        pointer_from_objref |>
        x -> convert(Ptr{T}, x) |>
        unsafe_load
end

# Should probably be `parse`
# TODO: untested!!!
# TODO: optimize!!!
function UInt256(x::AbstractString)
    n = length(x)
    @assert n <= 64
    if n <= 32
        to_unsigned((to_byte_typle(Base.parse(UInt128, x, base = 16))...,
                     to_byte_tuple(zero(UInt128))...))
    else
        to_unsigned(
                    (to_byte_tuple(Base.parse(UInt128, String(x[end - 31:end]),      base = 16))...,
                     to_byte_tuple(Base.parse(UInt128, String(x[       1:end - 32]), base = 16))...) )
    end
end

ext: jl


text:
function print_result(model, pkgs, stats, csr, adj, trans)
    if termination_status(model) == OPTIMAL
        for p = 1:pkgs
            println("Package $p")
            for k = 1:csr
                output = false
                for i = 1:stats, j = 1:stats
                    if adj[i, j] == 1 && trans[k][i, j] == 1
                        if value(model[:x][i, j, k, p]) == 1.0
                            if !output
                                println("\tCrowd Shipper $k")
                                output = true
                            end
                            println("\t\t$i -> $j")
                        end
                    end
                end
            end
        end
    else
        println("Model infeasible!")
        println()
    end
end

function print_element(model::Model, key::Symbol)
    element = model[key]
    for e in element
        println(e)
    end
end

max_stars_repo_name: FilippoRanza/crowd-shipping-models

max_issues_repo_name: FilippoRanza/crowd-shipping-models

max_forks_repo_name: FilippoRanza/crowd-shipping-models


text:
using FedDCD
include("./PrimalMethods.jl")
include("./DualMethods.jl")

# Softmax for mnist
fileTrain = "data/mnist.scale"
fileTest = "data/mnist.scale.t"
λ = 1e-2
participationRate = 0.3
localLr = 1e-2
numRounds = 100
# - FedAvg
μ = 0.0

RunFedAvgAndProx(
    fileTrain,
    fileTest,
    λ,
    μ,
    participationRate,
    localLr,
    numRounds,
    "results/exp1/FedAvg_logReg_mnist_lambda1e-2.txt"
)

# - FedProx
μ = 1e-4

RunFedAvgAndProx(
    fileTrain,
    fileTest,
    λ,
    μ,
    participationRate,
    localLr,
    numRounds,
    "results/exp1/FedProx_logReg_mnist_lambda1e-2.txt"
)

# - Scaffold
RunScaffold(
    fileTrain,
    fileTest,
    λ,
    participationRate,
    localLr,
    numRounds,
    "results/exp1/Scaffold_logReg_mnist_lambda1e-2.txt"
)

# Softmax for rcv1
fileTrain = "data/rcv1_train.binary"
fileTest = "data/rcv1_train.binary"
λ = 1e-3
participationRate = 0.3
localLr = 1e-1
numRounds = 100
# - FedAvg
μ = 0.0

RunFedAvgAndProx(
    fileTrain,
    fileTest,
    λ,
    μ,
    participationRate,
    localLr,
    numRounds,
    "results/exp1/FedAvg_logReg_rcv1_lambda1e-2.txt"
)

# - FedProx
μ = 1e-4

RunFedAvgAndProx(
    fileTrain,
    fileTest,
    λ,
    μ,
    participationRate,
    localLr,
    numRounds,
    "results/exp1/FedProx_logReg_rcv1_lambda1e-2.txt"
)

# - Scaffold
RunScaffold(
    fileTrain,
    fileTest,
    λ,
    participationRate,
    localLr,
    numRounds,
    "results/exp1/Scaffold_logReg_rcv1_lambda1e-2.txt"
)

# Softmax + MLP for mnist
fileTrain = "data/mnist.scale"
fileTest = "data/mnist.scale.t"
λ = 1e-2
participationRate = 0.3
localLr = 1e-2
numRounds = 100
# - FedAvg
μ = 0.0

# - FedAvg
RunFedAvgAndProxNN(
    fileTrain,
    fileTest,
    λ,
    μ,
    participationRate,
    localLr,
    numRounds,
    "results/exp1/FedAvg_MLP_mnist_lambda1e-2.txt"
)

# - FedProx

μ = 1e-4
RunFedAvgAndProxNN(
    fileTrain,
    fileTest,
    λ,
    μ,
    participationRate,
    localLr,
    numRounds,
    "results/exp1/FedProx_MLP_mnist_lambda1e-2.txt"
)

# - Scaffold
RunScaffoldNN(
    fileTrain,
    fileTest,
    λ,
    participationRate,
    localLr,
    numRounds,
    "results/exp1/Scaffold_MLP_mnist_lambda1e-2.txt"
)

# Softmax for mnist
fileTrain = "data/mnist.scale"
fileTest = "data/mnist.scale.t"
λ = 1e-2
participationRate = 0.1
localLr = 1e-3
numRounds = 500
# - FedAvg
μ = 0.0

RunFedAvgAndProx(
    fileTrain,
    fileTest,
    λ,
    μ,
    participationRate,
    localLr,
    numRounds,
    "results/exp2/FedAvg_logReg_mnist_lambda1e-2_tau1e-1.txt"
)

# - FedProx
μ = 1e-4

RunFedAvgAndProx(
    fileTrain,
    fileTest,
    λ,
    μ,
    participationRate,
    localLr,
    numRounds,
    "results/exp2/FedProx_logReg_mnist_lambda1e-2_tau1e-1.txt"
)

# - Scaffold
RunScaffold(
    fileTrain,
    fileTest,
    λ,
    participationRate,
    localLr,
    numRounds,
    "results/exp2/Scaffold_logReg_mnist_lambda1e-2_tau1e-1.txt"
)

# Softmax for rcv1
fileTrain = "data/rcv1_train.binary"
fileTest = "data/rcv1_train.binary"
λ = 1e-3
participationRate = 0.1
localLr = 1e-1
numRounds = 500
# - FedAvg
μ = 0.0

RunFedAvgAndProx(
    fileTrain,
    fileTest,
    λ,
    μ,
    participationRate,
    localLr,
    numRounds,
    "results/exp2/FedAvg_logReg_rcv1_lambda1e-2_tau1e-1.txt"
)

# - FedProx
μ = 1e-4

RunFedAvgAndProx(
    fileTrain,
    fileTest,
    λ,
    μ,
    participationRate,
    localLr,
    numRounds,
    "results/exp2/FedProx_logReg_rcv1_lambda1e-2_tau1e-1.txt"
)

# - Scaffold
RunScaffold(
    fileTrain,
    fileTest,
    λ,
    participationRate,
    localLr,
    numRounds,
    "results/exp2/Scaffold_logReg_rcv1_lambda1e-2_tau1e-1.txt"
)

# Softmax + MLP for mnist
fileTrain = "data/mnist.scale"
fileTest = "data/mnist.scale.t"
λ = 1e-2
participationRate = 0.1
localLr = 1e-2
numRounds = 500
# - FedAvg
μ = 0.0

# - FedAvg
RunFedAvgAndProxNN(
    fileTrain,
    fileTest,
    λ,
    μ,
    participationRate,
    localLr,
    numRounds,
    "results/exp2/FedAvg_MLP_mnist_lambda1e-2_tau1e-1.txt"
)

# - FedProx

μ = 1e-4
RunFedAvgAndProxNN(
    fileTrain,
    fileTest,
    λ,
    μ,
    participationRate,
    localLr,
    numRounds,
    "results/exp2/FedProx_MLP_mnist_lambda1e-2_tau1e-1.txt"
)

# - Scaffold
RunScaffoldNN(
    fileTrain,
    fileTest,
    λ,
    participationRate,
    localLr,
    numRounds,
    "results/exp2/Scaffold_MLP_mnist_lambda1e-2_tau1e-1.txt"
)

#######################################################################################################
##################################
#     Participation rate = 0.05
##################################
# Softmax for mnist
fileTrain = "data/mnist.scale"
fileTest = "data/mnist.scale.t"
λ = 1e-2
participationRate = 0.05
localLr = 1e-3
numRounds = 1000
# - FedAvg
μ = 0.0

RunFedAvgAndProx(
    fileTrain,
    fileTest,
    λ,
    μ,
    participationRate,
    localLr,
    numRounds,
    "results/exp2/FedAvg_logReg_mnist_lambda1e-2_tau5e-2.txt"
)

# - FedProx
μ = 1e-4

RunFedAvgAndProx(
    fileTrain,
    fileTest,
    λ,
    μ,
    participationRate,
    localLr,
    numRounds,
    "results/exp2/FedProx_logReg_mnist_lambda1e-2_tau5e-2.txt"
)

# - Scaffold
RunScaffold(
    fileTrain,
    fileTest,
    λ,
    participationRate,
    localLr,
    numRounds,
    "results/exp2/Scaffold_logReg_mnist_lambda1e-2_tau5e-2.txt"
)

# Softmax for rcv1
fileTrain = "data/rcv1_train.binary"
fileTest = "data/rcv1_train.binary"
λ = 1e-3
participationRate = 0.05
localLr = 1e-1
numRounds = 1000
# - FedAvg
μ = 0.0

RunFedAvgAndProx(
    fileTrain,
    fileTest,
    λ,
    μ,
    participationRate,
    localLr,
    numRounds,
    "results/exp2/FedAvg_logReg_rcv1_lambda1e-2_tau5e-2.txt"
)

# - FedProx
μ = 1e-4

RunFedAvgAndProx(
    fileTrain,
    fileTest,
    λ,
    μ,
    participationRate,
    localLr,
    numRounds,
    "results/exp2/FedProx_logReg_rcv1_lambda1e-2_tau5e-2.txt"
)

# - Scaffold
RunScaffold(
    fileTrain,
    fileTest,
    λ,
    participationRate,
    localLr,
    numRounds,
    "results/exp2/Scaffold_logReg_rcv1_lambda1e-2_tau5e-2.txt"
)

# Softmax + MLP for mnist
fileTrain = "data/mnist.scale"
fileTest = "data/mnist.scale.t"
λ = 1e-2
participationRate = 0.05
localLr = 1e-2
numRounds = 1000
# - FedAvg
μ = 0.0

# - FedAvg
RunFedAvgAndProxNN(
    fileTrain,
    fileTest,
    λ,
    μ,
    participationRate,
    localLr,
    numRounds,
    "results/exp2/FedAvg_MLP_mnist_lambda1e-2_tau5e-2.txt"
)

# - FedProx

μ = 1e-4
RunFedAvgAndProxNN(
    fileTrain,
    fileTest,
    λ,
    μ,
    participationRate,
    localLr,
    numRounds,
    "results/exp2/FedProx_MLP_mnist_lambda1e-2_tau5e-2.txt"
)

# - Scaffold
RunScaffoldNN(
    fileTrain,
    fileTest,
    λ,
    participationRate,
    localLr,
    numRounds,
    "results/exp2/Scaffold_MLP_mnist_lambda1e-2_tau5e-2.txt"
)

# Softmax for mnist
fileTrain = "data/mnist.scale"
fileTest = "data/mnist.scale.t"
λ = 1e-2
participationRate = 0.3
localLr = 1e-2
numRounds = 100

# Experiment 2
#   rcv1
#       τ = 0.3
RunFedDCD(
    "data/rcv1_train.binary",
    "data/rcv1_train.binary",
    1e-3,
    0.3,
    0.1,
    100,
    "results/exp2/FedDCD_rcv1_lambda1e-3_tau3e-1.txt"
)
RunAccFedDCD(
    "data/rcv1_train.binary",
    "data/rcv1_train.binary",
    1e-3,
    0.3,
    0.1,
    100,
    "results/exp2/AccFedDCD_rcv1_lambda1e-3_tau3e-1.txt"
)
#       τ = 0.1
RunFedDCD(
    "data/rcv1_train.binary",
    "data/rcv1_train.binary",
    1e-3,
    0.1,
    0.1,
    500,
    "results/exp2/FedDCD_rcv1_lambda1e-3_tau1e-1.txt"
)
RunAccFedDCD(
    "data/rcv1_train.binary",
    "data/rcv1_train.binary",
    1e-3,
    0.1,
    0.1,
    500,
    "results/exp2/AccFedDCD_rcv1_lambda1e-3_tau1e-1.txt"
)
#       τ = 0.05
RunFedDCD(
    "data/rcv1_train.binary",
    "data/rcv1_train.binary",
    1e-3,
    0.05,
    0.1,
    1000,
    "results/exp2/FedDCD_rcv1_lambda1e-3_tau5e-2.txt"
)
RunAccFedDCD(
    "data/rcv1_train.binary",
    "data/rcv1_train.binary",
    1e-3,
    0.05,
    0.1,
    1000,
    "results/exp2/AccFedDCD_rcv1_lambda1e-3_tau5e-2.txt"
)

#   mnist
#       τ = 0.3
RunFedDCDNN(
    "data/mnist.scale",
    "data/mnist.scale.t",
    1e-2,
    0.3,
    0.1,
    100,
    "results/exp2/FedDCD_mnist_lambda1e-2_tau3e-1.txt"
)
RunAccFedDCDNN(
    "data/mnist.scale",
    "data/mnist.scale.t",
    1e-2,
    0.3,
    0.1,
    100,
    "results/exp2/AccFedDCD_mnist_lambda1e-2_tau3e-1.txt"
)
#       τ = 0.1
RunFedDCDNN(
    "data/mnist.scale",
    "data/mnist.scale.t",
    1e-2,
    0.1,
    0.1,
    500,
    "results/exp2/FedDCD_mnist_lambda1e-2_tau1e-1.txt"
)
RunAccFedDCDNN(
    "data/mnist.scale",
    "data/mnist.scale.t",
    1e-2,
    0.1,
    0.1,
    500,
    "results/exp2/AccFedDCD_mnist_lambda1e-2_tau1e-1.txt"
)
#       τ = 0.05
RunFedDCDNN(
    "data/mnist.scale",
    "data/mnist.scale.t",
    1e-2,
    0.05,
    0.1,
    1000,
    "results/exp2/FedDCD_mnist_lambda1e-2_tau5e-2.txt"
)
RunAccFedDCDNN(
    "data/mnist.scale",
    "data/mnist.scale.t",
    1e-2,
    0.05,
    0.1,
    1000,
    "results/exp2/AccFedDCD_mnist_lambda1e-2_tau5e-2.txt"
)


text:
include("CrossEntropyLog.jl")
include("jacobian.jl")
using Base.Test

X = rand(Float32, (5, 3))
t = rand(Float32, (5, 3))
l = CrossEntropyLog{Float32}()

# Test grad. wrt. input X
target!(l, t)
Dfwd = jacobian_fwd(l, X)
Dbwd = jacobian_bwd(l, X)

@test_approx_eq Dfwd Dbwd


text:
module QJuliaPipeNFCG

using QJuliaInterface
using QJuliaEnums
using QJuliaBlas
using QJuliaReduce
using QJuliaSolvers

using LinearAlgebra
using Printf

using Plots

norm2    = QJuliaReduce.gnorm2
rdot     = QJuliaReduce.reDotProduct
cpy      = QJuliaBlas.cpy

@enum NormType begin
	  preconditioned_norm
	  unpreconditioned_norm
	  natural_norm
  end

@enum TruncType begin
      trunc_type_standard # <- bad with the preconditioner
	  trunc_type_notay    # <- significantly better with a preconditioner
	end

@inline function MatPrecon(out::AbstractArray, inp::AbstractArray, outSloppy::AbstractArray, inpSloppy::AbstractArray, K::Function)

    if pointer_from_objref(out) == pointer_from_objref(inp); return; end #nothing to do

    outSloppy .=@. 0.0
    cpy(inpSloppy, inp)       #noop for the alias refs
    K(outSloppy, inpSloppy)
    cpy(out, outSloppy)       #noop for the alias refs

  end

  function solver(x::AbstractArray, b::AbstractArray, Mat::Any, MatSloppy::Any, param::QJuliaSolvers.QJuliaSolverParam_qj, K::Function, extra_args...)

    is_preconditioned = param.inv_type_precondition != QJuliaEnums.QJULIA_INVALID_INVERTER

    solver_name = "PipeNFCG2"

    normcheck  = unpreconditioned_norm
    truncstrat = trunc_type_notay #trunc_type_standard

    mmax = param.nKrylov

    println("Running ", solver_name ," solver with truncation length ", mmax)

    if is_preconditioned == true; println("Preconditioner  :: ", param.inv_type_precondition); end

    if (param.maxiter == 0)
      if param.use_init_guess == false; x .=@. 0.0; end
      return
    end #if param.maxiter == 0

    mixed = (param.dtype_sloppy != param.dtype)

    if mixed == true; error("Running mixed precision solver.");end

    local r = zeros(param.dtype_sloppy, length(x))
    local u = zeros(param.dtype_sloppy, length(x))
    local w = zeros(param.dtype_sloppy, length(x))
    local m = zeros(param.dtype_sloppy, length(x))
    local n = zeros(param.dtype_sloppy, length(x))
    # set of vectors
    local p  = zero(Matrix{param.dtype_sloppy}(undef, length(x), mmax+1))
    local s  = zero(Matrix{param.dtype_sloppy}(undef, length(x), mmax+1))
    local q  = zero(Matrix{param.dtype_sloppy}(undef, length(x), mmax+1))
    local z  = zero(Matrix{param.dtype_sloppy}(undef, length(x), mmax+1))
    #
    local beta  = zeros(Float64, mmax+1)
    local eta   = zeros(Float64, mmax+1)

    dp = 0.0; norm2B = norm(b)

    Precond(out, inp) = MatPrecon(out, inp, out, inp, K)

    # Compute cycle initial residual
    if param.use_init_guess == true
      Mat(r, x)
      r .=@. b - r
    else
      r .=@. b
    end
    # Initial stage
    Precond(u, r)
    Mat(w, u)
    #
    gamma   = dot(r, u)
    delta   = dot(w, u)
    eta[1]  = delta; alpha = gamma/delta
    #
    dp  = normcheck == preconditioned_norm || normcheck == unpreconditioned_norm ? norm(u) : abs(gamma)

    # m = B(w)
    Precond(m, w)
    # n = Am
    Mat(n, m)
    p[:,1] .=@. u
    s[:,1] .=@. w
    q[:,1] .=@. m
    z[:,1] .=@. n
    # update x, r, z, w as zero iteration
    x .=@. x + alpha*p[:,1]
    r .=@. r - alpha*s[:,1]
    u .=@. u - alpha*q[:,1]
    w .=@. w - alpha*z[:,1]
    #
    stop = dp*dp*param.tol*param.tol

    println(solver_name," : Initial residual = ", dp / norm2B)

    gamma   = dot(r, u)
    delta   = dot(w, u)
    beta[1] = dot(s[:,1], u)

    k = 1; converged = false;
    mi = truncstrat == trunc_type_notay ? 1 : mmax

    while (k < param.maxiter && converged == false)

      idx = k % (mmax + 1) + 1
      #m = u + B(w-r)
      n .=@. w - r
      Precond(m, n)
      m .=@. u + m
      #n = Am
      Mat(n, m)

      # finish all global comms here
      eta[idx] = 0.0
      j = 0
      for i in max(0,k-mi):(k-1)
        kdx = (i % (mmax+1)) + 1; j += 1
	beta[j] /= -eta[kdx]
	eta[idx] -= ((beta[j])*(beta[j])) * eta[kdx]
      end
      #
      eta[idx] += delta
      if(eta[idx] <= 0.)
        println("Restart due to square root breakdown or exact zero of eta at it = ", k)
        break
      else
        alpha = gamma/eta[idx]
      end

      # project out stored search directions
      p[:,idx] .=@. u
      s[:,idx] .=@. w
      q[:,idx] .=@. m
      z[:,idx] .=@. n

      j = 0
      for i in max(0,k-mi):(k-1)
        kdx = (i % (mmax+1)) + 1; j += 1
	#
	p[:,idx] .=@. p[:,idx] + beta[j] * p[:,kdx]
	s[:,idx] .=@. s[:,idx] + beta[j] * s[:,kdx]
	q[:,idx] .=@. q[:,idx] + beta[j] * q[:,kdx]
	z[:,idx] .=@. z[:,idx] + beta[j] * z[:,kdx]
      end

      # Update x, r, z, w
      x .=@. x + alpha*p[:,idx]
      r .=@. r - alpha*s[:,idx]
      u .=@. u - alpha*q[:,idx]
      w .=@. w - alpha*z[:,idx]

      mi = truncstrat == trunc_type_notay ? ((k) % mmax)+1 : mmax

      gamma = dot(r, u)
      delta = dot(w, u)
      j = 0
      for i in max(0,k-mi+1):k
        kdx = (i % (mmax+1)) + 1; j += 1
	beta[j] = dot(s[:,kdx], u)
      end

      dp = normcheck == preconditioned_norm ? norm(u) : (normcheck == unpreconditioned_norm ? norm(r) : sqrt(gamma))
      # Check for convergence
      @printf("%s: %d iteration, iter residual: %1.15e \n", solver_name, k, dp/norm2B)

      k += 1
    end # while context

    return
  end # solver context

end # module cobtext

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
module RegisterCore

using CenterIndexedArrays
using ImageCore, ImageFiltering
using Requires

import Base: +, -, *, /

export
    # types
    MismatchArray,
    NumDenom,
    ColonFun,
    PreprocessSNF,
    # functions
    highpass,
    indmin_mismatch,
    maxshift,
    mismatcharrays,
    ratio,
    separate,
    paddedview,
    trimmedview

"""
The major functions/types exported by this module are:

- `NumDenom` and `MismatchArray`: packed pair representation of
  `(num,denom)` mismatch data
- `separate`: splits a `NumDenom` array into its component `num,denom` arrays
- `indmin_mismatch`: find the location of the minimum mismatch
- `highpass`: highpass filter an image before performing registration
- `PreprocessSNF`: shot-noise standardization and filtering

"""
RegisterCore

"""
    x = NumDenom(num, denom)

`NumDenom{T}` is an object containing a `(num,denom)` pair.
`x.num` is `num` and `x.denom` is `denom`.

Algebraically, `NumDenom` objects act like 2-vectors, and can be added and multiplied
by scalars:

    nd1 + nd2 = NumDenom(nd1.num + nd2.num, nd1.denom + nd2.denom)
    2*nd      = NumDenom(2*nd.num, 2*nd.denom)

Note that this is *not* what you'd get from normal arithmetic with ratios, where, e.g.,
`2*nd` would be expected to produce `NumDenom(2*nd.num, nd.denom)`.
The reason for calling these `*` and `+` is for use in `Interpolations.jl`, because it
allows interpolation to be performed on "both arrays" at once without
recomputing the interpolation coefficients. See the documentation for information
about how this is used for performing aperturered mismatch computations.

As a consequence, there is no `convert(Float64, nd::NumDenom)` method, because the algebra
above breaks any pretense that `NumDenom` numbers are somehow equivalent to ratios.
If you want to convert to a ratio, see [`ratio`](@ref).
"""
struct NumDenom{T<:Number}
    num::T
    denom::T
end
NumDenom(n::Gray, d::Gray) = NumDenom(gray(n), gray(d))
NumDenom(n::Gray, d) = NumDenom(gray(n), d)
NumDenom(n, d::Gray) = NumDenom(n, gray(d))
NumDenom(n, d) = NumDenom(promote(n, d)...)

(+)(p1::NumDenom, p2::NumDenom) = NumDenom(p1.num+p2.num, p1.denom+p2.denom)
(-)(p1::NumDenom, p2::NumDenom) = NumDenom(p1.num-p2.num, p1.denom-p2.denom)
(*)(n::Number, p::NumDenom) = NumDenom(n*p.num, n*p.denom)
(*)(p::NumDenom, n::Number) = n*p
(/)(p::NumDenom, n::Number) = NumDenom(p.num/n, p.denom/n)
Base.oneunit(::Type{NumDenom{T}}) where {T} = NumDenom{T}(oneunit(T),oneunit(T))
Base.oneunit(p::NumDenom) = oneunit(typeof(p))
Base.zero(::Type{NumDenom{T}}) where {T} = NumDenom(zero(T),zero(T))
Base.zero(p::NumDenom) = zero(typeof(p))
Base.promote_rule(::Type{NumDenom{T1}}, ::Type{T2}) where {T1,T2<:Number} = NumDenom{promote_type(T1,T2)}
Base.promote_rule(::Type{NumDenom{T1}}, ::Type{NumDenom{T2}}) where {T1,T2} = NumDenom{promote_type(T1,T2)}
Base.eltype(::Type{NumDenom{T}}) where {T} = T
Base.convert(::Type{NumDenom{T}}, p::NumDenom{T}) where {T} = p
Base.convert(::Type{NumDenom{T}}, p::NumDenom) where {T} = NumDenom{T}(p.num, p.denom)

Base.round(::Type{NumDenom{T}}, p::NumDenom) where T = NumDenom{T}(round(T, p.num), round(T, p.denom))

Base.convert(::Type{F}, p::NumDenom) where F<:AbstractFloat = error("`convert($F, ::NumDenom)` is deliberately not defined, see `?NumDenom`.")

function Base.show(io::IO, p::NumDenom)
    print(io, "NumDenom(")
    show(io, p.num)
    print(io, ",")
    show(io, p.denom)
    print(io, ")")
    return
end

const MismatchArray{ND<:NumDenom,N,A} = CenterIndexedArray{ND,N,A}

"""
    mxs = maxshift(D)

Return the `maxshift` value used to compute the mismatch array `D`.
"""
maxshift(A::MismatchArray) = A.halfsize

"""
`numdenom = MismatchArray(num, denom)` packs the array-pair
`(num,denom)` into a single `MismatchArray`.  This is useful
preparation for interpolation.
"""
function (::Type{M})(num::AbstractArray, denom::AbstractArray) where M<:MismatchArray
    size(num) == size(denom) || throw(DimensionMismatch("num and denom must have the same size"))
    T = promote_type(eltype(num), eltype(denom))
    numdenom = CenterIndexedArray{NumDenom{T}}(undef, size(num))
    _packnd!(numdenom, num, denom)
end

function _packnd!(numdenom::AbstractArray, num::AbstractArray, denom::AbstractArray)
    Rnd, Rnum, Rdenom = eachindex(numdenom), eachindex(num), eachindex(denom)
    if Rnum == Rdenom
        for (Idest, Isrc) in zip(Rnd, Rnum)
            @inbounds numdenom[Idest] = NumDenom(num[Isrc], denom[Isrc])
        end
    elseif Rnd == Rnum
        for (Inum, Idenom) in zip(Rnum, Rdenom)
            @inbounds numdenom[Inum] = NumDenom(num[Inum], denom[Idenom])
        end
    else
        for (Ind, Inum, Idenom) in zip(Rnd, Rnum, Rdenom)
            @inbounds numdenom[Ind] = NumDenom(num[Inum], denom[Idenom])
        end
    end
    numdenom
end

function _packnd!(numdenom::CenterIndexedArray, num::CenterIndexedArray, denom::CenterIndexedArray)
    @simd for I in eachindex(num)
        @inbounds numdenom[I] = NumDenom(num[I], denom[I])
    end
    numdenom
end

# The next are mostly used just for testing
"""
`mms = mismatcharrays(nums, denoms)` packs array-of-arrays num/denom pairs as an array-of-MismatchArrays.

`mms = mismatcharrays(nums, denom)`, for `denom` a single array, uses the same `denom` array for all `nums`.
"""
function mismatcharrays(nums::AbstractArray{A}, denom::AbstractArray{T}) where {A<:AbstractArray,T<:Number}
    first = true
    local mms
    for i in eachindex(nums)
        num = nums[i]
        mm = MismatchArray(num, denom)
        if first
            mms = Array{typeof(mm)}(undef, size(nums))
            first = false
        end
        mms[i] = mm
    end
    mms
end

function mismatcharrays(nums::AbstractArray{A1}, denoms::AbstractArray{A2}) where {A1<:AbstractArray,A2<:AbstractArray}
    size(nums) == size(denoms) || throw(DimensionMismatch("nums and denoms arrays must have the same number of apertures"))
    first = true
    local mms
    for i in eachindex(nums, denoms)
        mm = MismatchArray(nums[i], denoms[i])
        if first
            mms = Array{typeof(mm)}(undef, size(nums))
            first = false
        end
        mms[i] = mm
    end
    mms
end

"""
`num, denom = separate(mm)` splits an `AbstractArray{NumDenom}` into separate
numerator and denominator arrays.
"""
function separate(data::AbstractArray{NumDenom{T}}) where T
    num = Array{T}(undef, size(data))
    denom = similar(num)
    for I in eachindex(data)
        nd = data[I]
        num[I] = nd.num
        denom[I] = nd.denom
    end
    num, denom
end

function separate(mm::MismatchArray)
    num, denom = separate(mm.data)
    CenterIndexedArray(num), CenterIndexedArray(denom)
end

function separate(mma::AbstractArray{M}) where M<:MismatchArray
    T = eltype(eltype(M))
    nums = Array{CenterIndexedArray{T,ndims(M)}}(undef, size(mma))
    denoms = similar(nums)
    for (i,mm) in enumerate(mma)
        nums[i], denoms[i] = separate(mm)
    end
    nums, denoms
end

"""
    r = ratio(nd::NumDenom, thresh, fillval=NaN)

Return `nd.num/nd.denom`, unless `nd.denom < thresh`, in which case return `fillval` converted
to the same type as the ratio.
Choosing a `thresh` of zero will always return the ratio.
"""
@inline function ratio(nd::NumDenom{T}, thresh, fillval=convert(T,NaN)) where {T}
    r = nd.num/nd.denom
    return nd.denom < thresh ? oftype(r, fillval) : r
end
ratio(r::Real, thresh, fillval=NaN) = r

(::Type{M})(::Type{T}, dims::Dims) where {M<:MismatchArray,T} = CenterIndexedArray{NumDenom{T}}(undef, dims)
(::Type{M})(::Type{T}, dims::Integer...) where {M<:MismatchArray,T} = CenterIndexedArray{NumDenom{T}}(undef, dims)

function Base.copyto!(M::MismatchArray, nd::Tuple{AbstractArray, AbstractArray})
    num, denom = nd
    size(M) == size(num) == size(denom) || error("all sizes must match")
    for (IM, Ind) in zip(eachindex(M), eachindex(num))
        M[IM] = NumDenom(num[Ind], denom[Ind])
    end
    M
end

#### Utility functions ####

"""
`index = indmin_mismatch(numdenom, thresh)` returns the location of
the minimum value of what is effectively `num./denom`.  However, it
considers only those points for which `denom .> thresh`; moreover, it
will never choose an edge point.  `index` is a CartesianIndex into the
arrays.
"""
function indmin_mismatch(numdenom::MismatchArray{NumDenom{T},N}, thresh::Real) where {T,N}
    imin = CartesianIndex(ntuple(d->0, Val(N)))
    rmin = typemax(T)
    threshT = convert(T, thresh)
    @inbounds for I in CartesianIndices(map(trimedges, axes(numdenom)))
        nd = numdenom[I]
        if nd.denom > threshT
            r = nd.num/nd.denom
            if r < rmin
                imin = I
                rmin = r
            end
        end
    end
    return imin
end

function indmin_mismatch(r::CenterIndexedArray{T,N}) where {T<:Number,N}
    imin = CartesianIndex(ntuple(d->0, Val(N)))
    rmin = typemax(T)
    @inbounds for I in CartesianIndices(map(trimedges, axes(r)))
        rval = r[I]
        if rval < rmin
            imin = I
            rmin = rval
        end
    end
    return imin
end

trimedges(r::AbstractUnitRange) = first(r)+1:last(r)-1

### Miscellaneous

"""
`datahp = highpass([T], data, sigma)` returns a highpass-filtered
version of `data`, with all negative values truncated at 0.  The
highpass is computed by subtracting a lowpass-filtered version of
data, using Gaussian filtering of width `sigma`.  As it is based on
`Image.jl`'s Gaussian filter, it gracefully handles `NaN` values.

If you do not wish to highpass-filter along a particular axis, put
`Inf` into the corresponding slot in `sigma`.

You may optionally specify the element type of the result, which for
`Integer` or `FixedPoint` inputs defaults to `Float32`.
"""
function highpass(::Type{T}, data::AbstractArray, sigma) where T
    if any(isinf, sigma)
        datahp = convert(Array{T,ndims(data)}, data)
    else
        datahp = data - imfilter(T, data, KernelFactors.IIRGaussian(T, (sigma...,)), NA())
    end
    datahp[datahp .< 0] .= 0  # truncate anything below 0
    datahp
end
highpass(data::AbstractArray{T}, sigma) where {T<:AbstractFloat} = highpass(T, data, sigma)
highpass(data::AbstractArray, sigma) = highpass(Float32, data, sigma)

"""
`pp = PreprocessSNF(bias, sigmalp, sigmahp)` constructs an object that
can be used to pre-process an image as `pp(img)`. The "SNF" part of
the name means "shot-noise filtered," meaning that this preprocessor
is specifically designed for situations in which you are dominated by
shot noise (i.e., from photon-counting statistics).

The processing is of the form
```
    imgout = bandpass(√max(0,img-bias))
```
i.e., the image is bias-subtracted, square-root transformed (to turn
shot noise into constant variance), and then band-pass filtered using
Gaussian filters of width `sigmalp` (for the low-pass) and `sigmahp`
(for the high-pass).  You can pass `sigmalp=zeros(n)` to skip low-pass
filtering, and `sigmahp=fill(Inf, n)` to skip high-pass filtering.
"""
mutable struct PreprocessSNF  # Shot-noise filtered
    bias::Float32
    sigmalp::Vector{Float32}
    sigmahp::Vector{Float32}
end
# PreprocessSNF(bias::T, sigmalp, sigmahp) = PreprocessSNF{T}(bias, T[sigmalp...], T[sigmahp...])

function preprocess(pp::PreprocessSNF, A::AbstractArray)
    Af = sqrt_subtract_bias(A, pp.bias)
    imfilter(highpass(Af, pp.sigmahp), KernelFactors.IIRGaussian((pp.sigmalp...,)), NA())
end
(pp::PreprocessSNF)(A::AbstractArray) = preprocess(pp, A)
# For SubArrays, extend to the parent along any non-sliced
# dimension. That way, we keep any information from padding.
function (pp::PreprocessSNF)(A::SubArray)
    Bpad = preprocess(pp, paddedview(A))
    trimmedview(Bpad, A)
end
# ImageMeta method is defined under @require in __init__

function sqrt_subtract_bias(A, bias)
#    T = typeof(sqrt(one(promote_type(eltype(A), typeof(bias)))))
    T = Float32
    out = Array{T}(undef, size(A))
    for I in eachindex(A)
        @inbounds out[I] = sqrt(max(zero(T), convert(T, A[I]) - bias))
    end
    out
end

"""
`Apad = paddedview(A)`, for a SubArray `A`, returns a SubArray that
extends to the full parent along any non-sliced dimensions of the
parent. See also [`trimmedview`](@ref).
"""
paddedview(A::SubArray) = _paddedview(A, (), (), A.indices...)
_paddedview(A::SubArray{T,N,P,I}, newindexes, newsize) where {T,N,P,I} =
    SubArray(A.parent, newindexes)
@inline function _paddedview(A, newindexes, newsize, index, indexes...)
    d = length(newindexes)+1
    _paddedview(A, (newindexes..., pdindex(A.parent, d, index)), pdsize(A.parent, newsize, d, index), indexes...)
end
pdindex(A, d, i::Base.Slice) = i
pdindex(A, d, i::Real) = i
pdindex(A, d, i::UnitRange) = 1:size(A,d)
pdindex(A, d, i) = error("Cannot pad with an index of type ", typeof(i))

pdsize(A, newsize, d, i::Base.Slice) = tuple(newsize..., size(A,d))
pdsize(A, newsize, d, i::Real) = newsize
pdsize(A, newsize, d, i::UnitRange) = tuple(newsize..., size(A,d))

"""
`B = trimmedview(Bpad, A::SubArray)` returns a SubArray `B` with
`axes(B) = axes(A)`. `Bpad` must have the same size as `paddedview(A)`.
"""
function trimmedview(Bpad, A::SubArray)
    ndims(Bpad) == ndims(A) || throw(DimensionMismatch("dimensions $(ndims(Bpad)) and $(ndims(A)) of Bpad and A must match"))
    _trimmedview(Bpad, A.parent, 1, (), A.indices...)
end
_trimmedview(Bpad, P, d, newindexes) = view(Bpad, newindexes...)
@inline _trimmedview(Bpad, P, d, newindexes, index::Real, indexes...) =
    _trimmedview(Bpad, P, d+1, newindexes, indexes...)
@inline function _trimmedview(Bpad, P, d, newindexes, index, indexes...)
    dB = length(newindexes)+1
    Bsz = size(Bpad, dB)
    Psz = size(P, d)
    Bsz == Psz || throw(DimensionMismatch("dimension $dB of Bpad has size $Bsz, should have size $Psz"))
    _trimmedview(Bpad, P, d+1, (newindexes..., index), indexes...)
end

# For faster and type-stable slicing
struct ColonFun end
ColonFun(::Int) = Colon()

function __init__()
    @require ImageMetadata="bc367c6b-8a6b-528e-b4bd-a4b897500b49" begin
        (pp::PreprocessSNF)(A::ImageMetadata.ImageMeta) = ImageMetadata.shareproperties(A, pp(ImageMetadata.arraydata(A)))
    end
end

include("deprecations.jl")

end


text:
int_rules_5_1_1 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 5.1.1*(a+b*arcsin(c*x))^n =#
    @apply_utils Antiderivative((~(a') + ~(b') * asin(~(c') * ~x)) ^ ~(n'), ~x) => ~x * (~a + ~b * asin(~c * ~x)) ^ ~n - ~b * ~c * ~n * Antiderivative((~x * (~a + ~b * asin(~c * ~x)) ^ (~n - 1)) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && GtQ(~n, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * acos(~(c') * ~x)) ^ ~(n'), ~x) => ~x * (~a + ~b * acos(~c * ~x)) ^ ~n + ~b * ~c * ~n * Antiderivative((~x * (~a + ~b * acos(~c * ~x)) ^ (~n - 1)) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && GtQ(~n, 0)
    @apply_utils Antiderivative((~(a') + ~(b') * asin(~(c') * ~x)) ^ ~n, ~x) => (sqrt(1 - (~c) ^ 2 * (~x) ^ 2) * (~a + ~b * asin(~c * ~x)) ^ (~n + 1)) / (~b * ~c * (~n + 1)) + (~c / (~b * (~n + 1))) * Antiderivative((~x * (~a + ~b * asin(~c * ~x)) ^ (~n + 1)) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && LtQ(~n, -1)
    @apply_utils Antiderivative((~(a') + ~(b') * acos(~(c') * ~x)) ^ ~n, ~x) => (-(sqrt(1 - (~c) ^ 2 * (~x) ^ 2)) * (~a + ~b * acos(~c * ~x)) ^ (~n + 1)) / (~b * ~c * (~n + 1)) - (~c / (~b * (~n + 1))) * Antiderivative((~x * (~a + ~b * acos(~c * ~x)) ^ (~n + 1)) / sqrt(1 - (~c) ^ 2 * (~x) ^ 2), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && LtQ(~n, -1)
    @apply_utils Antiderivative((~(a') + ~(b') * asin(~(c') * ~x)) ^ ~n, ~x) => (1 / (~b * ~c)) * Subst(Antiderivative((~x) ^ ~n * cos(-(~a) / ~b + ~x / ~b), ~x), ~x, ~a + ~b * asin(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
    @apply_utils Antiderivative((~(a') + ~(b') * acos(~(c') * ~x)) ^ ~n, ~x) => (-1 / (~b * ~c)) * Subst(Antiderivative((~x) ^ ~n * sin(-(~a) / ~b + ~x / ~b), ~x), ~x, ~a + ~b * acos(~c * ~x)) <-- FreeQ([~a, ~b, ~c, ~n], ~x)
end

max_stars_repo_path: src/rules/5 Inverse trig functions/5.1 Inverse sine/5.1.1 (a+b arcsin(c x))^n.jl


text:
import Base.Random: rand

export propsel, propsel!

@doc """propsel(p::Population)

Create a new population through proportional selection on the given
population. Employs the algorithm describe in the Lipowski and
Lipowska paper below.

References:

https://en.wikipedia.org/wiki/Fitness_proportionate_selection

A. Lipowski and D. Lipowska, “Roulette-wheel selection via stochastic
acceptance,” Phys. A Stat. Mech. its Appl., vol. 391, no. 6, pp. 2193–2196,
2012.
"""
function propsel(p::Population)
  np = Population(p)
  propsel!(np)
  return np
end

function propsel(mp::MetaPopulation)
  np = MetaPopulation(mp)
  propsel!(np)
  return np
end

@doc """propsel!(p::Population)

Conduct proportional selection in-place.
"""
function propsel!(p::Population)
  fs = popfits(p)
  fmax = maximum(fs)
  if fmax == 0
    # all elements have fitness zero
    return
  end

  n = popsize(p)
  selected = zeros(Int64, n)
  k = 0
  while k < n
    i = rand(1:n)
    w = fs[i] / fmax
    if rand() < w
      selected[k + 1] = i
      k += 1
    end
  end

  p.genotypes[:] = [Genotype(p.genotypes[selected[i]]) for i = 1:n]
end

function propsel!(mp::MetaPopulation)
  for p = mp.populations
    propsel!(p)
  end
end


text:
function  grids(nr,nx,rbottom)

#   constructs an equi-spaced grid in radius and x=cos(theta)
#   This particular version uses the shell thickness as a length
#   scale. This means that the radius at the top boundary is not r = 1;
#
#   nr       -  number of grid points in radius
#   nx       -  number of grid point in cos(theta)
#   rbottom   -  radius of lower layer  (dimensionless)
#
#   output
#   r        - radial grid (excludes boundaries)
#   x        - x grid (excludes poles and equator)
#   theta    - colatitude (radian)

   rcmb = 1.0;
   dr = (rcmb-rbottom) / nr;
   r = (rbottom+dr/2:dr:rcmb-dr/2);     # exclude boundaries

#  avoid grid point on equator
   if (rem(nx,2) ==1)
      dx = 2/(nx-1);
   else
       dx = 2/nx;
   end

   x = (-1 + dx/2 : dx : 1 - dx/2);
   theta = acos.(x);   # colatitude (in radian)

   return r, x, theta

end

max_stars_repo_name: bbuffett/MacWaves

max_issues_repo_name: bbuffett/MacWaves

max_forks_repo_name: bbuffett/MacWaves


text:
module Classification
export Optimizer, OptimizerContext, GDContext, GradientDescent,
       StochasticGradientDescent, SGDContext, LogisticClassifier,
       fit!, add_designmatrix!, designmatrix, predict, σ, logcrossentropy,
       mse, predict, sigmoid!, logcrossentropy!, NesterovGradientDescent,
       NAGDContext, partition, addvalidationset!, stopearly!, postprocess!,
       score, accuracy, NeuralNet, weave!, addLayer!, addSigmoidLayer!,
       sigmoidactivation!, FullyConnected, forwardpropagation, backpropagation,
       feedforward, minibatch, traintestsplit, onehot, onecold, evaluate,
       addSoftmaxLayer!

include("optimizer.jl")
using .Optimizers
include("gd.jl")
import .GD: GDContext, GradientDescent, fit!
include("sgd.jl")
import .SGD: SGDContext, StochasticGradientDescent, fit!
using .SGD: partition
include("nag.jl")
import .NAGD: NAGDContext, NesterovGradientDescent, fit!
include("classifier.jl")
using .Classifiers: LogisticClassifier, designmatrix, add_designmatrix!,  predict, score, traintestsplit
import .Classifiers: fit!
include("neuralnet.jl")
using .NeuralNetwork

include("gd_net.jl")
include("sgd_net.jl")
include("onehot.jl")
end


text:
module DIIS
using LinearAlgebra
import Base.push!
import Base.length

"""
    DIISManager{T1<:AbstractFloat,
                T2<:AbstractFloat}

# Fields

    vecs::Array{Array{T1},1} vectors to extrapolate from.
    errs::Array{Array{T2},1} error vectors to build B matrix from.
    max_vec::Int64           max number of vectors to hold
"""
struct DIISManager{T1<:AbstractFloat,
                  T2 <: AbstractFloat }# where T <: AbstractFloat
    vecs::Array{Array{T1},1}
    errs::Array{Array{T2},1}
    max_vec::Int64
end

struct CROPManager{T<:AbstractFloat}
    Wopt::Array{Array{T}}
    Topt::Array{Array{T}}
    Waux::Array{Array{T}}
    Taux::Array{Array{T}}
end

function DIISManager{T1,T2}(;size=6) where { T1 <: AbstractFloat,
                                             T2 <: AbstractFloat }
    vecs = Array{Array{T1}}(undef,0)
    errs = Array{Array{T2}}(undef,0)
    DIISManager{T1,T2}(vecs,errs,size)
end

function CROPManager{T}(;size=6) where T
    Wopt = Array{Array{T}}(undef,0)
    Topt = Array{Array{T}}(undef,0)
    Waux = Array{Array{T}}(undef,0)
    Taux = Array{Array{T}}(undef,0)
    CROPManager{T}(Wopt,Topt,Waux,Taux)
end

function extrap(M::CROPManager{T}) where T
end

function length(M::DIISManager{T1,T2}) where { T1 <: AbstractFloat,
                                               T2 <: AbstractFloat }
    length(M.vecs)
end

function push!(M::DIISManager{T1,T2}, V::Array, E::Array) where { T1 <: AbstractFloat,
                                                                  T2 <: AbstractFloat }
    if length(M)+1 > M.max_vec
        norms = norm.(M.errs)
        idx = findmax(norms)[2]
        deleteat!(M.vecs,idx)
        deleteat!(M.errs,idx)
    end
    push!(M.vecs,convert(Array{T1},deepcopy(V)))
    push!(M.errs,convert(Array{T2},deepcopy(E)))
end

"""
    extrapolate(M::DIISManager{T1,T2}; add_res=false) where { T1 <: AbstractFloat,
                                                              T2 <: AbstractFloat }

Takes current state of M and produces an optimal (within subspace) trial vector.

# kwargs
    add_res=false      Do add estimated residual to trial vector?
"""
function extrapolate(M::DIISManager{T1,T2};add_res=false) where { T1 <: AbstractFloat,
                                                    T2 <: AbstractFloat }
    diis_size = length(M)
    B = ones(T1,diis_size+1,diis_size+1)*1
    B[end,end] = 0
    for (n1, e1) in enumerate(M.errs[1:end])
        for (n2, e2) in enumerate(M.errs[1:end])
            B[n1,n2] = sum(e1 .* e2) #(e1,e2)
        end
    end
    E = size(B,1)
    resid = zeros(T1,diis_size+1)
    resid[end] = 1
    ci = svd(B)\resid
    out = zeros(T1,size(M.vecs[1]))
    add_res ? outw = zeros(T1,size(M.vecs[1])) : nothing
    for num in 1:diis_size
        out += ci[num]*M.vecs[num]
        add_res ? outw += ci[num]*M.errs[num] : nothing
    end
    add_res ? out += outw : nothing
    out
end
end #module


text:
module GPIO

    """
    Mode of the GPIO
    """
    module Mode

        @enum ModeCode::UInt8 begin
            PULL_UP_INPUT_MODE = 0
            INPUT_BINDING_MODE = 1
            PUSH_PULL_OUTPUT_MODE = 2
            PWM_OUTPUT_MODE = 3
            OPEN_DRAIN_OUTPUT_MODE = 4
        end

    end

end

"""
    NexGPIO(nexSerial, name; pid=pid, cid=cid)

A Nextion control to General Purpose Input/Output (GPIO).

!!! note

    Only for enhanced model.
"""
struct NexGPIO <: AbstractNexObject
    _nid::NexID

    function NexGPIO(nexSerial::T, name::Name; pid=PageID(), cid=ComponentID()) where {T <: AbstractNexSerial}
        nid = NexID(nexSerial, name, pid, cid)
        new(nid)
    end
end


text:
"""
    replace_missing(a::AbstractGeoArray, newmissingval)
    replace_missing(a::AbstractGeoStack, newmissingval)

Replace missing values in the array or stack with a new missing value,
also updating the `missingval` field/s.

A `GeoArray` containing a newly allocated `Array` is always returned,
even when the missing value matches the current value.
"""
replace_missing(a::DiskGeoArray, args...) = replace_missing(GeoArray(a), args...)
function replace_missing(a::MemGeoArray, newmissingval=missing)
    newdata = if ismissing(missingval(a))
        if newmissingval === missing
            copy(parent(a))
        else
            collect(Missings.replace(parent(a), newmissingval))
        end
    else
        replace(parent(a), missingval(a) => newmissingval)
    end
    rebuild(a; data=newdata, missingval=newmissingval)
end
function replace_missing(stack::AbstractGeoStack, newmissingval=missing)
    rebuild(stack, map(a -> replace_missing(a, newmissingval, values(stack))))
end

"""
    boolmask(A::AbstractArray, [missingval])

Create a mask array of `Bool` values, from any AbstractArray. For `AbstractGeoArray`
the default `missingval` is `missingval(A)`, for all other `AbstractArray`s
it is `missing`.

The array returned from calling `boolmask` on a `AbstractGeoArray` is a
[`GeoArray`](@ref) with the same size and fields as the oridingl array
"""
function boolmask end
boolmask(A::AbstractGeoArray) =
    rebuild(A; data=boolmask(A, missingval(A)), missingval=false, name=:Bool_mask)
boolmask(A::AbstractArray, missingval::Missing=missing) = (a -> !ismissing(a)).(parent(A))
boolmask(A::AbstractArray, missingval) =
    if isnan(missingval)
        (a -> !isnan(a)).(parent(A))
    else
        (a -> !isapprox(a, missingval)).(parent(A))
    end

"""
    missingmask(A::AbstractArray, [missingval])

Create a mask array of `missing` or `true` values, from any AbstractArray.
For `AbstractGeoArray` the default `missingval` is `missingval(A)`,
for all other `AbstractArray`s it is `missing`.

The array returned from calling `boolmask` on a `AbstractGeoArray` is a
[`GeoArray`](@ref) with the same size and fields as the oridingl array
"""
function missingmask end
missingmask(A::AbstractGeoArray) =
    rebuild(A; data=missingmask(A, missingval(A)), missingval=missing, name=:Missing_mask)
missingmask(A::AbstractArray, missingval::Missing=missing) =
    (a -> ismissing(a) ? missing : true).(parent(A))
missingmask(A::AbstractArray, missingval) =
    if isnan(missingval)
        (a -> isnan(a) ? missing : true).(parent(A))
    else
        (a -> isapprox(a, missingval) ? missing : true).(parent(A))
    end


text:
using OrdinaryDiffEq, Plots
gr()

#Half-life of Carbon-14 is 5,730 years.
C₁ = 5.730

#Setup
u₀ = 1.0
tspan = (0.0, 1.0)

#Define the problem
radioactivedecay(u,p,t) = -C₁*u

#Pass to solver
prob = ODEProblem(radioactivedecay,u₀,tspan)
sol = solve(prob,Tsit5())

#Plot
plot(sol,linewidth=2,title ="Carbon-14 half-life", xaxis = "Time in thousands of years", yaxis = "Percentage left", label = "Numerical Solution")
plot!(sol.t, t->exp(-C₁*t),lw=3,ls=:dash,label="Analytical Solution")

# Simple Harmonic Oscillator Problem
using OrdinaryDiffEq, Plots

#Parameters
ω = 1

#Initial Conditions
x₀ = [0.0]
dx₀ = [π/2]
tspan = (0.0, 2π)

ϕ = atan((dx₀[1]/ω)/x₀[1])
A = √(x₀[1]^2 + dx₀[1]^2)

#Define the problem
function harmonicoscillator(ddu,du,u,ω,t)
    ddu .= -ω^2 * u
end

#Pass to solvers
prob = SecondOrderODEProblem(harmonicoscillator, dx₀, x₀, tspan, ω)
sol = solve(prob, DPRKN6())

#Plot
plot(sol, vars=[2,1], linewidth=2, title ="Simple Harmonic Oscillator", xaxis = "Time", yaxis = "Elongation", label = ["x" "dx"])
plot!(t->A*cos(ω*t-ϕ), lw=3, ls=:dash, label="Analytical Solution x")
plot!(t->-A*ω*sin(ω*t-ϕ), lw=3, ls=:dash, label="Analytical Solution dx")

# Simple Pendulum Problem
using OrdinaryDiffEq, Plots

#Constants
const g = 9.81
L = 1.0

#Initial Conditions
u₀ = [0,π/2]
tspan = (0.0,6.3)

#Define the problem
function simplependulum(du,u,p,t)
    θ = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -(g/L)*sin(θ)
end

#Pass to solvers
prob = ODEProblem(simplependulum, u₀, tspan)
sol = solve(prob,Tsit5())

#Plot
plot(sol,linewidth=2,title ="Simple Pendulum Problem", xaxis = "Time", yaxis = "Height", label = ["\\theta" "d\\theta"])

p = plot(sol,vars = (1,2), xlims = (-9,9), title = "Phase Space Plot", xaxis = "Velocity", yaxis = "Position", leg=false)
function phase_plot(prob, u0, p, tspan=2pi)
    _prob = ODEProblem(prob.f,u0,(0.0,tspan))
    sol = solve(_prob,Vern9()) # Use Vern9 solver for higher accuracy
    plot!(p,sol,vars = (1,2), xlims = nothing, ylims = nothing)
end
for i in -4pi:pi/2:4π
    for j in -4pi:pi/2:4π
        phase_plot(prob, [j,i], p)
    end
end
plot(p,xlims = (-9,9))

#Double Pendulum Problem
using OrdinaryDiffEq, Plots

#Constants and setup
const m₁, m₂, L₁, L₂ = 1, 2, 1, 2
initial = [0, π/3, 0, 3pi/5]
tspan = (0.,50.)

#Convenience function for transforming from polar to Cartesian coordinates
function polar2cart(sol;dt=0.02,l1=L₁,l2=L₂,vars=(2,4))
    u = sol.t[1]:dt:sol.t[end]

    p1 = l1*map(x->x[vars[1]], sol.(u))
    p2 = l2*map(y->y[vars[2]], sol.(u))

    x1 = l1*sin.(p1)
    y1 = l1*-cos.(p1)
    (u, (x1 + l2*sin.(p2),
     y1 - l2*cos.(p2)))
end

#Define the Problem
function double_pendulum(xdot,x,p,t)
    xdot[1]=x[2]
    xdot[2]=-((g*(2*m₁+m₂)*sin(x[1])+m₂*(g*sin(x[1]-2*x[3])+2*(L₂*x[4]^2+L₁*x[2]^2*cos(x[1]-x[3]))*sin(x[1]-x[3])))/(2*L₁*(m₁+m₂-m₂*cos(x[1]-x[3])^2)))
    xdot[3]=x[4]
    xdot[4]=(((m₁+m₂)*(L₁*x[2]^2+g*cos(x[1]))+L₂*m₂*x[4]^2*cos(x[1]-x[3]))*sin(x[1]-x[3]))/(L₂*(m₁+m₂-m₂*cos(x[1]-x[3])^2))
end

#Pass to Solvers
double_pendulum_problem = ODEProblem(double_pendulum, initial, tspan)
sol = solve(double_pendulum_problem, Vern7(), abs_tol=1e-10, dt=0.05);

#Obtain coordinates in Cartesian Geometry
ts, ps = polar2cart(sol, l1=L₁, l2=L₂, dt=0.01)
plot(ps...)

#Constants and setup
using OrdinaryDiffEq
initial2 = [0.01, 0.005, 0.01, 0.01]
tspan2 = (0.,500.)

#Define the problem
function double_pendulum_hamiltonian(udot,u,p,t)
    α  = u[1]
    lα = u[2]
    β  = u[3]
    lβ = u[4]
    udot .=
    [2(lα-(1+cos(β))lβ)/(3-cos(2β)),
    -2sin(α) - sin(α+β),
    2(-(1+cos(β))lα + (3+2cos(β))lβ)/(3-cos(2β)),
    -sin(α+β) - 2sin(β)*(((lα-lβ)lβ)/(3-cos(2β))) + 2sin(2β)*((lα^2 - 2(1+cos(β))lα*lβ + (3+2cos(β))lβ^2)/(3-cos(2β))^2)]
end

# Construct a ContiunousCallback
condition(u,t,integrator) = u[1]
affect!(integrator) = nothing
cb = ContinuousCallback(condition,affect!,nothing,
                        save_positions = (true,false))

# Construct Problem
poincare = ODEProblem(double_pendulum_hamiltonian, initial2, tspan2)
sol2 = solve(poincare, Vern9(), save_everystep = false, save_start=false, save_end=false, callback=cb, abstol=1e-16, reltol=1e-16,)

function poincare_map(prob, u₀, p; callback=cb)
    _prob = ODEProblem(prob.f, u₀, prob.tspan)
    sol = solve(_prob, Vern9(), save_everystep = false, save_start=false, save_end=false, callback=cb, abstol=1e-16, reltol=1e-16)
    scatter!(p, sol, vars=(3,4), markersize = 3, msw=0)
end

lβrange = -0.02:0.0025:0.02
p = scatter(sol2, vars=(3,4), leg=false, markersize = 3, msw=0)
for lβ in lβrange
    poincare_map(poincare, [0.01, 0.01, 0.01, lβ], p)
end
plot(p, xlabel="\\beta", ylabel="l_\\beta", ylims=(0, 0.03))

using OrdinaryDiffEq, Plots

#Setup
initial = [0.,0.1,0.5,0]
tspan = (0,100.)

#Remember, V is the potential of the system and T is the Total Kinetic Energy, thus E will
#the total energy of the system.
V(x,y) = 1//2 * (x^2 + y^2 + 2x^2*y - 2//3 * y^3)
E(x,y,dx,dy) = V(x,y) + 1//2 * (dx^2 + dy^2);

#Define the function
function Hénon_Heiles(du,u,p,t)
    x  = u[1]
    y  = u[2]
    dx = u[3]
    dy = u[4]
    du[1] = dx
    du[2] = dy
    du[3] = -x - 2x*y
    du[4] = y^2 - y -x^2
end

#Pass to solvers
prob = ODEProblem(Hénon_Heiles, initial, tspan)
sol = solve(prob, Vern9(), abs_tol=1e-16, rel_tol=1e-16);

# Plot the orbit
plot(sol, vars=(1,2), title = "The orbit of the Hénon-Heiles system", xaxis = "x", yaxis = "y", leg=false)

#Optional Sanity check - what do you think this returns and why?
@show sol.retcode

#Plot -
plot(sol, vars=(1,3), title = "Phase space for the Hénon-Heiles system", xaxis = "Position", yaxis = "Velocity")
plot!(sol, vars=(2,4), leg = false)

#We map the Total energies during the time intervals of the solution (sol.u here) to a new vector
#pass it to the plotter a bit more conveniently
energy = map(x->E(x...), sol.u)

#We use @show here to easily spot erratic behaviour in our system by seeing if the loss in energy was too great.
@show ΔE = energy[1]-energy[end]

#Plot
plot(sol.t, energy .- energy[1], title = "Change in Energy over Time", xaxis = "Time in iterations", yaxis = "Change in Energy")

function HH_acceleration!(dv,v,u,p,t)
    x,y  = u
    dx,dy = dv
    dv[1] = -x - 2x*y
    dv[2] = y^2 - y -x^2
end
initial_positions = [0.0,0.1]
initial_velocities = [0.5,0.0]
prob = SecondOrderODEProblem(HH_acceleration!,initial_velocities,initial_positions,tspan)
sol2 = solve(prob, KahanLi8(), dt=1/10);

# Plot the orbit
plot(sol2, vars=(3,4), title = "The orbit of the Hénon-Heiles system", xaxis = "x", yaxis = "y", leg=false)

plot(sol2, vars=(3,1), title = "Phase space for the Hénon-Heiles system", xaxis = "Position", yaxis = "Velocity")
plot!(sol2, vars=(4,2), leg = false)

energy = map(x->E(x[3], x[4], x[1], x[2]), sol2.u)
#We use @show here to easily spot erratic behaviour in our system by seeing if the loss in energy was too great.
@show ΔE = energy[1]-energy[end]

#Plot
plot(sol2.t, energy .- energy[1], title = "Change in Energy over Time", xaxis = "Time in iterations", yaxis = "Change in Energy")

sol3 = solve(prob, DPRKN6());
energy = map(x->E(x[3], x[4], x[1], x[2]), sol3.u)
@show ΔE = energy[1]-energy[end]
gr()
plot(sol3.t, energy .- energy[1], title = "Change in Energy over Time", xaxis = "Time in iterations", yaxis = "Change in Energy")

using DiffEqTutorials
DiffEqTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
using UnPack
using Test

@testset "Dict and NamedTuple" begin
    ###########################
    # Packing and unpacking @unpack, @pack!
    ##########################
    # Example with dict:
    d = Dict{Symbol,Any}(:a=>5.0,:b=>2,:c=>"Hi!")
    @unpack a, c = d
    @test a == 5.0 #true
    @test c == "Hi!" #true

    d = Dict("a"=>5.0,"b"=>2,"c"=>"Hi!")
    @unpack a, c = d
    @test a == 5.0 #true
    @test c == "Hi!" #true

    # Example with named tuple
    @eval d = (a=5.0, b=2, c="Hi!")
    @unpack a, c = d
    @test a == 5.0 #true
    @test c == "Hi!" #true
end

# having struct-defs inside a testset seems to be problematic in some julia version
mutable struct PropertyExample
    a
    last_set_property
end
@testset "getproperty" begin
    Base.getproperty(::PropertyExample, name::Symbol) = String(name)
    Base.setproperty!(d::PropertyExample, name::Symbol, value) =
        setfield!(d, :last_set_property, (name, value))
    Base.propertynames(::PropertyExample) = (:A, :B, :C)

    let d = PropertyExample(:should_be_ignored, nothing)
        @unpack a, b = d
        @test a == "a"
        @test b == "b"

        a = "a value"
        @pack! d = a
        @test getfield(d, :last_set_property) == (:a, "a value")
    end

    # TODO add test with non String string
end

mutable struct A; a; b; c; end
@testset "type" begin
    # Example with type:

    d = A(4,7.0,"Hi!")
    @unpack a, c = d
    @test a == 4 #true
    @test c == "Hi!" #true
end

@testset "Packing" begin
    # Example with dict:
    a = 5.0
    c = "Hi!"
    d = Dict{Symbol,Any}()
    @pack! d = a, c
    @test d==Dict{Symbol,Any}(:a=>5.0,:c=>"Hi!")

    d = Dict{String,Any}()
    @pack! d = a, c
    @test d==Dict{String,Any}("a"=>5.0,"c"=>"Hi!")

    # Example with type:
    a = 99
    c = "HaHa"
    d = A(4,7.0,"Hi")
    @pack! d = a, c
    @test d.a == 99
    @test d.c == "HaHa"
end

mutable struct UP1
    aUP1
    bUP1
end
struct UP2
    aUP2
    bUP2
end

@testset "old tests" begin
    uu = UP1(1,2)
    @test_throws ErrorException @unpack cUP1 = uu
    @test_throws ErrorException @unpack aUP1, cUP1 = uu

    aUP1, bUP1 = 0, 0
    @unpack aUP1 = uu
    @test aUP1==1
    @test bUP1==0
    aUP1, bUP1 = 0, 0
    @unpack aUP1, bUP1 = uu
    @test aUP1==1
    @test bUP1==2

    vv = uu
    aUP1 = 99
    @pack! uu = aUP1
    @test uu==vv
    @test uu.aUP1==99

    uu = UP2(1,2)
    @test_throws ErrorException @unpack cUP2 = uu
    @test_throws ErrorException @unpack aUP2, cUP2 = uu

    aUP1 = 99
    @test_throws ErrorException @pack! uu = aUP1
end

struct UP3
    a::Float64
    b::Int
end

@testset "inference" begin
    function f(u::UP3)
        @unpack a,b = u
        a,b
    end
    @inferred f(UP3(1,2))
end


text:
function hcompress(A::HMat)
    if A.blockType == LOWRANK
        if size(A.UMat,2)>A.MAXRANK
            (QU,RU) = qr(A.UMat);
            (QV,RV) = qr(A.VMat);
            (Utmp,Stmp,Vtmp) = svdtrunc(RU*RV',A.EPS,A.MAXRANK);
            if length(Stmp) > 0
                A.UMat = QU*(Utmp.*sqrt(Stmp)');
                A.VMat = QV*(Vtmp.*sqrt(Stmp)');
            else
                A.UMat = QU*Utmp;
                A.VMat = QV*Vtmp;
            end
        end
    elseif A.blockType == HMAT
        for C in A.childHMat
            hcompress(C);
        end
    end
end

function hh2l(A::HMat)
    Rcol = randn(A.width,A.MAXRANK+5);
    Rrow = randn(A.height,A.MAXRANK+5);
    ARcol = hmatvec(A,Rcol);
    ATRrow = hmatTvec(A,Rrow);
    (Qcol,) = qr(ARcol);
    (Qrow,) = qr(ATRrow);
    (Utmp,Stmp,Vtmp) = svdtrunc(pinv(Rrow'*Qcol)*Rrow'*ARcol*pinv(Qrow'*Rcol),A.EPS,A.MAXRANK);
    if length(Stmp) > 0
        UMat = Qcol*(Utmp.*sqrt(Stmp)');
        VMat = Qrow*(Vtmp.*sqrt(Stmp)');
    else
        UMat = Qcol*Utmp;
        VMat = Qrow*Vtmp;
    end
    return UMat,VMat;
end

max_stars_repo_name: YingzhouLi/HMat

max_issues_repo_name: YingzhouLi/HMat


text:
module fresnelltools

using LinearAlgebra

using ..matrixcore

export FresnellBoundrary, FresnellSlab, Grating, ContinousBorder, ThreeLayerSystem, Interrogator, BorderInterrogator

function FresnellBoundrary(n1::Number, n2::Number)::Matrix{Number}
    # Equation from Saleh & Teich 3.ed. eq.6.2-8 & eq. 6.2-9
    # Returns the scattering matrix for the relevant border

    FresnellBoundrary(n1, n2, 0)[1]
end

function FresnellBoundrary(n1::Number, n2::Number, θ1::Number)::Tuple{Matrix{Number}, Matrix{Number}, Number}
    # Equation from Saleh & Teich 3.ed. eq.6.2-8 & eq. 6.2-9
    # It is similar to the other FresnellBoundrary function, but takes incident angle into account
    # Returns the TE scattering matrix, the TM scattering matrix and the outgoing transmitted angle

    # The imaginary coefficient 1e-14im has no physical reason to exist, and is only added for stability in the program.
    # It is assumed to be needed because of type instability in the asin function i Julia.
    # See https://github.com/JuliaLang/julia/issues/24296
    if ((typeof(n1) <: Complex) & (typeof(n2) <: Real))
        complex_modifier = 1e-14im
    else
        complex_modifier = 0im
    end
    θ2 = asin(complex_modifier + sin(θ1) * n1 / n2)

    r_te_12 = (n1 * cos(θ1) - n2 * cos(θ2)) / (n1 * cos(θ1) + n2 * cos(θ2))
    r_te_21 = (n2 * cos(θ2) - n1 * cos(θ1)) / (n2 * cos(θ2) + n1 * cos(θ1))

    r_tm_12 = (n1 * sec(θ1) - n2 * sec(θ2)) / (n1 * sec(θ1) + n2 * sec(θ2))
    r_tm_21 = (n2 * sec(θ2) - n1 * sec(θ1)) / (n2 * sec(θ2) + n1 * sec(θ1))

    t_te_12 = 1 + r_te_12
    t_te_21 = 1 + r_te_21

    t_tm_12 = (1 + r_tm_12) * cos(θ1) / cos(θ2)
    t_tm_21 = (1 + r_tm_21) * cos(θ2) / cos(θ1)

    Ste = [
        t_te_12 r_te_21;
        r_te_12 t_te_21
    ]
    Stm = [
        t_tm_12 r_tm_21;
        r_tm_12 t_tm_21
    ]

    return Ste, Stm, θ2
end

function FresnellSlab(n::Number, k_0::Number, d::Number, θ::Number)::Matrix{Number}
    # Equation from Saleh & Teich 3.ed. eq.7.1-4
    delay = ℯ^(-1im * n * k_0 * d * cos(θ))
    S = [
        delay 0;
        0     delay
    ]

    if abs(S[1, 1]) < 1e-5
        S /= abs(S[1, 1])
        S *= 1e-5
    end

    S
end

function Grating(n_spechial::Number, n_normal::Number, d_spechial::Real, d_normal::Real, layers::Int, k_0::Number)::Matrix{Number}
    # A utility for creating the scattering matrix of a periodic grating

    border_normal_spechial = FresnellBoundrary(n_normal, n_spechial)
    bulk_spechial = FresnellSlab(n_spechial, k_0, d_spechial, 0)
    border_spechial_normal = FresnellBoundrary(n_spechial, n_normal)
    bulk_normal = FresnellSlab(n_normal, k_0, d_normal, 0)

    cell = CascadeScattering([border_normal_spechial, bulk_spechial, border_spechial_normal, bulk_normal])
    Repeating(cell, layers)
end

function Grating(n_spechial::Number, n_normal::Number, d_spechial::Real, d_normal::Real, layers::Int, k_0::Number, θ1::Real)::Tuple{Matrix{Number}, Matrix{Number}}
    # A utility similar to Grating, but it takes an angle and assumes an infinitely wide grating
    # Returns the TE and TM scattering matrices for the grating
    # Unlike the angleded FresnellBoundrary tool, no angle is returned, this is because the grating assumes normal material on both sides

    n_s_te, n_s_tm, θ2 = FresnellBoundrary(n_normal, n_spechial, θ1)
    s_n_te, s_n_tm, _  = FresnellBoundrary(n_spechial, n_normal, θ2)

    n_b = FresnellSlab(n_normal, k_0, d_normal, θ1)
    s_b = FresnellSlab(n_spechial, k_0, d_spechial, θ2)

    cell_te = CascadeScattering([n_s_te, s_b, s_n_te, n_b])
    cell_tm = CascadeScattering([n_s_tm, s_b, s_n_tm, n_b])

    Repeating(cell_te, layers), Repeating(cell_tm, layers)
end

function ContinousBorder(n_from::Number, n_to::Number, d::Real, k_0::Number, stepps::Int)::Matrix{Number}
    # A utility for creating a semi-continous change in refractive index
    # The gradient is approximated as linear

    components = Array{Matrix{Number}}(undef, 2 * stepps)

    n_stepp = (n_to - n_from) / stepps
    d_stepp = d / stepps
    for i in 1:stepps
        n1 = n_from + (i - 1) * n_stepp
        n2 = n1 + n_stepp

        components[2 * i - 1] = FresnellBoundrary(n1, n2)
        components[2 * i] = FresnellSlab(n2, k_0, d_stepp, 0)
    end

    CascadeScattering(components)
end

function ThreeLayerSystem(n_1::Number, n_bulk, n_3::Number, λ, θ, d)
    interface1_te, interface1_tm, θ2 = FresnellBoundrary(n_1, n_bulk(λ), θ)
    bulk = FresnellSlab(n_bulk(λ), 2*π / λ, d, θ2)
    interface2_te, interface2_tm, _ = FresnellBoundrary(n_bulk(λ), n_3, θ2)
    te_system = CascadeScattering([interface1_te, bulk, interface2_te])
    tm_system = CascadeScattering([interface1_tm, bulk, interface2_tm])

    te_system, tm_system
end

function Interrogator(layers::Vector{Function}, distances::Vector{Float64}, step::Float64)::Tuple{Vector{Float64}, Vector{Number}}
    # A small utility providing a default lambda for the Interrogator tool
    Power(Up, Um) = abs(Up + Um)^2
    Interrogator(layers, distances, step, Power)
end

function Interrogator(layers::Vector{Function}, distances::Vector{Float64}, step::Float64, expression::Function)::Tuple{Vector{Float64}, Vector{Number}}
    # A tool for inspecting the field inside a multilayer system

    total_system = Array{Matrix{Number}}(undef, length(layers))
    for i in 1:1:length(layers)
        total_system[i] = layers[i](distances[i])
    end
    U0_plus = 1
    U0_minus = (CascadeScattering(total_system) * [U0_plus, 0])[2]

    maxlen = ceil(Int, length(distances) * 2 + sum(distances) / step)
    d_values = Vector{Number}(undef, maxlen)
    u_values = Vector{Number}(undef, maxlen)

    d_values[1] = 0.0
    u_values[1] = expression(U0_plus, U0_minus)

    m_accumulated = I
    m_local = I
    m_partial = I
    d_accumulated = 0.0
    i = 1
    for (layer, distance) in zip(layers, distances)
        m_local = m_accumulated

        m_partial = StoM(layer(step))
        j = 0
        for _ in 0:step:distance
            i += 1
            j += 1
            d_accumulated += distance > step ? step : distance
            if (j % 1000) == 0
                m_local = StoM(layer(step * j)) * m_accumulated
            else
                m_local = m_partial * m_local
            end

            U_plus, U_minus = m_local * [U0_plus, U0_minus]
            measure = expression(U_plus, U_minus)
            d_values[i] = d_accumulated
            u_values[i] = measure
        end

        d_rest = distance % step
        if (d_rest > (step / 1000)) & (step > distance)
            m_partial = StoM(layer(d_rest))
            i += 1
            d_accumulated += d_rest
            m_local = m_partial * m_local

            (U_plus, U_minus) = m_local * [U0_plus, U0_minus]

            measure = expression(U_plus, U_minus)
            d_values[i] = d_accumulated
            u_values[i] = measure
        end

        m_accumulated = StoM(layer(distance)) * m_accumulated
    end

    u_values[1:i], d_values[1:i]
end

function BorderInterrogator(layers::Vector{Matrix{Number}}, distances::Vector{Float64}, expression::Function)::Tuple{Vector{Number}, Vector{Number}, Vector{Real}}
    # A tool for calculating the field at the borders in a structure independent of the Interrogator function
    # Returns the forward propagating fields, the backward propagating fields and the distances

    U0p = 1
    U0n = (CascadeScattering(layers) * [U0p, 0])[2]

    Upv = Vector{Number}()
    Unv = Vector{Number}()
    Dv = Vector{Float64}()

    push!(Upv, U0p)
    push!(Unv, U0n)
    push!(Dv, 0.0)

    accumulated_m = I
    accumulated_d = 0.0
    for (layer, distance) in zip(layers, distances)
        accumulated_m = StoM(layer) * accumulated_m
        accumulated_d += distance

        Up, Un = accumulated_m * [U0p, U0n]

        push!(Upv, Up)
        push!(Unv, Un)
        push!(Dv, accumulated_d)
    end

    return expression.(Upv), expression.(Unv), Dv
end

end # fresnelltools


text:
include("invHessianEstimationASym.jl")
include("invHessianEstimationSymReduced.jl")

########################
#
# util functions
#
########################
function _mul(S::Symmetric, θ::SparseIterate)
    p = size(S, 1)

    out = zeros(p, p)
    for ci=1:p
        for ri=1:p
            for k=1:p
                @inbounds out[ri, ci] += S[ri, k] * θ[(ci-1)*p + k]
            end
        end
    end
    out
end

function _mul(θ::SparseIterate, S::Symmetric)
    p = size(S, 1)

    out = zeros(p, p)
    for ci=1:p
        for ri=1:p
            for k=1:p
                @inbounds out[ri, ci] += θ[(k-1)*p + ri] * S[k, ci]
            end
        end
    end
    out
end


text:
#
# OIFITS.jl --
#
# Support for OI-FITS (optical interferometry data format) in Julia.
#
#------------------------------------------------------------------------------
#
# This file is part of OIFITS.jl which is licensed under the MIT "Expat"
# License:
#
# Copyright (C) 2015-2019: Éric Thiébaut.
#
#------------------------------------------------------------------------------

isdefined(Base, :__precompile__) && __precompile__(true)

module OIFITS

using Compat
using Compat: @debug, @error, @info, @warn

import Base: getindex, setindex!, haskey, keys, show

const Name = (isdefined(Core, :ASCIIString) ? ASCIIString : String)

include("oidata.jl")
include("fix-fitsio.jl")
include("misc.jl")
include("oifile.jl")
include("oiformat1.jl")
include("oiformat2.jl")
include("oipost.jl") # must be *after* oifile.jl and all oiformat*.jl
include("deprecations.jl")

end # module