text:
############################
#          SOLVE           #
############################
function solve!(prob::Problem, solver::ProjectedNewtonSolver)
    to = solver.stats[:timer]
    for i = 1:solver.opts.n_steps
        @timeit to "newton step" V_ = newton_step!(prob, solver)
        @timeit to "copy" begin
            copyto!(prob.X, V_.X)
            copyto!(prob.U, V_.U)
        end

        record_iteration!(prob,solver)
        solver.stats[:c_max][end] <= solver.opts.feasibility_tolerance ? break : nothing
    end

    return solver
end

function record_iteration!(prob::Problem, solver::ProjectedNewtonSolver)
    J = cost(prob)
    c_max = max_violation(prob)

    solver.stats[:iterations] += 1
    push!(solver.stats[:cost],J)
    push!(solver.stats[:c_max],c_max)
end

cost(prob::Problem, V::Union{PrimalDual,PrimalDualVars}) = cost(prob.obj, V.X, V.U, get_dt_traj(prob,V.U))

############################
#       CONSTRAINTS        #
############################
function dynamics_constraints!(prob::Problem, solver::DirectSolver, V=solver.V)
    N = prob.N
    X,U,dt = V.X, V.U, get_dt_traj(prob, V.U)
    solver.fVal[1] .= V.X[1] - prob.x0
    for k = 1:N-1
         evaluate!(solver.fVal[k+1], prob.model, X[k], U[k], dt[k])
         solver.fVal[k+1] .-= X[k+1]
     end
 end

function dynamics_jacobian!(prob::Problem, solver::ProjectedNewtonSolver, V=solver.V)
    n,m,N = size(prob)
    X,U, dt = V.X, V.U, get_dt_traj(prob,V.U)
    solver.∇F[1].xx .= Diagonal(I,n)
    solver.Y[1:n,1:n] .= Diagonal(I,n)
    part = (x=1:n, u =n .+ (1:m), x1=n+m .+ (1:n))
    p = num_constraints(prob)
    off1 = n
    off2 = 0
    for k = 1:N-1
        jacobian!(solver.∇F[k+1], prob.model, X[k], U[k], dt[k])
        solver.Y[off1 .+ part.x, off2 .+ part.x] .= solver.∇F[k+1].xx
        solver.Y[off1 .+ part.x, off2 .+ part.u] .= solver.∇F[k+1].xu
        solver.Y[off1 .+ part.x, off2 .+ part.x1] .= -Diagonal(I,n)
        off1 += n + p[k]
        off2 += n+m
    end
end

function update_constraints!(prob::Problem, solver::DirectSolver, V=solver.V)
    n,m,N = size(prob)
    for k = 1:N-1
        evaluate!(solver.C[k], prob.constraints[k], V.X[k], V.U[k])
    end
    evaluate!(solver.C[N], prob.constraints[N], V.X[N])
end

function active_set!(prob::Problem, solver::ProjectedNewtonSolver)
    n,m,N = size(prob)
    P = sum(num_constraints(prob)) + n*N
    a0 = copy(solver.a)
    for k = 1:N
        active_set!(solver.active_set[k], solver.C[k], solver.opts.active_set_tolerance)
    end
    if solver.opts.verbose && a0 != solver.a
        println("active set changed")
    end
end

function active_set!(a::AbstractVector{Bool}, c::AbstractArray{T}, tol::T=0.0) where T
    a0 = copy(a)
    equality, inequality = c.parts[:equality], c.parts[:inequality]
    a[equality] .= true
    a[inequality] .= c.inequality .>= -tol
end

######################################
#       CONSTRAINT JACBOBIANS        #
######################################
function constraint_jacobian!(prob::Problem, ∇C::Vector, X, U)
    n,m,N = size(prob)
    for k = 1:N-1
        jacobian!(∇C[k], prob.constraints[k], X[k], U[k])
    end
    jacobian!(∇C[N], prob.constraints[N], X[N])
end
constraint_jacobian!(prob::Problem, solver::DirectSolver, V=solver.V) =
    constraint_jacobian!(prob, solver.∇C, V.X, V.U)
# constraint_jacobian!(prob::Problem, Y::KKTJacobian, V=solver.V) =
#     constraint_jacobian!(prob, Y.∇C, V.X, V.U)

function active_constraints(prob::Problem, solver::ProjectedNewtonSolver)
    n,m,N = size(prob)
    a = solver.a.duals
    # return view(solver.Y, a, :), view(solver.y, a)
    return solver.Y.blocks[a,:], solver.y[a]
end

############################
#      COST EXPANSION      #
############################
function cost_expansion!(prob::Problem, solver::ProjectedNewtonSolver, V=solver.V) where T
    n,m,N = size(prob)
    NN = N*n + (N-1)*m
    H = solver.H
    g = solver.g
    dt = get_dt_traj(prob,V.U)

    part = (x=1:n, u=n .+ (1:m), z=1:n+m)
    part2 = (xx=(part.x, part.x), uu=(part.u, part.u), ux=(part.u, part.x), xu=(part.x, part.u))
    off = 0
    for k = 1:N-1
        # H[off .+ part.x, off .+ part.x] = Q[k].xx
        # H[off .+ part.x, off .+ part.u] = Q[k].ux'
        # H[off .+ part.u, off .+ part.x] = Q[k].ux
        # H[off .+ part.u, off .+ part.u] = Q[k].uu
        hess = PartedMatrix(view(H, off .+ part.z, off .+ part.z), part2)
        grad = PartedVector(view(g, off .+ part.z), part)
        hessian!(hess, prob.obj[k], V.X[k], V.U[k],dt[k])
        gradient!(grad, prob.obj[k], V.X[k], V.U[k],dt[k])

        off += n+m
    end

    hess = PartedMatrix(view(H, off .+ part.x, off .+ part.x), part2)
    grad = PartedVector(view(g, off .+ part.x), part)
    hessian!(hess, prob.obj[N], V.X[N])
    gradient!(grad, prob.obj[N], V.X[N])
end

######################
#     FUNCTIONS      #
######################
function update!(prob::Problem, solver::ProjectedNewtonSolver, V=solver.V, active_set=true)
    to = solver.stats[:timer]
    @timeit to "dynamics"       dynamics_constraints!(prob, solver, V)
    @timeit to "constraints"    update_constraints!(prob, solver, V)
    @timeit to "dyn jacobian"   dynamics_jacobian!(prob, solver, V)
    @timeit to "con jacobian"   constraint_jacobian!(prob, solver, V)
    @timeit to "cost expansion" cost_expansion!(prob, solver, V)
    if active_set
        @timeit to "active set" active_set!(prob, solver)
    end
end

function max_violation(solver::DirectSolver{T}) where T
    c_max = 0.0
    C = solver.C
    N = length(C)
    for k = 1:N
        if length(C[k].equality) > 0
            c_max = max(norm(C[k].equality,Inf), c_max)
        end
        if length(C[k].inequality) > 0
            c_max = max(pos(maximum(C[k].inequality)), c_max)
        end
        c_max = max(norm(solver.fVal[k], Inf), c_max)
    end
    return c_max
end

function calc_violations(solver::ProjectedNewtonSolver{T}) where T
    C = solver.C
    N = length(C)
    v = [zero(c) for c in C]
    v = zeros(N)
    for k = 1:N
        if length(C[k].equality) > 0
            v[k] = norm(C[k].equality,Inf)
        end
        if length(C[k].inequality) > 0
            v[k] = max(pos(maximum(C[k].inequality)), v[k])
        end
    end
    return v
end

function projection_solve!(prob, solver, V=solver.V, active_set_update=true)
    eps_feasible = solver.opts.feasibility_tolerance
    viol = norm(solver.y[solver.a.duals], Inf)
    max_projection_iters = 10

    count = 0
    while count < max_projection_iters && viol > eps_feasible
        viol = _projection_solve!(prob, solver, V, active_set_update)
        count += 1
    end
    return viol
end

function _projection_solve!(prob::Problem, solver::ProjectedNewtonSolver,
        V=solver.V, active_set_update=true)
    to = solver.stats[:timer]

    Z = primals(V)
    λ = duals(V)
    a = solver.a.duals
    max_refinements = 10
    convergence_rate_threshold = 1.1
    ρ = 1e-2

    # cost_expansion!(prob, solver, V)
    H = Diagonal(solver.H)

    @timeit to "dynamics"     dynamics_constraints!(prob, solver, V)
    @timeit to "constraints"  update_constraints!(prob, solver, V)
    @timeit to "dyn jacobian" dynamics_jacobian!(prob, solver, V)
    @timeit to "con jacobian" constraint_jacobian!(prob, solver, V)
    if active_set_update
        @timeit to "active set" active_set!(prob, solver)
    end
    Y,y = active_constraints(prob, solver)
    viol0 = norm(y,Inf)
    if solver.opts.verbose
        println("feas0: $viol0")
    end

    HinvY = H\Y'
    S = Symmetric(Y*HinvY)
    @timeit to "cholesky" Sreg = cholesky(S + ρ*I)
    viol_prev = viol0
    count = 0
    while count < max_refinements
        @timeit to "linesearch" viol = _projection_linesearch!(prob, solver, V, (S,Sreg), HinvY)
        convergence_rate = log10(viol)/log10(viol_prev)
        viol_prev = viol
        count += 1

        if solver.opts.verbose
            println("conv rate: $convergence_rate")
        end

        if convergence_rate < convergence_rate_threshold ||
                       viol < solver.opts.feasibility_tolerance
            break
        end
    end

    solver.stats[:S] = Sreg
    return viol_prev
end

function _projection_linesearch!(prob::Problem, solver::ProjectedNewtonSolver,
        V, S, HinvY)
    to = solver.stats[:timer]

    a = solver.a.duals
    y = solver.y[a]
    viol0 = norm(y,Inf)
    viol = Inf
    ρ = 1e-4

    Z = primals(V)
    V_ = copy(V)
    Z_ = primals(V_)
    α = 1.0
    ϕ = 0.5
    count = 1
    while true
        @timeit to "dual solve" δλ = reg_solve(S[1],y,S[2],1e-8,25)
        @timeit to "primal solve" δZ = -HinvY*δλ
        Z_ .= Z + α*δZ

        @timeit to "dynamics" dynamics_constraints!(prob, solver, V_)
        @timeit to "constraints" update_constraints!(prob, solver, V_)
        y = solver.y[a]
        viol = norm(y,Inf)

        if solver.opts.verbose
            println("feas: $viol")
        end
        if viol < viol0 || count > 10
            break
        else
            count += a
            α *= ϕ
        end
    end
    copyto!(Z,Z_)
    return viol
end

reg_solve(A, b, reg::Real, tol=1e-10, max_iters=10) = reg_solve(A, b, A + reg*I, tol, max_iters)
function reg_solve(A, b, B, tol=1e-10, max_iters=10)
    x = B\b
    count = 0
    while count < max_iters
        r = b - A*x
        # println("r_norm = $(norm(r))")

        if norm(r) < tol
            break
        else
            x += B\r
            count += 1
        end
    end
    # println("iters = $count")

    return x
end

function projection!(prob::Problem, solver::ProjectedNewtonSolver, V=solver.V, active_set_update=true)
    Z = primals(V)
    eps_feasible = solver.opts.feasibility_tolerance
    count = 0
    # cost_expansion!(prob, solver, V)
    H = Diagonal(solver.H)
    while true
        dynamics_constraints!(prob, solver, V)
        update_constraints!(prob, solver, V)
        dynamics_jacobian!(prob, solver, V)
        constraint_jacobian!(prob, solver, V)
        if active_set_update
            active_set!(prob, solver)
        end
        Y,y = active_constraints(prob, solver)
        HinvY = H\Y'

        viol = norm(y,Inf)
        if solver.opts.verbose
            println("feas: ", viol)
        end
        if viol < eps_feasible || count > 10
            break
        else
            δZ = -HinvY*((Y*HinvY)\y)
            Z .+= δZ
            count += 1
        end
    end
end

function primaldual_projection!(prob::Problem, solver::ProjectedNewtonSolver, V=solver.V, active_set_update=true)
    Z = primals(V)
    λ = duals(V)
    a = solver.a.duals
    eps_feasible = solver.opts.feasibility_tolerance
    count = 0
    ρ = 1e-6

    # cost_expansion!(prob, solver, V)
    H = Diagonal(solver.H)

    while true
        dynamics_constraints!(prob, solver, V)
        update_constraints!(prob, solver, V)
        dynamics_jacobian!(prob, solver, V)
        constraint_jacobian!(prob, solver, V)
        if active_set_update
            active_set!(prob, solver)
        end
        Y,y = active_constraints(prob, solver)
        HinvY = H\Y'

        viol = norm(y,Inf)
        if solver.opts.verbose
            println("feas: ", viol)
        end
        if viol < eps_feasible || count > 10
            if count == 0
                solver.stats[:S] = cholesky(Symmetric(Y*HinvY) + ρ*I)
            end
            break
        else
            λa = view(λ,a)

            S = cholesky(Symmetric(Y*HinvY) + ρ*I)
            δλ = S\y
            δZ = -HinvY*δλ

            Z .+= δZ
            λa .+= δλ
            count += 1

            solver.stats[:S] = S
        end
    end
end

function multiplier_projection!(prob::Problem, solver::ProjectedNewtonSolver, V=solver.V)
    g = solver.g
    a = solver.a.duals
    Y,y = active_constraints(prob, solver)
    λ = duals(V)[a]

    res0 = g + Y'λ
    δλ = -(Y*Y')\(Y*res0)
    λ_ = λ + δλ
    res = g + Y'*λ_
    copyto!(view(duals(V),a), λ_)
    res = norm(residual(prob, solver, V))
    return res, δλ
end

function solveKKT(prob::Problem, solver::ProjectedNewtonSolver, V=solver.V)
    a = solver.a
    δV = zero(V.V)
    λ = duals(V)[a.duals]
    Y,y = active_constraints(prob, solver)
    H,g = solver.H, solver.g
    Pa = length(y)
    A = [H Y'; Y zeros(Pa,Pa)]
    b = [g + Y'λ; y]
    δV[a] = -A\b
    return δV
end

using SuiteSparse
function solveKKT_Shur(prob::Problem, solver::ProjectedNewtonSolver, Hinv, V=solver.V)
    a = solver.a
    δV = zero(V.V)
    λ = duals(V)[a.duals]
    Y,y = active_constraints(prob, solver)
    g = solver.g
    r = g + Y'λ

    YHinv = Y*Hinv
    L = solver.stats[:S]::SuiteSparse.CHOLMOD.Factor{Float64}
    δλ = L\(y-YHinv*r)
    δz = -Hinv*(r+Y'δλ)

    δV[solver.parts.primals] .= δz
    δV[solver.parts.duals[solver.a.duals]] .= δλ
    return δV, δλ
end

function residual(prob::Problem, solver::ProjectedNewtonSolver, V=solver.V)
    a = solver.a
    Y,y = active_constraints(prob, solver)
    λ = duals(V)[a.duals]
    g = solver.g
    res = [g + Y'λ; y]
end

function line_search(prob::Problem, solver::ProjectedNewtonSolver, δV)
    α = 1.0
    s = 0.01
    J0 = cost(prob, solver.V)
    update!(prob, solver)
    res0 = norm(residual(prob, solver))
    count = 0
    solver.opts.verbose ? println("res0: $res0") : nothing
    while count < 10
        V_ = solver.V + α*δV

        # Calculate residual
        projection!(prob, solver, V_)
        # projection_solve!(prob, solver, V_)

        cost_expansion!(prob, solver, V_)
        res, = multiplier_projection!(prob, solver, V_)
        J = cost(prob, V_)

        # Calculate max violation
        viol = max_violation(solver)

        if solver.opts.verbose
            println("cost: $J \t residual: $res \t feas: $viol")
        end
        if res < (1-α*s)*res0
            solver.opts.verbose ? println("α: $α") : nothing
            return V_
        end
        count += 1
        α /= 2
    end
    return solver.V
end

function newton_step!(prob::Problem, solver::ProjectedNewtonSolver)
    V = solver.V
    verbose = solver.opts.verbose

    # Initial stats
    update!(prob, solver)
    J0 = cost(prob, V)
    res0 = norm(residual(prob, solver))
    viol0 = max_violation(solver)
    Hinv = inv(Diagonal(solver.H))

    # Projection
    verbose ? println("\nProjection:") : nothing
    # primaldual_projection!(prob, solver)
    projection_solve!(prob, solver)

    if solver.opts.solve_type == :feasible
        return solver.V
    end

    cost_expansion!(prob, solver)
    multiplier_projection!(prob, solver)

    # Solve KKT
    J1 = cost(prob, V)
    res1 = norm(residual(prob, solver))
    viol1 = max_violation(solver)
    δV, = solveKKT_Shur(prob, solver, Hinv)

    # Line Search
    verbose ? println("\nLine Search") : nothing
    V_ = line_search(prob, solver, δV)
    J_ = cost(prob, V_)
    res_ = norm(residual(prob, solver, V_))
    viol_ = max_violation(solver)

    # Print Stats
    if verbose
        println("\nStats")
        println("cost: $J0 → $J1 → $J_")
        println("res: $res0 → $res1 → $res_")
        println("viol: $viol0 → $viol1 → $viol_")
    end

    return V_
end

function buildL(L::KKTFactors,y_part)
    Pa = sum(y_part)
    S = PseudoBlockArray(zeros(Pa,Pa),y_part,y_part)
    N = length(y_part) ÷ 2

    S[Block(1,1)] = L.G[end].L
    S[Block(2,1)] = L.F[1]
    S[Block(2,2)] = L.G[1].L
    S[Block(3,1)] = L.L[1]
    S[Block(3,2)] = L.M[1]
    S[Block(3,3)] = L.H[1].L
    for k = 2:N-1
        S[Block(2k,2k-2)] = L.E[k]
        S[Block(2k,2k-1)] = L.F[k]
        S[Block(2k,2k  )] = L.G[k].L
        S[Block(2k+1,2k-2)] = L.K[k]
        S[Block(2k+1,2k-1)] = L.L[k]
        S[Block(2k+1,2k-0)] = L.M[k]
        S[Block(2k+1,2k+1)] = L.H[k].L
    end
    S[Block(2N,2N-2)] = L.E[N]
    S[Block(2N,2N-1)] = L.F[N]
    S[Block(2N,2N-0)] = L.G[N].L
    return S
end
buildL(solver::SequentialNewtonSolver) = buildL(solver.L, dual_partition(solver))

function buildY!(Y,solver::SequentialNewtonSolver)
    n,m,N = size(solver)

    # Get Active Jacobians
    xi,ui = 1:n, n .+ (1:m)
    a = solver.active_set
    C = [solver.∇C[k][a[k],1:n] for k = 1:N]
    D = [solver.∇C[k][a[k],n+1:n+m] for k = 1:N-1]

    Y[Block(1,1)] = Diagonal(I,n)
    for k = 1:N-1
        Y[Block(2k,2k-1)] = solver.∇F[k+1].xx
        Y[Block(2k,2k-0)] = solver.∇F[k+1].xu
        Y[Block(2k,2k+1)] = -Diagonal(I,n)
        Y[Block(2k+1,2k-1)] = C[k]
        Y[Block(2k+1,2k-0)] = D[k]
    end
    Y[Block(2N,2N-1)] = C[N]
    return Y
end
function buildY(solver::SequentialNewtonSolver)
    y_part = dual_partition(solver)
    z_part = repeat([n,m],N-1)
    push!(z_part,n)
    NN = sum(z_part)
    Pa = sum(y_part)

    Y = PseudoBlockArray(zeros(Pa,NN),y_part,z_part)
    buildY!(Y,solver)
    return Y
end

function buildS!(S,solver::SequentialNewtonSolver,C,D)
    n,m,N = size(solver)
    Qinv = solver.Qinv
    Rinv = solver.Rinv
    ∇F = view(solver.∇F,2:N)
    p = sum.(solver.active_set)
    pcum = insert!(cumsum(p), 1, 0)

    dinds = 1:n
    cinds = n .+ (1:p[1])
    S[dinds,dinds] = Qinv[1]
    S[dinds,n .+ dinds] = Qinv[1]*∇F[1].xx'
    S[dinds,n .+ cinds] = Qinv[1]*C[1]'

    for k = 1:N-1
        off = pcum[k] + k*n
        dinds = off .+ (1:n)
        cinds = off + n .+ (1:p[k])
        S[dinds,dinds] = ∇F[k].xx*Qinv[k]*∇F[k].xx' + ∇F[k].xu*Rinv[k]*∇F[k].xu' + Qinv[k+1]
        S[dinds,cinds] = ∇F[k].xx*Qinv[k]*C[k]'     + ∇F[k].xu*Rinv[k]*D[k]'
        S[cinds,cinds] = C[k]*Qinv[k]*C[k]'         + D[k]*Rinv[k]*D[k]'
        if k < N-1
            S[dinds,dinds .+ (p[k] + n)] = -Qinv[k+1]*∇F[k+1].xx'
            S[dinds, (off + p[k] + 2n) .+ (1:p[k+1])] = -Qinv[k+1]*C[k+1]'
        else
            S[dinds, (off + p[k] + n) .+ (1:p[k+1])] = -Qinv[k+1]*C[k+1]'
        end
    end
    off = pcum[N] + N*n
    cinds = off .+ (1:p[N])
    S[cinds,cinds] = C[N]*Qinv[N]*C[N]'
    return S
    return Symmetric(S)
end
function buildS(solver::SequentialNewtonSolver)
    y_part = dual_partition(solver)

    Pa = sum(y_part)
    S = zeros(Pa,Pa)

    C = [solver.∇C[k][a[k],1:n] for k = 1:N]
    D = [solver.∇C[k][a[k],n+1:n+m] for k = 1:N-1]

    buildS!(solver,C,D)
    return S
end

function Sinds(solver::SequentialNewtonSolver)
    n,m,N = size(solver)
    y_part = dual_partition(solver)
    inds = create_partition(Tuple(y_part))

    A = [(0:0,0:0) for k = 1:N]
    B = [(0:0,0:0) for k = 1:N]
    C = [(0:0,0:0) for k = 1:N]
    D = [(0:0,0:0) for k = 1:N]
    E = [(0:0,0:0) for k = 1:N]

    n,m,N = size(solver)
    p = sum.(solver.active_set)
    pcum = insert!(cumsum(p), 1, 0)

    dinds = 1:n
    cinds = n .+ (1:p[1])
    A[1] = (dinds,dinds)
    D[1] = (dinds,n .+ dinds)
    E[1] = (dinds,n .+ cinds)

    for k = 1:N-1
        off = pcum[k] + k*n
        dinds = off .+ (1:n)
        cinds = off + n .+ (1:p[k])
        A[k+1] = (dinds,dinds)
        B[k+1] = (dinds,cinds)
        C[k+1] = (cinds,cinds)
        if k < N-1
            D[k+1] = (dinds, dinds .+ (p[k] + n))
            E[k+1] = (dinds, (off + p[k] + 2n) .+ (1:p[k+1]))
        else
            E[k+1] = (dinds, (off + p[k] + n) .+ (1:p[k+1]))
        end
    end
    off = pcum[N] + N*n
    cinds = off .+ (1:p[N])
    C[1,1] = (cinds,cinds)
    S_part = (A=A,B=B,C=C,D=D,E=E)
end

function buildS!(S, solver, A, B, C, D, inds::NamedTuple)
    N = length(solver.Q)
    Qinv = solver.Qinv
    Rinv = solver.Rinv

    S[inds.A[1]...] = Qinv[1]
    S[inds.D[1]...] = Qinv[1]*A[1]'
    S[inds.E[1]...] = Qinv[1]*C[1]'

    for k = 1:N-1
        S[inds.A[k+1]...] = A[k]*Qinv[k]*A[k]' + B[k]*Rinv[k]*B[k]' + Qinv[k+1]
        S[inds.B[k+1]...] = A[k]*Qinv[k]*C[k]'     + B[k]*Rinv[k]*D[k]'
        S[inds.C[k+1]...] = C[k]*Qinv[k]*C[k]'         + D[k]*Rinv[k]*D[k]'
        if k < N-1
            S[inds.D[k+1]...] = -Qinv[k+1]*A[k+1]'
            S[inds.E[k+1]...] = -Qinv[k+1]*C[k+1]'
        else
            S[inds.E[k+1]...] = -Qinv[k+1]*C[k+1]'
        end
    end
    S[inds.C[1]...] = C[N]*Qinv[N]*C[N]'
    Symmetric(S)
    return nothing
    # return Symmetric(S)
end

function buildShurCompliment(prob::Problem, solver::SequentialNewtonSolver)
    n,m,N = size(prob)

    Qinv = solver.Qinv
    Rinv = solver.Rinv

    A = [F.xx for F in solver.∇F[2:end]]  # First jacobian is for initial condition
    B = [F.xu for F in solver.∇F[2:end]]
    C = [Array(F.x[a,:]) for (F,a) in zip(solver.∇C, solver.active_set)]
    D = [Array(F.u[a,:]) for (F,a) in zip(solver.∇C, solver.active_set)]

    P = num_active_constraints(solver)
    S = spzeros(P,P)

    _buildShurCompliment!(S, prob, solver, Qinv, Rinv, A, B, C, D)

    return Symmetric(S)
end

function _buildShurCompliment!(S, prob::Problem, solver::SequentialNewtonSolver, Qinv, Rinv, A, B, C, D)
    n,m,N = size(prob)

    p = sum.(solver.active_set)
    pcum = insert!(cumsum(p), 1, 0)

    dinds = 1:n
    cinds = n .+ (1:p[1])
    S[dinds,dinds] = Qinv[1]
    S[dinds,n .+ dinds] = Qinv[1]*A[1]'
    S[dinds,n .+ cinds] = Qinv[1]*C[1]'

    for k = 1:N-1
        off = pcum[k] + k*n
        dinds = off .+ (1:n)
        cinds = off + n .+ (1:p[k])
        S[dinds,dinds] = A[k]*Qinv[k]*A[k]' + B[k]*Rinv[k]*B[k]' + Qinv[k+1]
        S[dinds,cinds] = A[k]*Qinv[k]*C[k]' + B[k]*Rinv[k]*D[k]'
        S[cinds,cinds] = C[k]*Qinv[k]*C[k]' + D[k]*Rinv[k]*D[k]'
        if k < N-1
            S[dinds,dinds .+ (p[k] + n)] = -Qinv[k+1]*A[k+1]'
            S[dinds, (off + p[k] + 2n) .+ (1:p[k+1])] = -Qinv[k+1]*C[k+1]'
        else
            S[dinds, (off + p[k] + n) .+ (1:p[k+1])] = -Qinv[k+1]*C[k+1]'
        end
    end
    off = pcum[N] + N*n
    cinds = off .+ (1:p[N])
    S[cinds,cinds] = C[N]*Qinv[N]*C[N]'

end


text:
# Functions relating to the decomposition method

# Aliases
const DecBelief = Dict{Int64, SparseCat}
const DecState = Dict{Int64, SingleOCState}
const DecObs = Dict{Int64, SingleOCObs}

"""
    POMDPs.value(policy::AlphaVectorPolicy, b::SparseCat)
Return the value of the belief according to the QMDP policy
"""
function POMDPs.value(policy::AlphaVectorPolicy, b::SparseCat)
    val = zeros(length(actions(policy.pomdp)))
    for (i,s) in enumerate(b.vals)
        val += value(policy, s)*b.probs[i]
        # si = stateindex(pomdp, s)
        # val += policy.alphas[si, :]*b.p[i]
    end
    return val
end

# function POMDPs.value(policy::SARSOPPolicy, b::SparseCat)
#     pomdp = policy.pomdp
#     actions = policy.alphas.alpha_actions
#     vectors = alphas(policy)
#     n = size(vectors, 2)
#     utilities = zeros(n_actions(pomdp), n)
#     for i=1:n
#         act = actions[i] + 1
#         for (j, s) in enumerate(b.it)
#             s_itps, w_itps = interpolate_state(s, pomdp)
#             for (k, s_itp) in enumerate(s_itps)
#                 si = stateindex(pomdp, s_itp)
#                 utilities[act, i] += vectors[si,i]*b.p[j]*w_itps[k]
#             end
#         end
#     end
#     return maximum(utilities,2)
# end

"""
    POMDPs.value(policy::AlphaVectorPolicy, s::SingleOCState)
return the value of a given state using grid interpolations
"""
function POMDPs.value(policy::AlphaVectorPolicy, s::SingleOCState)
    # println("COUCOU")
    pomdp = policy.pomdp
    val = zeros(length(actions(pomdp)))

    ego_grid = RectangleGrid(get_X_grid(pomdp), get_V_grid(pomdp)) #XXX must not allSingleOCate the grid at each function call, find better implementation
    ego_indices, ego_weights = interpolants(ego_grid, [s.ego.posG.x, s.ego.v])

    if !off_the_grid(pomdp, s.ped)
        ped_grid = RectangleGrid(get_Y_grid(pomdp), get_V_ped_grid(pomdp))
        ped_indices, ped_weights = interpolants(ped_grid, [s.ped.posG.y, s.ped.v])
    else
        ped_indices = Int64[0]
        ped_weights = Float64[1.0]
    end

    sum_weight = 0.
    for i=1:length(ego_indices)
        for j=1:length(ped_indices)
            xg, vg = ind2x(ego_grid, ego_indices[i])
            ego = xv_to_state(pomdp, xg, vg)
            if ped_indices[1] == 0
                ped = get_off_the_grid(pomdp)
            else
                yg, v_pedg = ind2x(ped_grid, ped_indices[j])
                ped = yv_to_state(pomdp, yg, v_pedg)
            end
            state = SingleOCState(collision_checker(ego, ped, pomdp.ego_type, pomdp.ped_type),ego,ped)
            si = stateindex(pomdp, state)
            sum_weight += ego_weights[i]*ped_weights[j]
            @show size([a[si] for a in policy.alphas]*ego_weights[i]*ped_weights[j])
            val .+= [a[si] for a in policy.alphas]*ego_weights[i]*ped_weights[j]
        end
    end
    val /= sum_weight
    return val
end

function POMDPs.action(policy::AlphaVectorPolicy, b::SparseCat)
    alphas = policy.alphas' #size |A|x|S|
    util = zeros(n_actions(policy.pomdp)) # size |A|
    for i=1:n_actions(policy.pomdp)
        res = 0.0
        for (j,s) in enumerate(b.vals)
            si = stateindex(policy.pomdp, s)
            res += alphas[i][si]*b.probs[j]
        end
        util[i] = res
    end
    ihi = argmax(util)
    return policy.action_map[ihi]
end

# function POMDPs.action(policy::SARSOPPolicy, b::SingleOCDistribution)
#     pomdp = policy.pomdp
#     val = zeros(n_actions(pomdp))
#     actions = policy.alphas.alpha_actions
#     vectors = alphas(policy)
#     n = size(vectors, 2)
#     utilities = zeros(n)
#     for i=1:n
#         res = 0.
#         for (j, s) in enumerate(b.it)
#             si = stateindex(pomdp, s)
#             res += vectors[si, i]*b.p[j]
#         end
#         utilities[i] = res
#     end
#     a = actions[argmax(utilities)] + 1
#     return policy.action_map[a]
# end

# function POMDPs.action(policy::SARSOPPolicy, b::SingleOCDistribution)
#     max_val, ind_max = fargmax(value(policy, b))
#     return policy.action_map[ind_max]
# end

# function POMDPs.action(policy::Policy, b::Dict{Int64, B}) where B
#     val = fuse_value_min(policy, b)
#     max_ind = argmax(val)
#     return policy.action_map[max_ind]
# end

"""
    fuse_value(policy::Policy, b::Dict{Int64,SparseCat})
Fuse the value of different belief by summing them
"""
function fuse_value(policy::Policy, b::Dict{Int64,B}) where B
    val = zeros(n_actions(policy.pomdp))
    for id in keys(b)
        val_inc = value(policy, b[id])
        val += val_inc
        # println("For pedestrian ", id, "value is", val_inc)
    end
    return val
end

"""
    fuse_value_min(policy::Policy, b::Dict{Int64,SparseCat})
Fuse the value of different belief by taking the min of the two for each action
"""
function fuse_value_min(policy::Policy, b::Dict{Int64, B}) where B
    val = zeros(n_actions(policy.pomdp))
    for i=1:n_actions(policy.pomdp)
        min_val = +Inf
        for id in keys(b)
            v_id = value(policy, b[id])
            # println("For pedestrian ", id, "value is", v_id)
            v_id = v_id[i]
            if v_id < min_val
                min_val = v_id
            end
        end
        val[i] = min_val
    end
    return val
end

"""
    generate_scene(pomdp::SingleOCPOMDP, scene::Scene, a::SingleOCAction, rng::AbstractRNG)
generate the next scene using the pomdp's transition model
"""
function generate_scene(pomdp::SingleOCPOMDP, scene::Scene, a::SingleOCAction, rng::AbstractRNG)
    states = scene_to_states(scene, pomdp)
    statesp = Dict{Int64, SingleOCState}()
    for id in keys(states)
        statesp[id] = generate_s(pomdp, states[id], a, rng)
    end
    return states_to_scene(statesp, pomdp)
end

"""
    generate_scene(pomdp::SingleOCPOMDP, scene::Scene, a::SingleOCAction, rng::AbstractRNG)
generate a vector of states
"""
function generate_states(pomdp::SingleOCPOMDP, states::Dict{Int64, SingleOCState}, a::SingleOCAction, rng::AbstractRNG)
    statesp = Dict{Int64, SingleOCState}()
    for id in keys(states)
        statesp[id] = gen(DDNNout(:sp), pomdp, states[id], a, rng)
    end
    return statesp
end

"""
POMDPs.gen(::DDNNode{:o}, pomdp::SingleOCPOMDP, states::Dict{Int64, SingleOCState}, a::SingleOCAction, statesp::Dict{Int64, SingleOCState}, rng::AbstractRNG)
generate an observation from a vector of states
"""
function POMDPs.gen(::DDNNode{:o}, pomdp::SingleOCPOMDP, states::Dict{Int64, SingleOCState}, a::SingleOCAction, statesp::Dict{Int64, SingleOCState}, rng::AbstractRNG)
    o = Dict{Int64, SingleOCObs}()
    for id in keys(statesp)
        o[id] = gen(DDNNode(:o), pomdp, states[id], a, statesp[id], rng)
    end
    return o
end

# function POMDPs.update(up::SingleOCUpdater, bold::Dict{Int64, SingleOCDistribution}, a::SingleOCAction, o::Dict{Int64, SingleOCObs})
#     bnew = Dict{Int64, SingleOCDistribution}()
#     #TODO check for new comers
#     for id in keys(bold)
#         bnew[id] = update(up, bold[id], a, o[id])
#     end
#     return bnew
# end
#
# function POMDPs.update(up::Updater, bold::Dict{Int64, SingleOCDistribution}, a::Float64, o::Dict{Int64, SingleOCObs})
#     bnew = Dict{Int64, SingleOCDistribution}()
#     #TODO check for new comers
#     for id in keys(bold)
#         bnew[id] = update(up, bold[id], SingleOCAction(a), o[id])
#     end
#     return bnew
# end

# function POMDPs.update{U, B}(up::U, bold::Dict{Int64, B}, a::Float64, o::Dict{Int64, SingleOCObs})
#     bnew = Dict{Int64, B}()
#     #TODO check for new comers
#     for id in keys(bold)
#         bnew[id] = update(up, bold[id], SingleOCAction(a), o[id])
#     end
#     return bnew
# end

#### HELPERS ######################################################################################

"""
return the closest states and the assSingleOCiated weights using multi linear interpolation
"""
function interpolate_state(s::SingleOCState, pomdp::SingleOCPOMDP)
    ego_grid = RectangleGrid(get_X_grid(pomdp), get_V_grid(pomdp)) #XXX must not allSingleOCate the grid at each function call, find better implementation
    ego_indices, ego_weights = interpolants(ego_grid, [s.ego.posG.x, s.ego.v])

    if !isinf(s.ped.v)
        ped_grid = RectangleGrid(get_Y_grid(pomdp), get_V_ped_grid(pomdp))
        ped_indices, ped_weights = interpolants(ped_grid, [s.ped.posG.y, s.ped.v])
    else
        ped_indices = Int64[0]
        ped_weights = Float64[1.0]
    end

    states = SingleOCState[] #todo preallSingleOCate
    weights = Float64[]
    for i=1:length(ego_indices)
        for j=1:length(ped_indices)
            xg, vg = ind2x(ego_grid, ego_indices[i])
            ego = xv_to_state(pomdp, xg, vg)
            if ped_indices[1] == 0
                ped = get_off_the_grid(pomdp)
            else
                yg, v_pedg = ind2x(ped_grid, ped_indices[j])
                ped = yv_to_state(pomdp, yg, v_pedg)
            end
            state = SingleOCState(collision_checker(ego, ped, pomdp.ego_type, pomdp.ped_type),ego,ped)
            push!(states, state)
            push!(weights, ego_weights[i]*ped_weights[j])
        end
    end
    normalize!(weights, 1)
    return states, weights
end

"""
    scene_to_states(scene::Scene, pomdp::SingleOCPOMDP)
Convert a scene to a state representation, for each pedestrian present in the scene, a state object is created
"""
function scene_to_states(scene::Scene, pomdp::SingleOCPOMDP)
    ego = scene.entities[1].state
    @assert scene.entities[1].id == 1
    states = Dict{Int64, SingleOCState}()
    if scene.n == 1
        states[2] = SingleOCState(false, ego, get_off_the_grid(pomdp))
    end
    for i=1:scene.n - 1
        ped = scene.entities[i+1]
        states[ped.id] = SingleOCState(collision_checker(ego, ped.state, pomdp.ego_type, pomdp.ped_type), ego, ped.state)
    end
    return states
end

"""
    states_to_scene(states::Dict{Int64, SingleOCState}, pomdp::SingleOCPOMDP)
"""
function states_to_scene(states::Dict{Int64, SingleOCState}, pomdp::SingleOCPOMDP)
    scene = Scene()
    for id in keys(states)
        if scene.n == 0
            ego = Vehicle(states[id].ego, pomdp.ego_type, 1)
            push!(scene, ego)
        end
        ped = Vehicle(states[id].ped, pomdp.ped_type, id)
        push!(scene, ped)
    end
    return scene
end

"""
    rand_scene(rng::AbstractRNG, b::Dict{Int64, SparseCat}, env::CrosswalkEnv)
Sample a scene from a belief representation
"""
function rand_scene(rng::AbstractRNG, b::Dict{Int64, SparseCat}, pomdp::SingleOCPOMDP)
    scene = Scene()
    for id in keys(b)
        s = rand(rng, b[id])
        if scene.n == 0
            ego = Vehicle(s.ego, pomdp.ego_type, 1)
            push!(scene, ego)
        end
        ped = Vehicle(s.ped, pomdp.ped_type, id)
        push!(scene, ped)
    end
    return scene
end

"""
    initialize_two_beliefs(pomdp::SingleOCPOMDP)
Initialize a belief considering two pedestrians
"""
function initialize_two_beliefs(pomdp::SingleOCPOMDP)
    b = Dict{Int64, SparseCat}()
    b[2] = initialstate_distribution(pomdp)
    b[3] = initialstate_distribution(pomdp)
    return b
end

ext: jl

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
data_type = Float32

using Mocha
srand(123)

N = 1000
d = 4
lambda = 0.1
sinkhorn_iter = 10

n_hidden_layer = 2
n_hidden_unit = 64

test_wasserstein = true

X = rand(data_type, d, N)
Y_onehot = zeros(data_type, 1, N)
Y = zeros(data_type, d, N)
for i = 1:N
  Y[indmax(X[:,i]),i] = 1
  Y_onehot[i] = indmax(X[:,i])-1
end

# ground metric
M = zeros(d, d)
for i = 1:d
  for j = 1:d
    M[i,j] = abs(i-j)
  end
end

backend = DefaultBackend()
init(backend)

data_layer = MemoryDataLayer(batch_size=50,
    data=Array[X, test_wasserstein ? Y : Y_onehot])

hidden_layers = map(1:n_hidden_layer) do i
  InnerProductLayer(name="ip$i", output_dim=n_hidden_unit, neuron=Neurons.ReLU(),
      tops=[symbol("ip$i")], bottoms=[i == 1 ? :data : symbol("ip$(i-1)")])
end
pred_layer = InnerProductLayer(name="pred", output_dim=d,
    tops=[:pred], bottoms=[symbol("ip$n_hidden_layer")])
if test_wasserstein
  softmax_layer = SoftmaxLayer(bottoms=[:pred], tops=[:softpred])
  loss_layer = WassersteinLossLayer(bottoms=[:softpred, :label],
      M = M, lambda = lambda, sinkhorn_iter = sinkhorn_iter)
else
  loss_layer = SoftmaxLossLayer(bottoms=[:pred, :label])
end

common_layers = [hidden_layers..., pred_layer]
if test_wasserstein
  common_layers = [common_layers..., softmax_layer]
end
net = Net("Training", backend, [data_layer, common_layers..., loss_layer])
method = SGD()
params = make_solver_parameters(method, max_iter=10000,
                                lr_policy=LRPolicy.Fixed(0.01), mom_policy=MomPolicy.Fixed(0.9))
solver = Solver(method, params)
add_coffee_break(solver, TrainingSummary(), every_n_iter=100)

# for simplicity, we just use the training data as test...
data_layer_test = MemoryDataLayer(batch_size=100, data=Array[X, Y_onehot])
acc_layer = AccuracyLayer(bottoms=[:softpred, :label])
test_layers = [data_layer_test, common_layers..., acc_layer]
if !test_wasserstein
  softmax_layer = SoftmaxLayer(bottoms=[:pred], tops=[:softpred])
  test_layers = [test_layers..., softmax_layer]
end
test_net = Net("Testing", backend, test_layers)
add_coffee_break(solver, ValidationPerformance(test_net), every_n_iter=500)

solve(solver, net)

shutdown(backend)


text:
# Linear Models
@testset "Linear Models" begin
    @testset "Balanced Panel Data" begin
        # Balanced Panel Data
        data = dataset("Ecdat", "Crime")
        # Pooling
        model = fit(EconometricModel, @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris), data)
        @test sprint(show, model) ==
            "Continuous Response Model\nNumber of observations: 630\nNull Loglikelihood: 1633.30\nLoglikelihood: 1643.30\nR-squared: 0.0312\nLR Test: 19.99 ∼ χ²(3) ⟹  Pr > χ² = 0.0002\nFormula: CRMRTE ~ 1 + PrbConv + AvgSen + PrbPris\nVariance Covariance Estimator: OIM\n──────────────────────────────────────────────────────────────────────────────────────\n                    PE           SE      t-value  Pr > |t|         2.50%        97.50%\n──────────────────────────────────────────────────────────────────────────────────────\n(Intercept)   0.0186413    0.00431066    4.32447    <1e-04   0.0101762     0.0271064\nPrbConv      -0.00116473   0.000422062  -2.75961    0.0060  -0.00199356   -0.0003359\nAvgSen        0.000236185  0.000268216   0.88058    0.3789  -0.000290526   0.000762897\nPrbPris       0.0273394    0.00817642    3.34369    0.0009   0.0112829     0.0433959\n──────────────────────────────────────────────────────────────────────────────────────"
        @test nullloglikelihood(model) ≈ 1633.305 rtol = 1e-3
        @test loglikelihood(model) ≈ 1643.304 rtol = 1e-3
        @test dof(model) == 5
        # This corresponds to BIC from R's lm, Stata uses rank(vcov(model))
        @test bic(model) ≈ -3254.379 rtol = 1e-3
        @test coef(model) ≈ [0.0186413, -0.0011647, 0.0002362, 0.0273394] rtol = 1e-5
        @test vcov(model) ≈ [
            1.858177e-05 -1.620480e-07 -6.464872e-07 -2.860941e-05
            -1.620480e-07 1.781365e-07 -1.715589e-09 1.286486e-07
            -6.464872e-07 -1.715589e-09 7.193971e-08 8.182148e-09
            -2.860941e-05 1.286486e-07 8.182148e-09 6.685387e-05
            ] rtol = 1e-6
        @test confint(model) ≈ [
            0.0101761950 0.0271063941
            -0.0019935581 -0.0003358998
            -0.0002905263 0.0007628970
            0.0112828681 0.0433959412
            ] rtol = 1e-6
        # Between
        model = fit(
            BetweenEstimator,
            @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris),
            data,
            panel = :County,
            )
        @test sprint(show, model) ==
            "Between Estimator\nCounty with 630 groups\nBalanced groups with size 7\nNumber of observations: 90\nNull Loglikelihood: 239.56\nLoglikelihood: 244.40\nR-squared: 0.1029\nWald: 3.29 ∼ F(3, 86) ⟹ Pr > F = 0.0246\nFormula: CRMRTE ~ 1 + PrbConv + AvgSen + PrbPris\nVariance Covariance Estimator: OIM\n───────────────────────────────────────────────────────────────────────────────────\n                    PE          SE       t-value  Pr > |t|        2.50%      97.50%\n───────────────────────────────────────────────────────────────────────────────────\n(Intercept)  -0.00464339   0.0186571   -0.24888     0.8040  -0.0417325   0.0324457\nPrbConv      -0.00354113   0.0018904   -1.87322     0.0644  -0.00729911  0.00021685\nAvgSen        0.000848049  0.00116477   0.728086    0.4685  -0.00146743  0.00316353\nPrbPris       0.0730299    0.0332057    2.19932     0.0305   0.0070191   0.139041\n───────────────────────────────────────────────────────────────────────────────────"
        @test nullloglikelihood(model) ≈ 239.5638 rtol = 1e-3
        @test loglikelihood(model) ≈ 244.4479 rtol = 1e-3
        @test dof(model) == 5
        @test coef(model) ≈ [-0.00464339, -0.00354113, 0.00084805, 0.07302994] rtol = 1e-6
        @test vcov(model) ≈ [
            3.480883e-04 -5.985740e-06 -1.382250e-05 -5.104706e-04
            -5.985740e-06 3.573593e-06 -1.212431e-08 8.538922e-06
            -1.382250e-05 -1.212431e-08 1.356678e-06 3.953835e-06
            -5.104706e-04 8.538922e-06 3.953835e-06 1.102622e-03
            ] rtol = 1e-6
        # Using Stata's value, R's plm has confint incorrect
        @test confint(model) ≈ [
            -0.0417325 0.0324457
            -0.0072991 0.0002169
            -0.0014674 0.0031635
            0.0070191 0.1390408
            ] rtol = 1e-5
        # Within
        model = fit(
            EconometricModel,
            @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris + absorb(County)),
            data,
            )
        @test sprint(show, model) ==
            "Continuous Response Model\nNumber of observations: 630\nNull Loglikelihood: 1633.30\nLoglikelihood: 2273.95\nR-squared: 0.8707\nWald: 0.16 ∼ F(3, 537) ⟹ Pr > F = 0.9258\nFormula: CRMRTE ~ 1 + PrbConv + AvgSen + PrbPris + absorb(County)\nVariance Covariance Estimator: OIM\n──────────────────────────────────────────────────────────────────────────────────────\n                   PE           SE        t-value  Pr > |t|         2.50%       97.50%\n──────────────────────────────────────────────────────────────────────────────────────\n(Intercept)   0.0314527   0.00204638   15.3699       <1e-43   0.0274328    0.0354726\nPrbConv       6.65981e-6  0.0001985     0.0335507    0.9732  -0.000383272  0.000396592\nAvgSen        7.83181e-5  0.000127904   0.612318     0.5406  -0.000172936  0.000329572\nPrbPris      -0.0013419   0.00405182   -0.331185     0.7406  -0.00930126   0.00661746\n──────────────────────────────────────────────────────────────────────────────────────"
        model = fit(
            EconometricModel,
            @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris + absorb(County + Year)),
            data,
            )
        @test sprint(show, model) ==
            "Continuous Response Model\nNumber of observations: 630\nNull Loglikelihood: 1633.30\nLoglikelihood: 2290.25\nR-squared: 0.8775\nWald: 0.23 ∼ F(3, 531) ⟹ Pr > F = 0.8767\nFormula: CRMRTE ~ 1 + PrbConv + AvgSen + PrbPris + absorb(County + Year)\nVariance Covariance Estimator: OIM\n──────────────────────────────────────────────────────────────────────────────────────\n                    PE           SE       t-value  Pr > |t|         2.50%       97.50%\n──────────────────────────────────────────────────────────────────────────────────────\n(Intercept)   0.0319651    0.00206378   15.4886      <1e-44   0.027911     0.0360193\nPrbConv       7.90362e-5   0.000196089   0.403062    0.6871  -0.00030617   0.000464242\nAvgSen       -9.4788e-5    0.000134924  -0.702531    0.4827  -0.000359838  0.000170262\nPrbPris       0.000979533  0.0040432     0.242267    0.8087  -0.00696309   0.00892216\n──────────────────────────────────────────────────────────────────────────────────────"
        # Random Effects
        model = fit(
            RandomEffectsEstimator,
            @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris),
            data,
            panel = :County,
            time = :Year,
            )
        @test sprint(show, model) ==
            "One-way Random Effect Model\nLongitudinal dataset: County, Year\nBalanced dataset with 90 panels of length 7\nindividual error component: 0.0162\nidiosyncratic error component: 0.0071\nρ: 0.8399\nNumber of observations: 630\nNull Loglikelihood: 2225.79\nLoglikelihood: 2226.01\nR-squared: 0.0007\nWald: 0.15 ∼ F(3, 626) ⟹ Pr > F = 0.9291\nFormula: CRMRTE ~ PrbConv + AvgSen + PrbPris\nVariance Covariance Estimator: OIM\n──────────────────────────────────────────────────────────────────────────────────────\n                    PE           SE       t-value  Pr > |t|         2.50%       97.50%\n──────────────────────────────────────────────────────────────────────────────────────\n(Intercept)   0.0309293    0.00266706   11.5968      <1e-27   0.0256918    0.0361668\nPrbConv      -3.96634e-5   0.000198471  -0.199845    0.8417  -0.000429413  0.000350086\nAvgSen        8.2428e-5    0.000127818   0.644887    0.5192  -0.000168575  0.000333431\nPrbPris      -0.000123327  0.00404272   -0.030506    0.9757  -0.00806226   0.00781561\n──────────────────────────────────────────────────────────────────────────────────────"
        # IV Pooling
        model = fit(EconometricModel, @formula(CRMRTE ~ PrbConv + (AvgSen ~ PrbPris)), data)
        @test sprint(show, model) ==
            "Continuous Response Model\nNumber of observations: 630\nNull Loglikelihood: 1633.30\nLoglikelihood: -611.40\nR-squared: NaN\nLR Test: -4489.40 ∼ χ²(2) ⟹  Pr > χ² = 1.0000\nFormula: CRMRTE ~ 1 + PrbConv + (AvgSen ~ PrbPris)\nVariance Covariance Estimator: OIM\n─────────────────────────────────────────────────────────────────────────────────\n                   PE          SE        t-value  Pr > |t|       2.50%     97.50%\n─────────────────────────────────────────────────────────────────────────────────\n(Intercept)   2.17878     23.0244      0.0946291    0.9246  -43.0357    47.3932\nPrbConv       0.00456765   0.0638118   0.0715801    0.9430   -0.120743   0.129879\nAvgSen       -0.240139     2.57602    -0.093221     0.9258   -5.29883    4.81855\n─────────────────────────────────────────────────────────────────────────────────"
        # IV Between
        model = fit(
            BetweenEstimator,
            @formula(CRMRTE ~ PrbConv + (AvgSen ~ PrbPris)),
            data,
            panel = :County,
            )
        @test sprint(show, model) ==
            "Between Estimator\nCounty with 630 groups\nBalanced groups with size 7\nNumber of observations: 90\nNull Loglikelihood: 239.56\nLoglikelihood: 161.00\nR-squared: NaN\nWald: 0.73 ∼ F(2, 86) ⟹ Pr > F = 0.4844\nFormula: CRMRTE ~ 1 + PrbConv + (AvgSen ~ PrbPris)\nVariance Covariance Estimator: OIM\n─────────────────────────────────────────────────────────────────────────────────\n                   PE          SE       t-value  Pr > |t|       2.50%      97.50%\n─────────────────────────────────────────────────────────────────────────────────\n(Intercept)   0.250667    0.255826     0.979834    0.3299  -0.257898   0.759233\nPrbConv      -0.00331719  0.00481737  -0.688589    0.4929  -0.0128938  0.00625942\nAvgSen       -0.0242107   0.0286331   -0.845549    0.4002  -0.0811314  0.03271\n─────────────────────────────────────────────────────────────────────────────────"
        # IV Within
        model = fit(
            EconometricModel,
            @formula(CRMRTE ~ PrbConv + (AvgSen ~ PrbPris) + absorb(County)),
            data,
            )
        @test sprint(show, model) ==
            "Continuous Response Model\nNumber of observations: 630\nNull Loglikelihood: 1633.30\nLoglikelihood: 2245.89\nR-squared: NaN\nWald: 0.04 ∼ F(2, 537) ⟹ Pr > F = 0.9583\nFormula: CRMRTE ~ 1 + PrbConv + (AvgSen ~ PrbPris) + absorb(County)\nVariance Covariance Estimator: OIM\n──────────────────────────────────────────────────────────────────────────────────────\n                    PE           SE       t-value  Pr > |t|         2.50%       97.50%\n──────────────────────────────────────────────────────────────────────────────────────\n(Intercept)   0.0389703    0.025508      1.52777     0.1272  -0.0111374    0.089078\nPrbConv       2.46646e-5   0.000215805   0.114291    0.9090  -0.000399261  0.000448591\nAvgSen       -0.000826364  0.00285293   -0.289654    0.7722  -0.00643064   0.00477791\n──────────────────────────────────────────────────────────────────────────────────────"
        model = fit(
            EconometricModel,
            @formula(CRMRTE ~ PrbConv + (AvgSen ~ PrbPris) + absorb(County + Year)),
            data,
            )
        @test sprint(show, model) ==
            "Continuous Response Model\nNumber of observations: 630\nNull Loglikelihood: 1633.30\nLoglikelihood: 2280.31\nR-squared: NaN\nWald: 0.09 ∼ F(2, 531) ⟹ Pr > F = 0.9116\nFormula: CRMRTE ~ 1 + PrbConv + (AvgSen ~ PrbPris) + absorb(County + Year)\nVariance Covariance Estimator: OIM\n────────────────────────────────────────────────────────────────────────────────────\n                   PE           SE      t-value  Pr > |t|         2.50%       97.50%\n────────────────────────────────────────────────────────────────────────────────────\n(Intercept)  0.027409     0.0208187    1.31656     0.1886  -0.0134881    0.0683061\nPrbConv      6.00404e-5   0.000214966  0.279301    0.7801  -0.000362248  0.000482329\nAvgSen       0.000462028  0.0023309    0.198219    0.8429  -0.00411688   0.00504094\n────────────────────────────────────────────────────────────────────────────────────"
        # IV Random
        model = fit(
            RandomEffectsEstimator,
            @formula(CRMRTE ~ PrbConv + (AvgSen ~ PrbPris)),
            data,
            panel = :County,
            time = :Year,
            )
        @test sprint(show, model) ==
            "One-way Random Effect Model\nLongitudinal dataset: County, Year\nBalanced dataset with 90 panels of length 7\nindividual error component: 0.0413\nidiosyncratic error component: 0.0074\nρ: 0.9691\nNumber of observations: 630\nNull Loglikelihood: 2268.01\nLoglikelihood: 2248.34\nR-squared: NaN\nWald: 0.03 ∼ F(2, 626) ⟹ Pr > F = 0.9671\nFormula: CRMRTE ~ PrbConv + (AvgSen ~ PrbPris)\nVariance Covariance Estimator: OIM\n───────────────────────────────────────────────────────────────────────────────────────\n                    PE           SE        t-value  Pr > |t|         2.50%       97.50%\n───────────────────────────────────────────────────────────────────────────────────────\n(Intercept)   0.037747     0.0241684     1.56183      0.1188  -0.00971405   0.085208\nPrbConv       1.39581e-5   0.000200244   0.0697054    0.9445  -0.000379273  0.000407189\nAvgSen       -0.000688924  0.00266514   -0.258494     0.7961  -0.00592262   0.00454478\n───────────────────────────────────────────────────────────────────────────────────────"
    end
    @testset "Heteroscedasticity Consistent Variance Covariance Estimators" begin
        data = dataset("Ecdat", "Crime")
        # Pooling
        model = fit(EconometricModel, @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris), data)
        @test vcov(model) ≈ vcov(model, OIM)
        @test vcov(model, HC0) ≈ [
            1.51708E-05 3.58518E-07 -7.30083E-07 -2.03240E-05
            3.58518E-07 2.85571E-07 -3.58783E-08 -4.13434E-07
            -7.30083E-07 -3.58783E-08 9.36232E-08 -1.46277E-07
            -2.03240E-05 -4.13434E-07 -1.46277E-07 5.17390E-05
            ] rtol = 1e-6
        @test vcov(model, HC1) ≈ [
            1.52677E-05 3.60808E-07 -7.34748E-07 -2.04539E-05
            3.60808E-07 2.87396E-07 -3.61076E-08 -4.16076E-07
            -7.34748E-07 -3.61076E-08 9.42214E-08 -1.47212E-07
            -2.04539E-05 -4.16076E-07 -1.47212E-07 5.20696E-05
            ] rtol = 1e-6
        @test vcov(model, HC2) ≈ [
            1.56208E-05 9.05631E-08 -7.36000E-07 -2.08761E-05
            9.05631E-08 6.78782E-07 -5.03354E-08 -3.93296E-08
            -7.36000E-07 -5.03354E-08 9.59715E-08 -1.59383E-07
            -2.08761E-05 -3.93296E-08 -1.59383E-07 5.27586E-05
            ] rtol = 1e-6
        @test vcov(model, HC3) ≈ [
            1.66839E-05 -9.88963E-07 -7.13420E-07 -2.23048E-05
            -9.88963E-07 2.16675E-06 -1.04165E-07 1.50077E-06
            -7.13420E-07 -1.04165E-07 9.98145E-08 -2.14382E-07
            -2.23048E-05 1.50077E-06 -2.14382E-07 5.50463E-05
            ] rtol = 1e-6
        @test vcov(model, HC4) ≈ [
            3.21597E-05 -2.14119E-05 -4.18897E-09 -4.43245E-05
            -2.14119E-05 2.97197E-05 -1.09486E-06 3.08355E-05
            -4.18897E-09 -1.09486E-06 1.38192E-07 -1.26498E-06
            -4.43245E-05 3.08355E-05 -1.26498E-06 8.67924E-05
            ] rtol = 1e-6
        dv = fit(
            EconometricModel,
            @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris + County),
            data,
            contrasts = Dict(:County => DummyCoding()),
            )
        fe = fit(
            EconometricModel,
            @formula(CRMRTE ~ PrbConv + AvgSen + PrbPris + absorb(County)),
            data,
            )
        @test vcov(fe)[2:end, 2:end] ≈ vcov(dv)[2:4, 2:4]
        @test vcov(fe, HC1)[2:end, 2:end] ≈ vcov(dv, HC1)[2:4, 2:4]
    end
end


text:
abstract type Node end
struct Document <:Node
	kind::String
  definitions::Array{Node,1}
  Document(d)=new("Document",d)
end
struct OperationDefinition <:Node
	kind::String
	operation::String
	name
	variableDefinitions
	directives
	selectionSet
	OperationDefinition(op,name,vardef,direc,selec)=new("OperationDefinition",op,name,vardef,direc,selec)
end
struct SelectionSet <:Node
	kind::String
	selections::Array{Node,1}
	SelectionSet(sel) = new("SelectionSet",sel)
end
struct Name <:Node
	kind::String
	value::String
	Name(val)= new("Name",val)
end
struct Argument <:Node
	kind::String# "Argument";
	name::Name
	value #[ Literal | Variable | Reference ];
	Argument(n::Name,x)= new("Argument",n,x)
end
struct VariableDefinition <:Node
	kind::String
	variable
	tipe
	defaultValue
	VariableDefinition(va,tipe,defval)=new("VariableDefinition",va,tipe,defval)
end
struct Variable <:Node
	kind::String# "Variable";
	name
	Variable(n)= new("Variable",n)
end
struct FragmentSpread <:Node
	kind::String
	name
	directives
	FragmentSpread(n,d)=new("FragmentSpread",n,d)
end
struct InlineFragment <:Node
	kind::String
	typeCondition
	directives
  selectionSet
  InlineFragment(typecon,direc,selset)=InlineFragment("InlineFragment",typecon,direc,selset)
end
struct FragmentDefinition <:Node
	kind::String
	name
	typeCondition
	directives
	selectionSet
	FragmentDefinition(n,t,d,s)=new("FragmentDefinition",n,t,d,s)
end
struct IntValue <:Node
	kind::String
	value::Int64
	IntValue(v)=new("IntValue",v)
end
struct FloatValue <:Node
	kind::String
	value::Float64
	FloatValue(v)=new("FloatValue",v)
end
struct StringValue <:Node
	kind::String
	value::String
	StringValue(v)=new("StringValue",v)
end
struct BooleanValue <:Node
	kind::String
	value::Bool
	BooleanValue(v)=new("BooleanValue",v)
end
struct NullValue <:Node
	kind::String
	NullValue()=new("NullValue")
end
struct EnumValue <:Node
	kind::String
	value
	EnumValue(v)=new("EnumValue",v)
end
struct List_ <:Node
	kind::String
	value_s
	List_(vals)=new("List",vals)
end
struct Object_ <:Node
	kind::String
	fields
	Object_(f)=new("Object",f)
end
struct Object_Field <:Node
	kind::String
	name
	value
	Object_Field(n,v)=new("Object_Field",n,v)
end
struct Directive <:Node
	kind::String
	name
	arguments
	Directive(n,a) = new("Directive",n,a)
end
struct ListType <:Node
    kind::String
    tipe
    ListType(t)=new("ListType",t)
end
struct NonNullType <:Node
    kind::String
    tipe
    NonNullType(t)=new("NonNullType",t)
end
struct NamedType <:Node
   kind::String
   name
   NamedType(n)=new("NamedType",n)
end
struct SchemaDefinition <:Node
	kind::String
  directives
  operationTypes::Array{Node,1}
  SchemaDefinition(d,o)=new("SchemaDefinition",d,o)
end
struct OperationTypeDefinition <:Node
  kind::String
  operation
  tipe
  OperationTypeDefinition(o,t)=new("OperationTypeDefinition",o,t)
end
struct ScalarTypeDefinition <:Node
  kind::String
  name
  directives
  ScalarTypeDefinition(n,d)=new("ScalarTypeDefinition",n,d)
end
struct ObjectTypeDefinition <:Node
	kind::String
	name
	interfaces::Array{Node,1}
	directives
	fields::Array{Node,1}
	ObjectTypeDefinition(n,i,d,f)=new("ObjectTypeDefinition",n,i,d,f)
end
struct FieldDefinition <:Node
	kind::String# "Field";
	name
	arguments
	tipe
  directives
	FieldDefinition(n,arg,ti,dir)= new("FieldDefinition",n,arg,ti,dir)
end
struct Field <:Node
  kind::String
  alias
  name
  arguments
  directives
  selectionSet
  Field(a,n,arg,d,s)=new("Field",a,n,arg,d,s)
end
struct InputValueDefinition <:Node
  kind::String
  name
  tipe
  defaultValue
  directives
  InputValueDefinition(n,t,dfv,dir)=new("InputValueDefinition",n,t,dfv,dir)
end
struct InterfaceTypeDefinition <:Node
  kind::String
  name
  directives
  fields::Array{Node,1}
  InterfaceTypeDefinition(n,d,f)=new("InterfaceTypeDefinition",n,d,f)
end
struct UnionTypeDefinition <:Node
	kind::String
  name
  directives
  tipes::Array{Node,1}
  UnionTypeDefinition(n,d,t)=new("UnionTypeDefinition",n,d,t)
end
struct EnumTypeDefinition <:Node
  kind::String
  name
  directives
  _values::Array{Node,1}
  EnumTypeDefinition(n,d,v)=new("EnumTypeDefinition",n,d,v)
end
struct EnumValueDefinition <:Node
  kind::String
  name
  directives
  EnumValueDefinition(n,d)=new("EnumValueDefinition",n,d)
end
struct InputObjectTypeDefinition <:Node
  kind::String
  name
  directives
  fields
  InputObjectTypeDefinition(n,d,f)=new("InputObjectTypeDefinition",n,d,f)
end
struct TypeExtensionDefinition <:Node
  kind::String
  definition
  TypeExtensionDefinition(d)=new("TypeExtensionDefinition",d)
end
struct DirectiveDefinition <:Node
  kind::String
  name
  arguments
  locations
  DirectiveDefinition()=new("DirectiveDefinition")
end
struct Lexer
    token
    advance
    back
    function Lexer(str)
        lexers =Tokensgraphql(str)

        indice = 1
        fin =  length(lexers)

        function advance()
            if indice<=fin
                indice = indice+1
            end
        end

        function back()
            if indice>1
                indice = indice-1
            end
        end

        function current_token()
            return lexers[indice]
        end
        new(current_token,advance,back)
    end
end

"""
  Determines if the next token is of a given kind
"""
function peek(lexer::Lexer, kind)
  return lexer.token().kind === kind
end

"""
  If the next token is of the given kind, return true after advancing
  the lexer. Otherwise, do not change the parser state and return false.
"""
function next_token(lexer::Lexer, kind)
	check= lexer.token().kind === kind
	if (check)
		lexer.advance()
	end
	return check
end

"""
  If the next token is of the given kind, return that token after advancing
  the lexer. Otherwise, do not change the parser state and throw an error.
"""
function expect(lexer, kind)
	token = lexer.token()

	if (token.kind === kind)
		lexer.advance()
		return token
	end
	return throw(GraphQLError("{\"errors\":[{\"locations\": [{\"column\": $(token.startpos[2]),\"line\": $(token.startpos[1])}],\"message\": \"Syntax Error GraphQL request ($(token.startpos[2]):$(token.startpos[1])) Expected $kind, found $token \"}]}"))
end

"""
  If the next token is a keyword with the given value, return that token after
  advancing the lexer. Otherwise, do not change the parser state and return
  false.
"""
function expectKeyword(lexer, value::String)
  token = lexer.token()
  if ((token.kind === Tokens.NAME) & (token.val == value))
    lexer.advance()
    return token
  end
  return throw(GraphQLError("Expected $value, found $token in: line $(token.startpos[1]) , col $(token.startpos[2])"))
end

"""
  Helper function for creating an error when an unexpected lexed token
  is encountered.
"""
function unexpected(token)
  return "{\"errors\":[{\"locations\": [{\"column\": $(token.startpos[2]),\"line\": $(token.startpos[1])}],\"message\": \"Syntax Error GraphQL request ($(token.startpos[2]):$(token.startpos[1])) Unexpected character $(token) \"}]}"
end
function unexpected(lexer::Lexer)
  token =lexer.token()
  return "{\"errors\":[{\"locations\": [{\"column\": $(token.startpos[2]),\"line\": $(token.startpos[1])}],\"message\": \"Syntax Error GraphQL request ($(token.startpos[2]):$(token.startpos[1])) Unexpected character $(token) \"}]}"

end

"""
  Returns a possibly empty list of parse nodes, determined by
  the parseFn. This list begins with a lex token of openKind
  and ends with a lex token of closeKind. Advances the parser
  to the next lex token after the closing token.
"""
function an_y(lexer::Lexer, openKind, parseFn::Function, closeKind)
  expect(lexer, openKind)
  nodes = []
  while (!next_token(lexer, closeKind))
    push!(nodes,parseFn(lexer))
  end
  return nodes
end

"""
  Returns a non-empty list of parse nodes, determined by
  the parseFn. This list begins with a lex token of openKind
  and ends with a lex token of closeKind. Advances the parser
  to the next lex token after the closing token.
"""
function many(lexer::Lexer, openKind, parseFn::Function, closeKind)
 	expect(lexer, openKind)
 	nodes = [ parseFn(lexer) ]
 	while (!next_token(lexer, closeKind))
 		push!(nodes,parseFn(lexer))
 	end
 	return nodes
 end
"""
    Parse(str::String)

Given a GraphQL source, parses it into a AST.
Throws Diana.GraphQLError if a syntax error is encountered.
 """
function Parse(str::String)
	return parseDocument(Lexer(str))
end

 function parseDocument(lexer::Lexer)
 	start_token = lexer.token()
 	definitions = []
 	while true
 		push!(definitions,parseDefinition(lexer))
 		if(next_token(lexer, Tokens.ENDMARKER))
 			break
 		end
 	end
 	return Document(definitions)
 end

"""
  Definition :
    - OperationDefinition
    - FragmentDefinition
    - TypeSystemDefinition
"""
function parseDefinition(lexer::Lexer)
  if (peek(lexer, Tokens.LBRACE))
    return parseOperationDefinition(lexer)
  end

  if (peek(lexer, Tokens.NAME))
      # Note: subscription is an experimental non-spec addition.
      if((lexer.token().val)=="query")
        return parseOperationDefinition(lexer)
      end
      if((lexer.token().val)=="mutation")
        return parseOperationDefinition(lexer)
      end
      if((lexer.token().val)=="subscription")
        return parseOperationDefinition(lexer)
      end

      if((lexer.token().val)=="fragment")
      	return parseFragmentDefinition(lexer)
      end

      # Note: the Type System IDL is an experimental non-spec addition.
      if((lexer.token().val)=="schema")
        return parseTypeSystemDefinition(lexer) end
      if((lexer.token().val)=="scalar")
      return parseTypeSystemDefinition(lexer)
       end
      if((lexer.token().val)=="type")
      return parseTypeSystemDefinition(lexer)
       end
      if((lexer.token().val)=="interface")
      return parseTypeSystemDefinition(lexer)
       end
      if((lexer.token().val)=="union")
      return parseTypeSystemDefinition(lexer)
       end
      if((lexer.token().val)=="enum")
      return parseTypeSystemDefinition(lexer)
       end
      if((lexer.token().val)=="input")
      return parseTypeSystemDefinition(lexer)
       end
      if((lexer.token().val)=="extend")
      return parseTypeSystemDefinition(lexer)
       end
      if((lexer.token().val)=="directive")
      	return parseTypeSystemDefinition(lexer)
      end
    end

  return throw(GraphQLError(unexpected(lexer)))
end

"""
	Implements the parsing rules in the Operations section.

  OperationDefinition :
   - SelectionSet
   - OperationType Name? VariableDefinitions? Directives? SelectionSet
"""
function parseOperationDefinition(lexer::Lexer)
  start_token = lexer.token()
  if (peek(lexer, Tokens.LBRACE))
    pSS=parseSelectionSet(lexer)
  	return OperationDefinition("query",nothing,nothing,[],pSS)
  end

  operation = parseOperationType(lexer)
  name = ""
  if (peek(lexer, Tokens.NAME))
    name = parseName(lexer)
  end
  pVD=parseVariableDefinitions(lexer)
  pD=parseDirectives(lexer)
  pS=parseSelectionSet(lexer)
  return OperationDefinition(operation,name,pVD,pD,pS)
end

"""
  SelectionSet : { Selection+ }
"""
function parseSelectionSet(lexer::Lexer)
  start_token= lexer.token()
  return SelectionSet(many(lexer, Tokens.LBRACE, parseSelection,Tokens.RBRACE))
end

"""
  Selection :
    - Field
    - FragmentSpread
    - InlineFragment
"""
function parseSelection(lexer::Lexer)
  return peek(lexer, Tokens.SPREAD) ? parseFragment(lexer) : parseField(lexer)
end

"""
  Field : Alias? Name Arguments? Directives? SelectionSet?

  Alias : Name : NAME
"""
function parseField(lexer::Lexer)
  start_token = lexer.token()

  nameOrAlias = parseName(lexer)
  alias=""
  name=""
  if (next_token(lexer, Tokens.COLON))
    alias = nameOrAlias
    name = parseName(lexer)
  else
    alias = nothing
    name = nameOrAlias
  end
  pA= parseArguments(lexer)
  pD= parseDirectives(lexer)
  pSS=(peek(lexer, Tokens.LBRACE) ? parseSelectionSet(lexer) : nothing)
  return Field(alias,name,pA,pD,pSS#=,loc(lexer, start_token) =#)
end

"""
  Converts a name lex token into a name parse node.
"""
function parseName(lexer::Lexer)
  token = expect(lexer, Tokens.NAME)
  return Name(token.val)
end

"""
 Arguments : ( Argument+ )
"""
function parseArguments(lexer::Lexer)
  return peek(lexer, Tokens.LPAREN) ? many(lexer, Tokens.LPAREN, parseArgument, Tokens.RPAREN) : []
end

"""
  Argument : Name : Value
"""
function parseArgument(lexer::Lexer)
  start_token = lexer.token()
  pN=parseName(lexer)
  return Argument(pN,(expect(lexer, Tokens.COLON).val, parseValueLiteral(lexer, false)))
end

"""
  Given a string containing a GraphQL value (ex. `[42]`), parse the AST for
  that value.
  Throws Diana.GraphQLError if a syntax error is encountered.

  This is useful within tools that operate upon GraphQL Values directly and
  in isolation of complete GraphQL documents.

  Consider providing the results to the utility function: valueFromAST().
"""
function parseValue(source::String)   ##export
  lexer = Lexer(source)
  value = parseValueLiteral(lexer, false)
  expect(lexer, Tokens.ENDMARKER)
  return value
end

"""
  Given a string containing a GraphQL Type (ex. `[Int!]`), parse the AST for
  that type.
  Throws Diana.GraphQLError if a syntax error is encountered.

  This is useful within tools that operate upon GraphQL Types directly and
  in isolation of complete GraphQL documents.

  Consider providing the results to the utility function: typeFromAST().
"""
function parseType(source::String) ##export
  sourceObj = source
  lexer = createLexer(sourceObj)
  tipe = parseTypeReference(lexer)
  expect(lexer, Tokens.ENDMARKER)
  return tipe
end

"""
 Implements the parsing rules in the Document section.

  OperationType : one of query mutation subscription
"""
function parseOperationType(lexer::Lexer)
	operationToken = expect(lexer, Tokens.NAME)
	if(operationToken.val =="query")
		return "query"
	end
	if(operationToken.val =="mutation")
		return "mutation"
	end
	# Note: subscription is an experimental non-spec addition.
	if(operationToken.val =="subscription")
		return "subscription"
	end

	return throw(GraphQLError(unexpected(operationToken)))
end

"""
  VariableDefinitions : ( VariableDefinition+ )
"""
function parseVariableDefinitions(lexer::Lexer)
  return peek(lexer, Tokens.LPAREN) ? many(lexer,Tokens.LPAREN, parseVariableDefinition,Tokens.RPAREN) : []
end

"""
  VariableDefinition : Variable : Type DefaultValue?
"""
function parseVariableDefinition(lexer::Lexer)
  start_token = lexer.token()
  return VariableDefinition(parseVariable(lexer),((expect(lexer, Tokens.COLON)).val, parseTypeReference(lexer)),next_token(lexer, Tokens.EQUALS) ? parseValueLiteral(lexer, true) : nothing)
end

"""
  Variable : \$ Name
"""
function parseVariable(lexer::Lexer)
  start_token = lexer.token()
  expect(lexer, Tokens.DOLLAR)
  return Variable(parseName(lexer))
end

"""
Implements the parsing rules in the Fragments section.

  Corresponds to both FragmentSpread and InlineFragment in the spec.

  FragmentSpread : ... FragmentName Directives?

  InlineFragment : ... TypeCondition? Directives? SelectionSet
"""
function parseFragment(lexer::Lexer)
  start_token = lexer.token()
  expect(lexer, Tokens.SPREAD)
  if ((peek(lexer, Tokens.NAME)) & (lexer.token().val !== "on"))
  	return FragmentSpread(parseFragmentName(lexer),parseDirectives(lexer))
  end

  typeCondition = nothing

  if (lexer.token().val == "on")
    lexer.advance()
    typeCondition = parseNamedType(lexer)
  end
  pD=parseDirectives(lexer)
  pS=parseSelectionSet(lexer)
  return InlineFragment(typeCondition,pD,pS)
end

"""
  FragmentDefinition :
    - fragment FragmentName on TypeCondition Directives? SelectionSet

  TypeCondition : NamedType
"""
function parseFragmentDefinition(lexer::Lexer)
  start_token = lexer.token()
  expectKeyword(lexer,"fragment")
  pFN= parseFragmentName(lexer)
  v=((expectKeyword(lexer,"on").val), parseNamedType(lexer))
  pD=parseDirectives(lexer)
  pS=parseSelectionSet(lexer)
  return FragmentDefinition(pFN,v,pD,pS)
end

"""
  FragmentName : Name but not `on`
"""
function parseFragmentName(lexer::Lexer)
  if (lexer.token().val == "on")
    return throw(GraphQLError(unexpected(lexer)))
  end
  return parseName(lexer)
end

"""
Implements the parsing rules in the Values section.

  Value[Const] :
    - [~Const] Variable
    - IntValue
    - FloatValue
    - StringValue
    - BooleanValue
    - NullValue
    - EnumValue
    - ListValue[?Const]
    - ObjectValue[?Const]

  BooleanValue : one of `true` `false`

  NullValue : `null`

  EnumValue : Name but not `true`, `false` or `null`
"""
function parseValueLiteral(lexer::Lexer, isConst::Bool)
  token = lexer.token()
    if(token.kind==Tokens.LSQUARE)
      return parseList(lexer, isConst)
    end

    if(token.kind==Tokens.LBRACE)
      return parseObject(lexer, isConst)
    end

    if(token.kind==Tokens.INT)
      lexer.advance()
      return IntValue(token.val)
    end

    if(token.kind==Tokens.FLOAT)
      lexer.advance()
      return FloatValue(token.val)
    end

    if(token.kind==Tokens.STRING)
      lexer.advance()
      return StringValue(token.val)
    end

    if(token.kind==Tokens.NAME)
      if((token.val == "true" )| (token.val == "false"))
        lexer.advance()
        return BooleanValue(token.val == "true")
      elseif(token.val == "null")
        lexer.advance()
        return NullValue(#=loc(lexer, token)=#)
      end

      lexer.advance()
      return EnumValue(token.val)
    end

    if(token.kind==Tokens.DOLLAR)
      if (!isConst)
        return parseVariable(lexer)
      end
      #break
    end

  return throw(GraphQLError(unexpected(lexer)))
end

function parseConstValue(lexer::Lexer) ##export
  return parseValueLiteral(lexer, true)
end

function parseValueValue(lexer::Lexer)
  return parseValueLiteral(lexer, false)
end

"""
  ListValue[Const] :
    - [ ]
    - [ Value[?Const]+ ]
"""
function parseList(lexer::Lexer, isConst::Bool)
  start_token = lexer.token()
  item = isConst ? parseConstValue : parseValueValue
  return List_(an_y(lexer, Tokens.LSQUARE, item, Tokens.RSQUARE))
end

"""
  ObjectValue[Const] :
    - { }
    - { ObjectField[?Const]+ }
"""
function parseObject(lexer::Lexer, isConst::Bool)
  start_token = lexer.token()
  expect(lexer, Tokens.LBRACE)
  fields = []
  while (!next_token(lexer, Tokens.RBRACE))
    push!(fields,parseObjectField(lexer, isConst))
  end
  return Object_(fields)
end

"""
  ObjectField[Const] : Name : Value[?Const]
"""
function parseObjectField(lexer::Lexer, isConst::Bool)
  start_token = lexer.token()
  return Object_Field(parseName(lexer),((expect(lexer, Tokens.COLON)).val, parseValueLiteral(lexer, isConst)))
end

"""
Implements the parsing rules in the Directives section.

  Directives : Directive+
"""
function parseDirectives(lexer::Lexer)
  directives = []
  while (peek(lexer,Tokens.AT))
    push!(directives,parseDirective(lexer))
  end
  return directives
end

"""
  Directive : @ Name Arguments?
"""
function parseDirective(lexer::Lexer)
  start_token = lexer.token()
  expect(lexer, Tokens.AT)
  return Directive(parseName(lexer),parseArguments(lexer))
end

"""
Implements the parsing rules in the Types section.

  Type :
    - NamedType
    - ListType
    - NonNullType
"""
function parseTypeReference(lexer::Lexer) ##export
  start_token = lexer.token()
  tipe=""
  if(next_token(lexer, Tokens.LSQUARE))
    tipe = parseTypeReference(lexer)
    expect(lexer, Tokens.RSQUARE)
    tipe = ListType(tipe)
  else
    tipe = parseNamedType(lexer)
  end
  if (next_token(lexer, Tokens.BANG))
  	return NonNullType(tipe)
  end
  return tipe
end

"""
  NamedType : Name
"""
function parseNamedType(lexer::Lexer)  ##export
  start_token = lexer.token()
  return NamedType(parseName(lexer))
end

"""
Implements the parsing rules in the Type Definition section.

  TypeSystemDefinition :
    - SchemaDefinition
    - TypeDefinition
    - TypeExtensionDefinition
    - DirectiveDefinition

  TypeDefinition :
    - ScalarTypeDefinition
    - ObjectTypeDefinition
    - InterfaceTypeDefinition
    - UnionTypeDefinition
    - EnumTypeDefinition
    - InputObjectTypeDefinition
"""
function parseTypeSystemDefinition(lexer::Lexer)
  if (peek(lexer, Tokens.NAME))
      if(lexer.token().val=="schema") return parseSchemaDefinition(lexer) end
      if(lexer.token().val=="scalar") return parseScalarTypeDefinition(lexer) end
      if(lexer.token().val=="type") return parseObjectTypeDefinition(lexer) end
      if(lexer.token().val=="interface") return parseInterfaceTypeDefinition(lexer) end
      if(lexer.token().val=="union") return parseUnionTypeDefinition(lexer) end
      if(lexer.token().val=="enum") return parseEnumTypeDefinition(lexer) end
      if(lexer.token().val=="input") return parseInputObjectTypeDefinition(lexer) end
      if(lexer.token().val=="extend") return parseTypeExtensionDefinition(lexer) end
      if(lexer.token().val=="directive") return parseDirectiveDefinition(lexer) end
  end

  return throw(GraphQLError(unexpected(lexer)))
end

"""
  SchemaDefinition : schema Directives? { OperationTypeDefinition+ }

  OperationTypeDefinition : OperationType : NamedType
"""
function parseSchemaDefinition(lexer::Lexer)
  start_token = lexer.token()
  expectKeyword(lexer,"schema")
  directives = parseDirectives(lexer)
  operationTypes = many(lexer,Tokens.LBRACE,parseOperationTypeDefinition,Tokens.RBRACE)
  return SchemaDefinition(directives,operationTypes)
end

function parseOperationTypeDefinition(lexer::Lexer)
  start_token = lexer.token()
  operation = parseOperationType(lexer)
  expect(lexer, Tokens.COLON)
  tipe = parseNamedType(lexer)
  return OperationTypeDefinition(operation,tipe)
end

"""
  ScalarTypeDefinition : scalar Name Directives?
"""
function parseScalarTypeDefinition(lexer::Lexer)
  start_token = lexer.token()
  expectKeyword(lexer,"scalar")
  name = parseName(lexer)
  directives = parseDirectives(lexer)
  return ScalarTypeDefinition(name,directives)
end

"""
  ObjectTypeDefinition :
    - type Name ImplementsInterfaces? Directives? { FieldDefinition+ }
"""
function parseObjectTypeDefinition(lexer::Lexer)
  start_token = lexer.token()
  expectKeyword(lexer, "type")
  name = parseName(lexer)
  interfaces = parseImplementsInterfaces(lexer)
  directives = parseDirectives(lexer)
  fields = an_y(lexer,Tokens.LBRACE,parseFieldDefinition,Tokens.RBRACE)
  return ObjectTypeDefinition(name,interfaces,directives,fields)
end

"""
  ImplementsInterfaces : implements NamedType+
"""
function parseImplementsInterfaces(lexer::Lexer)
  tipes = []

  try
  expectKeyword(lexer,"implements")
    while true
    	push!(tipes,parseNamedType(lexer))
        if(!peek(lexer, Tokens.NAME))
        	break
        end
    end
  return tipes
  catch
    return []
  end
end

"""
  FieldDefinition : Name ArgumentsDefinition? : Type Directives?
"""
function parseFieldDefinition(lexer::Lexer)
  start_token = lexer.token()
  name = parseName(lexer)
  args = parseArgumentDefs(lexer)
  expect(lexer, Tokens.COLON)
  tipe = parseTypeReference(lexer)
  directives = parseDirectives(lexer)
  return FieldDefinition(name,args,tipe,directives)
end

"""
  ArgumentsDefinition : ( InputValueDefinition+ )
"""
function parseArgumentDefs(lexer::Lexer)
  if (!peek(lexer, Tokens.LPAREN))
    return []
  end
  return many(lexer, Tokens.LPAREN, parseInputValueDef, Tokens.RPAREN)
end

"""
  InputValueDefinition : Name : Type DefaultValue? Directives?
"""
function parseInputValueDef(lexer::Lexer)
  start_token = lexer.token()
  name = parseName(lexer)
  expect(lexer, Tokens.COLON)
  tipe = parseTypeReference(lexer)
  defaultValue = nothing
  if(next_token(lexer, Tokens.EQUALS))
    defaultValue = parseConstValue(lexer)
  end
  directives = parseDirectives(lexer)
  return InputValueDefinition(name,tipe,defaultValue,directives)
end

"""
  InterfaceTypeDefinition : interface Name Directives? { FieldDefinition+ }
"""
function parseInterfaceTypeDefinition(lexer::Lexer)
  start_token = lexer.token()
  expectKeyword(lexer,"interface")
  name = parseName(lexer)
  directives = parseDirectives(lexer)
  fields = an_y(lexer,Tokens.LBRACE,parseFieldDefinition,Tokens.RBRACE)

  return InterfaceTypeDefinition(name,directives,fields)
end

"""
  UnionTypeDefinition : union Name Directives? = UnionMembers
"""
function parseUnionTypeDefinition(lexer::Lexer)
  start_token = lexer.token()
  expectKeyword(lexer, "union")
  name = parseName(lexer)
  directives = parseDirectives(lexer)
  expect(lexer, Tokens.EQUALS)
  tipes = parseUnionMembers(lexer)
  return UnionTypeDefinition(name,directives,tipes,#=loc(lexer, start_token)=#)
end

"""
	UnionMembers :
    - `|`? NamedType
    - UnionMembers | NamedType
"""
function parseUnionMembers(lexer::Lexer)
  # Optional leading pipe
  next_token(lexer, Tokens.PIPE)
  members = []
  while true
  	push!(members,parseNamedType(lexer))
  	if(!next_token(lexer, Tokens.PIPE))
  	    break
  	end
  end

  return members
end

"""
  EnumTypeDefinition : enum Name Directives? { EnumValueDefinition+ }
"""
function parseEnumTypeDefinition(lexer::Lexer)
  start_token = lexer.token()
  expectKeyword(lexer, "enum")
  name = parseName(lexer)
  directives = parseDirectives(lexer)
  _values = many(lexer,Tokens.LBRACE,parseEnumValueDefinition,Tokens.RBRACE)
  return EnumTypeDefinition(name,directives,_values,#=loc(lexer, start_token),=#)
end

"""
  EnumValueDefinition : EnumValue Directives?

  EnumValue : Name
"""
function parseEnumValueDefinition(lexer::Lexer)
  start_token = lexer.token()
  name = parseName(lexer)
  directives = parseDirectives(lexer)
  return EnumValueDefinition(name,directives#=,loc(lexer, start_token),=#)
end

"""
  InputObjectTypeDefinition : input Name Directives? { InputValueDefinition+ }
"""
function parseInputObjectTypeDefinition(lexer::Lexer)
  start_token = lexer.token()
  expectKeyword(lexer, "input")
  name = parseName(lexer)
  directives = parseDirectives(lexer)
  fields = an_y(lexer,Tokens.LBRACE,parseInputValueDef,Tokens.RBRACE)
  return InputObjectTypeDefinition(name,directives,fields)
end

"""
  TypeExtensionDefinition : extend ObjectTypeDefinition
"""
function parseTypeExtensionDefinition(lexer::Lexer)
  start_token = lexer.token()
  expectKeyword(lexer, "extend")
  definition = parseObjectTypeDefinition(lexer)
  TypeExtensionDefinition
  return TypeExtensionDefinition(definition)
end

"""
  DirectiveDefinition :
    - directive @ Name ArgumentsDefinition? on DirectiveLocations
"""
function parseDirectiveDefinition(lexer::Lexer)
  start_token = lexer.token()
  expectKeyword(lexer, "directive")
  expect(lexer, Tokens.AT)
  name = parseName(lexer)
  args = parseArgumentDefs(lexer)
  expectKeyword(lexer, "on")
  locations = parseDirectiveLocations(lexer)
  return DirectiveDefinition(name,args,locations)
end

"""
  DirectiveLocations :
    - `|`? Name
    - DirectiveLocations | Name
"""
function parseDirectiveLocations(lexer::Lexer)
  # Optional leading pipe
  next_token(lexer, Tokens.PIPE)
  locations = []
  while true
  	push!(locations,parseName(lexer))
    if(!next_token(lexer, Tokens.PIPE))
      break
    end
  end
  return locations
end

function Base.show(io::IO, x::Node)
  t = typeof(x)
  print(io, "\n\e[33m < \e[37m")
    for f in fieldnames(t)
      elemento = getfield(x, f)
      v= true
      try
        v = length(elemento)>0
      catch
        v=true
      end
      if ((elemento!==nothing) & v)
        if(f == :kind)
          print(io,"Node :: ",elemento)
        elseif(f == :value) | (f == :operation)
          print(io," ,",f ," : \e[92m",elemento,"\e[37m")
        else
          print(io," ,",f ," : ",elemento)
        end
      end
    end

  print(io,"\e[33m > \e[37m")
end


text:
if abspath(PROGRAM_FILE) == @__FILE__
    a::Int64 = 2
    @assert(~(a) == -3)
    -1
    +1
end


text:
module cohort
    using Random
    include("./fixed_parameters.jl")
    using .fixed_parameters
    include("./human.jl")
    using .human
    include("./class.jl")
    using .class

    mutable struct Cohort
        quarantine::Bool
        warning::Bool
        quarantine_time::Int64
        crew::Array
        houses::Array
        classes::Array

        house_time::Any
        set_Quarantine::Any
        warn::Any

        function Cohort(num_ex_students::Int64, num_ex_teachers::Int64, intern_random_infection_rate::Float64, extern_random_infection_rate::Float64, teacher_break_rules::Float64, student_break_rules::Float64, houses::Array, class_info::Array)
            my_crew = [Array{Any}(undef, num_ex_students), Array{Any}(undef, num_ex_teachers)]
            for i in 1:num_ex_students
                spawnpoint = (rand(d:0.01:lund_size[1]-d), rand(cwh:0.01:lund_size[2]-hwh))
                my_crew[1][i] = Human(spawnpoint, student_break_rules, false, extern_random_infection_rate)
            end
            for i in 1:num_ex_teachers
                spawnpoint = (rand(d:0.01:lund_size[1]-d), rand(cwh:0.01:lund_size[2]-hwh))
                my_crew[2][i] = Human(spawnpoint, teacher_break_rules, false, extern_random_infection_rate)
            end

            for i in houses
                for j in 1:2
                    for k in i.crew[j]
                        push!(my_crew[j], k)
                    end
                end
            end

            my_classes = Array{Any}(undef, length(class_info))
            Random.shuffle!(my_crew[1])
            Random.shuffle!(my_crew[2])
            start = [0,0]
            for (index, info) in enumerate(class_info)
                my_classes[index] = Class(info[1], info[2])
                my_classes[index].pick(my_classes[index], my_crew, start)
                start[1] += my_classes[index].num_students
                start[2] += my_classes[index].num_teachers
            end

            new(false, false, 0, my_crew, houses, my_classes, house_time, set_Quarantine, warn)
        end
    end

    function house_time(self::Cohort, fakewalls::Array)
        for i in self.houses
            i.house_time(i, fakewalls)
        end
        for i in self.crew
            for j in i
                if !j.intern
                    j.position = (nothing, nothing)
                end
            end
        end
    end

    function set_Quarantine(self::Cohort)
        self.quarantine = true
        for i in self.houses
            i.quarantine = true
        end
    end

    function warn(self::Cohort)
        self.warning = true
        for i in self.houses
            i.warning = true
        end
        for i in self.classes
            i.warning = true
        end
    end

    export Cohort

end


text:
function test_pooling_layer(backend::Backend, pooling::PoolingFunction, has_padding::Bool, n_input, T, eps)
  println("-- Testing Pooling($(typeof(pooling))) $(has_padding? "with padding":"") on $(typeof(backend)){$T}...")
  println("    > Setup")

  if has_padding
    input_w   = 16
    input_h   = 10
    padding   = (2,2)
  else
    input_w   = 18
    input_h   = 12
    padding   = (0, 0)
  end
  input_chann = 3
  input_num   = 24
  kernel_w    = 3
  kernel_h    = 4
  stride_w    = 1
  stride_h    = 2

  dims = [rand(12:16, 4) for i =1:n_input]
  dims[1] = [input_w, input_h, input_chann, input_num]
  if isa(backend, AbstractGPUBackend)
    for i = 1:n_input
      # cuDNN pooling have different behavior on the boundary.
      # For mean pooling, cuDNN restricts the pooling size when it runs beyond the boundary.
      # Thus at the boundary, the mean values are larger than Mocha's default behavior,
      # because in Mocha's CPU implementation, the mean will count zeros beyond the the
      # boundary. I think this minor difference at the boundary is not super important.
      # So instead of messing around the implementation of GPU pooling (e.g. add explicit
      # padding when the pooling could run beyond the boundary), I use a workaround
      # in the unit test to make sure that the pooling always ends exactly at the boundary
      # in our unit tests.
      dims[i][1] -= (dims[i][1] + 2*padding[1] - kernel_w) % stride_w
      dims[i][2] -= (dims[i][2] + 2*padding[2] - kernel_h) % stride_h
    end
  end

  layer = PoolingLayer(kernel=(kernel_w,kernel_h), stride=(stride_w,stride_h), pad=padding,
      tops=Array{Symbol}(n_input), bottoms=Array{Symbol}(n_input), pooling=pooling)

  input = [rand(T, dims[i]...) for i = 1:n_input]
  inputs = Blob[make_blob(backend, x) for x in input]
  diffs = Blob[make_blob(backend, x) for x in input]

  state = setup(backend, layer, inputs, diffs)

  println("    > Forward")
  forward(backend, state, inputs)

  payloads = Array{Any}(n_input)
  for i = 1:n_input
    expected_output, payloads[i] = pooling_forward(state, i, input[i])
    got_output = similar(expected_output)
    copy!(got_output, state.blobs[i])
    @test all(-eps .< expected_output - got_output .< eps)
  end

  println("    > Backward")
  top_diff = [rand(T, size(state.blobs[i])) for i = 1:n_input]
  for i = 1:n_input
    copy!(state.blobs_diff[i], top_diff[i])
  end

  backward(backend, state, inputs, diffs)

  for i = 1:n_input
    expected_grad = pooling_backward(state, i, input[i], top_diff[i], payloads[i])
    got_grad = similar(expected_grad)
    copy!(got_grad, diffs[i])
    @test all(-eps .< expected_grad - got_grad .< eps)
  end

  shutdown(backend, state)
end

function pooling_forward(state, i, input::Array)
  width, height, channels, num = size(input)
  pooled_width = get_width(state.blobs[i])
  pooled_height = get_height(state.blobs[i])
  kernel_size = state.layer.kernel[1] * state.layer.kernel[2]

  output = zeros(eltype(input), pooled_width, pooled_height, channels, num)
  if isa(state.layer.pooling, Pooling.Max)
    mask = similar(output, Int)
  end

  for n = 1:num
    for c = 1:channels
      for ph = 1:pooled_height
        for pw = 1:pooled_width
          hstart = (ph-1)*state.layer.stride[2] - state.layer.pad[2] + 1
          wstart = (pw-1)*state.layer.stride[1] - state.layer.pad[1] + 1
          hend = min(hstart + state.layer.kernel[2] - 1, height)
          wend = min(wstart + state.layer.kernel[1] - 1, width)
          hstart = max(1, hstart)
          wstart = max(1, wstart)

          region = view(input, wstart:wend, hstart:hend, c, n)
          if isa(state.layer.pooling, Pooling.Max)
            index = indmax(region)
            mask[pw, ph, c, n] = index # note this is local index in region
            output[pw, ph, c, n] = region[index]
          elseif isa(state.layer.pooling, Pooling.Mean)
            output[pw, ph, c, n] = sum(region) / kernel_size
          else
            error("Unknown pooling $(state.layer.pooling)")
          end
        end
      end
    end
  end

  if isa(state.layer.pooling, Pooling.Max)
    return (output, mask)
  else
    return (output, nothing)
  end
end

function pooling_backward(state, i, input::Array, diff::Array, payload::Any)
  width, height, channels, num = size(input)
  pooled_width = get_width(state.blobs[i])
  pooled_height = get_height(state.blobs[i])
  kernel_size = state.layer.kernel[1] * state.layer.kernel[2]

  gradient = zeros(eltype(input), size(input))

  for n = 1:num
    for c = 1:channels
      for ph = 1:pooled_height
        for pw = 1:pooled_width
          hstart = (ph-1)*state.layer.stride[2] - state.layer.pad[2] + 1
          wstart = (pw-1)*state.layer.stride[1] - state.layer.pad[1] + 1
          hend = min(hstart + state.layer.kernel[2] - 1, height)
          wend = min(wstart + state.layer.kernel[1] - 1, width)
          hstart = max(1, hstart)
          wstart = max(1, wstart)

          region = view(gradient, wstart:wend, hstart:hend, c, n)
          if isa(state.layer.pooling, Pooling.Max)
            index = payload[pw, ph, c, n]
            region[index] += diff[pw, ph, c, n]
          elseif isa(state.layer.pooling, Pooling.Mean)
            region[:] += diff[pw, ph, c, n] / kernel_size
          else
            error("Unknown pooling $(state.layer.pooling)")
          end
        end
      end
    end
  end

  return gradient
end

function test_pooling_layer(backend::Backend, n_input, T, eps)
  test_pooling_layer(backend, Pooling.Max(), false, n_input, T, eps)
  test_pooling_layer(backend, Pooling.Mean(), false, n_input, T, eps)
  test_pooling_layer(backend, Pooling.Max(), true, n_input, T, eps)
  test_pooling_layer(backend, Pooling.Mean(), true, n_input, T, eps)
end

function test_pooling_layer(backend::Backend)
  test_pooling_layer(backend, 4, Float64, 1e-7)
  test_pooling_layer(backend, 2, Float32, 1e-3)
end

if test_cpu
  test_pooling_layer(backend_cpu)
end
if test_gpu
  test_pooling_layer(backend_gpu)
end


text:
using LightGraphs
export sosextractcycle, extractatomic, pushforward

#function dist(X, Y)
#    maximum(x -> abs(x[1] - x[2]) / max(x[1], x[2]), zip(X, Y))
#end

function findapprox(a, x, tol)
    _a(i) = dot(a[i][1], x)
    for i in eachindex(a)
        #if dist # isapprox(a[i][1], x, rtol=tol) # rtol is too restrictive since [0.999, 0.001] should be the same as [0.999, 0.002] but rtol says no
        if isapprox(_a(i), 1., rtol=tol)
            return i
        end
    end
    return 0
end

function pushapprox!(a, x, id, tol)
    i = findapprox(a, x, tol)
    if iszero(i)
        push!(a, (x, id))
        id
    else
        a[i][2]
    end
end

function project_atom(atom::WeightedDiracMeasure)
    nx = norm(atom.center, 2)
    return WeightedDiracMeasure(atom.center / nx, atom.weight * nx)
end
function pushforward(A::AbstractMatrix, atom::WeightedDiracMeasure)
    return project_atom(WeightedDiracMeasure(A * atom.center, atom.weight))
end
function pushforward(s::AbstractDiscreteSwitchedSystem, t,
                     atom::WeightedDiracMeasure)
    return pushforward(dynamicfort(s, t), atom)
end
function pushforward(s::AbstractDiscreteSwitchedSystem, t,
                     atom::WeightedDiracMeasure, γ)
    return pushforward(dynamicfort(s, t) / γ, atom)
end
function pushforward(s::AbstractDiscreteSwitchedSystem, t,
                     atomic::AtomicMeasure, args...)
    vars = variables(s, target(s, t))
    atoms = [pushforward(s, t, atom, args...) for atom in atomic.atoms]
    return AtomicMeasure(vars, atoms)
end

function extractatomic(s::AbstractDiscreteSwitchedSystem, d, t, ranktol,
                       dual=getlyap(s, d).dual)
    atoms = extractatoms(dual[t], ranktol)
    if atoms === nothing
        return nothing
    else
        return AtomicMeasure(atoms.variables, project_atom.(atoms.atoms))
    end
end

function extractstates(args...)
    atoms = extractatomic(args...)
    if atoms === nothing
        return MultivariateMoments.WeightedDiracMeasure{Float64}[]
    else
        atoms.atoms
    end
end

function startstates(s::AbstractDiscreteSwitchedSystem, edgestates, G, B, disttol)
    a = Tuple{Vector{Float64}, Int}[]
    for (t, states) in edgestates
        for atom in states
            x = atom.center
            i = pushapprox!(a, x, length(G)+1, disttol)
            σ = symbol(s, t)
            if i == length(G)+1
                push!(G, [(σ, t)])
                push!(B, x)
            else
                push!(G[i], (σ, t))
            end
        end
    end
    a
end

# Atom extraction

"""
    sosextractcycle(s::AbstractDiscreteSwitchedSystem, dual, d::Integer;
                    ranktols=1e-5, disttols=1e-5)

Extract cycles of high growth rate from atomic occupation measures given by the
infeasibility certificates of highest growth rate computed by
[`soslyap`](@ref). The method is detailed in [LJP17].

* [LJP17] B. Legat, R. M. Jungers, and P. A. Parrilo.
[Certifying unstability of Switched Systems using Sum of Squares Programming](https://arxiv.org/abs/1710.01814),
arXiv preprint arXiv:1710.01814, **2017**.
"""
function sosextractcycle(s::AbstractDiscreteSwitchedSystem, dual, d::Integer;
                         ranktols=1e-5, disttols=1e-5)
    smp = nothing
    for ranktol in ranktols
        # This part is the more costly since it does atom extraction
        # It is run only once for all disttols which is nice
        edgestates = map(u -> begin
                             map(t -> (t, extractstates(s, d, t, ranktol, dual)),
                                 out_transitions(s, u))
                         end,
                         states(s))

        for disttol in disttols
            G = Vector{Tuple{Int, transitiontype(s)}}[] # G[u] = list of edges (σ, v) going out of u
            B = Vector{Float64}[] # B[u] = values of the state at node u
            GB = map(u -> startstates(s, edgestates[u], G, B, disttol), states(s))

            # Try grouping with different distances

            g = DiGraph(length(G))
            w = fill(Inf, nv(g), nv(g))
            ij2t = Dict{NTuple{2, Int}, transitiontype(s)}()
            for u in eachindex(G)
                x = B[u]
                for (σ, t) in G[u]
                    atom = pushforward(s, t, WeightedDiracMeasure(x, one(eltype(x))))
                    y = atom.center
                    j = findapprox(GB[target(s, t)], y, disttol)
                    if !iszero(j)
                        v = GB[target(s, t)][j][2]
                        add_edge!(g, u, v)
                        # we use a minimum cycle arithmetic mean algo but we need the maximum cycle geometric mean
                        w[u, v] = -log(atom.weight)
                        ij2t[(u, v)] = t
                    end
                end
            end

            c, λmin = karp_minimum_cycle_mean(g, w)

            if !isempty(c)
                period = map(i -> ij2t[(c[i], c[(i % length(c)) + 1])], eachindex(c))

                newsmp = periodicswitching(s, period)
                notifyperiodic!(s, newsmp)
                if smp === nothing || isbetter(newsmp, smp)
                    smp = newsmp
                end
            end
        end
    end

    if smp !== nothing
        updatesmp!(s, smp)
    end

    return smp
end
function sosextractcycle(s::AbstractDiscreteSwitchedSystem, d::Integer;
                         tol=1e-5, ranktols=tol, disttols=tol,
                         kws...)::Union{Nothing, periodicswitchingtype(s)}
    lyap = getlyap(s, d; tol=tol, kws...)
    sosextractcycle(s, lyap.dual, d; ranktols=ranktols, disttols=disttols)
end


text:
@testset "Vector" begin
    x = randn(2)
    syms = (:a, :b)

    # Constructors
    test_rrule(LArray{syms}, x)
    constructor = @SLArray (2, ) syms
    test_rrule(constructor, x)

    # `getproperty`
    lx = LArray{syms}(x)
    test_rrule(getproperty, lx, first(syms))
    test_rrule(getproperty, lx, last(syms))

    slx = constructor(x)
    test_rrule(getproperty, slx, first(syms))
    test_rrule(getproperty, slx, last(syms))
end


text:
# Autogenerated wrapper script for ranger_jll for armv7l-linux-musleabihf-cxx03
export libranger

JLLWrappers.@generate_wrapper_header("ranger")
JLLWrappers.@declare_library_product(libranger, "libranger.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libranger,
        "lib/libranger.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Compare

function main()
    ids = []
    open("../input/input.txt") do file
        for line in eachline(file)
            id = map(x -> x[1], split(line, ""))
            for other in ids
                diff = other .!= id
                if sum(diff) == 1
                    # Print what's in common.
                    println(join(id[.!diff]))
                    # No labeled break in Julia yet.
                    # See also: https://github.com/JuliaLang/julia/issues/5334
                    @goto done
                end
            end
            push!(ids, id)
        end
        @label done
    end
end

main()

end

max_stars_repo_name: tjpalmer/advent

max_issues_repo_name: tjpalmer/advent

max_forks_repo_name: tjpalmer/advent


text:
export sequence_of_three, alternate_items

function sequence_of_three(start)
    out = []
    push!(out, start, start+1, start+2)
    return out
end

function alternate_items(nrepeats::Integer, item1, item2)
    out = []
    for _ = 1:nrepeats           # _ is often used for a variable you'll not use again
        push!(out, item1, item2)
    end
    return out
end