text:
# ------------------------------------
# muscl
# ------------------------------------
function muscl(Qbase, QbaseU, QbaseD, QbaseL, QbaseR, QbaseF, QbaseB, QconU, QconD, QconL, QconR, QconF, QconB,
                cellxmax, cellymax, cellzmax, nval, icell, specific_heat_ratio, volume, nodes)

    #=
    Q system
    VecA system is sama
    Q is defined on(near) boundary

            nodes[2,3]  QU[2,3]
                |--------------------|
                |       QD[2,3]      |
                |                    |
        QL[2,2] | QR[2,2]    QL[3,2] | QR[3,2]
                |      cell(2,2)     |
                |                    |
                |       QU[2,2]      |
                |--------------------|
            nodes[2,2]  QD[2,2]   nodes[3,2]

    =#
    for l in 1:nval
        for k in 1+icell-1:cellzmax-icell+1
            for j in 1+icell-1:cellymax-icell+1
                for i in 1+icell:cellxmax-icell+1
                    # ------------------------------
                    # L
                    # ------------------------------

                    lim = 0.25 * (((nodes[i-2,j+1,k+1,1] - nodes[i-1,j+1,k+1,1])^2
                                 + (nodes[i-2,j+1,k+1,2] - nodes[i-1,j+1,k+1,2])^2
                                 + (nodes[i-2,j+1,k+1,3] - nodes[i-1,j+1,k+1,3])^2)^0.5
                                 +((nodes[i-2,  j,k+1,1] - nodes[i-1,  j,k+1,1])^2
                                 + (nodes[i-2,  j,k+1,2] - nodes[i-1,  j,k+1,2])^2
                                 + (nodes[i-2,  j,k+1,3] - nodes[i-1,  j,k+1,3])^2)^0.5
                                 +((nodes[i-2,j+1,  k,1] - nodes[i-1,j+1,  k,1])^2
                                 + (nodes[i-2,j+1,  k,2] - nodes[i-1,j+1,  k,2])^2
                                 + (nodes[i-2,j+1,  k,3] - nodes[i-1,j+1,  k,3])^2)^0.5
                                 +((nodes[i-2,  j,  k,1] - nodes[i-1,  j,  k,1])^2
                                 + (nodes[i-2,  j,  k,2] - nodes[i-1,  j,  k,2])^2
                                 + (nodes[i-2,  j,  k,3] - nodes[i-1,  j,  k,3])^2)^0.5)

                    li  = 0.25 * (((nodes[i-1,j+1,k+1,1] - nodes[  i,j+1,k+1,1])^2
                                 + (nodes[i-1,j+1,k+1,2] - nodes[  i,j+1,k+1,2])^2
                                 + (nodes[i-1,j+1,k+1,3] - nodes[  i,j+1,k+1,3])^2)^0.5
                                 +((nodes[i-1,  j,k+1,1] - nodes[  i,  j,k+1,1])^2
                                 + (nodes[i-1,  j,k+1,2] - nodes[  i,  j,k+1,2])^2
                                 + (nodes[i-1,  j,k+1,3] - nodes[  i,  j,k+1,3])^2)^0.5
                                 +((nodes[i-1,j+1,  k,1] - nodes[  i,j+1,  k,1])^2
                                 + (nodes[i-1,j+1,  k,2] - nodes[  i,j+1,  k,2])^2
                                 + (nodes[i-1,j+1,  k,3] - nodes[  i,j+1,  k,3])^2)^0.5
                                 +((nodes[i-1,  j,  k,1] - nodes[  i,  j,  k,1])^2
                                 + (nodes[i-1,  j,  k,2] - nodes[  i,  j,  k,2])^2
                                 + (nodes[i-1,  j,  k,3] - nodes[  i,  j,  k,3])^2)^0.5)

                    lip = 0.25 * (((nodes[  i,j+1,k+1,1] - nodes[i+1,j+1,k+1,1])^2
                                 + (nodes[  i,j+1,k+1,2] - nodes[i+1,j+1,k+1,2])^2
                                 + (nodes[  i,j+1,k+1,3] - nodes[i+1,j+1,k+1,3])^2)^0.5
                                 +((nodes[  i,  j,k+1,1] - nodes[i+1,  j,k+1,1])^2
                                 + (nodes[  i,  j,k+1,2] - nodes[i+1,  j,k+1,2])^2
                                 + (nodes[  i,  j,k+1,3] - nodes[i+1,  j,k+1,3])^2)^0.5
                                 +((nodes[  i,j+1,  k,1] - nodes[i+1,j+1,  k,1])^2
                                 + (nodes[  i,j+1,  k,2] - nodes[i+1,j+1,  k,2])^2
                                 + (nodes[  i,j+1,  k,3] - nodes[i+1,j+1,  k,3])^2)^0.5
                                 +((nodes[  i,  j,  k,1] - nodes[i+1,  j,  k,1])^2
                                 + (nodes[  i,  j,  k,2] - nodes[i+1,  j,  k,2])^2
                                 + (nodes[  i,  j,  k,3] - nodes[i+1,  j,  k,3])^2)^0.5)

                    dip = 2*(Qbase[i,j,k,l]   - Qbase[i-1,j,k,l]) / (li  + lip)
                    dim = 2*(Qbase[i-1,j,k,l] - Qbase[i-2,j,k,l]) / (lim + li)
                    QbaseL[i,j,k,l] = Qbase[i-1,j,k,l] + 0.5*li*minmod(dim, dip)

                    # ------------------------------
                    # R
                    # ------------------------------

                    lim = li
                    li  = lip
                    lip = 0.25 * (((nodes[i+1,j+1,k+1,1] - nodes[i+2,j+1,k+1,1])^2
                                 + (nodes[i+1,j+1,k+1,2] - nodes[i+2,j+1,k+1,2])^2
                                 + (nodes[i+1,j+1,k+1,3] - nodes[i+2,j+1,k+1,3])^2)^0.5
                                 +((nodes[i+1,  j,k+1,1] - nodes[i+2,  j,k+1,1])^2
                                 + (nodes[i+1,  j,k+1,2] - nodes[i+2,  j,k+1,2])^2
                                 + (nodes[i+1,  j,k+1,3] - nodes[i+2,  j,k+1,3])^2)^0.5
                                 +((nodes[i+1,j+1,  k,1] - nodes[i+2,j+1,  k,1])^2
                                 + (nodes[i+1,j+1,  k,2] - nodes[i+2,j+1,  k,2])^2
                                 + (nodes[i+1,j+1,  k,3] - nodes[i+2,j+1,  k,3])^2)^0.5
                                 +((nodes[i+1,  j,  k,1] - nodes[i+2,  j,  k,1])^2
                                 + (nodes[i+1,  j,  k,2] - nodes[i+2,  j,  k,2])^2
                                 + (nodes[i+1,  j,  k,3] - nodes[i+2,  j,  k,3])^2)^0.5)

                    dip = 2*(Qbase[i+1,j,k,l] - Qbase[i,j,k,l]) / (li  + lip)
                    dim = 2*(Qbase[i,j,k,l]   - Qbase[i-1,j,k,l]) / (lim + li)
                    QbaseR[i,j,k,l] = Qbase[i,j,k,l]  - 0.5*li*minmod(dip, dim)
                end
            end
        end
    end

    for l in 1:nval
        for k in 1+icell-1:cellzmax-icell+1
            for j in 1+icell:cellymax-icell + 1
                for i in 1+icell-1:cellxmax-icell+1
                    # ------------------------------
                    # D
                    # ------------------------------

                    lim = 0.25 * (((nodes[i+1,j-2,k+1,1] - nodes[i+1,j-1,k+1,1])^2
                                 + (nodes[i+1,j-2,k+1,2] - nodes[i+1,j-1,k+1,2])^2
                                 + (nodes[i+1,j-2,k+1,3] - nodes[i+1,j-1,k+1,3])^2)^0.5
                                 +((nodes[  i,j-2,k+1,1] - nodes[  i,j-1,k+1,1])^2
                                 + (nodes[  i,j-2,k+1,2] - nodes[  i,j-1,k+1,2])^2
                                 + (nodes[  i,j-2,k+1,3] - nodes[  i,j-1,k+1,3])^2)^0.5
                                 +((nodes[i+1,j-2,  k,1] - nodes[i+1,j-1,  k,1])^2
                                 + (nodes[i+1,j-2,  k,2] - nodes[i+1,j-1,  k,2])^2
                                 + (nodes[i+1,j-2,  k,3] - nodes[i+1,j-1,  k,3])^2)^0.5
                                 +((nodes[  i,j-2,  k,1] - nodes[  i,j-1,  k,1])^2
                                 + (nodes[  i,j-2,  k,2] - nodes[  i,j-1,  k,2])^2
                                 + (nodes[  i,j-2,  k,3] - nodes[  i,j-1,  k,3])^2)^0.5)

                    li  = 0.25 * (((nodes[i+1,j-1,k+1,1] - nodes[i+1,  j,k+1,1])^2
                                 + (nodes[i+1,j-1,k+1,2] - nodes[i+1,  j,k+1,2])^2
                                 + (nodes[i+1,j-1,k+1,3] - nodes[i+1,  j,k+1,3])^2)^0.5
                                 +((nodes[  i,j-1,k+1,1] - nodes[  i,  j,k+1,1])^2
                                 + (nodes[  i,j-1,k+1,2] - nodes[  i,  j,k+1,2])^2
                                 + (nodes[  i,j-1,k+1,3] - nodes[  i,  j,k+1,3])^2)^0.5
                                 +((nodes[i+1,j-1,  k,1] - nodes[i+1,  j,  k,1])^2
                                 + (nodes[i+1,j-1,  k,2] - nodes[i+1,  j,  k,2])^2
                                 + (nodes[i+1,j-1,  k,3] - nodes[i+1,  j,  k,3])^2)^0.5
                                 +((nodes[  i,j-1,  k,1] - nodes[  i,  j,  k,1])^2
                                 + (nodes[  i,j-1,  k,2] - nodes[  i,  j,  k,2])^2
                                 + (nodes[  i,j-1,  k,3] - nodes[  i,  j,  k,3])^2)^0.5)

                    lip = 0.25 * (((nodes[i+1,  j,k+1,1] - nodes[i+1,j+1,k+1,1])^2
                                 + (nodes[i+1,  j,k+1,2] - nodes[i+1,j+1,k+1,2])^2
                                 + (nodes[i+1,  j,k+1,3] - nodes[i+1,j+1,k+1,3])^2)^0.5
                                 +((nodes[  i,  j,k+1,1] - nodes[  i,j+1,k+1,1])^2
                                 + (nodes[  i,  j,k+1,2] - nodes[  i,j+1,k+1,2])^2
                                 + (nodes[  i,  j,k+1,3] - nodes[  i,j+1,k+1,3])^2)^0.5
                                 +((nodes[i+1,  j,  k,1] - nodes[i+1,j+1,  k,1])^2
                                 + (nodes[i+1,  j,  k,2] - nodes[i+1,j+1,  k,2])^2
                                 + (nodes[i+1,  j,  k,3] - nodes[i+1,j+1,  k,3])^2)^0.5
                                 +((nodes[  i,  j,  k,1] - nodes[  i,j+1,  k,1])^2
                                 + (nodes[  i,  j,  k,2] - nodes[  i,j+1,  k,2])^2
                                 + (nodes[  i,  j,  k,3] - nodes[  i,j+1,  k,3])^2)^0.5)

                    dip = 2*(Qbase[i,j,k,l]   - Qbase[i,j-1,k,l]) / (li  + lip)
                    dim = 2*(Qbase[i,j-1,k,l] - Qbase[i,j-2,k,l]) / (lim + li)
                    QbaseD[i,j,k,l] = Qbase[i,j-1,k,l] + 0.5*li*minmod(dim, dip)

                    # ------------------------------
                    # U
                    # ------------------------------

                    lim = li
                    li  = lip
                    lip = 0.25 * (((nodes[i+1,j+1,k+1,1] - nodes[i+1,j+2,k+1,1])^2
                                 + (nodes[i+1,j+1,k+1,2] - nodes[i+1,j+2,k+1,2])^2
                                 + (nodes[i+1,j+1,k+1,3] - nodes[i+1,j+2,k+1,3])^2)^0.5
                                 +((nodes[  i,j+1,k+1,1] - nodes[  i,j+2,k+1,1])^2
                                 + (nodes[  i,j+1,k+1,2] - nodes[  i,j+2,k+1,2])^2
                                 + (nodes[  i,j+1,k+1,3] - nodes[  i,j+2,k+1,3])^2)^0.5
                                 +((nodes[i+1,j+1,  k,1] - nodes[i+1,j+2,  k,1])^2
                                 + (nodes[i+1,j+1,  k,2] - nodes[i+1,j+2,  k,2])^2
                                 + (nodes[i+1,j+1,  k,3] - nodes[i+1,j+2,  k,3])^2)^0.5
                                 +((nodes[  i,j+1,  k,1] - nodes[  i,j+2,  k,1])^2
                                 + (nodes[  i,j+1,  k,2] - nodes[  i,j+2,  k,2])^2
                                 + (nodes[  i,j+1,  k,3] - nodes[  i,j+2,  k,3])^2)^0.5)

                    dip = 2*(Qbase[i,j+1,k,l] - Qbase[i,j,k,l])   / (li  + lip)
                    dim = 2*(Qbase[i,j,k,l]   - Qbase[i,j-1,k,l]) / (lim + li)
                    QbaseU[i,j,k,l] = Qbase[i,j,k,l]  - 0.5*li*minmod(dip, dim)

                    #=if isequal(QbaseD[i,j,k,l], NaN) == true
                        println("  ")
                        println(i)
                        println(j)
                        println(k)
                        println(dip)
                        println(dim)
                    end
                    =#
                end
            end
        end
    end

    for l in 1:nval
        for k in 1+icell:cellzmax-icell+1
            for j in 1+icell-1:cellymax-icell + 1
                for i in 1+icell-1:cellxmax-icell+1
                    # ------------------------------
                    # B
                    # ------------------------------

                    lim = 0.25 * (((nodes[i+1,j+1,k-2,1] - nodes[i+1,j+1,k-1,1])^2
                                 + (nodes[i+1,j+1,k-2,2] - nodes[i+1,j+1,k-1,2])^2
                                 + (nodes[i+1,j+1,k-2,3] - nodes[i+1,j+1,k-1,3])^2)^0.5
                                 +((nodes[  i,j+1,k-2,1] - nodes[  i,j+1,k-1,1])^2
                                 + (nodes[  i,j+1,k-2,2] - nodes[  i,j+1,k-1,2])^2
                                 + (nodes[  i,j+1,k-2,3] - nodes[  i,j+1,k-1,3])^2)^0.5
                                 +((nodes[i+1,  j,k-2,1] - nodes[i+1,  j,k-1,1])^2
                                 + (nodes[i+1,  j,k-2,2] - nodes[i+1,  j,k-1,2])^2
                                 + (nodes[i+1,  j,k-2,3] - nodes[i+1,  j,k-1,3])^2)^0.5
                                 +((nodes[  i,  j,k-2,1] - nodes[  i,  j,k-1,1])^2
                                 + (nodes[  i,  j,k-2,2] - nodes[  i,  j,k-1,2])^2
                                 + (nodes[  i,  j,k-2,3] - nodes[  i,  j,k-1,3])^2)^0.5)

                    li  = 0.25 * (((nodes[i+1,j+1,k-1,1] - nodes[i+1,j+1,  k,1])^2
                                 + (nodes[i+1,j+1,k-1,2] - nodes[i+1,j+1,  k,2])^2
                                 + (nodes[i+1,j+1,k-1,3] - nodes[i+1,j+1,  k,3])^2)^0.5
                                 +((nodes[  i,j+1,k-1,1] - nodes[  i,j+1,  k,1])^2
                                 + (nodes[  i,j+1,k-1,2] - nodes[  i,j+1,  k,2])^2
                                 + (nodes[  i,j+1,k-1,3] - nodes[  i,j+1,  k,3])^2)^0.5
                                 +((nodes[i+1,  j,k-1,1] - nodes[i+1,  j,  k,1])^2
                                 + (nodes[i+1,  j,k-1,2] - nodes[i+1,  j,  k,2])^2
                                 + (nodes[i+1,  j,k-1,3] - nodes[i+1,  j,  k,3])^2)^0.5
                                 +((nodes[  i,  j,k-1,1] - nodes[  i,  j,  k,1])^2
                                 + (nodes[  i,  j,k-1,2] - nodes[  i,  j,  k,2])^2
                                 + (nodes[  i,  j,k-1,3] - nodes[  i,  j,  k,3])^2)^0.5)

                    lip = 0.25 * (((nodes[i+1,j+1,  k,1] - nodes[i+1,j+1,k+1,1])^2
                                 + (nodes[i+1,j+1,  k,2] - nodes[i+1,j+1,k+1,2])^2
                                 + (nodes[i+1,j+1,  k,3] - nodes[i+1,j+1,k+1,3])^2)^0.5
                                 +((nodes[  i,j+1,  k,1] - nodes[  i,j+1,k+1,1])^2
                                 + (nodes[  i,j+1,  k,2] - nodes[  i,j+1,k+1,2])^2
                                 + (nodes[  i,j+1,  k,3] - nodes[  i,j+1,k+1,3])^2)^0.5
                                 +((nodes[i+1,  j,  k,1] - nodes[i+1,  j,k+1,1])^2
                                 + (nodes[i+1,  j,  k,2] - nodes[i+1,  j,k+1,2])^2
                                 + (nodes[i+1,  j,  k,3] - nodes[i+1,  j,k+1,3])^2)^0.5
                                 +((nodes[  i,  j,  k,1] - nodes[  i,  j,k+1,1])^2
                                 + (nodes[  i,  j,  k,2] - nodes[  i,  j,k+1,2])^2
                                 + (nodes[  i,  j,  k,3] - nodes[  i,  j,k+1,3])^2)^0.5)

                    dip = 2*(Qbase[i,j,k,l]   - Qbase[i,j,k-1,l]) / (li  + lip)
                    dim = 2*(Qbase[i,j,k-1,l] - Qbase[i,j,k-2,l]) / (lim + li)
                    QbaseB[i,j,k,l] = Qbase[i,j,k-1,l] + 0.5*li*minmod(dim, dip)

                    # ------------------------------
                    # F
                    # ------------------------------

                    lim = li
                    li  = lip
                    lip = 0.25 * (((nodes[i+1,j+1,k+1,1] - nodes[i+1,j+1,k+2,1])^2
                                 + (nodes[i+1,j+1,k+1,2] - nodes[i+1,j+1,k+2,2])^2
                                 + (nodes[i+1,j+1,k+1,3] - nodes[i+1,j+1,k+2,3])^2)^0.5
                                 +((nodes[  i,j+1,k+1,1] - nodes[  i,j+1,k+2,1])^2
                                 + (nodes[  i,j+1,k+1,2] - nodes[  i,j+1,k+2,2])^2
                                 + (nodes[  i,j+1,k+1,3] - nodes[  i,j+1,k+2,3])^2)^0.5
                                 +((nodes[i+1,  j,k+1,1] - nodes[i+1,  j,k+2,1])^2
                                 + (nodes[i+1,  j,k+1,2] - nodes[i+1,  j,k+2,2])^2
                                 + (nodes[i+1,  j,k+1,3] - nodes[i+1,  j,k+2,3])^2)^0.5
                                 +((nodes[  i,  j,k+1,1] - nodes[  i,  j,k+2,1])^2
                                 + (nodes[  i,  j,k+1,2] - nodes[  i,  j,k+2,2])^2
                                 + (nodes[  i,  j,k+1,3] - nodes[  i,  j,k+2,3])^2)^0.5)

                    dip = 2*(Qbase[i,j,k+1,l] - Qbase[i,j,k,l])   / (li  + lip)
                    dim = 2*(Qbase[i,j,k,l]   - Qbase[i,j,k-1,l]) / (lim + li)
                    QbaseF[i,j,k,l] = Qbase[i,j,k,l]  - 0.5*li*minmod(dip, dim)
                end
            end
        end
    end

    QconU = base_to_conservative(QbaseU, QconU, cellxmax, cellymax, cellzmax, specific_heat_ratio)
    QconD = base_to_conservative(QbaseD, QconD, cellxmax, cellymax, cellzmax, specific_heat_ratio)
    QconL = base_to_conservative(QbaseL, QconL, cellxmax, cellymax, cellzmax, specific_heat_ratio)
    QconR = base_to_conservative(QbaseR, QconR, cellxmax, cellymax, cellzmax, specific_heat_ratio)
    QconF = base_to_conservative(QbaseF, QconF, cellxmax, cellymax, cellzmax, specific_heat_ratio)
    QconB = base_to_conservative(QbaseB, QconB, cellxmax, cellymax, cellzmax, specific_heat_ratio)

    #=
    println("")
    println(Qbase[49,50,2])
    println(QbaseL[50,50,2])
    println(Qbase[50,50,2])
    =#

    return QbaseU, QbaseD, QbaseL, QbaseR, QbaseF, QbaseB, QconU, QconD, QconL, QconR, QconF, QconB
end

function minmod(a,b)
    ans = sign(a) * max(0, min(a*sign(b), b*sign(a)))
    return ans
end


text:
# Test functionality of domain constants
path = joinpath(dirname(pathof(PDDL)), "..", "test", "constants")

domain = load_domain(joinpath(path, "domain.pddl"))
@test domain.name == Symbol("taxi")
problem = load_problem(joinpath(path, "problem.pddl"))
@test problem.name == Symbol("taxi-problem")

# Check that constants are loaded correctly
@test domain.constants == @pddl("red", "green", "yellow", "blue", "intaxi")

# Check that types of constants only resolve when domain is provided
state = init_state(problem)
@test state[pddl"(pasloc red)"] == false
@test state[domain, pddl"(pasloc red)"] == true

# Execute plan and check that it succeeds
plan = @pddl(
    "(move loc6 loc5 west)",
    "(move loc5 loc0 north)",
    "(pickup loc0 red)",
    "(move loc0 loc1 east)",
    "(move loc1 loc6 south)",
    "(move loc6 loc11 south)",
    "(move loc11 loc12 east)",
    "(move loc12 loc13 east)",
    "(move loc13 loc14 east)",
    "(move loc14 loc9 north)",
    "(move loc9 loc4 north)",
    "(dropoff loc4 green)"
)
state = execute(plan, state, domain)
@test satisfy(problem.goal, state)[1] == true

ext: jl


text:
"An abstract type for Runge-Kutta solvers."
abstract type AbstractRungeKuttaSolver <: AbstractInitialValueSolver end

"An abstract type for the solution of an [`AbstractInitialValueProblem`](@ref) obtained with an [`AbstractRungeKuttaSolver`](@ref)."
abstract type AbstractRungeKuttaSolution <: AbstractInitialValueSolution end

"An abstract type for the cache of an [`AbstractRungeKuttaSolver`](@ref)."
abstract type AbstractRungeKuttaCache end

"An abstract type for the parameters of an [`AbstractRungeKuttaSolver`](@ref)."
abstract type AbstractRungeKuttaParameters end

"An abstract type for the Butcher tableau of an [`AbstractRungeKuttaSolver`](@ref)."
abstract type AbstractButcherTableau <: AbstractRungeKuttaParameters end

"An abstract type for the parameters of an adaptive [`AbstractRungeKuttaSolver`](@ref)."
abstract type AbstractAdaptiveParameters <: AbstractRungeKuttaParameters end


text:
"""
A (abstract) `Polytope` describes the combinatorial structure used to discretize a surface from a set of given points (vertices)
linked together with edges, faces, 3D cells and so on. It does not include any location information, and therefore cannot
be represented as geometry. It is usually represented with a Hasse diagram, which provides a way to read off the structure
of a polytope as a partially ordered set.
Only when a polytope is associated with vertex data can we visually represent a possible _realization_ of this particular
polytope. Such a realization is called a geometric polytope.
"""
struct Polytope{N}
    graph::SimpleDiGraph{Int}
    nverts::Int
    start_indices::Vector{Int}
end

const Polygon = Polytope{2}
const Polyhedron = Polytope{3}

Polytope(g::SimpleDiGraph, k::Integer, nverts::Integer, start_indices::AbstractVector{<:Integer}) = Polytope{k}(g, nverts, start_indices)

"""
Build a `Polytope` from the structure provided in `indices`.
The nth argument passed as `indices` describes connectivity from the (n-1)th face to the nth face.

## Example

```julia
julia> Polytope(4,
    [(1, 2), (2, 3), (3, 1), (2, 3), (3, 4), (4, 2)], # edges (segments)
    [(1, 2, 3), (4, 5, 6)] # 2 triangle faces from edges 1, 2, 3 and 4, 5, 6 respectively
)

julia> Polytope(4,
    [(1, 2, 3, 1), (2, 3, 4, 2)], # edges (closed chains)
    [(1,), (2,)]) # faces (made from a closed chain -> PolySurface with no inner chains)
```
"""
function Polytope(connectivity::Vector{<:AbstractVector})
    faces_starts = start_indices(connectivity)
    nnodes = last(faces_starts) + length(last(connectivity))
    g = SimpleDiGraph(nnodes)
    add_faces!.(Ref(g), faces_starts, vcat(0, faces_starts[1:end-1]), connectivity)
    Polytope(g, length(connectivity), first(faces_starts), faces_starts)
end

Polytope(connectivity::AbstractVector...) = Polytope(collect(connectivity))

Base.:(==)(x::Polytope{K1}, y::Polytope{K2}) where {K1,K2} = false
Base.:(==)(x::Polytope{K}, y::Polytope{K}) where {K} = x.graph == y.graph && nverts(x) == nverts(y) && start_indices(x) == start_indices(y)

function add_faces!(g::SimpleDiGraph, src_start::Integer, dst_start::Integer, faces)
    for (src, dst_edge) ∈ enumerate(faces)
        for dst ∈ dst_edge
            add_edge!(g, src + src_start, dst + dst_start)
        end
    end
    g
end

start_indices(connectivity) = cumsum(maximum(maximum.(faces)) for faces ∈ connectivity)
start_indices(p::Polytope) = p.start_indices

nverts(p::Polytope) = p.nverts

rank(::Type{<:Polytope{N}}) where {N} = N
rank(p::Polytope) = rank(typeof(p))

function polytope_range(p::Polytope, k::Integer)
    @assert k ≠ rank(p)
    1:(start_indices(p)[k+1])
end

faces(p::Polytope, k) = k == rank(p) ? p : Polytope(p.graph[polytope_range(p, k)], k, nverts(p), start_indices(p)[1:k])

facets(p::Polytope) = faces(p, rank(p) - 1)

"""
Vertex associated with a specific position in space.
"""
struct Vertex{D,T}
    location::SVector{D,T}
end

Vertex(location::AbstractVector) = Vertex(SVector{length(location)}(location))
Vertex(coords::Number...) = Vertex(collect(coords))

Base.length(::Type{<:Vertex{D}}) where {D} = D

struct GeometricPolytope{D,P<:Polytope,V<:AbstractVector{<:Vertex{D}}}
    polytope::P
    vertices::V
    function GeometricPolytope(p::Polytope, vertices::AbstractVector{<:Vertex})
        @assert nverts(p) == length(vertices) "$(nverts(p)) vertices required for $p, found $(length(vertices))"
        new{length(eltype(vertices)),typeof(p),typeof(vertices)}(p, vertices)
    end
end

embeddim(::Type{<:GeometricPolytope{D}}) where {D} = D
embeddim(gp::GeometricPolytope) = embeddim(typeof(gp))


text:
using JuSOT
using Base.Test

# write your own tests here
@test 1 == 1

@test test_function1() == 1
@test test_ut1() == 1


text:
################### Filter and Threshold #####################
# pf = Filter(p)
a_f(r,u,v) = r^2*(∇(v)⊙∇(u))+v⊙u
function Filter(pvec,r,flag_f::Bool,P,Pf,Qf,dΩ,dΓ)
    if (flag_f)
        ph = FEFunction(P,pvec)
        op = AffineFEOperator(Pf,Qf) do u, v
            ∫( a_f(r,u,v))dΩ, ∫( v*ph )dΩ,∫( 0*v )dΓ
          end
        pf = solve(op)
        return get_free_values(pf)
    else
        return pvec
    end
end

# Threshold function
Threshold(β,η,flag_t::Bool,pf) = flag_t==false ? pf : ((tanh(β*η)+tanh(β*(pf-η)))/(tanh(β*η)+tanh(β*(1.0-η))))

max_stars_repo_name: Tiamat-Tech/hpinn

max_issues_repo_name: louhz/hpinn

max_forks_repo_name: louhz/hpinn


text:
module IndirectCoStateInit

using Reexport
@reexport using Heuristics
using StaticArrays

# Includes
include("CoStateInitializer.jl")
include("HeuristicsCoStateInitializer.jl")
include("FSSCoStateInitializer.jl")

# Exports
export FSSCoStateInitializer
export initialize!
export GetInitializedCostates

end


text:
push!(LOAD_PATH, joinpath(@__DIR__, ".."))

using BenchmarkTools
using CUDA
using Oceananigans
using Benchmarks

# Utility functions for generating tracer lists

function active_tracers(n)
    n == 0 && return []
    n == 1 && return [:b]
    n == 2 && return [:T, :S]
    throw(ArgumentError("Can't have more than 2 active tracers!"))
end

passive_tracers(n) = [Symbol("C" * string(m)) for m in 1:n]

tracer_list(n_active, n_passive) =
    Tuple(vcat(active_tracers(n_active), passive_tracers(n_passive)))

function buoyancy(n_active)
    n_active == 0 && return nothing
    n_active == 1 && return BuoyancyTracer()
    n_active == 2 && return SeawaterBuoyancy()
    throw(ArgumentError("Can't have more than 2 active tracers!"))
end

# Benchmark function

function benchmark_tracers(Arch, N, n_tracers)
    n_active, n_passive = n_tracers
    grid = RegularRectilinearGrid(size=(N, N, N), extent=(1, 1, 1))
    model = IncompressibleModel(architecture=Arch(), grid=grid, buoyancy=buoyancy(n_active),
                                tracers=tracer_list(n_active, n_passive))

    time_step!(model, 1) # warmup

    trial = @benchmark begin
        @sync_gpu time_step!($model, 1)
    end samples=10

    return trial
end

# Benchmark parameters

Architectures = has_cuda() ? [CPU, GPU] : [CPU]
Ns = [128]

# Each test case specifies (number of active tracers, number of passive tracers)
tracers = [(0, 0), (0, 1), (0, 2), (1, 0), (2, 0), (2, 3), (2, 5), (2, 10)]

# Run benchmarks

print_system_info()
suite = run_benchmarks(benchmark_tracers; Architectures, Ns, tracers)

df = benchmarks_dataframe(suite)
sort!(df, [:Architectures, :tracers, :Ns], by=(string, string, identity))
benchmarks_pretty_table(df, title="Arbitrary tracers benchmarks")

if GPU in Architectures
    df_Δ = gpu_speedups_suite(suite) |> speedups_dataframe
    sort!(df_Δ, [:tracers, :Ns])
    benchmarks_pretty_table(df_Δ, title="Arbitrary tracers CPU to GPU speedup")
end

for Arch in Architectures
    suite_arch = speedups_suite(suite[@tagged Arch], base_case=(Arch, Ns[1], (0, 0)))
    df_arch = speedups_dataframe(suite_arch, slowdown=true)
    sort!(df_arch, [:tracers, :Ns], by=(string, identity))
    benchmarks_pretty_table(df_arch, title="Arbitrary tracers relative performance ($Arch)")
end


text:
function POMDPs.requirements_info(solver::AbstractMCTSSolver, problem::Union{POMDP,MDP})
    if @implemented initialstate(::typeof(problem), ::MersenneTwister)
        requirements_info(solver, problem, initialstate(problem, MersenneTwister(1)))
    elseif statetype(typeof(problem)) <: Number
        s = one(statetype(typeof(problem)))
        requirements_info(solver, problem, s)
    else
        println("""
            Since MCTS is an online solver, most of the computation occurs in `action(policy, state)`. In order to view the requirements for this function, please, supply a state as the third argument to `requirements_info`, e.g.

                @requirements_info $(typeof(solver))() $(typeof(problem))() $(statetype(typeof(problem)))()

                """)
    end
end

function POMDPs.requirements_info(solver::AbstractMCTSSolver, problem::Union{POMDP,MDP}, s)
    policy = solve(solver, problem)
    requirements_info(policy, s)
end

function POMDPs.requirements_info(policy::AbstractMCTSPlanner, s)
    if !isequal(deepcopy(s), s)
        @warn("""
             isequal(deepcopy(s), s) returned false. Is isequal() defined correctly?

             For MCTS to work correctly, you must define isequal(::$(typeof(s)), ::$(typeof(s))) (see https://docs.julialang.org/en/stable/stdlib/collections/#Associative-Collections-1, https://github.com/andrewcooke/AutoHashEquals.jl#background, also consider using StaticArrays). This warning was thrown because isequal($(deepcopy(s)), $s) returned false.

             Note: isequal() should also be defined correctly for actions, but no warning will be issued.
             """)
    end
    if hash(deepcopy(s)) != hash(s)
        @warn("""
             hash(deepcopy(s)) was not equal to hash(s). Is hash() defined correctly?

             For MCTS to work correctly, you must define hash(::$(typeof(s)), ::UInt) (see https://docs.julialang.org/en/stable/stdlib/collections/#Associative-Collections-1, https://github.com/andrewcooke/AutoHashEquals.jl#background, also consider using StaticArrays). This warning was thrown because hash($(deepcopy(s))) != hash($s).

             Note: hash() should also be defined correctly for actions, but no warning will be issued.
             """)
    end
    @show_requirements action(policy, s)
end


text:
using FilePaths
using Base.Test

@testset "FilePaths" begin

include("mode.jl")
include("path.jl")

end


text:
"""
Throw an `AssertionError` if conditions like `op(exp1, exp2)` are `false`, where `op` is a conditional infix operator.

# Examples

```
julia> a = 3; b = 4;
julia> @assert_op a == b
ERROR: AssertionError: 3 == 4

julia> @assert_op a + 3 > b + 4
ERROR: AssertionError: 6 > 8
```
"""
macro assert_op(expr)
    assert_op(expr)
end

function assert_op(expr::Expr)
    # Only special case expressions of the form `expr1 == expr2`
    if length(expr.args) == 3 && expr.head == :call
        return assert_op(expr.args[1], expr.args[2], expr.args[3])
    else
        return :(@assert $(expr))
    end
end

function assert_op(op, exp1, exp2)
    return :(
        if !$op($(esc(exp1)), $(esc(exp2)))
            val1 = $(esc(exp1))
            val2 = $(esc(exp2))
            op_str = $(esc(op))
            throw(AssertionError("$val1 $op_str $val2"))
        end
    )
end


text:
# / Add this offset to convert Modified Julian Date to Julian Date.
const MJD_TO_JD = 2400000.5
# / Modified Julian Date of UNIX epoch (Jan. 1, 1970).
const UNIX_DT = DateTime(1970, 1, 1, 0, 0, 0)
const UNIX_JD = datetime2julian(UNIX_DT)

# / Seconds per half week.
const HALFWEEK = 302400
# / Seconds per whole week.
const FULLWEEK = 604800
# / Seconds per day.
const SEC_PER_DAY = 86400
# / Days per second.
const DAY_PER_SEC = 1.0 / SEC_PER_DAY

# / Milliseconds in a second.
const MS_PER_SEC = 1000
# / Seconds per millisecond.
const SEC_PER_MS = 1.0 / MS_PER_SEC

# / Milliseconds in a day.
const MS_PER_DAY = MS_PER_SEC * SEC_PER_DAY
# / Days per milliseconds.
const DAY_PER_MS = 1.0 / MS_PER_DAY

# system-specific constants

# GPS -------------------------------------------
# / 'Julian day' of GPS epoch (Jan. 6, 1980).
const GPS_EPOCH_DT = DateTime(1980, 1, 6, 0, 0, 0)
const GPS_EPOCH_JD = datetime2julian(GPS_EPOCH_DT)

# / Weeks per GPS Epoch
const GPS_WEEK_PER_EPOCH = 1024

# / Zcounts in a  day.
const ZCOUNT_PER_DAY = 57600
# / Days in a Zcount
const DAY_PER_ZCOUNT = 1.0 / ZCOUNT_PER_DAY
# / Zcounts in a week.
const ZCOUNT_PER_WEEK = 403200
# / Weeks in a Zcount.
const WEEK_PER_ZCOUNT = 1.0 / ZCOUNT_PER_WEEK
# / Z-counts per minute.
const ZCOUNT_PER_MINUTE = 40
# / Z-counts per hour.
const ZCOUNT_PER_HOUR = 2400

# GAL -------------------------------------------
# / 'Julian day' of GAL epoch (Aug 22 1999)
const GAL_EPOCH_DT = DateTime(1999, 8, 22, 0, 0, 0)
const GAL_EPOCH_JD = datetime2julian(GAL_EPOCH_DT)

const GAL_WEEK_PER_EPOCH = 4096

# QZS -------------------------------------------
# / 'Julian day' of QZS epoch (Jan. 1, 1980).
const QZS_EPOCH_DT = GPS_EPOCH_DT
const QZS_EPOCH_JD = GPS_EPOCH_JD
const QZS_WEEK_PER_EPOCH = 65535

# BDS -------------------------------------------
# / 'Julian day' of BDS epoch (Jan. 1, 2006).
const BDS_EPOCH_DT = DateTime(2006, 1, 1, 0, 0, 0)
const BDS_EPOCH_JD = datetime2julian(BDS_EPOCH_DT) # 2453736.5

const BDS_WEEK_PER_EPOCH = 8192

# IRN -------------------------------------------
# / 'Julian day' of IRN epoch (Aug 22, 1999).
const IRN_EPOCH_DT = DateTime(1999, 8, 22, 0, 0, 0)
const IRN_EPOCH_JD = datetime2julian(IRN_EPOCH_DT) # 2451412.5

const IRN_WEEK_PER_EPOCH = 1024


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Mortar3D.jl/blob/master/LICENSE

using Base.Test
using TimerOutputs
const to = TimerOutput()

test_files = String[]
push!(test_files, "test_calculate_normals.jl")
push!(test_files, "test_calculate_projections.jl")
push!(test_files, "test_calculate_polygon_clip.jl")

@testset "Mortar3D.jl" begin
    for fn in test_files
        timeit(to, fn) do
            include(fn)
        end
    end
end
println()
println("Test statistics:")
println(to)


text:
# code to generate R benchmark
# using RCall
# R"if(!require(pscl)){install.packages(\"pscl\");library(pscl)}"
# R"library(pscl)"
# R"data(\"bioChemists\", package = \"pscl\")"
# R"n <- dim(bioChemists)[1]"
# R"set.seed(13)"
# R"bioChemists$offpos = runif(n, min=0, max=3)"
# R"bioChemists$offzero = runif(n, min=0, max=3)"
# bioChemists=rcopy(R"bioChemists")
# writetable(joinpath(testdir,"data","bioChemists.csv"),bioChemists)

bioChemists=CSV.read(joinpath(testdir,"data","bioChemists.csv"))
bioChemists.marMarried=bioChemists.mar .== "Married"
bioChemists.femWomen=bioChemists.fem .== "Women"
bioChemists.art = convert(Array{Union{Float64, Missings.Missing},1}, bioChemists.art)

X=convert(Array{Float64,2},bioChemists[!,[:femWomen,:marMarried,:kid5,:phd,:ment]])
Xwconst=[ones(size(X,1)) X]
y=convert(Array{Float64,1},bioChemists.art)
const ixpartial = 50:60

###########################################################
# Xpos == Xzero
###########################################################
@testset "hurdle with Xpos == Xzero" begin

## logit-poisson
# R"fm_hp1 <- hurdle(art ~ fem + mar + kid5 + phd + ment, data = bioChemists)"
# print(R"summary(fm_hp1)")
# coefsR1=vec(rcopy(R"coef(fm_hp1, matrix = TRUE)"))
# yhatR1=vec(rcopy(R"predict(fm_hp1)"))
# yhatR1partial=vec(rcopy(R"predict(fm_hp1, newdata = bioChemists[ixpartial,])"))
# writetable(joinpath(testdir,"data","hurdle_coefsR1.csv"),DataFrame(coefsR=coefsR1))
# writetable(joinpath(testdir,"data","hurdle_yhatR1.csv"),DataFrame(yhatR1=yhatR1))
# writetable(joinpath(testdir,"data","hurdle_yhatR1partial.csv"),DataFrame(yhatR1partial=yhatR1partial))

coefsR1=vec(convert(Matrix{Float64},CSV.read(joinpath(testdir,"data","hurdle_coefsR1.csv"))))
yhatR1=vec(convert(Matrix{Float64},CSV.read(joinpath(testdir,"data","hurdle_yhatR1.csv"))))
yhatR1partial=vec(convert(Matrix{Float64},CSV.read(joinpath(testdir,"data","hurdle_yhatR1partial.csv"))))

# simple hurdle with GLM underlying
hurdlefit = fit(Hurdle,GeneralizedLinearModel,Xwconst,y)
showres = IOBuffer()
show(showres, hurdlefit)
showstr = String(take!(copy(showres)))
@test occursin("Hurdle regression",showstr)
@test occursin("Positive part coefficients",showstr)
@test occursin("Zero part coefficients",showstr)

coefsJ=vcat(coef(hurdlefit)...)
@test coefsJ ≈ coefsR1 rtol=1e-6
yhatJ=predict(hurdlefit, Xwconst)
@test yhatJ ≈ yhatR1 rtol=1e-6
yhatJpartial=predict(hurdlefit, Xwconst[ixpartial,:])
@test yhatJpartial ≈ yhatR1partial rtol=1e-6

# same simple hurdle through UNREGULATED lasso path
hurdleglrfit = fit(Hurdle,GammaLassoPath,X,y;λ=[0.0, 0.01])
coefsJ=vcat(coef(hurdleglrfit;select=MinAICc())...)
@test coefsJ ≈ coefsR1 rtol=1e-4
# rdist(coefsJ,coefsR1)
yhatJ = predict(hurdleglrfit, X; select=MinAICc())
@test yhatJ ≈ yhatR1 rtol=1e-4
yhatJpartial=predict(hurdleglrfit, X[ixpartial,:]; select=MinAICc())
@test yhatJpartial ≈ yhatR1partial rtol=1e-4
yhatJ = predict(hurdleglrfit, X; select=AllSeg())
@test size(yhatJ) == (size(X,1),2)

# regulated lasso path
hurdleglrfit = fit(Hurdle,GammaLassoPath,X,y; γ=0.0)
coefsJ=vcat(coef(hurdleglrfit;select=MinAICc())...)
@test coefsJ ≈ coefsR1 rtol=0.25
# rdist(coefsJ,coefsR1)
yhatJ = predict(hurdleglrfit, X; select=MinAICc())
@test yhatJ ≈ yhatR1 rtol=0.05
yhatJpartial=predict(hurdleglrfit, X[ixpartial,:]; select=MinAICc())
@test yhatJpartial ≈ yhatR1partial rtol=0.05

coefsJpos, coefsJzero = coef(hurdleglrfit;select=AllSeg())
@test size(coefsJpos,1) == 6
@test size(coefsJzero,1) == 6

showres = IOBuffer()
show(showres, hurdleglrfit)
showstr = String(take!(copy(showres)))
@test occursin("Hurdle regression",showstr)
@test occursin("Positive part regularization path",showstr)
@test occursin("Zero part regularization path",showstr)

# this one throws an error because we did not specify the same λ vector for both submodels so they have different lengths
@test_throws AssertionError predict(hurdleglrfit, X; select=AllSeg())

Random.seed!(1)
coefsJCVmin=vcat(coef(hurdleglrfit, MinCVKfold{MinCVmse}(5))...)
@test coefsJCVmin ≈ coefsR1 rtol=0.30
# rdist(coefsJCVmin,coefsR1)

# try with SharedArray
Xs = convert(SharedArray,X)
hurdleglrfitShared = fit(Hurdle,GammaLassoPath,Xs,y; γ=0.0)
coefsJShared=vcat(coef(hurdleglrfitShared;select=MinAICc())...)
@test coefsJShared ≈ coefsJ
yhatJ = predict(hurdleglrfitShared, X; select=MinAICc())
@test yhatJ ≈ yhatR1 rtol=0.05
yhatJpartial=predict(hurdleglrfitShared, X[ixpartial,:]; select=MinAICc())
@test yhatJpartial ≈ yhatR1partial rtol=0.05

# regulated gamma lasso path
hurdleglrfit = fit(Hurdle,GammaLassoPath,X,y; γ=10.0)
coefsJ=vcat(coef(hurdleglrfit;select=MinAICc())...)
@test coefsJ ≈ coefsR1 rtol=0.11
# rdist(coefsJ,coefsR1)
yhatJ = predict(hurdleglrfit, X; select=MinAICc())
@test yhatJ ≈ yhatR1 rtol=0.05
yhatJpartial=predict(hurdleglrfit, X[ixpartial,:]; select=MinAICc())
@test yhatJpartial ≈ yhatR1partial rtol=0.05

# using DataFrames formula interface
hurdlefit = fit(Hurdle,GeneralizedLinearModel,@formula(art ~ femWomen + marMarried + kid5 + phd + ment), bioChemists)
coefsJ=vcat(coef(hurdlefit)...)
@test coefsJ ≈ coefsR1 rtol=1e-6
yhatJ=predict(hurdlefit, Xwconst)
@test yhatJ ≈ yhatR1 rtol=1e-6
yhatJpartial=predict(hurdlefit, Xwconst[ixpartial,:])
@test yhatJpartial ≈ yhatR1partial rtol=1e-6

end

###########################################################
# Xpos ≠ Xzero
###########################################################
@testset "hurdle with Xpos ≠ Xzero" begin

# regulated gamma lasso path with different Xpos and Xzero
# R"fm_hp2 <- hurdle(art ~ fem + mar + kid5 | phd + ment, data = bioChemists)"
# print(R"summary(fm_hp2)")
# coefsR2=vec(rcopy(R"coef(fm_hp2, matrix = TRUE)"))
# yhatR2=vec(rcopy(R"predict(fm_hp2)"))
# yhatR2partial=vec(rcopy(R"predict(fm_hp2, newdata = bioChemists[ixpartial,])"))
# writetable(joinpath(testdir,"data","hurdle_coefsR2.csv"),DataFrame(coefsR=coefsR2))
# writetable(joinpath(testdir,"data","hurdle_yhatR2.csv"),DataFrame(yhatR2=yhatR2))
# writetable(joinpath(testdir,"data","hurdle_yhatR2partial.csv"),DataFrame(yhatR2partial=yhatR2partial))

coefsR2=vec(convert(Matrix{Float64},CSV.read(joinpath(testdir,"data","hurdle_coefsR2.csv"))))
yhatR2=vec(convert(Matrix{Float64},CSV.read(joinpath(testdir,"data","hurdle_yhatR2.csv"))))
yhatR2partial=vec(convert(Matrix{Float64},CSV.read(joinpath(testdir,"data","hurdle_yhatR2partial.csv"))))

Xpos = X[:,1:3]
Xzero = X[:,4:5]
Xzerowconst=[ones(size(X,1)) Xzero]
Xposwconst=[ones(size(X,1)) Xpos]
# ixpos = y.>0
# ypos = y[ixpos]
# countmap(ypos)

# simple hurdle with GLM underlying
hurdlefit = fit(Hurdle,GeneralizedLinearModel,Xzerowconst,y; Xpos=Xposwconst)
coefsJ=vcat(coef(hurdlefit)...)
@test coefsJ ≈ coefsR2 rtol=1e-5
yhatJ=predict(hurdlefit, Xzerowconst; Xpos=Xposwconst)
@test yhatJ ≈ yhatR2 rtol=1e-6
yhatJpartial=predict(hurdlefit, Xzerowconst[ixpartial,:]; Xpos=Xposwconst[ixpartial,:])
@test yhatJpartial ≈ yhatR2partial rtol=1e-6

# same simple hurdle through UNREGULATED lasso path
hurdleglrfit = fit(Hurdle,GammaLassoPath,Xzero,y; Xpos=Xpos, λ=[0.0,0.0001])
coefsJ=vcat(coef(hurdleglrfit;select=MinAICc())...)
@test coefsJ ≈ coefsR2 rtol=1e-4
# rdist(coefsJ,coefsR2)
yhatJ=predict(hurdleglrfit, Xzero; Xpos=Xpos, select=MinAICc())
@test yhatJ ≈ yhatR2 rtol=1e-4
yhatJpartial=predict(hurdleglrfit, Xzero[ixpartial,:]; Xpos=Xpos[ixpartial,:], select=MinAICc())
@test yhatJpartial ≈ yhatR2partial rtol=1e-4
yhatJ = predict(hurdleglrfit, Xzero; Xpos=Xpos, select=AllSeg())
@test size(yhatJ) == (size(X,1),2)

# regulated lasso path
hurdleglrfit = fit(Hurdle,GammaLassoPath,Xzero,y; Xpos=Xpos, γ=0.0)
coefsJ=vcat(coef(hurdleglrfit;select=MinAICc())...)
@test coefsJ ≈ coefsR2 rtol=0.10
# rdist(coefsJ,coefsR2)
yhatJ = predict(hurdleglrfit, Xzero; Xpos=Xpos, select=MinAICc())
@test yhatJ ≈ yhatR2 rtol=0.05
yhatJpartial=predict(hurdleglrfit, Xzero[ixpartial,:]; Xpos=Xpos[ixpartial,:], select=MinAICc())
@test yhatJpartial ≈ yhatR2partial rtol=0.05

coefsJpos, coefsJzero = coef(hurdleglrfit;select=AllSeg())
@test size(coefsJpos,1) == 4
@test size(coefsJzero,1) == 3

# try with SharedArray
Xzeros = convert(SharedArray,Xzero)
Xposs = convert(SharedArray,Xpos)
hurdleglrfitShared = fit(Hurdle,GammaLassoPath,Xzeros,y; Xpos=Xposs, γ=0.0)
coefsJShared=vcat(coef(hurdleglrfitShared;select=MinAICc())...)
@test coefsJShared ≈ coefsJ
yhatJ = predict(hurdleglrfitShared, Xzeros; Xpos=Xposs, select=MinAICc())
@test yhatJ ≈ yhatR2 rtol=0.05
yhatJpartial=predict(hurdleglrfitShared, Xzeros[ixpartial,:]; Xpos=Xposs[ixpartial,:], select=MinAICc())
@test yhatJpartial ≈ yhatR2partial rtol=0.05

# regulated gamma lasso path
hurdleglrfit = fit(Hurdle,GammaLassoPath,Xzero,y; Xpos=Xpos, γ=10.0)
coefsJ=vcat(coef(hurdleglrfit;select=MinAICc())...)
@test coefsJ ≈ coefsR2 rtol=0.15
# rdist(coefsJ,coefsR2)
yhatJ = predict(hurdleglrfit, Xzero; Xpos=Xpos, select=MinAICc())
@test yhatJ ≈ yhatR2 rtol=0.05
yhatJpartial=predict(hurdleglrfit, Xzero[ixpartial,:]; Xpos=Xpos[ixpartial,:], select=MinAICc())
@test yhatJpartial ≈ yhatR2partial rtol=0.05

# using DataFrames formula interface
hurdlefit = fit(Hurdle,GeneralizedLinearModel,@formula(art ~ phd + ment), bioChemists; fpos = @formula(art ~ femWomen + marMarried + kid5))
coefsJ=vcat(coef(hurdlefit)...)
@test coefsJ ≈ coefsR2 rtol=1e-5
# rdist(coefsJ,coefsR2)
yhatJ=predict(hurdlefit, Xzerowconst; Xpos=Xposwconst)
@test yhatJ ≈ yhatR2 rtol=1e-6
yhatJpartial=predict(hurdlefit, Xzerowconst[ixpartial,:]; Xpos=Xposwconst[ixpartial,:])
@test yhatJpartial ≈ yhatR2partial rtol=1e-6

end

###########################################################
# Xpos ≠ Xzero
###########################################################
@testset "hurdle with Xpos ≠ Xzero AND offset specified" begin

# regulated gamma lasso path with different Xpos and Xzero and an offset
# R"fm_hp3 <- hurdle(art ~ fem + mar + kid5 + offset(offpos) | phd + ment + offset(offzero), data = bioChemists)"
# print(R"summary(fm_hp3)")
# coefsR3=vec(rcopy(R"coef(fm_hp3, matrix = TRUE)"))
# yhatR3=vec(rcopy(R"predict(fm_hp3)"))
# yhatR3partial=vec(rcopy(R"predict(fm_hp3, newdata = bioChemists[ixpartial,])"))
# writetable(joinpath(testdir,"data","hurdle_coefsR3.csv"),DataFrame(coefsR=coefsR3))
# writetable(joinpath(testdir,"data","hurdle_yhatR3.csv"),DataFrame(yhatR3=yhatR3))
# writetable(joinpath(testdir,"data","hurdle_yhatR3partial.csv"),DataFrame(yhatR3partial=yhatR3partial))

coefsR3=vec(convert(Matrix{Float64},CSV.read(joinpath(testdir,"data","hurdle_coefsR3.csv"))))
yhatR3=vec(convert(Matrix{Float64},CSV.read(joinpath(testdir,"data","hurdle_yhatR3.csv"))))
yhatR3partial=vec(convert(Matrix{Float64},CSV.read(joinpath(testdir,"data","hurdle_yhatR3partial.csv"))))

offpos = convert(Vector{Float64},bioChemists.offpos)
offzero = convert(Vector{Float64},bioChemists.offzero)
Xpos = X[:,1:3]
Xzero = X[:,4:5]
Xzerowconst=[ones(size(X,1)) Xzero]
Xposwconst=[ones(size(X,1)) Xpos]
# ixpos = y.>0
# ypos = y[ixpos]
# countmap(ypos)

# simple hurdle with GLM underlying
hurdlefit = fit(Hurdle,GeneralizedLinearModel,Xzerowconst,y; Xpos=Xposwconst, offsetzero=offzero, offsetpos=offpos)
coefsJ=vcat(coef(hurdlefit)...)
@test coefsJ ≈ coefsR3 rtol=1e-5
yhatJ=predict(hurdlefit, Xzerowconst; Xpos=Xposwconst, offsetzero=offzero, offsetpos=offpos)
@test yhatJ ≈ yhatR3 rtol=1e-5
yhatJpartial=predict(hurdlefit, Xzerowconst[ixpartial,:]; Xpos=Xposwconst[ixpartial,:], offsetzero=offzero[ixpartial], offsetpos=offpos[ixpartial])
@test yhatJpartial ≈ yhatR3partial rtol=1e-6

# same simple hurdle through UNREGULATED lasso path
hurdleglrfit = fit(Hurdle,GammaLassoPath,Xzero,y; Xpos=Xpos, offsetzero=offzero, offsetpos=offpos, λ=[0.0,0.0001])
coefsJ=vcat(coef(hurdleglrfit;select=MinAICc())...)
@test coefsJ ≈ coefsR3 rtol=1e-4
# rdist(coefsJ,coefsR3)
yhatJ=predict(hurdleglrfit, Xzero; Xpos=Xpos, offsetzero=offzero, offsetpos=offpos, select=MinAICc())
@test yhatJ ≈ yhatR3 rtol=1e-4
yhatJpartial=predict(hurdleglrfit, Xzero[ixpartial,:]; Xpos=Xpos[ixpartial,:], offsetzero=offzero[ixpartial], offsetpos=offpos[ixpartial], select=MinAICc())
@test yhatJpartial ≈ yhatR3partial rtol=1e-4
yhatJ = predict(hurdleglrfit, Xzero; Xpos=Xpos, offsetzero=offzero, offsetpos=offpos, select=AllSeg())
@test size(yhatJ) == (size(X,1),2)

# regulated lasso path
hurdleglrfit = fit(Hurdle,GammaLassoPath,Xzero,y; Xpos=Xpos, offsetzero=offzero, offsetpos=offpos, γ=0.0)
coefsJ=vcat(coef(hurdleglrfit;select=MinAICc())...)
@test coefsJ ≈ coefsR3 rtol=0.10
# rdist(coefsJ,coefsR3)
yhatJ = predict(hurdleglrfit, Xzero; Xpos=Xpos, offsetzero=offzero, offsetpos=offpos, select=MinAICc())
@test yhatJ ≈ yhatR3 rtol=0.05
yhatJpartial=predict(hurdleglrfit, Xzero[ixpartial,:]; Xpos=Xpos[ixpartial,:], offsetzero=offzero[ixpartial], offsetpos=offpos[ixpartial], select=MinAICc())
@test yhatJpartial ≈ yhatR3partial rtol=0.05

coefsJpos, coefsJzero = coef(hurdleglrfit;select=AllSeg())
@test size(coefsJpos,1) == 4
@test size(coefsJzero,1) == 3

# try with SharedArray
Xzeros = convert(SharedArray,Xzero)
Xposs = convert(SharedArray,Xpos)
offzeros = convert(SharedArray,offzero)
offposs = convert(SharedArray,offpos)
hurdleglrfitShared = fit(Hurdle,GammaLassoPath,Xzeros,y; Xpos=Xposs, offsetzero=offzeros, offsetpos=offposs, γ=0.0)
coefsJShared=vcat(coef(hurdleglrfitShared;select=MinAICc())...)
@test coefsJShared ≈ coefsJ
yhatJ = predict(hurdleglrfitShared, Xzeros; Xpos=Xposs, offsetzero=offzeros, offsetpos=offposs, select=MinAICc())
@test yhatJ ≈ yhatR3 rtol=0.05
yhatJpartial=predict(hurdleglrfitShared, Xzeros[ixpartial,:]; Xpos=Xposs[ixpartial,:], offsetzero=offzeros[ixpartial], offsetpos=offposs[ixpartial], select=MinAICc())
@test yhatJpartial ≈ yhatR3partial rtol=0.05

# regulated gamma lasso path
hurdleglrfit = fit(Hurdle,GammaLassoPath,Xzero,y; Xpos=Xpos, offsetzero=offzero, offsetpos=offpos, γ=10.0)
coefsJ=vcat(coef(hurdleglrfit;select=MinAICc())...)
@test coefsJ ≈ coefsR3 rtol=0.15
# rdist(coefsJ,coefsR3)
yhatJ = predict(hurdleglrfit, Xzero; Xpos=Xpos, offsetzero=offzero, offsetpos=offpos, select=MinAICc())
@test yhatJ ≈ yhatR3 rtol=0.05
yhatJpartial=predict(hurdleglrfit, Xzero[ixpartial,:]; Xpos=Xpos[ixpartial,:], offsetzero=offzero[ixpartial], offsetpos=offpos[ixpartial], select=MinAICc())
@test yhatJpartial ≈ yhatR3partial rtol=0.05

# using DataFrames formula interface
hurdlefit = fit(Hurdle,GeneralizedLinearModel,@formula(art ~ phd + ment), bioChemists; fpos = @formula(art ~ femWomen + marMarried + kid5), offsetzero=offzero, offsetpos=offpos)
coefsJ=vcat(coef(hurdlefit)...)
@test coefsJ ≈ coefsR3 rtol=1e-5
# rdist(coefsJ,coefsR3)
yhatJ=predict(hurdlefit, Xzerowconst; Xpos=Xposwconst, offsetzero=offzero, offsetpos=offpos)
@test yhatJ ≈ yhatR3 rtol=1e-5
yhatJpartial=predict(hurdlefit, Xzerowconst[ixpartial,:]; Xpos=Xposwconst[ixpartial,:], offsetzero=offzero[ixpartial], offsetpos=offpos[ixpartial])
@test yhatJpartial ≈ yhatR3partial rtol=1e-5

end

###########################################################
# degenrate cases
###########################################################
@testset "hurdle degenerate cases" begin

# degenerate positive counts data case 1
include(joinpath(testdir,"data","degenerate_hurdle_1.jl"))
hurdle = @test_logs (:warn, r"failed to fit truncated counts model to positive") fit(Hurdle,GammaLassoPath,X,y; showwarnings=true)
coefsJ=vcat(coef(hurdle;select=MinAICc())...)
@test vec(coefsJ) ≈ [0.0, 0.0, -6.04112, 0.675767] rtol=1e-4
coefsJpos, coefsJzero = coef(hurdle;select=AllSeg())
@test size(coefsJpos,1) == 2
@test coefsJpos == zero(coefsJpos)
@test size(coefsJzero,1) == 2

ydeg = zero(y)
ydeg[1] = 2
ydeg[2] = 3
hurdleglm = @test_logs (:warn, r"failed to fit truncated counts model to positive") fit(Hurdle,GeneralizedLinearModel,[ones(size(X,1)) X],ydeg; showwarnings=true)

# degenerate positive counts data case 1 without >1
y0or1 = deepcopy(y)
y0or1[y.>1] .= 1
hurdle = @test_logs (:warn, r"ypos has no elements larger than 1") fit(Hurdle,GammaLassoPath,X,y0or1; showwarnings=true)
coefs0or1=vcat(coef(hurdle;select=MinAICc())...)
@test coefs0or1 == coefsJ
coefs0or1pos, coefs0or1zero = coef(hurdle;select=AllSeg())
@test coefs0or1pos == coefsJpos
@test coefs0or1zero == coefsJzero

@info("Testing hurdle degenerate cases. The following warnings about step-halving are expected ...")

# degenerate positive counts data case 1 without zeros
y0or1 = deepcopy(y)
y0or1[y.==0] .= 1
y0or1[1] = 3.0

hurdle = @test_logs (:warn, r"I\(y\) is all ones") fit(Hurdle,GammaLassoPath,X,y0or1; verbose=true, showwarnings=true)
coefs0or1=vcat(coef(hurdle;select=MinAICc())...)
@test vec(coefs0or1) ≈ [-10.2903, 0.820934, 0.0, 0.0] rtol=1e-4
coefs0or1pos, coefs0or1zero = coef(hurdle;select=AllSeg())
@test all(iszero,coefs0or1zero)

hurdleglm = @test_logs (:warn, r"I\(y\) is all ones") fit(Hurdle,GeneralizedLinearModel,[ones(size(X,1)) X],y0or1; verbose=true, showwarnings=true)

# degenerate positive counts data case 1 with only zeros
y0or1 = zero(y)
@test_throws ErrorException fit(Hurdle,GammaLassoPath,X,y0or1)
@test_logs (:warn, r"I\(y\) is all zeros") @test_throws ErrorException fit(Hurdle,GammaLassoPath,X,y0or1; verbose=true, showwarnings=true)

# degenerate positive counts data case 2
include(joinpath(testdir,"data","degenerate_hurdle_2.jl"))
hurdle = @test_logs (:warn, r"not enough variation in ypos") fit(Hurdle,GammaLassoPath,X,y; showwarnings=true)
coefsJ=vcat(coef(hurdle;select=MinAICc())...)
@test vec(coefsJ) ≈ [0.0,0.0,-5.30128195796556,0.1854148891565171] rtol=1e-4
coefsJpos, coefsJzero = coef(hurdle;select=AllSeg())
@test size(coefsJpos,1) == 2
@test coefsJpos == zero(coefsJpos)
@test size(coefsJzero,1) == 2

# degenerate positive counts data case 3
include(joinpath(testdir,"data","degenerate_hurdle_3.jl"))
hurdle = @test_logs (:warn, r"ypos has no elements larger than 1") fit(Hurdle,GammaLassoPath,X,y; showwarnings=true)
coefsJ=vcat(coef(hurdle;select=MinAICc())...)
# @test vec(coefsJ) ≈ [0.0, 0.0, -4.54363, 0.000458273] rtol=1e-4
@test length(coefsJ) == 4
coefsJpos, coefsJzero = coef(hurdle;select=AllSeg())
Matrix( coefsJzero)
@test size(coefsJpos,1) == 2
@test coefsJpos == zero(coefsJpos)
@test size(coefsJzero,1) == 2

include(joinpath(testdir,"data","degenerate_hurdle_4.jl"))
hurdle = @test_nowarn fit(Hurdle,GammaLassoPath,X,y; offset=offset)
# TODO: figure out why next test behaves differently when running tests manually
# as include("runtests.jl") or as ]test HurdleDMR
# hurdle = @test_logs (:warn, r"lambertw") (:warn, r"lambertw") fit(Hurdle,GammaLassoPath,X,y; offset=offset, showwarnings=true)

# this test case used to give numerical headaches to devresid(PositivePoisson(),...)
include(joinpath(testdir,"data","degenerate_hurdle_5.jl"))
hurdle = fit(Hurdle,GammaLassoPath,X,y; Xpos=Xpos, offset=offset)
@test !(ismissing(hurdle.mpos))
@test !(ismissing(hurdle.mzero))
@test typeof(hurdle.mpos) <: GammaLassoPath
@test typeof(hurdle.mzero) <: GammaLassoPath

end


text:
export DCT, IDCT

abstract type CosineTransform{N,C,T1,T2} <: LinearOperator end

"""
`DCT([domainType=Float64::Type,] dim_in::Tuple)`

`DCT(dim_in...)`

`DCT(x::AbstractArray)`

Creates a `LinearOperator` which, when multiplied with an array `x::AbstractArray{N}`, returns the `N`-dimensional Inverse Discrete Cosine Transform of `x`.

```julia
julia> DCT(Complex{Float64},(10,10))
ℱc  ℂ^(10, 10) -> ℂ^(10, 10)

julia> DCT(10,10)
ℱc  ℝ^(10, 10) -> ℂ^(10, 10)

julia> A = DCT(ones(3))
ℱc  ℝ^3 -> ℝ^3

julia> A*ones(3)
3-element Array{Float64,1}:
 1.73205
 0.0
 0.0

```

"""
struct DCT{N,
	   C<:RealOrComplex,
	   T1<:AbstractFFTs.Plan,
	   T2<:AbstractFFTs.Plan} <: CosineTransform{N,C,T1,T2}
	dim_in::NTuple{N,Int}
	A::T1
	At::T2
end

"""
`IDCT([domainType=Float64::Type,] dim_in::Tuple)`

`IDCT(dim_in...)`

`IDCT(x::AbstractArray)`

Creates a `LinearOperator` which, when multiplied with an array `x::AbstractArray{N}`, returns the `N`-dimensional Discrete Cosine Transform of `x`.

```julia
julia> IDCT(Complex{Float64},(10,10))
ℱc⁻¹  ℂ^(10, 10) -> ℂ^(10, 10)

julia> IDCT(10,10)
ℱc⁻¹  ℝ^(10, 10) -> ℂ^(10, 10)

julia> A = IDCT(ones(3))
ℱc⁻¹  ℝ^3 -> ℝ^3

julia> A*[1.;0.;0.]
3-element Array{Float64,1}:
 0.57735
 0.57735
 0.57735

```

"""
struct IDCT{N,
	    C<:RealOrComplex,
	    T1<:AbstractFFTs.Plan,
	    T2<:AbstractFFTs.Plan} <: CosineTransform{N,C,T1,T2}
	dim_in::NTuple{N,Int}
	A::T1
	At::T2
end

# Constructors
#standard constructor
DCT(T::Type,dim_in::NTuple{N,Int}) where {N} = DCT(zeros(T,dim_in))

DCT(dim_in::NTuple{N,Int}) where {N} = DCT(zeros(dim_in))
DCT(dim_in::Vararg{Int64}) = DCT(dim_in)
DCT(T::Type, dim_in::Vararg{Int64}) = DCT(T,dim_in)

function DCT(x::AbstractArray{C,N}) where {N,C<:RealOrComplex}
	A,At = plan_dct(x), plan_idct(x)
	DCT{N,C,typeof(A),typeof(At)}(size(x),A,At)
end

#standard constructor
IDCT(T::Type,dim_in::NTuple{N,Int}) where {N} = IDCT(zeros(T,dim_in))

IDCT(dim_in::NTuple{N,Int}) where {N} = IDCT(zeros(dim_in))
IDCT(dim_in::Vararg{Int64}) = IDCT(dim_in)
IDCT(T::Type, dim_in::Vararg{Int64}) = IDCT(T,dim_in)

function IDCT(x::AbstractArray{C,N}) where {N,C<:RealOrComplex}
	A,At = plan_idct(x), plan_dct(x)
	IDCT{N,C,typeof(A),typeof(At)}(size(x),A,At)
end

# Mappings

function mul!(y::AbstractArray{C,N},A::DCT{N,C,T1,T2},b::AbstractArray{C,N}) where {N,C,T1,T2}
	mul!(y,A.A,b)
end

function mul!(y::AbstractArray{C,N},A::AdjointOperator{DCT{N,C,T1,T2}},b::AbstractArray{C,N}) where {N,C,T1,T2}
	y .= A.A.At*b
end

function mul!(y::AbstractArray{C,N},A::IDCT{N,C,T1,T2},b::AbstractArray{C,N}) where {N,C,T1,T2}
	y .= A.A*b
end

function mul!(y::AbstractArray{C,N},A::AdjointOperator{IDCT{N,C,T1,T2}},b::AbstractArray{C,N}) where {N,C,T1,T2}
	mul!(y,A.A.At,b)
end

# Properties

size(L::CosineTransform) = (L.dim_in,L.dim_in)

fun_name(A::DCT)  = "ℱc"
fun_name(A::IDCT) = "ℱc⁻¹"

domainType(L::CosineTransform{N,C,T1,T2}) where {N,C,T1,T2} = C
codomainType(L::CosineTransform{N,C,T1,T2}) where {N,C,T1,T2} = C

is_AcA_diagonal(L::CosineTransform)     = true
is_AAc_diagonal(L::CosineTransform)     = true
is_orthogonal(  L::CosineTransform)     = true
is_invertible(L::CosineTransform)       = true
is_full_row_rank(L::CosineTransform)    = true
is_full_column_rank(L::CosineTransform) = true

diag_AcA(L::CosineTransform) = 1.
diag_AAc(L::CosineTransform) = 1.


text:
abstract type AbstractRenewableFormulation <: AbstractDeviceFormulation end
abstract type AbstractRenewableDispatchFormulation <: AbstractRenewableFormulation end
struct RenewableFullDispatch <: AbstractRenewableDispatchFormulation end
struct RenewableConstantPowerFactor <: AbstractRenewableDispatchFormulation end

########################### renewable generation variables #################################
function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    ::PSIContainer,
) where {T <: ActivePowerVariable, U <: PSY.RenewableGen}
    return AddVariableSpec(;
        variable_name = make_variable_name(T, U),
        binary = false,
        expression_name = :nodal_balance_active,
        lb_value_func = x -> 0.0,
        ub_value_func = x -> PSY.get_max_active_power(x),
    )
end

function AddVariableSpec(
    ::Type{T},
    ::Type{U},
    ::PSIContainer,
) where {T <: ReactivePowerVariable, U <: PSY.RenewableGen}
    return AddVariableSpec(;
        variable_name = make_variable_name(T, U),
        binary = false,
        expression_name = :nodal_balance_reactive,
    )
end

####################################### Reactive Power constraint_infos #########################
function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ReactivePowerVariable},
    ::Type{T},
    ::Type{<:AbstractDeviceFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.RenewableGen}
    return DeviceRangeConstraintSpec(;
        range_constraint_spec = RangeConstraintSpec(;
            constraint_name = make_constraint_name(
                RangeConstraint,
                ReactivePowerVariable,
                T,
            ),
            variable_name = make_variable_name(ReactivePowerVariable, T),
            limits_func = x -> PSY.get_reactive_power_limits(x),
            constraint_func = device_range!,
            constraint_struct = DeviceRangeConstraintInfo,
        ),
    )
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ReactivePowerVariable},
    ::Type{T},
    ::Type{<:RenewableConstantPowerFactor},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.RenewableGen}
    return DeviceRangeConstraintSpec(;
        custom_psi_container_func = custom_reactive_power_constraints!,
    )
end

function custom_reactive_power_constraints!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::Type{RenewableConstantPowerFactor},
) where {T <: PSY.RenewableGen}
    names = [PSY.get_name(d) for d in devices]
    time_steps = model_time_steps(psi_container)
    p_var = get_variable(psi_container, ACTIVE_POWER, T)
    q_var = get_variable(psi_container, REACTIVE_POWER, T)
    constraint_val = JuMPConstraintArray(undef, names, time_steps)
    assign_constraint!(psi_container, REACTIVE_RANGE, T, constraint_val)
    for t in time_steps, d in devices
        name = PSY.get_name(d)
        pf = sin(acos(PSY.get_power_factor(d)))
        constraint_val[name, t] =
            JuMP.@constraint(psi_container.JuMPmodel, q_var[name, t] == p_var[name, t] * pf)
    end
    return
end

function DeviceRangeConstraintSpec(
    ::Type{<:RangeConstraint},
    ::Type{ActivePowerVariable},
    ::Type{T},
    ::Type{<:AbstractRenewableDispatchFormulation},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward},
    use_parameters::Bool,
    use_forecasts::Bool,
) where {T <: PSY.RenewableGen}
    if !use_parameters && !use_forecasts
        return DeviceRangeConstraintSpec(;
            range_constraint_spec = RangeConstraintSpec(;
                constraint_name = make_constraint_name(
                    RangeConstraint,
                    ActivePowerVariable,
                    T,
                ),
                variable_name = make_variable_name(ActivePowerVariable, T),
                limits_func = x -> (min = 0.0, max = PSY.get_active_power(x)),
                constraint_func = device_range!,
                constraint_struct = DeviceRangeConstraintInfo,
            ),
        )
    end

    return DeviceRangeConstraintSpec(;
        timeseries_range_constraint_spec = TimeSeriesConstraintSpec(;
            constraint_name = make_constraint_name(RangeConstraint, ActivePowerVariable, T),
            variable_name = make_variable_name(ActivePowerVariable, T),
            parameter_name = use_parameters ? ACTIVE_POWER : nothing,
            forecast_label = "max_active_power",
            multiplier_func = x -> PSY.get_max_active_power(x),
            constraint_func = use_parameters ? device_timeseries_param_ub! :
                              device_timeseries_ub!,
        ),
    )
end

########################## Addition to the nodal balances ##################################

function NodalExpressionSpec(
    ::Type{T},
    ::Type{<:PM.AbstractPowerModel},
    use_forecasts::Bool,
) where {T <: PSY.RenewableGen}
    return NodalExpressionSpec(
        "max_active_power",
        REACTIVE_POWER,
        use_forecasts ? x -> PSY.get_max_reactive_power(x) : x -> PSY.get_reactive_power(x),
        1.0,
        T,
    )
end

function NodalExpressionSpec(
    ::Type{T},
    ::Type{<:PM.AbstractActivePowerModel},
    use_forecasts::Bool,
) where {T <: PSY.RenewableGen}
    return NodalExpressionSpec(
        "max_active_power",
        ACTIVE_POWER,
        use_forecasts ? x -> PSY.get_max_active_power(x) : x -> PSY.get_active_power(x),
        1.0,
        T,
    )
end

##################################### renewable generation cost ############################
function AddCostSpec(
    ::Type{T},
    ::Type{U},
    ::PSIContainer,
) where {T <: PSY.RenewableDispatch, U <: AbstractRenewableDispatchFormulation}
    # TODO: remove once cost_function is required
    cost_function = x -> isnothing(x) ? 1.0 : PSY.get_variable(x)
    return AddCostSpec(;
        variable_type = ActivePowerVariable,
        component_type = T,
        variable_cost = cost_function,
        multiplier = OBJECTIVE_FUNCTION_NEGATIVE,
    )
end


text:
## Cassette Powered Intervention
## =============================

@context ChangeCtx

function Cassette.execute(ctx::ChangeCtx, x::RandVar, ω::Ω)
  if ctx.metadata.id === x.id
    return ctx.metadata.x(ω)
  else
    return Cassette.RecurseInstead()
  end
end

"Causal Intervention: Set `θold` to `θnew` in `x`"
function change(θold::RandVar, θnew::RandVar, x::RandVar{T}) where T
  f = ω -> Cassette.overdub(ChangeCtx(metadata = (id = θold.id, x = θnew)), x, ω)
  RandVar{T}(f)
end

"Change where `θconst` is not a randvar, but constant"
change(θold, θconst, x) = change(θold, constant(θconst), x)


text:
using Anamnesis
using BenchmarkTools
using Random, LinearAlgebra

dtype = ComplexF64

callcount = [0, 0]
f = (x::AbstractVector{<:Number}) -> (callcount[1] += 1; x⋅x + 1)
sf = Scribe(f, IdDict())

x = rand(dtype, rand(64:256))
ξ = rand(dtype)
y = randstring(rand(0:31))

sf(x) == f(x)

# sg(ξ, y=y) == g(ξ, y=y)
# sf(x) == f(x)
# sg(ξ, y=y) == g(ξ, y=y)

ex = @macroexpand @anamnesis begin
        f2([1,2]) + h2(3)
        z = f2([1,2]) + h2(3)
        string(g2(1, y="fire"), g2(2, y="walk"), h2("with"), h2("me"))
        A = [h2("with") g2(1,y="fire")
             g2(2,y="walk") h2("me")]
end f2 g2


text:
using WIOD_brexit
using Test

@testset "WIOD_brexit.jl" begin
    # Write your tests here.

    #@test round(raw_long.value[1]) == 12924.0 ---- does not work, need to use something without external data

end


text:
module ZeroMeanFESpacesTests

using Test
using Gridap.Geometry
using Gridap.Fields
using Gridap.FESpaces
using Gridap.CellData
using Gridap.ReferenceFEs

domain = (0,1,0,1)
partition = (4,4)
model = CartesianDiscreteModel(domain,partition)

order = 2

trian = get_triangulation(model)
degree = order
dΩ = Measure(trian,degree)

_V = FESpace(model,ReferenceFE(lagrangian,Float64,order);conformity=:L2)

V = ZeroMeanFESpace(_V,dΩ)

matvecdata = ([],[],[])
matdata = ([],[],[])
vecdata = ([],[])
test_single_field_fe_space(V,matvecdata,matdata,vecdata)

U = TrialFESpace(V)
test_single_field_fe_space(U,matvecdata,matdata,vecdata)
@test isa(U,ZeroMeanFESpace)

fun(x) = sin(4*pi*(x[1]+x[2]^2)) + 3
uh = interpolate(fun, U)

mean1 = sum(∫(uh)*dΩ)

tol = 1.0e-10
@test abs(mean1) < tol

V = FESpace(model,ReferenceFE(lagrangian,Float64,order);conformity=:L2,constraint=:zeromean)
@test isa(V,ZeroMeanFESpace)

end # module


text:
using ReachabilityModels
using Test

@testset "ReachabilityModels.jl" begin
    # Write your own tests here.
end

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
"""
    moon98(date1, date2)

Approximate geocentric position and velocity of the Moon.

n.b. Not IAU-endorsed and without canonical status.

### Given ###

- `date1`: TT date part A (Notes 1,4)
- `date2`: TT date part B (Notes 1,4)

### Returned ###

- `pv`: Moon p,v, GCRS (AU, AU/d, Note 5)

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways, among
   others:

   | `date1`   |     `date2` | Method      |
   |:----------|:------------|:------------|
   | 2450123.7 |         0.0 | JD          |
   | 2451545.0 |     -1421.3 | J2000       |
   | 2400000.5 |     50123.2 | MJD         |
   | 2450123.5 |         0.2 | date & time |

   The JD method is the most natural and convenient to use in cases
   where the loss of several decimal digits of resolution is
   acceptable.  The J2000 method is best matched to the way the
   argument is handled internally and will deliver the optimum
   resolution.  The MJD method and the date & time methods are both
   good compromises between resolution and convenience.  The limited
   accuracy of the present algorithm is such that any of the methods
   is satisfactory.

2. This function is a full implementation of the algorithm
   published by Meeus (see reference) except that the light-time
   correction to the Moon's mean longitude has been omitted.

3. Comparisons with ELP/MPP02 over the interval 1950-2100 gave RMS
   errors of 2.9 arcsec in geocentric direction, 6.1 km in position
   and 36 mm/s in velocity.  The worst case errors were 18.3 arcsec
   in geocentric direction, 31.7 km in position and 172 mm/s in
   velocity.

4. The original algorithm is expressed in terms of "dynamical time",
   which can either be TDB or TT without any significant change in
   accuracy.  UT cannot be used without incurring significant errors
   (30 arcsec in the present era) due to the Moon's 0.5 arcsec/sec
   movement.

5. The result is with respect to the GCRS (the same as J2000.0 mean
   equator and equinox to within 23 mas).

6. Velocity is obtained by a complete analytical differentiation
   of the Meeus model.

7. The Meeus algorithm generates position and velocity in mean
   ecliptic coordinates of date, which the present function then
   rotates into GCRS.  Because the ecliptic system is precessing,
   there is a coupling between this spin (about 1.4 degrees per
   century) and the Moon position that produces a small velocity
   contribution.  In the present function this effect is neglected
   as it corresponds to a maximum difference of less than 3 mm/s and
   increases the RMS error by only 0.4%.

### References ###

- Meeus, J., Astronomical Algorithms, 2nd edition, Willmann-Bell,
  1998, p337.

- Simon, J.L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
  Francou, G. & Laskar, J., Astron.Astrophys., 1994, 282, 663

### Called ###

- [`s2pv`](@ref): spherical coordinates to pv-vector
- [`pfw06`](@ref): bias-precession F-W angles, IAU 2006
- [`ir`](@ref): initialize r-matrix to identity
- [`rz`](@ref): rotate around Z-axis
- [`rx`](@ref): rotate around X-axis
- [`rxpv`](@ref): product of r-matrix and pv-vector
"""
function moon98(date1, date2)
    pv = Array{Cdouble}(undef, 3, 2)
    ccall((:eraMoon98, liberfa), Cvoid,
          (Cdouble, Cdouble, Ptr{Cdouble}),
          date1, date2, pv)
    return cmatrix_to_array(pv)
end


text:
@doc doc"""
Enumtype including `LET_*`
"""
@enum(LoopElementType,
      LET_Cut,    # [1 1; 1 1]
      LET_FMLink, # [1 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 1]
      LET_AFLink, # [0 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
      LET_Vertex, # [0 0 0 0; 0 1 1 0; 0 1 1 0; 0 0 0 0]
      LET_Cross,  # [1 0 0 0; 0 0 1 0; 0 1 0 0; 0 0 0 1]
     )

@doc """
Loop element depicted as

```
|
o
|
```

or matrix

```
1 1 |+>
1 1 |->
```
""" LET_Cut

@doc """
Loop element depicted as

```
|  |
|--|
|  |
```

or matrix

```
1 0 0 0 |++>
0 0 0 0 |+->
0 0 0 0 |-+>
0 0 0 1 |-->
```
""" LET_FMLink

@doc """
Loop element depicted as

```
|  |
|~~|
|  |
```

or matrix

```
0 0 0 0 |++>
0 1 0 0 |+->
0 0 1 0 |-+>
0 0 0 0 |-->
```
""" LET_AFLink

@doc """
Loop element depicted as

```
|  |
~~~~
~~~~
|  |
```

or matrix

```
0 0 0 0 |++>
0 1 1 0 |+->
0 1 1 0 |-+>
0 0 0 0 |-->
```
""" LET_Vertex

@doc raw"""
Loop element depicted as

```
|  |
 \/
 /\
|  |
```

or matrix

```
1 0 0 0 |++>
0 0 1 0 |+->
0 1 0 0 |-+>
0 0 0 1 |-->
```
""" LET_Cross

@doc doc"""
(Imaginary-temporary and spatial) local operator as a perturbation with assigned loop element.
# Fields
- `let_type` : assigned loop element
- `isdiagonal` : operator is diagonal or not
    - in other words, two states connecting this perturbation are equivalent to
      each other or not.
- `time` : imaginary time ($\tau/\beta \in [0,1)$) which this perturbation acts on.
- `space` : spin or bond which this perturbation acts on.
            denotes `space` spin if space <= nspins
                or  `space - nspins` bond otherwise.
- `subspace` : subspin(s) index
- `bottom_id` :: index of node of union find assigned to a loop
- `top_id` :: index of node of union find assigned to the other loop

"""
mutable struct LocalLoopOperator
    let_type :: LoopElementType
    isdiagonal :: Bool
    time :: Float64
    space :: Int
    subspace :: Tuple{Int,Int}
    bottom_id :: Int
    top_id :: Int
end
function LocalLoopOperator(let_type::LoopElementType, time::Real, space::Int, subspace::Tuple{Int,Int})
    LocalLoopOperator(let_type, true, time, space, subspace, 0,0)
end


text:
importall Bukdu

type HTMLController <: ApplicationController
end

index(::HTMLController) = render(HTML, "
hello
")

Router() do
    get("/", HTMLController, index)
end

using Base.Test
conn = (Router)(get, "/")
@test 200 == conn.status
@test "
hello
" == conn.resp_body

logs = []
before(render, HTML) do t
    push!(logs, "b $t")
end
after(render, HTML) do t
    push!(logs, "a $t")
end

conn = (Router)(get, "/")
@test ["b
hello
", "a
hello
"] == logs

layout(::Layout, body) = """
$body
"""
show(::HTMLController) = render(HTML/Layout, "
hello
")
Router() do
    get("/say", HTMLController, show)
end

before(render, HTML/Layout) do t
    push!(logs, "bl $t")
end
after(render, HTML/Layout) do t
    push!(logs, "al $t")
end

empty!(logs)

conn = (Router)(get, "/say")
@test 200 == conn.status
@test "

hello

" == conn.resp_body

@test ["bl
hello
","b
hello
","a
hello
","al
hello
"] == logs


text:
for (frange, oneto) in ((mrange, OneToMRange),(srange ,OneToSRange))
    @testset "$oneto" begin
        let r = oneto(-5)
            @test isempty(r)
            @test length(r) == 0
            @test size(r) == (0,)
        end
        let r = oneto(3)
            @test !isempty(r)
            @test last(r) == r.stop
            @test_throws ErrorException r.notfield
            @test length(r) == 3
            @test size(r) == (3,)
            @test first(r) == 1
            @test last(r) == 3
            @test minimum(r) == 1
            @test maximum(r) == 3
            @test r[2] == 2
            @test r[2:3] == 2:3
            @test_throws BoundsError r[4]
            @test_throws BoundsError r[0]
            @test broadcast(+, r, 1) == 2:4
            @test 2*r == frange(2, step=2, stop=6)
            @test r + r == frange(2, step=2, stop=6)
            k = 0
            for i in r
                @test i == (k += 1)
            end
            @test intersect(r, oneto(2)) == oneto(2)
            @test intersect(r, OneTo(2)) == oneto(2)
            @test intersect(OneTo(2), r) == oneto(2)
            @test intersect(r, 0:5) == 1:3
            @test intersect(r, 2) == intersect(2, r) == frange(2, 2)
            @test findall(in(r), r) == findall(in(frange(1, length(r))), r) ==
                findall(in(r), frange(1, length(r))) == frange(1, length(r))
            io = IOBuffer()
            show(io, r)
            str = String(take!(io))
            @test str == "$(oneto)(3)"
            @test in(1, r) == true
        end

        @test oneto{Int}(oneto(10)) == oneto(10)
        @test oneto{UInt}(oneto(10)) == oneto(UInt(10))
        @test intersect(oneto(2), OneTo(2)) == oneto(2)
        @test intersect(OneTo(2), oneto(2)) == oneto(2)

        @test issubset(oneto(2), OneTo(2)) == true
        @test issubset(OneTo(2), oneto(2)) == true

        let r = oneto(7)
            @test findall(in(frange(2, (length(r) - 1))), r) == frange(2, (length(r) - 1))
            @test findall(in(r), frange(2, (length(r) - 1))) == frange(1, (length(r) - 2))
        end
        @test convert(oneto, 1:2) == oneto{Int}(2)
        @test_throws ArgumentError("first element must be 1, got 2") convert(oneto, 2:3)
        @test_throws ArgumentError("step must be 1, got 2") convert(oneto, 1:2:5)
        @test oneto(1:2) == oneto{Int}(2)
        @test oneto(1:1:2) == oneto{Int}(2)
        @test oneto{Int32}(1:2) == oneto{Int32}(2)
        @test oneto(Int32(1):Int32(2)) == oneto{Int32}(2)
        @test oneto{Int16}(3.0) == oneto{Int16}(3)
        @test_throws InexactError(:Int16, Int16, 3.2) oneto{Int16}(3.2)
    end
end


text:
abstract type AbstractScalarQuantity{T,N} <: AbstractArray{T,N} end

struct ScalarQuantity end

struct ScalarField{N,T,D<:AbstractArray{T,N},G<:AbstractAxisGrid} <: AbstractPICDataStructure{T,N,G}
    data::D
    grid::G
    name::String
end

struct ScalarVariable{T,D<:AbstractVector{T},G<:ParticlePositions} <: AbstractPICDataStructure{T,1,G}
    data::D
    grid::G
    name::String
end

ScalarField(data, grid; name="") = ScalarField(data, grid, name)
ScalarVariable(data, grid; name="") = ScalarVariable(data, grid, name)

@inline scalarness(::Type{<:ScalarField}) = ScalarQuantity()
@inline scalarness(::Type{<:ScalarVariable}) = ScalarQuantity()

# Indexing
@propagate_inbounds Base.getindex(::ScalarQuantity, f, i) = unwrapdata(f)[i]
@propagate_inbounds Base.setindex!(::ScalarQuantity, f, v, i) = unwrapdata(f)[i] = v

Base.eltype(::ScalarQuantity, f::Type{<:AbstractPICDataStructure{T}}) where T = T

# Broadcasting

Base.BroadcastStyle(::ScalarQuantity, ::Type{<:AbstractPICDataStructure}) = Broadcast.ArrayStyle{AbstractScalarQuantity}()

function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{AbstractScalarQuantity}}, ::Type{ElType}) where ElType
    # Scan the inputs for the AbstractPICDataStructure:
    f = find_field(bc)
    @debug "Building datastructure similar to $(typeof(f)) with eltype $ElType"
    grid = getdomain(f)
    # Keep the same grid for the output
    data = similar(unwrapdata(f), ElType, axes(bc))
    @debug "Output data type: $(typeof(data))"
    newstruct(f, data, grid)
end

function Base.similar(::ScalarQuantity, f, ::Type{S}, dims::Dims) where S
    # @debug "similar ScalarQuantity"
    grid = getdomain(f)
    data = similar(unwrapdata(f), S, dims)
    newstruct(f, data, grid)
end

Base.getproperty(::ScalarQuantity, f::AbstractPICDataStructure, key::Symbol) = getfield(f, key)
Base.propertynames(::ScalarQuantity, f::AbstractPICDataStructure) = fieldnames(typeof(f))


text:
@testset "Conversions" begin
    @test isequal(@data([1, 2, missing]),
                  convert(DataArray, @pdata([1, 2, missing])))

    # Test vector() and matrix() conversion tools
    dv = @data ones(5)
    @test isa(convert(Vector{Float64}, dv), Vector{Float64})
    dv[1] = missing
    # Should raise errors:
    # vector(dv)
    # convert(Vector{Float64}, dv)

    dm = @data ones(3, 3)
    @test isa(convert(Matrix{Float64}, dm), Matrix{Float64})
    dm[1, 1] = missing
    # Should raise errors:
    # matrix(dm)
    # convert(Matrix{Float64}, dm)

    a = DataArray(Any,2)
    convert(DataArray{Integer}, a)
    a[1] = 2
    convert(DataArray{Integer}, a)

    @test convert(DataArray, [1, missing]) isa DataVector{Int}
    @test isequal(convert(DataArray, [1, missing]), [1, missing])
    @test convert(DataArray{Int}, [1, missing]) isa DataVector{Int}
    @test isequal(convert(DataArray{Int}, [1, missing]), [1, missing])
    @test convert(DataArray{Any}, [1, missing]) isa DataVector{Any}
    @test isequal(convert(DataArray{Any}, [1, missing]), [1, missing])
    @test convert(DataArray{Int, 1}, [1, missing]) isa DataVector{Int}
    @test isequal(convert(DataArray{Int, 1}, [1, missing]), [1, missing])
    @test convert(DataArray{Any, 1}, [1, missing]) isa DataVector{Any}
    @test isequal(convert(DataArray{Any, 1}, [1, missing]), [1, missing])

    @test convert(Array, @data [1, missing]) isa Vector{Union{Int, Missing}}
    @test isequal(convert(Array, @data [1, missing]), [1, missing])
    @test_throws MethodError convert(Array{Int}, @data [1, missing])
    @test convert(Array{Union{Int, Missing}}, @data [1, missing]) isa Vector{Union{Int, Missing}}
    @test isequal(convert(Array{Union{Int, Missing}}, @data [1, missing]), [1, missing])
    @test convert(Array{Any}, @data [1, missing]) isa Vector{Any}
    @test isequal(convert(Array{Any}, @data [1, missing]), [1, missing])
    @test convert(Array{Union{Int, Missing}, 1}, @data [1, missing]) isa Vector{Union{Int, Missing}}
    @test isequal(convert(Array{Union{Int, Missing}}, @data [1, missing]), [1, missing])
    @test convert(Array{Any, 1}, @data [1, missing]) isa Vector{Any}
    @test isequal(convert(Array{Any, 1}, @data [1, missing]), [1, missing])
end


text:
module PDSQLiteTests

using Test
import ..AKCalc

@testset "Basic PlayerData entry" begin
  source = AKCalc.PDSQLite()
  base = AKCalc.OperatorBase("test_op_1", "Test Operator", 1)
  gd = AKCalc.GDMem()
  append!(gd, [base])
  op1 = AKCalc.Operator(base)
  @test append!(source, [op1]) == nothing
  @test AKCalc.operators(source, gd) == [op1]
end

end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    AbstractSpatialObject{T,N}

Spatial object in a `N`-dimensional space with coordinates of type `T`.
"""
abstract type AbstractSpatialObject{T,N} end

"""
    domain(object)

Return underlying domain of spatial `object`.
"""
domain(object) = object.domain

"""
    ndims(object)

Return the number of dimensions of domain underlying `object`.
"""
Base.ndims(::AbstractSpatialObject{T,N}) where {N,T} = N

"""
    coordtype(object)

Return the coordinate type of `object`.
"""
coordtype(::AbstractSpatialObject{T,N}) where {N,T} = T

"""
    coordnames(object)

Return the name of the coordinates in `object`.
"""
coordnames(object::AbstractSpatialObject{T,N}) where {N,T} = ntuple(i -> Symbol(:x,i), N)

"""
    npoints(object)

Return the number of points in `object`.
"""
npoints(object::AbstractSpatialObject) = npoints(domain(object))

"""
    coordinates(object, location)

Return the coordinates of the `location` in the `object`.
"""
function coordinates(object::AbstractSpatialObject{T,N},
                     location::Int) where {N,T}
  x = MVector{N,T}(undef)
  coordinates!(x, object, location)
  x
end

"""
    coordinates(object, locations)

Return the coordinates of `locations` in the `object`.
"""
function coordinates(object::AbstractSpatialObject{T,N},
                     locations::AbstractVector{Int}) where {N,T}
  X = Matrix{T}(undef, N, length(locations))
  coordinates!(X, object, locations)
  X
end

"""
    coordinates(object)

Return the coordinates of all locations in `object`.
"""
coordinates(object::AbstractSpatialObject) = coordinates(object, 1:npoints(object))

"""
    coordinates!(buff, object, locations)

Non-allocating version of [`coordinates`](@ref)
"""
function coordinates!(buff::AbstractMatrix, object::AbstractSpatialObject,
                      locations::AbstractVector{Int})
  for j in 1:length(locations)
    coordinates!(view(buff,:,j), object, locations[j])
  end
end

"""
    coordinates!(buff, object, location)

Non-allocating version of [`coordinates`](@ref).
"""
coordinates!(buff::AbstractVector, object::AbstractSpatialObject, location::Int) =
  coordinates!(buff, domain(object), location)