text:
using Test

##
model_pendulum, obj_uncon_pendulum = TrajectoryOptimization.Dynamics.pendulum!

u_min_pendulum = -4
u_max_pendulum = 2
x_min_pendulum = [-1;-3]
x_max_pendulum = [20; 20]

# -Constrained objective
obj_con_pendulum = ConstrainedObjective(obj_uncon_pendulum, u_min=u_min_pendulum, u_max=u_max_pendulum, x_min=x_min_pendulum, x_max=x_max_pendulum)
obj_con_pendulum.tf = 3.0
model = model_pendulum
obj = obj_con_pendulum
obj.Qf[:,:] = 0.0*Matrix(I,2,2)
# Solver
intergration = :rk4
N = 25
solver1 = Solver(model,obj,integration=intergration,N=N)

# control and state constraints and indices
p,pI,pE = get_num_constraints(solver1)
pu = 2
ps = 4
pu_idx = Array(1:2)
ps_idx = Array(3:6)
pN = model.n

X0 = line_trajectory(solver1)
U0 = ones(solver1.model.m,solver1.N)

solver1.opts.constraint_tolerance = 1e-2
solver1.opts.active_set_flag = false
results2, stat2 = solve(solver1,U0)
max_violation(results2)
cost(solver1,results2)
plot(to_array(results2.X)')
plot(to_array(results2.U)')
plot(to_array(results2.λ)')

results2.λN

solver1.opts.active_set_flag = true
solver1.opts.active_constraint_tolerance = 0.0
update_constraints!(results2,solver1)
update_jacobians!(results2, solver1)
Δv = _backwardpass_active_set!(results2,solver1)

results2.Kλ
results2.dλ
results2.λN
# J = forwardpass!(results2, solver1, Δv)
rollout_active_set!(results2,solver1,0.01)
# results2.X .= deepcopy(results2.X_)
# results2.U .= deepcopy(results2.U_)
update_constraints!(results2,solver1,results2.X_,results2.U_)
cost(solver1,results2,results2.X_,results2.U_)

results2.X_[end]
_cost(solver1,results2,results2.X_,results2.U_)
cost_constraints(solver1, results2)
plot(to_array(results2.X_)')
plot(to_array(results2.U_)')
plot!(to_array(results2.λ)')
results2.λN
max_violation(results2)

a = 1
function _backwardpass_active_set!(res::SolverVectorResults,solver::Solver)
    n,m,N = get_sizes(solver)
    m̄,mm = get_num_controls(solver)
    pN = n
    W = solver.obj.Q; R = solver.obj.R; Wf = solver.obj.Qf; xf = solver.obj.xf

    solver.state.minimum_time ? R_minimum_time = solver.opts.R_minimum_time : nothing
    solver.state.infeasible ? R_infeasible = solver.opts.R_infeasible*Matrix(I,n,n) : nothing

    dt = solver.dt

    # pull out values from results
    X = res.X; U = res.U; K = res.K; d = res.d; S = res.S; s = res.s
    L = res.L; l = res.l

    # Useful linear indices
    idx = Array(1:2*(n+mm))
    Idx = reshape(Array(1:(n+mm)^2),n+mm,n+mm)

    x_idx = idx[1:n]
    u_idx = idx[n+1:n+m]
    ū_idx = idx[n+1:n+mm]

    S = [k != N ? zeros(n+ps,n+ps) : zeros(n+pN,n+pN) for k = 1:N]
    s = [k != N ? zeros(n+ps) : zeros(n+pN) for k = 1:N]

    KK = res.Kλ
    dd = res.dλ
    KK[:] = [k != N-1 ? zeros(pu+ps,n) : zeros(pu+pN,n) for k = 1:N-1]
    dd[:] = [k != N-1 ? zeros(pu+ps) : zeros(pu+pN) for k = 1:N-1]

    # Boundary Conditions
    C = res.C; Iμ = res.Iμ; λ = res.λ; CxN = res.Cx_N

    Sxx = Wf + CxN'*res.IμN*CxN
    Sxλ = CxN'
    Sλλ = zeros(pN,pN)

    Sx = Wf*(X[N] - xf) + CxN'*(res.IμN*res.CN + res.λN)
    Sλ = res.CN

    S[N][1:n,1:n] = Sxx
    S[N][1:n,n+1:n+pN] = Sxλ
    S[N][n+1:n+pN,1:n] = Sxλ'
    S[N][n+1:n+pN,n+1:n+pN] = Sλλ

    s[N][1:n] = Sx
    s[N][n+1:n+pN] = Sλ

    # Initialize expected change in cost-to-go
    Δv = zeros(2)

    # Backward pass
    k = N-1
    while k >= 1
        solver.state.minimum_time ? dt = U[k][m̄]^2 : nothing

        x = X[k]
        u = U[k][1:m]

        ℓ1 = ℓ(x,u,W,R,xf)
        ℓx = W*(x - xf)
        ℓu = R*u
        ℓxx = W
        ℓuu = R
        ℓux = zeros(m,n)

        # Assemble expansion

        Lx = dt*ℓx
        Lu = dt*ℓu
        Lxx = dt*ℓxx
        Luu = dt*ℓuu
        Lux = dt*ℓux

        # Minimum time expansion components
        if solver.state.minimum_time
            h = U[k][m̄]

            l[n+m̄] = 2*h*(ℓ1 + R_minimum_time)

            tmp = 2*h*ℓu
            L[u_idx,m̄] = tmp
            L[m̄,u_idx] = tmp'
            L[m̄,m̄] = 2*(ℓ1 + R_minimum_time)

            L[m̄,x_idx] = 2*h*ℓx'
        end

        # Infeasible expansion components
        if solver.state.infeasible
            l[n+m̄+1:n+mm] = R_infeasible*U[k][m̄+1:m̄+n]
            L[n+m̄+1:n+mm,n+m̄+1:n+mm] = R_infeasible
        end

        # Final expansion terms
        if solver.state.minimum_time || solver.state.infeasible
            l[u_idx] = Lu

            L[u_idx,u_idx] = Luu
            L[u_idx,x_idx] = Lux

            Lu = l[n+1:n+mm]
            Luu = L[n+1:n+mm,n+1:n+mm]
            Lux = L[n+1:n+mm,1:n]
        end

        # Constraints
        Cx, Cu = res.Cx[k][ps_idx,:], res.Cu[k][pu_idx,:]
        Lx += Cx'*(Iμ[k][ps_idx,ps_idx]*C[k][ps_idx] + λ[k][ps_idx])
        Lu += Cu'*(Iμ[k][pu_idx,pu_idx]*C[k][pu_idx] + λ[k][pu_idx])
        Lxx += Cx'*Iμ[k][ps_idx,ps_idx]*Cx
        Luu += Cu'*Iμ[k][pu_idx,pu_idx]*Cu
        # Lux += Cu'*Iμ[k]*Cx

        # Compute gradients of the dynamics
        fdx, fdu = res.fdx[k], res.fdu[k]

        # Gradients and Hessians of Taylor Series Expansion of Q
        Qx = Lx + fdx'*s[k+1][1:n]
        Qu = Lu + fdu'*s[k+1][1:n]
        Qxx = Lxx + fdx'*S[k+1][1:n,1:n]*fdx
        Quu = Luu + fdu'*S[k+1][1:n,1:n]*fdu
        Qux = Lux + fdu'*S[k+1][1:n,1:n]*fdx

        if solver.opts.bp_reg_type == :state
            Quu_reg = Quu + res.ρ[1]*fdu'*fdu
            Qux_reg = Qux + res.ρ[1]*fdu'*fdx
        elseif solver.opts.bp_reg_type == :control
            Quu_reg = Quu + res.ρ[1]*I
            Qux_reg = Qux
        end

        # Regularization
        if !isposdef(Hermitian(Array(Quu_reg)))  # need to wrap Array since isposdef doesn't work for static arrays
            # increase regularization
            regularization_update!(res,solver,:increase)

            Sxx = S[N][1:n,1:n]
            Sxλ = S[N][1:n,n+1:n+pN]
            Sλλ = S[N][n+1:n+pN,n+1:n+pN]

            Sx = s[N][1:n]
            Sλ = s[N][n+1:n+pN]
            # reset backward pass
            k = N-1
            Δv[1] = 0.
            Δv[2] = 0.
            continue
        end

        # linearized dynamics for cost-to-go
        S̄xx = fdx'*Sxx*fdx
        S̄xλ = fdx'*Sxλ
        S̄xu = fdx'*Sxx*fdu
        S̄λλ = Sλλ
        S̄λu = Sxλ'*fdu
        S̄uu = fdu'*Sxx*fdu

        S̄x = fdx'*Sx
        S̄λ = Sλ
        S̄u = fdu'*Sx

        # get active set constraints and Jacobians
        active_set_idx_1u = res.active_set[k][pu_idx] # active set indices at index k
        active_set_idx_1s = res.active_set[k][ps_idx] # active set indices at index k

        k != N-1 ? active_set_idx_2 = res.active_set[k+1][ps_idx] : active_set_idx_2 = ones(Bool,pN) # active set indices at index k+1

        c_active_set_1u = res.C[k][pu_idx][active_set_idx_1u] # active set constraints at index k
        c_active_set_1s = res.C[k][ps_idx][active_set_idx_1s] # active set constraints at index k
        k != N-1 ? c_active_set_2 = S̄λ[active_set_idx_2] : c_active_set_2 = S̄λ # active set constraints at index k+1
        c_set_1u = res.C[k][pu_idx]
        c_set_1s = res.C[k][ps_idx]
        c_set_2 = S̄λ

        cx_active_set_1 = res.Cx[k][ps_idx,:][active_set_idx_1s,:] # active set Jacobians at index k
        cx_active_set_2 = (S̄xλ[:,active_set_idx_2])' # active set Jacobians at index k+1

        cx_set_1 = res.Cx[k][ps_idx,:]
        cx_set_2 = S̄xλ'

        cu_active_set_1 = res.Cu[k][pu_idx,:][active_set_idx_1u,:] # active set Jacobians at index k
        cu_active_set_2 = S̄λu[active_set_idx_2,:] # active set Jacobians at index k+1
        cu_set_1 = res.Cu[k][pu_idx,:]
        cu_set_2 = S̄λu

        # calculate control and multiplier gains
        tmp1 = [cu_active_set_1;cu_active_set_2]
        tmp2 = tmp1*(Quu_reg\tmp1')
        tmp3 = [zeros(length(c_active_set_1u),n);cx_active_set_2]
        tmp4 = [c_active_set_1u;c_active_set_2]
        tmp5 = [cu_set_1; cu_set_2]

        k != N-1 ? gains_active_idx = [res.active_set[k][pu_idx]; res.active_set[k+1][ps_idx]] : gains_active_idx = [res.active_set[k][pu_idx]; ones(Bool,pN)]

        if sum(gains_active_idx) != 0
            println("Iteration $k")
            println(tmp1)
            println(tmp2)
            println(tmp3)
            println(tmp4)
            # println(inv(tmp2))
            println("Cu")
            println(res.Cu[k])
            println("Cx")
            println(res.Cx[k])
            println("S̄λu:")
            println(S̄λu)
            println("Cx_N")
            println(res.Cx_N)
            println("Quu")
            println(Quu)
            println(Quu_reg)
            println(rank(tmp1))
            println(rank(tmp2))
            println(cond(tmp2))

            # println(res.Cx_N*fdu*inv(Quu_reg)*fdu'*res.Cx_N')
            iter = 0
            res.ρ[1] = 0.
            while rank(tmp2) < size(tmp2,1)
                tmp2 += res.ρ[1]*Matrix(I,size(tmp2,1),size(tmp2,1))
                regularization_update!(res,solver,:increase)

                iter += 1
                iter == 100 ? error("couldn't fix tmp2") : nothing
            end
            iter = 0

            KK[k][gains_active_idx,:] = -tmp2\(tmp1*(Quu_reg\Qux_reg) - tmp3)
            dd[k][gains_active_idx] = -tmp2\(tmp1*(Quu_reg\Qu) - tmp4)

            K[k] = -Quu_reg\(Qux_reg + tmp5'*KK[k])
            d[k] = -Quu_reg\(Qu + tmp5'*dd[k])

            println("KK: $(KK[k])")
            println("dd: $(dd[k])")
        end

        k != N-1 ? KKs = KK[k][ps_idx,:] : KKs = KK[k][Array(pu_idx[end]+1:pu_idx[end]+pN),:]
        KKu = KK[k][pu_idx,:]
        k != N-1 ? dds = dd[k][ps_idx] : dds = dd[k][Array(pu_idx[end]+1:pu_idx[end]+pN)]
        ddu = dd[k][pu_idx]

        Ŝxx = S̄xx + S̄xu*K[k] + S̄xλ*KKs + K[k]'*S̄xu' + K[k]'*S̄uu*K[k] + K[k]'*S̄λu'*KKs + KKs'*S̄xλ' + KKs'*S̄λu*K[k] + KKs'*S̄λλ*KKs
        Ŝx = S̄xu*d[k] + S̄xλ*dds + K[k]'*S̄uu*d[k] + KKs'*S̄λu*d[k] + K[k]'*S̄λu'*dds + KKs'*S̄λλ*dds + S̄x + K[k]'*S̄u + KKs'*S̄λ

        Δv[1] += S̄u'*d[k] + S̄λ'*dds
        Δv[2] += 0.5*d[k]'*S̄uu*d[k] + 0.5*d[k]'*S̄λu'*dds + 0.5*dds'*S̄λu*d[k] + 0.5*dds'*S̄λλ*dds

        Lλ = c_set_1s
        Lϕ = c_set_1u
        Luϕ = cu_set_1'
        Lxλ = cx_set_1'

        L̄xx = Lxx + Lux'*K[k] + K[k]'*Luu*K[k] + K[k]'*Luϕ*KKu + KKu'*Luϕ'*K[k] + K[k]'*Lux
        L̄xλ = Lxλ
        L̄x = Lux'*d[k] + K[k]'*Luu*d[k] + K[k]'*Luϕ*ddu + KKu'*Luϕ'*d[k] + Lx + K[k]'*Lu + KKu'*Lϕ
        L̄λ = Lλ

        Sxx = L̄xx + Ŝxx
        Sxλ = L̄xλ
        Sλλ = zeros(ps,ps)

        Sx = L̄x + Ŝx
        Sλ = L̄λ

        S[k][1:n,1:n] = Sxx
        S[k][1:n,n+1:n+ps] = Sxλ
        S[k][n+1:n+ps,1:n] = Sxλ'
        S[k][n+1:n+ps,n+1:n+ps] = Sλλ

        s[k][1:n] = Sx
        s[k][n+1:n+ps] = Sλ

        Δv[1] += Lu'*d[k] + Lϕ'*ddu
        Δv[2] += 0.5*d[k]'*Luu*d[k] + 0.5*d[k]'*Luϕ*ddu + 0.5*ddu'*Luϕ'*d[k]

        k = k - 1;
    end

    # decrease regularization after backward pass
    regularization_update!(res,solver,:decrease)

    return Δv
end

function rollout_active_set!(res::SolverVectorResults,solver::Solver,alpha::Float64)
    n,m,N = get_sizes(solver)
    m̄,mm = get_num_controls(solver)
    dt = solver.dt

    X = res.X; U = res.U; K = res.K; d = res.d; X_ = res.X_; U_ = res.U_

    X_[1] = solver.obj.x0;

    if solver.control_integration == :foh
        b = res.b
        du = alpha*d[1]
        U_[1] = U[1] + du
        dv = zero(du)
    end

    if solver.opts.active_set_flag
        λ = res.λ
        Kλ = res.Kλ
        Mλ = res.Mλ
        dλ = res.dλ

        δλ = alpha*dλ[1]
        λ[1][pu_idx] += δλ[pu_idx]
        λ[2][ps_idx] += δλ[ps_idx]
    end

    for k = 2:N
        δx = X_[k-1] - X[k-1]

        if solver.control_integration == :foh
            dv = K[k]*δx + b[k]*du + alpha*d[k]
            U_[k] = U[k] + dv
            solver.state.minimum_time ? dt = U_[k-1][m̄]^2 : nothing
            solver.fd(X_[k], X_[k-1], U_[k-1][1:m], U_[k][1:m], dt)
            du = dv
        else
            U_[k-1] = U[k-1] + K[k-1]*δx + alpha*d[k-1]
            solver.state.minimum_time ? dt = U_[k-1][m̄]^2 : nothing
            solver.fd(X_[k], X_[k-1], U_[k-1][1:m], dt)
        end

        # if solver.opts.active_set_flag
        #     δλ = Kλ[k-1]*δx + alpha*dλ[k-1]
        #     if k != N
        #         λ[k-1][pu_idx] += δλ[pu_idx]
        #         λ[k][ps_idx] += δλ[ps_idx]
        #     else
        #         λ[k-1][pu_idx] += δλ[pu_idx]
        #         res.λN[1:pN] += δλ[Array(pu_idx[end]+1:pu_idx[end]+pN)]
        #     end
        # end
        solver.state.infeasible ? X_[k] += U_[k-1][m̄.+(1:n)] : nothing

        # Check that rollout has not diverged
        if ~(norm(X_[k],Inf) < solver.opts.max_state_value && norm(U_[k-1],Inf) < solver.opts.max_control_value)
            return false
        end
    end

    # Calculate state derivatives and midpoints
    if solver.control_integration == :foh
        calculate_derivatives!(res, solver, X_, U_)
        calculate_midpoints!(res, solver, X_, U_)
    end

    # Update constraints
    update_constraints!(res,solver,X_,U_)

    return true
end

a = 1

# function _backwardpass_active_set!(res::SolverVectorResults,solver::Solver)
#     n,m,N = get_sizes(solver)
#     m̄,mm = get_num_controls(solver)
#
#     W = solver.obj.Q; R = solver.obj.R; Wf = solver.obj.Qf; xf = solver.obj.xf
#
#     solver.state.minimum_time ? R_minimum_time = solver.opts.R_minimum_time : nothing
#     solver.state.infeasible ? R_infeasible = solver.opts.R_infeasible*Matrix(I,n,n) : nothing
#
#     dt = solver.dt
#
#     # pull out values from results
#     X = res.X; U = res.U; K = res.K; d = res.d; S = res.S; s = res.s
#     L = res.L; l = res.l
#
#     # Useful linear indices
#     idx = Array(1:2*(n+mm))
#     Idx = reshape(Array(1:(n+mm)^2),n+mm,n+mm)
#
#     x_idx = idx[1:n]
#     u_idx = idx[n+1:n+m]
#     ū_idx = idx[n+1:n+mm]
#
#     # Boundary Conditions
#     S[N] = Wf
#     s[N] = Wf*(X[N] - xf)
#
#     # Initialize expected change in cost-to-go
#     Δv = zeros(2)
#
#     # Terminal constraints
#     if res isa ConstrainedIterResults
#         C = res.C; Iμ = res.Iμ; λ = res.λ
#         CxN = res.Cx_N
#         S[N] += CxN'*res.IμN*CxN
#         s[N] += CxN'*(res.IμN*res.CN + res.λN)
#     end
#
#     # Backward pass
#     k = N-1
#     while k >= 1
#         solver.state.minimum_time ? dt = U[k][m̄]^2 : nothing
#
#         x = X[k]
#         u = U[k][1:m]
#
#         ℓ1 = ℓ(x,u,W,R,xf)
#         ℓx = W*(x - xf)
#         ℓu = R*u
#         ℓxx = W
#         ℓuu = R
#         ℓux = zeros(m,n)
#
#         # Assemble expansion
#
#         Lx = dt*ℓx
#         Lu = dt*ℓu
#         Lxx = dt*ℓxx
#         Luu = dt*ℓuu
#         Lux = dt*ℓux
#
#         # Minimum time expansion components
#         if solver.state.minimum_time
#             h = U[k][m̄]
#
#             l[n+m̄] = 2*h*(ℓ1 + R_minimum_time)
#
#             tmp = 2*h*ℓu
#             L[u_idx,m̄] = tmp
#             L[m̄,u_idx] = tmp'
#             L[m̄,m̄] = 2*(ℓ1 + R_minimum_time)
#
#             L[m̄,x_idx] = 2*h*ℓx'
#         end
#
#         # Infeasible expansion components
#         if solver.state.infeasible
#             l[n+m̄+1:n+mm] = R_infeasible*U[k][m̄+1:m̄+n]
#             L[n+m̄+1:n+mm,n+m̄+1:n+mm] = R_infeasible
#         end
#
#         # Final expansion terms
#         if solver.state.minimum_time || solver.state.infeasible
#             l[u_idx] = Lu
#
#             L[u_idx,u_idx] = Luu
#             L[u_idx,x_idx] = Lux
#
#             Lu = view(l,n+1:n+mm)
#             Luu = view(L,n+1:n+mm,n+1:n+mm)
#             Lux = view(L,n+1:n+mm,1:n)
#         end
#
#         # Compute gradients of the dynamics
#         fdx, fdu = res.fdx[k], res.fdu[k]
#
#         # Gradients and Hessians of Taylor Series Expansion of Q
#         Qx = Lx + fdx'*s[k+1]
#         Qu = Lu + fdu'*s[k+1]
#         Qxx = Lxx + fdx'*S[k+1]*fdx
#         Quu = Luu + fdu'*S[k+1]*fdu
#         Qux = Lux + fdu'*S[k+1]*fdx
#
#         # Constraints
#         if res isa ConstrainedIterResults
#             Cx, Cu = res.Cx[k], res.Cu[k]
#             Qx .+= Cx'*(Iμ[k]*C[k] + λ[k])
#             Qu .+= Cu'*(Iμ[k]*C[k] + λ[k])
#             Qxx .+= Cx'*Iμ[k]*Cx
#             Quu .+= Cu'*Iμ[k]*Cu
#             Qux .+= Cu'*Iμ[k]*Cx
#         end
#
#         if solver.opts.bp_reg_type == :state
#             Quu_reg = Quu + res.ρ[1]*fdu'*fdu
#             Qux_reg = Qux + res.ρ[1]*fdu'*fdx
#         elseif solver.opts.bp_reg_type == :control
#             Quu_reg = Quu + res.ρ[1]*I
#             Qux_reg = Qux
#         end
#
#         # Regularization
#         if !isposdef(Hermitian(Array(Quu_reg)))  # need to wrap Array since isposdef doesn't work for static arrays
#             # increase regularization
#             regularization_update!(res,solver,:increase)
#
#             # reset backward pass
#             k = N-1
#             Δv[1] = 0.
#             Δv[2] = 0.
#             continue
#         end
#
#         if k == N-1
#             g = res.Cx_N*fdu
#             h = res.Cx_N*fdx
#
#             tmp = g*(Quu\g')
#             res.Kλ[N] = -tmp\(g*(Quu\Qux) - h)
#             res.Mλ[N] = -tmp\(g*(Quu\res.Cu[k]'))
#             res.dλ[N] = -tmp\(g*(Quu\Qu) - res.CN)
#
#             cx_active = res.Cx[k][res.active_set[k],:]
#             cu_active = res.Cu[k][res.active_set[k],:]
#             c_active = res.C[k][res.active_set[k]]
#
#             tmp = cu_active*(Quu\(cu_active' + (g'*res.Mλ[N])[:,res.active_set[k]]))
#             res.Kλ[k][res.active_set[k],:] = -tmp\(cu_active*(Quu\(Qux + g'*res.Kλ[N])) - cx_active)
#             res.dλ[k][res.active_set[k]] = -tmp\(cu_active*(Quu\(Qu + g'*res.dλ[N])) - c_active)
#
#             res.K[k] = -Quu\(Qux + g'*res.Kλ[N] + (res.Cu[k]' + g'*res.Mλ[N])*res.Kλ[k])
#             res.d[k] = -Quu\(Qu + g'*res.dλ[N] + (res.Cu[k]' + g'*res.Mλ[N])*res.dλ[k])
#         else
#             cx_active = res.Cx[k][res.active_set[k],:]
#             cu_active = res.Cu[k][res.active_set[k],:]
#             c_active = res.C[k][res.active_set[k]]
#
#             tmp = cu_active*(Quu\cu_active')
#             res.Kλ[k][res.active_set[k],:] = -tmp\(cu_active*(Quu\Qux) - cx_active)
#             res.dλ[k][res.active_set[k]] = -tmp\(cu_active*(Quu\Qu) - c_active)
#
#             res.K[k] = -Quu\(Qux + res.Cu[k]'*res.Kλ[k])
#             res.d[k] = -Quu\(Qu + res.Cu[k]'*res.dλ[k])
#         end
#
#         # Compute gains
#         # K[k] = -Quu_reg\Qux_reg
#         # d[k] = -Quu_reg\Qu
#
#         # Calculate cost-to-go (using unregularized Quu and Qux)
#         s[k] = Qx + K[k]'*Quu*d[k] + K[k]'*Qu + Qux'*d[k]
#         S[k] = Qxx + K[k]'*Quu*K[k] + K[k]'*Qux + Qux'*K[k]
#         S[k] = 0.5*(S[k] + S[k]')
#
#         # calculated change is cost-to-go over entire trajectory
#         Δv[1] += d[k]'*Qu
#         Δv[2] += 0.5*d[k]'*Quu*d[k]
#
#         k = k - 1;
#     end
#
#     # decrease regularization after backward pass
#     regularization_update!(res,solver,:decrease)
#
#     return Δv
# end

a = 1
# function _backwardpass_active_set!(res::SolverVectorResults,solver::Solver)
#     n,m,N = get_sizes(solver)
#     m̄,mm = get_num_controls(solver)
#     p,pI,pE = get_num_constraints(solver)
#     pN = n
#
#     W = solver.obj.Q; R = solver.obj.R; Wf = solver.obj.Qf; xf = solver.obj.xf
#
#     solver.state.minimum_time ? R_minimum_time = solver.opts.R_minimum_time : nothing
#     solver.state.infeasible ? R_infeasible = solver.opts.R_infeasible*Matrix(I,n,n) : nothing
#
#     dt = solver.dt
#
#     # pull out values from results
#     X = res.X; U = res.U; K = res.K; d = res.d;
#     L = res.L; l = res.l
#
#     # Useful linear indices
#     idx = Array(1:2*(n+mm))
#     Idx = reshape(Array(1:(n+mm)^2),n+mm,n+mm)
#
#     x_idx = idx[1:n]
#     u_idx = idx[n+1:n+m]
#     ū_idx = idx[n+1:n+mm]
#
#
#     S = [k != N ? zeros(n+p,n+p) : zeros(n+pN,n+pN) for k = 1:N]
#     s = [k != N ? zeros(n+p) : zeros(n+pN) for k = 1:N]
#
#     # Boundary Conditions
#     C = res.C; Iμ = res.Iμ; λ = res.λ; CxN = res.Cx_N
#
#     Sxx = Wf + CxN'*res.IμN*CxN
#     Sxλ = CxN'
#     Sλλ = zeros(pN,pN)
#
#     Sx = Wf*(X[N] - xf) + CxN'*(res.IμN*res.CN + res.λN)
#     Sλ = res.CN
#
#     S[N][1:n,1:n] = Sxx
#     S[N][1:n,n+1:n+pN] = Sxλ
#     S[N][n+1:n+pN,1:n] = Sxλ'
#     S[N][n+1:n+pN,n+1:n+pN] = Sλλ
#
#     s[N][1:n] = Sx
#     s[N][n+1:n+pN] = Sλ
#
#     # Initialize expected change in cost-to-go
#     Δv = zeros(2)
#
#     # Backward pass
#     k = N-1
#     while k >= 1
#         println("k: $k")
#         solver.state.minimum_time ? dt = U[k][m̄]^2 : nothing
#
#         x = X[k]
#         u = U[k][1:m]
#
#         ℓ1 = ℓ(x,u,W,R,xf)
#         ℓx = W*(x - xf)
#         ℓu = R*u
#         ℓxx = W
#         ℓuu = R
#         ℓux = zeros(m,n)
#
#         # Assemble expansion
#
#         Lx = dt*ℓx
#         Lu = dt*ℓu
#         Lxx = dt*ℓxx
#         Luu = dt*ℓuu
#         Lux = dt*ℓux
#
#         # Minimum time expansion components
#         if solver.state.minimum_time
#             h = U[k][m̄]
#
#             l[n+m̄] = 2*h*(ℓ1 + R_minimum_time)
#
#             tmp = 2*h*ℓu
#             L[u_idx,m̄] = tmp
#             L[m̄,u_idx] = tmp'
#             L[m̄,m̄] = 2*(ℓ1 + R_minimum_time)
#
#             L[m̄,x_idx] = 2*h*ℓx'
#         end
#
#         # Infeasible expansion components
#         if solver.state.infeasible
#             l[n+m̄+1:n+mm] = R_infeasible*U[k][m̄+1:m̄+n]
#             L[n+m̄+1:n+mm,n+m̄+1:n+mm] = R_infeasible
#         end
#
#         # Final expansion terms
#         if solver.state.minimum_time || solver.state.infeasible
#             l[u_idx] = Lu
#
#             L[u_idx,u_idx] = Luu
#             L[u_idx,x_idx] = Lux
#
#             Lu = l[n+1:n+mm]
#             Luu = L[n+1:n+mm,n+1:n+mm]
#             Lux = L[n+1:n+mm,1:n]
#         end
#
#         # Constraints
#         Cx, Cu = res.Cx[k], res.Cu[k]
#         Lx += Cx'*(Iμ[k]*C[k] + λ[k])
#         Lu += Cu'*(Iμ[k]*C[k] + λ[k])
#         Lxx += Cx'*Iμ[k]*Cx
#         Luu += Cu'*Iμ[k]*Cu
#         Lux += Cu'*Iμ[k]*Cx
#
#         # multiplier expansion terms
#         Lλ = res.C[k]
#         Lxλ = res.Cx[k]'
#         Lλλ = zeros(p,p)
#         Lλu = res.Cu[k]
#
#         # Compute gradients of the dynamics
#         fdx, fdu = res.fdx[k], res.fdu[k]
#
#         # Gradients and Hessians of Taylor Series Expansion of Q
#         Qx = Lx + fdx'*s[k+1][1:n]
#         Qu = Lu + fdu'*s[k+1][1:n]
#         Qxx = Lxx + fdx'*S[k+1][1:n,1:n]*fdx
#         Quu = Luu + fdu'*S[k+1][1:n,1:n]*fdu
#         Qux = Lux + fdu'*S[k+1][1:n,1:n]*fdx
#
#         # regularization
#         Quu_reg = Quu + res.ρ[1]*I
#         Qux_reg = Qux
#         Luu_reg = Luu + res.ρ[1]*I
#         Lux_reg = Lux
#
#         # Regularization
#         if !isposdef(Hermitian(Array(Quu_reg)))  # need to wrap Array since isposdef doesn't work for static arrays
#             # increase regularization
#             regularization_update!(res,solver,:increase)
#
#             Sxx = S[N][1:n,1:n]
#             Sxλ = S[N][1:n,n+1:n+pN]
#             Sλλ = S[N][n+1:n+pN,n+1:n+pN]
#
#             Sx = s[N][1:n]
#             Sλ = s[N][n+1:n+pN]
#
#             # reset backward pass
#             k = N-1
#             Δv[1] = 0.
#             Δv[2] = 0.
#             continue
#         end
#
#         # linearized dynamics for cost-to-go
#         S̄xx = fdx'*Sxx*fdx
#         S̄xλ = fdx'*Sxλ
#         S̄xu = fdx'*Sxx*fdu
#         S̄λλ = Sλλ
#         S̄λu = Sxλ'*fdu
#         S̄uu = fdu'*Sxx*fdu
#
#         S̄x = fdx'*Sx
#         S̄λ = Sλ
#         S̄u = fdu'*Sx
#
#         # get active set constraints and Jacobians
#         active_set_idx_1 = res.active_set[k] # active set indices at index k
#         k != N-1 ? active_set_idx_2 = res.active_set[k+1] : active_set_idx_2 = ones(Bool,pN) # active set indices at index k+1
#
#         c_active_set_1 = res.C[k][active_set_idx_1] # active set constraints at index k
#         c_active_set_2 = S̄λ[active_set_idx_2] # active set constraints at index k+1
#         c_set_1 = res.C[k]
#         c_set_2 = S̄λ
#
#         println("S̄λ: $S̄λ")
#         k != N-1 ? println("C: $(res.C[k+1])") : println("C: $(res.CN)")
#         println("active set: $(active_set_idx_2)")
#
#
#         cx_active_set_1 = res.Cx[k][active_set_idx_1,:] # active set Jacobians at index k
#         cx_active_set_2 = (S̄xλ[:,active_set_idx_2])' # active set Jacobians at index k+1
#
#         cx_set_1 = res.Cx[k]
#         cx_set_2 = S̄xλ'
#
#         cu_active_set_1 = res.Cu[k][active_set_idx_1,:] # active set Jacobians at index k
#         cu_active_set_2 = S̄λu[active_set_idx_2,:] # active set Jacobians at index k+1
#         cu_set_1 = res.Cu[k]
#         cu_set_2 = S̄λu
#
#         if sum(active_set_idx_2) != 0
#             # calculate multiplier gains
#             tmp = cu_active_set_2*(Quu_reg\cu_active_set_2')
#
#             res.Kλ[k+1][active_set_idx_2,:] = -tmp\(cu_active_set_2*(Quu_reg\Qux_reg) - cx_active_set_2)
#             res.Mλ[k+1][active_set_idx_2,active_set_idx_1] = -tmp\(cu_active_set_2*(Quu_reg\cu_active_set_1'))
#             res.dλ[k+1][active_set_idx_2] = -tmp\(cu_active_set_2*(Quu_reg\Qu) - c_active_set_2)
#
#             # reset inactive gains to zero (this may be unnecessary)
#             res.Kλ[k+1][.!active_set_idx_2,:] .= 0.
#             res.Mλ[k+1][.!active_set_idx_2,.!active_set_idx_1] .= 0.
#             res.dλ[k+1][.!active_set_idx_2] .= 0.
#         end
#
#         if k != 1
#             # calculate control gains
#             res.K[k] = -Quu_reg\(Qux_reg + cu_set_2'*res.Kλ[k+1])
#             res.M[k] = -Quu_reg\(cu_set_1' + cu_set_2'*res.Mλ[k+1])
#             res.d[k] = -Quu_reg\(Qu + cu_set_2'*res.dλ[k+1])
#         elseif k == 1
#             if sum(res.active_set[1]) != 0
#                 # calculate final multiplier gains
#                 tmp = cu_active_set_1*(Quu_reg\(cu_active_set_1' + cu_active_set_2'*res.Mλ[k+1][active_set_idx_2,active_set_idx_1]))
#                 res.Kλ[1][active_set_idx_1,:] = -tmp\(cu_active_set_1*(Quu_reg\(Qux_reg + cu_active_set_2'*res.Kλ[k+1][active_set_idx_2,:])) - cx_active_set_1)
#                 res.dλ[1][active_set_idx_1] = -tmp\(cu_active_set_1*(Quu_reg\(Qu + cu_active_set_2'*res.dλ[k+1][active_set_idx_2])) - c_active_set_1)
#
#                 # reset inactive gains to zero (this may be unnecessary)
#                 res.Kλ[1][.!active_set_idx_1,:] .= 0.
#                 res.dλ[1][.!active_set_idx_1] .= 0.
#             end
#
#             # calculate final control gains
#             res.K[1] = -Quu_reg\(Qux_reg + cu_set_1'*res.Kλ[1] + cu_set_2'*res.Kλ[k+1] + cu_set_2'*res.Mλ[k+1]*res.Kλ[1])
#             res.d[1] = -Quu_reg\(Qu + cu_set_1'*res.dλ[1] + cu_set_2'*res.dλ[k+1] + cu_set_2'*res.Mλ[k+1]*res.dλ[1])
#         end
#
#         # plug in multiplier equation
#         Ŝxx = S̄xx + S̄xλ*res.Kλ[k+1] + res.Kλ[k+1]'*S̄xλ' + res.Kλ[k+1]'*S̄λλ*res.Kλ[k+1]
#         Ŝxλ = S̄xλ*res.Mλ[k+1] + res.Kλ[k+1]'*S̄λλ*res.Mλ[k+1]
#         Ŝxu = res.Kλ[k+1]'*S̄λu + S̄xu
#         Ŝλλ = res.Mλ[k+1]'*S̄λλ*res.Mλ[k+1]
#         Ŝλu = res.Mλ[k+1]'*S̄λu
#         Ŝuu = S̄uu
#
#         Ŝx = S̄xλ*res.dλ[k+1] + res.Kλ[k+1]'*S̄λλ*res.dλ[k+1] + S̄x + res.Kλ[k+1]'*S̄λ
#         Ŝλ = res.Mλ[k+1]'*S̄λλ*res.dλ[k+1] + res.Mλ[k+1]'*S̄λ
#         Ŝu = S̄λu'*res.dλ[k+1] + S̄u
#
#         Δv[1] += res.dλ[k+1]'*S̄λ
#         Δv[2] += 0.5*res.dλ[k+1]'*S̄λλ*res.dλ[k+1]
#         println("Δv: $Δv")
#
#         # calculate cost-to-go
#         Q̂xx = Lxx + Ŝxx
#         Q̂xλ = Lxλ + Ŝxλ
#         Q̂ux = Lux_reg + Ŝxu'
#         Q̂uu = Luu_reg + Ŝuu
#         Q̂λu = Lλu + Ŝλu
#         Q̂λλ = Lλλ + Ŝλλ
#         Q̂x = Lx + Ŝx
#         Q̂u = Lu + Ŝu
#         Q̂λ = Lλ + Ŝλ
#
#         Sxx = Q̂xx + Q̂ux'*res.K[k] + res.K[k]'*Q̂ux + res.K[k]'*Q̂uu*res.K[k]
#         Sxλ = Q̂xλ + Q̂ux'*res.M[k] + res.K[k]'*Q̂λu' + res.K[k]'*Q̂uu*res.M[k]
#         Sλλ = Q̂λλ + res.M[k]'*Q̂λu' + res.M[k]'*Q̂uu*res.M[k]
#
#         Sx = Q̂x + Q̂ux'*res.d[k] + res.K[k]'*Q̂uu*res.d[k] + res.K[k]'*Q̂u
#         Sλ = Q̂λ + Q̂λu*res.d[k] + res.M[k]'*Q̂uu*res.d[k] + res.M[k]'*Q̂u
#
#         S[k][1:n,1:n] = Sxx
#         S[k][1:n,n+1:n+p] = Sxλ
#         S[k][n+1:n+p,1:n] = Sxλ'
#         S[k][n+1:n+p,n+1:n+p] = Sλλ
#
#         s[k][1:n] = Sx
#         s[k][n+1:n+p] = Sλ
#
#         Δv[1] += res.d[k]'*Q̂u
#         Δv[2] += 0.5*res.d[k]'*Q̂uu*res.d[k]
#
#         if k == 1
#             Δv[1] += res.dλ[1]'*Sλ
#             Δv[2] += 0.5*res.dλ[1]'*Sλλ*res.dλ[1]
#         end
#         println("Δv: $Δv")
#         k = k - 1;
#     end
#
#     # decrease regularization after backward pass
#     regularization_update!(res,solver,:decrease)
#
#     return Δv
# end
#
# function rollout!(res::SolverVectorResults,solver::Solver,alpha::Float64)
#     n,m,N = get_sizes(solver)
#     m̄,mm = get_num_controls(solver)
#     dt = solver.dt
#
#     X = res.X; U = res.U; K = res.K; d = res.d; X_ = res.X_; U_ = res.U_
#
#     X_[1] = solver.obj.x0;
#
#     if solver.control_integration == :foh
#         b = res.b
#         du = alpha*d[1]
#         U_[1] = U[1] + du
#         dv = zero(du)
#     end
#
#     if solver.opts.active_set_flag
#         λ = res.λ
#         Kλ = res.Kλ
#         Mλ = res.Mλ
#         dλ = res.dλ
#
#         M = res.M
#
#         δλ = alpha*dλ[1]
#         λ[1] += δλ
#     end
#
#     for k = 2:N
#         δx = X_[k-1] - X[k-1]
#
#         if solver.control_integration == :foh
#             dv = K[k]*δx + b[k]*du + alpha*d[k]
#             U_[k] = U[k] + dv
#             solver.state.minimum_time ? dt = U_[k-1][m̄]^2 : nothing
#             solver.fd(X_[k], X_[k-1], U_[k-1][1:m], U_[k][1:m], dt)
#             du = dv
#         elseif solver.opts.active_set_flag
#             U_[k-1] = U[k-1] + K[k-1]*δx + M[k-1]*δλ + alpha*d[k-1]
#             δλ = Kλ[k]*δx + Mλ[k]*δλ + alpha*dλ[k]
#             k != N ? λ[k] += δλ : res.λN += δλ
#             solver.state.minimum_time ? dt = U_[k-1][m̄]^2 : nothing
#             solver.fd(X_[k], X_[k-1], U_[k-1][1:m], dt)
#         else
#             U_[k-1] = U[k-1] + K[k-1]*δx + alpha*d[k-1]
#             solver.state.minimum_time ? dt = U_[k-1][m̄]^2 : nothing
#             solver.fd(X_[k], X_[k-1], U_[k-1][1:m], dt)
#         end
#
#         solver.state.infeasible ? X_[k] += U_[k-1][m̄.+(1:n)] : nothing
#
#         # Check that rollout has not diverged
#         if ~(norm(X_[k],Inf) < solver.opts.max_state_value && norm(U_[k-1],Inf) < solver.opts.max_control_value)
#             return false
#         end
#     end
#
#     # Calculate state derivatives and midpoints
#     if solver.control_integration == :foh
#         calculate_derivatives!(res, solver, X_, U_)
#         calculate_midpoints!(res, solver, X_, U_)
#     end
#
#     # Update constraints
#     update_constraints!(res,solver,X_,U_)
#
#     return true
# end

# function rollout_active_set!(res::SolverVectorResults,solver::Solver,alpha::Float64)
#     n,m,N = get_sizes(solver)
#     m̄,mm = get_num_controls(solver)
#     dt = solver.dt
#
#     X = res.X; U = res.U; K = res.K; d = res.d; X_ = res.X_; U_ = res.U_
#
#     X_[1] = solver.obj.x0;
#
#     if solver.control_integration == :foh
#         b = res.b
#         du = alpha*d[1]
#         U_[1] = U[1] + du
#         dv = zero(du)
#     end
#
#     if solver.opts.active_set_flag
#         λ = res.λ
#         Kλ = res.Kλ
#         Mλ = res.Mλ
#         dλ = res.dλ
#
#         δλ = alpha*dλ[1]
#         λ[1] += δλ
#     end
#
#     for k = 2:N
#         δx = X_[k-1] - X[k-1]
#
#         if solver.control_integration == :foh
#             dv = K[k]*δx + b[k]*du + alpha*d[k]
#             U_[k] = U[k] + dv
#             solver.state.minimum_time ? dt = U_[k-1][m̄]^2 : nothing
#             solver.fd(X_[k], X_[k-1], U_[k-1][1:m], U_[k][1:m], dt)
#             du = dv
#         else
#             U_[k-1] = U[k-1] + K[k-1]*δx + alpha*d[k-1]
#             solver.state.minimum_time ? dt = U_[k-1][m̄]^2 : nothing
#             solver.fd(X_[k], X_[k-1], U_[k-1][1:m], dt)
#         end
#
#         if solver.opts.active_set_flag
#             k != N ? δλ = Kλ[k]*δx + alpha*dλ[k] : nothing
#             k != N ? λ[k] += δλ : res.λN .+= Kλ[k]*δx + Mλ[k]*δλ + alpha*dλ[k]
#         end
#         solver.state.infeasible ? X_[k] += U_[k-1][m̄.+(1:n)] : nothing
#
#         # Check that rollout has not diverged
#         if ~(norm(X_[k],Inf) < solver.opts.max_state_value && norm(U_[k-1],Inf) < solver.opts.max_control_value)
#             return false
#         end
#     end
#
#     # Calculate state derivatives and midpoints
#     if solver.control_integration == :foh
#         calculate_derivatives!(res, solver, X_, U_)
#         calculate_midpoints!(res, solver, X_, U_)
#     end
#
#     # Update constraints
#     update_constraints!(res,solver,X_,U_)
#
#     return true
# end


text:
module TimeSeriesController
using DataFrames, Query, Dates
using ..Stats, ..Histogram, ..FFT, ..Correlogram
using JSON

export call
min_time(times, time_category) = findmin(times)
max_time(times, time_category) = findmax(times)

function group_dataframe(df, grouping_type)
  if grouping_type == "min"
    return df |>
    @map({time=_.time, value=Stats.min(convert(Array{Float64,1}, _.value)) }) |>
    DataFrame
  elseif grouping_type == "max"
    return df |>
    @map({time=_.time, value=Stats.max(convert(Array{Float64,1}, _.value)) }) |>
    DataFrame
  elseif grouping_type == "average"
    return df |>
    @map({time=_.time, value=Stats.mean(convert(Array{Float64,1}, _.value)) }) |>
    @mutate(value=Float64(_.value)) |>
    DataFrame
  elseif grouping_type == "sum"
    return df |>
    @map({time=_.time, value=Stats.sum(convert(Array{Float64,1}, _.value)) }) |>
    DataFrame
  else
    return df |>
    @map({time=_.time, value=Stats.count(convert(Array{Float64,1}, _.value)) }) |>
    DataFrame
  end
end

function format_date_time(time, time_category)
  if time_category == "date"
    return Dates.format(time, dateformat"yyyy-mm-dd")
  elseif time_category == "year"
    return Dates.format(time, dateformat"yyyy")
  elseif time_category == "year_month"
    return Dates.format(time, dateformat"yyyy-mm")
  elseif time_category == "date_hour"
    return Dates.format(time, dateformat"yyyy-mm-dd HH")
  elseif time_category == "date_hour_minute"
    return Dates.format(time, dateformat"yyyy-mm-dd HH:MM")
  else
    return Dates.format(time, dateformat"yyyy-mm-dd HH:MM:SS")
  end
end

function create_records_with_time(records, time_category, grouping_type, missing_complement)
  new_records = hcat(records...)
  ts = DateTime.(new_records[1, :], dateformat"y-m-d H:M:S")
  times = format_date_time.(ts, time_category)
  min, min_index = min_time(ts, time_category)
  max, max_index = max_time(ts, time_category)

  # grouping
  df = DataFrame(time=times, value=convert(Array{Float64,1}, new_records[2, :]))
  grouping_df = df |>
  @groupby(_.time) |>
  @map({time=getindex(_, [1, 1])[1].time, value=convert(Array{Float64,1}, get.(_, :value, [0]))})

  new_records = Matrix{Union{Float64, String}}(group_dataframe(grouping_df, grouping_type))

  # missing value complement
  existed_times = new_records[:, 1]
  if time_category == "date"
    interval = Day(1)
  elseif time_category == "year"
    interval = Year(1)
  elseif time_category == "year_month"
    interval = Month(1)
  elseif time_category == "date_hour"
    interval = Hour(1)
  elseif time_category == "date_hour_minute"
    interval = Minute(1)
  else
    interval = Second(1)
  end
  for time = format_date_time.(min:interval:max, time_category)
    if !(time in existed_times)
      new_records = vcat(new_records, [time 0])
    end
  end

  return sortslices(rotr90(new_records), dims=2)
end

function create_records_without_time(records)
  new_records = hcat(records...)
  permutedims(hcat(new_records[1, :], convert(Array{Float64,1}, new_records[2, :])))
end

function create_records(records, time_category, grouping_type, missing_complement)
  if time_category != "number"
    create_records_with_time(records, time_category, grouping_type, missing_complement)
  else
    create_records_without_time(records)
  end
end

function hist(values, time_category, grouping_type)
  values = convert(Array{Float64, 1}, values)
  Histogram.fit(values, 8)
end

function correlogram(values)
  values = convert(Array{Float64,1}, values)
  Correlogram.call(values)
end

function fft(values, time_category)
  arr = convert(Array{Float64, 1}, values)
  FFT.fft(arr)
end

"""
call(message)

Computes histgram, correlogram, fft for time series.

# Arguments
- `message::Dict`: the arguments dict
  - `"records"`: the array of array of rows
  - `"time_category"`: time category (date, year, year_month, date_hour, date_hour_minute, secound, number)
  - `"grouping_type"`: grouping policy for duplicated time (min, max, average, sum, count)

# Return
the json string of result object

# Examples
```jldoctest
julia> TimeSeriesController.call(Dict("records" => [["2021-04-01", 2.0], ["2021-04-02", 1.0], ["2021-04-03", 0.0]], "time_category" => "date", "grouping_type" => "average"))
"{\"correlogram\":[[0.0,1.0],[1.0,0.0],[2.0,-0.5]],\"fft\":[[\"1/6\",1.0],[\"1/2\",0.5773502691896257]],\"hist\":[[0.0,1.0],[0.5,0.0],[1.0,1.0],[1.5,0.0],[2.0,1.0]]}"

```
"""
function call(message)
  _records = message["records"]
  time_category = message["time_category"]
  grouping_type = message["grouping_type"]
  missing_complement = "zero"

  records = create_records(_records, time_category, grouping_type, missing_complement)
  values = records[2, :]

  dict = Dict(:hist => hist(values, time_category, grouping_type),
              :correlogram => correlogram(values),
              :fft => fft(values, time_category)) |> JSON.json
end
end

ext: jl


text:
#        field symbol   name        is array
fields = [ (:X,       "positions",    true),
           (:P,       "momenta",      true),
           (:M,       "masses",       true),
           (:Z,       "numbers",      true),
           (:cell,    "cell",         false),
           (:calc,    "calculator",   false),
           (:cons,    "constraint",   false),
           (:pbc,     "pbc",          false),
           ]

for (S, name, isarray) in fields
   set_name = Meta.parse("set_$(name)!")
   get_name = Meta.parse("$name")
   if isarray
      @eval begin
         function $(get_name)(at::Atoms)
            return copy(at.$S)
         end
         function $(set_name)(at::Atoms, Q)
            if length(at.$S) != length(Q)
               at.$S = copy(Q)
            else
               at.$S .= Q
            end
            return at
         end
      end
   else
      function $(get_name)(at::Atoms)
         return at.$S
      end
      function $(set_name)(at::Atoms, Q)
         at.$S = Q
         return at
      end
   end
end


text:
abstract type AbstractDiffusion <: AbstractPhysicsComponent end

Base.@kwdef struct ConstantViscosity{FT} <: AbstractDiffusion
    μ::FT
    ν::FT
    κ::FT
end

@inline function calc_diffusive_flux_argument!(grad, ::Nothing, _...)
    grad.∇ρ = 0
    grad.∇u = @SVector [0, 0, 0]
    grad.∇θ = 0

    return nothing
end

@inline function calc_diffusive_flux_argument!(grad, diff::ConstantViscosity, state::Vars, aux::Vars, t::Real)
    ρ = state.ρ
    ρu = state.ρu
    ρθ = state.ρθ

    u = ρu / ρ
    θ = ρθ / ρ

    grad.∇ρ = ρ
    grad.∇u = u
    grad.∇θ = θ

    return nothing
end

@inline function calc_diffusive_flux!(gradflux, ::Nothing, _...)
    gradflux.μ∇ρ = @SVector [0, 0, 0]
    gradflux.ν∇u = @SMatrix zeros(3, 3)
    gradflux.κ∇θ = @SVector [0, 0, 0]

    return nothing
end

@inline function calc_diffusive_flux!(gradflux, diff::ConstantViscosity, grad::Grad, state::Vars, aux::Vars, t::Real)
    μ = diff.μ * I
    ν = diff.ν * I
    κ = diff.κ * I

    gradflux.μ∇ρ = -μ * grad.∇ρ
    gradflux.ν∇u = -ν * grad.∇u
    gradflux.κ∇θ = -κ * grad.∇θ

    return nothing
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
struct GTIData <: JAXTAMData
    mission        :: Mission
    instrument     :: Symbol
    obsid          :: String
    e_range        :: Tuple{Float64,Float64}
    bin_time       :: Real
    gti_index      :: Int
    gti_start_time :: Real
    counts         :: Array
    times          :: StepRangeLen
    group          :: Int
end

"""
    _lc_filter_gtis(binned_times::StepRangeLen, binned_counts::Array{Int,1}, gtis::DataFrames.DataFrame, mission::Symbol, instrument::Symbol, obsid::String; min_gti_sec=16)

Splits the lightcurve (count) data into GTIs

First, removes GTIs under `min_gti_sec`, then puts the lightcurve data
into a `Dict{Int,GTIData}`, with the key as the `index` of the GTI
"""
function _lc_filter_gtis(mission::Mission, instrument::Symbol, obsid::String, e_range::Tuple{Float64,Float64},
        binned_times::StepRangeLen, binned_counts::Array{Int,1}, gtis::DataFrames.DataFrame;
        min_gti_sec=16
    )
    gti_data = Dict{Int,GTIData}()

    gti_groups = gtis[:group]
    gti_times  = hcat(gtis[:start], gtis[:stop])

    bin_time = binned_times[2] - binned_times[1]

    gti_mask = (gtis[:stop] .- gtis[:start]) .>= min_gti_sec # Exclude GTIs under `min_gti_sec`
    gti_times[gti_mask.==false, :] .= -1 # Set excluded GTIs to -1, GTIs aren't just discarded so that their absolute index can be kept track of

    excluded_gti_count = count(gti_mask.==false)

    @info "               -> prelim. excluded $excluded_gti_count GTIs under $(min_gti_sec)s"

    if excluded_gti_count == length(gti_groups)
        # All GTIs have been excluded
        throw(JAXTAMError("All GTIs under $(min_gti_sec)s, none left", :gtis, BoundsError(Any[], (1,))))
    end

    # Do this as the first GTI being zero screws with the initial start index
    # in the upcoming `for gtis` loop
    if gti_times[1, 1] == 0
        gti_times[1, 1] = eps()
    end

    # Dodgy way to convert a matrix into an array of arrays
    # so each GTI is stored as an array of [start; finish]
    # and each of those GTI arrays is an array itself
    # makes life a bit easier for the following `for gti in gtis` loop
    gti_times = [gti_times[x, :] for x in 1:size(gti_times, 1)]

    @info "               -> sorting GTIs"

    # Offset if binned_times are a subarray of the lightcurve
    if typeof(binned_times) <: SubArray
        subarray_offset = binned_times.offset1
    else
        subarray_offset = 0
    end

    for (i, gti) in enumerate(gti_times) # For each GTI, store the selected times and count rate within that GTI
        if gti[1] == -1 # Short GTIs have their start time set to -1, skipped here
            continue
        end

        start = ceil(Int, gti[1]/bin_time)+1  - subarray_offset
        stop  = floor(Int, gti[2]/bin_time)-1 - subarray_offset

        # Subtract GTI start time from all times, so all start from t=0
        array_times = binned_times[start:stop].-gti[1]
        range_times = array_times[1]:bin_time:array_times[end]

        gti_data[Int(i)] = GTIData(mission, instrument, obsid, e_range, bin_time, i, start,
            binned_counts[start:stop], range_times, gti_groups[i])
    end

    total_counts = sum(binned_counts)[1]
    gti_counts   = sum([sum(gti.counts) for gti in values(gti_data)])
    count_delta  = gti_counts-total_counts
    delta_prcnt  = round(count_delta/total_counts*100, digits=2)

    @info "               -> original counts: $total_counts | remaining: $gti_counts | delta: $count_delta ($delta_prcnt %)"

    if delta_prcnt > 10
        @warn "Count delta > 10% of total counts"
    end

    return gti_data
end

"""
    _gtis(lc::BinnedData)

Calls `_lc_filter_gtis` using `BinnedData` input
"""
function _gtis(lc::BinnedData)
    gti_data = _lc_filter_gtis(lc.mission, lc.instrument, lc.obsid, lc.e_range, lc.times, lc.counts, lc.gtis)

    return gti_data
end

"""
    _gtis_save(gtis, gti_dir::String)

Splits up the GTI data into a `_meta.feather` file containing non-array variables for each GTI (index, start/stop times, etc...)
and multiple `_gti.feather` files for each GTI containing the counts and times
"""
function _gtis_save(mission::Mission, obs_row::DataFrames.DataFrameRow{DataFrames.DataFrame},
        instrument::Union{String,Symbol}, e_range::Tuple,
        gtis_dir::String, gtis::Dict{Int64,JAXTAM.GTIData}; log=true
    )
    @info "Saving $instrument gti data"
    mkpath(gtis_dir)
    log_entries = Dict{Union{Symbol,Int},String}()

    gti_indecies = [k for k in keys(gtis)]
    gti_starts   = [t.gti_start_time for t in values(gtis)]
    gti_example  = gtis[gti_indecies[1]]
    gti_groups   = [o.group for o in values(gtis)]

    gti_basename = string("$(gti_example.instrument)_lc_$(gti_example.bin_time)_gti")
    gtis_meta    = DataFrame(mission=string(_mission_name(gti_example.mission)), instrument=String(gti_example.instrument),
        obsid=gti_example.obsid, bin_time=gti_example.bin_time, indecies=gti_indecies, starts=gti_starts, groups=gti_groups)

    Feather.write(joinpath(gtis_dir, "$(gti_basename)_meta.feather"), gtis_meta)
    log_entries[:path_meta] = joinpath(gtis_dir, "$(gti_basename)_meta.feather")

    for index in gti_indecies
        gti_savepath = joinpath(gtis_dir, "$(gti_basename)_$(index).feather")
        Feather.write(gti_savepath, DataFrame(counts=gtis[index].counts, times=gtis[index].times))
        log_entries[index] = gti_savepath
    end

    if log
        _log_add(mission, obs_row,
            Dict("data" =>
                Dict(:gtis =>
                    Dict(e_range =>
                        Dict(gti_example.bin_time =>
                            Dict(instrument =>
                                log_entries
                            )
                        )
                    )
                )
            )
        )
    end
    @info "Saving $instrument gti data finished"
end

"""
    _gtis_load(gti_dir, instrument, bin_time)

Loads and parses the `_meta` and `_gti` files, puts into a `GTIData` constructor, returns `Dict{Int,GTIData}`
"""
function _gtis_load(gti_dir::String, instrument, e_range::Tuple{Float64,Float64}, bin_time)
    @info "Loading $instrument gti data"
    bin_time = float(bin_time)

    gti_basename  = string("$(instrument)_lc_$(bin_time)_gti")
    gti_meta_path = joinpath(gti_dir, "$(gti_basename)_meta.feather")

    gti_meta = Feather.read(gti_meta_path)

    gti_data = Dict{Int,GTIData}()

    for gti_row in DataFrames.eachrow(gti_meta)
        gti_mission = _mission_symbol_to_type(Symbol(gti_row[:mission]))
        gti_inst    = Symbol(gti_row[:instrument])
        gti_obsid   = gti_row[:obsid]
        gti_bin_t   = gti_row[:bin_time]
        gti_idx     = gti_row[:indecies]
        gti_starts  = gti_row[:starts]
        current_gti = Feather.read(joinpath(gti_dir, "$(gti_basename)_$(gti_idx).feather"))
        gti_counts  = current_gti[:counts]
        gti_times   = current_gti[:times]
        gti_times   = gti_times[1]:gti_bin_t:gti_times[end] # Convert Array to Step Range
        gti_group   = gti_row[:groups]

        gti_data[gti_idx] = GTIData(gti_mission, gti_inst, gti_obsid, e_range, gti_bin_t, gti_idx,
            gti_starts, gti_counts, gti_times, gti_group)
    end

    @info "Loading $instrument gti data finished"

    return gti_data
end

"""
    gtis(mission_name::Symbol, obs_row::DataFrames.DataFrame, bin_time::Number;  =false)

Handles file management, checks to see if GTI files exist already and loads them, if files do not
exist then the `_gits` function is ran, then the data is saved

TODO: FIX E_RANGE PROPAGATION WHEN LCURVE_DATA IS GIVEN
"""
function gtis(mission::Mission, obs_row::DataFrames.DataFrameRow{DataFrames.DataFrame}, bin_time::Number; overwrite=false,
        e_range::Tuple{Float64,Float64}=_mission_good_e_range(mission),
        lcurve_data::Dict{Symbol,BinnedData}=Dict{Symbol,BinnedData}()
    )
    obsid       = obs_row[:obsid]
    instruments = _mission_instruments(mission)
    gti_files   = _log_query(mission, obs_row, "data", :gtis, e_range, bin_time)
    gtis_path   = abspath(_obs_path_local(mission, obs_row; kind=:jaxtam), "JAXTAM", _log_query_path(; category=:data, kind=:gtis, e_range=e_range, bin_time=bin_time))

    gtis = Dict{Symbol,Dict{Int64,JAXTAM.GTIData}}()
    for instrument in instruments
        if ismissing(gti_files) || !haskey(gti_files, instrument) || overwrite
            @info "Missing gti files for $instrument"

            if !haskey(lcurve_data, instrument)
                lcurve_data = lcurve(mission, obs_row, bin_time, e_range=e_range)
            end

            @assert _recursive_first(lcurve_data).e_range == e_range "lcurve energy range does not match input: $(_recursive_first(gtis_data).e_range) != $e_range"

            @info "Selecting GTIs for $instrument"
            gtis_data = _gtis(lcurve_data[Symbol(instrument)])

            _gtis_save(mission, obs_row, instrument, e_range, gtis_path, gtis_data)

            gtis[instrument] = gtis_data
        else
            gtis[instrument] = _gtis_load(gtis_path, instrument, e_range, bin_time)
        end
    end

    if ismissing(_log_query(mission, obs_row, "meta", :countrates, e_range))
        countrates = Dict{Symbol,Float64}()
        for instrument in instruments
            countrates[instrument] = mean([sum(gti[2].counts)/gti[2].times[end] for gti in gtis[instrument]])
        end
        _log_add(mission, obs_row, Dict("meta" => Dict(:countrates => Dict(e_range => countrates))))
    end

    return gtis
end

function gtis(mission::Mission, obsid::String, bin_time::Number; overwrite=false)
    obs_row = master_query(mission, :obsid, obsid)

    return gtis(mission, obs_row, bin_time, overwrite=overwrite)
end


text:
module NestedFixedPoint

using StocasEstimatorsBase, Optim
import StocasEstimatorsBase: fit, inner_iterations, outer_iterations
export fit_nfxp, fit, NFXP

include("types.jl")
include("logit.jl")
include("fit.jl")
end # module


text:
__precompile__()

module IMGedu

using Images, Colors, ImageView

export dilation_,
       erosion_,
       filter_gauss_localsigma,
       filter_gauss,
       filter_mean,
       filter_median,
       filter_sobel,
       grad_morf,
       median_threshold,
       morph_thinning,
       otsu_threshold_

# Functions
include("dilation_.jl")
include("erosion_.jl")
include("filter_gauss_localsigma.jl")
include("filter_gauss.jl")
include("filter_mean.jl")
include("filter_median.jl")
include("filter_sobel.jl")
include("grad_morf.jl")
include("median_threshold.jl")
include("morph_thinning.jl")
include("otsu_threshold_.jl")

end


text:
"""
    bin_bounds(minx::T, binwidth::T, maxx::T) where T<:Number
    ->  NTuple(minx, minx + i∙binwidth, maxx) where i = 1:(no_of_bins + 1)
"""
function bin_bounds(minx::T, binwidth::T, maxx::T) where T<:Number
    nbins = Int(round((maxx - minx) / binwidth))
    bounds = [minx + i * binwidth for i = 0:nbins]
    bounds[end] = maxx
    NTuple{length(bounds), eltype(bounds)}(bounds)
end
"""
    normalize_binwise(no_of_binbounds, binbounds, x)
    --> Float64 ∈ [0, 1], in steps
"""
function normalize_binwise(no_of_binbounds, binbounds, x)::Float64
    if x < first(binbounds)|| x > last(binbounds) || isnan(x) || isinf(x)
        NaN
    else
        bin = min(searchsortedlast(binbounds, x), no_of_binbounds - 1)
        (bin - 1)/ (no_of_binbounds - 2)
    end
end

_nancolor(colorscheme::ColorSchemes.ColorScheme) = _nancolor(colorscheme.colors)
_nancolor(colorscheme) = color_with_lumin(rotate_hue(last(colorscheme), 170°), 3 * lumin(last(colorscheme)))

"""
AbstractColorLegend. For subtypes, define functions:
    constructor including keyword arguments
    (l::legendtype)(invalue) -> color
    show
    bin_bounds
    draw_legend
"""
abstract type AbstractColorLegend end

"""
    show(io::IO, m::MIME"image/svg+xml", l<:AbstractColorLegend)
"""
show

"""
    bin_bounds(l<:AbstractColorLegend )
    ->  NTuple(minx, minx + i∙binwidth, maxx) where i = 1:(no_of_bins + 1)
"""
bin_bounds

"""
    draw_legend(p::Point, l<:AbstractColorLegend; include_nan_color = true, max_vert_height = 7EM)
        -> text values (the function also draws the legend)

Input:
    p                is the upper left position for the legend
    l                is the legend definition, which might include a .name field other than :Legend
    include_nan_color keyword argument. Default is 'true', which draws half-height boxes for values
                     outside l.maxlegend and l.minlegend
    max_vert_height  The legend will be truncated at the middle if the stack of boxes is taller. Typical figure height

The legend is vertically aligned.
"""
draw_legend

# Below: BinLegend. Other legend types in separate files, similar structure.

"""
    BinLegend(;maxlegend = missing,
                    binwidth = missing,
                    noofbins = missing,
                    binbounds = NTuple{0, Missing}()::NTuple,
                    minlegend = 0 * maxlegend,
                    colorscheme = leonardo,
                    nan_color = color_with_lumin(rotate_hue(last(colorscheme.colors), 170°), 3 * lumin(last(colorscheme.colors))),
                    name = :BinLegend::Symbol)

Subtype of AbstractColorLegend, suitable for 2 < bins < 100, but can take more. For scalar quantities or numbers, not vectors.

Call an object of this type to convert a quantity or other number to a color. More than ~6 bins look pretty but is not easy to read.

If called with values outside the defined range, 'nan_color' is returned.

Required input, all keyword arguments:
    maxlegend         top of last color bin - values above will take 'nan_color'
    One of these three:
    1) binwidth       Width of each color bin, e.g. (10N - 0N) / 3 = 3.33N. Use the same number type as in input!
    2) noofbins       Bins are sized equally in the closed interval minlengend, maxlegend.
    3) binbounds      Sorted from small to large. You can get a propsal by calling 'bin_bounds'.

Other input:
    minlegend         Bottom of first bin. Default is zero. Values below are 'nan_color'
    colorscheme,      Example: MechanicalSketch.ColorSchemes.Paired_6 (ref. github ColorSchemes.jl)
    nan_color         Output for input ∉ [minlegend, maxlegend] or ∈{inf, NaN, missing}. Default is nearly white,
                      but a transparent color can be useful
    name              keyword argument, must be a Symbol. Default is :Legend. This is used by 'draw_legend'.
"""
struct BinLegend{T, N} <: AbstractColorLegend
    binbounds::SVector{N, T}
    colorscheme::ColorSchemes.ColorScheme
    nan_color::ColorSchemes.Colorant
    name::Symbol
    function BinLegend(::Val{N}, binbounds::NTuple{N, T}, colorscheme, nan_color, name) where {T, N}
         sabibo = SVector{N, T}(binbounds)
         new{T, N}(sabibo, colorscheme, nan_color, name)
    end
 end
 # Outer constructor (don't use directly)
 function BinLegend(binbounds::NTuple{N, T}, colorscheme, nan_color, name) where {T, N}
    @assert length(binbounds) == N
    @assert isbitstype(T)
    @assert issorted(binbounds)
    @assert typeof(nan_color) == typeof(get(colorscheme, 0.5, :clamp))
    BinLegend(Val{N}(), binbounds, colorscheme, nan_color, name)
 end
 # Main, easygoing constructor
 function BinLegend(;maxlegend = missing,
                    binwidth = missing,
                    noofbins = missing,
                    binbounds = NTuple{0, Missing}()::NTuple,
                    minlegend = 0 * maxlegend,
                    colorscheme = ColSchemeNoMiddle,
                    nan_color = _nancolor(colorscheme),
                    name = :BinLegend::Symbol)
    # Check that one of three is missing
    @assert ismissing(maxlegend) + isempty(binbounds) == 1 "Provide either 'maxlegend' or 'binbounds'"
    @assert ismissing(binwidth) + ismissing(noofbins) + isempty(binbounds) == 2 "Provide one of 'binwidth', 'noofbins' or 'binbounds'"
    @assert typeof(maxlegend) == typeof(minlegend)
    bibo = if isempty(binbounds)
               if ismissing(binwidth)
                   binwidth = (maxlegend - minlegend) / noofbins
                end
                @assert typeof(binwidth) == typeof(minlegend)
                bin_bounds(minlegend, binwidth, maxlegend)
            else
                binbounds
            end
    # Pack colors as a ColorScheme
    colscheme = if colorscheme isa ColorSchemes.ColorScheme
        colorscheme
    else
        ColorSchemes.ColorScheme(colorscheme, "Na", "Unknown")
    end
    BinLegend(bibo, colscheme, nan_color, name)
end
# Innermost call of BinLegend type, don't use directly, trait dispatch below.
function (l::BinLegend{T, N})(::Val{N}, q::T) where {T <: Number, N}
    # scale values in defined legend range to defined steps in interval [0, 1]
    normalizedbin = normalize_binwise(N, bin_bounds(l), q)
    # return color
    if isnan(normalizedbin)
        l.nan_color
    else
        get(l.colorscheme, normalizedbin, :clamp)
    end
end
"""
Call the BinLegend type with a number or quantity, just like you would call a function, return a color
"""
(l::BinLegend{T, N})(q::T) where {T, N} = l(Val{N}(), q)

# Obligatory extended functions for this legend type
"""
    bin_bounds(l::BinLegend)
    ->  NTuple(minx, minx + i∙binwidth, maxx) where i = 1:(no_of_bins + 1)
"""
bin_bounds(l::BinLegend) = l.binbounds

function show(io::IO, mime::MIME"image/svg+xml", l::BinLegend)
    cols = [l.nan_color, l.(l.binbounds)[1:(end - 1)]..., l.nan_color]
    colnoalpha = color_without_alpha.(cols)
    buf = IOBuffer()
    show(buf, mime, colnoalpha)
    svgdocstr = String(take!(buf))
    invalid = """<svg xmlns="http://www.w3.org/2000/svg" version="1.1"""
    valid = "<svg "
    svgelem  = valid * split(svgdocstr, invalid)[2]
    htmlstart = "<!DOCTYPE html>"
    vmin = first(l.binbounds)
    vmax = last(l.binbounds)
    htmlpara = """

$(l.name)

Out-of-bounds , $(vmin) to $(vmax) in $(length(l.binbounds) - 1) bins, Out-of-bounds

        """
    htmlend = ""
    htmldoc = htmlstart * htmlpara * svgelem * htmlend
    write(io, htmldoc)
    flush(io)
end

"""
draw_background(pt_upleft, width, height;
                         background_hue = color_with_lumin(PALETTE[8], 0.1),
                         background_opacity = 0.2)
Background for legend
"""
function draw_background(pt_upleft, width::T, height::T;
                         background_hue = color_with_lumin(PALETTE[8], 0.1),
                         background_opacity = 0.2) where T
    @layer begin
        sethue(background_hue)
        setopacity(background_opacity)
        if T <: Quantity
            rect(pt_upleft, width, height; action = :fill)
        else
            rect(pt_upleft, width, height, :fill)
        end
    end
end
"""
    draw_legend(p::Point, l::BinLegendVector{T, N, M, U}; include_nan_color = true,
        max_vert_height = 10EM,
        background_hue = color_with_lumin(PALETTE[8], 0.1),
        background_opacity = 0.5)
        -> (upper left point, lower right point)

Input:
    p                is the upper left position for the legend
    l                is the legend definition, which might include a .name field other than :Legend
    include_nan_color keyword argument. Default is 'true', which draws half-height boxes for values
                     outside l.maxlegend and l.minlegend
    max_vert_height  The legend will be truncated at the middle if the stack of boxes is taller. Typical figure height

The legend is vertically aligned.
"""
function draw_legend(p::Point, l::BinLegend; include_nan_color = true,
                    max_vert_height = 10EM,
                    background_hue = color_with_lumin(PALETTE[8], 0.1),
                    background_opacity = 0.2)
    # Edge of bin values
    binb_sinking = Vector{Any}(reverse(bin_bounds(l)))
    # In bin values
    vals_sinking = reverse(bin_bounds(l))[2:end]
    # Can we fit the entire legend, or is it too tall?
    max_vert_height_pix = scale_to_pt(max_vert_height)
    estimate_full_height = 0.0
    while true
        estimate_full_height = (length(vals_sinking) + 2 + include_nan_color) * row_height()
        estimate_full_height < max_vert_height_pix && break
        # Truncate at middle
        vals_sinking = [vals_sinking[1:div(end, 2) - 1 ] ; [NaN * vals_sinking[1]] ;vals_sinking[div(end, 2) + 2: end] ]
        binb_sinking = [binb_sinking[1:div(end, 2) - 1 ] ; [".."] ;binb_sinking[div(end, 2) + 2: end] ]
    end

    # Positions
    pos_table = p + Point(0, 1) * row_height()
    box_tl  = pos_table + Point(0.5,  0.65    ) * row_height()
    colwidth = max(5.0EM, column_widths(;l.name => binb_sinking)[1])
    fullwidth = pos_table[1] + colwidth - p[1]
    fullheight = pos_table[2] + estimate_full_height - p[2]
    draw_background(p, fullwidth, fullheight; background_hue, background_opacity)
    # Plot
    ptul, ptbr = _draw_legend_boxes(box_tl, l, include_nan_color, vals_sinking)
    text_table_fixed_columns(pos_table, colwidth ; l.name => binb_sinking)

    p, p + Point(fullwidth, fullheight)
end

"""
    _draw_legend_boxes(pt_tl::Point, legend, include_nan_color, vals_sinking)
    --> pt_topleft, pt_br
"""
function _draw_legend_boxes(pt_tl::Point, legend, include_nan_color, vals_sinking)
    gsave()
    pt_topleft = pt_tl + include_nan_color .* (0.0, 0.5 * row_height())
    nancol = legend.nan_color
    if include_nan_color
        # Draw nan_color above, half height
        box_fill_outline(pt_tl +(0.0, 0.5 * row_height()), nancol;
                          height = 0.5 * row_height())
    end
    for v in vals_sinking
        colo = legend(v)
        pt_tl += (0.0, row_height())
        box_fill_outline(pt_tl, colo)
    end
    if include_nan_color
        # Draw nan_color below, half height
        pt_tl += (0.0, 0.5 * row_height())
        box_fill_outline(pt_tl + (0.0, 0.5 * row_height()), nancol;
                          height = 0.5 * row_height())
    end
    pt_br = pt_tl + (EM, row_height())
    grestore()
    pt_topleft, pt_br
end


text:
# ---------------------------------------------------------------------------- #
#
#   compJacob.jl
#
#   Several functions to compute jacobians of an element for both 2D and 3D
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

"""
    compJacob( master::Master2D, nodes::Matrix{Float64},
               ∂ξ∂x::Matrix{Float64}, jcw::Vector{Float64},
               ∂x∂ξ::Matrix{Float64})

Returns Jacobian of a 2D cell.
"""
function compJacob!( master::Master2D, nodes::Matrix{Float64},
                     ∂ξ∂x::Matrix{Float64}, jcw::Matrix{Float64},
                     ∂x∂ξ::Matrix{Float64} )

  ∂x∂ξ[:,1] = master.∇ϕ[:,:,1]' * nodes[:,1] # ∂x∂ξ
  ∂x∂ξ[:,2] = master.∇ϕ[:,:,2]' * nodes[:,1] # ∂x∂η
  ∂x∂ξ[:,3] = master.∇ϕ[:,:,1]' * nodes[:,2] # ∂y∂ξ
  ∂x∂ξ[:,4] = master.∇ϕ[:,:,2]' * nodes[:,2] # ∂y∂η

  jac = ∂x∂ξ[:,1].*∂x∂ξ[:,4] - ∂x∂ξ[:,2].*∂x∂ξ[:,3]

  for ii in 1:length(jac)
    jcw[ii,ii] = master.gwts[ii] * jac[ii]
  end

  ∂ξ∂x[:,1] =  1./jac .* ∂x∂ξ[:,4]
  ∂ξ∂x[:,2] = -1./jac .* ∂x∂ξ[:,2]
  ∂ξ∂x[:,3] = -1./jac .* ∂x∂ξ[:,3]
  ∂ξ∂x[:,4] =  1./jac .* ∂x∂ξ[:,1]

end

"""
    getderbfel( master::Master2D, ∂ξ∂x::Matrix{Float64} )

Returns ∇ϕ in the global coordinate system for 2D meshes.
"""

function getderbfel( master::Master2D, ∂ξ∂x::Matrix{Float64} )

  ∇ϕc = fill(0.0, size(master.∇ϕ))

  ∇ϕc[:,:,1] = master.∇ϕ[:,:,1] * diagm( ∂ξ∂x[:,1] ) + master.∇ϕ[:,:,2] * diagm( ∂ξ∂x[:,3] )
  ∇ϕc[:,:,2] = master.∇ϕ[:,:,1] * diagm( ∂ξ∂x[:,2] ) + master.∇ϕ[:,:,2] * diagm( ∂ξ∂x[:,4] )

  return ∇ϕc

end

"""
    compJacobFace( mesh::Mesh2D, master::Master2D, el::Int64, face::Int64 )

Returns Jacobian on on the `face` in element `el` for a 2D mesh.
"""
function compJacobFace( mesh::Mesh2D, master::Master2D, el::Int64, face::Int64 )

  indnod = 1
  rotdir = false

  if mesh.t2f[el,face] < 0
    # face DOES NOT follow counter-clockwise rotation
    indnod = 2
    rotdir = true
  end

  nod    = master.perm[:,face,indnod]

  ∂x₁∂ξ₁ = master.∇ϕ1d' * mesh.nodes[nod,1,el]
  ∂x₂∂ξ₁ = master.∇ϕ1d' * mesh.nodes[nod,2,el]

  p1d  = master.ϕ1d'  * mesh.nodes[nod,:,el]

  jac  = sqrt.( ∂x₁∂ξ₁.^2 + ∂x₂∂ξ₁.^2 )
  jcw  = master.gwts1d .* jac

  normal = [ ∂x₂∂ξ₁ -∂x₁∂ξ₁ ] ./ [jac jac]

  if rotdir
    normal *= -1.
  end

  tangent = [ normal[:,2] -normal[:,1] ]

  return (master.ϕ1d, p1d, nod, normal, jcw)

end

"""
    compJacob( master::Master2D, ∇ϕ::Array{Float64,3}, gwts::Vector{Float64},
               nodes::Matrix{Float64} )

Returns Jacobian of a 2D mesh.
"""
function compJacob( master::Master2D, nodes::Matrix{Float64} )

  ∂x∂ξ = master.∇ϕ[:,:,1]' * nodes[:,1]
  ∂x∂η = master.∇ϕ[:,:,2]' * nodes[:,1]
  ∂y∂ξ = master.∇ϕ[:,:,1]' * nodes[:,2]
  ∂y∂η = master.∇ϕ[:,:,2]' * nodes[:,2]

  jac = ∂x∂ξ.*∂y∂η - ∂x∂η.*∂y∂ξ
  jcw = diagm( master.gwts ) * jac

  ∂ξ∂x =  1./jac .* ∂y∂η
  ∂η∂x = -1./jac .* ∂y∂ξ
  ∂ξ∂y = -1./jac .* ∂x∂η
  ∂η∂y =  1./jac .* ∂x∂ξ

  ∂ξ∂x_vec = fill(0.0, size(∂x∂ξ,1), 4)
  ∂ξ∂x_vec[:,1] = ∂ξ∂x
  ∂ξ∂x_vec[:,2] = ∂ξ∂y
  ∂ξ∂x_vec[:,3] = ∂η∂x
  ∂ξ∂x_vec[:,4] = ∂η∂y

  return (jcw, ∂ξ∂x_vec)

end

"""
    compJacob!( mesh::Mesh3D, master::Master3D )

Computes Jacobians for 3D meshes and saves them in `mesh`.
"""
function compJacob!( mesh::Mesh3D, master::Master3D )

  # http://www.csun.edu/~lcaretto/me692/Coordinate%20transformations.pdf

  ∂x∂ξ = fill( 0.0, size(master.∇ϕ,2), size(mesh.nodes,3), 3, 3 )
  ∂ξ∂x2 = fill( 0.0, size(master.∇ϕ,2), size(mesh.nodes,3), 3, 3 )

  ∂x∂ξ[:,:,1,1] = master.∇ϕ[:,:,1]' * mesh.nodes[:,1,:]
  ∂x∂ξ[:,:,1,2] = master.∇ϕ[:,:,2]' * mesh.nodes[:,1,:]
  ∂x∂ξ[:,:,1,3] = master.∇ϕ[:,:,3]' * mesh.nodes[:,1,:]
  ∂x∂ξ[:,:,2,1] = master.∇ϕ[:,:,1]' * mesh.nodes[:,2,:]
  ∂x∂ξ[:,:,2,2] = master.∇ϕ[:,:,2]' * mesh.nodes[:,2,:]
  ∂x∂ξ[:,:,2,3] = master.∇ϕ[:,:,3]' * mesh.nodes[:,2,:]
  ∂x∂ξ[:,:,3,1] = master.∇ϕ[:,:,1]' * mesh.nodes[:,3,:]
  ∂x∂ξ[:,:,3,2] = master.∇ϕ[:,:,2]' * mesh.nodes[:,3,:]
  ∂x∂ξ[:,:,3,3] = master.∇ϕ[:,:,3]' * mesh.nodes[:,3,:]

  jac = ∂x∂ξ[:,:,1,1].*∂x∂ξ[:,:,2,2].*∂x∂ξ[:,:,3,3] + ∂x∂ξ[:,:,2,1].*∂x∂ξ[:,:,3,2].*∂x∂ξ[:,:,1,3] + ∂x∂ξ[:,:,3,1].*∂x∂ξ[:,:,1,2].*∂x∂ξ[:,:,2,3] -
        ∂x∂ξ[:,:,3,1].*∂x∂ξ[:,:,2,2].*∂x∂ξ[:,:,1,3] - ∂x∂ξ[:,:,2,1].*∂x∂ξ[:,:,1,2].*∂x∂ξ[:,:,3,3] - ∂x∂ξ[:,:,1,1].*∂x∂ξ[:,:,3,2].*∂x∂ξ[:,:,2,3]
  jcw = diagm( master.gwts ) * jac

  for ii in 1:3, jj in 1:3
    iip1 = ii + 1
    jjp1 = jj + 1
    iip2 = ii + 2
    jjp2 = jj + 2
    if jjp1 != 3
      jjp1 = rem(jjp1,3)
    end
    if jjp2 != 3
      jjp2 = rem(jjp2,3)
    end
    if iip1 != 3
      iip1 = rem(iip1,3)
    end
    if iip2 != 3
      iip2 = rem(iip2,3)
    end
    ∂ξ∂x2[:,:,ii,jj] = 1./jac .* ( ∂x∂ξ[:,:,jjp1,iip1].*∂x∂ξ[:,:,jjp2,iip2] - ∂x∂ξ[:,:,jjp1,iip2].*∂x∂ξ[:,:,jjp2,iip1] )
  end

  ∂ξ∂x = fill(0.0, size(master.∇ϕ,2), size(mesh.nodes,3), 9)
  ∂ξ∂x[:,:,1] = ∂ξ∂x2[:,:,1,1]#∂ξ₁∂x₁
  ∂ξ∂x[:,:,2] = ∂ξ∂x2[:,:,1,2]#∂ξ₁∂x₂
  ∂ξ∂x[:,:,3] = ∂ξ∂x2[:,:,1,3]#∂ξ₁∂x₃

  ∂ξ∂x[:,:,4] = ∂ξ∂x2[:,:,2,1]#∂ξ₂∂x₁
  ∂ξ∂x[:,:,5] = ∂ξ∂x2[:,:,2,2]#∂ξ₂∂x₂
  ∂ξ∂x[:,:,6] = ∂ξ∂x2[:,:,2,3]#∂ξ₂∂x₃

  ∂ξ∂x[:,:,7] = ∂ξ∂x2[:,:,3,1]#∂ξ₃∂x₁
  ∂ξ∂x[:,:,8] = ∂ξ∂x2[:,:,3,2]#∂ξ₃∂x₂
  ∂ξ∂x[:,:,9] = ∂ξ∂x2[:,:,3,3]#∂ξ₃∂x₃

  mesh.jcw  = jcw
  mesh.∂ξ∂x = ∂ξ∂x

end

"""
    getderbfel( master::Master3D, ∂ξ∂x::Matrix{Float64} )

Returns ∇ϕ in the global coordinate system for 3D meshes.
"""

function getderbfel( master::Master3D, ∂ξ∂x::Matrix{Float64} )

  ∇ϕc = fill(0.0, size(master.∇ϕ))

  ∇ϕc[:,:,1] = master.∇ϕ[:,:,1] * diagm( ∂ξ∂x[:,1] ) + master.∇ϕ[:,:,2] * diagm( ∂ξ∂x[:,4] ) + master.∇ϕ[:,:,3] * diagm( ∂ξ∂x[:,7] )
  ∇ϕc[:,:,2] = master.∇ϕ[:,:,1] * diagm( ∂ξ∂x[:,2] ) + master.∇ϕ[:,:,2] * diagm( ∂ξ∂x[:,5] ) + master.∇ϕ[:,:,3] * diagm( ∂ξ∂x[:,8] )
  ∇ϕc[:,:,3] = master.∇ϕ[:,:,1] * diagm( ∂ξ∂x[:,3] ) + master.∇ϕ[:,:,2] * diagm( ∂ξ∂x[:,6] ) + master.∇ϕ[:,:,3] * diagm( ∂ξ∂x[:,9] )

  return ∇ϕc

end

"""
    compJacobFace( mesh::Mesh3D, master::Master3D, el::Int64, face::Int64 )

Returns Jacobian on the `face` in element `el` for a 3D mesh.
"""
function compJacobFace( mesh::Mesh3D, master::Master3D, el::Int64, face::Int64 )

  nod  = master.perm[ :, face, abs.(mesh.t2f[el,face+4]) ]
  #nod  = master.perm[ :, 1, 1 ]

  p2d  = master.ϕ2D'  * mesh.nodes[nod,:,el]

  ∂x₁∂ξ₁ = master.∇ϕ2D[:,:,1]' * mesh.nodes[nod,1,el]
  ∂x₁∂ξ₂ = master.∇ϕ2D[:,:,2]' * mesh.nodes[nod,1,el]
  ∂x₂∂ξ₁ = master.∇ϕ2D[:,:,1]' * mesh.nodes[nod,2,el]
  ∂x₂∂ξ₂ = master.∇ϕ2D[:,:,2]' * mesh.nodes[nod,2,el]
  ∂x₃∂ξ₁ = master.∇ϕ2D[:,:,1]' * mesh.nodes[nod,3,el]
  ∂x₃∂ξ₂ = master.∇ϕ2D[:,:,2]' * mesh.nodes[nod,3,el]

  # cross product to find normal vector, normal = ∂x∂ξ₁ × ∂x∂ξ₂
  normal = hcat( ( ∂x₂∂ξ₁ .* ∂x₃∂ξ₂ - ∂x₂∂ξ₂ .* ∂x₃∂ξ₁ ),
                -( ∂x₁∂ξ₁ .* ∂x₃∂ξ₂ - ∂x₁∂ξ₂ .* ∂x₃∂ξ₁ ),
                 ( ∂x₁∂ξ₁ .* ∂x₂∂ξ₂ - ∂x₁∂ξ₂ .* ∂x₂∂ξ₁ ) )
  # normalize the normal vector
  jac      = sqrt.( normal[:,1].^2 + normal[:,2].^2 + normal[:,3].^2 )
  normal ./= jac * [1,1,1]'

  if mesh.t2f[el,face+4] < 0
    normal *= -1.0
  end

  jcw = master.gwts2D .* jac

  return (master.ϕ2D, p2d, nod, normal, jcw)

end

"""
    compJacob( master::Master3D, nodes::Matrix{Float64},
               ∂ξ∂x::Matrix{Float64}, jcw::Vector{Float64},
               ∂x∂ξ::Matrix{Float64})

Returns Jacobian of a 3D cell.
"""
function compJacob!( master::Master3D, nodes::Matrix{Float64},
                     ∂ξ∂x::Matrix{Float64}, jcw::Matrix{Float64},
                     ∂x∂ξ::Matrix{Float64} )

  # http://www.csun.edu/~lcaretto/me692/Coordinate%20transformations.pdf

  ∂x∂ξ[:,1] = master.∇ϕ[:,:,1]' * nodes[:,1] # 1 1 - 1
  ∂x∂ξ[:,2] = master.∇ϕ[:,:,2]' * nodes[:,1] # 1 2 - 2
  ∂x∂ξ[:,3] = master.∇ϕ[:,:,3]' * nodes[:,1] # 1 3 - 3
  ∂x∂ξ[:,4] = master.∇ϕ[:,:,1]' * nodes[:,2] # 2 1 - 4
  ∂x∂ξ[:,5] = master.∇ϕ[:,:,2]' * nodes[:,2] # 2 2 - 5
  ∂x∂ξ[:,6] = master.∇ϕ[:,:,3]' * nodes[:,2] # 2 3 - 6
  ∂x∂ξ[:,7] = master.∇ϕ[:,:,1]' * nodes[:,3] # 3 1 - 7
  ∂x∂ξ[:,8] = master.∇ϕ[:,:,2]' * nodes[:,3] # 3 2 - 8
  ∂x∂ξ[:,9] = master.∇ϕ[:,:,3]' * nodes[:,3] # 3 3 - 9

  jac = ∂x∂ξ[:,1].*∂x∂ξ[:,5].*∂x∂ξ[:,9] + ∂x∂ξ[:,4].*∂x∂ξ[:,8].*∂x∂ξ[:,3] + ∂x∂ξ[:,7].*∂x∂ξ[:,2].*∂x∂ξ[:,6] -
        ∂x∂ξ[:,7].*∂x∂ξ[:,5].*∂x∂ξ[:,3] - ∂x∂ξ[:,4].*∂x∂ξ[:,2].*∂x∂ξ[:,9] - ∂x∂ξ[:,1].*∂x∂ξ[:,8].*∂x∂ξ[:,6]

  for ii in 1:length(jac)
    jcw[ii,ii] = master.gwts[ii] * jac[ii]
  end

  for ii in 1:3, jj in 1:3
    kk = 3 * (ii - 1) + jj

    iip1 = ii + 1
    jjp1 = jj + 1

    iip2 = ii + 2
    jjp2 = jj + 2

    if jjp1 != 3
      jjp1 = rem(jjp1,3)
    end
    if jjp2 != 3
      jjp2 = rem(jjp2,3)
    end
    if iip1 != 3
      iip1 = rem(iip1,3)
    end
    if iip2 != 3
      iip2 = rem(iip2,3)
    end

    kk11 = 3 * (jjp1 - 1) + iip1
    kk22 = 3 * (jjp2 - 1) + iip2
    kk12 = 3 * (jjp1 - 1) + iip2
    kk21 = 3 * (jjp2 - 1) + iip1

    ∂ξ∂x[:,kk] = 1./jac .* ( ∂x∂ξ[:,kk11].*∂x∂ξ[:,kk22] - ∂x∂ξ[:,kk12].*∂x∂ξ[:,kk21] )
  end

end


text:
__precompile__()
module Camoforge

using Colors
using Images
using PyCall

export get_image
export most_common_colors

include("image_utils/get_image.jl")
include("image_utils/most_common_colors.jl")

end


text:
############################
### Administration Users ###
############################

immutable CreateUserOption
	userName::String
	email::String
	sourceID::Int64
	fullName::Nullable{String}
	loginName::Nullable{String}
	password::Nullable{String}
	sendNotify::Nullable{Bool}
end

function CreateUserOption(userName::String,email::String; sourceID::Int64 = 1, kwargs...)
	kwargsDict = Dict{Symbol,Any}(k=>v for (k,v) in kwargs)

	args = [userName, email, sourceID]

	for field in fieldnames(CreateUserOption)[4:end]
		push!(args,get(kwargsDict,field, fieldtype(CreateUserOption,field)()))
	end

	CreateUserOption(args...)
end

function adminCreateUser(c::Client,userOpt::CreateUserOption)

end


text:
"""
The concept is that, I want to use my account to pay for something.
I can just pick one for payment. If the balance is less than the payment, it will switch to next one.
If the handler cannot deal with the event or the responsibility, the handler dispatch to the next.
"""

abstract type Account end
struct NullAccount <: Account end
can_pay(acc::Account, amount) = acc.balance >= amount

function pay(acc::Account, amount)
    if can_pay(acc, amount)
        acc.balance -= amount
    else
        pay(acc.successor, amount)
    end
end

mutable struct Bank <: Account
    balance::Float64
    successor::Account
end
Bank(balance) = Bank(balance, NullAccount())

mutable struct Paypal <: Account
    balance::Float64
    successor::Account
end
Paypal(balance) = Paypal(balance, NullAccount())

mutable struct Bitcoin <: Account
    balance::Float64
    successor::Account
end
Bitcoin(balance) = Bitcoin(balance, NullAccount())


text:
#using MyExample
#using Test

#import MyExample
#MyExample.my_f(2,1)

using SafeTestsets
@safetestset "My f tests" begin include("My_f_tests.jl") end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

module CoDa

using Tables
using TableOperations
using Distributions
using StatsBase
using StaticArrays
using UnicodePlots
using LinearAlgebra
using FillArrays
using AxisArrays
using Statistics
using Random
using Printf

import Tables
import ScientificTypes as ST
import TableTransforms as TT
import TableTransforms: Transform, Stateless
import TableTransforms: assertions, isrevertible
import TableTransforms: apply, revert, reapply
import Distances: Metric, result_type
import Base: +, -, *, /, ==
import Base: adjoint, inv
import LinearAlgebra: norm, ⋅
import Random: rand

include("compositions.jl")
include("codaarrays.jl")
include("scitypes.jl")
include("distances.jl")
include("transforms.jl")
include("covariances.jl")
include("matrices.jl")

export
  # compositions
  Composition,
  parts, components,
  distance, norm, ⋅,

  # arrays
  CoDaArray,
  compose,

  # distances
  Aitchison,

  # transforms
  Closure,
  Remainder,
  LogRatio,
  ALR, CLR, ILR,
  alr, alrinv,
  clr, clrinv,
  ilr, ilrinv,

  # covariances
  variation,
  alrcov,
  clrcov,
  lrarray,

  # matrices
  Diagonal, I,
  JMatrix, J,
  FMatrix, F,
  GMatrix, G,
  HMatrix, H
end


text:
# ---
# title: Higher-order functions
# description: a set of helper opeartions that allow applying any function
# ---

# These operations are useful to perform an operation that is not explicitly
# defined in Augmentor.

using Augmentor
using Random
using Statistics: mean

Random.seed!(1337)

DecreaseContrast = MapFun(pixel -> pixel / 2)
IncreaseBrightness = AggregateThenMapFun(img -> mean(img),
                                         (pixel, M) -> pixel + M / 5)

img_in = testpattern(RGB, ratio=0.5)
img_out = augment(img_in, DecreaseContrast |> IncreaseBrightness)

# ## References

#md # ```@docs
#md # MapFun
#md # AggregateThenMapFun
#md # ```


text:
#IDM.jl
#A separate file for all the IDM stuff

#############
##IDM Model##
#############

immutable IDMParam
	a::Float64 #max  comfy acceleration
	b::Float64 #max comfy brake speed
	T::Float64 #desired safety time headway
	v0::Float64 #desired speed
	s0::Float64 #minimum headway (e.g. if x is less than this then you crashed)
	del::Float64 #'accel exponent'
end #IDMParam
==(a::IDMParam,b::IDMParam) = (a.a==b.a) && (a.b==b.b) &&(a.T == b.T)&&(a.v0==b.v0)&&(a.s0==b.s0)&&(a.del==b.del)
Base.hash(a::IDMParam,h::UInt64=zero(UInt64)) = hash(a.a,hash(a.b,hash(a.T,hash(a.v0,hash(a.s0,hash(a.del,h))))))

+(a::IDMParam, b::IDMParam) = IDMParam(a.a+b.a, a.b+b.b, a.T+b.T, a.v0+b.v0, a.s0+b.s0, a.del+b.del)
-(a::IDMParam, b::IDMParam) = a+(-1.0*b)
*(a::Real, p::IDMParam) = IDMParam(a*p.a, a*p.b, a*p.T, a*p.v0, a*p.s0, a*p.del)
# .*(p::IDMParam,a::Float64) = a*p
.*(p::IDMParam, v::Vector{Float64}) = IDMParam(v[1]*p.a, v[2]*p.b, v[3]*p.T, v[4]*p.v0, v[5]*p.s0, v[6]*p.del)
.^(p::IDMParam, n::Real) = IDMParam(p.a^n, p.b^n, p.T^n, p.v0^n, p.s0^n, p.del^n)
sqrt(p::IDMParam) = IDMParam(sqrt(p.a), sqrt(p.b), sqrt(p.T), sqrt(p.v0), sqrt(p.s0), sqrt(p.del))
zero(::Type{IDMParam}) = IDMParam(0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
nan(::Type{IDMParam}) = IDMParam(NaN, NaN, NaN, NaN, NaN, NaN)
abs(p::IDMParam) = IDMParam(abs(p.a), abs(p.b), abs(p.T), abs(p.v0), abs(p.s0), abs(p.del))
max(a::IDMParam, b::IDMParam) = IDMParam(max(a.a,b.a), max(a.b,b.b), max(a.T,b.T), max(a.v0,b.v0), max(a.s0,b.s0), max(a.del,b.del))

function get_idm_s_star(p::IDMParam, v::Float64, dv::Float64)
    return p.s0 + max(0.,v*p.T+v*dv/(2*sqrt(p.a*p.b)))
end

# dv is positive if vehicles are closing
function get_idm_dv(p::IDMParam,dt::Float64,v::Float64,dv::Float64,s::Float64)
	s_ = get_idm_s_star(p, v, dv)
    dvdt = p.a*(1.-(v/p.v0)^p.del - (s_/s)^2)
  if (dvdt) > 1.01*p.a
	  print(p.a)
	  print(" || ")
	  print((v/p.v0)^p.del)
	  print(" || ")
	  println(dvdt)
  end
	#dvdt = min(max(dvdt,-p.b),p.a)
	return (dvdt*dt)::Float64
end #get_idm_dv

# these aren't used anymore, but need to be there for the tests
IDMParam(a::Float64,b::Float64,T::Float64,v0::Float64,s0::Float64;del::Float64=4.) = IDMParam(a,b,T,v0,s0,del)
function build_cautious_idm(v0::Float64,s0::Float64)
	T = 2.
	a = 1.
	b = 1.
	return IDMParam(a,b,T,v0,s0)
end

function build_aggressive_idm(v0::Float64,s0::Float64)
	T = 0.8
	a = 2.
	b = 2.
	return IDMParam(a,b,T,v0,s0)
end

function build_normal_idm(v0::Float64,s0::Float64)
	T = 1.4
	a = 1.5
	b = 1.5
	return IDMParam(a,b,T,v0,s0)
end

#TODO: use Enum or something to avoid misspelling errors
function IDMParam(s::AbstractString,v0::Float64,s0::Float64)
	if lowercase(s) == "cautious"
		return build_cautious_idm(v0,s0)
	elseif lowercase(s) == "normal"
		return build_normal_idm(v0,s0)
	elseif lowercase(s) == "aggressive"
		return build_aggressive_idm(v0,s0)
	else
		error("No such idm phenotype: \"$(s)\"")
	end
end

max_stars_repo_name: zhcao/HighwayExit

max_issues_repo_name: zhcao/HighwayExit

max_forks_repo_name: zhcao/HighwayExit


text:
```Write three functions that compute the sum of the numbers in a list: using a
for-loop, a while-loop and recursion. (Subject to availability of these
constructs in your language of choice.```

list = [1,2,3,4,5,6,7.0]

s = 0.0
for i in x
    global s += i
end
s

s = 0.0
c = 1
while c <= size(x,1)
    global s += x[c]
    global c += 1
end
s

function countdown(n)
    if n ≤ 0
        println("Blastoff!")
    else
        print(n, " ")
        countdown(n-1)
    end
end

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
using Test
using TreeTools

println("##### iterators #####")

begin
	nwk = "(((A1:1,A2:1,B:1,C:1)i_2:0.5,(D1:0,D2:0)i_3:1.5)i_1:1,(E:1.5,(F:1.0,(G:0.,H:0.)i_6:0.5)i_5:1)i_4:1)i_r"
	local nodes = sort(["A1", "A2", "B", "C", "i_2", "D1", "D2", "i_3", "i_1", "E", "F", "G", "H", "i_6", "i_5", "i_4", "i_r"])
	local leaves = sort(["A1", "A2", "B", "C", "D1", "D2", "E", "F", "G", "H"])

	@testset "1" begin
		t = node2tree(TreeTools.parse_newick(nwk))
		@test sort([x.label for x in POT(t)]) == nodes
		@test sort([x.label for x in POTleaves(t)]) == leaves
	end
end

max_stars_repo_name: PierreBarrat/TreeTools

max_issues_repo_name: PierreBarrat/TreeTools

max_forks_repo_name: PierreBarrat/TreeTools


text:
#=
learner.jl:

Author: Peter Wolf (opus111@gmail.com)

Port of the FastAI V2 Learner API to Julia

Methods for handling the training loop

The original source is here

https://github.com/fastai/fastai2/blob/master/fastai2/learner.py

The documentation is copied from here

https://github.com/fastai/fastai2/blob/master/docs/learner.html

=#

"""
Basic class handling tweaks of the training loop by changing a [Learner](@ref) in various events

The training loop is defined in [Learner](@ref) a bit below and consists in a minimal set of instructions: looping through the data we:

compute the output of the model from the input
calculate a loss between this output and the desired target
compute the gradients of this loss with respect to all the model parameters
update the parameters accordingly

Any tweak of this training loop is defined in a Callback to avoid over-complicating the code of the training loop, and to make it easy to mix and match different techniques (since they'll be defined in different callbacks).

A callback can implement the following methods:

before_fit
after_fit
after_cancel_fit

before_epoch
after_epoch
after_cancel_epoch

before_epoch_train
after_epoch_train
after_cancel_epoch_train

before_batch_train
batch_train_loss
after_batch_train
after_cancel_batch_train

before_epoch_validate
after_epoch_validate
after_cancel_epoch_validate

before_batch_validate
after_batch_validate
batch_validate_loss
after_cancel_batch_validate

By default handling of these events do nothing.  Special behavior is implemented by overriding these methods

"""
abstract type AbstractCallback end

struct CancelFitException <: Exception end
struct CancelEpochTrainException <: Exception end
struct CancelBatchTrainException <: Exception end
struct CancelEpochValidateException <: Exception end
struct CancelBatchValidateException <: Exception end

"""
An `AbstractLearner` groups together a model, train and validate data,
  optimizer, loss function and callbacks.

[Callbacks](@ref) are used for every tweak of the training loop.
Callbacks receive epoch, batch and loss information which they may pass on to [Metrics](@ref).

# Require interface
- `data_bunch(l::AbstractLearner)`
- `data_bunch!(l::AbstractLearner, data_bunch)`
- `model(l::AbstractLearner)`
- `model!(l::AbstractLearner, model)`
- `loss(l::AbstractLearner)`
- `loss!(l::AbstractLearner, loss)`
- `opt(l::AbstractLearner)`
- `opt!(l::AbstractLearner, opt)`
- `add_cb!(learner::AbstractLearner, cb::AbstractCallback)`
- `cbs(learner::AbstractLearner)`
- `fit!(learner::AbstractLearner, epoch_count)`
"""
abstract type AbstractLearner end

"""
    Learner(data_bunch, model; opt = Flux.ADAM(), loss = Flux.mse)

A `Learner` is the standard grouping of a data bunch, model, optimizer, and loss.
"""
mutable struct Learner <: AbstractLearner
    cbs:: Array{AbstractCallback}
    db::DataBunch
    model
    opt
    loss
end
Learner(data_bunch, model; opt=Flux.ADAM(), loss=Flux.mse) = Learner([],data_bunch,model,opt,loss)

"""
Get the data bunch for `l`.
"""
data_bunch(l::Learner) = l.db
"""
Set the data bunch for `l` to `data_bunch`.
"""
data_bunch!(l::Learner, data_bunch) = l.db = data_bunch

"""
Get the model for `l`.
"""
model(l::Learner) = l.model
"""
Set the model for `l` to `model`.
"""
model!(l::Learner, model) = l.model = model

"""
Get the loss for `l`.
"""
loss(l::Learner) = l.loss
"""
Set the loss for `l` to `loss`.
"""
loss!(l::Learner,loss) = l.loss = loss

"""
Get the optimizer for `l`.
"""
opt(l::Learner) = l.opt
"""
Set the optimizer for `l` to `opt`.
"""
opt!(l::Learner,opt) = l.opt=opt

"""
Add `cb` to the list of callbacks for `learner`.
"""
add_cb!(learner::Learner,cb::AbstractCallback) = push!(learner.cbs,cb)
"""
Get the list of callbacks for `learner`.
"""
cbs(learner::Learner) = learner.cbs

# pass event to all callbacks
_cbs_before_fit(learner::Learner, epoch_count, batch_size) =  for cb in learner.cbs before_fit(cb,learner,epoch_count,batch_size) end
_cbs_after_fit(learner::Learner) =  for cb in learner.cbs after_fit(cb,learner) end
_cbs_after_cancel_fit(learner::Learner) =  for cb in learner.cbs after_cancel_fit(cb,learner) end

_cbs_before_epoch(learner::Learner, epoch) =  for cb in learner.cbs before_epoch(cb,learner,epoch) end
_cbs_after_epoch(learner::Learner, epoch) =  for cb in learner.cbs after_epoch(cb,learner,epoch) end
_cbs_after_cancel_epoch(learner::Learner, epoch) =  for cb in learner.cbs after_cancel_epoch(cb,learner,epoch) end

_cbs_before_epoch_train(learner::Learner, epoch) =  for cb in learner.cbs before_epoch_train(cb,learner,epoch) end
_cbs_after_epoch_train(learner::Learner, epoch) =  for cb in learner.cbs after_epoch_train(cb,learner,epoch) end
_cbs_after_cancel_epoch_train(learner::Learner, epoch) =  for cb in learner.cbs after_cancel_epoch_train(cb,learner,epoch) end

_cbs_before_batch_train(learner::Learner, epoch, batch) =  for cb in learner.cbs before_batch_train(cb,learner, epoch, batch) end
_cbs_batch_train_loss(learner::Learner, epoch, batch, loss) =  for cb in learner.cbs batch_train_loss(cb,learner, epoch, batch, loss) end
_cbs_after_batch_train(learner::Learner, batch, epoch) =  for cb in learner.cbs after_batch_train(cb,learner, epoch, batch) end
_cbs_after_cancel_batch_train(learner::Learner, batch, epoch) =  for cb in learner.cbs after_cancel_batch_train(cb,learner, epoch, batch) end

_cbs_before_epoch_validate(learner::Learner, epoch) =  for cb in learner.cbs before_epoch_validate(cb,learner,epoch) end
_cbs_after_epoch_validate(learner::Learner, epoch) =  for cb in learner.cbs after_epoch_validate(cb,learner,epoch) end
_cbs_after_cancel_epoch_validate(learner::Learner, epoch) =  for cb in learner.cbs after_cancel_epoch_validate(cb,learner,epoch) end

_cbs_before_batch_validate(learner::Learner, epoch, batch) =  for cb in learner.cbs before_batch_validate(cb,learner, epoch, batch) end
_cbs_after_batch_validate(learner::Learner, epoch, batch) =  for cb in learner.cbs after_batch_validate(cb,learner, epoch, batch) end
_cbs_batch_validate_loss(learner::Learner, epoch, batch, loss) =  for cb in learner.cbs batch_validate_loss(cb,learner, epoch, batch, loss) end
_cbs_after_cancel_batch_validate(learner::Learner, epoch, batch) =  for cb in learner.cbs after_cancel_batch_validate(cb,learner, epoch, batch) end

function _do_batch_train(learner::Learner, batch_data, epoch, batch, ps)
    try
        _loss(xy) = learner.loss(learner.model(xy[1]),xy[2])
        _cbs_before_batch_train(learner,epoch, batch)
        loss, back = Zygote.pullback(ps) do
            sum(_loss.(batch_data))
        end
        _cbs_batch_train_loss(learner, epoch, batch, loss)
        gs = back(1)
        update!(learner.opt, ps, gs)
    #catch CancelBatchTrainException
    #    _cbs_after_cancel_batch_train(learner,batch_index,epoch)
    finally
        _cbs_after_batch_train(learner,batch,epoch)
    end
end

function _do_epoch_train(learner::Learner, epoch, ps)
    try
        data = learner|> data_bunch |> train
        _cbs_before_epoch_train(learner,epoch)
        for (batch,batch_data) in enumerate(data)
            _do_batch_train(learner,batch_data,epoch,batch,ps)
        end
    #catch CancelEpochTrainException
    #    _cbs_after_cancel_epoch_train(learner,epoch)
    finally
        _cbs_after_epoch_train(learner,epoch)
    end
end

function _do_batch_validate(learner::Learner, batch_data, epoch, batch)
    try
        _loss(xy) = learner.loss(learner.model(xy[1]),xy[2])
        _cbs_before_batch_validate(learner,epoch,batch)
        loss = sum(_loss.(batch_data))
        _cbs_batch_validate_loss(learner,epoch,batch,loss)
    #catch CancelBatchValidateException
    #    _cbs_after_cancel_batch_validate(learner,epoch,batch)
    finally
        _cbs_after_batch_validate(learner,epoch,batch)
    end
end

function _do_epoch_validate(learner::Learner, epoch)
    try
        _cbs_before_epoch_validate(learner,epoch)
        data = learner|> data_bunch |> valid
        _cbs_before_epoch_validate(learner,epoch)
        for (batch,batch_data) in enumerate(data)
            _do_batch_validate(learner,batch_data,epoch,batch)
        end
    #catch CancelEpochValidateException
    #    _cbs_after_cancel_epoch_validate(learner,epoch)
    finally
        _cbs_after_epoch_validate(learner,epoch)
    end
end

"""
    fit(learner::Learner, epoch_count)

Fit [`model(learner)`](@ref model(::Learner)) for `epoch_count` epochs invoking callbacks for `learner`.
"""
function fit!(learner::Learner, epoch_count)
    try
        batch_size = train(data_bunch(learner)).batchsize
        ps = params(learner.model)
        #@infiltrate
        _cbs_before_fit(learner,epoch_count,batch_size)
        for epoch in 1:epoch_count
            try
                _cbs_before_epoch(learner,epoch)
                _do_epoch_train(learner,epoch,ps)
                _do_epoch_validate(learner,epoch)
            #catch CancelEpochException
            #    _cbs_after_cancel_epoch(learner,epoch)
            finally
                _cbs_after_epoch(learner,epoch)
            end
        end
    #catch CancelFitException
    #    _cbs_after_cancel_fit(learner)
    finally
        _cbs_after_fit(learner)
    end
end

"""
    implements_learner(T::DataType)

Test if a type implements the [`AbstractLearner`](@ref) interface.
"""
function implements_learner(T::DataType)
    return hasmethod(model,(T,)) &&
        hasmethod(model!,(T,Any)) &&
        hasmethod(loss,(T,)) &&
        hasmethod(loss!,(T,Any)) &&
        hasmethod(opt,(T,)) &&
        hasmethod(opt!,(T,Any)) &&
        hasmethod(cbs,(T,)) &&
        hasmethod(add_cb!,(T,AbstractCallback))
end

@assert implements_learner(Learner)


text:
### make U-Net ###

function make_model_unet(N)
    layers = Any[]
    maxDepth = minimum(round.(Int,log2.(N)).-1)
    numLayers = min(maxDepth, 3)

    H = zeros(Int, numLayers, length(N))
    H[1,:] .= N
    for l=1:numLayers-1
      H[l+1,:] .= ceil.(Int, vec(H[l,:])./2)
    end
    needCrop = Int.(isodd.(H))

    interp = :nearest #trilinear

    inChan = 1
    outChan = 64
    for l=1:numLayers
      push!(layers,UNetConvBlock(inChan,outChan,kernel=(3,3,3), stride=(2,2,2), pad=1))
      push!(layers,UNetConvBlock(outChan,outChan,kernel=(3,3,3), pad=1))
      inChan = outChan
      outChan *= 2
    end
    outChan = inChan
    push!(layers, Upsample(interp,scale=(2,2,2)))
    if sum(needCrop[:,end])>0
      push!(layers, Crop( (0,needCrop[end,1],0,needCrop[end,2],0,needCrop[end,3]) ) )
    end
    push!(layers, UNetConvBlock(inChan,outChan, kernel=(3,3,3), pad=1))
    push!(layers, UNetConvBlock(inChan, outChan, kernel=(1,1,1), pad=0))
    outChan = outChan ÷ 2
    for l=1:(numLayers-1)
      push!(layers, Upsample(interp,scale=(2,2,2)))
      if sum(needCrop[:,end-l])>0
        push!(layers, Crop( (0,needCrop[end-l,1],0,needCrop[end-l,2],0,needCrop[end-l,3]) ) )
      end
      push!(layers, UNetConvBlock(inChan,outChan,kernel=(3,3,3), pad=1))
      push!(layers, UNetConvBlock(outChan,outChan,kernel=(1,1,1), pad=0))
      inChan = outChan
      outChan = outChan ÷ 2
    end
    push!(layers, Conv((1, 1, 1), 64=>1, pad = 0, relu ;init=_random_normal))

    return Chain(layers...)
  end