text:
"""
Γ^{4}，有效势
"""

module FlowEquation

using ..Basics
using ..Triangulated
using ..Refine
using ..Fermi
using ..Fermi.Surface
using ..Fermi.Patch

export Gamma4
export ECGamma4
export pi_αβ_plus_ec, pi_αβ_minus_ec
export all_bubble_ec_mt

export dl_ec_mt

export TFGamma4
export pi_αβ_plus_tf, pi_αβ_minus_tf
export all_bubble_tf_mt

export dl_tf_mt

export pi_αβ_minus_tf_ult, pi_αβ_plus_tf_ult
export all_bubble_tf_ult_mt

export dl_tf_mix_ult_mt

export fliter_away_surface
export refine_to_surface
export engpeak_to_surface
export refine_list_triangle

export TFGamma4_refine_ltris_mt
export TFGamma4_addition_ltris_mt

"""
最后一个动量可以由动量守恒确定
"""
struct Gamma4{T <: Fermi.Abstract2DModel, P <: Basics.AbstractTriangle}
    model :: T
    λ_0 :: Float64
    V :: Array{Float64, 7}
    k4tab :: Array{Int64, 7}
    #patch允许不同能带取不同的位置，但是为了方便必须有一样的patch数目
    #对相同的布里渊区，patch切分的方法相同，所以相同的patchnum和ltris意味着相同的lpats
    patchnum :: Int64
    patches :: Vector{Vector{Basics.Point2D}}
    ltris :: Vector{P}
    lpats :: Vector{Int64}
    ladjs :: Union{Nothing, TYPE_LADJS{P}}
    #Vector{Tuple{Union{Missing, P}, Union{Missing, P}, Union{Missing, P}}}
    ltris_pat :: Union{Nothing, Vector{Vector{P}}}
end

"""
能量截断的Γ4
"""
function ECGamma4(
    model::Union{QuadrateSystem, TriangularSystem},
    λ_0::Float64, patchnum::Int64, splitnum::Int64
    )
    mpats::Vector{Vector{Basics.Point2D}} = []
    for midx in 1:1:model.bandnum
        push!(mpats, patches_under_vonhove(
            model.brillouin, model.dispersion[midx], patchnum
        ))
    end
    #
    V = zeros(
        model.bandnum, model.bandnum, model.bandnum, model.bandnum,
        patchnum, patchnum, patchnum
    )
    #
    k4tab = Array{Int64, 7}(undef,
    model.bandnum, model.bandnum, model.bandnum, model.bandnum,
    patchnum, patchnum, patchnum
    )
    #
    find_algo = isa(model, QuadrateSystem) ? find_patch_index_squa :
    find_patch_index_hexa
    for idxs in CartesianIndices(k4tab)
        m1, m2, m3, m4, k1, k2, k3 = Tuple(idxs)
        k1p = mpats[m1][k1]
        k2p = mpats[m2][k2]
        k3p = mpats[m3][k3]
        k4p = model.kadd(k1p, k2p)
        k4p = model.kadd(k4p, -k3p)
        k4tab[m1, m2, m3, m4, k1, k2, k3] =
        find_algo(k4p, model.brillouin, patchnum)
    end
    #
    split_algo = isa(model, QuadrateSystem) ? split_square :
    split_hexagon
    ltris, ladjs = split_algo(model.brillouin, splitnum)
    #
    lpats = group_ltris_into_patches_mt(ltris, model.brillouin, patchnum)
    return Gamma4(
        model,
        λ_0,
        V,
        k4tab,
        patchnum,
        mpats,
        ltris, lpats, ladjs, nothing
    )
end

function TFGamma4(
    model::Union{QuadrateSystem, TriangularSystem},
    λ_0::Float64, patchnum::Int64, splitnum::Int64
    )
    mpats::Vector{Vector{Basics.Point2D}} = []
    for midx in 1:1:model.bandnum
        push!(mpats, patches_under_vonhove(
            model.brillouin, model.dispersion[midx], patchnum
        ))
    end
    #
    V = zeros(
        model.bandnum, model.bandnum, model.bandnum, model.bandnum,
        patchnum, patchnum, patchnum
    )
    #
    k4tab = Array{Int64, 7}(undef,
    model.bandnum, model.bandnum, model.bandnum, model.bandnum,
    patchnum, patchnum, patchnum
    )
    #
    find_algo = isa(model, QuadrateSystem) ? find_patch_index_squa :
    find_patch_index_hexa
    for idxs in CartesianIndices(k4tab)
        m1, m2, m3, m4, k1, k2, k3 = Tuple(idxs)
        k1p = mpats[m1][k1]
        k2p = mpats[m2][k2]
        k3p = mpats[m3][k3]
        k4p = model.kadd(k1p, k2p)
        k4p = model.kadd(k4p, -k3p)
        k4tab[m1, m2, m3, m4, k1, k2, k3] =
        find_algo(k4p, model.brillouin, patchnum)
    end
    #
    split_algo = isa(model, QuadrateSystem) ? split_square :
    split_hexagon
    ltris, ladjs = split_algo(model.brillouin, splitnum)
    lpats = group_ltris_into_patches_mt(ltris, model.brillouin, patchnum)
    #
    #因为现在所有的能带都必须有同一个lpats，所以只要每一个patch中有哪些tri
    ltris_pat = Vector{typeof(ltris)}(undef, patchnum)
    for idx in 1:1:patchnum
        ltris_pat[idx] = []
    end
    #
    for (tri, pat) in zip(ltris, lpats)
        push!(ltris_pat[pat], tri)
    end
    #
    return Gamma4(
        model,
        λ_0,
        V,
        k4tab,
        patchnum,
        mpats,
        ltris, lpats, nothing, ltris_pat
    )
end

include("bubble.jl")

include("ec_bubble.jl")

"""
能量截断的导数
"""
function dl_ec_mt(Γ4::Gamma4, bubb_pp::T1, bubb_fs::T2,
    bubb_nfs::T3, bubb_ex::T4, bubb_nex::T5) where {
        T1<:Bubble, T2<:Bubble,
        T3<:Bubble, T4<:Bubble, T5<:Bubble}
    #
    sys = Γ4.model
    dl_val = zeros(size(Γ4.V))
    #所有的自由度
    Threads.@threads for idxs in CartesianIndices(dl_val)
        b1, b2, b3, b4, i1, i2, i3 = Tuple(idxs)
        i4 = Γ4.k4tab[b1, b2, b3, b4, i1, i2, i3]
        #需要进行的积分
        value = 0.
        place_holder = Array{Int8, 3}(undef, sys.bandnum, sys.bandnum, Γ4.patchnum)
        for intidxs in CartesianIndices(place_holder)
            α, β, i_n = Tuple(intidxs)
            pi_min_αβ_n_qpp = bubb_pp.V[α, β, b1, b2, i_n, i1, i2]
            pi_plu_αβ_n_qfs = bubb_fs.V[α, β, b2, b3, i_n, i2, i3]
            pi_plu_αβ_n_nqfs = bubb_nfs.V[α, β, b2, b3, i_n, i2, i3]
            pi_plu_αβ_n_qex = bubb_ex.V[α, β, b1, b3, i_n, i1, i3]
            pi_plu_αβ_n_nqex = bubb_nex.V[α, β, b1, b3, i_n, i1, i3]
            #
            value += Γ4.V[b2, b1, α, β, i2, i1, i_n] *
            Γ4.V[b3, b4, α, β, i3, i4, i_n] * pi_min_αβ_n_qpp
            #
            value += Γ4.V[b1, b2, α, β, i1, i2, i_n] *
            Γ4.V[b4, b3, α, β, i4, i3, i_n] * pi_min_αβ_n_qpp
            #
            value += 2*Γ4.V[α, b4, b1, β, i_n, i4, i1] *
            Γ4.V[α, b2, b3, β, i_n, i2, i3] * pi_plu_αβ_n_qfs
            #
            value += 2*Γ4.V[α, b1, b4, β, i_n, i1, i4] *
            Γ4.V[α, b3, b2, β, i_n, i3, i2] * pi_plu_αβ_n_nqfs
            #
            value -= Γ4.V[b4, α, b1, β, i4, i_n, i1] *
            Γ4.V[α, b2, b3, β, i_n, i2, i3] * pi_plu_αβ_n_qfs
            #
            value -= Γ4.V[b1, α, b4, β, i1, i_n, i4] *
            Γ4.V[α, b3, b2, β, i_n, i3, i2] * pi_plu_αβ_n_nqfs
            #
            value -= Γ4.V[α, b4, b1, β, i_n, i4, i1] *
            Γ4.V[b2, α, b3, β, i2, i_n, i3] * pi_plu_αβ_n_qfs
            #
            value -= Γ4.V[α, b1, b4, β, i_n, i1, i4] *
            Γ4.V[b3, α, b2, β, i3, i_n, i2] * pi_plu_αβ_n_nqfs
            #
            value -= Γ4.V[b3, α, b1, β, i3, i_n, i1] *
            Γ4.V[b2, α, b4, β, i2, i_n, i4] * pi_plu_αβ_n_qex
            #
            value -= Γ4.V[b1, α, b3, β, i1, i_n, i3] *
            Γ4.V[b4, α, b2, β, i4, i_n, i2] * pi_plu_αβ_n_nqex
            #
        end#
        ##总的负号
        value = -value
        ##太大的时候就线暂停
        if abs(value) > 1e32
            value = 0
        end
        #计算完成
        dl_val[b1, b2, b3, b4, i1, i2, i3] = value
    end
    return dl_val
end

include("tf_bubble.jl")

"""
温度流的导数
"""
function dl_tf_mt(Γ4::Gamma4, bubb_pp::T1, bubb_fs::T2, bubb_ex::T3
    ) where {T1<:Bubble, T2<:Bubble, T3<:Bubble}
    #
    sys = Γ4.model
    dl_val = zeros(size(Γ4.V))
    #所有的自由度
    Threads.@threads for idxs in CartesianIndices(dl_val)
        b1, b2, b3, b4, i1, i2, i3 = Tuple(idxs)
        i4 = Γ4.k4tab[b1, b2, b3, b4, i1, i2, i3]
        #需要进行的积分
        value = 0.
        place_holder = Array{Int8, 3}(undef, sys.bandnum, sys.bandnum, Γ4.patchnum)
        for intidxs in CartesianIndices(place_holder)
            α, β, i_n = Tuple(intidxs)
            pi_min_αβ_n_qpp = bubb_pp.V[α, β, b1, b2, i_n, i1, i2]
            pi_plu_αβ_n_qfs = bubb_fs.V[α, β, b2, b3, i_n, i2, i3]
            pi_plu_αβ_n_qex = bubb_ex.V[α, β, b1, b3, i_n, i1, i3]
            #
            value += Γ4.V[b2, b1, α, β, i2, i1, i_n] *
                Γ4.V[b3, b4, α, β, i3, i4, i_n] * pi_min_αβ_n_qpp
            value += 2*Γ4.V[α, b4, b1, β, i_n, i4, i1] *
                Γ4.V[α, b2, b3, β, i_n, i2, i3] * pi_plu_αβ_n_qfs
            value -= Γ4.V[b4, α, b1, β, i4, i_n, i1] *
                Γ4.V[α, b2, b3, β, i_n, i2, i3] * pi_plu_αβ_n_qfs
            value -= Γ4.V[α, b4, b1, β, i_n, i4, i1] *
                Γ4.V[b2, α, b3, β, i2, i_n, i3] * pi_plu_αβ_n_qfs
            value -= Γ4.V[b3, α, b1, β, i3, i_n, i1] *
                Γ4.V[b2, α, b4, β, i2, i_n, i4] * pi_plu_αβ_n_qex
        end#结束一个微分数值的计算
        #
        value = -value
        if abs(value) > 1.e32
            value = 0.
        end
        #计算完成
        dl_val[b1, b2, b3, b4, i1, i2, i3] = value
    end# 结束对所有能带和动量的循环
    return dl_val
end

include("tf_bubble_ult.jl")

"""
温度流的导数, 将极低温的贡献和常态温度的混合
"""
function dl_tf_mix_ult_mt(Γ4::Gamma4, bubb_pp_com::T1, bubb_fs_com::T2,
    bubb_ex_com::T2, bubb_pp_ult::T1, bubb_fs_ult::T2, bubb_ex_ult::T2
    ) where {T1<:Bubble, T2<:Bubble}
    #
    sys = Γ4.model
    dl_val = zeros(size(Γ4.V))
    #附加的ult的贡献
    #小三角型的面积
    triarea = area(Γ4.ltris[1])
    #贡献的半高宽
    halfarea = 4 * (Γ4.λ_0*exp(-bubb_pp_com.lval))^2
    coef = bubb_pp_com.lval > 12 ? 1.0 : 0.0
    #max(0., 1 - (halfarea / triarea))#1 / (1 + exp(halfarea/triarea))
    #所有的自由度
    Threads.@threads for idxs in CartesianIndices(dl_val)
        b1, b2, b3, b4, i1, i2, i3 = Tuple(idxs)
        i4 = Γ4.k4tab[b1, b2, b3, b4, i1, i2, i3]
        #需要进行的积分
        value = 0.
        place_holder = Array{Int8, 3}(undef, sys.bandnum, sys.bandnum, Γ4.patchnum)
        for intidxs in CartesianIndices(place_holder)
            α, β, i_n = Tuple(intidxs)
            pi_min_αβ_n_qpp = bubb_pp_com.V[α, β, b1, b2, i_n, i1, i2]
            pi_plu_αβ_n_qfs = bubb_fs_com.V[α, β, b2, b3, i_n, i2, i3]
            pi_plu_αβ_n_qex = bubb_ex_com.V[α, β, b1, b3, i_n, i1, i3]
            #增加上ult的贡献
            pi_min_αβ_n_qpp += coef * bubb_pp_ult.V[α, β, b1, b2, i_n, i1, i2]
            pi_plu_αβ_n_qfs += coef * bubb_fs_ult.V[α, β, b2, b3, i_n, i2, i3]
            pi_plu_αβ_n_qex += coef * bubb_ex_ult.V[α, β, b1, b3, i_n, i1, i3]
            #
            value += Γ4.V[b2, b1, α, β, i2, i1, i_n] *
                Γ4.V[b3, b4, α, β, i3, i4, i_n] * pi_min_αβ_n_qpp
            value += 2*Γ4.V[α, b4, b1, β, i_n, i4, i1] *
                Γ4.V[α, b2, b3, β, i_n, i2, i3] * pi_plu_αβ_n_qfs
            value -= Γ4.V[b4, α, b1, β, i4, i_n, i1] *
                Γ4.V[α, b2, b3, β, i_n, i2, i3] * pi_plu_αβ_n_qfs
            value -= Γ4.V[α, b4, b1, β, i_n, i4, i1] *
                Γ4.V[b2, α, b3, β, i2, i_n, i3] * pi_plu_αβ_n_qfs
            value -= Γ4.V[b3, α, b1, β, i3, i_n, i1] *
                Γ4.V[b2, α, b4, β, i2, i_n, i4] * pi_plu_αβ_n_qex
        end#结束一个微分数值的计算
        #
        value = -value
        if abs(value) > 1.e32
            value = 0.
        end
        #计算完成
        dl_val[b1, b2, b3, b4, i1, i2, i3] = value
    end# 结束对所有能带和动量的循环
    return dl_val
end

"""
取掉离费米面太远的
"""
function fliter_away_surface(Γ4::Gamma4{T, P}, lval) where {T, P}
    lamb = Γ4.λ_0 * exp(-lval)
    resltris::Vector{P} = []
    reslpats::Vector{Int64} = []
    cert = 1#length(Γ4.ltris) > 2000000 ? 5 : 15
    for (tri, pat) in zip(Γ4.ltris, Γ4.lpats)
        engs = [Γ4.model.dispersion[bidx](
            tri.center.x, tri.center.y) / lamb for bidx in 1:1:Γ4.model.bandnum
        ]
        #如果还有可能有贡献
        if minimum(abs.(engs)) < cert #&& minimum(abs.(engs2)) > cert # && pat == 1# && minimum(abs.(engs)) > 0.1
            push!(resltris, tri)
            push!(reslpats, pat)
        end
    end
    ltris_pat = Vector{typeof(resltris)}(undef, Γ4.patchnum)
    for idx in 1:1:Γ4.patchnum
        ltris_pat[idx] = []
    end
    #
    for (tri, pat) in zip(resltris, reslpats)
        push!(ltris_pat[pat], tri)
    end
    #
    return Gamma4(
        Γ4.model,
        Γ4.λ_0,
        Γ4.V,
        Γ4.k4tab,
        Γ4.patchnum,
        Γ4.patches,
        resltris, reslpats, nothing, ltris_pat
    )
end

"""
将Γ4中的ltris进行重新的分割
"""
function TFGamma4_refine_ltris_mt(Γ4::Gamma4{T, P}, lval; maxnum=2500000) where {T, P}
    #
    if !isnothing(Γ4.ladjs)
        @warn "会丢失ladjs的信息"
    end
    if length(Γ4.ltris) >= maxnum
        @warn "已经超过最大数量"
        return fliter_away_surface(Γ4, lval)
    end
    #
    lamb = Γ4.λ_0 * exp(-lval)
    max2suf = Vector{Float64}(undef, Γ4.model.bandnum)
    for bidx in 1:1:length(max2suf)
        disp = Γ4.model.dispersion[bidx]
        suf = const_energy_triangle(Γ4.ltris, 0., disp)
        engs = zeros(length(suf))
        for (idx, tri) in enumerate(suf)
            center = tri.center
            engs[idx] = abs(disp(center.x, center.y)) / lamb
        end
        max2suf[bidx] = maximum(engs)
    end
    #如果已经足够贴近费米面，直接返回原来的
    if maximum(max2suf) < 0.0001
        return fliter_away_surface(Γ4, lval)
    end
    @info string(length(Γ4.ltris))*" -> "*string(4*length(Γ4.ltris))
    #
    #将所有的nltris进行重新的切分，切分一定是切成4个
    spltris = Vector{P}(undef, 4*length(Γ4.ltris))
    for (idx, tri) in enumerate(Γ4.ltris)
        spltris[4*idx-3:4*idx] = split_triangle(tri)
    end
    #去掉离得太远的
    resltris::Vector{P} = []
    cert = length(spltris) > 1000000 ? 5 : 15
    for tri in spltris
        engs = [Γ4.model.dispersion[bidx](
            tri.center.x, tri.center.y) / lamb for bidx in 1:1:Γ4.model.bandnum
        ]
        #如果还有可能有贡献
        if minimum(abs.(engs)) < cert
            push!(resltris, tri)
        end
    end
    #重新计算一下所属的patch，可能会发生变化
    reslpats = group_ltris_into_patches_mt(resltris, Γ4.model.brillouin, Γ4.patchnum)
    #重新计算一下每个patch包含的三角形
    ltris_pat = Vector{typeof(resltris)}(undef, Γ4.patchnum)
    for idx in 1:1:Γ4.patchnum
        ltris_pat[idx] = []
    end
    #
    for (tri, pat) in zip(resltris, reslpats)
        push!(ltris_pat[pat], tri)
    end
    #
    return Gamma4(
        Γ4.model,
        Γ4.λ_0,
        Γ4.V,
        Γ4.k4tab,
        Γ4.patchnum,
        Γ4.patches,
        resltris, reslpats, nothing, ltris_pat
    )
end

"""
将Γ4中的ltris，靠近费米面的切割的更加细致
"""
function TFGamma4_reweight_ltris_mt(Γ4::Gamma4{T, P}, lval
    ; maxnum=5000000) where {T, P}
    #
    @warn "会被移除"
    truncated = false
    if !isnothing(Γ4.ladjs)
        @warn "会丢失ladjs的信息"
    end
    if length(Γ4.ltris) >= maxnum
        @warn "已经超过最大数量"
        truncated = true
    end
    lamb = Γ4.λ_0 * exp(-lval)
    #切分离费米面近的
    ltris = Γ4.ltris
    newltris::Vector{P} = []
    for tri in ltris
        add_tri = false
        for didx in 1:1:Γ4.model.bandnum
            disp = Γ4.model.dispersion[didx]
            ver1 = tri.vertex[1]
            sgn1 = sign(disp(ver1.x, ver1.y))
            ver2 = tri.vertex[2]
            sgn2 = sign(disp(ver2.x, ver2.y))
            ver3 = tri.vertex[3]
            sgn3 = sign(disp(ver3.x, ver3.y))
            #
            engs = abs(disp(tri.center.x, tri.center.y)) / lamb
            #不在费米面上
            if sgn1 == sgn2 && sgn1 == sgn3
                #如果没有超过最大，或者还有明显贡献
                if !truncated || engs < 3
                    add_tri = true
                end
                #继续循环下一个能带看是否要refine
                continue
            end
            #如果没有超过最大，而且还不够精细
            if !truncated && engs > 0.0001
                ftris = split_triangle(tri)
                push!(newltris, ftris[1])
                push!(newltris, ftris[2])
                push!(newltris, ftris[3])
                push!(newltris, ftris[4])
                #如果在某个能带上已经切掉，就不再计算别的能带了
                add_tri = false
                break
            #费米面上的tri一定会被保留
            else
                add_tri = true
                #继续循环下一个能带看是否要refine
                continue
            end
        end
        if add_tri
            push!(newltris, tri)
        end
    end
    #
    refltris = newltris
    #重新计算一下所属的patch，可能会发生变化
    reflpats = group_ltris_into_patches_mt(refltris, Γ4.model.brillouin, Γ4.patchnum)
    #重新计算一下每个patch包含的三角形
    ltris_pat = Vector{typeof(refltris)}(undef, Γ4.patchnum)
    for idx in 1:1:Γ4.patchnum
        ltris_pat[idx] = []
    end
    #
    for (tri, pat) in zip(refltris, reflpats)
        push!(ltris_pat[pat], tri)
    end
    #
    return Gamma4(
        Γ4.model,
        Γ4.λ_0,
        Γ4.V,
        Γ4.k4tab,
        Γ4.patchnum,
        Γ4.patches,
        refltris, reflpats, nothing, ltris_pat
    )
end

"""
将等能面上面的一串三角形变成更接近的
"""
function refine_to_surface(ltris::Vector{P}, eng, disp) where P
    newltris::Vector{P} = []
    for tri in ltris
        ftris = split_triangle(tri)
        for ftri in ftris
            if @onsurface ftri disp eng
                push!(newltris, ftri)
            end
        end
    end
    return newltris
end

"""
将等能面上面的一串三角形变成更接近的
"""
function refine_to_surface(Γ4::Gamma4{T, P}) where {T, P}
    #
    refltris::Vector{P} = []
    for tri in Γ4.ltris
        ftris = split_triangle(tri)
        for ftri in ftris
            isonsurface = false
            for didx in 1:1:Γ4.model.bandnum
                disp = Γ4.model.dispersion[didx]
                isonsurface = isonsurface || @onsurface ftri disp 0.
            end
            if isonsurface
                push!(refltris, ftri)
            end
        end
    end
    #重新计算一下所属的patch，可能会发生变化
    reflpats = group_ltris_into_patches_mt(refltris, Γ4.model.brillouin, Γ4.patchnum)
    #重新计算一下每个patch包含的三角形
    ltris_pat = Vector{typeof(refltris)}(undef, Γ4.patchnum)
    for idx in 1:1:Γ4.patchnum
        ltris_pat[idx] = []
    end
    #
    for (tri, pat) in zip(refltris, reflpats)
        push!(ltris_pat[pat], tri)
    end
    #
    return Gamma4(
        Γ4.model,
        Γ4.λ_0,
        Γ4.V,
        Γ4.k4tab,
        Γ4.patchnum,
        Γ4.patches,
        refltris, reflpats, nothing, ltris_pat
    )
end

"""
获取现在所有ltris中最大的能量和最小能量
"""
function engpeak_to_surface(Γ4::Gamma4{T, P}) where {T, P}
    maxi = 0.
    mini = 100.
    for tri in Γ4.ltris
        #先找到距离哪个能带最近
        bmini = 100.
        for didx in 1:1:Γ4.model.bandnum
            disp = Γ4.model.dispersion[didx]
            eng = abs(disp(tri.center.x, tri.center.y))
            if eng < bmini
                bmini = eng
            end
        end
        #然后修改最大最小
        if bmini > maxi
            maxi = bmini
        end
        if bmini < mini
            mini = bmini
        end
    end
    return maxi, mini
end

"""
将等能面上面的一串三角形变成更细的
"""
function refine_list_triangle(ltris::Vector{P}, itertime) where P
    nowltris = ltris
    for idx in 1:1:itertime
        newltris = Vector{P}(undef, 4*length(nowltris))
        for (idx, tri) in enumerate(nowltris)
            newltris[4*idx-3:4*idx] = split_triangle(tri)
        end
        nowltris = newltris
    end
    return nowltris
end

function TFGamma4_addition_ltris_mt(Γ4::Gamma4{T, P}, blval
    ; maxnum=5000000) where {T, P}
    @warn "会被移除"
    if !isnothing(Γ4.ladjs)
        @warn "会丢失ladjs的信息"
    end
    if length(Γ4.ltris) >= maxnum
        @warn "已经超过最大数量"
    end
    lamb = Γ4.λ_0 * exp(-blval)
    #找到合适的有贡献的位置
    #println(eng)
    away::Vector{P} = []
    inner::Vector{P} = []
    for tri in Γ4.ltris
        isaway = true
        for bidx in 1:1:Γ4.model.bandnum
            disp = Γ4.model.dispersion[bidx]
            teng = disp(tri.center.x, tri.center.y)
            if abs(teng) <= lamb
                isaway = false
                break
            end
        end
        if isaway
            push!(away, tri)
        else
            push!(inner, tri)
        end
    end
    #将有贡献的位置拆分的更细致
    iterf = log(maxnum - length(away)) - log(length(inner))
    iterf = iterf / log(4)
    iteri = floor(Int64, iterf)
    inner = refine_list_triangle(inner, iteri)
    #
    refltris = vcat(inner, away)
    #重新计算一下所属的patch，可能会发生变化
    reflpats = group_ltris_into_patches_mt(refltris, Γ4.model.brillouin, Γ4.patchnum)
    #重新计算一下每个patch包含的三角形
    ltris_pat = Vector{typeof(refltris)}(undef, Γ4.patchnum)
    for idx in 1:1:Γ4.patchnum
        ltris_pat[idx] = []
    end
    #
    for (tri, pat) in zip(refltris, reflpats)
        push!(ltris_pat[pat], tri)
    end
    #
    return Gamma4(
        Γ4.model,
        Γ4.λ_0,
        Γ4.V,
        Γ4.k4tab,
        Γ4.patchnum,
        Γ4.patches,
        refltris, reflpats, nothing, ltris_pat
    )
end

end # end module


text:
using Flux, DiffEqFlux, OrdinaryDiffEq, DiffEqParamEstim, Plots, Optim, Dates

using DiffEqBase, Flux, DiffResults, DiffEqSensitivity, ForwardDiff
using Flux.Tracker: @grad
using DiffEqSensitivity: adjoint_sensitivities_u0
######################################
mutable struct saver
    losses::Array{Float64,1}
    times::Array{Dates.Time,1}
    count_epochs::Int128
end
function saver(n_epochs)
    losses = zeros(n_epochs)
    times = fill(Dates.Time(Dates.now()),n_epochs)
    count_epochs = 0
    return saver(losses,times,count_epochs)
end
function update_saver(saver, loss_i, time_i)
    epoch_i = saver.count_epochs
    saver.losses[epoch_i] = loss_i
    saver.times[epoch_i] = time_i
end
####################################################### Observation ###############################################################
# Start conditions for the two species in the system
u0 = Float32[2.; 0.]
# Number of evaluations of the neural ODE. It relates to the numbers of layers of the neural net (depth of network).
datasize = 30
# Time span in which of evaluation will be and actual timepoints of evaluations
tspan = (0.0f0, 1.5f0)
t = range(tspan[1], tspan[2], length = datasize)
# The true ODE (with the true parameters) which the neural net should learn
function trueODEfunc(du, u, p, t)
  true_A = [-0.1 2.0; -2.0 -0.1]
  du .= ((u.^3)'true_A)'
end
# Construction of the ODEProblem and solving the ODEProblem with Tsit5 solver
prob = ODEProblem(trueODEfunc, u0, tspan)
ode_data = Array(solve(prob,Tsit5(),saveat=t))
################################################### Building a neural ODE ##########################################################
# Derivative is modeled by a neural net. Chain concatinates the functions ode function and two dense layers.
dudt = Chain(x -> x.^3,
       Dense(2,50,tanh),
       Dense(50,2))
# Parameters of the model which are to be learnt. They are: W1 (2x50), b1 (50), W2 (50x2), b2 (2)
ps = Flux.params(dudt)
#build node
n_ode = x->neural_ode(dudt, x, tspan, Tsit5(), saveat=t, reltol=1e-7, abstol=1e-9)
# to set by user
init_train, number_sets, number_epochs_loss1, number_epochs_loss2 = 8000, 15, 1000, 20
# for saving
sa = saver(init_train+number_sets*(number_epochs_loss1+number_epochs_loss2))
# define collocation loss function with callback
function node_two_stage_function(model, x, tspan, saveat, ode_data,
            args...; kwargs...)
  dudt_(du,u,p,t) = du .= model(u)
  prob_fly = ODEProblem(dudt_,x,tspan)
  two_stage_method(prob_fly, saveat, ode_data)
end
loss_n_ode = node_two_stage_function(dudt, u0, tspan, t, ode_data, Tsit5(), reltol=1e-7, abstol=1e-9)
two_stage_loss_fct()=loss_n_ode.cost_function(ps)
cb1 = function ()
    sa.count_epochs = sa.count_epochs +  1
    update_saver(sa, Tracker.data(two_stage_loss_fct()),Dates.Time(Dates.now()))
    #println("\"",Tracker.data(two_stage_loss_fct()),"\" \"",Dates.Time(Dates.now()),"\";")
end
#define L2 loss function with callback
L2_loss_fct() = sum(abs2,ode_data .- n_ode(u0))
cb2 = function ()
    sa.count_epochs = sa.count_epochs +  1
    update_saver(sa, Tracker.data(L2_loss_fct()),Dates.Time(Dates.now()))
    #println("\"",Tracker.data(L2_loss_fct()),"\" \"",Dates.Time(Dates.now()),"\";")
end
#training call
data0 = Iterators.repeated((), init_train)
data1 = Iterators.repeated((), number_epochs_loss1)
data2 = Iterators.repeated((), number_epochs_loss2)
opt1 = ADAM(0.1)
opt2 = ADAM(0.1)
Flux.train!(two_stage_loss_fct, ps, data0, opt1, cb = cb1)
@time for i in 1:number_sets
    Flux.train!(L2_loss_fct, ps, data2, opt2, cb = cb2)
    Flux.train!(two_stage_loss_fct, ps, data1, opt1, cb = cb1)
end
# Call n_ode to get first prediction and to show startpoint for training.
pred = n_ode(u0)
scatter(t, ode_data[1,:], label="data")
    scatter!(t, Flux.data(pred[1,:]), label="prediction")
    scatter!(t, ode_data[2,:], label="data")
    scatter!(t, Flux.data(pred[2,:]), label="prediction")

#savefig("sogood.png")

using BSON: @save
@save "DiffEqFlux.jl/test/results/mix/model_mix_8000_15_200_20_epochs.bson" dudt

ext: jl


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "libevent"
version = v"2.1.12"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/libevent/libevent/releases/download/release-$(version)-stable/libevent-$(version)-stable.tar.gz",
                  "92e6de1be9ec176428fd2367677e61ceffc2ee1cb119035037a27d346b0403bb"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libevent-*/
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
# FIXME: Name is `libevent-2-1-7.dll` but `parse_dl_name_version` strips the trailing `-7`
products = [
    LibraryProduct(["libevent", "libevent-2-1"], :libevent)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("OpenSSL_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies, julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
include("magpoly.jl")

using .magpoly
const mp = magpoly
using DelimitedFiles, Random, Distributions, LinearAlgebra

spins_coupling = 1.0
#spins_coupling = 0.0
alpha = 0.5
inv_temps = [1.0, 0.9, 0.8, 0.74, 0.65, 0.62, 0.6, 0.55, 0.5, 0.4]
n_temps = length(inv_temps)
n_mono = 300
n_strides = 200
stride = 200

polymers = Array{mp.Magnetic_polymer}(undef, n_temps)
trajs = Array{mp.MC_data}(undef, n_temps)
for i_temp in 1:n_temps
    polymers[i_temp] = mp.Magnetic_polymer(n_mono, inv_temps[i_temp], spins_coupling, alpha)
    trajs[i_temp] = mp.MC_data(n_strides*stride)
    if isfile("simulation_data/final_config_$(n_mono)_$(inv_temps[1]).txt")
        mp.initialize_poly!(polymers[i_temp],"simulation_data/final_config_$(n_mono)_$(inv_temps[1]).txt")
    else
        mp.initialize_poly!(polymers[i_temp])
    end
end

fields = zeros(n_mono)
mp.set_fields!(polymers, fields)
mp.set_coupling!(polymers, 0.7)
println(polymers[1].J)

println("ok")


text:
function randomAssignments(dataset::Dataset; seed::Union{Nothing, Int}=nothing)::Matrix{Int}
    """ Generates random assignments of n nodes to q groups

    Parameters
    ----------
    dataset : Dataset
        Dataset representing an observed graph.

    Returns
    -------
    x : Array{Int64,2}
        Matrix of assignments of nodes to groups
    """
    n = dataset.n
    q = dataset.n_communities
    if seed != nothing
        Random.seed!(seed)
    end

    assign = rand(1:q, n)

    x = zeros(Int, n, q)
    for i=1:n
        x[i,assign[i]] = 1
    end
    @assert sum(x) == n
    return x
end

function findImprovingRelocation(dataset::Dataset, w::Matrix{Float64}, x::Matrix{Int},
                                 reeval_w::Bool; best_imp::Bool=true, time_limit::Float64=400.0)
    """
    Parameters
    ----------
    dataset : Dataset
        Dataset representing an observed graph.
    x : Matrix{Int}
        Matrix of assignments of nodes to groups.
    w : Matrix{Float64}
        Matrix of probabilities of the SBM.
    reeval_w : Bool
        Whether to re-evaluate w to evaluate every possible move.
    best_imp : Bool
        Whether to take the best possible improving move (true) or to take the first improving move (false).
    """
    if best_imp == false
        throw("not implemented")
    end

    start = time(); # Start counting time

    n = dataset.n
    q = dataset.n_communities

    # Calculate objective value
    best_obj = calculateObjective(dataset, w, x)
    best_solution = (copy(w), copy(x))
    improved = false

    # loop through all nodes
    for i=1:n
        # Get current assignment of node i
        current_assign = argmax(x[i,:])
        for g=1:(q-1) # Find the best relocate move for node i

            # Break if time limit is exceeded
            availableTime = time_limit - (time() - start)
            if availableTime <= 0
                (w, x) = best_solution
                return improved, best_obj, w, x
            end

            # Set new assignment for node i
            new_assign = (current_assign + g)
            if new_assign > q
                new_assign = new_assign % q
            end

            # Reassign move
            x[i,current_assign] = 0
            x[i,new_assign] = 1

            # Evaluate move
            if reeval_w ##### Evaluate move on corresponding optimal w
                w = optimalProbMatrix(dataset, x)
            end
            #####
            obj_value = calculateObjective(dataset, w, x)
            if obj_value < best_obj
                best_obj = obj_value
                best_solution = (copy(w), copy(x))
                improved = true
            end

            # Undo reassign move
            x[i,current_assign] = 1
            x[i,new_assign] = 0
        end
    end

    # Reassign using the best move
    ##### Update matrix of probabilities w
    (w, x) = best_solution
    return improved, best_obj, w, x
end

function localSearchAssignments(estimator::SBMEstimator, dataset::Dataset, w::Union{Nothing,Matrix{Float64}}, x::Matrix{Int}; time_limit::Float64=400.0)
    """ Local Search in the space of the assignments x,
    by looking for the best relocate move.
    The value of w (probability matrix) is fixed and
    does not get updated inside this function

    Parameters
    ----------
    estimator : SBMEstimator
        Optimization method specifications.
    dataset : Dataset
        Dataset representing an observed graph.
    w : Union{Nothing,Matrix{Float64}}
        Matrix of probabilities. If nothing, then the optimal w will be calculated for the given x
    x : Matrix{Int}
        Matrix of assignments of nodes to groups.
    time_limit : Float64
        Time limit in seconds.

    Returns
    -------
    x : Matrix{Int}
        Matrix of assignments of nodes to groups (local optimum).
    """

    # Check if a matrix w was given
    if w == nothing
        w = optimalProbMatrix(dataset, x)
    end

    start = time(); # Start counting time

    improved = true
    while improved
        improved = false
        # Break if time limit is exceeded
        availableTime = time_limit - (time() - start)
        if availableTime <= 0
            break
        end
        improved, best_obj, w, x = findImprovingRelocation(dataset, w, x, false, best_imp=true, time_limit=availableTime)
    end
    return x
end


text:
"""
    nc_write(data::AbstractArray{T}, f::AbstractString, dims::Vector{NcDim}, attrib = Dict();
        varname = "x", compress = 1, kwargs..., overwrite = false, mode = "c") where {T<:Real}
    nc_write!(data::AbstractArray{T}, f::AbstractString, dims::Vector{NcDim}, attrib = Dict();
        varname = "x", compress = 1, kwargs...) where {T<:Real}
    nc_write!(data::AbstractArray{T}, f::AbstractString, dims::Vector{<:AbstractString}, attrib = Dict();
        varname = "x", compress = 1, kwargs...) where {T<:Real}

# Arguments

- `type`: which type to save? Julia variable types.

- `kwargs`: Parameters passed to `ncvar_def`, and further to [`NCDatasets.defVar`].
   For examples:

    + fillvalue: A value filled in the NetCDF file to indicate missing data. It
       will be stored in the _FillValue attribute.

    + chunksizes: Vector integers setting the chunk size. The total size of a
       chunk must be less than 4 GiB.

    + shuffle: If true, the shuffle filter is activated which can improve the
       compression ratio.

    + checksum: The checksum method can be :fletcher32 or :nochecksum
       (checksumming is disabled, which is the default)

    + attrib: An iterable of attribute name and attribute value pairs, for
       example a Dict, DataStructures.OrderedDict or simply a vector of pairs
       (see example below)

- `mode`(not used): one of "r", "c", "w", see [NCDatasets.nc_open()] for details

# Examples:
```julia
dims = [
    NcDim("lon", lon, Dict("longname" => "Longitude", "units" => "degrees east"))
    NcDim("lat", lat, Dict("longname" => "Latitude", "units" => "degrees north"))
    NcDim("time", 1:ntime)
]
nc_write(data, f, dims; varname = "x", opt...)

seealso ncvar_def
```
"""
function nc_write(data::AbstractArray{T}, f::AbstractString, dims::Vector{NcDim}, attrib = Dict();
    varname = "x", compress = 1, overwrite = false, mode = "c", kwargs...) where {T<:Real}

    # check whether variable defined
    if !check_file(f) || overwrite
        if isfile(f)
            rm(f)
        end

        ds = nc_open(f, mode)
        ncdim_def(ds, dims)

        dimnames = names(dims)
        ncvar_def(ds, varname, data, dimnames, attrib; compress = compress, kwargs...)
        close(ds)
    else
        println("[file exist]: $(basename(f))")
    end
end

function nc_write!(data::AbstractArray{T}, f::AbstractString, dims::Vector{NcDim}, attrib = Dict();
    varname = "x", compress = 1, kwargs...) where {T<:Real}

    mode = check_file(f) ? "a" : "c";
    ds = nc_open(f, mode)
    ncvar_def(ds, varname, data, dims, attrib; compress = compress, kwargs...)
    close(ds)
end

function nc_write!(data::AbstractArray{T}, f::AbstractString, dims::Vector{<:AbstractString}, attrib = Dict();
    varname = "x", compress = 1, kwargs...) where {T<:Real}

    mode = check_file(f) ? "a" : "c";
    ds = nc_open(f, mode)
    ncvar_def(ds, varname, data, dims, attrib; compress = compress, kwargs...)
    close(ds)
end


text:
using Test
using LinearAlgebra
using OrdinaryDiffEq
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear, prob_ode_lotkavoltera, prob_ode_fitzhughnagumo
using Plots

prob = prob_ode_lotkavoltera
prob = ODEFilters.remake_prob_with_jac(prob)

@testset "Smoothing for small dt and large q" begin
    dt = 1e-4
    q = 4
    @test solve(
        prob, EK0(order=q, smooth=true, diffusionmodel=:dynamic),
        adaptive=false,
        dt=dt,
    ) isa DiffEqBase.AbstractODESolution
end

@testset "Smooth vs. non-smooth" begin
    q = 3
    dt = 1e-2

    sol_nonsmooth = solve(prob, EK0(order=q, smooth=false), dense=false, adaptive=false, dt=dt);
    sol_smooth = solve(prob, EK0(order=q, smooth=true), adaptive=false, dt=dt);

    @test sol_nonsmooth.t ≈ sol_smooth.t
    @test sol_nonsmooth[end] == sol_smooth[end]
    @test sol_nonsmooth[end-1] != sol_smooth[end-1]

    plot(sol_smooth, label="smooth"); plot!(sol_nonsmooth, label="nonsmooth")

    sol_best = solve(prob, Tsit5(), abstol=1e-12, reltol=1e-12);
    sol_best_u = sol_best.(sol_smooth.t);

    nonsmooth_errors = norm.(sol_nonsmooth.u - sol_best_u, 2)
    smooth_errors = norm.(sol_smooth.u - sol_best_u, 2)

    @test 2*maximum(nonsmooth_errors) > maximum(smooth_errors)
    @test 2*sum(nonsmooth_errors) > sum(smooth_errors)

    # Previously we compared the smooth and non-smooth dense output, but this
    # does not work anymore since dense output requires smoothing!
end


text:
"""
   function lagdiff(x::AbstractVector, d=1)

Returns `out[i] = x[i] - x[i-d]` with out padded with missings so that `length(out)==length(x)`.
"""
function lagdiff(x::AbstractVector, d=1)
  if d>0
    dx = vcat(fill(missing, d), x[(d+1):end] .- x[1:(end-d)])
  elseif d<0
    dx = vcat(-x[(d+1):end] .+ x[1:(end-d)], fill(missing,d))
  else
    dx = zero(x)
  end
  return(dx)
end

"""
    function smooth(x::AbstractVector; w=pdf(Normal(), range(-3, 3, length=7)))

Returns the running weighted mean of x*w. When `n=length(w)`, then
`out[i] = sum(x[i - n÷2 + j-1] * w[j] for j in 1:n)/sum(w)`.
"""
function smooth(x::AbstractVector;  w=pdf(Normal(), range(-3, 3, length=7)))
  sx = Vector{Union{Missing, Float64}}(undef, length(x))
  n = length(w)
  if n % 2 != 1
    error("only odd length windows allowed")
  end
  shift = -(n÷2):(n÷2)
  for i in 1:length(sx)
    s = findfirst((i .+ shift) .> 0)
    l = findlast((i .+ shift) .< length(x))
    sx[i] = sum(w[s:l].*x[i .+ shift[s:l]])./sum(w[s:l])
    if isnan(sx[i])
      sx[i] = 0
    end
  end
  return(sx)
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

import FMIImport: fmi2Real, fmi2StatusOK, fmi2False

import FMIImport: fmi2SetTime, fmi2SetContinuousStates, fmi2SetReal
import FMIImport: fmi2CompletedIntegratorStep, fmi2NewDiscreteStates!, fmi2EventInfo, fmi2GetEventIndicators!
import FMIImport: fmi2GetReal, fmi2GetDerivatives, fmi2GetContinuousStates
import FMIImport: fmi2EnterContinuousTimeMode, fmi2EnterEventMode

using FMIImport: fmi2ComponentStateContinuousTimeMode, fmi2ComponentStateEventMode

mutable struct FMU2ComponentShadow

    # compatibility
    t::Float64
    realComponent
    jacobianUpdate!
    compAddr
    log::Bool

    eventInfo::fmi2EventInfo
    state
    fmu

    jac_ẋy_x::Matrix{fmi2Real}
    jac_ẋy_u::Matrix{fmi2Real}
    jac_x::Array{fmi2Real}
    jac_u::Union{Array{fmi2Real}, Nothing}
    jac_t::fmi2Real

    # derivative vector caching
    t_ẋ_cache::Array{Float64, 1}
    ẋ_cache::Array{Array{Float64, 1}, 1}
    ẋ_interp

    t_y_cache::Array{Float64, 1}
    y_cache::Array{Array{Float64, 1}, 1}
    y_interp

    t_x_cache::Array{Float64, 1}
    x_cache::Array{Array{Float64, 1}, 1}
    x_interp

    t_jac_ẋy_x_cache::Array{Float64, 1}
    jac_ẋy_x_cache::Array{Matrix{Float64}, 1}
    jac_ẋy_x_interp

    t_jac_ẋy_u_cache::Array{Float64, 1}
    jac_ẋy_u_cache::Array{Matrix{Float64}, 1}
    jac_ẋy_u_interp

    function FMU2ComponentShadow()
        inst = new()
        inst.log = false
        inst.realComponent = nothing

        inst.compAddr = C_NULL

        inst.jacobianUpdate! = jacobianUpdate!

        inst.jac_x = Array{fmi2Real, 1}()
        inst.jac_u = nothing
        inst.jac_t = -1.0
        inst.jac_ẋy_x = zeros(fmi2Real, 0, 0)
        inst.jac_ẋy_u = zeros(fmi2Real, 0, 0)

        reset!(inst)

        return inst
    end
end

function setComponent!(comp::FMU2ComponentShadow, realComponent::Union{FMU2Component, Nothing})

    comp.realComponent = realComponent

    if realComponent != nothing
        comp.eventInfo = realComponent.eventInfo
        comp.state = realComponent.state
        comp.fmu = realComponent.fmu
        comp.t = realComponent.t
        comp.compAddr = realComponent.compAddr
    end
end

function jacobianUpdate!(mat, comp::FMU2ComponentShadow, rdx, rx)

    if rx == comp.fmu.modelDescription.stateValueReferences

        if comp.jac_ẋy_x_interp !== nothing
            @assert length(comp.jac_ẋy_x_interp) > 0 "empty interpolation polynominal"

            b = bounds(comp.jac_ẋy_x_interp.itp)[1]
            t = comp.t
            t = max(b[1], min(b[2], t))
            mat[:] = comp.jac_ẋy_x_interp(t)

        else
            comp.realComponent.jacobianUpdate!(mat, comp.realComponent, rdx, rx)

            if isLogging(comp)
                if length(comp.t_jac_ẋy_x_cache) == 0 || comp.t_jac_ẋy_x_cache[end] < comp.t
                    push!(comp.t_jac_ẋy_x_cache, comp.t)
                    push!(comp.jac_ẋy_x_cache, copy(mat))
                end
            end
        end

    elseif rx == comp.fmu.modelDescription.inputValueReferences

        if comp.jac_ẋy_u_interp !== nothing
            @assert length(comp.jac_ẋy_u_interp) > 0 "empty interpolation polynominal"

            b = bounds(comp.jac_ẋy_u_interp.itp)[1]
            t = comp.t
            t = max(b[1], min(b[2], t))
            mat[:] = comp.jac_ẋy_u_interp(t)

        else
            comp.realComponent.jacobianUpdate!(mat, comp.realComponent, rdx, rx)

            if isLogging(comp)
                if length(comp.t_jac_ẋy_u_cache) == 0 || comp.t_jac_ẋy_u_cache[end] < comp.t
                    push!(comp.t_jac_ẋy_u_cache, comp.t)
                    push!(comp.jac_ẋy_u_cache, copy(mat))
                end
            end
        end

    else
        @assert false "shadow jacobianUpdate!, unknown rx = $(rx)"
    end

end

function reset!(componentShadow::FMU2ComponentShadow)
    componentShadow.t_y_cache = []
    componentShadow.t_x_cache = []
    componentShadow.t_ẋ_cache = []
    componentShadow.t_jac_ẋy_x_cache = []
    componentShadow.t_jac_ẋy_u_cache = []

    componentShadow.y_cache = []
    componentShadow.x_cache = []
    componentShadow.ẋ_cache = []
    componentShadow.jac_ẋy_x_cache = []
    componentShadow.jac_ẋy_u_cache = []

    componentShadow.y_interp = nothing
    componentShadow.x_interp = nothing
    componentShadow.ẋ_interp = nothing
    componentShadow.jac_ẋy_x_interp = nothing
    componentShadow.jac_ẋy_u_interp = nothing
end

function prepare!(componentShadow::FMU2ComponentShadow)

    # if length(componentShadow.ẋ_cache) == 0 || length(componentShadow.jac_ẋy_x_cache) == 0
    #     return nothing
    # end

    if length(componentShadow.ẋ_cache) > 0
        componentShadow.ẋ_interp = LinearInterpolation(componentShadow.t_ẋ_cache, componentShadow.ẋ_cache)
        @debug "[ẋ_interp]"
    end

    if length(componentShadow.x_cache) > 0
        componentShadow.x_interp = LinearInterpolation(componentShadow.t_x_cache, componentShadow.x_cache)
        @debug "[x_interp]"
    end

    if length(componentShadow.y_cache) > 0
        componentShadow.y_interp = LinearInterpolation(componentShadow.t_y_cache, componentShadow.y_cache)
        @debug "[y_interp]"
    end

    if length(componentShadow.jac_ẋy_x_cache) > 0
        componentShadow.jac_ẋy_x_interp = LinearInterpolation(componentShadow.t_jac_ẋy_x_cache, componentShadow.jac_ẋy_x_cache)
        @debug "[jac_ẋy_x_interp]"
    end

    if length(componentShadow.jac_ẋy_u_cache) > 0
        componentShadow.jac_ẋy_u_interp = LinearInterpolation(componentShadow.t_jac_ẋy_u_cache, componentShadow.jac_ẋy_u_cache)
        @debug "[jac_ẋy_u_interp]"
    end

    # if componentShadow.realComponent != nothing
    #     if componentShadow.fmu.executionConfig.freeInstance
    #         fmi2FreeInstance!(componentShadow.realComponent)
    #     end
    #     componentShadow.realComponent = nothing
    # end
end

function isPrepared(componentShadow::FMU2ComponentShadow)
    return (componentShadow.ẋ_interp !== nothing)
end

function isLogging(componentShadow::FMU2ComponentShadow)
    return componentShadow.log
end

function fmi2SetTime(comp::FMU2ComponentShadow, t)
    # if isPrepared(comp)
    #     comp.t = t
    #     return fmi2StatusOK
    # else
    comp.t = t
    status = fmi2SetTime(comp.realComponent, t)
    comp.state = comp.realComponent.state
    return status
    # end
end

function fmi2EnterContinuousTimeMode(comp::FMU2ComponentShadow)
    # if isPrepared(comp)
    #     comp.state = fmi2ComponentStateContinuousTimeMode
    #     return fmi2StatusOK
    # else
    status = fmi2EnterContinuousTimeMode(comp.realComponent)
    comp.state = comp.realComponent.state
    return status
    # end
end

function fmi2EnterEventMode(comp::FMU2ComponentShadow)
    # if isPrepared(comp)
    #     comp.state = fmi2ComponentStateEventMode
    #     return fmi2StatusOK
    # else
        status = fmi2EnterEventMode(comp.realComponent)
        comp.state = comp.realComponent.state
        return status
    # end
end

function fmi2SetContinuousStates(comp::FMU2ComponentShadow, x)
    # if isPrepared(comp)
    #     # do nothing
    #     return fmi2StatusOK
    # else
        return fmi2SetContinuousStates(comp.realComponent, x)
    # end
end

function fmi2SetReal(comp::FMU2ComponentShadow, vrs, vals)
    # if isPrepared(comp)
    #     # do nothing
    #     return fmi2StatusOK
    # else

    return fmi2SetReal(comp.realComponent, vrs, vals)
    # end
end

function fmi2CompletedIntegratorStep(comp::FMU2ComponentShadow, val)
    # if isPrepared(comp)
    #     status = fmi2StatusOK
    #     enterEventMode = fmi2False
    #     terminateSimulation = fmi2False
    #     return (status, enterEventMode, terminateSimulation)
    # else
    #if comp.realComponent !== nothing
        status = fmi2CompletedIntegratorStep(comp.realComponent, val)
        comp.state = comp.realComponent.state
        return status
    #end

    #return (fmi2StatusOK, fmi2False, fmi2False)
end

function fmi2NewDiscreteStates!(comp::FMU2ComponentShadow, eventInfo)

    # if isPrepared(comp)
    #     eventInfo.newDiscreteStatesNeeded = fmi2False
    #     eventInfo.terminateSimulation = fmi2False
    #     eventInfo.nominalsOfContinuousStatesChanged = fmi2False
    #     eventInfo.valuesOfContinuousStatesChanged = fmi2False
    #     eventInfo.nextEventTimeDefined = fmi2False
    #     eventInfo.nextEventTime = 0.0

    #     return fmi2StatusOK
    # else
        status = fmi2NewDiscreteStates!(comp.realComponent, eventInfo)
        comp.state = comp.realComponent.state
        return status
    #end
end

function fmi2GetEventIndicators!(comp::FMU2ComponentShadow, out)
    # if isPrepared(comp)
    #     for i in 1:length(out)
    #         out[i] = 1.0
    #     end
    #     return fmi2StatusOK
    # else
        return fmi2GetEventIndicators!(comp.realComponent, out)
    #end
end

function fmi2GetReal(comp::FMU2ComponentShadow, getValueReferences)

    y = nothing

    if (comp.y_interp !== nothing)

        b = bounds(comp.y_interp.itp)[1]
        t = comp.t
        t = max(b[1], min(b[2], t))
        y = comp.y_interp(t)

    else

        y = fmi2GetReal(comp.realComponent, getValueReferences)

        if isLogging(comp)
            if length(comp.t_y_cache) == 0 || comp.t_y_cache[end] < comp.t
                push!(comp.y_cache, collect(ForwardDiff.value(e) for e in y) )
                push!(comp.t_y_cache, comp.t)
            end
        end
    end

    return y
end

function fmi2GetDerivatives(comp::FMU2ComponentShadow)

    ẋ = nothing

    if (comp.ẋ_interp !== nothing)

        b = bounds(comp.ẋ_interp.itp)[1]
        t = comp.t
        t = max(b[1], min(b[2], t))
        ẋ = comp.ẋ_interp(t)

    else

        ẋ = fmi2GetDerivatives(comp.realComponent)

        if isLogging(comp)
            if length(comp.t_ẋ_cache) == 0 || comp.t_ẋ_cache[end] < comp.t
                push!(comp.ẋ_cache, collect(ForwardDiff.value(e) for e in ẋ) )
                push!(comp.t_ẋ_cache, comp.t)
            end
        end

    end

    #@info "derivatives, prepared=$(isPrepared(comp)), ẋ=$(ẋ), length=$(length(comp.t_ẋ_cache))"

    return ẋ
end

function fmi2GetContinuousStates(comp::FMU2ComponentShadow)

    x = nothing

    if (comp.x_interp !== nothing)

        b = bounds(comp.x_interp.itp)[1]
        t = comp.t
        t = max(b[1], min(b[2], t))
        x = comp.x_interp(t)

    else

        x = fmi2GetContinuousStates(comp.realComponent)

        if isLogging(comp)
            if length(comp.t_x_cache) == 0 || comp.t_x_cache[end] < comp.t
                push!(comp.x_cache, x)
                push!(comp.t_x_cache, comp.t)
            end
        end

    end

    return x
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct PacketCaptureStorageLocation <: SwaggerModel
    storageId::Any # spec type: Union{ Nothing, String } # spec name: storageId
    storagePath::Any # spec type: Union{ Nothing, String } # spec name: storagePath
    filePath::Any # spec type: Union{ Nothing, String } # spec name: filePath

    function PacketCaptureStorageLocation(;storageId=nothing, storagePath=nothing, filePath=nothing)
        o = new()
        validate_property(PacketCaptureStorageLocation, Symbol("storageId"), storageId)
        setfield!(o, Symbol("storageId"), storageId)
        validate_property(PacketCaptureStorageLocation, Symbol("storagePath"), storagePath)
        setfield!(o, Symbol("storagePath"), storagePath)
        validate_property(PacketCaptureStorageLocation, Symbol("filePath"), filePath)
        setfield!(o, Symbol("filePath"), filePath)
        o
    end
end # type PacketCaptureStorageLocation

const _property_map_PacketCaptureStorageLocation = Dict{Symbol,Symbol}(Symbol("storageId")=>Symbol("storageId"), Symbol("storagePath")=>Symbol("storagePath"), Symbol("filePath")=>Symbol("filePath"))
const _property_types_PacketCaptureStorageLocation = Dict{Symbol,String}(Symbol("storageId")=>"String", Symbol("storagePath")=>"String", Symbol("filePath")=>"String")
Base.propertynames(::Type{ PacketCaptureStorageLocation }) = collect(keys(_property_map_PacketCaptureStorageLocation))
Swagger.property_type(::Type{ PacketCaptureStorageLocation }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_PacketCaptureStorageLocation[name]))}
Swagger.field_name(::Type{ PacketCaptureStorageLocation }, property_name::Symbol) =  _property_map_PacketCaptureStorageLocation[property_name]

function check_required(o::PacketCaptureStorageLocation)
    true
end

function validate_property(::Type{ PacketCaptureStorageLocation }, name::Symbol, val)
end


text:
#==============================================================================#
# Field initialization
#==============================================================================#

"""
    fieldinit!(sys; seed = -1, scale = 1.0)
    fieldinit!(sys, fields; field_type = :omega, interpolate = true)
    fieldinit!(sys, fpath; load_cell = :true, interpolate = true)

Initialize fields for all monomer types in the system.
Fields can be initialized randomly with no provided data,
loaded from a field dictionary, or loaded from file.
"""
function fieldinit!(sys::FieldSystem; seed::Integer = -1, scale::Real = 0.1)
    if seed >= 0
        Random.seed!(seed)
    end
    for (mid, omega) in sys.fields
        omega .= scale * randn(sys.dims)
    end
    residuals!(sys)
    return nothing
end

function fieldinit!(sys::FieldSystem, fields::Dict; field_type::Symbol = :omega, interpolate::Bool = true)
    if field_type == :omega
        dict = sys.fields
    else
        dict = sys.density
    end

    for (mid, arr) in fields
        if haskey(sys.monomers, mid)
            if size(arr) == sys.dims
                dict[mid] = arr
            elseif interpolate
                itp = interpolate_grid(arr, sys.dims)
                dict[mid] = itp
            else
                error("Input field dimensions $(size(arr)) does not match system grid $(sys.dims).")
            end
        else
            @warn "Monomer with ID = $(mid) not present in system -- skipping."
        end
    end

    # Update based on fields
    if field_type == :rho
        uniform_fields!(sys)
    end
    residuals!(sys)

    return nothing
end

function fieldinit!(sys::FieldSystem, fpath::AbstractString; load_cell::Bool = true, interpolate::Bool = true)
    data = loadfields(fpath)
    if load_cell
        cell = UnitCell(data[:dim], data[:crystal], data[:cell_params])
        add_cell!(sys, cell)
    end
    fieldinit!(sys, data[:fields]; field_type = data[:field_type], interpolate = interpolate)
    return nothing
end

#==============================================================================#
# Field IO
#==============================================================================#

"""
    savefields(sys, fpath; type = :field/:density)

Save fields from the system to a JSON file.
Argument `type` specifies whether density or omega fields are saved.
"""
function savefields(sys::FieldSystem, fpath::AbstractString; field_type::Symbol = :omega)
    d = Dict()
    d[:dim] = ndims(sys.dims)
    d[:grid] = sys.dims
    d[:mids] = sort(collect(keys(sys.monomers)))

    d[:date] = Dates.format(Dates.now(), "Y-m-d H:M:S")
    d[:crystal] = sys.cell.crystal
    d[:cell_params] = sys.cell.params
    d[:field_type] = field_type

    fields = field_type == :rho ? sys.density : sys.fields
    d[:fields] = typeof(fields)()
    for (mid, arr) in fields
        d[:fields][mid] = arr
    end

    json_string = JSON.json(d, 4)#JSON.json(d, 4)
    open(fpath, "w") do f
        write(f, json_string)
    end

    return nothing
end

"""
    loadfields(fpath)

Load fields and system data from a specified JSON file.
"""
function loadfields(fpath::AbstractString)
    raw = JSON.parsefile(fpath)

    d = Dict()
    d[:dim] = Int(raw["dim"])
    d[:grid] = Tuple(convert(Vector{Int}, raw["grid"]))

    d[:mids] = Tuple(convert(Vector{Int}, raw["mids"]))
    d[:crystal] = Symbol(raw["crystal"])
    d[:cell_params] = convert(Vector{Float64}, raw["cell_params"])
    d[:field_type] = Symbol(raw["field_type"])

    mids = d[:mids]
    grid = d[:grid]
    d[:fields] = Dict(mid => zeros(grid) for mid in mids)
    for mid in mids
        arr = d[:fields][mid]
        vecs = raw["fields"][string(mid)]
        for k = 1:grid[end]
            for j = 1:grid[2]
                arr[:,j,k] = vecs[k][j]
            end
        end
    end

    return d
end


text:
# using Flux
using Pkg

Pkg.activate(".")

using JuliaRL

using Random
using ProgressMeter

function test_run(α=0.5/32, ϵ=0.01, γ=1.0, tilings=32, tiles=2; n_episodes=1000, seed=1002349)

    # ϵ = 0.1
    # α = 0.5/tilings
    size_env_state = 2
    num_actions = 3

    lu = Learning.LinearRL.WatkinsQ(α)
    policy = EpsilonGreedyQPolicy(ϵ, 1:3)
    agent = LinearAgents.TileCoderAgent(lu, size_env_state, num_actions, tilings, tiles, γ, policy)

    rng = Random.MersenneTwister(seed)
    env = MountainCar(rng)
    cumulative_reward_array = zeros(Int64, n_episodes)
    # println("Here")
    @showprogress 0.1 "Episode: " for episode = 1:n_episodes
        # for episode = 1:n_episodes
        terminal = false
        num_steps = 0
        cumulative_reward = 0
        _, state = start!(env; rng=rng)
        action = start!(agent, state; rng=rng)
        while !terminal

            # print(num_steps, "\r")
            _, state_prime, reward, terminal = step!(env, action; rng=rng)

            action = step!(agent, state_prime, reward, terminal; rng=rng)
            # println(action)
            # println(reward)
            num_steps += 1
            cumulative_reward += reward

        end
        cumulative_reward_array[episode] = cumulative_reward
    end
    return cumulative_reward_array
end

max_stars_repo_name: mkschleg/JuliaRL

max_issues_repo_name: mkschleg/JuliaRL

max_forks_repo_name: mkschleg/JuliaRL


text:
using Documenter, Images, ImageDraw

include("logo.jl")
generate_logo()

makedocs(sitename="ImageDraw.jl",
            format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"))
deploydocs(repo = "github.com/JuliaImages/ImageDraw.jl.git")


text:
export CosCharFuncPricer, makeCosCharFuncPricer, priceEuropean

struct CosCharFuncPricer{T}
    τ::T
    a::T
    b::T
    uk::Array{T,1}
    phi::Array{T,1}
    pi::T
end

function makeCosCharFuncPricer(
    cf::CharFunc{MAINT,CR},
    τ::T,
    m::Int,
    l::Int,
) where {T,CR,MAINT}
    p = model(cf)
    c1, c2, c4 = computeCumulants(p, τ)
    c2 = c2 + sqrt(abs(c4))
    a = c1 - l * sqrt(abs(c2))
    b = c1 + l * sqrt(abs(c2))
    # println("a ",a," b ",b)
    piHigh = const_pi(cf)
    z = @. (1:m) * piHigh / (b - a)
    phiz = map(z -> evaluateCharFunc(cf, z, τ), z)
    phi = @. real(phiz) * cos(-z * a) - imag(phiz) * sin(-z * a)
    uk = Vector{typeof(piHigh)}(undef, m)
    return CosCharFuncPricer(τ, a, b, uk, phi, piHigh)
end

function makeCosCharFuncPricer(
    cf::CharFunc{MAINT,CR},
    τ::T,
    l::Int;
    tol::T = 1e-8,
) where {T,CR,MAINT}
    p = model(cf)
    c1, c2, c4 = computeCumulants(p, τ)
    c2 = c2 + sqrt(abs(c4))
    a = c1 - l * sqrt(abs(c2))
    b = c1 + l * sqrt(abs(c2))
    c0 = cinf(p, τ)
    lWinv = Float64(2 * c0 / (tol * (b - a)))
    lW = lambertW(lWinv)
    m = ceil(Int, lW / Float64(c0) / Base.pi * Float64(b - a))
    m = max(32, m)
    # println("a ",a," b ",b, " m ",m, " lwi ",lWinv, " c0 ",c0)
    piHigh = const_pi(cf)
    z = @. (1:m) * piHigh / (b - a)
    phiz = map(z -> evaluateCharFunc(cf, z, τ), z)
    phi = @. real(phiz) * cos(-z * a) - imag(phiz) * sin(-z * a)
    uk = Vector{typeof(piHigh)}(undef, m)
    return CosCharFuncPricer(τ, a, b, uk, phi, piHigh)
end

#we adopt here the alternative formula of LeFloch "More Robust Pricing of European Options Based on Fourier Cosine Series Expansions"
function priceEuropean(
    p::CosCharFuncPricer{T},
    isCall::Bool,
    strike::T,
    forward::T,
    τ::T,
    discountDf::T,
) where {T}
    if τ != p.τ
        throw(DomainError(τ, string("maturity is different from pricer maturity ", p.τ)))
    end
    local pricePut
    x = log(forward / strike)
    if x >= -p.a && x >= p.b
        pricePut = 0
    elseif x <= p.a || x <= -p.b
        pricePut = discountDf * (strike - forward)
    else
        uk = p.uk
        a = p.a
        b = p.b
        ea = exp(a)
        f = forward
        piHigh = p.pi
        logStrike = log(strike / f)
        estrike = strike / f
        coeff = (-(estrike - ea) + estrike * (logStrike - a))*2 / (b - a)
        uk0 = coeff

        @inbounds for i = 1:length(uk)
            z = i * piHigh / (b - a)
            kPid = (logStrike - a) * z
            sk, ck = sincos(kPid)
            chi = (ck * estrike - 1 * ea + z * (sk * estrike)) / (1 + z^2)
            psi = sk / z
            coeff = (-chi + estrike * psi)*2 / (b - a)
            uk[i] = coeff
        end
        sumPut = uk0 / 2
        @inbounds for k = 1:length(uk)
            sumPut += p.phi[k] * p.uk[k]
        end
        pricePut = discountDf * f * sumPut
    end
    if isCall
        return pricePut + discountDf * (forward - strike)
    end
    return pricePut
end


text:
using Test
using FlightMechanics

@testset "atmosphere" begin include("atmosphere.jl") end
@testset "coordinates" begin include("coordinates.jl") end
@testset "anemometry" begin include("anemometry.jl") end
@testset "mechanics" begin include("mechanics.jl") end
@testset "flight mechanics" begin include("flight_mechanics.jl") end
@testset "inputs" begin include("models/inputs.jl") end
@testset "fcs" begin include("models/fcs.jl") end
@testset "pfm" begin include("models/point_forces_moments.jl") end
@testset "mass" begin include("models/mass.jl") end
@testset "attitude" begin include("models/attitude.jl") end
@testset "position" begin include("models/position.jl") end
@testset "aerodynamics" begin include("models/aerodynamics.jl") end
@testset "aerostate" begin include("models/aero_state.jl") end
@testset "dynamic system" begin include("models/dynamic_system.jl") end
@testset "trimmer" begin include("models/trimmer.jl") end

@testset "ac: c310" begin include("aircrafts/c310.jl") end
@testset "ac: f16" begin include("aircrafts/f16.jl") end


text:
function maxpartsum_in(a)
    sum = 0
    mins = 0
    answer = 0
    k = 0
    m = 0
    min_p = 0
    for i in 1:lastindex(a)
        sum = sum+a[i]
        curr = sum - mins

        if curr > answer
            answer = curr
            k = min_p +1
            m = i
        end

        if sum < mins
            mins = sum
            min_p = i
        end
    end

    return k, m
end


text:
# Rocket.jl message passing specific extensions

using Rocket

mutable struct LimitStackSchedulerProps
    soft_depth :: Int
    hard_depth :: Int
end

struct LimitStackScheduler <: Rocket.AbstractScheduler
    soft_limit :: Int
    hard_limit :: Int
    props      :: LimitStackSchedulerProps
end

LimitStackScheduler(soft_limit::Int)                  = LimitStackScheduler(soft_limit, typemax(Int) - 1)
LimitStackScheduler(soft_limit::Int, hard_limit::Int) = LimitStackScheduler(soft_limit, hard_limit, LimitStackSchedulerProps(0, 0))

get_soft_limit(scheduler::LimitStackScheduler) = scheduler.soft_limit
get_hard_limit(scheduler::LimitStackScheduler) = scheduler.hard_limit

function increase_depth!(scheduler::LimitStackScheduler)
    scheduler.props.soft_depth = scheduler.props.soft_depth + 1
    scheduler.props.hard_depth = scheduler.props.hard_depth + 1
end

function decrease_depth!(scheduler::LimitStackScheduler)
    scheduler.props.soft_depth = scheduler.props.soft_depth - 1
    scheduler.props.hard_depth = scheduler.props.hard_depth - 1
end

get_soft_depth(scheduler::LimitStackScheduler)     = scheduler.props.soft_depth
set_soft_depth!(scheduler::LimitStackScheduler, v) = scheduler.props.soft_depth = v

get_hard_depth(scheduler::LimitStackScheduler)     = scheduler.props.hard_depth
set_hard_depth!(scheduler::LimitStackScheduler, v) = scheduler.props.hard_depth = v

Base.show(io::IO, scheduler::LimitStackScheduler) = print(io, "LimitStackScheduler(soft_limit = $(get_soft_limit(scheduler)), hard_limit = $(get_hard_limit(scheduler)))")

Base.similar(scheduler::LimitStackScheduler) = LimitStackScheduler(get_soft_limit(scheduler), get_hard_limit(scheduler))

Rocket.makeinstance(::Type, scheduler::LimitStackScheduler) = scheduler

Rocket.instancetype(::Type, ::Type{ <: LimitStackScheduler }) = LimitStackScheduler

macro limitstack(no_limit_cb, limit_cb)
    output = quote
        increase_depth!(instance)

        if get_hard_depth(instance) >= get_hard_limit(instance)
            error("Hard limit in LimitStackScheduler exceeded")
        end

        result = if get_soft_depth(instance) < get_soft_limit(instance)
            begin
                $no_limit_cb
            end
        else
            previous_soft_depth = get_soft_depth(instance)
            set_soft_depth!(instance, 0)
            r = begin
                $limit_cb
            end
            set_soft_depth!(instance, previous_soft_depth)
            r
        end

        decrease_depth!(instance)

        result
    end
    return esc(output)
end

struct LimitStackSubscription <: Teardown
    instance     :: LimitStackScheduler
    subscription :: Teardown
end

Rocket.as_teardown(::Type{ <: LimitStackSubscription }) = UnsubscribableTeardownLogic()

function Rocket.on_unsubscribe!(l::LimitStackSubscription)
    instance     = l.instance
    subscription = l.subscription
    @limitstack Rocket.unsubscribe!(subscription) begin
        @sync @async Rocket.unsubscribe!(subscription)
    end
    return nothing
end

function Rocket.scheduled_subscription!(source, actor, instance::LimitStackScheduler)
    subscription = @limitstack Rocket.on_subscribe!(source, actor, instance) begin
        condition  = Base.Condition()
        @async begin
            try
                subscription = Rocket.on_subscribe!(source, actor, instance)
                notify(condition, subscription)
            catch exception
                notify(condition, exception, error = true)
            end
        end
        wait(condition)
    end
    return LimitStackSubscription(instance, subscription)
end

function Rocket.scheduled_next!(actor, value, instance::LimitStackScheduler)
    @limitstack Rocket.on_next!(actor, value) begin
        @sync @async Rocket.scheduled_next!(actor, value, instance)
    end
end

function Rocket.scheduled_error!(actor, err, instance::LimitStackScheduler)
    @limitstack Rocket.on_error!(actor, err) begin
        @sync @async Rocket.scheduled_error!(actor, err, instance)
    end
end

function Rocket.scheduled_complete!(actor, instance::LimitStackScheduler)
    @limitstack Rocket.on_complete!(actor) begin
        @sync @async Rocket.scheduled_complete!(actor, instance)
    end
end


text:
module RawMetaMacros

import Docile.Formats: metamacro, @META_str

function metamacro(::META"raw_meta"raw, body, mod, obj)
    body
end

function metamacro(::META"nested_meta", body, mod, obj)
    body
end

"!!raw_meta(!!undefined())"
raw_metamacro = ()

"!!nested_meta(!!var(metamacro_type:nestable))"
nested_metamacro = ()

end


text:
"""
    NormalizeRow(dict, cols)

Normalizes the values of a row present in `TabularItem` for the columns
specified in `cols` using `dict`, which contains the column names as
dictionary keys and the mean and standard deviation tuple present as
dictionary values.

## Example

```julia
using DataAugmentation

cols = [:col1, :col2, :col3]
row = (; zip(cols, [1, 2, 3])...)
item = TabularItem(row, cols)
normdict = Dict(:col1 => (1, 1), :col2 => (2, 2))

tfm = NormalizeRow(normdict, [:col1, :col2])
apply(tfm, item)
```
"""
struct NormalizeRow{T, S} <: Transform
    dict::T
    cols::S
end

function apply(tfm::NormalizeRow, item::TabularItem; randstate=nothing)
    x = NamedTuple(Iterators.map(item.columns, item.data) do col, val
        if col in tfm.cols
            colmean, colstd = tfm.dict[col]
            val = (val - colmean)/colstd
        end
        (col, val)
    end)
    TabularItem(x, item.columns)
end

"""
    FillMissing(dict, cols)

Fills the missing values of a row present in `TabularItem` for the columns
specified in `cols` using `dict`, which contains the column names as
dictionary keys and the value to fill the column with present as
dictionary values.

## Example

```julia
using DataAugmentation

cols = [:col1, :col2, :col3]
row = (; zip(cols, [1, 2, 3])...)
item = TabularItem(row, cols)
fmdict = Dict(:col1 => 100, :col2 => 100)

tfm = FillMissing(fmdict, [:col1, :col2])
apply(tfm, item)
```
"""
struct FillMissing{T, S} <: Transform
    dict::T
    cols::S
end

function apply(tfm::FillMissing, item::TabularItem; randstate=nothing)
    x = NamedTuple(Iterators.map(item.columns, item.data) do col, val
        if col in tfm.cols && ismissing(val)
            val = tfm.dict[col]
        end
        (col, val)
    end)
    TabularItem(x, item.columns)
end

"""
    Categorify(dict, cols)

Label encodes the values of a row present in `TabularItem` for the
columns specified in `cols` using `dict`, which contains the column
names as dictionary keys and the unique values of column present
as dictionary values.

if there are any `missing` values in the values to be transformed,
they are replaced by 1.

## Example

```julia
using DataAugmentation

cols = [:col1, :col2, :col3]
row = (; zip(cols, ["cat", 2, 3])...)
item = TabularItem(row, cols)
catdict = Dict(:col1 => ["dog", "cat"])

tfm = Categorify(catdict, [:col1])
apply(tfm, item)
```
"""
struct Categorify{T, S} <: Transform
    dict::T
    cols::S
    function Categorify{T, S}(dict::T, cols::S) where {T, S}
        for (col, vals) in dict
            if any(ismissing, vals)
                dict[col] = filter(!ismissing, vals)
                @warn "There is a missing value present for category '$col' which will be removed from Categorify dict"
            end
        end
        new{T, S}(dict, cols)
    end
end

Categorify(dict::T, cols::S) where {T, S} = Categorify{T, S}(dict, cols)

function apply(tfm::Categorify, item::TabularItem; randstate=nothing)
    x = NamedTuple(Iterators.map(item.columns, item.data) do col, val
        if col in tfm.cols
            val = ismissing(val) ? 1 : findfirst(val .== tfm.dict[col]) + 1
        end
        (col, val)
    end)
    TabularItem(x, item.columns)
end


text:
using ArgParse
include("Common.jl")

function parse(edge_feats_list, sentences, weights, fpath_out)
	info("Start parsing ...")
	f_out = open(fpath_out, "w")

	for (i, sent) in enumerate(sentences)
		n = length(sent)
		edge_feats = edge_feats_list[i]
		edge_scores = getedgescores(edge_feats, weights, n)
		heads = eisner(edge_scores, n)

		for (j, token) in enumerate(sent)
			token_array = [token.id,
							token.form,
							token.lemma,
							token.cpostag,
							token.postag,
							token.feats,
							heads[j], # predicted head
							token.deprel,
							"\n"]
			write(f_out, join(token_array, "\t") )
		end
		write(f_out, "\n")
	end
end

function getoptions(args)
	s = ArgParseSettings("Dependency parser in Julia.")
	@add_arg_table s begin
		"input"
			help = "input file"
			arg_type = String
			required = true
		"model"
			help = "model file"
			arg_type = String
			required = true
		"output"
			help = "parsed output file"
			arg_type = String
			required = true
	end
	parsed_args = parse_args(args, s)

	@assert isreadable(parsed_args["input"]) # check if training file readable
	@assert isreadable(parsed_args["model"]) # check if model file readable

	return parsed_args["input"], parsed_args["model"], parsed_args["output"]
end

function main(args)
	info("PARSING")
	const (fpath_in, fpath_model, fpath_out) = getoptions(args)
	info("Input file:\t$fpath_in")
	info("Model file:\t$fpath_model")
	info("Output file:\t$fpath_out")
	info("-------------------")

	const sentences = readconlldata(fpath_in)
	info("-------------------")

	info("Loading model ...")
	const model = open(deserialize, fpath_model)
	info("-------------------")

	const edge_feats_list = createfeatures(sentences, model.featdict)
	info("-------------------")

	parse(edge_feats_list, sentences, model.weights, fpath_out)

	info("... All DONE!")
end

main(ARGS)


text:
# ---
# title: 1328. Break a Palindrome
# id: problem1328
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/break-a-palindrome/description/>
# hidden: true
# ---
#
# Given a palindromic string `palindrome`, replace **exactly one** character by
# any lowercase English letter so that the string becomes the lexicographically
# smallest possible string that **isn 't** a palindrome.
#
# After doing so, return the final string.  If there is no way to do so, return
# the empty string.
#
#
#
# **Example 1:**
#
#
#
#     Input: palindrome = "abccba"
#     Output: "aaccba"
#
#
# **Example 2:**
#
#
#
#     Input: palindrome = "a"
#     Output: ""
#
#
#
#
# **Constraints:**
#
#   * `1 <= palindrome.length <= 1000`
#   * `palindrome` consists of only lowercase English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# -----------------------------------------------------------
#
# Construction of multidimensional DG derivative matrices
# based on the 1-D derivative matrix
#
# -----------------------------------------------------------

# Efficiency criticality: HIGH
# Matrix computations are only performed once,
# but this can be the main bottleneck if not done right

# Accuracy criticality: HIGH
# Critical for accurate PDE evolution

function D_matrix(d::Int, k::Int, n::Int, srefVD::Array{NTuple{3, CartesianIndex{D}}, 1},
        srefDV::Dict{NTuple{3, CartesianIndex{D}}, Int}; scheme = "sparse") where D
    D_matrix(Val(d), k, n, srefVD, srefDV, Val(Symbol(scheme)))
end
function D_matrix(::Val{d}, k::Int, n::Int, srefVD::Array{NTuple{3, CartesianIndex{D}}, 1},
        srefDV::Dict{NTuple{3, CartesianIndex{D}}, Int}, scheme::Val{Scheme}) where {d, D, Scheme}

    # len = length(srefVD[:,1])
    len = length(srefVD)
    @assert length(srefVD[:,1]) == length(srefVD)
    V2D_1D = V2Dref(Val(1),k,n,Val(:sparse)) # Why not 'scheme'?
    D2V_1D = D2Vref(Val(1),k,n,Val(:sparse))
    I = Int[]; J = Int[]; V = Float64[]

    # 1-dimensional derivative matrix - this will give all coefficient info
    Dmat_1D = periodic_DLF_matrix(k, n)

    for j in 1:len
        lcm = srefVD[j]
        l = lcm[1][d]
        c = lcm[2][d]
        m = lcm[3][d]
        j_1D = D2V_1D[(CartesianIndex(l),CartesianIndex(c),CartesianIndex(m))]
        derivs = Dmat_1D[:, j_1D]::SparseVector{Float64,Int}

        for i_1D in derivs.nzind
            lcm_1D = V2D_1D[i_1D]
            # using tuple constructor independent of Julia's metaprogramming:
            # make_cartesian_index(d, arr1, arr2) takes arr::CartesianIndex{1}
            # and makes a new CartesianIndex{D} using arr2::CartesianIndex{D}
            # with the dth value replaced by arr1[1]
            level2 = make_cartesian_index(d, lcm_1D[1], lcm[1])
            cutoff(scheme, level2, n) && continue
            cell2 = make_cartesian_index(d, lcm_1D[2], lcm[2])
            mode2 = make_cartesian_index(d, lcm_1D[3], lcm[3])
            i = srefDV[(level2, cell2, mode2)]
            push!(I, i)
            push!(J, j)
            push!(V, derivs[i_1D])
        end
    end
    # dropzeros! does not seem helpful for this matrix:
    return sparse(I, J, V, len, len, +)
end

function D_matrix(D::Int, d::Int, k::Int, n::Int; scheme="sparse")
    VD = V2Dref(Val(D), k, n, Val(Symbol(scheme)))
    DV = D2Vref(Val(D), k, n, Val(Symbol(scheme)))
    return D_matrix(d, k, n, VD, DV; scheme=scheme)
end

function grad_matrix(D::Int, k::Int, n::Int; scheme="sparse")
    return [D_matrix(D, d, k, n; scheme=scheme) for d in 1:D]
end

function laplacian_matrix(D::Int, k::Int, n::Int; scheme="sparse")
    len = get_size(Val(D), k, n, Val(Symbol(scheme)))
    lap = spzeros(len, len)
    for i in 1:D
        D_op = D_matrix(D, i, k, n; scheme=scheme)
        lap += D_op * D_op
    end
    return lap
end


text:
struct ComplementBridge{T, B<:MOIBC.SetMapBridge{T}, F, S} <: MOIBC.SetMapBridge{T, MOI.AbstractVectorSet, S, MOI.AbstractFunction, MOI.AbstractFunction}
    con_idx::CI
end

function MOI.supports_constraint(
    ::Type{<:ComplementBridge{T, B}},
    ::Type{F},
    ::Type{<:CS.ComplementSet{CF,S}}
) where {T, B, F<:MOI.AbstractFunction, CF, S}
    if CF <: MOI.ScalarAffineFunction
        is_supported = MOI.supports_constraint(B, CF, S)
    else
        is_supported = MOI.supports_constraint(B, MOIU.scalar_type(CF), S)
    end
    !is_supported && return false
    S <: AbstractBoolSet && return true

    concrete_B = MOIBC.concrete_bridge_type(B, MOI.ScalarAffineFunction{T}, S)
    supported = supports_concreteB(concrete_B)
    return supported
end

function MOIBC.concrete_bridge_type(
    ::Type{<:ComplementBridge{T,B}},
    G::Type{<:MOI.AbstractFunction},
    ::Type{<:CS.ComplementSet{F,S}},
) where {T,B,F,S}
    if S <: AbstractBoolSet
        concrete_B = B
    else
        concrete_B = MOIBC.concrete_bridge_type(B, MOI.ScalarAffineFunction{T}, S)
    end
    return ComplementBridge{T,concrete_B,F,S}
end

function MOIB.added_constraint_types(
    ::Type{<:ComplementBridge{T,B,F,S}}
) where {T,B,F,S}
    if S <: AbstractBoolSet
        added_constraints = added_constraint_types(B, S)
    else
        added_constraints = MOIB.added_constraint_types(B)
    end
    return [(F, CS.ComplementSet{F,added_constraints[1][2]})]
end

function MOIB.added_constrained_variable_types(::Type{<:ComplementBridge{T,B}}) where {T,B}
    return MOIB.added_constrained_variable_types(B)
end

function MOIBC.bridge_constraint(::Type{<:ComplementBridge{T, B, F, S}}, model, fct, set) where {T, B, F, S}
    new_fct = map_function(B, fct, set.set)
    new_inner_set = MOIBC.map_set(B, set.set)
    new_set = CS.ComplementSet{F}(new_inner_set)
    if (new_fct isa SAF) && F <: SAF
        new_fct = get_saf(new_fct)
    end
    return ComplementBridge{T,B,F,S}(MOI.add_constraint(model, new_fct, new_set))
end

function map_function(
    bridge::Type{<:ComplementBridge{T, B}},
    fct,
    set::ComplementSet{F,S}
) where {T,B,F,S<:AbstractBoolSet}
    return map_function(B, fct, set.set)
end

function MOIBC.map_function(
    bridge::Type{<:ComplementBridge{T, B}},
    fct,
) where {T,B}
    return MOIBC.map_function(B, fct)
end

function MOIBC.map_set(
    bridge::Type{<:ComplementBridge{T, B}},
    set::ComplementSet{F},
) where {T,B,F}
    mapped_set = MOIBC.map_set(B, set.set)
    return ComplementSet{F}(mapped_set)
end


text:
using CUDA
using CuCountMap
using StatsBase: countmap

v = rand(Int16, 100_000_000);

using BenchmarkTools

@benchmark gpu_countmap = cucountmap($v)

# BenchmarkTools.Trial:
#   memory estimate:  4.17 MiB
#   allocs estimate:  190
#   --------------
#   minimum time:     41.275 ms (0.00% GC)
#   median time:      44.494 ms (0.00% GC)
#   mean time:        52.756 ms (0.48% GC)
#   maximum time:     297.796 ms (3.98% GC)
#   --------------
#   samples:          95
#   evals/sample:     1

@benchmark cpu_countmap = countmap($v)
# BenchmarkTools.Trial:
#   memory estimate:  4.17 MiB
#   allocs estimate:  37
#   --------------
#   minimum time:     132.618 ms (0.00% GC)
#   median time:      134.176 ms (0.00% GC)
#   mean time:        134.874 ms (0.19% GC)
#   maximum time:     145.168 ms (6.76% GC)
#   --------------
#   samples:          38
#   evals/sample:     1

cuv = CUDA.cu(v)
@benchmark gpu_countmap2 = countmap(cuv)

# BenchmarkTools.Trial:
#   memory estimate:  4.17 MiB
#   allocs estimate:  97
#   --------------
#   minimum time:     5.472 ms (0.00% GC)
#   median time:      5.768 ms (0.00% GC)
#   mean time:        6.125 ms (3.91% GC)
#   maximum time:     201.707 ms (96.90% GC)
#   --------------
#   samples:          816
#   evals/sample:     1


text:
using Documenter
using TaylorIntegration

makedocs(
    modules  = [TaylorIntegration],
    format   = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"),
    clean = false,
    sitename = "TaylorIntegration.jl",
    authors  = "Jorge A. Pérez-Hernández and Luis Benet",
    pages    = [
        "Home" => "index.md",
        "Background" => [
            "taylor_method.md",
            "lyapunov_spectrum.md",
            "jet_transport.md"
            ],
        "Examples" => [
            "simple_example.md",
            "kepler.md",
            "lorenz_lyapunov.md",
            "pendulum.md",
            "root_finding.md",
            "common.md"
            ],
        "Optimizing: `@taylorize`" => "taylorize.md",
        "API"  => "api.md",
        ]
)

deploydocs(
    repo   = "github.com/PerezHz/TaylorIntegration.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing,
    push_preview = true
)


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

module UnitCommitment
    include("log.jl")
    include("dotdict.jl")
    include("instance.jl")
    include("screening.jl")
    include("components.jl")
    include("variables.jl")
    include("constraints.jl")
    include("formulation.jl")
    #include("model.jl")
    include("model2.jl")
    include("sensitivity.jl")
    include("validate.jl")
    include("convert.jl")
    include("initcond.jl")
end


text:
using Documenter, Simulate

makedocs(
    modules = [Simulate],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "Simulate.jl",
    authors  = "Paul Bayer",
    pages = [
        "Home" => "index.md",
        "News" => "news.md",
        "Manual" => [
            "manual/intro.md",
            "manual/approach.md",
            "manual/parallel.md",
            "manual/timer.md",
            "manual/usage.md",
            "Performance" => [
                "performance/performance.md",
                "performance/par_perf.md",
                "performance/benchmarks.md"]],
        "Examples" => [
            "examples/examples.md",
            "examples/tabletennis.md",
            "examples/singleserver.md",
            "examples/postoffice/postoffice.md",
            "examples/dicegame/dicegame.md",
            "examples/house_heating/house_heating.md"],
        "Internals" => "manual/internals.md",
        "Troubleshooting" => "manual/troubleshooting.md",
        "Release notes" => "manual/history.md"
    ]
)

deploydocs(
    repo   = "github.com/pbayer/Simulate.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing,
    devbranch = "master",
    devurl = "dev",
    versions = ["stable" => "v^", "v#.#", "dev" => "dev"]
)


text:
function loadWBMeta(data::String)
    ## load data in fastest possible way
    ## data will be in cache afterwards
    ## without manual interference, file will be on disk
    if data == "countries"
        return loadCountryMeta()
    elseif data == "indicators"
        return loadIndicatorMeta()
    else
        error("Metadata is either countries or indicators")
    end
end

function loadCountryMeta()
    ## load data from cache, disk or internet
    ## data will be in cache afterwards
    if country_cache == false
        filename = joinpath(Pkg.dir("WorldBankDataTd"),
                            "data/countryMeta.csv")
        if isfile(filename) # load from disk
            countryData = readtable(filename)
            set_country_cache(countryData)
            return countryData
        else
            # download data, store to disk and load to cache
            return getWBMeta("countries")
        end

    else # return data from cache
        return country_cache
    end
end

function loadIndicatorMeta()
    ## load data from cache, disk or internet
    ## data will be in cache afterwards
    if indicator_cache == false
        filename = joinpath(Pkg.dir("WorldBankDataTd"),
                            "data/indicatorMeta.csv")
        if isfile(filename) # load from disk
            indicatData = readtable(filename)
            set_indicator_cache(indicatData)
            return indicatData
        else
            # download data, store to disk and load to cache
            return getWBMeta("indicators")
        end

    else # return data from cache
        return indicator_cache
    end
end

function getWBMeta(data::String)
    ## download data, save to disk and load in cache
    if data == "countries"
        countryData = download_countries()
        filename = joinpath(Pkg.dir("WorldBankDataTd"),
                            "data/countryMeta.csv")
        writetable(filename, countryData)
        set_country_cache(countryData)
        return countryData
    elseif data == "indicators"
        indicatData = download_indicators()
        filename = joinpath(Pkg.dir("WorldBankDataTd"),
                            "data/indicatorMeta.csv")
        writetable(filename, indicatData)
        set_indicator_cache(indicatData)
        return indicatData
    else
        error("Metadata is either countries or indicators")
    end
end

function parse_indicator(json::Array{Any,1})
    indicator = ASCIIString[]
    name = UTF8String[]
    description = UTF8String[]
    source_database = UTF8String[]
    source_databaseId = Int[]
    source_organization = UTF8String[]

    for d in json[2]
        append!(indicator,[d["id"]])
        append!(name,[d["name"]])
        append!(description,[d["sourceNote"]])
        append!(source_database, [d["source"]["value"]])
        push!(source_databaseId, int(d["source"]["id"]))
        append!(source_organization, [d["sourceOrganization"]])
    end

    DataFrame(indicator = indicator, name = name,
              description = description,
              source_database = source_database,
              source_databaseId = source_databaseId,
              source_organization = source_organization)
end

function parse_country(json::Array{Any,1})

    ## never missing
    iso3c = ASCIIString[] # no NAs
    iso2c = ASCIIString[] # no NAs
    name = UTF8String[] # no NAs

    capital = DataArray(UTF8String, 0)
    region = DataArray(UTF8String, 0)               # make DataArray? aggr.
    regionId = DataArray(UTF8String, 0)             # make DataArray -> "NA"
    income = DataArray(UTF8String, 0)               # make DataArray? aggr.
    incomeId = DataArray(UTF8String, 0)             # make DataArray
    lending = DataArray(UTF8String, 0)              # make DataArray? aggr.
    lendingId = DataArray(UTF8String, 0)            # make DataArray

    longitude = DataArray(Float64, 0)            # make it DataArray(Float)
    latitude = DataArray(Float64, 0)             # make it DataArray(Float)

    for d in json[2]
        ## never missing values
        append!(iso3c,[d["id"]])
        append!(iso2c,[d["iso2Code"]])
        append!(name,[d["name"]])

        ## append for DataArrays? "" and "NA" to NA
        clean_and_push!(capital,d["capitalCity"])
        clean_and_push!(region,d["region"]["value"])
        clean_and_push!(income,d["incomeLevel"]["value"])
        clean_and_push!(lending,d["lendingType"]["value"])
        clean_and_push!(regionId,d["region"]["id"])
        clean_and_push!(incomeId,d["incomeLevel"]["id"])
        clean_and_push!(lendingId,d["lendingType"]["id"])

        ## numeric values or Void -> NA
        push!(longitude, tofloat(d["longitude"]))
        push!(latitude, tofloat(d["latitude"]))
    end

    ## longitude = [tofloat(i) for i in longitude] # comprehension forces
                                        # NA to NaN if convert(Float,
                                        # NA) is defined
    ## latitude = [tofloat(i) for i in latitude]

    DataFrame(iso3c = iso3c, iso2c = iso2c, name = name,
              region = region, regionId = regionId,
              capital = capital, longitude = longitude,
              latitude = latitude, income = income,
              incomeId = incomeId, lending = lending,
              lendingId = lendingId)
end

function download_indicators()
    dat = download_parse_json("http://api.worldbank.org/indicators?per_page=25000&format=json")

    parse_indicator(dat)
end

function download_countries()
    dat = download_parse_json("http://api.worldbank.org/countries/all?per_page=25000&format=json")

    parse_country(dat)
end

country_cache = false
indicator_cache = false

function set_country_cache(df::AbstractDataFrame)
    global country_cache = df
    if any(isna(country_cache[:iso2c])) # the iso2c code for North Africa is NA
        country_cache[:iso2c][convert(DataArray{Bool,1}, isna(country_cache[:iso2c]))]="NA"
    end
end

function set_indicator_cache(df::AbstractDataFrame)
    global indicator_cache = df
end


text:
#! /usr/bin/env julia

write(STDOUT, "I'm writing some stuff to STDOUT\n")

f = open("stream_test.txt", "w")

write(f, "I'm writing some stuff to a file\n")

close(f)

write(STDOUT, "I should have just wrote to a file.\n")

max_stars_repo_name: sjkelly/JuliaPlayground

max_issues_repo_name: sjkelly/JuliaPlayground

max_forks_repo_name: sjkelly/JuliaPlayground


text:
using CommutativeRings
import CommutativeRings: ord2expo, expo2ord, cbin, varnames

const Z = ZZ{Int}

@testset "index calculations" begin

    @test_throws ArgumentError cbin(0, 0)
    @test cbin(0, 1) == (0, 0)
    @test cbin(10, 1) == (10, 10)
    @test cbin(10, 9) == (10, 10)
    @test cbin(100, 9) == (11, 55)

    @test_throws ArgumentError ord2expo(0, 1) == [7654320]
    @test ord2expo(7654321, 1) == [7654320]
    @test expo2ord([1234]) == 1235
    @test ord2expo(1, 3) == [0, 0, 0]
    @test expo2ord([0,0,0,0]) == 1
    @test ord2expo(binomial(20+4-1, 4)+1, 4) == [0,0,0,20]
    @test expo2ord([0,0,10]) == binomial(10 + 3 - 1, 3)+1

    @test expo2ord.(ord2expo.(1:10000, 2)) == 1:10000
    @test expo2ord.(ord2expo.(1:10000, 5)) == 1:10000
end

@testset "constructors" begin

    @test Z[] == Vector{Z}(undef, 0)
    @test Z[:x] <: UnivariatePolynomial
    @test Z[:x, :y] <: MultivariatePolynomial
    @test Z[[:x],[:y]] <: MultivariatePolynomial
end

@testset "varnames and generators $P " for P in (Z[:x, :y, :z], Z[[:x], [:y,:z]])
    x, y, z = monom.(P, [[1,0,0], [0,1,0], [0,0,1]])
    @test [x, y, z] == generators(P)
    @test varnames(P) == [:x, :y, :z]
    @test varnames(P(0)) == [:x, :y, :z]

    p = 3x^3*y + 4x^2*y^2 + 5x*y*z^2
    @test deg(p) == 4
    @test multideg(p) == [3,1,0]
    @test LC(p) == 3
    @test LM(p) == x^3*y
    @test LT(p) == 3*x^3*y

    @test basetype(p) == Z
    @test copy(p) == p
    q = copy(p)
    q.coeff[1] = 100
    @test p.coeff[1] == 5
    @test q.coeff[1] == 100
end

@testset "addition" for P in (Z[:x, :y], Z[[:x], [:y]])

    @test zero(P) + one(P) == one(P)
    @test one(P) + one(P) == 2*one(P)
    @test -one(P) + one(P) == zero(P)
    @test -one(P) == (-1) * one(P)
    @test one(P) - one(P) == zero(P)

    x, y = generators(P)
    a = 4*x^3 + 3
    b = y + 5*x^3
    @test 5a - 4b == -4y + 15
    @test a + b == 5*x^3 + 5*y + 3
    @test deg(zero(P)) == -1
    @test deg(one(P)) == 0
    @test deg(x + 1) == 1
    @test deg(x*y + 1) == 2
end

@testset "multiplication" for P in (Z[:x, :y], Z[[:x], [:y]])

    @test zero(P) * zero(P) == zero(P)
    @test zero(P) * one(P) == zero(P)
    @test one(P) * one(P) == one(P)
    x, y = generators(P)
    @test 0*x == zero(P)
    @test x*1 === x
    @test (x + y)^2 == x^2 + 2x*y + y^2
    @test (x + y) * (x - y) == x^2 - y^2
    xy = x * 5 + 2y^2 + x*y
    @test one(P) * xy == xy
    @test xy * zero(P) == zero(P)
    @test xy^2 == 25x^2 + x^2*y*10 + 20x*y^2 + x^2*y^2 + 4x*y^3 + 4y^4
    @test xy * x == 5x^2 + 2x*y^2 + x^2*y

    R = (ZZ/12)[:x,:y]
    xr = monom(R, [1, 0])
    @test (xr*3) * 4 == zero(R)
end

@testset "division" begin
    P = ZZ{BigInt}[:x, :y]
    x, y = generators(P)
    f = x^2 - y
    g = x^3 - x
    G = [f, g]
    h = f * (x^2+y^2) + g * (x + y) + 1
    a, r, d = pdivrem(h, G)
    @test d == 1
    @test sum(a .* G) + r == h

    @test rem(x + y, one(P)) == 0
    s, r = divrem(x + y, one(P))
    @test iszero(r) && s == x + y
end

@testset "Gröbner base" for P in (Z[:x, :y], Z[[:x], [:y]])
    x, y = generators(P)
    @test groebnerbase([x, x]) == [x]
    @test hash(0*x) == hash(0)
    @test hash(x^0) == hash(1)
    f = x^2 - y
    g = x^3 - x
    @test groebnerbase([f, g]) == [x^2 - y, x*y - x, y^2 - y]
end

@testset "blocked order" begin
    P = Z[[:t], [:x,:y]]
    t, x, y = generators(P)
    @test (t + 1)^2 == t^2 + 2*t + 1
    @test (y + 1)^2 == y^2 + 2*y + 1
    @test (x + t)^2 == x^2 + 2*x*t + t^2
    @test (x+y)(1, 2, 3) == 5
    @test (t*x^2*y^3 + 1)(1, 2, 3) == 1*4*27 + 1
end

@testset "conversions" for P in (Z[:t, :x, :y], Z[[:t], [:x,:y]])
    t, x, y = generators(P)
    @test generators(P) == [t, x, y]

    Q = Z[:x]
    xq = monom(Q)
    @test P(0*xq) isa P
    @test P(0*xq) == zero(P)
    @test P(5*xq^0) == 5
    @test P(xq) == monom(P, [0, 1, 0])

    R = Z[:z]
    zr = monom(R)
    @test P(0*zr) isa P
    @test P(0*zr) == zero(P)
    @test P(5*zr^0) == 5
    @test_throws ArgumentError P(zr)

    Q = Z[:y, :x]
    yq, xq = monom(Q, [1,0]), monom(Q, [0,1])
    @test P(0*xq) isa P
    @test P(0*yq) == zero(P)
    @test P(5*xq^0) == 5
    @test P(xq) == monom(P, [0, 1, 0])
    @test P(xq + yq*xq) == x + y*x

    Q = Z[:z, :y, :x]
    zq, yq, xq = monom(Q, [1,0,0]), monom(Q, [0,1,0]), monom(Q, [0,0,1])
    @test P(0*xq) isa P
    @test P(0*yq) == zero(P)
    @test P(5*xq^0) == 5
    @test P(xq) == monom(P, [0, 1, 0])
    @test P(xq + yq*xq) == x + y*x
    @test_throws ArgumentError P(zq + xq + yq)
end

@testset "extension and elimination" begin
    Z = ZZ/7583
    P = Z[:x, :y]
    x, y = generators(P)
    A = [x^2-y^2, y^3-2x*y-y^2+2x, x*y^2-3x*y+2x]
    B = [x*y]
    Q = lextend(P, :t)
    ida = Ideal(A)
    idb = Ideal(B)
    t, = generators(Q)
    idcQ = Ideal([Q.(A) .* t; Q.(B) .* (1-t)])
    idc = Ideal{P}(P.([x for x in idcQ.base if multideg(x)[1] == 0]))
    @test idc == intersect(ida, idb)
end

@testset "derivatives" begin
    P = Z[:x, :y]
    x, y = generators(P)
    @test derive(x, (1, 0)) == x^0
    @test derive(x, (0, 1)) == 0
    @test derive(x*y^10, (1, 1)) == 10*y^9
    @test derive((x+2y)^4, (0, 0)) == (x+2y)^4
    @test_throws ArgumentError derive(x*y, (-1, 0))
end


text:
function iscastreamfunction(init::AbstractDict,iroot::AbstractDict)

    for irun = 1 : init["nruns"]
        @info "$(Dates.now()) - Calculating Zonal-Averaged MERIDIONAL STREAMFUNCTION for RUN $irun"
        iscacalcpsi(init,iroot,irun,init["sealp"])
    end

end

function iscacalcpsi(
    init::AbstractDict,iroot::AbstractDict,
    irun::Integer,sealp::Real
)

    @info "$(Dates.now()) - Extracting MERIDIONAL WIND and PRESSURE data ..."
    inc = iscarawname(iroot,irun=irun); ids = Dataset(inc);
    vwind = ids["vcomp"][:]*1
    pfull = ids["pfull"][:]*1
    phsfc = ids["phalf"][end]*1
    pfull = vcat(0,pfull/phsfc)
    psfc  = ids["ps"][:]*1
    lat = ids["lat"][:]*1;
    close(ids)

    @info "$(Dates.now()) - Reshape / Permutedims data for calculation ..."
    nlon,nlat,npre,nt = size(vwind);
    vwind = permutedims(dropdims(mean(vwind,dims=1),dims=1),[1,3,2]);
    psfc  = dropdims(mean(psfc,dims=1),dims=1)
    vpsi = zeros(nlat,nt,npre); vii = zeros(npre+1); psiii = zeros(npre+1)

    @info "$(Dates.now()) - Performing Numerical Integration ..."
    for it = 1 : nt, ilat = 1 : nlat

        for ipre = 1 : npre; vii[ipre+1] = vwind[ilat,it,ipre] end
        psiii .= cumul_integrate(pfull*psfc[ilat,it],vii)
        for ipre = 1 : npre
            vpsi[ilat,it,ipre] = psiii[ipre+1] * 2 * pi * 6378e3 .* cosd.(lat[ilat]) / 9.81
        end

    end

    @info "$(Dates.now()) - Saving MERIDIONAL STREAMFUNCTION data for RUN $irun ..."
    vpsi = permutedims(vpsi,[1,3,2]); iscasavepsi(vpsi,init,iroot,irun)

end

function iscasavepsi(
    vpsi::Array{<:Real,3},
    init::AbstractDict,iroot::AbstractDict,irun::Integer
)

    imod,ipar,itime = iscainitialize(init,modID="cpre",parID="psi_v");
    inc = iscacalcname(ipar,iroot,irun=irun);
    if isfile(inc)
        @info "$(Dates.now()) - Stale NetCDF file $(inc) detected.  Overwriting ..."
        rm(inc);
    end

    @debug "$(Dates.now()) - Creating NetCDF file $(afnc) for MERIDIONAL STREAMFUNCTION for RUN $irun ..."
    ds = Dataset(inc,"c");

    ds.dim["longitude"] = length(imod["lon"])
    ds.dim["latitude"] = length(imod["lat"])
    ds.dim["time"] = Inf
    ds.dim["phalf"] = length(imod["phalf"])
    ds.dim["pfull"] = length(imod["pfull"])

    nclon = defVar(ds,"lon", Float64, ("longitude",), attrib = Dict(
        "long_name"      => "longitude",
        "units"          => "degrees_E",
        "cartesian_axis" => "X",
        "edges"          => "lonb",
    ))

    nclat = defVar(ds,"lat", Float64, ("latitude",), attrib = Dict(
        "long_name"      => "latitude",
        "units"          => "degrees_N",
        "cartesian_axis" => "Y",
        "edges"          => "latb",
    ))

    nctime = defVar(ds,"time", Float64, ("time",), attrib = Dict(
        "long_name"      => itime["ncattribs"]["long_name"],
        "units"          => itime["ncattribs"]["units"],
        "cartesian_axis" => itime["ncattribs"]["cartesian_axis"],
        "calendar_type"  => itime["ncattribs"]["calendar_type"],
        "calendar"       => itime["ncattribs"]["calendar"],
        "bounds"         => itime["ncattribs"]["bounds"],
    ))

    ncphalf = defVar(ds,"phalf", Float64, ("phalf",), attrib = Dict(
        "long_name"      => "approx half pressure level",
        "units"          => "hPa",
        "cartesian_axis" => "Z",
        "positive"       => "down",
    ))

    ncpfull = defVar(ds,"pfull", Float64, ("pfull",), attrib = Dict(
        "long_name"      => "approx full pressure level",
        "units"          => "hPa",
        "cartesian_axis" => "Z",
        "positive"       => "down",
    ))

    nclon[:] = imod["lon"]
    nclat[:] = imod["lat"]
    nctime[:] = itime["raw"]
    ncphalf[:] = imod["phalf"]
    ncpfull[:] = imod["pfull"]

    v = defVar(ds,"psi_v",Float32,("latitude","pfull","time"),attrib = Dict(
        "long_name"     => "Meridional Streamfunction",
        "units"         => "kg/s",
        "scale_factor"  => 1e9,
        "cell_methods"  => "time: mean",
        "time_avg_info" => "average_T1,average_T2,average_DT",
    ))

    v[:] = vpsi; close(ds)

    @info "$(Dates.now()) - Saved MERIDIONAL STREAMFUNCTION data for RUN $irun into NetCDF file."

end