text:
"""
```
@inline function solve_hjb(V::Union{Array{R},Array{S}}, a_lb::R, γ::R, ddeath::R, pam::R,
                           trans::R, ξ::R, τ_I::R, aggZ::Union{R,S}, w::Union{R,S},
                           r_b_vec::Union{Vector{R},Vector{S}}, y::Vector{U}, a::Vector{R},
                           b::Vector{R}, cost, util, deposit; permanent::Bool=false,
                           Vamin = 0.0, Vbmin = 1e-8) where {R<:Float64, S, U<:Number}
```

Solves Hamilton-Jacobi-Bellman equation.
"""
@inline function solve_hjb(V::Union{Array{R},Array{S}}, a_lb::R, γ::R, ddeath::R, pam::R,
                           trans::R, ξ::R, τ_I::R, aggZ::Union{R,S}, w::Union{R,S},
                           r_b_vec::Union{Vector{R},Vector{S}}, y::Vector{U}, a::Vector{R},
                           b::Vector{R}, cost, util, deposit; permanent::Bool=false,
                           Vamin = 0.0, Vbmin = 1e-8) where {R<:Float64, S, U<:Number}
    #----------------------------------------------------------------
    # HJB Equation
    #----------------------------------------------------------------
    I, J, N = size(V)
    c, s, d = similar(V), similar(V), similar(V)

    perm_c = permanent ? ddeath * pam - aggZ : ddeath * pam
    c0_c   = ((1-ξ) - τ_I) * w

    for i=1:I, j=1:J, n=1:N

        c0 = c0_c * real(y[n]) + b[i] * (r_b_vec[i] + perm_c) + trans

        # ---- Liquid Assets, Forward + Backward Difference ----
        VaF = (j==J) ? 0.0 : max((V[i,j+1,n] - V[i,j,n]) / (a[j+1]-a[j]), Vamin)
        VaB = (j==1) ? 0.0 : max((V[i,j,n] - V[i,j-1,n]) / (a[j]-a[j-1]), Vamin)

        # ---- Illiquid Assets, Forward + Backward Difference ----
        VbF = (i==I) ? 0.0 : max((V[i+1,j,n] - V[i,j,n]) / (b[i+1]-b[i]), Vbmin)
        VbB = (i==1) ? 0.0 : max((V[i,j,n] - V[i-1,j,n]) / (b[i]-b[i-1]), Vbmin)

        # Decisions conditional on a particular direction of derivative
        cF = (i==I) ? 0.0 : VbF ^ (-1 / γ)
        cB = (i==1) ? c0  : VbB ^ (-1 / γ)

        sF = (i==I) ? 0.0 : c0 - cF
        sB = (i==1) ? 0.0 : c0 - cB

        if (i==1 && j > 1) VbB = util(cB) end

        #----------------------------------------------------------------
        # Consumption & Savings Decision
        #----------------------------------------------------------------
        Hc0 =                  util(c0)
        HcF = (i==I) ? -1e12 : util(cF) + VbF * sF
        HcB =                  util(cB) + VbB * sB

        validF = (sF > 0)
        validB = (sB < 0)

        IcF = validF * max(!validB, (HcF >= HcB)) * (HcF >= Hc0)
        IcB = validB * max(!validF, (HcB >= HcF)) * (HcB >= Hc0)
        Ic0 = 1 - IcF - IcB

        c[i,j,n] = IcF * cF + IcB * cB + Ic0 * c0
        s[i,j,n] = IcF * sF + IcB * sB

        #----------------------------------------------------------------
        # Deposit Decision
        #----------------------------------------------------------------
        dFB = (i == 1 || j == J) ? 0.0 : deposit(VaF, VbB, max(a[j], a_lb))
        dBF = (i == I || j == 1) ? 0.0 : deposit(VaB, VbF, max(a[j], a_lb))
        dBB = (j == 1)           ? 0.0 : deposit(VaB, VbB, max(a[j], a_lb))

        HdFB = (i == 1 || j == J) ? -1e12 : VaF * dFB - VbB * (dFB + cost(dFB, a[j]))
        HdBB = (j == 1)           ? -1e12 : VaB * dBB - VbB * (dBB + cost(dBB, a[j]))
        HdBF = (i == I || j == 1) ? -1e12 : VaB * dBF - VbF * (dBF + cost(dBF, a[j]))

        validFB = (HdFB > 0) * (dFB >  0)
        validBB = (HdBB > 0) * (dBB <= 0) * (dBB >  -cost(dBB, a[j]))
        validBF = (HdBF > 0) *              (dBF <= -cost(dBF, a[j]))

        IcFB = max(!validBF, HdFB >= HdBF) * validFB * max(!validBB, HdFB >= HdBB)
        IcBF = max(!validFB, HdBF >= HdFB) * validBF * max(!validBB, HdBF >= HdBB)
        IcBB = max(!validFB, HdBB >= HdFB) * validBB * max(!validBF, HdBB >= HdBF)

        d[i,j,n] = IcFB * dFB + IcBF * dBF + IcBB * dBB
    end
    return c, s, d
end

"""
```
@inline function solve_hjb_lite(V::Union{Array{R},Array{S}}, a_lb::R, γ::R,
                           permanent::Bool, ddeath::R, pam::R, aggZ::Union{R,S}, ξ::R, τ_I::R,
                           w::Union{R,S}, trans::R, r_b::Union{R,S}, r_b_borr::Union{R,S},
                           y::U, a::Vector{R}, b::Vector{R}, cost, util, deposit;
                           Vamin = 0.0, Vbmin = 1e-8) where {R<:Float64, S, U<:Number}
```
Solves Hamilton-Jacobi-Bellman equation.
"""
@inline function solve_hjb_lite(V::Union{Array{R},Array{S}}, a_lb::R, γ::R,
                           permanent::Bool, ddeath::R, pam::R, aggZ::Union{R,S}, ξ::R, τ_I::R,
                           w::Union{R,S}, trans::R, r_b::Union{R,S}, r_b_borr::Union{R,S},
                           y::U, a::Vector{R}, b::Vector{R}, cost, util, deposit;
                           Vamin = 0.0, Vbmin = 1e-8) where {R<:Float64, S, U<:Number}
    #----------------------------------------------------------------
    # HJB Equation
    #----------------------------------------------------------------
    I, J = size(V)

    perm_c = permanent ? ddeath * pam - aggZ : ddeath * pam
    c0_c   = ((1-ξ) - τ_I) * w

    c, s, d = similar(V), similar(V), similar(V)
    for i=1:I, j=1:J
        c0 = c0_c * real(y) + b[i] * (((b[i] >= 0) ? r_b : r_b_borr) + perm_c) + trans

        # ---- Liquid Assets, Forward + Backward Difference ----
        VaF = (j==J) ? 0.0 : max((V[i,j+1] - V[i,j]) / (a[j+1]-a[j]), Vamin)
        VaB = (j==1) ? 0.0 : max((V[i,j] - V[i,j-1]) / (a[j]-a[j-1]), Vamin)

        # ---- Illiquid Assets, Forward + Backward Difference ----
        VbF = (i==I) ? 0.0 : max((V[i+1,j] - V[i,j]) / (b[i+1]-b[i]), Vbmin)
        VbB = (i==1) ? 0.0 : max((V[i,j] - V[i-1,j]) / (b[i]-b[i-1]), Vbmin)

        # Decisions conditional on a particular direction of derivative
        cF = (i==I) ? 0.0 : VbF ^ (-1 / γ)
        cB = (i==1) ? c0  : VbB ^ (-1 / γ)

        sF = (i==I) ? 0.0 : c0 - cF
        sB = (i==1) ? 0.0 : c0 - cB

        if (i==1 && j > 1) VbB = util(cB) end

        #----------------------------------------------------------------
        # Consumption & Savings Decision
        #----------------------------------------------------------------
        Hc0 =                  util(c0)
        HcF = (i==I) ? -1e12 : util(cF) + VbF * sF
        HcB =                  util(cB) + VbB * sB

        validF = (sF > 0)
        validB = (sB < 0)

        IcF = validF * max(!validB, (HcF >= HcB)) * (HcF >= Hc0)
        IcB = validB * max(!validF, (HcB >= HcF)) * (HcB >= Hc0)
        Ic0 = 1 - IcF - IcB

        c[i,j] = IcF * cF + IcB * cB + Ic0 * c0
        s[i,j] = IcF * sF + IcB * sB

        #----------------------------------------------------------------
        # Deposit Decision
        #----------------------------------------------------------------
        dFB = (i == 1 || j == J) ? 0.0 : deposit(VaF, VbB, max(a[j], a_lb))
        dBF = (i == I || j == 1) ? 0.0 : deposit(VaB, VbF, max(a[j], a_lb))
        dBB = (j == 1)           ? 0.0 : deposit(VaB, VbB, max(a[j], a_lb))

        HdFB = (i == 1 || j == J) ? -1e12 : VaF * dFB - VbB * (dFB + cost(dFB, a[j]))
        HdBB = (j == 1)           ? -1e12 : VaB * dBB - VbB * (dBB + cost(dBB, a[j]))
        HdBF = (i == I || j == 1) ? -1e12 : VaB * dBF - VbF * (dBF + cost(dBF, a[j]))

        validFB = (HdFB > 0) * (dFB >  0)
        validBB = (HdBB > 0) * (dBB <= 0) * (dBB >  -cost(dBB, a[j]))
        validBF = (HdBF > 0) *              (dBF <= -cost(dBF, a[j]))

        IcFB = max(!validBF, HdFB >= HdBF) * validFB * max(!validBB, HdFB >= HdBB)
        IcBF = max(!validFB, HdBF >= HdFB) * validBF * max(!validBB, HdBF >= HdBB)
        IcBB = max(!validFB, HdBB >= HdFB) * validBB * max(!validBF, HdBB >= HdBF)

        d[i,j] = IcFB * dFB + IcBF * dBF + dBB * IcBB
    end
    return c, s, d
end

@inline function solve_kfe(dab_g_tilde, b, g_inds, gg, ddeath, lambda, AT, indx, V_SS)
    my_inds = findall(x -> x != indx, 1:size(lambda,1))

    loc = findall(b .== 0)
    dab_g_tilde_inv       = vec(1.0 ./ dab_g_tilde)
    dab_g_small           = dab_g_tilde ./ dab_g_tilde[loc] * ddeath
    dab_g_small[loc]     .= 0.0
    death_process         = -ddeath * my_speye(length(gg))
    death_process[loc,:]  = vec(dab_g_small)
    gIntermediate = dab_g_tilde_inv .* ((AT * (dab_g_tilde .* gg)) +
                        sum([lambda[kk,indx] .* dab_g_tilde .* g_inds[:,kk] for kk in my_inds])) +
                                             death_process * gg
    cc = sum([lambda[indx,kk]*(V_SS[:,kk] - V_SS[:,indx]) for kk in my_inds])
    return gIntermediate, cc
end

@inline function add_diag(A::SparseMatrixCSC, c::AbstractFloat)
    for i=1:size(A,1)
        A[i,i] += c
    end
    return A
end

@inline function mul_diag(A::SparseMatrixCSC, c::AbstractFloat)
    for i=1:size(A,1)
        A[i,i] *= c
    end
    return A
end

"""
```
@inline function update_value_fn(A::SparseMatrixCSC, Vn::Union{Array{T,3},Array{R,3}},
                                 lambda::Matrix{R},   u::Union{Array{T,3},Array{R,3}},
                                 Delta::R, rrho::R, ddeath::R) where {R<:Float64, T<:Complex}
```
Update value function.
"""
@inline function update_value_fn(A::SparseMatrixCSC, Vn::Union{Array{T,3},Array{R,3}},
                                 lambda::Matrix{R},   u::Union{Array{T,3},Array{R,3}},
                                 Delta::R, rrho::R, ddeath::R) where {R<:Float64, T<:Complex}
    I, J, N = size(Vn)
    Vn_new = Array{Float64,3}(undef, I, J, N)
    for kk = 1:N
        indx_k = 1:N .!= kk
        Ak = A[1+(kk-1)*(I*J):kk*(I*J), 1+(kk-1)*(I*J):kk*(I*J)]
        Bk = add_diag(-Delta*Ak, 1 + Delta*(rrho + ddeath) - Delta*lambda[kk,kk])
        uk_stacked   = vec(u[:,:,kk])
        Vk_stacked   = vec(Vn[:,:,kk])
        Vkp_stacked  = sum(broadcast(*,  lambda[kk, indx_k]',
                                    reshape(Vn[:,:,indx_k], I*J, N-1)), dims=2)
        qk           = Delta .* uk_stacked + Vk_stacked + Delta .* Vkp_stacked
        Vn1k_stacked = Bk \ qk
        Vn_new[:,:,kk] = vec(Vn1k_stacked)
    end
    return Vn_new
end

@inline function set_grids(a, b, a_g, b_g, N)
    I   = length(b)
    J   = length(a)
    I_g = length(b_g)
    J_g = length(a_g)

    #b_grid   = permutedims(repeat(b, 1, J, N), [1 2 3])
    #a_grid   = permutedims(repeat(a, 1, I, N), [2 1 3])

    b_g_grid = permutedims(repeat(b_g, 1, J_g, N),  [1 2 3])
    a_g_grid = permutedims(repeat(a_g, 1, I_g, N),  [2 1 3])

    dbf_g_grid = Array{Float64}(undef, I_g, J_g, N)
    dbf_g_grid[1:I_g-1,:,:] = b_g_grid[2:I_g,:,:] - b_g_grid[1:I_g-1,:,:]
    dbf_g_grid[I_g,:,:] = dbf_g_grid[I_g-1,:,:]

    dbb_g_grid = Array{Float64}(undef, I_g,J_g,N)
    dbb_g_grid[2:I_g,:,:] = b_g_grid[2:I_g,:,:] - b_g_grid[1:I_g-1,:,:]
    dbb_g_grid[1,:,:] = dbb_g_grid[2,:,:]

    daf_g_grid = Array{Float64}(undef, I_g, J_g, N)
    daf_g_grid[:,1:J_g-1,:] = a_g_grid[:,2:J_g,:] - a_g_grid[:,1:J_g-1,:]
    daf_g_grid[:,J_g,:] = daf_g_grid[:,J_g-1,:]

    dab_g_grid = Array{Float64}(undef, I_g, J_g, N)
    dab_g_grid[:,2:J_g,:] = a_g_grid[:,2:J_g,:] - a_g_grid[:,1:J_g-1,:]
    dab_g_grid[:,1,:] = dab_g_grid[:,2,:]

    db_g_tilde       = 0.5*(dbb_g_grid[:,1,1] + dbf_g_grid[:,1,1])
    db_g_tilde[1]    = 0.5*dbf_g_grid[1,1,1]
    db_g_tilde[end]  = 0.5*dbb_g_grid[end,1,1]
    da_g_tilde       = 0.5*(dab_g_grid[1,:,1] + daf_g_grid[1,:,1])
    da_g_tilde[1]    = 0.5*daf_g_grid[1,1,1]
    da_g_tilde[end]  = 0.5*dab_g_grid[1,end,1]

    dab_g_tilde      = kron(da_g_tilde, db_g_tilde)
    dab_g_tilde_grid = reshape(repeat(dab_g_tilde, N, 1), I_g, J_g, N)

    return a_g_grid, b_g_grid, dab_g_tilde_grid, dab_g_tilde
end

"""
```
@inline function backward_difference(a_g::Vector{T}, b_g::Vector{T}) where {T<:Real}
```
Instantiates necessary difference vector.
"""
@inline function backward_difference(a_g::Vector{T}, b_g::Vector{T}) where {T<:Real}

    I_g, J_g = length(b_g), length(a_g)

    dbf_g = similar(b_g)
    dbf_g[1:I_g-1] = b_g[2:I_g] - b_g[1:I_g-1]
    dbf_g[I_g] = dbf_g[I_g-1]

    dbb_g = similar(b_g)
    dbb_g[2:I_g] = b_g[2:I_g] - b_g[1:I_g-1]
    dbb_g[1] = dbb_g[2]

    daf_g = similar(a_g)
    daf_g[1:J_g-1] = a_g[2:J_g] - a_g[1:J_g-1]
    daf_g[J_g] = daf_g[J_g-1]

    dab_g = similar(a_g)
    dab_g[2:J_g] = a_g[2:J_g] - a_g[1:J_g-1]
    dab_g[1] = dab_g[2]

    db_g_tilde       = 0.5*(dbb_g + dbf_g)
    db_g_tilde[1]    = 0.5*dbf_g[1]
    db_g_tilde[end]  = 0.5*dbb_g[end]

    da_g_tilde       = 0.5*(dab_g + daf_g)
    da_g_tilde[1]    = 0.5*daf_g[1]
    da_g_tilde[end]  = 0.5*dab_g[end]

    return kron(da_g_tilde, db_g_tilde)
end

"""
```
@inline function transition(ddeath::R, pam::R, ξ::R, w::S, a_lb::R, aggZ::S,
                            d::Array{U,3}, d_g::Array{U,3}, s::Array{U,3}, s_g::Array{U,3},
                            r_a::U, a::Array{R,1}, a_g::Array{R,1}, b::Array{R,1},
                            b_g::Array{R,1}, y::Vector{U}, cost, util, deposit;
                            permanent::Bool=false) where {R<:AbstractFloat, S<:Number, U<:Number}
```
Constructs matrix which captures transition dynamics for both HJB and KF equations.
"""
@inline function transition(ddeath::R, pam::R, ξ::R, w::S, a_lb::R, aggZ::S, d::Array{U,3},
                            d_g::Array{U,3}, s::Array{U,3}, s_g::Array{U,3}, r_a::U,
                            a::Array{R,1}, a_g::Array{R,1}, b::Array{R,1}, b_g::Array{R,1},
                            y::Vector{U}, cost; permanent::Bool=false) where {R<:AbstractFloat,
                                                                              S<:Number, U<:Number}
    I, J, N     = size(d)
    I_g, J_g, _ = size(d_g)

    perm_const = r_a + (permanent ? ddeath * pam - aggZ : ddeath * pam)
    ξw       = ξ * w

    χ, ζ = similar(d), similar(d)
    X, Z = similar(d), similar(d)

    # Compute drifts for HJB
    for i=1:I, j=1:J, n=1:N

        adrift = a[j] * perm_const + (ξw * y[n])
        bdrift = -d[i,j,n] - cost(d[i,j,n], a[j])

        χ[i,j,n] = (j==1) ? 0.0 : -(min(d[i,j,n], 0) + min(adrift, 0)) / (a[j] - a[j-1])
        ζ[i,j,n] = (j==J) ? 0.0 :  (max(d[i,j,n], 0) + max(adrift, 0)) / (a[j+1] - a[j])

        X[i,j,n] = (i==1) ? 0.0 : -(min(bdrift, 0) + min(s[i,j,n], 0)) / (b[i] - b[i-1])
        Z[i,j,n] = (i==I) ? 0.0 :  (max(bdrift, 0) + max(s[i,j,n], 0)) / (b[i+1] - b[i])
    end

    yyy = -vec(X .+ Z .+ χ .+ ζ)

    χ = reshape(χ, I*J, N)
    χ = circshift(χ, -I)
    ζ = reshape(ζ, I*J, N)
    ζ = circshift(ζ, I)

    X  = reshape(X, I*J, N)
    X  = circshift(X, -1)
    Z  = reshape(Z, I*J, N)
    Z  = circshift(Z, 1)

    A = spdiagm(0 => yyy, I => vec(ζ)[I+1:end], -I => vec(χ)[1:end-I],
                          1 => vec(Z)[2:end],   -1 => vec(X)[1:end-1])

    # TODO: is it possible that below (ddeathpam) is incorrect?
    χu, ζu = similar(d_g), similar(d_g)
    Xu, Zu = similar(d_g), similar(d_g)

    # Compute drifts for KFE
    for i=1:I_g, j=1:J_g, n=1:N
        audrift = d_g[i,j,n] + a_g[j] * (r_a + ddeath*pam) + ξw * y[n] -
            (permanent ? aggZ * a_g[j] : 0.0)

        budrift = s_g[i,j,n] - d_g[i,j,n] - cost(d_g[i,j,n], a_g[j]) -
            (permanent ? aggZ * b_g[i] : 0.0)

        χu[i,j,n] = (j == 1)   ? 0.0 : -min(audrift, 0) / (a_g[j] - a_g[j-1])
        ζu[i,j,n] = (j == J_g) ? 0.0 :  max(audrift, 0) / (a_g[j+1] - a_g[j])

        Xu[i,j,n] = (i == 1)   ? 0.0 : -min(budrift, 0) / (b_g[i] - b_g[i-1])
        Zu[i,j,n] = (i == I_g) ? 0.0 :  max(budrift, 0) / (b_g[i+1] - b_g[i])
    end
    yyyu  = -vec(χu .+ ζu .+ Xu .+ Zu)

    χu = reshape(χu,I_g*J_g,N)
    χu = circshift(χu,-I_g)
    ζu = reshape(ζu,I_g*J_g,N)
    ζu = circshift(ζu,I_g)

    Xu = reshape(Xu, I_g*J_g, N)
    Xu = circshift(Xu, -1)
    Zu = reshape(Zu, I_g*J_g, N)
    Zu = circshift(Zu, 1)

    AT = spdiagm(0 => yyyu, I_g => vec(ζu)[I_g+1:end], -I_g => vec(χu)[1:end-I_g],
                              1 => vec(Zu)[2:end],       -1 => vec(Xu)[1:end-1])
    return A, AT
end

@inline function transition_deriva_lite(permanent::Bool, ddeath::R, pam::R, ξ::R, w::S,
                                        a_lb::R, aggZ::S, d::Array{U,2}, d_g::Array{U,2},
                                        s::Array{U,2}, s_g::Array{U,2}, r_a::U, a::Array{R,1},
                                        a_g::Array{R,1}, b::Array{R,1}, b_g::Array{R,1},
                                        y::U, lambda_ii::R, cost) where {R<:AbstractFloat,
                                                                         S<:Number, U<:Number}
    I, J = size(d)
    I_g, J_g = size(d_g)

    perm_const = permanent ? ddeath * pam - aggZ : ddeath * pam

    χ, ζ = similar(d), similar(d)
    X, Z = similar(d), similar(d)

    for i=1:I, j=1:J
        adrift = a[j] * (r_a + perm_const) + (ξ * w * real(y))

        χ[i,j] = (j==1) ? 0.0 : -(min(d[i,j], 0) + min(adrift, 0)) / (a[j] - a[j-1])
        ζ[i,j] = (j==J) ? 0.0 :  (max(d[i,j], 0) + max(adrift, 0)) / (a[j+1] - a[j])

        X[i,j] = (i==1) ? 0.0 : -(min(-d[i,j] - cost(d[i,j], a[j]), 0) +
                                  min(s[i,j], 0)) / (b[i] - b[i-1])
        Z[i,j] = (i==I) ? 0.0 :  (max(-d[i,j] - cost(d[i,j], a[j]), 0) +
                                  max(s[i,j], 0)) / (b[i+1] - b[i])
    end
    yyy = -vec(X .+ Z .+ χ .+ ζ)

    χ = reshape(χ, I*J)
    χ = circshift(χ, -I)
    ζ = reshape(ζ, I*J)
    ζ = circshift(ζ, I)

    X  = reshape(X, I*J)
    X  = circshift(X, -1)
    Z  = reshape(Z, I*J)
    Z  = circshift(Z, 1)

    A = spdiagm(0 => yyy, I => vec(ζ)[I+1:end], -I => vec(χ)[1:end-I],
                          1 => vec(Z)[2:end],   -1 => vec(X)[1:end-1])

    χu, ζu = similar(d_g), similar(d_g)
    Xu, Zu = similar(d_g), similar(d_g)

    for i=1:I_g, j=1:J_g
        # Compute drifts for KFE -- is it possible that below (ddeathpam) is incorrect?
        audrift = d_g[i,j] + a_g[j] * (r_a + ddeath*pam) + ξ * w * y -
            (permanent ? aggZ * a_g[j] : 0.0)

        budrift = s_g[i,j] - d_g[i,j] - cost(d_g[i,j], max(a_g[j], a_lb)) -
            (permanent ? aggZ * b_g[i] : 0.0)

        χu[i,j] = (j==1)   ? 0.0 : -min(audrift, 0) / (a_g[j] - a_g[j-1])
        ζu[i,j] = (j==J_g) ? 0.0 :  max(audrift, 0) / (a_g[j+1] - a_g[j])

        Xu[i,j] = (i==1)   ? 0.0 : -min(budrift, 0) / (b_g[i] - b_g[i-1])
        Zu[i,j] = (i==I_g) ? 0.0 :  max(budrift, 0) / (b_g[i+1] - b_g[i])
    end
    yyyu  = -vec(χu .+ ζu .+ Xu .+ Zu) .+ lambda_ii

    χu = reshape(χu, I_g*J_g)
    χu = circshift(χu, -I_g)
    ζu = reshape(ζu, I_g*J_g)
    ζu = circshift(ζu, I_g)

    Xu = reshape(Xu, I_g*J_g)
    Xu = circshift(Xu, -1)
    Zu = reshape(Zu, I_g*J_g)
    Zu = circshift(Zu, 1)

    # AT is pre-transposed for efficiency
    AT = spdiagm(0 => yyyu, -I_g => vec(ζu)[I_g+1:end], I_g => vec(χu)[1:end-I_g],
                            -1 => vec(Zu)[2:end],          1 => vec(Xu)[1:end-1])
    return A, AT
end

"""
```
@inline function construct_problem_functions(γ::T, χ0::R, χ1::R, χ2::R,
                                             a_lb::R) where {T<:Number,R<:Number}
```
Initializes utility, adjustment cost, and deposit functions.
"""
@inline function construct_problem_functions(γ::T, χ0::R, χ1::R, χ2::R,
                                             a_lb::R) where {T<:Number,R<:Number}

    util = if γ == 1.0
        @inline util1(c::S) where {S<:Number} = log(c)
    else
        @inline util2(c::S) where {S<:Number} = 1.0 / (1.0 - γ) * (c ^ (1-γ) - 1.0)
    end

    @inline function deposit(Va, Vb, a)
        indx_plus  = ((Va / Vb - 1 - χ0) > 0)
        indx_minus = ((Va / Vb - 1 + χ0) < 0)
        return χ1 * (max(Va / Vb - 1 - χ0, 0)) ^ (1/χ2) * a * indx_plus +
            (-χ1) * (max(-(Va / Vb - 1) - χ0, 0)) ^ (1/χ2) * a * indx_minus
    end

    @inline function cost(d, a)
        d_scaled = abs(d / max(a, a_lb))
        return max(a, a_lb) * (χ0 * (d_scaled) + 1.0 / (1.0 + χ2) *
                               ((d_scaled)^(1+χ2) * χ1^(-χ2)))
    end
    return util, deposit, cost
end

"""
```
@inline function create_y_grid(y_size::Int64, ygrid_new::Int64)
```
Function initializes income grid.
"""
@inline function create_y_grid(y_size::Int64, ygrid_new::Int64)#, dataroot::String)

    y      = Array{Float64}(undef, 0)
    y_dist = Array{Float64}(undef, 0)
    λ      = Array{Float64}(undef, 0)

    if y_size == 2
        y  = [0.1, 1]

        # Transition probabilities
        λ1 = 0.5 # expected duration of unemployment is 2 quarters
        λ2 = (λ1 / (y[2] * .93 - y[1])) * (y[2] - y[2] * .93) # unemployment rate 7
        λ  = [-λ1 λ1; λ2 -λ2]

    elseif y_size == 30
        # Load grid and transition probabilities from external calibration
        #y      = load(dataroot * "income_grid_30.jld2", "y")
        #y_dist = load(dataroot * "income_grid_30.jld2", "y_dist")
        #λ      = load(dataroot * "income_transition_30.jld2", "lambda")
        y = [-4.0043487, -2.3712946, -2.2529091, -1.861575, -1.7594332, -1.751468,
             -1.7514113, -1.7434461, -1.6413042, -1.1315846, -0.61985494, -0.50146947,
             -0.1101354, -0.0079935576, -2.8361989e-5, 2.8362067e-5, 0.0079935649,
             0.11013543, 0.50146947, 0.619855, 1.1315847, 1.6413042, 1.7434461, 1.7514113,
             1.751468, 1.7594332, 1.8615751, 2.2529091, 2.3712946, 4.0043487]

        y_dist = [0.001696604, 0.0036905108, 0.054850618, 0.0034315242, 0.002345745,
                  0.0034002893, 0.0034002799, 0.0023457452, 0.0034315242, 0.0036905106,
                  0.11931293, 0.001696604, 0.11093998, 0.07583712, 0.10993017, 0.10992986,
                  0.075837123, 0.11093998, 0.001696604, 0.11931292, 0.0036905108, 0.0034315242,
                  0.002345745, 0.0034002893, 0.0034002799, 0.0023457452, 0.0034315242,
                  0.054850615, 0.0036905106, 0.001696604]

        λ = build_lambda(y_size)

        # Create grid of income shocks
        y = exp.(y)

    elseif y_size == 33
        # Load grid and transition probabilities from external calibration
        #y      = load(dataroot * "income_grid.jld2", "y")
        #y_dist = load(dataroot * "income_grid.jld2", "y_dist")
        #λ      = load(dataroot * "income_transition.jld2", "lambda")
        y = [-4.7415191, -3.9041665, -3.1769113, -2.8490016, -2.5749613, -2.1240338,
             -2.011649, -1.8925175, -1.6610012, -1.2843938, -1.2100737, -0.95648415,
             -0.68244382, -0.60812364, -0.23151632, -0.1191315, 0.0, 0.1191315, 0.23151632,
             0.60812364, 0.68244382, 0.95648415, 1.2100737, 1.2843938, 1.6610012, 1.8925175,
             2.011649, 2.1240338, 2.5749613, 2.8490016, 3.1769113, 3.9041665, 4.7415191]

        y_dist = [0.00026325521, 0.00073171922, 0.0015341331, 0.008947109, 0.0025562662,
                  0.0032773437, 0.024868536, 0.011062805, 0.0032773437, 0.052139733,
                  0.0025562662, 0.00026325521, 0.086878404, 0.0015341331, 0.11138526,
                  0.00073171922, 0.37598543, 0.00073171922, 0.11138526, 0.0015341331,
                  0.086878404, 0.00026325521, 0.0025562662, 0.052139733, 0.0032773437,
                  0.011062805, 0.024868536, 0.0032773437, 0.0025562662, 0.008947109,
                  0.0015341331, 0.00073171922, 0.00026325521]

        λ = build_lambda(y_size)
        # Create grid of income shocks
        y = exp.(y)
    end

    # Compute stationary income distribution
    y_dist = stat_dist(λ')   # stationary distribution
    y_mean = dot(y, y_dist)

    if     ygrid_new == 0
        y = y ./ y_mean
    elseif ygrid_new == 1
        y = 0.1907858 .* y ./ y_mean
    end

    y_mean = dot(y, y_dist)
    N = size(y, 2)
    return λ, y, y_mean, y_dist, N
end

"""
```
@inline function create_a_grid(agrid_new::Int64, J::Int64, J_g::Int64, amin::Int64, amax::Int64)
```
Function initializes illiquid asset grids.
"""
@inline function create_a_grid(agrid_new::Int64, J::Int64, J_g::Int64, amin::Int64, amax::Int64)
    a   = Array{Float64}(undef, J)
    a_g = Array{Float64}(undef, J_g)

    if agrid_new == 0

        a 		    = range(0, stop = 1, length = J) # set the grid
        coeff_power = 0.9
        power       = 8
        a 			= amax*((1 - coeff_power) * a + coeff_power * (a .^ power))

        a_g         = range(0, stop = 1, length = J_g)
        amax_g      = 150
        coeff_power = 0.9
        power       = 8
        a_g         = amax_g*((1-coeff_power) * a_g + coeff_power * (a_g .^ power))

    elseif agrid_new == 1

        agridparam = 0.15
        a_raw = range(0, stop = 1, length = J)
        a_raw = a_raw.^(1/agridparam)

        a = amin .+ (amax - amin) * a_raw

        for i = 1:9
            a[i] = (i-1) * a[10]/(10 - 1)
        end

        # Define new grid points for g here. For now, it is created the same.
        agridparam = 0.15

        a_raw = range(0, stop = 1, length = J_g)
        a_raw = a_raw.^(1/agridparam)
        a_g = amin .+ (amax - amin) * a_raw

        for i = 1:9
            a_g[i] = (i-1) * a_g[10]/(10 - 1)
        end

    elseif agrid_new == 2

        agridparam = 0.15

        a_raw = range(0, stop = 1, length = J)
        a_raw = a_raw.^(1/agridparam)
        a = amin .+ (amax - amin) * a_raw

        for i = 1:9
            a[i] = (i-1) * a[10]/(10 - 1)
        end

        # Define new grid points for g here. For now, it is created the same.
        agridparam = 0.21

        a_raw = range(0, stop = 1, length = J_g)
        a_raw = a_raw.^(1/agridparam)
        a_g = amin .+ (amax - amin) * a_raw

        for i = 1:9
            a_g[i] = (i-1) * a_g[10]/(10. - 1.)
        end

    end

    a_g_0pos = findall(!iszero, a_g .== 0)

    return a, a_g, a_g_0pos
end

"""
```
@inline function create_b_grid(bgrid_new::Int64, I::Int64)
```
Function to initialize liquid asset grid.
"""
@inline function create_b_grid(bgrid_new::Int64, I::Int64, I_g::Int64)
    I = 50 # TEMP
    I_neg = 10
    I_pos = 40

    bmin = Array{Float64}(undef, 0)
    bmax = Array{Float64}(undef, 0)
    b    = Array{Float64}(undef, 0)
    b_g  = Array{Float64}(undef, 0)

    if bgrid_new == 0

        # Lower bound
        b_    = -2		# lower bound on liquid assets

        # Preparations
        I_neg = Int(I / 10) # number of grid points on the negative part of the real line
        bmin  = b_ 		# lower bound
        bmax  = 80

        # Grid for positive part
        b_pos 			= range(0, stop = 1, length = Int(I-I_neg))
        coeff_power 	= 0.9
        power = 8
        b_pos 			= bmax*((1 - coeff_power) * b_pos + coeff_power * (b_pos.^power))

        # Grid for negative part
        b_neg 			= range(0, stop = 1, length = Int(I_neg+1))
        coeff_power 	= 0
        power           = 1
        b_neg 			= -bmin*((1 - coeff_power) * b_neg + coeff_power * (b_neg.^power)) .+ bmin
        b_neg           = b_neg[1:end-1]

        # Put together grid on liquid assets
        b = [b_neg;b_pos]

        ## Define new grid points for g here. For now, it is created the same.
        I_g_neg = 5
        bmax_g  = 70

        # Grid for positive part
        b_pos 			= range(0, stop = 1, length = Int(I_g-I_g_neg))
        coeff_power 	= 0.9
        power = 8
        b_pos 			= bmax_g*((1 - coeff_power) * b_pos + coeff_power * (b_pos.^power))

        # Grid for negative part
        b_neg 			= range(0, stop = 1, length = Int(I_g_neg+1))
        coeff_power 	= 0
        power = 1
        b_neg = -bmin*((1 - coeff_power) * b_neg + coeff_power * (b_neg.^power)) .+ bmin
        b_neg = b_neg[1:end-1]

        # Put together grid on liquid assets
        b_g = [b_neg;b_pos]

    elseif bgrid_new == 1

        # set parameters
        I_neg = I - I_pos
        bmin  = -1
        b0    = 0
        bmax  = 40
        bgridparam = 0.35
        bgridparam_neg = 0.4
        b_ = bmin

        # positive part
        bp_raw = range(0, stop = 1, length = I_pos)
        bp_raw = bp_raw .^ (1 / bgridparam)
        bp     = b0 .+ (bmax - b0) * bp_raw

        # negative part
        bn_raw = range(0, stop = 1, length = Int64(I_neg/2+1))
        bn_raw = bn_raw.^(1/bgridparam_neg)
        bn = bmin .+ (bmin/2 - bmin) * bn_raw
        bn = [bn; Vector{Float64}(undef, Int64(I_neg/2-1))]
        for i = Int64(I_neg/2+2) : I_neg
            bn[i] = b0 - (bn[Int64(I_neg)+2-i] - bn[1])
        end

        # Put everything together
        b = [bn; bp]

        ## Define new grid points for g here. For now, it is created the same.
        # set parameters
        I_neg = I_g - I_pos
        bmin  = -1
        b0    = 0
        bmax  = 40
        bgridparam = 0.35
        bgridparam_neg = 0.4
        b_ = bmin

        # positive part
        bp_raw = range(0, stop = 1, length = I_pos)
        bp_raw = bp_raw.^(1/bgridparam)
        bp = b0 .+ (bmax - b0) * bp_raw

        # negative part
        bn_raw = range(0, stop = 1, length = Int64(I_neg/2+1))
        bn_raw = bn_raw.^(1/bgridparam_neg)
        bn = bmin .+ (bmin/2 - bmin) * bn_raw
        bn = [bn;Vector{Float64}(undef,Int64(I_neg/2-1))]
        for i = Int64(I_neg/2+2):I_neg
            bn[i] = b0 - (bn[Int64(I_neg+2-i)] - bn[1])
        end

        # Put everything together
        b_g = [bn; bp]

    elseif bgrid_new == 2

        # set parameters
        bmin = -1
        b0   = 0
        bmax = 40
        bgridparam = 0.25
        bgridparam_neg = 0.4
        b_ = bmin

        # positive part
        bp_raw = range(0, stop = 1, length = I_pos)
        bp_raw = bp_raw.^(1/bgridparam)
        bp = b0 .+ (bmax - b0) * bp_raw

        # negative part
        bn_raw = range(0, stop = 1, length = Int(I_neg/2+1))
        bn_raw = bn_raw.^(1/bgridparam_neg)
        bn = bmin .+ (bmin/2 - bmin) * bn_raw
        bn = [bn;Vector{Float64}(undef,Int(I_neg/2-1))]

        for i = I_neg/2+2:I_neg
            bn[i] = b0 - (bn[I_neg+2-i] - bn[1])
        end

        # put everything together
        b = [bn;bp]

        ## Define new grid points for g here. For now, it is created the same.
        # set parameters
        I_pos = 40
        I_neg = I_g - I_pos
        bmin = -1
        b0 = 0
        bmax = 40
        bgridparam = 0.6
        bgridparam_neg = 0.4
        b_ = bmin

        # positive part
        bp_raw = range(0, stop = 1, length = I_pos)
        bp_raw = bp_raw.^(1 / bgridparam)
        bp = b0 .+ (bmax - b0) * bp_raw

        # negative part
        bn_raw = range(0, stop = 1, length = Int(I_neg/2+1))
        bn_raw = bn_raw.^(1/bgridparam_neg)
        bn = bmin .+ (bmin/2 - bmin) * bn_raw
        bn = [bn;Vector{Float64}(undef,Int(I_neg/2-1))]
        for i = Int(I_neg/2+2):I_neg
            bn[i] = b0 - (bn[Int(I_neg+2-i)] - bn[1])
        end

        # put everything together
        b_g = [bn;bp]
    end
    b_g_0pos = findall(b_g .== 0)

    return I_neg, bmin, bmax, b, b_g, b_g_0pos
end


text:
doc"""
  SigmoidNumbers.find_lsb(p)

  returns the exponent of the least significant bit of the posit.
"""
function find_lsb{N,ES}(p::Posit{N,ES})
  (sgn, exp, frc) = posit_components(p)

  exp - (64 - trailing_zeros(frc)) - (sgn && frc != zero(UInt64))

end

ext: jl

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
export
   expected_runtime

#measure decorating functions
function measure(qwe::QWSearch, state::QSearchState)
   measure(qwe, state.state)
end

function measure(qwe::QWSearch, state::QSearchState, vertices::Vector{Int})
   measure(qwe, state.state, vertices)
end

#Evolve decorating function
function evolve(qwe::QWSearch, state::QSearchState)
   evolve(qwe, state.state)
end

function evolve(qwe::QWSearch, state::QSearchState, runtime::T) where T<:Real
   evolve(qwe, state.state, runtime)
end

#Expected runtime
"""
    expected_runtime(runtime, probability)
    expected_runtime(state)

Returns the expected runtime needed for quantum walk, considering it as Bernoulli
process. It equals to `runtime/probability`. In the case of `state` provided the
measurement is made, penalty is included.
"""
function expected_runtime(runtime::Real, probability::Real)
   runtime/probability
end

function expected_runtime(qws::QWSearch, state::QSearchState)
   expected_runtime(state.runtime + state.penalty, sum(state.probability))
end


text:
using Distributions
using DelimitedFiles
using Printf
using Dates

# start timer
start_time = Dates.now()

# set initial parameters' values
population_size = 50
dimension_size = 3
donors_number = 1
receivers_number = 1
maximum_evaluations = 5000
bound = 1000

# other dependent parameters, no need to change
current_evaluation = population_size
x_upper_bound = bound
x_lower_bound = -bound
y_upper_bound = bound
y_lower_bound = -bound
z_upper_bound = bound
z_lower_bound = 0

# BUILDING = [20, 50, 200]     #b1
# BUILDING = [20, 50, 250]     #b2
# BUILDING = [20, 50, 300]     #b3
# BUILDING = [10, 50, 250]     #b4
# BUILDING = [30, 50, 250]     #b5
BUILDING = [50, 50, 250]     #b6

# users location data file path
user_data = string("users/UserLocations_", BUILDING[1], "_", BUILDING[2], "_", BUILDING[3], ".dat")
Users_Locations = readdlm(user_data)

# enable/disable store output into file mode
store_output = true
store_file_path = string("tests/output_", BUILDING[1], "_",BUILDING[2], "_", BUILDING[3], ".dat")

# ========================================================= #
#                       Functions                           #
# ========================================================= #

# generating the initial population
function generate_drones()
    population = zeros(population_size, dimension_size)
    i = 1
    while i <= population_size
        new_drone = generate_drone()
        if !drone_is_inside(new_drone)
            population[i, :] = new_drone
            i += 1
        end
    end
    return population
end

# generating just one individual
function generate_drone()
    drone = zeros(1, dimension_size)
    drone[1, 1] = x_lower_bound + rand() * (x_upper_bound - x_lower_bound)
    drone[1, 2] = y_lower_bound + rand() * (y_upper_bound - y_lower_bound)
    drone[1, 3] = z_lower_bound + rand() * (z_upper_bound - z_lower_bound)
    return drone
end

# calculating fitness of all individuals in population
function calculate_fitnesses(drones)
    fitnesses = zeros(population_size)
    for i = 1:population_size
        fitnesses[i] = fitness(drones[i, :])
    end
    return fitnesses
end

# calculation fitness for one individual
function fitness(drone)
    total_sum = 0
    for index = 1:(size(Users_Locations)[1])
        total_sum += fitness_per_user(Users_Locations[index, :], drone)
    end
    return total_sum
end

# calculating signal strength for one user in the building
function fitness_per_user(user, drone)
    dotProduct = 1
    up = 0
    dp = 0
    d2D = sqrt(drone[1] ^ 2)
    d3D = 0

    for i = 1:dimension_size
        d3D += (user[i] - drone[i]) ^ 2
        dotProduct += user[i] * drone[i]
        up += user[i] ^ 2
        dp += drone[i] ^ 2
    end

    d3D = sqrt(d3D)
    mag_mul = sqrt(up) * sqrt(dp)

    result = (20.0 * log10(d3D) + 20.0 * ( 0.301 ) + 32.4) + (14.0 + 15.0 * ((1.0 - dotProduct/mag_mul) ^ 2) ) + (0.5 * d2D)
    return result
end

# perform infection between two individuals
function perform_infection(x_k, x_m)
    j = rand(1:dimension_size)
    x = copy(x_k)
    x[j] = x_k[j] + (rand(Uniform(-1.0, 1.0)) * (x_k[j] - x_m[j]))
    return check_bounds(x)
end

# check if exceeded bounds
function check_bounds(drone)
    # check drone's x location
    if drone[1] > x_upper_bound
        drone[1] = x_upper_bound
    elseif drone[1] < x_lower_bound
        drone[1] = x_lower_bound
    end
    # check drone's y location
    if drone[2] > y_upper_bound
        drone[2] = y_upper_bound
    elseif drone[2] < y_lower_bound
        drone[2] = y_lower_bound
    end
    # check drone's z location
    if drone[3] > z_upper_bound
        drone[3] = z_upper_bound
    elseif drone[3] < z_lower_bound
        drone[3] = z_lower_bound
    end
    while drone_is_inside(drone)
        drone = generate_drone()
    end
    return drone
end

# get lists of indexes of doreceivers_numbers and recievers
function get_donors_and_receivers_indexes(fitnesses)
    donors = zeros(Int64, donors_number)
    receivers = zeros(Int64, receivers_number)
    sorted_indexes = sortperm(fitnesses)
    for i = 1:donors_number
        donors[i] = sorted_indexes[i]
    end
    reverse!(sorted_indexes)
    for i = 1:receivers_number
        receivers[i] = sorted_indexes[i]
    end
    return donors, receivers
end

# performing plasma tranfer from donor to receiver indvidual
function perform_plasma_transfer(receiver, donor)
    for j = 1:dimension_size
        receiver[j] += rand(Uniform(-1, 1)) * (receiver[j] - donor[j])
    end
    return check_bounds(receiver)
end

# updating donor's parameters
function update_donor(donor)
    for j = 1:dimension_size
        donor[j] += rand(Uniform(-1, 1)) * donor[j]
    end
    return check_bounds(donor)
end

# compare individual's fitness with global fitness value
function compare_with_best_fitness(x)
    if fitness(x) < best_fitness
        global best_fitness = fitness(x)
        best_index = argmin(fitnesses)
        best = population[best_index, :]
        if !store_output
            println("Best: $best_fitness \t - location: $best")
        end
    end
end

# check if the drone is inside the building or not
function drone_is_inside(drone)
    x_in = false
    y_in = false
    z_in = false
    if 0 <= drone[1] <= BUILDING[1]
        x_in = true
    end
    if 0 <= drone[2] <= BUILDING[2]
        y_in = true
    end
    if 0 <= drone[3] <= BUILDING[3]
        z_in = true
    end
    return (x_in && y_in && z_in)
end

# ========================================================= #
#                      Start of IPA                         #
# ========================================================= #

# generating initial population
population = generate_drones()

# calculating fitness of population
fitnesses = calculate_fitnesses(population)

# finding best individual fitness
best_fitness = minimum(fitnesses)

# print initial best fitness value and other statistics
if !store_output
    println("Initial best fitness value: $best_fitness")
    println("Number of parameters: $dimension_size")
    println("Population size: $population_size")
end

while current_evaluation < maximum_evaluations

    # start of infection phase
    for index = 1:population_size
        if current_evaluation < maximum_evaluations
            global current_evaluation += 1
            random_index = rand(1:population_size)
            while random_index == index
                random_index = rand(1:population_size)
            end
            current_individual = population[index, :]
            random_individual = population[random_index, :]
            infected_individual = perform_infection(current_individual, random_individual)
            fitness_of_infected = fitness(infected_individual)
            if fitness_of_infected < fitnesses[index]
                population[index, :] = copy(infected_individual)
                fitnesses[index] = fitness_of_infected
                compare_with_best_fitness(infected_individual)
            end
        else
            break # if exceeded maximum evaluation number
        end
    end

    # start of plasma transfering phase
    # generating dose_control and treatment_control vectors
    dose_control = ones(Int64, receivers_number)
    treatment_control = ones(Int64, receivers_number)

    # get indexes of both of donors and receivers
    donors_indexes, receivers_indexes = get_donors_and_receivers_indexes(fitnesses)

    for i = 1:receivers_number
        receiver_index = receivers_indexes[i]
        random_donor_index = donors_indexes[Int(rand(1:donors_number))]
        current_receiver = population[receiver_index, :]
        random_donor = population[random_donor_index, :]
        while treatment_control[i] == 1
            if current_evaluation < maximum_evaluations
                global current_evaluation += 1
                treated_individual = perform_plasma_transfer(current_receiver, random_donor)
                treated_fitness = fitness(treated_individual)
                if dose_control[i] == 1
                    if treated_fitness < fitnesses[random_donor_index]
                        dose_control[i] += 1
                        population[receiver_index, :] = copy(treated_individual)
                        fitnesses[receiver_index] = treated_fitness
                    else
                        population[receiver_index, :] = copy(random_donor)
                        fitnesses[receiver_index] = fitnesses[random_donor_index]
                        treatment_control[i] = 0
                    end
                else
                    if treated_fitness < fitnesses[receiver_index]
                        population[receiver_index, :] = copy(treated_individual)
                        fitnesses[receiver_index] = treated_fitness
                    else
                        treatment_control[i] = 0
                    end
                end
                compare_with_best_fitness(population[receiver_index, :])
            else
                break # if exceeded maximum evaluation number
            end
        end
    end

    # start of donors updating phase
    for i = 1:donors_number
        if current_evaluation < maximum_evaluations
            global current_evaluation += 1
            donor_index = donors_indexes[i]
            if (current_evaluation / maximum_evaluations) > rand()
                population[donor_index, :] = update_donor(population[donor_index, :])
            else
                population[donor_index, :] = generate_drone()
            end
            fitnesses[donor_index] = fitness(population[donor_index, :])
            compare_with_best_fitness(population[donor_index, :])
        else
            break # if exceeded maximum evaluation number
        end
    end
end

if !store_output
    # print elapsed time
    end_time = Dates.now()
    total_time = (Dates.value(end_time) - Dates.value(start_time)) / 1000
    @printf("Elapsed time: %.2f seconds\n", total_time)

    # print best fitness value in scientific notation
    @printf("Best fitness value: %.6e\n", best_fitness)
end

if store_output
    best_index = argmin(fitnesses)
    best = population[best_index, :]
    x = best[1]
    y = best[2]
    z = best[3]
    output_str = "$best_fitness $x $y $z\n"
    open(store_file_path, "a+") do file
        write(file, output_str)
    end
end


text:
# Load X and y variable
using JLD
data = load("logisticData.jld")
(X,y) = (data["X"],data["y"])

# Add bias variable, initialize w, set regularization and optimization parameters
(n,d) = size(X)
X = [ones(n,1) X]
d += 1
w = zeros(d,1)
lambda = 1
maxPasses = 500
progTol = 1e-4
verbose = true

## Run and time coordinate descent to minimize L2-regularization logistic loss

# Start timer
tic()

# Compute Lipschitz constant of 'f'
(eigenvalues,~) = eig(X'*X)
L = .25*maximum(eigenvalues) + lambda;

# Start running coordinate descent
w_old = copy(w);
for k in 1:maxPasses*d

    # Choose variable to update 'j'
    j = rand(1:d)

    # Compute partial derivative 'g_j'
    yXw = y.*(X*w);
    sigmoid = 1./(1+exp.(-yXw));
    g = -X'*(y.*(1-sigmoid)) + lambda*w;
    g_j = g[j];

    # Update variable
    w[j] -= (1/L)*g_j;

    # Check for lack of progress after each "pass"
    if mod(k,d) == 0
        yXw = y.*(X*w)
        f = sum(log.(1 + exp.(-y.*(X*w)))) + (lambda/2)norm(w)^2
        delta = norm(w-w_old,Inf);
        if verbose
            @printf("Passes = %d, function = %.4e, change = %.4f\n",k/d,f,delta);
        end
        if delta < progTol
            @printf("Parameters changed by less than progTol on pass\n");
            break;
        end
        w_old = copy(w);
    end

end

# End timer
toc()


text:
module CUSPARSE

using ..APIUtils

using ..CUDA
using ..CUDA: CUstream, cuComplex, cuDoubleComplex, libraryPropertyType, cudaDataType
using ..CUDA: libcusparse, unsafe_free!, @retry_reclaim

using CEnum

using Memoize

using LinearAlgebra
using LinearAlgebra: HermOrSym

using Adapt

using DataStructures

using SparseArrays

const SparseChar = Char

# core library
include("libcusparse_common.jl")
include("error.jl")
include("libcusparse.jl")
include("libcusparse_deprecated.jl")

include("array.jl")
include("util.jl")
include("types.jl")

# low-level wrappers
include("helpers.jl")
include("management.jl")
include("level1.jl")
include("level2.jl")
include("level3.jl")
include("preconditioners.jl")
include("conversions.jl")
include("generic.jl")

# high-level integrations
include("interfaces.jl")

# cache for created, but unused handles
const handle_cache_lock = ReentrantLock()
const idle_handles = DefaultDict{CuContext,Vector{cusparseHandle_t}}(()->cusparseHandle_t[])

function handle()
    state = CUDA.active_state()
    handle, stream = get!(task_local_storage(), (:CUSPARSE, state.context)) do
        new_handle = @lock handle_cache_lock begin
            if isempty(idle_handles[state.context])
                cusparseCreate()
            else
                pop!(idle_handles[state.context])
            end
        end

        finalizer(current_task()) do task
            @spinlock handle_cache_lock begin
                push!(idle_handles[state.context], new_handle)
            end
        end
        # TODO: cusparseDestroy to preserve memory, or at exit?

        cusparseSetStream(new_handle, state.stream)

        new_handle, state.stream
    end::Tuple{cusparseHandle_t,CuStream}

    if stream != state.stream
        cusparseSetStream(handle, state.stream)
        task_local_storage((:CUSPARSE, state.context), (handle, state.stream))
    end

    return handle
end

end


text:
# The linearized forward modeling operator `F` and its adjoint (aka Migration operator)
# are the building blocks of iterative optimization schemes.
# For the `PoissonExpt`, we have the
# functionality to compute `F*x` without explicitly storing the operator matrix
# (see `LinearMaps.jl`).
# The perturbed field `δψ` due to a perturbation in `σ` is given by
# ```math
# δψ=-A⁻¹(σ₀)A(δσ)ψ₀,
# ```
# where
# ```math
# ∇⋅(σ₀(x,z)∇) ψ₀(t)=A(σ₀)ψ₀(t)=f(t).
# ```
# Lets start a tutorial.

# ### Load some packages
using GeoPhyInv
using SparseArrays
using StatsBase
using LinearAlgebra
using Random
using ProgressMeter
using LinearAlgebra
using Test
using ForwardDiff
using Calculus
using LinearMaps

# ### Setting up the spatial and temporal grids

nx=21
nz=21
nt=4
nznx=nz*nx
mgrid=[range(-div(nz,2), step=1.0, length=nz), range(-div(nx,2), step=1.0, length=nx)]
tgrid=range(0.0,step=0.5, length=nt)
@info "Grids are all set."

# ### Allocating medium parameters

Qv=abs.(randn(nz,nx))
η=abs.(randn(nz,nx))
k=abs.(randn(nz,nx))
σ=abs.(randn(nz,nx))
p=randn(nz,nx,nt)

σobs=abs.(randn(nz,nx))
Qobs=abs.(randn(nz,nx))
nr=10 # number of abstract receivers
ACQ=sprandn(nr,nz*nx,0.6); # choose a random acquisition operator
@info "We are ready for the PoissonExpt."

# ### Create an `Expt`, and then extract a linear forward map out of it
paE=PoissonExpt(p, tgrid, mgrid, Qv, k, η, σ, ACQ, σobs=σobs, Qobs=Qobs,)
F=operator_Born(paE, σ); # extract the linearized forward operator from `Expt`
GIPh.Utils.test_linearmap(F) # finally do some tests on the linearmap

# ### Usage

δx=randn(size(F,2)) # random model pertubation
δd=F*δx # corresponding pertubation in data
@info string("Length of data: (nt*nr)=",length(δd))


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(MathOptInterface.optimize!),Optimizer{Incremental{GLPK.Optimizer}, Incremental{Ipopt.Optimizer}, DefaultExt}})   # time: 11.415445
    Base.precompile(Tuple{typeof(reform_epigraph_min!),GlobalOptimizer{Incremental{GLPK.Optimizer}, Incremental{Ipopt.Optimizer}, DefaultExt},ParsedProblem,BufferedNonlinearFunction{1, NS}})   # time: 1.4057246
    Base.precompile(Tuple{typeof(forward_pass!),Evaluator,BufferedNonlinearFunction{2, NS}})   # time: 0.881403
    Base.precompile(Tuple{Type{Optimizer}})   # time: 0.5019401
    Base.precompile(Tuple{typeof(MathOptInterface.empty!),Optimizer{Incremental{GLPK.Optimizer}, Incremental{Ipopt.Optimizer}, DefaultExt}})   # time: 0.0625456
    Base.precompile(Tuple{typeof(initialize!),RelaxCache{1, NS},DirectedTree})   # time: 0.0522575
    Base.precompile(Tuple{Type{BufferedNonlinearFunction},JuMP._FunctionStorage,MathOptInterface.NLPBoundsPair,Dict{Int64, Vector{Int64}},Vector{JuMP._Derivatives.Linearity},OperatorRegistry,Vector{Float64},NS})   # time: 0.030106
    Base.precompile(Tuple{typeof(MathOptInterface.is_empty),Optimizer{Incremental{GLPK.Optimizer}, Incremental{Ipopt.Optimizer}, DefaultExt}})   # time: 0.0280181
    Base.precompile(Tuple{typeof(initialize!),RelaxCache{2, NS},DirectedTree})   # time: 0.0236066
    Base.precompile(Tuple{typeof(_propagate_constraint!),Evaluator,BufferedNonlinearFunction{1, NS}})   # time: 0.0191075
    Base.precompile(Tuple{typeof(rprop!),Relax,Evaluator,BufferedNonlinearFunction{2, NS}})   # time: 0.0188975
    isdefined(EAGO, Symbol("#157#160")) && Base.precompile(Tuple{getfield(EAGO, Symbol("#157#160")),BufferedNonlinearFunction{2, NS}})   # time: 0.0069586
    Base.precompile(Tuple{typeof(relax!),GlobalOptimizer{Incremental{GLPK.Optimizer}, Incremental{Ipopt.Optimizer}, DefaultExt},BufferedNonlinearFunction{1, NS},Int64,Bool})   # time: 0.005456
    Base.precompile(Tuple{typeof(MathOptInterface.add_variable),Optimizer{Incremental{GLPK.Optimizer}, Incremental{Ipopt.Optimizer}, DefaultExt}})   # time: 0.0048171
    Base.precompile(Tuple{typeof(MathOptInterface.add_constraint),Optimizer{Incremental{GLPK.Optimizer}, Incremental{Ipopt.Optimizer}, DefaultExt},MathOptInterface.SingleVariable,MathOptInterface.LessThan{Float64}})   # time: 0.0040936
    Base.precompile(Tuple{typeof(MathOptInterface.add_constraint),Optimizer{Incremental{GLPK.Optimizer}, Incremental{Ipopt.Optimizer}, DefaultExt},MathOptInterface.SingleVariable,MathOptInterface.GreaterThan{Float64}})   # time: 0.0038925
    Base.precompile(Tuple{typeof(MathOptInterface.set),Optimizer{Incremental{GLPK.Optimizer}, Incremental{Ipopt.Optimizer}, DefaultExt},MathOptInterface.RawParameter,Int64})   # time: 0.0033517
    Base.precompile(Tuple{typeof(MathOptInterface.set),Optimizer{Incremental{GLPK.Optimizer}, Incremental{Ipopt.Optimizer}, DefaultExt},MathOptInterface.NLPBlock,MathOptInterface.NLPBlockData})   # time: 0.0032988
    Base.precompile(Tuple{typeof(MathOptInterface.set),Optimizer{Incremental{GLPK.Optimizer}, Incremental{Ipopt.Optimizer}, DefaultExt},MathOptInterface.ObjectiveFunction{MathOptInterface.SingleVariable},MathOptInterface.SingleVariable})   # time: 0.0022907
    Base.precompile(Tuple{Type{RelaxCache{1, NS}}})   # time: 0.0020124
    Base.precompile(Tuple{Type{RelaxCache{2, NS}}})   # time: 0.0017719
    isdefined(EAGO, Symbol("#133#135")) && Base.precompile(Tuple{getfield(EAGO, Symbol("#133#135")),BufferedNonlinearFunction{2, NS}})   # time: 0.0014442
    Base.precompile(Tuple{typeof(MathOptInterface.set),Optimizer{Incremental{GLPK.Optimizer}, Incremental{Ipopt.Optimizer}, DefaultExt},MathOptInterface.ObjectiveSense,MathOptInterface.OptimizationSense})   # time: 0.0012987
    return
end

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
"""
	prunenode!(node::TreeNode)

Prune node `node` by detaching it from its ancestor. Return pruned `node` and its previous
  ancestor.
"""
function prunenode!(node::TreeNode)
	if node.isroot
		@warn "Trying to prune root: no op."
		return node, TreeNode()
	end
	anc = node.anc
	for (i,c) in enumerate(anc.child)
		if c == node
			deleteat!(anc.child,i)
			break
		end
	end
	if isempty(anc.child)
		anc.isleaf = true
	end
	node.anc = nothing
	node.isroot = true
	return node, anc
end

"""
	prunenode(node::TreeNode)

Prune node `node` by detaching it from its ancestor. Return pruned `node` and previous root `r`. The tree defined by `node` is copied before the operation, and therefore not modified.
"""
function prunenode(node::TreeNode)
	if node.isroot
		@warn "Trying to prune root: no op."
		node_ = deepcopy(node)
		return node_, node_
	end
	node_ = deepcopy(node)
	r = node_findroot(node_)
	anc = node_.anc
	for (i,c) in enumerate(anc.child)
		if c == node_
			deleteat!(anc.child,i)
			break
		end
	end
	node_.anc = nothing
	node_.isroot = true
	return node_, r
end

"""
	prunesubtree!(tree, r::TreeNode)
	prunesubtree!(tree, labellist)

Prune and subtree corresponding to the MRCA of labels in `labellist`.
  Return the root of the subtree as well as its previous direct ancestor.
"""
function prunesubtree!(tree, labellist; clade_only=true)
	if clade_only && !isclade(labellist, tree)
		error("Can't prune non-clade $labellist")
	end
	r = lca(tree, labellist)
	return prunesubtree!(tree, r)
end

function prunesubtree!(tree, r::TreeNode; remove_singletons=true)
	if r.isroot
		@error "Trying to prune root"
	end

	a = r.anc
	delnode(n) = begin
		delete!(tree.lnodes, n.label)
		if n.isleaf
			delete!(tree.lleaves, n.label)
		end
	end
	map!(delnode, r)
	prunenode!(r)
	if remove_singletons
		remove_internal_singletons!(tree, ptau=true)
	end
	return r, a
end

"""
	remove_internal_singletons!(tree; ptau=true)

Remove nodes with one child. Root node is left as is.
If `ptau`, the length of branches above removed nodes is added to the branch length above their children.
"""
function remove_internal_singletons!(tree; ptau=true)
	root = tree.root
	for n in values(tree.lnodes)
		if length(n.child) == 1
			if !n.isroot
				delete_node!(n, ptau=ptau)
				delete!(tree.lnodes, n.label)
			end
		end
	end
	# If root itself is a singleton, delete its child and regraft onto it.
	if length(tree.root.child) == 1 && !tree.root.child[1].isleaf
		n = tree.root.child[1]
		delete_node!(n, ptau=ptau)
		delete!(tree.lnodes, n.label)
	end
	#
	#node2tree!(tree, root)
end

"""
	graftnode!(r::TreeNode, n::TreeNode ; tau=n.tau)

Graft `n` on `r`.
"""
function graftnode!(r::TreeNode, n::TreeNode ; tau=n.tau)
	if !n.isroot || n.anc != nothing
		@error "Trying to graft non-root node $(n)."
	end
	push!(r.child, n)
	r.isleaf = false
	n.anc = r
	n.isroot = false
	n.tau = tau
	return nothing
end

"""
	delete_node!(node::TreeNode; ptau=false)

Delete `node`. If it is an internal node, its children are regrafted on `node.anc`. Returns the new `node.anc`.
If `ptau`, branch length above `node` is added to the regrafted branch.
Otherwise, the regrafted branch's length is unchanged. Return modified `node.anc`.

Note that the previous node will still be in the dictionary `lnodes` and `lleaves` (if a leaf) and the print function will fail on the tree,
to fully remove from the tree and apply the print function use `delete_node!(t::Tree, label; ptau=false)`
"""
function delete_node!(node::TreeNode; ptau=false)
	if node.isroot
		@error "Cannot delete root node"
		error()
	end
	out = node.anc
	if node.isleaf
		prunenode!(node)
	else
		base_tau = node.tau
		for c in reverse(node.child)
			nc = prunenode!(c)[1]
			graftnode!(node.anc, nc, tau = (base_tau*ptau + nc.tau))
		end
		prunenode!(node)
	end
	return out
end

function delete_node!(t::Tree, label; ptau=false)
	delete_node!(t.lnodes[label]; ptau)
	delete!(t.lnodes, label)
	haskey(t.lleaves, label) && delete!(t.lleaves, label)
	remove_internal_singletons!(t; ptau)
	return nothing
end

function delete_null_branches!(node::TreeNode; threshold = 1e-10)
	if !node.isleaf
		if !ismissing(node.tau) && node.tau < threshold && !node.isroot
			nr = delete_node!(node)
			for c in nr.child
				delete_null_branches!(c, threshold=threshold)
			end
		else
			for c in node.child
				delete_null_branches!(c,threshold=threshold)
			end
		end
	elseif !ismissing(node.tau) && node.tau < threshold && !node.isroot
		node.tau = 0.
	end
	return node
end

"""
	delete_null_branches!(tree::Tree; threshold=1e-10)

Delete internal node with branch length smaller than `threshold`. Propagates recursively down the tree. For leaf nodes, set branch length to 0 if smaller than `threshold`.
- If `node` needs not be deleted, call `delete_null_branches!` on its children
- If need be, call `delete_null_branches!` on `node.anc.child`

Calling `delete_null_branches!(node::TreeNode; threshold = 1e-10)` on a node deletes null branches recursively from input `node`, nodes with null branchs will still be in
`lnodes` and `lleaves` (if a leaf) dictionaries. To print the resulting tree use delete_null_branches!(tree::Tree; threshold=1e-10) to remove all null branches from the tree,
or otherwise remove nodes with null branches from `lnodes` and `lleaves`.
"""
function delete_null_branches!(tree::Tree; threshold=1e-10)
	delete_null_branches!(tree.root, threshold=threshold)
	node2tree!(tree, tree.root)
	return nothing
end

function delete_branches!(f, n::TreeNode)
	if !n.isroot && f(n)
		if !n.isleaf
			if !ismissing(n.tau)
				for c in n.child
					if !ismissing(c.tau)
						c.tau += n.tau
					end
				end
			end
			nr = delete_node!(n)
			for c in nr.child
				delete_branches!(f, c)
			end
		else
			n.tau = ismissing(n.tau) ? missing : 0.
		end
	else
		for c in n.child
			delete_branches!(f, c)
		end
	end

	return n
end

"""
	delete_branches!(f, tree::Tree)

Delete branch above node `n` if `f(n)` returns `true` when called on node `n`. When called on `tree` propagates recursively down the tree.
Only when called on `tree` will nodes be additionally deleted from the `lnodes` and `lleaves` (if a leaf) dictionaries (necessary for printing resulting tree).
"""
function delete_branches!(f, tree::Tree)
	delete_branches!(f, tree.root)
	remove_internal_singletons!(tree)
	node2tree!(tree, tree.root)
	return nothing
end

#=
Reroot the tree to which `node` belongs at `node`.
- If `node.isroot`,
- Else if `newroot == nothing`, reroot the tree defined by `node` at `node`. Call `reroot!(node.anc; node)`.
- Else, call `reroot!(node.anc; node)`, then change the ancestor of `node` to be `newroot`.
=#
function reroot!(node::Union{TreeNode,Nothing}; newroot::Union{TreeNode, Nothing}=nothing)
	# Breaking cases
	if node.anc == nothing || node.isroot
		if !(node.anc == nothing && node.isroot)
			@warn "There was a problem with input tree: previous root node has an ancestor."
		elseif newroot != nothing
			i = findfirst(c->c.label==newroot.label, node.child)
			deleteat!(node.child, i)
			node.anc = newroot
			node.tau = newroot.tau
			node.isroot = false
		end
	else # Recursion
		if newroot == nothing
			if node.isleaf
				@warn "Rooting on a leaf node..."
			end
			node.isroot = true
			reroot!(node.anc, newroot=node)
			push!(node.child, node.anc)
			node.anc = nothing
			node.tau = missing
		else
			i = findfirst(c->c.label==newroot.label, node.child)
			deleteat!(node.child, i)
			reroot!(node.anc, newroot=node)
			push!(node.child, node.anc)
			node.anc = newroot
			node.tau = newroot.tau
		end
	end
end
"""
	reroot!(tree::Tree, node::AbstractString)

Reroot `tree` at `tree.lnodes[node]`.
"""
function reroot!(tree::Tree, node::AbstractString)
	reroot!(tree.lnodes[node])
	tree.root = tree.lnodes[node]

	return nothing
end

"""
	add_internal_singleton!(n::TreeNode, a::TreeNode, τ::Real)

Add internal singleton above `n` and below `a`, at heigh `τ` above `n`.
Return the singleton.
"""
function add_internal_singleton!(n::TreeNode, a::TreeNode, τ::Real, label; v = false)
	# Branch length above n and above the singleton
	nτ, sτ = n.tau >= τ ? (τ, n.tau - τ) : (n.tau, 0.)

	s = TreeNode(; tau = sτ, label)
	prunenode!(n)
	n.tau = nτ
	graftnode!(a, s)
	graftnode!(s, n)
	return s
end
function add_internal_singleton!(n::TreeNode, a::TreeNode, τ::Missing, label; v = false)
	@assert ismissing(n.tau)
	prunenode!(n)
	s = TreeNode(; label)
	graftnode!(a, s)
	graftnode!(s, n)
	return s
end


text:
using NBody

function main()
  nb = NBody.read_nbody_json()
  qs = NBody.quartiles(nb)

  s = """
  Quartiles
  r(1/4):   $(@sprintf("%.4g", qs[1]))
  r(2/4):   $(@sprintf("%.4g", qs[2]))
  r(3/4):   $(@sprintf("%.4g", qs[3]))
  r(1):     $(@sprintf("%.4g", qs[4]))

  """
  write(STDERR, s)
  # NBody.write_snapshot(nb)
end

main()


text:
using GalerkinSparseGrids
using Base.Test
using Cubature

#--------------------------------------
# Testing Differentiation
#--------------------------------------

print("Testing differentiation 1-D DG basis... ")

k=3
for l in 2:5
    frefVD = V2Dref(1, k, l);
    frefDV = D2Vref(1, k, l);
    D_op = D_matrix(1, k, l, frefVD, frefDV; scheme="full")
    vcoeffs = vcoeffs_DG(1, k, l, x->cos(2*pi*x[1]); scheme="full")
    dvcoeffs = *(D_op,vcoeffs)
    dict= V2D(1, k, l, vcoeffs; scheme="full")
    ddict = V2D(1, k, l, dvcoeffs; scheme="full")
    err = x->(reconstruct_DG(ddict,[x[1]])+2*pi*sin(2*pi*x[1]))^2
    @test hquadrature(err, 0, 1, abstol=1.0e-10, maxevals=500)[1] < 1/(1<<(k+l-2))
end

println("Test Passed.")

print("Testing differentiation 2-D full DG basis... ")

D = 2
k = 3
for l in 2:5
    D_op = D_matrix(D, 1, k, l; scheme="full")
    vcoeffs = vcoeffs_DG(D, k, l, x->cos(2*pi*x[1])*cos(2*pi*x[2]); scheme="full")
    dvcoeffs = *(D_op, vcoeffs)
    dict= V2D(D, k, l, vcoeffs; scheme="full")
    ddict = V2D(D, k, l, dvcoeffs; scheme="full")
    err = x->(reconstruct_DG(ddict,[x[1],x[2]])+2*pi*sin(2*pi*x[1])*cos(2*pi*x[2]))^2
    @test hcubature(err, [0,0], [1,1], abstol=1.0e-10, maxevals=500)[1] < 1/(1<<(k+l-2))
end

println("Test Passed.")


text:
# Parameter of neural network
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

export Parameter, get_params, get_grads, set_params!, par2vec, vec2par

mutable struct Parameter
    data
    grad
end

"""
    Class for trainable network parameters.

 *Fields:*

 - `Parameter.data`: weights

 - `Parameter.grad`: gradient

"""
Parameter(x) = Parameter(x, nothing)

# Size and length for parameter types
size(x::Parameter) = size(x.data)
length(x::Parameter) = length(x.data)

@Flux.functor Parameter

"""
    clear_grad!(NL::NeuralNetLayer)

or

    clear_grad!(P::AbstractArray{Parameter, 1})

 Set gradients of each `Parameter` in the network layer to `nothing`.
"""
function clear_grad!(P::AbstractArray{Parameter, 1})
    for j=1:length(P)
        P[j].grad = nothing
    end
end

function get_grads(p::Parameter)
    return Parameter(p.grad)
end

function get_grads(pvec::Array{Parameter, 1})
    g = Array{Parameter, 1}(undef, length(pvec))
    for i = 1:length(pvec)
        g[i] = get_grads(pvec[i])
    end
    return g
end

function set_params!(pold::Parameter, pnew::Parameter)
    pold.data = pnew.data
    pold.grad = pnew.grad
end

function set_params!(pold::Array{Parameter, 1}, pnew::Array{Parameter, 1})
    for i = 1:length(pold)
        set_params!(pold[i], pnew[i])
    end
end

## Algebraic utilities for parameters

function dot(p1::Parameter, p2::Parameter)
    return dot(p1.data, p2.data)
end

function norm(p::Parameter)
    return norm(p.data)
end

function +(p1::Parameter, p2::Parameter)
    return Parameter(p1.data+p2.data)
end

function +(p1::Parameter, p2::Float32)
    return Parameter(p1.data+p2)
end

function +(p1::Float32, p2::Parameter)
    return p2+p1
end

function -(p1::Parameter, p2::Parameter)
    return Parameter(p1.data-p2.data)
end

function -(p1::Parameter, p2::Float32)
    return Parameter(p1.data-p2)
end

function -(p1::Float32, p2::Parameter)
    return -(p2-p1)
end

function -(p::Parameter)
    return Parameter(-p.data)
end

function *(p1::Parameter, p2::Float32)
    return Parameter(p1.data*p2)
end

function *(p1::Float32, p2::Parameter)
    return p2*p1
end

function /(p1::Parameter, p2::Float32)
    return Parameter(p1.data/p2)
end

function /(p1::Float32, p2::Parameter)
    return Parameter(p1/p2.data)
end

# Shape manipulation

function par2vec(x::Parameter)
    return vec(x.data), size(x.data)
end

function vec2par(x::Array{Float32, 1}, s::NTuple{N, Int64}) where {N}
    return Parameter(reshape(x, s))
end

function par2vec(x::Array{Parameter, 1})
    v = Array{Float32, 1}(undef, 0)
    s = Array{Any, 1}(undef, 0)
    for i = 1:length(x)
        vi, si = par2vec(x[i])
        v = cat(v, vi; dims=1)
        s = cat(s, si; dims=1)
    end
    return v, s
end

function vec2par(x::Array{Float32, 1}, s::Array{Any, 1})
    xpar = Array{Parameter, 1}(undef, length(s))
    idx_i = 0
    for i = 1:length(s)
        xpar[i] = vec2par(x[idx_i+1:idx_i+prod(s[i])], s[i])
        idx_i += prod(s[i])
    end
    return xpar
end


text:
using JLD2, Base.Test
@test JLD2.Lookup3.hash(UInt8[]) == 0xdeadbeef
@test JLD2.Lookup3.hash(UInt8[], 1, 0, 0xdeadbeef) == 0xbd5b7dde
@test JLD2.Lookup3.hash(UInt8[23]) == 0xa209c931
@test JLD2.Lookup3.hash(UInt8[0]) == 0x8ba9414b
@test JLD2.Lookup3.hash(UInt8[23, 187]) == 0x8ba7a6c9
@test JLD2.Lookup3.hash(UInt8[0, 0]) == 0x62cd61b3
@test JLD2.Lookup3.hash(UInt8[23, 187, 98]) == 0xcebdf4f0
@test JLD2.Lookup3.hash(UInt8[0, 0, 0]) == 0x6bd0060f
@test JLD2.Lookup3.hash(UInt8[23, 187, 98, 217]) == 0x2c88bb51
@test JLD2.Lookup3.hash(UInt8[0, 0, 0, 0]) == 0x049396b8
large_buf = Vector{UInt8}(3093)
for i = 1:length(large_buf)
    large_buf[i] = ((i-1) * 3) % UInt8
end
@test JLD2.Lookup3.hash(large_buf) == 0x1bd2ee7b
@test JLD2.Lookup3.hash(fill!(large_buf, 0)) == 0x930c7afc


text:
module Mazes
export Maze, generate_btree!, generate_sidewinder!

using Graphs, Random

struct Maze
    dim::Tuple # TODO
    graph::SimpleGraph
end
function Maze(dim)
    graph = SimpleGraph(prod(dim))
    return Maze(dim, graph)
end
function Base.show(io::IO, maze::Maze)
    if length(maze.dim) == 2
        empty = "   "
        rows, cols = maze.dim
        println(io, "+" * "---+"^cols)
        for row = 1:rows
            sides = "|"
            floors = "+"
            for col = 1:cols
                # julia uses column major order
                cell = cols * (col - 1) + row
                right = cols * (col) + row
                down = cols * (col - 1) + row + 1
                sides = sides * empty * (
                    if has_edge(maze.graph, cell, right)
                        " "
                    else
                        "|"
                    end
                )
                floors = floors * (
                    if has_edge(maze.graph, cell, down)
                        "   "
                    else
                        "---"
                    end
                ) * "+"
            end
            println(io, sides)
            println(io, floors)
        end
    else
        print("Higher dimensions are not implemented.")
    end
end

function generate_btree!(maze)
    if length(maze.dim) == 2
        rows, cols = maze.dim
        for cell in 1:prod(maze.dim)
            row = mod(cell-1, rows) + 1
            col = fld(cell-1, rows) + 1
            if row + 1 > rows
                linked = cols * (col) + row
            elseif col + 1 > cols
                linked = cols * (col - 1) + row + 1
            else
                down, right = shuffle([0,1])
                linked = cols * (col + right - 1) + row + down
            end
            add_edge!(maze.graph, cell, linked)
        end
    else
        print("Higher dimensions are not implemented.")
    end
end

function generate_sidewinder!(maze) # BUG
    if length(maze.dim) == 2
        rows, cols = maze.dim
        run = []
        for cell in 1:prod(maze.dim)
            row = mod(cell-1, rows) + 1
            col = fld(cell-1, rows) + 1
            if row + 1 > rows
                linked = cols * (col) + row
            elseif col + 1 > cols
                linked = cols * (col - 1) + row + 1
            else
                append!(run, cell)
                forward = rand([true, false])
                if forward
                    linked = cols * (col - 1) + row + 1
                else
                    choice = rand(run)
                    add_edge!(maze.graph, choice, choice+1)
                    empty!(run)
                    continue
                end
            end
            add_edge!(maze.graph, cell, linked)
        end
    else
        print("Higher dimensions are note implemented.")
    end
end

end

max_stars_repo_name: njprov/mazes

max_issues_repo_name: njprov/mazes

max_forks_repo_name: njprov/mazes


text:
### A Pluto.jl notebook ###
# v0.14.6

using Markdown
using InteractiveUtils

# ╔═╡ 0d5cd670-bdf4-11eb-3dfd-b36df6a45072
begin
    import Pkg
    Pkg.activate(mktempdir())
    Pkg.add([
        Pkg.PackageSpec(name="Agents", version="4"),
        Pkg.PackageSpec(name="PlutoUI", version="0.7"),
        Pkg.PackageSpec(name="Plots", version="1"),
        Pkg.PackageSpec(name="DataFrames", version="1"),

    ])
    using Agents
	using PlutoUI
	using Plots
	using Random
	using DataFrames
	PlutoUI.TableOfContents()
end

# ╔═╡ af61b274-831b-4fac-bc84-ab365867ddb3
md"""

# Agent-based Modeling

Agent-based Modeling (ABM) is a simulation method where the autonomous agents interacting with the environment (space) and/or each other by a set of rules.

The most obvious example of ABM is non-player characters (NPCs) in computer games (e.g. GTA V).

ABM is able to model heterogeneously, i.e. it does not require the environment to be well stirred (as opposed to ODEs), continuous (as opposed to to PDEs), nor need the characteristics of each kind of agents to be identical (as opposed to SSAs).

This gives ABM more flexible to model individual behaviors. (e.g. traffic jam, disease spread, molecular interactions)
"""

# ╔═╡ 529ac59e-5cdc-42f9-a2c3-8a36c3ef551d
md"""
## Elements of ABM

To use `Agents.jl`, we need to define:

- The [**space**](https://juliadynamics.github.io/Agents.jl/stable/api/#Available-spaces) where the agents live
- The **agents** with self-defined properties. See [`@agents`](https://juliadynamics.github.io/Agents.jl/stable/api/#@agent-macro)
- The **model** to hold the `space`, the `agent`s, and other parameters (called `properties`)
- The stepping function `step!()` to tell how the model evolve.
"""

# ╔═╡ eac86881-943f-41e9-a085-64ef2ad1554b
md"""
⚠ *WARNING*

Currently (2021-05-26) `abm_plot()` and `abm_video()` [do not work](https://discourse.julialang.org/t/error-no-backend-available-glmakie-cairomakie-wglmakie-when-calling-abm-video-in-interactivedynamics-jl/61626). There might be something wrong between `Agents.jl` and `Makie.jl`.
"""

# ╔═╡ 9668dd2f-42a0-4a1a-80f5-152af1074154
md"""
## Could I do ABM from scratch?
"""

# ╔═╡ 21a958a2-b535-4735-8b90-27ff2580858a
md"""
## Resources

- [Documentation](https://juliadynamics.github.io/Agents.jl/stable/) of `Agents.jl`.
- [sir-julia](https://github.com/epirecipes/sir-julia) : Various implementations of the classical SIR model in Julia.

"""

# ╔═╡ 0d4fd70d-892f-45d9-97a8-5e31e516ae9f
md"""

## Example 1 : Schelling's segregation model

Taken from `Agents.jl` [tutorial](https://juliadynamics.github.io/Agents.jl/stable/examples/schelling/).

- Space  : 2D grid space with a Chebyshev metric. This leads to 8 neighboring positions per position (except at the edges of the grid).
- Agents : They belong to one of two groups (0 or 1).
- Model : Each position of the grid can be occupied by at most one agent.
- For each step
  - If an agent has at least 3 neighbors belonging to the same group, then it is happy.
  - If an agent is unhappy, it keeps moving to new locations until it is happy.

"""

# ╔═╡ 51ce431e-8edf-4caa-8454-03f46bf4d9c3
md"""
To define an agent type, we should make a mutable struct derived from `AbstractAgent`

With 2 mandatory fields
- `id::Int` . The identifier number of the agent.
- `pos` . For agents on a 2D grid, the position field should be a tuple of 2 integers.

On top of that, we could define other properties for the agents.

Alternatively, we can use [`@agent`](https://juliadynamics.github.io/Agents.jl/stable/api/#@agent-macro) macro to let `Agents.jl` set up the mandatory fields for us.

```julia
@agent SchellingAgent GridAgent{2} begin
    mood::Bool
    group::Int
end
```
"""

# ╔═╡ 27dc679f-863e-45e5-b159-0135b3666463
mutable struct SchellingAgent <: AbstractAgent
    id::Int             # The identifier number of the agent
    pos::NTuple{2, Int} # The x, y location of the agent on a 2D grid
    mood::Bool          # whether the agent is happy in its position. (true = happy)
    group::Int          # The group of the agent, determines mood as it interacts with neighbors
end

# ╔═╡ 2c3940e7-7f2d-46bf-a8df-1f4ee9af25fe
md"""
It is recommeded to write a function to make ABMs so that it will be easy to recreate the model and change its parameters.
"""

# ╔═╡ 144bf87e-aac4-419f-b385-6ae853bc64ea
function make_schelling(; numagents = 320,
					  griddims = (20, 20),
		  		      min_to_be_happy = 3,
				      seed = 125)
    space = GridSpace(griddims, periodic = false)
    properties = (min_to_be_happy = 3, )

	# For demos to be reproducible
	# You don't need this in production code
    rng = Random.MersenneTwister(seed)

    model = ABM(
        SchellingAgent, space;
        properties, rng, scheduler = Schedulers.randomly
    )

    # populate the model with agents, adding equal amount of the two types of agents
    # at random positions in the model
    for n in 1:numagents
        agent = SchellingAgent(n, (1, 1), false, n < numagents / 2 ? 1 : 2)
        add_agent_single!(agent, model)
    end
    return model
end

# ╔═╡ d2d1fc91-d5b3-45f0-8237-f7acb9cbe8e4
md"""
And we define a stepping function in the format

`agent_step!(agent, model)`
"""

# ╔═╡ dc52f1f4-ac80-40c0-a3fa-afeaf11352f3
function agent_step!(agent::SchellingAgent, model)
    minhappy = model.min_to_be_happy
    count_neighbors_same_group = 0
    # For each neighbor, get group and compare to current agent's group
    # and increment count_neighbors_same_group as appropriately.
    # Here `nearby_agents` (with default arguments) will provide an iterator
    # over the nearby agents one grid point away, which are at most 8.
    for neighbor in nearby_agents(agent, model)
        if agent.group == neighbor.group
            count_neighbors_same_group += 1
        end
    end
    # After counting the neighbors, decide whether or not to move the agent.
    # If count_neighbors_same_group is at least the min_to_be_happy, set the
    # mood to true. Otherwise, move the agent to a random position.
    if count_neighbors_same_group ≥ minhappy
        agent.mood = true
    else
        move_agent_single!(agent, model)
    end
    return
end

# ╔═╡ d1282e3d-3dd4-4bae-83d5-4c29d27f6c97
begin
	groupcolor(a) = a.group == 1 ? :blue : :orange
	groupmarker(a) = a.group == 1 ? :circle : :rect
	model = make_schelling(griddims = (50, 50), numagents = 1800)

	anim = @animate for i in 1:25
		step!(model, agent_step!)
		plotabm(model; ac = groupcolor, am = groupmarker, as = 5, aspect_ratio=:equal, xlims=(0.0, 50.5), ylims = (0.0, 50.5), title = "Schelling Step $i", size=(600, 600))
	end

	mp4(anim, fps = 5)
end

# ╔═╡ 02a2223b-904b-4980-bb28-35ece350c6d2
md"""
## Appendix
Running environment and some auxillary functions.

**You can use [this helper tool](https://fonsp.com/article-test-3/pkghelper.html) to generate these commands!**
"""

# ╔═╡ Cell order:
# ╠═af61b274-831b-4fac-bc84-ab365867ddb3
# ╠═529ac59e-5cdc-42f9-a2c3-8a36c3ef551d
# ╠═eac86881-943f-41e9-a085-64ef2ad1554b
# ╠═9668dd2f-42a0-4a1a-80f5-152af1074154
# ╠═21a958a2-b535-4735-8b90-27ff2580858a
# ╠═0d4fd70d-892f-45d9-97a8-5e31e516ae9f
# ╠═51ce431e-8edf-4caa-8454-03f46bf4d9c3
# ╠═27dc679f-863e-45e5-b159-0135b3666463
# ╟─2c3940e7-7f2d-46bf-a8df-1f4ee9af25fe
# ╠═144bf87e-aac4-419f-b385-6ae853bc64ea
# ╠═d2d1fc91-d5b3-45f0-8237-f7acb9cbe8e4
# ╠═dc52f1f4-ac80-40c0-a3fa-afeaf11352f3
# ╠═d1282e3d-3dd4-4bae-83d5-4c29d27f6c97
# ╠═02a2223b-904b-4980-bb28-35ece350c6d2
# ╠═0d5cd670-bdf4-11eb-3dfd-b36df6a45072

max_stars_repo_name: stevengogogo/pluto-notebooks

max_issues_repo_name: stevengogogo/pluto-notebooks

max_forks_repo_name: stevengogogo/pluto-notebooks


text:
@testset "784.letter-case-permutation.jl" begin
    @test letter_case_permutation("a1b2") == ["a1b2", "A1b2", "a1B2", "A1B2"]
    @test letter_case_permutation("3z4") == ["3z4", "3Z4"]
    @test letter_case_permutation("12345") == ["12345"]
end


text:
"""
Estimate the bias of a bootstrap sampling.

```julia
bs = bootstrap(mean, randn(20), BasicSampling(100))

bias(bs)
```

"""
bias(t0, t1::AbstractVector) = mean(t1) - t0

bias(bs::BootstrapSample) = [bias(t0, t1) for (t0, t1) in zip(original(bs), straps(bs))]

bias(bs::BootstrapSample, idx::Int) = bias(original(bs, idx), straps(bs, idx))

"""
Estimate the standard error of a bootstrap sampling.

```julia
bs = bootstrap(mean, randn(20), BasicSampling(100))

stderror(bs)
```

"""
stderror(t1::AbstractVector) = std(t1)

stderror(bs::BootstrapSample) = [stderror(t1) for t1 in straps(bs)]

stderror(bs::BootstrapSample, idx::Int) = stderror(straps(bs, idx))

original(bs::BootstrapSample) = bs.t0

original(bs::BootstrapSample, idx::Int) = original(bs)[idx]

straps(bs::BootstrapSample) = bs.t1

straps(bs::BootstrapSample, idx::Int) = straps(bs)[idx]

data(bs::BootstrapSample) = bs.data

model(bs::ParametricBootstrapSample) = bs.model

sampling(bs::BootstrapSample) = bs.sampling


text:
using Test

include("test-c.jl")
include("test-cxx.jl")
include("test-cxx-types.jl")
include("test-f90.jl")
include("test-haskell.jl")
include("test-rust.jl")


text:
struct SeqLogger <: AbstractLogger
    server_url::String
    headers::Dict{String,String}
    min_level::Logging.LogLevel
    event_properties::Ref{String}
    event_batch::Vector{String}
    batch_size::Int
end

"""
    SeqLogger(
        server_url::AbstractString;
        min_level::Logging.LogLevel=Logging.Info,
        api_key::AbstractString="",
        batch_size::Int=10,
        event_properties...
    )

Logger to post log events to a `Seq` log server.

### Inputs
- `server_url` -- `Seq` server url (e.g. `"http://localhost:5341"`)
- `min_level` -- (optional, `default=Logging.Info`) minimal log level to filter the log events
- `api_key` --  (optional, `default=""`) API-key string for registered Applications
- `batch_size` -- (optional, `default=10`) number of log events sent to `Seq` server in single post
- `event_properties` -- (optional) global log event properties

#### Global Log Event Properties
The `SeqLogger` constructor allows to add global log event properties to the logger using   keyword-arguments.
```julia
SeqLogger("http://localhost:5341"; App="DJSON", Env="PROD", Id="24e0d145-d385-424b-b6ec-081aa17d504a")
```

#### Local Log Event Properties
For each individual log event, additional log event properties can be added which
only apply to a single log event.
```julia
@info "Log additional user id {userId}" userId="1"
```
Note: This only works, if the [`Logging.current_logger`](@ref) is of type `SeqLogger` or "contains" a `SeqLogger`.
"""
function SeqLogger(
    server_url::AbstractString;
    min_level::Logging.LogLevel=Logging.Info,
    api_key::AbstractString="",
    batch_size::Int=10,
    event_properties...
)
    url_endpoint = joinurl(server_url, "api/events/raw")
    headers = Dict("Content-Type" => "application/vnd.serilog.clef")
    if !isempty(api_key)
        headers["X-Seq-api_key"] = api_key
    end
    event_properties_str = stringify(; event_properties...)
    return SeqLogger(
        url_endpoint,
        headers,
        min_level,
        Ref(event_properties_str),
        String[],
        batch_size
    )
end

Logging.shouldlog(logger::SeqLogger, arg...) = true
Logging.min_enabled_level(logger::SeqLogger) = logger.min_level
Logging.catch_exceptions(logger::SeqLogger) = false

"""
    Logging.with_logger(@nospecialize(f::Function), logger::SeqLogger)

Extends the function [`Logging.with_logger`](@ref) for `SeqLogger`s.

### Note
After running the function `f`, the `SeqLogger` needs to flush the log events to
make sure that the entire log batch is sent to the `Seq` server.
"""
function Logging.with_logger(@nospecialize(f::Function), seqLogger::SeqLogger)
    result = Base.CoreLogging.with_logstate(f, Base.CoreLogging.LogState(seqLogger))
    flush_events(seqLogger)
    return result
end

const l = ReentrantLock()

"""
    Logging.handle_message(logger::SeqLogger, args...; kwargs...)

Extends the function [`Logging.handle_message`](@ref) for `SeqLogger`s.

### Note
- If the event batch of the `logger` is "full", the log events are flushed.
- A `ReentrantLock` is used to make the `logger` thread-safe.
"""
function Logging.handle_message(logger::SeqLogger, args...; kwargs...)
    lock(l)
    try
        message_args = LoggingExtras.handle_message_args(args...; kwargs...)
        event_str = parse_event_str_from_args(logger, message_args)
        push!(logger.event_batch, event_str)
        if length(logger.event_batch) >= logger.batch_size
            flush_events(logger)
        end
    finally
        unlock(l)
    end
    return nothing
end

"""
    parse_event_str_from_args(logger::SeqLogger, message_args::NamedTuple)::String

Create a log event string from the named tuple `message_args`.
"""
function parse_event_str_from_args(logger::SeqLogger, message_args::NamedTuple)
    line_event_properties = stringify(; _file=message_args.file, _line=message_args.line)
    is_error_log_event = message_args.level == Logging.Error
    # Ignore "back_trace" log event if log level is Error (use special @x option)
    log_event_kwargs = [
        key => value for (key, value) in message_args.kwargs
        if !(is_error_log_event && key == :back_trace)
    ]
    kwarg_event_properties = stringify(; log_event_kwargs...)
    clean_log_msg = replace_invalid_character("$(message_args.message)")
    at_time = "\"@t\":\"$(now())\""
    at_msg = "\"@mt\":\"$(clean_log_msg)\""
    at_level = "\"@l\":\"$(to_seq_level(message_args.level))\""
    default_event_properties = [at_time, at_msg, at_level]
    if is_error_log_event
        back_trace = [value for (key, value) in message_args.kwargs if key == :back_trace]
        clean_exp_msg = isempty(back_trace) ? clean_log_msg : replace_invalid_character(back_trace[1])
        at_exception = "\"@x\":\"$(clean_exp_msg)\""
        push!(default_event_properties, at_exception)
    end
    additonal_event_properties = [
        event for event in (line_event_properties,
                            kwarg_event_properties,
                            logger.event_properties[]) if !isempty(event) ]
    event_str = join([default_event_properties..., additonal_event_properties...], ",")
    return "{$event_str}"
end

"""
    flush_events(logger::SeqLogger)

Post all log events contained in the event batch of `logger` to the `Seq` server and
clear the event batch afterwards.
"""
function flush_events(logger::SeqLogger)
    if !isempty(logger.event_batch)
        eventBatchJson = join(logger.event_batch, "\n")
        empty!(logger.event_batch)
        post_log_events(logger, eventBatchJson)
    end
    return nothing
end

"""
    flush_current_logger()

Post the events in the logger batch event for the logger for the current task,
or the global logger if none is attached to the task.

### Note
In the main moduel of `Atom`, the `current_logger` is `Atom.Progress.JunoProgressLogger()`.
Therefore, if you set `SeqLogger` as a [Logging.global_logger](@ref) in in `Atom`
use [`flush_global_logger`](@ref).
"""
function flush_current_logger()
    current_logger = Logging.current_logger()
    flush_events(current_logger)
    nothing
end

"""
    flush_global_logger()

Post the events in the logger batch event for the global logger.

### Note
If the logger is run with [`Logging.with_logger`](@ref), this is considered a
current logger [Logging.current_logger](@ref) and  [`flush_current_logger`](@ref).
needs to be used.
"""
function flush_global_logger()
    current_logger = Logging.global_logger()
    flush_events(current_logger)
    nothing
end

"""
    post_request(logger::SeqLogger, json_str::AbstractString)

Send POST request with body `json_str` to `Seq` server.
"""
function post_log_events(logger::SeqLogger, json_str::AbstractString)
    return HTTP.request("POST", logger.server_url, logger.headers, json_str)
end

"""
    event_property!(logger::SeqLogger; kwargs...)

Add one or more event properties to the list of global event properties in `logger`.

### Example
```julia
event_property!(logger, user="Me", userId=1)
```
### Note
If a new event property with identical name as an existing on is added with
`event_property!`, the existing property in `new_event_properties` is not
replaced, the new property is just added to `new_event_properties`.
However, this still works since on the `Seq` side the raw post events considers the last property key as the final one if more than one has the same key.
"""
function event_property!(logger::SeqLogger; kwargs...)
    new_event_properties = stringify(; kwargs...)
    if isempty(logger.event_properties[])
        logger.event_properties[] = new_event_properties
    else
        logger.event_properties[] = join([logger.event_properties[],
                                          new_event_properties], ",")
    end
    return nothing
end

"""
    run_with_logger(f::Function, logger::AbstractLogger, args...)

Helper function that applies advanced event logging to the execution of `f(args...)`.

In addition to normal log events (`@info`, ...), the logger catches exceptions
to create error log events. Afterwards, the exception will continue propagation
as if it had not been caught.

### Inputs
- `f` -- function to execute
- `args` -- function arguments

### Retunrs
Result from applying `f(args...)` including thrown exceptions.

### Notes
The function `run_with_logger` applies a new logger to the function call
`f(args...)`. All loggers applied on a "higher level" are replaced for the lifespan
of the function call.

### Example
```julia
logger = ConsoleLogger(stderr, Logging.Info)
run_with_logger(logger, -3) do number
    @info "Compute square root for negative number: \$number"
    sqrt(number)
end
```
"""
function run_with_logger(f::Function, logger::AbstractLogger, args...)
    with_logger(logger) do
        # If there is an error in the function call, catch the exception, create
        # a log event and rethrow the exception.
        try
            return f(args...)
        catch exception
            back_trace = sprint(showerror, exception, catch_backtrace())
            exception_message = sprint(showerror, exception)
            # SeqLogger extract the back_trace `parse_event_str_from_args`
            @error exception_message back_trace=back_trace
            # Loggers including a SeqLogger needs to flash the event batch before rethrowing
            flush_events(logger)
            rethrow()
        end
    end
end


text:
using Test, Tectonic, Tectonic.Biber

tryrun(f) = try; (f(); true); catch err; rethrow(err); end

mktempdir() do dir
    for f in ["test.tex", "bib.tex", "CITATION.bib"]
        cp(joinpath(@__DIR__, f), joinpath(dir, f))
    end
    cd(dir) do
        @testset "Tectonic" begin
            @testset "tectonic" begin
                @test isfile(binary())
                @test version() == v"0.8.0"
                # Compile a file, clean up afterwards.
                @test !isfile("test.pdf")
                @test tryrun() do
                    tectonic() do bin
                        run(`$bin test.tex`)
                    end
                end
                @test isfile("test.pdf")
            end
            @testset "biber" begin
                @test isfile(Biber.binary())
                @test Biber.version() == v"2.14.0"
                @test tryrun() do
                    biber() do _
                        tectonic() do bin
                            run(`$bin bib.tex`)
                        end
                    end
                end
                @test isfile("bib.pdf")
            end
        end
    end
end


text:
# This code is a julia translation of the following project:
# https://bitbucket.org/happyalu/mcep_alpha_calc

# mcepalpha computes appropriate α for a given sampling frequency.
function mcepalpha(fs::Real;
                   start::Float64=0.0,
                   stop::Float64=1.0,
                   step::Float64=0.001,
                   numpoints::Integer=1000)
    α_candidates = start:step:stop
    mel = melscale_vector(fs, numpoints)
    distances = [rms_distance(mel, warping_vector(α, numpoints)) for
                 α in α_candidates]
    return α_candidates[indmin(distances)]
end

function melscale_vector(fs::Real, len::Integer)
    step = (fs / 2.0) / len
    melscalev = 1000.0/log(2)*log(1 + step.*(1:len)./1000.0)
    return melscalev / melscalev[end]
end

function warping_vector(α::Float64, len::Integer)
    step = π / len
    ω = step .* (1:len)
    num = (1-α*α) * sin(ω)
    den = (1+α*α) * cos(ω) - 2*α
    warpfreq = atan(num./den)
    warpfreq[warpfreq .< 0] += π
    return warpfreq / warpfreq[end]
end

rms_distance(v1, v2) = sumabs2(v1 - v2) / length(v1)


text:
@testset "Constraints" begin
    T = 5
    nx = 2
    nu = 1
    nw = 0
    dim_x = [nx for t = 1:T]
    dim_u = [nu for t = 1:T-1]
    dim_w = [nw for t = 1:T]
    x = [rand(dim_x[t]) for t = 1:T]
    u = [[rand(dim_u[t]) for t = 1:T-1]..., zeros(0)]
    w = [rand(dim_w[t]) for t = 1:T]

    ct = (x, u, w) -> [-ones(nx) - x; x - ones(nx)]
    cT = (x, u, w) -> x

    cont = StageConstraint(ct, nx, nu, nw, [t for t = 1:T-1], :inequality)
    conT = StageConstraint(cT, nx, 0, nw, [T], :equality)

    cons = [cont, conT]
    nc = DirectTrajectoryOptimization.num_con(cons)
    nj = DirectTrajectoryOptimization.num_jac(cons)
    idx_c = DirectTrajectoryOptimization.constraint_indices(cons)
    idx_j = DirectTrajectoryOptimization.jacobian_indices(cons)
    c = zeros(nc)
    j = zeros(nj)
    cont.val(c[idx_c[1]], x[1], u[1], w[1])
    conT.val(c[idx_c[T]], x[T], u[T], w[T])

    DirectTrajectoryOptimization.eval_con!(c, idx_c, cons, x, u, w)
    # info = @benchmark DirectTrajectoryOptimization.eval_con!($c, $idx_c, $cons, $x, $u, $w)

    @test norm(c - vcat([ct(x[t], u[t], w[t]) for t = 1:T-1]..., cT(x[T], u[T], w[T]))) < 1.0e-8
    DirectTrajectoryOptimization.eval_jac!(j, idx_j, cons, x, u, w)
    # info = @benchmark DirectTrajectoryOptimization.eval_jac!($j, $idx_j, $cons, $x, $u)

    dct = [-I zeros(nx, nu); I zeros(nx, nu)]
    dcT = Diagonal(ones(nx))
    dc = cat([dct for t = 1:T-1]..., dcT, dims=(1,2))
    sp = DirectTrajectoryOptimization.sparsity_jacobian(cons, dim_x, dim_u)
    j_dense = zeros(nc, sum(dim_x) + sum(dim_u))
    for (i, v) in enumerate(sp)
        j_dense[v[1], v[2]] = j[i]
    end
    @test norm(j_dense - dc) < 1.0e-8
end


text:
include("TypeDef.jl")
include("initializations.jl")
# include("convolve.jl")
include("padding.jl")
# include("pooling.jl")
# include("layerForProp.jl")
# include("layerBackProp.jl")
include("layerUpdateParams.jl")
include("chain.jl")
include("img2col.jl")
include("unroll.jl")
include("dimensions.jl")
# include("img2colConvolve.jl")
# include("fastPooling.jl")
# include("NNConv.jl")

###
include("parallelLayerForProp.jl")
include("parallelLayerBackProp.jl")


text:
const test_checkpoint_inputfile = "input_vals_channel_dg.jl"
#const test_checkpoint_moddict = Dict{String, Any}("Flux_name" => "RoeFlux", "use_DG" => true, "new_fname" => "input_vals_channel_dg")

mutable struct TestCheckpointData <: AbstractCheckpointData
  a::Int
  b::Array{Float64, 1}
end

function test_checkpoint(mesh, sbp, eqn, opts)
  @testset "----- Testing Checkpointer -----" begin
    myrank = mesh.myrank
    rand!(eqn.q_vec)

    chkpointer = Checkpointer(myrank, 2, "abc")

    chkpoint_data = TestCheckpointData(42, rand(10))

    # check initialization
    @test ( chkpointer.ncheckpoints )== 2
    @test ( chkpointer.paths[2] )== joinpath(pwd(), "abc_checkpoint2")
    @test ( isdir(joinpath(pwd(), chkpointer.paths[2])) )== true
    for i=1:chkpointer.ncheckpoints
      @test ( chkpointer.status[i] )== Utils.CheckpointFree
      @test ( chkpointer.history[i] )== -1
    end

    @test ( countFreeCheckpoints(chkpointer) )== 2
    @test ( getNextFreeCheckpoint(chkpointer) )== 1

    # save a checkpoint
    @test ( saveNextFreeCheckpoint(chkpointer, mesh, sbp, eqn, opts, chkpoint_data) )== 1

    @test ( countFreeCheckpoints(chkpointer) )== 1
    @test ( getNextFreeCheckpoint(chkpointer) )== 2
    @test ( getOldestCheckpoint(chkpointer) )== 1
    @test ( getLastCheckpoint(chkpointer) )== 1
    @test ( chkpointer.status[1] )== Utils.CheckpointUsed
    @test ( chkpointer.history[1] )== 1

    q1 = copy(eqn.q_vec)
    scale!(eqn.q_vec, 2)
    q2 = copy(eqn.q_vec)
    @test ( saveNextFreeCheckpoint(chkpointer, mesh, sbp, eqn, opts, chkpoint_data) )== 2
    @test ( countFreeCheckpoints(chkpointer) )== 0
    @test ( chkpointer.status[2] )== Utils.CheckpointUsed

    @test ( getOldestCheckpoint(chkpointer) )== 1
    @test_throws Exception  saveNextFreeCheckpoint(chkpointer, mesh, sbp, eqn, opts, chkpoint_data)

    # try loading a checkpoint
    scale!(eqn.q_vec, 2)
    loadLastCheckpoint(chkpointer, mesh, sbp, eqn, opts)
    @test isapprox( norm(eqn.q_vec - q2), 0.0) atol=1e-14

    chkpointer2 = Checkpointer(opts, mesh.myrank)

    @test ( chkpointer.ncheckpoints )== chkpointer2.ncheckpoints
    @test isapprox( norm(chkpointer.history - chkpointer2.history), 0.0) atol=1e-14
    @test isapprox( norm(chkpointer.status - chkpointer2.status), 0.0) atol=1e-14
    for i=1:chkpointer.ncheckpoints
      @test ( chkpointer.paths[i] )== chkpointer2.paths[i]
    end

    chkpoint_data2 = readCheckpointData(chkpointer, getLastCheckpoint(chkpointer), mesh.myrank)

    @test ( chkpoint_data.a )== chkpoint_data2.a
    @test isapprox( norm(chkpoint_data.b - chkpoint_data2.b ), 0.0) atol=1e-14

    # corrupt the flagfile to force the loading of the earlier checkpoint
    Utils.deleteFlagFile(chkpointer, getLastCheckpoint(chkpointer))
    opts["most_recent_checkpoint"] = 1
    opts["most_recent_checkpoint_path"] = chkpointer.paths[1]

    chkpointer2 = Checkpointer(opts, mesh.myrank)
    @test ( countFreeCheckpoints(chkpointer2) )== 1
    @test ( getLastCheckpoint(chkpointer2) )== 1

    # test with more checkpoints
    chkpointer = Checkpointer(myrank, 5)
    saveNextFreeCheckpoint(chkpointer, mesh, sbp, eqn, opts, chkpoint_data)
    saveNextFreeCheckpoint(chkpointer, mesh, sbp, eqn, opts, chkpoint_data)
    saveNextFreeCheckpoint(chkpointer, mesh, sbp, eqn, opts, chkpoint_data)
    saveNextFreeCheckpoint(chkpointer, mesh, sbp, eqn, opts, chkpoint_data)

    freeCheckpoint(chkpointer, 2)
    @test ( countFreeCheckpoints(chkpointer) )== 2
    @test ( getNextFreeCheckpoint(chkpointer) )== 2
    @test ( getLastCheckpoint(chkpointer) )== 4
    @test ( getOldestCheckpoint(chkpointer) )== 1

    saveNextFreeCheckpoint(chkpointer, mesh, sbp, eqn, opts, chkpoint_data)
    @test ( countFreeCheckpoints(chkpointer) )== 1
    @test ( getNextFreeCheckpoint(chkpointer) )== 5
    @test ( getLastCheckpoint(chkpointer) )== 2
    @test ( getOldestCheckpoint(chkpointer) )== 1

  end

  return nothing
end

add_func2!(EulerTests, test_checkpoint, test_dg_inputfile, [TAG_SHORTTEST, TAG_CHECKPOINT])


text:
##
using Revise
using GLMakie
import Pkg
Pkg.activate(".")
##
import CSV
using DataFrames
using StaticArrays
import Statistics
using NearestNeighbors
using ProgressMeter
##
filepath = "./Data/STORM/desmin_alphaactinin_600nm.csv"
localization_df = DataFrame(CSV.File(filepath))
##
function subset_xy_bounding_box(localization_df::DataFrame, xlim, ylim)
    selection_criteria_x = (xlim[1] .< localization_df.x .< xlim[2])
    selection_criteria_y = (ylim[1] .< localization_df.y .< ylim[2])
    selection_criteria = selection_criteria_x .& selection_criteria_y
    return localization_df[selection_criteria, :]
end

function split_dataframe_per_probe(localization_df::DataFrame)
    probes = sort!(Vector(localization_df[:, :probe]))
    probe_ids = unique(probes)
    return [localization_df[localization_df.probe .== id, :] for id in probe_ids]
end

function pointcloud_from_dataframe(localization_df::DataFrame)
    point_cloud = @inbounds @views [SVector{3, Float32}(localization_df[i, [:x, :y, :z]]...) for i = 1:size(localization_df)[1]]
    return point_cloud
end
##
xlim, ylim = (25000, 30000), (10000, 15000)
df = subset_xy_bounding_box(localization_df, xlim, ylim)
df_list =  split_dataframe_per_probe(df)
point_clouds = pointcloud_from_dataframe.(df_list)
##
import StatsBase: fit, Histogram
function point_to_set_correlation(point_coordinates, set_coordinates, Δr)
    set_tree = KDTree(set_coordinates)
    _, point_to_set_dists = nn(set_tree, point_coordinates)
    r = collect(0.0:Δr:maximum(point_to_set_dists)+Δr)

    point_to_set_hist = fit(Histogram, point_to_set_dists, r, closed = :left)
    return r[1:end-1], point_to_set_hist.weights
##
import RecursiveArrayTools: VectorOfArray
using Clustering

point_matrix = reshape(reinterpret(Float32, point_clouds[1]), (3, length(point_clouds[1])))
##
point_clouds[1][1] == point_matrix[:, 1]

#size(point_matrix)
##
clustering_result = dbscan(point_matrix, 5, min_cluster_size = 2)
##
VectorOfArray isa AbstractMatrix
##
import StatsBase: fit, Histogram
bin_size = 20
point_set_dist = desmin_pt_alphaactinin_set_dist
distances = collect(0.0:bin_size:maximum(point_set_dist)+bin_size)
h = fit(Histogram, point_set_dist, distances, closed=:left)
##
let
    fig = Figure(resolution = (1200, 600))
    ax1 = GLMakie.Axis(fig[1, 1], title = "Point(desmin)-to-Set(α-actinin) correlation function: Δr = $(bin_size)nm; x ∈ $(Tuple(xlim./1000)) μm; y ∈ $(Tuple(ylim./1000)) μm", xlabel = "Distance (nm)", ylabel = "Correlation (a.u.)")
    #ax2 = GLMakie.Axis(fig[1, 2], yscale = log10, title = "Point(desmin)-to-Set(α-actinin) correlation function (Δr = $(bin_size)nm)", xlabel = "Distance (pixels)", ylabel = "Correlation (a.u.)")
    scatterlines!(ax1, distances[1:end-1], h.weights)
    ax1.xticks = LinearTicks(10)
    #scatterlines!(ax2, distances[1:end-1], h.weights)
    fig
end