text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

const BATCH_WINDOW_SMALL_MILLIS = 1
const BATCH_WINDOW_LARGE_MILLIS = 3

const DISPATCHER_THREADS = 3

mutable struct Dispatcher
    input_channel::Channel{Proto.Transaction_Client}
    direct_dispatch_channel::Channel{Proto.ProtoType}
    dispatch_channel::Channel{Proto.ProtoType}
    dispatch_timer::Optional{Vector{Timer}}
end

function Dispatcher(input_channel::Channel{Proto.Transaction_Client})
    direct_dispatch_channel = Channel{Proto.ProtoType}(10)
    dispatch_channel = Channel{Proto.ProtoType}(10)

    disp_timer = Vector{Timer}()
    for _ in 1:DISPATCHER_THREADS
        push!(disp_timer, batch_requests(dispatch_channel,input_channel))
    end

    Threads.@spawn process_direct_requests(direct_dispatch_channel, input_channel)

    return Dispatcher(input_channel, direct_dispatch_channel, dispatch_channel,disp_timer)
end

"""
function process_direct_requests(in_channel::Channel{Proto.ProtoType}, out_channel::Channel{Proto.Transaction_Client})
    This function process the incoming request directly to the server
"""
function process_direct_requests(in_channel::Channel{Proto.ProtoType}, out_channel::Channel{Proto.Transaction_Client})
    for tmp_res in in_channel
        client_res = TransactionRequestBuilder.client_msg([tmp_res])
        put!(out_channel, client_res)
    end
    @debug "process_direct_requests was closed"
end

"""
function batch_requests(in_channel::Channel{Proto.ProtoType}, out_channel::Channel{Proto.Transaction_Client})
    This function contains the whole logic for batching incomming requests. The inner runner function will be
    called every x ms and will send the collected request in one Transaction_Client message to teh server.
"""
function batch_requests(in_channel::Channel{Proto.ProtoType}, out_channel::Channel{Proto.Transaction_Client})
    time_to_run = BATCH_WINDOW_SMALL_MILLIS / 1000

    function runner(controller)
        @async sleeper(controller)
        answers = Proto.Transaction_Req[]
        try
            while controller.running
                yield()
                if isready(in_channel)
                    tmp_res = take!(in_channel)
                    push!(answers, tmp_res)
                end
            end
            if length(answers) > 0
                @async put!(out_channel, TransactionRequestBuilder.client_msg(answers))
            end
        catch ex
            throw(TypeDBClientException("batch_requests runner failure",ex))
        end
    end

    cb(timer) = (runner(Controller(true, time_to_run)))
    t = Timer(cb,time_to_run, interval = time_to_run)
    wait(t)

    return t
end

function close(dispatcher::Dispatcher)
    safe_close(dispatcher.direct_dispatch_channel)
    safe_close(dispatcher.dispatch_channel)
    safe_close(dispatcher.dispatch_timer)
end


text:
import Test: @test, @testset, @test_throws
import Luna: Output

@testset "Linear" begin
    # compare radial single pol, to radial linear pol at 45 degrees,
    # in capillary (non birefringent) these should be identical

    using Luna
    import LinearAlgebra: norm
    a = 13e-6
    gas = :Ar
    pres = 5
    τ = 30e-15
    λ0 = 800e-9
    energy = 1e-6
    grid = Grid.RealGrid(5e-2, 800e-9, (160e-9, 3000e-9), 1e-12)

    dens0 = PhysData.density(gas, pres)
    densityfun(z) = dens0
    responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
    energyfun, energyfunω = Fields.energyfuncs(grid)
    modes = (
         Capillary.MarcatiliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false),
    )
    inputs = Fields.GaussField(λ0=λ0, τfwhm=τ, energy=energy)
    Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs,
                                modes, :y; full=true)
    statsfun = Stats.collect_stats(grid, Eω,
                                Stats.ω0(grid),
                                Stats.peakintensity(grid, modes),
                                Stats.fwhm_r(grid, modes),
                                Stats.energy(grid, energyfunω))
    output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
    linop = LinearOps.make_const_linop(grid, modes, λ0)
    Luna.run(Eω, grid, linop, transform, FT, output, status_period=10)

    modes = (
        Capillary.MarcatiliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=0.0, loss=false),
        Capillary.MarcatiliMode(a, gas, pres, n=1, m=1, kind=:HE, ϕ=π/2, loss=false)
    )
    inf = (Fields.GaussField(λ0=λ0, τfwhm=τ, energy=energy/2.0),)
    # same field in each mode
    inputs = ((mode=1, fields=inf), (mode=2, fields=inf))
    Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs,
                                modes, :xy; full=true)
    statsfun = Stats.collect_stats(grid, Eω,
                                Stats.ω0(grid),
                                Stats.peakintensity(grid, modes, components=:xy),
                                Stats.fwhm_r(grid, modes, components=:xy),
                                Stats.energy(grid, energyfunω))
    outputp = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
    linop = LinearOps.make_const_linop(grid, modes, λ0)
    Luna.run(Eω, grid, linop, transform, FT, outputp, status_period=10)

    Iω = dropdims(abs2.(output.data["Eω"]), dims=2)
    Iωp = dropdims(sum(abs2.(outputp.data["Eω"]), dims=2), dims=2)

    @test norm(Iω - Iωp)/norm(Iω) < 1.07e-12
    @test all(output["stats"]["peakintensity"] .≈ outputp["stats"]["peakintensity"])
    @test all(output["stats"]["energy"] .≈ sum(outputp["stats"]["energy"], dims=1))
    @test all(output["stats"]["fwhm_r"] .≈ outputp["stats"]["fwhm_r"])
end

##
@testset "Linear, LP11" begin
    using Luna
    import LinearAlgebra: norm
    a = 125e-6
    gas = :Ar
    pres = 0.167
    flength = 0.1

    τ = 10e-15
    λ0 = 800e-9
    energy = 150e-6
    grid = Grid.RealGrid(5e-2, 800e-9, (160e-9, 3000e-9), 1e-12)

    densityfun = let dens0=PhysData.density(gas, pres)
        z -> dens0
    end
    responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
    energyfun, energyfunω = Fields.energyfuncs(grid)
    modes = (
        Capillary.MarcatiliMode(a, gas, pres, n=0, m=1, kind=:TM, ϕ=0.0, loss=false),
        Capillary.MarcatiliMode(a, gas, pres, n=2, m=1, kind=:HE, ϕ=π/4, loss=false),
    )
    inf = (Fields.GaussField(λ0=λ0, τfwhm=τ, energy=energy/2.0),)
    # same field in each mode
    inputs = ((mode=1, fields=inf), (mode=2, fields=inf))
    Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs,
                                modes, :xy; full=true)
    statsfun = Stats.collect_stats(grid, Eω,
                                Stats.ω0(grid),
                                Stats.peakintensity(grid, modes),
                                Stats.fwhm_r(grid, modes),
                                Stats.energy(grid, energyfunω))
    output = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
    linop = LinearOps.make_const_linop(grid, modes, λ0)
    Luna.run(Eω, grid, linop, transform, FT, output, status_period=10)

    modes = (
        Capillary.MarcatiliMode(a, gas, pres, n=0, m=1, kind=:TM, ϕ=0.0, loss=false),
        Capillary.MarcatiliMode(a, gas, pres, n=2, m=1, kind=:HE, ϕ=-π/4, loss=false),
    )
    inf = (Fields.GaussField(λ0=λ0, τfwhm=τ, energy=energy/2.0),)
    # same field in each mode
    inputs = ((mode=1, fields=inf), (mode=2, fields=inf))
    Eω, transform, FT = Luna.setup(grid, densityfun, responses, inputs,
                                modes, :xy; full=true)
    statsfun = Stats.collect_stats(grid, Eω,
                                Stats.ω0(grid),
                                Stats.peakintensity(grid, modes, components=:xy),
                                Stats.fwhm_r(grid, modes, components=:xy),
                                Stats.energy(grid, energyfunω))
    outputp = Output.MemoryOutput(0, grid.zmax, 201, statsfun)
    linop = LinearOps.make_const_linop(grid, modes, λ0)
    Luna.run(Eω, grid, linop, transform, FT, outputp, status_period=10)

    Iω = dropdims(sum(abs2.(output.data["Eω"]), dims=2), dims=2)
    Iωp = dropdims(sum(abs2.(outputp.data["Eω"]), dims=2), dims=2)

    @test norm(Iω - Iωp)/norm(Iω) < 1.07e-12
    @test all(output["stats"]["peakintensity"] .≈ outputp["stats"]["peakintensity"])
    @test all(sum(output["stats"]["energy"], dims=1) .≈ sum(outputp["stats"]["energy"], dims=1))
    @test all(output["stats"]["fwhm_r"] .≈ outputp["stats"]["fwhm_r"])
end

ext: jl


text:
export TabularCFRPolicy

struct InfoStateNode{M<:AbstractVector{Bool}}
    strategy::Vector{Float64}
    cumulative_regret::Vector{Float64}
    cumulative_strategy::Vector{Float64}
    mask::M
end

function InfoStateNode(mask)
    n = sum(mask)
    InfoStateNode(
        fill(1 / n, n),
        zeros(n),
        zeros(n),
        mask
    )
end

#####
# TabularCFRPolicy
#####

mutable struct TabularCFRPolicy{S,T,R<:AbstractRNG} <: AbstractCFRPolicy
    nodes::Dict{S,InfoStateNode}
    behavior_policy::QBasedPolicy{TabularLearner{S,T},WeightedExplorer{true,R}}
    is_reset_neg_regrets::Bool
    is_linear_averaging::Bool
    weighted_averaging_delay::Int
    is_alternating_update::Bool
    rng::R
    n_iteration::Int
end

(p::TabularCFRPolicy)(env::AbstractEnv) = p.behavior_policy(env)

RLBase.prob(p::TabularCFRPolicy, env::AbstractEnv) = prob(p.behavior_policy, env)

"""
    TabularCFRPolicy(;kwargs...)

Some useful papers while implementing this algorithm:

- [An Introduction to Counterfactual Regret Minimization](http://modelai.gettysburg.edu/2013/cfr/cfr.pdf)
- [MONTE CARLO SAMPLING AND REGRET MINIMIZATION FOR EQUILIBRIUM COMPUTATION AND DECISION-MAKING IN LARGE EXTENSIVE FORM GAMES](http://mlanctot.info/files/papers/PhD_Thesis_MarcLanctot.pdf)
- [Solving Large Imperfect Information Games Using CFR⁺](https://arxiv.org/pdf/1407.5042.pdf)
- [Revisiting CFR⁺ and Alternating Updates](https://arxiv.org/pdf/1810.11542v1.pdf)
- [Solving Imperfect-Information Games via Discounted Regret Minimization](https://arxiv.org/pdf/1809.04040.pdf)

# Keyword Arguments

- `is_alternating_update=true`: If `true`, we update the players alternatively.
- `is_reset_neg_regrets=true`: Whether to use **regret matching⁺**.
- `is_linear_averaging=true`
- `weighted_averaging_delay=0`. The averaging delay in number of iterations. Only valid when `is_linear_averaging` is set to `true`.
- `state_type=String`, the data type of information set.
- `rng=Random.GLOBAL_RNG`
"""
function TabularCFRPolicy(;
    is_reset_neg_regrets = true,
    is_linear_averaging = true,
    weighted_averaging_delay = 0,
    is_alternating_update = true,
    state_type = String,
    rng = Random.GLOBAL_RNG,
    n_iteration = 1,
)
    TabularCFRPolicy(
        Dict{state_type,InfoStateNode}(),
        TabularRandomPolicy(;
            rng = rng,
            table = Dict{state_type,Vector{Float64}}(),
            is_normalized = true,
        ),
        is_reset_neg_regrets,
        is_linear_averaging,
        weighted_averaging_delay,
        is_alternating_update,
        rng,
        n_iteration,
    )
end

"Update the `behavior_policy`"
function RLBase.update!(p::TabularCFRPolicy)
    for (k, v) in p.nodes
        s = sum(v.cumulative_strategy)
        if s != 0
            m = v.mask
            strategy = zeros(length(m))
            strategy[m] .= v.cumulative_strategy ./ s
            update!(p.behavior_policy, k => strategy)
        else
            # The TabularLearner will return uniform distribution by default.
            # So we do nothing here.
        end
    end
end

"Run one interation"
function RLBase.update!(p::TabularCFRPolicy, env::AbstractEnv)
    w = p.is_linear_averaging ? max(p.n_iteration - p.weighted_averaging_delay, 0) : 1
    if p.is_alternating_update
        for x in players(env)
            if x != chance_player(env)
                cfr!(p.nodes, env, x, w)
                regret_matching!(p)
            end
        end
    else
        cfr!(p.nodes, env, nothing, w)
        regret_matching!(p)
    end
    p.n_iteration += 1
end

"""
Symbol meanings:

π: reach prob
π′: new reach prob
π₋ᵢ: opponents' reach prob
p: player to update. `nothing` means simultaneous update.
w: weight
v: counterfactual value **before weighted by opponent's reaching probability**
V: a vector containing the `v` after taking each action with current information set. Used to calculate the **regret value**
"""
function cfr!(nodes, env, p, w, π = Dict(x => 1.0 for x in players(env)))
    if is_terminated(env)
        reward(env, p)
    else
        if current_player(env) == chance_player(env)
            v = 0.0
            for (a, pₐ) in zip(action_space(env), prob(env))
                if pₐ > 0
                    π′ = copy(π)
                    π′[current_player(env)] *= pₐ
                    v += pₐ * cfr!(nodes, child(env, a), p, w, π′)
                end
            end
            v
        else
            v = 0.0
            legal_actions = legal_action_space(env)
            M = legal_action_space_mask(env)
            node = get!(nodes, state(env), InfoStateNode(M))

            is_update = isnothing(p) || p == current_player(env)
            V = is_update ? Vector{Float64}(undef, length(legal_actions)) : nothing

            for (i, a) in enumerate(legal_actions)
                σᵢ = node.strategy[i]
                π′ = copy(π)
                π′[current_player(env)] *= σᵢ

                vₐ = cfr!(nodes, child(env, a), p, w, π′)
                is_update && (V[i] = vₐ)
                v += σᵢ * vₐ
            end

            if is_update
                πᵢ = π[current_player(env)]
                π₋ᵢ = reduce(*, (prob for (p, prob) in π if p != current_player(env)))
                node.cumulative_regret .+= π₋ᵢ .* (V .- v)
                node.cumulative_strategy .+= w .* πᵢ .* node.strategy
            end
            v
        end
    end
end

function regret_matching!(p::TabularCFRPolicy)
    for node in values(p.nodes)
        regret_matching!(node; is_reset_neg_regrets = p.is_reset_neg_regrets)
    end
end

regret_matching!(node::InfoStateNode; kwargs...) =
    regret_matching!(node.strategy, node.cumulative_regret; kwargs...)

function regret_matching!(strategy, cumulative_regret; is_reset_neg_regrets = true)
    if is_reset_neg_regrets
        for i in 1:length(cumulative_regret)
            if cumulative_regret[i] < 0
                cumulative_regret[i] = 0
            end
        end
    end
    s = mapreduce(x -> max(0, x), +, cumulative_regret)
    if s > 0
        strategy .= max.(0.0, cumulative_regret) ./ s
    else
        fill!(strategy, 1 / length(strategy))
    end
end


text:
"""
   initial_cube(sim, l=100u"nm"; species="electron")

Select a cube of lenght `l` of the given `species` at the start
of the simulation `sim`.
"""
function initial_cube(sim, l=100u"nm"; species="electron")
    file = sim[1]
    r = file["grid/$species"]

    findall(rᵢ->rᵢ[1] ∈ zero(l)..l &&
            rᵢ[2] ∈ -l/2..l/2 &&
            rᵢ[3] ∈ -l/2..l/2, r)
end

function indexof(apparitions, haystack, n)
    found = zeros(Int, n)
    k = 1
    @inbounds for i in haystack
        if apparitions[i]
            found[k] = i
            k += 1
        end
    end

    return found
end

"""
    track_particles(sim, idxs; species="electron")

Track the particles with the initial index `idx` and species given
by the `species` argument in the simulation `sim`.
"""
function track_particles(sim, idxs; species="electron")
    n_total = get_npart(sim[1], species)
    n_searched = length(idxs)
    apparitions = falses(n_total)
    apparitions[idxs] .= true

    map(eachindex(sim.files)) do i
        file = sim[i]
        @time ids, r = file["id/$species", "grid/$species"]
        @debug "Loaded $(file.name)"
        ids = Int.(ids.data)
        if issorted(ids) && length(ids) == length(idxs)
            box_ids = ids
        else
            @time box_ids = indexof(apparitions, ids, n_searched)
        end
        filter!(!iszero, box_ids)

        @debug "Done with $i"
        SVector{3}.(r[1][box_ids], r[2][box_ids], r[3][box_ids])
    end
end


text:
"""
This is the code used to randomly generate a tree to demonstrate the plotting functionality of PhyloTrees.jl in the README.
"""

using PhyloTrees, Plots, StatsBase

g = Tree()
addnode!(g)
for i = 1:100
  branch!(g, sample(collect(keys(g.nodes))), 10*rand())
end
plot(g, dpi=200)
png(joinpath(@__DIR__, "treeplot.png"))


text:
int main() {
    boolean b=false;
    if(b) return 0;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
module Misc

	using Pkg

	const DATASETS_FOLDER = joinpath(dirname(dirname(@__DIR__)), "datasets")
	const IMG_LENGTH = 784

end


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

using JuliaFEM
using JuliaFEM.Testing
using Logging
Logging.configure(level=DEBUG)

@testset "create and manipulate fields" begin

@testset "updating time dependent fields" begin
    f = Field(0.0 => 1.0)
    @test last(f).time == 0.0
    @test last(f).data == 1.0
    update!(f, 0.0 => 2.0)
    @test last(f).time == 0.0
    @test last(f).data == 2.0
    @test length(f) == 1
    update!(f, 1.0 => 3.0)
    @test last(f).time == 1.0
    @test last(f).data == 3.0
    @test length(f) == 2
end

@testset "updating time invariant fields" begin
    f = Field(1.0)
    @test f.data == 1.0
    update!(f, 2.0)
    @test f.data == 2.0
end

@testset "field defined using function" begin
    g(xi, t) = xi[1]*t
    f = Field(g)
    v = f([1.0], 2.0)
    @test isapprox(v, 2.0)
end

@testset "dictionary fields" begin
    f1 = Dict{Int64, Vector{Float64}}(1 => [0.0, 0.0], 2 => [0.0, 0.0])
    f2 = Dict{Int64, Vector{Float64}}(1 => [1.0, 1.0], 2 => [1.0, 1.0])
    f = Field(0.0 => f1, 1.0 => f2)
    debug("field = $f")
    @test isa(f, DVTV)
    @test isapprox(f(0.0)[1], [0.0, 0.0])
    @test isapprox(f(1.0)[2], [1.0, 1.0])

    f = Field(0.0 => f1)
    update!(f, 1.0 => f2)
    @test isa(f, DVTV)
    @test isapprox(f(0.0)[1], [0.0, 0.0])
    @test isapprox(f(1.0)[2], [1.0, 1.0])

    f = Field(f1)
    @test isapprox(f(0.0)[1], [0.0, 0.0])
    @test isapprox(f[1], [0.0, 0.0])

    f = Field(f1)
    @test isa(f, DVTI)
end

end


text:
#all kernel functions are normalized such that their support radius is h

@fastmath function pos(x::Float64)::Float64
	return (x > 0.0 ? x : 0.0)
end

"""
    spline23(h::Float64, r::Float64)::Float64

Returns ``w(r)``, the value of a 2d cubic spline ``w`` with support radius `h`.
Integrates to unity.

"""
@fastmath function spline23(h::Float64, r::Float64)::Float64
	x = r/h
	if (x < 0.5)
		# 40/7pi = 1.8189136353359467
		return 1.8189136353359467*(1.0 - 6.0*x^2 + 6.0*x^3)/h^2
	elseif (x < 1.0)
		# 80/7pi = 3.6378272706718935
		return 3.6378272706718935*(1.0 - x)^3/h^2
	end
	return 0.0
end

"""
    Dspline23(h::Float64, r::Float64)::Float64

Returns ``\\frac{\\text{d}w}{\\text{d}r}(r)``, the first derivative of a 2d cubic spline ``w`` with support radius `h`.
Integrates to unity.

"""
@fastmath function Dspline23(h::Float64, r::Float64)::Float64
	x = r/h
	if (x < 0.5)
		# 240/7pi = 10.91348181201568
		return -10.91348181201568*(2.0*x - 3.0*x^2)/h^3
	elseif (x < 1.0)
		return -10.91348181201568*(1.0 - x)^2/h^3
	end
	return 0.0
end

"""
    rDspline23(h::Float64, r::Float64)::Float64

Returns ``\\frac{1}{r}\\,\\frac{\\text{d}w}{\\text{d}r}(r)``, the reduced first derivative of a 2d cubic spline ``w`` with support radius `h`.
Integrates to unity.

"""
@fastmath function rDspline23(h::Float64, r::Float64)::Float64
	x = r/h
	if (x < 0.5)
		# 240/7pi = 10.91348181201568
		return -10.91348181201568*(2.0 - 3.0*x)/h^4
	elseif (x < 1.0)
		return -10.91348181201568*(1.0 - x)^2/(x*h^4)
	end
	return 0.0
end

"""
    spline24(h::Float64, r::Float64)::Float64

Returns ``w(r)``, the value of a 2d quartic spline ``w`` with support radius `h`.
Integrates to unity.

"""
@fastmath function spline24(h::Float64, r::Float64)::Float64
	x = r/h
	return 6.222175110452539*(pos(1.0 - x)^4 - 5*pos(0.6 - x)^4 + 10*pos(0.2 - x)^4)/h^2
end

"""
    Dspline24(h::Float64, r::Float64)::Float64

Returns ``\\frac{\\text{d}w}{\\text{d}r}(r)``, the first derivative of a 2d quartic spline ``w`` with support radius `h`.
Integrates to unity.

"""
@fastmath function Dspline24(h::Float64, r::Float64)::Float64
	x = r/h
	return -24.888700441810155*(pos(1.0 - x)^3 - 5*pos(0.6 - x)^3 + 10*pos(0.2 - x)^3)/h^3
end

"""
    rDspline24(h::Float64, r::Float64)::Float64

Returns ``\\frac{1}{r}\\,\\frac{\\text{d}w}{\\text{d}r}(r)``, the reduced first derivative of a 2d quartic spline ``w`` with support radius `h`.
Integrates to unity.

"""
@fastmath function rDspline24(h::Float64, r::Float64)::Float64
	x = r/h
	if x > 0.2
		return -24.888700441810155*(pos(1.0 - x)^3 - 5*pos(0.6 - x)^3)/(x*h^4)
	end
	return -24.888700441810155*(1.2 - 6.0*x^2)/h^4
end

"""
    wendland2(h::Float64, r::Float64)::Float64

Returns ``w(r)``, the value of a 2d quintic Wendland kernel ``w`` with support radius `h`.
Integrates to unity.

"""
@fastmath function wendland2(h::Float64, r::Float64)::Float64
	x = r/h
	if (x > 1.0)
		return 0.0
	end
	# 7/pi ≐ 2.228169203286535
	return 2.228169203286535*((1.0 - x)^4)*(1.0 + 4.0*x)/h^2
end

"""
    Dwendland2(h::Float64, r::Float64)::Float64

Returns ``\\frac{\\text{d}w}{\\text{d}r}(r)``, the first derivative of a 2d quintic Wendland kernel ``w`` with support radius `h`.
Integrates to unity.

"""
@fastmath function Dwendland2(h::Float64, r::Float64)::Float64
	x = r/h
	if (x > 1.0)
		return 0.0
	end
	# 140/pi ≐ 44.563384065730695
	return -44.563384065730695*x*((1.0 - x)^3)/h^3
end

"""
    rDwendland2(h::Float64, r::Float64)::Float64

Returns ``\\frac{1}{r}\\,\\frac{\\text{d}w}{\\text{d}r}(r)``, the reduced first derivative of a 2d quintic Wendland kernel ``w`` with support radius `h`.
Integrates to unity.

"""
@fastmath function rDwendland2(h::Float64, r::Float64)::Float64
	x = r/h
	if (x > 1.0)
		return 0.0
	end
	# 140/pi ≐ 44.563384065730695
	return -44.563384065730695*((1.0 - x)^3)/h^4
end

"""
    wendland3(h::Float64, r::Float64)::Float64

Returns ``w(r)``, the value of a 3d quintic Wendland kernel ``w`` with support radius `h`.
Integrates to unity.

"""
@fastmath function wendland3(h::Float64, r::Float64)::Float64
	x = r/h
	if (x > 1.0)
		return 0.0
	end
	# 21/2pi ≐ 3.3422538049298023
	return 3.3422538049298023*((1.0 - x)^4)*(1.0 + 4.0*x)/h^3
end

"""
    Dwendland3(h::Float64, r::Float64)::Float64

Returns ``\\frac{\\text{d}w}{\\text{d}r}(r)``, the first derivative of a 3d quintic Wendland kernel ``w`` with support radius `h`.
Integrates to unity.

"""
@fastmath function Dwendland3(h::Float64, r::Float64)::Float64
	x = r/h
	if (x > 1.0)
		return 0.0
	end
	# 210/pi ≐ 66.84507609859604
	return -66.84507609859604*x*((1.0 - x)^3)/h^4
end

"""
    rDwendland3(h::Float64, r::Float64)::Float64

Returns ``\\frac{1}{r}\\,\\frac{\\text{d}w}{\\text{d}r}(r)``, the reduced first derivative of a 3d quintic Wendland kernel ``w`` with support radius `h`.
Integrates to unity.

"""
@fastmath function rDwendland3(h::Float64, r::Float64)::Float64
	x = r/h
	if (x > 1.0)
		return 0.0
	end
	# 210/pi ≐ 66.84507609859604
	return -66.84507609859604*((1.0 - x)^3)/h^5
end

@fastmath function DDwendland3(h::Float64, r::Float64)::Float64
	x = r/h
	if (x > 1.0)
		return 0.0
	end
	# 210/pi ≐ 66.84507609859604
	return -66.84507609859604*((1.0 - 4.0*x)*(1.0 - x)^2)/h^5
end

@fastmath function spiky3(h::Float64, r::Float64)::Float64
	x = r/h
	#15/2pi = 2.3873241463784303
	return 2.3873241463784303*(1.0 - x)^2/h^3
end

@fastmath function Dspiky3(h::Float64, r::Float64)::Float64
	x = r/h
	#15/pi = 4.7746482927568605
	return -(x < 1.0)*2.3873241463784303*(1.0 - x)/h^4
end

@fastmath function rDspiky3(h::Float64, r::Float64)::Float64
	x = r/h
	#15/pi = 4.7746482927568605
	return -(x < 1.0)*2.3873241463784303*(1.0 - x)/((x + 1e-6)*h^5)
end

@fastmath function DDspiky3(h::Float64, r::Float64)::Float64
	x = r/h
	#15/pi = 4.7746482927568605
	return (x < 1.0)*2.3873241463784303/h^5
end


text:
# ---
# title: 856. Score of Parentheses
# id: problem856
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String, Stack
# link: <https://leetcode.com/problems/score-of-parentheses/description/>
# hidden: true
# ---
#
# Given a balanced parentheses string `S`, compute the score of the string based
# on the following rule:
#
#   * `()` has score 1
#   * `AB` has score `A + B`, where A and B are balanced parentheses strings.
#   * `(A)` has score `2 * A`, where A is a balanced parentheses string.
#
#
#
# **Example 1:**
#
#
#
#     Input: "()"
#     Output: 1
#
#
# **Example 2:**
#
#
#
#     Input: "(())"
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: "()()"
#     Output: 2
#
#
# **Example 4:**
#
#
#
#     Input: "(()(()))"
#     Output: 6
#
#
#
#
# **Note:**
#
#   1. `S` is a balanced parentheses string, containing only `(` and `)`.
#   2. `2 <= S.length <= 50`
#
#
## @lc code=start
using LeetCode

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


text:
import Base: getindex, eltype, convert, length
import Distances: evaluate, colwise
import NearestNeighbors: inrange

### State Typedefs
export AbstractState, State, Path, AbstractBitsState, BitsState
export statevec2mat, statemat2vec, changeprecision, dense

"The `supertype` for all non-`AbstractVector` states."
abstract AbstractState
"Type union of all states, encompassing `AbstractVector`s and `AbstractState`s."
typealias State Union{AbstractVector, AbstractState}
"Shorthand for a list (`Vector`) of states."
typealias Path{S<:State} Vector{S}
"The `supertype` for all non-`SVector` states that satisfy `isbits(S) == true`."
abstract AbstractBitsState <: AbstractState
"Type union encompassing all states that may satisfy `isbits(S) == true`."
typealias BitsState Union{SVector, FieldVector, AbstractBitsState}

"Reinterprets a list of `BitsState`s as a `Matrix` without copying/re-allocating memory."    # TODO: methods for nonbits
function statevec2mat{S<:BitsState}(V::Vector{S})
    isbits(S) ? reinterpret(eltype(S), V, (length(S), length(V))) : error("$S must satisfy isbits($S) = true")
end
statevec2mat{S<:BitsState}(V::S...) = statevec2mat(collect(V))
statevec2mat{N,S<:BitsState}(V::NTuple{N,S}) = statevec2mat(collect(V))
"Reinterprets a `Matrix` as a column-wise list of `BitsState`s without copying/re-allocating memory."
function statemat2vec{S<:BitsState,T}(::Type{S}, M::Matrix{T})
    if isbits(S)
        length(S) == size(M,1) ? reinterpret(S, M, (size(M,2),)) : error("$S size mismatch with Matrix")
    else
        error("$S must satisfy isbits($S) = true")
    end
end
"Convert `x` to a similar type with numerical precision `T <: AbstractFloat` (e.g. `Float64`, `Float32`, etc.)."
changeprecision{T<:AbstractFloat,S}(::Type{T}, x::S) = convert(changeprecision(T,S), x)
changeprecision{T<:AbstractFloat,S,N}(::Type{T}, ::Type{Array{S,N}}) = Array{changeprecision(T,S),N}
changeprecision{T<:AbstractFloat,S<:Real,N}(::Type{T}, ::Type{Array{S,N}}) = Array{T,N}
changeprecision{T<:AbstractFloat,N,S<:Real}(::Type{T}, ::Type{SVector{N,S}}) = SVector{N,T}
changeprecision{T<:AbstractFloat,N1,N2,S<:Real,L}(::Type{T}, ::Type{SMatrix{N1,N2,S,L}}) = SMatrix{N1,N2,T,L}
changeprecision{T<:AbstractFloat,S<:Real}(::Type{T}, ::Type{S}) = T
changeprecision{T<:AbstractFloat,S}(::Type{T}, ::Type{S}) = S

## SE2State
export SE2State

"State consisting of a position (x,y) ∈ R^2 and an angle θ ∈ R."
immutable SE2State{T} <: FieldVector{T}
    x::T
    y::T
    θ::T
end
changeprecision{T<:AbstractFloat,S<:Real}(::Type{T}, ::Type{SE2State{S}}) = SE2State{T}

### StateSpace Typedefs
export StateSpace, State2Workspace

"The `supertype` for all state spaces; the type parameter is the state type."
abstract StateSpace{S<:State}
"Encodes a transformation from the state space (dynamics) to the workspace (obstacles)."
abstract State2Workspace

### Metrics and QuasiMetrics
export QuasiMetric, PreMetric, Metric, ChoppedMetric, ChoppedQuasiMetric, ChoppedPreMetric

"A quasimetric satisfies positivity, positive definiteness, and the triangle inequality (no symmetry)."
abstract QuasiMetric <: PreMetric

## Fallbacks (should usually be extended)
colwise{S<:State}(d::PreMetric, v::S, W::Vector{S}) = eltype(S)[evaluate(d, v, w) for w in W]
colwise{S<:State}(d::PreMetric, W::Vector{S}, v::S) = eltype(S)[evaluate(d, w, v) for w in W]

## Chopped, Lower-Bounded Metrics
"""
Evaluates as `m(v, w) <= chopval ? m(v, w) : Inf` for points `v` and `w`.\n
`lowerbound` should be easier to evaluate than `m` and must satisfy `lowerbound(v, w) ≤ m(v, w)` for all `v`, `w`.
"""
type ChoppedMetric{M<:Metric,B<:Metric,T<:AbstractFloat} <: PreMetric
    m::M
    lowerbound::B
    chopval::T
end
"""
Evaluates as `m(v, w) <= chopval ? m(v, w) : Inf` for points `v` and `w`.\n
`lowerbound` should be easier to evaluate than `m` and must satisfy `lowerbound(v, w) ≤ m(v, w)` for all `v`, `w`.
"""
type ChoppedQuasiMetric{M<:QuasiMetric,B<:Metric,T<:AbstractFloat} <: PreMetric
    m::M
    lowerbound::B
    chopval::T
end
"Type union encompassing `ChoppedMetric` and `ChoppedQuasiMetric`."
typealias ChoppedPreMetric{M,B,T} Union{ChoppedMetric{M,B,T}, ChoppedQuasiMetric{M,B,T}}
function evaluate(clbm::ChoppedPreMetric, v::AbstractVector, w::AbstractVector)    # TODO: maybe remove `AbstractState`?
    lb = evaluate(clbm.lowerbound, v, w)
    lb > clbm.chopval && return oftype(clbm.chopval, Inf)   # Inf is a sentinel value more numerically robust than using
    d = evaluate(clbm.m, v, w)                              # chopval; technically breaks the metric inequality
    d <= clbm.chopval ? d : oftype(clbm.chopval, Inf)
end
for M in (:ChoppedMetric, :ChoppedQuasiMetric)
    @eval changeprecision{T<:AbstractFloat}(::Type{T}, x::$M) = $M(changeprecision(T, x.m),
                                                                   changeprecision(T, x.lowerbound),
                                                                   T(x.chopval))
end
typealias SymmetricDistance Union{Metric, ChoppedMetric}
typealias AsymmetricDistance Union{QuasiMetric, ChoppedQuasiMetric}

### Steering
export ControlInfo, NullControl, StepControl, DurationAndTargetControl
export ControlSequence, ZeroOrderHoldControl, TimestampedTrajectoryControl
export duration, splitcontrol

abstract ControlInfo
immutable NullControl <: ControlInfo end
immutable StepControl{T<:AbstractFloat,N} <: ControlInfo
    t::T
    u::SVector{N,T}
end
StepControl{T}(t::T, u::AbstractVector{T}) = StepControl(t, SVector(u))
immutable DurationAndTargetControl{T<:AbstractFloat,N} <: ControlInfo
    t::T
    x::SVector{N,T}
end
DurationAndTargetControl{T}(t::T, x::AbstractVector{T}) = DurationAndTargetControl(t, SVector(x))
typealias ControlSequence{C<:ControlInfo} Vector{C}
typealias ZeroOrderHoldControl{C<:StepControl} Vector{C}
typealias TimestampedTrajectoryControl{C<:DurationAndTargetControl} Vector{C}

duration(x::ControlInfo) = x.t
duration(x::ControlSequence) = sum(duration(s) for s in x)
function splitcontrol{T}(x::StepControl{T}, t::Real)
    if t < 0
        StepControl(T(0), x.u), x
    elseif t < x.t
        StepControl(T(t), x.u), StepControl(x.t - T(t), x.u)
    else
        x, StepControl(T(0), x.u)
    end
end
function splitcontrol{T,R<:Real}(x::StepControl{T}, t::AbstractVector{R})
    issorted(t) || error("Times should be sorted as input to splitcontrol.")
    d = duration(x)
    push!([StepControl(T(dt), x.u) for dt in diff(clamp(t, T(0), d))], StepControl(d-T(t[end]), x.u))
end
function splitcontrol{T<:AbstractFloat}(x0::ZeroOrderHoldControl, t::Union{T, AbstractVector{T}})
    issorted(t) || error("Times should be sorted as input to splitcontrol.")
    x = copy(x0)    # TODO: not the most memory efficient, but simple for now
    t = clamp(t, T(0), duration(x))
    splitu = [similar(x, 0) for i in 1:length(t)+1]
    d = T(0)
    uidx = 1
    for (i, ti) in enumerate(t)
        while ti > d + duration(x[uidx]) && uidx < length(x0)    # second clause necessary because of numerical error
            push!(splitu[i], x[uidx])
            d += duration(x[uidx])
            uidx += 1
        end
        u1, u2 = splitcontrol(x[uidx], ti - d)
        push!(splitu[i], u1)
        x[uidx] = u2
        d = ti
    end
    splitu[end] = x[uidx:end]
    splitu
end


text:
using DataFrames,CSV
#Zigui writes the code to remove missing value and duplicate rows.
function remove_missing_duplicate(file; delim = ',', header = 0,miss = missing)
    """
    Remove missing value and duplicate
    """
    dat = CSV.read(file, delim = delim, header = header) # no header
    dat = unique!(dat) #remove duplicate
    nCol = ncol(dat)
    dat = dropmissing(dat) #Drop the rows containing missing
    if ismissing(miss)
        return dat
    else
        for i in 1:nCol
            dat = dat[dat[:,i].!=miss,:] #Drop the rows containing missing value we defined
        end
    end
    return dat
end

#Testing
dat = remove_missing_duplicate("pedigree_clean.txt",miss = "0")

#Jiayi modifies and optimizes her original code. She deal with the question that same individuals having different parents.
# remove missing value # may need to add an if loop to allow different missing Type
    dropmissing!(dat)

# remove rows that same individuals having different parents
    SameInd = Array{Int,1}(undef, 0)
    Duplicate_index = dat[nonunique(dat, 1), 1]
    unique!(Duplicate_index)
    for ID in Duplicate_index
        IDindex = findall(x -> x == ID, dat[1])
        append!(SameInd, IDindex)
    end

    deleterows!(dat, SameInd)

#Tianjing writes the code to deal with the problem that son has different parents and problem that
#individual appears in both sire and dam

function remove_dup_son(dat)
    """
    remove duplicate son
    """
    dup_ele = nonunique(dat, 1)
    dup_index=findall(x -> x==true, dup_ele)
    deleterows!(dat, dup_index)

    return dat
end

function same_sire_dam(dat)
    """
    remove element in both dam and sire
    """

    nRow = nrow(dat)
    row_record=[]
    for i in 1:nRow
        dam=dat[i,2]
        if dam in dat[3] && dam != '0' #in sire
            push!(row_record, i)
        end
    end

    deleterows!(dat, row_record)

    return dat
end

#Testing
remove_dup_son(dat)
same_sire_dam(dat)


text:
""" [Soares2004](@cite) """
Soares_q_tot(::Type{FT}) where {FT} = z -> if z <= 1350.0
        FT(5.0e-3) - FT(3.7e-4) * z / 1000
    else
        FT(5.0e-3) - FT(3.7e-4) * FT(1.35) - FT(9.4e-4) * (z - 1350) / 1000
    end
""" [Soares2004](@cite) """
Soares_θ_liq_ice(::Type{FT}) where {FT} = z -> if z <= 1350.0
        FT(300.0)
    else
        FT(300) + 2 * (z - 1350) / 1000
    end
""" [Soares2004](@cite) """
Soares_u(::Type{FT}) where {FT} = z -> FT(0.01)

""" [Soares2004](@cite) """
Soares_tke(::Type{FT}) where {FT} = z -> if z <= 1600.0
        FT(0.1) * FT(1.46) * FT(1.46) * (1 - z / 1600)
    else
        FT(0)
    end


text:
function test_method(func)
    frame = zeros(101, 101)
    frame[77, 77] = 1.0
    # central fwhm is all NaN
    @test func(frame, (51, 51), 3) === NaN
    # centered on point is infinite, because bkg stddev is 0
    @test func(frame, (77, 77), 3) ≈ Inf
end

@testset "snr" begin
    test_method(snr)
end

@testset "significance" begin
    test_method(significance)

    X = randn(rng, 101, 101) .+ 10
    # expect small sample to matter ~3fwhm out
    @test abs(significance(X, (58, 51), 3)) < abs(snr(X, (58, 51), 3))
end

@testset "noise" begin
    data = ones(101, 101)

    @test noise(data, (51, 51), 10) === NaN
end

@testset "detectionmap - $func" for func in (snr, significance, noise)
    data = ones(101, 101)
    _map = detectionmap(func, data, 10)
    # center should be fill value
    @test _map[51, 51] ≈ 0
    # there should be some values for the rest
    @test maximum(_map) > 0
    @test abs(minimum(_map)) ≥ 0
end


text:
using Netpbm, Images, ColorTypes, FixedPointNumbers, FileIO
using FactCheck

facts("IO") do
    workdir = joinpath(tempdir(), "Images")
    isdir(workdir) && rm(workdir, recursive=true)
    mkdir(workdir)

    context("Gray pgm") do
        a = rand(2,3)
        aa = convert(Array{U8}, a)
        fn = File(format"PGMBinary", joinpath(workdir, "3by2.pgm"))
        save(fn, a)
        b = load(fn)
        @fact convert(Array, b) --> aa
        save(fn, aa)
        b = load(fn)
        @fact convert(Array, b) --> aa
        aaimg = Images.grayim(aa')
        b = load(fn)
        @fact b --> aaimg
        aa = convert(Array{UFixed16}, a)
        save(fn, aa)
        b = load(fn)
        @fact convert(Array, b) --> aa
    end

    context("Color ppm") do
        fn = File(format"PPMBinary", joinpath(workdir, "3by2.ppm"))
        img = Images.colorim(rand(3,3,2))
        img24 = convert(Images.Image{RGB24}, img)
        save(fn, img24)
        b = load(fn)
        imgrgb8 = convert(Images.Image{RGB{U8}}, img)
        @fact Images.data(imgrgb8) --> Images.data(b)

        bb = load(fn)
        @fact data(bb) --> data(imgrgb8)
    end

    context("Colormap usage") do
        datafloat = reshape(linspace(0.5, 1.5, 6), 2, 3)
        dataint = round(UInt8, 254*(datafloat .- 0.5) .+ 1)  # ranges from 1 to 255
        # build our colormap
        b = RGB(0,0,1)
        w = RGB(1,1,1)
        r = RGB(1,0,0)
        cmaprgb = Array(RGB{U8}, 255)
        f = linspace(0,1,128)
        cmaprgb[1:128] = [(1-x)*b + x*w for x in f]
        cmaprgb[129:end] = [(1-x)*w + x*r for x in f[2:end]]
        img = Images.ImageCmap(dataint, cmaprgb)
        save(File(format"PPMBinary", joinpath(workdir,"cmap.ppm")), img)
        cmaprgb = Array(RGB, 255) # poorly-typed cmap, issue #336
        cmaprgb[1:128] = [(1-x)*b + x*w for x in f]
        cmaprgb[129:end] = [(1-x)*w + x*r for x in f[2:end]]
        img = Images.ImageCmap(dataint, cmaprgb)
        save(File(format"PPMBinary", joinpath(workdir, "cmap.ppm")), img)
    end

    context("Clamping (Images issue #256)") do
        A = grayim(rand(2,3))
        A[1,1] = -0.4
        fn = File(format"PGMBinary", joinpath(workdir, "2by3.pgm"))
        @fact_throws InexactError save(fn, A)
        save(fn, A, mapi=mapinfo(Clamp, A))
        B = load(fn)
        A[1,1] = 0
        @fact B --> map(Gray{U8}, A)
    end
end
FactCheck.exitstatus()


text:
@inline digitsRequired(bitsOfPrecision) = ceil(Int, bitsOfPrecision*0.3010299956639812)

# rubric digit display counts
# values for stringsmall, stringcompact, string, stringlarge, stringall
#@enum EXTENT small=1 compact=2 medium=3 large=4 all=5
const small=1;
const compact=2;
const medium=3;
const large=4;
const all=5;

macro I16(x)
    begin quote
        ($x)%Int16
    end end
end

const midDigits = [ @I16(8), @I16(15), @I16(25), @I16(50), @I16(1200) ]
const radDigits = [ @I16(3), @I16(6),  @I16(12), @I16(25), @I16(64)   ]

const nExtents  = length(midDigits);  const maxDigits = 1200;

function set_midpoint_digits_shown(idx::Int, ndigits::Int)
    (1 <= idx     <= nExtents) || throw(ErrorException("invalid EXTENT index ($idx)"))
    (1 <= ndigits <= maxDigits)  || throw(ErrorException("midpoint does not support an $(ndigits) digit count"))
    @inbounds midDigits[ idx ] = @I16(ndigits)
    return ndigits
end
function set_radius_digits_shown(idx::Int, ndigits::Int)
    (1 <= idx     <= nExtents) || throw(ErrorException("invalid EXTENT index ($idx)"))
    (1 <= ndigits <= maxDigits)  || throw(ErrorException("radius does support an ($ndigits) digit count"))
    @inbounds midDigits[ idx ] = @I16(ndigits)
    return ndigits
end

set_midpoint_digits_shown(ndigits::Int) = midDigits[Int(medium)] = ndigits
set_radius_digits_shown(ndigits::Int) = radDigits[Int(medium)] = ndigits

function get_midpoint_digits_shown(idx::Int)
    (1 <= idx <= nExtents) || throw(ErrorException("invalid EXTENT index ($idx)"))
    @inbounds m = midDigits[ idx ]
    return m
end
function get_radius_digits_shown(idx::Int)
    (1 <= idx <= nExtents) || throw(ErrorException("invalid EXTENT index ($idx)"))
    @inbounds r = radDigits[ idx ]
    return r
end

get_midpoint_digits_shown() = midDigits[ Int(medium) ]
get_radius_digits_shown()   = radDigits[ Int(medium) ]

digits_offer_bits(ndigits::Int) = convert(Int, div(abs(ndigits)*log2(10), 1))

# prepare the scene for balanced in displaybles interfacing

const nonfinite_strings = ["NaN", "+Inf", "-Inf", "±Inf"];

function string_nonfinite{T<:ArbFloat}(x::T)::String
    global nonfinite_strings
    idx = isnan(x) + isinf(x)<<2 - ispositive(x)<<1 - isnegative(x)
    return nonfinite_strings[idx]
end

#=
function string{T<:ArbFloat}(x::T)::String
    s = if isfinite(x)
            stringmedium(x, ndigits)
        else
            string_nonfinite(x)
        end
    return s
end
=#

function string{T<:ArbFloat}(x::T, ndigits::Int)::String
    rdigits = min(ndigits, get_radius_digits_shown(medium))

    s = if isfinite(x)
            isexact(x) ? string_exact(x, ndigits) : string_inexact(x, ndigits, rdigits)
        else
            string_nonfinite(x)
        end
    return s
end
@inline string{T<:ArbFloat}(x::T, ndigits::Int16) = string(x, ndigits%Int)

function string{T<:ArbFloat}(x::T, mdigits::Int, rdigits::Int)::String
    s = if isfinite(x)
            isexact(x) ? string_exact(x, mdigits) : string_inexact(x, mdigits, rdigits)
        else
            string_nonfinite(x)
        end
    return s
end
@inline string{T<:ArbFloat}(x::T, mdigits::Int16, rdigits::Int16) = string(x, mdigits%Int, rdigits%Int)

function arb_string{P}(x::ArbFloat{P}, digs::Int, mode::Int)::String
    cstr = ccall(@libarb(arb_get_str), Ptr{UInt8}, (Ptr{ArbFloat}, Int, UInt), &x, digs, mode%UInt)
    s = unsafe_string(cstr)
    return s
end

function string_exact{T<:ArbFloat}(x::T, mdigits::Int)::String
    digs = Int(mdigits)
    cstr = ccall(@libarb(arb_get_str), Ptr{UInt8}, (Ptr{ArbFloat}, Int, UInt), &x, digs, 2%UInt)
    s = unsafe_string(cstr)
    return cleanup_numstring(s, isinteger(x))
end
@inline string_exact{T<:ArbFloat}(x::T, mdigits::Int16) = string_exact(x, mdigits%Int)

function string_inexact{T<:ArbFloat}(x::T, mdigits::Int, rdigits::Int)::String
    mid = string_exact(midpoint(x), mdigits)
    rad = string_exact(radius(x), rdigits)
    return string(mid, "±", rad)
end
@inline string_inexact{T<:ArbFloat}(x::T, mdigits::Int16, rdigits::Int16) = string_inexact(x, mdigits%Int, rdigits%Int)
@inline string_inexact{T<:ArbFloat}(x::T, mdigits::Int, rdigits::Int16) = string_inexact(x, mdigits, rdigits%Int)
@inline string_inexact{T<:ArbFloat}(x::T, mdigits::Int16, rdigits::Int) = string_inexact(x, mdigits%Int, rdigits)

function cleanup_numstring(numstr::String, isaInteger::Bool)::String
    # is there an exponent
    body,expn = ('e' in numstr) ? split(numstr, 'e') : (numstr, "")
    if !isaInteger
        body = rstrip(body, '0')
    elseif '.' in body
        body = string(split(numstr, '.')[1])
    end
    if expn != ""
        body = string(body,"e",expn)
    end
    if body[end]=='.'
        body = string(body,"0")
    end
    return body
end

#=
function string_pm{T<:ArbFloat}(x::T)
    P = precision(T)
    mdigs = min(digitsRequired(P),get_midpoint_digits_shown(medium))
    rdigs = get_radius_digits_shown(medium)
    s = isexact(x) ? string_exact(x, mdigs) : string_inexact(x, mdigs, rdigs)
    return s
end
=#

function stringsmall_pm{P}(x::ArbFloat{P})::String
    if isexact(x)
        return stringsmall(x)
    end
    digs = min(digitsRequired(P), get_midpoint_digits_shown(small))
    sm = string_exact(midpoint(x), digs)

    sr = try
            string(Float32(radius(x)))
         catch
            string(round(radius(x),23,2))
         end

    return string(sm,"±", sr)
end
stringsmall{P}(x::ArbFloat{P})::String =
    string_exact(x, min(get_midpoint_digits_shown(small),digitsRequired(P)))

function stringcompact_pm{P}(x::ArbFloat{P})::String
    if isexact(x)
        return stringcompact(x)
    end
    digs = min(digitsRequired(P), get_midpoint_digits_shown(compact))
    sm = string_exact(midpoint(x), digs)

    sr = try
            string(Float32(radius(x)))
         catch
            string(round(radius(x),23,2))
         end

    return string(sm,"±", sr)
end
stringcompact{P}(x::ArbFloat{P})::String =
    string_exact(x, min(get_midpoint_digits_shown(compact),digitsRequired(P)))

function stringmedium_pm{P}(x::ArbFloat{P})::String
    if isexact(x)
        return stringmedium(x)
    end
    digs = min(digitsRequired(P), get_midpoint_digits_shown(medium))
    sm = string_exact(midpoint(x), digs)
    sr = try
            string(Float64(radius(x)))
         catch
            string(round(radius(x),58,2))
         end

    return string(sm,"±", sr)
end
stringmedium{P}(x::ArbFloat{P})::String =
    string_exact(x, min(get_midpoint_digits_shown(medium),digitsRequired(P)))

@inline string_pm{P}(x::ArbFloat{P}) = stringmedium_pm(x)
@inline string{P}(x::ArbFloat{P}) = stringmedium(x)

function stringlarge_pm{P}(x::ArbFloat{P})::String
    if isexact(x)
        return stringlarge(x)
    end
    digs = min(digitsRequired(P), get_midpoint_digits_shown(large))
    sm = string_exact(midpoint(x), digs)
    sr = try
            string(Float64(radius(x)))
         catch
            string(round(radius(x),58,2))
         end

    return string(sm,"±", sr)
end
stringlarge{P}(x::ArbFloat{P})::String =
    string_exact(x, min(get_midpoint_digits_shown(large),digitsRequired(P)))

function stringall_pm{P}(x::ArbFloat{P})::String
    if isexact(x)
        return stringall(x)
    end
    sm = string_exact(midpoint(x), digitsRequired(P))
    sr = try
            string(Float64(radius(x)))
         catch
            string(round(radius(x),58,2))
         end

    return string(sm,"±", sr)
end
stringall{P}(x::ArbFloat{P})::String =
    string_exact(x, digitsRequired(P))


text:
# Autogenerated wrapper script for Qt5Declarative_jll for x86_64-apple-darwin
export libqt5qml, libqt5qmlmodels, libqt5qmlworkerscript, libqt5quick, libqt5quickparticles, libqt5quickshapes, libqt5quicktest, libqt5quickwidgets

using Qt5Base_jll
JLLWrappers.@generate_wrapper_header("Qt5Declarative")
JLLWrappers.@declare_library_product(libqt5qml, "@rpath/QtQml")
JLLWrappers.@declare_library_product(libqt5qmlmodels, "@rpath/QtQmlModels")
JLLWrappers.@declare_library_product(libqt5qmlworkerscript, "@rpath/QtQmlWorkerScript")
JLLWrappers.@declare_library_product(libqt5quick, "@rpath/QtQuick")
JLLWrappers.@declare_library_product(libqt5quickparticles, "@rpath/QtQuickParticles")
JLLWrappers.@declare_library_product(libqt5quickshapes, "@rpath/QtQuickShapes")
JLLWrappers.@declare_library_product(libqt5quicktest, "@rpath/QtQuickTest")
JLLWrappers.@declare_library_product(libqt5quickwidgets, "@rpath/QtQuickWidgets")
function __init__()
    JLLWrappers.@generate_init_header(Qt5Base_jll)
    JLLWrappers.@init_library_product(
        libqt5qml,
        "lib/QtQml.framework/Versions/5/QtQml",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5qmlmodels,
        "lib/QtQmlModels.framework/Versions/5/QtQmlModels",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5qmlworkerscript,
        "lib/QtQmlWorkerScript.framework/Versions/5/QtQmlWorkerScript",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quick,
        "lib/QtQuick.framework/Versions/5/QtQuick",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickparticles,
        "lib/QtQuickParticles.framework/Versions/5/QtQuickParticles",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickshapes,
        "lib/QtQuickShapes.framework/Versions/5/QtQuickShapes",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quicktest,
        "lib/QtQuickTest.framework/Versions/5/QtQuickTest",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt5quickwidgets,
        "lib/QtQuickWidgets.framework/Versions/5/QtQuickWidgets",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module Tools

export kw_qreg, kw_creg, kw_gate, qasm_id, qasm_int, qasm_f64, qasm_str, cmp_ast, cmp_exp

using ..Types
using MLStyle
using RBNF: Token

const kw_qreg = Token{:reserved}("qreg")
const kw_creg = Token{:reserved}("creg")
const kw_gate = Token{:reserved}("gate")

qasm_id(s) = qasm_id(string(s))
qasm_id(s::String) = Token{:id}(s)
qasm_int(n::Int) = Token{:int}(string(n))
qasm_f64(x) = Token{:float64}(string(Float64(x)))
qasm_str(x) = Token{:str}(string("\"", x, "\""))

Base.isapprox(lhs::ASTNode, rhs::ASTNode) = cmp_ast(lhs, rhs)
Base.isapprox(lhs::Token, rhs::Token) = cmp_ast(lhs, rhs)

function cmp_ast(lhs::Token{A}, rhs::Token{B}) where {A, B}
    A === B || return false
    return lhs.str == rhs.str
end

function cmp_ast(x::Vector, y::Vector)
    length(x) == length(y) || return false
    for (x_stmt, y_stmt) in zip(x, y)
        cmp_ast(x_stmt, y_stmt) || return false
    end
    return true
end

cmp_ast(lhs::ASTNode, rhs::ASTNode) = false
function cmp_ast(lhs::T, rhs::T) where {T <: ASTNode}
    for name in fieldnames(T)
        a = getfield(lhs, name)
        b = getfield(rhs, name)
        cmp_ast(a, b) || return false
    end
    return true
end

function cmp_ast(lhs::Instruction, rhs::Instruction)
    lhs.name == rhs.name || return false
    all(map(cmp_exp, lhs.cargs, rhs.cargs)) || return false
    return cmp_ast(lhs.qargs, rhs.qargs)
end

function cmp_ast(lhs::UGate, rhs::UGate)
    cmp_exp(lhs.y, rhs.y) || return false
    cmp_exp(lhs.z1, rhs.z1) || return false
    cmp_exp(lhs.z2, rhs.z2) || return false
    return true
end

cmp_ast(lhs::Neg, rhs::Neg) = cmp_exp(lhs.val, rhs.val)

function cmp_exp(lhs::Tuple, rhs::Tuple)
    length(lhs) == length(rhs) || return false
    return all(map(cmp_ast, lhs, rhs))
end

function cmp_exp(lhs, rhs)
    @switch (lhs, rhs) begin
        @case (Neg(a::Token{:float64}), b::Token{:float64}) || (b::Token{:float64}, Neg(a::Token{:float64}))
            startswith(b.str, '-') && a.str == b.str[2:end]
        @case _
            cmp_ast(lhs, rhs)
    end
end

cmp_ast(lhs, rhs) = lhs == rhs

@deprecate issimilar(x, y) cmp_ast(x, y)

end


text:
@testset "Hill climbing" begin
    values = sin.(0:0.01:2)
    _, idx = @inferred PhasedArray.hill_climbing(values, 10)
    @test idx == CartesianIndex(158)
    values_2D = values * values'
    _, indices = @inferred PhasedArray.hill_climbing(values_2D, 10, 10)
    @test indices == CartesianIndex(158, 158)
end

@testset "Estimate DOA" begin
    ant_pos = @SMatrix([-1 1 -1 1; -1 -1 1 1; 0 0 0 0])
    manifold = IdealManifold(1, 1 / 4 * ant_pos; c0 = 1)

    @test est_doa(manifold, a -> abs(SVector( 1,  1,  1,  1)' * a) / norm(a)) == Spherical(1.0, 0.0, π / 2)
    @test est_doa(manifold, a -> abs(SVector( 1,  1, -1, -1)' * a) / norm(a)) == Spherical(1.0, π / 2, 0.0)
    @test est_doa(manifold, a -> abs(SVector(-1,  1, -1,  1)' * a) / norm(a)) == Spherical(1.0, 0.0, 0.0)

    doa = Spherical(1, 137π / 180, 40π / 180)
    @test est_doa(manifold, a -> abs(get_steer_vec(manifold, doa)' * a) / norm(a)) == doa

    @test PhasedArray.sph2cart(est_doa(manifold, a -> abs(SVector(1, 1, 1, 1)' * a) / norm(a), init_az = 0.1, init_el = π / 2 - 0.1)) ≈ SVector(0, 0, 1)
    @test est_doa(manifold, a -> abs(SVector( 1,  1, -1, -1)' * a) / norm(a), init_az = π / 2 + 0.1, init_el = 0.1) == Spherical(1.0, π / 2, 0.0)
    @test est_doa(manifold, a -> abs(SVector(-1,  1, -1,  1)' * a) / norm(a), init_az = 0.1, init_el = 0.1) == Spherical(1.0, 0.0, 0.0)

    doa = Spherical(1, 137π / 180, 40π / 180)
    @test est_doa(manifold, a -> abs(get_steer_vec(manifold, doa)' * a) / norm(a), init_az = 137π / 180 + 0.1, init_el = 40π / 180 - 0.1) == doa

    @test est_doa_by_signal_subspace(manifold, SVector( 1,  1,  1,  1)) == Spherical(1.0, 0.0, π / 2)
    @test est_doa_by_signal_subspace(manifold, SVector( 1,  1, -1, -1)) == Spherical(1.0, π / 2, 0.0)
    @test est_doa_by_signal_subspace(manifold, SVector(-1,  1, -1,  1)) == Spherical(1.0, 0.0, 0.0)

    doa = Spherical(1, 137π / 180, 40π / 180)
    @test est_doa_by_signal_subspace(manifold, get_steer_vec(manifold, doa)) == doa

    steer_vec = SVector( 1,  1,  1,  1)
    noise_subspace = I - steer_vec * steer_vec' / (steer_vec' * steer_vec)
    @test est_doa_by_noise_subspace(manifold, noise_subspace) == Spherical(1.0, 0.0, π / 2)
    @test est_doa_by_music(manifold, noise_subspace) == Spherical(1.0, 0.0, π / 2)
    steer_vec = SVector( 1,  1, -1, -1)
    noise_subspace = I - steer_vec * steer_vec' / (steer_vec' * steer_vec)
    @test est_doa_by_noise_subspace(manifold, noise_subspace) == Spherical(1.0, π / 2, 0.0)
    @test est_doa_by_music(manifold, noise_subspace) == Spherical(1.0, π / 2, 0.0)

    doa = Spherical(1, 137π / 180, 40π / 180)
    steer_vec = get_steer_vec(manifold, doa)
    noise_subspace = I - steer_vec * steer_vec' / (steer_vec' * steer_vec)
    @test est_doa_by_noise_subspace(manifold, noise_subspace) == doa
    @test est_doa_by_music(manifold, noise_subspace) == doa
end


text:
using HorseML.LossFunction
using HorseML.Clustering
using HorseML.Clustering: fit!
using DataFrames
using CSV

@testset "LossFunction" begin
    sy, st = rand(), rand()
    y, t = rand(20), rand(20)
    my1, my2, my3 = fill!(Array{Float64}(undef, 1, 1), sy), y[:, :], rand(20, 2)
    t1, t2 = rand(10, 10, 2, 3), rand(10, 10, 2, 3)

    LOSSES = [
        mse, cee, mae, huber, logcosh_loss, poisson, hinge, smooth_hinge
    ]

    for f in LOSSES
        @test typeof(f(sy, st)) <: Number
        @test typeof(f(y, t)) == Float64
        @test typeof(f(y, t, reduction="sum")) == Float64
        @test size(f(y, t, reduction="none")) == (20, )
        @test_throws ArgumentError f(y, t, reduction="NaN")
        @test f(my1, st) == f(sy, st)
        @test f(my2, st) == f(y, st)
        @test_throws DimensionMismatch f(my3, st)
        @test f(my2, t) == f(y, t)
        @test_throws DimensionMismatch f(my3, t)
        @test_nowarn f(t1, t2)
    end
end


text:
#This file is the Free-Wake code

function rigidwake(λ,filename)
  pointfile = open(filename,"w")
  write(pointfile,"  TITLE = 'RIGID WAKE OF ROTOR'\n")
  ktip = 0.88 #Default Tip Vortex ScrollUp Position
  # dt = dψ/Ω
  Nrw = Int64(2.0*π/dψ)
  rwps = Array{Vector}(Nrw,Nb,NR)
  rbp = Array{Float64}(NR,Nb)
  ζ0 = Array{Float64}(Nrw,Nb,NR)
  # Rtip = Array{Float64}(NR)
  for k in 1:NR
    # Rtip[k] = ktip*R
    for j in 1:Nb
      write(pointfile,"  VARIABLES = 'X','Y','Z'\n")
      write(pointfile,"  ZONE I = $(Nrw)\tF = POINT\n")
      rbp[k,j] = ecut+R/Nb*(j-1)
      for i in 1:Nrw
        ψ = (i-1)*dψ
        ψk = ψ+(k-1)*2*π/NR
        ζ0[i,j,k] = -ψ
        rwps[i,j,k] = [rbp[k,j]*cos(ψk)+μ*R*ψk;
                      -rbp[k,j]*sin(ψk);
                      λ*R*ψ]
        write(pointfile,
              "\t$(rwps[i,j,k][1])\t$(rwps[i,j,k][2])\t$(rwps[i,j,k][3])\n")
      end
    end
  end
  close(pointfile)
  return rwps,ζ0
end

function Gama(Cl,V∞,chord)
  gama = abs(0.5*Cl*V∞*chord)
  return gama
end

function freewake(wp0,Nwp,ψ,ζ0,Cl,v_all)
  #wake parameters
  ktip = 0.88
  dt = dψ/Ω
  rcut = 2*R
  # Nwp = Int64(2*π/dψ)
  bp = Array{Vector}(NR,Nb+1)
  for k in 1:NR
    ψk = ψ+2*π/NR*(k-1)
    for i in 1:(Nb+1)
      bp[k,i] = [(ecut+R/NR*(i-1))*cos(ψk),(ecut+R/NR*(i-1))*sin(ψk),0]
    end
  end

  Rtip = Array{Float64}(NR)
  for k in 1:NR
    Rtip[k] = ktip*R
  end

  #begin :: solve the Γ of wakepoints
  gama_bw = Array{Float64}(NR,Nb)
  for k in 1:NR
    for i in 1:Nb
      gama_bw[k,i] = Gama(Cl[k,i],v_all[k,i][2],c[k,i])
    end
  end
  #end :: the Γ of wakeponits solved

  #begin :: solve the Γ of wakepoints
  gama_wp = Array{Float64}(NR,Nb+1)
  for k in 1:NR
    gama_wp[k,1] = 0-gama_bw[k,1]
    gama_wp[k,Nb+1] = gama_bw[k,Nb]-0.0
    for i in 2:Nb
      gama_wp[k,i] = gama_bw[k,i-1]-gama_bw[k,i]
    end
  end
  #end :: the Γ of wakepoints solved

  # gama_ = sum(gama_bw)/NR # average of whole gama
  # vgamaind = Array{Vector}(Nb,Nwp)

  #begin :: initialize rigid wake
  # wp0 = rigidwake(uniforminflow(ψ)[3],"RIGIDWAKE.PLT")[1]
  #end :: rigid wake initialization

  #begin :: calculate bound vortex induced influnce on wakepoints
  vgamaind = Array{Vector}(Nwp,Nb,NR)
  for k in 1:NR
    for l in 1:Nb
      for i in 1:Nwp
        vgind = [0.0,0.0,0.0]
        for m in 1:NR
          for n in 1:Nb
            vgind += bslinevnocore(gama_bw[m,n],bp[m,n],bp[m,n+1],wp[i,l,k])
          end
        end
        vgamaind[i,l,k] = vgind
      end
    end
  end
  #end :: bound vortex induced influnce solving

  #begin :: calculate wake vortex influnce on wakepoints
  vwakeind = Array{Vector}(Nwp,Nb,NR)
  for k in 1:NR
    for l in 1:Nb
      for i in 1:Nwp
        vwvind = [0.0,0.0,0.0]
        for m in 1:NR
          for n in 1:Nb
            for r in 1:(Nwp-1)
              t_wp = abs(ψ-ζ0[r,n,m])/Ω
              vwvind += bslinev(gama_wp[m,n],
                                wp[r,n,m],wp[r+1,n,m],wp[i,l,k],t_wp)
              # vwvind += bslinevnocore(gama_wp[m,n],
              #                   wp[r,n,m],wp[r+1,n,m],wp[i,l,k])
            end
          end
        end
        vwakeind[i,l,k] = vwvind
      end
    end
  end
  #end :: wake vortex influnce on wakepoints calculated

  #begin :: wakepoints move with its local velocity
  vwpsind = Array{Vector}(Nwp,Nb,NR)
  vwps = Array{Vector}(Nwp,Nb,NR)
  dwps = Array{Vector}(Nwp,Nb,NR)
  for k in 1:NR
    for i in 1:Nb
      for j in 1:Nwp
        vwpsind[j,i,k] = vwakeind[j,i,k]+vgamaind[j,i,k]
        vwps[j,i,k] = vwpsind[j,i,k]+v_for
        dwps[j,i,k] = vwps[j,i,k]*dt
        wp[j,i,k] = wp[j,i,k]+dwps[j,i,k]
        ζ0[j,i,k] = ζ0[j,i,k]+dψ
      end
    end
  end
  #end :: wakepoints moved

  return wp,gama_wp
end

function wpgen(ψ,ζ0,wp,Nwp)  #new wakepoints generate
  l = 1
  if norm(ψ)<=1e-3
    return wp,ζ0
  else
    wpsadd = Array{Vector}(l,Nb,NR)
    ζadd = Array{Float64}(l,Nb,NR)
    rbp = Array{Float64}(NR,Nb)
    for k in 1:NR
      ψk = ψ+2*pi/NR*(k-1)
      for i in 1:Nb
        for l in 1
          rbp[k,i] = ecut+R/Nb*(i-1)
          wpsadd[l,i,k] = [rbp[k,i]*cos(ψk),
                            rbp[k,i]*sin(ψk),
                            0.0]
          ζadd[l,i,k] = ψ
        end
      end
    end
    newwp = maprepend(wp,wpsadd,Nwp,Nb,NR)
    newζ0 = maprepend(ζ0,ζadd,Nwp,Nb,NR)
    return newwp,newζ0
  end
end

function vbe(ψ,ζ0,gama_wp,wp) #calculate velocity of blade elemnts(both sys and ro coordinations)
  bp = Array{Vector}(NR,Nb+1)
  for k in 1:NR
    ψk = ψ+2*π/NR*(k-1)
    for i in 1:(Nb+1)
      bp[k,i] = [(ecut+R/NR*(i-1))*cos(ψk),(ecut+R/NR*(i-1))*sin(ψk),0]
    end
  end
  vb_ind = Array{Vector}(NR,Nb)
  for k in 1:NR
    for l in 1:Nb
      vbind = [0.0,0.0,0.0]
      for m in 1:NR
        for n in 1:Nb
          for r in 1:(Nwp-1)
            t = abs(ψ-ζ0[r,n,m])/Ω
            vbind += bslinev(gama_wp[m,n],wp[r,n,m],wp[r+1,n,m],bp[k,l],t)
          end
        end
      end
      vb_ind[k,l] = vbind
    end
  end

  rbp = Array{Float64}(NR,Nb) #radius of balde element points
  v_ro = Array{Vector}(NR,Nb)
  vb_sys = Array{Vector}(NR,Nb)
  vb_ro = Array{Vector}(NR,Nb)
  for k in 1:NR
    ψk = ψ+2*π/NR*(k-1)
    for i in 1:Nb
      rbp[k,i] = ecut+R/Nb*(i-1)
      v_ro[k,i] = [Ω*rbp[k,i]*sin(ψk),
                    -Ω*rbp[k,i]*cos(ψk),
                    0.0]
      vb_sys[k,i] = v_ro[k,i]+v_for+vb_ind[k,i]
      vb_ro[k,i] = systoro(vb_sys[k,i],ψ)
    end
  end
  return vb_sys,vb_ro
end

max_stars_repo_name: dingzhw/RotorTrim

max_issues_repo_name: helicopterliu/RotorTrim

max_forks_repo_name: helicopterliu/RotorTrim


text:
module BBVis

using JuMP, MathProgBase, CPLEX, Gurobi, Gadfly

export NodeData, BBData, record_node_data,
       plot_mip_gap, plot_mip_gaps, plot_mip_gaps_compat, plot_times_compat, plot_progress

const sampling_rate = 10

immutable NodeData
    node::Int
    time::Float64
    incumbent_value::Float64
    bestbound::Float64
end

type BBData
    nodes::Vector{NodeData}
    name::String
end

function node_data_callback(cb, model::Model)
    node_data::Vector{NodeData} = model.ext[:bbvis].nodes
    nodes_explored  = MathProgBase.cbgetexplorednodes(cb)
    (nodes_explored ≤ 100) || (nodes_explored % sampling_rate == 0) || return
    start_time      = isempty(node_data) ? 0.0 : node_data[1].time
    current_time    = time()
    has_incumbent   = CPLEX.cbgetfeasibility(cb)
    incumbent_value = MathProgBase.cbgetobj(cb)
    bestbound       = MathProgBase.cbgetbestbound(cb)
    push!(node_data, NodeData(nodes_explored,
                              current_time - start_time,
                              has_incumbent == 1 ? incumbent_value : NaN,
                              bestbound))
    nothing
end

function record_node_data(model::JuMP.Model, name::String)
    model.ext[:bbvis] = BBData(NodeData[], name)
    addInfoCallback(model, cb -> node_data_callback(cb,model))
end

function plot_mip_gap(model::JuMP.Model)
    nodes = model.ext[:bbvis].nodes
    Gadfly.plot(layer(x=[n.node for n in nodes], y=[n.incumbent_value for n in nodes], Geom.line),
                layer(x=[n.node for n in nodes], y=[n.bestbound       for n in nodes], Geom.line),
                Guide.xlabel("Node number"), Guide.ylabel(""), Stat.step)
end

function plot_mip_gaps(models::JuMP.Model...)
    layers = Any[]
    for (it,model) in enumerate(models)
        _nodes = model.ext[:bbvis].nodes
        nodes = filter(x -> !isnan(x.incumbent_value), _nodes)
        name  = fill(model.ext[:bbvis].name, length(nodes))
        append!(layers, [layer(x=[n.node for n in nodes], y=[n.incumbent_value for n in nodes], color=[name], Geom.line, Stat.step),
                         layer(x=[n.node for n in nodes], y=[n.bestbound       for n in nodes], color=[name], Geom.line, Stat.step)])
    end
    return plot(layers..., Guide.xlabel("Node number"), Guide.ylabel(""), Stat.step)
end

function plot_times(models::JuMP.Model...)
    layers = Any[]
    for (it,model) in enumerate(models)
        _nodes = model.ext[:bbvis].nodes
        nodes = filter(x -> !isnan(x.incumbent_value), _nodes)
        name  = fill(model.ext[:bbvis].name, length(nodes))
        append!(layers, [layer(x=[n.time for n in  nodes], y=[n.incumbent_value for n in  nodes], color=[ name], Geom.line, Stat.step),
                         layer(x=[n.time for n in  nodes], y=[n.bestbound       for n in  nodes], color=[ name], Geom.line, Stat.step)])
    end
    return plot(layers..., Guide.xlabel("Time"), Guide.ylabel(""), Stat.step)
end

function plot_progress(models::JuMP.Model...)
    node_layers = Any[]
    for (it,model) in enumerate(models)
        nodes = model.ext[:bbvis].nodes
        name  = fill(model.ext[:bbvis].name, length(nodes))
        append!(node_layers, [layer(x=[n.node for n in nodes], y=[n.incumbent_value for n in nodes], color=[name], Stat.step, Geom.line),
                              layer(x=[n.node for n in nodes], y=[n.bestbound       for n in nodes], color=[name], Stat.step, Geom.line)])
    end

    time_layers = Any[]
    for (it,model) in enumerate(models)
        nodes = model.ext[:bbvis].nodes
        name  = fill(model.ext[:bbvis].name, length(nodes))
        append!(time_layers, [layer(x=[n.time for n in nodes], y=[n.incumbent_value for n in nodes], color=[name], Stat.step, Geom.line),
                              layer(x=[n.time for n in nodes], y=[n.bestbound       for n in nodes], color=[name], Stat.step, Geom.line)])
    end
    return vstack(plot(time_layers..., Guide.xlabel("Time (s)"),    Guide.ylabel(""), Stat.step),
                  plot(node_layers..., Guide.xlabel("Node number"), Guide.ylabel(""), Stat.step))
end

function plot_mip_gaps_compat(models::JuMP.Model...)
    max_nodes = 0
    for model in models
        max_nodes = max(max_nodes, model.ext[:bbvis].nodes[end].node)
    end
    plots = Any[]
    for (it,model) in enumerate(models)
        nodes = model.ext[:bbvis].nodes
        name  = model.ext[:bbvis].name
        # p = plot(layer(x=vcat([n.node for n in nodes],max_nodes), y=vcat([n.incumbent_value for n in nodes], nodes[end].incumbent_value), Stat.step, Geom.line),
        #          layer(x=vcat([n.node for n in nodes],max_nodes), y=vcat([n.bestbound       for n in nodes], nodes[end].bestbound),       Stat.step, Geom.line),
        #          Guide.xlabel("Node number"), Guide.ylabel(name), Stat.step)
        p = plot(layer(x=[n.node for n in nodes], y=[n.incumbent_value for n in nodes], Stat.step, Geom.line),
                 layer(x=[n.node for n in nodes], y=[n.bestbound       for n in nodes], Stat.step, Geom.line),
                 Guide.xlabel("Node number"), Guide.ylabel(name), Stat.step, Scale.x_continuous(maxvalue=max_nodes))
        push!(plots, p)
    end
    vstack(plots...)
end

function plot_times_compat(models::JuMP.Model...)
    max_nodes = 0
    for model in models
        max_nodes = max(max_nodes, model.ext[:bbvis].nodes[end].node)
    end
    plots = Any[]
    for (it,model) in enumerate(models)
        nodes = model.ext[:bbvis].nodes
        name  = model.ext[:bbvis].name
        p = plot(layer(x=[n.time for n in nodes], y=[n.incumbent_value for n in nodes], Stat.step, Geom.line),
                 layer(x=[n.time for n in nodes], y=[n.bestbound       for n in nodes], Stat.step, Geom.line),
                 Guide.xlabel("Node number"), Guide.ylabel(name), Stat.step, Scale.x_continuous(maxvalue=max_nodes))
        push!(plots, p)
    end
    vstack(plots...)
end

end


text:
# Autogenerated wrapper script for helpme_jll for aarch64-linux-gnu-cxx11
export libhelpme

using FFTW_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("helpme")
JLLWrappers.@declare_library_product(libhelpme, "libhelpme.so")
function __init__()
    JLLWrappers.@generate_init_header(FFTW_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libhelpme,
        "lib/libhelpme.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#include("../src/BayesianOptimization.jl")
using BayesianOptimization

function foo(a, b)
    return sin(a) + b
end

xmax, ymax, opt = optimize(foo, 1.0:5.0, -2.0:4.0; iter=5, restarts=2)

#plot(gp)


text:
"""
    PerlinNoise <: NeutralLandscapeMaker

    PerlinNoise(; kw...)

Create a Perlin noise neutral landscape model with values ranging 0-1.

# Keywords

- `periods::Tuple{Int,Int}=(1,1)`: the number of periods of Perlin noise across row and
    column dimensions for the first octave.
- `octaves::Int=1`: the number of octaves that will form the Perlin noise.
- `lacunarity::Int=2` : the rate at which the frequency of periods increases for each
    octive.
- `persistance::Real=0.5` : the rate at which the amplitude of periods decreases for each
    octive.
- `valley::Symbol=`:u`: the kind of valley bottom that will be mimicked: `:u` produces
    u-shaped valleys, `:v` produces v-shaped valleys, and `:-` produces flat bottomed
    valleys.

Note: This is a memory-intensive algorithm with some settings. Be careful using larger
prime numbers for `period` when also using a large array size, high lacuarity and/or many
octaves. Memory use scales with the lowest common multiple of `periods`.
"""
Base.@kwdef struct PerlinNoise{PT<:Real} <: NeutralLandscapeMaker
    periods::Tuple{Int,Int} = (1, 1)
    octaves::Int = 1
    lacunarity::Int = 2
    persistance::PT = 0.5
    valley::Symbol = :u
end

function _landscape!(A, alg::PerlinNoise)
    # nRow must equal nCol so determine the dimension of the smallest square
    dim = max(size(A)...)
    # Check the dim is a multiple of each octives maximum number of periods and
    # expand dim if needed
    rperiodsmax = alg.periods[1] * alg.lacunarity^(alg.octaves - 1)
    cperiodsmax = alg.periods[2] * alg.lacunarity^(alg.octaves - 1)
    periodsmultiple = lcm(rperiodsmax, cperiodsmax) # lowest common multiple
    if dim % periodsmultiple != 0
        dim = ceil(Int, dim / periodsmultiple) * periodsmultiple
    end

    # Generate the Perlin noise
    noise = zeros(eltype(A), (dim, dim))
    meshbuf = Array{eltype(A),3}(undef, dim, dim, 2)
    nbufs = ntuple(_->Array{eltype(A),2}(undef, dim, dim), 4)
    for octave in 0:alg.octaves-1
        octave_noise!(noise, meshbuf, nbufs, alg, octave, dim, dim)
    end

    # Randomly extract the desired array size
    noiseview = _view_from_square(noise, size(A)...)

    # Rescale the Perlin noise to mimic different kinds of valley bottoms
    return if alg.valley == :u
        A .= noiseview
    elseif alg.valley == :v
        A .= abs.(noiseview)
    elseif alg.valley == :-
        A .= noiseview.^2
    else
        error("$(alg.valley) not recognised for `valley` use `:u`, `:v` or `:-`")
    end
end

function octave_noise!(
    noise, mesh, (n11, n21, n12, n22), alg::PerlinNoise, octave, nrow, ncol
)
    f(t) = @fastmath 6 * t ^ 5 - 15 * t ^ 4 + 10 * t ^ 3 # Wut

    # Mesh
    rp, cp = alg.periods .* alg.lacunarity^(octave)
    delta = (rp / nrow, cp / ncol)
    ranges = range(0, rp-delta[1], length=nrow), range(0, cp-delta[2], length=ncol)
    _mesh!(mesh, ranges...)
    @fastmath mesh .%= 1

    # Gradients
    # This allocates, but the gradients size changes with octave so it needs to
    # be a very smart in-place `repeat!` or some kind of generator, and the improvement
    # may not be than large (~20%).
    angles = 2pi .* rand(rp + 1, cp + 1)
    @fastmath gradients = cat(cos.(angles), sin.(angles); dims=3)
    d = (nrow ÷ rp, ncol ÷ cp)
    grad = repeat(gradients, inner=[d[1], d[2], 1])

    g111 = @view grad[1:ncol,         1:ncol,         1]
    g211 = @view grad[end-nrow+1:end, 1:ncol,         1]
    g121 = @view grad[1:ncol,         end-ncol+1:end, 1]
    g221 = @view grad[end-nrow+1:end, end-ncol+1:end, 1]
    g112 = @view grad[1:ncol,         1:ncol,         2]
    g212 = @view grad[end-nrow+1:end, 1:ncol,         2]
    g122 = @view grad[1:ncol,         end-ncol+1:end, 2]
    g222 = @view grad[end-nrow+1:end, end-ncol+1:end, 2]

    # Ramps
    m1 = @view mesh[:, :, 1]
    m2 = @view mesh[:, :, 2]
    n11 .= ((m1     .+ m2     ) .* g111 .+ (m1     .+ m2     ) .* g112)
    n21 .= ((m1 .-1 .+ m2     ) .* g211 .+ (m1 .-1 .+ m2     ) .* g212)
    n12 .= ((m1     .+ m2 .- 1) .* g121 .+ (m1     .+ m2 .- 1) .* g122)
    n22 .= ((m1 .-1 .+ m2 .- 1) .* g221 .+ (m1 .-1 .+ m2 .- 1) .* g222)

    # Interpolation
    mesh .= f.(mesh)
    t1 = @view mesh[:, :, 1]
    t2 = @view mesh[:, :, 2]
    noise .+= sqrt(2) .*  (alg.persistance ^ octave) .*
        ((1 .- t2) .* (n11 .* (1 .- t1) .+ t1 .* n21) .+
               t2  .* (n12 .* (1 .- t1) .+ t1 .* n22))
    return noise
end

function _mesh!(A, x, y)
    for (i, ival) in enumerate(x), j in 1:length(y)
        A[i, j, 1] = ival
    end
    for i in 1:length(x), (j, jval) in enumerate(y)
        A[i, j, 2] = jval
    end
    return A
end

function _view_from_square(source, nrow, ncol)
    # Extract a portion of the array to match the dimensions
    dim = size(source, 1)
    startrow = rand(1:(dim - nrow + 1))
    startcol = rand(1:(dim - ncol + 1))
    return @view source[startrow:startrow + nrow - 1, startcol:startcol + ncol - 1]
end


text:
module FileCacheAdapter

using Serialization, Logging
using Nullables
using Genie

"""
    to_cache(key::Union{String,Symbol}, content::Any; dir::String = "") :: Nothing

Persists `content` onto the file system under the `key` key.
"""
function to_cache(key::Union{String,Symbol}, content::Any; dir::String = "") :: Nothing
  open(cache_path(string(key), dir = dir), "w") do io
    serialize(io, content)
  end

  nothing
end

"""
    from_cache(key::Union{String,Symbol}, expiration::Int; dir::String = "") :: Nullable

Retrieves from cache the object stored under the `key` key if the `expiration` delta (in seconds) is in the future.
"""
function from_cache(key::Union{String,Symbol}, expiration::Int; dir::String = "") :: Nullable
  file_path = cache_path(string(key), dir = dir)

  ( ! isfile(file_path) || stat(file_path).ctime + expiration < time() ) && return Nullable()

  Genie.config.log_cache && @info("Found file system cache for $key at $file_path")

  output = open(file_path) do io
    deserialize(io)
  end

  Nullable(output)
end

"""
    purge(key::Union{String,Symbol}) :: Nothing

Removes the cache data stored under the `key` key.
"""
function purge(key::Union{String,Symbol}; dir::String = "") :: Nothing
  rm(cache_path(string(key), dir = dir))

  nothing
end

"""
    purge_all(; dir::String = "") :: Nothing

Removes all cached data.
"""
function purge_all(; dir::String = "") :: Nothing
  rm(cache_path("", dir = dir), recursive = true)
  mkpath(cache_path("", dir = dir))
  open(joinpath(cache_path("", dir = dir), ".gitkeep"), "w") do f
    write(f, "")
  end

  nothing
end

"""
    cache_path(key::Union{String,Symbol}; dir::String = "") :: String

Computes the path to a cache `key` based on current cache settings.
"""
function cache_path(key::Union{String,Symbol}; dir::String = "") :: String
  path = joinpath(Genie.CACHE_PATH, dir)
  ! isdir(path) && mkpath(path)

  joinpath(path, string(key))
end

end


text:
function update_automerge_comment!(data::GitHubAutoMergeData,
                                   body::AbstractString)::Nothing
    if data.read_only
        @info "`read_only` mode; skipping updating automerge comment."
        return nothing
    end
    api = data.api
    repo = data.registry
    pr = data.pr
    auth = data.auth
    whoami = data.whoami

    my_comments = my_retry(() -> get_all_my_pull_request_comments(api,
                                                                  repo,
                                                                  pr;
                                                                  auth = auth,
                                                                  whoami = whoami))
    num_comments = length(my_comments)
    _body = string(strip(string(body,
                                "\n",
                                "<!---\n",
                                "this_is_the_single_automerge_comment\n",
                                "--->\n")))
    if num_comments > 1
        for i = 2:num_comments
            comment_to_delete = my_comments[i]
            try
                my_retry(() -> delete_comment!(api, repo,
                                               pr,
                                               comment_to_delete;
                                               auth = auth))
            catch ex
                @error("Ignoring error: ", exception=(ex, catch_backtrace()))
            end
        end
        comment_to_update = my_comments[1]
        if strip(comment_to_update.body) != _body
            my_retry(() -> edit_comment!(api, repo,
                                         pr,
                                         comment_to_update,
                                         _body;
                                         auth = auth))
        end
    elseif num_comments == 1
        comment_to_update = my_comments[1]
        if strip(comment_to_update.body) != _body
            my_retry(() -> edit_comment!(api, repo,
                                         pr,
                                         comment_to_update,
                                         _body;
                                         auth = auth))
        end
    else
        always_assert(num_comments < 1)
        my_retry(() -> post_comment!(api, repo,
                                     pr,
                                     _body;
                                     auth = auth))
    end
    return nothing
end


text:
module Calibration

using Distributions

    function monteCarlo(paramSpace;samples=100000)
        params = Dict{Symbol,Array}()

        keyList = collect(keys(paramSpace))

        for k in keyList
            p = paramSpace[k]
            minV = p[:lower]
            maxV = p[:upper]
            params[k] = collect(rand(Uniform(minV,maxV), samples))
        end

        return params

    end

    function latinHypercube(func,forcing,obs,param_bounds;samples=10000)

    end

# end module
end

max_stars_repo_name: jsta/hymod

max_issues_repo_name: jsta/hymod

max_forks_repo_name: jsta/hymod


text:
export hamiltonian_cycle

#check if it is safe to add vertex v to the path
function isSafe(v::T, G::SimpleDigraph{T}, path::Deque{T}) where {T}
    prev = last(path)

    #check if the added vertex is an out_neighbor of the previous vertex
    Nv = out_neighbors(G, prev)
    if length(Nv) == 0
        return false
    end
    if (!in(v, Nv))
        return false
    end

    #check if the vertex already exists in the path
    if (in(v, path))
        return false
    end

    return true
end

function hamCycle(G::SimpleDigraph{T}, path::Deque{T}) where {T}
    #if all vertices are included in the cycle
    if length(path) == NV(G)
        #check if last vertex is connected to first vertex in path
        Nv = out_neighbors(G, last(path))
        if in(first(path), Nv)
            return true
        else
            return false
        end
    end

    #try different vertices as the next candidate in the Hamiltonian cycle
    vlist = collect(G.V)
    for v in vlist
        if (isSafe(v, G, path))
            push!(path, v)
            #cycle
            if (hamCycle(G, path) == true)
                return true
            end
            pop!(path)
        end
    end
    return false
end

#check if a directed graph contains a hamiltonian cycle
function directed_ham_cycle(G::SimpleDigraph{T}) where {T}
    result = Deque{T}()
    vlist = collect(G.V)

    #check if there are any isolated vertices
    simpleG = simplify(G)
    if (!is_connected(simpleG))
        return result
    end

    for v in vlist
        push!(result, v)
        if (hamCycle(G, result))
            return result
        end
        result = Deque{T}()
    end

    return result
end

#export the result as an array
function hamiltonian_cycle(G::SimpleDigraph{T}) where {T}
    return collect(directed_ham_cycle(G))
end


text:
"""
Writes the 3 plots needed for the project.
"""
function scenario_plots(processed_sojurn_times, mean_system_job_totals, proportion_orbiting_jobs, scenario_number::Int64, lambda_range::Array{Float64}, mode::Int)
    println("Plotting data...")

    if mode == 1
        PyPlot.scatter(lambda_range, mean_system_job_totals)
        PyPlot.xticks(lambda_range)
        PyPlot.title("Mean Number of Items in the Total System for Scenario $scenario_number\n")
        PyPlot.xlabel("λ")
        PyPlot.ylabel("Mean Items in System")
        PyPlot.savefig("mean_items_scenario_$scenario_number.png")
        PyPlot.close()

        PyPlot.plot(lambda_range, proportion_orbiting_jobs)
        PyPlot.ylim(bottom = 0.0)
        PyPlot.xticks(lambda_range)
        PyPlot.title("Proportion of Jobs in Orbit at T=10e^7 for Scenario $scenario_number\n")
        PyPlot.xlabel("λ")
        PyPlot.ylabel("Jobs in Orbit / Total Jobs in System")
        PyPlot.savefig("proportion_orbiting_jobs_scenario_$scenario_number.png")
        PyPlot.close()

    elseif mode == 2
        for i in processed_sojurn_times
            PyPlot.plot(sort(i), (1:length(i))./length(i))
        end
        PyPlot.legend(["λ = $i" for i in lambda_range])
        PyPlot.title("Empirical Distribution of Sojurn Times For Scenario $scenario_number\n")
        PyPlot.xlabel("Recorded Sojurn Times (Seconds)")
        PyPlot.ylabel("Fraction of Data")
        PyPlot.savefig("histogram_scenario_$scenario_number.png")
        PyPlot.close()
        # PyPlot.hist(processed_sojurn_times)
        # PyPlot.yscale("log")
        # PyPlot.legend(["λ = $i" for i in lambda_range])
        # PyPlot.title("Empirical Distribution of Sojurn Times For Scenario $scenario_number\n")
        # PyPlot.xlabel("Recorded Sojurn Times (Seconds)")
        # PyPlot.ylabel("Frequency")
        # PyPlot.savefig("histogram_scenario_$scenario_number.png")
        # PyPlot.close()
    end
end


text:
using Documenter, IGA

makedocs(
    sitename = "IGA.jl",
    format = :html,
    pages = ["Home" => "index.md",
             "API Documentation" => "api.md"]
)

max_stars_repo_name: michakraus/iga-demos

max_issues_repo_name: michakraus/iga-demos

max_forks_repo_name: michakraus/iga-demos


text:
using Documenter, LimitOrderBook

makedocs(modules = [LimitOrderBook], sitename = "LimitOrderBook.jl")

deploydocs(
    repo = "github.com/JuliaLang/LimitOrderBook.jl.git",
)


text:
# paths to npm.cmd and nodejs.exe
const NPM_PATH = npm_cmd()
const NODE_DIR = dirname(nodejs_cmd().exec[1])

# default path and model file name
const MODEL_DIR = "_julia"
const MODEL_NAME = "model.jl"

"""
    heta_update(version::String)

To install or update heta-compiler from NPM.

Arguments:

- `version` : `heta compiler` version. If the value is not provided, `heta_update` installs the latest version of `heta compiler`
"""
heta_update() = run(`$NPM_PATH i -g heta-compiler --prefix $NODE_DIR`)
# install particular version
heta_update(version::String) = run(`$NPM_PATH i -g heta-compiler@$version --prefix $NODE_DIR`)
# install from GitHub's repository

"""
    heta_update_dev(branch::String = "master")

Installs heta-compiler from GitHub's repository <https://github.com/hetalang/heta-compiler>.

Arguments:

- `branch` : branch to install, default is "master".
"""
heta_update_dev(branch::String = "master") = run(`$NPM_PATH i -g https://github.com/hetalang/heta-compiler.git\#$branch --prefix $NODE_DIR`)

"""
    heta_build(
      target_dir::AbstractString;
      declaration::String = "platform",
      skip_export::Bool = false,
      units_check::Bool = false,
      log_mode::String = "error",
      debug::Bool = false,
      julia_only::Bool = false,
      dist_dir::String = "dist",
      meta_dir::String = "meta",
      source::String = "index.heta",
      type::String = "heta"
    )

Builds the model from Heta-based reactions

See `heta comiler` docs for details:
https://hetalang.github.io/#/heta-compiler/cli-references?id=running-build-with-cli-options

Arguments:

- `target_dir` : path to a Heta platform directory
- `declaration` : path to declaration file. Default is `"platform"`
- `skip_export` : if set to `true` no files will be created. Default is `false`
- `units_check` : if set to `true` units will be checked for the consistancy
- `log_mode` : log mode. Default is `"error"`
- `debug` : turn on debug mode. Default is `false`
- `julia_only` : export only julia-based model. Default is `false`
- `dist_dir` : directory path, where to write distributives to. Default is `"dist"`
- `meta_dir` : meta directory path. Default is `"meta"`
- `source` : path to the main heta module. Default is `"index.heta"`
- `type` : type of the source file. Default is `"heta"`
"""
function heta_build(
  target_dir::AbstractString;
  declaration::String = "platform",
  skip_export::Bool = false,
  units_check::Bool = false,
  log_mode::String = "error",
  debug::Bool = false,
  julia_only::Bool = false,
  dist_dir::String = "dist",
  meta_dir::String = "meta",
  source::String = "index.heta",
  type::String = "heta"
)
  # check if heta is installed
  heta_build_path = Sys.iswindows() ? "$NODE_DIR/node_modules/heta-compiler" : "$NODE_DIR/lib/node_modules/heta-compiler"

  !isdir(heta_build_path) && throw("Heta compiler is not installed. Run `heta_update()` to install it.")

  # convert to absolute path
  _target_dir = abspath(target_dir)

  # check if the dir contains src, index.heta, platform.json
  #isdir("$_target_dir/src") && throw("src directory not found in $_target_dir")
  #!isfile("$_target_dir/index.heta") && throw("index.heta file not found in $_target_dir")
  #!isfile("$_target_dir/platform.json") && throw("platform.json file not found in $_target_dir")

  # cmd options supported by heta-compiler
  options_array = String[]

  declaration != "platform" && push!(options_array, "--declaration", declaration)
  skip_export && push!(options_array, "--skip-export")
  units_check && push!(options_array, "--units-check")
  log_mode != "error" && push!(options_array, "--log-mode", log_mode)
  debug && push!(options_array, "--debug")
  julia_only && push!(options_array, "--julia-only")
  dist_dir != "dist" && push!(options_array, "--dist-dir", dist_dir)
  meta_dir != "meta" && push!(options_array, "--meta-dir", meta_dir)
  source != "index.heta" && push!(options_array, "--source", source)
  type != "heta" && push!(options_array, "--type", type)

  # build the dist
  #=
  if  Sys.iswindows()
    heta_cmd = "heta.cmd"
  else
    heta_cmd = "heta" # not tested on unix
  end
  =#

  run_build = run(ignorestatus(`$NODE_DIR/node $heta_build_path/bin/heta-build.js $options_array $_target_dir`))
  return run_build.exitcode
end

"""
    load_platform(
      target_dir::AbstractString;
      rm_out::Bool = true,
      julia_only::Bool = true,
      dist_dir::String = ".",
      source::String = "index.heta",
      type::String = "heta",
      kwargs...
    )

Converts heta model to Julia and outputs `Platform` type.

See `heta comiler` docs for details:
https://hetalang.github.io/#/heta-compiler/cli-references?id=running-build-with-cli-options

Arguments:

- `target_dir` : path to a Heta platform directory
- `rm_out` : should the file with Julia model be removed after the model is loaded. Default is `true`
- `dist_dir` : directory path, where to write distributives to. Default is `"."`
- kwargs : other arguments supported by `heta_build`

"""
function load_platform(
  target_dir::AbstractString;
  rm_out::Bool = true,
  dist_dir::String = ".",
  kwargs...
)
  # convert heta model to julia
  build_exitcode = heta_build(target_dir; julia_only = true, dist_dir = dist_dir, kwargs...)

  # check the exitcode (0 - success, 1 - failure)
  build_exitcode == 1 && throw("Compilation errors. Likely there is an error in the code of the model. See logs")

  #convert to absolute path
  _target_dir = abspath(target_dir)

  # load model to Main
  return load_jlplatform("$_target_dir/$MODEL_DIR/$MODEL_NAME"; rm_out)
end

"""
    load_jlplatform(
      model_jl::AbstractString;
      rm_out::Bool = false
    )

Loads prebuild julia model as part of `Platform`

Arguments:

- `model_jl` : path to Julia model file
- `rm_out` : should the file with Julia model be removed after the model is loaded. Default is `false`
"""
function load_jlplatform(
  model_jl::AbstractString;
  rm_out::Bool = false
)
  # load model to Main
  Base.include(Main, model_jl)

  # remove output after model load
  rm_out && rm(dirname("$model_jl"); force=true, recursive=true)

  # tmp fix to output model without task
  # (models, tasks,) = Base.invokelatest(Main.Platform)

  platform = Base.invokelatest(Platform, Main.__platform__...)

  return platform
end

# tmp solution to add model only
"""
    load_jlmodel(
      model_jl::AbstractString;
      rm_out::Bool = false
    )

Loads prebuild julia model without `Platform`

Arguments:

- `model_jl` : path to Julia model file
- `rm_out` : should the file with Julia model be removed after the model is loaded. Default is `false`
"""
function load_jlmodel(model_jl::AbstractString; rm_out::Bool = false)
  platform = load_jlplatform(model_jl; rm_out = rm_out)

  first_model = [values(platform.models)...][1]

  return first_model
end


text:
"Initialize LIKWIDs NUMA module."
function init_numa()
    ret = LibLikwid.numa_init()
    if ret == 0
        _numainfo[] = unsafe_load(LibLikwid.get_numaTopology())
        _build_jl_numa()
        numa_initialized[] = true
        return true
    end
    return false
end

"Close and finalize LIKWIDs NUMA module."
function finalize_numa()
    LibLikwid.numa_finalize()
    numa_initialized[] = false
    _numainfo[] = nothing
    return nothing
end

function _build_jl_numa()
    ni = _numainfo[]

    nnodes = ni.numberOfNodes
    _nodes = unsafe_wrap(Array, ni.nodes, nnodes)
    nodes = Vector{NumaNode}(undef, nnodes)
    for (i, n) in enumerate(_nodes)
        nodes[i] = NumaNode(
            n.id,
            n.totalMemory, # kB
            n.freeMemory, # kB
            n.numberOfProcessors,
            unsafe_wrap(Array, n.processors, n.numberOfProcessors),
            n.numberOfDistances,
            unsafe_wrap(Array, n.distances, n.numberOfDistances),
        )
    end

    numainfo[] = NumaTopology(nnodes, nodes)
    return nothing
end

"""
    get_numa_topology() -> NumaTopology
Get the NUMA topology of the machine.

Automatically initializes the topology, NUMA, and affinity modules,
i.e. calls [`LIKWID.init_topology`](@ref), [`LIKWID.init_numa`](@ref), and [`LIKWID.init_affinity`](@ref).
"""
function get_numa_topology()
    if !topo_initialized[]
        init_topology() || error("Couldn't init topology.")
    end
    if !numa_initialized[]
        init_numa() || error("Couldn't init numa.")
    end
    if !affinity_initialized[]
        init_affinity() || error("Couldn't init affinity.")
    end
    return numainfo[]
end


text:
using Documenter, CUDAdrv

makedocs(
    modules = [CUDAdrv],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "CUDAdrv.jl",
    pages = [
        "Home"    => "index.md",
        "Manual"  => [
            "man/usage.md"
        ],
        "Library" => [
            "lib/api.md",
        ]
    ],
    doctest = true
)


text:
module test_sort_bibliography

using Test
using Bibliography

# TODO: replace this sorting testbib with a full example bibliography like
#       https://github.com/plk/biblatex/tree/dev/bibtex/bib/biblatex
#       in order to have a more real life testing scenario (the results
#       could be adapted from
#       github.com/plk/biblatex/blob/dev/testfiles/19-alphabetic-prefixed.pdf
#       )

testbib = Bibliography.import_bibtex("../examples/test_sort.bib")

@testset "sort_bibliography! : Errors" begin
    @test_throws ArgumentError sort_bibliography!(testbib, :nnn)
    @test_throws MethodError sort_bibliography!(testbib, "nyt")
    @test_throws MethodError sort_bibliography!(testbib[testbib.keys[1]], :nyt)
end

result = [
    "1998BBaA1a"
    "1998BBbA2a"
    "2007JJaAbe"
    "2007JJbAbe"
    "2010JJBaAaa"
    "2010JJBbAaa"
    "2011JJBcAaa"
    "2018DJaMad"
    "2018DJaMae"
    "SaA3a2010"
]
@testset "sort_bibliography! : Ordering: :key" begin
    sort_bibliography!(testbib)
    @test testbib.keys == result
    sort_bibliography!(testbib, :key)
    @test testbib.keys == result
end

result = [
    "2007JJbAbe"
    "2007JJaAbe"
    "2018DJaMae"
    "2018DJaMad"
    "2011JJBcAaa"
    "2010JJBaAaa"
    "2010JJBbAaa"
    "SaA3a2010"
    "1998BBbA2a"
    "1998BBaA1a"
]
@testset "sort_bibliography! : Ordering: :nyt" begin
    sort_bibliography!(testbib, :nyt)
    @test testbib.keys == result
end

result = [
    "2007JJbAbe"
    "2007JJaAbe"
    "2018DJaMae"
    "2018DJaMad"
    "2011JJBcAaa"
    "2010JJBaAaa"
    "2010JJBbAaa"
    "SaA3a2010"
    "1998BBaA1a"
    "1998BBbA2a"
]
@testset "sort_bibliography! : Ordering: :nty" begin
    sort_bibliography!(testbib, :nty)
    @test testbib.keys == result
end

result = [
    "2007JJbAbe"
    "2007JJaAbe"
    "2018DJaMae"
    "2018DJaMad"
    "2011JJBcAaa"
    "2010JJBaAaa"
    "2010JJBbAaa"
    "SaA3a2010"
    "1998BBaA1a"
    "1998BBbA2a"
]
@testset "sort_bibliography! : Ordering: :y" begin
    sort_bibliography!(testbib, :nty)
    @test testbib.keys == result
end

end # module test_sort_bibliography


text:
module ShowGraphvizTests

include("examples.jl")
include("test_examples.jl")
include("test_doctest.jl")

end  # module ShowGraphvizTests


text:
###### Integrators #####

type EvolveArgs
  integ_method::Function
  soften_len
  dt
  dt_param
  time_end
  dt_output
  dt_stats
  init_out
  x_info
  dynamic_step
  exact_time
  final_snapshot
end

function evolve(nb, arg::EvolveArgs)
  dt = arg.dt
  x_info = arg.x_info
  dt_stats = arg.dt_stats
  dt_output = arg.dt_output
  integ_method = arg.integ_method

  current_time = 0
  step = 0

  t_stats = current_time + dt_stats
  t_checkp = current_time + dt_output
  t_end = current_time + arg.time_end

  if !arg.dynamic_step
    t_stats -= (1//2)*dt
    t_checkp -= (1//2)*dt
    t_end -= (1//2)*dt
  end

  init_energy!(nb)
  write_stats(nb, step, x_info)

  if arg.init_out; write_snapshot(nb); end

  while current_time < t_end
    # Handle dynamic step
    if arg.dynamic_step
      dt = arg.dt_param * collision_time_scale(nb)

      if  arg.exact_time && current_time + dt > t_checkp
        dt = nextfloat(t_checkp - current_time)
      end
    end

    # Call the integration method
    # nb.step = step
    integ_method(nb, dt)

    current_time += dt
    nb.time = current_time
    step += 1
    if current_time >= t_stats
      write_stats(nb, step, x_info)
      t_stats += dt_stats
    end
    if current_time >= t_checkp
      write_snapshot(nb)
      t_checkp += dt_output
    end
  end

  # Write final snapshot and stats
  if arg.final_snapshot
    write_snapshot(nb)
  end
  write_stats(nb, step, x_info)
end

function forward_step(nb::NBodySystem, dt)
  nb_accel(x) = accel(x, nb)
  old_accel = map(nb_accel, nb.bodies)
  for b in nb.bodies
    b.pos += b.vel*dt
  end
  for (i,b) in enumerate(nb.bodies)
    nb.bodies[i].vel += old_accel[i]*dt
  end
end

function leapfrog_step(nb::NBodySystem, dt)
  for b in nb.bodies
    b.vel += accel(b, nb) * 0.5 * dt
  end
  for b in nb.bodies
    b.pos += b.vel * dt
  end
  for b in nb.bodies
    b.vel += accel(b, nb) * 0.5 * dt
  end
end

function rk2_step(nb::NBodySystem, dt)
  old_pos = Array{Vector{Float64}}(nb.N)
  for (i, b) in enumerate(nb.bodies); old_pos[i] = b.pos; end

  half_vel = Array{Vector{Float64}}(nb.N)
  for (i, b) in enumerate(nb.bodies)
    half_vel[i] = b.vel + accel(b, nb) * 0.5 * dt
  end
  for b in nb.bodies
    b.pos += b.vel * 0.5 * dt
  end
  for b in nb.bodies
    b.vel += accel(b, nb) * dt
  end
  for (i, b) in enumerate(nb.bodies)
    b.pos = old_pos[i] + half_vel[i] * 0.5 * dt
  end
end

function rk4_step(nb::NBodySystem, dt)
  # old_pos = Array{Vector{Float64}}(nb.N)
  # for (i, b) in enumerate(nb.bodies); old_pos[i] = b.pos; end;
  old_pos = map(x -> x.pos, nb.bodies)

  nb_accel(x) = accel(x, nb)
  a0 = map(nb_accel, nb.bodies)
  for (i, b) in enumerate(nb.bodies)
    b.pos = old_pos[i] + b.vel*0.5*dt + a0[i]*0.125*dt*dt
  end

  a1 = map(nb_accel, nb.bodies)
  for (i, b) in enumerate(nb.bodies)
    b.pos = old_pos[i] + b.vel*dt + a1[i]*0.5*dt*dt
  end

  a2 = map(nb_accel, nb.bodies)
  for (i, b) in enumerate(nb.bodies)
    b.pos = old_pos[i] + b.vel*dt + (a0[i]+a1[i]*2)*(1//6)*dt*dt
  end

  for (i, b) in enumerate(nb.bodies)
    b.vel += (a0[i]+a1[i]*4+a2[i])*(1//6)*dt
  end
end

function hermite_step(nb::NBodySystem, dt)
  old_pos = map(x -> x.pos, nb.bodies)
  old_vel = map(x -> x.vel, nb.bodies)

  old_acc = map(x -> accel(x, nb), nb.bodies)
  old_jerk = map(x -> jerk(x, nb), nb.bodies)

  for (i, b) in enumerate(nb.bodies)
    b.pos += b.vel*dt + old_acc[i]*0.5*dt*dt + old_jerk[i]*(dt*dt*dt/6)
  end

  for (i, b) in enumerate(nb.bodies)
    b.vel += old_acc[i]*dt + old_jerk[i]*(dt*dt/2)
  end

  for (i, b) in enumerate(nb.bodies)
    b.vel = old_vel[i] + (old_acc[i] + accel(b, nb))*(dt/2) + (old_jerk[i] - jerk(b, nb))*(dt*dt/12)
  end

  for (i, b) in enumerate(nb.bodies)
    b.pos = old_pos[i] + (old_vel[i] + b.vel)*(dt/2) + (old_acc[i] - accel(b, nb))*(dt*dt/12)
  end
end


text:
using
  NetCDF,
  EarthConstants,
  JLD2

filename = "SO_met_100day.nc"

ncinfo(filename)

t = ncread(filename, "time") * day
t_hrs = t / hour

# Radiation
shortwave = ncread(filename, "sw")
longwave = ncread(filename, "lw")

latentheating = ncread(filename, "qlat")
sensibleheating = ncread(filename, "qsens")

stress_x = ncread(filename, "tx")
stress_y = ncread(filename, "ty")

precip = ncread(filename, "precip")

filename = "example_southern_ocean_forcing.jld2"
file = jldopen(filename, "a")

file["t"]               = t
file["shortwave"]       = shortwave
file["longwave"]        = longwave
file["stress_x"]        = stress_x
file["stress_y"]        = stress_y
file["latentheating"]   = latentheating
file["sensibleheating"] = sensibleheating
file["precip"]          = precip

close(file)

#=
# Make plots
zeroline(t) = plot(t, 0t, color="k", alpha=0.2, linewidth=1)

close("all")
fig, axs = subplots(nrows=3, sharex=true, figsize=(12, 12))

sca(axs[1])
plot(t_hrs, stress_east;  label="East")
plot(t_hrs, stress_north; label="West")
zeroline(t_hrs)

xlabel("Time (hours)", labelpad=12)
ylabel("Wind stress (N/m\$^2\$)", labelpad=12)

cornerspines(horizontal="top")
legend(fontsize=12, labelspacing=0.5)

sca(axs[2])
plot(t_hrs, heating_latent;   label="Latent heat flux")
plot(t_hrs, heating_sensible; label="Sensible heat flux")
plot(t_hrs, shortwave;        label="Short wave radiation")
plot(t_hrs, longwave;         label="Long wave radiation")
zeroline(t_hrs)

ylabel("Heat flux (W/m)", labelpad=24)
legend(fontsize=10, labelspacing=0.2)

sidespine()

sca(axs[3])
plot(t_hrs, 1e6*precip)

xlabel("Time (hours)")
ylabel("Precipitation (\$\\mu\$m/s)", labelpad=24)

cornerspines()

tightshow()
=#


text: print("This file is here to make this a Julia repo so it shows up on the Julia Hacktoberfest issue tracker")


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule UBPF_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("UBPF")
JLLWrappers.@generate_main_file("UBPF", UUID("502467ad-4a4a-57e4-9860-6b433130b33f"))
end  # module UBPF_jll


text:
function BigFractionFormat()
    return BigFractionFormat(())
end

function BigFractionFormat(arg0::NumberFormat)
    return BigFractionFormat((NumberFormat,), arg0)
end

function BigFractionFormat(arg0::NumberFormat, arg1::NumberFormat)
    return BigFractionFormat((NumberFormat, NumberFormat), arg0, arg1)
end

function format(obj::BigFractionFormat, arg0::BigFraction, arg1::StringBuffer, arg2::FieldPosition)
    return jcall(obj, "format", StringBuffer, (BigFraction, StringBuffer, FieldPosition), arg0, arg1, arg2)
end

function format(obj::BigFractionFormat, arg0::Object, arg1::StringBuffer, arg2::FieldPosition)
    return jcall(obj, "format", StringBuffer, (Object, StringBuffer, FieldPosition), arg0, arg1, arg2)
end

function format_big_fraction(::Type{BigFractionFormat}, arg0::BigFraction)
    return jcall(BigFractionFormat, "formatBigFraction", JString, (BigFraction,), arg0)
end

function get_available_locales(::Type{BigFractionFormat})
    return jcall(BigFractionFormat, "getAvailableLocales", Vector{Locale}, ())
end

function get_improper_instance(::Type{BigFractionFormat})
    return jcall(BigFractionFormat, "getImproperInstance", BigFractionFormat, ())
end

function get_improper_instance(::Type{BigFractionFormat}, arg0::Locale)
    return jcall(BigFractionFormat, "getImproperInstance", BigFractionFormat, (Locale,), arg0)
end

function get_proper_instance(::Type{BigFractionFormat})
    return jcall(BigFractionFormat, "getProperInstance", BigFractionFormat, ())
end

function get_proper_instance(::Type{BigFractionFormat}, arg0::Locale)
    return jcall(BigFractionFormat, "getProperInstance", BigFractionFormat, (Locale,), arg0)
end

function parse(obj::BigFractionFormat, arg0::JString)
    return jcall(obj, "parse", BigFraction, (JString,), arg0)
end

function parse(obj::BigFractionFormat, arg0::JString, arg1::ParsePosition)
    return jcall(obj, "parse", BigFraction, (JString, ParsePosition), arg0, arg1)
end