text:
import Krylov.KRYLOV_SOLVERS

function init_Ksolver(M, v, sp::SolverParams)
  kmethod = sp.kmethod
  if kmethod ∈ (:gpmr, :diom, :fom, :dqgmres, :gmres)
    return eval(KRYLOV_SOLVERS[kmethod])(M, v, sp.mem)
  end
  return eval(KRYLOV_SOLVERS[kmethod])(M, v)
end

ksolve!(
  KS::MinresSolver{T},
  K,
  rhs::AbstractVector{T},
  M;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T} = minres!(
  KS,
  K,
  rhs,
  M = M,
  verbose = verbose,
  atol = zero(T),
  rtol = zero(T),
  ratol = atol,
  rrtol = rtol,
)

ksolve!(
  KS::MinresQlpSolver{T},
  K,
  rhs::AbstractVector{T},
  M;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T} = minres_qlp!(KS, K, rhs, M = M, verbose = verbose, atol = atol, rtol = rtol)

ksolve!(
  KS::SymmlqSolver{T},
  K,
  rhs::AbstractVector{T},
  M;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T} = symmlq!(KS, K, rhs, M = M, verbose = verbose, atol = atol, rtol = rtol)

ksolve!(
  KS::CgSolver{T},
  K,
  rhs::AbstractVector{T},
  M;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T} = cg!(KS, K, rhs, M = M, verbose = verbose, atol = atol, rtol = rtol)

ksolve!(
  KS::CgLanczosSolver{T},
  K,
  rhs::AbstractVector{T},
  M;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T} = cg_lanczos!(KS, K, rhs, M = M, verbose = verbose, atol = atol, rtol = rtol)

ksolve!(
  KS::CrSolver{T},
  K,
  rhs::AbstractVector{T},
  M;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T} = cr!(KS, K, rhs, M = M, verbose = verbose, atol = atol, rtol = rtol)

ksolve!(
  KS::BilqSolver{T},
  K,
  rhs::AbstractVector{T},
  M;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T} = bilq!(KS, K, rhs, verbose = verbose, atol = atol, rtol = rtol)

ksolve!(
  KS::QmrSolver{T},
  K,
  rhs::AbstractVector{T},
  M;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T} = qmr!(KS, K, rhs, verbose = verbose, atol = atol, rtol = rtol)

ksolve!(
  KS::UsymlqSolver{T},
  K,
  rhs::AbstractVector{T},
  M;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T} = usymlq!(KS, K, rhs, rhs, verbose = verbose, atol = atol, rtol = rtol)

ksolve!(
  KS::UsymqrSolver{T},
  K,
  rhs::AbstractVector{T},
  M;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T} = usymqr!(KS, K, rhs, rhs, verbose = verbose, atol = atol, rtol = rtol)

ksolve!(
  KS::BicgstabSolver{T},
  K,
  rhs::AbstractVector{T},
  M;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T} = bicgstab!(KS, K, rhs, verbose = verbose, atol = atol, rtol = rtol)

ksolve!(
  KS::DiomSolver{T},
  K,
  rhs::AbstractVector{T},
  M;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T} = diom!(KS, K, rhs, verbose = verbose, atol = atol, rtol = rtol)

ksolve!(
  KS::FomSolver{T},
  K,
  rhs::AbstractVector{T},
  M;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T} = fom!(KS, K, rhs, verbose = verbose, atol = atol, rtol = rtol)

ksolve!(
  KS::DqgmresSolver{T},
  K,
  rhs::AbstractVector{T},
  M;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T} = dqgmres!(KS, K, rhs, verbose = verbose, atol = atol, rtol = rtol)

ksolve!(
  KS::GmresSolver{T},
  K,
  rhs::AbstractVector{T},
  M;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T} = gmres!(KS, K, rhs, verbose = verbose, atol = atol, rtol = rtol)

ksolve!(
  KS::TricgSolver{T},
  A,
  ξ1::AbstractVector{T},
  ξ2::AbstractVector{T},
  M,
  N;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
  gsp::Bool = false,
) where {T} =
  tricg!(KS, A, ξ1, ξ2, M = M, N = N, flip = true, verbose = verbose, atol = atol, rtol = rtol)

ksolve!(
  KS::TrimrSolver{T},
  A,
  ξ1::AbstractVector{T},
  ξ2::AbstractVector{T},
  M,
  N;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
  gsp::Bool = false,
) where {T} = trimr!(
  KS,
  A,
  ξ1,
  ξ2,
  M = M,
  N = (gsp == true ? I : N),
  τ = -one(T),
  ν = (gsp ? zero(T) : one(T)),
  verbose = verbose,
  atol = atol,
  rtol = rtol,
)

function ksolve!(
  KS::GpmrSolver{T},
  A,
  ξ1::AbstractVector{T},
  ξ2::AbstractVector{T},
  M,
  N;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
  gsp::Bool = false,
) where {T}
  sqrtδI = sqrt(N.λ) * I
  return gpmr!(
    KS,
    A,
    A',
    ξ1,
    ξ2,
    C = sqrt.(M),
    D = gsp ? I : sqrtδI,
    E = sqrt.(M),
    F = gsp ? I : sqrtδI,
    λ = -one(T),
    μ = gsp ? zero(T) : one(T),
    verbose = verbose,
    atol = atol,
    rtol = rtol,
  )
end

# gpmr solver for K3.5
function ksolve!(
  KS::GpmrSolver{T},
  A,
  ξ1::AbstractVector{T},
  ξ2::AbstractVector{T},
  M,
  N::AbstractLinearOperator{T};
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T}
  return gpmr!(
    KS,
    A,
    A',
    ξ1,
    ξ2,
    C = M,
    D = N,
    E = transpose(M),
    F = N,
    λ = -one(T),
    verbose = verbose,
    atol = atol,
    rtol = rtol,
  )
end

function ksolve!(
  KS::LslqSolver{T},
  A,
  ξ1::AbstractVector{T},
  M,
  δ::T;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T}
  return lslq!(
    KS,
    A,
    ξ1,
    M = M,
    N = δ > zero(T) ? (one(T) / δ) * I : I,
    verbose = verbose,
    atol = atol,
    btol = rtol,
    etol = zero(T),
    utol = zero(T),
    conlim = T(Inf),
    sqd = δ > zero(T),
  )
end

function ksolve!(
  KS::LsqrSolver{T},
  A,
  ξ1::AbstractVector{T},
  M,
  δ::T;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T}
  return lsqr!(
    KS,
    A,
    ξ1,
    M = M,
    N = δ > zero(T) ? (one(T) / δ) * I : I,
    verbose = verbose,
    axtol = atol,
    btol = rtol,
    # atol = atol,
    # rtol = rtol,
    etol = zero(T),
    conlim = T(Inf),
    sqd = δ > zero(T),
  )
end

function ksolve!(
  KS::LsmrSolver{T},
  A,
  ξ1::AbstractVector{T},
  M,
  δ::T;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T}
  return lsmr!(
    KS,
    A,
    ξ1,
    M = M,
    N = δ > zero(T) ? (one(T) / δ) * I : I,
    verbose = verbose,
    axtol = zero(T), # atol,
    btol = zero(T), # rtol,
    atol = atol,
    rtol = rtol,
    etol = zero(T),
    conlim = T(Inf),
    sqd = δ > zero(T),
  )
end

function ksolve!(
  KS::LnlqSolver{T},
  A,
  ξ2::AbstractVector{T},
  M,
  δ::T;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T}
  return lnlq!(
    KS,
    A,
    ξ2,
    N = M,
    M = δ > zero(T) ? (one(T) / δ) * I : I,
    verbose = verbose,
    atol = atol,
    rtol = rtol,
    sqd = δ > zero(T),
  )
end

function ksolve!(
  KS::CraigSolver{T},
  A,
  ξ2::AbstractVector{T},
  M,
  δ::T;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T}
  return craig!(
    KS,
    A,
    ξ2,
    N = M,
    M = δ > zero(T) ? (one(T) / δ) * I : I,
    λ = δ > zero(T) ? one(T) : zero(T),
    verbose = verbose,
    atol = atol,
    rtol = rtol,
    btol = zero(T),
    conlim = T(Inf),
  )
end

function ksolve!(
  KS::CraigmrSolver{T},
  A,
  ξ2::AbstractVector{T},
  M,
  δ::T;
  verbose::Integer = 0,
  atol::T = T(sqrt(eps(T))),
  rtol::T = T(sqrt(eps(T))),
) where {T}
  return craigmr!(
    KS,
    A,
    ξ2,
    N = M,
    M = δ > zero(T) ? (one(T) / δ) * I : I,
    verbose = verbose,
    atol = atol,
    rtol = rtol,
    sqd = true,
  )
end

function kscale!(rhs::AbstractVector{T}) where {T <: Real}
  rhsNorm = norm(rhs)
  if rhsNorm != zero(T)
    rhs ./= rhsNorm
  end
  return rhsNorm
end

function kunscale!(sol::AbstractVector{T}, rhsNorm::T) where {T <: Real}
  if rhsNorm != zero(T)
    sol .*= rhsNorm
  end
end

function update_kresiduals_history!(
  res::AbstractResiduals{T},
  K,
  sol::AbstractVector{T},
  rhs::AbstractVector{T},
) where {T <: Real}
  if typeof(res) <: ResidualsHistory
    mul!(res.KΔxy, K, sol) # krylov residuals
    res.Kres .= res.KΔxy .- rhs
  end
end

function update_kresiduals_history_K1struct!(
  res::AbstractResiduals{T},
  A,
  E,
  tmp,
  δ,
  sol::AbstractVector{T},
  rhs::AbstractVector{T},
  formul::Symbol,
) where {T <: Real}
  if typeof(res) <: ResidualsHistory
    mul!(tmp, A', sol)
    tmp ./= E
    @views mul!(res.KΔxy, A, tmp)
    if formul == :K1_1
      mul!(res.Kres, A, rhs ./ E, -one(T), zero(T))
    elseif formul == :K1_2
      res.Kres .= .-rhs
    end
    res.Kres .+= res.KΔxy .+ δ .* sol # residual computation
  end
end


text:
# make CompoundPeriod iterable

eltype(x::CompoundPeriod) = Period
length(x::CompoundPeriod) = length(x.periods)

eltype(x::ReverseCompoundPeriod) = Period
length(x::ReverseCompoundPeriod) = length(x.cperiod.periods)

function Base.iterate(x::CompoundPeriod)
    iszero(length(x)) && return nothing
    return (x.periods[1], 1)
end

function Base.iterate(x::CompoundPeriod, state::Int)
    state === length(x) && return nothing
    state += 1
    return (x.periods[state], state)
end

function Base.iterate(x::ReverseCompoundPeriod)
    n = length(x)
    iszero(n) && return nothing
    return (x.cperiod.periods[n], n)
end

function Base.iterate(x::ReverseCompoundPeriod, state::Int)
    state === 1 && return nothing
    state -= 1
    return (x.cperiod.periods[state], state)
end

# make Period iterable for interoperability with CompoundPeriod

eltype(x::Period) = typeof(x)
length(x::Period) = 1

Base.iterate(x::Period) = (x, 1)
Base.iterate(x::Period, state::Int) = nothing

ext: jl


text:
using Documenter, BioCore

makedocs(
    format = :html,
    sitename = "BioCore.jl",
    pages = [
        "Home" => "index.md",
        "Using file format specimens" => "testing.md"
    ],
    authors = "The BioJulia Organisation and other contributors."
)

deploydocs(
    repo = "github.com/BioJulia/BioCore.jl.git",
    julia = "1.0",
    osname = "linux",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
using Symbolics, SparseArrays, LinearAlgebra, BenchmarkTools

@Symbolics.variables x[1:31]

S =[x[1]  0     0     0     0     0     0     0     0     0     0     0     0     0
    0     x[2]  0     0     0     0     0     0     0     0     0     0     0     0
    0     0     x[3]  0     0     0     0     0     0     0     0     0     0     0
    0     0     0     x[4]  0     0     0     x[15] 0     0     0     0     0     0
    0     0     0     0     x[5]  0     x[16] 0     x[17] 0     0     0     0     0
    0     0     0     0     0     x[6]  0     0     0     x[18] 0     0     0     0
    0     0     0     0     x[16] 0     x[7]  0     0     0     x[19] 0     0     0
    0     0     0     x[15] 0     0     0     x[8]  0     0     0     0     0     0
    0     0     0     0     x[17] 0     0     0     x[9]  0     x[20] 0     0     0
    0     0     0     0     0     x[18] 0     0     0     x[10] 0     0     0     0
    0     0     0     0     0     0     x[19] 0     x[20] 0     x[11] 0     0     0
    0     0     0     0     0     0     0     0     0     0     0     x[12] 0     0
    0     0     0     0     0     0     0     0     0     0     0     0     x[13] 0
    0     0     0     0     0     0     0     0     0     0     0     0     0     x[14]]

F =[1.0 0 0 0 0 0 0 0 0 0 0 0 0 0
    0 1 0 0 0 0 0 0 0 0 0 0 0 0
    0 0 1 0 0 0 0 0 0 0 0 0 0 0
    0 0 0 1 0 0 0 0 0 0 0 0 0 0
    0 0 0 0 1 0 0 0 0 0 0 0 0 0
    0 0 0 0 0 1 0 0 0 0 0 0 0 0
    0 0 0 0 0 0 1 0 0 0 0 0 0 0
    0 0 0 0 0 0 0 1 0 0 0 0 0 0
    0 0 0 0 0 0 0 0 1 0 0 0 0 0
    0 0 0 0 0 0 0 0 0 1 0 0 0 0
    0 0 0 0 0 0 0 0 0 0 1 0 0 0]

A =[0  0  0  0  0  0  0  0  0  0  0     1     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[21] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[22] 0     0
    0  0  0  0  0  0  0  0  0  0  0     0     1     0
    0  0  0  0  0  0  0  0  0  0  0     0     x[23] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[24] 0
    0  0  0  0  0  0  0  0  0  0  0     0     x[25] 0
    0  0  0  0  0  0  0  0  0  0  0     0     0     1
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[26]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[27]
    0  0  0  0  0  0  0  0  0  0  0     0     0     x[28]
    0  0  0  0  0  0  0  0  0  0  0     0     0     0
    0  0  0  0  0  0  0  0  0  0  0     x[29] 0     0
    0  0  0  0  0  0  0  0  0  0  0     x[30] x[31] 0]

S = sparse(S)

#F
F = sparse(F)

#A
A = sparse(A)

par = [x[i] for i = 1:31]
invia = I + A + A^2 + A^3 + A^4
imp_cov_sym = F*invia*S*permutedims(invia)*permutedims(F)
imp_cov_sym = Array(imp_cov_sym)
imp_cov_sym = Symbolics.simplify.(imp_cov_sym)
imp_cov_sym = imp_cov_sym[tril(trues(size(F, 1), size(F, 1)))]

∇Σ_sym = Symbolics.sparsejacobian(imp_cov_sym, par)
∇²Σ_sym = Symbolics.sparsejacobian(vec(permutedims(∇Σ_sym)), par)
∇Σ_sym = Array(∇Σ_sym)
# H_array = [Symbolics.sparsejacobian(∇Σ_sym[i, :], par) for i = 1:66]
H_array = [Symbolics.sparsehessian(imp_cov_sym[i], par) for i = 1:66]

#= @variables x[1:31]
par = [x[i] for i = 1:31]
testarray = [x[1]; x[1]*x[2]; x[5]*x[9]]

J = Symbolics.sparsejacobian(testarray, par)
J
H = Symbolics.sparsejacobian(J[3, :], par)
H_array[2].nzval =#

#= function similar_sparse_float(S)
    S_tuple = findnz(S)
    S_tuple_new = (S_tuple[1], S_tuple[2], ones(size(S_tuple[3], 1)))
    new = sparse(S_tuple_new..., size(S)...)
end

H_pre = similar_sparse_float.(H_array) =#

nfact = 5
nitem = 30

## Model definition
nobs = nfact*nitem
nnod = nfact+nobs
@variables x[1:Int64(nobs + nobs)]

x = [x[i] for i = 1:size(x, 1)]
#F
Ind = collect(1:nobs)
Jnd = collect(1:nobs)
V = fill(1,nobs)
F = sparse(Ind, Jnd, V, nobs, nnod)

#S
Ind = collect(1:nnod)
Jnd = collect(1:nnod)
V = [x[1:nobs]; fill(1.0, nfact)]
S = sparse(Ind, Jnd, V, nnod, nnod)

#A
Ind = collect(1:nobs)
Jnd = vcat([fill(nobs+i, nitem) for i in 1:nfact]...)
V = x[(nobs+1):(2*nobs)]
A = sparse(Ind, Jnd, V, nnod, nnod)

invia = I + A
imp_cov_sym = F*invia*S*permutedims(invia)*permutedims(F)
imp_cov_sym = Array(imp_cov_sym)
imp_cov_sym = Symbolics.simplify.(imp_cov_sym)
imp_cov_sym = imp_cov_sym[tril(trues(size(F, 1), size(F, 1)))]

∇Σ_sym = Symbolics.sparsejacobian(imp_cov_sym, x)
#∇²Σ_sym = Symbolics.sparsejacobian(vec(permutedims(∇Σ_sym)), par)
∇Σ_sym = Array(∇Σ_sym)
# H_array = [Symbolics.sparsejacobian(∇Σ_sym[i, :], par) for i = 1:66]
@time H_array = [Symbolics.sparsehessian(imp_cov_sym[i], x) for i = 1:size(imp_cov_sym, 1)]

HT = zeros(Num,size(H_array[1])...)
@variables J[1:size(H_array, 1)]

for i in 1:size(H_array, 1)
    HT += J[i]*H_array[i]
end

using StatsBase
nobs = 80
npar = 100
nnd = Int64(nobs*(nobs+1)/2)

testsym = StatsBase.sample(imp_cov_sym, nnd)
testpar = StatsBase.sample(x, npar)

∇testsym = Symbolics.sparsejacobian(testsym, testpar)
∇testsym = Array(∇testsym)
#∇²Σ_sym = ModelingToolkit.jacobian(vec(permutedims(∇Σ_sym)), x)
@time H_array = [Symbolics.sparsejacobian(∇testsym[i, :], testpar) for i = 1:nnd]
# H_array = [sparse(∇²Σ_sym[((i-1)*31+1):i*31, :]) for i = 1:66]

#= jacobian_fun =
    eval(ModelingToolkit.build_function(
        ∇Σ_sym,
        x
    )[2]) =#

hessian_fun =
    eval(ModelingToolkit.build_function(
        ∇²Σ_sym,
        x
    )[2])

H = zeros(2046, 31)
@benchmark hessian_fun(H, randpar)
#H_array = [sparse(H[((i-1)*31+1):i*31, :]) for i = 1:66]

#model_ls_ana(randpar)

J = (-2*(grad_ls.s-model_ls_ana.imply.imp_cov)'*grad_ls.V)'

@variables J[1:nnd]

HT = zeros(Num, npar, npar)
for i in 1:nnd
    HT += J[i]*H_array[i]
end

@variables Jsym[1:66]

HT = zeros(Num, 31,31)
for i in 1:66
    HT += Jsym[i]*H_array[i]
end

HT = simplify.(HT)

hessian_fun =
    eval(ModelingToolkit.build_function(
        HT,
        x,Jsym
    )[2])

T = zeros(31,31)

@benchmark hessian_fun(T, randpar, J)

T

maximum(abs.(H0 - (initial_P + T)))

f(Hsum, x, J)

maximum(abs.(H0 - (initial_P+HT)))

vectorized_hessian_fun =
    eval(ModelingToolkit.build_function(
        H_array,
        x
    )[2])

vectorized_hessian_fun_alloc =
    eval(ModelingToolkit.build_function(
        H_array,
        x
    )[1])

function mypre(x)
    return initial_P
end

function precondprep_WLS(P, x, model)
    model.imply.gradient_fun(model.imply.∇Σ, x)
    H = 2*model.imply.∇Σ'*model.diff.functions[1].V*model.imply.∇Σ
    return H
end

A_ldiv_B!(pgr, P, gr) = copyto!(pgr, P \ gr)
dot(x, P, y) = dot(x, P*y)

D = duplication_matrix(observed.n_man)
S = inv(observed.obs_cov)
S = kron(S,S)
V = 0.5*(D'*S*D)

LinearAlgebra.ldiv!(P::Matrix{Float64}, b::Vector{Float64}) = P \ b
LinearAlgebra.ldiv!(x::Vector{Float64}, P::Matrix{Float64}, b::Vector{Float64}) = copyto!(x, P \ b)

########################## one factor model ###########################

using Symbolics, SparseArrays, LinearAlgebra, BenchmarkTools

nfact = 1
nitem = 5

## Model definition
nobs = nfact*nitem
nnod = nfact+nobs
@variables x[1:Int64(nobs + nobs)]

x = [x[i] for i = 1:size(x, 1)]
#F
Ind = collect(1:nobs)
Jnd = collect(1:nobs)
V = fill(1,nobs)
F = sparse(Ind, Jnd, V, nobs, nnod)

#S
Ind = collect(1:nnod)
Jnd = collect(1:nnod)
V = [x[1:nobs]; fill(1.0, nfact)]
S = sparse(Ind, Jnd, V, nnod, nnod)

#A
Ind = collect(1:nobs)
Jnd = vcat([fill(nobs+i, nitem) for i in 1:nfact]...)
V = x[(nobs+1):(2*nobs)]
A = sparse(Ind, Jnd, V, nnod, nnod)

S[6,6] = x[6]
A[1,6] = 1

par = [x[i] for i = 1:size(x, 1)]
invia = I + A + A^2 + A^3 + A^4
imp_cov_sym = F*invia*S*permutedims(invia)*permutedims(F)
imp_cov_sym = Array(imp_cov_sym)
imp_cov_sym = Symbolics.simplify.(imp_cov_sym)
imp_cov_sym = imp_cov_sym[tril(trues(size(F, 1), size(F, 1)))]
#imp_cov_sym = vec(imp_cov_sym)
∇Σ_sym = Symbolics.sparsejacobian(imp_cov_sym, par)
∇²Σ_sym = Symbolics.sparsejacobian(vec(permutedims(∇Σ_sym)), par)
∇Σ_sym = Array(∇Σ_sym)
# H_array = [Symbolics.sparsejacobian(∇Σ_sym[i, :], par) for i = 1:66]
H_array = [Symbolics.sparsehessian(imp_cov_sym[i], par) for i = 1:size(imp_cov_sym, 1)]

@Symbolics.variables x[1:3, 1:3]

a = reshape([x[i, j] for j in 1:3 for i in 1:3], (3,3))

a = [x[1] x[2]
    x[3] x[4]]

s = inv(a)

@Symbolics.variables λ[1:2], ω[1:4]

S =[ω[1]  0     0     0
    0     ω[2]  0     0
    0     0     ω[3]  0
    0     0     0     ω[4]]

F =[1.0 0 0 0
    0 1 0 0
    0 0 1 0]

A =[0  0  0  1
    0  0  0  λ[1]
    0  0  0  λ[2]
    0  0  0  0]

S = sparse(S)

#F
F = sparse(F)

#A
A = sparse(A)

invia = I + A + A^2
imp_cov_sym = F*invia*S*permutedims(invia)*permutedims(F)
imp_cov_sym = Array(imp_cov_sym)
imp_cov_sym = Symbolics.simplify.(imp_cov_sym)

Σ_inv = inv(imp_cov_sym)

@Symbolics.variables s[1:3, 1:3]

Σ_obs = reshape([s[i, j] for j in 1:3 for i in 1:3], (3,3))

F = log(det(imp_cov_sym)) + tr(Σ_inv*Σ_obs)

F = simplify(F)


text:
const src_dir = dirname(normpath(@__DIR__, ".."))
include(string(src_dir, "/walk_logic/run_sim!.jl"))
include(string(src_dir, "/walk_logic/one_step!.jl"))
include(string(src_dir, "/walk_logic/approach_wall!.jl"))

"""
    boundary_check(initxy, dx, dy) -> Tuple{Array{Float64,1}, String}

Confirm the next coordinates after initxy with unit displacement vectors dx, dy.

Resolve collisions with escape boundaries, sensors, and walls.
"""
function boundary_check(initxy::Array{Float64,1}, dx::Float64, dy::Float64)
    proposedxy, endinglayer = evaluate_proposed(initxy, dx, dy)
    initx, inity = initxy
    propx, propy = proposedxy
    endingstepsize = STEP_SIZE_DICT[ endinglayer ]

    if inwalls(propx, propy)
        return wallcases(initxy, dx, dy, endingstepsize)
    elseif !inescapebounds(propx, propy)
        return undef, "escape"
    elseif insensor(proposedxy)
        return undef, sensorcases(propx) # not sensorcases(initx) in walls walk_logic
    else
        return proposedxy, "no collision"
    end
end

"sort between the four identical sensors"
function sensorcases(propx::Float64)::String
    if inspawnsensorx(propx)
        return "spawn sensor"
    elseif insecondsensorx(propx)
        return "second sensor"
    elseif inthirdsensorx(propx)
        return "third sensor"
    elseif infourthsensorx(propx)
        return "fourth sensor"
    end
end

"No side or top wall collision"
function wallcases(initxy::Array{Float64,1},
                   dx::Float64, dy::Float64,
                   ending_step_size::Float64)
    return approach_wall!(initxy, dx, dy, ending_step_size), "no collision"
end

"""
Same as walk_logic
"""
function north!(proposedxy::Array{Float64,1},
                initialxy::Array{Float64,1},
                dx::Float64, dy::Float64,
                start::String, endup::String)
    initx, inity = initialxy
    propx, propy = proposedxy
    scale = STEP_SIZE_DICT[start]
    scale2 = STEP_SIZE_DICT[endup]
    slope = (propx - initx) / (propy - inity)           # Reciprocal of traditional slope (delta x over delta y)

    if start == "ppd" || endup == "ppd"
        yintersection = PPD_MAX_Y
    else
        yintersection = ENZYME_MAX_Y
    end

    xintersection = slope * (yintersection - inity) + initx         # x-component of the head of first segment, ending at the border crossing
                # Reciprocal slope lets us multiply the change in y to obtain the change in x
    ydistleft = scale * dy + (inity - yintersection)                # y-component legnth of the original second segment
    xdistleft = scale * dx + (initx - xintersection)                # x-component length of the original second segment
    # First divide the x- and y- component lengths of the original segment steplength, scale -> recover the unit vector length of second segment
    # Then multiply the x- and y-components of the unit vector length of  second semester by scale2 -> obtain correct length of the second segment in new layer
    # Lastly append the correct lengths of the new second segment to the head of the first segment vector
    proposedxy[1] = xintersection + xdistleft * scale2/scale
    proposedxy[2] = yintersection + ydistleft * scale2/scale
    return proposedxy
end

"""
    eastwest!(proposedxy, initialxy, dx, dy, start, endup, EastORWest) -> Array{Float64,1}

Corrects displacement vectors with head and tail in two different layers and
crossing the vertical boundary between water-enzyme, ppd-enzyme.
"""
function eastwest!(proposedxy::Array{Float64,1},
                    initialxy::Array{Float64,1},
                    dx::Float64, dy::Float64,
                    start::String, endup::String,
                    EastOrWest::String)
    initx, inity = initialxy
    propx, propy = proposedxy
    scale = STEP_SIZE_DICT[start]
    scale2 = STEP_SIZE_DICT[endup]
    slope = (propy - inity) / (propx - initx)           # Traditional slope (delta y over delta x)

    enzwater = start == "water" || endup == "water"
    if EastOrWest == "E"
        if enzwater
            xintersection = ENZYME_RIGHT_X
        else
            xintersection = PPD_RIGHT_X
        end
    else
        if enzwater
            xintersection = ENZYME_LEFT_X
        else
            xintersection = PPD_LEFT_X
        end
    end

    yintersection = slope * (xintersection - initx) + inity         # y-component of the head of first segment, ending at the border crossing
                # Multiply the slope and the change in y to obtain the change in x
    xdistleft = scale * dx + (initx - xintersection)                # x-component legnth of the original second segment
    ydistleft = scale * dy + (inity - yintersection)                # y-component legnth of the original second segment
    # Divide the x- and y- component lengths of the original segment steplength, scale -> recover the unit vector length of second segment
    # Multiply the x- and y-components of the unit vector length of  second semester by scale2 -> obtain correct length of the second segment in new layer
    # Append the correct lengths of the new second segment to the head of the first segment vector
    proposedxy[1] = xintersection + xdistleft * scale2/scale
    proposedxy[2] = yintersection + ydistleft * scale2/scale
    return proposedxy
end

"""
    evaluate_proposed(initxy, dx, dy) -> Tuple{Array{Float64,1},String}

Identify location of xy coordinate and compute correct steplengths/displacement vectors.
"""
function evaluate_proposed(initxy::Array{Float64,1}, dx::Float64, dy::Float64)
    adjustedxy, whichsensor = adjustxy(initxy)
    if inwater(adjustedxy)
        newxy, endinglayer = watercalc(adjustedxy, dx, dy)
    elseif inenz(adjustedxy)
        newxy, endinglayer = enzcalc(adjustedxy, dx, dy)
    elseif inppd(adjustedxy)
        newxy, endinglayer = ppdcalc(adjustedxy, dx, dy)
    else
        println("evaluate_proposed Error", initxy)
        return [-1.0, -1.0], "evaluate_proposed Error"
    end
    return revertxy(newxy, whichsensor), endinglayer
end

"Helper when xy is initially in water"
function watercalc(initxy::Array{Float64,1}, dx::Float64, dy::Float64)::Tuple{Array{Float64,1},String}
    initx, inity = initxy
    proposedxy = [initx + WATER_STEP_SIZE * dx, inity + WATER_STEP_SIZE * dy]

    if inwater(proposedxy)
        return proposedxy, "water"
    elseif inenz(proposedxy)
        directionOfEntry = whereoutsideenz(initxy) # tail in water
        if directionOfEntry == "whereoutsideenz Error"
            println(directionOfEntry, " in watercalc ", initxy, proposedxy)
        elseif directionOfEntry == "N"
            confirmxy = north!(proposedxy, initxy, dx, dy, "water", "enz")
        elseif directionOfEntry == "NE"
            confirmxy = WATER_TO_ENZ_NORTHEAST
        elseif directionOfEntry == "NW"
            confirmxy = WATER_TO_ENZ_NORTHWEST
        else # E or W
            confirmxy = eastwest!(proposedxy, initxy, dx, dy, "water", "enz", directionOfEntry)
        end
        return confirmxy, "enz"
    elseif inppd(proposedxy)
        print("Water to ppd is impossible in Quinto", initxy, proposedxy)
    else # collision
        return proposedxy, "water"
    end
end

"Helper when xy is initially in the enzyme layer"
function enzcalc(initxy::Array{Float64,1}, dx::Float64, dy::Float64)::Tuple{Array{Float64,1},String}
    initx, inity = initxy
    proposedxy = [initx + ENZ_STEP_SIZE * dx, inity + ENZ_STEP_SIZE * dy]

    if inenz(proposedxy)
        return proposedxy, "enz"
    elseif inwater(proposedxy)
        directionOfEntry = whereoutsideenz(proposedxy) # head in water
        if directionOfEntry == "whereoutsideenz Error"
            println(directionOfEntry, " in enzcalc ", initxy, proposedxy)
        elseif directionOfEntry == "N"
            confirmxy = north!(proposedxy, initxy, dx, dy, "enz", "water")
        elseif directionOfEntry == "NE"
            confirmxy = ENZ_TO_WATER_NORTHEAST
        elseif directionOfEntry == "NW"
            confirmxy = ENZ_TO_WATER_NORTHWEST
        else # E or W
            confirmxy = eastwest!(proposedxy, initxy, dx, dy, "enz", "water", directionOfEntry)
        end
        return confirmxy, "water"
    elseif inppd(proposedxy)
        directionOfEntry = whereoutsideppd(initxy) # head in water
        if directionOfEntry == "whereoutsideppd Error"
            println(directionOfEntry, " in enzcalc ", initxy, proposedxy)
        elseif directionOfEntry == "N"
            confirmxy = north!(proposedxy, initxy, dx, dy, "enz", "ppd")
        elseif directionOfEntry == "NE"
            confirmxy = ENZ_TO_PPD_NORTHEAST
        elseif directionOfEntry == "NW"
            confirmxy = ENZ_TO_PPD_NORTHWEST
        else # E or W
            confirmxy = eastwest!(proposedxy, initxy, dx, dy, "enz", "ppd", directionOfEntry)
        end
        return confirmxy, "ppd"
    else # collision
        return proposedxy, "enz"
    end
end

"Helper when xy is initially in the ppd layer"
function ppdcalc(initxy::Array{Float64,1}, dx::Float64, dy::Float64)::Tuple{Array{Float64,1},String}
    initx, inity = initxy
    PPD_STEP_SIZE = STEP_SIZE_DICT["ppd"]
    proposedxy = [initx + PPD_STEP_SIZE * dx, inity + PPD_STEP_SIZE * dy] # no flow

    if inppd(proposedxy)
        return proposedxy, "ppd"
    elseif inenz(proposedxy)
        directionOfEntry = whereoutsideppd(proposedxy) # head in water
        if directionOfEntry == "whereoutsideppd Error"
            println(directionOfEntry, " in enzcalc ", initxy, proposedxy)
        elseif directionOfEntry == "N"
            confirmxy = north!(proposedxy, initxy, dx, dy, "ppd", "enz")
        elseif directionOfEntry == "NE"
            confirmxy = PPD_TO_ENZ_NORTHEAST
        elseif directionOfEntry == "NW"
            confirmxy = PPD_TO_ENZ_NORTHWEST
        else # E or W
            confirmxy = eastwest!(proposedxy, initxy, dx, dy, "ppd", "enz", directionOfEntry)
        end
        return confirmxy, "ppd"
    elseif inwater(proposedxy)
        print("PPD to water is impossible in Quinto", initxy, proposedxy)
    else # collision
        return proposedxy, "ppd"
    end
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
printstyled("  AH (Ad Hoc)\n", color=:light_green)
ah1_file  = path*"/SampleFiles/AH/ah1.f"
ah1_fstr  = path*"/SampleFiles/AH/ah1.*"
ahc_file  = path*"/SampleFiles/AH/lhz.ah"
ah_resp   = path*"/SampleFiles/AH/BRV.TSG.DS.lE21.resp"
ah2_file  = path*"/SampleFiles/AH/ah2.f"
ah2_fstr  = path*"/SampleFiles/AH/ah2.*"

printstyled("    v1\n", color=:light_green)
redirect_stdout(out) do
  S = verified_read_data("ah1", ah1_file, v=3)
  S = verified_read_data("ah1", ah1_file, full=true)
  S = verified_read_data("ah1", ah1_fstr, full=true, v=3, strict=false)
  @test S.n == 3

  S = read_data("ah1", ah1_fstr, full=true, v=3, strict=true)
  @test S.n == 4
  @test S.fs[1] == 4.0
  @test isapprox(S.gain[1], 64200.121094)
  @test isapprox(S.loc[1].lat, 35.599899)
  @test isapprox(S.loc[1].lon, -85.568802)
  @test isapprox(S.loc[1].el, 481.0)
  @test any([occursin("gdsn_tape",s) for s in S.notes[1]])
  @test any([occursin("demeaned",s) for s in S.notes[1]])
  @test length(S.resp[1].p) == 24
  @test length(S.resp[1].z) == 7
  @test u2d(S.t[1][1,2]*μs) == DateTime("1984-04-20T06:42:00.12")
  @test length(S.x[1]) == 720
  @test eltype(S.x[1]) == Float32
  @test isapprox(S.x[1][1:4], [-731.41247559, -724.41247559, -622.41247559, -470.4125061])

  C = verified_read_data("ah1", ahc_file, v=1, full=true)[1]

  # Station
  @test isapprox(C.loc.lat, 36.5416984)
  @test isapprox(C.loc.lon, 138.2088928)
  @test isapprox(C.loc.el, 422.0)
  @test isapprox(C.gain, 1.178e8)
  @test length(C.resp.p) == 10
  @test isapprox(C.resp.p[1], -0.0123f0 + 0.0123f0im)
  @test isapprox(C.resp.p[2], -0.0123f0 - 0.0123f0im)
  @test length(C.resp.z) == 3

  # Data
  @test length(C.x) == 1079
  @test eltype(C.x) == Float32
  @test C.fs == 1.0
  @test u2d(C.t[1,2]*μs) == DateTime("1990-05-12T04:49:54.49")

  # Event
  @test isapprox(C.misc["ev_lat"], 49.037)
  @test isapprox(C.misc["ev_lon"], 141.847)
  @test isapprox(C.misc["ev_dep"], 606.0)
  @test string(u2d(C.misc["ot"]*μs)) == "1990-05-12T04:50:08.7"
  @test startswith(C.misc["data_comment"], "Streckeisen STS-1V/VBB Seismometer")
  @test startswith(C.misc["event_comment"], "null")

  C = verified_read_data("ah1", ah_resp, full=true, vl=true)[1]
  @test isapprox(C.loc.lat, 53.058060)
  @test isapprox(C.loc.lon, 70.282799)
  @test isapprox(C.loc.el, 300.0)
  @test isapprox(C.gain, 0.05)
  @test isapprox(C.resp.a0, 40.009960)
  @test length(C.resp.p) == 7
  @test isapprox(C.resp.p[1], -0.1342653f0 + 0.1168836f0im)
  @test length(C.resp.z) == 4
  @test startswith(C.misc["data_comment"], "DS response in counts/nm")
  @test startswith(C.misc["event_comment"], "Calibration_for_hg_TSG")
  @test any([occursin("brv2ah: ahtedit",s) for s in C.notes])
  @test any([occursin("demeaned",s) for s in C.notes])
  @test any([occursin("modhead",s) for s in C.notes])
  @test any([occursin("ahtedit",s) for s in C.notes])
end
printstyled("      append existing channel\n", color=:light_green)
test_chan_ext(ah1_file, "ah1", "nu.RSN..IPZ", 4.0, 1, 451291110190001)
test_chan_ext(ah1_file, "ah1", "nu.RSC..IPZ", 4.0, 1, 451291320120000)

printstyled("    v2\n", color=:light_green)
redirect_stdout(out) do
  S = verified_read_data("ah2", ah2_file, v=3)
  S = verified_read_data("ah2", ah2_file, v=3, full=true)
  S = verified_read_data("ah2", ah2_fstr, v=3, full=true, vl=true, strict=false)
  @test S.n == 1

  S = read_data("ah2", ah2_fstr, v=3, full=true, strict=true, vl=true)
  @test S.n == 4
  @test S.fs[1] == 4.0
  @test isapprox(S.gain[1], 64200.121094)
  @test isapprox(S.loc[1].lat, 35.599899)
  @test isapprox(S.loc[1].lon, -85.568802)
  @test isapprox(S.loc[1].el, 481.0)
  @test any([occursin("gdsn_tape",s) for s in S.notes[1]])
  @test any([occursin("demeaned",s) for s in S.notes[1]])
  @test length(S.resp[1].p) == 24
  @test length(S.resp[1].z) == 7
  @test u2d(S.t[1][1,2]*μs) == DateTime("1984-04-20T06:42:00.12")
  @test length(S.x[1]) == 720
  @test eltype(S.x[1]) == Float32
  @test isapprox(S.x[1][1:4], [-731.41247559, -724.41247559, -622.41247559, -470.4125061])
end
printstyled("      append existing channel\n", color=:light_green)
test_chan_ext(ah2_file, "ah2", "nu.RS..IP", 4.0, 1, 451291320120000)

printstyled("      custom user attributes\n", color=:light_green)
tmp_ah = "tmp.ah"
io = open(ah2_file, "r")
buf = read(io, 10636)
skip(io, 4)
append!(buf, reinterpret(UInt8, [bswap(Int32(1))]))
append!(buf, reinterpret(UInt8, [bswap(Int32(6))]))
append!(buf, codeunits("PEBKAC\0\0"))
append!(buf, reinterpret(UInt8, [bswap(Int32(41))]))
append!(buf, codeunits("problem exists between keyboard and chair\0\0\0"))
append!(buf, read(io))
close(io)

ah_out = open(tmp_ah, "w")
write(ah_out, buf)
close(ah_out)

S = read_data("ah2", tmp_ah, strict=true)
@test haskey(S.misc[4], "PEBKAC")
@test S.misc[4]["PEBKAC"] == "problem exists between keyboard and chair"
safe_rm(tmp_ah)


text:
using JLD2
using Random
using LinearAlgebra
using Statistics
using CSV
using DataFrames
using FreqTables
using Distributions

function q1()
	cd("C:/Users/jakek/Documents/econ_phd/structural_modeling_course/forked_class/ProblemSets/PS1-julia-intro")
	### 1.
	# a.
	Random.seed!(1234)
	A = 15 .*rand(10,7) .-5
	B = sqrt(15) .*randn(10,7) .-2
	C = [A[1:5, 1:5] B[1:5,6:7]]
	D = [max(elem,0) for elem in A]

	# b.
	length(A)
	# c.
	length(unique(D))
	# d
	E = reshape(B, (1, length(B)))'
	E_alt = vec(B)
	# e
	F = cat(A,B; dims=3)
	# f
	F = permutedims(F, [3,1,2])
	# g
	G = kron(B,C)
	# error opccurs when we try kron(F,G) due to differences in dimensions.
	# h
	@save "matrixpractice.jld" A B C D E F G
	# i
	@save "firstmatrix.jld" A B C D
	# j
	C_csv = convert(DataFrame,C)
	CSV.write("Cmatrix.csv",C_csv)
	# k
	D_dat = convert(DataFrame,D)
	CSV.write("Dmatrix.dat",D_dat, delim="\t")
	return(A,B,C,D)
end

A,B,C,D = q1()

## 2.
# a.
function q2(A,B,C,D)
	AB = [A[i,j]*A[i,j] for i=1:size(A,1),j=1:size(A,2)]
	AB2 = A .*B
	# b.
	Cprime = Float64[]
	for i in 1:length(C)
		if abs(C[i])<5
		append!(Cprime, C[i])
		end
	end
	Cprime2 = C[abs.(C) .<5]
	Cprime == Cprime2

	# c
	X = Array{Float64}(undef, 15169, 6,5)
	for t in 1:size(X,3)
	X[1:15169,1,t] .= 1
	X[1:15169,2,t] .= rand(Bernoulli(0.75*(6-t)/5),15169)
	X[1:15169,3,t] .= vec((5*(t-1)) .*randn(15169,1) .+(15+t-1))
	X[1:15169,4,t] .= vec(1/MathConstants.e .*randn(15169,1) .+(MathConstants.pi*(6-t)/3))
	X[1:15169,5,t] .= vec(rand(Binomial(20,0.6),15169))
	X[1:15169,6,t] .= vec(rand(Binomial(20,0.5),15169))
	end

	# d.
	beta = [(0.25*t +0.75 , log(t),-t^(1/2),exp(t)-exp(t+1),t,t/3) for t in 1:5]
	beta = reduce(hcat, getindex.(beta,i) for i in eachindex(beta[1]))'

	# e
	Y = [ X[n,:,t]'*beta[:,t] + randn(1,1)[1]*0.36 for t in 1:5, n in 1:15169]'
end

q2(A,B,C,D)

function q3()
	## 3.
	# a.
	data = CSV.read("nlsw88.csv")
	@save "nlsw88.jld" data
	# b.
	mean(data.never_married)
	mean(data.collgrad)
	# c
	prop(freqtable(data.race))
	# d
	summarystats=describe(data)
	# e
	freqtable(data.industry, data.occupation)
	# f split apply combine part.
	subset = select(data, :wage, :occupation, :industry)
	combine(groupby(subset, All(:industry, :occupation)), :wage => mean)
end
q3()

function q4()
	# 4.
	# a
	@load "firstmatrix.jld"
	# b., c,d,e
	function matrixops(x,y)
		if  size(x) != size(y)
			error("inputs must have the same size")
		end
		# creates 3 outputs: element by element product, product, and sum
		return(x .*y, x'*y , x .+y)
	end
	E1, E2, E3 = matrixops(A,B)
	E1
	E2
	E3
	# f
	try
		matrixops(C,D)
	catch err
		print(err)
	end
	# g
	@load "nlsw88.jld"
	ttl_exp = convert(Array,data.ttl_exp)

	wage = convert(Array,data.wage)
	E1, E2, E3 = matrixops(ttl_exp, wage)
	E1
	E2
	E3
end

q4()


text:
# Computes derivatives over the timestep.
using PyPlot
using Limbdark
import Limbdark: Transit_Struct
# include("../../../src/integrate_lightcurve.jl")

# Test it out:

include("../../../test/loglinspace.jl")

t1 = -1.5; t2 = 1.5; nt = 10000; dt = 0.3
t = zeros(nt)
t .= linearspace(t1,t2,nt)
# The following compares different tolerances and maxdepths:
r = 0.1; b0 = 0.5; u_n = [0.3,0.3]; nu = length(u_n)
favg0 = zeros(nt,5+nu)
favg1 = zeros(nt,5+nu)
neval1= zeros(Int64,nt)
depthmax1= zeros(Int64,nt)
trans = transit_init(r,b0,u_n,true)
param = [0.0,1.0,b0]   # [t_0,v,b_0]

# Now, carry out speed test:

function compute_lightcurve!(trans::Transit_Struct{T},param::Array{T,1},t::Array{T,1},favg0::Array{T,2},nt::Int64) where {T <: Real}
@inbounds for i=1:nt
  trans.b = sqrt(param[3]^2+(param[2]*(t[i]-param[1]))^2)
  favg0[i,1] =  transit_poly!(trans)-1
  favg0[i,2] =  trans.dfdrb[1]
  favg0[i,3] =  trans.dfdrb[2]/trans.b*param[2]^2*(param[1]-t[i])
  favg0[i,4] =  trans.dfdrb[2]*param[2]/trans.b*(t[i]-param[1])^2
  favg0[i,5] =  trans.dfdrb[2]*param[3]/trans.b
  favg0[i,6:5+trans.n] = trans.dfdu
end
return
end

compute_lightcurve!(trans,param,t,favg0,nt)
@time compute_lightcurve!(trans,param,t,favg0,nt)

integrate_lightcurve!(trans,param,t,dt,favg1,nt,1e-10,30,neval1,depthmax1)

# Now plot the results:
fig,axes = subplots(2,3,figsize=(12,8))
if VERSION >= v"0.7"
  ax = axes[1]
  ax.plot(t,favg0[:,1].+1,linestyle="--",linewidth=2,label=L"${F}$")
  ax.plot(t,favg1[:,1].+1,linewidth=2,label=L"$\overline{{F}}$")
  ax.plot([-dt/2,-dt/2,-dt/2,dt/2,dt/2,dt/2],[0.9925,0.9935,0.993,0.993,0.9925,0.9935],color="b")
  ax.annotate(L"$\Delta t$",xy=[-0.08,0.991])
  ax.legend(loc="upper center")
  ax.axis([-1.5,1.5,0.988,1.002])
  ax.set_xlabel("Time")
  ax.set_ylabel("Flux")
  ax = axes[2]
  ax.plot(t,favg0[:,2]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial r}$")
  ax.plot(t,favg1[:,2]*100.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial r}$")
  ax.legend(loc="upper center")
  ax.axis([-1.5,1.5,-30,2])
  ax.set_xlabel("Time")
  ax.set_ylabel("Flux derivative [pct]")
  ax = axes[3]
  ax.plot(t,favg0[:,3]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial t_0}$")
  ax.plot(t,favg1[:,3]*100.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial t_0}$")
  ax.set_xlabel("Time")
  ax.set_ylabel("Flux derivative [pct]")
  ax.legend(loc="lower left")
  ax = axes[4]
  ax.plot(t,favg0[:,4]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial v}$")
  ax.plot(t,favg1[:,4]*100.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial v}$")
  ax.set_xlabel("Time")
  ax.set_ylabel("Flux derivative [pct]")
  ax.axis([-1.5,1.5,-1,6])
  ax.legend(loc="upper center")
  ax = axes[5]
  ax.plot(t,favg0[:,5]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial b_0}$")
  ax.plot(t,favg1[:,5]*100.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial b_0}$")
  ax.set_xlabel("Time")
  ax.set_ylabel("Flux derivative [pct]")
  ax.legend(loc="upper center")
  ax = axes[6]
  ax.plot(t,favg0[:,6]*1000.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial u_1}$")
  ax.plot(t,favg1[:,6]*1000.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial u_1}$")
  ax.plot(t,favg0[:,7]*1000.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial u_2}$")
  ax.plot(t,favg1[:,7]*1000.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial u_2}$")
  ax.set_xlabel("Time")
  ax.set_ylabel("Flux derivative [ppt]")
  ax.legend(loc="upper center")
  #b = sqrt.(param[3]^2+(param[2]*(t-param[1])).^2)
  #plot(dt^2*(d2fdb2.*(param[2]*(t-param[1])./b).^2+dfdb.*param[2]^2./b)/24.0)
else
  # Now plot the results:
  fig,axes = subplots(2,3,figsize=(12,8))
  ax = axes[1]
  ax[:plot](t,favg0[:,1].+1,linestyle="--",linewidth=2,label=L"${F}$")
  ax[:plot](t,favg1[:,1].+1,linewidth=2,label=L"$\overline{{F}}$")
  ax[:plot]([-dt/2,-dt/2,-dt/2,dt/2,dt/2,dt/2],[0.9925,0.9935,0.993,0.993,0.9925,0.9935],color="b")
  ax[:annotate](L"$\Delta t$",xy=[-0.08,0.991])
  ax[:legend](loc="upper center")
  ax[:axis]([-1.5,1.5,0.988,1.002])
  ax[:set_xlabel]("Time")
  ax[:set_ylabel]("Flux")
  ax = axes[2]
  ax[:plot](t,favg0[:,2]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial r}$")
  ax[:plot](t,favg1[:,2]*100.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial r}$")
  ax[:legend](loc="upper center")
  ax[:axis]([-1.5,1.5,-30,2])
  ax[:set_xlabel]("Time")
  ax[:set_ylabel]("Flux derivative [pct]")
  ax = axes[3]
  ax[:plot](t,favg0[:,3]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial t_0}$")
  ax[:plot](t,favg1[:,3]*100.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial t_0}$")
  ax[:set_xlabel]("Time")
  ax[:set_ylabel]("Flux derivative [pct]")
  ax[:legend](loc="lower left")
  ax = axes[4]
  ax[:plot](t,favg0[:,4]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial v}$")
  ax[:plot](t,favg1[:,4]*100.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial v}$")
  ax[:set_xlabel]("Time")
  ax[:set_ylabel]("Flux derivative [pct]")
  ax[:axis]([-1.5,1.5,-1,6])
  ax[:legend](loc="upper center")
  ax = axes[5]
  ax[:plot](t,favg0[:,5]*100.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial b_0}$")
  ax[:plot](t,favg1[:,5]*100.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial b_0}$")
  ax[:set_xlabel]("Time")
  ax[:set_ylabel]("Flux derivative [pct]")
  ax[:legend](loc="upper center")
  ax = axes[6]
  ax[:plot](t,favg0[:,6]*1000.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial u_1}$")
  ax[:plot](t,favg1[:,6]*1000.,linewidth=2,label=L"$\frac{\partial\overline{{F}}}{\partial u_1}$")
  ax[:plot](t,favg0[:,7]*1000.,linestyle="--",linewidth=2,label=L"$\frac{\partial{F}}{\partial u_2}$")
  ax[:plot](t,favg1[:,7]*1000.,linewidth=2,label=L"$\frac{\partial\overline{F}}{\partial u_2}$")
  ax[:set_xlabel]("Time")
  ax[:set_ylabel]("Flux derivative [ppt]")
  ax[:legend](loc="upper center")
end
savefig("integrate_transit_gradient.pdf", bbox_inches="tight")

max_issues_repo_name: rodluger/limbdark


text:
module UtilityFunctions

#Julia Libraries/Modules
using GSL

using TypesParticles, TypesBasis

export gaussprod1D,getgaussoverlap
export flattenbasisfunc
export boysfunction

@doc raw"""
TODO write description

"""
function gaussprod1D(npx::Float64,np_alpha::Float64,
                         mpx::Float64,mp_alpha::Float64)
    p = np_alpha + mp_alpha;
    q = (np_alpha * mp_alpha) / p;
    #ratio = (np_alpha*npx + mp_alpha*mpx)/ p;
    diff = npx-mpx;
    kernel = exp(-q*diff^2);
    return kernel
end #gaussprod1D

@doc raw"""

TODO write description

"""
function getgaussoverlap(np::Int,mp::Int,
                             nbasis::GaussOrbitals,mbasis::GaussOrbitals)
    np_alpha = nbasis.alphas[np];
    mp_alpha = mbasis.alphas[mp];
    spatialx =  gaussprod1D(nbasis.ax,np_alpha,mbasis.ax,mp_alpha);
    spatialy =  gaussprod1D(nbasis.ay,np_alpha,mbasis.ay,mp_alpha);
    spatialz =  gaussprod1D(nbasis.az,np_alpha,mbasis.az,mp_alpha);
    sprod = spatialx*spatialy*spatialz;
    term1 = sqrt(pi/(np_alpha+mp_alpha))^3; #relabel based on equations
    term2 = nbasis.norms[np]*mbasis.norms[mp]; #rekabek based on equations
    primoverlap = sprod*term1*term2;
    return primoverlap
end #getgaussoverlap

@doc """
flattenbasisfunction(natoms::Int, number_basis_functions::Int,
                     basis_functions::Array{GaussOrbitals})

description: This utility function takes a 2-dimensional array with elements having data-type
{GaussOrbitals} and flattens it to a 1-dimensional of data-type {Gaussorbitals} that has as
length equal to number of electrons times number of basis functions.

input:

output:

module requirements: TypesParticles.jl, TypesBasis.jl

"""
function flattenbasisfunc(natoms::Int,numbasisfunc::Int,
                              basisfunc::Array{GaussOrbitals})
    totalbasissize = natoms*numbasisfunc;
    flatbasisfunc = Array{GaussOrbitals,1}(undef,totalbasissize);
    ai = 1;
    for a=1:natoms
        tmpbf = basisfunc[a,:]
        for i=1:numbasisfunc
            flatbasisfunc[ai] = tmpbf[i];
            ai += 1;
        end
    end
    return flatbasisfunc
end #flattenbasisfunc

@doc raw"""
function boysfunction(x::Float64;n::Int=0)

description: This is an implementation for integration of coulomb basis intergral using
boys functions for Gaussian,  that are atomic-centered, basis functions. Requires use of gamma and incomplete gamma functions obtained through GSL library..

input:
output:

module requirements: This function needs access to the gamma and incomplete gamma functions. For this function call I use the GSL (GNU Scientific Library) interface which has sf_gamma and complimentary gamma function sf_gamma_inc_P

"""
function boysfunction(x::Float64;n::Int=0)
    if x == 0.00e0
        out = 1.00e0/(2.00e0*n+1.00e0);
    else
        nhalf = n + 0.5e0;
        gamma = GSL.sf_gamma(nhalf);
        igamma = GSL.sf_gamma_inc_P(nhalf,x);
        out = gamma*igamma/(2*x^(n+0.5e0));
    end
    return out

end #boysfunction

end #module


text:
Pixels = BitArray;
mutable struct Checkerboard
    w::Int64
    h::Int64
    pixelsize::Int64
    nowisalpha::Bool
    pixelsalpha::Pixels
    pixelsbeta::Pixels
end
Checkerboard(width, height) = begin
    offsetx = 5
    offsety = 15
    pixelsize = 10
    boardx = (width-offsetx) ÷ pixelsize
    boardy = (height-offsety) ÷ pixelsize
    pxl = rand(Bool, (boardx, boardy))
    Checkerboard(
        boardx,
        boardy,
        pixelsize,
        true,
        pxl,
        copy(pxl)
    )
end

pixels(board::Checkerboard) = board.nowisalpha ? board.pixelsalpha : board.pixelsbeta
oldpixels(board::Checkerboard) = board.nowisalpha ? board.pixelsbeta : board.pixelsalpha

neighborindices = filter(x -> Tuple(x) != (0,0), CartesianIndex(-1,-1):CartesianIndex(1,1))

function valid_neighbor(x::CartesianIndex,w,h)
    return x[1] > 1 &&
        x[1] < w-1 &&
        x[2] > 1 &&
        x[2] < h-1
end

function num_neighbors(pixels::Pixels, i::CartesianIndex, w, h)
    return count(ni -> valid_neighbor(ni+i,w,h) && pixels[ni+i], neighborindices)
end

function evolve_tick(board::Checkerboard)
    board.nowisalpha = !board.nowisalpha
    newpixels = pixels(board)
    oldp = oldpixels(board)
    for i in CartesianIndices(oldp)
        neigh = num_neighbors(oldp, i, board.w, board.h)
        if neigh == 3
            newpixels[i] = true
        elseif neigh < 2 || neigh > 3
            newpixels[i] = false
        else
            newpixels[i] = oldp[i]
        end
    end
    return nothing
end

function draw_item(board::Checkerboard, canvas::Gtk.GtkCanvas, scene::Scene)
    ctx = Gtk.getgc(canvas)
    pixs = pixels(board)
    for i in CartesianIndices(pixs)
        pixs[i] ? set_source_rgba(ctx,0,0,0,1) : set_source_rgba(ctx,1,1,1,0)
        pixelsize = board.pixelsize
        x = i[1]
        y = i[2]
        rectangle(ctx, 5+((x-1)*pixelsize), 15+((y-1)*pixelsize), pixelsize, pixelsize)
        fill(ctx)
    end
    return nothing
end


text:
"""
    Instance

This module is the second step for getting started with Vulkan. It's purpose is
to create a window with a Vulkan instance.

See https://vulkan-tutorial.com/en/Drawing_a_triangle/Setup/Instance
"""
module Instance

export main

using GLFW
using VulkanCore
using VulkanCore.LibVulkan

"""
    main()

The application's entry point.
"""
main() = (
    buildApp()
    |> initWindow
    |> initVulkan
    |> mainLoop
    |> cleanup
)

"""
The application's state. Anything the main app needs to run will be included
here.
"""
Base.@kwdef mutable struct App
    window::Union{Some{GLFW.Window}, Nothing}
    instance::Ref{VkInstance}
end

"""
Construct a new application instance.
"""
function buildApp()
    App(;
        window = nothing,
        instance = Ref(C_NULL)
    )
end

"""
    initWindow(app::App) :: App

Create the application's main window. Configure to support Vulkan instead of
OpenGL.
"""
function initWindow(app::App) :: App
    # NO_API disables GLFW's attempt to load an OpenGL context
    GLFW.WindowHint(GLFW.CLIENT_API, GLFW.NO_API)

    # Resizing takes extra management, so don't support it for now
    GLFW.WindowHint(GLFW.RESIZABLE, false)

    app.window = GLFW.CreateWindow(800, 600, "Vulkan") |> Some

    app
end

"""
    initVulkan(app::App) :: App

Initialize Vulkan resources for this application.
"""
function initVulkan(app::App) :: App

    # Create an AppInfo struct instance on the heap.
    #
    # Q: Why not just on the stack then Ref it on demand?
    # A: Because this struct is not passed directly into a function. Instead,
    #    the pointer to this struct is saved in the InstanceCreateInfo struct.
    #    This means the location must not change during the subsequent calls.
    appInfoRef = begin
        sType = VK_STRUCTURE_TYPE_APPLICATION_INFO
        pNext = C_NULL
        pApplicationName = pointer("Vulkan Demo")
        applicationVersion = VK_MAKE_VERSION(1, 0, 0)
        pEngineName = pointer("No Engine")
        engineVersion = VK_MAKE_VERSION(1, 0, 0)
        apiVersion = VK_API_VERSION_1_2
        VkApplicationInfo(
            sType,
            pNext,
            pApplicationName,
            applicationVersion,
            pEngineName,
            engineVersion,
            apiVersion
        ) |> Ref
    end

    # Note the usage of `GC.@preserve` to ensure that AppInfo doesn't get GC'd
    # or moved between the creation of the InstanceCreateInfo struct and the
    # call to vkCreateInstance.
    result = GC.@preserve appInfoRef begin
        createInfo = begin
            sType = VK_STRUCTURE_TYPE_INSTANCE_CREATE_INFO
            pNext = C_NULL
            flags = UInt32(0)
            pApplicationInfo =
                Base.unsafe_convert(Ptr{VkApplicationInfo}, appInfoRef)
            enabledExtensionCount = UInt32(0)
            ppEnabledExtensionNames =
                GLFW.GetRequiredInstanceExtensions(Ref(enabledExtensionCount))
            enabledLayerCount = UInt32(0)
            ppEnabledLayerNames = C_NULL
            VkInstanceCreateInfo(
                sType,
                pNext,
                flags,
                pApplicationInfo,
                enabledLayerCount,
                ppEnabledLayerNames,
                enabledExtensionCount,
                ppEnabledExtensionNames
            )
        end

        vkCreateInstance(Ref(createInfo), C_NULL, app.instance)
    end

    if result == VK_SUCCESS
        @info "successfully created vulkan instance"
    else
        error("Unable to create a vulkan instance")
    end

    # Get supported extensions

    # First get the number of extensions
    # Note the usage of ref to get mutable data under the pointer
    refCount = Ref(Cuint(0))
    vkEnumerateInstanceExtensionProperties(C_NULL, refCount, C_NULL)

    # Now get the actual extensions
    extensions = Vector{VkExtensionProperties}(undef, refCount[])
    vkEnumerateInstanceExtensionProperties(
        C_NULL, refCount, pointer(extensions))

    to_string(name) = rstrip(name |> collect |> String, '\0')
    names = [e.extensionName |> to_string for e in extensions]
    @info "supported extensions $names"

    app
end

"""
    mainLoop(app::App) :: App

The application's main loop. This function blocks until the window is closed
by clicking the 'x' or similarly exited.
"""
function mainLoop(app::App) :: App
    window = something(app.window)

    while !GLFW.WindowShouldClose(window)
        GLFW.PollEvents()
    end

    app
end

"""
    cleanup(app::App) :: App

Cleanup any and all graphical resources.
"""
function cleanup(app::App) :: App

    # the second argument is for custom allocator callbacks
    vkDestroyInstance(app.instance[], C_NULL)

    app.window |> something |> GLFW.DestroyWindow
    app.window = nothing
    app
end

end # End Instance Module


text:
module InteractModels
# Use the README as the module docs
@doc let
    path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    read(path, String)
end InteractModels

using Interact,
      Reexport

@reexport using ModelParameters

export InteractModel, ui

include("interactive.jl")

end


text:
type Matrix{T} <: ElementalMatrix{T}
    obj::Ptr{Void}
end

for (elty, ext) in ((:ElInt, :i),
                    (:Float32, :s),
                    (:Float64, :d),
                    (:Complex64, :c),
                    (:Complex128, :z))
    @eval begin
        function Matrix(::Type{$elty})
            obj = Ref{Ptr{Void}}(0)
            err = ccall(($(string("ElMatrixCreate_", ext)), libEl), Cuint,
                (Ref{Ptr{Void}},), obj)
            err == 0 || throw(ElError(err))
            return Matrix{$elty}(obj[])
        end

        function getindex(A::Matrix{$elty}, i::Integer, j::Integer)
            x = Ref{$elty}(0)
            err = ccall(($(string("ElMatrixGet_", ext)), libEl), Cuint,
                (Ptr{Void}, ElInt, ElInt, Ref{$elty}),
                A.obj, i - 1, j - 1, x)
            err == 0 || throw(ElError(err))
            return x[]
        end

        function resize!(A::Matrix{$elty}, i::Integer, j::Integer = 1) # to mimic vector behavior
            err = ccall(($(string("ElMatrixResize_", ext)), libEl), Cuint,
                (Ptr{Void}, ElInt, ElInt),
                A.obj, i, j)
            err == 0 || throw(ElError(err))
            return A
        end

        function height(A::Matrix{$elty})
            rs = Ref{ElInt}(0)
            err = ccall(($(string("ElMatrixHeight_", ext)), libEl), Cuint,
                (Ptr{Void}, Ref{ElInt}),
            A.obj, rs)
            err == 0 || throw(ElError(err))
            return rs[]
        end

        function width(A::Matrix{$elty})
            rs = Ref{ElInt}(0)
            err = ccall(($(string("ElMatrixWidth_", ext)), libEl), Cuint,
                (Ptr{Void}, Ref{ElInt}),
            A.obj, rs)
            err == 0 || throw(ElError(err))
            return rs[]
        end

        function lockedBuffer(A::Matrix{$elty})
            rp = Ref{Ptr{$elty}}(0)
            err = ccall(($(string("ElMatrixLockedBuffer_", ext)), libEl), Cuint,
                (Ptr{Void}, Ref{Ptr{$elty}}),
                A.obj, rp)
            err == 0 || throw(ElError(err))
            return rp[]
        end

        function setindex!(A::Matrix{$elty}, x::Number, i::Integer, j::Integer)
            err = ccall(($(string("ElMatrixSet_", ext)), libEl), Cuint,
                (Ptr{Void}, ElInt, ElInt, $elty),
                A.obj, i - 1, j - 1, x)
            err == 0 || throw(ElError(err))
            return A
        end
    end
end

Matrix() = Matrix(Float64)

# Julia convenience

pointer(A::Matrix) = lockedBuffer(A)

function similar{T}(::Matrix, ::Type{T}, sz::Dims)
    A = Matrix(T)
    resize!(A, sz...)
    return A
end

countnz(A::Matrix) = length(A)


text:
function curve_fit(model::Function, xpts, ydata, p0)
	# assumes model(xpts, params...) = ydata + noise
	# minimizes F(p) = sum(ydata - f(xdata)).^2 using leastsq()
	# returns p, f(p), g(p) where
	#   p - best fit parameters
	#   f(p) - vector of residuals
	#   g(p) - estimated Jacobian at p

	# construct the cost function
	f(p) = model(xpts, p) - ydata

	# construct Jacobian function
	g = Calculus.jacobian(f)

	results = levenberg_marquardt(f, g, p0)
	p = results.minimum
	return p, f(p), g(p)
end

estimate_errors(p, residuals, J) = estimate_errors(p, residuals, J, .95)

function estimate_errors(p, residuals, J, alpha)
	# estimate_errors(p, residuals, J, alpha) computes (1-alpha) error estimates for the parameters from leastsq
	#   p - parameters
	#   residuals - vector of residuals
	#   J - Jacobian
	#   alpha - compute alpha percent confidence interval, (e.g. alpha=0.95 for 95% CI)

	# mean square error is: standard square error / degrees of freedom
	n, p = size(J)
	mse = sse(residuals)/(n-p)

	# compute the covariance matrix from the QR decomposition
	Q,R = qr(J)
	Rinv = inv(R)
	covar = Rinv*Rinv*mse

	# then the standard errors are given by the sqrt of the diagonal
	std_error = sqrt(diag(covar))

	# scale by quantile of the student-t distribution
	dist = TDist(n-p)
	std_error *= quantile(dist, alpha)
end


text:
module EMLturbine
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...streams.EMLstreams
	import EMLtypes.length
	include("turbine/HidraulicTurbine.jl")
	include("turbine/HidraulicTurbineGenerator.jl")
end


text:
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_band.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_config.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_dense.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_direct.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_iterative.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_linearsolver.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function SUNLinSolGetType(S::SUNLinearSolver)
    ccall((:SUNLinSolGetType, libsundials_sundials), SUNLinearSolver_Type, (SUNLinearSolver,), S)
end

function SUNLinSolSetATimes(S::SUNLinearSolver, A_data, ATimes::ATimesFn)
    ccall((:SUNLinSolSetATimes, libsundials_sundials), Cint, (SUNLinearSolver, Ptr{Void}, ATimesFn), S, A_data, ATimes)
end

function SUNLinSolSetPreconditioner(S::SUNLinearSolver, P_data, Pset::PSetupFn, Psol::PSolveFn)
    ccall((:SUNLinSolSetPreconditioner, libsundials_sundials), Cint, (SUNLinearSolver, Ptr{Void}, PSetupFn, PSolveFn), S, P_data, Pset, Psol)
end

function __SUNLinSolSetScalingVectors(S::SUNLinearSolver, s1::N_Vector, s2::N_Vector)
    ccall((:SUNLinSolSetScalingVectors, libsundials_sundials), Cint, (SUNLinearSolver, N_Vector, N_Vector), S, s1, s2)
end

function SUNLinSolSetScalingVectors(S, s1, s2)
    __s1 = convert(NVector, s1)
    __s2 = convert(NVector, s2)
    __SUNLinSolSetScalingVectors(S, convert(N_Vector, __s1), convert(N_Vector, __s2))
end

function SUNLinSolInitialize(S::SUNLinearSolver)
    ccall((:SUNLinSolInitialize, libsundials_sundials), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolSetup(S::SUNLinearSolver, A::SUNMatrix)
    ccall((:SUNLinSolSetup, libsundials_sundials), Cint, (SUNLinearSolver, SUNMatrix), S, A)
end

function __SUNLinSolSolve(S::SUNLinearSolver, A::SUNMatrix, x::N_Vector, b::N_Vector, tol::realtype)
    ccall((:SUNLinSolSolve, libsundials_sundials), Cint, (SUNLinearSolver, SUNMatrix, N_Vector, N_Vector, realtype), S, A, x, b, tol)
end

function SUNLinSolSolve(S, A, x, b, tol)
    __x = convert(NVector, x)
    __b = convert(NVector, b)
    __SUNLinSolSolve(S, A, convert(N_Vector, __x), convert(N_Vector, __b), tol)
end

function SUNLinSolNumIters(S::SUNLinearSolver)
    ccall((:SUNLinSolNumIters, libsundials_sundials), Cint, (SUNLinearSolver,), S)
end

function SUNLinSolResNorm(S::SUNLinearSolver)
    ccall((:SUNLinSolResNorm, libsundials_sundials), realtype, (SUNLinearSolver,), S)
end

function SUNLinSolResid(S::SUNLinearSolver)
    ccall((:SUNLinSolResid, libsundials_sundials), N_Vector, (SUNLinearSolver,), S)
end

function SUNLinSolLastFlag(S::SUNLinearSolver)
    ccall((:SUNLinSolLastFlag, libsundials_sundials), Clong, (SUNLinearSolver,), S)
end

function SUNLinSolSpace(S::SUNLinearSolver, lenrwLS, leniwLS)
    ccall((:SUNLinSolSpace, libsundials_sundials), Cint, (SUNLinearSolver, Ptr{Clong}, Ptr{Clong}), S, lenrwLS, leniwLS)
end

function SUNLinSolFree(S::SUNLinearSolver)
    ccall((:SUNLinSolFree, libsundials_sundials), Cint, (SUNLinearSolver,), S)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_matrix.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function SUNMatGetID(A::SUNMatrix)
    ccall((:SUNMatGetID, libsundials_sundials), SUNMatrix_ID, (SUNMatrix,), A)
end

function SUNMatClone(A::SUNMatrix)
    ccall((:SUNMatClone, libsundials_sundials), SUNMatrix, (SUNMatrix,), A)
end

function SUNMatDestroy(A::SUNMatrix)
    ccall((:SUNMatDestroy, libsundials_sundials), Void, (SUNMatrix,), A)
end

function SUNMatZero(A::SUNMatrix)
    ccall((:SUNMatZero, libsundials_sundials), Cint, (SUNMatrix,), A)
end

function SUNMatCopy(A::SUNMatrix, B::SUNMatrix)
    ccall((:SUNMatCopy, libsundials_sundials), Cint, (SUNMatrix, SUNMatrix), A, B)
end

function SUNMatScaleAdd(c::realtype, A::SUNMatrix, B::SUNMatrix)
    ccall((:SUNMatScaleAdd, libsundials_sundials), Cint, (realtype, SUNMatrix, SUNMatrix), c, A, B)
end

function SUNMatScaleAddI(c::realtype, A::SUNMatrix)
    ccall((:SUNMatScaleAddI, libsundials_sundials), Cint, (realtype, SUNMatrix), c, A)
end

function __SUNMatMatvec(A::SUNMatrix, x::N_Vector, y::N_Vector)
    ccall((:SUNMatMatvec, libsundials_sundials), Cint, (SUNMatrix, N_Vector, N_Vector), A, x, y)
end

function SUNMatMatvec(A, x, y)
    __x = convert(NVector, x)
    __y = convert(NVector, y)
    __SUNMatMatvec(A, convert(N_Vector, __x), convert(N_Vector, __y))
end

function SUNMatSpace(A::SUNMatrix, lenrw, leniw)
    ccall((:SUNMatSpace, libsundials_sundials), Cint, (SUNMatrix, Ptr{Clong}, Ptr{Clong}), A, lenrw, leniw)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_nvector.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __N_VGetVectorID(w::N_Vector)
    ccall((:N_VGetVectorID, libsundials_sundials), N_Vector_ID, (N_Vector,), w)
end

function N_VGetVectorID(w)
    __w = convert(NVector, w)
    __N_VGetVectorID(convert(N_Vector, __w))
end

function __N_VClone(w::N_Vector)
    ccall((:N_VClone, libsundials_sundials), N_Vector, (N_Vector,), w)
end

function N_VClone(w)
    __w = convert(NVector, w)
    __N_VClone(convert(N_Vector, __w))
end

function __N_VCloneEmpty(w::N_Vector)
    ccall((:N_VCloneEmpty, libsundials_sundials), N_Vector, (N_Vector,), w)
end

function N_VCloneEmpty(w)
    __w = convert(NVector, w)
    __N_VCloneEmpty(convert(N_Vector, __w))
end

function __N_VDestroy(v::N_Vector)
    ccall((:N_VDestroy, libsundials_sundials), Void, (N_Vector,), v)
end

function N_VDestroy(v)
    __v = convert(NVector, v)
    __N_VDestroy(convert(N_Vector, __v))
end

function __N_VSpace(v::N_Vector, lrw, liw)
    ccall((:N_VSpace, libsundials_sundials), Void, (N_Vector, Ptr{sunindextype}, Ptr{sunindextype}), v, lrw, liw)
end

function N_VSpace(v, lrw, liw)
    __v = convert(NVector, v)
    __N_VSpace(convert(N_Vector, __v), lrw, liw)
end

function __N_VGetArrayPointer(v::N_Vector)
    ccall((:N_VGetArrayPointer, libsundials_sundials), Ptr{realtype}, (N_Vector,), v)
end

function N_VGetArrayPointer(v)
    __v = convert(NVector, v)
    __N_VGetArrayPointer(convert(N_Vector, __v))
end

function __N_VSetArrayPointer(v_data, v::N_Vector)
    ccall((:N_VSetArrayPointer, libsundials_sundials), Void, (Ptr{realtype}, N_Vector), v_data, v)
end

function N_VSetArrayPointer(v_data, v)
    __v = convert(NVector, v)
    __N_VSetArrayPointer(v_data, convert(N_Vector, __v))
end

function __N_VLinearSum(a::realtype, x::N_Vector, b::realtype, y::N_Vector, z::N_Vector)
    ccall((:N_VLinearSum, libsundials_sundials), Void, (realtype, N_Vector, realtype, N_Vector, N_Vector), a, x, b, y, z)
end

function N_VLinearSum(a, x, b, y, z)
    __x = convert(NVector, x)
    __y = convert(NVector, y)
    __z = convert(NVector, z)
    __N_VLinearSum(a, convert(N_Vector, __x), b, convert(N_Vector, __y), convert(N_Vector, __z))
end

function __N_VConst(c::realtype, z::N_Vector)
    ccall((:N_VConst, libsundials_sundials), Void, (realtype, N_Vector), c, z)
end

function N_VConst(c, z)
    __z = convert(NVector, z)
    __N_VConst(c, convert(N_Vector, __z))
end

function __N_VProd(x::N_Vector, y::N_Vector, z::N_Vector)
    ccall((:N_VProd, libsundials_sundials), Void, (N_Vector, N_Vector, N_Vector), x, y, z)
end

function N_VProd(x, y, z)
    __x = convert(NVector, x)
    __y = convert(NVector, y)
    __z = convert(NVector, z)
    __N_VProd(convert(N_Vector, __x), convert(N_Vector, __y), convert(N_Vector, __z))
end

function __N_VDiv(x::N_Vector, y::N_Vector, z::N_Vector)
    ccall((:N_VDiv, libsundials_sundials), Void, (N_Vector, N_Vector, N_Vector), x, y, z)
end

function N_VDiv(x, y, z)
    __x = convert(NVector, x)
    __y = convert(NVector, y)
    __z = convert(NVector, z)
    __N_VDiv(convert(N_Vector, __x), convert(N_Vector, __y), convert(N_Vector, __z))
end

function __N_VScale(c::realtype, x::N_Vector, z::N_Vector)
    ccall((:N_VScale, libsundials_sundials), Void, (realtype, N_Vector, N_Vector), c, x, z)
end

function N_VScale(c, x, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    __N_VScale(c, convert(N_Vector, __x), convert(N_Vector, __z))
end

function __N_VAbs(x::N_Vector, z::N_Vector)
    ccall((:N_VAbs, libsundials_sundials), Void, (N_Vector, N_Vector), x, z)
end

function N_VAbs(x, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    __N_VAbs(convert(N_Vector, __x), convert(N_Vector, __z))
end

function __N_VInv(x::N_Vector, z::N_Vector)
    ccall((:N_VInv, libsundials_sundials), Void, (N_Vector, N_Vector), x, z)
end

function N_VInv(x, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    __N_VInv(convert(N_Vector, __x), convert(N_Vector, __z))
end

function __N_VAddConst(x::N_Vector, b::realtype, z::N_Vector)
    ccall((:N_VAddConst, libsundials_sundials), Void, (N_Vector, realtype, N_Vector), x, b, z)
end

function N_VAddConst(x, b, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    __N_VAddConst(convert(N_Vector, __x), b, convert(N_Vector, __z))
end

function __N_VDotProd(x::N_Vector, y::N_Vector)
    ccall((:N_VDotProd, libsundials_sundials), realtype, (N_Vector, N_Vector), x, y)
end

function N_VDotProd(x, y)
    __x = convert(NVector, x)
    __y = convert(NVector, y)
    __N_VDotProd(convert(N_Vector, __x), convert(N_Vector, __y))
end

function __N_VMaxNorm(x::N_Vector)
    ccall((:N_VMaxNorm, libsundials_sundials), realtype, (N_Vector,), x)
end

function N_VMaxNorm(x)
    __x = convert(NVector, x)
    __N_VMaxNorm(convert(N_Vector, __x))
end

function __N_VWrmsNorm(x::N_Vector, w::N_Vector)
    ccall((:N_VWrmsNorm, libsundials_sundials), realtype, (N_Vector, N_Vector), x, w)
end

function N_VWrmsNorm(x, w)
    __x = convert(NVector, x)
    __w = convert(NVector, w)
    __N_VWrmsNorm(convert(N_Vector, __x), convert(N_Vector, __w))
end

function __N_VWrmsNormMask(x::N_Vector, w::N_Vector, id::N_Vector)
    ccall((:N_VWrmsNormMask, libsundials_sundials), realtype, (N_Vector, N_Vector, N_Vector), x, w, id)
end

function N_VWrmsNormMask(x, w, id)
    __x = convert(NVector, x)
    __w = convert(NVector, w)
    __id = convert(NVector, id)
    __N_VWrmsNormMask(convert(N_Vector, __x), convert(N_Vector, __w), convert(N_Vector, __id))
end

function __N_VMin(x::N_Vector)
    ccall((:N_VMin, libsundials_sundials), realtype, (N_Vector,), x)
end

function N_VMin(x)
    __x = convert(NVector, x)
    __N_VMin(convert(N_Vector, __x))
end

function __N_VWL2Norm(x::N_Vector, w::N_Vector)
    ccall((:N_VWL2Norm, libsundials_sundials), realtype, (N_Vector, N_Vector), x, w)
end

function N_VWL2Norm(x, w)
    __x = convert(NVector, x)
    __w = convert(NVector, w)
    __N_VWL2Norm(convert(N_Vector, __x), convert(N_Vector, __w))
end

function __N_VL1Norm(x::N_Vector)
    ccall((:N_VL1Norm, libsundials_sundials), realtype, (N_Vector,), x)
end

function N_VL1Norm(x)
    __x = convert(NVector, x)
    __N_VL1Norm(convert(N_Vector, __x))
end

function __N_VCompare(c::realtype, x::N_Vector, z::N_Vector)
    ccall((:N_VCompare, libsundials_sundials), Void, (realtype, N_Vector, N_Vector), c, x, z)
end

function N_VCompare(c, x, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    __N_VCompare(c, convert(N_Vector, __x), convert(N_Vector, __z))
end

function __N_VInvTest(x::N_Vector, z::N_Vector)
    ccall((:N_VInvTest, libsundials_sundials), Cint, (N_Vector, N_Vector), x, z)
end

function N_VInvTest(x, z)
    __x = convert(NVector, x)
    __z = convert(NVector, z)
    __N_VInvTest(convert(N_Vector, __x), convert(N_Vector, __z))
end

function __N_VConstrMask(c::N_Vector, x::N_Vector, m::N_Vector)
    ccall((:N_VConstrMask, libsundials_sundials), Cint, (N_Vector, N_Vector, N_Vector), c, x, m)
end

function N_VConstrMask(c, x, m)
    __c = convert(NVector, c)
    __x = convert(NVector, x)
    __m = convert(NVector, m)
    __N_VConstrMask(convert(N_Vector, __c), convert(N_Vector, __x), convert(N_Vector, __m))
end

function __N_VMinQuotient(num::N_Vector, denom::N_Vector)
    ccall((:N_VMinQuotient, libsundials_sundials), realtype, (N_Vector, N_Vector), num, denom)
end

function N_VMinQuotient(num, denom)
    __num = convert(NVector, num)
    __denom = convert(NVector, denom)
    __N_VMinQuotient(convert(N_Vector, __num), convert(N_Vector, __denom))
end

function __N_VCloneEmptyVectorArray(count::Cint, w::N_Vector)
    ccall((:N_VCloneEmptyVectorArray, libsundials_sundials), Ptr{N_Vector}, (Cint, N_Vector), count, w)
end

function N_VCloneEmptyVectorArray(count, w)
    __w = convert(NVector, w)
    __N_VCloneEmptyVectorArray(convert(Cint, count), convert(N_Vector, __w))
end

function __N_VCloneVectorArray(count::Cint, w::N_Vector)
    ccall((:N_VCloneVectorArray, libsundials_sundials), Ptr{N_Vector}, (Cint, N_Vector), count, w)
end

function N_VCloneVectorArray(count, w)
    __w = convert(NVector, w)
    __N_VCloneVectorArray(convert(Cint, count), convert(N_Vector, __w))
end

function __N_VDestroyVectorArray(vs, count::Cint)
    ccall((:N_VDestroyVectorArray, libsundials_sundials), Void, (Ptr{N_Vector}, Cint), vs, count)
end

function N_VDestroyVectorArray(vs, count)
    __N_VDestroyVectorArray(vs, convert(Cint, count))
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_pcg.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_spbcgs.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_spfgmr.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_spgmr.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_sptfqmr.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_types.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_fconfig.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_fnvector.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_math.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __SUNRpowerI(base::realtype, exponent::Cint)
    ccall((:SUNRpowerI, libsundials_sundials), realtype, (realtype, Cint), base, exponent)
end

function SUNRpowerI(base, exponent)
    __SUNRpowerI(base, convert(Cint, exponent))
end

function SUNRpowerR(base::realtype, exponent::realtype)
    ccall((:SUNRpowerR, libsundials_sundials), realtype, (realtype, realtype), base, exponent)
end

function SUNRsqrt(x::realtype)
    ccall((:SUNRsqrt, libsundials_sundials), realtype, (realtype,), x)
end

function SUNRabs(x::realtype)
    ccall((:SUNRabs, libsundials_sundials), realtype, (realtype,), x)
end

function SUNRexp(x::realtype)
    ccall((:SUNRexp, libsundials_sundials), realtype, (realtype,), x)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_sparse.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/sundials/sundials_version.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function __SUNDIALSGetVersion(version, len::Cint)
    ccall((:SUNDIALSGetVersion, libsundials_sundials), Cint, (Cstring, Cint), version, len)
end

function SUNDIALSGetVersion(version, len)
    __SUNDIALSGetVersion(version, convert(Cint, len))
end

function __SUNDIALSGetVersionNumber(major, minor, patch, label, len::Cint)
    ccall((:SUNDIALSGetVersionNumber, libsundials_sundials), Cint, (Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cstring, Cint), major, minor, patch, label, len)
end

function SUNDIALSGetVersionNumber(major, minor, patch, label, len)
    __SUNDIALSGetVersionNumber(major, minor, patch, label, convert(Cint, len))
end


text:
# Check to make sure types are imported properly
@test MultiSelectMenu <: TerminalMenus.AbstractMenu

# Invalid Menu Params
@test_throws ErrorException MultiSelectMenu(["one"])
@test_throws ErrorException MultiSelectMenu(["one", "two", "three"], pagesize=1)

# Constructor
@test MultiSelectMenu(["one", "two", "three"]).pagesize == 3
@test MultiSelectMenu(string.(1:30), pagesize=-1).pagesize == 30
@test MultiSelectMenu(string.(1:4), pagesize=10).pagesize == 4
@test MultiSelectMenu(string.(1:100)).pagesize == 10

multi_menu = MultiSelectMenu(string.(1:20))
@test TerminalMenus.options(multi_menu) == string.(1:20)
@test TerminalMenus.header(multi_menu) == "[press: d=done, a=all, n=none]"

# Output
TerminalMenus.config() # Use default chars
CONFIG = TerminalMenus.CONFIG

multi_menu = MultiSelectMenu(string.(1:10))
buf = IOBuffer()
TerminalMenus.writeLine(buf, multi_menu, 1, true)
@test String(take!(buf)) == string(CONFIG[:cursor], " ", CONFIG[:unchecked], " 1")
TerminalMenus.config(cursor='+')
TerminalMenus.writeLine(buf, multi_menu, 1, true)
@test String(take!(buf)) == string("+ ", CONFIG[:unchecked], " 1")
TerminalMenus.config(charset=:unicode)
TerminalMenus.writeLine(buf, multi_menu, 1, true)
@test String(take!(buf)) == string(CONFIG[:cursor], " ", CONFIG[:unchecked], " 1")

# Test SDTIN
multi_menu = MultiSelectMenu(string.(1:10))
@test simulateInput(Set([1,2]), multi_menu, :enter, :down, :enter, 'd')

max_stars_repo_name: djsegal/julia-fork

max_issues_repo_name: djsegal/julia-fork

max_forks_repo_name: djsegal/julia-fork


text:
function PlotVortices(xi,yi,nplus,nminus,Lx,Ly,R)

n = nplus+nminus;

xp = xi(1:nplus);
yp = yi(1:nplus);
xn = xi(nplus+1:n);
yn = yi(nplus+1:n);

plot(xp,yp,'.r','Markersize',10)
hold on
plot(xn,yn,'.b','Markersize',10)
xlim([-Lx/2 Lx/2])
ylim([-Ly/2 Ly/2])
set(gca,'plotboxaspectratio',[Lx/Ly 1 1])
if nargin == 7
    theta= 0:pi/100:2*pi;
    u = R*cos(theta);
    v = R*sin(theta);
    plot(u,v,'-k','Linewidth',1)
end


text:
function companion_matrix(c::Vector{T}) where T
    n=length(c)-1

    if n==0
        zeros(T,0,0)
    else
        A=zeros(T,n,n)
        for k=1:n
            A[k,end]=-c[k]/c[end]
        end
        for k=2:n
            A[k,k-1]=one(T)
        end
        A
    end
end

# if isdir(Pkg.dir("AMVW"))
#     using AMVW
#     function complexroots(cfs::Vector)
#         c=chop(cfs,10eps())
#
#         # Only use special routine for large roots
#         if length(c)≥70
#             Main.AMVW.rootsAMVW(c)
#         else
#             hesseneigvals(companion_matrix(c))
#         end
#     end
# else
complexroots(cfs::Vector{T}) where {T<:Union{Float64,ComplexF64}} =
    hesseneigvals(companion_matrix(chop(cfs,10eps())))
# end

function complexroots(cfs::Vector{T}) where T<:Union{BigFloat,Complex{BigFloat}}
    a = Fun(Taylor(Circle(BigFloat)),cfs)
    ap = a'
    rts = Array{Complex{BigFloat}}(complexroots(Vector{ComplexF64}(cfs)))
    # Do 3 Newton steps
    for _ = 1:3
        rts .-= a.(rts)./ap.(rts)
    end
    rts
end

complexroots(neg::Vector, pos::Vector) =
    complexroots([reverse(chop(neg,10eps()), dims=1);pos])
complexroots(f::Fun{Laurent{DD,RR}}) where {DD,RR} =
    mappoint.(Ref(Circle()), Ref(domain(f)),
        complexroots(f.coefficients[2:2:end],f.coefficients[1:2:end]))
complexroots(f::Fun{Taylor{DD,RR}}) where {DD,RR} =
    mappoint.(Ref(Circle()), Ref(domain(f)), complexroots(f.coefficients))

function roots(f::Fun{Laurent{DD,RR}}) where {DD,RR}
    irts=filter!(z->in(z,Circle()),complexroots(Fun(Laurent(Circle()),f.coefficients)))
    if length(irts)==0
        Complex{Float64}[]
    else
        rts=fromcanonical.(f, tocanonical.(Ref(Circle()), irts))
        if isa(domain(f),PeriodicSegment)
            sort!(real(rts))  # Make type safe?
        else
            rts
        end
    end
end

roots(f::Fun{Fourier{D,R}}) where {D,R} = roots(Fun(f,Laurent))


text:
# ---
# title: Dopamine\_IQN\_Atari(breakout)
# cover: assets/Dopamine_IQN_Atari_breakout_evaluating_avg_score.svg
# description: Use the IQN to play the atari game breakout.
# date: 2021-05-22
# author: "[Jun Tian](https://github.com/findmyway)"
# ---

# This experiment tries to use the same config in [google/dopamine](https://github.com/google/dopamine/blob/master/dopamine/agents/implicit_quantile/configs/implicit_quantile.gin) to run the atari games with IQN, except the following two major differences:

# - We use the `BSpline(Linear())` instead of `cv2.INTER_AREA` method to resize the image.
# - `ADAM` in Flux.jl do not support setting `epsilon`. (This should be a minor issue.)

# On a machine with a Nvidia 2080Ti GPU card, the training speed of this
# experiment is about **138 steps/sec**. The testing speed about **695
# steps/sec**. For comparison, the training speed of dopamine is about **79
# steps/sec**.

# Following are some basic stats. The evaluation result seems to be aligned with
# the result reported in
# [dopamine](https://github.com/google/dopamine/blob/2a7d91d283/baselines/data/breakout.json).

# Average reward per episode in evaluation mode:
# ![](assets/Dopamine_IQN_Atari_breakout_evaluating_avg_score.svg)

# Average episode length in evaluation mode:
# ![](assets/Dopamine_IQN_Atari_breakout_evaluating_avg_length.svg)

# Average episode length in training mode:
# ![](assets/Dopamine_IQN_Atari_breakout_training_episode_length.svg)

# Training loss per updated:
# ![](assets/Dopamine_IQN_Atari_breakout_training_loss.svg)

# Reward per episode in training mode:
# ![](assets/Dopamine_IQN_Atari_breakout_training_reward.svg)

#+ tangle=true
using ReinforcementLearning
using ArcadeLearningEnvironment
using CUDA
using Flux
using Flux.Losses: huber_loss
using Dates
using Random
using Setfield
using Statistics
using Logging
using TensorBoardLogger

#+ tangle=false
## START TODO: move into a common file
using ImageTransformations: imresize!

"""
    ResizeImage(img::Array{T, N})
    ResizeImage(dims::Int...) -> ResizeImage(Float32, dims...)
    ResizeImage(T::Type{<:Number}, dims::Int...)

By default the `BSpline(Linear())`` method is used to resize the `state` field
of an observation to size of `img` (or `dims`). In some other packages, people
use the
[`cv2.INTER_AREA`](https://github.com/google/dopamine/blob/2a7d91d2831ca28cea0d3b0f4d5c7a7107e846ab/dopamine/discrete_domains/atari_lib.py#L511-L513),
which is not supported in `ImageTransformations.jl` yet.
"""
struct ResizeImage{T,N}
    img::Array{T,N}
end

ResizeImage(dims::Int...) = ResizeImage(Float32, dims...)
ResizeImage(T::Type{<:Number}, dims::Int...) = ResizeImage(Array{T}(undef, dims))

function (p::ResizeImage)(state::AbstractArray)
    imresize!(p.img, state)
    p.img
end

function atari_env_factory(
    name,
    state_size,
    n_frames,
    max_episode_steps = 100_000;
    seed = nothing,
    repeat_action_probability = 0.25,
    n_replica = nothing,
)
    init(seed) =
        RewardTransformedEnv(
            StateCachedEnv(
                StateTransformedEnv(
                    AtariEnv(;
                        name = string(name),
                        grayscale_obs = true,
                        noop_max = 30,
                        frame_skip = 4,
                        terminal_on_life_loss = false,
                        repeat_action_probability = repeat_action_probability,
                        max_num_frames_per_episode = n_frames * max_episode_steps,
                        color_averaging = false,
                        full_action_space = false,
                        seed = seed,
                    );
                    state_mapping=Chain(
                        ResizeImage(state_size...),
                        StackFrames(state_size..., n_frames)
                    ),
                    state_space_mapping= _ -> Space(fill(0..256, state_size..., n_frames))
                )
            );
            reward_mapping = r -> clamp(r, -1, 1)
        )

    if isnothing(n_replica)
        init(seed)
    else
        envs = [
            init(isnothing(seed) ? nothing : hash(seed + i))
            for i in 1:n_replica
        ]
        states = Flux.batch(state.(envs))
        rewards = reward.(envs)
        terminals = is_terminated.(envs)
        A = Space([action_space(x) for x in envs])
        S = Space(fill(0..255, size(states)))
        MultiThreadEnv(envs, states, rewards, terminals, A, S, nothing)
    end
end

"Total reward per episode before reward reshaping"
Base.@kwdef mutable struct TotalOriginalRewardPerEpisode <: AbstractHook
    rewards::Vector{Float64} = Float64[]
    reward::Float64 = 0.0
end

function (hook::TotalOriginalRewardPerEpisode)(
    ::PostActStage,
    agent,
    env::RewardTransformedEnv,
)
    hook.reward += reward(env.env)
end

function (hook::TotalOriginalRewardPerEpisode)(::PostEpisodeStage, agent, env)
    push!(hook.rewards, hook.reward)
    hook.reward = 0
end

"Total reward of each inner env per episode before reward reshaping"
struct TotalBatchOriginalRewardPerEpisode <: AbstractHook
    rewards::Vector{Vector{Float64}}
    reward::Vector{Float64}
end

function TotalBatchOriginalRewardPerEpisode(batch_size::Int)
    TotalBatchOriginalRewardPerEpisode([Float64[] for _ in 1:batch_size], zeros(batch_size))
end

function (hook::TotalBatchOriginalRewardPerEpisode)(
    ::PostActStage,
    agent,
    env::MultiThreadEnv{<:RewardTransformedEnv},
)
    for (i, e) in enumerate(env.envs)
        hook.reward[i] += reward(e.env)
        if is_terminated(e)
            push!(hook.rewards[i], hook.reward[i])
            hook.reward[i] = 0.0
        end
    end
end

## END TODO: move into a common file

#+ tangle=true
function RL.Experiment(
    ::Val{:Dopamine},
    ::Val{:IQN},
    ::Val{:Atari},
    name::AbstractString;
    save_dir = nothing,
    seed = nothing,
)
    rng = Random.GLOBAL_RNG
    Random.seed!(rng, seed)
    device_rng = CUDA.functional() ? CUDA.CURAND.RNG() : rng
    Random.seed!(device_rng, isnothing(seed) ? nothing : hash(seed + 1))

    if isnothing(save_dir)
        t = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
        save_dir = joinpath(pwd(), "checkpoints", "dopamine_IQN_atari_$(name)_$(t)")
    end

    lg = TBLogger(joinpath(save_dir, "tb_log"), min_level = Logging.Info)

    N_FRAMES = 4
    STATE_SIZE = (84, 84)

    env = atari_env_factory(name, STATE_SIZE, N_FRAMES; seed = isnothing(seed) ? nothing : hash(seed + 2))
    N_ACTIONS = length(action_space(env))
    Nₑₘ = 64

    init = glorot_uniform(rng)

    create_model() =
        ImplicitQuantileNet(
            ψ = Chain(
                x -> x ./ 255,
                CrossCor((8, 8), N_FRAMES => 32, relu; stride = 4, pad = 2, init = init),
                CrossCor((4, 4), 32 => 64, relu; stride = 2, pad = 2, init = init),
                CrossCor((3, 3), 64 => 64, relu; stride = 1, pad = 1, init = init),
                x -> reshape(x, :, size(x)[end]),
            ),
            ϕ = Dense(Nₑₘ, 11 * 11 * 64, relu; init = init),
            header = Chain(
                Dense(11 * 11 * 64, 512, relu; init = init),
                Dense(512, N_ACTIONS; init = init),
            ),
        ) |> gpu

    agent = Agent(
        policy = QBasedPolicy(
            learner = IQNLearner(
                approximator = NeuralNetworkApproximator(
                    model = create_model(),
                    optimizer = ADAM(0.00005),  # epsilon is not set here
                ),
                target_approximator = NeuralNetworkApproximator(model = create_model()),
                κ = 1.0f0,
                N = 64,
                N′ = 64,
                Nₑₘ = Nₑₘ,
                K = 32,
                γ = 0.99f0,
                stack_size = 4,
                batch_size = 32,
                update_horizon = 3,
                min_replay_history = 20_000,
                update_freq = 4,
                target_update_freq = 8_000,
                default_priority = 1.0f2,  # only valid when trajectory is CircularArrayPSARTTrajectory
                rng = rng,
                device_rng = device_rng,
            ),
            explorer = EpsilonGreedyExplorer(
                ϵ_init = 1.0,
                ϵ_stable = 0.01,
                decay_steps = 250_000,
                kind = :linear,
                rng = rng,
            ),
        ),
        trajectory = CircularArraySARTTrajectory(
            capacity = haskey(ENV, "CI") : 1_000 : 1_000_000,
            state = Matrix{Float32} => STATE_SIZE,
        ),
    )

    EVALUATION_FREQ = 250_000
    MAX_EPISODE_STEPS_EVAL = 27_000

    total_reward_per_episode = TotalOriginalRewardPerEpisode()
    steps_per_episode = StepsPerEpisode()
    hook = ComposedHook(
        total_reward_per_episode,
        steps_per_episode,
        DoEveryNStep() do t, agent, env
            with_logger(lg) do
                @info "training" loss = agent.policy.learner.loss
            end
        end,
        DoEveryNEpisode() do t, agent, env
            with_logger(lg) do
                @info "training" reward = total_reward_per_episode.rewards[end] episode_length =
                    steps_per_episode.steps[end] log_step_increment = 0
            end
        end,
        DoEveryNStep(;n=EVALUATION_FREQ) do t, agent, env
            @info "evaluating agent at $t step..."
            p = agent.policy
            p = @set p.explorer = EpsilonGreedyExplorer(0.001; rng = rng)  # set evaluation epsilon
            h = ComposedHook(TotalOriginalRewardPerEpisode(), StepsPerEpisode())
            s = @elapsed run(
                p,
                atari_env_factory(
                    name,
                    STATE_SIZE,
                    N_FRAMES,
                    MAX_EPISODE_STEPS_EVAL;
                    seed = isnothing(seed) ? nothing : hash(seed + t)
                ),
                StopAfterStep(125_000; is_show_progress = false),
                h,
            )

            avg_score = mean(h[1].rewards[1:end-1])
            avg_length = mean(h[2].steps[1:end-1])

            @info "finished evaluating agent in $s seconds" avg_length = avg_length avg_score = avg_score
            with_logger(lg) do
                @info "evaluating" avg_length = avg_length avg_score = avg_score log_step_increment = 0
            end
        end,
    )

    stop_condition = StopAfterStep(
        haskey(ENV, "CI") ? 10_000 : 50_000_000,
        is_show_progress=!haskey(ENV, "CI")
    )
    Experiment(agent, env, stop_condition, hook, "# IQN <-> Atari($name)")
end

#+ tangle=false
ex = E`Dopamine_IQN_Atari(breakout)`
run(ex)


text:
# Define shortcuts for Tcl/Tk methods.
#
# Usage:
#
#     using Tcl
#     using Tcl.ShortNames
#

module ShortNames

using ...Tcl

export
    cget,
    configure,
    getinterp,
    getparent,
    getpath,
    grid,
    list,
    pack,
    place,

    Button,
    Canvas,
    Checkbutton,
    Entry,
    Frame,
    Label,
    Labelframe,
    Listbox,
    Menu,
    Menubutton,
    Message,
    Panedwindow,
    Radiobutton,
    Scale,
    Scrollbar,
    Spinbox,
    TButton,
    TCheckbutton,
    TCombobox,
    TEntry,
    TFrame,
    TLabel,
    TLabelframe,
    TMenubutton,
    TNotebook,
    TPanedwindow,
    TProgressbar,
    TRadiobutton,
    TScale,
    TScrollbar,
    TSeparator,
    TSizegrip,
    TSpinbox,
    Text,
    Toplevel,
    Treeview

const cget          = Tcl.cget
const configure     = Tcl.configure
const getinterp     = Tcl.getinterp
const getparent     = Tcl.getparent
const getpath       = Tcl.getpath
const grid          = Tcl.grid
const list          = Tcl.list
const pack          = Tcl.pack
const place         = Tcl.place

# Use the same short names as the Tk class names given by `winfo class $w`.
const Button        = TkButton
const Canvas        = TkCanvas
const Checkbutton   = TkCheckbutton
const Entry         = TkEntry
const Frame         = TkFrame
const Label         = TkLabel
const Labelframe    = TkLabelframe
const Listbox       = TkListbox
const Menu          = TkMenu
const Menubutton    = TkMenubutton
const Message       = TkMessage
const Panedwindow   = TkPanedwindow
const Radiobutton   = TkRadiobutton
const Scale         = TkScale
const Scrollbar     = TkScrollbar
const Spinbox       = TkSpinbox
const TButton       = TtkButton
const TCheckbutton  = TtkCheckbutton
const TCombobox     = TtkCombobox
const TEntry        = TtkEntry
const TFrame        = TtkFrame
const TLabel        = TtkLabel
const TLabelframe   = TtkLabelframe
const TMenubutton   = TtkMenubutton
const TNotebook     = TtkNotebook
const TPanedwindow  = TtkPanedwindow
const TProgressbar  = TtkProgressbar
const TRadiobutton  = TtkRadiobutton
const TScale        = TtkScale
const TScrollbar    = TtkScrollbar
const TSeparator    = TtkSeparator
const TSizegrip     = TtkSizegrip
const TSpinbox      = TtkSpinbox
const Text          = TkText
const Toplevel      = TkToplevel
const Treeview      = TtkTreeview

end #module


text:
module XPT

import Datetime: DateTime, datetime
import Base: show, close
import DataArrays: DataArray, NA
import DataFrames: DataFrame
import StrPack: @struct, calcsize, unpack

export XPTFile, readdf, closexpt

include("util.jl")
include("XPTInfo.jl")
include("NameStr.jl")
include("DSInfo.jl")
include("IBM2Float64.jl")
include("XPTFile.jl")

end # module


text:
"""
    DevicePtr{T,A}
A memory address that refers to data of type `T` that is accessible from the GPU. It is the
on-device counterpart of `MtlDAdrv.MtlPtr`, additionally keeping track of the address space
`A` where the data resides (shared, global, constant, etc). This information is used to
provide optimized implementations of operations such as `unsafe_load` and `unsafe_store!.`
"""
DevicePtr

# constructors
###DevicePtr{T,A}(x::Union{Int,UInt,DevicePtr}) where {T,A<:AddressSpace} = Base.bitcast(DevicePtr{T,A}, x)
#DevicePtr{T}(ptr::MtlPtr{T}) where {T} = DevicePtr{T,AS.Generic}(ptr)
#DevicePtr(ptr::MtlPtr{T}) where {T} = DevicePtr{T,AS.Generic}(ptr)

## conversions
###Base.convert(::Type{DevicePtr{T,A}}, x::Union{Int,UInt}) where {T,A<:AddressSpace} = DevicePtr{T,A}(x)

# between host and device pointers
#Base.convert(::Type{MtlPtr{T}},  p::DevicePtr)  where {T}                 = Base.bitcast(MtlPtr{T}, p)
Base.convert(::Type{MtlBuffer{T}},   p::DevicePtr)  where {T}               = MtlBuffer{T}(Base.bitcast(MTL.MTLBuffer, p))
#Base.convert(::Type{DevicePtr{T,A}}, p::MtlPtr) where {T,A<:AddressSpace}   = Base.bitcast(DevicePtr{T,A}, p)
#Base.convert(::Type{DevicePtr{T}},   p::MtlPtr)   where {T}                 = Base.bitcast(DevicePtr{T,AS.Generic}, p)
#Base.convert(::Type{DevicePtr{T,A}}, p::MtlBuffer) where {T,A<:AddressSpace}   = Base.bitcast(DevicePtr{T,A}, p.handle)
#Base.convert(::Type{DevicePtr{T}},   p::MtlBuffer)   where {T}                 = Base.bitcast(DevicePtr{T,AS.Generic}, p.handle)

# between CPU pointers, for the purpose of working with `ccall`
Base.unsafe_convert(::Type{Ptr{T}}, x::DevicePtr{T}) where {T} = reinterpret(Ptr{T}, x)

# between device pointers
Base.convert(::Type{<:DevicePtr}, p::DevicePtr)                         = throw(ArgumentError("cannot convert between incompatible device pointer types"))
Base.convert(::Type{DevicePtr{T,A}}, p::DevicePtr{T,A})   where {T,A}   = p
Base.unsafe_convert(::Type{DevicePtr{T,A}}, p::DevicePtr) where {T,A}   = Base.bitcast(DevicePtr{T,A}, p)
## identical addrspaces
Base.convert(::Type{DevicePtr{T,A}}, p::DevicePtr{U,A}) where {T,U,A} = Base.unsafe_convert(DevicePtr{T,A}, p)
## convert to & from generic
Base.convert(::Type{DevicePtr{T,AS.Generic}}, p::DevicePtr)               where {T}     = Base.unsafe_convert(DevicePtr{T,AS.Generic}, p)
Base.convert(::Type{DevicePtr{T,A}}, p::DevicePtr{U,AS.Generic})          where {T,U,A} = Base.unsafe_convert(DevicePtr{T,A}, p)
Base.convert(::Type{DevicePtr{T,AS.Generic}}, p::DevicePtr{T,AS.Generic}) where {T}     = p  # avoid ambiguities
## unspecified, preserve source addrspace
Base.convert(::Type{DevicePtr{T}}, p::DevicePtr{U,A}) where {T,U,A} = Base.unsafe_convert(DevicePtr{T,A}, p)

## memory operations

# operand types supported by llvm.nvvm.ldg.global
# NOTE: CUDA 8.0 supports more caching modifiers, but those aren't supported by LLVM yet
const LDGTypes = Union{UInt8, UInt16, UInt32, UInt64,
                       Int8, Int16, Int32, Int64,
                       Float32, Float64}

# interface

export unsafe_cached_load

Base.unsafe_load(p::DevicePtr{T}, i::Integer=1, align::Val=Val(1)) where {T} =
    pointerref(p, i, align)

Base.unsafe_store!(p::DevicePtr{T}, x, i::Integer=1, align::Val=Val(1)) where {T} =
    pointerset(p, convert(T, x), i, align)

# NOTE: fall back to normal pointerref for unsupported types. we could be smarter here,
#       e.g. destruct/load/reconstruct, but that's too complicated for what it's worth.
unsafe_cached_load(p::DevicePtr, i::Integer=1, align::Val=Val(1)) =
    pointerref(p, i, align)


text:
"""ここで実行"""

using CPUTime
using Plots