text:
Base.:+(A::Point2D{Cartesian},B::Point2D{Cartesian}) = Point2D{Cartesian}(A.x + B.x , A.z + B.z)

Base.:+(A::Point2D{Polar},B::Point2D{Polar}) = Point2D{Polar}(A.x + B.x , A.z + B.z)

Base.:*(A::Point2D{Polar},b::Number) = Point2D{Polar}(A.x * b , A.z * b)

Base.:*(A::Point2D{Cartesian},b::Number) = Point2D{Cartesian}(A.x * b , A.z * b)

function particlesadvection!(Particles, Δt; multiplier=1)
    Threads.@threads for i in axes(Particles,1)
        @inbounds Particles[i].CCart.x += Particles[i].UCart.x * Δt * multiplier
        @inbounds Particles[i].CCart.z += Particles[i].UCart.z * Δt * multiplier
    end
end

function particle_coordinates(P::Vector{PINFO}; coordinates_system = "cartesian")
    np = length(P)
    xp = Vector{Float64}(undef, np)
    zp = similar(xp)

    if coordinates_system == "cartesian"
        @batch for i in 1:np
            @inbounds xp[i] = P[i].CCart.x
            @inbounds zp[i] = P[i].CCart.z
        end

    elseif coordinates_system == "polar"
        @batch for i in 1:np
            @inbounds xp[i] = P[i].CPolar.x
            @inbounds zp[i] = P[i].CPolar.z
        end

    end

    return xp, zp
end

function particle_coordinates(P::PINFO)
    xp = P.CCart.x
    zp = P.CCart.z
    return xp, zp
end

function particle_coordinates!(xp, zp, P::Vector{PINFO})
    Threads.@threads for i in eachindex(P)
        @inbounds xp[i] = P[i].CCart.x
        @inbounds zp[i] = P[i].CCart.z
    end
end

#= INTERPOLATION FROM INTEGRATION POINT TO PARTICLE + utilities ======================
    EXTRAPOLATION FROM INTEGRATION POINTS TO NODAL POSITIONS
    Local coordinates of triangular element
          (0,1)                   η
            |\                    ^
            | \                   |
            |  \                  |
            |   \                ------> ξ
            |    \                |
            |     \
    (0,1/2) |      \ (1/2,1/2)
            |       \
            |        \
            |         \
            |          \
            |           \
            |____________\
        (0,0) (1/2,0) (1,0)
    New triangular coordinates for triangular elements in ip framework:
        xi'  = 2xi-1/3   --> xi  = xi'/2 + 1/6
        eta' = 2eta-1/3  --> eta = eta'/2 + 1/6
    Value at nodal coordinates
        w_i_vertex = sum(w_i_ip * N_i_node)
=========================================================================================#
function Fij2particle(particle_fields, particle_info, particle_weights, gr, F)

    EL2NOD_P1, EL2NOD = gr.e2n_p1, gr.e2n

    # (1) IP to NODE ----------------------------------------------------------
    # # -- Nodal element coordinates
    # ξ = @SVector [0.0, 1.0, 0.0]
    # η = @SVector [0.0, 0.0, 1.0]
    # # -- Transformation to ip triangular element coordinates
    # ξv = @. 2*ξ - 1/3
    # ηv = @. 2*η - 1/3

    # -- Transformation to ip triangular element coordinates
    X = @SMatrix [
        -1/3    -1/3
        2-1/3   -1/3
        -1/3    2-1/3
    ]

    # -- Shape functions
    nn3 = Val(3)
    NN, _ = shape_functions_triangles(X,nn3)

    # -- Map from integration points to the 6 element nodes
    m, n = size(F)
    Fxx = Array{Float64,2}(undef, m, n)
    Fzz, Fxz, Fzx = similar(Fxx), similar(Fxx), similar(Fxx)
    @inbounds Threads.@threads for i in CartesianIndices(F)
         Fxx[i] = F[i][1,1]
         Fxz[i] = F[i][1,2]
         Fzx[i] = F[i][2,1]
         Fzz[i] = F[i][2,2]
    end

    Fxx_nodal = nodalvalue(Fxx, NN)
    Fzz_nodal = nodalvalue(Fzz, NN)
    Fxz_nodal = nodalvalue(Fxz, NN)
    Fzx_nodal = nodalvalue(Fzx, NN)

    # -- Smooth fields
    nnod = maximum(EL2NOD)
    w = Vector{Float64}(undef,nnod) # weight buffer
    dummy = similar(w) # output buffer
    Fxx_smooth = smoothfield(Fxx_nodal,EL2NOD,w,dummy)
    Fzz_smooth = smoothfield(Fzz_nodal,EL2NOD,w,dummy)
    Fxz_smooth = smoothfield(Fxz_nodal,EL2NOD,w,dummy)
    Fzx_smooth = smoothfield(Fzx_nodal,EL2NOD,w,dummy)

    kernelF2particle!(
        particle_fields,
        Fxx_smooth,
        Fzz_smooth,
        Fxz_smooth,
        Fzx_smooth,
        EL2NOD_P1,
        particle_info,
        particle_weights
    )

    return particle_fields
end

@inline function barycentric_particle(v::NamedTuple, pc::Point2D{T}) where {T}
    x1,x2,x3 = v.x[1], v.x[2], v.x[3]
    z1,z2,z3 = v.z[1], v.z[2], v.z[3]
    x,z = pc.x, pc.z
    b1 = @muladd ((z2 - z3) * (x - x3) + (x3 - x2) * (z - z3)) /
         ((z2 - z3) * (x1 - x3) + (x3 - x2) * (z1 - z3))
    b2 = @muladd ((z3 - z1) * (x - x3) + (x1 - x3) * (z - z3)) /
         ((z2 - z3) * (x1 - x3) + (x3 - x2) * (z1 - z3))
    b3 = 1 - b1 - b2
    @SVector [b1, b2, b3]
end

@inline dot_F2part(F, t, N, i1, i2, idx) =
    @muladd F[t, idx[1]][i1, i2]*N[1] + F[t, idx[2]][i1, i2]*N[2] + F[t, idx[3]][i1, i2]*N[3]

function kernel_Fij_particle!(particle_fields, particle_info, ipx, ipz, F, i)
    t = particle_info[i].t_parent
    idx = @SVector [4, 6, 5]
    v = (
        x=(view(ipx, t, idx)),
        z=(view(ipz, t, idx))
    ) # Corner ips := [4 6 5]

    N = barycentric_particle(v, particle_info[i].CPolar)

    particle_fields.Fxx[i] = dot_F2part(F, t, N, 1, 1, idx)
    particle_fields.Fzx[i] = dot_F2part(F, t, N, 2, 1, idx)
    particle_fields.Fxz[i] = dot_F2part(F, t, N, 1, 2, idx)
    particle_fields.Fzz[i] = dot_F2part(F, t, N, 2, 2, idx)

end

function F2particle(particle_fields, particle_info, ipx, ipz, F)
    @batch for ipart in 1:length(particle_info)
        kernel_Fij_particle!(particle_fields, particle_info, ipx, ipz, F, ipart)
    end
    particle_fields
end

function ip2node(EL2NOD, F)
    # -- Nodal element coordinates
    ξ = Float64.([0,1,0])
    η = Float64.([0,0,1])

    # -- Transformation to ip triangular element coordinates
    ξv = @. 2.0*ξ - 1.0/3.0
    ηv = @. 2.0*η - 1.0/3.0

    # -- Shape functions
    nn3 = Val(3)
    NN,_ = shape_functions_triangles([ξv ηv], nn3)
    # NN = [SVector{3}(NN[i]) for i=1:3]

    # -- Map from integration points to the 6 element nodes
    F_nodal = nodalvalue(F, NN)

    # -- Smooth fields
    nnod = maximum(EL2NOD)
    w = Vector{Float64}(undef,nnod) # weight buffer
    dummy = similar(w) # output buffer
    smoothfield(F_nodal, EL2NOD, w, dummy)
    # smoothfield(F_nodal, view(EL2NOD,1:6,:), area_el, w, dummy)
end

function Fijkernel!(FF, idx1, idx2, NN, EL2NOD, EL2NOD_P1, w, dummy,
                   particle_info, particle_weights, particle_fields)

    m, n = size(FF)
    F = Array{Float64,2}(undef, m, n)
    Threads.@threads for i in CartesianIndices(FF)
        @inbounds F[i] = FF[i][idx1, idx2]
    end

    # F = [FF[i][idx1,idx2] for i in CartesianIndices(FF)]
    F_nodal = nodalvalue(F, NN)
    F_smooth = smoothfield(F_nodal, view(EL2NOD,1:6,:), w, dummy)
    # smoothfield(F_nodal, EL2NOD, area_el, w, dummy)

    Threads.@threads for i in axes(particle_info,1)
        @inbounds particle_fields.Fxx[i] = node2particle(
            view(F_smooth, view(EL2NOD_P1,:,particle_info[i].t)),
            particle_weights[i].barycentric
        )
    end
end

function Fij2particle_thread!(particle_fields,particle_info,particle_weights,EL2NOD_P1,EL2NOD,
    area_el,F)

    # (1) IP to NODE ----------------------------------------------------------
    # -- Nodal element coordinates
    ξ = Float64.([0,1,0])
    η = Float64.([0,0,1])

    # -- Transformation to ip triangular element coordinates
    ξv = @. 2.0*ξ - 1.0/3.0
    ηv = @. 2.0*η - 1.0/3.0

    # -- Shape functions# -- Shape functions
    nn3 = Val(3)
    NN,_ = shape_functions_triangles([ξv ηv], nn3)

    # -- Map from integration points to the 6 element nodes
    w  = Vector{Float64}(undef,nnod) # weight buffer
    dummy = similar(w) # output buffer

    Fijkernel!(F, 1, 1, NN, EL2NOD, EL2NOD_P1, w, dummy, particle_info, particle_weights, particle_fields)
    Fijkernel!(F, 1, 2, NN, EL2NOD, EL2NOD_P1, w, dummy, particle_info, particle_weights, particle_fields)
    Fijkernel!(F, 2, 1, NN, EL2NOD, EL2NOD_P1, w, dummy, particle_info, particle_weights, particle_fields)
    Fijkernel!(F, 2, 2, NN, EL2NOD, EL2NOD_P1, w, dummy, particle_info, particle_weights, particle_fields)

end

"""
Interpolate smoothed nodal Fij to particles
"""
function kernelF2particle!(particle_fields,Fxx_smooth,Fzz_smooth,Fxz_smooth,Fzx_smooth, EL2NOD_P1,particle_info,particle_weights)
    Threads.@threads for i in axes(particle_info,1)
        # idx = view(EL2NOD_P1,:,particle_info[i].t)
        @inbounds particle_fields.Fxx[i] =
            node2particle(view(Fxx_smooth, view(EL2NOD_P1,:,particle_info[i].t)), particle_weights[i].barycentric)
        @inbounds particle_fields.Fzz[i] =
            node2particle(view(Fzz_smooth, view(EL2NOD_P1,:,particle_info[i].t)), particle_weights[i].barycentric)
        @inbounds particle_fields.Fxz[i] =
            node2particle(view(Fxz_smooth, view(EL2NOD_P1,:,particle_info[i].t)), particle_weights[i].barycentric)
        @inbounds particle_fields.Fzx[i] =
            node2particle(view(Fzx_smooth, view(EL2NOD_P1,:,particle_info[i].t)), particle_weights[i].barycentric)
    end
end

function smoothfield(F, e2n, w, dummy)
    fill!(dummy, 0.0)
    fill!(w, 0)
    for i in 1:6
        idx = view(e2n,i,:)
        @views w[idx] .+= 1
        @views dummy[idx] .+= view(F, :, i)
    end
    dummy./w
end

"""
Extrapolate field from integration point to nodal positions
    (*) A field @ integration points
    (*) NN linear shape functions
"""
@inline function nodalvalue(A, NN)
    A3 = takecornerips(A)
    Anodal = Array{Float64,2}(undef, size(A3,1), size(A,2))
    for u in CartesianIndices(A3)
        @inbounds Anodal[u] = mydot(NN[u.I[2]], view(A3, u.I[1],:))
    end
    averageedge!(Anodal)
    return Anodal
end

"""
Take values at the linear integration points
"""
@inline function takecornerips(A)
    m,n = size(A,1), Int(size(A,2)/2)
    B   = Array{Float64,2}(undef,m,n)
    for i in 1:m
        @inbounds B[i,1] = A[i,4]
        @inbounds B[i,2] = A[i,6]
        @inbounds B[i,3] = A[i,5]
    end
    B
end

takecornerips(A::Array{SArray{Tuple{2,2},Float64,2,4},2}) = @views A[:, [4, 6, 5]]

"""
Average field at nodes located on the edges of the element
"""
@inline function averageedge!(A)
    for i in axes(A,1)
        @inbounds A[i,4] = 0.5 * (A[i,1] + A[i,2])
        @inbounds A[i,5] = 0.5 * (A[i,3] + A[i,2])
        @inbounds A[i,6] = 0.5 * (A[i,1] + A[i,3])
    end
end

"""
node2particle(A::Vector{T}, barycentric::Vector{T})

    Linear interpolation from element nodes to arbitrary point within the element using barycentric coordinates (equivalent linear  FEM shape functions.)

        A := field defined at triangular vertices
        N := barycentric coordinates of i-th particle
"""
@inline node2particle(A_p1,N) = mydot(A_p1,N)

function velocities2ip(Ucartesian, EL2NOD)
    # =========== PREPARE INTEGRATION POINTS & DERIVATIVES wrt LOCAL COORDINATES
    nip = 6
    nnodel = 3
    nel = size(EL2NOD,2)
    ni, nn = Val(nip), Val(nnodel)
    N, dNds, _, w_ip  = _get_SF(ni,nn)
    Ux_ip = Array{Float64,2}(undef,nel,nip)
    Uz_ip = similar(Ux_ip)

    for iel in 1:nel
        Ux_el, Uz_el = getvelocity(Ucartesian, view(EL2NOD,:,iel))
        @inbounds for ip in 1:nip
            Ux_ip[iel, ip] = mydot(Ux_el[1:3], N[ip])
            Uz_ip[iel, ip] = mydot(Uz_el[1:3], N[ip])
        end

    end
    return Ux_ip, Uz_ip
end

function getvelocity(U,idx)
    Ux = @SVector [U[idx[1]].x, U[idx[2]].x, U[idx[3]].x, U[idx[4]].x, U[idx[5]].x, U[idx[6]].x]
    Uz = @SVector [U[idx[1]].z, U[idx[2]].z, U[idx[3]].z, U[idx[4]].z, U[idx[5]].z, U[idx[6]].z]
    return Ux, Uz
end

getvelocityX(U,idx) = @SVector [U[idx[1]].x, U[idx[2]].x, U[idx[3]].x]
getvelocityZ(U,idx) = @SVector [U[idx[1]].z, U[idx[2]].z, U[idx[3]].z]

function getvelocity(U)
    Ux = @inbounds [u.x for u in U]
    Uz = @inbounds [u.z for u in U]
    return Ux,Uz
end

# INTERPOLATE TEMPERATURE FROM NODES TO PARTICLES =======================================
@inline function initial_particle_temperature!(particle_fields, e2n_p1, T, particle_info, particle_weights)
    np = length(particle_info)
    Threads.@threads for i in 1:np
        # -- Temperature interpolation (weighted average)
        @inbounds particle_fields.T[i] = weightednode2particle(
            view(T, view(e2n_p1,:,particle_info[i].t)),
            particle_weights[i].barycentric
            )
    end

end

function interpolate_temperature!(
                        T0,
                        particle_fields,
                        gr,
                        ρ,
                        T,
                        particle_info,
                        particle_weights,
                        VarT,
                        nθ,
                        nr,
                        Δt,
                        ΔT,
                    )

    particle_fields = temperature2particle(
            T0,
            particle_fields,
            gr.e2n_p1,
            ρ,
            T,
            particle_info,
            particle_weights,
            VarT,
            nθ,
            nr,
            Δt,
        )

    ΔT_remaining = ΔTsubgrid2node(
            ΔT,
            particle_fields,
            particle_info,
            particle_weights,
            gr
        )

    corrected_particle_temperature!(
            particle_fields, ΔT_remaining
        )
end

# INTERPOLATE TEMPERATURE FROM NODES TO PARTICLES =======================================
@inline function temperature2particle(T0, particle_fields, e2n_p1, ρ, T, particle_info, particle_weights, VarT, nθ, nr, Δt)

    Cp, κ = VarT.Cp, VarT.κ

    # allocations
    T, ΔT_subgrid = particle_fields.T, particle_fields.ΔT_subgrid
    t = [p.t for p in particle_info]
    # nt = Threads.nthreads()
    # Tᵢ0 = Vector{Float64}(undef, nt)
    # Tᵢ = Vector{Float64}(undef, nt)
    # ρᵢ = Vector{Float64}(undef, nt)

    # constant variables
    Δθ = 2π*1.22/(nθ)
    Δr = 1/(nr-1)
    #=
        Gerya's book, pp. 143
        Δt_diff = Cp * ρₘ / (κₘ * (2/Δx² + Δz²) )
        0 ≤ d ≤ 1 , numerical diffusion coefficient
    =#
    multiplier = Δt/(Cp/(κ*(2/((Δθ/2)^2) + 2/((Δr/2)^2)))) # == Δt/Δt_diff
    d = 1

    # loop
    Threads.@threads for i in eachindex(particle_info)
        @inbounds T[i], ΔT_subgrid[i] =
            _temperature2particle(
                T[i],
                view(T0, view(e2n_p1,:,t)),
                particle_weights[i].barycentric,
                d,
                multiplier
            )
    end

    particle_fields.T.= T
    particle_fields.ΔT_subgrid .= ΔT_subgrid

    return particle_fields

end

function _temperature2particle(Ti, T0, barycentric, d, multiplier)
    # Old T
    @inbounds Tᵢ0 = weightednode2particle(
       T0,
       barycentric
    )

    ΔT_subgrid = (Tᵢ0 - Ti)*(1 - exp(-d*multiplier))
    Ti += ΔT_subgrid

    return Ti, ΔT_subgrid
end

function corrected_particle_temperature!(particle_fields, ΔT_remaining)
    Threads.@threads for i in eachindex(ΔT_remaining)
        @inbounds particle_fields.T[i] += ΔT_remaining[i]
    end
end

function F2ip(F, particle_fields, particle_info, particle_weights, nel)
    np = length(particle_info)
    weight = Array{Float64,2}(undef, np, 6)
    Threads.@threads for i in axes(particle_info,1)
        ## THIS IS TYPE UNSTABLE -> FIX ASAP
        @inbounds @views weight[i,:] .= 1.0./(particle_weights[i].ip_weights.^2) # inverse of the distance
    end

    F2ip_inner_tkernel(F, particle_info, particle_fields, weight, nel)

    # # weight times field
    # t = Vector{Int32}(undef, np)
    # wxx, wzz, wxz, wzx =  similar(weight), similar(weight), similar(weight), similar(weight)
    # Threads.@threads for i in axes(weight,1)
    #     for j in axes(weight,2)
    #         @inbounds w = weight[i,j]
    #         @inbounds wxx[i,j] = w*particle_fields.Fxx[i]
    #         @inbounds wzz[i,j] = w*particle_fields.Fzz[i]
    #         @inbounds wxz[i,j] = w*particle_fields.Fxz[i]
    #         @inbounds wzx[i,j] = w*particle_fields.Fzx[i]
    #     end
    #    @inbounds t[i] = Int32(particle_info[i].t_parent) # t-th element hosting particle i-ith
    # end

    # # Interpolate from particle to integration point
    # for j in 1:6
    #     accw = 1.0./accumarray(t,view(weight,:,j)) # weights sum
    #     Fxx  = accumarray(t,view(wxx,:,j)).*accw
    #     Fzz  = accumarray(t,view(wzz,:,j)).*accw
    #     Fxz  = accumarray(t,view(wxz,:,j)).*accw
    #     Fzx  = accumarray(t,view(wzx,:,j)).*accw
    #     Threads.@threads for i in 1:nel
    #         @inbounds F[i,j] = @SMatrix [Fxx[i] Fxz[i]; Fzx[i] Fzz[i]]
    #     end
    # end

    return F
end

function F2ip_inner_kernel(F, particle_info, particle_fields, weight, nel)
    Fxx = fill(0.0, nel, 6)
    Fzz = fill(0.0, nel, 6)
    Fzx = fill(0.0, nel, 6)
    Fxz = fill(0.0, nel, 6)
    for i in axes(weight,1)
        ω = ntuple(j->weight[i,j]::Float64, 6)
        @inbounds @views @. Fxx[particle_info[i].t_parent, :] += particle_fields.Fxx[i]* ω
        @inbounds @views @. Fxz[particle_info[i].t_parent, :] += particle_fields.Fxz[i]* ω
        @inbounds @views @. Fzx[particle_info[i].t_parent, :] += particle_fields.Fzx[i]* ω
        @inbounds @views @. Fzz[particle_info[i].t_parent, :] += particle_fields.Fzz[i]* ω
    end
    Threads.@threads for i in eachindex(Fxx)
        @inbounds F[i] = @SMatrix [Fxx[i] Fxz[i]; Fzx[i] Fzz[i]]
    end
    return F
end

function F2ip_inner_tkernel(F, particle_info, particle_fields, weight, nel)
    nt = Threads.nthreads()
    Fxx = [fill(0.0, nel, 6) for _ in 1:nt]
    Fzz = [fill(0.0, nel, 6) for _ in 1:nt]
    Fzx = [fill(0.0, nel, 6) for _ in 1:nt]
    Fxz = [fill(0.0, nel, 6) for _ in 1:nt]
    accw = [fill(0.0, nel, 6) for _ in 1:nt]

    # Acummarray on each thread
    Threads.@threads for i in axes(weight,1)
        # @inbounds ω = ntuple(j->weight[i,j], 6)
        @inbounds @views Fxx[Threads.threadid()][particle_info[i].t_parent, :]  .+= particle_fields.Fxx[i].* ntuple(j->weight[i,j], 6)
        @inbounds @views Fxz[Threads.threadid()][particle_info[i].t_parent, :]  .+= particle_fields.Fxz[i].* ntuple(j->weight[i,j], 6)
        @inbounds @views Fzx[Threads.threadid()][particle_info[i].t_parent, :]  .+= particle_fields.Fzx[i].* ntuple(j->weight[i,j], 6)
        @inbounds @views Fzz[Threads.threadid()][particle_info[i].t_parent, :]  .+= particle_fields.Fzz[i].* ntuple(j->weight[i,j], 6)
        @inbounds @views accw[Threads.threadid()][particle_info[i].t_parent, :] .+= ntuple(j->weight[i,j], 6)
    end

    # Acummulate on first array
    for i in 2:nt
        @inbounds Fxz[1] .+= Fxz[i]
        @inbounds Fxx[1] .+= Fxx[i]
        @inbounds Fzx[1] .+= Fzx[i]
        @inbounds Fzz[1] .+= Fzz[i]
        @inbounds accw[1] .+= accw[i]
    end

    # Store F tensor
    Threads.@threads for i in CartesianIndices(Fxx[1])
        @inbounds F[i] = @SMatrix [Fxx[1][i]/accw[1][i] Fxz[1][i]/accw[1][i]; Fzx[1][i]/accw[1][i] Fzz[1][i]/accw[1][i]]
    end
    return F
end

function F2ip_inner_tkernel_unsafe(F, particle_info, particle_fields, weight, nel)
    Fxx = fill(0.0, nel, 6)
    Fzz = fill(0.0, nel, 6)
    Fzx = fill(0.0, nel, 6)
    Fxz = fill(0.0, nel, 6)
    accw = fill(0.0, nel, 6)
    # Acummarray on each thread
    Threads.@threads for i in axes(weight,1)
        ω = ntuple(j->weight[i,j], 6)
        @inbounds @views @. Fxx[particle_info[i].t_parent, :] += particle_fields.Fxx[i]* ω
        @inbounds @views @. Fxz[particle_info[i].t_parent, :] += particle_fields.Fxz[i]* ω
        @inbounds @views @. Fzx[particle_info[i].t_parent, :] += particle_fields.Fzx[i]* ω
        @inbounds @views @. Fzz[particle_info[i].t_parent, :] += particle_fields.Fzz[i]* ω
        @inbounds @views @. accw[particle_info[i].t_parent, :] += ω
    end

    # Store F tensor
    Threads.@threads for i in CartesianIndices(Fxx)
        @inbounds F[i] = @SMatrix [Fxx[i]/accw[i] Fxz[i]/accw[i]; Fzx[i]/accw[i] Fzz[i]/accw[i]]
    end
    return F
end

#=
    Interpolate from particles to nodes:
        T_node = ∑T_particle_i*w_i / ∑w_i
    where
        w = 1/distance(node -> ip)^n
=#
function T2node(T, particle_fields, particle_info, particle_weights, gr, IDs)

    EL2NOD_P1 = gr.e2n_p1

    #=  Build arrays  =#
    t = [particle_info[i].t for i in axes(particle_info,1)]
    w = Array{Float64,2}(undef, size(particle_weights,1),size(particle_weights[1].node_weights,1))
    for i in axes(particle_weights,1)
        dummy = SVector{3}(particle_weights[i].barycentric)
        for j in 1:3
            # @inbounds w[i,j] = 1/(dummy[j]^2)
            @inbounds w[i,j] = dummy[j]
        end
    end

    # T2node_inner_tkernel(T, particle_info, particle_fields, w, EL2NOD_P1)
    # applybounds!(T, 1.0, 0.0)
    # T

    nodes = transpose(view(EL2NOD_P1,:,t))
    A = w.*particle_fields.T
    nn = Array(nodes)
    A1 = kernel2(nn, A, true)
    A2 = kernel2(nn, w, true)
    T .= vec(A1./A2)
    fixT!(T, 0.0, 1.0, IDs)
    # T[IDs .== "outter"] .= 0.0
    applybounds!(T, 1.0, 0.0)
    T
end

function T2node_inner_tkernel_unsafe(T, particle_info, particle_fields, weight, EL2NOD_P1)
    fill!(T, 0.0)
    accw = fill(0.0, length(T))
    # Acummarray on each thread
    Threads.@threads for i in axes(weight,1)
        ω = ntuple(j->weight[i,j], 3)
        @inbounds @views T[view(EL2NOD_P1,:,particle_info[i].t)] .+= particle_fields.T[i].*ω
        @inbounds @views accw[view(EL2NOD_P1,:,particle_info[i].t)] .+= ω
    end

    # Store T tensor
    Threads.@threads for i in eachindex(T)
        @inbounds T[i] /= accw[i]
    end
    return T
end

function T2node_inner_tkernel(T, particle_info, particle_fields, weight, EL2NOD_P1)
    fill!(T, 0.0)
    nt = Threads.nthreads()
    Tnt = [fill(0.0, length(T)) for _ in 1:nt-1]
    accw = [fill(0.0, length(T)) for _ in 1:nt]
    # Acummarray on each thread
    Threads.@threads for i in axes(weight,1)
        ω = ntuple(j->weight[i,j], 3)
        if Threads.threadid() == 1
            @inbounds @views T[view(EL2NOD_P1,:,particle_info[i].t)] .+= particle_fields.T[i].*ω
        else
            @inbounds @views Tnt[Threads.threadid()-1][view(EL2NOD_P1,:,particle_info[i].t)] .+= particle_fields.T[i].*ω
        end
        @inbounds @views accw[Threads.threadid()][view(EL2NOD_P1,:,particle_info[i].t)] .+= ω
    end

    # Acummulate on first array
    for i in 2:nt
        @inbounds T .+= Tnt[i-1]
        @inbounds accw[1] .+= accw[i]
    end

    # Store T tensor
    Threads.@threads for i in eachindex(T)
        @inbounds T[i] /= accw[1][i]
    end
    return T
end

function ΔTsubgrid2node(ΔT, particle_fields, particle_info, particle_weights, gr)

    EL2NOD_P1 = gr.e2n_p1

    #=  Build arrays  =#
    t = [p.t for p in particle_info]
    w = Matrix{Float64}(undef, size(particle_weights,1), size(particle_weights[1].node_weights,1))
    # t = Vector{Int64}(undef, size(particle_weights,1))
    A = similar(w)
    Threads.@threads for i in axes(particle_weights,1)
        for j in 1:3
            @inbounds ω = particle_weights[i].barycentric[j]
            @inbounds w[i,j] = ω
            @inbounds A[i,j] = ω*particle_fields.ΔT_subgrid[i]
        end
    end

    # nodes = transpose(view(EL2NOD_P1,:,t))
    # nn = Array(nodes)
    nn = view(EL2NOD_P1, :, t)'
    A1 = kernel2(nn, A, true) # accumarray -> ∑Tᵢ*ωᵢ
    A2 = kernel2(nn, w, true) # accumarray -> ∑ωᵢ

    ΔT_remaining = similar(A1)
    @turbo for i in axes(A1,1)
        ΔT_remaining[i] = ΔT[i] - A1[i]/A2[i]
    end

    nodalfield2partcle(ΔT_remaining, EL2NOD_P1, particle_weights, particle_info) # interpolate to particles

    # # ΔT_subgrid in element vertices node
    # ΔT_subgrid_vert = vec(A1./A2) # -> ∑Tᵢ*ωᵢ/∑ωᵢ
    # # Add ΔT_subgrid in center node
    # ΔT_subgrid_el = view(ΔT_subgrid_vert, view(EL2NOD,1:6,:)) # element T
    # ΔT_subgrid_bubble = vec(mean(ΔT_subgrid_el, dims = 1))
    # ΔT_subgrid = vcat(ΔT_subgrid_vert, ΔT_subgrid_bubble)

end

"""
nodalfield2partcle(field, e2n, particle_weights, particle_info)

    Interpolate field from node to corresponding particle
"""
function nodalfield2partcle(field, e2n, particle_weights, particle_info)
    # memory allocations
    np = length(particle_info)
    field_particle = Vector{Float64}(undef, np)
    t = [p.t for p in particle_info]
    barycentric = [p.barycentric for p in particle_weights]

    # interpolate from node to corresponding particle
    Threads.@threads for i in 1:np
        @inbounds field_particle[i] = weightednode2particle(
            (@views field[e2n[:, t[i]]]),
            barycentric[i]
        )
    end
    field_particle
end

"""
    applybounds!(A::AbstractArray{T}, upper::T, lower::T) where {T<:Real}
"""
function applybounds!(A::AbstractArray{T}, upper::T, lower::T) where {T<:Real}
    Threads.@threads for i in eachindex(A)
        @inbounds A[i] = ifelse(A[i] > upper, upper, A[i])
        @inbounds A[i] = ifelse(A[i] < lower, lower, A[i])
    end
end

"""
    applybounds(A::T, upper::T,lower::T) where {T<:Real}
"""
applybounds(A::T, upper::T, lower::T) where {T<:Real} = max(min(A, upper), lower)

ip2particle(ip_field,iw) = mydot(ip_field, iw) / sum(iw)

function kernel1(a, b, collapse=false)
    m,n = maximum(a),size(a,2)
    v = fill(0.0,m,n)
    for j in axes(v,2), _ in axes(b,1)
        @inbounds v[a[i,j],j] += b[i][j]
    end
    collapse == true ? sum(v,dims=2) : v
end

function kernel2(a, b, collapse=false)
    m, n = maximum(a), size(a,2)
    v = fill(0.0,m,n)
    Threads.@threads for j in axes(v,2)
        for i in axes(b,1)
            @inbounds v[a[i,j],j] += b[i,j]
        end
    end
    collapse == true ? sum(v,dims=2) : v
end

@inline function T2particle(particle_fields, e2n_p1, T, particle_info, particle_weights)
    @inbounds Threads.@threads for i in axes(particle_info,1)
        # -- Temperature interpolation (weighted average)
        particle_fields.T[i] = weightednode2particle(view(T, view(e2n_p1,:,particle_info[i].t)),
                                                     fw.(particle_weights[i].barycentric))
    end
   return particle_fields
end

@inline fw(x) = x^3 * exp(-3*x)

@inline function weightednode2particle(A, ω)
    a, b = 0.0, 0.0
    @turbo warn_check_args=false for i in axes(A,1)
        a += ω[i]*A[i]
        b += ω[i]
    end
    a/b
end

function T2particle_cubic!(particle_fields, gr, T, particle_info)
    Tp = quasicubic_interpolation(gr, particle_info, T)
    Threads.@threads for i in axes(Tp,1)
        @inbounds particle_fields.T[i] = max(min(Tp[i], 1.0), 0.0)
    end
end

function velocities2particle_cubic!(gr, particle_info, U)
    Ux0, Uz0 = getvelocity(U)

    Uxp = quasicubic_interpolation(gr, particle_info, Ux0)
    Uzp = quasicubic_interpolation(gr, particle_info, Uz0)

    Threads.@threads for i in axes(particle_info,1)
        @inbounds particle_info[i].UCart.x = Uxp[i]
        @inbounds particle_info[i].UCart.z = Uzp[i]
    end

end

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
"""
    euler(f,a,b,y0,M)
berisi program untuk mencari solusi persamaan differensial `y'=f(t,y)`
dengan masalah nilai awal `y(a) = y0` pada interval `[a, b]`. Program ini secara default
berisi 5 masukan, yaitu fungsi `f(t,y)`, titik ujung interval penyelesaian `[a,b]`, nilai awal
`y0`, dan jumlah sub-interval `M`.
# Example
```jl
julia> f(t,y) = (t-y)/2;

julia> sol = euler(f,0,3,1,6)
7×2 Array{Float64,2}:
 0.0  1.0
 0.5  0.75
 1.0  0.6875
 1.5  0.765625
 2.0  0.949219
 2.5  1.21191
 3.0  1.53394
```
return solusi masalah nilai awal `sol`.
"""
function euler(f,a,b,y0,M)
    M = Int(M)
    h = (b-a)/M;
    T = a:h:b;
    Y = Array{Float64}(undef,length(T),1)
    Y[1] = y0;
    # Mulai langkah iterasi euler
    for k = 1:M
        # Rumus iterasi euler
        Y[k+1] = Y[k]+h*f(T[k],Y[k])
    end
    sol = [T Y]
end

ext: jl


text:
using RobustPCA
using Base.Test

D = [0.462911    0.365901  0.00204357    0.692873    0.935861;
    0.0446199    0.108606   0.0664309   0.0736707    0.264429;
     0.320581    0.287788    0.073133    0.188872    0.526404;
     0.356266    0.197536 0.000718338    0.513795    0.370094;
     0.677814    0.011651    0.818047   0.0457694    0.471477]

A = [0.462911   0.365901  0.00204356   0.345428   0.623104;
     0.0446199  0.108606  0.0429271    0.0736707  0.183814;
     0.320581   0.203777  0.073133     0.188872   0.472217;
     0.30725    0.197536  0.000717701  0.201626   0.370094;
     0.234245   0.011651  0.103622     0.0457694  0.279032]

E = [0.0        0.0        0.0        0.347445  0.312757 ;
     0.0        0.0        0.0235038  0.0       0.0806151;
     0.0        0.0840109  0.0        0.0       0.0541868;
     0.0490157  0.0        6.5061e-7  0.312169  0.0      ;
     0.443569   0.0        0.714425   0.0       0.192445]

@time Â, Ê = inexact_alm_rpca(D, nonnegativeE=true, nonnegativeA=true)

@test_approx_eq_eps Â A 1.0e-6
@test_approx_eq_eps Ê E 1.0e-6
# reconstruction eror
@test norm(D - (Â + Ê))/norm(D) < 1.0e-7


text:
"""
sample!([rng, ]a, k)

Sample `k` element from array `a` without repetition and eventually excluding elements in `exclude`.

### Optional Arguments
- `exclude=()`: elements in `a` to exclude from sampling.

### Implementation Notes
Changes the order of the elements in `a`. For a non-mutating version, see [`sample`](@ref).
"""
function sample!(rng::AbstractRNG, a::AbstractArray, k::Integer; exclude = ())
    length(a) < k + length(exclude) && error("Array too short.")
    res = Vector{eltype(a)}()
    sizehint!(res, k)
    n = length(a)
    i = 1
    while length(res) < k
        r = rand(rng, 1:(n - i + 1))
        if !(a[r] in exclude)
            push!(res, a[r])
            a[r], a[n - i + 1] = a[n - i + 1], a[r]
            i += 1
        end
    end
    res
end

sample!(a::AbstractArray, k::Integer; exclude = ()) = sample!(getRNG(), a, k; exclude = exclude)

"""
    sample([rng,] r, k)

Sample `k` element from unit range `r` without repetition and eventually excluding elements in `exclude`.

### Optional Arguments
- `exclude=()`: elements in `a` to exclude from sampling.

### Implementation Notes
Unlike [`sample!`](@ref), does not produce side effects.
"""
sample(a::UnitRange, k::Integer; exclude = ()) = sample!(getRNG(), collect(a), k; exclude = exclude)

getRNG(seed::Integer = -1) = seed >= 0 ? MersenneTwister(seed) : Base.Random.GLOBAL_RNG


text:
module TestCaches

using Test
using RigidBodyDynamics
using RigidBodyDynamics.Contact
import Random

function randmech()
    rand_tree_mechanism(Float64,
        QuaternionFloating{Float64},
        [Revolute{Float64} for i = 1 : 5]...,
        [Fixed{Float64} for i = 1 : 5]...,
        [Prismatic{Float64} for i = 1 : 5]...,
        [Planar{Float64} for i = 1 : 5]...,
        [SPQuatFloating{Float64} for i = 1:2]...,
        [SinCosRevolute{Float64} for i = 1:2]...
    )
end

function cachetest(cache, eltypefun)
    x64 = @inferred cache[Float64]
    @test eltypefun(x64) == Float64
    @test cache[Float64] === x64
    @test @allocated(cache[Float64]) == 0

    x32 = @inferred cache[Float32]
    @test eltypefun(x32) == Float32
    @test cache[Float32] === x32
    @test @allocated(cache[Float32]) == 0

    @test cache[Float64] === x64
    @test @allocated(cache[Float64]) == 0
end

@testset "StateCache" begin
    Random.seed!(1)
    mechanism = randmech()
    cache = StateCache(mechanism)
    cachetest(cache, RigidBodyDynamics.state_vector_eltype)
end

@testset "DynamicsResultCache" begin
    @testset "Basic mechanism" begin
        Random.seed!(2)
        mechanism = randmech()
        cache = DynamicsResultCache(mechanism)
        cachetest(cache, result -> eltype(result.v̇))
    end

    @testset "Mechanism with contact points (Issue #483)" begin
        Random.seed!(3)
        mechanism = randmech()
        contactmodel = SoftContactModel(hunt_crossley_hertz(k = 500e3),
            ViscoelasticCoulombModel(0.8, 20e3, 100.))
        body = rand(bodies(mechanism))
        add_contact_point!(body,
            ContactPoint(Point3D(default_frame(body), 0.0, 0.0, 0.0), contactmodel))
        cache = DynamicsResultCache(mechanism)
        cachetest(cache, result -> eltype(result.v̇))
    end
end

@testset "SegmentedVectorCache" begin
    Random.seed!(3)
    mechanism = randmech()
    state = MechanismState(mechanism)
    cache = SegmentedVectorCache(RigidBodyDynamics.ranges(velocity(state)))
    cachetest(cache, eltype)
end

end # module


text:
using Test
using ClimateMachine
using ClimateMachine.Mesh.Grids
using ClimateMachine.Mesh.Topologies
using ClimateMachine.MPIStateArrays
using ImperoHooks
import ImperoHooks: launch_volume_gradient!, launch_interface_gradient!
using Impero, Printf, MPI, LinearAlgebra, Statistics, GaussQuadrature
include(pwd() * "/test/gradient_test_utils.jl")

ClimateMachine.init()
const ArrayType = ClimateMachine.array_type()
const mpicomm = MPI.COMM_WORLD
const FT = Float64
Ω = Circle(-1,1) × Circle(-1,1) × Circle(-1,1)
dims = ndims(Ω)

ClimateMachine.gpu_allowscalar(true)

if 2 == ndims(Ω)
    grid = DiscontinuousSpectralElementGrid(Ω, elements = (1,1), polynomialorder = (4,4), array = ArrayType)
else
    grid = DiscontinuousSpectralElementGrid(Ω, elements = (3,2,3), polynomialorder = (3,3,3), array = ArrayType)
end

x, y, z = coordinates(grid)
nrealelem = size(x)[2] # fix this later to only depend on grid intrinsically
ijksize = prod(polynomialorders(grid) .+ 1)
device = array_device(x)
dim = ndims(Ω)
N = round(Int, size(x)[1]^(1/dim)) - 1
dependencies = nothing

# Initialize State s
Q  = MPIStateArray{FT}(mpicomm, ArrayType, ijksize, nrealelem, 1)
∇Q = MPIStateArray{FT}(mpicomm, ArrayType, ijksize, nrealelem, 3)
exact_∇Q = copy(∇Q)
cartesian_∇Q = copy(∇Q) .* 0.0
##
event = launch_volume_gradient!(∇Q, Q, grid)
wait(event)

## Test Block 1: Volume Test, gradient
a = 1
b = 1
c = 1
@. Q.realdata[:,:, 1] = a * sin(π*x) + b * sin(π*y) + c * sin(π*z)
@. exact_∇Q.realdata[:,:, 1] = a * π*cos(π*x)
@. exact_∇Q.realdata[:,:, 2] = b * π*cos(π*y)
@. exact_∇Q.realdata[:,:, 3] = c * π*cos(π*z)

∇!(cartesian_∇Q, Q, grid)

event = launch_volume_gradient!(∇Q, Q, grid)
wait(event)
tol = eps(1e5)
L∞(x) = maximum(abs.(x))
println(L∞(∇Q - exact_∇Q))
@testset "Gradient Test" begin
    @test L∞(∇Q - cartesian_∇Q) < tol
end

tol = 0.68 # small element sizes with p or h refinement this will get smaller
@testset "Exact Gradient Test" begin
    @test L∞(∇Q - exact_∇Q) < tol
end

## Test Block 2: Interface Test

event = launch_interface_gradient!(∇Q, Q, grid,)
wait(event)

tol = eps(1000.0)
@testset "Gradient Interface Test" begin
    @test L∞(∇Q) < tol
end

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
# initial.jl

import Requests: post, json

const RequireToInitial = "You need to initial"

request(app::App, verb::AbstractString, pair::Pair) = request(app, verb, Dict(pair))
function request(app::App, verb::AbstractString, param::Dict)
    if isempty(app.url)
        Dict("typ"=>"symbol", "value"=>RequireToInitial)
    else
        json(post("$(app.url)$(verb)"; json = param))
    end
end

"""
    initial(url)::QueryResult

Create an endpoint by connecting to the AppConsole Server
"""
initial(url::AbstractString) = initial(App(url))
function initial(app::App)
    verb = "/initial"
    param = Dict()
    dict = request(app, verb, param)
    result = QueryResult(dict, app, verb, param)
    set_endpoint(result)
    result
end

"""
    set_endpoint(result)

Set the global variable ``current_endpoint`` using the query result
"""
function set_endpoint(result::QueryResult)
    global current_endpoint
    current_endpoint = result.app
end

"""
    valueof(result)

Return the value of result
"""
function valueof(result::QueryResult)
    return result.info.value
end


text:
mutable struct RhoComparison <: AbstractModel
    params::NamedTuple{(:env_alpha,:patch_aplha,:beta,:reset,:bias), NTuple{5,Float64}}
    limits::NTuple{5,Float64}
end

function init(::Type{WaltonForaging.RhoComparison})
    return RhoComparison(
        (env_alpha = rand(0.0:0.0001:0.1), #=alpha_environment = [0,0.1] learning rate of environment=#
            patch_aplha = rand(0.0:0.0001:0.2), #=alpha_patch = [0, 0.2] learning rate of the patch=#
            beta = rand(2:0.0001:10), #=beta = [2, 10] noise parameter or inverse temperature=#
            reset = rand(0.0:0.0001:0.5), #=reset = [0, 0.5] patch env rew rate at the entrance of a trials, because it has decayed in the last trial=#
            bias = 0.1 #=bias = [2, 10] bias towards staying in the softmax=#
        ),
        (1.0, 1.0, Inf64, 1.0, Inf64)
    )
end

Distributions.params(m::RhoComparison) = collect(m.params)

function updateparams!(m::RhoComparison,v::Vector{Float64})
    (length(v) == 5) || error("Incorrect number of params for $(typeof(m))")
    env_alpha,patch_aplha,beta,reset,bias = tuple(v...)
    m.params = (;env_alpha,patch_aplha,beta,reset,bias)
end

function Pstay(m::RhoComparison, env_rew_rate::Float64, patch_rew_rate::Float64)
    env_alpha, patch_alpha, beta, reset_val, bias = tuple(params(m)...)
    return exp(bias + beta*patch_rew_rate)/(exp(bias + beta*patch_rew_rate) + exp(beta*env_rew_rate))
end

function Pleave(m::RhoComparison, env_rew_rate::Float64, patch_rew_rate::Float64)
    env_alpha, patch_alpha, beta, reset_val, bias = tuple(params(m)...)
    # return 1 - exp(bias + beta*patch_rew_rate)/exp(bias + beta*patch_rew_rate) + exp(beta*env_rew_rate)
    return exp(beta*env_rew_rate)/(exp(beta*env_rew_rate) + exp(bias + beta*patch_rew_rate))
end

function Poutcome(m::RhoComparison, env_rew_rate::Float64, patch_rew_rate::Float64, leave::Int64)
    if leave == 1
        Pleave(m, env_rew_rate, patch_rew_rate)
    elseif leave == 0
        Pstay(m, env_rew_rate, patch_rew_rate)
    else
        error("Wrong argument for Poutcome calculation: $leave")
    end
end

function Likelihood(m::RhoComparison,df)
    # if I start using subdf for patches conditions I will need:
    # an initial env reward rate before splitting per session (maybe a separate variable to feed to loglikelihood)
    # how do I minimize the same patch reset value for the whole session?

    env_alpha, patch_alpha, beta, reset_val, bias = tuple(params(m)...)
    final_env_rho = sum(df.REWARD/df[end,:TIME])
    transform!(df,[:REWARD,:BIN,:ENV_INITIALVALUE,:NEWSESSION] => ((r,b,i,n) -> env_get_rew_rate(r,b,i,n,env_alpha)) => :EnvRewRate)
    transform!(df, [:REWARD,:BIN,:NEWTRIAL] => ((r,b,n) -> patch_get_rew_rate(r, b, reset_val, n, patch_alpha)) => :PatchRewRate)
    transform!(df, [:EnvRewRate,:PatchRewRate,:LEAVE] => ByRow((e,p,l) -> Poutcome(m,e,p,l)) => :Probability)
    return filter(r -> r.SIDE == "travel", df).Probability
end

function MixedModels.simulate(m::RhoComparison,df)
    transform(df, :Probability => (x-> convert(Float64, x < rand())) => :SimulatedModel)
end


text:
function DiffEqBase.unwrap_cache(integrator::SDEIntegrator, is_stiff)
  alg   = integrator.alg
  cache = integrator.cache
  iscomp = alg isa StochasticCompositeAlgorithm
  if !iscomp
    return cache
  elseif alg.choice_function isa AutoSwitch
    num = is_stiff ? 2 : 1
    return cache.caches[num]
  else
    return cache.caches[integrator.cache.current]
  end
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related with string searching.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "String search" begin
    str = """
        Test 1 😅 \e[38;5;231;48;5;243mTest 2 😅 \e[38;5;201;48;5;243mTest\e[0m
        Test 1 😅 \e[38;5;231;48;5;243mTest 2 😅 \e[38;5;201;48;5;243mTest\e[0m
        """

    expected = [
        (11, 9),
        (35, 9)
    ]
    search_result = string_search(str, r"Test 2 😅")
    @test search_result == expected
end

@testset "String search by line" begin
    str = """
        Test 1 😅 \e[38;5;231;48;5;243mTest 2 😅 \e[38;5;201;48;5;243mTest\e[0m
        Test 1 😅 \e[38;5;231;48;5;243mTest 2 😅 \e[38;5;201;48;5;243mTest\e[0m
        Test 1 😅 \e[38;5;231;48;5;243mTest 2 😅 \e[38;5;201;48;5;243mTest\e[0m
        Test 1 😅 \e[38;5;231;48;5;243mTest 2 😅 \e[38;5;201;48;5;243mTest\e[0m
        Test 1 😅 \e[38;5;231;48;5;243mTest 2 😅 \e[38;5;201;48;5;243mTest\e[0m
        """

    search_result = string_search_per_line(str, r"Test 2 😅")

    @test search_result[1] == [(11, 9)]
    @test search_result[2] == [(11, 9)]
    @test search_result[3] == [(11, 9)]
    @test search_result[4] == [(11, 9)]
    @test search_result[5] == [(11, 9)]
    @test length(search_result) == 5
end


text:
# Line drawing algorithm

"""
    drawline(func::Function, ci1::CartesianIndex{2}, ci2::CartesianIndex{2}, eachstep::Bool=false)

At each step along the line from `ci1` to `ci2` call `func` with a single argument, the row, column coordinates
of the pixel as a `Tuple{Int, Int}`.

If `eachstep=true` then `func` is called once when the row changes and once when the column changes.  When
`eachstep=false` then `func` is called less frequently and both row and col may change between calls.
"""
function drawline(
    func::Function,
    ci1::CartesianIndex{2},
    ci2::CartesianIndex{2},
    eachstep::Bool = false,
)
    r, c = ci1.I
    dr, dc = -abs(ci2.I[1] - ci1.I[1]), abs(ci2.I[2] - ci1.I[2])
    sr, sc = ci1.I[1] < ci2.I[1] ? 1 : -1, ci1.I[2] < ci2.I[2] ? 1 : -1
    err = dc + dr
    func(ci1.I)
    while !((c == ci2.I[2]) && (r == ci2.I[1]))
        e2 = 2 * err
        if e2 >= dr
            err += dr
            c += sc
            if eachstep
                func((r, c))
            end
        end
        if e2 <= dc
            err += dc
            r += sr
            if eachstep
                func((r, c))
            end
        end
        if !eachstep
            func((r, c))
        end
    end
end

"""
    drawray(
        func::Function,  # returns Bool
        ci1::CartesianIndex{2}, # starting point
        ci2::CartesianIndex{2}, # towards this point
        bounds::CartesianIndices{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}},  # stay within
        eachstep::Bool=false
    )
    function drawray(
        func::Function,
        ci1::CartesianIndex{2},
        ci2::CartesianIndex{2},
        eachstep::Bool = false,
    )

    Draw a ray starting at `ci1` towards `ci2` and call `func` with a single argument, the row, column coordinates
of the pixel as a `Tuple{Int, Int}`.  `func` must return a `Bool` - true to continue drawing the ray and false to
stop drawing the ray.  The first version checks the bounds each step and the second doesn't.

If `eachstep=true` then `func` is called once when the row changes and once when the column changes.  When
`eachstep=false` then `func` is called less frequently and both row and col may change between calls.
"""

function drawray(
    func::Function,
    ci1::CartesianIndex{2},
    ci2::CartesianIndex{2},
    bounds::CartesianIndices{2, Tuple{Vararg{OrdinalRange{Int64, Int64}, 2}}},
    eachstep::Bool = false,
)
    r, c = ci1.I
    dr, dc = -abs(ci2.I[1] - ci1.I[1]), abs(ci2.I[2] - ci1.I[2])
    sr, sc = ci1.I[1] < ci2.I[1] ? 1 : -1, ci1.I[2] < ci2.I[2] ? 1 : -1
    err = dc + dr
    inb(r,c) = ((r in bounds.indices[1]) && (c in bounds.indices[2]))
    if !(inb(r,c) && func((r,c)))
        return nothing
    end
    while true
        e2 = 2 * err
        if e2 >= dr
            err += dr
            c += sc
            if eachstep && !(inb(r,c) && func((r, c)))
                return nothing
            end
        end
        if e2 <= dc
            err += dc
            r += sr
            if eachstep && !(inb(r,c) && func((r, c)))
                return nothing
            end
        end
        if (!eachstep) && !(inb(r,c) && func((r, c)))
            return nothing
        end

    end
end

function drawray(
    func::Function,
    ci1::CartesianIndex{2},
    ci2::CartesianIndex{2},
    eachstep::Bool = false,
)
    r, c = ci1.I
    dr, dc = -abs(ci2.I[1] - ci1.I[1]), abs(ci2.I[2] - ci1.I[2])
    sr, sc = ci1.I[1] < ci2.I[1] ? 1 : -1, ci1.I[2] < ci2.I[2] ? 1 : -1
    err = dc + dr
    if !func((r,c))
        return nothing
    end
    while true
        e2 = 2 * err
        if e2 >= dr
            err += dr
            c += sc
            if eachstep && !func((r, c))
                return nothing
            end
        end
        if e2 <= dc
            err += dc
            r += sr
            if eachstep && !func((r, c))
                return nothing
            end
        end
        if (!eachstep) && !func((r, c))
            return nothing
        end

    end
end

max_stars_repo_name: NicholasWMRitchie/NeXLSpectrum

max_issues_repo_name: NicholasWMRitchie/NeXLSpectrum

max_forks_repo_name: NicholasWMRitchie/NeXLSpectrum


text:
using Maple

maps = filter(x -> endswith(x, ".bin"), readdir("."))

for mapFile in maps
  println("loading map: $(mapFile)")
  mapleMap = loadMap(mapFile)
  correctPackage = mapFile[1:end-4]
  println("map loaded")
  println("package: $(mapleMap.package)")
  println("correct package: $(correctPackage)")
  mapleMap.package = correctPackage
  println("package changed to $(mapleMap.package)")
  encodeMap(mapleMap, mapFile)
  println("saved fixed map: $(mapFile)")
end

max_stars_repo_name: npdmfixup/Fuji

max_issues_repo_name: npdmfixup/Fuji

max_forks_repo_name: npdmfixup/Fuji


text:
using AlgebraicDynamics.UWDDynam
using AlgebraicDynamics.DWDDynam

using OrdinaryDiffEq

using Base.Iterators

"""    trajectory(r::DiscreteResourceSharer, u0::AbstractVector, p, nsteps::Int; dt::Int = 1)

Evolves the resouce sharer `r` for `nsteps` times with step size `dt`, initial condition `u0`, and parameters `p`.
"""
function trajectory(r::DiscreteResourceSharer, u0::AbstractVector, p, T::Int; dt::Int= 1)
  prob = DiscreteProblem(r, u0, (0, T), p)
  sol = solve(problem, FunctionMap(); dt = dt)
  return sol.xs
end

"""    trajectory(m::DiscreteMachine, u0::AbstractVector, xs::AbstractVector, p, nsteps::Int; dt::Int = 1)

Evolves the machine `m` for `nsteps` times with step size `dt`, initial condition `u0`, and parameters `p`. Any inputs to `m` are determied by `xs`. If `m` has no inputs then you can omit `xs`.
"""
function trajectory(m::DiscreteMachine, u0::AbstractVector, xs,  p, T::Int; dt::Int= 1)
  prob = DiscreteProblem(m, u0, xs, (0, T), p)
  sol = solve(prob, FunctionMap(); dt = dt)
  return sol.xs
end


text:
module JLD2DebugTools

using Reexport
@reexport using JLD2

let
    jlnames = names(JLD2; all=true)
    filter!(n -> !startswith(string(n), '#'), jlnames)
    for name in jlnames
        @eval using JLD2: $name
    end
end

export describe, committed_datatypes

function Base.show(io::IO, roffset::RelOffset)
    print(io, "RelOffset($(Int(roffset.offset)))")
end

include("datasets.jl")
############################################################################################
##                                      describe                                          ##
############################################################################################

describe(f::JLDFile, name::AbstractString) =
    describe(f.root_group, name)

function describe(g::Group, name::AbstractString)
    f = g.f
    f.n_times_opened == 0 && throw(ArgumentError("file is closed"))

    (g, name) = pathize(g, name, false)

    roffset = lookup_offset(g, name)
    if roffset == UNDEFINED_ADDRESS
        haskey(g.unwritten_child_groups, name) && return g.unwritten_child_groups[name]
        throw(KeyError(name))
    end

    if isgroup(f, roffset)
        let loaded_groups = f.loaded_groups
            get!(()->load_group(f, roffset), loaded_groups, roffset)
        end
    else
        describe_dataset(f, roffset, name)
    end
end

function describe_dataset(f::JLDFile, offset::RelOffset, name)
    io = f.io
    seek(io, fileoffset(f, offset))
    cio = begin_checksum_read(io)
    sz = read_obj_start(cio)
    pmax = position(cio) + sz

    # Messages
    dataspace = ReadDataspace()
    attrs = EMPTY_READ_ATTRIBUTES
    datatype_class::UInt8 = 0
    datatype_offset::Int64 = 0
    data_offset::Int64 = 0
    data_length::Int = -1
    chunked_storage::Bool = false
    filter_id::UInt16 = 0

    storage_type = 0xFF
    while position(cio) <= pmax-4
        msg = jlread(cio, HeaderMessage)
        endpos = position(cio) + msg.size
        if msg.msg_type == HM_DATASPACE
            dataspace = read_dataspace_message(cio)
        elseif msg.msg_type == HM_DATATYPE
            datatype_class, datatype_offset = read_datatype_message(cio, f, (msg.flags & 2) == 2)
        elseif msg.msg_type == HM_FILL_VALUE
            (jlread(cio, UInt8) == 3 && jlread(cio, UInt8) == 0x09) || throw(UnsupportedFeatureException())
        elseif msg.msg_type == HM_DATA_LAYOUT
            jlread(cio, UInt8) == 4 || throw(UnsupportedVersionException())
            storage_type = jlread(cio, UInt8)
            if storage_type == LC_COMPACT_STORAGE
                data_length = jlread(cio, UInt16)
                data_offset = position(cio)
            elseif storage_type == LC_CONTIGUOUS_STORAGE
                data_offset = fileoffset(f, jlread(cio, RelOffset))
                data_length = jlread(cio, Length)
            elseif storage_type == LC_CHUNKED_STORAGE
                # TODO: validate this
                flags = jlread(cio, UInt8)
                dimensionality = jlread(cio, UInt8)
                dimensionality_size = jlread(cio, UInt8)
                skip(cio, Int(dimensionality)*Int(dimensionality_size))

                chunk_indexing_type = jlread(cio, UInt8)
                chunk_indexing_type == 1 || throw(UnsupportedFeatureException("Unknown chunk indexing type"))
                data_length = jlread(cio, Length)
                jlread(cio, UInt32)
                data_offset = fileoffset(f, jlread(cio, RelOffset))
                chunked_storage = true
            else
                throw(UnsupportedFeatureException("Unknown data layout"))
            end
        elseif msg.msg_type == HM_FILTER_PIPELINE
            version = jlread(cio, UInt8)
            version == 2 || throw(UnsupportedVersionException("Filter Pipeline Message version $version is not implemented"))
            nfilters = jlread(cio, UInt8)
            nfilters == 1 || throw(UnsupportedFeatureException())
            filter_id = jlread(cio, UInt16)
            issupported_filter(filter_id) || throw(UnsupportedFeatureException("Unknown Compression Filter $filter_id"))
        elseif msg.msg_type == HM_ATTRIBUTE
            if attrs === EMPTY_READ_ATTRIBUTES
                attrs = ReadAttribute[read_attribute(cio, f)]
            else
                push!(attrs, read_attribute(cio, f))
            end
        elseif (msg.flags & 2^3) != 0
            throw(UnsupportedFeatureException())
        end
        seek(cio, endpos)
    end
    seek(cio, pmax)

    filter_id != 0 && !chunked_storage && throw(InvalidDataException("Compressed data must be chunked"))

    # Checksum
    end_checksum(cio) == jlread(io, UInt32) || throw(InvalidDataException("Invalid Checksum"))

    # TODO verify that data length matches
    #val = read_data(f, dataspace, datatype_class, datatype_offset, data_offset, data_length,
     #               filter_id, offset, attrs)
    #val

    #println(dataspace)
    #println(datatype_class)
    #println(datatype_offset)
    #println(data_length)
    #println(attrs)

    if datatype_class == typemax(UInt8) # Committed datatype
        typestr = stringify_committed_datatype(f, f.datatype_locations[h5offset(f, datatype_offset)])
    else
        seek(io, datatype_offset)
        @read_datatype io datatype_class dt begin
            rr = jltype(f, dt)
            typestr = string(typeof(rr).parameters[1])
            #seek(io, data_offset)
            #read_dataspace = (dataspace, NULL_REFERENCE, data_length, filter_id)
            #data = read_data(f, rr, read_dataspace, attrs)
        end
    end
    println("Dataset \"$name\" at position: $offset")
    storage_type_str = Dict(0x00 => "Compact Storage",
                            0x01 => "Contiguous Storage",
                            0x02 => "Chunked Storage")[storage_type]
    dspace = Dict(0x00 => "Scalar", 0x01 => "Simple", 0x02 => "Null")[dataspace.dataspace_type]

    println("Storage type: $storage_type_str")
    println("Dataspace: $dspace")

    ndims, ndim_offset = get_ndims_offset(f, dataspace, attrs)
    seek(io, ndim_offset)
    dims = array_dims(io, ndims)
    println("Dimensions: $dims")

    if ndims > 0
        println("Type signature: Array{$typestr, $(Int(ndims))}")
    else
        println("Type signature: $typestr")
    end
    if filter_id != 0
        println("Compressed with $(ID_TO_DECOMPRESSOR[filter_id][2])")
    end
    println("Bytes on disk: $data_length")
    #= if datatype_class != typemax(UInt8) # basic datatype
        print("Value")
        print(stdout, data)
    end =#
    return nothing
end

function array_dims(io::IO, ndims)
    ds = tuple(reverse!(read!(io, Vector{Int64}(undef, ndims)))...)
end

export committed_datatypes
function committed_datatypes(f::JLDFile)
    if !haskey(f, "_types")
        println("No committed types")
    end
    for (n, (offset,cdt)) in enumerate(f.datatype_locations)
        # skip Datatype
        #= if cdt.index == 1
            continue
        end =#
        print("$n) ")
        typestr = stringify_committed_datatype(f, cdt; showfields=true)
        #println("$n)\n$typestr")
        println(typestr)
    end
end

function stringify_committed_datatype(f, cdt; showfields=false)
    io = f.io
    dt, attrs = read_committed_datatype(f, cdt)
    written_type_str = ""
    julia_type_str = ""
    for attr in attrs
        if !(attr.name == :julia_type || attr.name == :written_type)
            continue
        end
        datatype_offset = attr.datatype_offset
        data_offset = attr.data_offset
        rr = jltype(f, f.datatype_locations[h5offset(f, datatype_offset)])
        @assert attr.dataspace.dataspace_type == DS_SCALAR
        seek(io, data_offset)
        str = jlconvert_string(rr, f, io.curptr)
        if attr.name == :written_type
            written_type_str = str
        elseif attr.name ==:julia_type
            julia_type_str = str
        end
    end

#=     end
    attr = first(attrs)
    datatype_offset = attr.datatype_offset
    data_offset = attr.data_offset

    rr = jltype(f, f.datatype_locations[h5offset(f, datatype_offset)])
    @assert attr.dataspace.dataspace_type == DS_SCALAR

    seek(io, data_offset)
    typestr = jlconvert_string(rr, f, io.curptr)
 =#
    #@show typestr
    if !showfields ||
        dt isa BasicDatatype ||
        dt isa VariableLengthDatatype ||
        startswith(julia_type_str, r"Tuple|Union") ||
        julia_type_str == "DataType"
        return julia_type_str
    end
    typestr = "struct "*julia_type_str
    if !isempty(written_type_str)
        typestr = "julia type: $typestr\n"*
        "written type:\n"*"struct $written_type_str"
    end
    field_datatypes = read_field_datatypes(f, attrs)
    #rodr = reconstruct_odr(f, dt, field_datatypes)
    #types = typeof(rodr).parameters[2].parameters
    # Get the type and ODR information for each field
    #types = Vector{Any}(undef, length(dt.names))
    #h5types = Vector{Any}(undef, length(dt.names))
    for i = 1:length(dt.names)
        if !isempty(field_datatypes) && (ref = field_datatypes[i]) != NULL_REFERENCE
            #dtrr = jltype(f, f.datatype_locations[ref])
            fieldtype = stringify_committed_datatype(f, f.datatype_locations[ref])
        else
            # These are normal julia types
            #@show dt.members[i]
            dtrr = jltype(f, dt.members[i])
            fieldtype = string(typeof(dtrr).parameters[1])
            #field = stringify_committed_datatype(f, dt.members[i])

        end
        #types[i], h5types[i] = typeof(dtrr).parameters
        #println(dt.names[i],"::",field)
        typestr *= "\n\t"*"$(dt.names[i])::$(fieldtype)"
    end
    #if length(dt.names) > 0

    typestr = typestr*"\n"*"end"
    #end
    #@show field_datatypes
    #@show rodr
    #@show types

#=     structexpr =Expr(:struct, false, reconname,
                   Expr(:block, Any[ Expr(Symbol("::"), dt.names[i], types[i]) for i = 1:length(types) ]...,
                        # suppress default constructors, plus a bogus `new()` call to make sure
                        # ninitialized is zero.
                        Expr(:if, false, Expr(:call, :new))))
    @info structexpr
 =#
    return typestr
end

function stringify_object(f, offset)
    io = f.io
    seek(io, fileoffset(f, offset))
    cio = begin_checksum_read(io)
    sz = read_obj_start(cio)
    pmax = position(cio) + sz

    # Messages
    dataspace = ReadDataspace()
    attrs = EMPTY_READ_ATTRIBUTES
    datatype_class::UInt8 = 0
    datatype_offset::Int64 = 0
    data_offset::Int64 = 0
    data_length::Int = -1
    chunked_storage::Bool = false
    filter_id::UInt16 = 0
    while position(cio) <= pmax-4
        msg = jlread(cio, HeaderMessage)
        endpos = position(cio) + msg.size
        if msg.msg_type == HM_DATASPACE
            dataspace = read_dataspace_message(cio)
        elseif msg.msg_type == HM_DATATYPE
            datatype_class, datatype_offset = read_datatype_message(cio, f, (msg.flags & 2) == 2)
        elseif msg.msg_type == HM_FILL_VALUE
            (jlread(cio, UInt8) == 3 && jlread(cio, UInt8) == 0x09) || throw(UnsupportedFeatureException())
        elseif msg.msg_type == HM_DATA_LAYOUT
            jlread(cio, UInt8) == 4 || throw(UnsupportedVersionException())
            storage_type = jlread(cio, UInt8)
            if storage_type == LC_COMPACT_STORAGE
                data_length = jlread(cio, UInt16)
                data_offset = position(cio)
            elseif storage_type == LC_CONTIGUOUS_STORAGE
                data_offset = fileoffset(f, jlread(cio, RelOffset))
                data_length = jlread(cio, Length)
            elseif storage_type == LC_CHUNKED_STORAGE
                # TODO: validate this
                flags = jlread(cio, UInt8)
                dimensionality = jlread(cio, UInt8)
                dimensionality_size = jlread(cio, UInt8)
                skip(cio, Int(dimensionality)*Int(dimensionality_size))

                chunk_indexing_type = jlread(cio, UInt8)
                chunk_indexing_type == 1 || throw(UnsupportedFeatureException("Unknown chunk indexing type"))
                data_length = jlread(cio, Length)
                jlread(cio, UInt32)
                data_offset = fileoffset(f, jlread(cio, RelOffset))
                chunked_storage = true
            else
                throw(UnsupportedFeatureException("Unknown data layout"))
            end
        elseif msg.msg_type == HM_FILTER_PIPELINE
            version = jlread(cio, UInt8)
            version == 2 || throw(UnsupportedVersionException("Filter Pipeline Message version $version is not implemented"))
            nfilters = jlread(cio, UInt8)
            nfilters == 1 || throw(UnsupportedFeatureException())
            filter_id = jlread(cio, UInt16)
            issupported_filter(filter_id) || throw(UnsupportedFeatureException("Unknown Compression Filter $filter_id"))
        elseif msg.msg_type == HM_ATTRIBUTE
            if attrs === EMPTY_READ_ATTRIBUTES
                attrs = ReadAttribute[read_attribute(cio, f)]
            else
                push!(attrs, read_attribute(cio, f))
            end
        elseif (msg.flags & 2^3) != 0
            throw(UnsupportedFeatureException())
        end
        seek(cio, endpos)
    end
    seek(cio, pmax)

    filter_id != 0 && !chunked_storage && throw(InvalidDataException("Compressed data must be chunked"))

    # Checksum
    end_checksum(cio) == jlread(io, UInt32) || throw(InvalidDataException("Invalid Checksum"))

    # TODO verify that data length matches
    #val = read_data(f, dataspace, datatype_class, datatype_offset, data_offset, data_length,
    #                filter_id, offset, attrs)
    # An excerpt from the commented function

    if datatype_class == typemax(UInt8) # Committed datatype
        rr = jltype(f, f.datatype_locations[h5offset(f, datatype_offset)])
        seek(io, data_offset)
        typestr = jlconvert_string(rr, f, io.curptr)
    else
        seek(io, datatype_offset)
        @read_datatype io datatype_class dt begin
            rr = jltype(f, dt)
            seek(io, data_offset)
            read_dataspace = (dataspace, NULL_REFERENCE, data_length, filter_id)
            res = read_data(f, rr, read_dataspace, nothing)
            string(res)
        end
    end
end

function typestring_from_refs(f::JLDFile, ptr::Ptr)
    # Test for a potential null pointer indicating an empty array
    isinit = jlunsafe_load(convert(Ptr{UInt32}, ptr)) != 0
    if isinit
        refs = jlconvert(ReadRepresentation{RelOffset, Vlen{RelOffset}}(), f, ptr, NULL_REFERENCE)
        #println("datatypes at refs $(Int.(getproperty.(refs,:offset)))")
        params =  Any[let
            # If the reference is to a committed datatype, read the datatype
            nulldt = CommittedDatatype(UNDEFINED_ADDRESS, 0)
            cdt = get(f.datatype_locations, ref, nulldt)
            res = if cdt !== nulldt
                stringify_committed_datatype(f, cdt)
            else
                stringify_object(f, ref)
            end
            if startswith(res, r"Core|Main")
                res = last(split(res, "."; limit=2))
            end
            res
        end for ref in refs]
        return params
    end
    return []
end

function jlconvert_string(rr::ReadRepresentation{T,DataTypeODR()},
                        f::JLDFile,
                        ptr::Ptr) where T
    mypath = String(jlconvert(ReadRepresentation{UInt8,Vlen{UInt8}}(), f, ptr, NULL_REFERENCE))
    params = typestring_from_refs(f, ptr+odr_sizeof(Vlen{UInt8}))
    if startswith(mypath, r"Core|Main")
        mypath = last(split(mypath, "."; limit=2))
    end

    if !isempty(params)
        return mypath*"{"*interleave(params, ",")*"}"
    else
        return mypath
    end
end

function jlconvert_string(::ReadRepresentation{Union, UnionTypeODR()}, f::JLDFile,
    ptr::Ptr)#, header_offset::RelOffset)
    # Skip union type description in the beginning
    ptr += odr_sizeof(Vlen{String})
    # Reconstruct a Union by reading a list of DataTypes and UnionAlls
    # Lookup of RelOffsets is taken from jlconvert of DataTypes
    datatypes = typestring_from_refs(f, ptr)
    unionalls = typestring_from_refs(f, ptr+odr_sizeof(Vlen{RelOffset}))
    "Union{"*interleave(vcat(datatypes, unionalls), ",")*"}"
    #= v = Union{datatypes..., unionalls...}
    track_weakref!(f, header_offset, v)
    v =#
end

function interleave(strs, int)
    if isempty(strs)
        return ""
    elseif length(strs) == 1
        return only(strs)
    end
    prod(str*int*" " for str in strs[1:end-1]) * last(strs)
end
end


text:
module Epicrop

using DataFrames, Dates, Interpolations

include("helpers.jl")

export hlipmodel
export bacterialblight
export brownspot
export leafblast
export sheathblight
export tungro

#re-exports from DataFrames.jl
export DataFrame

"""
hlipmodel(
    wth,
    emergence,
    onset,
    duration,
    rhlim,
    rainlim,
    H0,
    I0,
    RcA,
    RcT,
    RcOpt,
    p,
    i,
    Sx,
    a,
    RRS,
    RRG)

Runs a healthy-latent-infectious-postinfectious (HLIP) model using weather data and optimal curve values for respective crop diseases.

## Keywords

- `wth` a data frame of weather on a daily time-step containing data with the following field names.
    | Field | value |
    |-------|-------------|
    |YYYYMMDD | Date as Year Month Day, YYYY-MM-DD, (ISO8601) |
    |DOY |  Consecutive day of year, commonly called "Julian date" |
    |TEMP | Mean daily temperature (°C) |
    |RHUM | Mean daily relative humidity (%) |
    |RAIN | Mean daily rainfall (mm) |
- `emergence`: expected date of plant emergence entered as a `Dates.Date` object. From Table 1 Savary _et al._ 2012.
- `onset` expected number of days until the onset of disease after emergence date. From Table 1 Savary _et al._ 2012.
- `duration`: simulation duration (growing season length). From Table 1 Savary _et al._ 2012.
- `rhlim`: threshold to decide whether leaves are wet or not (usually 90%). From Table 1 Savary _et al._ 2012.
- `rainlim`: threshold to decide whether leaves are wet or not. From Table 1 Savary _et al._ 2012.
- `H0`: initial number of plant's healthy sites. From Table 1 Savary _et al._ 2012.
- `I0`: initial number of infective sites. From Table 1 Savary _et al._ 2012.
- `RcA`: crop age modifier for *Rc* (the basic infection rate corrected for removals). From Table 1 Savary _et al._ 2012.
- `RcT`: temperature modifier for *Rc* (the basic infection rate corrected for removals). From Table 1 Savary _et al._ 2012.
- `RcOpt`: potential basic infection rate corrected for removals. From Table 1 Savary _et al._ 2012.
- `i`: duration of infectious period. From Table 1 Savary _et al._ 2012.
- `p`: duration of latent period. From Table 1 Savary _et al._ 2012.
- `Sx`: maximum number of sites. From Table 1 Savary _et al._ 2012.
- `a`: aggregation coefficient. From Table 1 Savary _et al._ 2012.
- `RRS`: relative rate of physiological senescence. From Table 1 Savary _et al._ 2012.
- `RRG`: relative rate of growth. From Table 1 Savary _et al._ 2012.

## Returns

A `DataFrame` with the model's output. Latitude and longitude are included for mapping purposes if they are present in the input weather data.
"""

function hlipmodel(;
    wth,
    emergence,
    onset,
    duration,
    rhlim,
    rainlim,
    H0,
    I0,
    RcA,
    RcT,
    RcOpt,
    p,
    i,
    Sx,
    a,
    RRS,
    RRG,
)

    if !(typeof(emergence) == Dates.Date)
        throw(DomainError(emergence, "emergence must be a Date object"))
    end

    final_day = emergence + Dates.Day(duration - 1)
    season = Base.collect(emergence:Dates.Day(1):final_day)

    if !(emergence >= wth[1, "YYYYMMDD"] ||
        final_day > Base.findmax(wth[:, "YYYYMMDD"])[1])
        throw(DomainError(wth, "incomplete weather data or dates do not align"))
    end

    if (H0 < 0)
        throw(DomainError(H0,
            "H0 cannot be < 0, check your initial number of healthy sites"))
    end

    if (I0 < 0)
        throw(DomainError(I0,
            "I0 cannot be < 0, check your initial number of infective sites"))
    end

    season_wth = wth[Base.in(season - Dates.Day(1)).(wth.YYYYMMDD), :]

    return(
        _hliploop(
            season = season,
            season_wth = season_wth,
            onset = onset,
            duration = duration,
            rhlim = rhlim,
            rainlim = rainlim,
            H0 = H0,
            I0 = I0,
            RcA = RcA,
            RcT = RcT,
            RcOpt = RcOpt,
            p = p,
            i = i,
            Sx = Sx,
            a = a,
            RRS = RRS,
            RRG = RRG
        )
    )
end

function _hliploop(;
    season,
    season_wth,
    onset,
    duration,
    rhlim,
    rainlim,
    H0,
    I0,
    RcA,
    RcT,
    RcOpt,
    p,
    i,
    Sx,
    a,
    RRS,
    RRG)

    # Create `infday` for use in for loop below.
    infday = 0

    # Create vectors of preallocated output variables.
    cofr = Base.zeros(duration)
    rc = Base.zeros(duration)
    RHCoef = Base.zeros(duration)
    latency = Base.zeros(duration)
    infectious = Base.zeros(duration)
    intensity = Base.zeros(duration)
    rsenesced = Base.zeros(duration)
    rgrowth = Base.zeros(duration)
    rtransfer = Base.zeros(duration)
    infection = Base.zeros(duration)
    diseased = Base.zeros(duration)
    senesced = Base.zeros(duration)
    removed = Base.zeros(duration)
    now_infectious = Base.zeros(duration)
    now_latent = Base.zeros(duration)
    sites = Base.zeros(duration)
    total_sites = Base.zeros(duration)
    Rc_age = _fn_rc(RcA, 1:duration)
    Rc_temp = _fn_rc(RcT, season_wth[!, :TEMP])

    for d in 1:duration
        d_1 = d - 1

        # Start of the state calculations.
        if d == 1
            # Start crop growth.
            sites[d] = H0
            rsenesced[d] = RRS * sites[d]
        else
            if d > i
                removed_today = infectious[infday+1]
            else
                removed_today = 0
            end

            sites[d] = sites[d_1] + rgrowth[d_1] - infection[d_1] - rsenesced[d_1]
            rsenesced[d] = removed_today + RRS * sites[d]
            senesced[d] = senesced[d_1] + rsenesced[d_1]

            latency[d] = infection[d_1]
            latday = d - p
            latday = Base.max(1, latday)
            now_latent[d] = Base.sum(latency[latday:d])

            infectious[d] = rtransfer[d_1]
            infday = d - i
            infday = Base.max(1, infday)
            now_infectious[d] = Base.sum(infectious[infday:d])
        end

        if (season_wth[!, :RHUM][d] >= rhlim || season_wth[!, :RAIN][d] >= rainlim)
            RHCoef[d] = 1
        end

        rc[d] = RcOpt * (Rc_age[d] * Rc_temp[d] * RHCoef[d])
        diseased[d] = sum(infectious) + now_latent[d] + removed[d]
        removed[d] = sum(infectious) - now_infectious[d]

        cofr[d] = 1 - (diseased[d] / (sites[d] + diseased[d]))

        # Initialise disease.
        if d == onset
            infection[d] = I0
        elseif d > onset
            infection[d] = now_infectious[d] * rc[d] * (cofr[d]^a)
        else
            infection[d] = 0
        end

        if d >= p
            rtransfer[d] = latency[latday]
        else
            rtransfer[d] = 0
        end

        total_sites[d] = diseased[d] + sites[d]
        rgrowth[d] = RRG * sites[d] * (1 - (total_sites[d] / Sx))
        intensity[d] = (diseased[d] - removed[d]) / (total_sites[d] - removed[d])
    end

    res = DataFrames.DataFrame(
        simday = 1:duration,
        dates = season,
        sites = sites,
        latent = now_latent,
        infectious = now_infectious,
        removed = removed,
        senseced = senesced,
        rateinf = infection,
        rtransfer = rtransfer,
        rgrowth = rgrowth,
        rsenesced = rsenesced,
        diseased = diseased,
        intensity = intensity,
        audpc = _audpc(intensity)
    )

    if hasproperty(season_wth, "LAT") && hasproperty(season_wth, "LON")
        res = DataFrames.DataFrame(
            insertcols!(res, :lat => season_wth[1, "LAT"], :lon => season_wth[1, "LON"])
        )
    end

    return res
end

function _fn_rc(Rc, x)
    itp = Interpolations.LinearInterpolation(Rc[:, 1], Rc[:, 2], extrapolation_bc = 0)
    x = itp.(x)
    return x
end

function _audpc(intensity)
    n = length(intensity) - 1
    intvec = Base.zeros(n)
    out = 0.0

    for i in 1:n
        intvec[i] = (intensity[i] + intensity[i + 1]) / 2
        out = sum(intvec)
    end

    return out
end

end # module


text:
"""
    module Structs

Contains all the structs in ModiaMath
"""
module TypesAndStructs

# Types
export AbstractSimulationModel
"""
    type ModiaMath.AbstractSimulationModel

Struct that is used as simulation model (has field: simulationState)
"""
abstract type AbstractSimulationModel end

export AbstractComponentWithVariables
"""
    type ModiaMath.AbstractComponentWithVariables

Struct that contains ModiaMath.AbstractVariables as field or as field
in a sub-struct.
"""
abstract type AbstractComponentWithVariables end

export AbstractComponentInternal
" The internal part of a component (has at least fields \"name\" and \"within\") "
abstract type AbstractComponentInternal end

export AbstractVariable
"""
    type ModiaMath.AbstractVariable <: ModiaMath.AbstractComponentWithVariables

A Variable used as element of the DAE model description and is
included in the result (if no residue)
"""
abstract type AbstractVariable <: AbstractComponentWithVariables end

export AbstractRealVariable
"""
    ModiaMath.AbstractRealVariable <: ModiaMath.AbstractVariable

A real [`ModiaMath.AbstractVariable`](@ref) (either scalar or array)
"""
abstract type AbstractRealVariable <: AbstractVariable end

# Use by Logger only
export Logger
"""
    mutable struct Logger - Log model evaluations
"""
mutable struct Logger
    log::Bool           # = true, if logging

    # log categories
    statistics::Bool
    progress::Bool
    infos::Bool
    warnings::Bool
    events::Bool

    Logger() = new(false, true, true, true, true, true)
end

export SimulationStatistics
"""
    mutable struct SimulationStatistics - Collect statistics of the last simulation run.

The following data is stored in this structure:

- `structureOfDAE`: Structure of DAE
- `cpuTimeInitialization`: CPU-time for initialization
- `cpuTimeIntegration`: CPU-time for integration
- `startTime`: start time of the integration
- `stopTime`: stop time of the integration
- `interval`: communication interval of the integration
- `tolerance`: relative tolerance used for the integration
- `nEquations`: number of equations (length of y and of yp)
- `nConstraints`: number of constraint equations
- `nResults`: number of time points stored in result data structure
- `nSteps`: number of (successful) steps
- `nResidues`: number of calls to compute residues (includes residue calls for Jacobian)
- `nZeroCrossing`: number of calls to compute zero crossings
- `nJac`: number of calls to compute Jacobian
- `nTimeEvents`: number of time events
- `nRestartEvents`: number of events with integrator restart
- `nErrTestFails`: number of fails of error tests
- `h0`: stepsize used at the first step
- `hMin`: minimum integration stepsize
- `hMax`: maximum integration stepsize
- `orderMax`: maximum integration order
- `sparseSolver` = true: if sparse solver used, otherwise dense solver
- `nGroups`: if sparseSolver, number of column groups to compute Jacobian
  (e.g. if nEquations=100, nGroups=5, then 5+1=6 model evaluations are needed
  to compute the Jacobian numerically, instead of 101 model evaluations without
  taking the sparseness structure into account).
"""
mutable struct SimulationStatistics
    structureOfDAE::Any
    cpuTimeInitialization::Float64
    cpuTimeIntegration::Float64
    startTime::Float64
    stopTime::Float64
    interval::Float64
    tolerance::Float64
    nEquations::Int
    nConstraints::Union{Int,Missing}
    nResults::Int
    nSteps::Int
    nResidues::Int
    nZeroCrossings::Int
    nJac::Int
    nTimeEvents::Int
    nStateEvents::Int
    nRestartEvents::Int
    nErrTestFails::Int
    h0::Float64
    hMin::Float64
    hMax::Float64
    orderMax::Int
    sparseSolver::Bool
    nGroups::Int

    SimulationStatistics(structureOfDAE, nEquations::Int, nConstraints::Union{Int,Missing}, sparseSolver::Bool, nGroups::Int) =
    new(structureOfDAE, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, nEquations, nConstraints, 0, 0, 0, 0, 0, 0, 0, 0, 0, floatmax(Float64),
        floatmax(Float64), 0.0, 0, sparseSolver, nGroups)
end

end


text:
 module Meshsub

export meshsub

using Printf: @sprintf
using Triangulate: TriangulateIO, triangulate
using ..Sheets: SV2, RWGSheet

"""
    meshsub(;points, seglist, area,
            segmarkers=Int[], holes=Array{Float64}(undef,2,0),  switches="") -> sheet::RWGSheet

Triangulate polygonal region(s) and create a sheet object with the triangulation fields properly initialized.

#  Required Input Keyword Arguments:

- `points::Matrix{<:Real}` with size(points,1) == 2. Input list of boundary segment vertex coordinates.
- `seglist::Matrix{<:Integer}` with size(seglist,1) == 2. List of boundary segments for the
               regions(s) to be triangulated.  seglist[:,n] contains the point
               indices for the n'th boundary segment.
- `area::Real` Desired maximum area of any triangle in the generated mesh.
- `ntri::Int`  Desired number of triangles in the generated mesh.

#  Optional Input Keyword Arguments:

- `segmarkers::Vector{<:Integer}` with length(segmarkers) == size(seglist,2). Each boundary
               polygon is numbered starting with 1. These boundary numbers are entered
               in segmarkers to denote which boundary a particular segment belongs to.
               Required if more than a single boundary is specified (i.e. for multiply connected regions).
- `holes::Matrix{<:Real} with size(holes,1) == 2.  Each column `n` contains the coordinates
               of a point within the n'th *hole* (i.e., region which is not to be triangulated). There
               must be a boundary polygon defined for the boundary of each hole.  Required if more
               than a single boundary is specified (i.e., for multiply connected regions).
- `switches::String` A set of switches to pass to `Triangulation.triangulate` instead of the default
               `"pDa\$(area)q20.0Qe"`.  If you pass your own switches, be sure to include `"e"`.

#  Return value:
    sheet      A variable of type RWGSheet with fields ρ, e1, e2, fe, and fv properly initialized.

"""
function meshsub(;points::Matrix{<:Real}, seglist::Matrix{<:Integer},
                 segmarkers=Int[],
                 holes=Array{Float64}(undef,2,0),
                 area::Real=0.0,
                 ntri::Int,
                 switches::String="")::RWGSheet

    size(points,1) == 2 || throw(ArgumentError("points first dimension must be 2"))
    size(seglist,1) == 2 || throw(ArgumentError("seglist first dimension must be 2"))
    isempty(segmarkers) || (length(segmarkers) == size(seglist,2)) ||
        throw(ArgumentError("wrong length(segmarkers)"))
    isempty(holes) || (size(holes,1) == 2) ||
        throw(ArgumentError("holes first dimension must be 2"))
    minangstr = "30.0" # Angle quality requirement
    if isempty(switches)
        iter = 14 # Number of triangulation iterations to reach desired ntri
        astr = @sprintf("%.14f", area)
        switches = "pDa$(astr)q$(minangstr)Qe"
    else
        iter = 1
    end

    # Set up for call to triangulate:
    triin = TriangulateIO()
    triin.pointlist=Matrix{Cdouble}(points)
    triin.segmentlist=Matrix{Cint}(seglist)
    isempty(segmarkers) || (triin.segmentmarkerlist = Vector{Cint}(segmarkers))
    isempty(holes) || (triin.holelist = Matrix{Cdouble}(holes))

    areaold = area
    triout = deepcopy(triin) # Establish scope
    for k in 1:iter
        # Perform triangulation:
        (triout, vorout)=triangulate(switches, triin)
        ntrinew = size(triout.trianglelist,2)
        (0.9ntri ≤ ntrinew ≤ 1.1ntri) && break
        correction = (ntrinew/ntri)^(1/3)
        areanew = areaold * correction
        astr = @sprintf("%.14f", areanew)
        switches = "pDa$(astr)q$(minangstr)Qe"
        areaold = areanew
    end

    # Set up the sheet data structures:
    sh = RWGSheet()
    plist = triout.pointlist
    sh.ρ = [SV2(plist[:,i]) for i in 1:size(plist, 2)]
    sh.e1 = triout.edgelist[1,:]
    sh.e2 = triout.edgelist[2,:]
    sh.fv = copy(triout.trianglelist)
    # Set up sh.fe so that for any triangle, edge 1 connects nodes 1 and 3,
    # edge 2 connects nodes 1 and 2, and edge 3 connects nodes 2 and 3:
    sh.fe = similar(sh.fv) # Preallocation
    for tri in 1:size(sh.fv,2)
        n1, n2, n3 = sh.fv[:,tri]
        for (side, m1, m2) in [(1,n2,n3), (2,n1,n3), (3,n1,n2)]
                found = false
            for i in 1:length(sh.e1)
                if (sh.e1[i] == m1 && sh.e2[i] == m2) || (sh.e1[i] == m2 && sh.e2[i] == m1)
                    sh.fe[side,tri] = i
                    found = true
                    break
                end
            end
            !found && error("triangle edge not found connecting vertices $m1 and $m2 for face $tri")
        end
    end
    sh.fr = zeros(Float64, size(sh.fv,2)) # Face resistance vector
    return sh
end # function

end # module


text:
"""
    get_authors_titles(dir::String) -> Tuple{Vector{String}, Vector{String}}

Extract the authors and titles from all PDFs found recursively in the passed directory.

# Arguments
- `dir::String`: $(DIR_ARGUMENT)

# Returns
- `Tuple{Vector{String}, Vector{String}}`: the authors and titles

# Throws
- [`DirectoryDoesNotExist`](@ref): $(DIRECTORY_DOES_NOT_EXIST_EXCEPTION)
- Exceptions from: [`get_author_title`](@ref)

# Example
```jldoctest; output = false
using PDFHighlights

path_to_pdf_dir = joinpath(pathof(PDFHighlights) |> dirname |> dirname, "test", "pdf")

get_authors_titles(path_to_pdf_dir) == (["Pavel Sobolev"], ["A dummy PDF for tests"])

# output

true
```
"""
function get_authors_titles(dir::String)::Tuple{Vector{String}, Vector{String}}

    authors = String[]
    titles = String[]

    if isdir(dir)
        for (root, dirs, files) in walkdir(dir), file in files
            if endswith(file, ".pdf")
                author, title = get_author_title(joinpath(root, file))
                push!(authors, author)
                push!(titles, title)
            end
        end
    else
        throw(DirectoryDoesNotExist(dir))
    end

    return authors, titles

end


text:
export tap_on_subscribe, TapBeforeSubscription, TapAfterSubscription

import Base: show

"""
    TapBeforeSubscription

One of the strategies for `tap_on_subscribe` operator. With `TapBeforeSubscription` tap callback will be called before actual subscription.

See also: [`tap_on_subscribe`](@ref), [`TapAfterSubscription`](@ref)
"""
struct TapBeforeSubscription end

"""
    TapAfterSubscription

One of the strategies for `tap_on_subscribe` operator. With `TapBeforeSubscription` tap callback will be called after actual subscription.

See also: [`tap_on_subscribe`](@ref), [`TapBeforeSubscription`](@ref)
"""
struct TapAfterSubscription end

"""
    tap_on_subscribe(tapFn::F, strategy::S = TapBeforeSubscription()) where { F <: Function }

Creates a tap operator, which performs a side effect on the subscription on the source Observable, but return an Observable that is identical to the source.

# Arguments
- `tapFn::Function`: side-effect tap function with `() -> Nothing` signature
- `strategy`: (optional), specifies the order of a side-effect and an actual subscription, uses `TapBeforeSubscription` by default

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source = from([ 1, 2, 3 ])
subscribe!(source |> tap_on_subscribe(() -> println("Someone subscribed")), logger())
;

# output

Someone subscribed
[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Data: 3
[LogActor] Completed
```

```jldoctest
using Rocket

source = from([ 1, 2, 3 ])
subscribe!(source |> tap_on_subscribe(() -> println("Someone subscribed"), TapAfterSubscription()), logger())
;

# output

[LogActor] Data: 1
[LogActor] Data: 2
[LogActor] Data: 3
[LogActor] Completed
Someone subscribed
```

See also: [`TapBeforeSubscription`](@ref), [`TapAfterSubscription`](@ref), [`tap`](@ref), [`tap_on_unsubscribe`](@ref), [`tap_on_complete`](@ref), [`logger`](@ref)
"""
tap_on_subscribe(tapFn::F, strategy::S = TapBeforeSubscription()) where { F <: Function, S } = TapOnSubscribeOperator{F, S}(tapFn, strategy)

struct TapOnSubscribeOperator{F, S} <: InferableOperator
    tapFn    :: F
    strategy :: S
end

operator_right(::TapOnSubscribeOperator, ::Type{L}) where L = L

function on_call!(::Type{L}, ::Type{L}, operator::TapOnSubscribeOperator{F, S}, source) where { L, F, S }
    return proxy(L, source, TapOnSubscribeProxy{F, S}(operator.tapFn, operator.strategy))
end

struct TapOnSubscribeProxy{F, S} <: SourceProxy
    tapFn    :: F
    strategy :: S
end

source_proxy!(::Type{L}, proxy::TapOnSubscribeProxy{F, T}, source::S) where { L, S, F, T } = TapOnSubscribeSource{L, S, F, T}(proxy.tapFn, proxy.strategy, source)

@subscribable struct TapOnSubscribeSource{L, S, F, T} <: Subscribable{L}
    tapFn    :: F
    strategy :: T
    source   :: S
end

on_subscribe!(source::TapOnSubscribeSource, actor) = __on_subscribe_with_tap(source.strategy, source, actor)

function __on_subscribe_with_tap(::TapBeforeSubscription, source::TapOnSubscribeSource, actor)
    source.tapFn()
    return subscribe!(source.source, actor)
end

function __on_subscribe_with_tap(::TapAfterSubscription, source::TapOnSubscribeSource, actor)
    subscription = subscribe!(source.source, actor)
    source.tapFn()
    return subscription
end

Base.show(io::IO, ::TapOnSubscribeOperator)          = print(io, "TapOnSubscribeOperator()")
Base.show(io::IO, ::TapOnSubscribeProxy)             = print(io, "TapOnSubscribeProxy()")
Base.show(io::IO, ::TapOnSubscribeSource{L}) where L = print(io, "TapOnSubscribeSource($L)")


text:
# include all different thermo files here
#TODO Automate it to include all jl files in the Thermo folder
include("Thermo/IdealGasThermo.jl")
# Metaprogramming loop to generate and export all mass and mean functions
property_names=((:cv,"Heat capacity at constant volume (cv)"),
                (:cp,"Heat capacity at constant pressure (cp)"),
                (:s,"entropy (s)"),
                (:s0,"reference entropy (s0)"),
                (:h,"enthalpy (h)"),
                (:a,"helmholz free energy (a)"),
                (:g,"gibbs free energy (g)"),
                (:u,"internal energy"))
for (phi, doc_name) in property_names
    # decide if dimless factor is R or RT
    if phi in (:cv, :cp, :s, :s0)
        dimmless = Symbol(:_R)
        RRT =:(R)
    else
        dimmless = Symbol(:_RT)
        RRT =:(R*T)
    end
    #define the 5 functions for each quantity
    cal_phi_dimless = Symbol(:cal_,phi,dimmless)
    cal_phi = Symbol(:cal_,phi)
    cal_phimass =  Symbol(:cal_,phi,:mass)
    cal_phi_mean = Symbol(cal_phi,:_mean)
    cal_phimass_mean = Symbol(cal_phimass,:_mean)

    @eval begin
        # Dispatches the call from a solution object to it's thermo object
        $cal_phi_dimless(gas::Arrhenius.Solution, T, p, X)=$cal_phi_dimless(gas,gas.thermo, T, p, X)
        export $cal_phi_dimless
        """
            $($cal_phi)(Solution, T, p, X)

        calculates the molar $($doc_name) for each species
        """
        function $cal_phi(gas::Arrhenius.Solution,thermo::Thermo,
                           T::Real, p::Real, X::AbstractArray)
            return $cal_phi_dimless(gas,thermo, T, p, X) * $RRT
        end
        $cal_phi(gas::Arrhenius.Solution, T, p, X)=$cal_phi(gas,gas.thermo, T, p, X)
        export $cal_phi
        """
            $($cal_phi_mean)(Solution, T, p, X)

        calculates the mean mole based $($doc_name) of the mixture
        """
        function $cal_phi_mean(gas::Arrhenius.Solution,thermo::Thermo,
                                T::Real, p::Real, X::AbstractArray)
            return dot(X,$cal_phi_dimless(gas,thermo, T, p, X)) * $RRT
        end
        $cal_phi_mean(gas::Arrhenius.Solution, T, p, X)=$cal_phi_mean(gas,gas.thermo, T, p, X)
        export $cal_phi_mean
        """
            $($cal_phimass)(Solution, T, p, X)

        calculates the partial mass based $($doc_name) for each species
        """
        function $cal_phimass(gas::Arrhenius.Solution,thermo::Thermo,
                              T::Real, p::Real, X::AbstractArray)
            return $cal_phi(gas,thermo, T, p, X) ./gas.MW
        end
        $cal_phimass(gas::Arrhenius.Solution, T, p, X)=$cal_phimass(gas,gas.thermo, T, p, X)
        export $cal_phimass
        """
            $($cal_phimass_mean)(Solution, T, p, X)

        calculates the mean mass based $($doc_name) of the mixture
        """
        function $cal_phimass_mean(gas::Arrhenius.Solution,thermo::Thermo,
                                    T::Real, p::Real, X::AbstractArray)
            return $cal_phi_mean(gas,thermo, T, p, X) / dot(X,gas.MW)
        end
        $cal_phimass_mean(gas::Arrhenius.Solution, T, p, X)=$cal_phimass_mean(gas,gas.thermo, T, p, X)
        export $cal_phimass_mean
    end
end

# ------ Deprecated ----------#

"get specific of heat capacity"
function get_cp(gas, T, X, mean_MW)
    cp_T = [1.0, T, T^2, T^3, T^4]
    if T <= 1000.0
        cp = @view(gas.thermo.nasa_low[:, 1:5]) * cp_T
    else
        cp = @view(gas.thermo.nasa_high[:, 1:5]) * cp_T
    end
    # TODO: not sure if inplace operation will be an issue for AD
    if !gas.thermo.isTcommon
        ind_correction = @. (T > 1000.0) & (T < gas.thermo.Trange[:, 2])
        cp[ind_correction] .=
            @view(gas.thermo.nasa_low[ind_correction, 1:5]) * cp_T
    end
    cp_mole = dot(cp, X) * R
    cp_mass = cp_mole / mean_MW
    return cp_mole, cp_mass
end
export get_cp

"get specific of heat capacity"
function get_cv(cp_mole, cp_mass, mean_MW)
    cv_mole = cp_mole - R
    cv_mass = cv_mole / mean_MW
    return cv_mole, cv_mass
end
export get_cv

"get enthaphy (H) per mole"
function get_H(gas, T, Y, X)
    H_T = [1.0, T / 2.0, T^2 / 3.0, T^3 / 4.0, T^4 / 5.0, 1.0 / T]
    if T <= 1000.0
        h_mole = @view(gas.thermo.nasa_low[:, 1:6]) * H_T * R * T
    else
        h_mole = @view(gas.thermo.nasa_high[:, 1:6]) * H_T * R * T
    end
    if !gas.thermo.isTcommon
        ind_correction = @. (T > 1000.0) & (T < gas.thermo.Trange[:, 2])
        h_mole[ind_correction] .=
            @view(gas.thermo.nasa_low[ind_correction, 1:6]) * H_T * R * T
    end
    # H_mole = dot(h_mole, X)
    return h_mole
end
export get_H

"get enthaphy (H) per mass"
function H_mass_func(gas, h_mole, Y)
    return dot(h_mole ./ gas.MW, Y)
end
export H_mass_func

"get enthaphy (U) per mole"
function get_U(h_mole, T)
    u_mole = h_mole .- (R * T)
    return u_mole
end
export get_U

"get enthaphy (U) per mass"
function U_mass_func(gas, u_mole, Y)
    return dot(u_mole ./ gas.MW, Y)
end
export U_mass_func

"get entropy (S)"
function get_S(gas, T, P, X)
    S_T = [log(T), T, T^2 / 2.0, T^3 / 3.0, T^4 / 4.0, 1.0]
    if T <= 1000.0
        S0 = @view(gas.thermo.nasa_low[:, [1, 2, 3, 4, 5, 7]]) * S_T * R
    else
        S0 = @view(gas.thermo.nasa_high[:, [1, 2, 3, 4, 5, 7]]) * S_T * R
    end
    if !gas.thermo.isTcommon
        ind_correction = @. (T > 1000.0) & (T < gas.thermo.Trange[:, 2])
        S0[ind_correction] .=
            @view(gas.thermo.nasa_low[ind_correction, [1, 2, 3, 4, 5, 7]]) *
            S_T * R
    end
    # _X = @. S0 - R * log(clamp(X, 1.e-30, Inf))
    # s_mole = _X .- R * (P / one_atm)
    # S_mole = dot(s_mole, X)
    return S0
end
export get_S

"get entropy (S) per unit mass"
function S_mass_func(gas, s_mole, Y)
    return dot(s_mole ./ gas.MW, Y)
end
export S_mass_func


text:
module SpringCollab2020UnderwaterSwitchController

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/UnderwaterSwitchController" UnderwaterSwitchController(x::Integer, y::Integer, flag::String="underwater_switch")

const placements = Ahorn.PlacementDict(
    "Underwater Switch Controller (Spring Collab 2020)" => Ahorn.EntityPlacement(
        UnderwaterSwitchController
    )
)

function Ahorn.selection(entity::UnderwaterSwitchController)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::UnderwaterSwitchController, room::Maple.Room) = Ahorn.drawSprite(ctx, "ahorn/SpringCollab2020/underwater_switch", 0, 0)

end


text:
"""
    DynamicsOutputs

Infrastructure for saving quantities during trajectories.

Defines all available options that can be used within the `output` tuple along with
the functions that perform the saving operation.
"""
module DynamicsOutputs

using RecursiveArrayTools: ArrayPartition
using DiffEqCallbacks: SavedValues, SavingCallback
using ComponentArrays: ComponentVector
using RingPolymerArrays: get_centroid

using NQCModels: NQCModels
using NQCDynamics:
    Estimators,
    DynamicsUtils,
    AbstractSimulation

using ..DynamicsUtils:
    get_positions,
    get_velocities,
    get_quantum_subsystem

struct EnsembleSaver{N,F,S<:AbstractSimulation}
    function_names::NTuple{N, Symbol}
    output_functions::F
    sim::S
end

function EnsembleSaver(function_names::NTuple{N, Symbol}, sim) where {N}
    output_functions = tuple([getfield(DynamicsOutputs, f) for f in function_names]...)
    EnsembleSaver(function_names, output_functions, sim)
end

function output_template(output::EnsembleSaver, savepoints, u0)
    sol = (t=savepoints, u=[u0 for _ in savepoints])
    out = output(sol, 0)[1]
    for i in eachindex(out)
        for j in eachindex(out[i])
            out[i][j] = zero(out[i][j])
        end
    end
    return out
end

function (output::EnsembleSaver)(sol, _)
    out = [Any[] for _ in 1:length(sol.t)]
    for (i, (t, u)) in enumerate(zip(sol.t, sol.u))
        sizehint!(out[i], length(output.function_names)+1)
        push!(out[i], t)
        evaluate_output!(out[i], u, t, output.sim, output.output_functions...)
    end
    return (out, false)
end

"""
    create_saving_callback(quantities::NTuple{N, Symbol}; saveat=[]) where {N}

Get the `SavingCallback` that will populate `saved_values` with the result obtained
by evaluating the functions provided in `function_names`.
"""
function create_saving_callback(function_names::NTuple{N, Symbol}; saveat=[]) where {N}
    saved_values = SavedValues(Float64, Vector{Any})
    saving_function = OutputSaver(function_names)
    SavingCallback(saving_function, saved_values; saveat=saveat), saved_values
end

struct OutputSaver{N,F}
    function_names::NTuple{N, Symbol}
    output_functions::F
end

"""
    OutputSaver(function_names::NTuple{N, Symbol}) where {N}

Used to obtain the outputs for all functions given in `function_names`.
"""
function OutputSaver(function_names::NTuple{N, Symbol}) where {N}
    output_functions = tuple([getfield(DynamicsOutputs, f) for f in function_names]...)
    OutputSaver(function_names, output_functions)
end

"""
    (output::OutputSaver)(u, t, integrator)

Evaluates every function listed in `output.function_names` and returns all the results.
"""
function (output::OutputSaver)(u, t, integrator)
    out = Any[]
    sizehint!(out, length(output.function_names))
    evaluate_output!(out, u, t, integrator.p, output.output_functions...)
    return out
end

"""
Used to recursively evaluate every function for the [`OutputSaver`](@ref).

See here for a description of why it is written like this:
https://stackoverflow.com/questions/55840333/type-stability-for-lists-of-closures
"""
function evaluate_output!(out, u, t, sim, f::F, output_functions...) where {F}
    push!(out, f(u, t, sim))
    return evaluate_output!(out, u, t, sim, output_functions...)
end
evaluate_output!(out, u, t, sim) = out

export force
export velocity
export position
export centroid_position
export potential
export hamiltonian
export kinetic
export u
export quantum_subsystem
export state
export noise
export population
export adiabatic_population
export friction

"Get the force"
force(u, t, sim) = -copy(sim.calculator.derivative)

"Get the velocity"
velocity(u, t, sim) = copy(get_velocities(u))

"Get the velocity of the ring polymer centroid"
centroid_velocity(u, t, sim) = get_centroid(get_velocities(u))

"Get the position"
position(u, t, sim) = copy(get_positions(u))

"Get the position of the ring polymer centroid"
centroid_position(u, t, sim) = get_centroid(get_positions(u))

"Evaluate the potential from the model"
potential(u, t, sim) = DynamicsUtils.classical_potential_energy(sim, u)

"Evaluate the classical Hamiltonian"
hamiltonian(u, t, sim) = DynamicsUtils.classical_hamiltonian(sim, u)

"Evaluate the classical kinetic energy"
kinetic(u, t, sim) = DynamicsUtils.classical_kinetic_energy(sim, get_velocities(u))

"Get all the dynamics variables. This is the default"
u(u, t, sim) = copy(u)
u(u::ArrayPartition, t, sim) = ComponentVector(v=copy(get_velocities(u)), r=copy(get_positions(u)))

"""
Get the quantum subsystem of the dynamics variables.
Requires that `DynamicsUtils.get_quantum_subsystem` is implemented for your chosen method.
"""
quantum_subsystem(u, t, sim) = copy(get_quantum_subsystem(u))

mapping_position(u, t, sim) = copy(DynamicsUtils.get_mapping_positions(u))
mapping_momentum(u, t, sim) = copy(DynamicsUtils.get_mapping_momenta(u))

"""
Get the currently occupied state from the dynamics variables.
Requires that the dynamics variable has a field `state`.
Currently this is for surface hopping methods only.
Use [`population`](@ref) or [`adiabatic_population`](@ref) for most other methods.
"""
state(u, t, sim) = copy(u.state)

"Evaluate the diabatic population"
population(u, t, sim) = Estimators.diabatic_population(sim, u)
total_population(u, t, sim) = sum(Estimators.diabatic_population(sim, u))

"Evaluate the adiabatic population"
adiabatic_population(u, t, sim) = Estimators.adiabatic_population(sim, u)

end # module


text:
using Test

using YAXArrays.YAXTools: PickAxisArray
@testset "PickAxisArray" begin
    p = reshape(1:24, 2, 3, 4)
    a = PickAxisArray(p, [:, true, false, true])
    @test a[1, 5, 1] == p[:, 1, 1]
    @test a[1, 6, 2] == p[:, 1, 2]
    @test a[3, 7, 4] == p[:, 3, 4]

    @test @view(a[1, 5, 1]) == p[:, 1, 1]
    @test @view(a[1, 6, 2]) == p[:, 1, 2]
    @test @view(a[3, 7, 4]) == p[:, 3, 4]

    @test eltype(a) <: Vector{Int}

    a = PickAxisArray(p, [true, false, :, :])
    @test a[1] == p[1, :, :]
    @test view(a, 2) == p[2, :, :]

    @test eltype(a) <: Array{Int,2}

    using OffsetArrays
    p = OffsetArray(reshape(1:48, 4, 3, 4), 0:3, 1:3, 1:4)
    a = PickAxisArray(p, [(1, 1), false, :, :])

    @test a[1, 5] == p[0:2, :, :]
    @test a[2, 6] == p[1:3, :, :]

    @test view(a, 1, 4) == p[0:2, :, :]
    @test view(a, 2, 5) == p[1:3, :, :]

    @test eltype(a) <: Array{Int,3}

    p = OffsetArray(reshape(1:60, 3, 4, 5), 1:3, 0:3, 0:4)
    a = PickAxisArray(p, [false, :, (1, 1), (1, 1)])
    @test a[5, 1, 1] == p[:, 0:2, 0:2]
    @test a[6, 1, 2] == p[:, 0:2, 1:3]
    @test a[6, 2, 3] == p[:, 1:3, 2:4]

    @test view(a, 5, 1, 1) == p[:, 0:2, 0:2]
    @test view(a, 6, 1, 2) == p[:, 0:2, 1:3]
    @test view(a, 6, 2, 3) == p[:, 1:3, 2:4]

    @test eltype(a) <: Array{Int,3}

    p = OffsetArray(reshape(1:60, 3, 4, 5), 1:3, 0:3, 0:4)
    a = PickAxisArray(p, [false, :, (1, 1), (1, 1)], [2, 3, 1])
    @test a[5, 1, 1] == permutedims(p[:, 0:2, 0:2], [2, 3, 1])
    @test a[6, 1, 2] == permutedims(p[:, 0:2, 1:3], [2, 3, 1])
    @test a[6, 2, 3] == permutedims(p[:, 1:3, 2:4], [2, 3, 1])

    @test view(a, 5, 1, 1) == permutedims(p[:, 0:2, 0:2], [2, 3, 1])
    @test view(a, 6, 1, 2) == permutedims(p[:, 0:2, 1:3], [2, 3, 1])
    @test view(a, 6, 2, 3) == permutedims(p[:, 1:3, 2:4], [2, 3, 1])

    @test eltype(a) <: Array{Int,3}
end

max_forks_repo_name: meggart/YAXArrays


text:
module MCIntegration
include("utility/utility.jl")

include("montecarlo.jl")
export montecarlo, Configuration, Diagram, FermiK, BoseK, Tau, TauPair

end


text:
"""
    PoissonLikelihood()

Poisson likelihood with rate as exponential of samples from GP `f`. This is to be used if
we assume that the uncertainity associated with the data follows a Poisson distribution.
"""
struct PoissonLikelihood end

(l::PoissonLikelihood)(f::Real) = Poisson(exp(f))

(l::PoissonLikelihood)(fs::AbstractVector{<:Real}) = Product(Poisson.(exp.(fs)))