text:
export Polygon

#= Polygon (for a base shape) =#

# Assume the followings for the polygon represented by Polygon:
# - The polygon is convex.
# - The vertices are listed in the counter-clockwise order around the origin.
mutable struct Polygon{K,K2} <: Shape{2,4}  # K2 = 2K
    v::SMatrix{2,K,Float64,K2}  # vertices
    n::SMatrix{2,K,Float64,K2}  # direction normals to edges
    Polygon{K,K2}(v,n) where {K,K2} = new(v,n)  # suppress default outer constructor
end

function Polygon(v::SMatrix{2,K,<:Real}) where {K}
    # Sort the vertices in the counter-clockwise direction
    w = v .- mean(v, dims=Val(2))  # v in center-of-mass coordinates
    ϕ = mod.(atan.(w[2,:], w[1,:]), 2π)  # SVector{K}: angle of vertices between 0 and 2π; `%` does not work for negative angle
    if !issorted(ϕ)
        # Do this only when ϕ is not sorted, because the following uses allocations.
        ind = MVector{K}(sortperm(ϕ))  # sortperm(::SVector) currently returns Vector, not MVector
        v = v[:,ind]  # SVector{K}: sorted v
    end

    # Calculate the increases in angle between neighboring edges.
    ∆v = hcat(diff(v, dims=Val(2)), SMatrix{2,1}(v[:,1]-v[:,end]))  # SMatrix{2,K}: edge directions
    ∆z = ∆v[1,:] + im * ∆v[2,:]  # SVector{K}: edge directions as complex numbers
    icurr = ntuple(identity, Val(K-1))
    inext = ntuple(x->x+1, Val(K-1))
    ∆ϕ = angle.(∆z[SVector(inext)] ./ ∆z[SVector(icurr)])  # angle returns value between -π and π

    # Check all the angle increases are positive.  If they aren't, the polygon is not convex.
    all(∆ϕ .> 0) || throw("v = $v should represent vertices of convex polygon.")

    n = [∆v[2,:] -∆v[1,:]]'  # SMatrix{2,K}; outward normal directions to edges
    n = n ./ hypot.(n[1,:], n[2,:])'  # normalize

    return Polygon{K,2K}(v,n)
end

Polygon(v::AbstractMatrix{<:Real}) = (K = size(v,2); Polygon(SMatrix{2,K}(v)))

# Regular polygon
function Polygon{K}(c::SVector{2,<:Real},
                    r::Real,  # distance between center and each vertex
                    θ::Real=0.0  # angle from +y-direction towards first vertex
                    ) where {K}
    ∆θ = 2π / K

    θs = π/2 + θ .+ ∆θ .* SVector(ntuple(k->k-1, Val(K)))  # SVector{K}: angles of vertices
    v = c .+ r .* [cos.(θs) sin.(θs)]'  # SMatrix{2,K}: locations of vertices

    return Polygon(v)
end

Polygon{K}(c::AbstractVector{<:Real},  # [x, y]: center of regular polygon
           r::Real,  # radius: distance from center to vertices
           θ::Real=0.0  # angle of first vertex
           ) where {K} =
   Polygon{K}(SVector{2}(c), r, θ)

Base.:(==)(s1::Polygon, s2::Polygon) = s1.v==s2.v && s1.n==s2.n  # assume sorted v
Base.isapprox(s1::Polygon, s2::Polygon) = s1.v≈s2.v && s1.n≈s2.n  # assume sorted v
Base.hash(s::Polygon, h::UInt) = hash(s.v, hash(s.n, hash(:Polygon, h)))

function level(x::SVector{2,<:Real}, s::Polygon)
    c = mean(s.v, dims=Val(2))  # center of mass

    d = sum(s.n .* (x .- c), dims=Val(1))
    r = sum(s.n .* (s.v .- c), dims=Val(1))
    @assert all(r .> 0)

    return maximum(d ./ r) - 1.0
end

function surfpt_nearby(x::SVector{2,<:Real}, s::Polygon{K}) where {K}
    # Calculate the signed distances from x to edge lines.
    ∆xe = sum(s.n .* (x .- s.v), dims=Val(1))[1,:]  # SVector{K}: values of equations of edge lines
    abs∆xe = abs.(∆xe)  # SVector{K}

    # Determine if x is outside of edges, inclusive.
    sz = abs.((-)(bounds(s)...))  # SVector{2}
    onbnd = abs∆xe .≤ Base.rtoldefault(Float64) * max(sz.data...)  # SVector{K}
    isout = (∆xe.>0) .| onbnd  # SVector{K}

    # For x inside the polygon, it is easy to find the closest surface point: we can simply
    # pick the closest edge and find its point closest to x.
    # For x outside the polygon, there are many cases depending on how many edges x lies
    # outside of.  However, x that is sufficiently close to the polygon should be outside of
    # at most two edges.  (Outside two edges near the vertices, and one edge elsewhere.)
    # Therefore, we will consider only cout ≤ 2 below.
    cout = count(isout)
    if cout == 2  # x is outside two edges
        # We could choose to find ind corresponding to the two nonnegative ∆xe, but such an
        # operation leads to allocations because it does not create an SVector.  Instead,
        # find the closest vertex directly, which is the surface point we are looking for
        # for cout = 2.
        ∆xv = x .- s.v
        l∆xv = hypot.(∆xv[1,:], ∆xv[2,:])
        imin = argmin(l∆xv)
        surf = s.v[:,imin]
        imin₋₁ = mod1(imin-1,K)

        if onbnd[imin] && onbnd[imin₋₁]  # x is very close to vertex imin
            nout = s.n[:,imin] + s.n[:,imin₋₁]
        else
            nout = x - s.v[:,imin]
        end
        nout = normalize(nout)
    else  # cout ≤ 1 or cout ≥ 3
        cout ≤ 1 || @warn "x = $x is outside $cout edges: too far from polygon with vertices $(s.v); " *
                            "result could be inaccurate."
        # Choose the closest edge to x.
        # If cout = 0, all ∆xe are negative, so the largest ∆xe is the smallest in magnitude
        # and corresponds to the closest edge.
        # If cout = 1, all but one ∆xe are negative, so again the largest ∆xe is the only
        # nonnegative one and corresponds to the edge outside which x lies.
        # Even for cout = 3, this seems to be a reasonable choice from a simple geometric
        # argument.
        imax = argmax(∆xe)
        vmax, nmax = s.v[:,imax], s.n[:,imax]

        ∆x = (nmax⋅(vmax-x)) .* nmax
        surf = x + ∆x
        nout = nmax
    end

    return surf, nout
end

translate(s::Polygon, ∆::SVector{2,<:Real}) = (s2 = deepcopy(s); s2.v = s2.v .+ ∆; s2)

function bounds(s::Polygon)
    l = minimum(s.v, dims=Val(2))[:,1]
    u = maximum(s.v, dims=Val(2))[:,1]

    return (l, u)
end


text:
module AirSeaFluxes

using Dierckx
using LinearAlgebra

include("simpleflux.jl")
include("bulkformulae.jl")
include("holtslag.jl")
include("kpp.jl")
include("AOGCM1D_helpers.jl")
include("AOGCM1D.jl")
include("calc_pCO2.jl")

export simpleflux, bulkformulae, AOGCM1D, calc_pco2

end # module

ext: jl


text:
using Documenter, ShiftedArrays

DocMeta.setdocmeta!(ShiftedArrays, :DocTestSetup, :(using ShiftedArrays); recursive=true)

makedocs(
    # options
    modules = [ShiftedArrays],
    sitename = "ShiftedArrays.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
    ),
    pages = Any[
        "Introduction" => "index.md",
        "API" => "api.md",
    ],
    strict = true,
)

# Deploy built documentation from Travis.
# =======================================

deploydocs(
    # options
    repo = "github.com/JuliaArrays/ShiftedArrays.jl.git",
    push_preview = true,
)


text:
function generate_reaction_string(reaction_string::String)

    # split the reaction into its components -
	component_array = split(reaction_string,',');

    # reaction name -
	reaction_name = String.(component_array[1]);
	reactant_phrase = String.(component_array[2]);
	product_phrase = String.(component_array[3]);

    # reaction string
    return "$(reactant_phrase) -> $(product_phrase)"
end

function generate_reaction_strings(reactions::Array{String,1}; expand::Bool = false)::Array{String,1}

    # initialize -
    reaction_string_array = Array{String,1}()

    # should we expand the reversible reactions?
    reactions_to_process = reactions;
    if (expand == true)
        reactions_to_process = expand_reversible_reactions(reactions);
    end

    # first: let's discover the species list -
	for reaction_string ∈ reactions_to_process

        # rstring -
        rstring = generate_reaction_string(reaction_string)
        push!(reaction_string_array, rstring)
    end

    # return -
    return reaction_string_array
end


text:
# ## -------------------------------------------------------------------
# # MSE per method
# let

#     for (iJR, Data) in [
#             (ChK.iJR904, ChK.KayserData),
#             (ChN.iJR904, ChN.NanchenData),
#             (ChF.iJR904, ChF.FolsomData),
#             (ChH.iJR904, ChH.HeerdenData),
#         ]

#         model = iJR.load_model("max_model")
#         datfile = iJR.procdir("dat.bson")
#         DAT = UJL.load_data(datfile; verbose = false)

#         EXCH_FLX_IDERS = DAT[:EXCH_FLX_IDERS]
#         EXPS = DAT[:EXPS]

#         p = plot(;xlabel = "experiment", ylabel = "MSE")
#         for method in [:ME_MAX_POL]
#             MSEs = []
#             for exp in EXPS

#                 sum = 0.0
#                 N = 0
#                 glc_flx = DAT[:exp, :flx, "GLC", exp]
#                 for ider in EXCH_FLX_IDERS
#                     model_val = DAT[method, :flx, ider, exp]
#                     exp_val = DAT[:exp, :flx, ider, exp]
#                     sum += (model_val/glc_flx - exp_val/glc_flx)^2
#                     N += 1
#                 end
#                 push!(MSEs, sum / N)
#             end
#             scatter!(p, EXPS, MSEs; color = method_colors[method],
#                 label = string(method), m = 8, alpha = 0.8,
#                 legend = :topleft
#             )
#             plot!(p, EXPS, MSEs; color = method_colors[method],
#                 label = "", ls = :dash, alpha = 0.8
#             )
#         end
#         src = nameof(Data)
#         mysavefig(p, "MSE_per_method"; src)
#     end
# end

## -------------------------------------------------------------------
let
    me_method = :ME_MAX_POL
    fba_method = :FBA_Z_FIX_MIN_COST

    p0 = plot(;xlabel = "politope box volume (log)", ylabel = "MSE")
    all_p = deepcopy(p0)
    zeroth = 1e-2
    for (iJR, Data) in [
                (ChK.iJR904, ChK.KayserData),
                (ChN.iJR904, ChN.NanchenData),
                (ChF.iJR904, ChF.FolsomData),
                (ChH.iJR904, ChH.HeerdenData),
            ]

        src = nameof(Data)

        model = iJR.load_model("max_model")
        datfile = iJR.procdir("dat.bson")
        DAT = UJL.load_data(datfile; verbose = false)

        EXCH_FLX_IDERS = ["GLC", "AC"]
        EXPS = DAT[:EXPS]
        rxns_map = iJR.load_rxns_map()

        exglcidx = ChU.rxnindex(model, iJR.EX_GLC_IDER)
        biomidx = ChU.rxnindex(model, iJR.BIOMASS_IDER)

        model_ids = [rxns_map[id] for id in EXCH_FLX_IDERS]
        model_idxs = [ChU.rxnindex(model, model_id) for model_id in model_ids]

        # Reduce Pol
        p = deepcopy(p0)
        vols = []
        me_MSEs = []
        fba_MSEs = []

        # ---------------------------------------------------
        # ME_MAX_POL
        for exp in EXPS
            # box vol
            cgD_X = -Data.ciD_X(:GLC, exp)
            D = Data.val(:D, exp)

            ChU.lb!(model, exglcidx, cgD_X)
            ChU.bounds!(model, biomidx, D, D)
            vol = box_vol(model, model_idxs) .+ zeroth
            isnan(vol) && continue

            push!(vols, vol)

            # MSE
            exp_glc = DAT[:exp, :flx, "GLC", exp]
            exp_flxs = DAT[:exp, :flx, EXCH_FLX_IDERS, exp]

            # ME
            me_flxs = DAT[me_method, :flx, EXCH_FLX_IDERS, exp]
            MSEsum = sum(((exp_flxs ./ exp_glc) .- (me_flxs ./ exp_glc)).^2)
            MSE = MSEsum / length(EXCH_FLX_IDERS)
            push!(me_MSEs, MSE)

            any(MSE .> 0.3) && let
                println(); @info("Outlayer", exp, src); println()
            end

            # FBA
            fba_flxs = DAT[fba_method, :flx, EXCH_FLX_IDERS, exp]
            MSEsum = sum(((exp_flxs ./ exp_glc) .- (fba_flxs ./ exp_glc)).^2)
            MSE = MSEsum / length(EXCH_FLX_IDERS)
            push!(fba_MSEs, MSE)

        end
        sidxs = sortperm(vols)

        marker = source_markers[Data]
        for _p in [p, all_p]

            scatter!(_p, vols[sidxs], me_MSEs[sidxs];
                label = "", m = (8, marker),
                color = method_colors[me_method]
            )

            scatter!(_p, vols[sidxs], fba_MSEs[sidxs];
                label = "", m = (8, marker),
                color = method_colors[fba_method]
            )
        end

        mysavefig(p, "pol_box_vs_MSE"; src)
    end
    plot!(all_p; xscale = :log10)
    mysavefig(all_p, "pol_box_vs_MSE")
end

## -------------------------------------------------------------------
# TODO: work better inter project comunication
# MSE FBA vs ME
let

    DAT_FILE_PREFFIX =  "maxent_ep_dat"
    me_method = :ME_MAX_POL
    fba_method = :FBA_Z_FIX_MIN_COST

    p0 = plot(;xlabel = "experiment", ylabel = "MSE [FBA - ME]",)
    for (iJR, Data) in [
        (ChK.iJR904, ChK.KayserData),
        (ChN.iJR904, ChN.NanchenData),
        (ChF.iJR904, ChF.FolsomData),
        (ChH.iJR904, ChH.HeerdenData),
    ]
        src = nameof(Data)
        marker = (8, source_markers[Data])

        datfile = iJR.procdir("dat.bson")
        DAT = UJL.load_data(datfile; verbose = false)

        # Preload FBA Data
        LP_DAT = ChE.load_LP_DAT(src)

        EXPS = DAT[:EXPS]
        MSEs = []
        p = deepcopy(p0)
        for exp in EXPS

            # norm
            exp_glc = DAT[:exp, :flx, "GLC", exp]

            # load ME dat
            ME_DAT = UJL.mysavename(DAT_FILE_PREFFIX, "jls";
                method = me_method, exp
            ) |> iJR.procdir |> deserialize

            # common rxns
            me_model = ME_DAT[:model]
            fba_model = LP_DAT[fba_method, :model, exp]
            RXNS = intersect(me_model.rxns, fba_model.rxns)

            # ME flxs
            epouts = ME_DAT[:epouts]
            exp_beta = maximum(keys(epouts))
            epout = epouts[exp_beta]
            me_flxs = ChU.av.([me_model], [epout], RXNS)
            @show length(me_flxs)

            # FBA
            fbaout = LP_DAT[fba_method, :fbaout, exp]
            fba_flxs = ChU.av.([fba_model], [fbaout], RXNS)
            @show length(fba_flxs)
            println()

            # MSE
            MSEsum = sum(((me_flxs ./ exp_glc) .- (fba_flxs ./ exp_glc)).^2)
            MSE = MSEsum / length(fba_flxs)
            push!(MSEs, MSE)

        end
        scatter!(p, EXPS, MSEs; label = "",
            alpha = 0.8, color = :black, m = 6
        )
        plot!(p, EXPS, MSEs; label = "",
            alpha = 0.6, ls = :dash, lw = 3, color = :black
        )
        mysavefig(p, "MSE_fba_vs_me"; src)
    end
end

## -------------------------------------------------------------------
# # MSE per ider
# let
#     p = plot(;xlabel = "experiment", ylabel = "MSE")
#     for method in ALL_METHODS
#         MSEs = []

#         for ider in EXCH_FLX_IDERS
#             sum = 0.0
#             N = 0
#             for exp in EXPS
#                 glc_flx = DAT[method, :Fd, :flx, "GLC", exp]
#                 model_val = DAT[method, :ep, :flx, ider, exp]
#                 exp_val = DAT[method, :Fd, :flx, ider, exp]
#                 sum += (model_val/glc_flx - exp_val/glc_flx)^2
#                 N += 1
#             end
#             push!(MSEs, sum / N)
#         end

#         scatter!(p, EXCH_FLX_IDERS, MSEs; color = method_colors[method],
#             label = string(method), m = 8, alpha = 0.8,
#             legend = :topleft
#         )
#         plot!(p, EXCH_FLX_IDERS, MSEs; color = method_colors[method],
#             label = "", ls = :dash, alpha = 0.8
#         )
#     end
#     mysavefig(p, "MSE_per_ider")
# end

# ## -------------------------------------------------------------------
# # MSE per beta
# let
#     method = ME_MAX_POL

#     ps = Plots.Plot[]
#     for exp in EXPS
#         p = plot(;title = string("exp: ", exp),
#             xlabel = "beta", ylabel = "MSE"
#         )

#         datfile = INDEX[method, :DFILE, exp]
#         dat = deserialize(datfile)
#         epouts = dat[:epouts]
#         betas = epouts |> keys |> collect |> sort
#         exp_beta = maximum(betas) # dat[:exp_beta]
#         model = dat[:model]

#         MSEs = []
#         for beta in betas
#             epout = epouts[beta]

#             sum = 0.0
#             N = 0

#             glc_flx = Fd.uval(:GLC, exp)
#             for ider in EXCH_FLX_IDERS

#                 model_met = Fd_mets_map[ider]
#                 model_exch = Fd_rxns_map[ider]
#                 model_exchi = ChU.rxnindex(model, model_exch)

#                 model_flx = ChU.av(model, epout, model_exchi)
#                 exp_flx = Fd.uval(ider, exp)

#                 sum += (model_flx/glc_flx - exp_flx/glc_flx)^2
#                 N += 1
#             end

#             push!(MSEs, sum / N)
#         end

#         scatter!(p, first.(betas), MSEs; color = :black,
#             label = "", m = 8, alpha = 0.8
#         )
#         plot!(p, first.(betas), MSEs; color = :black,
#             label = "", ls = :dash, alpha = 0.8
#         )
#         vline!(p, [first(exp_beta)]; color = :black,
#             label = "", ls = :dot, lw = 3, alpha = 0.9
#         )
#         push!(ps, p)
#     end
#     mysavefig(ps, "MSE_vs_beta")
# end


text:
# Realization and dynamic range scaling
using RSDeltaSigmaPort
using RSDeltaSigmaPort.EasyPlot #set, cons
import Printf: @sprintf
j=im

function stringify(vA, digits=6)
	#ctx=IOContext(stdout, :compact=>true)
	#IOContext does not support # of digits
	#_s(v) = @sprintf("%.6f", v) #Cannot parameterize format specifier
	return [round(v, digits=6) for v in vA]
end

#==
===============================================================================#
println("\n*** Modulator realization and scaling")
OSR=42
order=5; opt=1

NTF = synthesizeNTF(order, OSR, opt=opt)
plot = plotNTF(NTF, OSR)
plot.title = string(order, "th-Order Modulator")
#saveimage(:png, "dsdemo3_ntf.png", plot, AR=2/1, width=900)
displaygui(plot)

@info("Realizing NTF..."); flush(stdout); flush(stderr)
#-------------------------------------------------------------------------------
a,g,b,c = realizeNTF(NTF)
b = [b[1] zeros(1,length(b)-1)] #Use a single feed-in for the input

println("\nUnscaled modulator")
println("   DAC feedback coefficients = ", join(stringify(a), ", "))
println("   resonator feedback coefficients = ", join(stringify(g), ", "))

println()
@info("Calculating the state maxima...")
#-------------------------------------------------------------------------------
ABCD = stuffABCD(a,g,b,c)
u = range(0, stop=0.6, length=30)
N = 10^4
plot = plotStateMaxima(u, ABCD, N=N)
displaygui(plot)

println()
@info("Calculating the scaled coefficients...")

#-------------------------------------------------------------------------------
ABCDs, umax = scaleABCD(ABCD, N=N)
as, gs, bs, cs = mapABCD(ABCDs)
println("\nScaled modulator, umax=", @sprintf("%.2f", umax))
println("   DAC feedback coefficients = ", join(stringify(as), ", "))
println("   resonator feedback coefficients = ", join(stringify(gs), ", "))
println("   interstage coefficients = ", join(stringify(cs), ", "))
println("   feed-in coefficients = ", join(stringify(bs), ", "))

println()
@info("Calculating the state maxima...")
#-------------------------------------------------------------------------------
u = range(0, stop=umax, length=30)
N = 10^4

plot = plotStateMaxima(u, ABCDs, N=N)
displaygui(plot)

:END_OF_DEMO


text:
type SysexEvent <: TrackEvent
    dT::Int
    data::Array{UInt8,1}
end

function issysexevent(b::UInt8)
    b == 0xF0
end

function readsysexevent(dT::UInt8, f::IO)
    data = UInt8[]
    read(f, UInt8) # Eat the SYSEX that's on top of f
    datalength = readvariablelength(f)
    b = read(f, UInt8)
    while isdatabyte(b)
        push!(data, b)
        b = read(f, UInt8)
    end

    if b != 0xF7
        error("Invalid sysex event, did not end with 0xF7")
    end
    # The last byte of sysex event is F7. We leave that out of the data, and write it back when we write the event.

    if length(data) + 1 != datalength
        error("Invalid sysex event. Expected $(datalength) bytes, received $(length(data) + 1)")
    end
    SysexEvent(dT, data)
end

function writeevent(f::IO, event::SysexEvent)
    writevariablelength(f, event.dT)
    write(f, SYSEX)
    writevariablelength(f, length(event.data) + 1) # +1 for the ending F7
    write(f, event.data)
    write(f, 0xF7)
end

export SysexEvent


text:
alpine1(x) = sum(abs(x .* sin(x) + 0.1 .* x))
alpine2(x) = prod(sqrt(x) .* sin(x))
chung_reynolds(x) = sum(x .^ 2) ^ 2
qing(x) = sum((x .^ 2 - (1:length(x))) .^ 2)


text:
using Literate

EXAMPLEDIR = joinpath(@__DIR__, "src", "examples")
GENERATEDDIR = joinpath(@__DIR__, "src", "examples", "generated")

for example in filter!(x -> endswith(x, ".jl"), readdir(EXAMPLEDIR))
    input = abspath(joinpath(EXAMPLEDIR, example))
    Literate.markdown(input, GENERATEDDIR)
end


text:
##############################################################################################
module DataAssimilationBenchmarks
##############################################################################################
##############################################################################################
# imports and exports
include("methods/DeSolvers.jl")
include("methods/EnsembleKalmanSchemes.jl")
include("models/L96.jl")
include("models/IEEE39bus.jl")
include("experiments/GenerateTimeSeries.jl")
include("experiments/FilterExps.jl")
include("experiments/SmootherExps.jl")
include("experiments/SingleExperimentDriver.jl")
using .DeSolvers
using .EnsembleKalmanSchemes
using .L96
using .IEEE39bus
using .GenerateTimeSeries
using .FilterExps
using .SingleExperimentDriver
export DeSolvers, EnsembleKalmanSchemes, L96, IEEE39bus, GenerateTimeSeries, FilterExps,
       SingleExperimentDriver

##############################################################################################
##############################################################################################

function Info()
    print("  _____        _                         ")
    printstyled("_",color=9)
    print("           ")
    printstyled("_",color=2)
    print(" _       _   ")
    printstyled("_              \n",color=13)
    print(" |  __ \\      | |          /\\           ")
    printstyled("(_)",color=9)
    printstyled("         (_)",color=2)
    print(" |     | | ")
    printstyled("(_)             \n",color=13)
    print(" | |  | | __ _| |_ __ _   /  \\   ___ ___ _ _ __ ___  _| | __ _| |_ _  ___  _ __   \n")
    print(" | |  | |/ _` | __/ _` | / /\\ \\ / __/ __| | '_ ` _ \\| | |/ _` | __| |/ _ \\| '_ \\  \n")
    print(" | |__| | (_| | || (_| |/ ____ \\\\__ \\__ \\ | | | | | | | | (_| | |_| | (_) | | | | \n")
    print(" |_____/ \\__,_|\\__\\__,_/_/    \\_\\___/___/_|_| |_| |_|_|_|\\__,_|\\__|_|\\___/|_| |_| \n")
    print("\n")
    print("  ____                  _                          _         ")
    printstyled(" _ ", color=12)
    print("_\n")
    print(" |  _ \\                | |                        | |        ")
    printstyled("(_)",color=12)
    print(" |                \n")
    print(" | |_) | ___ _ __   ___| |__  _ __ ___   __ _ _ __| | _____   _| |                \n")
    print(" |  _ < / _ \\ '_ \\ / __| '_ \\| '_ ` _ \\ / _` | '__| |/ / __| | | |                \n")
    print(" | |_) |  __/ | | | (__| | | | | | | | | (_| | |  |   <\\__ \\_| | |                \n")
    print(" |____/ \\___|_| |_|\\___|_| |_|_| |_| |_|\\__,_|_|  |_|\\_\\___(_) |_|                \n")
    print("                                                            _/ |                  \n")
    print("                                                           |__/                   \n")

    print("\n")
    printstyled(" Welcome to DataAssimilationBenchmarks!\n", bold=true)
    print(" Version 0.10, Copyright 2021 Colin James Grudzien (cgrudz@mailbox.org)\n")
    print(" Licensed under the Apache License, Version 2.0 \n")
    print(" https://github.com/cgrudz/DataAssimilationBenchmarks/blob/master/LICENSE.md\n")
    print("\n")
    printstyled(" Description\n", bold=true)
    print(" This is my personal data assimilation benchmark research code with an emphasis on testing and validation\n")
    print(" of ensemble-based filters and smoothers in chaotic toy models. DataAssimilationBenchmarks is a wrapper\n")
    print(" library including the core numerical solvers for ordinary and stochastic differential  equations,\n")
    print(" solvers for data assimilation routines and the core process model code for running twin experiments\n")
    print(" with benchmark models. These methods can be run stand-alone in other programs by calling these\n")
    print(" functions from the DeSolvers, EnsembleKalmanSchemes and L96 sub-modules from this library.\n")
    print(" Future solvers and models will be added as sub-modules in the methods and models directories respectively.\n")
    print("\n")
    print(" In order to get the full functionality of this package you you will need to install the dev version.\n")
    print(" This provides the access to edit all of the outer-loop routines for setting up twin experiments. \n")
    print(" These routines are defined in the modules in the \"experiments\" directory.  The \"slurm_submit_scripts\"\n")
    print(" directory includes routines for parallel submission of experiments in Slurm.  Data processing scripts\n")
    print(" and visualization scripts (written in Python with Matplotlib and Seaborn) are included in the \"analysis\"\n")
    print(" directory. \n")
    print(" \n")
    print(" Instructions on how to install the dev version of this package are included in the README.md\n")

    nothing

end

##############################################################################################

end

max_stars_repo_name: sukhreens/DataAssimilationBenchmarks

max_issues_repo_name: sukhreens/DataAssimilationBenchmarks

max_forks_repo_name: sukhreens/DataAssimilationBenchmarks


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using RadMod2D

include("./models2D.jl")
include("./plot2D.jl")

using CairoMakie
CairoMakie.activate!(type = "svg")
filetype = ".svg"
# CairoMakie.activate!(type = "png")
# filetype = ".png"

function fig_pinhole_therm2D_eps(eps_lab)
    # two rectangles connected with pinhole
    m = model_two_rectangles_with_holes(0.02)
    # blocking with tilewalk
    vfmat = zeros(Float64, m.nelem, m.nelem)
    existing_vf!(m, vfmat)
    n = 30
    dx, dy = create_tiles(m, n)
    @time t_occ = check_tile_occupation(m, dx, dy, n)
    @time blocking_vf_with_tiles!(m, vfmat, dx, dy, n, t_occ)
    @time calculating_vf!(m, vfmat, normit = true)
    # solve Qp
    temp = zeros(m.nelem,1)
    set_bc_part!(m, temp, 1:2, 400)
    set_bc_part!(m, temp, 3, 600)
    set_bc_part!(m, temp, 4:5, 400)
    # set_bc_part!(m, temp, 1:5, 600)
    set_bc_part!(m, temp, 6:10, 400)
    epsilon = zeros(m.nelem,1)
    set_bc_part!(m, epsilon, 1:10, 1.0)
    set_bc_part!(m, epsilon, 8, eps_lab)
    @time Qp, G = tempsolver(m, vfmat, temp, epsilon)
    Qp_parts = [sum(Qp[m.elem2par[i].first:m.elem2par[i].last,1]) for i = 1:m.npar]
    # rect2_min = minimum(Qp[m.elem2par[6].first:m.elem2par[10].last,1])
    # rect2_max = maximum(Qp[m.elem2par[6].first:m.elem2par[10].last,1])
    # println("Qp of rect2 between :", rect2_min, " and ", rect2_max)
    area = [m.elem[i].area for i = m.elem2par[1].first:m.elem2par[end].last]
    qp_area = Qp[:] ./ area[:]
    qp_area_kw = qp_area ./ 1000
    rect2_min = minimum(qp_area_kw[m.elem2par[6].first:m.elem2par[10].last,1])
    rect2_max = maximum(qp_area_kw[m.elem2par[6].first:m.elem2par[10].last,1])
    println("qp_area_kw of rect2 between :", rect2_min, " and ", rect2_max)
    # modifying qw for displaying only incoming heat flux on rect2
    qp_area_kw_mod = qp_area_kw
    qp_area_kw_mod[qp_area_kw.>=rect2_max] .= 0.0
    qp_area_kw_mod[qp_area_kw.<=rect2_min] .= rect2_min
    qp_area_kw_mod .*= (-1)
    # 2D plot
    fig = Figure(resolution = (270, 400), font = "Arial", fontsize = 10)
    ax = fig[1, 1] = Axis(fig)
    linewidth = 1.0
    setup_axis(ax, linewidth)
    ax.xlabel = "X in m"
    ax.ylabel = "Y in m"
    ax.aspect = DataAspect()
    cbar = plot_model_with_value(fig, ax, m, qp_area_kw_mod, "Incoming heat flux density in kW/m", cmap = :viridis, linewidth = linewidth*1.5, showcbar = true)
    # cropping here at 0 and only displaying incoming heat flux - also making all positive
    # cbar.ticks = -50:20:50
    # cbar.limits = (-10,0.1)
    cbar.limits = (0, 0.2)
    cbar.ticks = 0.0:0.05:0.2
    cbar.vertical = false
    cbar.width = 120
    cbar.height = 12
    cbar.labelpadding = 1.0 # default 5.0
    cbar.ticklabelpad = 1.0 # default 3.0
    cbar.ticksize = 2.0 # default 6.0
    cbar.tellwidth = false # erlaubt es das colorbar und plot andere groeßen zueinander haben
    # display(fig)
    eps_save = replace(string(eps_lab),"." => "n")
    save("fig_c2_therm2D_eps"*eps_save*"_cb"*filetype, fig)
    # for i = 1:m.nnodes
    #     println(m.nodes[i].x, ";", m.nodes[i].y)
    # end
    # for i = 1:m.nelem
    #     println(m.elem[i].node1, ";", m.elem[i].node2, ";", m.elem[i].com, ";", m.elem[i].nvec, ";", m.elem[i].area)
    # end
    # for i = 1:m.nelem
    #     println(temp[i,:], ";", epsilon[i,:], ";", qp_area_kw_mod[i,:])
    # end
end

# fig_pinhole_therm2D_eps(0.5)
fig_pinhole_therm2D_eps(1.0)

# include("./test/run_pinhole.jl")


text:
@testset "ExponentialLikelihood" begin
    lik = ExponentialLikelihood()
    test_interface(lik, SqExponentialKernel(), rand(10))
end


text:
"""
    semidiscretize(u0func, semi::AbstractSemidiscretization, tspan)

Apply the semidiscretization `semi` to the initial data given by `u0func`
and return an `ODEProblem` with time span `tspan`.
"""
function semidiscretize(u0func, semi::AbstractSemidiscretization, tspan)
    u0 = compute_coefficients(u0func, semi.derivative)
    ode = ODEProblem(semi, u0, tspan)
end

function evaluate_coefficients(u, semi::AbstractSemidiscretization)
    evaluate_coefficients(u, semi.derivative)
end

function integrate(func, u, semi::AbstractSemidiscretization)
    integrate(func, u, semi.derivative)
end

struct ScalarIntegralQuantities{T} <: FieldVector{2,T}
    mass::T
    energy::T
end

function FilterCallback(filter::AbstractFilter)
    condition = (t, u, integrator) -> true
    affect!(integrator) = filter(integrator.u)

    DiscreteCallback(condition, affect!, save_positions=(false,false))
end


text:
############################################################################################
for iter in eachindex(objectives)
    @testset "Likelihood estimate, all models" begin
        ## Resampling methods
        for _resample in resamplemethods
            _obj = deepcopy(objectives[iter])
            ## Sample data given model and create new objective
            dat, lat = simulate(_rng, _obj.model; Nsamples = 500)
            _tagged = Tagged(_obj.model, :latent)
            fill!(_obj.model, _tagged, (; latent = lat))
            _obj = Objective(_obj.model, dat, _tagged)

            ## Define PF default tuning parameter
            pfdefault = ParticleFilterDefault(;
                coverage = 2.00,
                threshold = 0.75,
                referencing = Marginal(),
                resampling = _resample
            )
            ## Initialize kernel and check if it can be run
            pfkernel = ParticleFilter(
                _rng,
                _obj,
                pfdefault
            )
            _val, _diag = propose(_rng, pfkernel, _obj)
            ## Check Marginal PF likelihood
            _type = infer(_rng, BaytesFilters.AbstractDiagnostics, pfkernel, _obj.model, _obj.data)
            diags = Vector{_type}(undef, 50)
            for idx in eachindex(diags)
                _, diags[idx] = propose(_rng, pfkernel, _obj)
            end
            ℓobjective_approx = [diags[idx].base.ℓobjective for idx in eachindex(diags)]
            _, ℓobjective_exact = filter_forward(_obj)
            @test ℓobjective_exact ≈ mean(ℓobjective_approx) atol = 50.0 #!NOTE: Initial distribution slightly different in forward filter for HSMM, rest should be atol ~ 1.0
            #Check for results
            results(diags, pfkernel, 2, [.1, .2, .5, .8, .9])
        end
    end
end

############################################################################################
# Number of Particle estimate
for iter in eachindex(objectives)
    @testset "Number Particle estimation, all models" begin
        ## Resampling methods
        for _resample in resamplemethods
            ## Referencing methods
            for _reference in referencemethods
                _obj = deepcopy(objectives[iter])
                ## Sample data given model and create new objective
                dat, lat = simulate(_rng, _obj.model; Nsamples = 500)
                _tagged = Tagged(_obj.model, :latent)
                fill!(_obj.model, _tagged, (; latent = lat))
                _obj = Objective(_obj.model, dat, _tagged)
                ## Define PF default tuning parameter
                pfdefault = ParticleFilterDefault(;
                    referencing = _reference,
                    resampling = _resample
                )
                ## Check if we can initiate from Constructor
                constructor = ParticleFilterConstructor(:latent, pfdefault)
                ## Check number of particles
                _ℓobjectivevariance = 0.50
                _maxiter = 10
                _margin = 1.0
                _NparticlesEstimate = estimate_Nparticles(
                    _rng, constructor, _obj, _ℓobjectivevariance;
                    Nchains = 4,
                    margin = _margin,
                    itermax = _maxiter,
                    mincoverage = 0.1,
                    printoutput = false
                )
                _variance = _NparticlesEstimate.ℓobjective_variance[_NparticlesEstimate.itermax]
                _boundaries = (_variance - _margin, _variance + _margin)
                if _NparticlesEstimate.itermax < _maxiter
                    @test _boundaries[begin] < _variance < _boundaries[end]
                end
            end
        end
    end
end


text:
using Metrics

@testset "Regression" begin
    y_true = [1, 1, 0, 0, 1]
    y_pred = [0.9, 0.9, 0.1, 0.2, 0.88]

    @test mae(y_pred, y_true) == 0.12399999999999997
    @test mse(y_pred, y_true) == 0.01688
    @test male(y_pred, y_true) == 13.70276763512634
    @test msle(y_pred, y_true) == 464.8426856611348

    @test r2_score(y_pred, y_true) == 0.9296666666666666
    @test adjusted_r2_score(y_pred, y_true, 3) == 0.7186666666666666
end


text:
# ------------------------------------------------------------------------------------------
# ## 函数
#
# 话题：
# 1. 如何声明一个函数
# 2. Julia 中的鸭子类型（Duck-typing）
# 3. 原地修改与非原地修改函数（Mutating vs. non-mutating functions）
# 4. 高阶函数
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### 如何声明一个函数
#
# Julia 提供了几种不同的声明函数的方法。
#
# #### 方法一：用 `function` 和 `end` 关键字
# ------------------------------------------------------------------------------------------

function sayhi(name)
    println("Hi $name, it's great to see you!")
end

function f(x)
    x^2
end

# ------------------------------------------------------------------------------------------
# 我们可以这样调用函数：
# ------------------------------------------------------------------------------------------

sayhi("C-3PO")

f(42)

# ------------------------------------------------------------------------------------------
# #### 方法二：用 `=`
# 等价的，我们可以只用一行声明函数，拿上面两个函数为例
# ------------------------------------------------------------------------------------------

sayhi2(name) = println("Hi $name, it's great to see you!")

f2(x) = x^2

sayhi2("R2D2")

f2(42)

# ------------------------------------------------------------------------------------------
# #### 方法三：匿名函数
# 最后，我们可以声明“匿名”函数
# ------------------------------------------------------------------------------------------

sayhi3 = name -> println("Hi $name, it's great to see you!")

f3 = x -> x^2

sayhi3("Chewbacca")

f3(42)

# ------------------------------------------------------------------------------------------
# ### Julia 中的鸭子类型（Duck-typing）
# > *“假如它嘎嘎像个鸭鸭，那它就是个鸭鸭”*
#
# 在未指定函数参数的类型时，不管输入参数是什么类型，只要函数体可以正常执行，
# 那么 Julia 函数就可以应用在这种类型的参数上
#
#
#
# 例如，函数 `sayhi` 可以接受这个写作整数的小电视角色的名字……
# ------------------------------------------------------------------------------------------

sayhi(55595472)

# ------------------------------------------------------------------------------------------
# 函数 `f` 可以接受一个矩阵
# ------------------------------------------------------------------------------------------

A = rand(3, 3)
A

f(A)

# ------------------------------------------------------------------------------------------
# 函数 `f` 也可以接受字符串比如"hi"，因为 `*` 在接受字符串时用作字符串拼接。
# ------------------------------------------------------------------------------------------

f("hi")

# ------------------------------------------------------------------------------------------
# 然而，`f` 函数无法接受一个向量（vector）。
# 不同于 `A^2` 是已经定义好的矩阵乘方，当 `v` 是向量时的 `v^2` 是未定义的，会导致错误。
# ------------------------------------------------------------------------------------------

v = rand(3)

f(v)

# ------------------------------------------------------------------------------------------
# ### 原地修改与非原地修改函数（Mutating vs. non-mutating functions）
#
# 约定：函数名以 `!` 结束的原地修改传入的变量，函数名非以 `!` 结束的则不会改变传入的变量。
#
# 例如下面的 `sort` 和 `sort!`
# ------------------------------------------------------------------------------------------

v = [3, 5, 2]

sort(v)

v

# ------------------------------------------------------------------------------------------
# `sort(v)`会返回一个新的、将 `v` 排序后的数组，`v` 本身不会改变
#
# 而当执行 `sort!(v)` 后，`v` 将会被原地排序修改
# ------------------------------------------------------------------------------------------

sort!(v)

v

# ------------------------------------------------------------------------------------------
# ### 一些高阶函数
#
# #### `map`
#
# `map` 是 Julia 中的高阶函数。“高阶”指的是它会*接受一个函数*作为它的传入参数。
# `map` 将传入的函数作用于传入的数据结构的每一个元素。例如执行
#
# ```julia
# map(f, [1, 2, 3])
# ```
# 将得到一个输出数组，数组的每一个元素为 `f` 作用在数组 `[1, 2, 3]` 的对应元素
# ```julia
# [f(1), f(2), f(3)]
# ```
# ------------------------------------------------------------------------------------------

map(f, [1, 2, 3])

# ------------------------------------------------------------------------------------------
# 这里我们把向量 `[1, 2, 3]` 的每一个元素都平方了，而不是对向量 `[1, 2, 3]` 本身进行平方。
#
# `map` 可以接受匿名函数作为参数，像这样
# ------------------------------------------------------------------------------------------

x -> x^3

# ------------------------------------------------------------------------------------------
# 通过
# ------------------------------------------------------------------------------------------

map(x -> x^3, [1, 2, 3])

# ------------------------------------------------------------------------------------------
# 这就完成了对数组 `[1, 2, 3]` 所有元素立方的计算！
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### 广播（broadcast）
#
# 广播（`broadcast`）和 `map` 一样也是高阶函数。
# 广播是一种广义化的 `map`，所以 `map` 能做的它都能做。并不止如此，广播的调用方法和 `map` 一样。
# ------------------------------------------------------------------------------------------

broadcast(f, [1, 2, 3])

# ------------------------------------------------------------------------------------------
# 我们又一次通过函数 `f` 把向量 `[1, 2, 3]` 的每一个元素都平方了
# ——这次是用“广播”函数 `f` 的方式！
#
# 广播（`broadcast`）有个语法糖：调用要广播的函数时，在函数名和输入变量之间加一个 `.`。
#
# 比如：
# ```julia
# broadcast(f, [1, 2, 3])
# ```
# 就等价于
# ```julia
# f.([1, 2, 3])
# ```
# ------------------------------------------------------------------------------------------

f.([1, 2, 3])

# ------------------------------------------------------------------------------------------
# 再一次强调 `map` 和广播方式与如下直接调用函数的区别。
#
# ```julia
# f([1, 2, 3])
# ```
# 我们可以求向量中每个元素的平方，但是不能求向量的平方！
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# 为了再讲清楚这一点，我们再看一下以下调用的区别
#
# ```julia
# f(A)
# ```
# 与
# ```julia
# f.(A)
# ```
# 其中 `A` 是个矩阵：
# ------------------------------------------------------------------------------------------

A = [i + 3*j for j in 0:2, i in 1:3]

f(A)

# ------------------------------------------------------------------------------------------
# 和之前一样，对于矩阵 `A` 来说
# ```
# f(A) = A^2 = A * A
# ```
#
# 另一方面，
# ------------------------------------------------------------------------------------------

B = f.(A)

# ------------------------------------------------------------------------------------------
# 结果包含了矩阵 `A` 中每个元素的平方
#
# 广播的点语法让比较复杂的元素间运算的表示更自然、更贴近数学表示。
#
# 例如：
# ------------------------------------------------------------------------------------------

A .+ 2 .* f.(A) ./ A

# ------------------------------------------------------------------------------------------
# 而不是
# ------------------------------------------------------------------------------------------

broadcast(x -> x + 2 * f(x) / x, A)

# ------------------------------------------------------------------------------------------
# 这两种方法本质上完全一样
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### 练习
#
# #### 6.1
# 声明一个函数 `add_one` 用来给它的输入参数加1。
# ------------------------------------------------------------------------------------------

function add_one(n)
    # 请在此编写你的代码
    n
end

@assert add_one(1) == 2

@assert add_one(11) == 12

# ------------------------------------------------------------------------------------------
# #### 6.2
# 使用函数`map`或者`broadcast`给矩阵`A`的每一个元素加1并把它赋值给`A1`。
# ------------------------------------------------------------------------------------------

# 初始化 A。请勿更改
A = [i + 3*j for j in 0:2, i in 1:3]

# 请在此编写你的代码
A1 = map(x->x, A)

@assert A1 == [2 3 4; 5 6 7; 8 9 10]

# ------------------------------------------------------------------------------------------
# #### 6.3
# 使用广播的点语法给矩阵`A`的每一个元素加1并把它赋值给`A2`。
# ------------------------------------------------------------------------------------------

# 请在此编写你的代码
A2 = A1

@assert A2 == [3 4 5; 6 7 8;9 10 11]

# ------------------------------------------------------------------------------------------
# 请在完成练习运行 `@assert` 块以检测结果是否正确，若无报错即为正确。
# ------------------------------------------------------------------------------------------

max_stars_repo_name: NahsiN/JuliaTutorials

max_issues_repo_name: NahsiN/JuliaTutorials

max_forks_repo_name: NahsiN/JuliaTutorials


text:
function toString(b::SymbolicType)
    #@show b
    b = Basic(b)
    if b.ptr == C_NULL
       error("Trying to print an uninitialized SymEngine Basic variable.")
    end
    a = ccall((:basic_str_julia, libsymengine), Cstring, (Ref{Basic}, ), b)
    string = unsafe_string(a)
    @show string
    ccall((:basic_str_free, libsymengine), Nothing, (Cstring, ), a)
    @show string
    return string
end

Base.show(io::IO, b::SymbolicType) = print(io, toString(b))

" show symengine logo "
mutable struct AsciiArt x end
function ascii_art()
    out = ccall((:ascii_art_str, libsymengine),  Ptr{UInt8},  ())
    AsciiArt(unsafe_string(out))
end

Base.show(io::IO, x::AsciiArt) = print(io, x.x)


text:
"""
Fill volume with spheres of certain radius.

In brief, the model generates a box with spheres in hexagonal packaging, removes the ones outside the volume.
The positions are perturbed by the noise term and finally the system is left to relax by a simple particle model.

# Parameters
 - **f** (Function) Function that returns true if center of sphere is inside the volume
 - **box** (Array{Float,2}) Maximum box where to fill the volumes.
 - **r** (Number) Radius of the spheres

# Optional keyword parameters

 - **N** (Int) Maximum number of particles inside the volume. If NaN (default), there is not upper bound.
 - **noise** (Number) Noise ratio used to perturb the particles from the hexagonal lattice.
 - **platform** (String) Platform in which perform the relaxation step after the noise perturbation.

# Example
```Julia
julia> using AgentBasedModels;
julia> f(x,y,z) = sqrt(x^2+y^2+z^2) < 10;
julia> pos = fillVolumeSpheres(f,[[-10,-10,-10],[10,10,10]],1,noise=0.25);
```
![Figure](assets/FillVolumeSpheres.png)
Figure rendered with [Makie.jl](https://github.com/JuliaPlots/Makie.jl) using meshscatter function.
"""
function fillVolumeSpheres(f,box,rr;N=NaN,noise=0.1,platform="cpu")

    r = rr*0.9
    volumeX, volumeY, volumeZ = latticeCompactHexagonal(box,r,noiseRatio=noise)

    #Extrude
    ext = zeros(Bool,length(volumeX))
    for i in 1:length(volumeX)
        ext[i] = f(volumeX[i],volumeY[i],volumeZ[i])
    end
    volumeX = volumeX[ext]; volumeY = volumeY[ext]; volumeZ = volumeZ[ext];

    #Remove exceding cells
    if N != NaN
        l = length(volumeX)
        while length(volumeX) > N
            rem = rand(1:length(volumeX),l-N)
            nrem = [i for i in 1:l if !(i in rem)]
            volumeX = volumeX[nrem]; volumeY = volumeY[nrem]; volumeZ = volumeZ[nrem];
            l = length(volumeX)
        end
    end

    #Model relax
    relax!(volumeX, volumeY, volumeZ, rr, noise, platform)

    #Remove exceding cells after relaxation
    if N != NaN
        l = length(volumeX)
        while length(volumeX) > N
            rem = rand(1:length(volumeX),l-N)
            nrem = [i for i in 1:l if !(i in rem)]
            volumeX = volumeX[nrem]; volumeY = volumeY[nrem]; volumeZ = volumeZ[nrem];
            l = length(volumeX)
        end
    end

    return volumeX, volumeY, volumeZ
end

function relax!(volumeX, volumeY, volumeZ, r, n, platform="cpu")

    #Community
    if platform == "cpu"
        model = MODEL_BASIC
    else
        model = MODEL_BASIC_GPU
    end

    com = Community(model,N=length(volumeX))
    com[:x] = volumeX
    com[:y] = volumeY
    com[:z] = volumeZ
    com[:D] = 5.
    com[:r] = r
    com[:n] = n
    com[:nn] = 1

    counter = 10

    c = model.evolve(com,dt=0.1,tMax=20.,tSaveStep=19.9)
    com = c[end]
    f = maximum(sqrt.(com[:fx].^2+com[:fy].^2+com[:fz].^2))
    while f > 0.001 && counter > 0
        c = model.evolve(com,dt=0.1,tMax=20.,tSaveStep=19.9)
        com = c[end]
        f = maximum(sqrt.(com[:fx].^2+com[:fy].^2+com[:fz].^2))

        counter -= 1
    end

    volumeX .= c[end][:x]
    volumeY .= c[end][:y]
    volumeZ .= c[end][:z]

    return nothing

end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
export Graph, erdos_renyi, planted_clique, get_dcs_test_instance

typealias TGraph{T} SparseMatrixCSC{T,Int64}
typealias Graph TGraph{Float64}
typealias GraphVec Array{Graph, 1}

function erdos_renyi(n::Integer, p::Real; is_directed::Bool=false)
  A = spones(sprand(n, n, p))
  if is_directed==false
    A = triu(A, 1)
    A = A + A'
  else
    for i=1:n
      A[i, i]=0.0
    end
  end
  return A
end

function erdos_renyi_fast(n::Integer, p::Real; is_directed::Bool=false)
  A = spzeros(Bool, n, n)
  if is_directed
    for i=1:n
      @simd for j=(i+1):n
        @inbounds A[i, j] = rand() < p
        @inbounds A[j, i] = rand() < p
      end
    end
  else
    for i=1:n
      B = sprandbool(1, n - i, p)
      # bi = find(B)
      @simd for j= 1 : (n-i)
        @inbounds A[i, j+i] = B[j]
        @inbounds A[j+i, i] = B[j]
      end

    end
  end
  return A
end

function planted_clique(n::Integer, p::Real, k::Integer; is_directed=false)
  A = erdos_renyi(n, p; is_directed = is_directed)
  S = randperm(n)[1:k;]
  for i=1:k
    for j=1:k
      if i != j
        A[S[i], S[j]] = 1
      end
    end
  end
  return A
end

function get_dcs_test_instance(n::Integer = 15, p::Real = 0.1, k::Integer = 0)
  if k==0 && VERSION.minor == 3
    k = max(5, iceil(sqrt(n)))
  elseif k==0
    k = max(5, ceil(Integer, sqrt(n)))
  end

  G = Graph[]
  push!(G,erdos_renyi(n, p))
  push!(G,planted_clique(n, p, k))

  return G
end


text:
# ---
# title: 789. Escape The Ghosts
# id: problem789
# author: Indigo
# date: 2021-01-29
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/escape-the-ghosts/description/>
# hidden: true
# ---
#
# You are playing a simplified PAC-MAN game on an infinite 2-D grid. You start
# at the point `[0, 0]`, and you are given a destination point `target =
# [xtarget, ytarget]`, which you are trying to get to. There are several ghosts
# on the map with their starting positions given as an array `ghosts`, where
# `ghosts[i] = [xi, yi]` represents the starting position of the `ith` ghost.
# All inputs are **integral coordinates**.
#
# Each turn, you and all the ghosts may independently choose to either **move 1
# unit** in any of the four cardinal directions: north, east, south, or west or
# **stay still**. All actions happen **simultaneously**.
#
# You escape if and only if you can reach the target **before** any ghost
# reaches you. If you reach any square (including the target) at the **same
# time** as a ghost, it **does not** count as an escape.
#
# Return `true` _if it is possible to escape, otherwise return_`false` _._
#
#
#
# **Example 1:**
#
#
#
#     Input: ghosts = [[1,0],[0,3]], target = [0,1]
#     Output: true
#     Explanation: You can reach the destination (0, 1) after 1 turn, while the ghosts located at (1, 0) and (0, 3) cannot catch up with you.
#
#
# **Example 2:**
#
#
#
#     Input: ghosts = [[1,0]], target = [2,0]
#     Output: false
#     Explanation: You need to reach the destination (2, 0), but the ghost at (1, 0) lies between you and the destination.
#
#
# **Example 3:**
#
#
#
#     Input: ghosts = [[2,0]], target = [1,0]
#     Output: false
#     Explanation: The ghost can reach the target at the same time as you.
#
#
# **Example 4:**
#
#
#
#     Input: ghosts = [[5,0],[-10,-2],[0,-5],[-2,-2],[-7,1]], target = [7,7]
#     Output: false
#
#
# **Example 5:**
#
#
#
#     Input: ghosts = [[-1,0],[0,1],[-1,0],[0,1],[-1,0]], target = [0,0]
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `1 <= ghosts.length <= 100`
#   * `ghosts[i].length == 2`
#   * `-104 <= xi, yi <= 104`
#   * There can be **multiple ghosts** in the same location.
#   * `target.length == 2`
#   * `-104 <= xtarget, ytarget <= 104`
#
#
## @lc code=start
using LeetCode

function escape_ghosts(ghosts::Vector{Vector{Int}}, target::Vector{Int})
    dist1 = sum(abs, target)
    dist2 = minimum(sum(abs, g .- target) for g in ghosts)
    return dist1 < dist2
end
## @lc code=end


text:
module VersionVigilante

import Pkg

include("compare_versions.jl")
include("main.jl")
include("parse_project.jl")
include("utils.jl")

end # module


text:
export nagents, AbstractAgent, ABM, AgentBasedModel,
random_activation, by_id, fastest, partial_activation, random_agent,
property_activation, pos_vertex

abstract type AbstractSpace end

"""
All agents must be a mutable subtype of `AbstractAgent`.
Your agent type **must have** at least the `id` field, and if there is a
space structure the `pos` field, (fields are expected in this order)
```julia
mutable struct MyAgent{P} <: AbstractAgent
    id::Int
    pos::P
end
```
Only for grid spaces, `pos` can be an `NTuple`. For arbitrary graph spaces
it must always be an integer (the graph node number).

Your agent type may have other additional fields relevant to your system,
for example variable quantities like "status" or other "counters".
"""
abstract type AbstractAgent end

function correct_pos_type(n, model)
    if typeof(model.space) <: GraphSpace
        return coord2vertex(n, model)
    elseif typeof(model.space) <: GridSpace
        return vertex2coord(n, model)
    end
end

SpaceType=Union{Nothing, AbstractSpace}
struct AgentBasedModel{A<:AbstractAgent, S<:SpaceType, F, P}
    agents::Dict{Int,A}
    space::S
    scheduler::F
    properties::P
end
const ABM = AgentBasedModel
agenttype(::ABM{A}) where {A} = A

"""
    AgentBasedModel(agent_type [, space]; scheduler, properties)
Create an agent based model from the given agent type,
and the `space` (from [`Space`](@ref)).
`ABM` is equivalent with `AgentBasedModel`.
The agents are stored in a dictionary `model.agents`, where the keys are the
agent IDs, while the values are the agents themselves.
It is recommended however to use [`id2agent`](@ref) to get an agent.

`space` can be omitted, in which it will equal to `nothing`.
This means that all agents are virtualy in one node and have no spatial structure.
If space is omitted, some functions that fascilitate agent-space interactions will not work.

Optionally provide a `scheduler` that creates the order with which agents
are activated in the model, and `properties`
for additional model-level properties.
This is accessed as `model.properties` for later use.
"""
function AgentBasedModel(
        ::Type{A}, space::S = nothing;
        scheduler::F = fastest, properties::P = nothing
        ) where {A<:AbstractAgent, S<:SpaceType, F, P}
    agents = Dict{Int, A}()
    return ABM{A, S, F, P}(agents, space, scheduler, properties)
end

function Base.show(io::IO, abm::ABM{A}) where {A}
    s = "AgentBasedModel with $(nagents(abm)) agents of type $(nameof(A))"
    if abm.space == nothing
        s*= "\n no space"
    else
        s*= "\n space: $(nameof(typeof(abm.space))) with $(nv(abm)) nodes and $(ne(abm)) edges"
    end
    s*= "\n scheduler: $(nameof(abm.scheduler))"
    print(io, s)
    if abm.properties ≠ nothing
        print(io, "\n properties: ", abm.properties)
    end
end

"""
    random_agent(model)
Return a random agent from the model.
"""
random_agent(model) = model.agents[rand(keys(model.agents))]

"""
    nagents(model::ABM)
Return the number of agents in the `model`.
"""
nagents(model::ABM) = length(model.agents)

####################################
# Schedulers
####################################
"""
    fastest
Activate all agents once per step in the order dictated by the agent's container,
which is arbitrary (the keys sequence of a dictionary).
This is the fastest way to activate all agents once per step.
"""
fastest(model) = keys(model.agents)

"""
    by_id
Activate agents at each step according to their id.
"""
function by_id(model::ABM)
  agent_ids = sort(collect(keys(model.agents)))
  return agent_ids
end

@deprecate as_added by_id

"""
    random_activation
Activate agents once per step in a random order.
Different random ordering is used at each different step.
"""
function random_activation(model::ABM)
  order = shuffle(collect(keys(model.agents)))
end

"""
    partial_activation(p)
At each step, activate only `p` percentage of randomly chosen agents.
"""
function partial_activation(p::Real)
    function partial(model::ABM{A, S, F, P}) where {A, S, F, P}
        ids = collect(keys(model.agents))
        return randsubseq(ids, p)
    end
    return partial
end

"""
    property_activation(property)
At each step, activate the agents in an order dictated by their `property`,
with agents with greater `property` acting first. `property` is a `Symbol`, which
just dictates which field the agents to compare.
"""
function property_activation(p::Symbol)
    function by_property(model::ABM{A, S, F, P}) where {A, S, F, P}
        ids = collect(keys(model.agents))
        properties = [getproperty(model.agents[id], p) for id in ids]
        s = sortperm(properties)
        return ids[s]
    end
end


text:
const ESC = "\u001B"

gotoANSI(x, y) = print("$ESC[$(y);$(x)H")

gotoANSI(3, 6)
println("Hello")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Documenter, AutomotiveInteraction

makedocs(
    modules = [AutomotiveInteraction],
    format = Documenter.HTML(),
    sitename="AutomotiveInteraction.jl"
    )

deploydocs(
    repo = "github.com/sisl/AutomotiveInteraction.jl.git"
    )

max_stars_repo_name: sisl/AutomotiveInteraction

max_issues_repo_name: sisl/AutomotiveInteraction

max_forks_repo_name: sisl/AutomotiveInteraction


text:
using ModelingToolkit, StaticArrays, LinearAlgebra
using DiffEqBase, SparseArrays
using Test
using NonlinearSolve
using ModelingToolkit: value

canonequal(a, b) = isequal(simplify(a), simplify(b))

# Define some variables
@parameters t σ ρ β
@variables x y z

function test_nlsys_inference(name, sys, vs, ps)
    @testset "NonlinearSystem construction: $name" begin
        @test Set(states(sys))     == Set(value.(vs))
        @test Set(parameters(sys)) == Set(value.(ps))
    end
end

# Define a nonlinear system
eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
@named ns = NonlinearSystem(eqs, [x,y,z], [σ,ρ,β], defaults = Dict(x => 2))
@test eval(toexpr(ns)) == ns
test_nlsys_inference("standard", ns, (x, y, z), (σ, ρ, β))
@test begin
    f = eval(generate_function(ns, [x,y,z], [σ,ρ,β])[2])
    du = [0.0, 0.0, 0.0]
    f(du, [1,2,3], [1,2,3])
    du ≈ [1, -3, -7]
end

# Now nonlinear system with only variables
@variables x y z
@parameters σ ρ β

# Define a nonlinear system
eqs = [0 ~ σ*(y-x),
       y ~ x*(ρ-z),
       β*z ~ x*y]
@named ns = NonlinearSystem(eqs, [x,y,z], [σ,ρ,β])
jac = calculate_jacobian(ns)
@testset "nlsys jacobian" begin
    @test canonequal(jac[1,1], σ * -1)
    @test canonequal(jac[1,2], σ)
    @test canonequal(jac[1,3], 0)
    @test canonequal(jac[2,1], ρ - z)
    @test canonequal(jac[2,2], -1)
    @test canonequal(jac[2,3], x * -1)
    @test canonequal(jac[3,1], y)
    @test canonequal(jac[3,2], x)
    @test canonequal(jac[3,3], -1 * β)
end
nlsys_func = generate_function(ns, [x,y,z], [σ,ρ,β])
jac_func = generate_jacobian(ns)
f = @eval eval(nlsys_func)

# Intermediate calculations
a = y - x
# Define a nonlinear system
eqs = [0 ~ σ*a,
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
@named ns = NonlinearSystem(eqs, [x,y,z], [σ,ρ,β])
nlsys_func = generate_function(ns, [x,y,z], [σ,ρ,β])
nf = NonlinearFunction(ns)
jac = calculate_jacobian(ns)

@test ModelingToolkit.jacobian_sparsity(ns).colptr == sparse(jac).colptr
@test ModelingToolkit.jacobian_sparsity(ns).rowval == sparse(jac).rowval

jac = generate_jacobian(ns)

prob = NonlinearProblem(ns,ones(3),ones(3))
sol = solve(prob,NewtonRaphson())
@test sol.u[1] ≈ sol.u[2]

@test_throws ArgumentError NonlinearProblem(ns,ones(4),ones(3))

@variables u F s a
eqs1 = [
        0 ~ σ*(y-x) + F,
        0 ~ x*(ρ-z)-u,
        0 ~ x*y - β*z,
        0 ~ x + y - z - u,
       ]

lorenz = name -> NonlinearSystem(eqs1, [x,y,z,u,F], [σ,ρ,β], name=name)
lorenz1 = lorenz(:lorenz1)
@test_throws ArgumentError NonlinearProblem(lorenz1, zeros(5))
lorenz2 = lorenz(:lorenz2)
@named connected = NonlinearSystem([s ~ a + lorenz1.x
                             lorenz2.y ~ s
                             lorenz1.F ~ lorenz2.u
                             lorenz2.F ~ lorenz1.u], [s, a], [], systems=[lorenz1,lorenz2])
@test_nowarn alias_elimination(connected)

# system promotion
using OrdinaryDiffEq
@variables t
D = Differential(t)
@named subsys = convert_system(ODESystem, lorenz1, t)
@named sys = ODESystem([D(subsys.x) ~ subsys.x + subsys.x], t, systems=[subsys])
sys = structural_simplify(sys)
u0 = [subsys.x => 1, subsys.z => 2.0]
prob = ODEProblem(sys, u0, (0, 1.0), [subsys.σ=>1,subsys.ρ=>2,subsys.β=>3])
sol = solve(prob, Rodas5())
@test sol[subsys.x] + sol[subsys.y] - sol[subsys.z] ≈ sol[subsys.u]
@test_throws ArgumentError convert_system(ODESystem, sys, t)

@parameters t σ ρ β
@variables x y z

# Define a nonlinear system
eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
@named ns = NonlinearSystem(eqs, [x,y,z], [σ,ρ,β])
np = NonlinearProblem(ns, [0,0,0], [1,2,3], jac=true, sparse=true)
@test calculate_jacobian(ns, sparse=true) isa SparseMatrixCSC

# issue #819
@testset "Combined system name collisions" begin
       function makesys(name)
           @parameters a
           @variables x f

           NonlinearSystem([0 ~ -a * x + f], [x,f], [a]; name)
       end

       function issue819()
           sys1 = makesys(:sys1)
           sys2 = makesys(:sys1)
           @test_throws ArgumentError NonlinearSystem([sys2.f ~ sys1.x, sys1.f ~ 0], [], [], systems = [sys1, sys2], name=:foo)
       end
       issue819()
end

# issue #1115
@testset "Extending a NonlinearSystem with no iv" begin
    @parameters a b
    @variables x y
    eqs1 = [
        0 ~ a * x
    ]
    eqs2 = [
        0 ~ b * y
    ]

    @named sys1 = NonlinearSystem(eqs1, [x], [a])
    @named sys2 = NonlinearSystem(eqs2, [y], [b])
    @named sys3 = extend(sys1, sys2)

    @test isequal(union(Set(parameters(sys1)), Set(parameters(sys2))), Set(parameters(sys3)))
    @test isequal(union(Set(states(sys1)), Set(states(sys2))), Set(states(sys3)))
    @test isequal(union(Set(equations(sys1)), Set(equations(sys2))), Set(equations(sys3)))
end

# observed variable handling
@variables t x(t) RHS(t)
@parameters τ
@named fol = NonlinearSystem([0  ~ (1 - x)/τ], [x], [τ]; observed=[RHS ~ (1 - x)/τ])
@test isequal(RHS, @nonamespace fol.RHS)
RHS2 = RHS
@unpack RHS = fol
@test isequal(RHS, RHS2)

# issue #1358
@variables t
@variables v1(t) v2(t) i1(t) i2(t)
eq = [
    v1 ~ sin(2pi*t)
    v1 - v2 ~ i1
    v2 ~ i2
    i1 ~ i2
]
@named sys = ODESystem(eq)
@test length(equations(structural_simplify(sys))) == 0

max_stars_repo_name: paulrschrater/ModelingToolsJulia

max_issues_repo_name: paulrschrater/ModelingToolsJulia

max_forks_repo_name: paulrschrater/ModelingToolsJulia


text:
function report(
    reportdir::String, invprob::InvestmentProblem{R,G1,G2,G3,I,T,P}
) where {R,G1,G2,G3,I,T,P}

    status = termination_status(invprob.model)

    if status != MOI.OPTIMAL
        status == MOI.OPTIMIZE_NOT_CALLED && error("Model has not yet been solved")
        error("Model was not solved to optimality: status is $(status)")
    end

    isdir(reportdir) || mkpath(reportdir)

    annualmarkets = DataFrame(
         scenario=String[], market=String[], demand=Float64[], price=Float64[])

    hourlymarkets = DataFrame(
        scenario=String[], market=String[], region=String[],
        period=String[], timestep=Int[], demand=Float64[], price=Float64[])

    dispatches = DataFrame(
        scenario=String[], market=String[], region=String[], gen=String[],
        period=String[], timestep=Int[], value=Float64[])

    for (scenario, isleaf) in scenarios(invprob)

        println("\n", scenario.name)
        oneoffdiscount, recurringdiscount = discount(scenario)
        periodweights = scenario.periodweights

        # Dispatch
        for market in [:energydispatch, :raisereserve, :lowerreserve],
            (r, regionname) in enumerate(invprob.regionnames),
            gentype in [:thermal, :variable, :storage]

            gentech = getfield(invprob.technologies, gentype)
            gennames = gentech.name

            if gentype == :storage
                gencapacities = gentech.maxdischarge
                dispatch = if market == :energydispatch
                    scenario.operations.storage.energydischarge .-
                    scenario.operations.storage.energycharge
                else
                    getfield(scenario.operations.storage, market)
                end
            else
                gencapacities = gentech.maxgen
                dispatch = getfield(getfield(scenario.operations, gentype), market)
            end

            marketname = market == :energydispatch ? "energy" : string(market)

            for (g, genname) in enumerate(gennames),
                (p, periodname) in enumerate(invprob.periodnames),
                t in 1:T

                    push!(dispatches,
                          (scenario=scenario.name, market=marketname,
                           region=regionname, gen=genname,
                           period=periodname, timestep=t,
                           value=value(dispatch[r,g,t,p])))

            end

        end

        # Capacity Market
        capmkt = scenario.markets.capacity
        capmkt_demand = value(capmkt.totalcontribution)
        capmkt_price =
            if capmkt_demand < capmkt.maxprice_capacity
                capmkt.maxprice
            elseif capmkt_demand < capmkt.midprice_capacity
                value(capmkt.seg1price)
            elseif capmkt_demand < capmkt.zeroprice_capacity
                value(capmkt.seg2price)
            else
                0
            end
        push!(annualmarkets, (
                scenario=scenario.name, market="capacity",
                demand=capmkt_demand, price=capmkt_price))

        # REC market
        recmkt = scenario.markets.rec
        recmkt_demand = recmkt.demand
        recmkt_price = dual(recmkt.marketclearing)
        push!(annualmarkets, (
            scenario=scenario.name, market="rec",
            demand=recmkt_demand, price=recmkt_price/recurringdiscount))

        # Hourly, regional market demand + prices
        for market in [:energy, :raisereserve, :lowerreserve],
            (r, regionname) in enumerate(invprob.regionnames)

            demand = getfield(scenario.markets, market).demand
            price = dual.(getfield(scenario.markets, market).marketclearing)

            for (p, periodname) in enumerate(invprob.periodnames), t in 1:T

                push!(hourlymarkets, (
                    scenario=scenario.name, market=string(market),
                    region=regionname, period=periodname, timestep=t,
                    demand=demand[r,t,p],
                    price=price[r,t,p]/periodweights[p]/recurringdiscount))

            end

        end

        if isleaf

            # TODO: Quantile plot of energy prices for each leaf scenario?

            # Save series of buildout progressions, by leaf scenario

            scen = scenario
            buildouts = DataFrame[]

            while !isnothing(scen)
                println(scen.name)
                push!(buildouts, investmentstates(scen))
                scen = scen.parent
            end

            CSV.write(joinpath(reportdir, "buildouts_$(scenario.name).csv"),
                      vcat(reverse(buildouts)...))

        end

    end

    dispatches.value[-1e-4 .< dispatches.value .< 0] .= 0.
    CSV.write(joinpath(reportdir, "dispatch.csv"), dispatches)
    CSV.write(joinpath(reportdir, "markets_hourly.csv"), hourlymarkets)
    CSV.write(joinpath(reportdir, "markets_annual.csv"), annualmarkets)

end

function investmentstates(scen::Scenario)

    result = DataFrame(
        scenario=String[], tech=String[], state=String[],
        count=Int[], capacity=Float64[])

    for techtype in [:thermal, :variable, :storage]

        techs = getfield(scen.investmentproblem.technologies, techtype)
        investments = getfield(scen.investments, techtype)

        for state in [:vesting, :holding, :building, :dispatching, :retired]

            counts = round.(Int, value.(getfield(investments, state)))

            for i in 1:length(techs.name)

                techname = techs.name[i]
                techcapacity = techtype == :storage ?
                    techs.maxdischarge[i] : techs.maxgen[i]

                count = sum(counts[:, i])
                push!(result, (scenario=scen.name, tech=techname,
                               state=string(state),
                               count=count,
                               capacity=count*techcapacity))

            end

        end

    end

    return result

end


text:
using EquationEditor
const EE = EquationEditor
const logger =  EE.WebSocketLogger(stderr, Base.CoreLogging.Debug, right_justify = 10)
const server, task = serveEE(logger)
EE.open_browser("chrome", url = "localhost:8000/interact.html")

max_stars_repo_name: hustf/EquationEditor

max_issues_repo_name: hustf/EquationEditor

max_forks_repo_name: hustf/EquationEditor


text:
function test_ilmm(rng, kernels, H, x_train, x_test, y_train, y_test)

    # Construct ILMM and equivalent GP against which to test.
    ilmm = ILMM(independent_mogp(map(GP, kernels)), H)
    n_ilmm = GP(LinearMixingModelKernel(kernels, H'))

    ilmmx = ilmm(x_train, 1e-6)
    n_ilmmx = n_ilmm(x_train, 1e-6)

    @test isapprox(mean(ilmmx), mean(n_ilmmx))
    @test isapprox(var(ilmmx), var(n_ilmmx))
    @test isapprox(cov(ilmmx), cov(n_ilmmx))
    @test isapprox(logpdf(ilmmx, y_train), logpdf(n_ilmmx, y_train))
    @test _is_approx(marginals(ilmmx), marginals(n_ilmmx))
    @test length(rand(rng, ilmmx)) == size(H, 1) * length(x_train.x)

    p_ilmmx = posterior(ilmmx, y_train)
    p_n_ilmmx = posterior(n_ilmmx, y_train)

    pi = p_ilmmx(x_test, 1e-6)
    pni = p_n_ilmmx(x_test, 1e-6)

    @test isapprox(mean(pi), mean(pni))
    @test isapprox(var(pi), var(pni))
    @test isapprox(logpdf(pi, y_test), logpdf(pni, y_test))
    @test _is_approx(marginals(pi), marginals(pni))
    @test length(rand(rng, pi)) == size(H, 1) * length(x_test.x)

    test_sampling_consistency(rng, ilmm, x_train)

    @test gradient(logpdf, ilmmx, y_train) isa Tuple
    @test gradient(logpdf, pi, y_test) isa Tuple

    @testset "primary_public_interface" begin
        test_finitegp_primary_and_secondary_public_interface(rng, ilmmx)
        test_finitegp_primary_and_secondary_public_interface(rng, pi)
    end
end

@testset "ilmm" begin
    rng = Random.seed!(04161999)
    x_train, x_test, y_train, y_test = generate_toy_data(rng)

    @testset "Full Rank, Dense H" begin
        H = rand(3, 3)
        kernels = [SEKernel(), Matern32Kernel(), Matern32Kernel()]
        test_ilmm(rng, kernels, H, x_train, x_test, y_train, y_test)
    end

    @testset "M Latent Processes" begin
        H = rand(3, 2)
        k1, k2 = SEKernel(), Matern32Kernel()
        test_ilmm(rng, [k1, k2], H, x_train, x_test, y_train, y_test)

        @testset "util" begin
            Σ = Diagonal(Fill(2, 3))
            @test noise_var(Σ) == 2

            y = rand(16)
            @test size(reshape_y(y, 8)) == (2, 8)
            @test size(reshape_y(y, 2)) == (8, 2)

            fs = independent_mogp([GP(Matern32Kernel())])
            H = rand(2, 1)
            x = MOInputIsotopicByOutputs(ColVecs(rand(2, 2)), 2)
            ilmm = ILMM(fs, H)
            ilmmx = ilmm(x, 0.1)
            @test (fs, H, 0.1, x.x) == unpack(ilmmx)

            @test get_latent_gp(ilmm) == fs
        end
    end

    @testset "1 Latent Processes" begin
        H = rand(3, 1)
        k = SEKernel()
        test_ilmm(rng, [k], H, x_train, x_test, y_train, y_test)
    end
end
@info "Ran ilmm tests."


text:
module BenchmarkProfiles

import Plots, NaNMath

export performance_ratios, performance_profile
export data_ratios, data_profile

include("performance_profiles.jl")
include("data_profiles.jl")

end


text:
function get_satellites(obj::EphemerisFile)
    return jcall(obj, "getSatellites", Map, ())
end


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 = "Xorg_libXpm"
version = v"3.5.13"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://www.x.org/archive/individual/lib/libXpm-$(version).tar.bz2",
                  "9cd1da57588b6cb71450eff2273ef6b657537a9ac4d02d0014228845b935ac25"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libXpm-*/
# We need a native xgettext
apk add gettext
./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 = [p for p in supported_platforms() if Sys.islinux(p) || Sys.isfreebsd(p)]

# The products that we will ensure are always built
products = [
    LibraryProduct("libXpm", :libXpm)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(name="Xorg_xorgproto_jll", uuid="c4d99508-4286-5418-9131-c86396af500b")),
    Dependency(PackageSpec(name="Xorg_libX11_jll", uuid="4f6342f7-b3d2-589e-9d20-edeb45f2b2bc")),
    # TOOD: verify Gettext is actually needed at runtime
    Dependency(PackageSpec(name="Gettext_jll", uuid="78b55507-aeef-58d4-861c-77aaff3498b1"), v"0.20.1"; compat="=0.20.1"),
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
######################################################################
# solver.jl
#
# defines solver and planner types and provides convenience functions
#
# action(planner, b) function provided in action.jl
######################################################################

struct DefaultPolicy <: Policy end
struct DefaultUpdater <: Updater end

"""
Implements anytime error minimization search (AEMS).
Specifically, this is AEMS2, which generally outperforms AEMS1.

Example with some optional arguments:

`AEMSSolver(max_iterations = 100, action_selector = :L)`

Fields:

- `max_iterations`
    Maximum node expansions per action.

- `max_time`
    Maximum time (in seconds) to spend per action.

- `updater`
    Updater used to propagate belief nodes. Defaults to DiscreteUpdater.

- `lower_bound`
    Subtype of `Policy`. Defaults to fixed action policy.

- `upper_bound`
    Subtype of `Policy`. Defaults to FIB policy.

- `root_manager`
    Allowable values are `:clear`, `:belief`, or `:user`.
    Defaults to `:clear`.

- `action_selector`
    Allowable values are `:U` or `:L`. Defaults to `:U`
"""
mutable struct AEMSSolver{U<:Updater, PL<:Policy, PU<:Policy} <: Solver
    max_iterations::Int
    max_time::Float64   # max time per action, in seconds

    updater::U

    lower_bound::PL
    upper_bound::PU

    root_manager::Symbol
    action_selector::Symbol
end
function AEMSSolver(;
                    max_iterations::Int = 1000,
                    max_time::Real = 1.0,
                    updater = DefaultUpdater(),
                    lower_bound = DefaultPolicy(),
                    upper_bound = DefaultPolicy(),
                    root_manager::Symbol = :clear,
                    action_selector::Symbol = :U
                   )

    @assert in(action_selector, (:U, :L))
    @assert in(root_manager, (:clear, :belief, :user))

    return AEMSSolver(max_iterations, float(max_time), updater, lower_bound, upper_bound, root_manager, action_selector)
end

# planner
struct AEMSPlanner{S <: AEMSSolver,
                   P <: POMDP,
                   U <: Updater,
                   PL <: Policy,
                   PU <: Policy,
                   GT <: Graph,
                   A,
                   O
                  } <: Policy

    solver::S           # contains solver parameters
    pomdp::P            # model
    updater::U
    G::GT
    lower_bound::PL      # lower bound
    upper_bound::PU      # upper bound
    action_list::Vector{A}
    obs_list::Vector{O}
    root_manager::Symbol
end
function AEMSPlanner(s, pomdp::POMDP, up, lb, ub)
    b0 = initialize_belief(up, initialstate(pomdp))
    bn_type = BeliefNode{typeof(b0)}
    G = Graph{bn_type}(discount(pomdp))
    a_list = ordered_actions(pomdp)
    o_list = ordered_observations(pomdp)
    rm = s.root_manager
    return AEMSPlanner(s, pomdp, up, G, lb, ub, a_list, o_list, rm)
end

"""
Clears existing action/observation graph.
"""
clear_graph!(planner::AEMSPlanner) = clear_graph!(planner.G)

"""
Returns tuple `(nb, na)`, numbers of belief and action nodes in graph
"""
count_nodes(planner::AEMSPlanner) = planner.G.nb, planner.G.na

function update_root(planner::AEMSPlanner, a, o)
    if planner.root_manager != :user
        error("User is trying to update root but (root_manager != :user).")
    end
    ai = actionindex(planner.pomdp, a)
    oi = obsindex(planner.pomdp, o)

    original_root = get_root(planner.G)

    an_ind = original_root.children.start + ai - 1
    an = get_an(planner, an_ind)

    new_root_ind = an.children.start + oi - 1
    planner.G.root_ind = new_root_ind
    return planner      # to prevent it from returning planner.G.root_ind
end

# convenience functions
get_bn(planner::AEMSPlanner, bn_idx::Int) = planner.G.belief_nodes[bn_idx]
get_an(planner::AEMSPlanner, an_idx::Int) = planner.G.action_nodes[an_idx]

# SOLVE
function solve(solver::AEMSSolver, pomdp::POMDP)
    # if no updater was given, default to discrete updater
    up = solver.updater
    if typeof(up) == DefaultUpdater
        up = DiscreteUpdater(pomdp)
    end

    # if no lower bound was given to solver, default to blind
    lb = solver.lower_bound
    if typeof(lb) == DefaultPolicy
        lb = solve(FixedActionSolver(), pomdp)
    end

    # if no upper bound was passed to solver, default to FIB
    ub = solver.upper_bound
    if typeof(ub) == DefaultPolicy
        ub = solve(FIBSolver(), pomdp)
    end

    AEMSPlanner(solver, pomdp, up, lb, ub)
end

updater(planner::AEMSPlanner) = planner.updater

function determine_root_node(planner::AEMSPlanner, b)
    if planner.root_manager == :clear
        clear_graph!(planner.G)
    end
    if planner.G.nb == 0
        L = value(planner.lower_bound, b)
        U = value(planner.upper_bound, b)
        bn_root = BeliefNode(b, L, U)
        add_node(planner.G, bn_root)
        return bn_root
    end

    if planner.root_manager == :belief
        bn_root = get_root(planner.G)

        # if this root has same belief, make it the root
        if bn_root.b == b
            return bn_root
        end

        # search through children to see if they might be root
        for an_idx in bn_root.children
            an = get_an(planner, an_idx)
            for bn_idx in an.children
                bn = get_bn(planner, bn_idx)
                if bn.b == b
                    planner.G.root_ind = bn.ind
                    return get_root(planner.G)
                end
            end
        end
        error("None of the child beliefs match input")
        # another option is to clear graph and start again
    end

    # if we get to here, then we must have :user
    return get_root(planner.G)
end


text:
# TODO: add more target distributions and make them iteratable
# TODO: Integrate with https://github.com/xukai92/VecTargets.jl to achieve goal noted
#       above.

# Dimension of testing distribution
const D = 5
# Tolerance ratio
const TRATIO = Int == Int64 ? 1 : 2
# Deterministic tolerance
const DETATOL = 1e-3 * D * TRATIO
# Random tolerance
const RNDATOL = 5e-2 * D * TRATIO * 2

# Hand-coded multivariate Gaussian

const gaussian_m = zeros(D)
const gaussian_s = ones(D)

struct Gaussian{Tm, Ts}
    m::Tm
    s::Ts
end

function ℓπ_gaussian(g::AbstractVecOrMat{T}, s) where {T}
    return .-(log(2 * T(pi)) .+ 2 .* log.(s) .+ abs2.(g) ./ s.^2) ./ 2
end

ℓπ_gaussian(m, s, x) = ℓπ_gaussian(m .- x, s)

function ∇ℓπ_gaussianl(m, s, x)
    g = m .- x
    v = ℓπ_gaussian(g, s)
    return v, g
end

function get_ℓπ(g::Gaussian)
    ℓπ(x::AbstractVector) = sum(ℓπ_gaussian(g.m, g.s, x))
    ℓπ(x::AbstractMatrix) = dropdims(sum(ℓπ_gaussian(g.m, g.s, x); dims=1); dims=1)
    return ℓπ
end

function get_∇ℓπ(g::Gaussian)
    function ∇ℓπ(x::AbstractVector)
        val, grad = ∇ℓπ_gaussianl(g.m, g.s, x)
        return sum(val), grad
    end
    function ∇ℓπ(x::AbstractMatrix)
        val, grad = ∇ℓπ_gaussianl(g.m, g.s, x)
        return dropdims(sum(val; dims=1); dims=1), grad
    end
    return ∇ℓπ
end

ℓπ = get_ℓπ(Gaussian(gaussian_m, gaussian_s))
∂ℓπ∂θ = get_∇ℓπ(Gaussian(gaussian_m, gaussian_s))

# For the Turing model
# @model gdemo() = begin
#     s ~ InverseGamma(2, 3)
#     m ~ Normal(0, sqrt(s))
#     1.5 ~ Normal(m, sqrt(s))
#     2.0 ~ Normal(m, sqrt(s))
#     return s, m
# end

using Distributions: logpdf, InverseGamma, Normal
using Bijectors: invlink, logpdf_with_trans

function invlink_gdemo(θ)
    s = invlink(InverseGamma(2, 3), θ[1])
    m = θ[2]
    return [s, m]
end

function ℓπ_gdemo(θ)
    s, m = invlink_gdemo(θ)
    logprior = logpdf_with_trans(InverseGamma(2, 3), s, true) + logpdf(Normal(0, sqrt(s)), m)
    loglikelihood = logpdf(Normal(m, sqrt(s)), 1.5) + logpdf(Normal(m, sqrt(s)), 2.0)
    return logprior + loglikelihood
end

using Distributions: MvNormal
import Turing

Turing.@model function mvntest(θ=missing, x=missing)
    θ ~ MvNormal(zeros(D), 2)
    x ~ Normal(sum(θ), 1)
    return θ, x
end

function get_primitives(x, modelgen)
    spl_prior = Turing.SampleFromPrior()
    function ℓπ(θ)
        vi = Turing.VarInfo(model)
        vi[spl_prior] = θ
        model(vi, spl_prior)
        Turing.getlogp(vi)
    end
    adbackend = Turing.Core.ForwardDiffAD{40}
    alg_ad = Turing.HMC{adbackend}(0.1, 1)
    model = modelgen(missing, x)
    vi = Turing.VarInfo(model)
    spl = Turing.Sampler(alg_ad, model)
    Turing.Core.link!(vi, spl)
    ∂ℓπ∂θ = θ -> Turing.Core.gradient_logp(adbackend(), θ, vi, model, spl)
    θ₀ = Turing.VarInfo(model)[Turing.SampleFromPrior()]
    return ℓπ, ∂ℓπ∂θ, θ₀
end

function rand_θ_given(x, modelgen, metric, κ; n_samples=20)
    ℓπ, ∂ℓπ∂θ, θ₀ = get_primitives(x, modelgen)
    h = Hamiltonian(metric, ℓπ, ∂ℓπ∂θ)
    samples, stats = sample(h, κ, θ₀, n_samples; verbose=false, progress=false)
    s = samples[end]
    return length(s) == 1 ? s[1] : s
end

# Test function
geweke_g(θ, x) = cat(θ, x; dims=1)

test_show(x) = test_show(s -> length(s) > 0, x)
function test_show(pred, x)
    io = IOBuffer(; append = true)
    show(io, x)
    s = read(io, String)
    @test pred(s)
end


text:
using BenchmarkTools
using Printf
using Transducers

#### Parsing

function parseEntry(s::AbstractString)::Array{Bool,1}
    out = Array{Bool, 1}(undef, length(s))
    for i = 1:length(s)
        out[i] = (s[i] == '#')
    end
    out
end

#### Reading

function readMap(path)
    readlines(path) |> Map(line -> parseEntry(line)) |> collect
end

#### Logic

function countTrees(map, slope)
    slopeX, slopeY = slope
    count = 0; x = 1; y = 1
    while y <= length(map)
        row = map[y]
        if row[((x - 1) % length(row)) + 1]
            count += 1
        end
        x += slopeX; y += slopeY
    end
    count
end

function multiCountTrees(map, slopes)
    foldxt(*, slopes |> Map(slope -> countTrees(map, slope)))
end

#### Test

exampleMap = readMap("src/advent/day3/example.txt")
testMap = readMap("src/advent/day3/input.txt")

######### Part 1

testSlope = (3, 1)
@assert countTrees(exampleMap, testSlope) == 7

@printf("part1 = %d\n", countTrees(testMap, testSlope))
@btime countTrees(testMap, testSlope);

######### Part 2

testSlopes = ((1, 1), (3, 1), (5, 1), (7, 1), (1, 2))

@assert multiCountTrees(exampleMap, testSlopes) == 336

@printf("part2 = %d\n", multiCountTrees(testMap, testSlopes))
@btime multiCountTrees(testMap, testSlopes);


text:
@doc raw"""
```
DiscreteShape(iterations::Int = 2000, smoothness::Int = 1000, albedo::Real = Inf, illumination_direction::Vector{T} where T <: Real = [Inf, Inf, Inf])
```
Contins `DiscreteShape` algorithm which ttempts to produce a heightmap from a
grayscale image by minimization of a set of Euler-Lagrange equations. This is done
discretely at each point in the image utilizing the second derivates of the
surface normals and their Fourier Transforms.
# Output
Returns a `DiscreteShape` functor which can be run on an image to attempt reconstruction
of the surface contianed in the image.
# Details
The algorithm attempts to minimise the brightness deviation ``ϵ`` through the
brightness constraint ``\epsilon_1`` and smoothness constraint ``\epsilon_2``
as defined bellow:
```math
\epsilon=\iint((\epsilon_1+\lambda\epsilon_2))dxdy=\iint((E(x,y)-R(p,q))^2+
\lambda(p_x^2+q_x^2+p_y^2+q_y^2))
```
This is minimised using the Euler-Lagrange equations defined as;
```math
\dfrac{\delta\epsilon}{\delta p}-\dfrac{\delta}{\delta x}\dfrac{\delta\epsilon}{
\delta p_x}-\dfrac{\delta}{\delta y}\dfrac{\delta\epsilon}{\delta p_y}=0
```
and
```math
\dfrac{\delta
\epsilon}{\delta q}-\dfrac{\delta}{\delta x}\dfrac{\delta\epsilon}{\delta q_x}-
\dfrac{\delta}{\delta y}\dfrac{\delta\epsilon}{\delta q_y}=0
```
which become:
```math
\begin{gathered}
-2(E-R)\dfrac{\delta R}{\delta p}-2\lambda p_{xx}-2\lambda p_{yy}=0\\
-2(E-R)\dfrac{\delta R}{\delta p}-2\lambda q_{xx}-2\lambda q_{yy}=0
\end{gathered}
```
which can be further simplified to give:
```math
\begin{gathered}
\nabla^2p=\dfrac{1}{\lambda}(R-E)\dfrac{\delta R}{\delta p}\\
\nabla^2q=\dfrac{1}{\lambda}(R-E)\dfrac{\delta R}{\delta q}
\end{gathered}
```
where ``\nabla^2p=p_{xx}+p_{yy}`` and ``\nabla^2q=q_{xx}+q_{yy}`` are Laplacians
of p and q.

However, for computation we are dealing with a discreate case of these equations
which can be defined as below:
```math
\begin{gathered}
p_{i,j}=\bar{p}_{i,j}+\dfrac{1}{4\lambda}(E-R)\dfrac{\delta R}{\delta p}\\q_{i,j}=
\bar{q}_{i,j}+\dfrac{1}{4\lambda}(E-R)\dfrac{\delta R}{\delta q}\\
\end{gathered}
```
where
```math
\bar{p}_{i,j}=\dfrac{p_{i+1,j}+p_{i-1,j}+p_{i,j+1}+p_{i,j-1}}{4}
```
and
```math
\bar{q}_{i,j}=\dfrac{q_{i+1,j}+q_{i-1,j}+q_{i,j+1}+q_{i,j-1}}{4}
```
Finally, the algorithm needs to enforce integrability on p and q and retrieve
the surface Z. This can be done by taking the Fast Fourier Transform of p and q
()``c_p(\omega_x,\omega_y)`` and ``c_q(\omega_x,\omega_y)``) and then using the
Inverse Fast Fourier Transform to recover Z and update p and q as per bellow:
```math
\begin{gathered}
p=\sum c_p(\omega_x,\omega_y)e^{j(\omega_xx+\omega_yy)}\\
q=\sum c_q(\omega_x,\omega_y)e^{j(\omega_xx+\omega_yy)}\\Z=\sum c(\omega_x,\omega_y)
e^{j(\omega_xx+\omega_yy)}\\
\end{gathered}
```
where
```math
c(\omega_x,\omega_y)=\dfrac{-j(\omega_xc_p(\omega_x,\omega_y)+\omega_yc_q(\omega_x
,\omega_y))}{\omega_x^2+\omega_y^2}
```
# Arguments
The arguments are described in more detail below.
## `albedo`
A `Real` that specifies the albedo (amount of light reflected) of the image. If
`albedo` is specified to must the `illumination_direction`. Defults to `Inf` which
will trigger the algorithm to run `estimate_img_properties`.
## `illumination_direction`
A `Vector{T} where T <: Real` that specifies the tilt value to be used by the
algorithm. The `illumination_direction` should be a vector of the form [x,y,z]
where x,y,z are int he range [0,1]. If `illumination_direction` is specified
to must the `albedo`. Defults to `[Inf, Inf, Inf]` which
will trigger the algorithm to run `estimate_img_properties`.
## `iterations`
An `Int` that specifies the number of iterations the algorithm is to perform. If
left unspecified a default value of 2000 is used.
## `smoothness`
An `Int` that specifies the strength of the smoothness constraint in the minimised
function. Defults to `1000`.
!!! note
    If `albedo` and `illumination_direction` are not defined (i.e. have defulted
    to `Inf`) they will be calculated at runtime using `estimate_img_properties`.
# Example
Compute the heightmap for a synthetic image generated by `generate_surface`.
```julia
using Images, Makie, ShapeFromShading

#generate synthetic image
img = generate_surface(SynthSphere(50), 1.0, [0.2,0,0.9])

#calculate the heightmap (using 500 iterations)
discreteShape = DiscreteShape(iterations = 500, albedo = 1.0, illumination_direction = [0.2,0,0.9])
Z,p,q = discreteShape(img)

#normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)

#display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:2
surface(r, r, Z)
```
# Reference
1. S. Elhabian, "Hands on Shape from Shading", Computer Vision and Image Processing, 2008.
"""
function (algorithm::DiscreteShape)(img::AbstractArray)
    ρ = algorithm.albedo
    I = algorithm.illumination_direction
    smoothness = algorithm.smoothness
    iterations = algorithm.iterations

    if ρ == Inf || I == [Inf, Inf, Inf]
        ρ,I,σ,τ = estimate_img_properties(img)
    end

    E = Array{Float64}(img)

    #initialize variables
    p = zeros(Complex{Float64},axes(E))
    q = zeros(Complex{Float64},axes(E))
    R = zeros(Complex{Float64},axes(E))
    Z = zeros(axes(E))
    return solve_EulerLagrange(ρ, I, iterations, p, q, R, smoothness, E, Z)
end


text:
using ReadWriteLocks
using Test

@testset "ReadWriteLock" begin

@testset "Constructor" begin
    @test ReadWriteLock() isa ReadWriteLock{ReentrantLock}
    @test ReadWriteLock(0, false, Threads.Mutex()) isa ReadWriteLock{Threads.Mutex}
end

@testset "Single-threaded tests" begin
    @testset "Initialization" begin
        rwlock = ReadWriteLock()

        @test rwlock.read_lock.rwlock == rwlock
        @test rwlock.write_lock.rwlock == rwlock
        @test rwlock.readers == 0
        @test rwlock.writer == false
        @test read_lock(rwlock) == rwlock.read_lock
        @test write_lock(rwlock) == rwlock.write_lock
    end
end  # Single-threaded

@testset "Two-threaded tests" begin
    @testset "Read locks" begin
        NUM_LOCKS = 10

        @testset "$NUM_LOCKS locks" begin
            rwlock = ReadWriteLock()
            rlock = read_lock(rwlock)

            @async begin
                @test rwlock.writer == false
                @test rwlock.readers == 0
                for i = 1:NUM_LOCKS
                    @test lock!(rlock) == nothing
                    # lock(rlock)
                    @test rwlock.readers == i
                end
            end

            sleep(1)

            @test rwlock.readers == NUM_LOCKS

            @async begin
                @test rwlock.writer == false
                @test rwlock.readers == NUM_LOCKS
                for i in NUM_LOCKS:-1:1
                    @test unlock!(rlock) == nothing
                    @test rwlock.readers == i - 1
                end
                @test rwlock.readers == 0
            end

            sleep(1)

            @test rwlock.readers == 0
        end
    end

    @testset "Write locks" begin
        @testset "two locks" begin
            rwlock = ReadWriteLock()
            wlock = write_lock(rwlock)

            c = Channel{Symbol}(1)
            put!(c, :pretest)

            @async begin
                @test rwlock.writer == false
                @test rwlock.readers == 0
                @test lock!(wlock) == nothing
                @test rwlock.writer == true
                @test rwlock.readers == 0
                @test take!(c) == :pretest
                put!(c, :prelock)
                @test lock!(wlock) == nothing

                # this code should never be reached
                @test take!(c) == :prelock
                put!(c, :postlock)
            end

            sleep(1)

            @test take!(c) == :prelock
            put!(c, :posttest)

            @test rwlock.writer == true
            @test rwlock.readers == 0
        end

        @testset "unlock" begin
            rwlock = ReadWriteLock()
            wlock = write_lock(rwlock)

            c = Channel{Symbol}(1)
            put!(c, :pretest)

            @async begin
                @test rwlock.writer == false
                @test rwlock.readers == 0
                @test lock!(wlock) == nothing
                @test rwlock.writer == true
                @test rwlock.readers == 0
                @test take!(c) == :pretest
                put!(c, :preunlock)
                @test unlock!(wlock) == nothing
                @test take!(c) == :preunlock
                put!(c, :postunlock)
            end

            sleep(1)

            @test take!(c) == :postunlock
            put!(c, :posttest)

            @test rwlock.writer == false
            @test rwlock.readers == 0
        end
    end

    @testset "read and write locks" begin
        @testset "write then read" begin
            rwlock = ReadWriteLock()
            wlock = write_lock(rwlock)
            rlock = read_lock(rwlock)

            c = Channel{Symbol}(1)
            put!(c, :pretest)

            @async begin
                @test rwlock.writer == false
                @test rwlock.readers == 0
                @test lock!(wlock) == nothing
                @test rwlock.writer == true
                @test rwlock.readers == 0
                @test take!(c) == :pretest
                put!(c, :prelock)
                @test lock!(rlock) == nothing

                # this code should never be reached
                @test take!(c) == :prelock
                put!(c, :postlock)
            end

            sleep(1)

            @test take!(c) == :prelock
            put!(c, :posttest)

            @test rwlock.writer == true
            @test rwlock.readers == 0
        end

        @testset "read then write" begin
            rwlock = ReadWriteLock()
            wlock = write_lock(rwlock)
            rlock = read_lock(rwlock)

            c = Channel{Symbol}(1)
            put!(c, :pretest)

            @async begin
                @test rwlock.writer == false
                @test rwlock.readers == 0
                @test lock!(rlock) == nothing
                @test rwlock.writer == false
                @test rwlock.readers == 1
                @test take!(c) == :pretest
                put!(c, :prelock)
                @test lock!(wlock) == nothing

                # this code should never be reached
                @test take!(c) == :prelock
                put!(c, :postlock)
            end

            sleep(1)

            @test take!(c) == :prelock
            put!(c, :posttest)

            @test rwlock.writer == false
            @test rwlock.readers == 1
        end
    end
end  # Two-threaded

end  # ReadWriteLocks


text:
"""
    EMAFn(x,m=1)

Equally weighted moving average (over the current and m-1 lagged values) for a matrix `x`.

"""
function EMAFn(x,m=1)

    (T,n) = (size(x,1),size(x,2))
    q = m - 1                               #m=2, q=1 (avg over [t-1,t])

    (q < 0) && error("q must be >= 0")

    y = fill(NaN,T,n)                       #y[t] = (x[t-1] + x[t])/2
    for t = 1:T
        vv     = max(1,t-q):t                 #with q=2; 1:1,1:2,1:3,2:4
        y[t,:] = sum(x[vv,:],dims=1)/m
    end

    return y

end

max_stars_repo_name: duzaichuan/FinancialEconometrics

max_issues_repo_name: duzaichuan/FinancialEconometrics

max_forks_repo_name: duzaichuan/FinancialEconometrics


text:
export Symp, isGeneric, nModes, ⊗, dsum, ⊕
export Omega, Ω

struct Symp{T} <: AbstractMatrix{T}
    S:: AbstractMatrix{T}
    function Symp{T}(S::AbstractMatrix{T}) where T
        if !(T<:Real) return error("not real") end
        n = LinearAlgebra.checksquare(S) ÷ 2
        Ω = cat(fill([0 1; -1 0], n)...; dims=(1,2))
        if (norm(S*Ω*S' - Ω) >= tolerance) return error("not symplectic") end
        return new{T}(S)
    end
end

# Define AbstractMatrix Interface `size`
Base.size(S::Symp, dim::Integer)::Integer = size(S.S, dim)
Base.size(S::Symp) = size(S.S)
# Define AbstractMatrix Interface `getindex`
function Base.getindex(S::Symp{T}, i::Integer, j::Integer)::T where {T}
    return S.S[i, j]
end
# Type Conversions
Symp(S::AbstractMatrix) = Symp{typeof(S).parameters[1]}(S)
Symp(S::Symp) = S
Base.Matrix(S::Symp)::AbstractMatrix = S.S
Base.convert(::Type{Symp}, S::AbstractMatrix) = Symp(S)
Base.convert(::Type{Symp}, S::Symp) = S
Base.convert(::Type{AbstractMatrix}, S::Symp) = Matrix(S)

# Fundamental Constructors
# From Number of Modes to Symplectic Form
Omega(n::Int)::Symp = cat(fill([0 1; -1 0], n)...; dims=(1,2))
Omega(S::Symp)::Symp = Omega(size(S)[1] ÷ 2)
Ω = Omega

# Single Mode Gaussian Operation

# Properties
isSquare = LinearAlgebra.checksquare
# isGeneric(m::AbstractMatrix)::Bool = all(x -> abs(x) > tolerance, m )
isGeneric(m::AbstractMatrix)::Bool = true
isGeneric(S::Symp)::Bool = isGeneric(S.S)
nModes(S::Symp)::Int = size(S)[1] ÷ 2       # Number of Modes

# Operators
Base.transpose(S::Symp)::Symp = Base.transpose( S.S )
Base.adjoint(S::Symp)::Symp = Base.adjoint( S.S )
# Unary Minus and Unary Plus
Base.:-(S::Symp)::Symp = Base.:-( S.S )
Base.:+(S::Symp)::Symp = Base.:+( S.S )
# Inverse
Base.inv(S::Symp)::Symp = - Ω(S) * S' * Ω(S)
Base.:*(S1::Symp, S2::Symp, Ss::Vararg{Symp})::Symp = *(S1.S, S2.S,  [s.S for s in Ss]... )
⊗ = Base.kron
# Direct Sum
dsum(M1::AbstractMatrix, M2::AbstractMatrix)::AbstractMatrix = cat(M1, M2; dims=(1,2))
dsum(M1::AbstractMatrix, M2::AbstractMatrix, Ms::Vararg{AbstractMatrix})::AbstractMatrix = dsum(dsum(M1, M2),  mapreduce(identity, dsum, Ms ))
dsum(S1::Symp, S2::Symp, Ss::Vararg{Symp})::Symp = dsum(S1.S, S2.S,  [s.S for s in Ss]... )
⊕ = dsum
Base.:^(S::Symp, p::Int)::Symp = p < 0 ? inv(S).S ^ (-p)  : S.S^p


text:
# Autogenerated wrapper script for cryptopp_jll for x86_64-apple-darwin
export libcryptopp

JLLWrappers.@generate_wrapper_header("cryptopp")
JLLWrappers.@declare_library_product(libcryptopp, "@rpath/libcryptopp.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcryptopp,
        "lib/libcryptopp.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()