text:
using LinearAlgebra: norm

#=
  Solvers for GLM
=#
abstract type AbstractSolver end

# Weights functions
@inline function W(distrib::AbstractDistribution, link::AbstractLink, mu::Array{T, 1}, eta::Array{T, 1}) where {T <: AbstractFloat}
  return ( (deta_dmu(link, mu, eta).^2) .* variance(distrib, mu)).^(-0.5)
end
function W(distrib::AbstractDistribution, link::AbstractLink, mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  return [( (deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-0.5) for i in 1:nBlocks]
end
function W(::Block1DParallel, distrib::AbstractDistribution, link::AbstractLink, mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  ret::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, nBlocks)
  @threads for i in 1:nBlocks
    ret[i] = ((deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-0.5)
  end
  return ret
end

#=
  GLM Solver uses (X'WX)^(-1) * X'Wy
=#
struct VanillaSolver <: AbstractSolver end

# Weights for the VanillaSolver
@inline function W(::VanillaSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{T, 1}, eta::Array{T, 1}) where {T <: AbstractFloat}
  return ((deta_dmu(link, mu, eta).^2) .* variance(distrib, mu)).^(-1)
end
function W(::VanillaSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  return [((deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-1) for i in 1:nBlocks]
end
function W(::Block1DParallel, ::VanillaSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  ret::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, nBlocks)
  @threads for i in 1:nBlocks
    ret[i] = ((deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-1)
  end
  return ret
end

#=
  Solver Using QR Decomposition
=#
struct QRSolver <: AbstractSolver end

# Weights for the QRSolver
@inline function W(::QRSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{T, 1}, eta::Array{T, 1}) where {T <: AbstractFloat}
  return ((deta_dmu(link, mu, eta).^2) .* variance(distrib, mu)).^(-0.5)
end
function W(::QRSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  return [((deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-0.5) for i in 1:nBlocks]
end
function W(::Block1DParallel, ::QRSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  ret::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, nBlocks)
  @threads for i in 1:nBlocks
    ret[i] = ((deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-0.5)
  end
  return ret
end

#=
  TODO:
  1. Unify the QRSolve and Vanilla solve interface as in
     the D implementation.
=#
function solve(::Type{<: VanillaSolver}, X::Array{T, 2}, y::Array{T, 1}) where {T <: AbstractFloat}
  ret::Array{T, 1} = Array{T, 1}(undef, 0)
  try
    ret = (X'*X)/(X'*y)
  catch
    ret = inv(X'*X)*(X'*y)
  end
  return ret
  #return inv(X'*X)*(X'*y)
end

import LinearAlgebra.LAPACK.geqrf!
import LinearAlgebra.LAPACK.orgqr!
import LinearAlgebra.BLAS.gemv
import LinearAlgebra.LAPACK.trtrs!

# import LinearAlgebra.UpperTriangular
# import LinearAlgebra.LAPACK.trtri!
# import LinearAlgebra.BLAS.trsm!

#=
  The QR Solver
=#
function getUpperR(QR::Array{T, 2}) where {T <: AbstractFloat}
  _, p = size(QR)
  R = zeros(T, (p, p))
  for i in 1:p
    for j in 1:i
      R[j, i] = QR[j, i]
    end
  end
  return R
end

function solve(::Type{<: QRSolver}, Q::Array{T, 2}, y::Array{T, 1}) where {T <: AbstractFloat}
  Q = copy(Q)
  _, p = size(Q)
  tau = zeros(T, p)
  geqrf!(Q, tau)
  R = getUpperR(Q)
  orgqr!(Q, tau, p)
  z = gemv('T', 1.0, Q, y)
  # z = Q' * y
  trtrs!('U', 'N', 'N', R, z)
  return z
end

#============================  Matrix Inverses ============================#
#= Matrix Inverse By LU Decomposition =#
import LinearAlgebra.LAPACK.getrf!
import LinearAlgebra.LAPACK.getri!

#=
  Inverse Types
=#
abstract type AbstractInverse end

#= GETRI Matrix Inverse =#
struct GETRIInverse <: AbstractInverse end
#= POTRI Matrix Inverse =#
struct POTRIInverse <: AbstractInverse end
#= SYTRF Matrix Inverse =#
struct SYTRFInverse <: AbstractInverse end
#= GESVD Matrix Inverse =#
struct GESVDInverse <: AbstractInverse end

#=
  Inverse of a square matrix using the GETRI function
  uses LU decomposition. The geqrf() does the QR decomposition
  and the getri() function carries out an inverse from the output
  of geqrf().
=#
function inv(::GETRIInverse, A::Array{T, 2}) where {T <: AbstractFloat}
  A = copy(A)
  p, _ = size(A)
  A, ipiv, info = getrf!(A)
  @assert(info == 0, "getrf! error in argument " * string(info))
  getri!(A, ipiv)
  return A
end

#= Matrix Inverse By Cholesky Decomposition For Positive Definite Matrices =#
import LinearAlgebra.LAPACK.potrf!
import LinearAlgebra.LAPACK.potri!

#=
  # Examples:
  x = rand(100, 10);
  x = x' * x
  sum(abs.(Base.inv(x) .- inv(POTRIInverse(), x)))
=#
function inv(::POTRIInverse, A::Array{T, 2}) where {T <: AbstractFloat}
  p, _ = size(A)
  A = copy(A)
  A, info = potrf!('U', A)
  @assert(info == 0, "potrf error in argument " * string(info))
  potri!('U', A)
  for i in 1:p
    for j in 1:i
      A[i, j] = A[j, i]
    end
  end
  return A
end

#= Matrix Inverse By LU Decomposition For Symmetric Indefinite Matrices =#
import LinearAlgebra.LAPACK.sytri!# (uplo, A, ipiv)
import LinearAlgebra.LAPACK.sytrf!# (uplo, A) -> (A, ipiv, info)

#=
  # Examples:
  x = rand(100, 10);
  x = x' * x
  sum(abs.(Base.inv(x) .- inv(SYTRFInverse(), x)))
=#
function inv(::SYTRFInverse, A::Array{T, 2}) where {T <: AbstractFloat}
  p, _ = size(A)
  A = copy(A)
  A, ipiv, info = sytrf!('U', A)
  sytri!('U', A, ipiv)
  for i in 1:p
    for j in 1:i
      A[i, j] = A[j, i]
    end
  end
  return A
end

#= Generalized Matrix Inverse By SVD  =#
import LinearAlgebra.LAPACK.gesvd! # (jobu, jobvt, A) -> (U, S, VT)

#= Matrix Vector Sweep Function =#
function multSweep(xw::Array{T, 2}, w::Array{T, 1})::Array{T, 2} where {T <: AbstractFloat}
  n, p = size(xw)
  for j in 1:p
    for i in 1:n
      xw[i, j] = xw[i, j] * w[i]
    end
  end
  return xw
end

#=
  # Examples:
  x = rand(100, 10);
  x = x' * x
  sum(abs.(Base.inv(x) .- inv(GESVDInverse(), x)))
=#
function inv(::GESVDInverse, A::Array{T, 2}) where {T <: AbstractFloat}
  p, _ = size(A)
  A = copy(A)
  U, S, Vt = gesvd!('A', 'A', A)
  S = [el > 1E-9 ? 1/el : el for el in S]
  return multSweep(Vt, S)' * U'
end

#============================== Linear Equation Solvers ==============================#
#= For Ax = b =#
#==============================  GESV Solver ==============================#

function calcXWX(xwx::Array{T, 2}, x::Array{T, 2}, w::Array{T, 1}, z::Array{T, 1}) where {T <: AbstractFloat}
  # n, p = size(x)
  xw = copy(x)
  xw = multSweep(xw, w)
  xwx = xw' * x
  xwz = xw' * z
  return xwx, xwz
end

function calcXWX(xwx::Array{T, 2}, x::Array{Array{T, 2}, 1}, w::Array{Array{T, 1}, 1}, z::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  p::Int64 = size(x[1])[2]
  nBlocks::Int64 = length(x)
  xwx = zeros(T, (p, p))
  xwz = zeros(T, (p, 1))
  for i in 1:nBlocks
    xw = copy(x[i])
    xw = multSweep(xw, w[i])
    xwx += xw' * x[i]
    xwz += xw' * z[i]
  end
  return xwx, xwz
end

function calcXWX(::Block1DParallel, xwx::Array{T, 2}, x::Array{Array{T, 2}, 1}, w::Array{Array{T, 1}, 1}, z::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  p::Int64 = size(x[1])[2]
  nBlocks::Int64 = length(x)
  xwx = [zeros(T, (p, p)) for i in 1:nthreads()]
  xwz = [zeros(T, (p, 1)) for i in 1:nthreads()]
  @threads for i in 1:nBlocks
    xw = copy(x[i])
    xw = multSweep(xw, w[i])
    xwx[threadid()] += xw' * x[i]
    xwz[threadid()] += xw' * z[i]
  end
  for i in 2:nthreads()
    xwx[1] .+= xwx[i]
    xwz[1] .+= xwz[i]
  end
  return xwx[1], xwz[1]
end

function calcXW(xwx::Array{T, 2}, x::Array{T, 2}, w::Array{T, 1}, z::Array{T, 1}) where {T <: AbstractFloat}
  # n, p = size(x)
  xw = copy(x)
  xw = multSweep(xw, w)
  return xw, z .* w
end

struct GESVSolver <: AbstractSolver end
# Weights for the GESVSolver
@inline function W(::GESVSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{T, 1}, eta::Array{T, 1}) where {T <: AbstractFloat}
  return ((deta_dmu(link, mu, eta).^2) .* variance(distrib, mu)).^(-1)
end
function W(::GESVSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  return [((deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-1) for i in 1:nBlocks]
end
function W(matrixType::Block1DParallel, ::GESVSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  ret::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, nBlocks)
  @threads for i in 1:nBlocks
    ret[i] = ((deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-1)
  end
  return ret
end

# Solver
import LinearAlgebra.LAPACK.gesv!
function solve!(::GESVSolver, R::Array{T, 2},
    xwx::Array{T, 2}, xw::Array{T, 2},
    x::Array{T, 2}, z::Array{T, 1}, w::Array{T, 1},
    coef::Array{T, 1}) where {T <: AbstractFloat}
  xwx, xwz = calcXWX(xwx, x, w, z)
  coef, _, ipiv = gesv!(copy(xwx), xwz)
  return (R, xwx, xw, x, z, w, coef)
end
function solve!(::GESVSolver, R::Array{T, 2},
  xwx::Array{T, 2}, xw::Array{Array{T, 2},1},
  x::Array{Array{T, 2}, 1}, z::Array{Array{T, 1}, 1},
  w::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}
  xwx, xwz = calcXWX(xwx, x, w, z)
  coef, _, ipiv = gesv!(copy(xwx), xwz)
  return (R, xwx, xw, x, z, w, coef)
end
function solve!(matrixType::Block1DParallel, ::GESVSolver, R::Array{T, 2},
  xwx::Array{T, 2}, xw::Array{Array{T, 2},1},
  x::Array{Array{T, 2}, 1}, z::Array{Array{T, 1}, 1},
  w::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}
  xwx, xwz = calcXWX(matrixType, xwx, x, w, z)
  coef, _, ipiv = gesv!(copy(xwx), xwz)
  return (R, xwx, xw, x, z, w, coef)
end
# Default Covariance Function
function cov(::AbstractSolver, invType::AbstractInverse, R::Array{T, 2}, xwx::Array{T, 2}, xw::Array{T, 2})::Array{T, 2} where {T <: AbstractFloat}
  return inv(invType, xwx)
end

function cov(::AbstractSolver, invType::AbstractInverse, R::Array{T, 2}, xwx::Array{T, 2})::Array{T, 2} where {T <: AbstractFloat}
  return inv(invType, xwx)
end

#==============================  POSV Solver ==============================#
#= Solver For Positive Definite Matrices =#
# Weights for posv!(uplo, A, B) -> (A, B)
import LinearAlgebra.LAPACK.posv! # (uplo, A, B) -> (A, B)

struct POSVSolver <: AbstractSolver end
# Weights for the POSVSolver
@inline function W(::POSVSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{T, 1}, eta::Array{T, 1}) where {T <: AbstractFloat}
  return ((deta_dmu(link, mu, eta).^2) .* variance(distrib, mu)).^(-1)
end
function W(::POSVSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  return [((deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-1) for i in 1:nBlocks]
end
function W(matrixType::Block1DParallel, ::POSVSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  ret::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, nBlocks)
  @threads for i in 1:nBlocks
    ret[i] = ((deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-1)
  end
  return ret
end
# Solver
function solve!(::POSVSolver, R::Array{T, 2},
  xwx::Array{T, 2}, xw::Array{T, 2},
  x::Array{T, 2}, z::Array{T, 1}, w::Array{T, 1},
  coef::Array{T, 1}) where {T <: AbstractFloat}
  xwx, xwz = calcXWX(xwx, x, w, z)
  _, coef = posv!('U', copy(xwx), xwz)
  return (R, xwx, xw, x, z, w, coef)
end
function solve!(::POSVSolver, R::Array{T, 2},
  xwx::Array{T, 2}, xw::Array{Array{T, 2},1},
  x::Array{Array{T, 2}, 1}, z::Array{Array{T, 1}, 1},
  w::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}
  xwx, xwz = calcXWX(xwx, x, w, z)
  _, coef = posv!('U', copy(xwx), xwz)
  return (R, xwx, xw, x, z, w, coef)
end
function solve!(matrixType::Block1DParallel, ::POSVSolver, R::Array{T, 2},
  xwx::Array{T, 2}, xw::Array{Array{T, 2},1},
  x::Array{Array{T, 2}, 1}, z::Array{Array{T, 1}, 1},
  w::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}
  xwx, xwz = calcXWX(matrixType, xwx, x, w, z)
  _, coef = posv!('U', copy(xwx), xwz)
  return (R, xwx, xw, x, z, w, coef)
end
#==============================  SYSV Solver ==============================#
#= Solver For Symmetric Indefinite Matrices =#
import LinearAlgebra.LAPACK.sysv! # (uplo, A, B) -> (B, A, ipiv)

struct SYSVSolver <: AbstractSolver end
# Weights for the SYSVSolver
@inline function W(::SYSVSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{T, 1}, eta::Array{T, 1}) where {T <: AbstractFloat}
  return ((deta_dmu(link, mu, eta).^2) .* variance(distrib, mu)).^(-1)
end
function W(::SYSVSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  return [((deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-1) for i in 1:nBlocks]
end
function W(matrixType::Block1DParallel, ::SYSVSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  ret::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, nBlocks)
  @threads for i in 1:nBlocks
    ret[i] = ((deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-1)
  end
  return ret
end
# Solver
function solve!(::SYSVSolver, R::Array{T, 2},
  xwx::Array{T, 2}, xw::Array{T, 2},
  x::Array{T, 2}, z::Array{T, 1}, w::Array{T, 1},
  coef::Array{T, 1}) where {T <: AbstractFloat}
  xwx, xwz = calcXWX(xwx, x, w, z)
  coef, _, _ = sysv!('U', copy(xwx), xwz)
  return (R, xwx, xw, x, z, w, coef)
end
function solve!(::SYSVSolver, R::Array{T, 2},
  xwx::Array{T, 2}, xw::Array{Array{T, 2},1},
  x::Array{Array{T, 2}, 1}, z::Array{Array{T, 1}, 1},
  w::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}
  xwx, xwz = calcXWX(xwx, x, w, z)
  coef, _, _ = sysv!('U', copy(xwx), xwz)
  return (R, xwx, xw, x, z, w, coef)
end
function solve!(matrixType::Block1DParallel, ::SYSVSolver, R::Array{T, 2},
  xwx::Array{T, 2}, xw::Array{Array{T, 2},1},
  x::Array{Array{T, 2}, 1}, z::Array{Array{T, 1}, 1},
  w::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}
  xwx, xwz = calcXWX(matrixType, xwx, x, w, z)
  coef, _, _ = sysv!('U', copy(xwx), xwz)
  return (R, xwx, xw, x, z, w, coef)
end
#============================== Least Squares Solvers ==============================#
#==============================  GELS Solver ==============================#
#= Solver For Symmetric Indefinite Matrices =#
import LinearAlgebra.LAPACK.gels! # (trans, A, B) -> (F, B, ssr)

struct GELSSolver <: AbstractSolver end
# Weights for the GELSSolver
@inline function W(::GELSSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{T, 1}, eta::Array{T, 1}) where {T <: AbstractFloat}
  return ((deta_dmu(link, mu, eta).^2) .* variance(distrib, mu)).^(-0.5)
end
function W(::GELSSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  return [((deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-0.5) for i in 1:nBlocks]
end
function W(matrixType::Block1DParallel, ::GELSSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  ret::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, nBlocks)
  @threads for i in 1:nBlocks
    ret[i] = ((deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-0.5)
  end
  return ret
end

# Solver
function solve!(::GELSSolver, R::Array{T, 2},
  xwx::Array{T, 2}, xw::Array{T, 2},
  x::Array{T, 2}, z::Array{T, 1}, w::Array{T, 1},
  coef::Array{T, 1}) where {T <: AbstractFloat}
  _, p = size(x)
  xw, zw = calcXW(xwx, x, w, z)
  F, coef, ssr = gels!('N', xw, zw)
  for j in 1:p
    for i in 1:j
      if i <= j
        R[i, j] = xw[i, j]
      end
    end
  end
  return (R, xwx, xw, x, z, w, coef)
end
function cov(::GELSSolver, invType::AbstractInverse, R::Array{T, 2}, xwx::Array{T, 2}, xw::Array{T, 2})::Array{T, 2} where {T <: AbstractFloat}
  xwx = R' * R
  return inv(invType, xwx)
end
#==============================  GELSY Solver ==============================#
import LinearAlgebra.LAPACK.gelsy! # (A, B, rcond) -> (B, rnk)
struct GELSYSolver <: AbstractSolver end
# Weights for the GELSYSolver
@inline function W(::GELSYSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{T, 1}, eta::Array{T, 1}) where {T <: AbstractFloat}
  return ((deta_dmu(link, mu, eta).^2) .* variance(distrib, mu)).^(-0.5)
end
function W(::GELSYSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{T, 1}, eta::Array{T, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  return [((deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-0.5) for i in 1:nBlocks]
end
function W(matrixType::Block1DParallel, ::GELSYSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{T, 1}, eta::Array{T, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  ret::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, nBlocks)
  @threads for i in 1:nBlocks
    ret[i] = ((deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-0.5)
  end
  return ret
end

# Solver
function solve!(::GELSYSolver, R::Array{T, 2},
  xwx::Array{T, 2}, xw::Array{T, 2},
  x::Array{T, 2}, z::Array{T, 1}, w::Array{T, 1},
  coef::Array{T, 1}) where {T <: AbstractFloat}
  _, p = size(x); rcond = T(0)
  xw, zw = calcXW(xwx, x, w, z)
  coef, rnk = gelsy!(copy(xw), zw, rcond)
  return (R, xwx, xw, x, z, w, coef)
end
function cov(::GELSYSolver, invType::AbstractInverse, R::Array{T, 2}, xwx::Array{T, 2}, xw::Array{T, 2})::Array{T, 2} where {T <: AbstractFloat}
  xwx = xw' * xw
  return inv(invType, xwx)
end
#==============================  GELSD Solver ==============================#
#= Solver For Symmetric Indefinite Matrices =#
import LinearAlgebra.LAPACK.gelsd! # (A, B, rcond) -> (B, rnk)

struct GELSDSolver <: AbstractSolver end
# Weights for the GELSDSolver
@inline function W(::GELSDSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{T, 1}, eta::Array{T, 1}) where {T <: AbstractFloat}
  return ((deta_dmu(link, mu, eta).^2) .* variance(distrib, mu)).^(-0.5)
end
function W(::GELSDSolver, distrib::AbstractDistribution, link::AbstractLink, mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  return [((deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-0.5) for i in 1:nBlocks]
end
# Solver
function solve!(::GELSDSolver, R::Array{T, 2},
  xwx::Array{T, 2}, xw::Array{T, 2},
  x::Array{T, 2}, z::Array{T, 1}, w::Array{T, 1},
  coef::Array{T, 1}) where {T <: AbstractFloat}
  _, p = size(x)
  rcond::T = T(0)
  xw, zw = calcXW(xwx, x, w, z)
  coef, _ = gelsd!(copy(xw), zw, rcond)
  return (R, xwx, xw, x, z, w, coef)
end
function cov(::GELSDSolver, invType::AbstractInverse, R::Array{T, 2}, xwx::Array{T, 2}, xw::Array{T, 2})::Array{T, 2} where {T <: AbstractFloat}
  xwx = xw' * xw
  return inv(invType, xwx)
end

#==============================  GRADIENT DESCENT SOLVERS ==============================#
#=======================================================================================#

function XWX(x::Array{T, 2}, w::Array{T, 1}) where {T <: AbstractFloat}
  xw = copy(x)
  xw = multSweep(xw, w)
  xwx = xw' * x
  return xwx
end

function XWX(x::Array{Array{T, 2}, 1}, w::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  p::Int64 = size(x[1])[2]
  nBlocks::Int64 = length(x)
  xwx = zeros(T, (p, p))
  for i in 1:nBlocks
    xw = copy(x[i])
    xw = multSweep(xw, w[i])
    xwx += xw' * x[i]
  end
  return xwx
end

function XWX(::Block1DParallel, x::Array{Array{T, 2}, 1}, w::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  p::Int64 = size(x[1])[2]
  nBlocks::Int64 = length(x)
  xwx = [zeros(T, (p, p)) for i in 1:nthreads()]
  @threads for i in 1:nBlocks
    xw = copy(x[i])
    xw = multSweep(xw, w[i])
    xwx[threadid()] += xw' * x[i]
  end
  for i in 2:nthreads()
    xwx[1] .+= xwx[i]
  end
  return xwx[1]
end

abstract type AbstractGradientDescentSolver end

# Weights for the GradientDescentSolver
@inline function W(::AbstractGradientDescentSolver,
            distrib::AbstractDistribution, link::AbstractLink,
            mu::Array{T, 1}, eta::Array{T, 1}) where {T <: AbstractFloat}
  return ((deta_dmu(link, mu, eta).^2) .* variance(distrib, mu)).^(-1)
end
function W(::AbstractGradientDescentSolver, distrib::AbstractDistribution,
            link::AbstractLink, mu::Array{Array{T, 1}, 1},
            eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  return [((deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-1) for i in 1:nBlocks]
end
function W(matrixType::Block1DParallel, ::AbstractGradientDescentSolver,
            distrib::AbstractDistribution, link::AbstractLink,
            mu::Array{Array{T, 1}, 1},
            eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}
  nBlocks::Int64 = length(mu)
  ret::Array{Array{T, 1}, 1} = Array{Array{T, 1}, 1}(undef, nBlocks)
  @threads for i in 1:nBlocks
    ret[i] = ((deta_dmu(link, mu[i], eta[i]).^2) .* variance(distrib, mu[i])).^(-1)
  end
  return ret
end

# GD Covariance Function
function cov(::AbstractGradientDescentSolver, invType::AbstractInverse, xwx::Array{T, 2})::Array{T, 2} where {T <: AbstractFloat}
  return inv(invType, xwx)
end

#==============================  GRADIENTS ==============================#

#=
  Internal gradient function
=#
function _gradient(distrib::AbstractDistribution, link::AbstractLink,
            y::Array{T, 1}, x::Array{T, 2}, mu::Array{T, 1},
            eta::Array{T, 1}) where {T <: AbstractFloat}
  n, p = size(x)
  grad = zeros(T, p)
  numer::Array{T, 1} = y .- mu
  X2::T = sum((numer.^2) ./ variance(distrib, mu))
  tmp::Array{T, 1} = numer ./(deta_dmu(link, mu, eta) .* variance(distrib, mu))
  for j in 1:p
    grad[j] += sum(tmp .* x[:, j])
  end
  return (X2 = X2, grad = grad)
end

#=
  Gradient function for regular data
=#
function gradient(::RegularData, distrib::AbstractDistribution, link::AbstractLink,
            y::Array{T, 1}, x::Array{T, 2}, mu::Array{T, 1},
            eta::Array{T, 1}) where {T <: AbstractFloat}
  n, p = size(x)
  tmp = _gradient(distrib, link, y, x, mu, eta)
  df = n - p
  @assert(df > 0, "Number of items not greater than number of parameters")
  phi::T = tmp.X2/df
  grad = tmp.grad ./phi
  return grad
end

#=
  Gradient function for blocked data
=#
function gradient(::Block1D, distrib::AbstractDistribution, link::AbstractLink,
  y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}

  nBlocks::Int64 = length(y); n::Int64 = 0
  p::Int64 = size(x[1])[2]; grad::Array{T, 1} = zeros(T, p)
  X2::T = T(0)

  for i in 1:nBlocks
    n += length(y[i])
    tmp = _gradient(distrib, link, y[i], x[i], mu[i], eta[i])
    # println("size(grad): ", size(grad))
    # println("size(tmp.grad): ", size(tmp.grad))
    grad .+= tmp.grad
    X2 += tmp.X2
  end

  df = n - p
  phi::T = X2/df
  @assert(df > 0, "Number of items not greater than number of parameters")
  grad ./= phi
  return grad
end

function gradient(::Block1DParallel, distrib::AbstractDistribution, link::AbstractLink,
  y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}) where {T <: AbstractFloat}

  nBlocks::Int64 = length(y); p::Int64 = size(x[1])[2]; n::Int64 = 0;
  # grad::Array{T, 1} = zeros(T, p); X2::T = T(0)

  nStore = [Int64(0) for i in 1:nthreads()]
  X2Store = [T(0) for i in 1:nthreads()]
  gradStore = [Base.zeros(T, p) for i in 1:nthreads()]

  distribStore = [copy(distrib) for i in 1:nthreads()]
  linkStore = [copy(link) for i in 1:nthreads()]

  @threads for i in 1:nBlocks
    tmp = _gradient(distribStore[threadid()], linkStore[threadid()], y[i], x[i], mu[i], eta[i])
    gradStore[threadid()] .+= tmp.grad
    X2Store[threadid()] += tmp.X2
    nStore[threadid()] += length(y[i])
  end

  n = Int64(0); grad = zeros(T, p); X2 = T(0)
  for i in 1:nthreads()
    # nStore[1] += nStore[i]
    # gradStore[1] .+= gradStore[i]
    # X2Store[1] += X2Store[i]
    n += nStore[i]
    grad .+= gradStore[i]
    X2 += X2Store[i]
  end

  # df = nStore[1] - p
  # phi::T = X2Store[1]/df
  df = n - p
  phi::T = X2/df
  @assert(df > 0, "Number of items not greater than number of parameters")

  # return gradStore[1]./phi
  grad ./= phi
  return grad
end

#==============================  Simple Gradient Descent Solver ==============================#

mutable struct GradientDescentSolver{T} <: AbstractGradientDescentSolver
  learningRate::T
  function GradientDescentSolver(learningRate::T = T(1E-6)) where {T <: AbstractFloat}
    return new{T}(learningRate)
  end
end
#= Copy constructor =#
function copy(solver::GradientDescentSolver)
  return GradientDescentSolver(solver.learningRate)
end

function solve!(dataType::RegularData, solver::GradientDescentSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{T, 1}, x::Array{T, 2}, mu::Array{T, 1},
  eta::Array{T, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  grad = gradient(dataType, distrib, link, y, x, mu, eta);
  coef .+= solver.learningRate .* grad;
  return solver, coef
end

function solve!(dataType::Block1D, solver::GradientDescentSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  grad = gradient(dataType, distrib, link, y, x, mu, eta);
  coef .+= solver.learningRate .* grad;
  return solver, coef
end

function solve!(dataType::Block1DParallel, solver::GradientDescentSolver,
  distrib::AbstractDistribution, link::AbstractLink,
  y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  grad = gradient(dataType, distrib, link, y, x, mu, eta);
  coef .+= solver.learningRate .* grad;
  return solver, coef
end

#==============================  Gradient Descent Momentum Solver ==============================#

mutable struct MomentumSolver{T} <: AbstractGradientDescentSolver
  learningRate::T
  momentum::T
  delta::Array{T, 1}
  function MomentumSolver(learningRate::T, momentum::T, p::Int64) where {T <: AbstractFloat}
    return new{T}(learningRate, momentum, zeros(T, p))
  end
  function MomentumSolver(learningRate::T, momentum::T, delta::Array{T, 1}) where {T <: AbstractFloat}
    return new{T}(learningRate, momentum, delta)
  end
end
#= Copy constructor =#
function copy(solver::MomentumSolver)
  return MomentumSolver(solver.learningRate, solver.momentum, copy(solver.delta))
end

function solve!(dataType::RegularData, solver::MomentumSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{T, 1}, x::Array{T, 2}, mu::Array{T, 1},
  eta::Array{T, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  grad = gradient(dataType, distrib, link, y, x, mu, eta)
  solver.delta .= (solver.momentum .* solver.delta) .+ solver.learningRate .* grad
  coef .+= solver.delta
  return solver, coef
end

function solve!(dataType::Block1D, solver::MomentumSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  grad = gradient(dataType, distrib, link, y, x, mu, eta)
  solver.delta .= (solver.momentum .* solver.delta) .+ solver.learningRate .* grad
  coef .+= solver.delta
  return solver, coef
end

function solve!(dataType::Block1DParallel, solver::MomentumSolver,
  distrib::AbstractDistribution, link::AbstractLink,
  y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  grad = gradient(dataType, distrib, link, y, x, mu, eta)
  solver.delta .= (solver.momentum .* solver.delta) .+ solver.learningRate .* grad
  coef .+= solver.delta
  return solver, coef
end

#==============================  Gradient Descent Nesterov Solver ==============================#
#= Default Modify/Unmodify coefficients for Nesterov =#
function NesterovModifier(obj::AbstractGradientDescentSolver, coef::Array{T, 1}) where {T <: AbstractFloat}
  return coef
end
function NesterovUnModifier(obj::AbstractGradientDescentSolver, coef::Array{T, 1}) where {T <: AbstractFloat}
  return coef
end

mutable struct NesterovSolver{T} <: AbstractGradientDescentSolver
  learningRate::T
  momentum::T
  delta::Array{T, 1}
  function NesterovSolver(learningRate::T, momentum::T, p::Int64) where {T <: AbstractFloat}
    return new{T}(learningRate, momentum, zeros(T, p))
  end
  function NesterovSolver(learningRate::T, momentum::T, delta::Array{T, 1}) where {T <: AbstractFloat}
    return new{T}(learningRate, momentum, delta)
  end
end
#= Copy constructor =#
function copy(solver::NesterovSolver)
  return NesterovSolver(solver.learningRate, solver.momentum, copy(solver.delta))
end

function solve!(dataType::RegularData, solver::NesterovSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{T, 1}, x::Array{T, 2}, mu::Array{T, 1},
  eta::Array{T, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  grad = gradient(dataType, distrib, link, y, x, mu, eta)
  solver.delta .= (solver.momentum .* solver.delta) .+ solver.learningRate .* grad
  coef .+= solver.delta
  return solver, coef
end

function solve!(dataType::Block1D, solver::NesterovSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  grad = gradient(dataType, distrib, link, y, x, mu, eta)
  solver.delta .= (solver.momentum .* solver.delta) .+ solver.learningRate .* grad
  coef .+= solver.delta
  return solver, coef
end

function solve!(dataType::Block1DParallel, solver::NesterovSolver,
  distrib::AbstractDistribution, link::AbstractLink,
  y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  grad = gradient(dataType, distrib, link, y, x, mu, eta)
  solver.delta .= (solver.momentum .* solver.delta) .+ solver.learningRate .* grad
  coef .+= solver.delta
  return solver, coef
end

function NesterovModifier(obj::NesterovSolver, coef::Array{T, 1}) where {T <: AbstractFloat}
  coef .-= obj.momentum .* obj.delta
  return coef
end
function NesterovUnModifier(obj::NesterovSolver, coef::Array{T, 1}) where {T <: AbstractFloat}
  coef .+= obj.momentum .* obj.delta
  return coef
end

#==============================  Gradient Descent Adagrad Solver ==============================#
#= Modify/Unmodify coefficients for Adagrad =#
mutable struct AdagradSolver{T} <: AbstractGradientDescentSolver
  learningRate::T
  G::Array{T, 1}
  epsilon::T
  function AdagradSolver(learningRate::T, p::Int64, epsilon::T) where {T <: AbstractFloat}
    return new{T}(learningRate, zeros(T, p), epsilon)
  end
  function AdagradSolver(learningRate::T, G::Array{T, 1}, epsilon::T) where {T <: AbstractFloat}
    return new{T}(learningRate, G, epsilon)
  end
end
#= Copy constructor =#
function copy(solver::AdagradSolver)
  return AdagradSolver(solver.learningRate, copy(solver.G), solver.epsilon)
end

function solve!(dataType::RegularData, solver::AdagradSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{T, 1}, x::Array{T, 2}, mu::Array{T, 1},
  eta::Array{T, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  grad = gradient(dataType, distrib, link, y, x, mu, eta)
  solver.G += coef.^2
  coef .+= (solver.learningRate .* grad)./((solver.G .+ solver.epsilon).^0.5)
  return solver, coef
end

function solve!(dataType::Block1D, solver::AdagradSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  grad = gradient(dataType, distrib, link, y, x, mu, eta)
  solver.G += coef.^2
  coef .+= (solver.learningRate .* grad)./((solver.G .+ solver.epsilon).^0.5)
  return solver, coef
end

function solve!(dataType::Block1DParallel, solver::AdagradSolver,
  distrib::AbstractDistribution, link::AbstractLink,
  y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  grad = gradient(dataType, distrib, link, y, x, mu, eta)
  solver.G += coef.^2
  coef .+= (solver.learningRate .* grad)./((solver.G .+ solver.epsilon).^0.5)
  return solver, coef
end

#==============================  Gradient Descent Adadelta Solver ==============================#
#= Modify/Unmodify coefficients for Adagrad =#
mutable struct AdadeltaSolver{T} <: AbstractGradientDescentSolver
  momentum::T
  G::Array{T, 1}
  B::Array{T, 1}
  epsilon::T
  function AdadeltaSolver(momentum::T, p::Int64, epsilon::T) where {T <: AbstractFloat}
    return new{T}(momentum, zeros(T, p), zeros(T, p), epsilon)
  end
  function AdadeltaSolver(momentum::T, G::Array{T, 1}, B::Array{T, 1}, epsilon::T) where {T <: AbstractFloat}
    return new{T}(momentum, G, B, epsilon)
  end
end
#= Copy constructor =#
function copy(solver::AdadeltaSolver)
  return AdadeltaSolver(solver.momentum, copy(solver.G), copy(solver.B), solver.epsilon)
end

macro AdadeltaUpdate()
  expr = quote
    grad = gradient(dataType, distrib, link, y, x, mu, eta)
    momentumM1 = (1 - solver.momentum)
    solver.G = (solver.momentum .* solver.G) .+ (momentumM1 .* (grad.^2))
    diff = (grad .* ((solver.B .+ solver.epsilon).^0.5)) ./((solver.G .+ solver.epsilon).^0.5)
    coef .+= diff
    solver.B .+= (solver.momentum .* solver.B) .+ (momentumM1 .* (diff.^2))
    return solver, coef
  end
  return esc(expr)
end

function solve!(dataType::RegularData, solver::AdadeltaSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{T, 1}, x::Array{T, 2}, mu::Array{T, 1},
  eta::Array{T, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  @AdadeltaUpdate
end

function solve!(dataType::Block1D, solver::AdadeltaSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  @AdadeltaUpdate
end

function solve!(dataType::Block1DParallel, solver::AdadeltaSolver,
  distrib::AbstractDistribution, link::AbstractLink,
  y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  @AdadeltaUpdate
end

#==============================  Gradient Descent RMSprop Solver ==============================#
#= Modify/Unmodify coefficients for Adagrad =#
mutable struct RMSpropSolver{T} <: AbstractGradientDescentSolver
  learningRate::T
  G::Array{T, 1}
  epsilon::T
  function RMSpropSolver(learningRate::T, p::Int64, epsilon::T) where {T <: AbstractFloat}
    return new{T}(learningRate, zeros(T, p), epsilon)
  end
  function RMSpropSolver(learningRate::T, G::Array{T, 1}, epsilon::T) where {T <: AbstractFloat}
    return new{T}(learningRate, G, epsilon)
  end
end
#= Copy constructor =#
function copy(solver::RMSpropSolver)
  return RMSpropSolver(solver.learningRate, copy(solver.G), solver.epsilon)
end

macro RMSpropUpdate()
  expr = quote
    grad = gradient(dataType, distrib, link, y, x, mu, eta)
    solver.G = (0.9 .* solver.G) .+ (0.1 .* (grad.^2))
    diff = (solver.learningRate .* grad) ./((solver.G .+ solver.epsilon).^0.5)
    coef .+= diff
    return solver, coef
  end
  return esc(expr)
end

function solve!(dataType::RegularData, solver::RMSpropSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{T, 1}, x::Array{T, 2}, mu::Array{T, 1},
  eta::Array{T, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  @RMSpropUpdate
end

function solve!(dataType::Block1D, solver::RMSpropSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  @RMSpropUpdate
end

function solve!(dataType::Block1DParallel, solver::RMSpropSolver,
  distrib::AbstractDistribution, link::AbstractLink,
  y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  @RMSpropUpdate
end

#==============================  Gradient Descent Adam Solver ==============================#
function iteration(obj::AbstractGradientDescentSolver, iter::Int64)
  return obj
end
#= Adam =#
mutable struct AdamSolver{T} <: AbstractGradientDescentSolver
  learningRate::T
  b1::T
  b2::T
  m::Array{T, 1}
  v::Array{T, 1}
  epsilon::T
  iter::Int64
  function AdamSolver(learningRate::T, b1::T, b2::T, p::Int64, epsilon::T) where {T <: AbstractFloat}
    return new{T}(learningRate, b1, b2, zeros(T, p), zeros(T, p), epsilon, Int64(1))
  end
  function AdamSolver(learningRate::T, b1::T, b2::T, m::Array{T, 1}, v::Array{T, 1}, epsilon::T, iter::Int64) where {T <: AbstractFloat}
    return new{T}(learningRate, b1, b2, m, v, epsilon, iter)
  end
end
#= Copy constructor =#
function copy(solver::AdamSolver)
  return AdamSolver(solver.learningRate, solver.b1, solver.b2, copy(solver.m), copy(solver.v), solver.epsilon, solver.iter)
end
function iteration(obj::AdamSolver, iter::Int64)
  obj.iter = iter
  return obj
end

macro AdamUpdate()
  expr = quote
    grad = gradient(dataType, distrib, link, y, x, mu, eta)
    solver.m = (solver.b1 .* solver.m) .+ ((1 - solver.b1) .* grad)
    solver.v = (solver.b2 .* solver.v) .+ ((1 - solver.b2) .* (grad.^2))
    mp = solver.m ./(1 - (solver.b1^solver.iter))
    vp = solver.v ./(1 - (solver.b2^solver.iter))
    diff = (solver.learningRate .* mp)./(solver.epsilon .+ (vp.^0.5))
    coef .+= diff
    return solver, coef
  end
  return esc(expr)
end

function solve!(dataType::RegularData, solver::AdamSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{T, 1}, x::Array{T, 2}, mu::Array{T, 1},
  eta::Array{T, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  @AdamUpdate
end

function solve!(dataType::Block1D, solver::AdamSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  @AdamUpdate
end

function solve!(dataType::Block1DParallel, solver::AdamSolver,
  distrib::AbstractDistribution, link::AbstractLink,
  y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  @AdamUpdate
end

#==============================  Gradient Descent AdaMax Solver ==============================#
#= AdaMax =#
mutable struct AdaMaxSolver{T} <: AbstractGradientDescentSolver
  learningRate::T
  b1::T
  b2::T
  m::Array{T, 1}
  v::T
  iter::Int64
  function AdaMaxSolver(learningRate::T, b1::T, b2::T, p::Int64) where {T <: AbstractFloat}
    return new{T}(learningRate, b1, b2, zeros(T, p), T(0), Int64(1))
  end
  function AdaMaxSolver(learningRate::T, b1::T, b2::T, m::Array{T, 1}, v::T, iter::Int64) where {T <: AbstractFloat}
    return new{T}(learningRate, b1, b2, m, v, iter)
  end
end
#= Copy constructor =#
function copy(solver::AdaMaxSolver)
  return AdaMaxSolver(solver.learningRate, solver.b1, solver.b2, copy(solver.m), solver.v, solver.iter)
end
function iteration(obj::AdaMaxSolver, iter::Int64)
  obj.iter = iter
  return obj
end

macro AdaMaxUpdate()
  expr = quote
    grad = gradient(dataType, distrib, link, y, x, mu, eta)

    absG = norm(grad)
    solver.v = (solver.b2 * solver.v) + ((1 - solver.b2) * absG)
    u = max(solver.b2 * solver.v, absG)
    solver.m = (solver.b1 .* solver.m) .+ ((1 - solver.b1) .* grad)
    mp = solver.m ./(1 - (solver.b1^solver.iter))
    diff = (solver.learningRate .* mp)./u

    coef .+= diff
    return solver, coef
  end
  return esc(expr)
end

function solve!(dataType::RegularData, solver::AdaMaxSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{T, 1}, x::Array{T, 2}, mu::Array{T, 1},
  eta::Array{T, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  @AdaMaxUpdate
end

function solve!(dataType::Block1D, solver::AdaMaxSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  @AdaMaxUpdate
end

function solve!(dataType::Block1DParallel, solver::AdaMaxSolver,
  distrib::AbstractDistribution, link::AbstractLink,
  y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  @AdaMaxUpdate
end

#==============================  Gradient Descent NAdam Solver ==============================#
#= NAdam =#
mutable struct NAdamSolver{T} <: AbstractGradientDescentSolver
  learningRate::T
  b1::T
  b2::T
  m::Array{T, 1}
  v::Array{T, 1}
  epsilon::T
  iter::Int64
  function NAdamSolver(learningRate::T, b1::T, b2::T, p::Int64, epsilon::T) where {T <: AbstractFloat}
    return new{T}(learningRate, b1, b2, zeros(T, p), zeros(T, p), epsilon, Int64(1))
  end
  function NAdamSolver(learningRate::T, b1::T, b2::T, m::Array{T, 1}, v::Array{T, 1}, epsilon::T, iter::Int64) where {T <: AbstractFloat}
    return new{T}(learningRate, b1, b2, m, v, epsilon, iter)
  end
end
#= Copy constructor =#
function copy(solver::NAdamSolver)
  return NAdamSolver(solver.learningRate, solver.b1, solver.b2, copy(solver.m), copy(solver.v), solver.epsilon, solver.iter)
end
function iteration(obj::NAdamSolver, iter::Int64)
  obj.iter = iter
  return obj
end

macro NAdamUpdate()
  expr = quote
    grad = gradient(dataType, distrib, link, y, x, mu, eta)

    solver.m .= (solver.b1 .* solver.m) .+ ((1 - solver.b1) .* grad)
    mp = solver.m ./(1 - (solver.b1^solver.iter))
    solver.v .= (solver.b2 .* solver.v) .+ ((1 - solver.b2) .* (grad.^2))
    vp = solver.v ./(1 - (solver.b2.^solver.iter))
    tmp1 = solver.learningRate ./((vp.^0.5) .+ solver.epsilon);
    tmp2 = (solver.b1 .* mp) .+ ((1 - solver.b1) .* grad) ./(1 - (solver.b1^solver.iter))

    coef .+= (tmp1 .* tmp2)

    return solver, coef
  end
  return esc(expr)
end

function solve!(dataType::RegularData, solver::NAdamSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{T, 1}, x::Array{T, 2}, mu::Array{T, 1},
  eta::Array{T, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  @NAdamUpdate
end

function solve!(dataType::Block1D, solver::NAdamSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  @NAdamUpdate
end

function solve!(dataType::Block1DParallel, solver::NAdamSolver,
  distrib::AbstractDistribution, link::AbstractLink,
  y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  @NAdamUpdate
end

#==============================  Gradient Descent AMSGrad Solver ==============================#
#= AMSGrad =#
mutable struct AMSGradSolver{T} <: AbstractGradientDescentSolver
  learningRate::T
  b1::T
  b2::T
  m::Array{T, 1}
  v::Array{T, 1}
  vp::Array{T, 1}
  epsilon::T
  iter::Int64
  function AMSGradSolver(learningRate::T, b1::T, b2::T, p::Int64, epsilon::T) where {T <: AbstractFloat}
    return new{T}(learningRate, b1, b2, zeros(T, p), zeros(T, p), zeros(T, p), epsilon, Int64(1))
  end
  function AMSGradSolver(learningRate::T, b1::T, b2::T, m::Array{T, 1}, v::Array{T, 1}, vp::Array{T, 1}, epsilon::T, iter::Int64) where {T <: AbstractFloat}
    return new{T}(learningRate, b1, b2, m, v, vp, epsilon, iter)
  end
end
#= Copy constructor =#
function copy(solver::AMSGradSolver)
  return AMSGradSolver(solver.learningRate, solver.b1, solver.b2, copy(solver.m), copy(solver.v), copy(solver.vp), solver.epsilon, solver.iter)
end
function iteration(obj::AMSGradSolver, iter::Int64)
  obj.iter = iter
  return obj
end

macro AMSGradUpdate()
  expr = quote
    grad = gradient(dataType, distrib, link, y, x, mu, eta)

    solver.m .= (solver.b1 .* solver.m) .+ ((1 - solver.b1) .* grad)
    solver.v .= (solver.b2 .* solver.v) .+ ((1 - solver.b2) .* (grad.^2))
    solver.vp .= max.(solver.vp, solver.v)
    diff = (solver.learningRate .* solver.m)./((solver.vp .^0.5) .+ solver.epsilon)

    coef .+= diff

    return solver, coef
  end
  return esc(expr)
end

function solve!(dataType::RegularData, solver::AMSGradSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{T, 1}, x::Array{T, 2}, mu::Array{T, 1},
  eta::Array{T, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  @AMSGradUpdate
end

function solve!(dataType::Block1D, solver::AMSGradSolver, distrib::AbstractDistribution,
  link::AbstractLink, y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  @AMSGradUpdate
end

function solve!(dataType::Block1DParallel, solver::AMSGradSolver,
  distrib::AbstractDistribution, link::AbstractLink,
  y::Array{Array{T, 1}, 1}, x::Array{Array{T, 2}, 1},
  mu::Array{Array{T, 1}, 1}, eta::Array{Array{T, 1}, 1}, coef::Array{T, 1}) where {T <: AbstractFloat}

  @AMSGradUpdate
end

#==============================  GRADIENT DESCENT INITIALIZERS ==============================#

# init type return.
InitDouble{T} = Tuple{Array{T, 1}, Array{T, 1}} # where {T <: AbstractFloat}

# Initializer for gradient descent solvers
function init!(::AbstractGradientDescentSolver,
               ::AbstractDistribution, y::Array{T},
               wts::Array{T, 1}) #= ::InitDouble{T} =# where {T <: AbstractFloat}
  y = y[:, 1]
  return y, wts
end

function init!(::AbstractGradientDescentSolver,
               ::AbstractDistribution, y::Array{Array{T, 2}, 1},
               wts::Array{Array{T, 1}, 1}) #= ::InitDouble{Array{T, 1}}=# where {T <: AbstractFloat}
  nBlocks::Int64 = length(y)
  if length(size(y[1])) == 2
    y = [y[i][:, 1] for i in 1:nBlocks]
  end
  return y, wts;
end

function init!(::AbstractGradientDescentSolver,
               ::Block1DParallel, ::AbstractDistribution,
               y::Array{Array{T, 2}, 1},
               wts::Array{Array{T, 1}, 1}) #=::InitDouble{Array{T, 1}} =# where {T <: AbstractFloat}
  nBlocks::Int64 = length(y)
  isMat::Bool = length(size(y[1])) == 2
  yNew = Array{Array{T, 1}, 1}(undef, nBlocks);
  @threads for i in 1:nBlocks
    yNew[i] = isMat ? y[i][:, 1] : y[i]
  end
  return yNew, wts;
end

#===================== BinomialDistribution Initializers ====================#
function init!(::AbstractGradientDescentSolver,
      ::BinomialDistribution, y::Array{T},
      wts::Array{T, 1}) #= ::InitDouble{T}=# where {T <: AbstractFloat}
  if size(y)[2] == 1
    y = y[:, 1]
  elseif size(y)[2] == 2
    n = size(y)[1]
    events = y[:, 1]
    N = events .+ y[:, 2]
    y = [if N[i] != 0; events[i]/N[i]; else T(0) end for i in 1:n]
    if length(wts) != 0
      wts .*= N
    else
      wts = N
    end
  else
    error("There was a problem")
  end
  return y, wts
end

function init!(::AbstractGradientDescentSolver,
               ::BinomialDistribution, y::Array{Array{T, 2}, 1},
               wts::Array{Array{T, 1}, 1}) #= ::InitDouble{Array{T, 1}} =# where {T <: AbstractFloat}
  nBlocks::Int64 = length(y)
  if size(y[1])[2] == 1
    y = [y[i][:, 1] for i in 1:nBlocks]
  elseif size(y[1])[2] == 2
    events = [y[i][:, 1] for i in 1:nBlocks]
    N = [events[i] .+ y[i][:, 2] for i in 1:nBlocks]
    yNew = Array{Array{T, 1}, 1}(undef, nBlocks)
    for i in 1:nBlocks
      n = size(y[i])[1]
      tmp = zeros(T, n)
      for j in 1:n
        tmp[j] = if N[i][j] != 0; events[i][j]/N[i][j]; else T(0) end
      end
      yNew[i] = tmp
    end
    y = yNew
    if length(wts) != 0
      [wts[i] .*= N[i] for i in 1:nBlocks]
    else
      wts = N
    end
  else
    error("There was a problem with the dimensions of y")
  end
  return y, wts
end

function init!(::AbstractGradientDescentSolver,
          ::Block1DParallel, ::BinomialDistribution,
          y::Array{Array{T, 2}, 1},
          wts::Array{Array{T, 1}, 1}) #= ::InitDouble{Array{T, 1}} =# where {T <: AbstractFloat}
  nBlocks::Int64 = length(y)
  yNew = Array{Array{T, 1}, 1}(undef, nBlocks)
  events = Array{Array{T, 1}, 1}(undef, nBlocks)
  N = Array{Array{T, 1}, 1}(undef, nBlocks)
  if size(y[1])[2] == 1
    @threads for i in 1:nBlocks
      yNew[i] = y[i][:, 1]
    end
    y = yNew
  elseif size(y[1])[2] == 2
    @threads for i in 1:nBlocks
      events[i] = y[i][:, 1]
    end
    @threads for i in 1:nBlocks
      N[i] = events[i] .+ y[i][:, 2]
    end
    @threads for i in 1:nBlocks
      n = size(y[i])[1]
      tmp = zeros(T, n)
      for j in 1:n
        tmp[j] = if N[i][j] != 0; events[i][j]/N[i][j]; else T(0) end
      end
      yNew[i] = tmp
    end
    y = yNew
    if length(wts) != 0
      @threads for i in 1:nBlocks
        wts[i] = wts[i] .* N[i]
      end
    else
      wts = N
    end
  else
    error("There was a problem with the dimensions of y")
  end
  return y, wts
end

max_stars_repo_name: dataPulverizer/glmSolver

max_issues_repo_name: dataPulverizer/glmSolver

max_forks_repo_name: dataPulverizer/glmSolver


text:
using Test, SafeTestsets

@safetestset "TopOptProblems.jl" begin
    include("problems.jl")
    include("metadata.jl")
end

ext: jl


text:
# compare-hyperparameters-with-kernels
#
# This script compares the quality of samples generated from different
# step size hyperparameter settings of a Stochastic Gradient Langevin Dynamics
# sampler for a Gaussian mixture model posterior target distribution.

include("experiment_utils.jl")

# set the discrepancy
@parsestringcli discrepancytype "k" "discrepancytype" "maternradial"

using Iterators: product

using SteinDistributions:
    SteinGMMPosterior,
    randgmm,
    numdimensions,
    numdatapoints,
    gradlogdensity
using SteinDiscrepancy:
    SteinGaussianKernel,
    SteinMaternRadialKernel,
    SteinGaussianPowerKernel,
    SteinInverseMultiquadricKernel,
    ksd
using SteinSamplers: runsgld

## Experiment settings
# Target size of sample to draw from sampler
n = 1000
# Run one complete experiment per random number generator seed
seeds = 1:1
# Sampler settings for SGLD
distname = "sgld-gmm"
# set the kernel
kernel = nothing
if discrepancytype == "maternradial"
    kernel = SteinMaternRadialKernel()
elseif discrepancytype == "gaussian"
    kernel = SteinGaussianKernel()
elseif discrepancytype == "gaussianpower"
    kernel = SteinGaussianPowerKernel(2.0)
elseif discrepancytype == "inversemultiquadric"
    kernel = SteinInverseMultiquadricKernel(0.5)
end

## Target distribution settings for SteinGMMPosterior
# Generate random dataset from GMM
truex = [0.0, 1.0] # Parameter vector used to generate dataset
L = 100 # Number of datapoints to generate
data_generation_seed = 1
y = @setseed data_generation_seed randgmm(truex, L)
# save the y values for later
save_json(
    Dict{Any, Any}("y" => y);
    dir="compare-hyperparameters-with-kernels-y",
    numsamples=L,
    x=string(truex)
)
# Create posterior distribution based on dataset
target = SteinGMMPosterior(y)
# define gradlogp
function gradlogp(x::Array{Float64,1})
    gradlogdensity(target, x)
end
# Fix mini-batch size
batchsize = 5
# Number of sweeps to make through the data
numsweeps = round(Int, ceil(n / floor(numdatapoints(target)/batchsize)))
# Hyperparameter to vary: step size epsilon
epsilons = [5e-2, 1e-2, 5e-3, 1e-3, 5e-4, 1e-4, 5e-5, 1e-5]
# Define SGLD sampler with step size epsilon and starting sample value x0
sampler(x0,epsilon) = runsgld(target, x0;
                              epsilonfunc=t -> epsilon,
                              numsweeps=numsweeps,
                              batchsize=batchsize)

## Draw and evaluate samples for different hyperparameter and seed combinations
d = numdimensions(target)
for (epsilon, seed) in product(epsilons, seeds)
    # Generate starting sample value x0
    x0 = @setseed seed randn(d)

    println("Generating sample for seed=$(seed), epsilon=$(epsilon)")
    # Run sampler
    X = []; numgrad = NaN
    @trycatchcontinue(
        begin
            (X, numgrad) = @setseed seed sampler(x0,epsilon)
            if any(isnan(X))
                println("[$(distname):seed=$(seed),epsilon=$(epsilon)] NaNs found. Skipping.")
                continue
            end
        end,
        println("[$(distname):seed=$(seed),epsilon=$(epsilon)]:")
    )

    # Compute Stein discrepancy for first n points in sample
    res = nothing
    @trycatchcontinue(
        begin
            res = ksd(points=X[1:n,:], gradlogdensity=gradlogp, kernel=kernel)
        end,
        println("[n=$(n)|sampler=$(distname)]:")
    )
    if res == nothing
        continue
    end
    # compute the coercive function
    kerneldiscrepancy = sqrt(res.discrepancy2)

    # Package results
    instance_data = Dict{Any, Any}(
        "distname" => distname,
        "discrepancytype" => discrepancytype,
        "n" => n,
        "d" => d,
        "seed" => seed,
        "epsilon" => epsilon,
        "X" => X[1:n,:],
        "q" => res.weights,
        "steindiscrepancy" => kerneldiscrepancy,
        "numgrad" => numgrad * (n / size(X, 1)),
        "solvetime" => res.solvetime,
    )

    # Save results
    save_json(
        instance_data;
        dir="compare-hyperparameters-with-kernels",
        distname=distname,
        discrepancytype=discrepancytype,
        n=n,
        epsilon=epsilon,
        d=d,
        seed=seed
    )
end


text:
module Day2

using Advent: save_input

const regex = r"([0-9]+)-([0-9]+) ([a-z]): ([a-z]+)$"

function parse_line(line)
    low, high, letter, password = match(regex, line).captures
    parse(Int, low), parse(Int, high), only(letter), password
end

function solve_a(filename)
    count(eachline(filename)) do line
        low, high, letter, password = parse_line(line)
        count(==(letter), password) in low:high
    end
end

function solve_b(filename)
    count(eachline(filename)) do line
        low, high, letter, password = parse_line(line)
        (password[low] == letter) ⊻ (password[high] == letter)
    end
end

solve(filename = save_input(2020, 2)) = solve_a(filename), solve_b(filename)

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
using Logging
using Random

"Return a random index to be filled from the garden mask."
function randomindex(mask::Matrix{Bool})::Int
    while true
        i = rand(1:length(mask))
        if mask[i]
            return i
        end
    end
end

"Swap to the elements corresponding to the two provided indices."
function swap!(garden::Matrix{Int}, i::Int, j::Int)
    t = garden[i]
    garden[i] = garden[j]
    garden[j] = t
    garden
end

"Return the neighbours to be filled of the cell at the given index."
function neighbours(garden::Matrix{Int}, idx::Int)::Vector{Int}
    m, n = size(garden)
    i, j = Tuple(CartesianIndices(garden)[idx])
    neighbourindices = [(i, j-1), (i, j+1), (i-1, j), (i+1, j)]
    # cells filled with 0 are not part of the garden
    [
        garden[k, l] for (k, l) in neighbourindices
        if 0 < k <= m && 0 < l <= n && garden[k, l] != 0
    ]
end

"Compute the cost difference when swapping the two provided indices."
function deltacost(garden::Matrix{Int}, costs::Matrix{Float64}, i::Int, j::Int)::Float64
    cost = 0
    for k in neighbours(garden, i)
        cost += costs[k, garden[j]] - costs[k, garden[i]]
    end
    for k in neighbours(garden, j)
        cost += costs[k, garden[i]] - costs[k, garden[j]]
    end
    cost
end

"Compute the total cost of a garden configuration."
function gardencost(garden::Matrix{Int}, mask::Matrix{Bool}, costs::Matrix{Float64})::Float64
    cost = 0
    for i = 1:length(mask)
        if mask[i] == 0
            continue
        end
        for k in neighbours(garden, i)
            cost += costs[k, garden[i]]
        end
    end
    cost / 2
end

"Update the garden using Metropolis-Hastings, using the inverse temperature beta."
function update!(
    garden::Matrix{Int},
    mask::Matrix{Bool},
    costs::Matrix{Float64},
    beta::Float64 = 5.0
)
    N = length(garden)
    i = randomindex(mask)
    j = randomindex(mask)
    while i == j
        j = randomindex(mask)
    end
    d = deltacost(garden, costs, i, j)
    @debug "cost difference $d"
    if rand() < exp(- beta * d)
        @debug "swapping indices $i and $j"
        swap!(garden, i, j)
        return d
    end
    return 0.0
end

"Fill the garden randomly with a predefined number for each plant."
function fillgardenrandomly!(garden::Matrix{Int}, mask::Matrix{Bool}, plants::DataFrame)
    cells = vcat([repeat([plant], count) for (plant, count) in eachrow(plants)]...)
    # fill the remaining slots with random plants
    diffcount = sum(mask) - length(cells)
    cells = vcat(cells, rand(cells, diffcount))
    garden[mask] = shuffle!(indexin(cells, plants.name))
    garden
end

"Update the garden for a given number of steps, and return the total cost over time."
function gardenevolution!(garden::Matrix{Int}, mask::Matrix{Bool}, costs::Matrix{Float64}; steps::Int = 10000, beta::Float64 = 5.0)
    gardencosts = [gardencost(garden, mask, costs)]
    for i = 1:steps
        update!(garden, mask, costs, beta)
        if mod(i, 1000) == 0
            append!(gardencosts, gardencost(garden, mask, costs))
        end
    end
    gardencosts
end

max_stars_repo_name: dlozeve/GardenOptim

max_issues_repo_name: dlozeve/GardenOptim

max_forks_repo_name: dlozeve/GardenOptim


text:
using MCMCStorage, Test, DocStringExtensions

using MCMCStorage.StanCSV: parse_variable_name, collapse_contiguous_dimensions,
   parse_schema

####
#### MCMCChains
####

@testset "schema, layout, views" begin
    named_dims = (a = (), b = (1, 2), c = (2, 3, 4))
    s = Chains.IndexSchema(named_dims)
    l = vcat([(:a, )],
             vec([(:b, Tuple(i)...) for i in CartesianIndices((1:1, 1:2))]),
             vec([(:c, Tuple(i)...) for i in CartesianIndices((1:2, 1:3, 1:4))]))
    @test Chains.labels(s) == l
    l = mapreduce(prod, +, named_dims)
    @test length(s) == l
    v = 1:l
    m = reshape(v, 1, :)
    @test v[s.a] == 1
    @test v[s.b] == reshape(2:3, 1, 2)
    @test v[s.c] == reshape(4:l, 2, 3, 4)
    @test v[s] == (a = v[s.a], b = v[s.b], c = v[s.c])
    @test m[:, s.a] == [1]
    @test m[:, s.b] == [reshape(2:3, 1, 2)]
    @test m[:, s.c] == [reshape(4:l, 2, 3, 4)]
    @test m[1, s] == (a = m[1, s.a], b = m[1, s.b], c = m[1, s.c])
    @test m[:, s] == [(a = m[1, s.a], b = m[1, s.b], c = m[1, s.c])]
end

@testset "chains" begin
    sample = Float64.(hcat(1:10, 2:2:20, 3:3:30))
    sch = Chains.IndexSchema((a = (), b = (2, )))
    chain = Chains.Chain(sch, sample; warmup = 3, is_ordered = true, thinning = 2)
    @test Chains.labels(sch) == Chains.labels(chain)
    @test Chains.sample_matrix(chain) == sample[4:end, :]
    @test Chains.sample_matrix(chain, Val(true)) == sample
    @test Chains.thinning(chain) == 2
    @test Chains.warmup(chain) == 3
    p = [(f = Float64(i); (a = i, b = [2*i, 3*i])) for i in 4:10]
    @test collect(chain[:, :]) == p
    @test collect(chain[:, :a]) == map(x -> x.a, p)
    @test collect(chain[:, :b]) == map(x -> x.b, p)
    c2 = vcat(chain, chain)
    @test Chains.sample_matrix(c2) == vcat(sample[4:end, :], sample[4:end, :])
    @test Chains.warmup(c2) == 0
    @test Chains.thinning(c2) == nothing
end

####
#### StanCSV
####

@testset "variable name parsing" begin
    @test parse_variable_name("a") == (:a => ())
    @test parse_variable_name("kappa.1.3") == (:kappa => (1, 3))
    @test parse_variable_name("stepsize__") == (:stepsize__ => ())
    @test_throws ArgumentError parse_variable_name("b.")
    @test_throws ArgumentError parse_variable_name("b.12.")
    @test_throws ArgumentError parse_variable_name("")
end

@testset "collapsing" begin
    v = [:a => (1, 1), :a => (2, 1), :a => (1, 2), :a => (2, 2)]
    @test collapse_contiguous_dimensions(v, 1) == (:a, (2, 2), 5)
    @test_throws ArgumentError collapse_contiguous_dimensions(v, 2)
    @test_throws ArgumentError collapse_contiguous_dimensions(v[1:(end-1)], 1)
    @test_throws ArgumentError collapse_contiguous_dimensions(v[[1, 2, 4]], 1)
end

"""
$(SIGNATURES)

Function that generates a vector of `name => index` pairs. For testing.
"""
function make_parsed_header(named_dims)
    v = Any[]
    for (name, dims) in pairs(named_dims)
        append!(v, collect(name => Tuple(ix) for ix in CartesianIndices(dims)))
    end
    v
end

@testset "parsing schema" begin
    nd = (a = (), b = (1, 2), c = (3, 4, 7), d = ())
    @test parse_schema(make_parsed_header(nd)) == Chains.IndexSchema(nd)
end

@testset "reading CSV data" begin

    # Print `n` comment lines to `io`.
    function comment!(io::IO, n = 1)
        for _ in 1:n
            println(io, "# comment")
        end
    end

    # generate test data
    io = IOBuffer()
    comment!(io, rand(1:5))
    println(io, "a, b.1, b.2, c.1.1, c.2.1, c.1.2, c.2.2")
    comment!(io, rand(1:3))
    for i in 1:10
        f = Float64(i)
        print(io, f, ", ", f + 1, " ,", f + 2) # space deliberately mix in spaces
        for j in 1:4
            print(io, ",",  i + j + 4.0) # c[:, :]
        end
        println(io)
        rand() < 0.5 && comment!(io) # intersperse with comments
    end
    comment!(io, rand(1:4))
    contents = String(take!(io))

    # read test data
    chain = StanCSV.read_chain(IOBuffer(contents))
    sch = Chains.schema(chain)
    @test sch == Chains.IndexSchema((a = (), b = (2, ), c = (2, 2)))
    @test chain[:, :a] == Float64.(1:10)
    @test chain[:, :b] == [[i + 1.0, i + 2.0] for i in 1:10]
    @test chain[:, :c] == [(i + 4.0) .+ [1 3; 2 4] for i in 1:10]
end

@testset "corner cases" begin
    @test_throws(ArgumentError("Fewer than 2 fields in line."),
                 StanCSV.read_chain(IOBuffer("a,b\n1\n2\n")))
    c1 = StanCSV.read_chain(IOBuffer("a,b\n1,2\n"))
    c2 = Chains.Chain(Chains.IndexSchema((a = (), b = ())), collect(reshape(1.0:2.0, 1, 2)))
    @test Chains.schema(c1) == Chains.schema(c2)
    @test Chains.sample_matrix(c1) == Chains.sample_matrix(c2)
end

@testset "filename discovery" begin
    mktempdir() do dir
        base = "foo_"
        id_files = [i => joinpath(dir, base * string(i) * ".csv") for i in 1:10]
        for (_, file) in id_files
            touch(file)
        end
        @test StanCSV.matching_files(joinpath(dir, base)) == id_files
    end
end


text:
function FieldAuxiliaryElements(arg0::FieldOrbit, arg1::jint)
    return FieldAuxiliaryElements((FieldOrbit, jint), arg0, arg1)
end

function get_alpha(obj::FieldAuxiliaryElements)
    return jcall(obj, "getAlpha", RealFieldElement, ())
end

function get_b(obj::FieldAuxiliaryElements)
    return jcall(obj, "getB", RealFieldElement, ())
end

function get_beta(obj::FieldAuxiliaryElements)
    return jcall(obj, "getBeta", RealFieldElement, ())
end

function get_c(obj::FieldAuxiliaryElements)
    return jcall(obj, "getC", RealFieldElement, ())
end

function get_date(obj::FieldAuxiliaryElements)
    return jcall(obj, "getDate", FieldAbsoluteDate, ())
end

function get_ecc(obj::FieldAuxiliaryElements)
    return jcall(obj, "getEcc", RealFieldElement, ())
end

function get_frame(obj::FieldAuxiliaryElements)
    return jcall(obj, "getFrame", Frame, ())
end

function get_gamma(obj::FieldAuxiliaryElements)
    return jcall(obj, "getGamma", RealFieldElement, ())
end

function get_h(obj::FieldAuxiliaryElements)
    return jcall(obj, "getH", RealFieldElement, ())
end

function get_k(obj::FieldAuxiliaryElements)
    return jcall(obj, "getK", RealFieldElement, ())
end

function get_keplerian_period(obj::FieldAuxiliaryElements)
    return jcall(obj, "getKeplerianPeriod", RealFieldElement, ())
end

function get_le(obj::FieldAuxiliaryElements)
    return jcall(obj, "getLe", RealFieldElement, ())
end

function get_lm(obj::FieldAuxiliaryElements)
    return jcall(obj, "getLM", RealFieldElement, ())
end

function get_lv(obj::FieldAuxiliaryElements)
    return jcall(obj, "getLv", RealFieldElement, ())
end

function get_mean_motion(obj::FieldAuxiliaryElements)
    return jcall(obj, "getMeanMotion", RealFieldElement, ())
end

function get_p(obj::FieldAuxiliaryElements)
    return jcall(obj, "getP", RealFieldElement, ())
end

function get_q(obj::FieldAuxiliaryElements)
    return jcall(obj, "getQ", RealFieldElement, ())
end

function get_retrograde_factor(obj::FieldAuxiliaryElements)
    return jcall(obj, "getRetrogradeFactor", jint, ())
end

function get_sma(obj::FieldAuxiliaryElements)
    return jcall(obj, "getSma", RealFieldElement, ())
end

function get_vector_f(obj::FieldAuxiliaryElements)
    return jcall(obj, "getVectorF", FieldVector3D, ())
end

function get_vector_g(obj::FieldAuxiliaryElements)
    return jcall(obj, "getVectorG", FieldVector3D, ())
end

function get_vector_w(obj::FieldAuxiliaryElements)
    return jcall(obj, "getVectorW", FieldVector3D, ())
end


text:
module Fuel

export total, correct_total

function fuel_for_weight(weight)
    fld(weight, 3) - 2
end

function total(weights)
    return sum(fuel_for_weight, weights)
end

function correct_fuel_for_weight(weight)
    fuel = fuel_for_weight(weight)
    fuel <= 0 ? 0 : fuel += correct_fuel_for_weight(fuel)
end

function correct_total(weights)
    return sum(correct_fuel_for_weight, weights)
end

end


text:
"""
The geometry consists of a “hook” in the form of a curved strip rigidly clamped
at one end and loaded with a unit in-plane shear along the width at the other
end. It has two circular segments that are connected at the tangent point. The
smaller segment has a mean radius of 0.3556 m (14 inches) and spans 60° from
the clamped end to the tangent point. The larger segment spans 150° from the
tangent point to the free end and has a mean radius of 1.1684 m (46 inches).
The hook is 0.0508 m (2 inches) thick and 0.508 m (20 inches) wide, modeled as
linear elastic with an elastic modulus of 22.77 MPa (3300 psi) and a Poisson's
ratio of 0.35. In most tests the shear force is applied through the use of a
distributing coupling constraint. The coupling constraint provides coupling
between a reference node on which the load is prescribed and the nodes located
on the free end. The distributed nodal loads on the free end are equivalent to
a uniformly distributed load of 8.7563 N/m (0.05 lb/in).

"""
module raasch_examples

using LinearAlgebra, Statistics
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FESetShellQ4Module: FESetShellQ4
using FinEtoolsFlexStructures.FEMMShellT3FFModulew
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!
using VisualStructures: plot_nodes, plot_midline, render, plot_space_box, plot_midsurface, space_aspectratio, save_to_json
using FinEtools.MeshExportModule.VTKWrite: vtkwrite

function _execute(input = "raasch_s4_1x9.inp", drilling_stiffness_scale = 1.0, visualize = true, nL = 9, nW = 1)
    # The physical quantities are provided in english units, dimensions in inches.10
    E = 3300.0;
    nu = 0.35;
    thickness  =  2.0;
    tolerance = thickness/2
    # Reference solution for the vertical deflection under the load. Obtained with a refined 20-node solid model.
    ref_sol = 5.022012648671993;
    R = 46.0;
    formul = FEMMShellT3FFModule
    # formul = FEMMShellT3DSGMTModule
    @show formul

    if input == ""
        fens, fes = T3block(210.0, 20.0, nL, nW)
        fens.xyz = xyz3(fens)
        for k in 1:count(fens)
            a, z = fens.xyz[k, :]
            if a <= 60.0
                fens.xyz[k, :] .= (14*sin(a*pi/180), 14*(1-cos(a*pi/180)), z)
            else
                a = a - 60.0 + 30.0
                xc = 14*sin(60*pi/180) + 46*cos(30*pi/180)
                yc = 14*(1-cos(60*pi/180)) - 46*sin(30*pi/180)
                fens.xyz[k, :] .= (xc-46*cos(a*pi/180), yc+46*sin(a*pi/180), z)
            end
        end
        input = "raasch_$(nL)x$(nW)"
    else
        output = import_ABAQUS(joinpath(dirname(@__FILE__()), input))
        fens = output["fens"]
        fes = output["fesets"][1]

        connected = findunconnnodes(fens, fes);
        fens, new_numbering = compactnodes(fens, connected);
        fes = renumberconn!(fes, new_numbering);

        fens, fes = Q4toT3(fens, fes)

        # fens, fes = T3refine(fens, fes)# .
        # fens, fes = T3refine(fens, fes)# .
        # fens, fes = T3refine(fens, fes)# .
    end

    @show count(fens), count(fes)

    # plots = cat(plot_space_box([[0 0 -R/2]; [R/2 R/2 R/2]]),
    #     plot_nodes(fens),
    #     plot_midsurface(fens, fes);
    # dims = 1)
    # pl = render(plots)

    mater = MatDeforElastIso(DeforModelRed3D, E, nu)

        # Report

    @info "Mesh: $input"

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = formul.make(IntegDomain(fes, TriRule(1), thickness), mater)
    femm.drilling_stiffness_scale = drilling_stiffness_scale
    stiffness = formul.stiffness
    associategeometry! = formul.associategeometry!

    # Construct the requisite fields, geometry and displacement
    # Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

    # Apply EBC's
    # Clamped end
    l1 = selectnode(fens; box = Float64[0 0 -Inf Inf -Inf Inf], inflate = tolerance)
    for i in [1,2,3,4,5,6]
        setebc!(dchi, l1, true, i)
    end

    applyebc!(dchi)
    numberdofs!(dchi);

    # Assemble the system matrix
    associategeometry!(femm, geom0)
    K = stiffness(femm, geom0, u0, Rfield0, dchi);

    # Load
    bfes = meshboundary(fes)
    l1 = selectelem(fens, bfes, box = [97.96152422706632 97.96152422706632 -16 -16 0 20], inflate = 1.0e-6)
    lfemm = FEMMBase(IntegDomain(subset(bfes, l1), GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[0, 0, 0.05, 0, 0, 0]);
    F = distribloads(lfemm, geom0, dchi, fi, 3);
    @assert  isapprox(sum(F), 1.0)

    # @infiltrate
    # Solve
    U = K\F
    scattersysvec!(dchi, U[:])
    nl = selectnode(fens; box = Float64[97.96152422706632 97.96152422706632 -16 -16 0 20], inflate = 1.0e-6)
    targetu =  mean(dchi.values[nl, 3])
    @info "Solution: $(round(targetu, digits=8)),  $(round(targetu/ref_sol, digits = 4)*100)%"

    # formul._resultant_check(femm, geom0, u0, Rfield0, dchi)

    # Generate a graphical display of displacements and rotations
    scalars = []
    for nc in 1:6
        push!(scalars, ("dchi$nc", deepcopy(dchi.values[:, nc])))
    end
    vectors = []
    push!(vectors, ("U", deepcopy(dchi.values[:, 1:3])))
    push!(vectors, ("UR", deepcopy(dchi.values[:, 4:6])))
    vtkwrite("$input-dchi.vtu", fens, fes; scalars = scalars, vectors = vectors)

    # Generate a graphical display of resultants
    function csys!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        cross3!(view(csmatout, :, 3), view(tangents, :, 1), view(tangents, :, 2))
        r = view(csmatout, :, 3)
        csmatout[:, 3] .= vec(r)/norm(vec(r))
        csmatout[:, 2] .= (0.0, 0.0, 1.0)
        cross3!(view(csmatout, :, 1), view(csmatout, :, 2), view(csmatout, :, 3))
        return csmatout
    end
    ocsys = CSys(3, 3, csys!)
    scalars = []
    for nc in 1:3
        fld = fieldfromintegpoints(femm, geom0, dchi, :moment, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
        push!(scalars, ("m$nc", fld.values))
    end
    vtkwrite("$(input)-m.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
    scalars = []
    for nc in 1:3
        fld = fieldfromintegpoints(femm, geom0, dchi, :membrane, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
        push!(scalars, ("n$nc", fld.values))
    end
    vtkwrite("$(input)-n.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])
    scalars = []
    for nc in 1:2
        fld = fieldfromintegpoints(femm, geom0, dchi, :shear, nc, outputcsys = ocsys)
            # fld = elemfieldfromintegpoints(femm, geom0, dchi, :moment, nc)
        push!(scalars, ("q$nc", fld.values))
    end
    vtkwrite("$(input)-q.vtu", fens, fes; scalars = scalars, vectors = [("u", dchi.values[:, 1:3])])

    # Visualization
    if visualize
        scattersysvec!(dchi, (R/2)/maximum(abs.(U)).*U)
        update_rotation_field!(Rfield0, dchi)
        plots = cat(plot_space_box([[0 0 -R]; [R R R]]),
            plot_nodes(fens),
            plot_midsurface(fens, fes; x = geom0.values, u = dchi.values[:, 1:3], R = Rfield0.values);
            dims = 1)
        pl = render(plots)
    end
    return targetu/ref_sol
end

function test_convergence()

    @info "Raasch hook"

    # for m in ["1x9", "3x18", "5x36", "10x72", "20x144"]
        # _execute("raasch_s4_" * m * ".inp", 1.0, false)
    for n in 1:7
        _execute("", 1.0, false, 9*2^(n-1), 1*2^(n-1))
    end
    return true
end

function test_dep_drilling_stiffness_scale()

    @info "Raasch hook"

    all_results = []
    all_drilling_stiffness_scale = [1000.0, 1.0, 0.1, 0.0001, 0.000001]
    for drilling_stiffness_scale in all_drilling_stiffness_scale
        results = Float64[]
        # for m in ["1x9", "3x18", "5x36", "10x72", "20x144"]
            # v = _execute("raasch_s4_" * m * ".inp", drilling_stiffness_scale, false)
        for n in 1:7
            v = _execute("", drilling_stiffness_scale, false, 9*2^(n-1), 1*2^(n-1))
                    push!(results, v)
        end
        push!(all_results, results)
    end
    return all_drilling_stiffness_scale, all_results
end

function allrun()
    println("#####################################################")
    println("# test_convergence ")
    test_convergence()
    println("#####################################################")
    println("# test_dep_drilling_stiffness_scale ")
    test_dep_drilling_stiffness_scale()
    return true
end # function allrun

@info "All examples may be executed with "
println("using .$(@__MODULE__); $(@__MODULE__).allrun()")

end # module
nothing


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md
using JuliaFEM
using JuliaFEM.Core
using JuliaFEM.API: Model

"""
Function for creating solver and all the necessary components
for the calculation
"""
function get_solver(model::Model, case_name::String)
    case = model.load_cases[case_name]

    # Create core elements
    core_elements = create_core_elements(model)

    # Add Neumann boundary conditions to elements
    add_neumann_bcs!(model, case, core_elements)

    # Create Dirichlet boundary conditions
    dirichlet_arr = create_dirichlet_bcs(model, case, core_elements)

    # Create solver
    solver = create_solver(model, case,
                            core_elements, dirichlet_arr)
	return solver
end

"""
"""
function create_solver(model, case, core_elements, dirichlet_arr)
    field_problem = JuliaFEM.Core.(case.problem)()
    all_elsets = model.elsets

    # Pushing all the defined element sets into the calculation
    for element_set in case.sets
        el_ids = all_elsets[element_set].elements
        for each in el_ids
            push!(field_problem, core_elements[each])
        end
    end

    # Creating the solver and pushing problems and
	solver = JuliaFEM.Core.(case.solver)()
	push!(solver, field_problem)
	push!(solver, dirichlet_arr...)
    return solver
end

"""
"""
function create_dirichlet_bcs(model, case, core_elements)
    dirichlet_arr = Any[]
    dirichlet_bcs = case.dirichlet_boundary_conditions
    elsets = model.elsets
    field_problem = JuliaFEM.Core.(case.problem)()
    for each in dirichlet_bcs
        set_name = each.set_name
        value = each.value
        problem = JuliaFEM.Core.DirichletProblem(
           JuliaFEM.Core.get_unknown_field_name(field_problem),
           JuliaFEM.Core.get_unknown_field_dimension(field_problem))
        set_for_bc = each.set_name
        set_ids = elsets[set_for_bc]
        bc = each.value
        for el_id in set_ids.elements
            core_element = core_elements[el_id]
            core_element[bc[1]] = bc[2]
            push!(problem, core_element)
        end
        push!(dirichlet_arr, problem)
    end
    return dirichlet_arr
end

"""
"""
function create_core_elements(model)
    core_elements = Dict()
    nodes = model.nodes
    all_elements = model.elements
    element_ids = keys(all_elements)
    for el_id in element_ids
        element = all_elements[el_id]
        el_type = element.element_type
        el_id = element.id
        mat = element.material
        conn = element.connectivity
        core_element = JuliaFEM.Core.(el_type)(conn)
        core_element["geometry"] = map(x->nodes[x].coords, conn)
        for each in keys(mat.scalar_data)
            core_element[each] = mat.scalar_data[each]
        end
        core_elements[el_id] = core_element
        model.elements[el_id].results = core_element
    end
    return core_elements
end

"""
"""
function add_neumann_bcs!(model, case, core_elements)
    neumann_bcs = case.neumann_boundary_conditions
    elsets = model.elsets
    for each in neumann_bcs
        set_for_bc = each.set_name
        set_ids = elsets[set_for_bc]
        bc = each.value
        for el_id in set_ids.elements
            core_element = core_elements[el_id]
            core_element[bc[1]] = bc[2]
        end
        if !(set_for_bc in case.sets)
            push!(case.sets, set_for_bc)
        end
    end
end

"""
"""
function solve!(model::Model, case_name::String, time::Float64)
    # Create solver
	solver = get_solver(model, case_name)

    # Solve problem at given time
    solver(time)
end


text:
"Abstract struct for thermal generation technologies"
abstract type
    ThermalGen <: Generator
end

""""
Data Structure for thermal generation technologies.
    The data structure contains all the information for technical and economical modeling.
    The data fields can be filled using named fields or directly.

    Examples

"""
struct ThermalDispatch <: ThermalGen
    name::String
    available::Bool
    bus::Bus
    tech::Union{TechThermal,Nothing}
    econ::Union{EconThermal,Nothing}
end

ThermalDispatch(; name = "init",
                status = false,
                bus = Bus(),
                tech = TechThermal(),
                econ = EconThermal()) = ThermalDispatch(name, status, bus, tech, econ)

""""
Data Structure for thermal generation technologies subjecto to seasonality constraints.
    The data structure contains all the information for technical and economical modeling and an extra field for a time series.
    The data fields can be filled using named fields or directly.

    Examples

"""
struct ThermalGenSeason <: ThermalGen
    name::String
    available::Bool
    bus::Bus
    tech::Union{TechThermal,Nothing}
    econ::Union{EconThermal,Nothing}
    scalingfactor::TimeSeries.TimeArray
end

ThermalGenSeason(; name = "init",
                status = false,
                bus = Bus(),
                tech = TechThermal(),
                econ = EconThermal(),
                scalingfactor = TimeSeries.TimeArray(Dates.today(),ones(1))) = ThermalGenSeason(name, status, bus, tech, econ, scalingfactor)


text:
@testset "Test generated structs" begin
    descriptor_file =
        joinpath(@__DIR__, "..", "src", "descriptors", "power_system_structs.json")
    existing_dir = joinpath(@__DIR__, "..", "src", "models", "generated")
    @test IS.test_generated_structs(descriptor_file, existing_dir)
end


text:
include("adjoint_lsssolve.jl")
"""
Implements the NILSAS algorithm due to Ni and Talnikar 2019.
"Adjoint sensitivity analysis on chaotic dynamical
systems by Non-Intrusive Least Squares Adjoint
Shadowing (NILSAS)"

Below, u -> F(u) is the time-1 primal map and u_trj refers
to a primal trajectory of length n.
Inputs:
    du_trj: dxdxn n-length timeseries of Jacobian matrices
    		along u_trj
    X: dxn n-length timeseries of perturbation vectors.
    		X[:,i] = dF/ds(u_trj[i,:]) which is in the
    		tangent space of u_trj[i+1,:].
    f: dxn n-length timeseries of the neutral CLV. By
    		assumption, the only neutral CLV is the RHS of
    		the primal ode: du/dt = f(u,s).
    		In case of maps, u_n = F(u_{n+1}), f is zeros.
    J:mxn n-length timeseries of m scalar functions J(u)
    		evaluated along u_trj.
    dJ: mxdxn n-length timeseries of DJ(u) evaluated along
    		u_trj
    s: array of parameters
    d_u: integer dimension of the unstable subspace. You can 		set a guess value to obtain correct sensitivities
    	as long as the guess value is greater than or equal
    	to the true d_u.
Outputs:
    vsh: dxn n-length timeseries of shadowing tangent
    	vectors that are bounded solutions of the
    	conventional tangent equation:
    		v[:,i+1] = du_trj[:,i]*v[:,i] + X[:,i]

    dJds: scalar-valued sensitivity of interest:
    	d/ds sum(J) ≈ d/ds <J, mu>, where mu is
    	the SRB measure.
"""
function lss(du_trj, X, f, s, d_u, g)
    d, n = size(X)
    lyap_exps = zeros(d_u)
    R = zeros(d_u,d_u,n)
    Q = zeros(d,d_u,n)
    v = zeros(d,1,n) #assume one parameter
    v[:,1,1] = zeros(d)
    A = qr!(rand(d,d_u))
    Q[:,:,1] = Array(A.Q)
    R[:,:,1] = A.R

    b = zeros(d_u,1,n) #1 is the # parameters
    pf_Q = zeros(d_u,1,n) #1 is the # parameters
    pf_v = zeros(n)

    ff = zeros(n)
    if all(isapprox.(f,0.))
    	ff = ones(n)
    else
    	[ff[i] = sum(x -> x^2, f[:,i]) for i = 1:n]
    end
    pf_Q[:,1,1] = Q[:,:,1]'*g[:,1]
    [Q[:,j,1] = Q[:,j,1] - (pf_Q[j,1,1]*
     		f[:,1]/ff[1]) for j=1:d_u]

    for i=2:n
        v[:,:,i] = du_trj[:,:,i-1]*v[:,:,i-1] +
    				X[:,i-1]
    	Q[:,:,i] = du_trj[:,:,i-1]*Q[:,:,i-1]
    	pf_Q[:,:,i] = Q[:,:,i]'*g[:,i]
    	#[Q[:,j,i] = Q[:,j,i] - pf_Q[j,1,i]*f[:,i]/
    	#				ff[i] for j=1:d_u]
    	A = qr!(Q[:,:,i])
        Q[:,:,i] = Array(A.Q)
        R[:,:,i] = A.R
    	pf_v[i] = dot(v[:,1,i],g[:,i])
    	v[:,:,i] = v[:,:,i] - (pf_v[i]/ff[i]*f[:,i])
        b[:,:,i] = (Q[:,:,i]')*v[:,:,i]
    	v[:,:,i] = v[:,:,i] - Q[:,:,i]*b[:,:,i]
	lyap_exps .+= log.(abs.(diag(R[:,:,i])))./n
    end

    b = reshape(collect(b),d_u, n)
    a = lsssolve(R,b,sum(pf_v),pf_Q[:,1,:])
    v = reshape(collect(v),d,n)
    vsh = zeros(d, n)
    xi = zeros(n)
    vsh[:,1] = Q[:,:,1]*reshape(a[:,1], d_u, 1) +
    			  v[:,1]
    for i = 2:n
    	vsh[:,i] = v[:,i] + Q[:,:,i]*reshape(a[:,i],d_u,1)
    end
    return vsh, xi
end
function compute_sens(vsh, xi, dJ, f)
    m, d, n = size(dJ)
    dJds = zeros(m)
    for i = 1:n
    	dJds .= dJds .+ dJ[:,:,i]*vsh[:,i]
    end
    return dJds
end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
struct Pendulum{I, T} <: Model{I, T}
    n::Int
    m::Int
    d::Int

    mass    # mass
    b       # friction
    lc      # length to center of mass
    g       # gravity
end

function f(model::Pendulum, x, u, w)
    @SVector [x[2],
              (u[1] / ((model.mass * model.lc * model.lc))
                - model.g * sin(x[1]) / model.lc
                - model.b * x[2] / (model.mass * model.lc * model.lc))]
end

function k(model::Pendulum, x)
    @SVector [model.lc * sin(x[1]),
              -1.0 * model.lc * cos(x[1])]
end

n, m, d = 2, 1, 2
model = Pendulum{Midpoint, FixedTime}(n, m, d, 1.0, 0.1, 0.5, 9.81)


text:
"""
Defintions common to all implementations of `Transport`.
"""
module Common

export Message,
  Transport,
  Connection,
  Connections,
  connectOne,
  connectAll,
  Messenger,
  Listener,
  connectTo,
  messengerTo,
  listenOn,
  sendTo,
  receivedMessages,
  sendMessage,
  connection,
  listener,
  call,
  gcall,
  callAny

using Base
using Logging
using Random

using ...Utils

"""
An abstract payload for sending data via a `Transport`
"""
abstract type Message end

"""
A `Transport` delivers messages sent to groups or individual recipients
"""
abstract type Transport end

"""
A `Connection` represents a conduit between a source and a destination: messages may be
sent to the destination from the source, and optionally they may also be received
from the destination to the source. Closing a connection prevents further message exchange
using the connection.
"""
abstract type Connection end

"""
A listener is typically a task that is awaiting incoming connections; closing a listener will
prevent further connections to that listener.
"""
Listener = Task

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

"""
Create a connection (if one is possible) that enables sending and receiving
messages to the recipient. The provided handler should take 1 argument,
a `Connection` for sending or receiving messages with the recipient.
"""
function connectTo(transport::Transport, recipient)::Connection end

"""
Listen for incoming connections and invoke the indicated function when they appear.
Returns a `Listener`. The provided handler should take 1 argument, a `Connection`;
a new handler will be created for each connection accepted by the listener. A listener
"""
function listenOn(handler::Function, transport::Transport, address)::Listener end

"""
Send a message on the indicated connection
"""
function sendTo(connection::Connection, message) end

"""
Receive a message over the indicated connection
"""
function receiveFrom(connection::Connection) end

"""
Close the connection, freeing resources associated with it--and stopping any tasks
moving data on the connection.
"""
function Base.close(connection::Connection) end

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

"""
A messenger handles asynchronous transmission / reception of messages over a `Connection`,
until closed. In addition to stopping the tasks for transmission / reception, closing
the messenger also closes the underlying connection.
"""
struct Messenger{T}
  resource::T
  inbound::Channel
  outbound::Channel
  errors::Channel
end

"""
Construct a messenger for the given resource
"""
function Messenger(connection)
  errors = Channel()
  messenger = Messenger(
    connection,
    Channel() do inbound
      messageReceiver(connection, inbound, errors)
    end,
    Channel() do outbound
      messageSender(connection, outbound, errors)
    end,
    errors,
  )
  messenger
end

"""
Continuously receive messages over the connection and places them in the
messenger's inbound channel, until the messenger is closed.
"""
function messageReceiver(connection::Connection, inbound::Channel, errors::Channel)
  try
    @debug "Preparing to receive messages"
    while true
      # Receive using the connection's transport
      @debug "Receiving message over transport connection: $connection"
      msg = receiveFrom(connection)
      put!(inbound, msg)
      @debug "Message $msg received over transport connection: $connection"
    end
  catch ex
    if isa(ex, InvalidStateException) && (ex.state == :closed)
      @debug "Inbound closed"
    else
      @error "Error receiving messages" exception = (ex, stacktrace(catch_backtrace()))
    end
  finally
    @debug "Finished receive messages"
  end
end

"""
Continously send any messages on the messenger's outbound channel,
until the messenger is closed
"""
function messageSender(connection::Connection, outbound::Channel, errors::Channel)
  try
    @debug "Preparing to deliver messages"
    for msg in outbound
      @debug "Delivering message $msg over transport connection"
      sendTo(connection, msg)
      @debug "Message $msg delivered"
    end
  catch ex
    if (isa(ex, InvalidStateException) && (ex.state == :closed)) || isa(ex, Base.IOError)
      @debug "Outbounbd closed"
    else
      @error "Error delivering messages" exception = (ex, stacktrace(catch_backtrace()))
    end
  finally
    @debug "Finished delivering messages"
  end
end

"""
Enqueue message for sending over the transport connection
"""
function sendMessage(messenger::Messenger, message)
  @debug "Enqueuing message to send"
  try
    put!(messenger.outbound, message)
  catch ex
    if (isa(ex, InvalidStateException) && (ex.state == :closed))
      @debug "Sending channel closed"
    else
      rethrow()
    end
  end
  @debug "Message enqueued"
end

"""
Return a channel for processing messages received on the Messenger
"""
function receivedMessages(messenger::Messenger)
  messenger.inbound
end

function Base.close(connection::Messenger)
  close(connection.resource)
  close(connection.outbound)
  close(connection.inbound)
  close(connection.errors)
end

"""
Connect to the recipient via the underlying transport and construct
a `Messenger` with the resulting `Connection`
"""
function messengerTo(transport, recipient)
  Messenger(connectTo(transport, recipient))
end

"""
Create a `Connection` to the recepient and a `Messenger` to exchange messages,
then invoke the handler with the `Messenger` as an argument. Close the messenger
as soon as the connection exits.
"""
function connection(handler::Function, transport::Transport, recipient)
  @debug "Beginning connection to $recipient"
  try
    messenger = messengerTo(transport, recipient)
    try
      handler(messenger)
    finally
      close(messenger)
    end
  finally
    @debug "Finished connection to $recipient"
  end
end

struct Connections
  transport::Transport
  messengers::Dict{Any,Messenger}
end

function Base.getindex(connections::Connections, address)
  get(connections.messengers, address, messengerTo(connections.transport, address))
end

Connections(transport::Transport) = Connections(transport, Dict{Any,Messenger}())

"""
Create a connection to the address and if successful retain it in `connections`
"""
function connectOne(connections::Connections, address, timeout)
  connections[address] = connectTo(connections.transport, address)
end

"""
Asynchronously connect to as many destination addresses as possible,
ignoring exceptions if encountered
"""
function connectAll(connections::Connections, addresses, timeout)
  @sync bounded(timeout) do
    for address in addresses
      @async bounded(timeout) do
        if !haskey(connections, address)
          try
            connectOne(connections, address, timeout)
          catch ex
            @error "Error connecting to all" exception = (ex, stacktrace(catch_backtrace()))
          end
        end
      end
    end
  end
  connections
end

"""
Create a listener on the underlying transport at the indicated address,
invoking the handler for each `Connection` discveroed. A listner should be `close`d
when no longer needed, to free up resources
"""
function listener(handler::Function, transport::Transport, address)
  @async begin
    @debug "Beginning to listen"
    listenOn(transport, address) do connection
      @debug "Beginning to handle a connection"
      messenger = Messenger(connection)
      try
        handler(messenger)
      catch ex
        if isa(ex, ListenerClosedException)
          @debug "Closing listener"
        else
          @error "Error while handling" exception = (ex, stacktrace(catch_backtrace()))
        end
      finally
        close(messenger)
        @debug "Finished handling a connection"
      end
    end
    @debug "Finished listening"
  end
end

"""
Exception sent to a listener to indicate that it has been closed
"""
struct ListenerClosedException <: Exception end

function Base.close(listener::Listener)
  schedule(listener, ListenerClosedException(), error = true)
end

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

"""
Send a message through the specified `Messenger`, and wait up to `timeout`
seconds for a response. If no response received in the alloted time, then
throw a `TimeoutException`; otherwise, return the reeived response.
"""
function call(messenger::Messenger, timeout, message)
  bounded(timeout) do
    sendMessage(messenger, message)
    response = take!(receivedMessages(messenger))
    return response
  end
end

"""
Send the message through each of the `Messenger`s, but only allow
`timeout` seconds for a response. Return all results received within
the alloted time, if any.
"""
function gcall(messengers, timeout, message)
  sz = length(messengers)
  responses = Channel()
  results = Vector()
  try
    bounded(timeout) do
      for messenger in messengers
        @async begin
          response = call(messenger, timeout, message)
          put!(responses, response)
        end
      end
      for response in responses
        push!(results, response)
        if length(results) >= sz
          break
        end
      end
    end
  catch ex
    if !isa(ex, TimeoutException)
      rethrow()
    end
  finally
    Base.close(responses)
  end
  results
end

"""
Exception thrown when there is no server that can respond to a client request
"""
struct NoAvailableServerException <: Exception end

"""
Call the messengers in order with the message, returning after the first
messenger that responds without an exception. No other messengers after the first
successful call will receive the message.
"""
function callAny(
  transport::Transport,
  messenger::Union{Messenger,Nothing},
  addresses,
  timeout,
  msg::Message,
)
  for address in shuffle(addresses)
    if messenger === nothing
      try
        messenger = connectTo(transport, address)
      catch ex
        @error "Error connecting to server" exception = (ex, stacktrace(catch_backtrace()))
        # let's try the next address
        continue
      end
    end
    try
      return (messenger, call(messenger, timeout, msg))
    catch ex
      @error "Error calling server" exception = (ex, stacktrace(catch_backtrace()))
      close(messenger)
      messenger = nothing
    end
  end
  # if we wind up here, we couldn't contact any server reliably
  throw(NoAvailableServerException())
end

end


text:
using Base.Test

@test istrue(true)
@test istrue(false) == false
@test istrue(NA) == false

@test isfalse(false)
@test isfalse(true) == false
@test isfalse(NA) == false

@test isna(true) == false
@test isna(false) == false
@test isna(NA) == true