text:
include("dependencies_for_runtests.jl")

using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVerticalDiffusivity

for closure in closures
    @eval begin
        using Oceananigans.TurbulenceClosures: $closure
    end
end

function closure_instantiation(closurename)
    closure = getproperty(TurbulenceClosures, closurename)()
    return true
end

function constant_isotropic_diffusivity_basic(T=Float64; ν=T(0.3), κ=T(0.7))
    closure = ScalarDiffusivity(T; κ=(T=κ, S=κ), ν=ν)
    return closure.ν == ν && closure.κ.T == κ
end

function anisotropic_diffusivity_convenience_kwarg(T=Float64; νh=T(0.3), κh=T(0.7))
    closure = HorizontalScalarDiffusivity(κ=(T=κh, S=κh), ν=νh)
    return closure.ν == νh && closure.κ.T == κh && closure.κ.T == κh
end

function run_constant_isotropic_diffusivity_fluxdiv_tests(FT=Float64; ν=FT(0.3), κ=FT(0.7))
          arch = CPU()
       closure = ScalarDiffusivity(FT, κ=(T=κ, S=κ), ν=ν)
          grid = RectilinearGrid(FT, size=(3, 1, 4), extent=(3, 1, 4))
    velocities = VelocityFields(grid)
       tracers = TracerFields((:T, :S), grid)
         clock = Clock(time=0.0)

    u, v, w = velocities
       T, S = tracers

    for k in 1:4
        interior(u)[:, 1, k] .= [0, -1/2, 0]
        interior(v)[:, 1, k] .= [0, -2,   0]
        interior(w)[:, 1, k] .= [0, -3,   0]
        interior(T)[:, 1, k] .= [0, -1,   0]
    end

    model_fields = merge(datatuple(velocities), datatuple(tracers))
    fill_halo_regions!(merge(velocities, tracers), arch, nothing, model_fields)

    U, C = datatuples(velocities, tracers)

    @test ∇_dot_qᶜ(2, 1, 3, grid, closure, C.T, Val(1), clock, nothing) == - 2κ
    @test ∂ⱼ_τ₁ⱼ(2, 1, 3, grid, closure, clock, U, nothing) == - 2ν
    @test ∂ⱼ_τ₂ⱼ(2, 1, 3, grid, closure, clock, U, nothing) == - 4ν
    @test ∂ⱼ_τ₃ⱼ(2, 1, 3, grid, closure, clock, U, nothing) == - 6ν

    return nothing
end

function anisotropic_diffusivity_fluxdiv(FT=Float64; νh=FT(0.3), κh=FT(0.7), νz=FT(0.1), κz=FT(0.5))
          arch = CPU()
      closureh = HorizontalScalarDiffusivity(FT, ν=νh, κ=(T=κh, S=κh))
      closurez = VerticalScalarDiffusivity(FT, ν=νz, κ=(T=κz, S=κz))
          grid = RectilinearGrid(arch, FT, size=(3, 1, 4), extent=(3, 1, 4))
           eos = LinearEquationOfState(FT)
      buoyancy = SeawaterBuoyancy(FT, gravitational_acceleration=1, equation_of_state=eos)
    velocities = VelocityFields(grid)
       tracers = TracerFields((:T, :S), grid)
         clock = Clock(time=0.0)

    u, v, w, T, S = merge(velocities, tracers)

    interior(u)[:, 1, 2] .= [0,  1, 0]
    interior(u)[:, 1, 3] .= [0, -1, 0]
    interior(u)[:, 1, 4] .= [0,  1, 0]

    interior(v)[:, 1, 2] .= [0,  1, 0]
    interior(v)[:, 1, 3] .= [0, -2, 0]
    interior(v)[:, 1, 4] .= [0,  1, 0]

    interior(w)[:, 1, 2] .= [0,  1, 0]
    interior(w)[:, 1, 3] .= [0, -3, 0]
    interior(w)[:, 1, 4] .= [0,  1, 0]

    interior(T)[:, 1, 2] .= [0,  1, 0]
    interior(T)[:, 1, 3] .= [0, -4, 0]
    interior(T)[:, 1, 4] .= [0,  1, 0]

    model_fields = merge(datatuple(velocities), datatuple(tracers))
    fill_halo_regions!(merge(velocities, tracers), arch, nothing, model_fields)

    U, C = datatuples(velocities, tracers)

    return (∇_dot_qᶜ(2, 1, 3, grid, closureh, C.T, Val(1), clock, nothing) == -  8κh &&
            ∇_dot_qᶜ(2, 1, 3, grid, closurez, C.T, Val(1), clock, nothing) == - 10κz &&
              ∂ⱼ_τ₁ⱼ(2, 1, 3, grid, closureh, clock, U, nothing) == - (2νh) &&
              ∂ⱼ_τ₁ⱼ(2, 1, 3, grid, closurez, clock, U, nothing) == - (4νz) &&
              ∂ⱼ_τ₂ⱼ(2, 1, 3, grid, closureh, clock, U, nothing) == - (4νh) &&
              ∂ⱼ_τ₂ⱼ(2, 1, 3, grid, closurez, clock, U, nothing) == - (6νz) &&
              ∂ⱼ_τ₃ⱼ(2, 1, 3, grid, closureh, clock, U, nothing) == - (6νh) &&
              ∂ⱼ_τ₃ⱼ(2, 1, 3, grid, closurez, clock, U, nothing) == - (8νz))
end

function time_step_with_variable_isotropic_diffusivity(arch)
    grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 2, 3))
    closure = ScalarDiffusivity(ν = (x, y, z, t) -> exp(z) * cos(x) * cos(y) * cos(t),
                                κ = (x, y, z, t) -> exp(z) * cos(x) * cos(y) * cos(t))

    model = NonhydrostaticModel(; grid, closure)
    time_step!(model, 1, euler=true)
    return true
end

function time_step_with_variable_anisotropic_diffusivity(arch)

    clov = VerticalScalarDiffusivity(ν = (x, y, z, t) -> exp(z) * cos(x) * cos(y) * cos(t),
                                     κ = (x, y, z, t) -> exp(z) * cos(x) * cos(y) * cos(t))
    cloh = HorizontalScalarDiffusivity(ν = (x, y, z, t) -> exp(z) * cos(x) * cos(y) * cos(t),
                                       κ = (x, y, z, t) -> exp(z) * cos(x) * cos(y) * cos(t))
    for clo in (clov, cloh)

        model = NonhydrostaticModel(grid=RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 2, 3)), closure=clo)

        time_step!(model, 1, euler=true)
    end
    return true
end

function time_step_with_tupled_closure(FT, arch)
    closure_tuple = (AnisotropicMinimumDissipation(FT), ScalarDiffusivity(FT))

    model = NonhydrostaticModel(closure=closure_tuple,
                                grid=RectilinearGrid(arch, FT, size=(1, 1, 1), extent=(1, 2, 3)))

    time_step!(model, 1, euler=true)

    return true
end

function run_time_step_with_catke_tests(arch, closure)
    grid = RectilinearGrid(arch, size=(1, 1, 1), extent=(1, 2, 3))
    buoyancy = BuoyancyTracer()

    # These shouldn't work (need :e in tracers)
    @test_throws ArgumentError HydrostaticFreeSurfaceModel(; grid, closure, buoyancy, tracers=:b)
    @test_throws ArgumentError HydrostaticFreeSurfaceModel(; grid, closure, buoyancy, tracers=(:b, :E))

    # CATKE isn't supported with NonhydrostaticModel (we don't diffuse vertical velocity)
    @test_throws ErrorException NonhydrostaticModel(; grid, closure, buoyancy, tracers=(:b, :e))

    model = HydrostaticFreeSurfaceModel(; grid, closure, buoyancy, tracers = (:b, :e))

    # Default boundary condition is Flux, Nothing... with CATKE this has to change.
    @test !(model.tracers.e.boundary_conditions.top.condition isa BoundaryCondition{Flux, Nothing})

    # Can we time-step?
    time_step!(model, 1)
    @test true

    # Once more for good measure
    time_step!(model, 1)
    @test true

    # Return model if we want to do more tests
    return model
end

function compute_closure_specific_diffusive_cfl(closurename)
    grid = RectilinearGrid(CPU(), size=(1, 1, 1), extent=(1, 2, 3))
    closure = getproperty(TurbulenceClosures, closurename)()

    model = NonhydrostaticModel(; grid, closure)
    dcfl = DiffusiveCFL(0.1)
    @test dcfl(model) isa Number

    tracerless_model = NonhydrostaticModel(; grid, closure, buoyancy=nothing, tracers=nothing)
    dcfl = DiffusiveCFL(0.2)
    @test dcfl(tracerless_model) isa Number

    return nothing
end

@testset "Turbulence closures" begin
    @info "Testing turbulence closures..."

    @testset "Closure instantiation" begin
        @info "  Testing closure instantiation..."
        for closure in closures
            @test closure_instantiation(closure)
        end
    end

    @testset "Constant isotropic diffusivity" begin
        @info "  Testing constant isotropic diffusivity..."
        for T in float_types
            @test constant_isotropic_diffusivity_basic(T)
            run_constant_isotropic_diffusivity_fluxdiv_tests(T)
        end
    end

    @testset "Constant anisotropic diffusivity" begin
        @info "  Testing constant anisotropic diffusivity..."
        for T in float_types
            @test anisotropic_diffusivity_convenience_kwarg(T)
            @test anisotropic_diffusivity_fluxdiv(T, νz=zero(T), νh=zero(T))
            @test anisotropic_diffusivity_fluxdiv(T)
        end
    end

    @testset "Time-stepping with variable diffusivities" begin
        @info "  Testing time-stepping with presribed variable diffusivities..."
        for arch in archs
            @test time_step_with_variable_isotropic_diffusivity(arch)
            @test time_step_with_variable_anisotropic_diffusivity(arch)
        end
    end

    @testset "Time-stepping with CATKE closure" begin
        @info "  Testing time-stepping with CATKE closure and closure tuples with CATKE..."
        for arch in archs
            warning = false

            @info "    Testing time-stepping CATKE by itself..."
            closure = CATKEVerticalDiffusivity(; warning)
            run_time_step_with_catke_tests(arch, closure)

            @info "    Testing time-stepping CATKE in a 2-tuple with ScalarDiffusivity..."
            closure = (CATKEVerticalDiffusivity(; warning), ScalarDiffusivity())
            model = run_time_step_with_catke_tests(arch, closure)
            @test first(model.closure) === closure[1]

            # Test that closure tuples with CATKE are correctly reordered
            @info "    Testing time-stepping CATKE in a 2-tuple with ScalarDiffusivity..."
            closure = (ScalarDiffusivity(), CATKEVerticalDiffusivity(; warning))
            model = run_time_step_with_catke_tests(arch, closure)
            @test first(model.closure) === closure[2]

            # These are slow to compile...
            @info "    Testing time-stepping CATKE in a 3-tuple..."
            closure = (ScalarDiffusivity(), CATKEVerticalDiffusivity(; warning), ScalarDiffusivity())
            model = run_time_step_with_catke_tests(arch, closure)
            @test first(model.closure) === closure[2]
        end
    end

    @testset "Closure tuples" begin
        @info "  Testing time-stepping with a tuple of closures..."
        for arch in archs
            for FT in float_types
                @test time_step_with_tupled_closure(FT, arch)
            end
        end
    end

    @testset "Diagnostics" begin
        @info "  Testing turbulence closure diagnostics..."
        for closure in closures
            compute_closure_specific_diffusive_cfl(closure)
        end
    end
end


text:
module Wikidata

    using HTTP, JSON

    struct WikidataEntity
        code::String
        dataDict::Dict
    end

    """
    expects an Wikidata Identifier for a specific item as a parameter, such as WikidataEntity('Q42') for Douglas Adams
    returns a WikidataEntity consisting of code::String and dataDict::Dict.
    # Example
    adams = Wikidata.WikidataEntity("Q42")
    println(Wikidata.label(adams))
    """
    function WikidataEntity(name::String)
        resp = HTTP.get("https://www.wikidata.org/wiki/Special:EntityData/$(name).json")
        #resp = HTTP.request("get", "www.google.com")
        str = String(resp.body)
        jobj = JSON.Parser.parse(str)["entities"][name]
        #descriptions = jobj["descriptions"]
        WikidataEntity(name, jobj)
    end

    """
    returns the english label (the title, the name) of the WikidataEntity
    # Example
    adams = Wikidata.WikidataEntity("Q42")
    println(Wikidata.label(adams))
    """
    function label(x::WikidataEntity)
        x.dataDict["labels"]["en"]["value"]
end
    """
    returns true, if the property for the given property identifier exists for this Entity
    # Example
    #load a Wikidata-Entity (Douglas Adams)
    adams = Wikidata.WikidataEntity("Q42")
    println(Wikidata.label(adams))

    #check if place of birth exists for this entity:
    if(Wikidata.hasproperty(adams, "P19"))
        placeofbirth = Wikidata.getproperty(adams, "P19")
        println(Wikidata.label(placeofbirth[1]))
    end
    """
    function hasproperty(x::WikidataEntity, property::String)
        return haskey(x.dataDict["claims"], property)
end
    """
        returns a list of properties of different type.
        Supported properties: item

        # Example
        #load a Wikidata-Entity (Douglas Adams)
        adams = Wikidata.WikidataEntity("Q42")
        println(Wikidata.label(adams))

        #check if place of birth exists for this entity:
        if(Wikidata.hasproperty(adams, "P19"))
            placeofbirth = Wikidata.getproperty(adams, "P19")
            println(Wikidata.label(placeofbirth[1]))
        end
    """
    function getproperty(x::WikidataEntity, property::String)
        if (!hasproperty(x, property))
            throw(ArgumentError("Entity does not have property $(property)"))
        end
        propertylist = x.dataDict["claims"][property]
        returnlist = Any[]
        for p in propertylist
            datatype = p["mainsnak"]["datavalue"]["type"]
            #item type
            if(datatype == "wikibase-entityid")
                itemcode = p["mainsnak"]["datavalue"]["value"]["id"]
                push!(returnlist, WikidataEntity(itemcode))
            #coordinates type
            elseif(datatype == "globecoordinate")
                latlontuple = (p["mainsnak"]["datavalue"]["value"]["latitude"], p["mainsnak"]["datavalue"]["value"]["longitude"])
                push!(returnlist, latlontuple)
            else
                throw("Unimplemented: Datatype $(datatype) is not implemented yet. Feel free to send a pull request to https://github.com/KimBue/Wikidata.jl" )

            end
        end
        return returnlist

end

end # module

ext: jl

max_stars_repo_name: KimBue/Wikidata


text:
p = (a=1, b="b")

# macro unp(x)
#     local nt = @eval $x
#     for k in keys(nt)
#         local v = nt[k]
#         @eval $k = $v
#     end
# end

function qt(x)
    return
end

function unp(x)
    exs = [Expr(:(=), k, x[k]) for k in keys(x)]
    eval.(exs)
end

function unq(x)
    exs = [Expr(:(=), k, QuoteNode(x[k])) for k in keys(x)]
    eval.(exs)
end


text:
function Decimal64(arg0::jdouble)
    return Decimal64((jdouble,), arg0)
end

function abs(obj::Decimal64)
    return jcall(obj, "abs", Decimal64, ())
end

function acos(obj::Decimal64)
    return jcall(obj, "acos", Decimal64, ())
end

function acosh(obj::Decimal64)
    return jcall(obj, "acosh", Decimal64, ())
end

function add(obj::Decimal64, arg0::Decimal64)
    return jcall(obj, "add", Decimal64, (Decimal64,), arg0)
end

function add(obj::Decimal64, arg0::jdouble)
    return jcall(obj, "add", Decimal64, (jdouble,), arg0)
end

function asin(obj::Decimal64)
    return jcall(obj, "asin", Decimal64, ())
end

function asinh(obj::Decimal64)
    return jcall(obj, "asinh", Decimal64, ())
end

function atan(obj::Decimal64)
    return jcall(obj, "atan", Decimal64, ())
end

function atan2(obj::Decimal64, arg0::Decimal64)
    return jcall(obj, "atan2", Decimal64, (Decimal64,), arg0)
end

function atanh(obj::Decimal64)
    return jcall(obj, "atanh", Decimal64, ())
end

function byte_value(obj::Decimal64)
    return jcall(obj, "byteValue", jbyte, ())
end

function cbrt(obj::Decimal64)
    return jcall(obj, "cbrt", Decimal64, ())
end

function ceil(obj::Decimal64)
    return jcall(obj, "ceil", Decimal64, ())
end

function compare_to(obj::Decimal64, arg0::Decimal64)
    return jcall(obj, "compareTo", jint, (Decimal64,), arg0)
end

function copy_sign(obj::Decimal64, arg0::Decimal64)
    return jcall(obj, "copySign", Decimal64, (Decimal64,), arg0)
end

function copy_sign(obj::Decimal64, arg0::jdouble)
    return jcall(obj, "copySign", Decimal64, (jdouble,), arg0)
end

function cos(obj::Decimal64)
    return jcall(obj, "cos", Decimal64, ())
end

function cosh(obj::Decimal64)
    return jcall(obj, "cosh", Decimal64, ())
end

function divide(obj::Decimal64, arg0::Decimal64)
    return jcall(obj, "divide", Decimal64, (Decimal64,), arg0)
end

function divide(obj::Decimal64, arg0::jdouble)
    return jcall(obj, "divide", Decimal64, (jdouble,), arg0)
end

function double_value(obj::Decimal64)
    return jcall(obj, "doubleValue", jdouble, ())
end

function equals(obj::Decimal64, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function exp(obj::Decimal64)
    return jcall(obj, "exp", Decimal64, ())
end

function expm1(obj::Decimal64)
    return jcall(obj, "expm1", Decimal64, ())
end

function float_value(obj::Decimal64)
    return jcall(obj, "floatValue", jfloat, ())
end

function floor(obj::Decimal64)
    return jcall(obj, "floor", Decimal64, ())
end

function get_exponent(obj::CalculusFieldElement)
    return jcall(obj, "getExponent", jint, ())
end

function get_field(obj::Decimal64)
    return jcall(obj, "getField", Field, ())
end

function get_real(obj::Decimal64)
    return jcall(obj, "getReal", jdouble, ())
end

function hash_code(obj::Decimal64)
    return jcall(obj, "hashCode", jint, ())
end

function hypot(obj::Decimal64, arg0::Decimal64)
    return jcall(obj, "hypot", Decimal64, (Decimal64,), arg0)
end

function int_value(obj::Decimal64)
    return jcall(obj, "intValue", jint, ())
end

function is_infinite(obj::Decimal64)
    return jcall(obj, "isInfinite", jboolean, ())
end

function is_na_n(obj::Decimal64)
    return jcall(obj, "isNaN", jboolean, ())
end

function linear_combination(obj::Decimal64, arg0::Decimal64, arg1::Decimal64, arg2::Decimal64, arg3::Decimal64)
    return jcall(obj, "linearCombination", Decimal64, (Decimal64, Decimal64, Decimal64, Decimal64), arg0, arg1, arg2, arg3)
end

function linear_combination(obj::Decimal64, arg0::Decimal64, arg1::Decimal64, arg2::Decimal64, arg3::Decimal64, arg4::Decimal64, arg5::Decimal64)
    return jcall(obj, "linearCombination", Decimal64, (Decimal64, Decimal64, Decimal64, Decimal64, Decimal64, Decimal64), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(obj::Decimal64, arg0::Decimal64, arg1::Decimal64, arg2::Decimal64, arg3::Decimal64, arg4::Decimal64, arg5::Decimal64, arg6::Decimal64, arg7::Decimal64)
    return jcall(obj, "linearCombination", Decimal64, (Decimal64, Decimal64, Decimal64, Decimal64, Decimal64, Decimal64, Decimal64, Decimal64), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function linear_combination(obj::Decimal64, arg0::Vector{Decimal64}, arg1::Vector{Decimal64})
    return jcall(obj, "linearCombination", Decimal64, (Vector{Decimal64}, Vector{Decimal64}), arg0, arg1)
end

function linear_combination(obj::Decimal64, arg0::Vector{jdouble}, arg1::Vector{Decimal64})
    return jcall(obj, "linearCombination", Decimal64, (Vector{jdouble}, Vector{Decimal64}), arg0, arg1)
end

function linear_combination(obj::Decimal64, arg0::jdouble, arg1::Decimal64, arg2::jdouble, arg3::Decimal64)
    return jcall(obj, "linearCombination", Decimal64, (jdouble, Decimal64, jdouble, Decimal64), arg0, arg1, arg2, arg3)
end

function linear_combination(obj::Decimal64, arg0::jdouble, arg1::Decimal64, arg2::jdouble, arg3::Decimal64, arg4::jdouble, arg5::Decimal64)
    return jcall(obj, "linearCombination", Decimal64, (jdouble, Decimal64, jdouble, Decimal64, jdouble, Decimal64), arg0, arg1, arg2, arg3, arg4, arg5)
end

function linear_combination(obj::Decimal64, arg0::jdouble, arg1::Decimal64, arg2::jdouble, arg3::Decimal64, arg4::jdouble, arg5::Decimal64, arg6::jdouble, arg7::Decimal64)
    return jcall(obj, "linearCombination", Decimal64, (jdouble, Decimal64, jdouble, Decimal64, jdouble, Decimal64, jdouble, Decimal64), arg0, arg1, arg2, arg3, arg4, arg5, arg6, arg7)
end

function log(obj::Decimal64)
    return jcall(obj, "log", Decimal64, ())
end

function log10(obj::Decimal64)
    return jcall(obj, "log10", Decimal64, ())
end

function log1p(obj::Decimal64)
    return jcall(obj, "log1p", Decimal64, ())
end

function long_value(obj::Decimal64)
    return jcall(obj, "longValue", jlong, ())
end

function multiply(obj::Decimal64, arg0::Decimal64)
    return jcall(obj, "multiply", Decimal64, (Decimal64,), arg0)
end

function multiply(obj::Decimal64, arg0::jdouble)
    return jcall(obj, "multiply", Decimal64, (jdouble,), arg0)
end

function multiply(obj::Decimal64, arg0::jint)
    return jcall(obj, "multiply", Decimal64, (jint,), arg0)
end

function negate(obj::Decimal64)
    return jcall(obj, "negate", Decimal64, ())
end

function new_instance(obj::Decimal64, arg0::jdouble)
    return jcall(obj, "newInstance", Decimal64, (jdouble,), arg0)
end

function pow(obj::Decimal64, arg0::Decimal64)
    return jcall(obj, "pow", Decimal64, (Decimal64,), arg0)
end

function pow(obj::Decimal64, arg0::jdouble)
    return jcall(obj, "pow", Decimal64, (jdouble,), arg0)
end

function pow(obj::Decimal64, arg0::jint)
    return jcall(obj, "pow", Decimal64, (jint,), arg0)
end

function reciprocal(obj::Decimal64)
    return jcall(obj, "reciprocal", Decimal64, ())
end

function remainder(obj::Decimal64, arg0::Decimal64)
    return jcall(obj, "remainder", Decimal64, (Decimal64,), arg0)
end

function remainder(obj::Decimal64, arg0::jdouble)
    return jcall(obj, "remainder", Decimal64, (jdouble,), arg0)
end

function rint(obj::Decimal64)
    return jcall(obj, "rint", Decimal64, ())
end

function root_n(obj::Decimal64, arg0::jint)
    return jcall(obj, "rootN", Decimal64, (jint,), arg0)
end

function round(obj::RealFieldElement)
    return jcall(obj, "round", jlong, ())
end

function scalb(obj::Decimal64, arg0::jint)
    return jcall(obj, "scalb", Decimal64, (jint,), arg0)
end

function short_value(obj::Decimal64)
    return jcall(obj, "shortValue", jshort, ())
end

function signum(obj::Decimal64)
    return jcall(obj, "signum", Decimal64, ())
end

function sin(obj::Decimal64)
    return jcall(obj, "sin", Decimal64, ())
end

function sin_cos(obj::CalculusFieldElement)
    return jcall(obj, "sinCos", FieldSinCos, ())
end

function sinh(obj::Decimal64)
    return jcall(obj, "sinh", Decimal64, ())
end

function sqrt(obj::Decimal64)
    return jcall(obj, "sqrt", Decimal64, ())
end

function subtract(obj::Decimal64, arg0::Decimal64)
    return jcall(obj, "subtract", Decimal64, (Decimal64,), arg0)
end

function subtract(obj::Decimal64, arg0::jdouble)
    return jcall(obj, "subtract", Decimal64, (jdouble,), arg0)
end

function tan(obj::Decimal64)
    return jcall(obj, "tan", Decimal64, ())
end

function tanh(obj::Decimal64)
    return jcall(obj, "tanh", Decimal64, ())
end

function to_degrees(obj::Decimal64)
    return jcall(obj, "toDegrees", Decimal64, ())
end

function to_radians(obj::Decimal64)
    return jcall(obj, "toRadians", Decimal64, ())
end

function to_string(obj::Decimal64)
    return jcall(obj, "toString", JString, ())
end


text:
import Knet: update!, train!

abstract type AbstractTrainable end
abstract type AbstractLayer <: AbstractTrainable end
# TODO: rename StaticLayer -> AbstractStaticLayer and StochasticLayer -> AbstractStochasticLayer
abstract type StaticLayer <: AbstractLayer end
abstract type StochasticLayer <: AbstractLayer end

"""
Initialise optimizer for each trainable parameters.
"""
function initoptim!(model::AbstractTrainable, otype; args...)
    for p in Knet.params(model)
        p.opt = otype(; args...)
    end
end

"""
Update all parameters in `ps` by back-propgating `losstape`.
"""
function update!(ps::Array, losstape::AutoGrad.Tape)
    for p in ps
        g = AutoGrad.grad(losstape, p)
        @assert p.opt != nothing "$p has no optimizer set up"
        update!(AutoGrad.value(p), g, p.opt)
    end
end

"""
Dynamic Adam - an Adam that can initialise first and second momentum during optimization.
This is used to deal dynamic models.
"""
mutable struct DynamicAdam
    adam::Knet.Adam
    isreset::Bool
end

DynamicAdam(isreset=true; kwargs...) = DynamicAdam(Knet.Adam(; kwargs...), isreset)

for T in (Array{Float32},Array{Float64},Knet.KnetArray{Float32},Knet.KnetArray{Float64}); @eval begin
    function update!(w::$T, g::$T, p::DynamicAdam)
        if !(p.adam.fstm === nothing) && length(g) > length(p.adam.fstm)
            F = eltype($T)
            if p.isreset
                p.adam.fstm = $T(zeros(F, size(g)...))
                p.adam.scndm = $T(zeros(F, size(g)...))
            else
                # Vector
                if length(size(g)) == 1
                    d_pad = length(g) - length(p.adam.fstm)
                    p.adam.fstm = vcat([p.adam.fstm, $T(zeros(F, d_pad))]...)
                    p.adam.scndm = vcat([p.adam.scndm, $T(zeros(F, d_pad))]...)
                # Matrix
                else
                    # Expand rows
                    if size(g, 1) > size(p.adam.fstm, 1)
                        d_pad = size(g, 1) - size(p.adam.fstm, 1)
                        p.adam.fstm = vcat([p.adam.fstm, $T(zeros(F, d_pad, size(g, 2)))]...)
                        p.adam.scndm = vcat([p.adam.scndm, $T(zeros(F, d_pad, size(g, 2)))]...)
                    # Expand cols
                    else
                        d_pad = size(g, 2) - size(p.adam.fstm, 2)
                        p.adam.fstm = hcat([p.adam.fstm, $T(zeros(F, size(g, 1), d_pad))]...)
                        p.adam.scndm = hcat([p.adam.scndm, $T(zeros(F, size(g, 1), d_pad))]...)
                    end
                end
            end
        end
        update!(w, g, p.adam)
    end
end; end

"""
Update all parameters of `model` by back-propgating `losstape`.
"""
update!(model::AbstractTrainable, losstape::AutoGrad.Tape) = update!(Knet.params(model), losstape)

"""
Compute the number of parameters.
"""
function numparams(model::AbstractTrainable)
    n = 0
    for p in Knet.params(model)
        n += prod(size(p))
    end
    return n
end

"""
Neural models.

The difference between a layer and a model is that for models, loss are defined.

All `NeuralModel` should be compatible with the following signature.
```julia
Base.show(nm::NeuralModel)
(nm::NeuralModel)(data, training::Val{true})    # loss
(nm::NeuralModel)(data, training::Val{false})   # evaluate metric
```
"""
abstract type NeuralModel <: AbstractTrainable end
Base.show(nm::NeuralModel) = error("Method not implemented")

"""
Train the model on a dataset.

For each batch in the dataset, do:
1. Compute the gradient of loss;
2. Update model parameters using the optimizers set for each parameters.

It returns batch averaged loss in the end.
"""
function train!(model::NeuralModel, dataloader; kwargs...)
    loss_list = []
    # Pop `:epoch` out from from `kwargs` if exist
    local epoch = nothing
    if :epoch in keys(kwargs)
        epoch = kwargs[:epoch]
        kwargs = Dict(kwargs...)
        # NOTE: the line below is uncommented because `:epoch` is useful for annealing
        #       so we need to deal with possible unwanted keyword agurments manually
        # pop!(kwargs, :epoch)
    end
    for (i, data_batch) in enumerate(dataloader)
        # Calculate the iteration till now if given the epoch
        if epoch != nothing
            kwargs[:iter] = length(dataloader) * (epoch - 1) + i
        end
        losstape = Knet.@diff model(data_batch, Val(:true); kwargs...)
        update!(model, losstape)
        push!(loss_list, Knet.value(losstape))
    end
    return mean(loss_list)
end

"""
Evaluate the model on a dataset.

It returns batch averaged metric value in the end.
"""
function evaluate(model::NeuralModel, dataloader; kargs...)
    loss_list = []
    for data_batch in dataloader
        push!(loss_list, model(data_batch, Val(:false); kargs...)[1])
    end
    return mean(loss_list)
end

export initoptim!, update!, DynamicAdam, numparams
export AbstractTrainable, AbstractLayer, StaticLayer, StochasticLayer
export NeuralModel, train!, evaluate

include("activations.jl")
export log1pexp, logexpm1, softplus, invsoftplus, leaky_relu
include("layers.jl")
export Dense, LazyDense, DynamicIn, DynamicOut, Chain
include("stochastic_layers.jl")
# TODO: somehow make below automatic
export GaussianDense, GaussianLazyDense, GaussianDynamicIn, GaussianDynamicOut
export GaussianLogVarDense, GaussianLogVarLazyDense, GaussianLogVarDynamicIn, GaussianLogVarDynamicOut
export KumaraswamyDense, KumaraswamyLazyDense, KumaraswamyDynamicIn, KumaraswamyDynamicOut
export BernoulliDense, BernoulliLazyDense, BernoulliDynamicIn, BernoulliDynamicOut
export BernoulliLogitDense, BernoulliLogitLazyDense, BernoulliLogitDynamicIn, BernoulliLogitDynamicOut
export GumbelBernoulliLogitDense, GumbelBernoulliLogitLazyDense, GumbelBernoulliLogitDynamicIn, GumbelBernoulliLogitDynamicOut
include("sbc.jl")
export MeanFieldSBC, StructuredSBC

include("rho.jl")
export Rho, updategrad!


text:
using Celeste
using Celeste.SDSSIO

@testset "main" begin

@testset "read_config" begin
    config, datasets = Celeste.read_config("configs/plain.yml")
    @test datasets["sdss"] == SDSSDataSet("data")

    config, datasets = Celeste.read_config("configs/notplain.yml")
    @test datasets["sdss"] == SDSSDataSet("data";
                                          dirlayout = :sdss,
                                          iostrategy = :masterrpc,
                                          compressed = true,
                                          slurp = true)
end

@testset "executable" begin
    # ensure data exists locally
    SampleData.fetch_sdss_data(3900, 6, 269, "all")

    Celeste.main(["sdss", "164.39", "164.41", "39.11", "39.13",
                  "--config", "configs/plain.yml"])
end

end


text:
# Julia program to implement Selection sort

function selection_sort(arr, n)
    for i in 1:n
        # Let the first element in this iteration be the minimum value
        min = i
        # Find the least element in the concerned part of the array
        for j in (i + 1):n
            if arr[min] > arr[j]
                min = j
            end
        end
        # Swap the least element with the beginning of the concerned part of the array
        arr[i], arr[min] = arr[min], arr[i]
    end
end

print("How many numbers do you want to sort? ")
n = readline()
n = parse(Int, n)
if (n <= 0)
    println("No numbers to sort!!!")
    exit()
end
arr = Int[]
print("Enter the numbers: ")
arr = [parse(Int, num) for num in split(readline())]
selection_sort(arr, n)
print("The numbers in sorted order is ")
print(arr)

"""
Time Complexity - O(n^2), where 'n' is the size of the array
Space Complexity - O(n)

SAMPLE INPUT AND OUTPUT

SAMPLE I
How many numbers do you want to sort? 5
Enter the numbers: 4 2 5 1 3
The numbers in sorted order is [1, 2, 3, 4, 5]

"""


text:
using TradingLogic
using Reactive

@testset "Working with signals"
  @testset "Change detection"
    s_inp = Reactive.Signal(5)
    s_chg = TradingLogic.schange(s_inp)
    push!(s_inp, 5)
    @test s_chg.value == false
    push!(s_inp, 8)
    @test s_chg.value == true
    push!(s_inp, 8)
    @test s_chg.value == false
  end
  @testset "Buffering signal history and SMA calculation"
    # simple vector test case
    nsma = 10
    vval = rand(100)
    vsma = fill(NaN, 100)
    [vsma[i] = mean(vval[(i-nsma+1):i]) for i = 10:100]

    s_inp = Reactive.Signal(vval[1])
    s_sma = Reactive.map(mean, Reactive.foldp(
                            TradingLogic.sighistbuffer!,
                            TradingLogic.initbuff(nsma, s_inp.value),
                            s_inp))
    vsma_sig = Array(Float64, 100)
    vsma_sig[1] = s_sma.value
    for i = 2:100
      push!(s_inp, vval[i])
      vsma_sig[i] = s_sma.value
    end
    @test all(isnan(vsma_sig[1:9])) == true
    @test vsma_sig[10:100] == roughly(vsma[10:100])
  end
end


text:
# array with no elements
a = []
# 1-dimensional array
b = [4, 5, 6]
println("b =",b)
println("first =",b[1])   # this is 4
println("last  =",b[end]) # this is 6
# 2-dimensional arrays
m  = [1 2 4; 3 4 5]
println("m=", m)
println("first =",m[1,1])     # this is 4
println("last  =",m[end,end]) # this is 5

max_stars_repo_name: elucian/julia

max_issues_repo_name: elucian/julia

max_forks_repo_name: elucian/julia


text:
#### visualize 2D.

# dot(u,x)+c = 0 to y = m*x + b.
function get2Dline(u::Vector{T}, c::T) where T

    m = -u[1]/u[2]
    b = c/u[2]

    return m, b
end

# traverse from the root towards the leaves.
# as we traverse, build up the boundary visualizations (t, y).
function getpartitionlines!(y_set::Vector{Vector{T}},
    t_set,
    node::BinaryNode{PartitionDataType{T}},
    level::Int,
    min_t, max_t, max_N_t::Int,
    centroid::Vector{T},
    max_dist::T) where T

    # draw line.
    m, b = get2Dline(node.data.hp.v, node.data.hp.c)
    t = LinRange(min_t, max_t, max_N_t)
    y = m .* t .+ b

    # prune according to distance from centroid.
    coordinates = collect( [t[n]; y[n]] for n = 1:length(t) )
    inds = findall(xx->norm(xx-centroid)<max_dist, coordinates)
    y = y[inds]
    t = t[inds]

    # prune according to constraints imposed by parents to current node.
    y_pruned, t_pruned = prunepartitionline(node, collect(y), collect(t))

    # store.
    push!(y_set, y_pruned)
    push!(t_set, t_pruned)
    # push!(y_set, y)
    # push!(t_set, t)

    # do not recurse at the level before leaf nodes, which is level 1.
    if level != 2

        # recurse.
        getpartitionlines!(y_set, t_set, node.left, level-1, min_t, max_t, max_N_t, centroid, max_dist)
        getpartitionlines!(y_set, t_set, node.right, level-1, min_t, max_t, max_N_t, centroid, max_dist)
    end

    return nothing
end

function prunepartitionline(node, y::Vector{T}, t::Vector{T}) where T

    #
    @assert length(y) == length(t)
    y_pruned = y
    t_pruned = t

    if isdefined(node, :parent)

        # node.
        c = node.parent.data.hp.c
        v = node.parent.data.hp.v

        # depending on whether current node was a left or right child of parent, use different constraints.
        constraint_func = xx->(dot(v,xx) < c)
        if node.parent.right == node

            constraint_func = xx->!(dot(v,xx) < c)
        end

        coordinates = collect( [t[n]; y[n]] for n = 1:length(t) )
        inds = findall(constraint_func, coordinates)
        y_pruned = y[inds]
        t_pruned = t[inds]

        y_pruned, t_pruned = prunepartitionline(node.parent, y_pruned, t_pruned)

    end

    return y_pruned, t_pruned
end

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
using StructJuMP, JuMP
using StructJuMPSolverInterface

include("select_solver.jl")
#an example model
scen = 2
m = StructuredModel(num_scenarios=scen)
@variable(m, x[1:2])
@NLobjective(m, Min, (x[1]+x[2])^2)
@NLconstraint(m, x[1] * x[2] == 10)

for i in 1:scen
    bl = StructuredModel(parent=m)
    @variable(bl, y)
    @NLconstraint(bl, x[2]^2 + x[1]*y ≤ 5)
    @NLobjective(bl, Min, (x[1]+x[2])*y)
end

structJuMPSolve(m)

getVarValue(m)


text:
#  ------------------------------------------------------------------------------------------
#  Copyright (c) Microsoft Corporation. All rights reserved.
#  Licensed under the MIT License (MIT). See LICENSE in the repo root for license information.
#  ------------------------------------------------------------------------------------------

try
    AzureClusterlessHPC == Main.TestCore.AzureClusterlessHPC
    print("Module already loaded.\n")
catch
    using AzureClusterlessHPC, PyCall, Test, SyntaxTree, Random
end

pool_id = "test_pool"
job_id = ["test_job_1", "test_job_2"]
task_id = Array{Dict, 1}(undef, 0)
push!(task_id, Dict("taskname" => "task_1", "pool" => 1)); push!(task_id, Dict("taskname" => "task_2", "pool" => 1))
num_tasks = 2
blobcontainer = "test_container"
output = [BlobFuture(blobcontainer, BlobRef("outfile_1")), BlobFuture(blobcontainer, BlobRef("outfile_2"))]
batch_client = [PyObject(nothing)]
blob_client = [PyObject(nothing)]

###################################################################################################
# Constructors

function test_vars(bctrl)
    @test typeof(bctrl) <: BatchController
    @test typeof(bctrl.pool_id) <: Union{String, Array}
    @test typeof(bctrl.job_id) <: Union{Array, String}
    @test typeof(bctrl.task_id) <: Array{Dict, 1}
    @test typeof(bctrl.num_tasks) <: Integer
    @test typeof(bctrl.output) <: Array{Any, 1}
    @test typeof(bctrl.blobcontainer) <: String
    @test typeof(bctrl.batch_client) <: Union{Array, Nothing}
    @test typeof(bctrl.blob_client) <: Union{Array, Nothing}
end

# Generic constructur
bctrl = BatchController(pool_id, job_id, task_id, num_tasks, output, blobcontainer, batch_client, blob_client)
test_vars(bctrl)

# Shortcut constructor
bctrl = BatchController(job_id, task_id, num_tasks, output)
test_vars(bctrl)

###################################################################################################
# Fetch

# Fetch for given task
output = [BlobFuture(blobcontainer, BlobRef("outfile_1")), BlobFuture(blobcontainer, BlobRef("outfile_2"))]
bctrl_empty = BatchController(pool_id, job_id, task_id, num_tasks, output, blobcontainer, [nothing], [nothing])

idx = 1
out = AzureClusterlessHPC.fetch(bctrl_empty, idx; destroy_blob=false, timeout=60, wait_for_completion=true)
@test typeof(out) == Tuple{}

# Fetch given task and overwrite input
output = [BlobFuture(blobcontainer, BlobRef("outfile_1")), BlobFuture(blobcontainer, BlobRef("outfile_2"))]
bctrl_empty = BatchController(pool_id, job_id, task_id, num_tasks, output, blobcontainer, [nothing], [nothing])
@test ~isnothing(output[1].blob)
@test ~isnothing(output[2].blob)

out = AzureClusterlessHPC.fetch!(bctrl_empty, idx; destroy_blob=false, timeout=60, wait_for_completion=true)
@test isnothing(out)
@test isnothing(output[1].blob)
@test ~isnothing(output[2].blob)

# Fetch all
output = [BlobFuture(blobcontainer, BlobRef("outfile_1")), BlobFuture(blobcontainer, BlobRef("outfile_2"))]
bctrl_empty = BatchController(pool_id, job_id, task_id, num_tasks, output, blobcontainer, [nothing], [nothing])
out = AzureClusterlessHPC.fetch(bctrl_empty; destroy_blob=false, timeout=60)
@test typeof(out) == Array{Any, 1}
@test length(out) == 2

# Fetch all tasks and overwrite input
output = [BlobFuture(blobcontainer, BlobRef("outfile_1")), BlobFuture(blobcontainer, BlobRef("outfile_2"))]
bctrl_empty = BatchController(pool_id, job_id, task_id, num_tasks, output, blobcontainer, [nothing], [nothing])
@test ~isnothing(output[1].blob)
@test ~isnothing(output[2].blob)

out = AzureClusterlessHPC.fetch!(bctrl_empty; destroy_blob=false, timeout=60)
@test isnothing(out[1])
@test isnothing(out[2])
@test isnothing(output[1].blob)
@test isnothing(output[2].blob)

# Fetch reduce
output = [BlobFuture(blobcontainer, BlobRef("outfile_1")), BlobFuture(blobcontainer, BlobRef("outfile_2"))]
bctrl_empty = BatchController(pool_id, job_id, task_id, num_tasks, output, blobcontainer, [nothing], [nothing])
AzureClusterlessHPC.fetchreduce(bctrl_empty)

# Fetch reduce overwrite input (non-tuple)
reduce = [9999]
output = [BlobFuture(blobcontainer, BlobRef("outfile_1")), BlobFuture(blobcontainer, BlobRef("outfile_2"))]
bctrl_empty = BatchController(pool_id, job_id, task_id, num_tasks, output, blobcontainer, [nothing], [nothing])
out = AzureClusterlessHPC.fetchreduce!(bctrl_empty, reduce)
@test out == reduce[]

# Fetch reduce overwrite input tuple
reduce = ([8888], [9999])
output = [BlobFuture(blobcontainer, BlobRef("outfile_1")), BlobFuture(blobcontainer, BlobRef("outfile_2"))]
bctrl_empty = BatchController(pool_id, job_id, task_id, num_tasks, output, blobcontainer, [nothing], [nothing])
out = AzureClusterlessHPC.fetchreduce!(bctrl_empty, reduce)
@test out == reduce


text:
using HorizonSideRobots

function n_size_chess_marking(r::Robot,n::Int)
    x, y = go_to_bottom_left_corner_return_coords_shift_of_start(r)

    all = show_path_after_moving(r, Nord)+1
    while isborder(r, Sud) == false
        move!(r, Sud)
    end

    side=Nord
    t=0
    while isborder(r, Ost) == false
        if n <= t < 2n
            for i in 1:n
                move!(r,side)
            end
        end

        direction_mark(r,side,n,all)

        t = t + 1
        if t > 2n
            t -= 2n
        end

        move!(r, Ost)
        side = HorizonSide(mod(Int(side) + 2, 4))
        show_path_after_moving(r, side)
        side = HorizonSide(mod(Int(side) + 2, 4))
    end
    if n <= t < 2n
        for i in 1:n
            move!(r, side)
        end
    end
    direction_mark(r, side, n, all)

    go_to_bottom_left_corner_return_coords_shift_of_start(r)
    move_n(r, Ost, x)
    move_n(r, Nord, y)
end

function direction_mark(r::Robot, side::HorizonSide, n::Int, all::Int)
    k = 0
    p = 0
    while isborder(r, side) == false
        p = 0
        for i in 1:n
            if isborder(r, side) == false
                putmarker!(r)
                move!(r, side)
                k = k + 1
                p += 1
            end
        end
        if (isborder(r, side) == true) && (p != n)
            putmarker!(r)
        end
        for j in 1:n
            if isborder(r, side) == false
                move!(r, side)
                k += 1
            end
        end
    end
    l = 2n + 1
    while all >= l
        if all == l
            move!(r, HorizonSide(mod(Int(side) + 2, 4)))
            if ismarker(r) == false
                move!(r, side)
                putmarker!(r)
            end
        end
        l= l + 3n
    end
    return k
end

function go_to_bottom_left_corner_return_coords_shift_of_start(r::Robot)
    x=0
    y=0
    while isborder(r, West) == false
        move!(r, West)
        x = x + 1
    end
    while isborder(r, Sud) == false
        move!(r, Sud)
        y = y + 1
    end
    return x, y
end

function show_path_after_moving(r::Robot, side::HorizonSide)
    n = 0
    while isborder(r, side) == false
        move!(r, side)
        n = n + 1
    end
    return n
end

function move_n(r::Robot, side::HorizonSide, n::Int)
    for i in 1:n
        move!(r, side)
    end
end


text:
# Copyright © 2021, UChicago Argonne, LLC
# All Rights Reserved.

# This source code is licensed under the BSD-style license found in the
# LICENSE file in the root directory of this source tree.

function get_fixed_nodes(ysol, fixlevel)
    if fixlevel < 0
        return nothing
    elseif fixlevel == 0
        return keys(ysol)
    else
        fixed_nodes = keys(ysol)
        for i in 1:fixlevel
            fixed_nodes = get_nonleaves(fixed_nodes)
        end
        return fixed_nodes
    end
end

function reset_variables(model, df_vars, bounds; unused_nodes=nothing, fixed_nodes=nothing)
    """

    """
    variables = JuMP.all_variables(model)

    # ##### (1) reset variables to original setting
    for row in eachrow(df_vars)
        var = variables[row.index]
        if JuMP.is_fixed(var)
            JuMP.unfix(var)
            if row.vartype in ["c"]
                JuMP.set_lower_bound(var, bounds["c_lb"])
                JuMP.set_lower_bound(var, bounds["c_ub"])
            end
        end
    end

    if isnothing(unused_nodes)
        unused_nodes = []
    end

    if isnothing(fixed_nodes)
        fixed_nodes = []
    end

    for row in eachrow(df_vars)
        var = variables[row.index]
        if row.vartype in ["y", "c"]
            if row.nodeid in unused_nodes
                JuMP.fix(var, 0.0, force=true)
            elseif row.nodeid in fixed_nodes
                JuMP.fix(var, row.sol, force=true)
            else
                if row.vartype in ["y"] && row.sol > 0.5
                    JuMP.set_start_value(var, 1.0)
                elseif row.vartype in ["y"] && row.sol <= 0.5
                    JuMP.set_start_value(var, 0.0)
                else
                    JuMP.set_start_value(var, row.sol)
                end
            end
        end
    end

    model
end

# function reset_variables_by_maxdist_fixlevel(model, df_vars, bounds, ysol, maxdist, fixlevel)
#     fixed_nodes = get_fixed_nodes(ysol, fixlevel)
#     nodes_ground = Set(1:maximum(df_vars.nodeid))
#     unused_nodes = get_unused_nodes(ysol, nodes_ground, maxdist)
#     reset_variables(model, df_vars, bounds; unused_nodes=nothing, fixed_nodes=fixed_nodes)
# end

function add_neighbor_constraints(model, df_vars, mindist, maxdist)
    variables = JuMP.all_variables(model)
    df_y_1 = filter(row -> row.vartype == "y" && row.sol > 0.5, df_vars)
    df_y_0 = filter(row -> row.vartype == "y" && row.sol <= 0.5, df_vars)
    active_nodes = unique(df_y_1.nodeid)
    @constraint(model, mindist <= sum(variables[row.index] for row in eachrow(df_y_0) if !(row.nodeid in active_nodes))
                + sum(1 - variables[i] for i in df_y_1.index) <= maxdist)

    return model
end

function set_restricted_MINLP(model, df_vars, bounds, ysol, mindist, maxdist, fixlevel)
    fixed_nodes = get_fixed_nodes(ysol, fixlevel)
    nodes_ground = Set(1:maximum(df_vars.nodeid))
    unused_nodes = get_unused_nodes(ysol, nodes_ground, maxdist)
        model = reset_variables(model, df_vars, bounds; unused_nodes=unused_nodes, fixed_nodes=fixed_nodes)
    model = add_neighbor_constraints(model, df_vars, mindist, maxdist)
    model
end

function initialize_STreCH_parameters(params)
    mindist = 0
    maxdist = params["INT_SOLU_DIST_INIT"]
    fixlevel = params["INT_FIXLEVEL_INIT"]
    nodelimit = params["INT_NODELIMIT_INIT"]
    mindist, maxdist, fixlevel, nodelimit
end

function update_STreCH_parameters(params, mindist, maxdist, fixlevel, nodelimit)
    if fixlevel < params["INT_FIXLEVEL_MAX"]
        fixlevel += 1
    else
        if maxdist < params["INT_SOLU_DIST_MAX"]
            fixlevel = params["INT_FIXLEVEL_INIT"]
            mindist = maxdist + 1
            maxdist = maxdist + params["INT_SOLU_DIST_INCR"]
        else
            if nodelimit < params["INT_NODELIMIT_MAX"]
                fixlevel = params["INT_FIXLEVEL_INIT"]
                mindist = 0
                maxdist = params["INT_SOLU_DIST_INIT"]
                nodelimit = 10 * nodelimit
            else
                return nothing, nothing, nothing, nothing
            end
        end
    end

    mindist, maxdist, fixlevel, nodelimit
end

function update_df_sol(df_sol, obs, df_vars, num_iter, time)
    ysol, csol = get_ysol_csol(df_vars)
    formula = get_formula(ysol, csol)[1]
    objval = compute_err_formula(obs, formula, "mse"; err2Inf=true)
    append!(df_sol, Dict("index"=>num_iter, "formula"=>formula, "objval"=>objval, "treesize"=>length(ysol), "time"=>time))
    df_sol
end

function STreCH(dir; datafile="data.txt", paramfile="params.txt")
    EPSILON = 1e-12
    params = read_parameters(dir * paramfile)
    df_sol = DataFrame(index=Int[], formula=String[], objval=Float64[], treesize=Int[], time=Float64[])

    io = open_logger(dir * "log.log", params["LOGLEVEL"])

    t_begin = time()
    num_iter = 1

    ### Parameters
    nodes = OrderedSet(1:(2^(params["INT_MAXDEPTH"]+1)-1))
    obs = DelimitedFiles.readdlm(dir * datafile)

    ### Solve an initial problem
    model, df_vars, bounds = get_MINLP_model(nodes, obs, params["OPERATORS"], params["FORMULATION"], params)
    model = reset_variables(model, df_vars, bounds, unused_nodes=8:1000)
    set_optimizer(model, SCIP.Optimizer)

    ### Set absgap, timelimit
    remaining_time = params["DBL_TIMELIMIT"] - (time() - t_begin)
    set_optimizer_attribute(model, "limits/time", remaining_time)
    set_optimizer_attribute(model, "display/verblevel", 2)

    ### Solve
    optimize!(model)

    ### Save the result
    objval_best, ysol_best, csol_best, df_vars_best = save_solution(dir, model, df_vars, obs, save_df_sol=false)
    df_sol = update_df_sol(df_sol, obs, df_vars_best, num_iter, time() - t_begin)
    CSV.write(dir * "df_sol.csv", df_sol)

    ### Improve
    mindist, maxdist, fixlevel, nodelimit = initialize_STreCH_parameters(params)
    while remaining_time > 0 && objval_best > EPSILON
        num_iter += 1
        @info "num_iter = $num_iter, fixlevel = $fixlevel, mindist=$mindist, maxdist=$maxdist, nodelimit=$nodelimit"
        model, df_vars, bounds = get_MINLP_model(nodes, obs, params["OPERATORS"], params["FORMULATION"], params)
        model = set_restricted_MINLP(model, df_vars_best, bounds, ysol_best, mindist, maxdist, fixlevel)

        ### Set the optimizer and its attributes
        set_optimizer(model, SCIP.Optimizer)
        remaining_time = params["DBL_TIMELIMIT"] - (time() - t_begin)
        absgap = objval_best < 1e-06 ? 0.0 : 0.99 * objval_best
        set_optimizer_attribute(model, "limits/time", remaining_time)
        set_optimizer_attribute(model, "limits/nodes", nodelimit)
        set_optimizer_attribute(model, "limits/absgap", absgap)
        set_optimizer_attribute(model, "display/verblevel", 2)

        ### Solve
        optimize!(model)
        termination_status = JuMP.termination_status(model)
        @info "num_iter = $num_iter , $termination_status"

        ### Read solution
        objval_new, ysol_new, csol_new, df_vars_new = save_solution(dir, model, df_vars, obs, save_df_sol=false)
        df_sol = update_df_sol(df_sol, obs, df_vars_new, num_iter, time() - t_begin)
        CSV.write(dir * "df_sol.csv", df_sol)

        ### Compare with the best solution and update paramters for the next iteration
        if objval_new < objval_best - EPSILON
            objval_best = objval_new
            ysol_best = ysol_new
            csol_best = csol_new
            df_vars_best = df_vars_new
            mindist, maxdist, fixlevel, nodelimit = initialize_STreCH_parameters(params)
        else
            """
            update parameters and keep the current solution
            """
            mindist, maxdist, fixlevel, nodelimit = update_STreCH_parameters(params, mindist, maxdist, fixlevel, nodelimit)
            if isnothing(mindist)
                break
            end
        end

        remaining_time = params["DBL_TIMELIMIT"] - (time() - t_begin)
    end

    # ### Remove files
    # Glob.glob(dir * "*.csv")
    # for file in Glob.glob(dir * "*.csv")
    #     if !occursin("df_sol.csv", file) && !occursin("df_vars.csv", file)
    #         rm(file)
    #     end
    # end

    close(io)
end

max_stars_repo_name: jongeunkim/STreCH

max_issues_repo_name: jongeunkim/STreCH

max_forks_repo_name: jongeunkim/STreCH


text:
export ang_to_au
export au_to_ang

export eV_to_au
export au_to_eV

export eV_per_ang_to_au
export au_to_eV_per_ang

export u_to_au
export au_to_u

export ps_inv_to_au
export au_to_ps_inv

export ps_to_au
export au_to_ps

export austrip_kwargs

auconvertstrip(u::Unitful.Units, x::Number) = ustrip(auconvert(u, x))

ang_to_au(x) = austrip(x*u"Å")
au_to_ang(x) = auconvertstrip(u"Å", x)

eV_to_au(x) = austrip(x*u"eV")
au_to_eV(x) = auconvertstrip(u"eV", x)

eV_per_ang_to_au(x) = austrip(x*u"eV/Å")
au_to_eV_per_ang(x) = auconvertstrip(u"eV/Å", x)

u_to_au(x) = austrip(x*u"u")
au_to_u(x) = auconvertstrip(u"u", x)

ps_inv_to_au(x) = austrip(x*u"ps^-1")
au_to_ps_inv(x) = auconvertstrip(u"ps^-1", x)

ps_to_au(x) = austrip(x*u"ps")
au_to_ps(x) = auconvertstrip(u"ps", x)

function austrip_kwargs(;kwargs...)
    stripped_kwargs = Dict{Symbol,Any}()
    for (key, val) in kwargs
        if val isa Unitful.Quantity
            stripped_kwargs[key] = austrip(val)
        else
            stripped_kwargs[key] = val
        end
    end
    stripped_kwargs
end


text:
using ArgCheck

using Test
using Random

include("checks.jl")
include("perf.jl")


text:
"""Definitions of various utility functions for metaprogramming"""
module Util

include("../docstrings.jl")

using Base.Meta
## AST manipulation utility functions

# useful shortcuts for nested macros
"""Add a dollar expression"""
dollar(v) = Expr(:$, v)
"Make a block expression from an array of exprs"
block(vs...) = Expr(:block, vs...)
"Add a & expression"
amp(v) = Expr(:&, v)

export dollar
export block
export amp

include("cleaning.jl")
export rmlines
export binarize!
export cleanast
export cleanblock
export interpolate_dollar
export remove_assertions

include("walks.jl")
export df_walk
export df_walk!
export bf_walk
export bf_walk!

include("fixpoint.jl")
export normalize
export normalize_nocycle

end


text:
#TODO : seed, copy
"""
{0,1,...,n-1}

Example usage:
discrete_obj.observation_space = Discrete(2)
"""
mutable struct Discrete <: AbstractSpace
    n::Int
    dtype::DataType

    Discrete(N::Int) = new(N, Int64)
end

sample(discrete_obj::Discrete) = rand(1:discrete_obj.n)

function contains(discrete_obj::Discrete, x::Union{Number, Array})
    as_int::Union{Number, Array, Nothing} = nothing
    try
        as_int = Int.(x)
    catch InexactError
        return false
    end
    return all(1 .<= as_int .<= discrete_obj.n)
end

Base.:(==)(discrete_obj::Discrete, other::Discrete) = discrete_obj.n == other.n

Base.length(discrete_obj::Discrete) = discrete_obj.n


text:
using Documenter, QuasiMonteCarlo

makedocs(;
    modules=[QuasiMonteCarlo],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/cscherrer/QuasiMonteCarlo.jl/blob/{commit}{path}#L{line}",
    sitename="QuasiMonteCarlo.jl",
    authors="Chad Scherrer",
    assets=String[],
)

deploydocs(;
    repo="github.com/cscherrer/QuasiMonteCarlo.jl",
)


text:
using Catalyst, LinearAlgebra, DiffEqJump, Test

@parameters t k[1:20]
@variables A(t) B(t) C(t) D(t)
rxs = [Reaction(k[1], nothing, [A]),            # 0 -> A
       Reaction(k[2], [B], nothing),            # B -> 0
       Reaction(k[3],[A],[C]),                  # A -> C
       Reaction(k[4], [C], [A,B]),              # C -> A + B
       Reaction(k[5], [C], [A], [1], [2]),      # C -> A + A
       Reaction(k[6], [A,B], [C]),              # A + B -> C
       Reaction(k[7], [B], [A], [2], [1]),      # 2B -> A
       Reaction(k[8], [A,B], [A,C]),            # A + B -> A + C
       Reaction(k[9], [A,B], [C,D]),            # A + B -> C + D
       Reaction(k[10], [A], [C,D], [2], [1,1]), # 2A -> C + D
       Reaction(k[11], [A], [A,B], [2], [1,1]), # 2A -> A + B
       Reaction(k[12], [A,B,C], [C,D], [1,3,4], [2, 3]),          # A+3B+4C -> 2C + 3D
       Reaction(k[13], [A,B], nothing, [3,1], nothing),           # 3A+B -> 0
       Reaction(k[14], nothing, [A], nothing, [2]),               # 0 -> 2A
       Reaction(k[15]*A/(2+A), [A], nothing; only_use_rate=true), # A -> 0 with custom rate
       Reaction(k[16], [A], [B]; only_use_rate=true),             # A -> B with custom rate.
       Reaction(k[17]*A*exp(B), [C], [D], [2], [1]),              # 2C -> D with non constant rate.
       Reaction(k[18]*B, nothing, [B], nothing, [2]),             # 0 -> 2B with non constant rate.
       Reaction(k[19]*t, [A], [B]),                                # A -> B with non constant rate.
       Reaction(k[20]*t*A, [B,C], [D],[2,1],[2])                  # 2A +B -> 2C with non constant rate.
  ]
@named rs = ReactionSystem(rxs,t,[A,B,C,D],k)
odesys = convert(ODESystem,rs)
sdesys = convert(SDESystem,rs)

# test show
io = IOBuffer()
show(io, rs)
str = String(take!(io))
@test count(isequal('\n'), str) < 30

# defaults test
def_p = [ki => float(i) for (i, ki) in enumerate(k)]
def_u0 = [A => 0.5, B => 1., C=> 1.5, D => 2.0]
defs = merge(Dict(def_p), Dict(def_u0))

@named rs = ReactionSystem(rxs,t,[A,B,C,D],k; defaults=defs)
odesys = convert(ODESystem,rs)
sdesys = convert(SDESystem,rs)
js = convert(JumpSystem,rs)
nlsys = convert(NonlinearSystem,rs)

@test ModelingToolkit.get_defaults(rs) ==
       ModelingToolkit.get_defaults(odesys) ==
       ModelingToolkit.get_defaults(sdesys) ==
       ModelingToolkit.get_defaults(js) ==
       ModelingToolkit.get_defaults(nlsys) ==
       defs

u0map = [A=>5.] # was 0.5
pmap = [k[1]=>5.] # was 1.
prob = ODEProblem(rs, u0map, (0,10.), pmap)
@test prob.p[1] == 5.
@test prob.u0[1] == 5.
u0 = [10., 11., 12., 13.]
ps = [float(x) for x in 100:119]
prob = ODEProblem(rs, u0, (0,10.), ps)
@test prob.p == ps
@test prob.u0 == u0

# hard coded ODE rhs
function oderhs(u,k,t)
       A = u[1]; B = u[2]; C = u[3]; D = u[4];
       du = zeros(eltype(u),4)
       du[1] = k[1] - k[3]*A + k[4]*C + 2*k[5]*C - k[6]*A*B + k[7]*B^2/2 - k[9]*A*B - k[10]*A^2 - k[11]*A^2/2 - k[12]*A*B^3*C^4/144 - 3*k[13]*A^3*B/6 + 2*k[14] - k[15]*A/(2+A) - k[16] - k[19]*t*A
       du[2] = -k[2]*B + k[4]*C - k[6]*A*B - k[7]*B^2 - k[8]*A*B - k[9]*A*B + k[11]*A^2/2 - 3*k[12]*A*B^3*C^4/144 - k[13]*A^3*B/6 + k[16] + 2*k[18]*B + k[19]*t*A - 2*k[20]*t*A*B^2*C
       du[3] = k[3]*A - k[4]*C - k[5]*C  + k[6]*A*B + k[8]*A*B + k[9]*A*B + k[10]*A^2/2 - 2*k[12]*A*B^3*C^4/144 - 2*k[17]*A*exp(B)*C^2/2 - k[20]*t*A*B^2*C
       du[4] = k[9]*A*B + k[10]*A^2/2 + 3*k[12]*A*B^3*C^4/144 + k[17]*A*exp(B)*C^2/2 + 2*k[20]*t*A*B^2*C
       du
end

# sde noise coefs
function sdenoise(u,k,t)
       A = u[1]; B = u[2]; C = u[3]; D = u[4];
       G = zeros(eltype(u),length(k),length(u))
       z = zero(eltype(u))

       G = [sqrt(k[1]) z  z z;
            z -sqrt(k[2]*B) z z;
            -sqrt(k[3]*A) z sqrt(k[3]*A) z;
            sqrt(k[4]*C) sqrt(k[4]*C) -sqrt(k[4]*C) z;
            2*sqrt(k[5]*C) z -sqrt(k[5]*C) z;
            -sqrt(k[6]*A*B) -sqrt(k[6]*A*B) sqrt(k[6]*A*B) z;
            sqrt(k[7]*B^2/2) -2*sqrt(k[7]*B^2/2) z z;
            z -sqrt(k[8]*A*B) sqrt(k[8]*A*B) z;
            -sqrt(k[9]*A*B) -sqrt(k[9]*A*B) sqrt(k[9]*A*B) sqrt(k[9]*A*B);
            -2*sqrt(k[10]*A^2/2) z sqrt(k[10]*A^2/2) sqrt(k[10]*A^2/2);
            -sqrt(k[11]*A^2/2) sqrt(k[11]*A^2/2) z z;
            -sqrt(k[12]*A*B^3*C^4/144) -3*sqrt(k[12]*A*B^3*C^4/144) -2*sqrt(k[12]*A*B^3*C^4/144) 3*sqrt(k[12]*A*B^3*C^4/144);
            -3*sqrt(k[13]*A^3*B/6) -sqrt(k[13]*A^3*B/6) z z;
            2*sqrt(k[14]) z z z;
            -sqrt(k[15]*A/(2+A)) z z z;
            -sqrt(k[16]) sqrt(k[16]) z z;
            z z -2*sqrt(k[17]*A*exp(B)*C^2/2) sqrt(k[17]*A*exp(B)*C^2/2);
            z 2*sqrt(k[18]*B) z z;
            -sqrt(k[19]*t*A) sqrt(k[19]*t*A) z z;
            z -2*sqrt(k[20]*t*A*B^2*C) -sqrt(k[20]*t*A*B^2*C) +2*sqrt(k[20]*t*A*B^2*C)]'
       return G
end

# test by evaluating drift and diffusion terms
p  = rand(length(k))
u  = rand(length(k))
t  = 0.
du = oderhs(u,p,t)
G  = sdenoise(u,p,t)
sdesys = convert(SDESystem,rs)
sf = SDEFunction{false}(sdesys, states(rs), parameters(rs))
du2 = sf.f(u,p,t)
@test norm(du-du2) < 100*eps()
G2 = sf.g(u,p,t)
@test norm(G-G2) < 100*eps()

# test conversion to NonlinearSystem
ns = convert(NonlinearSystem,rs)
fnl = eval(generate_function(ns)[2])
dunl = similar(du)
fnl(dunl,u,p)
@test norm(du-dunl) < 100*eps()

# tests the noise_scaling argument.
p  = rand(length(k)+1)
u  = rand(length(k))
t  = 0.
G  = p[21]*sdenoise(u,p,t)
@variables η
sdesys_noise_scaling = convert(SDESystem,rs;noise_scaling=η)
sf = SDEFunction{false}(sdesys_noise_scaling, states(rs), parameters(sdesys_noise_scaling))
G2 = sf.g(u,p,t)
@test norm(G-G2) < 100*eps()

# tests the noise_scaling vector argument.
p  = rand(length(k)+3)
u  = rand(length(k))
t  = 0.
G  = vcat(fill(p[21],8),fill(p[22],3),fill(p[23],9))' .* sdenoise(u,p,t)
@variables η[1:3]
sdesys_noise_scaling = convert(SDESystem,rs;noise_scaling=vcat(fill(η[1],8),fill(η[2],3),fill(η[3],9)))
sf = SDEFunction{false}(sdesys_noise_scaling, states(rs), parameters(sdesys_noise_scaling))
G2 = sf.g(u,p,t)
@test norm(G-G2) < 100*eps()

# tests using previous parameter for noise scaling
p  = rand(length(k))
u  = rand(length(k))
t  = 0.
G  = [p p p p]' .* sdenoise(u,p,t)
sdesys_noise_scaling = convert(SDESystem,rs;noise_scaling=k)
sf = SDEFunction{false}(sdesys_noise_scaling, states(rs), parameters(sdesys_noise_scaling))
G2 = sf.g(u,p,t)
@test norm(G-G2) < 100*eps()

# test with JumpSystem
js = convert(JumpSystem, rs)

midxs = 1:14
cidxs = 15:18
vidxs = 19:20
@test all(map(i -> typeof(equations(js)[i]) <: DiffEqJump.MassActionJump, midxs))
@test all(map(i -> typeof(equations(js)[i]) <: DiffEqJump.ConstantRateJump, cidxs))
@test all(map(i -> typeof(equations(js)[i]) <: DiffEqJump.VariableRateJump, vidxs))

pars = rand(length(k)); u0 = rand(1:10,4); ttt = rand();
jumps = Vector{Union{ConstantRateJump, MassActionJump, VariableRateJump}}(undef,length(rxs))

jumps[1] = MassActionJump(pars[1], Vector{Pair{Int,Int}}(), [1 => 1]);
jumps[2] = MassActionJump(pars[2], [2 => 1], [2 => -1]);
jumps[3] = MassActionJump(pars[3], [1 => 1], [1 => -1, 3 => 1]);
jumps[4] = MassActionJump(pars[4], [3 => 1], [1 => 1, 2 => 1, 3 => -1]);
jumps[5] = MassActionJump(pars[5], [3 => 1], [1 => 2, 3 => -1]);
jumps[6] = MassActionJump(pars[6], [1 => 1, 2 => 1], [1 => -1, 2 => -1, 3 => 1]);
jumps[7] = MassActionJump(pars[7], [2 => 2], [1 => 1, 2 => -2]);
jumps[8] = MassActionJump(pars[8], [1 => 1, 2 => 1], [2 => -1, 3 => 1]);
jumps[9] = MassActionJump(pars[9], [1 => 1, 2 => 1], [1 => -1, 2 => -1, 3 => 1, 4 => 1]);
jumps[10] = MassActionJump(pars[10], [1 => 2], [1 => -2, 3 => 1, 4 => 1]);
jumps[11] = MassActionJump(pars[11], [1 => 2], [1 => -1, 2 => 1]);
jumps[12] = MassActionJump(pars[12], [1 => 1, 2 => 3, 3 => 4], [1 => -1, 2 => -3, 3 => -2, 4 => 3]);
jumps[13] = MassActionJump(pars[13], [1 => 3, 2 => 1], [1 => -3, 2 => -1]);
jumps[14] = MassActionJump(pars[14], Vector{Pair{Int,Int}}(), [1 => 2]);

jumps[15] = ConstantRateJump((u,p,t) -> p[15]*u[1]/(2+u[1]), integrator -> (integrator.u[1] -= 1))
jumps[16] = ConstantRateJump((u,p,t) -> p[16], integrator -> (integrator.u[1] -= 1; integrator.u[2] += 1;))
jumps[17] = ConstantRateJump((u,p,t) -> p[17]*u[1]*exp(u[2])*binomial(u[3],2), integrator -> (integrator.u[3] -= 2; integrator.u[4] += 1))
jumps[18] = ConstantRateJump((u,p,t) -> p[18]*u[2], integrator -> (integrator.u[2] += 2))

jumps[19] = VariableRateJump((u,p,t) -> p[19]*u[1]*t, integrator -> (integrator.u[1] -= 1; integrator.u[2] += 1))
jumps[20] = VariableRateJump((u,p,t) -> p[20]*t*u[1]*binomial(u[2],2)*u[3], integrator -> (integrator.u[2] -= 2; integrator.u[3] -= 1; integrator.u[4] += 2))

statetoid = Dict(state => i for (i,state) in enumerate(states(js)))
jspmapper = ModelingToolkit.JumpSysMajParamMapper(js, pars)
symmaj = ModelingToolkit.assemble_maj(equations(js).x[1], statetoid, jspmapper)
maj    = MassActionJump(symmaj.param_mapper(pars), symmaj.reactant_stoch, symmaj.net_stoch, symmaj.param_mapper, scale_rates=false)
for i in midxs
  @test abs(jumps[i].scaled_rates - maj.scaled_rates[i]) < 100*eps()
  @test jumps[i].reactant_stoch == maj.reactant_stoch[i]
  @test jumps[i].net_stoch == maj.net_stoch[i]
end
for i in cidxs
  crj = ModelingToolkit.assemble_crj(js, equations(js)[i], statetoid)
  @test isapprox(crj.rate(u0,p,ttt), jumps[i].rate(u0,p,ttt))
  fake_integrator1 = (u=zeros(4),p=p,t=0); fake_integrator2 = deepcopy(fake_integrator1);
  crj.affect!(fake_integrator1);
  jumps[i].affect!(fake_integrator2);
  @test fake_integrator1 == fake_integrator2
end
for i in vidxs
  crj = ModelingToolkit.assemble_vrj(js, equations(js)[i], statetoid)
  @test isapprox(crj.rate(u0,p,ttt), jumps[i].rate(u0,p,ttt))
  fake_integrator1 = (u=zeros(4),p=p,t=0.); fake_integrator2 = deepcopy(fake_integrator1);
  crj.affect!(fake_integrator1); jumps[i].affect!(fake_integrator2);
  @test fake_integrator1 == fake_integrator2
end

# test for https://github.com/SciML/ModelingToolkit.jl/issues/436
@parameters t
@variables S(t) I(t)
rxs = [Reaction(1,[S],[I]), Reaction(1.1,[S],[I])]
@named rs = ReactionSystem(rxs, t, [S,I], [])
js = convert(JumpSystem, rs)
dprob = DiscreteProblem(js, [S => 1, I => 1], (0.0,10.0))
jprob = JumpProblem(js, dprob, Direct())
sol = solve(jprob, SSAStepper())

# test for https://github.com/SciML/ModelingToolkit.jl/issues/1042
jprob = JumpProblem(rs, dprob, Direct(), save_positions=(false,false))

@parameters k1 k2
@variables R(t)
rxs = [Reaction(k1*S, [S,I], [I], [2,3], [2]),
       Reaction(k2*R, [I], [R]) ]
@named rs = ReactionSystem(rxs, t, [S,I,R], [k1,k2])
@test isequal(oderatelaw(equations(rs)[1]), k1*S*S^2*I^3/(factorial(2)*factorial(3)))
@test_skip isequal(jumpratelaw(equations(eqs)[1]), k1*S*binomial(S,2)*binomial(I,3))
dep = Set()
ModelingToolkit.get_variables!(dep, rxs[2], Set(states(rs)))
dep2  = Set([R,I])
@test dep == dep2
dep = Set()
ModelingToolkit.modified_states!(dep, rxs[2], Set(states(rs)))
@test dep == Set([R,I])

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

@test isequal2(jumpratelaw(rxs[1]), k1*S*S*(S-1)*I*(I-1)*(I-2)/12)
@test isequal2(jumpratelaw(rxs[1]; combinatoric_ratelaw=false), k1*S*S*(S-1)*I*(I-1)*(I-2))
@test isequal2(oderatelaw(rxs[1]), k1*S*S^2*I^3/12)
@test isequal2(oderatelaw(rxs[1]; combinatoric_ratelaw=false), k1*S*S^2*I^3)

#test ODE scaling:
os = convert(ODESystem,rs)
@test isequal2(equations(os)[1].rhs, -2*k1*S*S^2*I^3/12)
os = convert(ODESystem,rs; combinatoric_ratelaws=false)
@test isequal2(equations(os)[1].rhs, -2*k1*S*S^2*I^3)

# test ConstantRateJump rate scaling
js = convert(JumpSystem,rs)
@test isequal2(equations(js)[1].rate, k1*S*S*(S-1)*I*(I-1)*(I-2)/12)
js = convert(JumpSystem,rs;combinatoric_ratelaws=false)
@test isequal2(equations(js)[1].rate, k1*S*S*(S-1)*I*(I-1)*(I-2))

# test MassActionJump rate scaling
rxs = [Reaction(k1, [S,I], [I], [2,3], [2]),
       Reaction(k2, [I], [R]) ]
@named rs = ReactionSystem(rxs, t, [S,I,R], [k1,k2])
js = convert(JumpSystem, rs)
@test isequal2(equations(js)[1].scaled_rates, k1/12)
js = convert(JumpSystem,rs; combinatoric_ratelaws=false)
@test isequal2(equations(js)[1].scaled_rates, k1)


text:
include("src/Modeling.jl")
include("src/Controls.jl")
using Plots
ENV["GKSwstype"]="nul" # Call this if using WSL

mutable struct Params
    mc::Float64
    mp::Float64
    g::Float64
    l::Float64
    function Params()
        return new(10., 1., 10., 2.)
    end
end

# Differential equation used in model of inverted pendulum
function xdot(x::Vector{Float64}, params::Params, force::Float64)
    return [
        x[2],
        (-params.mp * params.g * sin(x[3]) * cos(x[3]) + params.mp * params.l * sin(x[3]) * x[4] ^ 2 + force) /
        (params.mc + params.mp - params.mp * sin(x[3])^2),
        x[4],
        (-params.g * (params.mc + params.mp) * sin(x[3]) + params.mp * params.l * sin(x[3]) * cos(x[3]) * x[4] ^ 2 + force * cos(x[3])) /
        (-params.l * (params.mc + params.mp) + params.mp * params.l * cos(x[3]) ^ 2)
    ]
end

params = Params()

# Linearized State-Space model for making the controller
A = [
    0 1 0 0;
    0 0 (params.mp * params.g / params.mc) 0;
    0 0 0 1;
    0 0 (params.g * (params.mc + params.mp) / params.mc / params.l) 0
]

B = [0; 1 / params.mc; 0; -1 / params.mc / params.l]

C = [1 0 0 0]

ss = StateSpace(A, B, C)
controller = Controller(ss, [-1, -1.1, -1.2, -1.3])

# Setting up time constraints
dt = 0.01
t = 0:dt:20
x = [[0., 0., 0., 0.]] #Initial conditions of pendulum
u = [0.]

# Iterating through all the time steps
for i in t[2:end]
    push!(u, get_u(controller, last(x), 0.0))
    push!(x, rk4(xdot, last(x), dt, params, last(u)))
    if i == 1.0
        last(x)[4] = -2.0 # Instantaneous change in angular velocity, simulating flicking the ball
    end
end

# Making GIF

anim = @animate for i in 1:length(t)
    p1 = plot([x[i][1], x[i][1] + params.l * sin(x[i][3])],
        [0, params.l * cos(x[i][3])],
        leg = false,
        xlims = (-10, 1),
        ylims = (-1, 2.5),
        aspect_ratio = 1)
    scatter!(p1, [x[i][1]], [0], m = (:rect, 12))
    scatter!(p1, [x[i][1] + params.l * sin(x[i][3])], [params.l * cos(x[i][3])], m = (:circle, 8))
    p2 = plot(t, [j[1] for j in x], label = "Position", legend_position = :bottomright)
    scatter!(p2, [t[i]], [x[i][1]], label = "")
    p3 = plot(t, [j[3] for j in x], label = "Angle")
    scatter!(p3, [t[i]], [x[i][3]], label = "")
    plot(p1,p2, p3, layout = (3,1), size = (1200, 800))
end every 10;
gif(anim, "gifs/Inverted_pendulum.gif");


text:
module FastFibration

include("imports.jl")

include("graph_def.jl")
include("fiber_def.jl")
include("strongcomp_def.jl")

include("graph_prop.jl")
include("basic_ops.jl")
include("io_ops.jl")
include("topology.jl")
include("generation.jl")
include("fibration_functions.jl")

include("utils.jl")
end


text:
abstract type AbstractRequest end

struct RequestImmediate <: AbstractRequest end
Base.getindex(::RequestImmediate) = C.CEED_REQUEST_IMMEDIATE[]

struct RequestOrdered <: AbstractRequest end
Base.getindex(::RequestOrdered) = C.CEED_REQUEST_ORDERED[]

#=
# CeedRequest is not fully implemented in libCEED. When it is implemented, the
# following can be used as a starting point for the Julia interface.

struct Request <: AbstractRequest
    ref::RefValue{C.CeedRequest}
end

Request() = Request(Ref{C.CeedRequest}())

Base.wait(req::Request) = C.CeedRequestWait(req[])
=#

max_stars_repo_name: AdelekeBankole/libCEED

max_issues_repo_name: AdelekeBankole/libCEED

max_forks_repo_name: AdelekeBankole/libCEED


text:
matrixdict = Dict("hilb" => hilb, "hadamard" => hadamard,
                  "cauchy" => cauchy, "circul" => circul,
                  "dingdong" => dingdong, "frank" => frank,
                  "invhilb" => invhilb, "forsythe" => forsythe,
                  "magic" => magic, "grcar" => grcar,
                  "triw" => triw, "moler" => moler,
                  "pascal" => pascal, "kahan" => kahan,
                  "pei" => pei, "vand" => vand,
                  "invol" => invol, "chebspec" => chebspec,
                  "lotkin" => lotkin, "clement" => clement,
                  "fiedler" => fiedler, "minij" => minij,
                  "binomial" => binomialm, "tridiag" => tridiag,
                  "lehmer" => lehmer, "parter" => parter,
                  "chow" => chow, "randcorr" => randcorr,
                  "poisson" => poisson, "neumann" => neumann,
                  "rosser" => rosser, "sampling" => sampling,
                  "wilkinson" => wilkinson, "rando" => rando,
                  "randsvd" => randsvd, "rohess" => rohess,
                  "kms" => kms, "wathen" => wathen,
                  "oscillate" => oscillate, "toeplitz" => toeplitz,
                  "hankel" => hankel, "golub" => golub,
                  "companion" => companion,
                  "prolate" => prolate, "deriv2" => deriv2,
                  "shaw" => shaw, "wing" => wing,
                  "foxgood" => foxgood, "heat" => heat,
                  "baart" => baart, "phillips" => phillips,
                  "gravity" => gravity, "blur" => blur,
                  "spikes" => spikes, "ursell" => ursell,
                  "parallax" => parallax, "erdrey" => erdrey,
                  "gilbert" => gilbert, "smallworld" => smallworld
                  );

matrixclass = Dict("symmetric" => ["hilb", "cauchy", "circul", "dingdong",
                                   "invhilb", "moler", "pascal", "pei",
                                   "clement", "fiedler", "minij", "tridiag",
                                   "lehmer", "randcorr", "poisson", "wilkinson",
                                   "kms", "wathen", "oscillate", "prolate",
                                   "hankel"],

             "inverse" => ["hilb", "hadamard", "cauchy", "invhilb",
                           "forsythe", "magic", "triw", "moler", "pascal",
                           "kahan", "pei", "vand", "invol", "lotkin",
                           "clement", "fiedler", "minij", "tridiag",
                           "lehmer", "poisson", "kms" ],

             "ill-cond" => ["hilb", "cauchy", "frank", "invhilb",
                            "forsythe", "triw", "moler", "pascal",
                            "kahan","pei", "vand", "invol", "lotkin",
                            "tridiag", "rosser", "randsvd", "kms",
                            "oscillate", "prolate", "golub"],

             "pos-def" => ["hilb", "cauchy", "circul", "invhilb",
                           "moler", "pascal", "pei", "minij", "tridiag",
                           "lehmer", "poisson", "kms", "wathen", "oscillate"],

             "eigen" =>   ["hadamard", "circul", "dingdong", "frank",
                           "forsythe", "grcar", "pascal", "invol","chebspec",
                           "lotkin", "clement", "fiedler", "minij",
                           "tridiag", "parter", "chow", "poisson", "neumann",
                           "rosser", "sampling", "wilkinson","wathen",
                           "oscillate"],

             "sparse" => ["poisson", "neumann", "wathen", "blur", "erdrey", "gilbert",
                          "smallworld"],

             "random" => ["rosser", "rando", "randcorr", "randsvd", "rohess",
                          "wathen", "oscillate", "golub", "erdrey", "gilbert", "smallworld"],

             "regprob" => ["deriv2", "shaw", "wing", "foxgood", "heat",
                           "baart", "phillips", "gravity", "blur",
                           "spikes", "ursell", "parallax"],
              "graph" => ["erdrey", "gilbert", "smallworld"]
               );


text:
function writePhysicalNames(physicalNames::PhysicalNames)
    docCharArray::Array{Char, 1} = collect("\$PhysicalNames\n")
    docCharArray = vcat(docCharArray, collect(string(length(physicalNames.dimTags))), ['\n'])
    for i ∈ 1:length(physicalNames.dimTags)
        docCharArray = vcat(docCharArray, collect(string(physicalNames.dimTags[i][1])), [' '],
        collect(string(physicalNames.dimTags[i][2])), [' ','"'], collect(physicalNames.names[i]),
        ['"', '\n'])
    end
return vcat(docCharArray, collect("\$EndPhysicalNames\n"))
end

function writeNodes(nodes::Array{Float64, 2})
    @assert length(nodes[1, :]) == 3 "Nodes must have 3 columns or 3 dimensions."
    docCharArray::Array{Char, 1} = collect("\$Nodes\n")
    docCharArray = vcat(docCharArray, collect(string(length(nodes[:,1]))),['\n'])
    for nodeNo ∈ 1:length(nodes[:,1])
        docCharArray = vcat(docCharArray, [' '], collect(string(nodeNo)), [' '], collect(string(nodes[nodeNo,1])),
        [' '], collect(string(nodes[nodeNo,2])), [' '], collect(string(nodes[nodeNo,3])), ['\n'])
    end
    return vcat(docCharArray, collect("\$EndNodes\n"))
end

function writeElements(elementArray::Array{Element, 1})
    docCharArray::Array{Char, 1} = vcat(collect("\$Elements\n"), collect(string(length(elementArray))), ['\n'])
    for elementNo ∈ 1:length(elementArray)
        docCharArray = vcat(docCharArray, collect(string(elementNo)), [' '],
        collect(string(elementArray[elementNo].gmshElementType)), [' '],
        collect(string(length(elementArray[elementNo].attributes))), [' '])
        for attribNo ∈ 1:length(elementArray[elementNo].attributes)
            docCharArray = vcat(docCharArray, collect(string(elementArray[elementNo].attributes[attribNo])))
        end
        for nodeNo ∈ 1:length(elementArray[elementNo].nodeTags)
            docCharArray = vcat(docCharArray, [' '], collect(string(elementArray[elementNo].nodeTags[nodeNo])))
        end
        docCharArray = vcat(docCharArray, ['\n'])
    end
    return vcat(docCharArray, collect("\$EndElements\n"))
end
function writeGmshv2(fileName::String, physicalNames::PhysicalNames, nodes::Array{Float64, 2}, elementArray::Array{Element, 1})
    #Mesh File Header
    docCharArray::Array{Char, 1} = collect("\$MeshFormat\n2.2 0 8\n\$EndMeshFormat\n")
    #Mesh Physical Names
    docCharArray = vcat(docCharArray, writePhysicalNames(physicalNames),
    writeNodes(nodes), writeElements(elementArray))
    file = open(fileName, "w")
    write(file, String(docCharArray))
    close(file)
    return String(docCharArray)
end


text:
"""
    AST.record(::ASTMDP, sym::Symbol, val)

Recard an ASTMetric specified by `sym`.
"""
function record(mdp::ASTMDP, sym::Symbol, val)
    if mdp.params.debug
        push!(getproperty(mdp.metrics, sym), val)
    end
end

function record(mdp::ASTMDP; prob=1, logprob=exp(prob), miss_distance=Inf, reward=-Inf, event=false)
    AST.record(mdp, :prob, prob)
    AST.record(mdp, :logprob, logprob)
    AST.record(mdp, :miss_distance, miss_distance)
    AST.record(mdp, :reward, reward)
    AST.record(mdp, :event, event)
end


text:
#!/usr/bin/env julia

using VibronicToolkit

using ArgParse

s = ArgParseSettings()
s.autofix_names = true
@add_arg_table s begin
    "--conf"
        metavar = "FILE"
        help = "path to config file"
        required = true
    "--beta"
        metavar = "T"
        help = "reciprocal temperature"
        arg_type = Float64
        required = true
    "--basis-size"
        metavar = "N"
        help = "single-mode basis size"
        arg_type = Int
        required = true
    "--num-links"
        metavar = "P"
        help = "number of Trotter links"
        arg_type = Int
        required = true
end
c = parse_args(ARGS, s, as_symbols=true)

sys = read(c[:conf], DenseSystem)
beta = c[:beta]
basis_size = c[:basis_size]
P = c[:num_links]

trotter = Trotter(sys, beta, basis_size, P)

println("Z: $(trotter.Z)")
println("E: $(trotter.E)")
println("Cv: $(trotter.Cv)")


text:
using CSV, DataFrames, Distributed, Dates, LinearAlgebra, Distributions, DelimitedFiles, LambertW, SharedArrays

# Custom package
using Jevo

#=
In this script we are testing a variety of parameters to identify the range of
non-equilibrium ratios to which theoretical prediction are accurate.
=#

# Get date to append to output file
date = Dates.format(Dates.today(), "yyyy_mm_dd")

# Get number of workers as a script argument
if length(ARGS) == 1
    addprocs(parse(Int64, ARGS[1]))
elseif length(ARGS) > 1
    throw(ArgumentError("Only one command line argument (cores)."))
end

# Import packages needed for all workers
@everywhere  begin
    using Jevo
    using Distributions
    using DelimitedFiles
    using SharedArrays
    using LambertW
    using LinearAlgebra
end

# Parameters
@everywhere begin
    reps = 100
    steps = 10^7
    rho = [0, 0.1, 0.5, 1., 2, 3, 10]
    l_0 = 10
    N = 100
    nu = 1/N^2

    f0_eps = [(50, 2.5), (100, 2), (200, 1.5), (500, 1)]

    n = 20

    # Length scale
    l0(ϵ, n, f0=f0) = 1/ϵ * lambertw(2 * ϵ^2 * N * 10 * f0 * (n-1)/n^2* exp(5 * ϵ))
    l0(κ, l, ϵ, n, f0=f0) = 1/ϵ * lambertw(2 * ϵ^2 * N * l * f0 * (n-1)/n^2 * exp(5 * ϵ)/(1 + κ))

    # Find length cost to set length in equilibrium (derived from genetic load)
    #lambda(l_opt) = l0()^2/l_opt^2 * n^2 / (n-1) * 1/ϵ

    # Find length cost to set length in equilibrium (derived from genetic load modified by length dynamics)
    lambda_new(l_opt, ϵ, n, f0=f0) = l0(ϵ, n, f0)^3 / l_opt^2 * n^2 / (n-1) * 1/2

    #l_opt_static(κ, ϵ, n) = l0(ϵ, n) * sqrt(n^2 / ((n - 1) * ϵ * lambda(10)) * (1+κ))
    l_opt_dynamic(κ, ϵ, n, f0=f0) = l0(ϵ, n, f0)^(3/2) * sqrt(n^2 / ((n - 1) * 2 * lambda_new(10, ϵ, n, f0)) * (1+κ))
end

# Function to run one simulation
@everywhere function run(N, rho, nu, l_0, steps, f0, ϵ)
    # Initiate population
    pop = Jevo.mono_pop(N=100, l=l_0)
    Jevo.initiate!(pop)
    emat = ϵ * (ones(4, 4) - Matrix{Float64}(I, 4, 4))
    F = Jevo.fermi_fitness(f0=f0, fl=lambda_new(10, ϵ, n, f0)/(2N*l0(ϵ, n, f0)), epsilon=ϵ)
    for i in 1:steps
        Jevo.bp_substitution!(pop, emat, F)
        if rand() < rho/N
            Jevo.driver_mutation!(pop)
        end
        if rand() < nu
            Jevo.l_substitution!(pop, emat, F)
        end
        # Recover lost sites
        if length(pop.seqs) < 6
            Jevo.initiate!(pop)
        end
    end
    return Jevo.get_energy(pop, emat), length(pop.seqs)
end

for (_f0, ϵ) in f0_eps
    println("Running parameters f0=$_f0 and epsilon=$ϵ:")
    f0 = _f0 / 2N
    # Write Metadata
    open(date*"_sub_$(f0)_$(ϵ)_METADATA.txt", "a") do io
        write(io, "N=$N\n")
        write(io, "f0=$f0\n")
        write(io, "repetitions=$reps\n")
        write(io, "steps=$steps\n")
        write(io, "rho=$rho\n")
        write(io, "nu=$nu\n")
        write(io, "l_0=$l_0\n")
        write(io, "epsilon=$ϵ")
	write(io, "n=$n")
    end
    E_results = SharedArray{Float64, 2}(length(rho), reps)
    l_results = SharedArray{Float64, 2}(length(rho), reps)
    rho_list = SharedArray{Float64, 2}(length(rho), reps)
    # Run simulations and enjoy speed
    @sync @distributed for j in 1:reps
        for r in 1:length(rho)
            E, L = run(N, rho[r], nu, l_0, steps, f0, ϵ)
            E_results[r, j] = E
            l_results[r, j] = L
            rho_list[r, j] = rho[r]
        end
        println("Run $j done.")
    end
    df = DataFrame(gamma=[(E_results...)...], l=[(l_results...)...], rho=[(rho_list...)...])
    CSV.write(date * "_sub_$(f0)_$(ϵ)_results.csv", df)
end


text:
# Code written entirely by Eric Alcaide
# 8-puzzle example for SearchMethods.jl

using SearchMethods

# The goal is to order the 8-puzzle from an unordered state.
# 0 represents a blank space.
#  |-----------------|
#  |  0  |  1  |  2  |
#  |-----------------|
#  |  3  |  4  |  5  |
#  |-----------------|
#  |  6  |  7  |  8  |
#  |-----------------|

# Formalize the initial game state as an array
state = [3,1,2,4,0,5,6,7,8]

# Generate new nodes from a parent node.
function expand_node(node, tree)
	neighbors = []
	# Locate 0 -> determine moves
	loc = findin(node, 0)[1]
	# Check if moves are possible (up, down, left, right) + expand
	# expand up node
	if loc > 3
		up = copy(node)
		up[loc], up[loc-3] = up[loc-3], up[loc]
		push!(neighbors, up)
	end
	# expand down node
	if loc < 7
		down = copy(node)
		down[loc], down[loc+3] = down[loc+3], down[loc]
		push!(neighbors, down)
	end
	# expand left node
	if loc % 3 != 1
		left = copy(node)
		left[loc], left[loc-1] = left[loc-1], left[loc]
		push!(neighbors, left)
	end
	# expand right node
	if loc % 3 != 0
		right = copy(node)
		right[loc], right[loc+1] = right[loc+1], right[loc]
		push!(neighbors, right)
	end

	return neighbors
end

# Check if the puzzle is ordered
function goal_check(node)
	if node == [0,1,2,3,4,5,6,7,8] return true
	else return false
	end
end

# TESTS FOR UNINFORMED SEARCH ALGORITHMS
println("UNINFORMED METHODS")
# Sample test for the Iterative Deepening Search algorithm.
println("IDS Example")
@time println(SearchMethods.ids(state, false, goal_check, expand_node, 2))
# Sample test for the Breadth Frirst Search algorithm.
println("BFS Example")
@time println(SearchMethods.bfs(state, false, goal_check, expand_node))
# Sample test for the Depth Limited Search algorithm.
println("DLS Example")
@time println(SearchMethods.dls(state, false, goal_check, expand_node, 3))


text:
# A simple example to demonstrate the use of cross validation
#
# Here, we consider a simple model: using a mean vector to represent
# a set of samples. The goodness of the model is assessed in terms
# of the RMSE (root-mean-square-error) evaluated on the testing set
#

using MLBase

# functions

function compute_center(X::Matrix{Float64})
	c = vec(mean(X, 2))
	@printf("training on %d samples => (%.4f, %.4f)\n",
		size(X,2), c[1], c[2])
	return c
end

function compute_rmse(c::Vector{Float64}, X::Matrix{Float64})
	v = sqrt(mean(sum(abs2(X .- c),1)))
	@printf("RMSE on test set: %.6f\n\n", v)
	return v
end

# data

const n = 200
const data = [2., 3.] .+ randn(2, n)

# cross validation

(c, v, inds) = cross_validate(
	inds -> compute_center(data[:, inds]),     		# training function
	(c, inds) -> compute_rmse(c, data[:, inds]),    # evaluation function
	n,    			# total number of samples
	Kfold(n, 5), 	# cross validation plan: 5-fold
	Reverse)		# smaller score indicates better model

# display results

@printf("best model = (%.4f, %.4f), score = %.6f\n", c[1], c[2], v)


text:
using Plots,  LaTeXStrings

function plot_fitness(step=0.1, imagename="fitness-function.png")
    x = 0:step:1
    fitness(i) = 2/(i+1) - 1

    plot(x, fitness.(x), color = :orange, line = 3, guidefontsize=14, legend = false,
         xlabel = L"t_{et}(i)", ylabel = L"f(i)", title = "Função fitness")

    savefig(imagename)
end

function plot_peaks_surface(step=0.1, imagename="peaks-function.png")
    xs = -3:step:3
    ys = -3:step:3
    x_grid = [x for x = xs for y = ys]
    y_grid = [y for x = xs for y = ys]
    peaks(x, y) = 3*(1-x)^2 * exp(-x^2 -(y+1)^2) -10(x/5 - x^3 - y^5)*exp(-x^2-y^2)-1/3*exp(-(x+1)^2-y^2)

    plot(x_grid, y_grid, -1 * peaks.(x_grid, y_grid), guidefontsize=14, label = false,
        xlabel = L"w_i", ylabel = L"w_j", title = L"E(w_i, w_j)", cbar=false,
        camera=(75,45), st = :surface, color = :bluesreds)

    savefig(imagename)
end

max_stars_repo_name: Lucas-Muniz/MNEAT

max_issues_repo_name: Lucas-Muniz/MNEAT

max_forks_repo_name: Lucas-Muniz/MNEAT


text:
#############################
# Internal Structures
#############################

const MaybeTask = Union{Nothing, Task}

mutable struct LoggedQueuePool <: AbstractThreadPool
    inq  :: Channel{Task}
    outq :: Channel{Task}
    cnt  :: Threads.Atomic{Int}
    lck  :: ReentrantLock
    recs :: Vector{ThreadRecord}
    ind  :: Int
    t0   :: UInt64
    log  :: ThreadLog

    LoggedQueuePool(tids, handler::Function) = begin
        N = length(tids)
        pool = new(Channel{Task}(N), Channel{Task}(N), Threads.Atomic{Int}(0), ReentrantLock(), ThreadRecord[], 0, time_ns(), ThreadLog())
        Threads.@threads for tid in 1:Threads.nthreads()
            if tid in tids
                @async handler(pool)
            end
        end
        return pool
    end
end

function _default_log_handler(pool)
    t0 = time_ns()
    tid = Threads.threadid()
    for t in pool.inq
        lock(pool.lck)
        pool.ind += 1
        ind = pool.ind
        push!(pool.recs, (ind, tid, true, time_ns()))
        unlock(pool.lck)
        schedule(t)
        wait(t)
        lock(pool.lck)
        push!(pool.recs, (ind, tid, false, time_ns()))
        unlock(pool.lck)
        put!(pool.outq, t)
        Threads.atomic_sub!(pool.cnt, 1)
    end
end

#############################
# Constructors / Finalizer
#############################

"""
    LoggedQueuePool(init_thrd=1, nthrds=Threads.nthreads())

The main LoggedQueuePool object. Its API mimics that of a `Channel{Task}`, but each
submitted task is executed on a different thread.  If `allow_primary` is true,
the assigned thread might be the primary, which will interfere with future
thread management for the duration of any heavy-computational (blocking)
processes.  If it is false, all assigned threads will be off of the primary.
Each thread will only be allowed one Task at a time, but each thread will
backfill with the next queued Task immediately on completion of the previous,
without regard to how bust the other threads may be.

The API for the `LoggedQueuePool` is Identical to that for [`QueuePool`](@ref).
"""
function LoggedQueuePool(init_thrd::Integer=1, nthrds::Integer=Threads.nthreads())
    thrd0 = min(init_thrd, Threads.nthreads())
    thrd1 = min(thrd0+nthrds-1, Threads.nthreads())
    return LoggedQueuePool(thrd0:thrd1, _default_log_handler)
end

#############################
# LoggedQueuePool API
#############################

"""
    Base.put!(pool::LoggedQueuePool, t::Task)

Put the task `t` into the pool, blocking until the pool has
an available thread.
"""
function Base.put!(pool::LoggedQueuePool, t::Task)
    Threads.atomic_add!(pool.cnt, 1)
    put!(pool.inq, t)
end

"""
    Base.put!(pool::LoggedQueuePool, fn::Function, args...)
    Base.put!(fn::Function, pool::LoggedQueuePool, args...)

Creates a task that runs `fn(args...)` and adds it to the pool, blocking
until the pool has an available thread.
"""
Base.put!(pool::LoggedQueuePool, fn::Function, args...) = Base.put!(pool, Task(()->fn(args...)))
Base.put!(fn::Function, pool::LoggedQueuePool, args...) = Base.put!(pool, Task(()->fn(args...)))

"""
    Base.take!(pool::LoggedQueuePool) -> Task

Takes the next available completed task from the pool, blocking until a
task is available.
"""
Base.take!(pool::LoggedQueuePool) = fetch(take!(pool.outq))

"""
    Base.close(pool::LoggedQueuePool)

Shuts down the pool, closing the internal thread handlers.  It is safe
to issue this command after all Tasks have been submitted, regardless of
the Task completion status. If issued while the pool is still active, it
will `yield` until all tasks have been completed.
"""
function Base.close(pool::LoggedQueuePool)
    close(pool.inq)
    while pool.cnt[] > 0
        sleep(0.1)
        yield()
    end
    close(pool.outq)
    _recordstolog!(pool.log, pool.recs, pool. t0)
end

"""
    Base.iterate(pool::LoggedQueuePool[, state])

Iterates over the completed Tasks, grabbing the next one available
and ending when the pool has been `close`ed.
"""
Base.iterate(pool::LoggedQueuePool, state=nothing) = iterate(pool.outq, state)
Base.IteratorSize(::LoggedQueuePool) = Base.SizeUnknown()
Base.eltype(::LoggedQueuePool) = Task

"""
    ThreadPools.isactive(pool::LoggedQueuePool)

Returns `true` if there are queued Tasks anywhere in the pool, either
awaiting execution, executing, or waiting to be retrieved.
"""
isactive(pool::LoggedQueuePool) = isready(pool.inq) || isready(pool.outq) || pool.cnt[] > 0

#############################
# Result Iterator
#############################

"""
    results(pool::LoggedQueuePool) -> result iterator

Identical to [`results(pool::LoggedQueuePool)`](@ref).
"""
results(pool::LoggedQueuePool) = ResultIterator(pool)

#############################
# ThreadPool Interface
#############################

function tmap(pool::LoggedQueuePool, fn::Function, itr)
    N = length(itr)
    sizehint!(pool.recs, N)
    result = Vector{_detect_type(fn, itr)}(undef, N)
    _fn = (ind, x) -> (ind, fn(x))
    @async begin
        for (ind, item) in enumerate(itr)
            put!(pool, _fn, ind, item)
        end
    end
    for i in 1:N
        (ind, y) = fetch(take!(pool))
        @inbounds result[ind] = y
    end
    return result
end


text:
function get_auc(x::Vector{T}, y::Vector{T}) where {T}
  AUC=zero(T)
  for i in 2:length(x)
    dx = x[i] - x[i-1]
    dy = y[i] - y[i-1]
    AUC += ( (dx*y[i-1]) + (0.5*dx*dy) )
  end
  return AUC
end


text:
"""
    HistoryFunction(h, sol, integrator)

Wrap history function `h`, solution `sol`, and integrator `integrator` to create a common
interface for retrieving values at any time point with varying accuracy.

Before the initial time point of solution `sol` values are calculated by history function
`h`, for time points in the time span of `sol` interpolated values of `sol` are returned,
and after the final time point of `sol` an inter- or extrapolation of the current state
of integrator `integrator` is retrieved.
"""
struct HistoryFunction{F1,F2,F3<:ODEIntegrator} <: Function
    h::F1
    sol::F2
    integrator::F3
end

function (f::HistoryFunction)(p, t, ::Type{Val{deriv}}=Val{0}; idxs=nothing) where deriv
    integrator = f.integrator

    @inbounds if integrator.tdir * t < integrator.tdir * f.sol.t[1]
        if deriv == 0 && idxs === nothing
            return f.h(p, t)
        elseif idxs === nothing
            return f.h(p, t, Val{deriv})
        elseif deriv == 0
            return f.h(p, t; idxs = idxs)
        else
            return f.h(p, t, Val{deriv}; idxs = idxs)
        end
    elseif integrator.tdir * t <= integrator.tdir * f.sol.t[end] # Put equals back
        return f.sol.interp(t, idxs, Val{deriv}, f.integrator.p)
    end

    # set boolean that indicates that history function was evaluated at time point past the
    # final time point of the current solution
    # NOTE: does not interfere with usual use of isout since this integrator is only used
    # for inter- and extrapolation of future values and saving of the solution but does not
    # affect whether time steps are accepted
    integrator.isout = true

    # handle extrapolations at initial time point
    if integrator.t == integrator.sol.prob.tspan[1]
        return constant_extrapolant(t, integrator, idxs, Val{deriv})
    else
        return OrdinaryDiffEq.current_interpolant(t, integrator, idxs, Val{deriv})
    end
end

function (f::HistoryFunction)(val, p, t, ::Type{Val{deriv}}=Val{0}; idxs=nothing) where deriv
    integrator = f.integrator

    @inbounds if integrator.tdir * t < integrator.tdir * f.sol.t[1]
        if deriv == 0 && idxs === nothing
            return f.h(val, p, t)
        elseif idxs === nothing
            return f.h(val, p, t, Val{deriv})
        elseif deriv == 0
            return f.h(val, p, t; idxs = idxs)
        else
            return f.h(val, p, t, Val{deriv}; idxs = idxs)
        end
    elseif integrator.tdir * t <= integrator.tdir * f.sol.t[end] # Put equals back
        return f.sol.interp(val, t, idxs, Val{deriv}, f.integrator.p)
    end

    # set boolean that indicates that history function was evaluated at time point past the
    # final time point of the current solution
    # NOTE: does not interfere with usual use of isout since this integrator is only used
    # for inter- and extrapolation of future values and saving of the solution but does not
    # affect whether time steps are accepted
    integrator.isout = true

    # handle extrapolations at initial time point
    if integrator.t == integrator.sol.prob.tspan[1]
        return constant_extrapolant!(val, t, integrator, idxs, Val{deriv})
    else
        return OrdinaryDiffEq.current_interpolant!(val, t, f.integrator, idxs, Val{deriv})
    end
end


text:
module RandomWrapper

using JavaCall

const AbstractWell = @jimport org.hipparchus.random.AbstractWell
const BaseRandomGenerator = @jimport org.hipparchus.random.BaseRandomGenerator
const Class = @jimport java.lang.Class
const Collection = @jimport java.util.Collection
const CorrelatedRandomVectorGenerator = @jimport org.hipparchus.random.CorrelatedRandomVectorGenerator
const DoubleStream = @jimport java.util.stream.DoubleStream
const ForwardingRandomGenerator = @jimport org.hipparchus.random.ForwardingRandomGenerator
const GaussianRandomGenerator = @jimport org.hipparchus.random.GaussianRandomGenerator
const HaltonSequenceGenerator = @jimport org.hipparchus.random.HaltonSequenceGenerator
const ISAACRandom = @jimport org.hipparchus.random.ISAACRandom
const InputStream = @jimport java.io.InputStream
const IntRandomGenerator = @jimport org.hipparchus.random.IntRandomGenerator
const IntStream = @jimport java.util.stream.IntStream
const IntegerDistribution = @jimport org.hipparchus.distribution.IntegerDistribution
const JDKRandomGenerator = @jimport org.hipparchus.random.JDKRandomGenerator
const JString = @jimport java.lang.String
const LongStream = @jimport java.util.stream.LongStream
const MersenneTwister = @jimport org.hipparchus.random.MersenneTwister
const NormalizedRandomGenerator = @jimport org.hipparchus.random.NormalizedRandomGenerator
const Object = @jimport java.lang.Object
const Random = @jimport java.util.Random
const RandomAdaptor = @jimport org.hipparchus.random.RandomAdaptor
const RandomDataGenerator = @jimport org.hipparchus.random.RandomDataGenerator
const RandomGenerator = @jimport org.hipparchus.random.RandomGenerator
const RandomVectorGenerator = @jimport org.hipparchus.random.RandomVectorGenerator
const RealDistribution = @jimport org.hipparchus.distribution.RealDistribution
const RealMatrix = @jimport org.hipparchus.linear.RealMatrix
const SobolSequenceGenerator = @jimport org.hipparchus.random.SobolSequenceGenerator
const StableRandomGenerator = @jimport org.hipparchus.random.StableRandomGenerator
const SynchronizedRandomGenerator = @jimport org.hipparchus.random.SynchronizedRandomGenerator
const UncorrelatedRandomVectorGenerator = @jimport org.hipparchus.random.UncorrelatedRandomVectorGenerator
const UniformRandomGenerator = @jimport org.hipparchus.random.UniformRandomGenerator
const UnitSphereRandomVectorGenerator = @jimport org.hipparchus.random.UnitSphereRandomVectorGenerator
const Well1024a = @jimport org.hipparchus.random.Well1024a
const Well19937a = @jimport org.hipparchus.random.Well19937a
const Well19937c = @jimport org.hipparchus.random.Well19937c
const Well44497a = @jimport org.hipparchus.random.Well44497a
const Well44497b = @jimport org.hipparchus.random.Well44497b
const Well512a = @jimport org.hipparchus.random.Well512a

export AbstractWell
export BaseRandomGenerator
export CorrelatedRandomVectorGenerator
export ForwardingRandomGenerator
export GaussianRandomGenerator
export HaltonSequenceGenerator
export ISAACRandom
export IntRandomGenerator
export JDKRandomGenerator
export MersenneTwister
export NormalizedRandomGenerator
export RandomAdaptor
export RandomDataGenerator
export RandomGenerator
export RandomVectorGenerator
export SobolSequenceGenerator
export StableRandomGenerator
export SynchronizedRandomGenerator
export UncorrelatedRandomVectorGenerator
export UniformRandomGenerator
export UnitSphereRandomVectorGenerator
export Well1024a
export Well19937a
export Well19937c
export Well44497a
export Well44497b
export Well512a
export doubles
export equals
export get_class
export get_generator
export get_next_index
export get_rank
export get_root_matrix
export hash_code
export ints
export longs
export next_beta
export next_boolean
export next_bytes
export next_deviate
export next_deviates
export next_double
export next_exponential
export next_float
export next_gamma
export next_gaussian
export next_hex_string
export next_int
export next_log_normal
export next_long
export next_normal
export next_normalized_double
export next_permutation
export next_poisson
export next_sample
export next_sample_with_replacement
export next_uniform
export next_vector
export next_zipf
export notify
export notify_all
export of
export set_seed
export skip_to
export to_string
export wait

include("abstract_well.jl")
include("base_random_generator.jl")
include("correlated_random_vector_generator.jl")
include("forwarding_random_generator.jl")
include("gaussian_random_generator.jl")
include("halton_sequence_generator.jl")
include("int_random_generator.jl")
include("isaac_random.jl")
include("jdk_random_generator.jl")
include("mersenne_twister.jl")
include("normalized_random_generator.jl")
include("random_adaptor.jl")
include("random_data_generator.jl")
include("random_generator.jl")
include("random_vector_generator.jl")
include("sobol_sequence_generator.jl")
include("stable_random_generator.jl")
include("synchronized_random_generator.jl")
include("uncorrelated_random_vector_generator.jl")
include("uniform_random_generator.jl")
include("unit_sphere_random_vector_generator.jl")
include("well1024a.jl")
include("well19937a.jl")
include("well19937c.jl")
include("well44497a.jl")
include("well44497b.jl")
include("well512a.jl")

end